Extension of the Reduced Integration Scheme to Calculate the Direct Exchange Areas in 3 D Rectangular Enclosures with Nonscattering Media

The evaluation of direct exchange areas (DEA) in zonal method is themost important task due to the heavy computer cost of multiintegrals together with the existing of singularities. A technique of variable transformation to reduce the fold of integrals, which was developed originally by Erkku (1959) to calculate the DEAs of uniformly zonal dividing cylindrical system, was extended by Tian and Chiu (2003) for nonuniformly zonal dividing cylindrical system with large thermal gradients. In this paper, we further extend the reduced integration scheme (RIS) to calculate the DEAs in three-dimensional rectangular system.The detail deductions of six-, five-, and fourfold integrals to threefold ones are presented; the DEAs in a rectangular system with assumption of gray medium are computed by the Gaussian quadrature integration (GQI) and the RIS comparatively. The comparisons reveal that the RIS can provide remarkable higher accuracy and efficiency than GQI. More interestingly and practicably, the singularities of DEAs can be decomposed and weakened obviously by RIS.


Introduction
In recent decades, the zonal method [1] has been being an outstanding tool for thermal radiation calculations, whether for the developments of the numerical methods themselves [2][3][4][5][6][7][8][9][10] or the applications in thermal engineering [11][12][13][14][15].It is well known that the most heavy computation cost of zonal method comes from the direct exchange areas (DEAs), which include the radiative properties of the participating medium and boundaries, together with the geometrical relations.After the DEAs, the total exchange areas (TEAs) can be solved, then the solution of zonal method.Usually, the calculations of DEAs involve six-, five-, and fourfold integrals for threedimensional systems, as can be seen in the following contents.Further, the analytical or exact solutions of DEAs are difficult or impossible for most cases.Therefore, considerable efforts have been made to reduce or simplify the computation cost of DEAs till now.
Larsen and Howell [3] defined the exchange factor method to reduce the computer cost using the reciprocity and symmetry.Sika [4] defined the DEAs of all kinds, say, surface to surface, surface to volume, volume to volume in axisymmetric cylindrical geometries with gray gas, and reduced multifold (four, five, and six) integrals to single-or doublefold integrals, and simplified formulae with high efficiency.Ma [5] illustrated that the GZM (generalized zoning method) can treat the anisotropic scattering using scattering-exchange areas; the GZM was used in one-dimensional case and compared with existing solutions.Goyhénèche and Sacadura [6] presented the explicit matrix relation for total exchange areas; the integration order can be reduced; then the TEAs in anisotropically scattering medium and anisotropically reflecting walls were obtained, and the weighted sum of gray gases model (WSGG) could be adopted in the zonal method.Tian and Chiu [7] deduced the formulas to reduce the multifold to less fold integrations by variable transformation for nonuniform zones in cylindrical systems.Two years later, Tian and Chiu [8] proposed a new approach to compute the DEAs by hybridization of the finite volume method (FVM) with the midpoint integration scheme; the FVM is for adjacent and overlapping zones; the midpoint integration scheme is for nonadjacent zones.El Hitti et al. [9] used zone method plus the thin layer approximation to simulate the glass sheets under high temperature condition; they evaluate the DEAs by flux planes approximation.At the same time, El Hitti et al. [10] evaluated the TEAs from DEAs of nongray system by replating algorithm.
The applications of zonal method can be found in many thermal processes, whether for offline analysis of parameter design or for online control of operating parameter optimization.Chapman et al. [11] performed a parametric study for a batch reheating furnace; the DEAs were evaluated by Monte Carlo method and the WSGG (weighted sum of gray gas) was adopted to consider different species of gases.Jeong and Ha [12] simulated the continuously annealing furnace using zonal method for radiation with the measured gas temperature and imaginary surface and discussed the correction of the soot effects on radiative properties.Barr [13] reviewed the thermal model of steady-state pusher-type reheating furnace; more details like the treatment of combination of radiation and forced convection were mentioned also.Chen and Jaluria [14] used the parallelization for DEA computation in their research.Steinboeck et al. [15] presented the thermal model of a slab reheating furnace, in which the heat conduction in the slab was solved by the weighted residual approach and the radiation was solved by the net radiation method based on zone method.
Recently, we have made the mathematical modeling of thermal processes both in conventional reheating furnace and regenerative reheating furnace [16] using zonal method to consider the radiative heat transfer.After cautious analyzing and comparing of above-mentioned approaches, the reduced integration scheme (RIS), which was described in detail in [7], is extended to simplify the computation of DEAs in rectangular systems, for uniform or nonuniform zone arrangement.Although Tian and Chiu [7] have already pointed out that this technique can be easily extended to three-dimensional rectangular enclosures, we still would like to give the detail deduction and results.More interestingly, the believable evidences of accuracy and efficiency will be given, especially for the treatment of singularities.This paper is organized as follows.In Section 2, the definition of DEAs for different zone pairs and detailed formulae are given.In Section 3, the detailed description of the RIS for DEAs is described.In Section 4, verifications and comparisons for specified cases are made.Finally, the last section gives the conclusions.

Definitions of DEAs in Rectangular Enclosures
In fact, the definitions of DEAs of all kinds of zone pair, say, surface-surface, gas-gas (or volume-volume), and gassurface (or volume-surface), can be found in [17,18].Here we present them briefly for the subsequent deductions.
Figure 1 shows the radiative exchange between surface zones   and   , whose areas are   and   , respectively.Then the formula of DEA between them can be obtained with the guarantee of energy conservation for semitransparent medium of extinction coefficient .While, as indicated in the title, the media considered here are nonscattering, hereafter, only the absorption coefficient   is used.Consider Similarly, as shown in Figure 2, the formula of DEA between two gas zones   and   , whose volumes are   and   individually, is as follows: The formula of DEA between one gas zone   and one surface zone   , whose volume is   and whose area is   individually, is as follows (see Figure 3): Due to the feature of radiative heat transfer, there exists the reciprocity for DEAs and it implies The energy conservation of radiative heat transfer leads to the following relations for surface zone of area   and gas zone of volume   , respectively: According to the reciprocity, only half of the DEAs need to be evaluated.The energy conservation can be used to check the accuracy of computation.

Simplified Calculations of DEAs by Reduced Integration Scheme
3.1.Double-Fold Integral into Three Single-Fold Integrals.As described in [7], the original technique to transform the double-fold integral into three single-fold integrals was proposed by Erkku [19].Here we briefly give its principle.
There is a twofold integral; for example, where the assumptions are made that  and  are constants,  ≥  > 0, and  is a function of  − .
Set  =  +  − ; we get After the variable substitution, the integral can be executed in three regions as shown in Figure 4:    Please note that this variable transformation is available for all computation of DEAs in rectangular systems.Now the twofold integral ( 7) is transformed to three onefold integrals.As pointed out by Tian and Chiu [7], the computation time can be saved remarkably using (11) instead of (7).

Simplification of DEAs in Rectangular System
. Considering a rectangular enclosure as shown in Figure 5, where the numbers mean the zone number of surface zone and the numbers with circle mean the zone number of gas zone, the simplification of multifold integral of DEAs is given for gas-gas, gas-surface, and surface-surface, sequentially.Please note that the lengths of each zone (gas or surface) in three directions are not necessary to be equal; in another words, the zone dividing can be nonuniform.
(i) Gas-Gas Zones     .In Cartesian coordinates system, formula (2) can be rewritten as where Set and assume Δ  ≥ Δ  ; then one can get the following after integral transformation: This time, Set  =  1 −  2 ,  =  1 −  2 ,  =  1 −  2 and use the RIS which is described in Section 3.1; formula ( 14) can be expressed as where (they are used later again) If the zone dividing is uniform in all directions, the above integrals can be further simplified as where (ii) Gas-Surface Zones     .Taking the same meanings of symbols after formula (12), formula (3) becomes There exist six cases for the surface zone   according to its different locations when the gas zone   is fixed; they are as follows: (A) the surface zones are perpendicular to the -axis but away from the origin, as the zones numbered 13∼16 in Figure 5(c), for example; (B) the surface zones are perpendicular to the -axis but in the same plane with the origin, as the zones numbered 17∼20, which are opposite to the zones of case (A), in Figure 5(c), for example; (C) the surface zones are perpendicular to the -axis but away from the origin, as the zones numbered 5∼8 in Figure 5(b), for example; (D) the surface zones are perpendicular to the -axis but in the same plane with the origin, as the zones numbered 1∼4, which are opposite to the zones of case (C), in Figure 5(b), for example; (E) the surface zones are perpendicular to the -axis but away from the origin, as the zones numbered 21∼24 in Figure 5(b), for example; (F) the surface zones are perpendicular to the -axis but in the same plane with the origin, as the zones numbered 9∼12, which are opposite to the zones of case (E), in Figure 5(b), for example.
Consider the case (A); now the cosine function in (3) is Then Set After changing the order of integration, one can get where and the capital  means the distance along the -axis.
Using the RIS, formula (24) becomes (, )  (, )  ()    1 . ( If the gas zone can be projected onto the surface zone completely along one direction, or its projection can overlay the surface zone completely after parallel moving along one direction, formula (26) can be further simplified as As to cases (B)∼(F), the cosine functions and () are different, but the similar formulae as (26) and ( 27) can be obtained.
(iii) Surface-Surface Zones     .As shown in Figure 5, there are two cases for the pairs of surface-surface zones as follows: (A) one surface zone is perpendicular to the other; (B) one surface zone is parallel to the other.
For case (A), taking one surface zone on the - plane and another one on the - plane, for example, the cosine functions in formula (1) are Then formula (1) can be rewritten as where (30) Set After changing the order of integration, one gets Using the RIS, formula (32) becomes For other situations, one surface zone on the - plane or parallel to - plane, and another one on the - plane or parallel to - plane, for example, the deductions can be carried out similarly.
For case (B), taking two surface zones, which are perpendicular to the -axis, for example, the cosine functions in formula (1) Then formula (1) becomes where (37) Set After changing the order of integration, one gets where Using the RIS, formula (39) becomes If the lengths of both surface zones are equal in and directions, separately, formula (41) can be further simplified as The formulae of other situations for case (B) can be deduced in the same way.We state that the situations with shielding effect do not belong to the cases studied here.
Up to now we know that, if the zone grids along three directions are uniform, the integrals of gas-gas zones can be reduced from one sixfold integral to eight threefold integrals and those of gas-surface zones from one fivefold to four threefold integrals.For the case of parallel surface zones, one fourfold integral is reduced to four twofold integrals; for the case of perpendicular surface zones, one fourfold integral is reduced to two threefold integrals.If the quadrature point number  is chosen for the single-fold integral, the  6 and  3 points are needed for a sixfold integral and a threefold integral, separately.Therefore, the quadrature points needed for one DEA of gas-gas zones can be reduced from  6 to 8 3 , from  5 to 4 3 for one DEA of gas-surface zones, and from  4 to 2 2 (perpendicular) or 4 2 (parallel) for DEA of surface-surface zones.Once the integral regions are determined, the greater  is, the higher the accuracy is.On the other hand, more computational time can be saved with the increasing of , and more higher accuracy can be achieved as well.

Cases Study and Comparisons
Returning to Figure 5 again, we take the cubic enclosure with the sizes of 2 m × 2 m × 2 m and assume the medium is gray with   = 0.5 m −1 to check the accuracy and efficiency of the RIS for DEAs.The cube is divided into 8 gas zones (numbers with circle) and 24 surface zones uniformly.Therefore the side length of each gas zone and surface zone is 1 m.All the following comparisons are made on a personal computer with Intel(R) Core (TM)2 Duo CPU P8700 (2.53 GHz) and 3.0 GB (DDR2) memory.
Same as in [7], the Gaussian quadrature integration is adopted as the direct method for DEAs themselves or the reduced three-or twofold integrals for comparisons.Although there may be some other simple methods which can be used for DEAs or the reduced three-or twofold integrals, we will not do much more comparisons between the RIS and them.In the following tables and figures, the symbols GQI and RIS mean the results by Gaussian quadrature integration and reduced integration scheme, respectively.

Accuracy Validation and Comparisons.
Taking the number of integral steps  = 8 (the corresponding integral step length is 1/ = 0.125 m) in each dimension, the DEAs of gas zone  1 to all gas zones   ( = 1, . . ., 8) in the enclosure are computed by GQI and RIS separately, and the results are listed in Table 1.It can be seen that both GQI and RIS almost give the same results except for the DEA of  1  1 .The difference for  1  1 between by GQI and by RIS comes from the so called "self-irradiation" which exits in  1  1 .The "selfirradiation" is one kind of singularities, which happens in the calculation of DEAs of gas zone to itself.Those results of DEAs of gas zone  1 to all surface zones   ( = 1, . . ., 24) in the enclosure are listed in Table 2.It is clear that both GQI and RIS almost give the same results except for  1  1 ,  1  9 , and  1  17 .The differences between by GQI and by RIS are due to the singularity when surface zones are adjacent to the gas zone  1 .Those results of DEAs of surface zone  1 to all surface zones   ( = 1, . . ., 24) in the enclosure are listed in Table 3. Again, both GQI and RIS almost give the same results except for  1  9 and  1  17 .The differences between by GQI and by RIS are caused by the singularity that the surface zones  9 and  17 are just neighbor to  1 .
It is a common knowledge that the accuracy of integrals will increase with the number of integral step  (the corresponding integral step length is 1/).Here we choose  1  1 ,  1  1 , and  1  9 as examples to compare the DEAs by GQI and by RIS against  in each dimension.The results are illustrated in Figure 6.Here, in plot (a), the results by GQI are given only for  being up to 16 because one week or even more days may be needed when  is up to 32 on the present personal computer.Obviously, the results by GQI are more sensitive to , and the relative errors (that is, with respect to the final result of RIS at  = 32) of  1  1 (plot (a)) and  1  1 (plot (b)) at  = 16 are still greater than 20%, while the results of RIS almost go to convergence when  ≥ 8.Note that the singularities are involved in both  1  1 and  1  1 as mentioned above.Thus we know that the RIS can provide much higher accuracy than the GQI especially for those DEAs in which the singularities are included.Figure 7 depicts the differences of DEAs by GQI and RIS for various .Whether for DEAs of  1   ( = 1, . . ., 8) as shown in plot (a) or for DEAs of  1   ( = 1, . . ., 24) as shown in plot (b), the differences between GQI and RIS are almost invisible except for those DEAs who include the singularities.However, the differences of both for DEAs with singularities,  1  1 ,  1  1 ,  1  9 , and  1  17 , for example, are decreasing remarkably when  increases from 2 to 8.For higher values of , the compared differences are not given in Figure 7 since almost several weeks computer time is needed by GQI on our personal computer.
As mentioned in Section 2, the energy conservation of DEAs in one enclosure can be applied as the criterion to check the accuracy of numerical methods.Taking the gas zone  1 as the example, there is the exact value of 2 for formula (6); in other words, ∑ 8 =1  1   + ∑ 24 =1  1   = 2 m 2 .The results by  into several ones as described in Section 3.1; at the same time, the singularity is also decomposed and weakened.
As to the singularities that appeared in the computation of DEAs, there exit three cases.The first case is the DEA of one gas zone to itself, the second case is the DEA of one gas zone to its adjacent surface, and the third case is the DEA between two surfaces with an intersecting line.In the computations of all DEAs by RIS, the folds are always reduced.That means the possibilities of singularity can be reduced by RIS, just as said in [7] like "reduction of severity of singularities."

CPU Time Comparisons.
The best way to evaluate the efficiency of one approach is the comparison of CPU time.Again we first take the single DEA of  1  1 ,  1  1 , and  1  9 , respectively, to compare the CPU time cost by GQI and RIS as  increases from 2 to 64.The results are listed in Table 5. There, the reason for the absence of some values by GQI is just because almost several weeks are needed for multifold of DEAs.For gas-gas DEAs, the CPU time of GQI is about 10 4 times that of RIS; for gas-surface and surfacesurface DEAs, they are about 10 3 and 10 2 times, separately.This phenomenon reveals that the larger folds of the original integrals are, the more CPU time can be saved when using RIS.
Then taking the sum of DEAs of     (,  = 1, . . ., 8) and     ( = 1, . . ., 8;  = 1, . . ., 24) as examples, CPU time costs by both GQI and RIS for different  = 2, 4, 8 are listed in Table 6.As one can see, at least one hundred times of CPU time can be saved using RIS instead of GQI for the same .With the increasing of , even more CPU time can be saved when using RIS.

Conclusions
The reduced integration scheme is extended to calculate the direct exchange areas in three-dimensional rectangular system.The detailed reductions of multifold (six-, five-, and fourfold) integrals to threefold integrals are presented for nonuniform zonal dividing.A uniform zonal dividing case with gray medium is adopted to comparatively check the accuracy and efficiency of the reduced integration scheme.The standard Gaussian quadrature integration is also applied for the purpose of comparisons.The comparisons reveal that only 10 −4 times, 10 −3 times, and 10 −2 times CPU of Gaussian quadrature integration for DEAs of gas-gas with sixfold integrals, gas-surface with fivefold integrals, and surfacesurface with fourfold integrals are needed, respectively.More attractively, the reduced integration scheme can remarkably decompose and weaken the singularities, which exist in the DEAs of adjacent zones or the volume zones themselves, obviously.Distance between two elements, m :

Figure 3 :
Figure 3: Radiative exchange between gas zone   with volume   and surface zone   with area   .

( 9 )
Changing the order of integration  , we get − )  .

Figure 4 :
Figure 4: Regions of integration after the transformation.
Figure 1: Radiative exchange between surface zone   with area   and surface zone   with area   .
j Figure 2: Radiative exchange between gas zone   with volume   and gas zone   with volume   .

Table 4 :
Verifications and relative errors of conservation for ∑   = 2 m 2 with the change of .
.One can see that the relative errors of RIS with respect to the exact value are even lesser than 0.1% when  is only up to 4, while, for the same conditions, the relative errors of GQI are still greater than 1%.The reason may be mainly because of the eliminating of singularities of DEAs by GQI and RIS differently.In RIS, one integral region is always decomposed

Table 5 :
Comparison of single DEA CPU time in seconds for different .