- Email: [email protected]
Relating Eddington factors to flux limiters
J. Quant. Spectrosc. Radiat. Transfer Vol. 31, No. 2, pp. 149-160, 1984 Printed in Great Britain.
0022-4073/84 $3.00+ .00 Pergamon Press Ltd.
RELATING EDDINGTON FACTORS TO FLUX LIMITERSt C. D. LEVERMORE A-Division, Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, CA 94550, U.S.A.
(Reveived 6 May 1983) Abstract--Variable Eddington factors and flux-limiters have been introduced in the P-1 and diffusion equations, respectively, to handle sittiations when the underlying intensity is too anisotropic for the unmodified theories to remain valid. We present a derivation of a relation between the two for which a new approach to the diffusion approximation is used, Algebraic expressions for Eddington factors satisfying the moment conditions are not satisfactory for closing the P-1 equations but, by using the derived relation, yield acceptable flux-limited diffusion theories.
1. I N T R O D U C T I O N
In many problems involving transport theory, resource constraints may necessitate adoption of some approximate method to calculate the transport phenomena. Two such methods are the P-1 and diffusion approximations, both of which require a degree of isotropy in the underlying intensity. In order to handle situations with significant anisotropy, variable Eddington factors have been introduced for the P-1 approximation 1-4 and flux-limiters for the diffusion approximation. 5-8While it has been understood that these two ideas are related, no explicit relationship has been developed until now. Assuming coherent scattering and stationary matter, the radiation transport equation for the specific intensity of photons, I(1~, r, t) is 9'1° (1/c)8tl + ~ • g I + trI = (c/4x)trAB + as f g(I2 • ~ ' ) I(O') dO'. d4~z
(1)
Here t2, r, and t are the angular, spatial and temporal variables, c is the speed of light, aA(r, t) is the absorption coefficient, trs(r, t) is the scattering coefficient, a = aA + as is the total interaction coefficient, B(r, t) is the black body energy density, and g(f2 • l~', r, t) is the scattering angular redistribution function which is normalized by 1=
g(~
• ~ ' ) dQ = 2x
fl
g ( ~ ) d/~.
These equations may be considered to be either frequency dependent with the frequency parameters suppressed or gray equations. Integrating Eq. (1) over the solid angle dQ yields 8,E + V " F + ctrA ( E -
B ) = O,
(2)
where the energy density, E(r, t), and the energy flux, F(r, t) are defined in terms of the zeroth and first moments of the intensity by E(r, t) = (l/c) f4~ I(1"2, r, t) dO
(3)
tThis work was performed under the auspices of the U.S. Department of Energy at the Lawrence Livermore National Laboratory under contract number W-7405-ENG-48. 149
150
C.D. LEVERMORE
and F(r, t) = I OI(I'L r, t) dO. d4=
(4)
An equation for the flux is obtained by multiplying Eq. (1) by fl and then integrating over the solid angle dQ, giving (1/c)dtF + c F • P + (a - glas)F = 0,
(5)
where g~(r, t) is the first moment of the angular redistribution function gl = 2re
f
l
#g(/t) d# 1
and the radiation pressure tensor, P(r, t), is defined in terms of the second moment of the intensity by
P(r, t) = (l/c) f4, 1~12I(12, r, t) dO.
(6)
In general, the equation for the nth moment of the intensity is obtained by taking the nth moment of Eq. (1) and will always involve the divergence of the (n + 1)th moment of the intensity. If these moment equations are to be useful in simplifying the study of radiative transport problems, then some approximations are needed to truncate this infinite sequence of coupled equations. Usually, this means justifying the use of some expression for a given moment in terms of lower ones. The problem of finding such an approximation that is also physically consistent is called the moment closure problem. The P-1 equations are usually written in terms of the Eddington tensor, T(r, t), defined by P=E~"
(7)
and consist of Eq. (2) and Eq. (5) rewritten as (1/c) t,F + cV • ( ~ E ) + oTV = O,
(8)
where aT = a -- g~as is the transport coefficient and T is considered to be a known function of E, F, and other variables in the problem. The classical Eddington approximation 1~ is T=(1/3)I
(9)
as for an isotropic intensity. We shall see in the next section that in the choice of an expression for the Eddington tensor, we must consider certain moment conditons; in Section 3; specific expressions are given that. depend on E and F in a nontrivial way. The diffusion equation is Eq. (2) with the flux, F, considered to be a known function of E, F E , and other variables in the problem. The classical Eddington approximation 11 is F = - (c/3ar)VE
(10)
as for an isotropic intensity in steady state. It is clear from Eq. (4) that
IF[ _< cE,
(1 l)
which is not generally satisfied by Eq. (10). A diffusion equation is called flux-limited if
Relating Eddington factors to flux limiters
151
the expression for the flux used satisfies Eq. (11) identically. In Section 4, we derive flux-limited diffusion theories from the P-1 equations.
2. MOMENT CONDITIONS The information contained in various moments is not independent, but is constrained by the fact that the moments are derived from a nonnegative density on the unit sphere and as such must satisfy equalities and inequalities. To see this for the lower moments, it is useful to introduce the new dimensionless quantities ~(I2, r, t), and f(r, t) defined by
I =cE~,
F=cEf;
is called the normalized intensity and f is called the normalized flux. Equations (3), (4), and (6) may be rewritten in terms of these new quantities as
1 = f4. 4~(12, r, t) d12, f(r, t) = f4,12~(~2, r, t) dO,
T(r, t) =
f4.
I212 ~(I2, r, t) dO,
(12a)
(12b)
(12c)
Hence, f and "r are the first and second moments of the nonnegative unit density ¢i on the unit sphere. Using the nonnegativity of 4, the fact that If21 = 1, and Eq. (12a), we see from Eq. (12b) t h a t f m u s t satisfy Ill _< 1,
(13)
This is just a restatement of the flux-limiting condition, Eq. (11). In order to state constraints on T, we need two concepts from linear algebra. The first is the result that the trace of a given matrix lfl, which is denoted tr(Iffl), is equal to the sum of the diagonal elements of l('I and is alsoequal to the sum of the eigenvalues of 1VI. The second is that a given symmetric matrix M is called nonnegative when x • Mx >__0 for an arbitrary vector x. This property is denoted by M > 0 and is equivalent to saying that all of the eigenvalues of 1('I are nonnegative or that the determinant of each leading principal minor is nonnegative. ~2 Since f and T are the first and second moments of a nonnegative unit density on the unit sphere, they must satisfy the constraints
tr ('F) = l,
(14)
1"-ff>_ 0.
(15)
The first is obtained directly from eqs. (12c) [using the fact tr(1212) = 12 • I2 = 1] and (12a) while the second is a consequence of the identity x. ('r - ff)x -- f4~ [I2 .x - f. x]2q~(~) dr2, which follows from eq. (12) for any vector x and the nonnegativity of 4. In other words eq. (14) is a reflection of the fact the density is defined on the unit sphere and Eq. (15) reflects the nonnegativity of the density. Taken together, Eqs. (14) and (15) imply Eq. (13). The significance of the constraints of Eqs. (14) and (15) to the closure problem lies in
152
C.D. LEVERMORE
the fact that, if a vector f and a symmetric matrix T satisfy Eqs. (14) and (15), then they may be realized as the first and second moments of a nonnegative unit density on the sphere. The density q~(~) will be uniquely determined only when Ifl-- l, in which case = ff and ~b(g2) is the Dirac weight function at the point f on the unit sphere. For a proof of these statements and a discussion of more general ideas, see Ref. 4. The form of "F and the resulting conditions, Eqs. (14) and (15), will be simplified somewhat when it is known or assumed that the intensity is symmetric about a preferred direction. This will be the case, for example, if the problem being considered has a spherical or slab symmetry. Let the preferred direction be denoted by the unit vector n(r, t). Since cb(~) must be invariant under rotations that fix n, it follows that q~ is a function only of n. ~ and that, by Eqs. (4) and (6), the vector f and the matrix "F must also be invariant under such rotations. The only way a vector can be invariant under such rotations is if it is proportional to n. It follows immediately that f has the form
f=fn, w h e r e f = If]. Furthermore, the vector 1"n must also be invariant under such rotations and so n must be an eigenvector of T with some eigenvalue X, "Fn --~ ~(n.
By symmetry, the plane perpendicular to n must then be an eigenspace of T with the eigenvalue determined by the trace identity, Eq. (14), to be (1 - Z ) / 2 . Thus T may be written in the form = i(i
- ~()/2 + n n ( 3 ~
-
1)/2
(16)
where I is the identity matrix. The inequality constraint, Eq. (15), then reduces to the simple requirement fa_
(17)
The quantity x(r, t), which is known as the Eddington factor, can be expressed directly in terms of the normalized intensity by the formula x(r, t) = f4~ [fa • n(r, t)] 2 q~(fa, r, t) d~.
(18)
3. EDDINGTON FACTORS In choosing an expression for the Eddington tensor to close the P-1 equations there are three points of view that may be adopted regarding the condition expressed in Eq. (15). The first would be to ignore the condition entirely and assume it is not significant. The second would be to ignore it while making the choice but to use it as a check on the consistency of the solution. The third would be to pick an expression which satisfies the condition identically. Application of the last point of view will be discussed below. There are two basic approaches to picking T in such a way that Eq. (15) is satisfied. One is to choose an expression for T without any reference to the underlying density q~ such as Kershaw 4 suggested, z(f) = (1 + 2Ja)/3.
(19)
This has the advantage of being simple and the disadvantage of being entirely ad hoc, with no underlying conception of the nature of the density it is trying to describe. The other approach attempts to overcome that last criticism by assuming a specific form for the normalized density ~b(f2) as a function of f and then use Eq. (18) to compute x ( f ) . As long
Relating Eddington factors to flux limiters
153
as the assumed form for q~ is a nonnegative unit density, the resulting expression for Z will satisfy Eq. (17). This approach is embodied in the following four examples which are recorded here for future reference. Making the simplifying assumption that photons obey Maxwell-Boltzmann statistics, Minerbo 3 used entropy arguments to obtain
~ (12 ) = (1/ 4n )( Z /sinh Z ) e z" ' ~, where Z is a parameter related to f by f = coth Z - (l/Z). A s f r a n g e s from 0 to 1, Z will range from 0 to ~ . The resulting Eddington factor is Z = 1 - (2/Z)[coth Z - (l/Z)]
(20)
and was called the maximum entropy Eddington factor. Any intensity may be approximated as a linear function o f n • Q, so Minerbo 3 suggested one could take
f(1/4n)(1 + 3fn. f~) ~(12) = [max [0, (1/9rr)(n. 1"2 - 3 f + 2)/(1 _f)2]
if0 < f < (1/3) if(l/3) < f < 1.
Then the Eddington factor is
)(1/3) z(f) = ~(1/2)(1 _ f ) 2 + f 2
if0 < f < (1/3) if(l/3) < f < 1.
(21)
This choice of q~ is piecewise linear in n. f2 and reproduces the classical Eddington approximation for f < (1/3). Using a Chapman-Enskog approach, Levermore 7'8 derived the intensity O(O) = (1/47r)(R coth R - R n . 1'2)-% where R is a parameter related to f by f = coth R - (I/R).
(22)
The corresponding Eddington factor is Z = coth R [coth R - (I/R)].
(23)
This choice of • has the form of the angular density of a normal mode for the transport equation with isotropic scattering, l° The last example is applied when dealing with frequency integrated quantities and assumes that the density is isotropic in some inertial frame. Using the contravariance of the stress-energy tensor 9 we then find that X = (1 + 3//2)/(3 + f12),
(24a)
f = 4/~/(3 +/32),
(24b)
where tic is the speed of the Lorentz frame in which the density is isotropic. As f ranges from 0 to 1, fl will also range from 0 to 1. These equations are simple enough that one can eliminate fl and solve for Z explicitly as a function o f f :
154
C. D. LEVERMORE
Z : (1/3) + 2fl/(2 + x/4 - 3f2) = (3 + 4f2)/(5 + 2x//-4- 3f21.
(25)
In Fig. 1, these five Eddington factors are plotted for comparison. Note all have the value of 1/3 at f = 0 and as functions of f a r e increasing and concave up. In all of the above the tensor T depends algebraicly on E and F in a nontrivial way, so if they are used to close the moment equations (12) one obtains a nonlinear hyperbolic system. Since solutions may then develop jump discontinuities in finite times, 13the resulting system may not be an acceptable description of transport. Thus other settings will be examined in which these Eddington factors may be employed.
4. FLUX LIMITED DIFFUSION One use of the moment equations is to derive a diffusion equation for the energy density E. The usual procedure is to first drop the 0re term in Eq. (8) to obtain F = - (c/ar)V" ('rE)
(26)
and then substitute the result into Eq. (12a). This will give a closed equation for E provided T is independent o f f (and therefore F). But ifT is chosen to be consistent with Eq. (15) then it must depend on f. The problem is usually gotten around by assuming the spatial variations of T can be neglected and using Eq. (26) to arrive at an algebraic relation between f, E, and VE:
TX = f,
(27)
where X(r, t) is the dimensionless gradient defined by X = - (VE)/(arE).
The idea is to now solve Eq. (21) for f in terms of X and then substitute the result into Eq. (2). However, assume for the moment that T is of the form given by eq. (16) and X is of the form X = Xn; then X and f are related by f = zOOX. Ifm is the smallest value that g(f) takes on, then the largest value X m a y take on consistent with Eq. (13) is less 1/m. Thus the resulting expression for the flux will violate Eq. (13) when X gets too large and the diffusion theory is not flux-limited. One might think this state of affairs arose due to the approximation leading to Eq. (27), but on closer examination one finds the step leading to Eq. (26) is the cause. This is seen by considering the case of extreme streaming, when an unidirectional front of radiation is penetrating a cold absorbing media (B ~ E, as ~ aa). Mathematically this means F,.~cEn,
P ,,~ E nn,
where n is the unit vector in the direction of the flux. Comparing the corresponding terms in Eqs. (2) and (8) one finds (1/c)Ot F ~ O,En, cV • P ~
c(n. V E ) n ,~, ( V • F)n,
(a -- gla~)F ,-~ c a a E n ,~, c a a ( E - B)n. So assuming OtF can be dropped in Eq. (8) means OrE must be dropped in Eq. (2) if one is to be consistent and the resulting equations reduce to the single equation
0.4
0.6
t
0.8
l.Or '
0.2
[MINERBO (LINEAR) ] - - - - ~ ]
[MINERBO (STATISTICAL)~---
ILEVERMORE(LORENTZ)]
[KERSHAW
04
[LEVERMORE(CHAPMAN-ENSKOG) I
,
f
Fig. 1.
o16
'
0.8
'
i,o
R
X
=b
5~
o
o
m
rn
c. D. LEVERMORE
156
n . V E + (a -- g , a , ) E ~ 0
or equivalently X --- 1. This result can be seen directly from Eq. (27) by setting f = 1 and using Eq. (15). Thus such a theory can only handle steady-state extreme streaming and could not be flux-limited. Evidently a new approximation is needed if one hopes to derive a flux-limited diffusion theory. Before proceeding further, the Eqs. (2) and (8) will be rewritten in terms of the dimensionless quantities that appear naturally in Eqs. (14) and (15). They become OrE + c V • (fE)
+ c a A ( E -- B ) = O,
c3,(fE) + c V • ( T E ) + c a r f E = O.
(28a) (28b)
Using Eq. (22a) to eliminate 0rE from Eq. (22b) one obtains [(1/c)Orf + f" V f ] E + V " [('r
- if)E] + COarfE = O,
(29)
where e) is the effective albedo O) =
[O's(1 --
gl)E + aABI/(arE).
The key assumption of the following diffusion theory is motivated by the observation that in both the isotropic and extreme streaming cases the spatial and temporal variations of f and "r are small. Neglecting terms involving derivations in Eq. (29) gives an algebraic relation between T, f, E and IZE; (T - f f ) R = f,
(30)
where R(r, t) is the dimensionless gradient defined by R = - (1/~ar)(VE/E).
The idea is to choose "r as a function of f consistent with Eq. (15) and solve Eq. (30) for f in terms of R and substitute the result into Eq. (28a) to obtain a diffusion equation. This task is simplified by assuming T has the form given by Eq. (16), in which case R and f will be proportional and R = IRI is related to f by R =f/[z(f)
-Jc].
(31)
Assuming the righthand side is an increasing function o f f , one may then solve for f as a function of R. This mild assumption on Z is satisfied by all the expressions considered in Section 3. If in addition, Z satisfies the constraint of Eq. (17) and Z(0)= 0, then a s f increases from 0 to 1, R will increase from 0 to oo. When f is considered as a function of R, then as R increases from 0 to oo,fwill increase from 0 to 1 and the resulting diffusion theory will be flux-limited. Defining 2(R) by f =2(R)R
gives the flux F = - (e/ooar)2(R)VE;
(32)
2(R) is called the flux-limiter associated with the Eddington factor Z(f)- The resulting diffusion equation is a nonlinear parabolic equation provided co :/: 0 so solutions will not form discontinuities. The functions 2 and Z are related implicitly by
Relating Eddington factors to flux limiters
~(R) = zoo - f
157
(33)
or, equivalently, by z(f) = 2(R) + 2 ( R ) 2 R 2 so that physically desirable properties of one may be related to analytic properties of the other. For example: (C1) Since Z must satisfy Eq. (17), then 2 must satisfy 2(R) + 2(R)2R 2 < 1. (C2) If when f = 0 there is no preferred direction then the underlying intensity must be assumed isotropic and X(0)= 1/3. Therefore 2 ( 0 ) = 1/3. The resulting diffusion theory will approach isotropic diffusion theory in the limit of R small and ¢0 near 1 as it should. (C3) It is often required that zOO be an increasing function of f since by Eq. (18) the Eddington factor is seen to be a rough measure of the anistropy of ~(I2) in the direction of f. This will be the case if and only if ).(R) + 2 ( R ) Z R 2 is increasing as a function of R. In particular this says that a flux-limiter should not decrease from its value at R = 0 faster than quadratically in R, which is also required to be consistent with isotropic diffusion theory. (C4) Since flux-limiting should be an effect that increases as gradients get bigger and one moves away from the classical diffusion limit, it is often assumed 2(R) is a decreasing function of R. This is equivalent to assuming ZOO _ f 2 is decreasing as a function of f, which immediately implies that the righthandside of Eq. (31) is increasing and so is stronger than that previous assumption.
5. A P P L I C A T I O N S
If these relations are applied to the Eddington factors mentioned in Section 3, all of which satisfy the above conditions, then the related flux-limiters are found to be the following: for Eq. (19), 2(R) = 2/(3 + x/9 + 4R2);
(34)
for Eq. (20), ,~, = ( 1 / Z 2) - -
csch2Z,
where the parameter Z is related to R by R = [coth Z - ( 1 / Z ) ] / [ ( 1 / Z 2) - csch 2 Z]; for Eq. (21),
(35)
158
C. D. LEVERMORE
)'2/(3 + 9 x / ~ 12R2)_ 2(R) = [1/(1 + R + ~/I + 2R)
i f 0 < R _<(3/2), i f ( 3 / 2 ) < R < o c;
(36)
for Eq. (23), 2(R) = (1/R)[coth R -- (I/R)];
(37)
,~ = 3(1 -/~2)2/(3 +/~2)2,
(38)
for Eq. (24),
where
R = 4/3(3 +/32)/(1 -/~2)2.
These five flux-limiters are plotted for comparison in Fig. 2. It should be emphasized that the above theory gives a unified framework for relating flux-limiters and Eddington factors and does not single out a particular flux-limiter or Eddington factor as "best" for any problem. In practice one is often confronted with an ad hoc flux limiter, about which physical insight might be gained by knowing the related Eddington factor. For example, we consider the sum flux-limiter proposed by Wilson,5 viz. 2(R) = 1/(3 + R).
(39)
The related Eddington factor is z(f)=(1
-f+3fz)/3,
which decreases from Z = 1/3 a t f = 0 to a minimum of X = 11/36 a t f = 1/6 and so violates condition (C3). This indicates the resulting expression for the flux probability under estimates the actual flux for gradients that are not too big (R < 3). On the other hand, the cut-off flux limiter that is essentially found in Winslow6 is 2"~"
fl/3
(a)=~I/R
if0
(40)
and has as its related Eddington factor
z(f) = (1/3) +f2, which violates condition (C1) i f f > x f ~ . This indicates the resulting expression for the flux probability over-estimates the actual flux for moderately sized gradients (2 < R < 10). It also clearly over-estimates the flux for larger gradients but the error is not serious since one presumes f is actually near 1 in that case. The flux-limiter proposed by Levermore] '8 is given by Eq. (37) and therefore suffers none of the difficulties of the previous two examples. 6. CONCLUDING REMARKS In order to complete the derivation of diffusion theory from the P-1 equations, a word should be said about boundary conditions. For the P-1 equations, some scalar combination of E and F is usually specified at each boundary point. To obtain the boundary conditions for diffusion equation merely use Eq. (32) to replace F with an expression involving E and IZEin the P-1 boundary conditions. The derived relation between Eddington factors and flux limiters, Eq. (33), was a direct
Relating Eddington factors to flux limiters
I
!
I
159
o
2
h C F
h n.
0
iv
Q
h
tl
C
"
h
oZ
--
~q L~
Q
I
I
N
I
0
--
C)
160
C . D . LEVERMORE
consequence of the new derivation of the diffusion equation presented in Section 4. A different result would be obtained if one could find a different derivation yielding a different flux-limited diffusion equation from the P-1 equations. However, even if another relation could be found, it could not be simpler than the present result and the two results would have to agree within the order of approximation over their common domain of validity. Equation (30) can be thought of as a generalization of Fick's law which allows for the possibility that the flux, F, is not in the same direction as - VE. However, it follows from the condition given in Eq. (15) that the angle between the two vectors cannot be greater than 90 ~'. This is a reflection of the fact the theory assumes some redistribution of the radiation by the matter. In the limit of a cold absorber (~o-~0) the theory becomes singular and like any other diffusion theory, can be expected to perform poorly at times. Perhaps the most striking fact about the flux-limited diffusion equations is that one can employ flux limiters that arise from Eddington factors given by single algebraic functions and obtain a well behaved theory. Yet these same Eddington factors can lead to anomalous results (shocks) when used in the P-1 equations. More work is needed to make our understanding of this phenomenon more precise.
1. 2. 3. 4. 5. 6. 7, 8. 9. 10. 11. 12. 13.
REFERENCES G. C. Pomraning, JQSRT 9, 407 (1969). G. C. Pomraning, Nukleonik 6, 348 (1965). G. N. Minerbo, JQSRT 20, 541 (1978). D. S. Kershaw, "Flux Limiting Nature's Own Way", Lawrence Livermore National Laboratory, Livermore, CA, UCRL-78378 (1976). J. M. LeBlanc and J. R. Wilson, Appl. J. 161, 541 (1970). A. M. Winslow, Nucl. Sci. Engng 32, 101 (1968). C. D. Levermore, "A Chapman-Enskog Approach to Flux-Limited Diffusion Theory", Lawrence Livermore National Laboratory, Livermore, CA, UCID-18229 (1979). C. D. Levermore and G. C. Pomraning, Appl. J. 248, 321 (1981). G. C. Pomraning, The Equations of Radiation Hydrodynamics. Pergamon Press, Oxford (1973). K. M. Case and P. F. Zweifel, Linear Transport Theory. Addison Wesley, Reading, Mass. (1967). A. Eddington, The Internal Constitution of Stars, Dover, New York (1926). K. Hoffman and R. Kunze, Linear Algebra, p. 328. Prentice-Hall, Englewood Cliffs, New Jersey (1961). R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. 2, p. 488. Wiley, New York (1962).
0022-4073/84 $3.00+ .00 Pergamon Press Ltd.
RELATING EDDINGTON FACTORS TO FLUX LIMITERSt C. D. LEVERMORE A-Division, Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, CA 94550, U.S.A.
(Reveived 6 May 1983) Abstract--Variable Eddington factors and flux-limiters have been introduced in the P-1 and diffusion equations, respectively, to handle sittiations when the underlying intensity is too anisotropic for the unmodified theories to remain valid. We present a derivation of a relation between the two for which a new approach to the diffusion approximation is used, Algebraic expressions for Eddington factors satisfying the moment conditions are not satisfactory for closing the P-1 equations but, by using the derived relation, yield acceptable flux-limited diffusion theories.
1. I N T R O D U C T I O N
In many problems involving transport theory, resource constraints may necessitate adoption of some approximate method to calculate the transport phenomena. Two such methods are the P-1 and diffusion approximations, both of which require a degree of isotropy in the underlying intensity. In order to handle situations with significant anisotropy, variable Eddington factors have been introduced for the P-1 approximation 1-4 and flux-limiters for the diffusion approximation. 5-8While it has been understood that these two ideas are related, no explicit relationship has been developed until now. Assuming coherent scattering and stationary matter, the radiation transport equation for the specific intensity of photons, I(1~, r, t) is 9'1° (1/c)8tl + ~ • g I + trI = (c/4x)trAB + as f g(I2 • ~ ' ) I(O') dO'. d4~z
(1)
Here t2, r, and t are the angular, spatial and temporal variables, c is the speed of light, aA(r, t) is the absorption coefficient, trs(r, t) is the scattering coefficient, a = aA + as is the total interaction coefficient, B(r, t) is the black body energy density, and g(f2 • l~', r, t) is the scattering angular redistribution function which is normalized by 1=
g(~
• ~ ' ) dQ = 2x
fl
g ( ~ ) d/~.
These equations may be considered to be either frequency dependent with the frequency parameters suppressed or gray equations. Integrating Eq. (1) over the solid angle dQ yields 8,E + V " F + ctrA ( E -
B ) = O,
(2)
where the energy density, E(r, t), and the energy flux, F(r, t) are defined in terms of the zeroth and first moments of the intensity by E(r, t) = (l/c) f4~ I(1"2, r, t) dO
(3)
tThis work was performed under the auspices of the U.S. Department of Energy at the Lawrence Livermore National Laboratory under contract number W-7405-ENG-48. 149
150
C.D. LEVERMORE
and F(r, t) = I OI(I'L r, t) dO. d4=
(4)
An equation for the flux is obtained by multiplying Eq. (1) by fl and then integrating over the solid angle dQ, giving (1/c)dtF + c F • P + (a - glas)F = 0,
(5)
where g~(r, t) is the first moment of the angular redistribution function gl = 2re
f
l
#g(/t) d# 1
and the radiation pressure tensor, P(r, t), is defined in terms of the second moment of the intensity by
P(r, t) = (l/c) f4, 1~12I(12, r, t) dO.
(6)
In general, the equation for the nth moment of the intensity is obtained by taking the nth moment of Eq. (1) and will always involve the divergence of the (n + 1)th moment of the intensity. If these moment equations are to be useful in simplifying the study of radiative transport problems, then some approximations are needed to truncate this infinite sequence of coupled equations. Usually, this means justifying the use of some expression for a given moment in terms of lower ones. The problem of finding such an approximation that is also physically consistent is called the moment closure problem. The P-1 equations are usually written in terms of the Eddington tensor, T(r, t), defined by P=E~"
(7)
and consist of Eq. (2) and Eq. (5) rewritten as (1/c) t,F + cV • ( ~ E ) + oTV = O,
(8)
where aT = a -- g~as is the transport coefficient and T is considered to be a known function of E, F, and other variables in the problem. The classical Eddington approximation 1~ is T=(1/3)I
(9)
as for an isotropic intensity. We shall see in the next section that in the choice of an expression for the Eddington tensor, we must consider certain moment conditons; in Section 3; specific expressions are given that. depend on E and F in a nontrivial way. The diffusion equation is Eq. (2) with the flux, F, considered to be a known function of E, F E , and other variables in the problem. The classical Eddington approximation 11 is F = - (c/3ar)VE
(10)
as for an isotropic intensity in steady state. It is clear from Eq. (4) that
IF[ _< cE,
(1 l)
which is not generally satisfied by Eq. (10). A diffusion equation is called flux-limited if
Relating Eddington factors to flux limiters
151
the expression for the flux used satisfies Eq. (11) identically. In Section 4, we derive flux-limited diffusion theories from the P-1 equations.
2. MOMENT CONDITIONS The information contained in various moments is not independent, but is constrained by the fact that the moments are derived from a nonnegative density on the unit sphere and as such must satisfy equalities and inequalities. To see this for the lower moments, it is useful to introduce the new dimensionless quantities ~(I2, r, t), and f(r, t) defined by
I =cE~,
F=cEf;
is called the normalized intensity and f is called the normalized flux. Equations (3), (4), and (6) may be rewritten in terms of these new quantities as
1 = f4. 4~(12, r, t) d12, f(r, t) = f4,12~(~2, r, t) dO,
T(r, t) =
f4.
I212 ~(I2, r, t) dO,
(12a)
(12b)
(12c)
Hence, f and "r are the first and second moments of the nonnegative unit density ¢i on the unit sphere. Using the nonnegativity of 4, the fact that If21 = 1, and Eq. (12a), we see from Eq. (12b) t h a t f m u s t satisfy Ill _< 1,
(13)
This is just a restatement of the flux-limiting condition, Eq. (11). In order to state constraints on T, we need two concepts from linear algebra. The first is the result that the trace of a given matrix lfl, which is denoted tr(Iffl), is equal to the sum of the diagonal elements of l('I and is alsoequal to the sum of the eigenvalues of 1VI. The second is that a given symmetric matrix M is called nonnegative when x • Mx >__0 for an arbitrary vector x. This property is denoted by M > 0 and is equivalent to saying that all of the eigenvalues of 1('I are nonnegative or that the determinant of each leading principal minor is nonnegative. ~2 Since f and T are the first and second moments of a nonnegative unit density on the unit sphere, they must satisfy the constraints
tr ('F) = l,
(14)
1"-ff>_ 0.
(15)
The first is obtained directly from eqs. (12c) [using the fact tr(1212) = 12 • I2 = 1] and (12a) while the second is a consequence of the identity x. ('r - ff)x -- f4~ [I2 .x - f. x]2q~(~) dr2, which follows from eq. (12) for any vector x and the nonnegativity of 4. In other words eq. (14) is a reflection of the fact the density is defined on the unit sphere and Eq. (15) reflects the nonnegativity of the density. Taken together, Eqs. (14) and (15) imply Eq. (13). The significance of the constraints of Eqs. (14) and (15) to the closure problem lies in
152
C.D. LEVERMORE
the fact that, if a vector f and a symmetric matrix T satisfy Eqs. (14) and (15), then they may be realized as the first and second moments of a nonnegative unit density on the sphere. The density q~(~) will be uniquely determined only when Ifl-- l, in which case = ff and ~b(g2) is the Dirac weight function at the point f on the unit sphere. For a proof of these statements and a discussion of more general ideas, see Ref. 4. The form of "F and the resulting conditions, Eqs. (14) and (15), will be simplified somewhat when it is known or assumed that the intensity is symmetric about a preferred direction. This will be the case, for example, if the problem being considered has a spherical or slab symmetry. Let the preferred direction be denoted by the unit vector n(r, t). Since cb(~) must be invariant under rotations that fix n, it follows that q~ is a function only of n. ~ and that, by Eqs. (4) and (6), the vector f and the matrix "F must also be invariant under such rotations. The only way a vector can be invariant under such rotations is if it is proportional to n. It follows immediately that f has the form
f=fn, w h e r e f = If]. Furthermore, the vector 1"n must also be invariant under such rotations and so n must be an eigenvector of T with some eigenvalue X, "Fn --~ ~(n.
By symmetry, the plane perpendicular to n must then be an eigenspace of T with the eigenvalue determined by the trace identity, Eq. (14), to be (1 - Z ) / 2 . Thus T may be written in the form = i(i
- ~()/2 + n n ( 3 ~
-
1)/2
(16)
where I is the identity matrix. The inequality constraint, Eq. (15), then reduces to the simple requirement fa_
(17)
The quantity x(r, t), which is known as the Eddington factor, can be expressed directly in terms of the normalized intensity by the formula x(r, t) = f4~ [fa • n(r, t)] 2 q~(fa, r, t) d~.
(18)
3. EDDINGTON FACTORS In choosing an expression for the Eddington tensor to close the P-1 equations there are three points of view that may be adopted regarding the condition expressed in Eq. (15). The first would be to ignore the condition entirely and assume it is not significant. The second would be to ignore it while making the choice but to use it as a check on the consistency of the solution. The third would be to pick an expression which satisfies the condition identically. Application of the last point of view will be discussed below. There are two basic approaches to picking T in such a way that Eq. (15) is satisfied. One is to choose an expression for T without any reference to the underlying density q~ such as Kershaw 4 suggested, z(f) = (1 + 2Ja)/3.
(19)
This has the advantage of being simple and the disadvantage of being entirely ad hoc, with no underlying conception of the nature of the density it is trying to describe. The other approach attempts to overcome that last criticism by assuming a specific form for the normalized density ~b(f2) as a function of f and then use Eq. (18) to compute x ( f ) . As long
Relating Eddington factors to flux limiters
153
as the assumed form for q~ is a nonnegative unit density, the resulting expression for Z will satisfy Eq. (17). This approach is embodied in the following four examples which are recorded here for future reference. Making the simplifying assumption that photons obey Maxwell-Boltzmann statistics, Minerbo 3 used entropy arguments to obtain
~ (12 ) = (1/ 4n )( Z /sinh Z ) e z" ' ~, where Z is a parameter related to f by f = coth Z - (l/Z). A s f r a n g e s from 0 to 1, Z will range from 0 to ~ . The resulting Eddington factor is Z = 1 - (2/Z)[coth Z - (l/Z)]
(20)
and was called the maximum entropy Eddington factor. Any intensity may be approximated as a linear function o f n • Q, so Minerbo 3 suggested one could take
f(1/4n)(1 + 3fn. f~) ~(12) = [max [0, (1/9rr)(n. 1"2 - 3 f + 2)/(1 _f)2]
if0 < f < (1/3) if(l/3) < f < 1.
Then the Eddington factor is
)(1/3) z(f) = ~(1/2)(1 _ f ) 2 + f 2
if0 < f < (1/3) if(l/3) < f < 1.
(21)
This choice of q~ is piecewise linear in n. f2 and reproduces the classical Eddington approximation for f < (1/3). Using a Chapman-Enskog approach, Levermore 7'8 derived the intensity O(O) = (1/47r)(R coth R - R n . 1'2)-% where R is a parameter related to f by f = coth R - (I/R).
(22)
The corresponding Eddington factor is Z = coth R [coth R - (I/R)].
(23)
This choice of • has the form of the angular density of a normal mode for the transport equation with isotropic scattering, l° The last example is applied when dealing with frequency integrated quantities and assumes that the density is isotropic in some inertial frame. Using the contravariance of the stress-energy tensor 9 we then find that X = (1 + 3//2)/(3 + f12),
(24a)
f = 4/~/(3 +/32),
(24b)
where tic is the speed of the Lorentz frame in which the density is isotropic. As f ranges from 0 to 1, fl will also range from 0 to 1. These equations are simple enough that one can eliminate fl and solve for Z explicitly as a function o f f :
154
C. D. LEVERMORE
Z : (1/3) + 2fl/(2 + x/4 - 3f2) = (3 + 4f2)/(5 + 2x//-4- 3f21.
(25)
In Fig. 1, these five Eddington factors are plotted for comparison. Note all have the value of 1/3 at f = 0 and as functions of f a r e increasing and concave up. In all of the above the tensor T depends algebraicly on E and F in a nontrivial way, so if they are used to close the moment equations (12) one obtains a nonlinear hyperbolic system. Since solutions may then develop jump discontinuities in finite times, 13the resulting system may not be an acceptable description of transport. Thus other settings will be examined in which these Eddington factors may be employed.
4. FLUX LIMITED DIFFUSION One use of the moment equations is to derive a diffusion equation for the energy density E. The usual procedure is to first drop the 0re term in Eq. (8) to obtain F = - (c/ar)V" ('rE)
(26)
and then substitute the result into Eq. (12a). This will give a closed equation for E provided T is independent o f f (and therefore F). But ifT is chosen to be consistent with Eq. (15) then it must depend on f. The problem is usually gotten around by assuming the spatial variations of T can be neglected and using Eq. (26) to arrive at an algebraic relation between f, E, and VE:
TX = f,
(27)
where X(r, t) is the dimensionless gradient defined by X = - (VE)/(arE).
The idea is to now solve Eq. (21) for f in terms of X and then substitute the result into Eq. (2). However, assume for the moment that T is of the form given by eq. (16) and X is of the form X = Xn; then X and f are related by f = zOOX. Ifm is the smallest value that g(f) takes on, then the largest value X m a y take on consistent with Eq. (13) is less 1/m. Thus the resulting expression for the flux will violate Eq. (13) when X gets too large and the diffusion theory is not flux-limited. One might think this state of affairs arose due to the approximation leading to Eq. (27), but on closer examination one finds the step leading to Eq. (26) is the cause. This is seen by considering the case of extreme streaming, when an unidirectional front of radiation is penetrating a cold absorbing media (B ~ E, as ~ aa). Mathematically this means F,.~cEn,
P ,,~ E nn,
where n is the unit vector in the direction of the flux. Comparing the corresponding terms in Eqs. (2) and (8) one finds (1/c)Ot F ~ O,En, cV • P ~
c(n. V E ) n ,~, ( V • F)n,
(a -- gla~)F ,-~ c a a E n ,~, c a a ( E - B)n. So assuming OtF can be dropped in Eq. (8) means OrE must be dropped in Eq. (2) if one is to be consistent and the resulting equations reduce to the single equation
0.4
0.6
t
0.8
l.Or '
0.2
[MINERBO (LINEAR) ] - - - - ~ ]
[MINERBO (STATISTICAL)~---
ILEVERMORE(LORENTZ)]
[KERSHAW
04
[LEVERMORE(CHAPMAN-ENSKOG) I
,
f
Fig. 1.
o16
'
0.8
'
i,o
R
X
=b
5~
o
o
m
rn
c. D. LEVERMORE
156
n . V E + (a -- g , a , ) E ~ 0
or equivalently X --- 1. This result can be seen directly from Eq. (27) by setting f = 1 and using Eq. (15). Thus such a theory can only handle steady-state extreme streaming and could not be flux-limited. Evidently a new approximation is needed if one hopes to derive a flux-limited diffusion theory. Before proceeding further, the Eqs. (2) and (8) will be rewritten in terms of the dimensionless quantities that appear naturally in Eqs. (14) and (15). They become OrE + c V • (fE)
+ c a A ( E -- B ) = O,
c3,(fE) + c V • ( T E ) + c a r f E = O.
(28a) (28b)
Using Eq. (22a) to eliminate 0rE from Eq. (22b) one obtains [(1/c)Orf + f" V f ] E + V " [('r
- if)E] + COarfE = O,
(29)
where e) is the effective albedo O) =
[O's(1 --
gl)E + aABI/(arE).
The key assumption of the following diffusion theory is motivated by the observation that in both the isotropic and extreme streaming cases the spatial and temporal variations of f and "r are small. Neglecting terms involving derivations in Eq. (29) gives an algebraic relation between T, f, E and IZE; (T - f f ) R = f,
(30)
where R(r, t) is the dimensionless gradient defined by R = - (1/~ar)(VE/E).
The idea is to choose "r as a function of f consistent with Eq. (15) and solve Eq. (30) for f in terms of R and substitute the result into Eq. (28a) to obtain a diffusion equation. This task is simplified by assuming T has the form given by Eq. (16), in which case R and f will be proportional and R = IRI is related to f by R =f/[z(f)
-Jc].
(31)
Assuming the righthand side is an increasing function o f f , one may then solve for f as a function of R. This mild assumption on Z is satisfied by all the expressions considered in Section 3. If in addition, Z satisfies the constraint of Eq. (17) and Z(0)= 0, then a s f increases from 0 to 1, R will increase from 0 to oo. When f is considered as a function of R, then as R increases from 0 to oo,fwill increase from 0 to 1 and the resulting diffusion theory will be flux-limited. Defining 2(R) by f =2(R)R
gives the flux F = - (e/ooar)2(R)VE;
(32)
2(R) is called the flux-limiter associated with the Eddington factor Z(f)- The resulting diffusion equation is a nonlinear parabolic equation provided co :/: 0 so solutions will not form discontinuities. The functions 2 and Z are related implicitly by
Relating Eddington factors to flux limiters
~(R) = zoo - f
157
(33)
or, equivalently, by z(f) = 2(R) + 2 ( R ) 2 R 2 so that physically desirable properties of one may be related to analytic properties of the other. For example: (C1) Since Z must satisfy Eq. (17), then 2 must satisfy 2(R) + 2(R)2R 2 < 1. (C2) If when f = 0 there is no preferred direction then the underlying intensity must be assumed isotropic and X(0)= 1/3. Therefore 2 ( 0 ) = 1/3. The resulting diffusion theory will approach isotropic diffusion theory in the limit of R small and ¢0 near 1 as it should. (C3) It is often required that zOO be an increasing function of f since by Eq. (18) the Eddington factor is seen to be a rough measure of the anistropy of ~(I2) in the direction of f. This will be the case if and only if ).(R) + 2 ( R ) Z R 2 is increasing as a function of R. In particular this says that a flux-limiter should not decrease from its value at R = 0 faster than quadratically in R, which is also required to be consistent with isotropic diffusion theory. (C4) Since flux-limiting should be an effect that increases as gradients get bigger and one moves away from the classical diffusion limit, it is often assumed 2(R) is a decreasing function of R. This is equivalent to assuming ZOO _ f 2 is decreasing as a function of f, which immediately implies that the righthandside of Eq. (31) is increasing and so is stronger than that previous assumption.
5. A P P L I C A T I O N S
If these relations are applied to the Eddington factors mentioned in Section 3, all of which satisfy the above conditions, then the related flux-limiters are found to be the following: for Eq. (19), 2(R) = 2/(3 + x/9 + 4R2);
(34)
for Eq. (20), ,~, = ( 1 / Z 2) - -
csch2Z,
where the parameter Z is related to R by R = [coth Z - ( 1 / Z ) ] / [ ( 1 / Z 2) - csch 2 Z]; for Eq. (21),
(35)
158
C. D. LEVERMORE
)'2/(3 + 9 x / ~ 12R2)_ 2(R) = [1/(1 + R + ~/I + 2R)
i f 0 < R _<(3/2), i f ( 3 / 2 ) < R < o c;
(36)
for Eq. (23), 2(R) = (1/R)[coth R -- (I/R)];
(37)
,~ = 3(1 -/~2)2/(3 +/~2)2,
(38)
for Eq. (24),
where
R = 4/3(3 +/32)/(1 -/~2)2.
These five flux-limiters are plotted for comparison in Fig. 2. It should be emphasized that the above theory gives a unified framework for relating flux-limiters and Eddington factors and does not single out a particular flux-limiter or Eddington factor as "best" for any problem. In practice one is often confronted with an ad hoc flux limiter, about which physical insight might be gained by knowing the related Eddington factor. For example, we consider the sum flux-limiter proposed by Wilson,5 viz. 2(R) = 1/(3 + R).
(39)
The related Eddington factor is z(f)=(1
-f+3fz)/3,
which decreases from Z = 1/3 a t f = 0 to a minimum of X = 11/36 a t f = 1/6 and so violates condition (C3). This indicates the resulting expression for the flux probability under estimates the actual flux for gradients that are not too big (R < 3). On the other hand, the cut-off flux limiter that is essentially found in Winslow6 is 2"~"
fl/3
(a)=~I/R
if0
(40)
and has as its related Eddington factor
z(f) = (1/3) +f2, which violates condition (C1) i f f > x f ~ . This indicates the resulting expression for the flux probability over-estimates the actual flux for moderately sized gradients (2 < R < 10). It also clearly over-estimates the flux for larger gradients but the error is not serious since one presumes f is actually near 1 in that case. The flux-limiter proposed by Levermore] '8 is given by Eq. (37) and therefore suffers none of the difficulties of the previous two examples. 6. CONCLUDING REMARKS In order to complete the derivation of diffusion theory from the P-1 equations, a word should be said about boundary conditions. For the P-1 equations, some scalar combination of E and F is usually specified at each boundary point. To obtain the boundary conditions for diffusion equation merely use Eq. (32) to replace F with an expression involving E and IZEin the P-1 boundary conditions. The derived relation between Eddington factors and flux limiters, Eq. (33), was a direct
Relating Eddington factors to flux limiters
I
!
I
159
o
2
h C F
h n.
0
iv
Q
h
tl
C
"
h
oZ
--
~q L~
Q
I
I
N
I
0
--
C)
160
C . D . LEVERMORE
consequence of the new derivation of the diffusion equation presented in Section 4. A different result would be obtained if one could find a different derivation yielding a different flux-limited diffusion equation from the P-1 equations. However, even if another relation could be found, it could not be simpler than the present result and the two results would have to agree within the order of approximation over their common domain of validity. Equation (30) can be thought of as a generalization of Fick's law which allows for the possibility that the flux, F, is not in the same direction as - VE. However, it follows from the condition given in Eq. (15) that the angle between the two vectors cannot be greater than 90 ~'. This is a reflection of the fact the theory assumes some redistribution of the radiation by the matter. In the limit of a cold absorber (~o-~0) the theory becomes singular and like any other diffusion theory, can be expected to perform poorly at times. Perhaps the most striking fact about the flux-limited diffusion equations is that one can employ flux limiters that arise from Eddington factors given by single algebraic functions and obtain a well behaved theory. Yet these same Eddington factors can lead to anomalous results (shocks) when used in the P-1 equations. More work is needed to make our understanding of this phenomenon more precise.
1. 2. 3. 4. 5. 6. 7, 8. 9. 10. 11. 12. 13.
REFERENCES G. C. Pomraning, JQSRT 9, 407 (1969). G. C. Pomraning, Nukleonik 6, 348 (1965). G. N. Minerbo, JQSRT 20, 541 (1978). D. S. Kershaw, "Flux Limiting Nature's Own Way", Lawrence Livermore National Laboratory, Livermore, CA, UCRL-78378 (1976). J. M. LeBlanc and J. R. Wilson, Appl. J. 161, 541 (1970). A. M. Winslow, Nucl. Sci. Engng 32, 101 (1968). C. D. Levermore, "A Chapman-Enskog Approach to Flux-Limited Diffusion Theory", Lawrence Livermore National Laboratory, Livermore, CA, UCID-18229 (1979). C. D. Levermore and G. C. Pomraning, Appl. J. 248, 321 (1981). G. C. Pomraning, The Equations of Radiation Hydrodynamics. Pergamon Press, Oxford (1973). K. M. Case and P. F. Zweifel, Linear Transport Theory. Addison Wesley, Reading, Mass. (1967). A. Eddington, The Internal Constitution of Stars, Dover, New York (1926). K. Hoffman and R. Kunze, Linear Algebra, p. 328. Prentice-Hall, Englewood Cliffs, New Jersey (1961). R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. 2, p. 488. Wiley, New York (1962).