- Email: [email protected]
Coordinate systems for the radiative transfer equation in curvilinear media
J. Qmt. Spectrosc. Radiar. Transfer Vol.48, No. 4, PP. 4274, Printed in Great Britain. All rights reserved
1992
OC22-4073/92 S5.00 + 0.00 Copyright 0 1992 Pergamon Press Ltd
COORDINATE SYSTEMS FOR THE RADIATIVE TRANSFER EQUATION IN CURVILINEAR MEDIA PETER I3. JON@ Department of Mechanical Engineering, Auburn University, AL 36849-5341,U.S.A.
YILDIZ BAYAZITOGLU Mechanical Engineering and Materials Science Department, Rice University, Houston, TX 77251-1892, U.S.A. (Received 27 September 1991; received for publication 18 March 1992)
Abstract-Spatial/directional coordinate systems are presented on which general directional boundary conditions may be satisfied for the solution of the radiative transfer equation in curvilinear, participating media. The general directional boundary conditions described herein may be applied to curvilinear media without regard to spatial symmetry, although their application depends upon use of a local coordinate system, defining the intensity direction, which has a polar axis directed normal to any nearby surfaces. In the most commonly used spatial/directional coordinate systems for cylindrical and spherical media, the orientation of the directional polar axis indicates that the non-symmetric boundary conditions may be applied only to, respectively, parallel disk and spherical annulus boundary geometries. Application of the general directional boundary conditions to non-symmetric curvilinear media of these geometries, referred to only partially or obliquely in previous works, is demonstrated completely and consisely in the paper. Also, new spatial/directional coordinate systems are derived, based on reorientation of the local directional coordinate system in spatially cylindrical and spherical geometries. These new coordinate systems allow application of the general directional boundary conditions to other curvilinear geometries, such as a cylindrical annulus, cylindrical wedge, conical annulus, and spherical wedge. For each coordinate system, expressions for the pathlength derivative of radiation intensity are given in simple, conservative, and discrete ordinates form. The techniques illustrated in the derivation of these spatial/directional coordinate systems for non-symmetric curvilinear media may be applied to other, more complicated boundary geometries as well.
INTRODUCTION
Solution of the radiative transfer equation (RTE) in curvilinear media requires the specification of directional derivatives and directional boundary conditions. These directional terms do not appear in Cartesian media. Directional derivatives are contained in the pathlength derivative term of the RTE and are based on the particular dual coordinate system used to define spatial location and local directional orientation (spatial/directional coordinate systems or SDCSs). There are two commonly used curvilinear SDCSs, one in spatially spherical media and one in spatially cylindrical media, although others may be derived, as is shown here. Directional boundary conditions are required to solve the RTE in curvilinear media. These boundary conditions may be inferred from spatial symmetry, as in many previously presented RTE solutions in curvilinear media (e.g. spherical symmetry, axial symmetry), or from the application of general directional boundary conditions to specific cases of boundary geometry. In this paper, a fresh exposition of these general boundary conditions is presented. As a result, a necessary correlation between the specific spatial geometry of a participating medium and the directional coordinate system used to solve the RTE in that medium is demonstrated for general curvilinear cases. The specific spatial geometries suited to the commonly used SDCSs are defined, and new SDCSs are derived for other spatial geometries. The pathlength derivative terms in both the commonly used and the new SDCSs are presented in simple (analytically correct), conservative (useful in numerical solutions), and discrete ordinates fTo whom all correspondence
should be addressed. 427
PETER D. JONES and YILDIZ BAYAZITOGLU
428
(a numerical method) form. In this introduction, the general usage of curvilinear SDCSs in the RTE is discussed, general directional boundary conditions are stated, the necessary correlation between spatial geometry and directional coordinate systems is explained, and the commonly used SDCSs in spherical and cylindrical spatial geometries are presented in simple, conservative, and discrete ordinates form. In the following analysis, new SDCSs are derived and the spatial geometries for which they are suited are described. The introductory discussion sets the analytical framework for derivation of the new SDCSs. The intensity of thermal radiation, Z (or Z, if spectral problems are considered), is a quantity which varies both spatially and directionally, and is a function of up to five spatial/directional scalars (three spatial and two directional). Variations in the intensity are described by the RTE’ as follows: az/as = B(S - Z),
(1)
where S is a source function for emission and scattering, p is an extinction coefficient, and s is a directional pathlength. The pathlength derivative, aZ/as, may be expanded into coordinate form in curvilinear coordinates, leading to the expressions’
a dr -+--+--+--+--ao a dd a i%=dsar ds a4 for spherical spatial coordinates,
dl// a ds a*
dB ds
a ae
dw
a
d0
a
do
a
and
a dr a d4 a dx a -+--+rax+ds%++ds ZXar ds a$
WI
for cylindrical spatial coordinates, where (r, 4,1+9) or (r, 4, x) define the spatial location, and the locally defined polar (0) and azimuthal (0) angles define the absolute direction of s relative to the spatial location. The directional derivatives, d6/ds and dw/ds, describe the changes in orientation of the directional coordinate system with changes in spatial location. In Cartesian spatial coordinates, the directional coordinate system has a fixed orientation (de/ds and do/ds are both zero) and thus (3) For problems in spatially curvilinear geometries, as opposed to Cartesian geometries, it is necessary to explicitly solve the directional dependence of I. In order to solve the first-order RTE, one boundary condition must be specified for each scalar variable.? For the spatial variables, direction-specific boundary conditions are usually taken as the emitted and reflected intensity at a boundary surface, facing into the radiatively participating medium. Spatial boundary conditions at any particular surface are given only for the half-range of direction (0, o) which faces into the medium, away from the surface. The spatial boundary condition for the remaining half-range of (0, o) is taken at another surface, hence the spatial boundary conditions are not completely directionally independent. For the directional variables, in general cases where there is no spherical or axial symmetry to exploit, the only location-independent boundary conditions are those which arise from the self-symmetry of the directional coordinate system. For instance, the azimuthal angle has a range 0 Q w G 21c, and the boundary condition3
I@, 4, *, 8, w = 0) = Z(r, #J,$, 0, w = 2x1
(44
ensures continuity across the branch cut. The polar angle, with its range 0 < 8 < 72,is not bounded by an explicit boundary condition; however, we may make use of several conditions which lead tAlthough the source term, S, may be integral in I, many iterative numerical solutions of the RTE (including ordinates method) may be written as first order within a radiation source iteration loop.
the discrete
Coordinate systems for the RTE in curvilinear media
to the computation
429
of polar boundary values. For instance, d6’ z
= 0. 8=O.n
(W
Equation (4b) allows for the fact that for s aligned with the directional polar axis (wherever the polar axis may point relative to the spatial coordinates), spatial movement produces no reorientation of the polar axis. Equation (4b) should not be interpreted to imply that (JZ/ae) = 0 at 8 = 0, rc. We also have
-az 80
=
0,
e=O,n
where the arc traversed by do is proportional to sin 8. Frequently, enough of the spatial directional derivatives in Eqs. (2) go to zero at 8 = 0,~ [d+/ds = (l/r) sin 8 cos o, for instance], that it is possible to calculate boundary values of Z(r, 4, t,b, 8 = 0, w) and Z(r, 4, #, 8 = rr, o).~ These boundary values, with Eq. (4a), constitute the general directional boundary conditions of the RTE. In order to solve the RTE for the boundary conditions of Eqs. (4), it is necessary for Z to vary smoothly with direction. Discontinuities in the directional distribution of Zmay occur near surfaces, where the difference between incoming and outgoing intensities can be large. The value of Z can vary widely over a change in direction between angled slightly away from the surface (into the medium) and slightly toward the surface (out of the medium), especially near diffuse surfaces. Such abrupt changes can destabilize a numerical or mathematical solution algorithm. Discontinuities in the directional distribution of Z are avoided only when the directional polar axis (0 = 0) is set up as normal to any proximate surface. The sweep of the directional azimuthal angle (w) is then parallel to that surface, and the large variation of Z with direction appears only in the al/a0 (polar) term. Advantage may then be taken of the existence of boundary values of Z in both directions, 0 = 0,7c. The whole range of direction (0, o) may be split into forward and backward facing half-ranges (forward: 0 < 0 < 7~12, 0 < w < 271; backward: n/2 < 8 < K, 0 < o < 21~).In each half-range, Eqs. (4b) and (4c) may be applied to find the boundary value of Z on the positive or negative polar axis, respectively. Since the RTE as written in Eq. (1) is first order, it is not necessary to solve across the potential discontinuity at 0 = n/2. Equation (4a), which must operate on the whole (undivided) range of o, is used to complete the solution Thus, in the absence of simplifying spatial symmetry or additional boundary values, the curvilinear RTE may only be solved using the general boundary conditions of Eqs. (4), and then only if the RTE is
Fig. 1. (a) Geometry of the SphRad spatial/directional coordinate system. (b) Spherical annulus geometry, which is well suited for non-symmetric application of the SphRad spatial/directional coordinate system.
430
PETERD. JONESand YILDIZBAYAZITOGLLJ
in an SDCS which is set up so that the directional equitorial plane coincides with the spatial boundary surfaces (i.e. the directional polar axis is normal to the surface). The most commonly used SDCS in spatially spherical media is suited by the latter condition to the spatial geometry of a spherical annulus [surfaces of constant radius, Fig. l(b)]. The pathlength derivative for this SDCS is written as’ written
az i
az as
-=cosO,--+-sin8,cosw,---+
ar
r
az w
sin O,sino, r sin 4
az 1 sin 8 g _ sin 8, sin 0, az --rae, rtan+ z’ a* r
(5a)
where the subscript r denotes that the directional polar (0,) and azimuthal (w,) angles are based on a polar axis which is an extension of the spatial radial coordinate [Fig. l(a)]. Thus the boundary conditions of Eqs. (4) may be applied in a spherical annulus geometry, as the directional equitorial plane lies in the surface. This SDCS is here termed ‘spatially spherical, radial polar’, or SphRad, to distinguish it from SDCSs with other spatial geometry (cylindrical) or polar axis orientations (polar and azimuthal in spatially spherical geometry, axial and angular in spatially cylindrical geometry). Terminology for the other SDCSs described in this paper is given in Table 1. In the analysis which follows, analogs to Eq. (5a) will be derived for geometries defined in spatially spherical coordinates which have surfaces defined by constant spatial polar or spatial azimuthal angles, as opposed to surfaces defined by constant radius. Equation (5a) may also be written in conservative form (‘conservative’ refers to conservation of neutrons; the RTE is similar in form to the neutron transport equation) by multiplying Z by d V dR, where dV = r2 sin 4 dr d+ d$ (differential element of volume) and dR = sin I3de do (differential element of solid angle), carrying out the differentiations, and dividing by dV dR. The coordinate derivatives (dr/d.s, de,/&, etc.) are also taken within the partial derivatives, giving the result2~5~6
az i - = - cos 0,: (Zr*) as r*
+sinre~io~Or
-& (Z sin 4) + sinr~i~~wr $ -&&(Zsin*B,)-s-&(Zsinw,).
(5b)
Conservative form is useful in numerical solutions, where it helps to regulate errors due to a finite mesh size. Equations (5a) and (5b) can be shown to be equivalent. This equivalence is a necessary condition for the conservative form, and serves as both a derivation tool and check for the new SDCSs presented in the following analysis. Direct numerical solution of the RTE may be achieved using the discrete ordinates method,2,4.6-* which is essentially a finite difference representation using a quadrature of points and integration weights for the directional mesh. Using a finite volume/finite solid angle cell approach,4q8 the conservative RTE may be integrated over the finite cell 6V6R to obtain the finite cell discrete Table 1. Nomenclature identifying the spatial/directional coordinate systems @DC.%). Abbreviation CylAXl+
I
Description spatially cylindrical, directional polar axin aliied
with axial spatial coordinate
W-d
rpatially cylindrical, directional polar axia aligned with radial spatial coordinate
cr-e
spatially cylindrical, directional polar axis aligned with angular spatial coordinate
SphRad +
rpatially 8phuica1, directional polar axiS aligned with radial spatial coordinate
SphPal
rpatially spherical, diiectional polar axb aligned with polar spatial coordinate
SphAsi
rpatially rpherical, ditional polar axb~ aligned with azimuthal spatial coordinate
tCommonly used coordinate system.
Coordinate systems for the RTE in curvilinear media
431
ordinates form of the pathlength derivative, written dy M1; dI’ d(l = LIj(Zrz)6,( -cos 4) ~,(~)cos(~,,)w, is +f sj(r2)6k(Zsin~)6,(~)sin(e,,)w,cos(w,.)w,
w,
+:sj(r2)Bk(~)6r(Z)sin(e,,)w,sin(o,.)w,
-:6j(r2)8k(-COS~)6r(~)6,(Z0,)w,
-~6j(r2)Bk(sin~)6,(~)sin(e,,)w,6,(ZR,),
(W
where the mesh indicies are rj, &, $,, 8,,, co,,. . The nomenclature Sj(Zr2) implies a point-centered cell difference over a single index, 4+,,2,k,,,m,nrj+,,2 - h_. ,,2,k,,,m,nrf_,,2, and w,,, and w, are integration weights for integrals over cells of directional angles,‘j viz. f(e) sin 8 de =f(&,)w,, 6. f(w) do =fbJ~n. WI s n s wn The suggestion of Abu-Shumaysg is followed, making the integration weight for the quadrature point (4 n), Wnl,“,equal to w,,,w,,, where w, and w, are independent. In Eq. (5c), the coefficients 0, and Q occupy places analogous to those of sin* 8, and sin o, in Eq. (5b). The reason for the substitution is common to all discrete simulations of the conservative RTE in curvilinear media, namely that making a direct discretization of Eq. (5b) results in an expression where al/as # 0 for constant Z, except when very restrictive conditions are placed on the spatial mesh and directional quadrature.‘v3 The common solution to this difficulty is to substitute the discrete coefficients 0, and R, for the analytical elements in the directional terms of Eq. (5b), where 0, and Q are derived specifically to allow the condition al/as = 0 for constant Z to be met for the conservative form. For the SphRad SDCS, the coefficients are defined by2v3 a,(@,) = 2w,c0se,,,, o~I+.~= 0, ~,(Q)=w,COS~,,m
(54
QL,=o=O,
Of)
where the initial values lie on particular directional cell boundaries, and Eqs. (se) and (5f) then provide recursion equations for 0, and Q on the remaining cell boundaries (m + l/2 or n f l/2). The commonly used SDCS in spatially cylindrical media has the directional polar axis aligned with the spatial axial coordinate [Fig. 2(a)]. This leads to the expression
az
az
as’
I
az
sin8,coso,--+-sinB,sino,-+cose,---sin8,sinw,---, ar r 84
(a)
az
i
ax
r
az aox
, lb)
Fig. 2. (a) Geometry of the CylAxi spatial/directional coordinate system. (b) Parallel disk geometry, which is well suited for non-symmetric application of the CylAxi spatial/directional coordinate system.
(64
432
PETERD. JONFSand YILDIZ BAYAZITOGLU
where the subscript x describes the polar axis alignment. 2*5,6 This expression has the particular utility that d6,/ds is uniquely zero, as the polar axis (spatial axial coordinate) does not change its orientation with changes in s. Therefore, only one directional boundary condition, Eq. (4a), is required. This spatially cylindrical, axial polar (CylAxi) SDCS has been applied to problems with the spatial geometry of a cylindrical annulus, Fig. 3(b). ‘OJ’In this geometry, the sweep of o, near a surface passes from a direction facing into the medium to a direction facing out, leading to a discontinuity in al/&, which is difficult to bridge numerically, and which makes the general boundary condition, Eq. (4a), very difficult to apply. This difficulty is avoided if the problem is spatially axially symmetric. In that case, there is a further directional boundary condition: al/&~, = 0 for o, such that ds lies in the plane of spatial axial symmetry. The directional boundary values may then be split between forward and backward facing half-ranges, and the discontinuity at the half-range boundaries is avoided. Alternatively, in non-symmetric cases, the values for the polar directional distribution can sometimes be used to generate a line of azimuthal boundary values in each half range, dependent on the polar angle, which may be used as substitutes for Eq. (4a). This approach was successfully used in Refs. 10 and 11 to solve the RTE in cylindrical annular media, using a SDCS (CylAxi) which would otherwise be ill-suited for this geometry. The CylAxi SDCS is best suited to a geometry of facing, parallel disks [Fig. 2(b)], where the directignal polar axis is normal to any surface, and the variation of Z with o, is everywhere continuous. Equation (6a) is written in conservative form in a manner analogous to that for spatially spherical geometry, using d V = r dr d4 dx and dR = sin 9, do, dw,, and resulting in
az i = - sin 0, cos w, g (Zr) + i sin 8, sin 0, z + cos 8, g - b sin 0, & as r 84
-
x
(Isin 0,).
(6b)
Equation (6b) satisfies the necessary condition of equivalence to simple form, Eq. (6a). Equation (6b) may be further adapted to a finite volume/finite solid angle cell discrete ordinates form in a manner similar to that for Eq. (5~) with the result
,& ss 6y
dV dQ= 6jUr)h(4) ~I(x)sin(e,,)w,cos(~,.)w, +
dj(r)
b(Z) 6Win(kJw,
W+,k, (6~)
where
(b)
Fig. 3. (a) Geometry of the CylRad spatial/directional coordinate system. (b) Cylindrical annuius geometry, which is well suited for non-symmetric application of the CylRad spatial/directional coordinate system.
433
Coordinate systems for the RTE in curvilinear media
The foregoing introductory discussion has outlined the occurrence of directional derivatives in the RTE in curvilinear media, presented general boundary conditions for these terms, shown the form of the pathlength derivative in the SDCSs most commonly employed for spatially spherical and cylindrical geometries, with finite cell discrete ordinates expressions, and discussed the boundary condition and continuity difficulties associated with using these SDCSs in ill-suited spatial geometries. For practical solution of the RTE in non-symmetric curvilinear media, it is often advantageous to use an SDCS suited to the geometry (with the directional polar axis normal to any surface) in order to apply the general boundary conditions of Eqs. (4). However, to the knowledge of the authors, the only SDCSs for which the pathlength derivative in the RTE has been derived are what have been termed the ‘usual’ SDCSs, SphRad and CylAxi. In the following analysis, alternative SDCSs are proposed which are suited for particular spherical and cylindrical spatial geometries. Pathlength derivative expressions are derived for these SDCSs in simple, conservative, and finite cell discrete ordinates form. ANALYSIS Spatially cylindrical, radial polar
An SDCS more suitable for solution of the RTE in cylindrical annular geometries than the commonly used CylAxi is the spatially cylindrical, radial polar SDCS [CylRad; Fig. 3(a)], which is based on a directional polar axis aligned with the spatial radial direction. With the CylRad SDCS, the directional equitorial plane always coincides with the surfaces of a cylindrical annulus, allowing the boundary conditions of Eqs. (4) to be applied. CylRad allows boundary values to be determined for 8, = 0, II using Eqs. (4b, c), while the distribution of Z with o, is uninterrupted and Eq. (4a) may be applied. The CylRad SDCS may be derived as a transformation from CylAxi by expressing the direction of the differential pathlength ds for each SDCS in the local spatial coordinate system (i, &a) as follows: ds -= sin 8, cos o,i + sin 0, sin w,f$ + cos 8,?, IdsI (74 ds -= cos 8, i - sin 8, sin W,fj + sin 8, cos CO, i, IdsI where the subscript x denotes polar (0) and azimuthal (0) directional angles based on a polar axis aligned with the local spatial axial component (CylAxi), while the subscript r denotes a polar axis aligned with the local spatial radial component (CylRad). Equating the s derivatives of the spatial component terms in Eqs. (7a) and using the CylAxi definitions
(7’4 from Eq. (6a), the CylRad directional derivatives may be derived, resulting in
de,_ ds - -;
do, 1 sin 8, sin’ w, , -=
1 --cos~,sino,coso,. r
ds
Using the directional derivatives and spatial components, the simple form of the CylRad pathlength derivative may be assembled similarly to Eq. (2b), resulting in
az
az i ar r
-=cos&---sin&sinw,-+sinf+cos~-
as
az
az
w
ax 1 --
r
sin
az i 8, sin’ 0, -ae, - -r cos 8, sin 0, cos W,-;;
7
.
(74
Following the methods demonstrated for CylAxi, the conservative form for CylRad is found to be 1 az i a -=-cos&--(Zr)--sintI,sino,-+sintJ,coso,-
as
r
ar
r
az w -S&(Z
az ax sin2 e,) - i cos 8, -& (I sin 0, cos o,), I
(7e)
PETERD.
434
JONFS and YILDIZ BAYAZIT~CLU
where the necessary equivalence between Eqs. (7d) and (7e) is satisfied. Analytically, either Eq. (7d) or (7e) is correct; however, Eq. (7e) is generally more stable in numerical simulations. Again following the methods demonstrated for CylAxi, the finite cell discrete ordinates form of the pathlength derivative for CylRad is (7f)
6v d~~dvda=6,(Z~)~,(~)6,(x)cos(e~,~)W~W~ ss -
6j(rPk(1)
Uxbin(%Jw,
sMw..N,
+~6j(r2)bk(~)8,(l)sin(B,,)w,cos(o,.)w,
-
6j(r)
4A4)
4(X)
4S@,)sin2b4,.h
-
6j(r)
6k(+)
G,(x)cos(e,~)w~6~(zn~)~
where in order to satisfy the necessary condition that al/as = 0 for constant Z, we choose d,(@,)
=
Wrn
cos(eLm)@lo,=ll= 0,
sin2(w, ,) ’
W4) = 0, Q,.,
l/2 =
4
sin(~&os(qn).
Ug) V’h)
Equations (7g) are the recursion and intitial value for 0, in the manner of Eqs. (5e), (5f), and (6d). Equations (7h) indicate that, in order to satisfy the necessary condition, Q is not allowed to vary within the finite cell; however, R, is discontinuous at the cell boundaries n &-l/2, and must vary from cell to cell in order to maintain the correspondence between Eqs. (7e) and (7f). Note that the sin* (w,.) term in the fourth line of Eq. (7f) is canceled by the choice of 6,(0,), and that the combination of Eqs. (7f) and (7g) does not contain the sin2(o,J term. In order to use general directional boundary conditions, without recourse to the spatial symmetry of a medium, it is necessary to use an SDCS with a directional polar axis normal to the boundary surfaces. The CylRad SDCS has a directional polar axis which is normal to surfaces of constant radius, e.g. cylindrical tubes. Hence, the CylRad SDCS is appropriate for solution of the RTE in cylindrical annular media [Fig. 3(b)]. Spatially cylindrical, angular polar
The spatially cylindrical coordinate system may be used to define other geometries besides a cylindrical annulus. For instance, a medium bounded by an array of rectangular plates set at angles to each other and hinged at a common point [cylindrical wedge, Fig. 4(b)] is readily described in
(b)
Fig. 4. (a) Geometry of the CylAng spatial/directional coordinate system. (b) Cylindrical wedge geometry, which is well suited for non-symmetric application of the CylAng spatial/directional coordinate system.
435
Coordinate systems for the RTE in curvilinear media
cylindrical coordinates. For such a geometry, the boundary surfaces are defined by constant angular coordinates, and the appropriate SDCS has a directional polar axis aligned with the spatial angular direction [Fig. 4(a)]. The pathlength derivative for this spatially cylindrical, angular polar (CylAng) SDCS is derived using ds -= IdsI
sin
8, cos o4 i +
cose,$ - sin 8, sin o&
with the result expressed in simple form
ar
CYI
I
sin 8~~0s~ - +-case @ar r as=
ar
az
- -sin 8, sino "ax 4a4
1 az +; cose,$cosw~~-
sin cf+ -- 1 r sin 8,
@b)
and in conservative form
ar -=as
1 az a az sin e,+cos0,--(Zr)+- cos e - -sin 8, sin 0 r #ax & r $84 I
sin e In finite cell discrete ordinates form, the CylAng pathlength derivative dy ,,!dV Is
da = 6jW) b(4)
4(~)W&,h,
sin o++,). (8~)
is expressed
cos(oo,,h
+~~(~)~~(Z)GI(X)COS(~~,,)W,W, -:sj(r2)8~(~)6,(l)sin(e~,,)w,sin(o~,.)w, +6j(r)6,(~)6,(x)G,(lO,)cos(o,.)w,
In order to satisfy the necessary condition, we choose &(@,) = - w, sin(&,,), WV
= 0,
a,,,
1/2 =
w,
@61s,=n,2= 0,
(8e)
W+,.),
WI
where Eqs. (8e) are a recursion and initial value set and Eqs. (8f) define constant cell values for R,, as described in the previous section. In addition to CylAxi, CylRad, and CylAng, further SDCSs may be developed for spatially cylindrical geometries. These may be derived by following the techniques demonstrated above, as long as the directional polar axis is kept normal to the medium’s boundary surfaces. For instance, a cylindrical spiral surface may be defined by ds -= IdsI
cos e,(cosait+ sin ~(4) + sin 8,sin o,( -sin c& + cos a$) + sin 8, cos 0,i
(9)
where a is the angle between the spatial axial coordinate and the normal to the spiral surface. Using Eq. (9) and following the steps outlined above, a pathlength derivative expression may be derived for a medium bounded by this cylindrical spiral surface. Spatially spherical polar polar
The spatial spherical coordinate system may be applied to other geometries besides spherical annuli. For these other geometries, it is necessary to use a directional coordinate system with a polar axis normal to the boundary surfaces. For example, if the boundary surfaces in a spherical coordinate system are defined as surfaces of constant spatial polar angle, over a range of radius and spatial azimuthal angle, then a geometry of the annulus between cones of different apex angle
436
PETERD. JONESand YILDIZ BAYAZIT~GLU
Fig. 5. (a) Geometry of the SphPol spatial/directional coordinate system. (b) Conical annulus geometry, which is well suited for non-symmetric application of the SphPol spatial/directional coordinate system. is described for solution The SphRad system,
[Fig. S(b)]. In order to use this SDCS, spatially spherical polar polar [SphPol, Fig. 5(a)], of the RTE, it is necessary to derive a new expression for the pathlength derivative. and SphPol pathlength directions, relative to the (i, $, $) local spatial coordinate
are
sin 8, sin 0, sink II/,
ds
1
ds
^ sin + cos e,4 -
(ds=cosB~i+sin0,cosw,$ ldsl = sin 8, cos 0,:
+
e. cos wg
-
*9
sin 4
which, when combined with the SphRad directional derivatives from Eq. (Sa), sin 8, sin w, r tan4 ’
dw, -=ds
dB, _ 1 . sin%,, ds - -;
(lob)
result in the SphPol directional derivatives de, _= ds
-- 1 cos 8, cos wg r
(1Oc)
Using Eqs. (lOa) and (lOc), the simple form of the pathlength assembled, viz.
az
az
- =sin04c0so as
-
+ar
i
+-00se r
az
+a4
-
derivative in SphPol may be
sin 8, sin o4 az rsin4 ZJ
sin 8, sin’ o+ --cos 8, cos We az do. (104 tan I$ tan $J > > 4 The conservative form may be determined by multiplying Z by dV = r2 sin 4 dr d4 d$, dR = sin 8, de, dw,, and the directional derivatives, differentiating, and dividing by d V dR, with the result
az
I
- =-sin0#coso as r2 COSW+ a + ~-(zsin86c0se,)r sin 8, ae, __qz. 1 r sin e+ am,
z(Zrz)+
2
-.$ (1 sin 4) _ sinre~i~~w,
$
4 ar
sin2 w+
r
sin o,+) -fg$&(’
+z sin e4 tan C#Jae,
sin’ ( 8,)
sin o4 cos 0~~).
UW
437
Coordinate systems for the RTE in curvilinear media
Generation of the finite cell discrete ordinates form for SphPol is straightforward, complicated by the multiple directional partial derivatives in Eq. (10e). The result is
ss ar
bv mz
dV da = dj(Zr*> U-~0s 4) 4($>sW,,,)w, +
fsj
43
though
cos(w&w,
sin4) 4(ll/)COs(~,,)w,w,
- f sj(r2) U+) UOsin(~~,,)w,sin(q,.)w, + i6jCr2) sk(-cos4) 6!($)~~(Z@I~)COS(WI$,~)W~ - 46j(r2) Usin 4) h(JI>hAZ@24)sin2(w~,,)w, -
i6jCr2)
wm
d,(IR,,) 6k(-cos4) a,(+> SW, m ) (IOf)
- i 6j(r2) Msin 4) 4(II/)cos(B$,,)w, 6,(ZR,,), where we choose
4d%J = -w, sin(e+J, U@*,)
= w,
MT+)
= 0,
4 (% ) = 0,
cos(4, m) sin2(04,, ) ’
o,,I,,=,, @*&dl=
= 0,
wk)
0,
(1W (1Oi)
fi~,~+ ,,2= w, sin(04,.), Q24,n+ 1,2= w, side,
, )cos(w+, .).
(W)
Spatially spherical, azimuthal polar
Construction of an SDCS in spatially spherical media where the directional polar axis is oriented along the direction of the local spatial azimuthal coordinate results in the spatially spherical, azimuthal polar (SphAzi) SDCS, Fig. 6(a), which may be applied to a spherical wedge geometry, Fig. 6(b). The SphAzi pathlength direction, relative to the local directional coordinate system V, $, ~0, is ~0~8, ds -=sin8ticoswtii+sin8+sino+$ - ---$. (114 sin 4 IdsI The SphAzi directional derivatives may be solved using Eqs. (5a), (lob), and (llc),
leading to
(11’4 COS~,COS
we
sin w+
which results in expressions for the pathlength derivative in simple form,
az cos 8, az az i az -_= sin 8, cos w - + - sin 8,sin ws - - a4 r sin4 a$ +iYr r as
COS~+COS W+
sin
W$
tan 8,tan f#~ and, in conservative form, az i as =7sin8,coswti~(IrZ)+
sinre~i~~w’ $
(I sin 4) + z
sin wti 1 + COSW# -&(I sin e,cose,) r sin 8,( tan f#~ ) 1 COSTS sin We). +z +z cosw,)-+ r sin 8, aw, r tan 8,tan 4 aa,
$
+-
(114
PETER D. JONES and YILDIZ BAYAZIT~GLU
438
(b)
Fig. 6. (a) Geometry of the SphAzi spatial/directional coordinate system. (b) Spherical wedge geometry, which is well suited for non-symmetric application of the SphAzi spatial/directional coordinate system.
Generation
development
of the finite cell discrete ordinates form of the pathlength derivative is similar to the for other SDCSs, with the result
az 6v 6n~dVdR=S,(z~2)b~(-cos~)g,(~)sin(8l,~)w~cos(w~,~)w, ss
+ + 6j(r2) +f +
&(Z sin 4) 8,(ll/)sin(O,,,)w,
sin(o+)w,
~j(r2)~,(~)~(~)cos(~,,,)w,w,
i6j(r2)k%(sin
+ i 6j(r2) 4&n - f
6j(r2)
4)sin(o,,.)
4) 4(+)
+
[
6,(-cos
1 sin(8e,,) - sin(%.,)
fik(-COs4) d,($)
1
4)cos(o,,.)]
S,($) &(ZO,)w,
MQ,)
wm6,(ZR,,),
(1W
sln(%, 1 m
where we choose 6,(0,)
[
= W, c?y:Vm)) - sin(&,) .m
MA,)
= - w,sin(f+,.),
hW2,)
= w,cos(q,),
= 1, 0,IBe
QI,I,,=,!2 = 0, ~,,I,,=,=@
x,2= 03
(llf> (1w (1 lh)
Each of Eqs. (1 lf, g, h) is a recursion and initial value set. In order to use the SphAzi SDCS with the directional polar axis normal to the medium’s boundary surfaces, it should be applied to a
Coordinate systems for the RTE in curvilinear media
439
Table 2. Summary of spatial/directional coordinate systems (SACS) results for the SLXS defining geometry; an example of an application; and the pathlength derivative in simple analytic form, conservative analytic form, and finite cell discrete ordinates form.
Result
SphRad
BphPdlspMrl
EF--ry
Fig.2a
Fig.30
Fig.4a
Fig.la
Fig.Sa
Fig.Ba
application
Fig.2b
Fig.3b
Fig.lb
Fig.lb
Fig.Sb
Fig.6b
(N/ad),
simple
%.@a)
JW’ld)
Eg.(8b)
E&a)
Eq.(lOd) Eq.(llc)
(N/ad),
conae-rvative
EOb)
Eq.(‘lc)
Eq.(8c)
E&b)
Eq.(lOe)
EgW
EgOf)
Eg.(8d)
JQ.(Se)
Eq.(lOf) Eq.(llc)
/,v &,, (aI/ad)mdn, discrete ordinates
Eq.(lld)
spatial geometry such as plane surfaces set at an arbitrary angle, much as for CylAng, but with the spatial location defined in spherical rather than cylindrical coordinates [Fig. 6(b)]. The SphAzi SDCS was investigated in a preliminary manner by Jones,3 although the form of the pathlength derivative derived in that work cannot be developed into a conservative form which satisfies the necessary condition that al/as = 0 for constant I. This condition is met in the earlier work only if Z is multiplied by dV rather than dV dR. SUMMARY The radiative transfer equation in curvilinear coordinates requires directional boundary conditions for its solution. For particular problems, it may be possible to derive directional boundary conditions which result from geometrical symmetry. However, in general cases without useful symmetry, the only available boundary conditions are Eq. (4a) and the boundary values which may be developed using Eqs. (4b) and (4c). In order to apply these general boundary conditions, the coordinate system which defines the direction of the radiation intensity must be set up in such a way that its polar axis is normal to the spatial boundaries, especially when in proximity to the boundary surfaces. Therefore, for a particular local orientation of the directional coordinate axes, only very specific spatial boundaries may be defined for use with the general directional boundary conditions. In the commonly used spatial/directional coordinate systems, widely presented in the literature and here termed CylAxi (spatially cylindrical, directional polar axis aligned with the axial spatial coordinate) and SphRad (spatially spherical, directional polar axis aligned with the radial spatial coordinate), the general directional boundary conditions may only be applied in spatial geometries of facing parallel disks (CylAxi) and spherical annuli (SphRad). In the present work, spatial/directional coordinate systems are derived in which the general directional boundary conditions may be applied (and the radiative transfer equation solved) for other spatial geometries. Intensity pathlengths derivatives, which contain the directional derivatives, are derived for each new spatial/directional coordinate system in simple, conservative, and finite cell discrete ordinates form. These results are summarized in Table 2. Spatial/directional coordinate systems may also be developed for spatial geometries defined by a combination of spatial coordinates, as opposed to the geometries presented here, which are defined by a single coordinate held constant. Using the techniques outlined in the present work, intensity pathlength derivatives may be derived for such coordinate systems. The radiative transfer equation may then be solved in these geometries using the general directional boundary conditions, unaided by boundary conditions derived from geometrical symmetry. REFERENCES 1. M. N. &igik, Radiative Transfer and Interactions with Conduction and Convection, Wiley, New York, NY (1973) [reprint: Werbel & Peck, New York, NY (1985)]. 2. E. E. Lewis and W. F. Miller Jr., Computational Methods of Neutron Transport, Wiley, New York, NY (1984). QSRT
48/4-G
PETERD.
440
3. P. D. Jones, “Radiation
4. 5. 6.
7. 8. 9.
10. 11.
JONES and YILDIZ BAYAZITOGLU
and Convection Heat Transfer in Particle-Laden Fluid Flow,” Ph.D. Thesis, Rice University (1990). J. R. Tsai, M. N. ozivik, and F. Santarelli, JQSRT 42, 187 (1989). R. Viskanta and M. P. Mengiic, Prog. Energy Combust. Sci. 13, 97 (1987). B. G. Carlson and K. D. Lathrop, in Computing Methods in Reactor Physics, H. Greenspan, C. N. Kelber, and D. Okrent eds., Gordon & Breach, New York, NY (1968). W. A. Fiveland, J. Heat Transfer 106, 699 (1984). W. A. Fiveland, J. Thermophys. Heat Transfer 2, 309 (1988). I. K. Abu-Shumays, Nucl. Sci. Engng 64, 299 (1977). A. Yucel and M. L. Williams, in ASME Proceedings of the National Heat Transfer Conference, HTD-Vol. 96, 1, H. R. Jacobs ed., American Society of Mechanical Engineers, New York, NY (1988). A. S. Jamalludin and P. J. Smith, in ASME Proceedings of the National Heat Transfer Conference, HTD-Vol. %, 1, H. R. Jacobs ed., American Society of Mechanical Engineers, New York, NY (1988).
1992
OC22-4073/92 S5.00 + 0.00 Copyright 0 1992 Pergamon Press Ltd
COORDINATE SYSTEMS FOR THE RADIATIVE TRANSFER EQUATION IN CURVILINEAR MEDIA PETER I3. JON@ Department of Mechanical Engineering, Auburn University, AL 36849-5341,U.S.A.
YILDIZ BAYAZITOGLU Mechanical Engineering and Materials Science Department, Rice University, Houston, TX 77251-1892, U.S.A. (Received 27 September 1991; received for publication 18 March 1992)
Abstract-Spatial/directional coordinate systems are presented on which general directional boundary conditions may be satisfied for the solution of the radiative transfer equation in curvilinear, participating media. The general directional boundary conditions described herein may be applied to curvilinear media without regard to spatial symmetry, although their application depends upon use of a local coordinate system, defining the intensity direction, which has a polar axis directed normal to any nearby surfaces. In the most commonly used spatial/directional coordinate systems for cylindrical and spherical media, the orientation of the directional polar axis indicates that the non-symmetric boundary conditions may be applied only to, respectively, parallel disk and spherical annulus boundary geometries. Application of the general directional boundary conditions to non-symmetric curvilinear media of these geometries, referred to only partially or obliquely in previous works, is demonstrated completely and consisely in the paper. Also, new spatial/directional coordinate systems are derived, based on reorientation of the local directional coordinate system in spatially cylindrical and spherical geometries. These new coordinate systems allow application of the general directional boundary conditions to other curvilinear geometries, such as a cylindrical annulus, cylindrical wedge, conical annulus, and spherical wedge. For each coordinate system, expressions for the pathlength derivative of radiation intensity are given in simple, conservative, and discrete ordinates form. The techniques illustrated in the derivation of these spatial/directional coordinate systems for non-symmetric curvilinear media may be applied to other, more complicated boundary geometries as well.
INTRODUCTION
Solution of the radiative transfer equation (RTE) in curvilinear media requires the specification of directional derivatives and directional boundary conditions. These directional terms do not appear in Cartesian media. Directional derivatives are contained in the pathlength derivative term of the RTE and are based on the particular dual coordinate system used to define spatial location and local directional orientation (spatial/directional coordinate systems or SDCSs). There are two commonly used curvilinear SDCSs, one in spatially spherical media and one in spatially cylindrical media, although others may be derived, as is shown here. Directional boundary conditions are required to solve the RTE in curvilinear media. These boundary conditions may be inferred from spatial symmetry, as in many previously presented RTE solutions in curvilinear media (e.g. spherical symmetry, axial symmetry), or from the application of general directional boundary conditions to specific cases of boundary geometry. In this paper, a fresh exposition of these general boundary conditions is presented. As a result, a necessary correlation between the specific spatial geometry of a participating medium and the directional coordinate system used to solve the RTE in that medium is demonstrated for general curvilinear cases. The specific spatial geometries suited to the commonly used SDCSs are defined, and new SDCSs are derived for other spatial geometries. The pathlength derivative terms in both the commonly used and the new SDCSs are presented in simple (analytically correct), conservative (useful in numerical solutions), and discrete ordinates fTo whom all correspondence
should be addressed. 427
PETER D. JONES and YILDIZ BAYAZITOGLU
428
(a numerical method) form. In this introduction, the general usage of curvilinear SDCSs in the RTE is discussed, general directional boundary conditions are stated, the necessary correlation between spatial geometry and directional coordinate systems is explained, and the commonly used SDCSs in spherical and cylindrical spatial geometries are presented in simple, conservative, and discrete ordinates form. In the following analysis, new SDCSs are derived and the spatial geometries for which they are suited are described. The introductory discussion sets the analytical framework for derivation of the new SDCSs. The intensity of thermal radiation, Z (or Z, if spectral problems are considered), is a quantity which varies both spatially and directionally, and is a function of up to five spatial/directional scalars (three spatial and two directional). Variations in the intensity are described by the RTE’ as follows: az/as = B(S - Z),
(1)
where S is a source function for emission and scattering, p is an extinction coefficient, and s is a directional pathlength. The pathlength derivative, aZ/as, may be expanded into coordinate form in curvilinear coordinates, leading to the expressions’
a dr -+--+--+--+--ao a dd a i%=dsar ds a4 for spherical spatial coordinates,
dl// a ds a*
dB ds
a ae
dw
a
d0
a
do
a
and
a dr a d4 a dx a -+--+rax+ds%++ds ZXar ds a$
WI
for cylindrical spatial coordinates, where (r, 4,1+9) or (r, 4, x) define the spatial location, and the locally defined polar (0) and azimuthal (0) angles define the absolute direction of s relative to the spatial location. The directional derivatives, d6/ds and dw/ds, describe the changes in orientation of the directional coordinate system with changes in spatial location. In Cartesian spatial coordinates, the directional coordinate system has a fixed orientation (de/ds and do/ds are both zero) and thus (3) For problems in spatially curvilinear geometries, as opposed to Cartesian geometries, it is necessary to explicitly solve the directional dependence of I. In order to solve the first-order RTE, one boundary condition must be specified for each scalar variable.? For the spatial variables, direction-specific boundary conditions are usually taken as the emitted and reflected intensity at a boundary surface, facing into the radiatively participating medium. Spatial boundary conditions at any particular surface are given only for the half-range of direction (0, o) which faces into the medium, away from the surface. The spatial boundary condition for the remaining half-range of (0, o) is taken at another surface, hence the spatial boundary conditions are not completely directionally independent. For the directional variables, in general cases where there is no spherical or axial symmetry to exploit, the only location-independent boundary conditions are those which arise from the self-symmetry of the directional coordinate system. For instance, the azimuthal angle has a range 0 Q w G 21c, and the boundary condition3
I@, 4, *, 8, w = 0) = Z(r, #J,$, 0, w = 2x1
(44
ensures continuity across the branch cut. The polar angle, with its range 0 < 8 < 72,is not bounded by an explicit boundary condition; however, we may make use of several conditions which lead tAlthough the source term, S, may be integral in I, many iterative numerical solutions of the RTE (including ordinates method) may be written as first order within a radiation source iteration loop.
the discrete
Coordinate systems for the RTE in curvilinear media
to the computation
429
of polar boundary values. For instance, d6’ z
= 0. 8=O.n
(W
Equation (4b) allows for the fact that for s aligned with the directional polar axis (wherever the polar axis may point relative to the spatial coordinates), spatial movement produces no reorientation of the polar axis. Equation (4b) should not be interpreted to imply that (JZ/ae) = 0 at 8 = 0, rc. We also have
-az 80
=
0,
e=O,n
where the arc traversed by do is proportional to sin 8. Frequently, enough of the spatial directional derivatives in Eqs. (2) go to zero at 8 = 0,~ [d+/ds = (l/r) sin 8 cos o, for instance], that it is possible to calculate boundary values of Z(r, 4, t,b, 8 = 0, w) and Z(r, 4, #, 8 = rr, o).~ These boundary values, with Eq. (4a), constitute the general directional boundary conditions of the RTE. In order to solve the RTE for the boundary conditions of Eqs. (4), it is necessary for Z to vary smoothly with direction. Discontinuities in the directional distribution of Zmay occur near surfaces, where the difference between incoming and outgoing intensities can be large. The value of Z can vary widely over a change in direction between angled slightly away from the surface (into the medium) and slightly toward the surface (out of the medium), especially near diffuse surfaces. Such abrupt changes can destabilize a numerical or mathematical solution algorithm. Discontinuities in the directional distribution of Z are avoided only when the directional polar axis (0 = 0) is set up as normal to any proximate surface. The sweep of the directional azimuthal angle (w) is then parallel to that surface, and the large variation of Z with direction appears only in the al/a0 (polar) term. Advantage may then be taken of the existence of boundary values of Z in both directions, 0 = 0,7c. The whole range of direction (0, o) may be split into forward and backward facing half-ranges (forward: 0 < 0 < 7~12, 0 < w < 271; backward: n/2 < 8 < K, 0 < o < 21~).In each half-range, Eqs. (4b) and (4c) may be applied to find the boundary value of Z on the positive or negative polar axis, respectively. Since the RTE as written in Eq. (1) is first order, it is not necessary to solve across the potential discontinuity at 0 = n/2. Equation (4a), which must operate on the whole (undivided) range of o, is used to complete the solution Thus, in the absence of simplifying spatial symmetry or additional boundary values, the curvilinear RTE may only be solved using the general boundary conditions of Eqs. (4), and then only if the RTE is
Fig. 1. (a) Geometry of the SphRad spatial/directional coordinate system. (b) Spherical annulus geometry, which is well suited for non-symmetric application of the SphRad spatial/directional coordinate system.
430
PETERD. JONESand YILDIZBAYAZITOGLLJ
in an SDCS which is set up so that the directional equitorial plane coincides with the spatial boundary surfaces (i.e. the directional polar axis is normal to the surface). The most commonly used SDCS in spatially spherical media is suited by the latter condition to the spatial geometry of a spherical annulus [surfaces of constant radius, Fig. l(b)]. The pathlength derivative for this SDCS is written as’ written
az i
az as
-=cosO,--+-sin8,cosw,---+
ar
r
az w
sin O,sino, r sin 4
az 1 sin 8 g _ sin 8, sin 0, az --rae, rtan+ z’ a* r
(5a)
where the subscript r denotes that the directional polar (0,) and azimuthal (w,) angles are based on a polar axis which is an extension of the spatial radial coordinate [Fig. l(a)]. Thus the boundary conditions of Eqs. (4) may be applied in a spherical annulus geometry, as the directional equitorial plane lies in the surface. This SDCS is here termed ‘spatially spherical, radial polar’, or SphRad, to distinguish it from SDCSs with other spatial geometry (cylindrical) or polar axis orientations (polar and azimuthal in spatially spherical geometry, axial and angular in spatially cylindrical geometry). Terminology for the other SDCSs described in this paper is given in Table 1. In the analysis which follows, analogs to Eq. (5a) will be derived for geometries defined in spatially spherical coordinates which have surfaces defined by constant spatial polar or spatial azimuthal angles, as opposed to surfaces defined by constant radius. Equation (5a) may also be written in conservative form (‘conservative’ refers to conservation of neutrons; the RTE is similar in form to the neutron transport equation) by multiplying Z by d V dR, where dV = r2 sin 4 dr d+ d$ (differential element of volume) and dR = sin I3de do (differential element of solid angle), carrying out the differentiations, and dividing by dV dR. The coordinate derivatives (dr/d.s, de,/&, etc.) are also taken within the partial derivatives, giving the result2~5~6
az i - = - cos 0,: (Zr*) as r*
+sinre~io~Or
-& (Z sin 4) + sinr~i~~wr $ -&&(Zsin*B,)-s-&(Zsinw,).
(5b)
Conservative form is useful in numerical solutions, where it helps to regulate errors due to a finite mesh size. Equations (5a) and (5b) can be shown to be equivalent. This equivalence is a necessary condition for the conservative form, and serves as both a derivation tool and check for the new SDCSs presented in the following analysis. Direct numerical solution of the RTE may be achieved using the discrete ordinates method,2,4.6-* which is essentially a finite difference representation using a quadrature of points and integration weights for the directional mesh. Using a finite volume/finite solid angle cell approach,4q8 the conservative RTE may be integrated over the finite cell 6V6R to obtain the finite cell discrete Table 1. Nomenclature identifying the spatial/directional coordinate systems @DC.%). Abbreviation CylAXl+
I
Description spatially cylindrical, directional polar axin aliied
with axial spatial coordinate
W-d
rpatially cylindrical, directional polar axia aligned with radial spatial coordinate
cr-e
spatially cylindrical, directional polar axis aligned with angular spatial coordinate
SphRad +
rpatially 8phuica1, directional polar axiS aligned with radial spatial coordinate
SphPal
rpatially spherical, diiectional polar axb aligned with polar spatial coordinate
SphAsi
rpatially rpherical, ditional polar axb~ aligned with azimuthal spatial coordinate
tCommonly used coordinate system.
Coordinate systems for the RTE in curvilinear media
431
ordinates form of the pathlength derivative, written dy M1; dI’ d(l = LIj(Zrz)6,( -cos 4) ~,(~)cos(~,,)w, is +f sj(r2)6k(Zsin~)6,(~)sin(e,,)w,cos(w,.)w,
w,
+:sj(r2)Bk(~)6r(Z)sin(e,,)w,sin(o,.)w,
-:6j(r2)8k(-COS~)6r(~)6,(Z0,)w,
-~6j(r2)Bk(sin~)6,(~)sin(e,,)w,6,(ZR,),
(W
where the mesh indicies are rj, &, $,, 8,,, co,,. . The nomenclature Sj(Zr2) implies a point-centered cell difference over a single index, 4+,,2,k,,,m,nrj+,,2 - h_. ,,2,k,,,m,nrf_,,2, and w,,, and w, are integration weights for integrals over cells of directional angles,‘j viz. f(e) sin 8 de =f(&,)w,, 6. f(w) do =fbJ~n. WI s n s wn The suggestion of Abu-Shumaysg is followed, making the integration weight for the quadrature point (4 n), Wnl,“,equal to w,,,w,,, where w, and w, are independent. In Eq. (5c), the coefficients 0, and Q occupy places analogous to those of sin* 8, and sin o, in Eq. (5b). The reason for the substitution is common to all discrete simulations of the conservative RTE in curvilinear media, namely that making a direct discretization of Eq. (5b) results in an expression where al/as # 0 for constant Z, except when very restrictive conditions are placed on the spatial mesh and directional quadrature.‘v3 The common solution to this difficulty is to substitute the discrete coefficients 0, and R, for the analytical elements in the directional terms of Eq. (5b), where 0, and Q are derived specifically to allow the condition al/as = 0 for constant Z to be met for the conservative form. For the SphRad SDCS, the coefficients are defined by2v3 a,(@,) = 2w,c0se,,,, o~I+.~= 0, ~,(Q)=w,COS~,,m
(54
QL,=o=O,
Of)
where the initial values lie on particular directional cell boundaries, and Eqs. (se) and (5f) then provide recursion equations for 0, and Q on the remaining cell boundaries (m + l/2 or n f l/2). The commonly used SDCS in spatially cylindrical media has the directional polar axis aligned with the spatial axial coordinate [Fig. 2(a)]. This leads to the expression
az
az
as’
I
az
sin8,coso,--+-sinB,sino,-+cose,---sin8,sinw,---, ar r 84
(a)
az
i
ax
r
az aox
, lb)
Fig. 2. (a) Geometry of the CylAxi spatial/directional coordinate system. (b) Parallel disk geometry, which is well suited for non-symmetric application of the CylAxi spatial/directional coordinate system.
(64
432
PETERD. JONFSand YILDIZ BAYAZITOGLU
where the subscript x describes the polar axis alignment. 2*5,6 This expression has the particular utility that d6,/ds is uniquely zero, as the polar axis (spatial axial coordinate) does not change its orientation with changes in s. Therefore, only one directional boundary condition, Eq. (4a), is required. This spatially cylindrical, axial polar (CylAxi) SDCS has been applied to problems with the spatial geometry of a cylindrical annulus, Fig. 3(b). ‘OJ’In this geometry, the sweep of o, near a surface passes from a direction facing into the medium to a direction facing out, leading to a discontinuity in al/&, which is difficult to bridge numerically, and which makes the general boundary condition, Eq. (4a), very difficult to apply. This difficulty is avoided if the problem is spatially axially symmetric. In that case, there is a further directional boundary condition: al/&~, = 0 for o, such that ds lies in the plane of spatial axial symmetry. The directional boundary values may then be split between forward and backward facing half-ranges, and the discontinuity at the half-range boundaries is avoided. Alternatively, in non-symmetric cases, the values for the polar directional distribution can sometimes be used to generate a line of azimuthal boundary values in each half range, dependent on the polar angle, which may be used as substitutes for Eq. (4a). This approach was successfully used in Refs. 10 and 11 to solve the RTE in cylindrical annular media, using a SDCS (CylAxi) which would otherwise be ill-suited for this geometry. The CylAxi SDCS is best suited to a geometry of facing, parallel disks [Fig. 2(b)], where the directignal polar axis is normal to any surface, and the variation of Z with o, is everywhere continuous. Equation (6a) is written in conservative form in a manner analogous to that for spatially spherical geometry, using d V = r dr d4 dx and dR = sin 9, do, dw,, and resulting in
az i = - sin 0, cos w, g (Zr) + i sin 8, sin 0, z + cos 8, g - b sin 0, & as r 84
-
x
(Isin 0,).
(6b)
Equation (6b) satisfies the necessary condition of equivalence to simple form, Eq. (6a). Equation (6b) may be further adapted to a finite volume/finite solid angle cell discrete ordinates form in a manner similar to that for Eq. (5~) with the result
,& ss 6y
dV dQ= 6jUr)h(4) ~I(x)sin(e,,)w,cos(~,.)w, +
dj(r)
b(Z) 6Win(kJw,
W+,k, (6~)
where
(b)
Fig. 3. (a) Geometry of the CylRad spatial/directional coordinate system. (b) Cylindrical annuius geometry, which is well suited for non-symmetric application of the CylRad spatial/directional coordinate system.
433
Coordinate systems for the RTE in curvilinear media
The foregoing introductory discussion has outlined the occurrence of directional derivatives in the RTE in curvilinear media, presented general boundary conditions for these terms, shown the form of the pathlength derivative in the SDCSs most commonly employed for spatially spherical and cylindrical geometries, with finite cell discrete ordinates expressions, and discussed the boundary condition and continuity difficulties associated with using these SDCSs in ill-suited spatial geometries. For practical solution of the RTE in non-symmetric curvilinear media, it is often advantageous to use an SDCS suited to the geometry (with the directional polar axis normal to any surface) in order to apply the general boundary conditions of Eqs. (4). However, to the knowledge of the authors, the only SDCSs for which the pathlength derivative in the RTE has been derived are what have been termed the ‘usual’ SDCSs, SphRad and CylAxi. In the following analysis, alternative SDCSs are proposed which are suited for particular spherical and cylindrical spatial geometries. Pathlength derivative expressions are derived for these SDCSs in simple, conservative, and finite cell discrete ordinates form. ANALYSIS Spatially cylindrical, radial polar
An SDCS more suitable for solution of the RTE in cylindrical annular geometries than the commonly used CylAxi is the spatially cylindrical, radial polar SDCS [CylRad; Fig. 3(a)], which is based on a directional polar axis aligned with the spatial radial direction. With the CylRad SDCS, the directional equitorial plane always coincides with the surfaces of a cylindrical annulus, allowing the boundary conditions of Eqs. (4) to be applied. CylRad allows boundary values to be determined for 8, = 0, II using Eqs. (4b, c), while the distribution of Z with o, is uninterrupted and Eq. (4a) may be applied. The CylRad SDCS may be derived as a transformation from CylAxi by expressing the direction of the differential pathlength ds for each SDCS in the local spatial coordinate system (i, &a) as follows: ds -= sin 8, cos o,i + sin 0, sin w,f$ + cos 8,?, IdsI (74 ds -= cos 8, i - sin 8, sin W,fj + sin 8, cos CO, i, IdsI where the subscript x denotes polar (0) and azimuthal (0) directional angles based on a polar axis aligned with the local spatial axial component (CylAxi), while the subscript r denotes a polar axis aligned with the local spatial radial component (CylRad). Equating the s derivatives of the spatial component terms in Eqs. (7a) and using the CylAxi definitions
(7’4 from Eq. (6a), the CylRad directional derivatives may be derived, resulting in
de,_ ds - -;
do, 1 sin 8, sin’ w, , -=
1 --cos~,sino,coso,. r
ds
Using the directional derivatives and spatial components, the simple form of the CylRad pathlength derivative may be assembled similarly to Eq. (2b), resulting in
az
az i ar r
-=cos&---sin&sinw,-+sinf+cos~-
as
az
az
w
ax 1 --
r
sin
az i 8, sin’ 0, -ae, - -r cos 8, sin 0, cos W,-;;
7
.
(74
Following the methods demonstrated for CylAxi, the conservative form for CylRad is found to be 1 az i a -=-cos&--(Zr)--sintI,sino,-+sintJ,coso,-
as
r
ar
r
az w -S&(Z
az ax sin2 e,) - i cos 8, -& (I sin 0, cos o,), I
(7e)
PETERD.
434
JONFS and YILDIZ BAYAZIT~CLU
where the necessary equivalence between Eqs. (7d) and (7e) is satisfied. Analytically, either Eq. (7d) or (7e) is correct; however, Eq. (7e) is generally more stable in numerical simulations. Again following the methods demonstrated for CylAxi, the finite cell discrete ordinates form of the pathlength derivative for CylRad is (7f)
6v d~~dvda=6,(Z~)~,(~)6,(x)cos(e~,~)W~W~ ss -
6j(rPk(1)
Uxbin(%Jw,
sMw..N,
+~6j(r2)bk(~)8,(l)sin(B,,)w,cos(o,.)w,
-
6j(r)
4A4)
4(X)
4S@,)sin2b4,.h
-
6j(r)
6k(+)
G,(x)cos(e,~)w~6~(zn~)~
where in order to satisfy the necessary condition that al/as = 0 for constant Z, we choose d,(@,)
=
Wrn
cos(eLm)@lo,=ll= 0,
sin2(w, ,) ’
W4) = 0, Q,.,
l/2 =
4
sin(~&os(qn).
Ug) V’h)
Equations (7g) are the recursion and intitial value for 0, in the manner of Eqs. (5e), (5f), and (6d). Equations (7h) indicate that, in order to satisfy the necessary condition, Q is not allowed to vary within the finite cell; however, R, is discontinuous at the cell boundaries n &-l/2, and must vary from cell to cell in order to maintain the correspondence between Eqs. (7e) and (7f). Note that the sin* (w,.) term in the fourth line of Eq. (7f) is canceled by the choice of 6,(0,), and that the combination of Eqs. (7f) and (7g) does not contain the sin2(o,J term. In order to use general directional boundary conditions, without recourse to the spatial symmetry of a medium, it is necessary to use an SDCS with a directional polar axis normal to the boundary surfaces. The CylRad SDCS has a directional polar axis which is normal to surfaces of constant radius, e.g. cylindrical tubes. Hence, the CylRad SDCS is appropriate for solution of the RTE in cylindrical annular media [Fig. 3(b)]. Spatially cylindrical, angular polar
The spatially cylindrical coordinate system may be used to define other geometries besides a cylindrical annulus. For instance, a medium bounded by an array of rectangular plates set at angles to each other and hinged at a common point [cylindrical wedge, Fig. 4(b)] is readily described in
(b)
Fig. 4. (a) Geometry of the CylAng spatial/directional coordinate system. (b) Cylindrical wedge geometry, which is well suited for non-symmetric application of the CylAng spatial/directional coordinate system.
435
Coordinate systems for the RTE in curvilinear media
cylindrical coordinates. For such a geometry, the boundary surfaces are defined by constant angular coordinates, and the appropriate SDCS has a directional polar axis aligned with the spatial angular direction [Fig. 4(a)]. The pathlength derivative for this spatially cylindrical, angular polar (CylAng) SDCS is derived using ds -= IdsI
sin
8, cos o4 i +
cose,$ - sin 8, sin o&
with the result expressed in simple form
ar
CYI
I
sin 8~~0s~ - +-case @ar r as=
ar
az
- -sin 8, sino "ax 4a4
1 az +; cose,$cosw~~-
sin cf+ -- 1 r sin 8,
@b)
and in conservative form
ar -=as
1 az a az sin e,+cos0,--(Zr)+- cos e - -sin 8, sin 0 r #ax & r $84 I
sin e In finite cell discrete ordinates form, the CylAng pathlength derivative dy ,,!dV Is
da = 6jW) b(4)
4(~)W&,h,
sin o++,). (8~)
is expressed
cos(oo,,h
+~~(~)~~(Z)GI(X)COS(~~,,)W,W, -:sj(r2)8~(~)6,(l)sin(e~,,)w,sin(o~,.)w, +6j(r)6,(~)6,(x)G,(lO,)cos(o,.)w,
In order to satisfy the necessary condition, we choose &(@,) = - w, sin(&,,), WV
= 0,
a,,,
1/2 =
w,
@61s,=n,2= 0,
(8e)
W+,.),
WI
where Eqs. (8e) are a recursion and initial value set and Eqs. (8f) define constant cell values for R,, as described in the previous section. In addition to CylAxi, CylRad, and CylAng, further SDCSs may be developed for spatially cylindrical geometries. These may be derived by following the techniques demonstrated above, as long as the directional polar axis is kept normal to the medium’s boundary surfaces. For instance, a cylindrical spiral surface may be defined by ds -= IdsI
cos e,(cosait+ sin ~(4) + sin 8,sin o,( -sin c& + cos a$) + sin 8, cos 0,i
(9)
where a is the angle between the spatial axial coordinate and the normal to the spiral surface. Using Eq. (9) and following the steps outlined above, a pathlength derivative expression may be derived for a medium bounded by this cylindrical spiral surface. Spatially spherical polar polar
The spatial spherical coordinate system may be applied to other geometries besides spherical annuli. For these other geometries, it is necessary to use a directional coordinate system with a polar axis normal to the boundary surfaces. For example, if the boundary surfaces in a spherical coordinate system are defined as surfaces of constant spatial polar angle, over a range of radius and spatial azimuthal angle, then a geometry of the annulus between cones of different apex angle
436
PETERD. JONESand YILDIZ BAYAZIT~GLU
Fig. 5. (a) Geometry of the SphPol spatial/directional coordinate system. (b) Conical annulus geometry, which is well suited for non-symmetric application of the SphPol spatial/directional coordinate system. is described for solution The SphRad system,
[Fig. S(b)]. In order to use this SDCS, spatially spherical polar polar [SphPol, Fig. 5(a)], of the RTE, it is necessary to derive a new expression for the pathlength derivative. and SphPol pathlength directions, relative to the (i, $, $) local spatial coordinate
are
sin 8, sin 0, sink II/,
ds
1
ds
^ sin + cos e,4 -
(ds=cosB~i+sin0,cosw,$ ldsl = sin 8, cos 0,:
+
e. cos wg
-
*9
sin 4
which, when combined with the SphRad directional derivatives from Eq. (Sa), sin 8, sin w, r tan4 ’
dw, -=ds
dB, _ 1 . sin%,, ds - -;
(lob)
result in the SphPol directional derivatives de, _= ds
-- 1 cos 8, cos wg r
(1Oc)
Using Eqs. (lOa) and (lOc), the simple form of the pathlength assembled, viz.
az
az
- =sin04c0so as
-
+ar
i
+-00se r
az
+a4
-
derivative in SphPol may be
sin 8, sin o4 az rsin4 ZJ
sin 8, sin’ o+ --cos 8, cos We az do. (104 tan I$ tan $J > > 4 The conservative form may be determined by multiplying Z by dV = r2 sin 4 dr d4 d$, dR = sin 8, de, dw,, and the directional derivatives, differentiating, and dividing by d V dR, with the result
az
I
- =-sin0#coso as r2 COSW+ a + ~-(zsin86c0se,)r sin 8, ae, __qz. 1 r sin e+ am,
z(Zrz)+
2
-.$ (1 sin 4) _ sinre~i~~w,
$
4 ar
sin2 w+
r
sin o,+) -fg$&(’
+z sin e4 tan C#Jae,
sin’ ( 8,)
sin o4 cos 0~~).
UW
437
Coordinate systems for the RTE in curvilinear media
Generation of the finite cell discrete ordinates form for SphPol is straightforward, complicated by the multiple directional partial derivatives in Eq. (10e). The result is
ss ar
bv mz
dV da = dj(Zr*> U-~0s 4) 4($>sW,,,)w, +
fsj
43
though
cos(w&w,
sin4) 4(ll/)COs(~,,)w,w,
- f sj(r2) U+) UOsin(~~,,)w,sin(q,.)w, + i6jCr2) sk(-cos4) 6!($)~~(Z@I~)COS(WI$,~)W~ - 46j(r2) Usin 4) h(JI>hAZ@24)sin2(w~,,)w, -
i6jCr2)
wm
d,(IR,,) 6k(-cos4) a,(+> SW, m ) (IOf)
- i 6j(r2) Msin 4) 4(II/)cos(B$,,)w, 6,(ZR,,), where we choose
4d%J = -w, sin(e+J, U@*,)
= w,
MT+)
= 0,
4 (% ) = 0,
cos(4, m) sin2(04,, ) ’
o,,I,,=,, @*&dl=
= 0,
wk)
0,
(1W (1Oi)
fi~,~+ ,,2= w, sin(04,.), Q24,n+ 1,2= w, side,
, )cos(w+, .).
(W)
Spatially spherical, azimuthal polar
Construction of an SDCS in spatially spherical media where the directional polar axis is oriented along the direction of the local spatial azimuthal coordinate results in the spatially spherical, azimuthal polar (SphAzi) SDCS, Fig. 6(a), which may be applied to a spherical wedge geometry, Fig. 6(b). The SphAzi pathlength direction, relative to the local directional coordinate system V, $, ~0, is ~0~8, ds -=sin8ticoswtii+sin8+sino+$ - ---$. (114 sin 4 IdsI The SphAzi directional derivatives may be solved using Eqs. (5a), (lob), and (llc),
leading to
(11’4 COS~,COS
we
sin w+
which results in expressions for the pathlength derivative in simple form,
az cos 8, az az i az -_= sin 8, cos w - + - sin 8,sin ws - - a4 r sin4 a$ +iYr r as
COS~+COS W+
sin
W$
tan 8,tan f#~ and, in conservative form, az i as =7sin8,coswti~(IrZ)+
sinre~i~~w’ $
(I sin 4) + z
sin wti 1 + COSW# -&(I sin e,cose,) r sin 8,( tan f#~ ) 1 COSTS sin We). +z +z cosw,)-+ r sin 8, aw, r tan 8,tan 4 aa,
$
+-
(114
PETER D. JONES and YILDIZ BAYAZIT~GLU
438
(b)
Fig. 6. (a) Geometry of the SphAzi spatial/directional coordinate system. (b) Spherical wedge geometry, which is well suited for non-symmetric application of the SphAzi spatial/directional coordinate system.
Generation
development
of the finite cell discrete ordinates form of the pathlength derivative is similar to the for other SDCSs, with the result
az 6v 6n~dVdR=S,(z~2)b~(-cos~)g,(~)sin(8l,~)w~cos(w~,~)w, ss
+ + 6j(r2) +f +
&(Z sin 4) 8,(ll/)sin(O,,,)w,
sin(o+)w,
~j(r2)~,(~)~(~)cos(~,,,)w,w,
i6j(r2)k%(sin
+ i 6j(r2) 4&n - f
6j(r2)
4)sin(o,,.)
4) 4(+)
+
[
6,(-cos
1 sin(8e,,) - sin(%.,)
fik(-COs4) d,($)
1
4)cos(o,,.)]
S,($) &(ZO,)w,
MQ,)
wm6,(ZR,,),
(1W
sln(%, 1 m
where we choose 6,(0,)
[
= W, c?y:Vm)) - sin(&,) .m
MA,)
= - w,sin(f+,.),
hW2,)
= w,cos(q,),
= 1, 0,IBe
QI,I,,=,!2 = 0, ~,,I,,=,=@
x,2= 03
(llf> (1w (1 lh)
Each of Eqs. (1 lf, g, h) is a recursion and initial value set. In order to use the SphAzi SDCS with the directional polar axis normal to the medium’s boundary surfaces, it should be applied to a
Coordinate systems for the RTE in curvilinear media
439
Table 2. Summary of spatial/directional coordinate systems (SACS) results for the SLXS defining geometry; an example of an application; and the pathlength derivative in simple analytic form, conservative analytic form, and finite cell discrete ordinates form.
Result
SphRad
BphPdlspMrl
EF--ry
Fig.2a
Fig.30
Fig.4a
Fig.la
Fig.Sa
Fig.Ba
application
Fig.2b
Fig.3b
Fig.lb
Fig.lb
Fig.Sb
Fig.6b
(N/ad),
simple
%.@a)
JW’ld)
Eg.(8b)
E&a)
Eq.(lOd) Eq.(llc)
(N/ad),
conae-rvative
EOb)
Eq.(‘lc)
Eq.(8c)
E&b)
Eq.(lOe)
EgW
EgOf)
Eg.(8d)
JQ.(Se)
Eq.(lOf) Eq.(llc)
/,v &,, (aI/ad)mdn, discrete ordinates
Eq.(lld)
spatial geometry such as plane surfaces set at an arbitrary angle, much as for CylAng, but with the spatial location defined in spherical rather than cylindrical coordinates [Fig. 6(b)]. The SphAzi SDCS was investigated in a preliminary manner by Jones,3 although the form of the pathlength derivative derived in that work cannot be developed into a conservative form which satisfies the necessary condition that al/as = 0 for constant I. This condition is met in the earlier work only if Z is multiplied by dV rather than dV dR. SUMMARY The radiative transfer equation in curvilinear coordinates requires directional boundary conditions for its solution. For particular problems, it may be possible to derive directional boundary conditions which result from geometrical symmetry. However, in general cases without useful symmetry, the only available boundary conditions are Eq. (4a) and the boundary values which may be developed using Eqs. (4b) and (4c). In order to apply these general boundary conditions, the coordinate system which defines the direction of the radiation intensity must be set up in such a way that its polar axis is normal to the spatial boundaries, especially when in proximity to the boundary surfaces. Therefore, for a particular local orientation of the directional coordinate axes, only very specific spatial boundaries may be defined for use with the general directional boundary conditions. In the commonly used spatial/directional coordinate systems, widely presented in the literature and here termed CylAxi (spatially cylindrical, directional polar axis aligned with the axial spatial coordinate) and SphRad (spatially spherical, directional polar axis aligned with the radial spatial coordinate), the general directional boundary conditions may only be applied in spatial geometries of facing parallel disks (CylAxi) and spherical annuli (SphRad). In the present work, spatial/directional coordinate systems are derived in which the general directional boundary conditions may be applied (and the radiative transfer equation solved) for other spatial geometries. Intensity pathlengths derivatives, which contain the directional derivatives, are derived for each new spatial/directional coordinate system in simple, conservative, and finite cell discrete ordinates form. These results are summarized in Table 2. Spatial/directional coordinate systems may also be developed for spatial geometries defined by a combination of spatial coordinates, as opposed to the geometries presented here, which are defined by a single coordinate held constant. Using the techniques outlined in the present work, intensity pathlength derivatives may be derived for such coordinate systems. The radiative transfer equation may then be solved in these geometries using the general directional boundary conditions, unaided by boundary conditions derived from geometrical symmetry. REFERENCES 1. M. N. &igik, Radiative Transfer and Interactions with Conduction and Convection, Wiley, New York, NY (1973) [reprint: Werbel & Peck, New York, NY (1985)]. 2. E. E. Lewis and W. F. Miller Jr., Computational Methods of Neutron Transport, Wiley, New York, NY (1984). QSRT
48/4-G
PETERD.
440
3. P. D. Jones, “Radiation
4. 5. 6.
7. 8. 9.
10. 11.
JONES and YILDIZ BAYAZITOGLU
and Convection Heat Transfer in Particle-Laden Fluid Flow,” Ph.D. Thesis, Rice University (1990). J. R. Tsai, M. N. ozivik, and F. Santarelli, JQSRT 42, 187 (1989). R. Viskanta and M. P. Mengiic, Prog. Energy Combust. Sci. 13, 97 (1987). B. G. Carlson and K. D. Lathrop, in Computing Methods in Reactor Physics, H. Greenspan, C. N. Kelber, and D. Okrent eds., Gordon & Breach, New York, NY (1968). W. A. Fiveland, J. Heat Transfer 106, 699 (1984). W. A. Fiveland, J. Thermophys. Heat Transfer 2, 309 (1988). I. K. Abu-Shumays, Nucl. Sci. Engng 64, 299 (1977). A. Yucel and M. L. Williams, in ASME Proceedings of the National Heat Transfer Conference, HTD-Vol. 96, 1, H. R. Jacobs ed., American Society of Mechanical Engineers, New York, NY (1988). A. S. Jamalludin and P. J. Smith, in ASME Proceedings of the National Heat Transfer Conference, HTD-Vol. %, 1, H. R. Jacobs ed., American Society of Mechanical Engineers, New York, NY (1988).