- Email: [email protected]
A new discrete ordinate algorithm for computing radiative transfer in one-dimensional atmospheres
Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
www.elsevier.com/locate/jqsrt
A new discrete ordinate algorithm for computing radiative transfer in one-dimensional atmospheres Hong-Shun Lia;∗ , Gilles Flamantb , Ji-Dong Lua a
School of Energy and Power Engineering, Huazhong University of Science & Technology, Wuhan 430074, China b IMP-CNRS, BP 5, Odeillo, 66125 Font-Romeu C/edex, France Received 30 August 2002; accepted 27 November 2002
Abstract Based on the discrete ordinates scheme with (an) in5nitely small weight(s) (DOS+ISW) which was recently proposed by the present authors, a new discrete ordinate algorithm is developed for computing the radiative transfer in 1-D atmospheres. The biggest feature of the new algorithm is that it is simple and universally applicable. It is also suited to parallel computation. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Radiative transfer; Atmospheric radiation; Discrete ordinates method
1. Introduction The ability to solve the radiative transfer equation in a fast and accurate fashion is central to the atmospheric science and remote sensing. Up to now, many methods for computing radiative transfer in the atmosphere have been proposed, such as the Eddington approximation [1], the FN method [2,3], the discrete ordinates method (DOM) [4–7], the Monte Carlo method [8], the spherical harmonics method [9,10], to name but a few. As pointed out in Ref. [11], this area of modeling is still the subject of ongoing research because it is so computationally demanding. The purpose of this work is to propose an easy-to-use comprehensive discrete ordinate algorithm for computing radiative transfer in 1-D atmospheres.
∗
Corresponding author. Fax: +86-27-8754-0724. E-mail address: [email protected] (H.-S. Li).
0022-4073/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-4073(02)00381-3
408
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Table 1 An ISW1 + LSO S8 quadrature set Direction
w
1 2 3 4 5 6 7 8 9 10 ISW1
0.1422555 0.1422555 0.1422555 0.1422555 0.5773503 0.5773503 0.5773503 0.8040087 0.8040087 0.9795543 0.4358899
0.1422555 0.5773503 0.8040087 0.9795543 0.1422555 0.5773503 0.8040087 0.1422555 0.5773503 0.1422555 0
0.9795543 0.8040087 0.5773503 0.1422555 0.8040087 0.5773503 0.1422555 0.5773503 0.1422555 0.1422555 0.9
0.1712359 0.0992284 0.0992284 0.1712359 0.0992284 0.4617179 0.0992284 0.0992284 0.0992284 0.1712359 1 × 10−100
2. DOS + ISW algorithm Recently, we proposed a new discrete ordinate scheme, named DOS + ISW (the discrete ordinates scheme with (an) in5nitely small weight(s)), to solve collimated irradiation problems [12]. In this and the following sections, we shall apply the DOS + ISW to atmospheric radiation problems. 2.1. DOS + ISW The DOS + ISW [12] is introduced brieHy in this subsection. The construction of the DOS + ISW is very simple: at 5rst, one or more new discrete directions are directly added to an existing discrete ordinate quadrature set (say, LSO S8 , LSH S10 [13], SRAP15 [14], etc.), and the weight(s) associated with the new discrete direction(s) is (are) set to be in5nitely small. Then, the new discrete direction(s) and the original quadrature set make up a new discrete ordinate quadrature set. As an example, Table 1 lists a new discrete direction with the weight of 1 × 10−100 added to the LSO S8 (because of the symmetry, only the ten directions in the 5rst octant of the LSO S8 are given). In Table 1, , and are the direction cosines with respect to x-, y- and z-axis, respectively, and w angular weight. This new discrete direction and the 80 discrete directions of the LSO S8 form a new discrete ordinate quadrature set, which is referred to as ISWM + LSO S8 , here M is the total number of the discrete direction(s) with (an) in5nitely small weight(s). In Table 1, M = 1. This kind of scheme is referred to as the discrete ordinate scheme with (an) in5nitely small weight(s) (DOS + ISW). In Ref. [12], the feasibility of the DOS + ISW is veri5ed theoretically and tested with a standard test problem in Ref. [15]. Since the weights of the new directions are in5nitely small, the new discrete directions have no inHuence on the moment equations of the DOS + ISW. Consequently, we can add any number of the new discrete directions with in5nitely small weights to any existing discrete ordinate quadrature set, without exception, and the new discrete directions can be arbitrarily arranged [12].
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
409
2.2. How to apply the DOS + ISW to atmospheric radiation problems When solving an atmospheric radiation problem, we need not only to capture the solar beam, but also to calculate the radiation intensities in many speci5ed directions (assuming that we need to calculate the radiation intensities in M speci5ed directions). This kind of problem can be easily solved by using the DOS+ISW. At 5rst, add one new discrete directions with in5nite small weights to an existing discrete ordinate quadrature set, say, LSO S8 , and let it be coincident with the direction of the solar beam; then, add M new discrete directions with in5nite small weights to the above quadrature set, they are used to calculate the intensities in M speci5ed directions. The M + 1 new directions and the LSO S8 make up a new quadrature set, referred to as ISWM +1 + LSO S8 . As aforementioned, the new discrete directions can be arbitrarily arranged; therefore, no matter which direction the solar beam comes from, it can be always exactly captured by the ISWM +1 + LSO S8 , meanwhile, ISWM +1 + LSO S8 can be used to calculate the radiation intensities in any directions. 2.3. Discrete ordinate equations In a 3-D Cartesian coordinates system, for the N discrete directions in the original quadrature set, the discrete ordinate equations may be written as [16] i
@Ii @Ii @Ii + i + i + Ii = Si ; @x @y @z
i = 1; 2; : : : ; N;
(1a)
where I is the radiation intensity, distinction coeJcient, and S source function. Eq. (1a) is valid for a gray medium, or on a spectral basis, for a nongray medium. For the new discrete direction with an in5nitely small weight which is coincident with the direction of the solar beam, the corresponding discrete ordinate equations are formally the same as Eq. (1a). But for ease of presentation, a subscript “s” is employed in order to distinguish it from Eq. (1a), namely s
@Is @Is @Is + s + s + Is = Ss : @x @y @z
(1b)
Similarly, for the other M new discrete directions with in5nitely small weights, a subscript “p” is added in order to distinguish them from Eqs. (1a) and (1b) pi
@Ipi @Ipi @Ipi + pi + pi + Ipi = Spi ; @x @y @z
i = 1; 2; : : : ; M:
(1c)
The source function in Eq. (1a) may be written as N
Si = (1 − !)Ib +
M
! ! ! ws si Is + wj ji Ij + wpj ji Ipj ; 4 j=1 4 4 j=1
i = 1; 2; : : : ; N;
(2)
where is the scattering phase function. Because wpj (j = 1; : : : ; M ) are in5nitely small, and Ipj (j = 1; : : : ; M ) are not in5nitely large, thus the fourth term on the right-hand side can be
410
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
ignored, and Eq. (2) is simpli5ed to N
Si = (1 − !)Ib +
! ! ws si Is ; wj ji Ij + 4 j=1 4
i = 1; 2; : : : ; N:
(3a)
Similarly, source functions in Eqs. (1b) and (1c) can be obtained as N
! ! ws ss Is ; Ss = (1 − !)Ib + wj js Ij + 4 j=1 4
(3b)
N
! ! ws si Is ; Spi = (1 − !)Ib + wj ji Ij + 4 j=1 4
i = 1; 2; : : : ; M:
(3c)
Eqs. (3a)–(3c) show that the source functions of Eqs. (1a)–(1c) have nothing to do with Ipj (j = 1; 2; : : : ; M ). 2.4. Boundary conditions It is clear that the DOS + ISW is the same as the standard DOM except for the new discrete directions with in5nitely small weights (here the standard DOM means the DOM described in Ref. [16]). Therefore, the treatment of the boundary conditions is the same as the standard DOM except for the boundary condition for the new discrete direction being coincident with the solar beam. See Fig. 1, assuming that ss is the direction from which the solar beam impinges onto the top of the atmosphere. The intensity of the solar beam incident on the top of the atmosphere may be written as [16, pp. 572]. Isun ( ; s0 ) = q0 [s − ss ];
= 0;
(4)
where q0 is the total radiative heat Hux of the solar beam through a surface normal to ss . Most of previous investigators employed Eq. (4) to calculate the intensity of the solar beam. However, in the DOS + ISW, Eq. (4) is replaced by Eq. (5b) in the following. Since the solar beam is “represented” by the new discrete direction being coincident with s0 , q0 can be expressed as q0 = Is ( )ws ;
=0
(5a)
or Is ( ) = q0 =ws ;
= 0:
(5b)
Eq. (5b) is used as the boundary condition for Eq. (1b). It is seen that if ws tends to zero, then Is ( ) will approach in5nity, and the control-angle of Is ( ) will approach in5nitely small (note that the control-angle of a discrete direction is numerically equal to its weight [16]). This trend just conforms to the physical characteristics of collimated irradiation, since the intensity of a collimated beam is zero for all directions except for its propagating direction where it is in5nitely large [16, pp. 573].
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
411
solar beam (ss)
τ =0
cos -1 µs τ
τ = τ0
ground
Fig. 1. Radiative transfer in the atmosphere.
As Modest pointed out [16, pp. 574]: “Collimated irradiation could simply be treated as “strongly directional emission”. However, the discontinuity of intensity with direction causes problems with analytical as well as numerical solution techniques, thus warranting a separate approach for this type of problem”. In the DOS + ISW, the problems caused by the discontinuity of intensity with direction are solved with Eq. (5b). On the top of the atmosphere, the boundary conditions for the other discrete directions may be written as Ii = 0 Ipi = 0
for i ¿ 0;
(6a)
for pi ¿ 0:
(6b)
Assuming that the directional emissivity and bidirectional reHectivity of the ground are #(si ) and $(sj ; si ), respectively. The corresponding boundary conditions at the bottom of the atmosphere may be written as [16] Ii = #(si )Ib + $(sj ; si )wj Ij + $(ss ; si )ws Is + $(spj ; si )wpj Ipj for i ¡ 0; = 0 ; (7) j ¿0
pj ¿0
where Ib is the black body intensity at the temperature of the ground, si stands for an arbitrary discrete direction which propagates upwards. Considering that wpj (s) are in5nitely small, and that Ipj are not in5nitely large, Eq. (7) can be simpli5ed to $(sj ; si )wj Ij + $(ss ; si )ws Is for i ¡ 0; = 0 : (8a) Ii = #(si )Ib + j ¿0
412
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Fig. 2. Solution of the discrete ordinate equations.
Similarly, for the new discrete direction with in5nitely small weights Ipi = #(spi )Ib + $(sj ; spi )wj Ij + $(ss ; spi )ws Is for pi ¡ 0; = 0 :
(8b)
j ¿0
Eqs. (8a) and (8b) show that, like the source functions, the boundary conditions at the bottom of the atmosphere also have nothing to do with Ipj (j = 1; 2; : : : ; M ). 2.5. Numerical procedures As aforementioned, there is no diNerence between the DOS + ISW and the standard DOM except for the new discrete directions with in5nitely small weights; therefore, any numerical method for solving the discrete ordinate equations in the standard DOM can be directly used for the DOS+ISW, such as Fiveland’s method [17], Kumar and his coworkers’ method [18], etc. Of course, if only one or several new discrete directions are added to an existing discrete ordinate quadrature set, there is no problem to SIMULTANEOUSLY solve the N +M +1 discrete ordinate equations. However, if M is very large, it may be impossible to simultaneously solve the discrete ordinate equations, because the computer storage requirement will become amazingly prohibitive. Thus we have to develop new numerical procedures so as to remove this obstacle. As mentioned above, for Eqs. (1a) and (1b), their source functions and boundary conditions have nothing to do with Ipj (j=1; 2; : : : ; M ); therefore, the numerical procedures for solving Eqs. (1a)–(1c) can be divided into two steps as follows. Step A: Only Eqs. (1a) and (1b) are solved. This step is the same as that of the standard DOM, thus any suitable numerical method can be employed. In this work, Fiveland’s method [17] is employed. Unfortunately, Fiveland’s method [17] is conditionally stable; therefore it is modi5ed in order to make it be unconditionally stable. See Fig. 2, the 1-D atmosphere is divided into a certain number of layers, in each layer, all the parameters (except for the intensity) are considered
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
413
as constants. For an arbitrary discrete direction upwards, evaluating Eq. (1a) between two layers gives the following form: IB = In exp(−R n =2) + Sn [1 − exp(−R n =2)];
(9a)
In+1 = IB exp(−R n+1 =2) + Sn+1 [1 − exp(−R n+1 =2)];
(9b)
where I refer to the intensity in an arbitrary discrete direction whose direction cosine is , Sn and Sn+1 are determined with Eqs. (3a) or (3b) by using the intensities in the nth and (n + 1)th layers, respectively. For the discrete directions downwards, use the same method as above. The other numerical procedures are the same as that in Ref. [17], the reader may consult Ref. [17] for the details, which will not repeated here for brevity. Obviously, an iterative procedure is necessary, it continues until the convergence is established. In this work, it is speci5ed when RIjn =Ijn ¡ 1 × 10−6 ;
(10)
where Ijn is the intensity in the jth discrete direction at the nth layer, and RIjn represents the change in the intensity from one iteration to the next iteration. Step B: Eq. (1c) is solved using the method described in the Step A. Correspondingly, the source functions are determined using Eq. (3c), and the boundary conditions using Eqs. (6b) and (8b). However, not like in the Step A, the M new discrete ordinate equations are not solved simultaneously, since M may be a very large number. On the contrary, Eqs. (1c), (3c), (6b) and (8b) show that these M discrete ordinate equations are independent of each other, thus they can be INDEPENDENTLY solved one by one by using the same method as described in the Step A, and parallel computation can be adopted. Note that Eqs. (3c), (6b) and (8b) show that the source functions and boundary conditions of Eq. (1c) are dependent on only Is and Ii (i = 1; 2; : : : ; N ). Since the iteration in the Step A has been convergent, it is not necessary to iterate again in the Step B. In fact, Eqs. (1c), (3c), (6b) and (8b) show that Ipi (i = 1; 2; : : : ; M ) in all the layers can be directly obtained. Once the intensity 5led is determined, the other parameters, such as the heat Hux, the albedo, and so on, can be readily obtained. The computations of these parameters are formally the SAME as that in the standard DOM, namely [16] q+ ( ) = j wj Ij ( ) + s ws Is ( ) + pj wpj Ipj ( ) j ¿0
=
pj ¿0
j wj Ij ( ) + s ws Is ( );
j ¿0
q− ( ) =
Ij j w j +
j ¡0
q( ) = q+ ( ) + q− ( );
pj ¡0
Ipj pj wpj =
(11a)
Ij j w j ;
(11b)
j ¡0
(12)
where q( ), q+ ( ) and q− ( ) are the net heat Hux, the heat Hux in the positive and negative -directions, respectively. Note that wpj (j = 1; 2; : : : ; M ) are in5nitely small. The albedo and transmissivity may be written as A = q− (0)=q+ (0);
(13)
414
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Tr = q+ ( 0 )=q+ (0);
(14)
respectively. For ease of presentation, hereafter the algorithm described in this section is referred to as the DOS + ISW algorithm. 3. Test problems In order to verify the DOS + ISW algorithm, we have made extensive tests for a variety of cases. For illustrative purpose, we employ three test problems to show the performance of the DOS + ISW algorithm. 3.1. Test problem 1: isotropical scattering Test problem 1 is taken from Ref. [6]. In Ref. [6], Stamnes and Conklin employed three test problems to verify their matrix DOM (i.e. the well-known DISORT later on [7]). Here we use the third test problem in Ref. [6] to verify the DOS + ISW algorithm, because among the three test problems, the third one is the most diJcult one. Consider a semi-in5nite, isotropically scattering atmosphere with an exponentially varying single scattering albedo, i.e. !( ) = exp(− =B);
(15)
where B is an arbitrary constant. The solar beam is incident on the top of the atmosphere in the direction of s =0:9. The objective of this test problem is to calculate the albedo and the exit angular intensity distribution. Garcia and Siewert [2] presented results obtained with their FN method for this test case. In this work, the FN solutions [2] will be employed as the benchmark. 3.2. Test problems 2 and 3: anisotropical scattering For test problem 2, consider an 1-D homogeneous atmosphere, the relevant parameters are as follows: Haze L phase function, single scattering albedo !=0:9, 0 =1, 0 =1, q0 =, the reHectivity of the ground is equal to zero. Test problem 3 is the same as test problem 2 except for 0 = 0:5. These two test problems were proposed by the Radiation Commission of the International Association of Meteorology and Atmospheric Physics [19] (hereafter referred to as “standard procedures”), and have been used as test cases by many scholars [5,10,20]. 4. Results 4.1. Results for test problem 1 Test problem 1 is solved using an ISW12 + LSO S8 quadrature set, the 12 new directions with in5nitely small weights are as follows: s = 0:9;
(16a)
p1 = −0:05;
(16b)
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
415
Table 2 Comparison of exit angular intensity distributions and albedo values for test problem 1 −
B=1
B = 10
This work
Ref. [6]a
Ref. [2]
This work
Ref. [6]a
Ref. [2]
0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.69127 0.65166 0.57393 0.50822 0.45447 0.41036 0.37374 0.34295 0.31675 0.29422 0.27464
0.69608 0.65568 0.57612 0.50942 0.45512 0.41070 0.37388 0.34298 0.31670 0.29410 0.27448
0.69801 0.65651 0.57637 0.50955 0.45523 0.41078 0.37396 0.34304 0.31676 0.29416 0.27454
0.97594 0.98762 0.98204 0.95860 0.92786 0.89443 0.86057 0.82740 0.79551 0.76515 0.73644
0.98160 0.99256 0.98468 0.95998 0.92860 0.89486 0.86078 0.82748 0.79548 0.76505 0.73628
0.98407 0.99353 0.98484 0.96002 0.92864 0.89488 0.86082 0.82752 0.79553 0.76511 0.73634
Albedo
0.20514
0.20509
0.20517
0.46559
0.46629
0.46626
a
Using an 8-stream approximation.
pi = (i − 1) × (−0:1);
i = 1; : : : ; 11:
(16c)
Among them, the 5rst one is used to capture the solar beam, while the rest are used to calculate the exit angular intensity distribution in 11 speci5ed directions. In order to make a comparison with the results in Ref. [6], we also adopt the total optical thickness of the atmosphere total = 30, so as to mimic the semi-in5nite situation, and also divide the atmosphere into 100 adaptive nodal points by also using Eqs. (15) and (16) in Ref. [6]. Table 2 lists the results for B = 1 and 10 obtained with the above ISW12 + LSO S8 , as well as the results in Ref. [6] (using the 8-stream approximation) and the FN solutions [2]. Note that, in the case of one-dimension and isotropical scattering, the LSO S8 is equivalent to an 8-stream approximation in Ref. [6]. It is seen that the present computations are very close to the FN solutions [2], and the biggest deviation is only 0.98%. Further computations show that the precision of the DOS + ISW algorithm can be improved by using higher order discrete ordinate quadrature sets. For example, if using an ISW12 + SRAP45 [14] quadrature set, which is equivalent to a 90-stream approximation, the present computation can agree 4 –5 decimal places (not shown here for brevity). Fig. 3 shows the exit angular intensity distribution for B = 1, 10, 100, 1000, 106 and ∞ obtained with the above ISW12 + LSO S8 , the results are in good agreement with the FN solutions [2] for all the cases except for B = 106 (in the case of B = ∞, a 5ctitious boundary with a diNuse reHectivity of unity is set at the position of = total = 30, since there is no absorption in the range of ¿ 30), the biggest deviation in these cases is only 1%. Unfortunately, in the case of B=106 the biggest deviation is up to 4.79%. This is because that the scattering albedo varies rapidly around = 106 (see Eq. (15)). Consequently, adopting total = 30 is not enough to represent the semi-in5nite situation. Thus, we adopt total = 30 × 106 , divide the atmosphere into 500 layers, and employ an ISW12 + SRAP8 [14]
416
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421 2
2.2 B=∞ 6
B = 10
1.5
1.7
I* (-)
B = 1000
1.2
1 B = 100
0.7
0.5
B=1 B = 10
0.2
0 0
0.2
0.4
0.6
0.8
1
− Fig. 3. The exit angular intensity distribution for test problem 1.
Table 3 The albedo for test problem 1 B
!0 = 0:7
1 10 102 103 106 ∞ a
0.9
0.999
1
Thisa
Ref. [2]
Thisa
Ref. [2]
Thisa
Ref. [2]
Thisa
Ref. [2]
0.11891 0.19672 0.21649 0.21894 0.21922 0.21922
0.11885 0.19699 0.21679 0.21924 0.21952 0.21952
0.17261 0.33968 0.41460 0.42813 0.42982 0.42982
0.17262 0.34028 0.41532 0.42885 0.43054 0.43054
0.20478 0.46372 0.69165 0.83063 0.92164 0.92210
0.20483 0.46466 0.69334 0.83281 0.91737 0.91773
0.20514 0.46559 0.69689 0.84428 0.94566 0.99873
0.20517 0.46626 0.69861 0.84656 0.98355 1.00000
This work, using the LSO S8 .
quadrature set, which is equivalent to a 16-stream approximation, to perform another computation, the corresponding result is shown in Fig. 3, the biggest deviation is reduced to less than 1%. Table 3 gives the albedo obtained with the above ISW12 + LSO S8 . The present computations are in good agreement with the FN solutions [2], the biggest deviation is only 0.47% for all the cases except for !0 = 1 and B = 106 . In that case the relative deviation is up to 4%. However, if using the same improvement as above (i.e. ISW12 + SRAP8 , total = 30 × 106 , 500 layers), the relative deviation is reduced to only 0.42%
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
417
4.2. Results for test problems 2 and 3 In order to solve test problem 2, we employ two new quadrature set: an ISW12 + LSH S12 and an ISW12 + SRAP12 (in each octant, the LSH S12 [21] and SRAP12 [14] have 21 and 90 discrete directions, respectively). The twelve new discrete directions in the former quadrature set are the same as that in the later one, namely i = 1:2 − 0:2 × i; 6 = 1 × 10−200 12 = 1:
i = 1; : : : ; 5; 7; : : : ; 11;
(17a)
(or 6 = −1 × 10−200 );
(17b) (17c)
Of them the 12th new direction is used to capture the solar beam, the other 11 ones, to calculate the intensities in 11 speci5ed directions (the sixth one is used for the direction of = 0. From Eqs. (9a) and (9b), it is seen that it is impossible to perform numerical computation if = 0). In addition, in the numerical computation, the atmosphere is divided into 200 layers. It is well known that theoretically a phase function satis5es 1 (*; * ) d* = 1; (18) 4 4 where * is the solid angle. In the DOM, the term on the left-hand side in Eq. (18) is discretized as 1 1 (*; * ) d* ≈ ij wj : (19) 4 4 4 j More often than not, the term on the right-hand side of Eq. (19) is not equal to one in the case of anisotropic scattering, this implies that the conservation of energy is not respected during the scattering process; therefore, we adopt the renormalizing factor fi proposed in Ref. [22], namely ij wj : (20) fi = 4 j
Then fi is applied to Eqs. (3a) and (3b). The reader is suggested to consult Ref. [22] for the details of fi , which is not repeated here for brevity. Table 4 gives intensity distribution for test problem 2 obtained with the above two quadrature sets. Since test problem 2 represents an overhead sun, the intensity is independent of azimuth. In this adhoc case, the ISW12 + LSH S12 and ISW12 + SRAP12 quadrature sets are equivalent to a 12-stream and a 24-stream approximation in the matrix DOM [5,6], respectively. It is seen that the present computations are in good agreement with the standard procedures solutions, the biggest deviation for the result using the ISW12 + LSH S12 is 2.36%, and for the result using the ISW12 + SRAP12 , only 0.22% except for the two cases of = 0:5 and = 1:0, and = 1:0 and = 1:0. In these two cases, the relative deviations with respect to the standard procedures solution are about 0.4%, however the relative deviations with respect to the Spherical Harmonics solution [10] are only 0.06%. Similar results for test problem 3 and for , − ,0 = 0 are tabulated in Table 5, the results are also obtained with the above two quadrature sets (but 12 = 0:5). It is seen that the situation in Table 5 is similar to that in Table 4. Table 6 shows the net Huxes at several optical depths for test problem 3,
418
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Table 4 Intensities for test problem 2
This work LSH S12 [21]
SRAP12 [14]
Ref. [5] DOM
Ref. [10] Spherical harmonics
Standard procedures [19] Spherical harmonicsa
Successive scatteringa
0
−1:0 −0:8 −0:6 −0:4 −0:2 −0:0
2:7776 − 2 3:1542 − 2 3:9608 − 2 5:3814 − 2 6:6692 − 2 5:2839 − 2
2:7988 − 2 3:1403 − 2 3:9075 − 2 5:3407 − 2 6:6794 − 2 5:2025 − 2
2:7845 − 2 3:1435 − 2 3:9183 − 2 5:3588 − 2 6:7021 − 2 5:1775 − 2
2:7972 − 2 3:1448 − 2 3:9131 − 2 5:3511 − 2 6:6956 − 2 5:1800 − 2
2:7940 − 2 3:1450 − 2 3:9130 − 2 5:3470 − 2 6:6910 − 2
2:831 − 2 3:166 − 2 3:935 − 2 5:369 − 2 6:698 − 2 5:168 − 2
0.5
−1:0 −0:8 −0:6 −0:4 −0:2 0:0 0:2 0:4 0:6 0:8 1:0
1:3638 − 2 1:6090 − 2 2:2244 − 2 3:6945 − 2 6:7596 − 2 9:6504 − 2 9:2499 − 2 8:9834 − 2 1:1620 − 1 2:4174 − 1 2:2530 + 0
1:3758 − 2 1:6075 − 2 2:2061 − 2 3:6812 − 2 6:6660 − 2 9:4075 − 2 9:1500 − 2 8:8198 − 2 1:1519 − 1 2:3833 − 1 2:2491 + 0
1:3693 − 2 1:6090 − 2 2:2124 − 2 3:6897 − 2 6:6729 − 2 9:3993 − 2 9:1552 − 2 8:8285 − 2 1:1531 − 1 2:3851 − 1 2:2398 + 0
1:3752 − 2 1:6096 − 2 2:2098 − 2 3:6854 − 2 6:6697 − 2 9:3998 − 2 9:1592 − 2 8:8332 − 2 1:1534 − 1 2:3849 − 1 2:2484 + 0
1:3740 − 2 1:6100 − 2 2:2100 − 2 3:6830 − 2 6:6680 − 2 9:4000 − 2 9:1550 − 2 8:8340 − 2 1:1532 − 1 2:3850 − 1 2:2396 + 0
1:390 − 2 1:620 − 2 2:218 − 2 3:687 − 2 6:662 − 2 9:387 − 2 9:153 − 2 8:848 − 2 1:157 − 1 2:395 − 1 2:239 + 0
1.0
0:0 0:2 0:4 0:6 0:8 1:0
8:0126 − 2 1:2273 − 1 1:5165 − 1 2:0292 − 1 3:9607 − 1 2:9912 + 0
7:9897 − 2 1:2396 − 1 1:4794 − 1 2:0033 − 1 3:8652 − 1 2:9796 + 0
7:9911 − 2 1:2417 − 1 1:4822 − 1 2:0067 − 1 3:8696 − 1 2:9672 + 0
7:9400 − 2 1:2421 − 1 1:4827 − 1 2:0070 − 1 3:8692 − 1 2:9777 + 0
1:2413 − 1 1:4826 − 1 2:0068 − 1 3:8694 − 1 2:9669 + 0
7:971 − 2 1:247 − 1 1:491 − 1 2:019 − 1 3:893 − 1 2:975 + 0
a
Taken from Ref. [5].
both results obtained with the ISW12 + LSH S12 and ISW12 + SRAP12 agree very well with the standard procedures solutions. On a personal computer with a CPU of Intel Pentium III 800, the computing times for the results obtained with the ISW12 + LSH S12 and ISW12 + SRAP12 in Table 4 are 2 and 26 s, respectively, and for the results in Table 5, 2 and 26 s again. 4.3. EAect of the weights of the new directions In order to examine the eNect of the weights of the new discrete directions on the computational results, diNerent weight values are tested, such as 10−5 , 10−6 , 10−7 , 10−8 and so on. Numerical tests show that, for the three test problems, ws and wpi (i=1; 2; : : : ; M ) have no eNect on the computational results if they are less than 1 × 10−10 . This implies that in a practical computation we can use any weight values provided they are small enough.
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
419
Table 5 Intensities for test problem 3
This work LSH S12 [21]
SRAP12 [14]
Ref. [5] DOM
Ref. [10] Spherical harmonics
Standard procedures [19] Spherical harmonicsa
Successive scatteringa
0
−1:0 −0:8 −0:6 −0:4 −0:2 −0:0
2:2535 − 2 6:4226 − 2 1:5170 − 1 3:2408 − 1 6:4837 − 1 1:0324 + 0
2:2845 − 2 6:5001 − 2 1:5104 − 1 3:2944 − 1 6:5636 − 1 1:0317 + 0
2:2856 − 2 6:4983 − 2 1:5096 − 1 3:2932 − 1 6:5686 − 1 1:0319 + 0
2:2819 − 2 6:4998 − 2 1:5099 − 1 3:2933 − 1 6:5683 − 1 1:0320 + 0
2:2810 − 2 6:5000 − 2 1:5100 − 1 3:2920 − 1 6:5640 − 1
2:2780 − 2 6:4620 − 2 1:5010 − 1 3:2740 − 1 6:5220 − 1 1:0300 + 0
0.5
−1:0 −0:8 −0:6 −0:4 −0:2 0:0 0:2 0:4 0:6 0:8 1:0
9:3060 − 3 2:6517 − 2 6:8688 − 2 1:7496 − 1 4:5823 − 1 1:0237 + 0 1:8258 + 0 2:9408 + 0 2:2271 + 0 5:7242 − 1 4:7335 − 2
9:3520 − 3 2:6608 − 2 6:8148 − 2 1:7539 − 1 4:5197 − 1 1:0087 + 0 1:8271 + 0 2:8729 + 0 2:1837 + 0 5:6945 − 1 4:7633 − 2
9:3645 − 3 2:6604 − 2 6:8109 − 2 1:7521 − 1 4:5150 − 1 1:0088 + 0 1:8265 + 0 2:8728 + 0 2:1840 + 0 5:6966 + 0 4:7590 − 2
9:3470 − 3 2:6608 − 2 6:8117 − 2 1:7521 − 1 4:5147 − 1 1:0087 + 0 1:8263 + 0 2:8730 + 0 2:1840 + 0 5:6958 − 1 4:7508 − 2
9:3440 − 3 2:6610 − 2 6:8110 − 2 1:7520 − 1 4:5140 − 1 1:0090 + 0 1:8260 + 0 2:8730 + 0 2:1840 + 0 5:6960 − 1 4:7600 − 2
9:3400 − 3 2:6500 − 2 6:7840 − 2 1:7450 − 1 4:5010 − 1 1:0080 + 0 1:8250 + 0 2:8620 + 0 2:1780 + 0 5:7050 − 1 4:7720 − 2
1.0
0:0 0:2 0:4 0:6 0:8 1:0
5:4019 − 1 1:2978 + 0 2:5316 + 0 2:2401 + 0 7:2811 − 1 8:3202 − 2
5:2720 − 1 1:2909 + 0 2:4006 + 0 2:1531 + 0 7:1959 − 1 8:3838 − 2
5:2772 − 1 1:2916 + 0 2:4008 + 0 2:1537 + 0 7:2002 − 1 8:3735 − 2
5:2400 − 1 1:2914 + 0 2:4010 + 0 2:1537 + 0 7:1993 − 1 8:3758 − 2
1:2910 + 0 2:4010 + 0 2:1530 + 0 7:1990 − 1 8:3750 − 2
5:2310 − 1 1:2890 + 0 2:3920 + 0 2:1480 + 0 7:2060 − 1 8:4030 − 2
a
Taken from Ref. [5].
5. Summary From Tables 2, 4 and 5, it is seen that the matrix DOM solutions [5,6] are slightly better than the present computations. Usually the DOS + ISW algorithm needs higher order discrete ordinate quadrature sets to achieve the same precision as that of the matrix DOM [5,6]. But the DOS + ISW algorithm has the following proper features: (1) Like the standard DOM, the DOS + ISW algorithm may be carried out to any arbitrary order and accuracy. (2) It is very simple, and is unconditionally stable. In the matrix DOM [5–7], the radiation intensity is expressed as the Fourier series, the user needs to list a matrix of coeJcients, and to compute the eigenvalues and eigenvectors. Whereas in the DOS + ISW algorithm, all these procedures are not necessary.
420
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Table 6 Net Huxes for test problem 3
This work
0.0 0.05 0.1 0.2 0.5 0.75 1.0 a
LSH S12 [21]
SRAP12 [14]
Ref. [5] DOM
1.3546 1.3352 1.3156 1.2765 1.1662 1.0878 1.0245
1.3450 1.3254 1.3055 1.2661 1.1563 1.0788 1.0155
1.3453 1.3257 1.3059 1.2664 1.1568 1.0791 1.0159
Ref. [10] spherical harmonics
Standard procedures [19] Spherical harmonicsa
Successive scatteringa
1.3453 1.3257 1.3059 1.2664 1.1568 1.0791 1.0159
1.3452 1.3257 1.3058 1.2664 1.1567 1.0791 1.0158
1.3453 1.3257 1.3059 1.2664 1.1568 1.0791 1.0159
Taken from Ref. [5].
(3) The new discrete directions have no inHuence on the moment equations of the DOS + ISW; therefore the user can add any number of the new discrete directions with in5nitely small weights to any discrete ordinate quadrature set in being, without exception, and the new discrete directions can be arbitrarily arranged. Consequently, no matter which direction the solar beam comes from, it can be exactly captured by the DOS + ISW algorithm, and the DOS + ISW algorithm can predict the radiation intensities in any directions. (4) The weights of the new discrete directions have no eNect on the computational results if they are small enough. In practical computations, we suggest that the user adopt the weights of less than 1 × 10−100 . (5) The DOS + ISW algorithm is universally applicable to any cases where the standard DOM is applicable to, such as homogeneous or inhomogeneous atmospheres, isotropic or anisotropic scattering, ! ¡ 1 or ! ≡ 1. Moreover, the atmosphere can be divided into any number of layers if the computer storage permits. The DOS + ISW algorithm could be used as an alternative method for calculating the radiative transfer in 1-D atmospheres. Acknowledgements This work is supported by the National Natural Science Foundation of China (No. 59906003), the Key Program of China for International Cooperation of Science and Technology (No. 2001CB711203), and the Science Foundation of Huazhong University of Science and Technology. References [1] Shettle EP, Weinman JA. The transfer of solar irradiance through inhomogeneous turbid atmosphere evaluated by Eddington’s approximation. J Atmos Sci 1970;27:1048–55. [2] Garcia RDM, Siewert CE. Radiative transfer in inhomogeneous atmospheres—numerical results. JQSRT 1981;25: 277–83.
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
421
[3] Garcia RDM, Siewert CE. The FN method in atmospheric radiative transfer. Int J Eng Sci 1998;36:1623–49. [4] Liou KN. A numerical experiment on Chandrasekhar’s discrete-ordinate method for radiative transfer: applications to cloudy and hazy atmospheres. J Atmos Sci 1973;30:1303–26. [5] Stamnes K, Dale H. A new look at the discretem ordinate method for radiative transfer calculations in anisotropically scattering atmospheres. II: intensity computation. J Atmos Sci 1981;38:2696–706. [6] Stamnes K, Conklin P. A new multi-layer discrete ordinate approach to radiative transfer in vertically inhomogeneous atmospheres. JQSRT 1984;31:273–82. [7] Stamnes K, Tsay S-C, Wiscombe W, Jawaweera K. Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Appl Opt 1988;27:2502–9. [8] Collins DG, Blattner WG, Wells MB, Horak HG. Backward Monte Carlo calculations of the polarization characteristics of the radiation emerging from spherical-shell atmospheres. Appl Opt 1972;11:2684–96. [9] Dave JV. A direct solution to the spherical harmonics approximation to the radiative transfer equation for an arbitrary solar elevation. Part I: theory. J Atmos Sci 1975;32:790–8. [10] Benassi M, Garcia RDM, Siewert CE. A high-order spherical harmonics solution to the standard problem in radiative transfer. Astrophys J 1984;280:853–64. [11] Evans KF. The spherical harmonics discrete ordinates method for three-dimensional atmospheric radiative transfer. J Atmos Sci 1998;55:429–46. [12] Li H-S, Flamant G, Lu J-D. An alternative discrete ordinate scheme for collimated irradiation problems. Int Comm Heat Mass Transfer, 2002, accepted for publication. [13] Fiveland WA. The selection of discrete ordinate quadrature sets for anisotropic scattering. In: Fiveland WA, et al., editors. Fundamental of radiation heat transfer, HTD-vol. 160. New York: ASME, 1991. p. 89 –96. [14] Li B-W, Yao Q, Cao X-Y, Cen K-F. A new discrete ordinate quadrature scheme for 3-D radiative heat transfer. ASME J Heat Transfer 1998;120:514–8. [15] Crosbie AL, Schrenker RG. Multiple scattering in a two-dimensional rectangular medium exposed to collimated radiation. J Quant Spectrosc Radiat Transfer 1985;33:101–25. [16] Modest MF. Radiative heat transfer. New York: McGraw-Hill, 1993. p. 541–71. [17] Fiveland WA. Discrete ordinate methods for radiative heat transfer in isotropically and anisotropically scattering media. J Heat Transfer 1987;109:809–12. [18] Kumar S, Majumdar A, Tien CL. The diNerential discrete-ordinate method for solutions of the equation of radiative transfer. J Heat Transfer 1990;112:424–9. [19] Lenoble J. Standard procedures to compute atmospheric radiative transfer in a scattering atmosphere. Boulder, Colorado: National Center for Atmosphere Research, 1977. [20] Wang Hong-Qi, Zhao Gao-Xiang. Computation of radiation intensity in strongly anisotropic scattering atmospheres by using the discrete ordinates method. Sci China Ser B (Chinese version) 1989; (12): 1330 –9. [21] Balsara D. Fast and accurate discrete ordinates methods for multidimensional radiative transfer. Part I, basic methods. JQSRT 2001;69:671–707. [22] Liu LH, Ruan LM, Tan HP. On the discrete ordinates method for radiative heat transfer in anisotropically scattering media. Int J Heat Mass Transfer 2002;45:3259–62.
www.elsevier.com/locate/jqsrt
A new discrete ordinate algorithm for computing radiative transfer in one-dimensional atmospheres Hong-Shun Lia;∗ , Gilles Flamantb , Ji-Dong Lua a
School of Energy and Power Engineering, Huazhong University of Science & Technology, Wuhan 430074, China b IMP-CNRS, BP 5, Odeillo, 66125 Font-Romeu C/edex, France Received 30 August 2002; accepted 27 November 2002
Abstract Based on the discrete ordinates scheme with (an) in5nitely small weight(s) (DOS+ISW) which was recently proposed by the present authors, a new discrete ordinate algorithm is developed for computing the radiative transfer in 1-D atmospheres. The biggest feature of the new algorithm is that it is simple and universally applicable. It is also suited to parallel computation. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Radiative transfer; Atmospheric radiation; Discrete ordinates method
1. Introduction The ability to solve the radiative transfer equation in a fast and accurate fashion is central to the atmospheric science and remote sensing. Up to now, many methods for computing radiative transfer in the atmosphere have been proposed, such as the Eddington approximation [1], the FN method [2,3], the discrete ordinates method (DOM) [4–7], the Monte Carlo method [8], the spherical harmonics method [9,10], to name but a few. As pointed out in Ref. [11], this area of modeling is still the subject of ongoing research because it is so computationally demanding. The purpose of this work is to propose an easy-to-use comprehensive discrete ordinate algorithm for computing radiative transfer in 1-D atmospheres.
∗
Corresponding author. Fax: +86-27-8754-0724. E-mail address: [email protected] (H.-S. Li).
0022-4073/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-4073(02)00381-3
408
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Table 1 An ISW1 + LSO S8 quadrature set Direction
w
1 2 3 4 5 6 7 8 9 10 ISW1
0.1422555 0.1422555 0.1422555 0.1422555 0.5773503 0.5773503 0.5773503 0.8040087 0.8040087 0.9795543 0.4358899
0.1422555 0.5773503 0.8040087 0.9795543 0.1422555 0.5773503 0.8040087 0.1422555 0.5773503 0.1422555 0
0.9795543 0.8040087 0.5773503 0.1422555 0.8040087 0.5773503 0.1422555 0.5773503 0.1422555 0.1422555 0.9
0.1712359 0.0992284 0.0992284 0.1712359 0.0992284 0.4617179 0.0992284 0.0992284 0.0992284 0.1712359 1 × 10−100
2. DOS + ISW algorithm Recently, we proposed a new discrete ordinate scheme, named DOS + ISW (the discrete ordinates scheme with (an) in5nitely small weight(s)), to solve collimated irradiation problems [12]. In this and the following sections, we shall apply the DOS + ISW to atmospheric radiation problems. 2.1. DOS + ISW The DOS + ISW [12] is introduced brieHy in this subsection. The construction of the DOS + ISW is very simple: at 5rst, one or more new discrete directions are directly added to an existing discrete ordinate quadrature set (say, LSO S8 , LSH S10 [13], SRAP15 [14], etc.), and the weight(s) associated with the new discrete direction(s) is (are) set to be in5nitely small. Then, the new discrete direction(s) and the original quadrature set make up a new discrete ordinate quadrature set. As an example, Table 1 lists a new discrete direction with the weight of 1 × 10−100 added to the LSO S8 (because of the symmetry, only the ten directions in the 5rst octant of the LSO S8 are given). In Table 1, , and are the direction cosines with respect to x-, y- and z-axis, respectively, and w angular weight. This new discrete direction and the 80 discrete directions of the LSO S8 form a new discrete ordinate quadrature set, which is referred to as ISWM + LSO S8 , here M is the total number of the discrete direction(s) with (an) in5nitely small weight(s). In Table 1, M = 1. This kind of scheme is referred to as the discrete ordinate scheme with (an) in5nitely small weight(s) (DOS + ISW). In Ref. [12], the feasibility of the DOS + ISW is veri5ed theoretically and tested with a standard test problem in Ref. [15]. Since the weights of the new directions are in5nitely small, the new discrete directions have no inHuence on the moment equations of the DOS + ISW. Consequently, we can add any number of the new discrete directions with in5nitely small weights to any existing discrete ordinate quadrature set, without exception, and the new discrete directions can be arbitrarily arranged [12].
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
409
2.2. How to apply the DOS + ISW to atmospheric radiation problems When solving an atmospheric radiation problem, we need not only to capture the solar beam, but also to calculate the radiation intensities in many speci5ed directions (assuming that we need to calculate the radiation intensities in M speci5ed directions). This kind of problem can be easily solved by using the DOS+ISW. At 5rst, add one new discrete directions with in5nite small weights to an existing discrete ordinate quadrature set, say, LSO S8 , and let it be coincident with the direction of the solar beam; then, add M new discrete directions with in5nite small weights to the above quadrature set, they are used to calculate the intensities in M speci5ed directions. The M + 1 new directions and the LSO S8 make up a new quadrature set, referred to as ISWM +1 + LSO S8 . As aforementioned, the new discrete directions can be arbitrarily arranged; therefore, no matter which direction the solar beam comes from, it can be always exactly captured by the ISWM +1 + LSO S8 , meanwhile, ISWM +1 + LSO S8 can be used to calculate the radiation intensities in any directions. 2.3. Discrete ordinate equations In a 3-D Cartesian coordinates system, for the N discrete directions in the original quadrature set, the discrete ordinate equations may be written as [16] i
@Ii @Ii @Ii + i + i + Ii = Si ; @x @y @z
i = 1; 2; : : : ; N;
(1a)
where I is the radiation intensity, distinction coeJcient, and S source function. Eq. (1a) is valid for a gray medium, or on a spectral basis, for a nongray medium. For the new discrete direction with an in5nitely small weight which is coincident with the direction of the solar beam, the corresponding discrete ordinate equations are formally the same as Eq. (1a). But for ease of presentation, a subscript “s” is employed in order to distinguish it from Eq. (1a), namely s
@Is @Is @Is + s + s + Is = Ss : @x @y @z
(1b)
Similarly, for the other M new discrete directions with in5nitely small weights, a subscript “p” is added in order to distinguish them from Eqs. (1a) and (1b) pi
@Ipi @Ipi @Ipi + pi + pi + Ipi = Spi ; @x @y @z
i = 1; 2; : : : ; M:
(1c)
The source function in Eq. (1a) may be written as N
Si = (1 − !)Ib +
M
! ! ! ws si Is + wj ji Ij + wpj ji Ipj ; 4 j=1 4 4 j=1
i = 1; 2; : : : ; N;
(2)
where is the scattering phase function. Because wpj (j = 1; : : : ; M ) are in5nitely small, and Ipj (j = 1; : : : ; M ) are not in5nitely large, thus the fourth term on the right-hand side can be
410
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
ignored, and Eq. (2) is simpli5ed to N
Si = (1 − !)Ib +
! ! ws si Is ; wj ji Ij + 4 j=1 4
i = 1; 2; : : : ; N:
(3a)
Similarly, source functions in Eqs. (1b) and (1c) can be obtained as N
! ! ws ss Is ; Ss = (1 − !)Ib + wj js Ij + 4 j=1 4
(3b)
N
! ! ws si Is ; Spi = (1 − !)Ib + wj ji Ij + 4 j=1 4
i = 1; 2; : : : ; M:
(3c)
Eqs. (3a)–(3c) show that the source functions of Eqs. (1a)–(1c) have nothing to do with Ipj (j = 1; 2; : : : ; M ). 2.4. Boundary conditions It is clear that the DOS + ISW is the same as the standard DOM except for the new discrete directions with in5nitely small weights (here the standard DOM means the DOM described in Ref. [16]). Therefore, the treatment of the boundary conditions is the same as the standard DOM except for the boundary condition for the new discrete direction being coincident with the solar beam. See Fig. 1, assuming that ss is the direction from which the solar beam impinges onto the top of the atmosphere. The intensity of the solar beam incident on the top of the atmosphere may be written as [16, pp. 572]. Isun ( ; s0 ) = q0 [s − ss ];
= 0;
(4)
where q0 is the total radiative heat Hux of the solar beam through a surface normal to ss . Most of previous investigators employed Eq. (4) to calculate the intensity of the solar beam. However, in the DOS + ISW, Eq. (4) is replaced by Eq. (5b) in the following. Since the solar beam is “represented” by the new discrete direction being coincident with s0 , q0 can be expressed as q0 = Is ( )ws ;
=0
(5a)
or Is ( ) = q0 =ws ;
= 0:
(5b)
Eq. (5b) is used as the boundary condition for Eq. (1b). It is seen that if ws tends to zero, then Is ( ) will approach in5nity, and the control-angle of Is ( ) will approach in5nitely small (note that the control-angle of a discrete direction is numerically equal to its weight [16]). This trend just conforms to the physical characteristics of collimated irradiation, since the intensity of a collimated beam is zero for all directions except for its propagating direction where it is in5nitely large [16, pp. 573].
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
411
solar beam (ss)
τ =0
cos -1 µs τ
τ = τ0
ground
Fig. 1. Radiative transfer in the atmosphere.
As Modest pointed out [16, pp. 574]: “Collimated irradiation could simply be treated as “strongly directional emission”. However, the discontinuity of intensity with direction causes problems with analytical as well as numerical solution techniques, thus warranting a separate approach for this type of problem”. In the DOS + ISW, the problems caused by the discontinuity of intensity with direction are solved with Eq. (5b). On the top of the atmosphere, the boundary conditions for the other discrete directions may be written as Ii = 0 Ipi = 0
for i ¿ 0;
(6a)
for pi ¿ 0:
(6b)
Assuming that the directional emissivity and bidirectional reHectivity of the ground are #(si ) and $(sj ; si ), respectively. The corresponding boundary conditions at the bottom of the atmosphere may be written as [16] Ii = #(si )Ib + $(sj ; si )wj Ij + $(ss ; si )ws Is + $(spj ; si )wpj Ipj for i ¡ 0; = 0 ; (7) j ¿0
pj ¿0
where Ib is the black body intensity at the temperature of the ground, si stands for an arbitrary discrete direction which propagates upwards. Considering that wpj (s) are in5nitely small, and that Ipj are not in5nitely large, Eq. (7) can be simpli5ed to $(sj ; si )wj Ij + $(ss ; si )ws Is for i ¡ 0; = 0 : (8a) Ii = #(si )Ib + j ¿0
412
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Fig. 2. Solution of the discrete ordinate equations.
Similarly, for the new discrete direction with in5nitely small weights Ipi = #(spi )Ib + $(sj ; spi )wj Ij + $(ss ; spi )ws Is for pi ¡ 0; = 0 :
(8b)
j ¿0
Eqs. (8a) and (8b) show that, like the source functions, the boundary conditions at the bottom of the atmosphere also have nothing to do with Ipj (j = 1; 2; : : : ; M ). 2.5. Numerical procedures As aforementioned, there is no diNerence between the DOS + ISW and the standard DOM except for the new discrete directions with in5nitely small weights; therefore, any numerical method for solving the discrete ordinate equations in the standard DOM can be directly used for the DOS+ISW, such as Fiveland’s method [17], Kumar and his coworkers’ method [18], etc. Of course, if only one or several new discrete directions are added to an existing discrete ordinate quadrature set, there is no problem to SIMULTANEOUSLY solve the N +M +1 discrete ordinate equations. However, if M is very large, it may be impossible to simultaneously solve the discrete ordinate equations, because the computer storage requirement will become amazingly prohibitive. Thus we have to develop new numerical procedures so as to remove this obstacle. As mentioned above, for Eqs. (1a) and (1b), their source functions and boundary conditions have nothing to do with Ipj (j=1; 2; : : : ; M ); therefore, the numerical procedures for solving Eqs. (1a)–(1c) can be divided into two steps as follows. Step A: Only Eqs. (1a) and (1b) are solved. This step is the same as that of the standard DOM, thus any suitable numerical method can be employed. In this work, Fiveland’s method [17] is employed. Unfortunately, Fiveland’s method [17] is conditionally stable; therefore it is modi5ed in order to make it be unconditionally stable. See Fig. 2, the 1-D atmosphere is divided into a certain number of layers, in each layer, all the parameters (except for the intensity) are considered
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
413
as constants. For an arbitrary discrete direction upwards, evaluating Eq. (1a) between two layers gives the following form: IB = In exp(−R n =2) + Sn [1 − exp(−R n =2)];
(9a)
In+1 = IB exp(−R n+1 =2) + Sn+1 [1 − exp(−R n+1 =2)];
(9b)
where I refer to the intensity in an arbitrary discrete direction whose direction cosine is , Sn and Sn+1 are determined with Eqs. (3a) or (3b) by using the intensities in the nth and (n + 1)th layers, respectively. For the discrete directions downwards, use the same method as above. The other numerical procedures are the same as that in Ref. [17], the reader may consult Ref. [17] for the details, which will not repeated here for brevity. Obviously, an iterative procedure is necessary, it continues until the convergence is established. In this work, it is speci5ed when RIjn =Ijn ¡ 1 × 10−6 ;
(10)
where Ijn is the intensity in the jth discrete direction at the nth layer, and RIjn represents the change in the intensity from one iteration to the next iteration. Step B: Eq. (1c) is solved using the method described in the Step A. Correspondingly, the source functions are determined using Eq. (3c), and the boundary conditions using Eqs. (6b) and (8b). However, not like in the Step A, the M new discrete ordinate equations are not solved simultaneously, since M may be a very large number. On the contrary, Eqs. (1c), (3c), (6b) and (8b) show that these M discrete ordinate equations are independent of each other, thus they can be INDEPENDENTLY solved one by one by using the same method as described in the Step A, and parallel computation can be adopted. Note that Eqs. (3c), (6b) and (8b) show that the source functions and boundary conditions of Eq. (1c) are dependent on only Is and Ii (i = 1; 2; : : : ; N ). Since the iteration in the Step A has been convergent, it is not necessary to iterate again in the Step B. In fact, Eqs. (1c), (3c), (6b) and (8b) show that Ipi (i = 1; 2; : : : ; M ) in all the layers can be directly obtained. Once the intensity 5led is determined, the other parameters, such as the heat Hux, the albedo, and so on, can be readily obtained. The computations of these parameters are formally the SAME as that in the standard DOM, namely [16] q+ ( ) = j wj Ij ( ) + s ws Is ( ) + pj wpj Ipj ( ) j ¿0
=
pj ¿0
j wj Ij ( ) + s ws Is ( );
j ¿0
q− ( ) =
Ij j w j +
j ¡0
q( ) = q+ ( ) + q− ( );
pj ¡0
Ipj pj wpj =
(11a)
Ij j w j ;
(11b)
j ¡0
(12)
where q( ), q+ ( ) and q− ( ) are the net heat Hux, the heat Hux in the positive and negative -directions, respectively. Note that wpj (j = 1; 2; : : : ; M ) are in5nitely small. The albedo and transmissivity may be written as A = q− (0)=q+ (0);
(13)
414
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Tr = q+ ( 0 )=q+ (0);
(14)
respectively. For ease of presentation, hereafter the algorithm described in this section is referred to as the DOS + ISW algorithm. 3. Test problems In order to verify the DOS + ISW algorithm, we have made extensive tests for a variety of cases. For illustrative purpose, we employ three test problems to show the performance of the DOS + ISW algorithm. 3.1. Test problem 1: isotropical scattering Test problem 1 is taken from Ref. [6]. In Ref. [6], Stamnes and Conklin employed three test problems to verify their matrix DOM (i.e. the well-known DISORT later on [7]). Here we use the third test problem in Ref. [6] to verify the DOS + ISW algorithm, because among the three test problems, the third one is the most diJcult one. Consider a semi-in5nite, isotropically scattering atmosphere with an exponentially varying single scattering albedo, i.e. !( ) = exp(− =B);
(15)
where B is an arbitrary constant. The solar beam is incident on the top of the atmosphere in the direction of s =0:9. The objective of this test problem is to calculate the albedo and the exit angular intensity distribution. Garcia and Siewert [2] presented results obtained with their FN method for this test case. In this work, the FN solutions [2] will be employed as the benchmark. 3.2. Test problems 2 and 3: anisotropical scattering For test problem 2, consider an 1-D homogeneous atmosphere, the relevant parameters are as follows: Haze L phase function, single scattering albedo !=0:9, 0 =1, 0 =1, q0 =, the reHectivity of the ground is equal to zero. Test problem 3 is the same as test problem 2 except for 0 = 0:5. These two test problems were proposed by the Radiation Commission of the International Association of Meteorology and Atmospheric Physics [19] (hereafter referred to as “standard procedures”), and have been used as test cases by many scholars [5,10,20]. 4. Results 4.1. Results for test problem 1 Test problem 1 is solved using an ISW12 + LSO S8 quadrature set, the 12 new directions with in5nitely small weights are as follows: s = 0:9;
(16a)
p1 = −0:05;
(16b)
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
415
Table 2 Comparison of exit angular intensity distributions and albedo values for test problem 1 −
B=1
B = 10
This work
Ref. [6]a
Ref. [2]
This work
Ref. [6]a
Ref. [2]
0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.69127 0.65166 0.57393 0.50822 0.45447 0.41036 0.37374 0.34295 0.31675 0.29422 0.27464
0.69608 0.65568 0.57612 0.50942 0.45512 0.41070 0.37388 0.34298 0.31670 0.29410 0.27448
0.69801 0.65651 0.57637 0.50955 0.45523 0.41078 0.37396 0.34304 0.31676 0.29416 0.27454
0.97594 0.98762 0.98204 0.95860 0.92786 0.89443 0.86057 0.82740 0.79551 0.76515 0.73644
0.98160 0.99256 0.98468 0.95998 0.92860 0.89486 0.86078 0.82748 0.79548 0.76505 0.73628
0.98407 0.99353 0.98484 0.96002 0.92864 0.89488 0.86082 0.82752 0.79553 0.76511 0.73634
Albedo
0.20514
0.20509
0.20517
0.46559
0.46629
0.46626
a
Using an 8-stream approximation.
pi = (i − 1) × (−0:1);
i = 1; : : : ; 11:
(16c)
Among them, the 5rst one is used to capture the solar beam, while the rest are used to calculate the exit angular intensity distribution in 11 speci5ed directions. In order to make a comparison with the results in Ref. [6], we also adopt the total optical thickness of the atmosphere total = 30, so as to mimic the semi-in5nite situation, and also divide the atmosphere into 100 adaptive nodal points by also using Eqs. (15) and (16) in Ref. [6]. Table 2 lists the results for B = 1 and 10 obtained with the above ISW12 + LSO S8 , as well as the results in Ref. [6] (using the 8-stream approximation) and the FN solutions [2]. Note that, in the case of one-dimension and isotropical scattering, the LSO S8 is equivalent to an 8-stream approximation in Ref. [6]. It is seen that the present computations are very close to the FN solutions [2], and the biggest deviation is only 0.98%. Further computations show that the precision of the DOS + ISW algorithm can be improved by using higher order discrete ordinate quadrature sets. For example, if using an ISW12 + SRAP45 [14] quadrature set, which is equivalent to a 90-stream approximation, the present computation can agree 4 –5 decimal places (not shown here for brevity). Fig. 3 shows the exit angular intensity distribution for B = 1, 10, 100, 1000, 106 and ∞ obtained with the above ISW12 + LSO S8 , the results are in good agreement with the FN solutions [2] for all the cases except for B = 106 (in the case of B = ∞, a 5ctitious boundary with a diNuse reHectivity of unity is set at the position of = total = 30, since there is no absorption in the range of ¿ 30), the biggest deviation in these cases is only 1%. Unfortunately, in the case of B=106 the biggest deviation is up to 4.79%. This is because that the scattering albedo varies rapidly around = 106 (see Eq. (15)). Consequently, adopting total = 30 is not enough to represent the semi-in5nite situation. Thus, we adopt total = 30 × 106 , divide the atmosphere into 500 layers, and employ an ISW12 + SRAP8 [14]
416
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421 2
2.2 B=∞ 6
B = 10
1.5
1.7
I* (-)
B = 1000
1.2
1 B = 100
0.7
0.5
B=1 B = 10
0.2
0 0
0.2
0.4
0.6
0.8
1
− Fig. 3. The exit angular intensity distribution for test problem 1.
Table 3 The albedo for test problem 1 B
!0 = 0:7
1 10 102 103 106 ∞ a
0.9
0.999
1
Thisa
Ref. [2]
Thisa
Ref. [2]
Thisa
Ref. [2]
Thisa
Ref. [2]
0.11891 0.19672 0.21649 0.21894 0.21922 0.21922
0.11885 0.19699 0.21679 0.21924 0.21952 0.21952
0.17261 0.33968 0.41460 0.42813 0.42982 0.42982
0.17262 0.34028 0.41532 0.42885 0.43054 0.43054
0.20478 0.46372 0.69165 0.83063 0.92164 0.92210
0.20483 0.46466 0.69334 0.83281 0.91737 0.91773
0.20514 0.46559 0.69689 0.84428 0.94566 0.99873
0.20517 0.46626 0.69861 0.84656 0.98355 1.00000
This work, using the LSO S8 .
quadrature set, which is equivalent to a 16-stream approximation, to perform another computation, the corresponding result is shown in Fig. 3, the biggest deviation is reduced to less than 1%. Table 3 gives the albedo obtained with the above ISW12 + LSO S8 . The present computations are in good agreement with the FN solutions [2], the biggest deviation is only 0.47% for all the cases except for !0 = 1 and B = 106 . In that case the relative deviation is up to 4%. However, if using the same improvement as above (i.e. ISW12 + SRAP8 , total = 30 × 106 , 500 layers), the relative deviation is reduced to only 0.42%
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
417
4.2. Results for test problems 2 and 3 In order to solve test problem 2, we employ two new quadrature set: an ISW12 + LSH S12 and an ISW12 + SRAP12 (in each octant, the LSH S12 [21] and SRAP12 [14] have 21 and 90 discrete directions, respectively). The twelve new discrete directions in the former quadrature set are the same as that in the later one, namely i = 1:2 − 0:2 × i; 6 = 1 × 10−200 12 = 1:
i = 1; : : : ; 5; 7; : : : ; 11;
(17a)
(or 6 = −1 × 10−200 );
(17b) (17c)
Of them the 12th new direction is used to capture the solar beam, the other 11 ones, to calculate the intensities in 11 speci5ed directions (the sixth one is used for the direction of = 0. From Eqs. (9a) and (9b), it is seen that it is impossible to perform numerical computation if = 0). In addition, in the numerical computation, the atmosphere is divided into 200 layers. It is well known that theoretically a phase function satis5es 1 (*; * ) d* = 1; (18) 4 4 where * is the solid angle. In the DOM, the term on the left-hand side in Eq. (18) is discretized as 1 1 (*; * ) d* ≈ ij wj : (19) 4 4 4 j More often than not, the term on the right-hand side of Eq. (19) is not equal to one in the case of anisotropic scattering, this implies that the conservation of energy is not respected during the scattering process; therefore, we adopt the renormalizing factor fi proposed in Ref. [22], namely ij wj : (20) fi = 4 j
Then fi is applied to Eqs. (3a) and (3b). The reader is suggested to consult Ref. [22] for the details of fi , which is not repeated here for brevity. Table 4 gives intensity distribution for test problem 2 obtained with the above two quadrature sets. Since test problem 2 represents an overhead sun, the intensity is independent of azimuth. In this adhoc case, the ISW12 + LSH S12 and ISW12 + SRAP12 quadrature sets are equivalent to a 12-stream and a 24-stream approximation in the matrix DOM [5,6], respectively. It is seen that the present computations are in good agreement with the standard procedures solutions, the biggest deviation for the result using the ISW12 + LSH S12 is 2.36%, and for the result using the ISW12 + SRAP12 , only 0.22% except for the two cases of = 0:5 and = 1:0, and = 1:0 and = 1:0. In these two cases, the relative deviations with respect to the standard procedures solution are about 0.4%, however the relative deviations with respect to the Spherical Harmonics solution [10] are only 0.06%. Similar results for test problem 3 and for , − ,0 = 0 are tabulated in Table 5, the results are also obtained with the above two quadrature sets (but 12 = 0:5). It is seen that the situation in Table 5 is similar to that in Table 4. Table 6 shows the net Huxes at several optical depths for test problem 3,
418
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Table 4 Intensities for test problem 2
This work LSH S12 [21]
SRAP12 [14]
Ref. [5] DOM
Ref. [10] Spherical harmonics
Standard procedures [19] Spherical harmonicsa
Successive scatteringa
0
−1:0 −0:8 −0:6 −0:4 −0:2 −0:0
2:7776 − 2 3:1542 − 2 3:9608 − 2 5:3814 − 2 6:6692 − 2 5:2839 − 2
2:7988 − 2 3:1403 − 2 3:9075 − 2 5:3407 − 2 6:6794 − 2 5:2025 − 2
2:7845 − 2 3:1435 − 2 3:9183 − 2 5:3588 − 2 6:7021 − 2 5:1775 − 2
2:7972 − 2 3:1448 − 2 3:9131 − 2 5:3511 − 2 6:6956 − 2 5:1800 − 2
2:7940 − 2 3:1450 − 2 3:9130 − 2 5:3470 − 2 6:6910 − 2
2:831 − 2 3:166 − 2 3:935 − 2 5:369 − 2 6:698 − 2 5:168 − 2
0.5
−1:0 −0:8 −0:6 −0:4 −0:2 0:0 0:2 0:4 0:6 0:8 1:0
1:3638 − 2 1:6090 − 2 2:2244 − 2 3:6945 − 2 6:7596 − 2 9:6504 − 2 9:2499 − 2 8:9834 − 2 1:1620 − 1 2:4174 − 1 2:2530 + 0
1:3758 − 2 1:6075 − 2 2:2061 − 2 3:6812 − 2 6:6660 − 2 9:4075 − 2 9:1500 − 2 8:8198 − 2 1:1519 − 1 2:3833 − 1 2:2491 + 0
1:3693 − 2 1:6090 − 2 2:2124 − 2 3:6897 − 2 6:6729 − 2 9:3993 − 2 9:1552 − 2 8:8285 − 2 1:1531 − 1 2:3851 − 1 2:2398 + 0
1:3752 − 2 1:6096 − 2 2:2098 − 2 3:6854 − 2 6:6697 − 2 9:3998 − 2 9:1592 − 2 8:8332 − 2 1:1534 − 1 2:3849 − 1 2:2484 + 0
1:3740 − 2 1:6100 − 2 2:2100 − 2 3:6830 − 2 6:6680 − 2 9:4000 − 2 9:1550 − 2 8:8340 − 2 1:1532 − 1 2:3850 − 1 2:2396 + 0
1:390 − 2 1:620 − 2 2:218 − 2 3:687 − 2 6:662 − 2 9:387 − 2 9:153 − 2 8:848 − 2 1:157 − 1 2:395 − 1 2:239 + 0
1.0
0:0 0:2 0:4 0:6 0:8 1:0
8:0126 − 2 1:2273 − 1 1:5165 − 1 2:0292 − 1 3:9607 − 1 2:9912 + 0
7:9897 − 2 1:2396 − 1 1:4794 − 1 2:0033 − 1 3:8652 − 1 2:9796 + 0
7:9911 − 2 1:2417 − 1 1:4822 − 1 2:0067 − 1 3:8696 − 1 2:9672 + 0
7:9400 − 2 1:2421 − 1 1:4827 − 1 2:0070 − 1 3:8692 − 1 2:9777 + 0
1:2413 − 1 1:4826 − 1 2:0068 − 1 3:8694 − 1 2:9669 + 0
7:971 − 2 1:247 − 1 1:491 − 1 2:019 − 1 3:893 − 1 2:975 + 0
a
Taken from Ref. [5].
both results obtained with the ISW12 + LSH S12 and ISW12 + SRAP12 agree very well with the standard procedures solutions. On a personal computer with a CPU of Intel Pentium III 800, the computing times for the results obtained with the ISW12 + LSH S12 and ISW12 + SRAP12 in Table 4 are 2 and 26 s, respectively, and for the results in Table 5, 2 and 26 s again. 4.3. EAect of the weights of the new directions In order to examine the eNect of the weights of the new discrete directions on the computational results, diNerent weight values are tested, such as 10−5 , 10−6 , 10−7 , 10−8 and so on. Numerical tests show that, for the three test problems, ws and wpi (i=1; 2; : : : ; M ) have no eNect on the computational results if they are less than 1 × 10−10 . This implies that in a practical computation we can use any weight values provided they are small enough.
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
419
Table 5 Intensities for test problem 3
This work LSH S12 [21]
SRAP12 [14]
Ref. [5] DOM
Ref. [10] Spherical harmonics
Standard procedures [19] Spherical harmonicsa
Successive scatteringa
0
−1:0 −0:8 −0:6 −0:4 −0:2 −0:0
2:2535 − 2 6:4226 − 2 1:5170 − 1 3:2408 − 1 6:4837 − 1 1:0324 + 0
2:2845 − 2 6:5001 − 2 1:5104 − 1 3:2944 − 1 6:5636 − 1 1:0317 + 0
2:2856 − 2 6:4983 − 2 1:5096 − 1 3:2932 − 1 6:5686 − 1 1:0319 + 0
2:2819 − 2 6:4998 − 2 1:5099 − 1 3:2933 − 1 6:5683 − 1 1:0320 + 0
2:2810 − 2 6:5000 − 2 1:5100 − 1 3:2920 − 1 6:5640 − 1
2:2780 − 2 6:4620 − 2 1:5010 − 1 3:2740 − 1 6:5220 − 1 1:0300 + 0
0.5
−1:0 −0:8 −0:6 −0:4 −0:2 0:0 0:2 0:4 0:6 0:8 1:0
9:3060 − 3 2:6517 − 2 6:8688 − 2 1:7496 − 1 4:5823 − 1 1:0237 + 0 1:8258 + 0 2:9408 + 0 2:2271 + 0 5:7242 − 1 4:7335 − 2
9:3520 − 3 2:6608 − 2 6:8148 − 2 1:7539 − 1 4:5197 − 1 1:0087 + 0 1:8271 + 0 2:8729 + 0 2:1837 + 0 5:6945 − 1 4:7633 − 2
9:3645 − 3 2:6604 − 2 6:8109 − 2 1:7521 − 1 4:5150 − 1 1:0088 + 0 1:8265 + 0 2:8728 + 0 2:1840 + 0 5:6966 + 0 4:7590 − 2
9:3470 − 3 2:6608 − 2 6:8117 − 2 1:7521 − 1 4:5147 − 1 1:0087 + 0 1:8263 + 0 2:8730 + 0 2:1840 + 0 5:6958 − 1 4:7508 − 2
9:3440 − 3 2:6610 − 2 6:8110 − 2 1:7520 − 1 4:5140 − 1 1:0090 + 0 1:8260 + 0 2:8730 + 0 2:1840 + 0 5:6960 − 1 4:7600 − 2
9:3400 − 3 2:6500 − 2 6:7840 − 2 1:7450 − 1 4:5010 − 1 1:0080 + 0 1:8250 + 0 2:8620 + 0 2:1780 + 0 5:7050 − 1 4:7720 − 2
1.0
0:0 0:2 0:4 0:6 0:8 1:0
5:4019 − 1 1:2978 + 0 2:5316 + 0 2:2401 + 0 7:2811 − 1 8:3202 − 2
5:2720 − 1 1:2909 + 0 2:4006 + 0 2:1531 + 0 7:1959 − 1 8:3838 − 2
5:2772 − 1 1:2916 + 0 2:4008 + 0 2:1537 + 0 7:2002 − 1 8:3735 − 2
5:2400 − 1 1:2914 + 0 2:4010 + 0 2:1537 + 0 7:1993 − 1 8:3758 − 2
1:2910 + 0 2:4010 + 0 2:1530 + 0 7:1990 − 1 8:3750 − 2
5:2310 − 1 1:2890 + 0 2:3920 + 0 2:1480 + 0 7:2060 − 1 8:4030 − 2
a
Taken from Ref. [5].
5. Summary From Tables 2, 4 and 5, it is seen that the matrix DOM solutions [5,6] are slightly better than the present computations. Usually the DOS + ISW algorithm needs higher order discrete ordinate quadrature sets to achieve the same precision as that of the matrix DOM [5,6]. But the DOS + ISW algorithm has the following proper features: (1) Like the standard DOM, the DOS + ISW algorithm may be carried out to any arbitrary order and accuracy. (2) It is very simple, and is unconditionally stable. In the matrix DOM [5–7], the radiation intensity is expressed as the Fourier series, the user needs to list a matrix of coeJcients, and to compute the eigenvalues and eigenvectors. Whereas in the DOS + ISW algorithm, all these procedures are not necessary.
420
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
Table 6 Net Huxes for test problem 3
This work
0.0 0.05 0.1 0.2 0.5 0.75 1.0 a
LSH S12 [21]
SRAP12 [14]
Ref. [5] DOM
1.3546 1.3352 1.3156 1.2765 1.1662 1.0878 1.0245
1.3450 1.3254 1.3055 1.2661 1.1563 1.0788 1.0155
1.3453 1.3257 1.3059 1.2664 1.1568 1.0791 1.0159
Ref. [10] spherical harmonics
Standard procedures [19] Spherical harmonicsa
Successive scatteringa
1.3453 1.3257 1.3059 1.2664 1.1568 1.0791 1.0159
1.3452 1.3257 1.3058 1.2664 1.1567 1.0791 1.0158
1.3453 1.3257 1.3059 1.2664 1.1568 1.0791 1.0159
Taken from Ref. [5].
(3) The new discrete directions have no inHuence on the moment equations of the DOS + ISW; therefore the user can add any number of the new discrete directions with in5nitely small weights to any discrete ordinate quadrature set in being, without exception, and the new discrete directions can be arbitrarily arranged. Consequently, no matter which direction the solar beam comes from, it can be exactly captured by the DOS + ISW algorithm, and the DOS + ISW algorithm can predict the radiation intensities in any directions. (4) The weights of the new discrete directions have no eNect on the computational results if they are small enough. In practical computations, we suggest that the user adopt the weights of less than 1 × 10−100 . (5) The DOS + ISW algorithm is universally applicable to any cases where the standard DOM is applicable to, such as homogeneous or inhomogeneous atmospheres, isotropic or anisotropic scattering, ! ¡ 1 or ! ≡ 1. Moreover, the atmosphere can be divided into any number of layers if the computer storage permits. The DOS + ISW algorithm could be used as an alternative method for calculating the radiative transfer in 1-D atmospheres. Acknowledgements This work is supported by the National Natural Science Foundation of China (No. 59906003), the Key Program of China for International Cooperation of Science and Technology (No. 2001CB711203), and the Science Foundation of Huazhong University of Science and Technology. References [1] Shettle EP, Weinman JA. The transfer of solar irradiance through inhomogeneous turbid atmosphere evaluated by Eddington’s approximation. J Atmos Sci 1970;27:1048–55. [2] Garcia RDM, Siewert CE. Radiative transfer in inhomogeneous atmospheres—numerical results. JQSRT 1981;25: 277–83.
H.-S. Li et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 407 – 421
421
[3] Garcia RDM, Siewert CE. The FN method in atmospheric radiative transfer. Int J Eng Sci 1998;36:1623–49. [4] Liou KN. A numerical experiment on Chandrasekhar’s discrete-ordinate method for radiative transfer: applications to cloudy and hazy atmospheres. J Atmos Sci 1973;30:1303–26. [5] Stamnes K, Dale H. A new look at the discretem ordinate method for radiative transfer calculations in anisotropically scattering atmospheres. II: intensity computation. J Atmos Sci 1981;38:2696–706. [6] Stamnes K, Conklin P. A new multi-layer discrete ordinate approach to radiative transfer in vertically inhomogeneous atmospheres. JQSRT 1984;31:273–82. [7] Stamnes K, Tsay S-C, Wiscombe W, Jawaweera K. Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Appl Opt 1988;27:2502–9. [8] Collins DG, Blattner WG, Wells MB, Horak HG. Backward Monte Carlo calculations of the polarization characteristics of the radiation emerging from spherical-shell atmospheres. Appl Opt 1972;11:2684–96. [9] Dave JV. A direct solution to the spherical harmonics approximation to the radiative transfer equation for an arbitrary solar elevation. Part I: theory. J Atmos Sci 1975;32:790–8. [10] Benassi M, Garcia RDM, Siewert CE. A high-order spherical harmonics solution to the standard problem in radiative transfer. Astrophys J 1984;280:853–64. [11] Evans KF. The spherical harmonics discrete ordinates method for three-dimensional atmospheric radiative transfer. J Atmos Sci 1998;55:429–46. [12] Li H-S, Flamant G, Lu J-D. An alternative discrete ordinate scheme for collimated irradiation problems. Int Comm Heat Mass Transfer, 2002, accepted for publication. [13] Fiveland WA. The selection of discrete ordinate quadrature sets for anisotropic scattering. In: Fiveland WA, et al., editors. Fundamental of radiation heat transfer, HTD-vol. 160. New York: ASME, 1991. p. 89 –96. [14] Li B-W, Yao Q, Cao X-Y, Cen K-F. A new discrete ordinate quadrature scheme for 3-D radiative heat transfer. ASME J Heat Transfer 1998;120:514–8. [15] Crosbie AL, Schrenker RG. Multiple scattering in a two-dimensional rectangular medium exposed to collimated radiation. J Quant Spectrosc Radiat Transfer 1985;33:101–25. [16] Modest MF. Radiative heat transfer. New York: McGraw-Hill, 1993. p. 541–71. [17] Fiveland WA. Discrete ordinate methods for radiative heat transfer in isotropically and anisotropically scattering media. J Heat Transfer 1987;109:809–12. [18] Kumar S, Majumdar A, Tien CL. The diNerential discrete-ordinate method for solutions of the equation of radiative transfer. J Heat Transfer 1990;112:424–9. [19] Lenoble J. Standard procedures to compute atmospheric radiative transfer in a scattering atmosphere. Boulder, Colorado: National Center for Atmosphere Research, 1977. [20] Wang Hong-Qi, Zhao Gao-Xiang. Computation of radiation intensity in strongly anisotropic scattering atmospheres by using the discrete ordinates method. Sci China Ser B (Chinese version) 1989; (12): 1330 –9. [21] Balsara D. Fast and accurate discrete ordinates methods for multidimensional radiative transfer. Part I, basic methods. JQSRT 2001;69:671–707. [22] Liu LH, Ruan LM, Tan HP. On the discrete ordinates method for radiative heat transfer in anisotropically scattering media. Int J Heat Mass Transfer 2002;45:3259–62.