Cylindrical Photocatalytic Reactors. Radiation Absorption and

3094

Ind. Eng. Chem. Res. 1997, 36, 3094-3109

Cylindrical Photocatalytic Reactors. Radiation Absorption and Scattering Effects Produced by Suspended Fine Particles in an Annular Space Roberto L. Romero,† Orlando M. Alfano,‡ and Alberto E. Cassano*,‡ INTEC,§ Universidad Nacional del Litoral (U.N.L.) and CONICET, Gu¨ emes 3450, 3000 Santa Fe, Argentina

The radiation field for an annular-type photocatalytic reactor has been modeled, and the resulting integro-differential equation has been solved resorting to the discrete ordinate method. The spatial domain was described in terms of two cylindrical coordinates and the field of direction with two angular variables. The radiation entering through the inner reactor wall was obtained from the extended source with voluminal emission model. The system parameters for the radiative-transfer equation were taken from experimental measurements. Resorting to computational experiments, the effects produced by different concentrations of (1) transparent scattering centers and (2) radiation absorbing (catalytic) particles were investigated. In the second case, Aldrich and Degussa P 25 titanium dioxide particles were used. The following information is reported: (1) specific intensities as a function of position inside the reactor and of the angular distribution of directions; (2) incident radiation profiles and local volumetric rate of energy absorption profiles, both as a function of the spatial position inside the reactor. Introduction By the end of the 1970s and early 1980s, photocatalytic reactions strongly called the attention of the scientific community when they were proposed as a suitable means to promote the separation of hydrogen and oxygen from water using solar illumination. Almost at the same time when, in the middle 1980s, the practical feasibility of this process was severely questioned, a second attractive application was envisaged. It was made evident after the publication of several papers by Ollis’ group (Pruden and Ollis, 1983; Hsiao et al., 1983; Ollis et al., 1984) showing that an equivalent approach could be used for degrading organic compounds. Consequently, one of the concepts that has received an increasing degree of attention has been the use of titanium dioxide as a catalyst for the lightinduced photolysis of pollutants (Ollis and Al-Ekabi, 1993). Catalysis by illuminated titanium dioxide is the result of the interaction of the electrons and holes generated in the photoactivated semiconductor with the surrounding medium; thus, as a consequence of light absorption, electron-hole pairs are formed in the solid particle that can recombine or participate in reduction and oxidation reactions that have as a result the decomposition of contaminants. The development of a reliable knowledge base is still in its initial stages in which problems related to catalyst preparation and immobilization, catalyst chemical and mechanical stability, chemistry and kinetic networks of pollutant degradation, intrinsic reaction kinetics including the effects of irradiance level and substrate competition, photocatalytic reactor design, and process integration with other water treatment technologies represent some of the main weaknesses (National Research Council, 1991). * To whom correspondence should be addressed. Fax: +5442-55 0944. E-mail: [email protected]. † CONICET Research Staff. ‡ Professor at U.N.L. and CONICET Research Staff. § Instituto de Desarrollo Tecnolo ´ gico para la Industria Quı´mica. S0888-5885(96)00664-1 CCC: $14.00

Contrasting with the profuse literature regarding the chemistry of these reactions, a parallel research effort has not been applied to the engineering aspects of these processes. Design methods are very scarce, and in most of the cases, they rely on empirical or semiempirical approximations (Jakob et al., 1993). Polychromatic radiation transport in participating media (with absorption, scattering, and chemical reaction) is properly represented by a set (in theory, an infinite set) of integro-differential equations (Spadoni et al., 1978; Cassano et al., 1995). Only for very simplified and somewhat unrealistic situations is it possible to get analytical solutions. Special numerical techniques are necessary that can be adapted from general transport theory and particularly from neutron transport applications (Carlson, 1953, 1963; Duderstadt and Martin, 1979). Very recently, some rigorous approaches based on the fundamentals of radiation transport theory have been proposed, but they have been restricted to small-scale reactors (Cabrera et al., 1994; Alfano et al., 1995; Brandi et al., 1996) or mostly to theoretical computations (Pasquali et al., 1996). Among the most attractive configurations for photocatalytic applications, the cylindrical reactor has often been used in solar through collectors (Turchi and Mehos, 1992; Pacheco et al., 1993). If artificial light is to be used in combination with continuous flow reactors, the cylindrical geometry of an annular reactor is perhaps the simplest and most efficient configuration. Cylindrical geometries pose a very challenging problem in that radiation transport in curvilinear coordinate systems gives rise to additional complications in the numerical procedures. These complications are even magnified by the presence of an artificial radiation source that translates into an extremely concentrated boundary condition. As will be seen in the next sections, the solution techniques for these problems are very unconventional as compared with a large variety of chemical engineering applications. The detailed description of the radiation-field modeling in an annular photocatalytic reactor will be presented in this paper. © 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3095 Table 1. Adopted Reactor and Lamp Characteristics for the Theoretical Model system annular reactor Tempax glass walls lamp Hanovia LL-189a-10/1200

parameter

value

rin rou LR nominal input power arc length (LL) lamp radius (rL) characteristic emission photochemical power between 295 and 405 nm total emission range

3 cm 5 cm 30.48 cm 1200 W 30.48 cm 1.03 cm isotropic and voluminal 6.5 × 10-4 einstein s-1 220-1370 nm

Table 2. Characteristics of Solid Particles material A-Ba

Porasil latexb Aldrich TiO2 Degussa TiO2 a

nominal diameter, nm 37 × 2 × 103 150-200 30-90

103-75

×

〈σλ*〉Pλ, cm2 g-1

〈κλ*〉Pλ, cm2 g-1

CPm, g cm-3 × 10-3

265 24 464 35 980 53 502

0 0 2758 5316

3-6 0.1-1.0 0.05-0.5 0.05-0.5

103

Water Associates (1970). b Poly(vinyltoluene). Note: σ ) σ*CPm; κ ) κ*CPm.

Figure 1. Annular photoreactor.

Reactor Geometry and Model Parameters Computational experiments will be carried out with data corresponding to the preliminary design of a small pilot-plant reactor. Thus, theoretical predictions will be obtained employing realistic values of a practicable system. The reacting system is a cylinder of annular cross-sectional area. This annular space surrounds a cylindrical tubular lamp (the radiation source) having, in this case, the same length as the reactor (Figure 1). The lamp is a polychromatic source of the arc type having voluminal emission (Cassano et al., 1995). Emission by the lamp is mainly in the ultraviolet and visible ranges, although some near-infrared radiation is also produced. The lamp is refrigerated by water circulating in a cooling jacket that is made of transparent quartz; this system absorbs most of the unnecessary infrared radiation normally produced by the lamp. The inner reactor wall is made of Tempax glass that cuts off most of the emitted radiation below 295 nm; however, this wall is transparent to radiation in the wavelength range of interest (295 e λ [)] nm e 405). Above 405 nm, there is no absorption by titanium dioxide since its band-gap potential is never smaller than 3.0 eV. The gap between both rector walls must be necessarily small (no larger than 2-5 cm) because, under the normally employed catalyst concentrations, absorption by titanium dioxide is very strong and, consequently, after a few centimeters the reactor is opaque to radiation. The annular space is filled with a spatially uniform suspension of solid particles in water. The fluid phase

is assumed to be transparent to the employed radiation in the above-mentioned wavelength range; hence, absorption may be the result of interaction of photons with the solid only. Different types of particles will be investigated: (i) transparent particles made of rather large-size silica beads (37 e dp [)] µm e 75), (ii) transparent particles made of a suspension of smallsize latex spherical beads (dp = 2 µm), (iii) Aldrich titanium dioxide (99% Anatase) absorbing catalytic particles (150 e dp [)] nm e 200) and Degussa P 25 titanium dioxide absorbing catalytic particles (30 e dp [)] nm e 90). Concentrations of the solid suspensions vary according to the respective optical properties. In order to observe appreciable effects with rather large particle sizes and no absorption, the solid concentration must be relatively large [(3 - 6) × 10-3 g cm-3], reaching in this way significant values of the optical density. With highly absorbing solids, the concentration is usually small [(0.05 - 0.5) × 10-3 g cm-3]. At almost neutral pH, titanium dioxide tends to agglomerate (Martı´n et al., 1993). Depending upon the solid characteristics and the hydrodynamics of the surrounding fluid (stirring, recycling, ultrasonic mixing, etc.), this agglomeration may be more or less significant. Cabrera et al. (1996) reported from a wide variety of experiments diameters that range from 300 to 700 nm. This is an important consideration for the radiationfield model as will be seen below. Other details can be found in Tables 1 and 2. Radiation Field The theoretical grounds for radiation transport in heterogeneous participating media are fairly wellknown (Ozisik, 1973; Siegel and Howell, 1992). An extension to heterogeneous, participating and reacting media was first proposed in a theoretical paper presented by Spadoni et al. in 1978. However, up to just recently, most of the system parameters that are required for a realistic application of these theoretical approaches to water suspensions were unknown. A systematic development, including for the first time applications with experimental verifications of this rigorous theory, was recently published by Alfano et al. (1994, 1995) and Cabrera et al. (1994, 1995, 1996). Cartesian geometry was always used. Essentially one must write, for monochromatic light, a radiation balance along a given direction of propagation together with the appropriate boundary conditions

3096 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

(Cassano et al., 1995). This balance is generally made with the following assumptions: (1) scattering is multiple but independent (the last is a correct hypothesis for the fairly dilute suspensions usually employed in photocatalysis), and (2) scattering is elastic, which means that after scattering, only changes in direction occur; i.e., it excludes changes in energy (frequency) and is also valid for photocatalytic systems. Under these conditions, along its traveling trajectory through a heterogeneous participating medium, a ray may lose radiation due to (i) absorption and (ii) out-scattering and may gain energy by (iii) internal emission and (iv) inscattering (coming from multiple scattered rays in other points of the reaction space). In photocatalytic processes, which are normally carried out at almost ambient temperatures, internal emission is always neglected. Absorption, out-scattering, and in-scattering are modeled with the appropriate constitutive equations. A minimum of two parameters and a distribution function for scattering are needed. They are the linear absorption coefficient (κλ), the linear scattering coefficient (σλ), and the phase function for scattering [p(Ω′ f Ω)]. The resulting radiative-transfer equation (RTE) is

In the elementary volume of radiation absorption, for a single-photon absorption, energy is absorbed according to

eaλ (x,t) ) κλ(x,t)Gλ(x,t)

(5)

eaλ is the spectral (monochromatic) local volumetric rate of energy absorption (LVREA) (Irazoqui et al., 1973) or the spectral-absorbed incident radiation. Once the radiation absorption distribution is known for any position in the three-dimensional space (x) and any time (t), an integration of the monochromatic LVREA over the frequency interval is required to account for polychromatic emission:

ea(x,t) )

∫λλ

max

min

κλ(x,t) dλ

∫Ω)4πIλ(x,Ω,t) dΩ

(6)

In the absence of emission, for a cylindrical geometry and considering that x ) x(r,z,β) and Ω ) Ω(θ,χ), the RTE takes the following form (Lee, 1962; Duderstadt and Martin, 1979):

[

∂Iλ(r,z,β,θ,χ,t) sin χ ∂Iλ(r,z,β,θ,χ,t) + ∂r r ∂χ sin χ ∂Iλ(r,z,β,θ,χ,t) ∂Iλ(r,z,β,θ,χ,t) + + cos θ r ∂β ∂z [κλ(r,z,β,t) + σλ(r,z,β,t)]Iλ(r,z,β,θ,χ,t) ) σλ(r,z,β,t) p(Ω′fΩ)Iλ(r,z,β,Ω′,t) dΩ′ (7) Ω′)4π 4π

sin θ cos χ

]

∫

The spectral specific intensity gives the amount of energy per unit normal area, unit solid angle (Ω) about the direction of propagation (Ω), unit time, and unit wavelength (frequency) interval. In this equation, the phase function [p(Ω′fΩ)] is normalized according to

∫Ω′)4πp(Ω′fΩ) dΩ′ ) 1

1 4π

Iλ(r,LL,β,θ,χ,t) ) Iλ(r,0,β,θ,χ,t) ) 0 (2)

Equation 1 must be solved with the appropriate boundary conditions at the bounding surfaces of the reacting space; i.e., (1) top and bottom of the reactor, (2) inner wall of the reactor, and (3) outer wall of the reactor (Figure 1). For each direction of radiation propagation, the bundles of radiation rays have an associated given solid angle. Since at every material point in the reaction space we are interested in the incoming radiation from all physically permitted directions (actually, solid angles), an integration over all the arriving rays will be required. This gives rise to the incident radiation concept, defined by

Gλ(x,t) )

∫Ω)4πIλ(x,Ω,t) dΩ

(3)

In eq 3, an integration over all possible directions Ω(θ,φ) over the entire spherical space has been performed. For polychromatic radiation, an integration over the wavelength (frequency) range of interest must be performed, accounting for the overlapping wavelength regions of lamp emission, reactor wall transmission (absorption and reflection), and radiation absorbing reacting (or catalytic) species absorption coefficient:

G(x,t) )

∫λ

λmax min

dλ

∫Ω)4πIλ(x,Ω,t) dΩ

The boundary conditions in this cylindrical geometry can be stated as follows: (1) At the reactor top and bottom circular (annular) spaces, there are no radiation entrances:

(4)

(8)

(2) At the inner wall of the reactor, the incoming radiation is that produced by the lamp placed at the reactor center line:

Iλ(rin,z,β,θ,χ,t) ) f (z,β,θ,χ,t)

(9)

An emission model for the lamp should provide this information. In the case of the annular reactor, when using lamps made of a transparent wall, this boundary condition could also include contributions from the opposite reaction space back-scattering. In a welldesigned annular reactor, when the lamp surface is opaque (that is, the case of lamps with superficial and, hence, diffuse emission), contributions due to backscattering will be negligible. An expression for f(z,β,θ,χ,t) will be developed farther ahead. (3) At the reactor outer wall, we could have a partially transmitting-partially reflecting boundary:

Iλ(rou,z,β,θref,χref,t) ) F′λ(θ,χ,n)Iλ(rou,z,β,θ,χ,t) (10) In the case that all incoming radiation (boundary condition 9) is absorbedscertainly a good design criterionseq 10 can be further simplified. In this work, at the external radius, we will assume null reflection. Boundary Condition and System Parameters Emission Model for the Lamp. The three-dimensional source with the voluminal isotropic emission

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3097

They are the analytical expressions of the planar angles that, measured from that point, define the limiting rays that are tangent to the boundaries of the lamp volume (top, bottom, and lateral cylinder). From Irazoqui et al. (1973), we find the following:

limits for the independent variable θ θi(φ) ) tan-1

{

}

rin cos φ - [rin2(cos2 φ - 1) + rL2]1/2 Xi - z (17)

in eq 17, i ) 1 w Xi ) LL and i ) 2 w Xi ) 0 Figure 2. Geometry for the emission model (adapted from Irazoqui et al. (1973)).

model will be used (Irazoqui et al., 1973). This model was developed for tubular, arc-type UV lamps (Figure 2). Considering a spherical coordinate system to describe the ray trajectory, that is, Ω ) Ω(θ,φ) and applying the RTE to the lamp emission, the following equation results (Cassano et al., 1995):

Ioλ(θ,φ) ) jeλ∆FL(x,θ,φ)YR,λ(θ,φ)

(11)

∆FL(x,θ,φ) ) F2(x,θ,φ) - F1(x,θ,φ)

(12)

with

It must be remarked that ∆FL is a function of the position x at the reactor inner wall and the direction of the incoming radiation represented by a system of spherical coordinates (θ,φ). Reflection and absorption at the reactor wall are accounted for by means of the wall transmission coefficient. The value of jeλ must be related to the lamp output power Pλ,L. By definition, jeλ is the energy emitted per unit volume, unit solid angle of emission, and unit time; then,

jeλ )

Pλ,L 4π2rL2LL

[

]

(rin2 - rL2)1/2 rin

(18)

Absorption and Scattering Coefficients and the Phase Function. Equation 7 shows that two parameters (the absorption and the scattering coefficients) and the phase function are required. Although the phase function is normally approximated with a wavelengthindependent representation, the optical parameters have a known dependence with wavelength (Cabrera et al., 1996). Rigorously speaking, one should solve the RTE for each emission line of the lamp and afterwards add the results. A second approximation is to use wavelength-averaged properties. The required computer-processing time can be greatly reduced if wavelength-averaged properties are used (a sort-of pseudomonochromatic approximation). Extension to a real polychromatic solution does not involve additional theoretical complexities (Claria´ et al., 1988). The pseudomonochromatic approximation (Cassano et al., 1995) is obtained using the lamp output power as a weighting factor: λn)405

rin cos φ ( [rin2(cos2 φ - 1) + rL2]1/2 (14) F1,2 ) sin θ and

2[rin2(cos2 φ - 1) + rL2]1/2 (15) sin θ

Finally, the boundary condition results:

Ioλ(θ,φ) )

-φ1 ) φ2 ) cos-1

∑

(13)

From eqs 11 and 12, we must know the value of ∆FL:

∆FL ) F2 - F1 )

limits for the independent variable φ

[rin2(cos2 φ - 1) + rL2]1/2 YR,λ(θ,φ) sin θ 2π2rL2LL (16) Pλ,L

Angular Limits for Calculating the Boundary Condition (Irazoqui et al., 1973; Cassano et al., 1995 (Figure 2)). Equation 16 indicates that at every point x, contributions from different directions (θ,φ) must be known. The limiting values of the θ,φ system are defined by (1) the lamp geometry, (2) the reactor geometry, and (3) the relative positions of the lamp and the reactor. In our case, these limits must be known at every point on the surface of the inner reactor wall.

〈Ψ*λ〉Pλ )

λ1)295

Ψλi*Pλi

λn)405

∑

λ1)295

(19)

P λi

This approach was applied to every optical property employed in this work (κ* and σ*). This procedure does not introduce significant distortion in the results because the obtained average is a reasonable representation of the actual distribution in wavelengths, particularly for computational experiments. The exact approach should have been used if the results from this simulation were to be compared with experimental data. The scattering coefficient for silica beds (Porasil A-B) was calculated from measurements performed by Alfano et al. (1995). These particles are transparent in the interested wavelength range. Considering their size, the geometric optics approximation can be safely applied, and a specularly reflecting phase function will be used (Siegel and Howell, 1992):

p(cos φ) )

Fλ′(cos φ) Fλ

(20)

Spherical particles of smaller size, but still transparent, were also investigated. A suspension of latex particles made of poly(vinyltoluene) in water was stud-

3098 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

ied by Daniel and Incropera (1977). They have experimentally measured the “cross section” for scattering and the phase function. The calculated scattering coefficient (Cabrera et al., 1995) was averaged as indicated above, and the phase function was incorporated into the calculation after a precise fitting of the reported experimental values. Two types of titanium dioxide suspensions (at different concentrations) were also used. The absorption and scattering coefficients have been measured by Cabrera et al. (1996) using a recently developed technique. Once more, these values were averaged as before. In this case, considering the characteristics of the titanium dioxide particles, a diffuse phase function was used (Siegel and Howell, 1992; Cabrera et al., 1995):

8 p ) (sin φ - φ cos φ) 3π

(21)

Table 2 provides the information about the “lamp output power, wavelength-averaged” optical properties used in this work. Numerical Solution A revision on rigorous and approximate methods of solution has been presented by Ozisik (1973) and Santarelli (1983). Some of these methods were originally developed to investigate neutron transport in nuclear engineering problems; this transport is governed by equations having very similar integro-differential mathematical structure. Santarelli and co-workers (Spadoni et al., 1978; Pasquali et al., 1996) have used Monte Carlo simulation; instead of computing radiation intensities by application of eqs 7-10 and from these results calculating the LVREA, they have simulated the behavior of a statistically significant number of photons. The life of these photons is followed from their “birth” at the radiation source to their “fate” by absorption or by scattering out of the reactor boundaries. Among the approximate solutions, the discrete ordinate method (DOM) (Duderstadt and Martin, 1979) is recognized as the one that provides the most accurate results; besides, it is very powerful for tackling a great variety of problems with almost no restrictions. The DOM transforms the integro-differential equation (eq 7) into a system of algebraic equations that can be solved by machine computation. Cartesian geometries have been successfully used by Alfano et al. (1994, 1995), Cabrera et al. (1994, 1995), and Brandi et al. (1996). A triple-discretization procedure must be developed: (i) for the polychromatic wavelength range, (ii) for the spatial variable (x: r,z,β), and (iii) for the directional variable (Ω: θ,φ). As proposed in this work, the wavelength problem can also be approximated by using a single averaged lamp output power and averaged optical properties (absorption and scattering coefficients). To get the solution, essentially, a radiation balance considering all angular (directions) contributions must be solved for each cell in the spatial domain. Solving this balance, radiation intensities are so obtained for each direction. They define, in this way, the calculation path over the spatial computational mesh. When this radiation balance is so accomplished, special care is needed in order to include those terms that take into account the existing coupling between the angular variables that are used to characterize the flight of photons. Let us explain this concept. For the numerical

Figure 3. Schematic representation of the angular redistribution problem.

calculation, the own nature of the radiation propagation phenomenon (which depends on position in space and on the trajectory of the radiation beams) makes it necessary to generate spatial and angular cells. Regarding the spatial description, the annular system has symmetry in the β direction, and consequently, the spatial cells will be part of an r-z conventional mesh. Then, under these circumstances, the spatial cell for our balance will be a ring of rectangular (∆r × ∆z) cross section surrounding the tubular source. On the other hand, direction Ω is described in terms of the angular coordinates µ and η (Figure 1). Let us consider, for example, the case in which a photon moves from a point such as that defined by r1 to another defined by r2 with a direction Ω (Figure 3). When we move along this direction (the one that must be used to write the radiation balance), µ1 ) cos γr1 changes into µ2 ) cos γr2, while η ) cos γz stays constant. These changes in µ generate along the numerical calculation some sort of energy transfer from one angular cell to the neighboring angular cell. This behavior is usually recognized with the name “the angular redistribution problem”. This is a consequence of the employed curvilinear coordinate system for the spherical propagation of radiation in space. This very particular situation can be well understood looking at eq 7. In the cylindrical coordinate system, it can be seen that the directional derivative of the RTE is not only a function of the spatial coordinates (r,z,β) but of the angular ones as well. To account for this coupling of directions along the numerical calculation, some additional terms must be incorporated into the balance of the spatial cells. When, giving a point in space, the final integration over all possible directions is performed, all these terms cancel among themselves, giving rise to a result that is usually termed the “conservative form” of the radiation energy balance. This form must be achieved if a successful approximation with the discrete ordinate method is desired. On the other hand, it has been shown that if a conventional finite difference technique is applied to an equation such as eq 7, one obtains a “nonconservative” computational scheme (Duderstadt and Martin, 1979). The main conclusion derived from this work is the strong recommendation that the set of equations required for establishing the DOM computational scheme be obtained from a specific energy balance. This balance must be done in the cell of the computational mesh devised for each particular system under consideration, after the corresponding spatial and directional discretizations have been decided.

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3099

Figure 5. (a) Angular cell for direction Ωm and its location on the sphere of radius equal to unity. Mesh of angular cells on one octant of the sphere having radius equal to one. (b) For an S8 array and (c) for an S16 array. Table 3. Area and Volume Elements for Equation 23

Figure 4. Cross section of the spatial cell. For direction Ωm, a and d are contours for entrances and b and c for emergences.

Let us look now at a spatial cell having a volume Vi, j (Figure 4). The r-z section of the cell has the following characteristics: (1) its center is defined by ri, zj; (2) its sides are equal to (ri+1/2, ri-1/2) and (zj+1/2, zj-1/2); and (3) each one of these sides defines the corresponding areas of radiation transfer given by Ai-1/2, j, Ai+1/2, j, Bi, j-1/2, and Bi, j+1/2. In the figure, we have also identified points of entrance or emergence of radiation (a, b, c, and d) which are associated with the directional mesh. When the Ω space is made discrete, one gets a set of characteristic angular cells. Let us look at a typical cell (m) having a center characterized by direction Ωm. This cell has angular limits given by Ωm+1/2 and Ωm-1/2. These limits are defined by the direction cosines associated with this particular direction. Thus, we have associated a pair of direction cosines µm and ηm to every discrete direction in the computational scheme. Each angular cell (Figure 5a) affects the radiation field at the point under consideration according to a given weighting factor (wm). This weight is directly related to that part of the area corresponding to a sphere of directions having a unit radius associated with the Ωm trajectory and the limits Ωm+1/2 and Ωm-1/2. For a given wavelength λ, the energy balance on the cell Vi, jwm is given by j j µm(Ai+1/2, j I i+1/2, - Ai-1/2, j I i-1/2, )+ m m

ηm(Bi, j+1/2I i,mj+1/2 - Bi, j-1/2I i,mj-1/2) + in- and out-radiation fluxes corresponding to the spatial cell i, j i, j (Ai+1/2, j - Ai-1/2, j)(Rm+1/2I m+1/2 wm-1 - Rm-1/2I m-1/2 wm-1) +

angular redistribution term

βI i,mjVi, j absorption and out-scattering

σ )

M

∑I

4π n)1

i, j n pnmwnVi, j

(22)

in-scattering

In eq 22, pnm provides the phase function for radiation coming in from direction n that is incorporated into direction m [p(ΩnfΩm)]. Table 3 gives the values of A, B, and V for the cylindrical cell where the radiation balance is performed. In eq 22, Rm+1/2 and Rm-1/2 are effective areas for the angular fluxes corresponding to Ωm+1/2 and Ωm-1/2. They are originated in the angular redistribution due

coeff

value

Ai+1/2, j Bi, j+1/2 Vi, j

2πri+1/2∆zj 2π[(ri+1/2 + ri-1/2)2-1]∆ri 2π[(ri+1/2 + ri-1/2)2-1]∆ri∆zj

Table 4. Quadrant Definition and Contours for Radiation Entrances and Emergences quadrant

η

µ

incoming

outgoing

I II III IV

>0 >0