EN (546767), страница 2
Текст из файла (страница 2)
A matrixrepresentation is more convenient for the analytic transformation and implementation of solutions forthe computer. The solution of the homogeneous equationtt rrtrttd P(τ)− B τ0(23)L(τ0 ) = e L(0) ≡ P(τ0 ) L(0) := − BP(τ) ,dτconnects the radiation at the slab lower boundary with an expression on the top and is called the propropagator [7].4Eurotherm Conference No.
95: Computational Thermal Radiation in Participating Media IVIOP PublishingJournal of Physics: Conference Series 369 (2012) 012021doi:10.1088/1742-6596/369/1/012021The problem of the solution (22) is connected with negative and positive exponents in theexpression that leads to fast worsening of the system matrix conditions. In order to eliminate thiseffect, we used the scale transformation [4].The matrixexponent is represented in the form [4]tt rr ΓτB τ0−10(24)e = Ue U ,rttt twhere U is the matrix of eigenvectors of the matrix B and Γ = diag(Γ − , Γ + ) is the diagonal matrixttof eigenvalues sorted by ascending ordering so that Γ − = −Γ + .Consequently, the equation (22) can be rewritten asttrrrttu12 L + (0) eΓ− τ0 ut 11 eΓ− τ0 ut 12 L + (τ0 ) J − u11tr− −Γt τ t(25)r = r ,+ tt−Γ + τ0 t + 0L(τ)L(0)eueuuu−0− J+ 2122 2122rttt tr J− t τ0 Γτt −1 rt −1 u11 u12 −1where J ≡ r = S ∫ e U M ∆ (τ, µ 0 )d τ , U ≡ tt . u 21 u 22 0 J+ Note that expression (25) contains the exponent only with the positive power.
Expressing inrrequation (25) the emerging radiation L + (0), L − (τ0 ) from a layer through the incident radiationrrL − (0), L + (τ0 ) we find a smooth regular part of the solution in the form of the so-called scatterers [7]:rrrtt L − (0) F− R − T− L + (0) t r(26)r = r + t, L + (τ0 ) F+ T+ R + L − (τ0 ) ttrrtttttt −1 F− t J − R − T− t u11− eΓ− τ0 u12 t − u12eΓ− τ0 u11 t = H t twhere r = H r , t , H ≡ −Γt τ t .tt−Γ τ− u 22 u 21 e + 0 u 21 − e + 0 u 22 F+ J + T+ R + The obtained system (26) is the desired solution that makes definition of the radiance distributionof the reflected and transmitted radiations possible.
Note that the expression (26) is the rigorousanalytical solution of the boundary value problem for VRTE discretized by DOM. The transfer fromVRTE (1) to the equation system (26) is possible due to the singularity elimination (2) and to thechange of the scattering integral (discretization) by the Gaussian quadrature (18).The solution obtained in the form of scatterers (26) has a functional type and allows entering aphotometric concept of the layer radiance factor by the reflection or transmission. This immediatelygives the possibility to formulate a matrix-operator method to replace the two adjacent layers by theequivalent one described by the expression identical to (32), but with effective parameters [3].Consequently, the solution in the form of scatterers possesses the property of invariance.Equation (26) can be reorganized in a form similar to the propagator (23):t t t tt trr T+ − R − T−−1R + R − T−−1 rt −1 tt −1 L(0) + F .(27)L(τ0 ) = T− − T− R +This kind of transformation has been called the “star product” [7].
Since one knows the differentialequation (23) for the propagator, one is able to get the one-point problem with initial conditions of thematrix Riccati equation for the matrix elements of the scatterers [7], e.g. the reflection:ttttt t tt t t t tt t t t tt b1b2 d R− t t td R+= b 2 + b1R − + R − b1 + R − b 2 R − ;= − b 2 + b1R + − R + b1 − R + b 2 R + ; B ≡ tt , (28)dτdτ − b 2 − b1 tand similarly for the transmission T± .It is easy to see that (28) is an equivalent to the equations of Ambartsumyan-Chandrasekhar [8],derived from the principle of invariant embedding.
Thus, after eliminating the solution anisotropic partand the equation discretization we reach a unique analytical solution in the matrix form (26). In thissense the method of successive orders of scattering is an iterative method for finding the inversetmatrix A in (26), and the Monte Carlo is a stochastic method.5Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IVIOP PublishingJournal of Physics: Conference Series 369 (2012) 012021doi:10.1088/1742-6596/369/1/012021Implementing the suggested solution (26) as the algorithm includes the calculation of the followingcomponents: zeroes and weights of the quadrature formula for the VRTE discretization (18), sourcefunction (4), eigenvectors and values of the system matrix (24), and products of matrices in (26).The practical implementation of the indicated algorithm depends mainly on the size of matrices in(26). Note that (26) is obtained for the regular part, therefore, its size is determined mostly not by themedium optical parameters and the boundary conditions, but by the elimination of the solutionanisotropic part.
In general, M≈N≈K. However, in case of correct elimination of the anisotropic partsolution, i.e. the smooth regular part is near to the isotropic angular distribution, it is possible thatM<<N<<K. This improves considerably the algorithm performance.4. Solution anisotropic part eliminationNow let’s analyse and compare all the main known methods of the anisotropic part elimination fromthe point of view of matrix size M and N reduction and the enhancement of calculation efficiency ofthe solution (26).4.1. Elimination of direct non-scattering radiationFor the first time Eddington proposed the elimination of the anisotropic part by subtracting the directnon-scattered radiation on the basis of the Bouguer law, Milne developed the theory and it was fullycrystallized in the works of Chandrasekhar [8].
Thus, the expression for the direct non-scatteredradiation in the plane of reference coincided with the plane of incidence is given by:rrtrrTL a (τ, µ, ϕ) ≡ L 0 (τ, µ, ϕ) = e −τ / µ0 R(ϕ) L0 δ(ˆl − ˆl 0 ), L0 = E0 [1 p cos ϕ0 p sin ϕ0 q ] ,(29)r0where L is the vector of Stokes parameters of incident radiation, E0 is the normal to beam irradiance,p, q are degrees of linear and circular polarization of the beam, and ϕ0 is the azimuth of polarization.In this case the source function in the equation for the regular part takes the formrΛ −τ / µ0 tt∆ (τ, ˆl ) =(30)eR(ϕ) x (ˆl , ˆl0 ) .4πUnder these conditions the radiation scattered into small angles slightly differs from the directradiation, and this method becomes inefficient.
Strongly anisotropic scattering leads to a significantincrease in the number of terms in the expansion of the scattering matrix in a series of generalizedspherical functions K, which correspondingly increases the size of matrices N and M in the solution(26).4.2. Delta-M methodSome works were done to solve the anisotropy problem. Various algorithms based on the truncation ofscattering phase function were used. Delta-M method is the most efficient among them; it could begeneralized in the vector case [9]. As in the scalar one the scattering phase function is represented as asum of a delta-function and a smooth parttttx (ˆl , ˆl′) = 4πf 1δ(ˆl − ˆl′) + (1 − f ) x * (ˆl , ˆl′)(31)tHere f is a part of anisotropy scattering; 1 is an identity matrix.Equation (31) for Greek matrix (11) can be written down astttt tt(32)χk = f 1 + (1 − f )χ*k , χ*k = (χ k − f 1) (1 − f ) .This phase matrix transformation leads to the scale transformation of optical depth and singlescattering albedo(33)τ* = (1 − Λf )τ, Λ* = Λ (1 − f ) (1 − Λf ) .Thus, delta-M method could decrease K significantly.
However such an approach distorts initialboundary problem and leads to the error in a small sighting angle (intrinsically we do not consider thecoarsest fraction of aerosol) or to oscillations in the angular radiance distribution.6Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IVIOP PublishingJournal of Physics: Conference Series 369 (2012) 012021doi:10.1088/1742-6596/369/1/0120214.3. TMS and IMS-methodsTo eliminate the problems in VRTE solution mentioned above, Nakajima and Tanaka [10] returned tothe idea of determination of the anisotropic part in the solution on the base of approximate analyticalpresentation of angular distribution of Stokes parameters for the first and second order of scattering:rrrrr(34)L(τ, ˆl ) = L 0 (τ, ˆl ) + L1 (τ, ˆl ) + L 2 (τ, ˆl ) + L r (τ, ˆl ) .14444244443rLa ( τ, ˆl )Here the first two orders satisfy the equations:rr t∂L1 (τ, ˆl ) rΛ −τ µ0 t(35)µ+ L1 (τ, ˆl ) =eR(ϕ) L0 x (ˆl 0 , ˆl ) ,∂τ4πrr∂L 2 (τ, ˆl ) rΛ t tˆˆ tˆ′ ˆ′′′µ+ L 2 (τ, ˆl ) =χχ(36)R()x(l,l)R()L1 ( τ, l ) d l .∂τ4π ∫Equations (35) and (36) are solved analytically.
However the solution of equation (36) is a tripleintegral, and the run-time of its calculations can exceed the run-time of the initial problem solution(26). For this reason an approximation describing Stokes parameters distribution near the incidentdirection is generally used. The most efficient way is to use the small angle approximation, where thedispersion of scattered and non-scattered rays is not taken into account.
This approximation isequivalent to the change of µ in (35), (36) by µ0. The most known codes using both delta-M and TMSmethods are DISORT [11] in the scalar case and Pstar [12] in the vector one.4.4. Small angle modification of spherical harmonic methodrIt could be shown that the natural evolution of TMS-method is the inclusion in L a (τ, ˆl ) of all orders ofscattering at small angle approximation [6]. This approach is described in [6] on the basis of smallangle modification of spherical harmonics method (MSH).In this case the anisotropic part of angular distribution of Stokes parameters in CP-representation isexpanded with spherical functions with reference to the incident direction l̂ 0 :∞∞r2k + 1 t k r k(37)L a (τ, ˆl ) = ∑ ∑Ym (ν )f m (τ) exp(imψ ) ,m =−∞ k = 0 4πhere ν = (ˆl , ˆl ) , ψ is azimuth l̂ in the system toward l̂ .00For the strong anisotropic angular distribution its spectrum is a smooth function of harmonic indexk.
For this reason we can get the differential equation [6] for expansion coefficientsrr∂ f m (τ, k ) tt r(38)µ0+ (1 − Λxk )f m (τ, k ) = 0 ,∂τthe solution of which is simplerrtt(39)f m (τ, k ) = exp −(1 − Λxk )τ µ0 f m (0, k ) .()After the conversion in the solution to SP–representation MSH [6] could be presented in thecoordinate system about the normal to the layer boundary:MKtrtt rt2k + 1 t m(40)Π k (µ) Zk ( τ)Π km (µ 0 ) D c L0 ,L a (τ, ˆl, ˆl 0 ) = ∑ ∑ (2 − δm ,0 )φc (mϕ)∑c =1,2 m = 0k =0 4πtttttwhere D1 =diag(1 1 0 0), D 2 =diag(0 0 1 1), Zk ( τ) = exp −(1 − Λχ k )τ µ0 .{}By substituting expression (40) in (4), and after some transformations [6] we get an expression forthe source functionKt trt rttµ0 ∆ cm ( τ, µi ) = ∑ (2k + 1)(µi − µ 0 )Π mk (µi ) d k Zk (τ)Π mk (µ0 ) Dc L0k =mt tt ttt rtttt+ A mK +1 Π mK (µi ) d K +1Z K +1 (τ)Π mK +1 (µ 0 ) − Π mK +1 (µi ) d K Z K (τ)Π mK (µ 0 ) D c L0 ,7(41)Eurotherm Conference No.
95: Computational Thermal Radiation in Participating Media IVIOP PublishingJournal of Physics: Conference Series 369 (2012) 012021doi:10.1088/1742-6596/369/1/012021ttt1twhere d k = 1 − Λχk , A mk =(k 2 − m 2 )( k 2 − s 2 )δ rs .rskThis approach makes the regular part of the solution almost an isotropic function. So M can besmall, and N does not practically depend on K. We implemented this method of the solutionanisotropic part elimination as two codes: MDOM in the scalar case and MVDOM in the vector one.Both codes are realized by languages FORTRAN as well as Matlab [13].5. Numerical implementationFirst of all we note that the numerical implementation of the solution (26) is based on solving the basicproblems of linear algebra.