Computation of Radiant Heat Transfer on a Non-Orthogonal Mesh Using the Finite-Volume Method

This article was downloaded by: [University of Glasgow] On: 20 October 2012, At: 06:29 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/unhb20 COMPUTATION OF RADIANT HEAT TRANSFER ON A NONORTHOGONAL MESH USING THE FINITE-VOLUME METHOD E. H. Chui a & G. D. Raithby a a Department of Mechanical Engineering, University of Waterloo, Waterloo, Canada N2L, 3G1, Canada Version of record first published: 23 Mar 2007. To cite this article: E. H. Chui & G. D. Raithby (1993): COMPUTATION OF RADIANT HEAT TRANSFER ON A NONORTHOGONAL MESH USING THE FINITE-VOLUME METHOD, Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology, 23:3, 269-288 To link to this article: http://dx.doi.org/10.1080/10407799308914901 PLEASE SCROLL DOWN FOR ARTICLE Full terms and conditions of use: http://www.tandfonline.com/page/terms-and-conditions This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. The publisher does not give any warranty express or implied or make any representation that the contents will be complete or accurate or up to date. The accuracy of any instructions, formulae, and drug doses should be independently verified with primary sources. The publisher shall not be liable for any loss, actions, claims, proceedings, demand, or costs or damages whatsoever or howsoever caused arising directly or indirectly in connection with or arising out of the use of this material. Numerical Heat Transfer, Part B, vol. 23, pp. 269-288, 1993 COMPUTATION OF RADIANT HEAT TRANSFER ON A NONORTHOGONAL MESH USING THE FINITE-VOLUME METHOD E. H. Chui and G. D. Raithby Departmenr of Mechanical Engineering. Universiry of Mrerloo, ffiterloo, Camah, Canada N2L 3Gl The finite-volume method has been shown to effictively predicl mdinnt exchange in geometrically simple enclosures where the medium is gray, absorbing, emitting, ond scat- tering. Cariesinn and circular cylindrical meshes have always been used. The present article shows that the method applies equally well lo geometrically complex enclosures where nonorihogonal, boundnry-fitted meshes are used. This development pennits radian! heat transfer to be computed on the same mesh employed to solve the equations of fluid motion. INTRODUCTION To predict the fluid flow and convective heat transfer in the geometrically complex systems of practical interest, computational methods are normally implemented on non- orthogonal, boundary-fitted meshes. Many such problems also require prediction of the radiant heat transfer, where the medium affects the heat transfer through emission, absorption, and scattering. The radiative transfer equation (or RTE) governs this radiant exchange. Although it is desirable to solve this equation on the same computational mesh used for the fluid flow, there has been very little discussion in the literature on solution methods that use nonorthogonal meshes (see reviews by Howell [I] and Viskanta and Mengiiq [2]). The spherical harmonics method of Mengiiq and Viskanta [3, 4]), the discrete ordinates method of Fiveland [5], the zonal method of Larsen and Howell [6], the six flux method of Siddall and Selquk [7], and the discrete transfer method of Lockwood and Shah [8] are all presented using rectangular or cylindrical meshes. The finite-volume method (FVM), developed by the present authors, has also been applied only on Cartesian 1:9] and cylindrical [lo] meshes. The goal of this article is to demonstrate the application of this method on nonorthogonal, quadrilateral element meshes. This work was supported by a contract from the Combustion and Carbonization Research Laboratory, CANMET. of the Energy Mines and Resources Depanment of the Canadian Federal Government, under the scientific authority of P. Hughes. The authors are also grateful for the additional financial assistance received from the Natural Science and Engineering Research Council of Canada. Copyright O 1993 'Itrylor & Francis 10407790/93 $10.00 + .00 D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 270 E. H. CHUI AND G. D. RAITHBY NOMENCLATURE area of surface panel f. m2 radiant intensity, wl(m2 sr) inscattering intensity. Eqs. (2) and (7) blackbody intensity, a@/* intensity at boundary intensity at integration point if within w' intensity at node P within w' upstream intensity for integration point on surface panel f within w' absorption coefficient, llm number of control volumes in the x and y directions number of control volumes in the polar (8) and azimuthal (4) direc- tions unit normal to a surface of control volume radiant heat flux, w/m2 radiant energy within o' crossing surface panel f, W position vector intensity source function within w', Eq. (6) distance in direction of s, m unit vector in direction of intensity path length, m temperature, K volume. m3 weight on node n within a quadrilat- eral element based on bilinear a p proximation. Eq. (9) Cartesian coordinates local coordinates for a quadrilateral element. Fig. 2A Kronecker delta function emissivity polar angle measured from z axis extinction coefficient, K, + a,. I lm Stefan-Boltzmann constant scattering coefficient. Ilm azimuthal angle measured from x axis in x-y plane scattering phase function, Eq. (2) new variable for phase function. Eq. (8) scattering angle, between s and 8' in Eq. (2) scattering albedo, aJ(a, + KJ solid angle. sr Subscripts and Superscripts b f if I nb N, S, E. W NW, NE. SW. SE P 4 W refers to blackbody refers to panel f on the surface of a control volume refers to the integration point on surface panel f associated with discrete solid angle w1 refers to neighboring volumes north, south, east, west northwest, northeast, southwest, southeast refers to nodal point P interpolation location for intensity upstream of point if enclosure wall Several features make the FVM attractive for the Dresent Dumoses. The method . . was originally derived with nonorthogonal meshes in mind. It also guarantees global conservation of radiant energy, and captures the diffusion limit exactly for a strongly participating medium. The FVM has been shown to give accurate results for rectangular and cylifldrical enclosures, anisotropic scattering is easy to treat, and rapid convergence can be achieved over a large range of optical thickness [ll]. In the sections that follow, the discretization of the RTE on the nonorthogonal mesh is described, a new solution procedure is outlined, details related to special treat- ment at the boundaries are given, and the results of four test problems are presented. For simplicity, only two-dimensional problems are considered, and the medium is assumed to be gray. D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 RADIANT HEAT TRANSFER ON NONORTHOGONAL MESH 271 RADIATIVE TRANSFER EQUATION The transfer of radiant energy is governed by the radiative transfer equation (RTE) [12]: It describes the variation of radiant intensity I(r, s) at location r in direction s. The four terms on the right side account for absorption, outscattering, emission, and inscattering, respectively. The inscattering term is defined by where I(r, s ' ) is the incident intensity from direction s f and a is the scattering phase function between s ' and s. The principal goal here is to obtain a discrete representation of I both in space and direction using a finite-volume approach. SPATIAL GRID Figure 1A shows how an irregular geometry can be subdivided using a mesh made up of quadrilateral "elements." Each grid line intersection defines a node, denoted by a heavy dot. AU dependent variables are stored at the nodes, and the equation for each variable is obtained from its discretized conservation equation for a volume surrounding the node. The volume associated with an interior node, such as P i n Fig. IA, is the sum of the four cross-hatched quadrants shown in Fig. 1B. The outside surfaces of the quadrants (the dashed lines) lie along lines that joint the midpoints of opposite sides of the element. The eight surface panels of the control volume are denoted by f - 1, 2, . . . , 8. For panel f, if is the integration point lying at the center of the panel, A, is its surface area, and n, is the unit outward normal. For nodes adjacent to the boundary, such as P in Fig. lC, the volume extends right to the physical boundary, and neighbor- ing nodes lie on the boundary. For two-dimensional problems, the volume is assumed to have unit depth normal to the plane of Fig. 1. Except for minor details related to the volume adjacent to a surface, the grid system just described was proposed by Baliga and Patankar [13] and used to make fluid flow computations by Schneider and Raw [14]. ANGULAR DISCRETIZATION The directional dependence of radiant intensity at every spatial node is captured by subdividing direction into discrete, nonoverlapping, solid angles of sue w'. Total flexi- bility is allowed in choosing the individual size of o' provided ~ w ' - 4n. In practice, w' is defined by a range of polar angle 0 and of azimuthal angle 4. D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 E. H. CHUl AND G. D. RAITHBY Fig. 1 (A) Sample computational mesh; (B) an interior control volume formed from assembly of four quadrants; (C) a control volume adjacent to the boundary. FINITE-VOLUME APPROXIMATION Implementing the FVM, the RTE [Eq. (I.)] is integrated over control volume Vp and solid angle w' to yield the conservation constraint that the net radiant energy leaving through the surface of Vp within w' equals the net generation of radiant energy within Vp and w' by emission, absorption, and scattering. Approximating this generation term by nodal point values, the conservation balance becomes where Q; is the radiant energy within w' that crosses the surface panel f, 1; is the intensity within w' at node P, while I,,, and are, respectively, the nodal blackbody intensity and the inscattering term. Defining I; as the intensity at integration point if within w' , D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 RADIANT HEAT TRANSFER ON NONORTHOGONAL MESH 273 where s is a unit vector in the direction of radiation and the integral D; is evaluated analytically. Using a Taylor series expansion about integration if, Raithby and Chui [9] introduced the closure relation where I$ is the intensity at point uf (see Fig. 1B) that lies distance S upstream from point if. K - KO + u, is the extinction coefficient and variable R', defined in Eq. (6), embodies the emitting and scattering characteristics of the medium. The first term on the right side of Eq. (5) represents the remnant of intensity I:, reaching if, while the last two terms incorporate the net change of intensity due to emission, absorption, and scattering from uf to if. Retention of the last term is crucial for the FVM to capture the optically thick limit. The introduction of Eq. (5) generates three new unknowns: R;, ( ~ R ' I ~ S ) ~ and 1;. These must be related to nodal values before Eq. (3) can be solved. Anention is next turned to this matter. Calculation of R', and iaR'las), The variable R', is defined as where I, - (&IT) is the blackbody intensity and i' is the inscattering term evaluated by approximating Eq. (2) with where the value of 6 is obtained by exact or very accurate numerical integration of the phase function @. This formulation was discussed and tested by Chui et al. [lo]. Since R' is specified at nodal locations (such as nodes 1 to 4 in Fig. 2), R; at integration points can be obtained by bilinear interpolation of nodal values; using the 0, y coordinate system shown in Fig. 2, D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 274 E. H. CHUI AND G. D. RAITHBY A €3 Ng. 2 Determination of integration point intensities on edges of quadrants Q1 to Q4. where w, is the bilinear interpolation weight on node n for the integration point if specified by the local coordinates (flu, yJ. Differentiating Eq. (9), rewritten for general location @, y ) where ax/& and aylas are direction cosines in direction s in the global x, y coordinate system. For a given grid, w,, awnlax and awnlay are calculated analytically once and stored. Equations (9) and (10) are then used when needed to determine R:, (13R'las)~ at all integration points. D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 RADIANT HEAT TRANWER ON NONORTHOGONAL MESH 275 Calculation of & fd is the intensity at point uf upstream of integration point if. For example, points ul and u2 associated with integration points if - il and i2, respectively, are shown in Fig. 2A. To find uf, the ray in the center of w' at point if is traced backward until it intersects with the surface of the element. &can then be obtained by linear interpolation (e.g., I:, can be obtained in terms of I: and I: in Fig. 2A). This practice leads to an equation for I; that involves nodes downstream from node P; this complicates the solu- tion of the equation, may lead to "wiggles" in the intensity, and reflects a physically incorrect influence. By slightly sacrificing accuracy, this downstream dependence can be avoided. For example, for quadrant Ql in Fig. 2A and for the I direction shown, the following approximation is used: I:, - I:, - I:. For Q2 in Fig. 28, I:, would be linearly interpolated from 1: and I:, as would I!,, from 1: and I:. Algebraic Equation for 1; With R;, ( a ~ ' l a s ) ~ and 1; defined, Eqs. (4) and (5) can be substituted into (3) to obtain an algebraic equation relating I', to the intensities at neighboring nodes: where the summation of nb is over the "upstream" neighbors. For the I direction shown in Fig. lB, Eq. (11) takes the form Equation (1 1) for each node together with boundary conditions provide the equa- tion set for the nodal intensities in direction I. Similar equations are formed for other directions; the nodes that lie upstream of course change with I. SOLUTION PROCEDURE Ovewiew hovided scattering is not present, I', in Eq. (1 1) is independent of downstream intensities. If the solver visits the nodes in the correct order, all terms on the right side of the equation will be known so I', at each node can be found by direct substitution. This allows the solution to be obtained by moving from node to node in the optimal "march- ing order." For a Cartesian mesh it is immediately clear what the marching order must be, but for the irregular grid in Fig. 1A the order is no longer obvious. To find a solution on irregular grids, there are two options. The equations could be solved by repeatedly sweeping across the grid until the convergence is achieved, without regard to the optimal order. Alternatively a "marching-order map" could be created that gives the optimal order in which the nodes should be visited. Once the angular discreti- zation is fixed, a marching-order map can be constructed for each intensity direction, and this map can be stored for repeated use. Because the map is easy to generate, because it needs to be revised only when the spatial grid or angular discretization is D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 276 E. H. CHUl AND G. D. RAITHBY changed, and because it reduces solution costs, the creation of the map is the preferred alternative. To obtain a preview of how a marching-order map might appear, the reader is referred to Fig. 4B. For the U-shaped region, the optimal order in which the nodes should be visited to solve for the intensity in the particular direction s is V1, V2, . . . , V15. To create the map, the boundary elements are swept first to find the starting location. Once this is found, the entire marching sequence is derived by repeated appli- cation of a simple algorithm. Details are now provided. Marching-Order Within an Element The element shown in Fig. 2 contains four quadrants. Quadrant Q1 is part of the control volume for node 1, which requires the intensities I:, and which in turn require I:, and I,!,. These and values of I:, and are similarly required for the other control volumes associated with quadrants 42, Q3, and 44 . According to the discretiza- tion scheme already described, I:, -. I:, - I:, so the equation for I: is independent of I:, I: and I:. Clearly, node 1 is farthest "upstream." Similarly, the equation for I: depends on I: but not on 1: or I:, and I: depends on I: but not on I: or 1:. Therefore, nodes 2 or 4 are equally upstream, but not so far as node I. Node 3 is least upstream. The order in which the nodes in Fig. 2 should be visited is therefore 1-2or4-3, where 2or4 means the order may be 2-4 or 4-2. For a different direction s, the order may be different. The general algorithm for node ordering within an element is shown in Fig. 3. Four sectors, denoted by a, to p4, are defined by the quadrant boundaries. Once the sector, in which the intensity direction s falls, is determined, the table in Fig. 3 specifies the marching order. The special cases that arise for boundary elements are also shown. The stars denote locations where the intensity is known from boundary conditions. Generation of the Marching-Order Map The elements that make up the grid for a U-shaped region are shown in Fig. 4A. The node numbering is converted from the local numbering in Fig. 2 to a global number- ing in Fig. 4 that uses (i, J ] , where i numbers the node along one set of grid lines and j along the other. To find the starting location, the boundary elements are visited in any order until an interior node is found for which the only nodes farther upstream are boundary nodes that have known intensities from boundary conditions [e.g., in Fig. 4A, the intensities at nodes (1, I), (2, l) , (3, I), etc., are known from boundary conditions]. Applying the element algorithm in Fig. 3 to element el in Fig. 4A shows that node (3, 2) is farthest upstream among the unknown intensities in that element. But from element e2, node (2, 2) is farther upstream than (3, 2) and from e3, there is no unknown intensity that lies still farther upstream. Node (2, 2) is therefore a suitable starting location. This is denoted in Fig. 4B as V1. With the farthest upstream node known, the map is generated by systematically finding which nodes are next farthest upstream. Applying the element algorithm to the elements that share the starting node (2, 2) shows that either (2, 3) or (3, 2) is next. D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 BOUNDARY B C Fig. 3 Variation of marching order with the sector p in which a given intensily direction s lies. Superscript * denotes boundary intensities that are known. D ow nl oa de d by [U niv ers ity of G las go w] at 06 :29 20 O cto be r 2 01 2 278 E. H. CHUl AND G. D. RAITHBY Fig. 4 (A) Determination of the star