- Email: [email protected]
Numerical solution of the unsteady radiative transfer equation
J. Quant. Spectrosc. Radiat. Transfer. Vol. 11, pp. 949-962. Pergamon Press 1971. Printed in Great Britain
NUMERICAL SOLUTION OF THE UNSTEADY RADIATIVE TRANSFER EQUATION B. C. HANKIN, J. G. HILL and A. G. P. WARHAM Mathematical Physics Department, A.W.R.E., Aldermaston, Berks. Abstract--In problems involving the transport of radiation through matter, it is necessary to solve the radiative transfer equation numerically. The scope for analytic solutions is restricted because of the complexity of the radiation-material interaction problem; this is particularly severe if the material is compressible. A study is made of various difference schemes within the framework of the discrete SN method tl~ with a view to obtaining a well-behaved solution. In addition to the usual requirements of consistency, stability and accuracy, there are further desirable features. Firstly, it is preferable to keep in the difference scheme, the photon conservation property as expressed by the basic transport equation. Secondly, the solution should be non-negative, not only because such values are unphysical, but also because of their effect on the radiation-material coupling. Such a coupling is usually non-linear, and negative solutions can lead to instabilities in the numerical solution. Finally, it is preferable that the method should lead to simple calculations. Results are given for an example involving the radiative heating of a semi-infinite slab. It is shown that the diffusion limit imposes a particularly severe test, and the method used by Lathrop in the computer code DTF-IV ~6) is favoured. A second example demonstrates a disturbing oscillatory behaviour arising from Diamond differencing which is basic to the DTF-IV method. A modification to this method designed to overcome this problem is examined. 1. I N T R O D U C T I O N
THE TRANSFER of radiation through matter is described by a linear Boltzmann equation when the photon-photon interaction is neglected. In practical problems the source function and mean free path often depend on material temperature, so that this equation is coupled with an equation describing the radiation-material interaction. This coupling can be very non-linear, reducing considerably the scope for analytic solutions, so that numerical methods have to be used. An added complexity arises if the material moves hydrodynamically. The present paper examines variations of the finite difference method known as the discrete SN method31) This method is chosen because it fulfills the usual requirements for a difference scheme of consistency and stability, and embodies the photon conservation property which is fundamental to the original transport equation. Much work has been done on the discrete SN method in the context of neutron transport. Such work has been aimed broadly at getting an accurate non-negative scheme. Our experience on radiation problems has shown that it is highly desirable that photon specific intensities are computed non-negative. Clearly negative intensities are unphysical but it is their consequences that can be more serious: the non-linear coupling between the transport equation and the photon-material interaction equation can lead to instabilities. F u r t h e r G R A N T (2) has pointed out that when, for example, mean free paths have a strong temperature dependence, it is necessary to ensure that the optically thick or diffusion limit is adequately treated. 949
950
B . C . HANKIN, J. G . HILL a n d A . G . P. WARHAM
There are also difficulties in calculating unsteady transport or multi-space dimensional transport in optically thin situations. These difficulties are common to both neutron and photon studies. 2. THE RADIATION-MATERIAL INTERACTION PROBLEM Let us consider the equation describing the plane unsteady transport of grey radiation in a purely absorbing medium in which local thermodynamic equilibrium exists. The transfer equation is c - xc3I/~t + p~3I/~3z + 2 - 11 = 2 - XB (1) where I = I(z, lz, t), the specific intensity at position z in direction cos- ~/~ at time t, c 2 B p T
= = = = =
velocity of light, 2(p, T) grey absorption mean free path, (ac/4rOT 4, the black body emission function, material density, material temperature
and ac/4 = Stefan's constant. If we integrate equation (1) over/~ we get the radiation energy equation ~3E/Ot + ~F/~z = 4rc2-1 [ B - (c/47t)E]
(2)
where +1
E = (2g/c) f
I d/~, the radiation energy density
-1
and
+
F = 2~ /
/~I d/~, the radiative flux.
1/ -1
The right hand side of equation (2) gives the difference between the emissions and absorptions. By equating this to the material cooling rate we get the following equation for the radiation-material interaction, - t3Em/8t = 4rc2- I ~ B - (c/4rc)E]
(3)
where E m = E,,(p, T), material energy density. No hydrodynamic effects are included in equation (3). The interaction problem consists of a solution of equations (1) and (3) with appropriate boundary conditions. E m is related to p and T by means of an equation of state. 2 is also a known function of p and T. The dependence of E m and 2 on p and T may exist in tabular form. 3. THE DISCRETE SN METHOD The discrete SN form for equation (1) is written (11 -- I°)/cA t + # j - 112([i -- I i - 1)/AZ + ~i---t1/2] =
~.Z-11/2Bi- 1/2
(4)
N u m e r i c a l s o l u t i o n o f the u n s t e a d y r a d i a t i v e transfer e q u a t i o n
951
where At, Az are the time and space mesh sizes respectively and 11, i0, I+ and Ii_ 1 are component intensities representing the flow of radiation energy across the space-time mesh as illustrated in Fig. 1. i is the mean intensity in the cell in the discrete direction #j_ 1/2 and determines the total absorption in the cell. 2+_ 1/2, B~_ 1/2 are average values of 2 and B in the time interval.
t I i I At
Ii-i
z
I°
FIG. 1. D e f i n i t i o n o f cell intensities.
If we multiply equation (4) by a weighting factor ogg_1/2 and sum over all j we get the following difference form of equation (2), the radiation energy equation : (E 1 -
E°)/At + (F i - Fi_ 1)/Az = 41z}tz_lu2[Bi_
1/2 -
(c/47r)E]
I5)
where E ° = (4n/c) ~ coj_ 1/2/°, similarly E 1 and E, i F i - 1 = 4~ E (Dj_ 1/21gj_ 1 / 2 1 / _ 1 and similarly Fi.
(6)
J
Equations (6) provide quadrature formulae for E and F depending on the choice of 09j_ 1/2, P j - 1/2. The more difficult problem of curved geometry is readily included formally in the discrete SN method. Considerable effort has been devoted to the choice of 09j_ 1/2, P j - 1/2 in the neutron context31) Integration over any space region for any time interval is readily applied to equation (5) thus ensuring the conservation property. Equation (4) introduces five unknown intensities so we require four more equations to complete the solution. Two equations are provided by the boundary conditions in the z and t directions so that two more equations are needed. These are the interpolation or
952
B.C. HANKIN, J. G. HILLand A. G. P. WARHAra
auxiliary relations which may be written in the following simple form :
i = ali+(1--a)li_l~ bI 1 + ( 1 - b ) I ° J
(7)
w h e r e 0 < a , b < 1. If we take a and b equal to ½ we complete the Diamond difference scheme. The method of solving equations (4) and (7) is described by CARLSON.(1)
4. THE NEGATIVITY PROBLEM For /~_ 1/2 positive we may regard 1i_ 1, I° as known (initially from the boundary conditions). If we eliminate i and 11 from equations (4) and (7) we get for the remaining unknown 11, (1 + aAz/bltj_ 1/2cAt + aAz/pj_
1/2~.i_ 1/2)Ii =
[1 - (1 - a)Az/bl~j_ 1/2cAt
- (1 - a)Az/#j_ 1/221-1/2]Ii- 1 + (Az/bpj_ 1/2cAt)I ° + (Az/t~j_ 1/22i- 1/2)Bi- 1/2.
(8)
Similarly for the unknown I1 we get (1 + bpj_ 1/2cAt/aAz + bcAt/2 i_ 1/2)11 = [1 - (1 - b)llj_ 1/2cAt/aAz - (1 - b)cAt/2 i_ 1/2110
+ (l~j_ 1/2cAt/aAz)I i_ 1 + (cAt/2i- 1/2)Bi- 1/2-
(9)
Clearly 11 and 11 are guaranteed non-negative if all the coefficients in equations (8) and (9) are positive, i is always guaranteed non-negative. Hence it is sufficient for (1 - a) < (Az/bltj_ 1~2cAt q- Az/#j_
i/2~.i_ 1/2)- 1
(1 -- b) < (Itj_ 1/2cAt/aAz + cat~2 I_ 1/2)- 1.
(10) (1 i)
The only values of a and b which satisfy inequalities (10) and (11) for all values are a = b = 1. Referring to equation (7) this corresponds to a forward difference or Step scheme. It is important to emphasize that inequalities (10) and (11) are sufficient but they are not in general necessary to give non-negativity. If we consider Diamond differencing in optically thin situations (AT,/J.i_l/2 small) relations (10) and (11) lead to
Hi_ 1/2cAt/Az _> 1
(12)
Az/laj_l/2CAt > 1.
(13)
and
Relations (12) and (13) are mutually exclusive except for the special case of the equality. Negative intensities will have a tendency therefore to arise if the factor !ui_ ~/2cAt/Az departs significantly from unity. In radiation-material interaction problems the timestep is often limited by the convergence in the solution of the interaction equation [equation (3)]. Hence inequality (12) tends to be violated giving rise to oscillations on the spatial mesh.
Numerical solution of the unsteady radiative transfer equation
953
It is clear also that negative values will tend to appear with Diamond differencing in optically thick problems (Az/2i_ 1/2 large). This problem has to be faced in radiative transfer studies. 5. SOLUTIONS TO THE NEGATIVITY PROBLEM The difficulty described in the previous section was recognized early in the development of the SN method and techniques for overcoming this were produced31~ In this Section these and more recent methods are described and a preliminary assessment is given in Section 6.
5.1 Forward difference or Step method This method was introduced in Section 4 above. The interpolation relations of equation (7) become ] = Ii = ll- 1 =I1
tgj_l{2 > O]
1~2- 1/2 < 0 i'"
(14)
For stability reasons the Solution of the SN equations proceeds for z increasing if/as _ 1/2 is positive and z decreasing for/22_ 1/2 negative. As shown in Section 4 above, the Step method guarantees non-negative intensities whatever the physical parameter values or mesh size. However, its truncation error is generally larger than that for D i a m o n d differencing. For this reason one of the earliest solutions ~1) to the negativity problem was a scheme combining the D i a m o n d and Step methods. Here a test is made to see whether the intensity (or angular flux in the neutron problem) is negative. If it is, the intensity is recalculated using a Step relation. 5.2 Woods method Many solutions have been attempted using a characteristics approach. The first simple application within the discrete S N framework is due to WOODSt3~ and has been described by LATHROP34J We give below a derivation of our own which follows WENDROFF.ts~ Wendroff's scheme does not inhibit negative intensities. We first rewrite equation (1) in characteristic form,
d l / d s + 2 - lI = 2 - 1 B
(15)
where dt/ds = c-1, dz/ds = p and s is the characteristic coordinate. If we assume 2, B constant within a space-time mesh we m a y write for the solution of equation (15) on the mesh I = I o exp( - s/2) + B[1 - e x p ( - s/2)]
(16)
where Io = intensity at s = 0. Let us now consider equation (4), the discrete S N form of the transport equation. We may write this as follows
( I k - - l k - 1)/As + ).7-11/21 = 27-11/2 Bi- 1/2
(17)
954
B.C. HANKIN, J. G. HILL and A. G. P. WARHAM
where
Ik- 1 = (I°/cat + #j- 1/2Ii - 1/Az)/As ] Ik = (I1/cAt + l~j- 1/2Ii/Az)/As
t"
(18)
and
As = cAtAz/(#j_ 1~2cAt + Az) We regard equation (17) as a difference form of equation (15) which implies that lk_ 1,
Ik lie on the same characteristic direction. Accordingly we use equation (16) to relate Ik_ 1 and lk, Ik = Ik- 1 exp(-- As/2i_ 1/2) + Bi-
1/2[ 1 - -
e x p ( - As~2 i_ 1/2)].
(19)
We are treating here the positive p direction. Equation (19) is decomposed by Woods into a pair of equations. The decomposition however depends on the mesh parameters. For cAt < Az/pj_ 1/2 we have
(Az/pj_ 1/2)11 = [(Az/pj_ 1/2 - cAt) I° + cAtli- 1] e x p ( - As~21_ 1/2) q- (Az//2j_ 1/2)Bi_ 1/211
-
e x p ( - As/,~ i_ 1/2)]
(20)
cAtI i = cAtI ° e x p ( - As/21_ 1/2) + cAtBi - x/2[ 1 - exp( - As/2i_ 1/2)]
(21)
and for cat > Az/p i_ 1/2 we have
(Az/ktj_ 1/2)11 = (Az/pj_ 1/2)Ii_ ~ e x p ( - As/2i_ 1/2) + (Az/p;_ 1/2)Bi- 1/211 - exp( - As~21_ 1/2)]
(22)
cAtI i = [(Az/pj_ 1/2)10+ (cAt - Az/#j_ 1/2)Ii_ 1] exp( - As~21_ 1/2) + cAtB~_ 1/211 - exp( - As~2 i_ 1/2)].
(23)
These relations are best understood by considering the characteristic rays on a space-time mesh. They guarantee non-negative intensities since all the coefficients are non-negative. The interpolation relations of equations (20), (21) or (22), (23) are more complicated than the Diamond or Step schemes. They require a test on the mesh parameter ratio and involve the evaluation of exponential functions. Consequently the method is computationally slower than the Diamond or Step schemes. Further, generalization to more dimensions is difficult. In the limits of cArp j_ ~/2/Az large or small the Woods equations give a truncation error which is no better than for Step differencing. This has led us to develop an improved scheme in which a linear variation of the intensity is implied along the edges of the mesh. Subsequently it was discovered that this method was one used by Lathrop in an early version of his code T W O T R A N . The resulting equations are even more complicated than the Woods equations. 5.3 D T F - I V method LATHROP (6) has proposed a simple solution to the negativity problem which is readily generalize d to multidimensional problems. In the computer code DTF-IV, written for
Numerical solutionof the unsteadyradiative transfer equation
955
one dimensional steady neutron transport, he employs Diamond differencing but tests to see whether fluxes are negative. Such values are set zero and the other cell fluxes are recomputed using the basic difference form of the transport equation. This ensures that the conservation property is not lost. This process is repeated as soon as negative values appear. It is clear that in the directions (space or time) where negative fluxes appear under Diamond differencing, interpolation relations of the form of equation (7) are abandoned. 5.4 T W O T R A N method In an alternative approach LATHROP (7) has attempted to keep interpolation relations like equation (7), choosing a, b as close to ½ as possible guaranteeing non-negative values of the flux. By considering limits on a and b imposed by inequalities (10) and (11) he chooses a = Max(½, h)
(24)
b = Max(½, l) where 1 - h = (2Az//~_ 1~2cAt + Az/l~ j_
1/22i_ 1/2)- 1
and
/ . 1 -- 1
(25)
(2/~j_uzcAt/Az + cAt~2 i _ 1/2)- 1
Non-negativity is ensured by guaranteeing that the coefficients in equations (8) and (9) are non-negative. It is important to note that the DTF-IV method can accept negative coefficients provided the resulting intensity is non-negative. Lathrop is dealing with two dimensional steady neutron transport in the computer code TWOTRAN, but the results apply in a straightforward manner to the time dependent problem. 5.5 Geometric mean method We suggest the following method which to our knowledge has not been considered before. We replace the interpolation relations, equation (7), by *he following geometric mean relations
]2 = iili_ 1~ i l i 0 f"
(26)
Substitution of equation (26) in equation (4), the difference form of the transport equation, leads to quadratic equations for the unknown intensities. The quadratic equation has in general one positive and one negative root so that by taking the positive value we ensure a non-negative scheme. 6. A PRELIMINARY ASSESSMENT The various methods described in Section 5 have been analysed in some detail and compared on simple test problems. Only the conclusions are given here.
956
B.C. HANKIN,J. G. HILLand A. G. P. WARHAM
The forward difference or Step method is not recommended because its accuracy generally does not compare well with other methods. However, a hybrid of Diamond and Step differencing (see Section 5.1) is a successful compromise. The methods based on characteristics such as that suggested by Woods are computationally slower than the other methods. They do not seem to gain significantly in accuracy so that any advantage, in convergence schemes for the interaction equation for example, is not apparent. Lathrop's methods in the computer codes DTF-IV and T W O T R A N are attractive since they are reasonably accurate and computationally simple. There is a tendency to oscillatory behaviour in the DTF-IV method which stems from the basic oscillatory behaviour of the Diamond scheme. This is a function of the type of problem and will be discussed later in relation to the radiative transfer problems. The geometric mean method compares unfavourably with the DTF-IV and T W O T R A N methods. It is computationally slower requiring the evaluation of square roots in solving the quadratic equations. Further there is a tendency for strong oscillations arising from a zero or small intensity boundary condition. This problem has been eased with a hybrid Diamond-Geometric Mean scheme. The conclusions given here agree well with those given recently by LATHROP. t4) In the remainder of this paper we shall discuss the methods within the context of radiative transfer problems. 7. RADIATIVE TRANSFER PROBLEMS In radiation-material interaction problems it is not always possible to produce spatial zones that are optically thin. With mean free paths a strong function of temperature and possibly the material in motion, zones may change their opacity considerably with time. It is essential therefore to ensure good behaviour in the diffusion limit for which the SN method was not originally developed. 7.1 Radiative heating o f a semi-infinite s l a b ~ r e s u l t s
Results are given below on the radiative heating of a semi-infinite purely absorbing slab of uniform density and constant mean free path. The equations involved in the solution of this problem are given in Section 2 with the mean free path 2 a constant. Further, an ideal gas equation of state is assumed for the material. Hence E m = pc~T
where cv = specific heat at constant volume. At time zero, radiation from a black-body field of temperature To is allowed to fall on to the slab. We non-dimensionalize the problem by taking z* = z/;~
and t* = ct/2.
Numerical solution of the unsteady radiative transfer equation
957
z* is the distance in mean free paths and t* is the distance in mean free paths travelled by a photon in time t. It is readily shown that the following parameter defines this problem, e = aTg/(pe v + aTe); where e is the fraction of the total energy in the radiation field in equilibrium at temperature To. For the results below e = 0-62826. In Fig. 2 the non-dimensionalized temperature T* ( = T/To) is plotted against z* for t* = 10. Results using the D T F - I V and Step methods are displayed along with an accurate Monte Carlo calculation given by CAMPBELL and NELSON.ts) The D T F - I V method shows considerably better agreement than the Step method. The results for Woods, T W O T R A N and Geometric Mean methods are not significantly different from the D T F - I V values. Az* = 0-5 and At* -- 1.0 for the difference methods.
I'0
0"e
"--']
0"6 T °
i-'] 0-4
"1 0"2
I
I
Z
~ - - - - - ~ - . . . . ~--
3
4
5
6
Zo FIG. 2. T e m p e r a t u r e d i s t r i b u t i o n for the case e = 0-62826, t* = 10. - . . . . M o n t e C a r l o c a l c u l a t i o n of CAMPBELL a n d NELSONtS); D T F - I V method~6); - - Step m e t h o d . Az* = 0-5, At* = 1.0 for the D T F - I V a n d Step calculations.
In Fig. 3 we have a diffusion situation since t* = 10 4. The D T F - I V , T W O T R A N , Geometric Mean and Step methods are compared with an equilibrium diffusion calculation which is essentially exact, given by CAMPBELL.(9) The T W O T R A N and Step methods are indistinguishable. The D T F - I V method is the only one which shows good agreement. Az* = At* = 10 for the non-diffusion methods. For these examples D i a m o n d differencing leads to an immediate collapse of the calculation because of large negative intensities at the front of the temperature wave. 7.2 Radiative heating o f a semi-infinite slab--discussion It is clear in the example above that the most testing situation arises when spatial meshes become optically thick. The only method which survives this test is that of DTF-IV. It is instructive to examine the solution of the transport equation in the diffusion limit in order to understand the failure of some methods.
958
B.C. HANKIN,J. G. HILLand A. G. P. WARHAM I0
~i__________
o.a
____ __~..____t____~________~ L
- .l..L~.a
,_ - ",.... "'b
0-6
t_._
I
i... _ -"1 ,i,
I.. _ _
i
'
I i
L...~ i
t
i
i
20
40
60
80
i
i
t I I00
L
120
t
f
i
140
160
180
Z*
FIG. 3. Temperature distribution for the case e = 0.62826, t* = 104.- - - - Equilibrium diffusion calculation of CAMPBELL~9~;- DTF-IV methodt6~;- - - Step and TWOTRANtT~ methods; - - - - - Geometric Mean method. Az* = At* = 10-0 for the DTF-IV, Step, TWOTRAN and Geometric Mean methods.
In this limit retardation effects are negligible and there is close equilibrium between the radiation a n d material. Hence we m a y m a k e the following a p p r o x i m a t i o n s : 2c-ldl/~t
<< ~
.
(27)
Substitution of a p p r o x i m a t i o n s (27) in equation (1) yields I _~ B - / ~ d B / d r
(28)
where dr
=
dz/2.
The radiative flux is given by [see definition after equation (2)], F -~ (47r/3)dB/dr.
(29)
It is i m p o r t a n t that F is c o m p u t e d correctly by the difference m e t h o d s since this determines the flow of radiation from one cell to the next. In understanding the failure of certain schemes it is pertinent to note that the leading term of equation (28) disappears in the flux integration. It is necessary for the difference m e t h o d s to have at least this property. Essentially we require that the interface intensity Ii should tend to some estimate of the source B on the interface which is independent of
~j- 1/2" If we consider D i a m o n d differencing which is basic to the D T F - I V m e t h o d we find that this condition is fulfilled. However, D i a m o n d relations p r o d u c e negative intensities at the front of the temperature wave so that it is necessary to examine the m e t h o d when D i a m o n d differencing is a b a n d o n e d . Here intensities are set zero and thus the condition required above is satisfied.
Numerical solution of the unsteady radiative transfer equation
959
The T W O T R A N method gives however
1, ___ 8 , - , / 2 ~- Bi+x/2
>
#s-i/2 <
~"
(30)
Here the radiative flux is directly proportional to B t _ x / 2 - B , + l / : and is too large by an amount which increases with optical thickness. T W O T R A N calculations with Az* = 1 show the temperature wave close to where it should be: however, reducing the space and hence time step by a factor 10 increases the computational effort by a factor 100. As an alternative approach the T W O T R A N interpolation relations [equation (7) with equation (24)] could be modified to ensure the correct behaviour in the diffusion limit. For the Step method, approximations (30) apply also so that the flux is predicted too large and the temperature wave travels too quickly. The same conclusion follows for Woods method from an examination of equations (20), (21), (22) and (23); in the diffusion limit the exponentials tend to zero and lead to approximations (30). It has proved possible to modify the Woods method to give the correct diffusion limit, but as previously mentioned, the equations are more complicated. The behaviour of the Geometric Mean method is more difficult to understand. It may be shown that in the optically thick cell limit for/~j_ 1/2 > 0 I, is given by It = B,_ 1/211 + 2/~s_ x/2Arl-_lx/2{1 - (B,_ 1/2/1,_ 1)2}] + O(Azf_2x/2) where Ar t_ ~/z = optical thickness of the cell i - 1. Whereas for #s_ 1/2 < 0 we have 1, = B, + x/2[1 - 2# s- 1/2Az?+xl/z{1 - ( B , + x/z/I, + x)2}] + O(Az~-21/2). It appears therefore that the quadratic terms have a dominating influence otherwise the temperature wave would travel too quickly, which is not in accord with the results above (Fig. 3). 8. OSCILLATORY BEHAVIOUR IN THE DTF-IV METHOD As a result of the conclusions reached in Section 6 and 7 the method of the code DTF-IV is favoured. However, in some problems a large scale oscillatory behaviour has been observed. Typically, radiation from a time dependent source falls onto a transparent region. The photon transit time of this region is comparable with the time scale of the source changes so that retardation effects are important. To demonstrate this effect a simple test problem is given below. 8.1 Test p r o b l e m - results An optically thin slab of constant mean free path is emitting radiation uniformly at temperature To in space and time. On one face radiation falls from a black-body field also at temperature To. At time zero the black-body field is removed and the radiative flux emerging from the other face is observed. The flux is given by the following analytic formula. F* = 1 t* < 2L*~ 1 - 2[E3(2L*)-(2L*/t*)2Ea(t*)] t* > 2L*J
(31)
960
B.C. HANKIN, J. G. HILL and A. G. P. WARHAM
where F* = 4F/ac T'~o, non-dimensionalized flux, acl4 = Stephan's constant, t* = ct/2 L* = L/2
L = half thickness of the slab, and E3(X ) =
third order exponential integral of argument X.
For a time equal to the photon transit time of the slab ( = 2L/c) the removal of the blackbody field cannot affect the flux emerging fr.om the opposite face. For times subsequent to this the flux will fall significantly if the slab is optically thin. From equation (31) we see that the asymptotic value is given by
F*, = 1-2E3(2L* ).
(32)
If the slab is optically thick we get the expected value of unity for F*, since the effect of the opposite face cannot be felt. In the example below L* = 10-2, which gives Fa* = 0.0381, AL* = 10-3 i.e. ten spatial cells and At* = 10 -4. In Fig. 4 results are given for the DTF-IV method under the S 2 and $4 approximations. We see large oscillations of wavelength typically equal to the photon transit time for a cell. The analytic solution is also displayed along with $4 calculations using the Step method and a modified DTF-IV method which is described in Section 8.2 below : these latter $4 calculations are indistinguishable and hence only one curve is seen. The T W O T R A N method gives results, not shown here, very close to those of the Step and modified DTF-IV methods.
1.6
/"x. ,~
1.4
1.2
1.0
?",
"
/
L
,ll!~.~ In i.'I !/\', i i!iY~i:rlll ~.~/ : II 9, i ' ~ 'ii. "' ~"; \i'/ ,, ~",..,; '~'~'"..\ ,:..\\
0"8
\,
/ ", ', ~.
\\
'4.,.
F° 0"6
',
0"4
02
I 001
I 0"02
] 0"03
I 0"04
.........
I 0.05
I 0"06
I 0.07'
t"
FIG. 4. Emergent flux v. time, L* = 0.0l. Analytic solution; - - . - - DTF-IV method <6), S 2 calculation; - - - DTF-IV method, $4 calculation; . . . . Step and modified DTF-IV methods. AL* = 0.001, At* = 0.0001 for the DTF-IV. Step and modified DTF-IV calculations.
Numerical solution of the unsteady radiative transfer equation
961
8.2 Test problem- discussion For this problem no negative intensities tend to arise in the DTF-IV method, so that the oscillatory behaviour is due to Diamond differencing. It appears that the black-body boundary condition is communicated instantly across the slab and, by way of compensation, Diamond differencing results in oscillations about the true solution. An analysis for a single cell (in the half slab) shows the emergent intensity for direction #j_ 1/2 has one maximum at a time approximately equal to the time for a photon travelling in this direction to cross half the cell, and is largely independent of timestep. There is therefore a ray effect arising from the discrete set of#j_ 1/2- For an $4 calculation one might expect two maxima corresponding to the two discrete directions. However, the direction with the largest pj._ 1/2 dominates and only one maximum occurs. A two cell calculation shows two maxima for both $2 and S 4 calculations. The ten cell calculation in Fig. 4 shows ten maxima in the $2 calculation and ten main maxima in the $4 calculation with smaller maxima at late times (t* > 0.02); the effect of two discrete/~j_ 1/2 in the S 4 calculation is noticeable here. Also from Fig. 4 we see that the flux in the $2 calculation begins to fall at a time later than predicted. This is a ray effect since the transit time along the discrete direction is increased by a factor x/3 (pj_ 1/2 = - l/x/3) over that for the forward direction. At early times (t* < 0.02) the mean value o f F * is appreciably larger than unity due to the choice of ogj_ 1/2, #~- 1/2. In order to get the correct forward flux from an isotropic field (in equilibrium I = acT~/4n) we require from equation (6) the quantity 4 ~ o9~_1/2Pj-1/2 equal to unity when summed over positive p j_ 1/2. For the S 2 and $4 calculations this quantity equals 1.1547 and 1.0328 respectively. As mentioned earlier the time interval between the maxima is the order of the photon transit time for a spatial zone. In these calculations AL*/At* = 10 so that many timesteps are involved in defining these oscillations. For AL*/At* approximately unity the time interval would be of the order of the timestep. It is clear from our studies that by using more spatial cells and a higher angular approximation a better numerical solution is obtained. However, the obvious penalty is an increase in calculation time. The Step and T W O T R A N results show no oscillations and considering the low angular approximation ($4) the agreement is quite good. It can be shown that the T W O T R A N method is using Diamond differencing in time and very closely Step differencing in space. Hence, since the timestep is small, there is little difference between the Step and T W O T R A N results. As previously noted the Step method has generally a large truncation error and the T W O T R A N method is incorrect in the diffusion limit. The possibility of extending the T W O T R A N method to cover optically thick cell problems has also been noted. An analysis of the negativity problem (see Section 4) for one dimensional steady transport has suggested a modification to the DTF-IV method which may also ease the oscillation problem. A method is currently being developed which is based on Diamond differencing but tests to see whether interface intensities in space and time are less than some source value for the mesh. If they are, they are set equal to this source value and the other intensities recomputed following the principle of the DTF-IV method. In the spatial direction this source is found by interpolating on to the interface assuming it is linear with optical depth. Currently, the source value in the time direction is taken equal to the value used in that timestep: the possibility of using an extrapolated value is being considered. This modified DTF-IV
962
B.C. HANKIN,J. G. HILL and A. G. P. WARHAM
scheme has been tested on a number of problems and gives results close to those of the normal DTF-IV scheme except for the example of the present section. For the radiative slab heating problem in the diffusion limit (see Section 7), it is essential that the test source in the spatial direction is an interpolated value as described above, otherwise the speed of the thermal wave is incorrect. For the example of the present Section this method shows smooth results but more recent calculations, not reported here, are not so encouraging. CONCLUSIONS
The present work has examined variations on the discrete SN method for solving the radiative transfer equation in radiation-material interaction problems. A preliminary assessment indicates that the method used by Lathrop for one dimensional steady neutron transport in the computer code DTF-IV, fulfills the main requirements of energy conservation, non-negativity of intensities, reasonable accuracy and computational simplicity. The behaviour of the methods in the diffusion limit is important in this assessment. It is suggested that the failure of the TWOTRAN method in this limit can be overcome. The work has been extended in an attempt to understand the oscillatory behaviour of the DTF-IV method in certain problems. A simple example shows this phenomenon arising from Diamond differencing alone. Some preliminary calculations using Diamond differencing and a more stringent application of the testing procedure in the DTF-IV method are encouraging. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9.
B. G. CARLSON,Meth. ComputationalPhys. 1, 1 (1963). I. P. GRANT, Mon. Not. R. Astr. Soc. 125, 417 (1963). D. WooDs, Private Communication. (1968). K. D. LATHROP,J. comp. Phys. 4, 475 (1969). B. WENDROFF,J. comp. Phys. 4, 211 (1969). K. D. LATHROP,LA-3373 (1965). K. D. LATHROP, LA-4058 (1969). P. M. CAMPBELLand R. NELSON, UCRL-7838 (1964). P. M. CAMPnELL,UCRL-12411 (1965).
NUMERICAL SOLUTION OF THE UNSTEADY RADIATIVE TRANSFER EQUATION B. C. HANKIN, J. G. HILL and A. G. P. WARHAM Mathematical Physics Department, A.W.R.E., Aldermaston, Berks. Abstract--In problems involving the transport of radiation through matter, it is necessary to solve the radiative transfer equation numerically. The scope for analytic solutions is restricted because of the complexity of the radiation-material interaction problem; this is particularly severe if the material is compressible. A study is made of various difference schemes within the framework of the discrete SN method tl~ with a view to obtaining a well-behaved solution. In addition to the usual requirements of consistency, stability and accuracy, there are further desirable features. Firstly, it is preferable to keep in the difference scheme, the photon conservation property as expressed by the basic transport equation. Secondly, the solution should be non-negative, not only because such values are unphysical, but also because of their effect on the radiation-material coupling. Such a coupling is usually non-linear, and negative solutions can lead to instabilities in the numerical solution. Finally, it is preferable that the method should lead to simple calculations. Results are given for an example involving the radiative heating of a semi-infinite slab. It is shown that the diffusion limit imposes a particularly severe test, and the method used by Lathrop in the computer code DTF-IV ~6) is favoured. A second example demonstrates a disturbing oscillatory behaviour arising from Diamond differencing which is basic to the DTF-IV method. A modification to this method designed to overcome this problem is examined. 1. I N T R O D U C T I O N
THE TRANSFER of radiation through matter is described by a linear Boltzmann equation when the photon-photon interaction is neglected. In practical problems the source function and mean free path often depend on material temperature, so that this equation is coupled with an equation describing the radiation-material interaction. This coupling can be very non-linear, reducing considerably the scope for analytic solutions, so that numerical methods have to be used. An added complexity arises if the material moves hydrodynamically. The present paper examines variations of the finite difference method known as the discrete SN method31) This method is chosen because it fulfills the usual requirements for a difference scheme of consistency and stability, and embodies the photon conservation property which is fundamental to the original transport equation. Much work has been done on the discrete SN method in the context of neutron transport. Such work has been aimed broadly at getting an accurate non-negative scheme. Our experience on radiation problems has shown that it is highly desirable that photon specific intensities are computed non-negative. Clearly negative intensities are unphysical but it is their consequences that can be more serious: the non-linear coupling between the transport equation and the photon-material interaction equation can lead to instabilities. F u r t h e r G R A N T (2) has pointed out that when, for example, mean free paths have a strong temperature dependence, it is necessary to ensure that the optically thick or diffusion limit is adequately treated. 949
950
B . C . HANKIN, J. G . HILL a n d A . G . P. WARHAM
There are also difficulties in calculating unsteady transport or multi-space dimensional transport in optically thin situations. These difficulties are common to both neutron and photon studies. 2. THE RADIATION-MATERIAL INTERACTION PROBLEM Let us consider the equation describing the plane unsteady transport of grey radiation in a purely absorbing medium in which local thermodynamic equilibrium exists. The transfer equation is c - xc3I/~t + p~3I/~3z + 2 - 11 = 2 - XB (1) where I = I(z, lz, t), the specific intensity at position z in direction cos- ~/~ at time t, c 2 B p T
= = = = =
velocity of light, 2(p, T) grey absorption mean free path, (ac/4rOT 4, the black body emission function, material density, material temperature
and ac/4 = Stefan's constant. If we integrate equation (1) over/~ we get the radiation energy equation ~3E/Ot + ~F/~z = 4rc2-1 [ B - (c/47t)E]
(2)
where +1
E = (2g/c) f
I d/~, the radiation energy density
-1
and
+
F = 2~ /
/~I d/~, the radiative flux.
1/ -1
The right hand side of equation (2) gives the difference between the emissions and absorptions. By equating this to the material cooling rate we get the following equation for the radiation-material interaction, - t3Em/8t = 4rc2- I ~ B - (c/4rc)E]
(3)
where E m = E,,(p, T), material energy density. No hydrodynamic effects are included in equation (3). The interaction problem consists of a solution of equations (1) and (3) with appropriate boundary conditions. E m is related to p and T by means of an equation of state. 2 is also a known function of p and T. The dependence of E m and 2 on p and T may exist in tabular form. 3. THE DISCRETE SN METHOD The discrete SN form for equation (1) is written (11 -- I°)/cA t + # j - 112([i -- I i - 1)/AZ + ~i---t1/2] =
~.Z-11/2Bi- 1/2
(4)
N u m e r i c a l s o l u t i o n o f the u n s t e a d y r a d i a t i v e transfer e q u a t i o n
951
where At, Az are the time and space mesh sizes respectively and 11, i0, I+ and Ii_ 1 are component intensities representing the flow of radiation energy across the space-time mesh as illustrated in Fig. 1. i is the mean intensity in the cell in the discrete direction #j_ 1/2 and determines the total absorption in the cell. 2+_ 1/2, B~_ 1/2 are average values of 2 and B in the time interval.
t I i I At
Ii-i
z
I°
FIG. 1. D e f i n i t i o n o f cell intensities.
If we multiply equation (4) by a weighting factor ogg_1/2 and sum over all j we get the following difference form of equation (2), the radiation energy equation : (E 1 -
E°)/At + (F i - Fi_ 1)/Az = 41z}tz_lu2[Bi_
1/2 -
(c/47r)E]
I5)
where E ° = (4n/c) ~ coj_ 1/2/°, similarly E 1 and E, i F i - 1 = 4~ E (Dj_ 1/21gj_ 1 / 2 1 / _ 1 and similarly Fi.
(6)
J
Equations (6) provide quadrature formulae for E and F depending on the choice of 09j_ 1/2, P j - 1/2. The more difficult problem of curved geometry is readily included formally in the discrete SN method. Considerable effort has been devoted to the choice of 09j_ 1/2, P j - 1/2 in the neutron context31) Integration over any space region for any time interval is readily applied to equation (5) thus ensuring the conservation property. Equation (4) introduces five unknown intensities so we require four more equations to complete the solution. Two equations are provided by the boundary conditions in the z and t directions so that two more equations are needed. These are the interpolation or
952
B.C. HANKIN, J. G. HILLand A. G. P. WARHAra
auxiliary relations which may be written in the following simple form :
i = ali+(1--a)li_l~ bI 1 + ( 1 - b ) I ° J
(7)
w h e r e 0 < a , b < 1. If we take a and b equal to ½ we complete the Diamond difference scheme. The method of solving equations (4) and (7) is described by CARLSON.(1)
4. THE NEGATIVITY PROBLEM For /~_ 1/2 positive we may regard 1i_ 1, I° as known (initially from the boundary conditions). If we eliminate i and 11 from equations (4) and (7) we get for the remaining unknown 11, (1 + aAz/bltj_ 1/2cAt + aAz/pj_
1/2~.i_ 1/2)Ii =
[1 - (1 - a)Az/bl~j_ 1/2cAt
- (1 - a)Az/#j_ 1/221-1/2]Ii- 1 + (Az/bpj_ 1/2cAt)I ° + (Az/t~j_ 1/22i- 1/2)Bi- 1/2.
(8)
Similarly for the unknown I1 we get (1 + bpj_ 1/2cAt/aAz + bcAt/2 i_ 1/2)11 = [1 - (1 - b)llj_ 1/2cAt/aAz - (1 - b)cAt/2 i_ 1/2110
+ (l~j_ 1/2cAt/aAz)I i_ 1 + (cAt/2i- 1/2)Bi- 1/2-
(9)
Clearly 11 and 11 are guaranteed non-negative if all the coefficients in equations (8) and (9) are positive, i is always guaranteed non-negative. Hence it is sufficient for (1 - a) < (Az/bltj_ 1~2cAt q- Az/#j_
i/2~.i_ 1/2)- 1
(1 -- b) < (Itj_ 1/2cAt/aAz + cat~2 I_ 1/2)- 1.
(10) (1 i)
The only values of a and b which satisfy inequalities (10) and (11) for all values are a = b = 1. Referring to equation (7) this corresponds to a forward difference or Step scheme. It is important to emphasize that inequalities (10) and (11) are sufficient but they are not in general necessary to give non-negativity. If we consider Diamond differencing in optically thin situations (AT,/J.i_l/2 small) relations (10) and (11) lead to
Hi_ 1/2cAt/Az _> 1
(12)
Az/laj_l/2CAt > 1.
(13)
and
Relations (12) and (13) are mutually exclusive except for the special case of the equality. Negative intensities will have a tendency therefore to arise if the factor !ui_ ~/2cAt/Az departs significantly from unity. In radiation-material interaction problems the timestep is often limited by the convergence in the solution of the interaction equation [equation (3)]. Hence inequality (12) tends to be violated giving rise to oscillations on the spatial mesh.
Numerical solution of the unsteady radiative transfer equation
953
It is clear also that negative values will tend to appear with Diamond differencing in optically thick problems (Az/2i_ 1/2 large). This problem has to be faced in radiative transfer studies. 5. SOLUTIONS TO THE NEGATIVITY PROBLEM The difficulty described in the previous section was recognized early in the development of the SN method and techniques for overcoming this were produced31~ In this Section these and more recent methods are described and a preliminary assessment is given in Section 6.
5.1 Forward difference or Step method This method was introduced in Section 4 above. The interpolation relations of equation (7) become ] = Ii = ll- 1 =I1
tgj_l{2 > O]
1~2- 1/2 < 0 i'"
(14)
For stability reasons the Solution of the SN equations proceeds for z increasing if/as _ 1/2 is positive and z decreasing for/22_ 1/2 negative. As shown in Section 4 above, the Step method guarantees non-negative intensities whatever the physical parameter values or mesh size. However, its truncation error is generally larger than that for D i a m o n d differencing. For this reason one of the earliest solutions ~1) to the negativity problem was a scheme combining the D i a m o n d and Step methods. Here a test is made to see whether the intensity (or angular flux in the neutron problem) is negative. If it is, the intensity is recalculated using a Step relation. 5.2 Woods method Many solutions have been attempted using a characteristics approach. The first simple application within the discrete S N framework is due to WOODSt3~ and has been described by LATHROP34J We give below a derivation of our own which follows WENDROFF.ts~ Wendroff's scheme does not inhibit negative intensities. We first rewrite equation (1) in characteristic form,
d l / d s + 2 - lI = 2 - 1 B
(15)
where dt/ds = c-1, dz/ds = p and s is the characteristic coordinate. If we assume 2, B constant within a space-time mesh we m a y write for the solution of equation (15) on the mesh I = I o exp( - s/2) + B[1 - e x p ( - s/2)]
(16)
where Io = intensity at s = 0. Let us now consider equation (4), the discrete S N form of the transport equation. We may write this as follows
( I k - - l k - 1)/As + ).7-11/21 = 27-11/2 Bi- 1/2
(17)
954
B.C. HANKIN, J. G. HILL and A. G. P. WARHAM
where
Ik- 1 = (I°/cat + #j- 1/2Ii - 1/Az)/As ] Ik = (I1/cAt + l~j- 1/2Ii/Az)/As
t"
(18)
and
As = cAtAz/(#j_ 1~2cAt + Az) We regard equation (17) as a difference form of equation (15) which implies that lk_ 1,
Ik lie on the same characteristic direction. Accordingly we use equation (16) to relate Ik_ 1 and lk, Ik = Ik- 1 exp(-- As/2i_ 1/2) + Bi-
1/2[ 1 - -
e x p ( - As~2 i_ 1/2)].
(19)
We are treating here the positive p direction. Equation (19) is decomposed by Woods into a pair of equations. The decomposition however depends on the mesh parameters. For cAt < Az/pj_ 1/2 we have
(Az/pj_ 1/2)11 = [(Az/pj_ 1/2 - cAt) I° + cAtli- 1] e x p ( - As~21_ 1/2) q- (Az//2j_ 1/2)Bi_ 1/211
-
e x p ( - As/,~ i_ 1/2)]
(20)
cAtI i = cAtI ° e x p ( - As/21_ 1/2) + cAtBi - x/2[ 1 - exp( - As/2i_ 1/2)]
(21)
and for cat > Az/p i_ 1/2 we have
(Az/ktj_ 1/2)11 = (Az/pj_ 1/2)Ii_ ~ e x p ( - As/2i_ 1/2) + (Az/p;_ 1/2)Bi- 1/211 - exp( - As~21_ 1/2)]
(22)
cAtI i = [(Az/pj_ 1/2)10+ (cAt - Az/#j_ 1/2)Ii_ 1] exp( - As~21_ 1/2) + cAtB~_ 1/211 - exp( - As~2 i_ 1/2)].
(23)
These relations are best understood by considering the characteristic rays on a space-time mesh. They guarantee non-negative intensities since all the coefficients are non-negative. The interpolation relations of equations (20), (21) or (22), (23) are more complicated than the Diamond or Step schemes. They require a test on the mesh parameter ratio and involve the evaluation of exponential functions. Consequently the method is computationally slower than the Diamond or Step schemes. Further, generalization to more dimensions is difficult. In the limits of cArp j_ ~/2/Az large or small the Woods equations give a truncation error which is no better than for Step differencing. This has led us to develop an improved scheme in which a linear variation of the intensity is implied along the edges of the mesh. Subsequently it was discovered that this method was one used by Lathrop in an early version of his code T W O T R A N . The resulting equations are even more complicated than the Woods equations. 5.3 D T F - I V method LATHROP (6) has proposed a simple solution to the negativity problem which is readily generalize d to multidimensional problems. In the computer code DTF-IV, written for
Numerical solutionof the unsteadyradiative transfer equation
955
one dimensional steady neutron transport, he employs Diamond differencing but tests to see whether fluxes are negative. Such values are set zero and the other cell fluxes are recomputed using the basic difference form of the transport equation. This ensures that the conservation property is not lost. This process is repeated as soon as negative values appear. It is clear that in the directions (space or time) where negative fluxes appear under Diamond differencing, interpolation relations of the form of equation (7) are abandoned. 5.4 T W O T R A N method In an alternative approach LATHROP (7) has attempted to keep interpolation relations like equation (7), choosing a, b as close to ½ as possible guaranteeing non-negative values of the flux. By considering limits on a and b imposed by inequalities (10) and (11) he chooses a = Max(½, h)
(24)
b = Max(½, l) where 1 - h = (2Az//~_ 1~2cAt + Az/l~ j_
1/22i_ 1/2)- 1
and
/ . 1 -- 1
(25)
(2/~j_uzcAt/Az + cAt~2 i _ 1/2)- 1
Non-negativity is ensured by guaranteeing that the coefficients in equations (8) and (9) are non-negative. It is important to note that the DTF-IV method can accept negative coefficients provided the resulting intensity is non-negative. Lathrop is dealing with two dimensional steady neutron transport in the computer code TWOTRAN, but the results apply in a straightforward manner to the time dependent problem. 5.5 Geometric mean method We suggest the following method which to our knowledge has not been considered before. We replace the interpolation relations, equation (7), by *he following geometric mean relations
]2 = iili_ 1~ i l i 0 f"
(26)
Substitution of equation (26) in equation (4), the difference form of the transport equation, leads to quadratic equations for the unknown intensities. The quadratic equation has in general one positive and one negative root so that by taking the positive value we ensure a non-negative scheme. 6. A PRELIMINARY ASSESSMENT The various methods described in Section 5 have been analysed in some detail and compared on simple test problems. Only the conclusions are given here.
956
B.C. HANKIN,J. G. HILLand A. G. P. WARHAM
The forward difference or Step method is not recommended because its accuracy generally does not compare well with other methods. However, a hybrid of Diamond and Step differencing (see Section 5.1) is a successful compromise. The methods based on characteristics such as that suggested by Woods are computationally slower than the other methods. They do not seem to gain significantly in accuracy so that any advantage, in convergence schemes for the interaction equation for example, is not apparent. Lathrop's methods in the computer codes DTF-IV and T W O T R A N are attractive since they are reasonably accurate and computationally simple. There is a tendency to oscillatory behaviour in the DTF-IV method which stems from the basic oscillatory behaviour of the Diamond scheme. This is a function of the type of problem and will be discussed later in relation to the radiative transfer problems. The geometric mean method compares unfavourably with the DTF-IV and T W O T R A N methods. It is computationally slower requiring the evaluation of square roots in solving the quadratic equations. Further there is a tendency for strong oscillations arising from a zero or small intensity boundary condition. This problem has been eased with a hybrid Diamond-Geometric Mean scheme. The conclusions given here agree well with those given recently by LATHROP. t4) In the remainder of this paper we shall discuss the methods within the context of radiative transfer problems. 7. RADIATIVE TRANSFER PROBLEMS In radiation-material interaction problems it is not always possible to produce spatial zones that are optically thin. With mean free paths a strong function of temperature and possibly the material in motion, zones may change their opacity considerably with time. It is essential therefore to ensure good behaviour in the diffusion limit for which the SN method was not originally developed. 7.1 Radiative heating o f a semi-infinite s l a b ~ r e s u l t s
Results are given below on the radiative heating of a semi-infinite purely absorbing slab of uniform density and constant mean free path. The equations involved in the solution of this problem are given in Section 2 with the mean free path 2 a constant. Further, an ideal gas equation of state is assumed for the material. Hence E m = pc~T
where cv = specific heat at constant volume. At time zero, radiation from a black-body field of temperature To is allowed to fall on to the slab. We non-dimensionalize the problem by taking z* = z/;~
and t* = ct/2.
Numerical solution of the unsteady radiative transfer equation
957
z* is the distance in mean free paths and t* is the distance in mean free paths travelled by a photon in time t. It is readily shown that the following parameter defines this problem, e = aTg/(pe v + aTe); where e is the fraction of the total energy in the radiation field in equilibrium at temperature To. For the results below e = 0-62826. In Fig. 2 the non-dimensionalized temperature T* ( = T/To) is plotted against z* for t* = 10. Results using the D T F - I V and Step methods are displayed along with an accurate Monte Carlo calculation given by CAMPBELL and NELSON.ts) The D T F - I V method shows considerably better agreement than the Step method. The results for Woods, T W O T R A N and Geometric Mean methods are not significantly different from the D T F - I V values. Az* = 0-5 and At* -- 1.0 for the difference methods.
I'0
0"e
"--']
0"6 T °
i-'] 0-4
"1 0"2
I
I
Z
~ - - - - - ~ - . . . . ~--
3
4
5
6
Zo FIG. 2. T e m p e r a t u r e d i s t r i b u t i o n for the case e = 0-62826, t* = 10. - . . . . M o n t e C a r l o c a l c u l a t i o n of CAMPBELL a n d NELSONtS); D T F - I V method~6); - - Step m e t h o d . Az* = 0-5, At* = 1.0 for the D T F - I V a n d Step calculations.
In Fig. 3 we have a diffusion situation since t* = 10 4. The D T F - I V , T W O T R A N , Geometric Mean and Step methods are compared with an equilibrium diffusion calculation which is essentially exact, given by CAMPBELL.(9) The T W O T R A N and Step methods are indistinguishable. The D T F - I V method is the only one which shows good agreement. Az* = At* = 10 for the non-diffusion methods. For these examples D i a m o n d differencing leads to an immediate collapse of the calculation because of large negative intensities at the front of the temperature wave. 7.2 Radiative heating o f a semi-infinite slab--discussion It is clear in the example above that the most testing situation arises when spatial meshes become optically thick. The only method which survives this test is that of DTF-IV. It is instructive to examine the solution of the transport equation in the diffusion limit in order to understand the failure of some methods.
958
B.C. HANKIN,J. G. HILLand A. G. P. WARHAM I0
~i__________
o.a
____ __~..____t____~________~ L
- .l..L~.a
,_ - ",.... "'b
0-6
t_._
I
i... _ -"1 ,i,
I.. _ _
i
'
I i
L...~ i
t
i
i
20
40
60
80
i
i
t I I00
L
120
t
f
i
140
160
180
Z*
FIG. 3. Temperature distribution for the case e = 0.62826, t* = 104.- - - - Equilibrium diffusion calculation of CAMPBELL~9~;- DTF-IV methodt6~;- - - Step and TWOTRANtT~ methods; - - - - - Geometric Mean method. Az* = At* = 10-0 for the DTF-IV, Step, TWOTRAN and Geometric Mean methods.
In this limit retardation effects are negligible and there is close equilibrium between the radiation a n d material. Hence we m a y m a k e the following a p p r o x i m a t i o n s : 2c-ldl/~t
<< ~
.
(27)
Substitution of a p p r o x i m a t i o n s (27) in equation (1) yields I _~ B - / ~ d B / d r
(28)
where dr
=
dz/2.
The radiative flux is given by [see definition after equation (2)], F -~ (47r/3)dB/dr.
(29)
It is i m p o r t a n t that F is c o m p u t e d correctly by the difference m e t h o d s since this determines the flow of radiation from one cell to the next. In understanding the failure of certain schemes it is pertinent to note that the leading term of equation (28) disappears in the flux integration. It is necessary for the difference m e t h o d s to have at least this property. Essentially we require that the interface intensity Ii should tend to some estimate of the source B on the interface which is independent of
~j- 1/2" If we consider D i a m o n d differencing which is basic to the D T F - I V m e t h o d we find that this condition is fulfilled. However, D i a m o n d relations p r o d u c e negative intensities at the front of the temperature wave so that it is necessary to examine the m e t h o d when D i a m o n d differencing is a b a n d o n e d . Here intensities are set zero and thus the condition required above is satisfied.
Numerical solution of the unsteady radiative transfer equation
959
The T W O T R A N method gives however
1, ___ 8 , - , / 2 ~- Bi+x/2
>
#s-i/2 <
~"
(30)
Here the radiative flux is directly proportional to B t _ x / 2 - B , + l / : and is too large by an amount which increases with optical thickness. T W O T R A N calculations with Az* = 1 show the temperature wave close to where it should be: however, reducing the space and hence time step by a factor 10 increases the computational effort by a factor 100. As an alternative approach the T W O T R A N interpolation relations [equation (7) with equation (24)] could be modified to ensure the correct behaviour in the diffusion limit. For the Step method, approximations (30) apply also so that the flux is predicted too large and the temperature wave travels too quickly. The same conclusion follows for Woods method from an examination of equations (20), (21), (22) and (23); in the diffusion limit the exponentials tend to zero and lead to approximations (30). It has proved possible to modify the Woods method to give the correct diffusion limit, but as previously mentioned, the equations are more complicated. The behaviour of the Geometric Mean method is more difficult to understand. It may be shown that in the optically thick cell limit for/~j_ 1/2 > 0 I, is given by It = B,_ 1/211 + 2/~s_ x/2Arl-_lx/2{1 - (B,_ 1/2/1,_ 1)2}] + O(Azf_2x/2) where Ar t_ ~/z = optical thickness of the cell i - 1. Whereas for #s_ 1/2 < 0 we have 1, = B, + x/2[1 - 2# s- 1/2Az?+xl/z{1 - ( B , + x/z/I, + x)2}] + O(Az~-21/2). It appears therefore that the quadratic terms have a dominating influence otherwise the temperature wave would travel too quickly, which is not in accord with the results above (Fig. 3). 8. OSCILLATORY BEHAVIOUR IN THE DTF-IV METHOD As a result of the conclusions reached in Section 6 and 7 the method of the code DTF-IV is favoured. However, in some problems a large scale oscillatory behaviour has been observed. Typically, radiation from a time dependent source falls onto a transparent region. The photon transit time of this region is comparable with the time scale of the source changes so that retardation effects are important. To demonstrate this effect a simple test problem is given below. 8.1 Test p r o b l e m - results An optically thin slab of constant mean free path is emitting radiation uniformly at temperature To in space and time. On one face radiation falls from a black-body field also at temperature To. At time zero the black-body field is removed and the radiative flux emerging from the other face is observed. The flux is given by the following analytic formula. F* = 1 t* < 2L*~ 1 - 2[E3(2L*)-(2L*/t*)2Ea(t*)] t* > 2L*J
(31)
960
B.C. HANKIN, J. G. HILL and A. G. P. WARHAM
where F* = 4F/ac T'~o, non-dimensionalized flux, acl4 = Stephan's constant, t* = ct/2 L* = L/2
L = half thickness of the slab, and E3(X ) =
third order exponential integral of argument X.
For a time equal to the photon transit time of the slab ( = 2L/c) the removal of the blackbody field cannot affect the flux emerging fr.om the opposite face. For times subsequent to this the flux will fall significantly if the slab is optically thin. From equation (31) we see that the asymptotic value is given by
F*, = 1-2E3(2L* ).
(32)
If the slab is optically thick we get the expected value of unity for F*, since the effect of the opposite face cannot be felt. In the example below L* = 10-2, which gives Fa* = 0.0381, AL* = 10-3 i.e. ten spatial cells and At* = 10 -4. In Fig. 4 results are given for the DTF-IV method under the S 2 and $4 approximations. We see large oscillations of wavelength typically equal to the photon transit time for a cell. The analytic solution is also displayed along with $4 calculations using the Step method and a modified DTF-IV method which is described in Section 8.2 below : these latter $4 calculations are indistinguishable and hence only one curve is seen. The T W O T R A N method gives results, not shown here, very close to those of the Step and modified DTF-IV methods.
1.6
/"x. ,~
1.4
1.2
1.0
?",
"
/
L
,ll!~.~ In i.'I !/\', i i!iY~i:rlll ~.~/ : II 9, i ' ~ 'ii. "' ~"; \i'/ ,, ~",..,; '~'~'"..\ ,:..\\
0"8
\,
/ ", ', ~.
\\
'4.,.
F° 0"6
',
0"4
02
I 001
I 0"02
] 0"03
I 0"04
.........
I 0.05
I 0"06
I 0.07'
t"
FIG. 4. Emergent flux v. time, L* = 0.0l. Analytic solution; - - . - - DTF-IV method <6), S 2 calculation; - - - DTF-IV method, $4 calculation; . . . . Step and modified DTF-IV methods. AL* = 0.001, At* = 0.0001 for the DTF-IV. Step and modified DTF-IV calculations.
Numerical solution of the unsteady radiative transfer equation
961
8.2 Test problem- discussion For this problem no negative intensities tend to arise in the DTF-IV method, so that the oscillatory behaviour is due to Diamond differencing. It appears that the black-body boundary condition is communicated instantly across the slab and, by way of compensation, Diamond differencing results in oscillations about the true solution. An analysis for a single cell (in the half slab) shows the emergent intensity for direction #j_ 1/2 has one maximum at a time approximately equal to the time for a photon travelling in this direction to cross half the cell, and is largely independent of timestep. There is therefore a ray effect arising from the discrete set of#j_ 1/2- For an $4 calculation one might expect two maxima corresponding to the two discrete directions. However, the direction with the largest pj._ 1/2 dominates and only one maximum occurs. A two cell calculation shows two maxima for both $2 and S 4 calculations. The ten cell calculation in Fig. 4 shows ten maxima in the $2 calculation and ten main maxima in the $4 calculation with smaller maxima at late times (t* > 0.02); the effect of two discrete/~j_ 1/2 in the S 4 calculation is noticeable here. Also from Fig. 4 we see that the flux in the $2 calculation begins to fall at a time later than predicted. This is a ray effect since the transit time along the discrete direction is increased by a factor x/3 (pj_ 1/2 = - l/x/3) over that for the forward direction. At early times (t* < 0.02) the mean value o f F * is appreciably larger than unity due to the choice of ogj_ 1/2, #~- 1/2. In order to get the correct forward flux from an isotropic field (in equilibrium I = acT~/4n) we require from equation (6) the quantity 4 ~ o9~_1/2Pj-1/2 equal to unity when summed over positive p j_ 1/2. For the S 2 and $4 calculations this quantity equals 1.1547 and 1.0328 respectively. As mentioned earlier the time interval between the maxima is the order of the photon transit time for a spatial zone. In these calculations AL*/At* = 10 so that many timesteps are involved in defining these oscillations. For AL*/At* approximately unity the time interval would be of the order of the timestep. It is clear from our studies that by using more spatial cells and a higher angular approximation a better numerical solution is obtained. However, the obvious penalty is an increase in calculation time. The Step and T W O T R A N results show no oscillations and considering the low angular approximation ($4) the agreement is quite good. It can be shown that the T W O T R A N method is using Diamond differencing in time and very closely Step differencing in space. Hence, since the timestep is small, there is little difference between the Step and T W O T R A N results. As previously noted the Step method has generally a large truncation error and the T W O T R A N method is incorrect in the diffusion limit. The possibility of extending the T W O T R A N method to cover optically thick cell problems has also been noted. An analysis of the negativity problem (see Section 4) for one dimensional steady transport has suggested a modification to the DTF-IV method which may also ease the oscillation problem. A method is currently being developed which is based on Diamond differencing but tests to see whether interface intensities in space and time are less than some source value for the mesh. If they are, they are set equal to this source value and the other intensities recomputed following the principle of the DTF-IV method. In the spatial direction this source is found by interpolating on to the interface assuming it is linear with optical depth. Currently, the source value in the time direction is taken equal to the value used in that timestep: the possibility of using an extrapolated value is being considered. This modified DTF-IV
962
B.C. HANKIN,J. G. HILL and A. G. P. WARHAM
scheme has been tested on a number of problems and gives results close to those of the normal DTF-IV scheme except for the example of the present section. For the radiative slab heating problem in the diffusion limit (see Section 7), it is essential that the test source in the spatial direction is an interpolated value as described above, otherwise the speed of the thermal wave is incorrect. For the example of the present Section this method shows smooth results but more recent calculations, not reported here, are not so encouraging. CONCLUSIONS
The present work has examined variations on the discrete SN method for solving the radiative transfer equation in radiation-material interaction problems. A preliminary assessment indicates that the method used by Lathrop for one dimensional steady neutron transport in the computer code DTF-IV, fulfills the main requirements of energy conservation, non-negativity of intensities, reasonable accuracy and computational simplicity. The behaviour of the methods in the diffusion limit is important in this assessment. It is suggested that the failure of the TWOTRAN method in this limit can be overcome. The work has been extended in an attempt to understand the oscillatory behaviour of the DTF-IV method in certain problems. A simple example shows this phenomenon arising from Diamond differencing alone. Some preliminary calculations using Diamond differencing and a more stringent application of the testing procedure in the DTF-IV method are encouraging. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9.
B. G. CARLSON,Meth. ComputationalPhys. 1, 1 (1963). I. P. GRANT, Mon. Not. R. Astr. Soc. 125, 417 (1963). D. WooDs, Private Communication. (1968). K. D. LATHROP,J. comp. Phys. 4, 475 (1969). B. WENDROFF,J. comp. Phys. 4, 211 (1969). K. D. LATHROP,LA-3373 (1965). K. D. LATHROP, LA-4058 (1969). P. M. CAMPBELLand R. NELSON, UCRL-7838 (1964). P. M. CAMPnELL,UCRL-12411 (1965).