- Email: [email protected]
The Iterative-DRESOR method to solve radiative transfer in a plane-parallel, anisotropic scattering medium with specular-diffuse boundaries
ARTICLE IN PRESS Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
The Iterative-DRESOR method to solve radiative transfer in a plane-parallel, anisotropic scattering medium with specular-diffuse boundaries Zhi-feng Huang, Qiang Cheng, Huai-chun Zhou State Key Laboratory of Coal Combustion, Huazhong University of Science and Technology, Wuhan 430074, Hubei, PR China
a r t i c l e i n f o
abstract
Article history: Received 28 October 2008 Received in revised form 12 February 2009 Accepted 11 March 2009
Using the intensity with high directional resolution obtained by the Basic-DRESOR method as an initial guess, which is substituted into the integrated radiative transfer equation (IRTE), an iterative algorithm is proposed, called the Iterative-DRESOR method. This method can reduce the error levels of the intensity from several percent using the Basic-DRESOR method to a level of less than 1.0 106 with acceptable computation costs. The method is also validated against the exact heat flux in literature in some cases. It further clarifies some uncertain results for the reflectance in a pure, linearly anisotropic scattering medium with specular-diffuse boundaries. The directional distributions of intensity are obviously influenced by the reflecting modes of the boundary, especially in the zone near the boundary. The reflecting mode of an emitting boundary has little effect on the transmittance or reflectance. The reflecting mode of a non-emitting boundary also has little effect on the transmittance, but it obviously influences the reflectance. The difference between the reflectance for specular and diffuse boundaries increases at first, and then decreases, as the optical thickness of the medium increases. The difference will decrease as the scattering albedo of the medium increases, and it is negligible when the medium is pure scattering. The effect of the scattering phase function of the medium on the difference can also not be ignored. The Iterative-DRESOR method is expected to strengthen the capability of the Monte Carlo method to produce accurate results and to validate the results of other methods to solve RTE. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Radiative transfer equation Plane-parallel DRESOR method Specular-diffuse boundaries Iterative algorithm Anisotropic scattering
1. Introduction As a fundamental problem, radiative heat transfer in absorbing, emitting, and anisotropic scattering, plane-parallel media has been investigated by many authors and with different techniques. Busbridge and Orchard [1,2] proposed an exact method for two different modes of anisotropically scattering phase function. The SKN (Synthetic Kernel) approximation is proposed by Altac [3] for solving radiative transfer problems in linearly anisotropically scattering, homogeneous and inhomogeneous, participating media. Sarma et al. [4] used the discrete transfer method to analyze radiative heat transfer problems in a participating medium subjected to collimated radiation. Pontaza and Reddy [5] employed least-squares-based finite element formulations for the numerical solution of the radiative transfer equation in
Corresponding author. Tel.: +86 27 87557815; fax: +86 27 87540249.
E-mail address: [email protected] (H.-c. Zhou). 0022-4073/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2009.03.007
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
Nomenclature Er Emax I i j M0 n^ N0 N1 Ni q Rid Rsd r s s^ S T
difference in reflectance by specular and diffuse boundaries maximum relative difference as defined in Eq. (9) radiation intensity, W/(m2sr) discrete grids index discrete angles index number of total discrete grids normal direction of boundary number of total discrete angles number of energy bundles iterative times heat flux, W/(m2) an intermediate value as defined in Eq. 10. DRESOR value, 1/m3, or 1/m2 position vector, m distance, m unit vector along a given direction radiative source function, W/(m2sr) temperature, K
y
ss r o o1,o2 t e O F
1073
directional angles scattering coefficient, 1/m reflectivity of boundary scattering albedo parameters of the scattering phase function optical thickness emissivity of boundary angular space, sr scattering phase function
Superscripts n d s
iterative number diffuse reflection specular reflection
Subscripts w b 1 2
wall blackbody left boundary right boundary
Greek letters
b
k
attenuation coefficient, 1/m absorption coefficient, 1/m
absorbing, emitting, anisotropically scattering media, allowing for spatially varying absorption and scattering coefficients. In these papers the boundaries considered were all transparent or diffusely reflecting. But in applications some surfaces are smooth and cannot be treated as diffuse reflectors [6], for example, the ceramic surfaces in combustion chambers of highspeed aircraft; also the effect of reflecting modes of boundaries on the layer equilibrium temperatures of an absorbing–emitting layer exposed to incident radiation was analyzed. An integral theory was presented by Machali [7,8] to consider the isotropy and anisotropy of scattering and the effects of diffuse and specular reflectivity of the bounding surfaces in a plane-parallel medium. Maruyama used the radiation element method [9] to solve the problem with a specular-diffuse boundary, and both collimated and diffuse incident irradiation were analyzed. But for a planeparallel, purely linearly anisotropic scattering medium with emitting, absorbing, and specular-diffuse boundaries the reflectance got [8,10] showed obvious differences; and the exact solution for this case has not been confirmed and needs further study. Furthermore, most of these methods cannot efficiently calculate radiative intensity with high directional resolution, which provides more information for inverse radiative transfer problems. This plays a key role in radiative image processing in some industrial applications, for example, monitoring of temperature distributions in industrial combustion furnaces [11]. A new method called the DRESOR method, based on the Monte Carlo method, can directly calculate intensity with high directional resolution. It has been applied to the isotropic scattering media [12,13], and then extended to solve the problem with anisotropic scattering media [14]. The results calculated by the DRESOR method (it is hereinafter called the Basic-DRESOR method so as to distinguish it from the present DRESOR method described in this paper) agreed well with those calculated by some other methods and the integral formulation of RTE. However, the relative errors for the intensity got by the Basic-DRESOR method are several percent, and it is hard to reduce them significantly, due to the inherent feature belonging to all the Monte Carlo methods. Iterative techniques were usually used to improve the accuracy of numerical methods. The one-dimensional radiative transfer equation can be transformed into an integral equation for the source function and the resulting equation can be solved by iteration [15]. Analytical studies on the convergence of iterative schemes [16] for the radiative transfer equation show that the convergence can be very slow if the single scattering albedo is close to unity. Sutton and Ozisik [17] improved the convergence significantly by using the P1 approximation of the spherical harmonics method for the initial guess, and focused on the dimensionless radiative heat flux (reflectance and transmittance). The convergence criterion adopted in [17] was that the difference of the source term between two successive iterations was less than a specific value, such as 1.0 103. If directional quantities such as the angular distribution of radiation intensity are needed, it might be desirable to use a larger number of quadrature points [17].
ARTICLE IN PRESS 1074
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
In this paper, the iterative algorithm is directly applied to the integral form of the radiative transfer equation, and the intensity got by the Basic-DRESOR method is used as an initial guess. The description of the Iterative-DRESOR method, together with the treatment method for the specularly reflecting boundary, which has not been clarified before in the Basic-DRESOR method, will be described in detail in Section 2. In Section 3, the method was validated by the exact reflectance in literature for an anisotropic scattering medium. The convergence of the Iterative-DRESOR method will be demonstrated; and the radiative flux for isotropic and linearly anisotropic scattering media got by the present method will be compared with those in literature. Furthermore, the effect of the reflecting modes, especially the specular boundary, on the intensity and the radiative flux will be examined thoroughly. Finally some conclusions will be drawn. 2. The Iterative DRESOR method 2.1. The Basic-DRESOR method The radiative transfer equation (RTE) for an emitting, absorbing and anisotropic scattering medium is [18] Z ss Iðr; s^ i ÞFðs^ i ; s^ Þ dOi s^ rIðr; s^ Þ ¼ kIb ðr; s^ Þ bIðr; s^ Þ þ 4p 4p and its integral form for the medium is given [18] as Z s Z s Z Iðr; s^ Þ ¼ Iw ðr; s^ Þ exp b ds00 þ Sðr0 ; s^ Þ exp 0
0
s0
b ds00 b ds0
(1)
(2)
0
where s ¼ |rrw|, s0 ¼ |rr0 |, and the direction of integration is switched to go along s00 (from point r toward the wall). The source function Sðr0 ; s^ Þ is [18] Z o Iðr0 ; s^ i ÞFðs^ i ; s^ Þ dOi (3) Sðr0 ; s^ Þ ¼ ð1 oÞIb ðr0 Þ þ 4p 4p generally, for opaque surfaces with arbitrary surface properties, the boundary condition is given [18] as Z rðrw ; s^ 0 ; s^ ÞIðrw ; s^ 0 Þjn^ s^ 0 j dO0 Iw ðrw ; s^ Þ ¼ ðrw ÞIb ðrw Þ þ ^ s^ 0 o0 n
(4)
In the previous works of the authors [11–14], the DRESOR values, Rsd , were introduced to quantitatively account the integral terms in Eqs. (3) and (4), and the method was called the DRESOR method. Then, the formula for Iðr; s^ Þ by the BasicDRESOR method was given as Z 1 pðrw ÞIb ðrw Þ þ Rsd ðr0w ; rw ; s^ Þ½pðr0w ÞIb ðr0w Þ dA0 Iðr; s^ Þ ¼ p W Z s Z Rsd ðr00 ; rw ; s^ Þ½4pbð1 oÞIb ðr00 Þ dV 00 exp b ds00 þ 0 V Z Z s 1 s 0 0 0 ^ 0 4pbð1 oÞIb ðr Þ þ Rd ðrw ; r ; sÞ½pðrw ÞIb ðr0w Þ dA0 þ 0 4p W Z s0 Z þ Rsd ðr00 ; r0 ; s^ Þ½4pbð1 oÞIb ðr00 Þ dV 00 exp b ds00 ds0 (5) V
0
From Eq. (5) it is obvious that, for an emitting, absorbing and scattering medium bounded by walls with arbitrary 0 properties, once all the distributions of kðrÞ, ss ðrÞ, ðrw Þ, Fðs^ i ; s^ Þ and rðrw ; s^ ; s^ Þ are given, all the DRESOR values can be calculated. Then Iðr; s^ Þ, the angular distribution of radiative intensity at any point r, can be expressed as a function of the three-dimensional distribution of Ib ðr0 Þ inside the medium and the two-dimensional distribution of Ib ðr0w Þ on the wall surfaces. For the calculation of the DRESOR values, the tracking procedure of every bundle is based on the path-length method following the general procedure of the Radiative Energy Absorption Distribution (READ) method, a Monte Carlo method [19]. In tracking the energy bundles, the ratio of energy scattered by the medium and/or reflected by the boundary surfaces is recorded in the DRESOR method, whereas the ratio of energy absorbed by the medium and boundary is recorded in the READ method. Once the intensity with high directional resolution is obtained, it is easy to deduce all the associated quantities, for example, the radiative flux. The description for the Basic-DRESOR method can be found in the literature [11–14], and it will not be repeated here again. 2.2. The Iterative-DRESOR method In the Iterative-DRESOR method, the intensity calculated using the Basic-DRESOR method is taken as an initial guess for an iterative algorithm, marked as Ið0Þ ðr; s^ Þ, and it will be substituted into Eqs. (2)–(4), then a new intensity with the same high directional resolution will be obtained.
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1075
Let IðnÞ ðr; s^ Þ be the intensity after the nth iteration, n ¼ 1,2,y. At first, it is substituted into the RHS of Eqs. (3) and (4) to ðrw ; s^ Þ as below: obtain an updated source function Sðnþ1Þ ðr; s^ Þ and boundary intensity Iðnþ1Þ w Z o IðnÞ ðr0 ; s^ i ÞFðs^ i ; s^ Þ dOi (6) Sðnþ1Þ ðr0 ; s^ Þ ¼ ð1 oÞIb ðr0 Þ þ 4p 4p Iðnþ1Þ ðrw ; s^ Þ ¼ ðrw ÞIb ðrw Þ þ w
Z ^ s^ 0 o0 n
rðrw ; s^ 0 ; s^ ÞIðnÞ ðrw ; s^ 0 Þjn^ s^ 0 j dO0
(7)
Then an updated intensity Iðnþ1Þ ðr; s^ Þ can be calculated by substituting Sðnþ1Þ ðr; s^ Þ and Iðnþ1Þ ðrw ; s^ Þ into the RHS of Eq. (2). w Z ^ Iðnþ1Þ ðr; s^ Þ ¼ Iðnþ1Þ ðr; s Þ exp w
s
b ds00 þ 0
Z
s
Z Sðnþ1Þ ðr0 ; s^ Þ exp
0
s0
b ds00 b ds0
(8)
0
The iteration will end when the maximum updating ratio for the intensity is smaller than a pre-determined value, such as 1.0 106. Emax ¼ maxfjIðnþ1Þ ðr; s^ Þ IðnÞ ðr; s^ Þj=IðnÞ ðr; s^ Þgo1:0 106
(9)
r;s^
If the condition of Eq. (9) is satisfied, it can be confirmed that the iterative algorithm is converged. This means not only that the solution of the RTE with high directional resolution is obtained, but also the accuracy of the solution reaches a very high level with relative errors less than 1.0 106. Therefore, the solution for the intensity can be regarded as a numerically exact one.
2.3. Treatment for the specularly reflecting boundary For the specular reflector the energy bundle will be reflected only at a specific angle equal to the incident angle on the opposite side. This causes the treatment for the specular reflecting boundary to be different from that for the diffuse reflecting boundary as described [14] and as discussed below. Before calculating the DRESOR values Rsd ðr; rw ; s^ j Þ for the specular reflecting boundary, Rid ðr; rw ; s^ i Þ, an intermediate variable is used to record the ratio of the energy incident from a unit solid angle around the direction s^ i to a unit area around rw to that emitted from a unit area or a unit volume around the point r. As the calculation begins, it is set to zero. N1 energy bundles are tracked, and the initial energy is set as E0 ¼ 1.0, experiencing gradual reduction as it is absorbed along its traveling path. During the tracking process, it is renewed by Rid;new ðr; rw ; s^ i Þ ¼ Rid;old ðr; rw ; s^ i Þ þ
E0 ðs^ i Þ N1
(10)
where E0 ðs^ i Þ means the residual energy possessed by the bundle arriving at the boundary from direction s^ i . For the specular reflecting boundary, the DRESOR values can be calculated as Rsd ðr; rw ; s^ j Þ ¼ rsrw Rid ðr; rw ; s^ M0 jþ1 Þ p=ðdO cos yj Þ
(11)
where dO is the solid angle for the discrete angle. For the diffuse reflecting boundary Rsd ðr; rw ; s^ j Þ ¼
M0 X
rdrw Rid ðr; rw ; s^ i Þ
(12)
i¼1
The other parts of treatment for the specular reflecting boundary are the same as those for the diffuse reflecting boundary.
Table 1 Calculation cases and conditions. Case Parameters studied 1
Reflectance q~ ¼ 1 qð0Þ=ðsT 41 Þ
2
Heat flux q~ (in units of 2e1sT41)
3
Radiative intensity, dimensionless reflectance and transmittance
Boundary 1 T1 ¼ 1000 K; e1 ¼ 1.0 T1 ¼ 1000 K; specularly reflecting T1 ¼ 1000 K; diffusely or specularly reflecting
Medium
t ¼ 1.0–10.0; o ¼ 1.0; F(y) ¼ 1+o1P1(cos y)+o2P2(cos y) t ¼ 0.01–5.0; o ¼ 1.0; F(y) ¼ 1+x cos y, x ¼ 1, 0, 1
t ¼ 0.01–5.0; o ¼ 0.0, 0.5, or 1.0; F(y) ¼ 1+x cos y, x ¼ 1, 0, 1
Boundary 2 T2 ¼ 0 K; e2 ¼ 1.0 T2 ¼ 500 K; e2 ¼ 0.8; diffusely reflecting T2 ¼ 0 K; diffusely or specularly reflecting
References Busbridge and Orchard [1,2] Machali [7,8]; Tan et al. [10].
ARTICLE IN PRESS 1076
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
3. Results and discussion 3.1. Objective, cases and conditions of calculation The Iterative DRESOR method is applied to plane-parallel and isotropic/anisotropic scattering media. Three cases are studied in the paper, as listed in Table 1. In Case 1, two boundaries are black with ¼ 1:0 and with two different temperatures T1 ¼ 1000 K, T2 ¼ 0. The medium is pure scattering with o ¼ 1:0. Nonlinearly scattering phase function F(y) ¼ 1+o1P1(cos y)+o2P2(cos y) is considered. The reflectance will be calculated. This case is used to validate the accuracy of the Iterative-DRESOR method. In Case 2, the medium is also pure scattering, and two boundaries are diffusely emitting, specularly reflecting and diffusely reflecting with temperatures T1 ¼ 1000 K, T2 ¼ 500 K. The linearly anisotropic scattering phase function (F(y) ¼ 1+x cos y , x ¼ 1, 0, 1) is considered. When x ¼ 0 the phase function turns into isotropic scattering. The emissivity of Boundary 2 is 0.8 and it is diffuse reflecting with rd2 ¼ 0:2. The emissivity of boundary 1 is variable, and it has specular reflectance rs1 ¼ 1 1 . In this case, there is an exact solution of the heat flux for the isotropic scattering, which can be used to verify the present method. There is not an exact solution for the anisotropic scattering, and some uncertainty in the solution will be clarified.
Table 2 The time taken to calculate the intensity with the phase function F1 for different optical thickness in Case 1.
t0
N0
1 2 3 4 5 6 7 8 9 10
M0
100 100 200 200 300 300 400 400 500 500
Calculation time by BasicDRESOR method (s)
180 180 180 180 180 180 180 180 180 180
32 50 99 107 194 210 311 326 438 453
Iterative number and calculation time by Iterative-DRESOR method after convergence Iterative number
Calculation time (s)
22 41 61 88 117 150 178 226 263 309
14 27 126 183 498 639 1291 1689 2903 3421
Table 3 Reflectance for the phase function F(y) ¼ 1+o1P1(cos y)+o2P2(cos y) for various optical thicknesses in Case 1.
t0
1 2 3 4 5 6 7 8 9 10
F0 (o1 ¼ 0, o2 ¼ 0) Exact
DRESOR
0.4466 0.6099 0.6984 0.7540 0.7923 0.8203 0.8417 0.8585 0.8721 0.8833
0.43646 0.58316 0.67332 0.72334 0.76488 0.78600 0.80974 0.82406 0.83922 0.84600
F1 (o1 ¼ 0, o2 ¼ 0.5) Present work 2.270% 4.384% 3.591% 4.066% 3.461% 4.181% 3.797% 4.012% 3.770% 4.223%
0.44659 0.60995 0.69838 0.75408 0.79241 0.82041 0.84177 0.85859 0.87219 0.88341
0.002% 0.008% 0.003% 0.011% 0.014% 0.013% 0.008% 0.010% 0.010% 0.012%
F2 (o1 ¼ 1.0, o2 ¼ 0.5) 1 2 3 4 5 6 7 8 9 10
0.3580 0.5157 0.6104 0.6739 0.7197 0.7541 0.7810 0.8026 0.8204 0.8352
0.34864 0.49419 0.58729 0.64481 0.69381 0.72562 0.75477 0.77006 0.79215 0.80330
Exact
DRESOR
0.4468 0.6101 0.6985 0.7541 0.7924 0.8204 0.8417 0.8585 0.8721 0.8833
0.43639 0.58686 0.68043 0.72624 0.76270 0.78789 0.81242 0.82399 0.84087 0.85132
Present work 2.330% 3.516% 2.587% 3.694% 3.748% 3.963% 3.479% 4.020% 3.581% 3.621%
0.44682 0.61013 0.69851 0.75417 0.79247 0.82046 0.84180 0.85862 0.87221 0.88343
0.004% 0.005% 0.001% 0.009% 0.009% 0.007% 0.012% 0.014% 0.013% 0.015%
3.579% 4.833% 3.781% 4.187% 4.202% 4.281% 3.998% 4.201% 4.008% 4.065%
0.30202 0.44900 0.54371 0.61046 0.66015 0.69860 0.72924 0.75423 0.77500 0.79254
0.007% 0.000% 0.002% 0.010% 0.008% 0.014% 0.005% 0.017% 0.013% 0.018%
F3 (o1 ¼ 1.5, o2 ¼ 0.5) 2.615% 4.171% 3.786% 4.317% 3.597% 3.777% 3.359% 4.054% 3.443% 3.819
0.35803 0.51572 0.61039 0.67398 0.71972 0.75420 0.78114 0.80276 0.82050 0.83532
0.008% 0.004% 0.002% 0.012% 0.003% 0.013% 0.018% 0.020% 0.012% 0.014%
0.3020 0.4490 0.5437 0.6104 0.6601 0.6985 0.7292 0.7541 0.7749 0.7924
0.29119 0.42730 0.52314 0.58484 0.63236 0.66860 0.70005 0.72242 0.74384 0.76019
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1077
In Case 3, the boundary conditions are the same as those in Case 2 except the temperature of Boundary 2 is set to 0, that is T2 ¼ 0.0. The medium is absorbing, scattering, and non-emitting with scattering phase function F(y) ¼ 1+x cos y , x ¼ 1, 0, 1. This case is used to compare the effects of different reflecting modes on the radiative intensity, the transmittance and the reflectance.
0.1 0.01
Emax
1E-3 1E-4 1E-5 1E-6 0
2
4
6
8
10
12
14
16
18
20
22
Ni 0.1 0.01
Emax
1E-3 1E-4 1E-5 1E-6 0
30
60
90
120 150 180 210 240 270 300 Ni
Fig. 1. The maximum error in intensity with phase function F1, (a) t0 ¼ 1.0 and (b) t0 ¼ 10.0 in Case 1.
Table 4 The time taken to calculate the phase function F(y) ¼ 1cos y with rs2 ¼ rd1 ¼ 0, rd2 ¼ 0:2, e2 ¼ 0.8 and different e1 in Case 2.
t0
N0
M0 Calculation time by Basic-DRESOR method (s) Iterative numbers and calculation time by Iterative-DRESOR method after convergence Iterative numbers
e1 ¼ 0.2 0.01 10 180 21 0.1 20 180 46 0.5 50 180 131 1.0 100 180 280 2.0 100 180 326 5.0 300 180 1081
Calculation time (s)
e1 ¼ 0.7
e1 ¼ 1.0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
17 31 89 190 200 678
2 5 16 43 53 270
8 14 29 46 89 233
6 10 20 30 57 153
4 8 16 25 49 139
0.3 1 7 30 60 1029
0.3 0.8 5 20 39 711
0.2 0.6 4 16 32 630
ARTICLE IN PRESS 1078
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
3.2. Validation of the Iterative-DRESOR method The Iterative-DRESOR method is validated in Case 1. Participating media with different optical thicknesses are subjected to a diffuse irradiation flux from Boundary 1. The medium is divided into N0 discrete grids and M0 discrete angles in the angle range of [0, p]. As the optical thickness increases, N0 will also increases to maintain accurate calculation. The number of energy bundles used in the DRESOR method is 100,000. The discrete number and times consumed with different optical thicknesses are listed in Table 2. The reflectance q~ ¼ 1 qð0Þ=ðsT 4s1 Þ is compared with the exact solution obtained by Busbridge and Orchard [1,2]; and the calculation results from the present method are compared with the exact solutions, as listed in Table 3. The comparisons are made for four different phase functions with different parameters o1 and o2 , and F0 is isotropic scattering and F1 is Rayleigh scattering phase functions. The calculations in this paper are all carried out on a Pentium IV, 3.0 GHz PC. As shown in Table 2, the new method does not consume much more time compared with the Basic-DRESOR method, especially when the optical thickness is not large, where only a few seconds are needed for the iterations. The iteration number increases as the optical thickness increasing, and the time consumed increases rapidly. But as shown in Fig. 1, Emax decreases continuously as the iteration time increases. The curve displays as an exponential curve after several iterations. As the optical thickness is 10.0, the calculation time by the Iterative-DRESOR method is 3421 s, being 7.55 times than 453 s by the Basic-DRESOR method. Compared with the reduction of calculation errors from the level of 1:0 102 to that of 1:0 106 , the consumption in time is a worthy cost. If the accuracy of the results is not required to be very high, the time taken can be reduced. For example, if the accuracy needed is 1:0 103 in Fig. 1(b), the iteration number needed is only 12 and the time consumed for iteration is only 136 s. The results mentioned above are done for a pure scattering medium. If the scattering albedo decreases, the iteration number will decrease quickly, and the time used can be further reduced. The exact reflectance in Case 1 from literature [1,2] and those calculated by the Basic-DRESOR method and the IterativeDRESOR method and their errors are listed in Table 3 for the four phase functions. Compared with the exact solutions, the accuracy of the Iterative-DRESOR method is significantly improved. The errors are all less than 0.02%. The comparison shows that the present method is very accurate without much increase in computing time.
Table 5 The dimensionless heat flux for isotropic scattering with rs2 ¼ rd1 ¼ 0, rd2 ¼ 0.2, e2 ¼ 0.8 and different e1 in Case 2.
t0
0.01 0.1 0.5 1.0 2.0 5.0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
Machali [8]
Present work
Machali [8]
Present work
Machali [8]
Present work
Heat flux
Heat flux
Relative errors (%)
Heat flux
Heat flux
Relative errors (%)
Heat flux
Heat flux
Relative errors (%)
0.44559 0.43851 0.41229 0.38554 0.34261 0.25772
0.445607 0.438524 0.412300 0.385547 0.342593 0.257644
0.004 0.003 0.002 0.002 0.005 0.029
0.39661 0.37798 0.31827 0.26854 0.20591 0.12165
0.396620 0.377987 0.318273 0.268539 0.205887 0.121600
0.003 0.002 0.001 0.000 0.011 0.041
0.37208 0.34928 0.28067 0.22788 0.16660 0.09254
0.372087 0.349286 0.280674 0.227880 0.166584 0.092497
0.001 0.002 0.001 0.000 0.010 0.046
Table 6 The dimensionless heat fluxes for linear anisotropic scattering with rs2 ¼ rd1 ¼ 0, rd2 ¼ 0.2, e2 ¼ 0.8 and different e1 in Case 2.
t0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
Machali [8]
Tan et al. [10]
DRESOR
Machali [8]
Tan et al. [10]
DRESOR
Machali [8]
Tan et al. [10]
DRESOR
0.44580 0.44055 0.42084 0.39918 0.36008 0.27062
0.445854 0.440978 0.422416 0.402770 0.369657 0.297461
0.445818 0.440585 0.421571 0.402085 0.369615 0.298730
0.39720 0.38335 0.33683 0.29299 0.22909 0.13191
0.397343 0.384319 0.339727 0.299160 0.242866 0.155902
0.397208 0.383397 0.338376 0.298464 0.243298 0.157344
0.37282 0.35586 0.30168 0.25370 0.18888 0.10124
0.372997 0.356935 0.304539 0.259759 0.201875 0.121452
0.372826 0.355916 0.303382 0.259410 0.202588 0.122814
x ¼ 1.0 0.01 0.44538 0.1 0.43649 0.5 0.40443 1.0 0.37421 2.0 0.33032 5.0 0.25047
0.445319 0.436010 0.402241 0.369110 0.318571 0.226868
0.445395 0.436483 0.403428 0.370316 0.319252 0.226488
0.39602 0.37277 0.30234 0.25005 0.19105 0.11614
0.395867 0.371612 0.298426 0.242584 0.177949 0.099421
0.396034 0.372727 0.300424 0.244068 0.178447 0.099085
0.37134 0.34295 0.26315 0.20903 0.15274 0.08791
0.371146 0.341579 0.259136 0.201901 0.141122 0.074488
0.371350 0.342898 0.261129 0.203184 0.141445 0.074181
x ¼ 1.0 0.01 0.1 0.5 1.0 2.0 5.0
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1079
3.3. Radiative heat transfer with specular-diffuse boundaries Consider a pure scattering plane-parallel medium in Case 2 with isotropic scattering. For different optical thicknesses, the discrete number of the medium, the calculation time by the Basic-DRESOR and the Iterative-DRESOR methods, and the iterative numbers, are all listed in Table 4. It can be seen that when the optical thickness is very small (o1.0), the time taken for the iterations in the Iterative-DRESOR method is negligible compared with that taken for the Basic-DRESOR method, but the accuracy is considerably improved. The dimensionless heat flux q~ (in units of 2e1sT41) with isotropic scattering is given in Table 5. Results calculated by the Iterative-DRESOR method agree well with those given by Machali [7] for the isotropic scattering, and the biggest error is less than 0.05%. The dimensionless heat flux from literature [8,10] and obtained by the Iterative-DRESOR method for two linear anisotropic scattering phase functions with different e1 are given in Table 6. When the optical thickness is bigger than 1.0, the results given by Machali [8] are different from those calculated by Tan et al. [10] and the present method. The results obtained by the Iterative-DRESOR method, a numerically exact method with relative errors of 1.0 106, agree very well with those of Tan et al. [10]. So, the method in [8] should be carefully checked when it is applied to large optical thickness media with anisotropic scattering. 3.4. Effect of specular reflecting boundary on the intensity The effects of different reflecting modes of boundaries on radiative intensity with high directional resolution are examined in Case 3, where N0 ¼ 100, and M0 ¼ 180. The intensity with high directional resolution calculated by the Iterative-DRESOR method and their relative errors are all shown in Figs. 2–4. i, j denote the discrete numbers of grid and angle, respectively. In these cases the medium is pure scattering with linear forward scattering and the optical thickness
180 160 140 2.2E3
3.1E3
120
2.8E3 2.5E3 3.5E3 3.8E3 4.1E3 4.8E3 4.4E3 5.1E3 5.4E3 6.1E3 5.7E3
j
100 80 60
6.4E3
40 20 0
10
20
30
40
50
60
70
80
90 100
i
1.0x10-6 9.0x10-7
error
8.0x10-7 7.0x10-7 6.0x10-7 5.0x10-7 4.0x10-7 3.0x10 10 20 30 40 50 60 70 80 i 90 100
180 160 140 120 100 80 60
j
-7
20
40
Fig. 2. (a) Intensity and (b) its relative errors for linear forward scattering with rd1 ¼ rs2 ¼ 0, rs1 ¼ 0:8 and rd2 ¼ 0:2 in Case 3.
ARTICLE IN PRESS 1080
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
t ¼ 1.0. The emissivity of the two boundaries are e1 ¼ 0.2, and e2 ¼ 0.8, respectively. The main difference lies in the reflecting modes of the two boundaries. In Fig. 2, Boundary 1 is a specularly reflecting surface and Boundary 2 is a diffusely reflecting surface with rs1 ¼ 0:8, rd2 ¼ 0:2. In Fig. 3, the two boundaries are diffusely reflecting surfaces with rd1 ¼ 0:8, rd2 ¼ 0:2. In Fig. 4, Boundary 1 is diffusely reflecting and Boundary 2 is a specularly reflecting surface with rd1 ¼ 0:8, rs2 ¼ 0:2. From these figures, it can be seen firstly that the intensity obtained by the Iterative-DRESOR method is very accurate, and the relative errors are less than 1.0 106. Secondly, the difference in intensity caused by different reflecting modes is also obvious, especially in the zone near the specular reflecting boundary. In order to show the difference in detail, the directional intensity at the two boundaries is plotted in Fig. 5(a) and (b). As shown in Fig. 5(a), the intensity at Boundary 1 around 90 degrees in the forward directions is larger with specular reflecting than that with a diffuse reflecting boundary. In contrast, as shown in Fig. 5(b), the intensity at Boundary 2 in the backward direction around 901 with specular reflecting boundary is smaller than that with a diffuse reflecting boundary. This varying distribution of directional intensity also appears in the angle range of [p/2, p] if Boundary 2 is specular reflecting; but if Boundary 2 is diffuse reflecting the directional intensity in the range will become more uniform. 3.5. Effect of specular reflecting boundary surfaces on the flux Researchers are usually more concerned about radiative flux. The effects of specular reflecting boundary on the reflectance and transmittance in Case 3 are shown in Table 7, where the medium is linearly forward scattering with o ¼ 0.5 and t ¼ 1.0. It can be seen that the reflecting mode of Boundary 1, the emitting boundary, has little effect on the reflectance and transmittance. The differences of reflectance and transmittance are all less than 1.52% with different reflecting mode of Boundary 1. The reflecting mode of Boundary 2, the non-emitting boundary, also has little effect on the transmittance, and the differences are all less than 1.0%. It should be stated that the effect of different reflecting modes of Boundary 2 on the
180 160 140 3.3E3
2.0E3
2.8E3
120
2.5E32.3E3
3.0E3
100
3.5E3 3.8E3 4.3E34.0E3 4.5E3 5.0E3
j
4.8E3
80 60
5.8E3
6.3E3
40
5.5E3 5.3E3
6.0E3
20 0
10
20
30
40
50
60
70
80
90 100
i
1.0x10-6 9.0x10-7
error
8.0x10-7 7.0x10-7 6.0x10-7 5.0x10-7
-7
3.0x10
10 20 30 40 50 60 70 80 i 90 100 40 20
180 160 140 120 100 80 60
j
4.0x10-7
Fig. 3. (a) Intensity and (b) its relative errors for linear forward scattering with rs1 ¼ rs2 ¼ 0, rd1 ¼ 0.8 and rd2 ¼ 0:2 in Case 3.
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1081
reflectance is obvious, and the differences are all bigger than 4% when r2 ¼ 0.8 and bigger than 2% when r2 ¼ 0.2. If the boundary’s reflectivity is equal to zero, there will be no reflectance, and consequently, no difference. The difference of reflectance between specular reflecting and diffuse reflecting of Boundary 2 is hereinafter expressed as Er for brevity. As shown in Table 7, as the reflectivity of Boundary 2 increases, Er increases. In addition, there are some other factors affecting the difference, such as the optical thickness, the scattering albedo, and the scattering phase function of the medium, which will be given below. Er variation with the optical thickness is shown in Fig. 6 with scattering albedo o ¼ 0.5, and the emissivities of the two boundaries being e1 ¼ 1.0, and e2 ¼ 0.2. As shown in Fig. 6, the difference increases as the optical thickness increases, but once the optical thickness has reached a specific value, the difference decreases as the optical thickness increases. As shown in Fig. 5(b), the directional intensity around 01 is larger than that in other directions, and less energy will be absorbed after it is reflected by the specular reflecting Boundary 2, leading to a larger reflectance than that of a diffuse reflecting Boundary 2. If the optical thickness of the medium increases further, less energy reflected by Boundary 2 can reach Boundary 1, and Er becomes smaller and smaller. The variation of Er with o ¼ 1.0, 0.5, and 0.0 is listed in Table 8, where the medium is isotropically scattering with t ¼ 1.0, e1 ¼ 1.0, and e2 ¼ 0.2. Er increases from 0.08% to 25.25% as the scattering albedo decreases from 1.0 to 0.0. When the scattering albedo becomes larger, less energy emitted from Boundary 1 can reach Boundary 2, and Er will decrease. If the medium is pure scattering, Er is negligible. The influence of phase function on Er is shown in Table 9. It can be seen that the effect of the phase function on the reflectance cannot be ignored, not only for the same reflecting condition of the boundary, but also between different reflecting modes of the boundary. It has been stated that Monte Carlo method can be used to solve tedious and difficult problems in radiative transfer, and because the method is known to produce accurate results if enough sample packets are used, the Monte Carlo method is
180 160 3E3
140
2.3E32.0E3
120
3.0E3
100 j 80 6.3E3
2.5E3
4.0E33.8E33.5E3 4.3E3 4.8E3
4.5E3
60
2.8E3
5.8E35.5E3 6.0E3
5.3E3
5.0E3
40 20 0
10
20
30
40
50
60
70
80
90 100
i
1.0x10-6 9.0x10-7
error
8.0x10-7 7.0x10-7 6.0x10-7 5.0x10-7 4.0x10-7 3.0x10 10 20 30 40 50 60 70 80 i 90 100
180 160 140 120 100 80 60
j
-7
20
40
Fig. 4. (a) Intensity and (b) its relative errors for linear forward scattering with rs1 ¼ rd2 ¼ 0, rd1 ¼ 0:8 and rs2 ¼ 0:2 in Case 3.
ARTICLE IN PRESS 1082
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
often used to validate the results of other methods [20]. The Iterative-DRESOR method, as a development of Monte Carlo method, will strengthen the capability of the Monte Carlo method in this function. In order to expand this methodology to more complex realistic cases, especially the spectral problems, such as FSK (full-spectrum k-distribution) method [21] and the LBL (line-by-line) method [22,23] can be adopted in the DRESOR method in the future work. 4. Conclusion The Iterative-DRESOR method is an iterative algorithm applied to the integral form of the radiative transfer equation using an initial guess of intensity with high directional resolution calculated by the Basic-DRESOR method. The treatment of the specular reflecting boundary in tracking the energy bundles is proposed, which has not been clarified before in the Basic-DRESOR method. For a plane-parallel, anisotropic scattering medium, bounded with black boundaries with different specularly reflecting
diffusely reflecting 90
Intensity at boundary1
8000
120
60
7000 6000 30
150
5000 4000 3000 2000
180
0 specularly reflecting
diffusely reflecting 90
Intensity at boundary2
5000
120
60
4000 3000
30
150
2000 1000 0 180
0
Fig. 5. (a) Intensity at Boundary 1 with different reflecting modes of Boundary 1, and (b) at Boundary 2 with different reflecting modes of Boundary 2 in Case 3.
Table 7 The effect of different reflecting modes of the two boundaries on dimensionless radiative flux with linearly forward scattering medium t ¼ 1.0, o ¼ 0.5 in Case 3. Conditions
r1 ¼ 0.0, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.8 r1 ¼ 0.8, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.0
r1 ¼ 0.0, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.8 r1 ¼ 0.8, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.0 a
Two diffuse boundaries
Specular Boundary 1, diffuse Boundary 2 a
Reflectance
Reflectance
Relative Error (%)
0.19480 0.16217 0.04611 0.09688 0.02608 0.02044
0.16215 0.04614 0.09706 0.02629 0.02075
0.01 0.07 0.19 0.81 1.52
Transmittance
Transmittance
Relative errora (%)
0.07359 0.07654 0.08718 0.28407 0.30634 0.36801
0.07657 0.08705 0.28348 0.30394 0.36447
0.04 0.15 0.21 0.78 0.96
Relative error to the corresponding data for two diffuse boundaries.
Diffuse Boundary 1, specular Boundary 2 Reflectance 0.20344 0.16965 0.04855 0.09900 0.02669
Relative Errora (%) 4.44 4.61 5.29 2.19 2.34
Transmittance
Relative errora (%)
0.07289 0.07599 0.08707 0.28348 0.30623
0.95 0.72 0.13 0.21 0.04
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1083
5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0
1
2
3
4
5
Fig. 6. The effect of different optical thickness with linear forward scattering o ¼ 0.5, e1 ¼ 1.0, and e2 ¼ 0.2 in Case 3.
Table 8 The effect of different scattering albedo on reflectance with isotropic scattering medium t ¼ 1.0, e1 ¼ 1.0, and e2 ¼ 0.2 in Case 3.
Reflectance with diffuse Boundary 2 Reflectance with specular Boundary 2 Er
o ¼ 1.0
o ¼ 0.5
o ¼ 0.0
0.82785 0.82847 0.08%
0.21846 0.22629 3.58%
0.03846 0.04817 25.25%
Table 9 The effect of different phase functions on reflectance with t ¼ 1.0, o ¼ 0.5, e1 ¼ 1.0, and e2 ¼ 0.2 in Case 3.
Reflectance with diffuse Boundary 2 Reflectance with specular Boundary 2 Er
Forward scattering
Isotropic scattering
Backward scattering
0.19480 0.20344 4.44%
0.21846 0.22629 3.58%
0.24016 0.24728 2.96%
temperatures, exact solutions for the reflectance are used to validate the accuracy of the new method, and the convergence of the new method is demonstrated. For media with isotropic and linearly anisotropic scattering, the radiative flux obtained by the present method is compared with those in the literature. Furthermore, the effects of the different reflecting boundaries, especially the specular reflecting boundary, on the radiative intensity and flux are examined thoroughly, and this paper focuses on the difference of reflectance between specular and diffuse reflecting boundaries. The main conclusions that are drawn are listed below. (1) The Iterative-DRESOR method can give intensity with high directional resolution at a high accuracy level with relative errors less than 1.0 106 with acceptable computation costs even for an optical thickness of 10.0. For optical thicknesses less than 1.0, the increase in the computation cost is negligible, but the accuracy is improved significantly. (2) The reflectance in a pure, anisotropic scattering medium with black boundaries obtained by the present method is validated by the exact data in the literature, and agreement at a high level is confirmed. (3) For a pure, anisotropic scattering medium with specular and diffuse reflecting boundaries, the numerically-exact results for the reflectance obtained by the present method clarified some uncertainty over the correct solution in the literature. (4) The directional distributions of intensity are obviously influenced by the reflecting modes of the boundaries, especially in the zone near the boundary, which may reflect specularly or diffusely. (5) The reflecting mode of the emitting boundary has little effect on the transmittance or reflectance. The reflecting mode of the non-emitting boundary also has little effect on the transmittance, but it obviously influences the reflectance. (6) The difference between the reflectance by the specular and diffuse boundaries increases at first, and then decreases, as the optical thickness of the medium increases. The difference will decrease as the scattering albedo of the medium
ARTICLE IN PRESS 1084
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
increases, and it is negligible when the medium is pure scattering. Also the effect of the scattering phase function of the medium on the difference cannot be ignored. From the conclusions above, we can expect that the Iterative-DRESOR method can strengthen the capability of the Monte Carlo method to produce accurate results and to validate the results of other methods.
Acknowledgment The study has been supported by the National Natural Science Foundation of China (No. 50636010, No. 50721005). References [1] Busbridge IW, Orchard SE. Reflection and transmission of light by a thick atmosphere according to phase function: 1+x cos y. The Astrophysical Journal 1967;149:655–64. [2] Orchard SE. Reflection and transmission of light by thick atmospheres of pure scatterers with a phase function: 1+w1P1(cos y)+w2P2(cos y). The Astrophysical Journal 1967;149:665–74. [3] Altac Z. The SKN approximation for solving radiative transfer problems in absorbing, emitting, and linearly anisotropically scattering plane-parallel medium: part 2. Journal of Heat Transfer 2002;124(4):685–95. [4] Sarma D, Mishra SC, Mahanta P. Analysis of collimated radiation in participating media using the discrete transfer method. Journal of Quantitative Spectroscopy and Radiative Transfer 2005;96(1):123–35. [5] Pontaza JP, Reddy JN. Least-squares finite element formulations for one-dimensional radiative transfer. Journal of Quantitative Spectroscopy and Radiative Transfer 2005;95(3):387–406. [6] Siegel R, Spuckler CM. Effects of refractive index and diffuse or specular boundaries on a radiating isothermal layer. Journal of Heat Transfer 1994;116(3):787–90. [7] Machali HF. Radiative transfer in participating media under conditions of radiative equilibrium. Journal of Quantitative Spectroscopy and Radiative Transfer 1995;53(2):201–10. [8] Machali HF, Madkour MA. Radiative transfer in a participating slab with anisotropic scattering and general boundary conditions. Journal of Quantitative Spectroscopy and Radiative Transfer 1995;54(5):803–13. [9] Maruyama S. Radiative heat transfer in anisotropic scattering media with specular boundary subjected to collimated irradiation. International Journal of Heat and Mass Transfer 1998;41(18):2847–56. [10] Tan H-P, Yi H-L, Wang PY, et al. Ray tracing method for transient coupled heat transfer in an anisotropic scattering layer. International Journal of Heat and Mass Transfer 2004;47(19–20):4045–59. [11] Zhou H-C, Lou C, Cheng Q, et al. Experimental investigations on visualization of three-dimensional temperature distributions in a large-scale pulverized-coal-fired boiler furnace. Proceedings of the Combustion Institute 2005;30(1):1699–706. [12] Zhou H-C, Chen D-L, Cheng Q. A new way to calculate radiative intensity and solve radiative transfer equation through using the Monte Carlo method. Journal of Quantitative Spectroscopy and Radiative Transfer 2004;83(3–4):459–81. [13] Cheng Q, Zhou H-C. The DRESOR method for a collimated irradiation on an isotropically scattering layer. Journal of Heat Transfer, Transactions of the ASME 2007;129(5):634–45. [14] Zhou H-C, Cheng Q, Huang Z-F, et al. The influence of anisotropic scattering on the radiative intensity in a gray, plane-parallel medium calculated by the DRESOR method. Journal of Quantitative Spectroscopy and Radiative Transfer 2007;104(1):99–115. [15] Ozisik MN. Radiative transfer and interactions with conduction and convection. New York: Wiley-Interscience; 1973. [16] Pao CV. An iterative method for radiative transfer is a slab with specularly reflecting boundary. Applied Mathematics and Computation 1975;1(4):353–64. [17] Sutton WH, Ozisik MN. An iterative solution for anisotropic radiative transfer in a slab. Journal of Heat Transfer 1979;101:695–8. [18] Modest MF. Radiative heat transfer. 2nd ed. San Diego: Academic Press; 2003. [19] Yang W-J, Taniguchi H, Kudo K. Radiative heat transfer by the Monte Carlo method. In: Hartnett JP, Irvines FT, editors. Advances in heat transfer, vol. 27. San Diego: Academic Press; 1995. p. 1–215. [20] Howell JR. The Monte Carlo method in radiative heat transfer. Journal of Heat Transfer—Transactions of the ASME 1998;120(3):547–60. [21] Wang L, Yang J, Modest MF, Haworth DC. Application of the full-spectrum k-distribution method to photon Monte Carlo solvers. Journal of Quantitative Spectroscopy and Radiative Transfer 2006;104(2):297–304. [22] Tang KC, Brewster MQ. Analysis of molecular gas radiation: real gas property effects. Journal of Thermophysics and Heat Transfer 1999;13(4):460–6. [23] Wang A, Modest MF. Spectral Monte Carlo models for nongray radiation analyses in inhomogeneous participating media. Journal of Quantitative Spectroscopy and Radiative Transfer 2007;50(19–20):3877–89.
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
The Iterative-DRESOR method to solve radiative transfer in a plane-parallel, anisotropic scattering medium with specular-diffuse boundaries Zhi-feng Huang, Qiang Cheng, Huai-chun Zhou State Key Laboratory of Coal Combustion, Huazhong University of Science and Technology, Wuhan 430074, Hubei, PR China
a r t i c l e i n f o
abstract
Article history: Received 28 October 2008 Received in revised form 12 February 2009 Accepted 11 March 2009
Using the intensity with high directional resolution obtained by the Basic-DRESOR method as an initial guess, which is substituted into the integrated radiative transfer equation (IRTE), an iterative algorithm is proposed, called the Iterative-DRESOR method. This method can reduce the error levels of the intensity from several percent using the Basic-DRESOR method to a level of less than 1.0 106 with acceptable computation costs. The method is also validated against the exact heat flux in literature in some cases. It further clarifies some uncertain results for the reflectance in a pure, linearly anisotropic scattering medium with specular-diffuse boundaries. The directional distributions of intensity are obviously influenced by the reflecting modes of the boundary, especially in the zone near the boundary. The reflecting mode of an emitting boundary has little effect on the transmittance or reflectance. The reflecting mode of a non-emitting boundary also has little effect on the transmittance, but it obviously influences the reflectance. The difference between the reflectance for specular and diffuse boundaries increases at first, and then decreases, as the optical thickness of the medium increases. The difference will decrease as the scattering albedo of the medium increases, and it is negligible when the medium is pure scattering. The effect of the scattering phase function of the medium on the difference can also not be ignored. The Iterative-DRESOR method is expected to strengthen the capability of the Monte Carlo method to produce accurate results and to validate the results of other methods to solve RTE. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Radiative transfer equation Plane-parallel DRESOR method Specular-diffuse boundaries Iterative algorithm Anisotropic scattering
1. Introduction As a fundamental problem, radiative heat transfer in absorbing, emitting, and anisotropic scattering, plane-parallel media has been investigated by many authors and with different techniques. Busbridge and Orchard [1,2] proposed an exact method for two different modes of anisotropically scattering phase function. The SKN (Synthetic Kernel) approximation is proposed by Altac [3] for solving radiative transfer problems in linearly anisotropically scattering, homogeneous and inhomogeneous, participating media. Sarma et al. [4] used the discrete transfer method to analyze radiative heat transfer problems in a participating medium subjected to collimated radiation. Pontaza and Reddy [5] employed least-squares-based finite element formulations for the numerical solution of the radiative transfer equation in
Corresponding author. Tel.: +86 27 87557815; fax: +86 27 87540249.
E-mail address: [email protected] (H.-c. Zhou). 0022-4073/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2009.03.007
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
Nomenclature Er Emax I i j M0 n^ N0 N1 Ni q Rid Rsd r s s^ S T
difference in reflectance by specular and diffuse boundaries maximum relative difference as defined in Eq. (9) radiation intensity, W/(m2sr) discrete grids index discrete angles index number of total discrete grids normal direction of boundary number of total discrete angles number of energy bundles iterative times heat flux, W/(m2) an intermediate value as defined in Eq. 10. DRESOR value, 1/m3, or 1/m2 position vector, m distance, m unit vector along a given direction radiative source function, W/(m2sr) temperature, K
y
ss r o o1,o2 t e O F
1073
directional angles scattering coefficient, 1/m reflectivity of boundary scattering albedo parameters of the scattering phase function optical thickness emissivity of boundary angular space, sr scattering phase function
Superscripts n d s
iterative number diffuse reflection specular reflection
Subscripts w b 1 2
wall blackbody left boundary right boundary
Greek letters
b
k
attenuation coefficient, 1/m absorption coefficient, 1/m
absorbing, emitting, anisotropically scattering media, allowing for spatially varying absorption and scattering coefficients. In these papers the boundaries considered were all transparent or diffusely reflecting. But in applications some surfaces are smooth and cannot be treated as diffuse reflectors [6], for example, the ceramic surfaces in combustion chambers of highspeed aircraft; also the effect of reflecting modes of boundaries on the layer equilibrium temperatures of an absorbing–emitting layer exposed to incident radiation was analyzed. An integral theory was presented by Machali [7,8] to consider the isotropy and anisotropy of scattering and the effects of diffuse and specular reflectivity of the bounding surfaces in a plane-parallel medium. Maruyama used the radiation element method [9] to solve the problem with a specular-diffuse boundary, and both collimated and diffuse incident irradiation were analyzed. But for a planeparallel, purely linearly anisotropic scattering medium with emitting, absorbing, and specular-diffuse boundaries the reflectance got [8,10] showed obvious differences; and the exact solution for this case has not been confirmed and needs further study. Furthermore, most of these methods cannot efficiently calculate radiative intensity with high directional resolution, which provides more information for inverse radiative transfer problems. This plays a key role in radiative image processing in some industrial applications, for example, monitoring of temperature distributions in industrial combustion furnaces [11]. A new method called the DRESOR method, based on the Monte Carlo method, can directly calculate intensity with high directional resolution. It has been applied to the isotropic scattering media [12,13], and then extended to solve the problem with anisotropic scattering media [14]. The results calculated by the DRESOR method (it is hereinafter called the Basic-DRESOR method so as to distinguish it from the present DRESOR method described in this paper) agreed well with those calculated by some other methods and the integral formulation of RTE. However, the relative errors for the intensity got by the Basic-DRESOR method are several percent, and it is hard to reduce them significantly, due to the inherent feature belonging to all the Monte Carlo methods. Iterative techniques were usually used to improve the accuracy of numerical methods. The one-dimensional radiative transfer equation can be transformed into an integral equation for the source function and the resulting equation can be solved by iteration [15]. Analytical studies on the convergence of iterative schemes [16] for the radiative transfer equation show that the convergence can be very slow if the single scattering albedo is close to unity. Sutton and Ozisik [17] improved the convergence significantly by using the P1 approximation of the spherical harmonics method for the initial guess, and focused on the dimensionless radiative heat flux (reflectance and transmittance). The convergence criterion adopted in [17] was that the difference of the source term between two successive iterations was less than a specific value, such as 1.0 103. If directional quantities such as the angular distribution of radiation intensity are needed, it might be desirable to use a larger number of quadrature points [17].
ARTICLE IN PRESS 1074
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
In this paper, the iterative algorithm is directly applied to the integral form of the radiative transfer equation, and the intensity got by the Basic-DRESOR method is used as an initial guess. The description of the Iterative-DRESOR method, together with the treatment method for the specularly reflecting boundary, which has not been clarified before in the Basic-DRESOR method, will be described in detail in Section 2. In Section 3, the method was validated by the exact reflectance in literature for an anisotropic scattering medium. The convergence of the Iterative-DRESOR method will be demonstrated; and the radiative flux for isotropic and linearly anisotropic scattering media got by the present method will be compared with those in literature. Furthermore, the effect of the reflecting modes, especially the specular boundary, on the intensity and the radiative flux will be examined thoroughly. Finally some conclusions will be drawn. 2. The Iterative DRESOR method 2.1. The Basic-DRESOR method The radiative transfer equation (RTE) for an emitting, absorbing and anisotropic scattering medium is [18] Z ss Iðr; s^ i ÞFðs^ i ; s^ Þ dOi s^ rIðr; s^ Þ ¼ kIb ðr; s^ Þ bIðr; s^ Þ þ 4p 4p and its integral form for the medium is given [18] as Z s Z s Z Iðr; s^ Þ ¼ Iw ðr; s^ Þ exp b ds00 þ Sðr0 ; s^ Þ exp 0
0
s0
b ds00 b ds0
(1)
(2)
0
where s ¼ |rrw|, s0 ¼ |rr0 |, and the direction of integration is switched to go along s00 (from point r toward the wall). The source function Sðr0 ; s^ Þ is [18] Z o Iðr0 ; s^ i ÞFðs^ i ; s^ Þ dOi (3) Sðr0 ; s^ Þ ¼ ð1 oÞIb ðr0 Þ þ 4p 4p generally, for opaque surfaces with arbitrary surface properties, the boundary condition is given [18] as Z rðrw ; s^ 0 ; s^ ÞIðrw ; s^ 0 Þjn^ s^ 0 j dO0 Iw ðrw ; s^ Þ ¼ ðrw ÞIb ðrw Þ þ ^ s^ 0 o0 n
(4)
In the previous works of the authors [11–14], the DRESOR values, Rsd , were introduced to quantitatively account the integral terms in Eqs. (3) and (4), and the method was called the DRESOR method. Then, the formula for Iðr; s^ Þ by the BasicDRESOR method was given as Z 1 pðrw ÞIb ðrw Þ þ Rsd ðr0w ; rw ; s^ Þ½pðr0w ÞIb ðr0w Þ dA0 Iðr; s^ Þ ¼ p W Z s Z Rsd ðr00 ; rw ; s^ Þ½4pbð1 oÞIb ðr00 Þ dV 00 exp b ds00 þ 0 V Z Z s 1 s 0 0 0 ^ 0 4pbð1 oÞIb ðr Þ þ Rd ðrw ; r ; sÞ½pðrw ÞIb ðr0w Þ dA0 þ 0 4p W Z s0 Z þ Rsd ðr00 ; r0 ; s^ Þ½4pbð1 oÞIb ðr00 Þ dV 00 exp b ds00 ds0 (5) V
0
From Eq. (5) it is obvious that, for an emitting, absorbing and scattering medium bounded by walls with arbitrary 0 properties, once all the distributions of kðrÞ, ss ðrÞ, ðrw Þ, Fðs^ i ; s^ Þ and rðrw ; s^ ; s^ Þ are given, all the DRESOR values can be calculated. Then Iðr; s^ Þ, the angular distribution of radiative intensity at any point r, can be expressed as a function of the three-dimensional distribution of Ib ðr0 Þ inside the medium and the two-dimensional distribution of Ib ðr0w Þ on the wall surfaces. For the calculation of the DRESOR values, the tracking procedure of every bundle is based on the path-length method following the general procedure of the Radiative Energy Absorption Distribution (READ) method, a Monte Carlo method [19]. In tracking the energy bundles, the ratio of energy scattered by the medium and/or reflected by the boundary surfaces is recorded in the DRESOR method, whereas the ratio of energy absorbed by the medium and boundary is recorded in the READ method. Once the intensity with high directional resolution is obtained, it is easy to deduce all the associated quantities, for example, the radiative flux. The description for the Basic-DRESOR method can be found in the literature [11–14], and it will not be repeated here again. 2.2. The Iterative-DRESOR method In the Iterative-DRESOR method, the intensity calculated using the Basic-DRESOR method is taken as an initial guess for an iterative algorithm, marked as Ið0Þ ðr; s^ Þ, and it will be substituted into Eqs. (2)–(4), then a new intensity with the same high directional resolution will be obtained.
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1075
Let IðnÞ ðr; s^ Þ be the intensity after the nth iteration, n ¼ 1,2,y. At first, it is substituted into the RHS of Eqs. (3) and (4) to ðrw ; s^ Þ as below: obtain an updated source function Sðnþ1Þ ðr; s^ Þ and boundary intensity Iðnþ1Þ w Z o IðnÞ ðr0 ; s^ i ÞFðs^ i ; s^ Þ dOi (6) Sðnþ1Þ ðr0 ; s^ Þ ¼ ð1 oÞIb ðr0 Þ þ 4p 4p Iðnþ1Þ ðrw ; s^ Þ ¼ ðrw ÞIb ðrw Þ þ w
Z ^ s^ 0 o0 n
rðrw ; s^ 0 ; s^ ÞIðnÞ ðrw ; s^ 0 Þjn^ s^ 0 j dO0
(7)
Then an updated intensity Iðnþ1Þ ðr; s^ Þ can be calculated by substituting Sðnþ1Þ ðr; s^ Þ and Iðnþ1Þ ðrw ; s^ Þ into the RHS of Eq. (2). w Z ^ Iðnþ1Þ ðr; s^ Þ ¼ Iðnþ1Þ ðr; s Þ exp w
s
b ds00 þ 0
Z
s
Z Sðnþ1Þ ðr0 ; s^ Þ exp
0
s0
b ds00 b ds0
(8)
0
The iteration will end when the maximum updating ratio for the intensity is smaller than a pre-determined value, such as 1.0 106. Emax ¼ maxfjIðnþ1Þ ðr; s^ Þ IðnÞ ðr; s^ Þj=IðnÞ ðr; s^ Þgo1:0 106
(9)
r;s^
If the condition of Eq. (9) is satisfied, it can be confirmed that the iterative algorithm is converged. This means not only that the solution of the RTE with high directional resolution is obtained, but also the accuracy of the solution reaches a very high level with relative errors less than 1.0 106. Therefore, the solution for the intensity can be regarded as a numerically exact one.
2.3. Treatment for the specularly reflecting boundary For the specular reflector the energy bundle will be reflected only at a specific angle equal to the incident angle on the opposite side. This causes the treatment for the specular reflecting boundary to be different from that for the diffuse reflecting boundary as described [14] and as discussed below. Before calculating the DRESOR values Rsd ðr; rw ; s^ j Þ for the specular reflecting boundary, Rid ðr; rw ; s^ i Þ, an intermediate variable is used to record the ratio of the energy incident from a unit solid angle around the direction s^ i to a unit area around rw to that emitted from a unit area or a unit volume around the point r. As the calculation begins, it is set to zero. N1 energy bundles are tracked, and the initial energy is set as E0 ¼ 1.0, experiencing gradual reduction as it is absorbed along its traveling path. During the tracking process, it is renewed by Rid;new ðr; rw ; s^ i Þ ¼ Rid;old ðr; rw ; s^ i Þ þ
E0 ðs^ i Þ N1
(10)
where E0 ðs^ i Þ means the residual energy possessed by the bundle arriving at the boundary from direction s^ i . For the specular reflecting boundary, the DRESOR values can be calculated as Rsd ðr; rw ; s^ j Þ ¼ rsrw Rid ðr; rw ; s^ M0 jþ1 Þ p=ðdO cos yj Þ
(11)
where dO is the solid angle for the discrete angle. For the diffuse reflecting boundary Rsd ðr; rw ; s^ j Þ ¼
M0 X
rdrw Rid ðr; rw ; s^ i Þ
(12)
i¼1
The other parts of treatment for the specular reflecting boundary are the same as those for the diffuse reflecting boundary.
Table 1 Calculation cases and conditions. Case Parameters studied 1
Reflectance q~ ¼ 1 qð0Þ=ðsT 41 Þ
2
Heat flux q~ (in units of 2e1sT41)
3
Radiative intensity, dimensionless reflectance and transmittance
Boundary 1 T1 ¼ 1000 K; e1 ¼ 1.0 T1 ¼ 1000 K; specularly reflecting T1 ¼ 1000 K; diffusely or specularly reflecting
Medium
t ¼ 1.0–10.0; o ¼ 1.0; F(y) ¼ 1+o1P1(cos y)+o2P2(cos y) t ¼ 0.01–5.0; o ¼ 1.0; F(y) ¼ 1+x cos y, x ¼ 1, 0, 1
t ¼ 0.01–5.0; o ¼ 0.0, 0.5, or 1.0; F(y) ¼ 1+x cos y, x ¼ 1, 0, 1
Boundary 2 T2 ¼ 0 K; e2 ¼ 1.0 T2 ¼ 500 K; e2 ¼ 0.8; diffusely reflecting T2 ¼ 0 K; diffusely or specularly reflecting
References Busbridge and Orchard [1,2] Machali [7,8]; Tan et al. [10].
ARTICLE IN PRESS 1076
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
3. Results and discussion 3.1. Objective, cases and conditions of calculation The Iterative DRESOR method is applied to plane-parallel and isotropic/anisotropic scattering media. Three cases are studied in the paper, as listed in Table 1. In Case 1, two boundaries are black with ¼ 1:0 and with two different temperatures T1 ¼ 1000 K, T2 ¼ 0. The medium is pure scattering with o ¼ 1:0. Nonlinearly scattering phase function F(y) ¼ 1+o1P1(cos y)+o2P2(cos y) is considered. The reflectance will be calculated. This case is used to validate the accuracy of the Iterative-DRESOR method. In Case 2, the medium is also pure scattering, and two boundaries are diffusely emitting, specularly reflecting and diffusely reflecting with temperatures T1 ¼ 1000 K, T2 ¼ 500 K. The linearly anisotropic scattering phase function (F(y) ¼ 1+x cos y , x ¼ 1, 0, 1) is considered. When x ¼ 0 the phase function turns into isotropic scattering. The emissivity of Boundary 2 is 0.8 and it is diffuse reflecting with rd2 ¼ 0:2. The emissivity of boundary 1 is variable, and it has specular reflectance rs1 ¼ 1 1 . In this case, there is an exact solution of the heat flux for the isotropic scattering, which can be used to verify the present method. There is not an exact solution for the anisotropic scattering, and some uncertainty in the solution will be clarified.
Table 2 The time taken to calculate the intensity with the phase function F1 for different optical thickness in Case 1.
t0
N0
1 2 3 4 5 6 7 8 9 10
M0
100 100 200 200 300 300 400 400 500 500
Calculation time by BasicDRESOR method (s)
180 180 180 180 180 180 180 180 180 180
32 50 99 107 194 210 311 326 438 453
Iterative number and calculation time by Iterative-DRESOR method after convergence Iterative number
Calculation time (s)
22 41 61 88 117 150 178 226 263 309
14 27 126 183 498 639 1291 1689 2903 3421
Table 3 Reflectance for the phase function F(y) ¼ 1+o1P1(cos y)+o2P2(cos y) for various optical thicknesses in Case 1.
t0
1 2 3 4 5 6 7 8 9 10
F0 (o1 ¼ 0, o2 ¼ 0) Exact
DRESOR
0.4466 0.6099 0.6984 0.7540 0.7923 0.8203 0.8417 0.8585 0.8721 0.8833
0.43646 0.58316 0.67332 0.72334 0.76488 0.78600 0.80974 0.82406 0.83922 0.84600
F1 (o1 ¼ 0, o2 ¼ 0.5) Present work 2.270% 4.384% 3.591% 4.066% 3.461% 4.181% 3.797% 4.012% 3.770% 4.223%
0.44659 0.60995 0.69838 0.75408 0.79241 0.82041 0.84177 0.85859 0.87219 0.88341
0.002% 0.008% 0.003% 0.011% 0.014% 0.013% 0.008% 0.010% 0.010% 0.012%
F2 (o1 ¼ 1.0, o2 ¼ 0.5) 1 2 3 4 5 6 7 8 9 10
0.3580 0.5157 0.6104 0.6739 0.7197 0.7541 0.7810 0.8026 0.8204 0.8352
0.34864 0.49419 0.58729 0.64481 0.69381 0.72562 0.75477 0.77006 0.79215 0.80330
Exact
DRESOR
0.4468 0.6101 0.6985 0.7541 0.7924 0.8204 0.8417 0.8585 0.8721 0.8833
0.43639 0.58686 0.68043 0.72624 0.76270 0.78789 0.81242 0.82399 0.84087 0.85132
Present work 2.330% 3.516% 2.587% 3.694% 3.748% 3.963% 3.479% 4.020% 3.581% 3.621%
0.44682 0.61013 0.69851 0.75417 0.79247 0.82046 0.84180 0.85862 0.87221 0.88343
0.004% 0.005% 0.001% 0.009% 0.009% 0.007% 0.012% 0.014% 0.013% 0.015%
3.579% 4.833% 3.781% 4.187% 4.202% 4.281% 3.998% 4.201% 4.008% 4.065%
0.30202 0.44900 0.54371 0.61046 0.66015 0.69860 0.72924 0.75423 0.77500 0.79254
0.007% 0.000% 0.002% 0.010% 0.008% 0.014% 0.005% 0.017% 0.013% 0.018%
F3 (o1 ¼ 1.5, o2 ¼ 0.5) 2.615% 4.171% 3.786% 4.317% 3.597% 3.777% 3.359% 4.054% 3.443% 3.819
0.35803 0.51572 0.61039 0.67398 0.71972 0.75420 0.78114 0.80276 0.82050 0.83532
0.008% 0.004% 0.002% 0.012% 0.003% 0.013% 0.018% 0.020% 0.012% 0.014%
0.3020 0.4490 0.5437 0.6104 0.6601 0.6985 0.7292 0.7541 0.7749 0.7924
0.29119 0.42730 0.52314 0.58484 0.63236 0.66860 0.70005 0.72242 0.74384 0.76019
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1077
In Case 3, the boundary conditions are the same as those in Case 2 except the temperature of Boundary 2 is set to 0, that is T2 ¼ 0.0. The medium is absorbing, scattering, and non-emitting with scattering phase function F(y) ¼ 1+x cos y , x ¼ 1, 0, 1. This case is used to compare the effects of different reflecting modes on the radiative intensity, the transmittance and the reflectance.
0.1 0.01
Emax
1E-3 1E-4 1E-5 1E-6 0
2
4
6
8
10
12
14
16
18
20
22
Ni 0.1 0.01
Emax
1E-3 1E-4 1E-5 1E-6 0
30
60
90
120 150 180 210 240 270 300 Ni
Fig. 1. The maximum error in intensity with phase function F1, (a) t0 ¼ 1.0 and (b) t0 ¼ 10.0 in Case 1.
Table 4 The time taken to calculate the phase function F(y) ¼ 1cos y with rs2 ¼ rd1 ¼ 0, rd2 ¼ 0:2, e2 ¼ 0.8 and different e1 in Case 2.
t0
N0
M0 Calculation time by Basic-DRESOR method (s) Iterative numbers and calculation time by Iterative-DRESOR method after convergence Iterative numbers
e1 ¼ 0.2 0.01 10 180 21 0.1 20 180 46 0.5 50 180 131 1.0 100 180 280 2.0 100 180 326 5.0 300 180 1081
Calculation time (s)
e1 ¼ 0.7
e1 ¼ 1.0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
17 31 89 190 200 678
2 5 16 43 53 270
8 14 29 46 89 233
6 10 20 30 57 153
4 8 16 25 49 139
0.3 1 7 30 60 1029
0.3 0.8 5 20 39 711
0.2 0.6 4 16 32 630
ARTICLE IN PRESS 1078
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
3.2. Validation of the Iterative-DRESOR method The Iterative-DRESOR method is validated in Case 1. Participating media with different optical thicknesses are subjected to a diffuse irradiation flux from Boundary 1. The medium is divided into N0 discrete grids and M0 discrete angles in the angle range of [0, p]. As the optical thickness increases, N0 will also increases to maintain accurate calculation. The number of energy bundles used in the DRESOR method is 100,000. The discrete number and times consumed with different optical thicknesses are listed in Table 2. The reflectance q~ ¼ 1 qð0Þ=ðsT 4s1 Þ is compared with the exact solution obtained by Busbridge and Orchard [1,2]; and the calculation results from the present method are compared with the exact solutions, as listed in Table 3. The comparisons are made for four different phase functions with different parameters o1 and o2 , and F0 is isotropic scattering and F1 is Rayleigh scattering phase functions. The calculations in this paper are all carried out on a Pentium IV, 3.0 GHz PC. As shown in Table 2, the new method does not consume much more time compared with the Basic-DRESOR method, especially when the optical thickness is not large, where only a few seconds are needed for the iterations. The iteration number increases as the optical thickness increasing, and the time consumed increases rapidly. But as shown in Fig. 1, Emax decreases continuously as the iteration time increases. The curve displays as an exponential curve after several iterations. As the optical thickness is 10.0, the calculation time by the Iterative-DRESOR method is 3421 s, being 7.55 times than 453 s by the Basic-DRESOR method. Compared with the reduction of calculation errors from the level of 1:0 102 to that of 1:0 106 , the consumption in time is a worthy cost. If the accuracy of the results is not required to be very high, the time taken can be reduced. For example, if the accuracy needed is 1:0 103 in Fig. 1(b), the iteration number needed is only 12 and the time consumed for iteration is only 136 s. The results mentioned above are done for a pure scattering medium. If the scattering albedo decreases, the iteration number will decrease quickly, and the time used can be further reduced. The exact reflectance in Case 1 from literature [1,2] and those calculated by the Basic-DRESOR method and the IterativeDRESOR method and their errors are listed in Table 3 for the four phase functions. Compared with the exact solutions, the accuracy of the Iterative-DRESOR method is significantly improved. The errors are all less than 0.02%. The comparison shows that the present method is very accurate without much increase in computing time.
Table 5 The dimensionless heat flux for isotropic scattering with rs2 ¼ rd1 ¼ 0, rd2 ¼ 0.2, e2 ¼ 0.8 and different e1 in Case 2.
t0
0.01 0.1 0.5 1.0 2.0 5.0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
Machali [8]
Present work
Machali [8]
Present work
Machali [8]
Present work
Heat flux
Heat flux
Relative errors (%)
Heat flux
Heat flux
Relative errors (%)
Heat flux
Heat flux
Relative errors (%)
0.44559 0.43851 0.41229 0.38554 0.34261 0.25772
0.445607 0.438524 0.412300 0.385547 0.342593 0.257644
0.004 0.003 0.002 0.002 0.005 0.029
0.39661 0.37798 0.31827 0.26854 0.20591 0.12165
0.396620 0.377987 0.318273 0.268539 0.205887 0.121600
0.003 0.002 0.001 0.000 0.011 0.041
0.37208 0.34928 0.28067 0.22788 0.16660 0.09254
0.372087 0.349286 0.280674 0.227880 0.166584 0.092497
0.001 0.002 0.001 0.000 0.010 0.046
Table 6 The dimensionless heat fluxes for linear anisotropic scattering with rs2 ¼ rd1 ¼ 0, rd2 ¼ 0.2, e2 ¼ 0.8 and different e1 in Case 2.
t0
e1 ¼ 0.2
e1 ¼ 0.7
e1 ¼ 1.0
Machali [8]
Tan et al. [10]
DRESOR
Machali [8]
Tan et al. [10]
DRESOR
Machali [8]
Tan et al. [10]
DRESOR
0.44580 0.44055 0.42084 0.39918 0.36008 0.27062
0.445854 0.440978 0.422416 0.402770 0.369657 0.297461
0.445818 0.440585 0.421571 0.402085 0.369615 0.298730
0.39720 0.38335 0.33683 0.29299 0.22909 0.13191
0.397343 0.384319 0.339727 0.299160 0.242866 0.155902
0.397208 0.383397 0.338376 0.298464 0.243298 0.157344
0.37282 0.35586 0.30168 0.25370 0.18888 0.10124
0.372997 0.356935 0.304539 0.259759 0.201875 0.121452
0.372826 0.355916 0.303382 0.259410 0.202588 0.122814
x ¼ 1.0 0.01 0.44538 0.1 0.43649 0.5 0.40443 1.0 0.37421 2.0 0.33032 5.0 0.25047
0.445319 0.436010 0.402241 0.369110 0.318571 0.226868
0.445395 0.436483 0.403428 0.370316 0.319252 0.226488
0.39602 0.37277 0.30234 0.25005 0.19105 0.11614
0.395867 0.371612 0.298426 0.242584 0.177949 0.099421
0.396034 0.372727 0.300424 0.244068 0.178447 0.099085
0.37134 0.34295 0.26315 0.20903 0.15274 0.08791
0.371146 0.341579 0.259136 0.201901 0.141122 0.074488
0.371350 0.342898 0.261129 0.203184 0.141445 0.074181
x ¼ 1.0 0.01 0.1 0.5 1.0 2.0 5.0
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1079
3.3. Radiative heat transfer with specular-diffuse boundaries Consider a pure scattering plane-parallel medium in Case 2 with isotropic scattering. For different optical thicknesses, the discrete number of the medium, the calculation time by the Basic-DRESOR and the Iterative-DRESOR methods, and the iterative numbers, are all listed in Table 4. It can be seen that when the optical thickness is very small (o1.0), the time taken for the iterations in the Iterative-DRESOR method is negligible compared with that taken for the Basic-DRESOR method, but the accuracy is considerably improved. The dimensionless heat flux q~ (in units of 2e1sT41) with isotropic scattering is given in Table 5. Results calculated by the Iterative-DRESOR method agree well with those given by Machali [7] for the isotropic scattering, and the biggest error is less than 0.05%. The dimensionless heat flux from literature [8,10] and obtained by the Iterative-DRESOR method for two linear anisotropic scattering phase functions with different e1 are given in Table 6. When the optical thickness is bigger than 1.0, the results given by Machali [8] are different from those calculated by Tan et al. [10] and the present method. The results obtained by the Iterative-DRESOR method, a numerically exact method with relative errors of 1.0 106, agree very well with those of Tan et al. [10]. So, the method in [8] should be carefully checked when it is applied to large optical thickness media with anisotropic scattering. 3.4. Effect of specular reflecting boundary on the intensity The effects of different reflecting modes of boundaries on radiative intensity with high directional resolution are examined in Case 3, where N0 ¼ 100, and M0 ¼ 180. The intensity with high directional resolution calculated by the Iterative-DRESOR method and their relative errors are all shown in Figs. 2–4. i, j denote the discrete numbers of grid and angle, respectively. In these cases the medium is pure scattering with linear forward scattering and the optical thickness
180 160 140 2.2E3
3.1E3
120
2.8E3 2.5E3 3.5E3 3.8E3 4.1E3 4.8E3 4.4E3 5.1E3 5.4E3 6.1E3 5.7E3
j
100 80 60
6.4E3
40 20 0
10
20
30
40
50
60
70
80
90 100
i
1.0x10-6 9.0x10-7
error
8.0x10-7 7.0x10-7 6.0x10-7 5.0x10-7 4.0x10-7 3.0x10 10 20 30 40 50 60 70 80 i 90 100
180 160 140 120 100 80 60
j
-7
20
40
Fig. 2. (a) Intensity and (b) its relative errors for linear forward scattering with rd1 ¼ rs2 ¼ 0, rs1 ¼ 0:8 and rd2 ¼ 0:2 in Case 3.
ARTICLE IN PRESS 1080
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
t ¼ 1.0. The emissivity of the two boundaries are e1 ¼ 0.2, and e2 ¼ 0.8, respectively. The main difference lies in the reflecting modes of the two boundaries. In Fig. 2, Boundary 1 is a specularly reflecting surface and Boundary 2 is a diffusely reflecting surface with rs1 ¼ 0:8, rd2 ¼ 0:2. In Fig. 3, the two boundaries are diffusely reflecting surfaces with rd1 ¼ 0:8, rd2 ¼ 0:2. In Fig. 4, Boundary 1 is diffusely reflecting and Boundary 2 is a specularly reflecting surface with rd1 ¼ 0:8, rs2 ¼ 0:2. From these figures, it can be seen firstly that the intensity obtained by the Iterative-DRESOR method is very accurate, and the relative errors are less than 1.0 106. Secondly, the difference in intensity caused by different reflecting modes is also obvious, especially in the zone near the specular reflecting boundary. In order to show the difference in detail, the directional intensity at the two boundaries is plotted in Fig. 5(a) and (b). As shown in Fig. 5(a), the intensity at Boundary 1 around 90 degrees in the forward directions is larger with specular reflecting than that with a diffuse reflecting boundary. In contrast, as shown in Fig. 5(b), the intensity at Boundary 2 in the backward direction around 901 with specular reflecting boundary is smaller than that with a diffuse reflecting boundary. This varying distribution of directional intensity also appears in the angle range of [p/2, p] if Boundary 2 is specular reflecting; but if Boundary 2 is diffuse reflecting the directional intensity in the range will become more uniform. 3.5. Effect of specular reflecting boundary surfaces on the flux Researchers are usually more concerned about radiative flux. The effects of specular reflecting boundary on the reflectance and transmittance in Case 3 are shown in Table 7, where the medium is linearly forward scattering with o ¼ 0.5 and t ¼ 1.0. It can be seen that the reflecting mode of Boundary 1, the emitting boundary, has little effect on the reflectance and transmittance. The differences of reflectance and transmittance are all less than 1.52% with different reflecting mode of Boundary 1. The reflecting mode of Boundary 2, the non-emitting boundary, also has little effect on the transmittance, and the differences are all less than 1.0%. It should be stated that the effect of different reflecting modes of Boundary 2 on the
180 160 140 3.3E3
2.0E3
2.8E3
120
2.5E32.3E3
3.0E3
100
3.5E3 3.8E3 4.3E34.0E3 4.5E3 5.0E3
j
4.8E3
80 60
5.8E3
6.3E3
40
5.5E3 5.3E3
6.0E3
20 0
10
20
30
40
50
60
70
80
90 100
i
1.0x10-6 9.0x10-7
error
8.0x10-7 7.0x10-7 6.0x10-7 5.0x10-7
-7
3.0x10
10 20 30 40 50 60 70 80 i 90 100 40 20
180 160 140 120 100 80 60
j
4.0x10-7
Fig. 3. (a) Intensity and (b) its relative errors for linear forward scattering with rs1 ¼ rs2 ¼ 0, rd1 ¼ 0.8 and rd2 ¼ 0:2 in Case 3.
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1081
reflectance is obvious, and the differences are all bigger than 4% when r2 ¼ 0.8 and bigger than 2% when r2 ¼ 0.2. If the boundary’s reflectivity is equal to zero, there will be no reflectance, and consequently, no difference. The difference of reflectance between specular reflecting and diffuse reflecting of Boundary 2 is hereinafter expressed as Er for brevity. As shown in Table 7, as the reflectivity of Boundary 2 increases, Er increases. In addition, there are some other factors affecting the difference, such as the optical thickness, the scattering albedo, and the scattering phase function of the medium, which will be given below. Er variation with the optical thickness is shown in Fig. 6 with scattering albedo o ¼ 0.5, and the emissivities of the two boundaries being e1 ¼ 1.0, and e2 ¼ 0.2. As shown in Fig. 6, the difference increases as the optical thickness increases, but once the optical thickness has reached a specific value, the difference decreases as the optical thickness increases. As shown in Fig. 5(b), the directional intensity around 01 is larger than that in other directions, and less energy will be absorbed after it is reflected by the specular reflecting Boundary 2, leading to a larger reflectance than that of a diffuse reflecting Boundary 2. If the optical thickness of the medium increases further, less energy reflected by Boundary 2 can reach Boundary 1, and Er becomes smaller and smaller. The variation of Er with o ¼ 1.0, 0.5, and 0.0 is listed in Table 8, where the medium is isotropically scattering with t ¼ 1.0, e1 ¼ 1.0, and e2 ¼ 0.2. Er increases from 0.08% to 25.25% as the scattering albedo decreases from 1.0 to 0.0. When the scattering albedo becomes larger, less energy emitted from Boundary 1 can reach Boundary 2, and Er will decrease. If the medium is pure scattering, Er is negligible. The influence of phase function on Er is shown in Table 9. It can be seen that the effect of the phase function on the reflectance cannot be ignored, not only for the same reflecting condition of the boundary, but also between different reflecting modes of the boundary. It has been stated that Monte Carlo method can be used to solve tedious and difficult problems in radiative transfer, and because the method is known to produce accurate results if enough sample packets are used, the Monte Carlo method is
180 160 3E3
140
2.3E32.0E3
120
3.0E3
100 j 80 6.3E3
2.5E3
4.0E33.8E33.5E3 4.3E3 4.8E3
4.5E3
60
2.8E3
5.8E35.5E3 6.0E3
5.3E3
5.0E3
40 20 0
10
20
30
40
50
60
70
80
90 100
i
1.0x10-6 9.0x10-7
error
8.0x10-7 7.0x10-7 6.0x10-7 5.0x10-7 4.0x10-7 3.0x10 10 20 30 40 50 60 70 80 i 90 100
180 160 140 120 100 80 60
j
-7
20
40
Fig. 4. (a) Intensity and (b) its relative errors for linear forward scattering with rs1 ¼ rd2 ¼ 0, rd1 ¼ 0:8 and rs2 ¼ 0:2 in Case 3.
ARTICLE IN PRESS 1082
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
often used to validate the results of other methods [20]. The Iterative-DRESOR method, as a development of Monte Carlo method, will strengthen the capability of the Monte Carlo method in this function. In order to expand this methodology to more complex realistic cases, especially the spectral problems, such as FSK (full-spectrum k-distribution) method [21] and the LBL (line-by-line) method [22,23] can be adopted in the DRESOR method in the future work. 4. Conclusion The Iterative-DRESOR method is an iterative algorithm applied to the integral form of the radiative transfer equation using an initial guess of intensity with high directional resolution calculated by the Basic-DRESOR method. The treatment of the specular reflecting boundary in tracking the energy bundles is proposed, which has not been clarified before in the Basic-DRESOR method. For a plane-parallel, anisotropic scattering medium, bounded with black boundaries with different specularly reflecting
diffusely reflecting 90
Intensity at boundary1
8000
120
60
7000 6000 30
150
5000 4000 3000 2000
180
0 specularly reflecting
diffusely reflecting 90
Intensity at boundary2
5000
120
60
4000 3000
30
150
2000 1000 0 180
0
Fig. 5. (a) Intensity at Boundary 1 with different reflecting modes of Boundary 1, and (b) at Boundary 2 with different reflecting modes of Boundary 2 in Case 3.
Table 7 The effect of different reflecting modes of the two boundaries on dimensionless radiative flux with linearly forward scattering medium t ¼ 1.0, o ¼ 0.5 in Case 3. Conditions
r1 ¼ 0.0, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.8 r1 ¼ 0.8, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.0
r1 ¼ 0.0, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.8 r1 ¼ 0.8, r2 ¼ 0.8 r1 ¼ 0.2, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.2 r1 ¼ 0.8, r2 ¼ 0.0 a
Two diffuse boundaries
Specular Boundary 1, diffuse Boundary 2 a
Reflectance
Reflectance
Relative Error (%)
0.19480 0.16217 0.04611 0.09688 0.02608 0.02044
0.16215 0.04614 0.09706 0.02629 0.02075
0.01 0.07 0.19 0.81 1.52
Transmittance
Transmittance
Relative errora (%)
0.07359 0.07654 0.08718 0.28407 0.30634 0.36801
0.07657 0.08705 0.28348 0.30394 0.36447
0.04 0.15 0.21 0.78 0.96
Relative error to the corresponding data for two diffuse boundaries.
Diffuse Boundary 1, specular Boundary 2 Reflectance 0.20344 0.16965 0.04855 0.09900 0.02669
Relative Errora (%) 4.44 4.61 5.29 2.19 2.34
Transmittance
Relative errora (%)
0.07289 0.07599 0.08707 0.28348 0.30623
0.95 0.72 0.13 0.21 0.04
ARTICLE IN PRESS Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
1083
5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0
1
2
3
4
5
Fig. 6. The effect of different optical thickness with linear forward scattering o ¼ 0.5, e1 ¼ 1.0, and e2 ¼ 0.2 in Case 3.
Table 8 The effect of different scattering albedo on reflectance with isotropic scattering medium t ¼ 1.0, e1 ¼ 1.0, and e2 ¼ 0.2 in Case 3.
Reflectance with diffuse Boundary 2 Reflectance with specular Boundary 2 Er
o ¼ 1.0
o ¼ 0.5
o ¼ 0.0
0.82785 0.82847 0.08%
0.21846 0.22629 3.58%
0.03846 0.04817 25.25%
Table 9 The effect of different phase functions on reflectance with t ¼ 1.0, o ¼ 0.5, e1 ¼ 1.0, and e2 ¼ 0.2 in Case 3.
Reflectance with diffuse Boundary 2 Reflectance with specular Boundary 2 Er
Forward scattering
Isotropic scattering
Backward scattering
0.19480 0.20344 4.44%
0.21846 0.22629 3.58%
0.24016 0.24728 2.96%
temperatures, exact solutions for the reflectance are used to validate the accuracy of the new method, and the convergence of the new method is demonstrated. For media with isotropic and linearly anisotropic scattering, the radiative flux obtained by the present method is compared with those in the literature. Furthermore, the effects of the different reflecting boundaries, especially the specular reflecting boundary, on the radiative intensity and flux are examined thoroughly, and this paper focuses on the difference of reflectance between specular and diffuse reflecting boundaries. The main conclusions that are drawn are listed below. (1) The Iterative-DRESOR method can give intensity with high directional resolution at a high accuracy level with relative errors less than 1.0 106 with acceptable computation costs even for an optical thickness of 10.0. For optical thicknesses less than 1.0, the increase in the computation cost is negligible, but the accuracy is improved significantly. (2) The reflectance in a pure, anisotropic scattering medium with black boundaries obtained by the present method is validated by the exact data in the literature, and agreement at a high level is confirmed. (3) For a pure, anisotropic scattering medium with specular and diffuse reflecting boundaries, the numerically-exact results for the reflectance obtained by the present method clarified some uncertainty over the correct solution in the literature. (4) The directional distributions of intensity are obviously influenced by the reflecting modes of the boundaries, especially in the zone near the boundary, which may reflect specularly or diffusely. (5) The reflecting mode of the emitting boundary has little effect on the transmittance or reflectance. The reflecting mode of the non-emitting boundary also has little effect on the transmittance, but it obviously influences the reflectance. (6) The difference between the reflectance by the specular and diffuse boundaries increases at first, and then decreases, as the optical thickness of the medium increases. The difference will decrease as the scattering albedo of the medium
ARTICLE IN PRESS 1084
Z.-f. Huang et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 110 (2009) 1072–1084
increases, and it is negligible when the medium is pure scattering. Also the effect of the scattering phase function of the medium on the difference cannot be ignored. From the conclusions above, we can expect that the Iterative-DRESOR method can strengthen the capability of the Monte Carlo method to produce accurate results and to validate the results of other methods.
Acknowledgment The study has been supported by the National Natural Science Foundation of China (No. 50636010, No. 50721005). References [1] Busbridge IW, Orchard SE. Reflection and transmission of light by a thick atmosphere according to phase function: 1+x cos y. The Astrophysical Journal 1967;149:655–64. [2] Orchard SE. Reflection and transmission of light by thick atmospheres of pure scatterers with a phase function: 1+w1P1(cos y)+w2P2(cos y). The Astrophysical Journal 1967;149:665–74. [3] Altac Z. The SKN approximation for solving radiative transfer problems in absorbing, emitting, and linearly anisotropically scattering plane-parallel medium: part 2. Journal of Heat Transfer 2002;124(4):685–95. [4] Sarma D, Mishra SC, Mahanta P. Analysis of collimated radiation in participating media using the discrete transfer method. Journal of Quantitative Spectroscopy and Radiative Transfer 2005;96(1):123–35. [5] Pontaza JP, Reddy JN. Least-squares finite element formulations for one-dimensional radiative transfer. Journal of Quantitative Spectroscopy and Radiative Transfer 2005;95(3):387–406. [6] Siegel R, Spuckler CM. Effects of refractive index and diffuse or specular boundaries on a radiating isothermal layer. Journal of Heat Transfer 1994;116(3):787–90. [7] Machali HF. Radiative transfer in participating media under conditions of radiative equilibrium. Journal of Quantitative Spectroscopy and Radiative Transfer 1995;53(2):201–10. [8] Machali HF, Madkour MA. Radiative transfer in a participating slab with anisotropic scattering and general boundary conditions. Journal of Quantitative Spectroscopy and Radiative Transfer 1995;54(5):803–13. [9] Maruyama S. Radiative heat transfer in anisotropic scattering media with specular boundary subjected to collimated irradiation. International Journal of Heat and Mass Transfer 1998;41(18):2847–56. [10] Tan H-P, Yi H-L, Wang PY, et al. Ray tracing method for transient coupled heat transfer in an anisotropic scattering layer. International Journal of Heat and Mass Transfer 2004;47(19–20):4045–59. [11] Zhou H-C, Lou C, Cheng Q, et al. Experimental investigations on visualization of three-dimensional temperature distributions in a large-scale pulverized-coal-fired boiler furnace. Proceedings of the Combustion Institute 2005;30(1):1699–706. [12] Zhou H-C, Chen D-L, Cheng Q. A new way to calculate radiative intensity and solve radiative transfer equation through using the Monte Carlo method. Journal of Quantitative Spectroscopy and Radiative Transfer 2004;83(3–4):459–81. [13] Cheng Q, Zhou H-C. The DRESOR method for a collimated irradiation on an isotropically scattering layer. Journal of Heat Transfer, Transactions of the ASME 2007;129(5):634–45. [14] Zhou H-C, Cheng Q, Huang Z-F, et al. The influence of anisotropic scattering on the radiative intensity in a gray, plane-parallel medium calculated by the DRESOR method. Journal of Quantitative Spectroscopy and Radiative Transfer 2007;104(1):99–115. [15] Ozisik MN. Radiative transfer and interactions with conduction and convection. New York: Wiley-Interscience; 1973. [16] Pao CV. An iterative method for radiative transfer is a slab with specularly reflecting boundary. Applied Mathematics and Computation 1975;1(4):353–64. [17] Sutton WH, Ozisik MN. An iterative solution for anisotropic radiative transfer in a slab. Journal of Heat Transfer 1979;101:695–8. [18] Modest MF. Radiative heat transfer. 2nd ed. San Diego: Academic Press; 2003. [19] Yang W-J, Taniguchi H, Kudo K. Radiative heat transfer by the Monte Carlo method. In: Hartnett JP, Irvines FT, editors. Advances in heat transfer, vol. 27. San Diego: Academic Press; 1995. p. 1–215. [20] Howell JR. The Monte Carlo method in radiative heat transfer. Journal of Heat Transfer—Transactions of the ASME 1998;120(3):547–60. [21] Wang L, Yang J, Modest MF, Haworth DC. Application of the full-spectrum k-distribution method to photon Monte Carlo solvers. Journal of Quantitative Spectroscopy and Radiative Transfer 2006;104(2):297–304. [22] Tang KC, Brewster MQ. Analysis of molecular gas radiation: real gas property effects. Journal of Thermophysics and Heat Transfer 1999;13(4):460–6. [23] Wang A, Modest MF. Spectral Monte Carlo models for nongray radiation analyses in inhomogeneous participating media. Journal of Quantitative Spectroscopy and Radiative Transfer 2007;50(19–20):3877–89.