Research Article

# 3D Inverse Boundary Design Problem of Conduction-radiation Heat Transfer

ABSTRACT

A numerical analysis of a combined heat transfer problem was considered in this study. The inverse boundary design problem of combined conduction-radiation was solved. The aim of this study was to find the strength of heaters in a 3D cubic enclosure to produce desired temperature and heat flux distribution on the design surface. Conjugate gradients method was chosen to perform the iterative search procedure for obtaining the optimal solution. The direct problem was solved using the finite volume method for both energy and radiative transfer equations. Effects of conduction-radiation parameter, extinction coefficient, single scattering albedo and aspect ratio of the enclosure were studied. The results showed that the combination of a direct method and an inverse one could lead to precise estimation of the heater power distribution. Various problems of the pure radiation, pure conduction and combined radiation-conduction were solved and it was observed that estimated heat flux on the design surface was so close to the desired value that the maximum root mean square error was less than 2%.

 Services Related Articles in ASCI Similar Articles in this Journal Search in Google Scholar View Citation Report Citation Science Alert

Received: June 13, 2011; Accepted: January 02, 2012; Published: February 22, 2012

INTRODUCTION

Coupled conduction and radiation transport processes find applications in many engineering problems, such as glass processing, cooling of electronic devices and fire protection. When, the medium is participating and/or the boundaries have large values of emissivity, the radiative effects become more important. In this case, the energy equation provides the local temperature which determines the blackbody intensity in the Radiative Transfer Equation (RTE). Furthermore, solving the RTE gives the divergence of the radiative flux that is a source term in the energy equation (Mezrhab et al., 2008).

Not all the heat transfer problems have the analytical solution. This fact stems from the complexity of the governing equations and/or the boundaries and their conditions. Various heat transfer problems have been solved by numerical techniques in recent years (Abdelkader et al., 2007; Esmaeilzadeh et al., 2008).

In this study two main numerical techniques are used: the Finite Volume Method (FVM) for the direct problem and the conjugate gradients method for the inverse problem. The former is used extensively in analyzing the heat transfer and the fluid dynamics. Kwaire et al. (2009) solved a three phase flow a non-Newtonian fluid by FVM. Xingwei and Chaoying (2011) applied FVM to investigate the aerodynamic characteristics of a forward flight flapping airfoil. The modeling of the flow field around an elliptic airfoil is carried out by Benazza et al. (2007) using a finite volume method. The numerical study of an MHD problem in the steel solidification process is done by Mechighel and Kadja (2007). Sabeur-Bendehina et al. (2005) used the optimization techniques as well as finite volume method to find the location and the length of the partitions in a natural convective flow inside an undulated enclosure. Moraveji et al. (2011) numerically studied the heat and mass transfer in the fluidized bed dryer by coupling of the finite volume with semi implicit method for pressure linked equations (SIMPLE) method.

While, many works in thermal design are dedicated to 2D problems, the effect of 3D combined conduction-radiation is never taken into consideration. In this study an attempt was made to solve three dimensional inverse boundary design problem. The problem is solved for various cases and effects of conduction-radiation parameter, extinction coefficient, single scattering albedo and aspect ratio of the enclosure are studied.

FORMULATION AND NUMERICAL MODEL

Direct problem: The physical domain is shown in Fig. 1a. The cubic enclosure of unit length (L = 1 m) is filled with a grey, absorbing, emitting and scattering gas.

 (1)

It is assumed that thermal conductivity k is independent of temperature. The divergence of radiative heat flux is given by:

 (2)

where, κ is the absorption coefficient, Ib is the blackbody intensity and G is the incident radiation. The dimensionless parameters for this problem are introduced as:

 (3)

So, the dimensionless energy equation is expressed as:

 (4)

Here, ω is the single scattering albedo and NCR is the conduction-radiation parameter. The derivatives are with respect to dimensionless length. To find the divergence of radiative heat flux the RTE must be solved.

RTE for a grey, absorbing, emitting and scattering medium in the direction s can be written as:

 (5)

where, the source function can be defined as:

 (6)

 Fig. 1(a-c): (a) Physical domain, (b) A typical three-dimensional control volume and (c) A typical control solid angle

In Eq. 5 and 6, r and s are the position vector and the unit vector describing the radiative intensity direction, respectively. I is the radiative intensity; β, κ and σ are extinction, absorption and scattering coefficient, respectively. Φ is scattering phase function which is dependent on both in-coming (s’) and out-going (s) intensities.

The discretization procedure of RTE is according to that of Chai et al. (1994) and is not repeated here. Integrating the RTE (Eq. 5) over 3D control volume (Fig. 1b) and over a control angle (Fig. 1c), the final discretized equation can be derived as (Talukdar et al., 2005):

 (7)

The coefficients and the source term can be found by Talukdar et al. (2005).

The radiative flux at enclosure walls can be calculated as:

 (8)

Where:

 (9)

Hw represents the irradiation and by definition is the integration of all incoming directional intensities.

The total heat flux qtot is the summation of both conductive and radiative heat fluxes. In dimensionless form, total heat flux on any wall can be formulated as:

 (10)

 Fig. 2: Arbitrary enclosure geometry

Inverse problem: The general formulation of the inverse problem at hand may be described by referring to Fig. 2 which is a schematic of a 3D enclosure with arbitrary shape. Suppose that it is desired to have a specified temperature and heat flux distribution on Γd (i.e., the design surface). An inverse design problem may be defined as determining the required distribution of heater strengths on Γh (heater surface) such that it provides the desired conditions on the design surface. Therefore, in the inverse problem the boundary condition is unknown on the heater surface, while the design surface is constrained by two simultaneous boundary conditions. The boundary conditions on the other surfaces of the enclosure are all known. In order to solve these types of problem, the heater surface is discretized into NH subelements whose corresponding heat fluxes are represented by a vector as:

 (11)

The desired heat flux q*d on the design surface Γd is also assumed to be discretized into ND known components as the following vector:

 (12)

To estimate the unknown heat flux vector q*h, the temperature of design surface is set to be the actual boundary condition of the direct conduction radiation problem. Then by assuming an initial guess for the heat flux at the heater surface the FVM is invoked to calculate heat flux distribution on the design surface. The calculated value of q*c are then compared with the desired value of the heat flux on the design surface and if necessary, the heat flux distribution on the heater surface is varied. This procedure is repeated until sufficient proximity between the calculated and the design values of the heat flux are reached. This can be facilitated by defining an objective function F (q*h) as:

 (13)

The problem is then to find the vector q*h which minimizes function F. The values of q*c,m are components of the calculated heat flux as a function of q*h. The optimal point of function F corresponds to the solution of the problem.

The computational algorithm of this minimization procedure consists of two main modules; the direct conduction-radiation problem computation and the search modules.

For the first module the finite volume method is used to solve the energy equation. The central differencing scheme (Spalding, 1972) is chosen to discretize the diffusive terms. The directional intensities are also obtained by the finite volume method for radiative transfer and they are used to calculate the radiant heat fluxes. The intensities at integration points are approximated by nodal intensities with adopting the step spatial difference scheme (f = 1) to get the positive intensities.

To begin the numerical procedure for the direct problem, the radiant heat fluxes are estimated using FVM. In this step the temperature field which is used to calculate the T* term is from the previous iteration. Then to update the temperature field, the energy equation is solved. The temperature of the adiabatic walls is calculated from Eq. 10 which is a 4th degree equation of T*w (due to presence of qw term) and is solved using the Newton-Raphson method (Deuflhard, 2004). The solution is terminated when the steady state condition is reached and the following condition is satisfied:

 (14)

where, n is the iteration number and ζ is the independent parameter which can stand for I* and T*. For the inverse module, CGM is used as the basic search method to minimize function S. The complete procedure of conjugate gradients algorithm is explained by Ozisik and Orlande (2000). In what follows a brief description concerning the manner of computing the sensitivity matrix and updating the unknown vector used in the CGM is given. In order to establish CGM, we need to determine sensitivity coefficients of the calculated heat fluxes with respect to the heat fluxes of the heater elements. These sensitivity coefficients form sensitivity matrix Z and are defined by:

 (15)

where, ND is total number of elements on the Γd and NH is element numbers on Γh. A boundary value problem is developed for computation of sensitivity coefficients by derivation from Eq. 1 with respect to unknown flux q*h,n. This is done by solving the direct problem (Eq. 4, 7) with the following boundary conditions (Pourshaghaghy et al., 2006; Payan et al., 2009):

 (16)

Note that sensitivity coefficient problem is linear with respect to inverse procedure. Therefore, these coefficients are found prior to start of inverse computations. For this reason, it is sufficient to solve Eq. 4, 7 one time for each n (1≤n≤NH) and then the value of computed q*tot on each element (i) is obtained by Eq. 10 and if element (i) coincides on the design surface mth element, the value of resulted q*tot is assigned directly as Zm,n (m = 1, 2,…, ND).

With known sensitivity values, we are able to calculate q*h for the next step using iterative procedure of CGM:

 (17)

where, k represents the step number. The process of obtaining λk, the search step size and dk, the direction of descent, is fully explained by Ozisik and Orlande (2000).

Inverse analysis procedure: As said before, the strength of heat fluxes on heater surface is unknown variable in this problem. All the material properties as well as other boundary conditions are specified. To start the problem an initial guess for the unknown vector q0h is made. A brief description of the inverse algorithm is as follows:

 • The sensitivity matrix is obtained and stored • Knowing the heat flux on the heater surface, the direct problem is solved. The numerical procedure to solve the direct problem is already explained. In this step only the temperature (and not heat flux) condition is imposed on the design surface • Dimensionless total heat flux on the design surface qc*k is calculated and stopping criteria is checked. If ||∇Fk|| = ||-2ZT (q*d-qc*k)||5x10-6 the solution is terminated; otherwise we proceed to the next step • The heat flux of the next step qhk+1 is computed. k is replaced by k+1 and we return to step 2

RESULTS AND DISCUSSION

Some preliminary efforts are dedicated to find the optimal mesh size for a balance between accuracy and convergence rate. It is decided that, as also suggested by Talukdar et al. (2008), the domain is subdivided into 30x30x30 uniform control volumes in x, y and z directions except for the case of β = 5 where a 50x50x50 uniform grid is used. Solid angle is discretized to 12x12 uniform control angles in θ and n directions.

In all studied cases, walls of enclosure are assumed black.

The physical domain of the problem is a cubic enclosure. The configuration of heater and design elements are depicted in Fig. 3a. Ten equal sized heaters are located on the north, west and east walls and numbered by P in Fig. 3a, while south surface is divided into 10 design surfaces of equal size (numbered by M in Fig. 3a). The unnumbered top and bottom faces are assumed adiabatic.

The inverse problem is to find heat fluxes on the heater surfaces when two different thermal conditions of temperature (T* = 0.5) and heat flux (q*tot = -0.5) are imposed on all design elements.

 Fig. 3(a-c): (a) Configuration of heater and design elements, (b) and (c) Results of validation check

Validation check: Usually inverse algorithms are checked by the results of the same but direct problem. This idea is brought in this study. For the configuration of Fig. 3a with N = 0.01, ω = 0 and τ = 1 when q*h = 0.5 was applied on all the heater elements, the direct problem is solved and calculated heat flux on the design surface q*c is displayed in Fig. 3b. Then this value of q*c is used as the desired value in the inverse problem and we have searched for the unknown heat flux on the heater surface q*h. The estimated heater power is shown in Fig. 3c and it is in good agreement with the input data. This comparison corroborates the accuracy of the employed numerical technique.

Effect of conduction-radiation parameter: The problem is solved for five different values of NCR and distribution of heat flux on the heater surfaces is exhibited in Fig. 4. In all of these cases, radiative scattering is assumed negligible. Two of these five values of NCR are actually extreme cases of pure conduction and pure radiation and it helps to understand when each of these modes of heat transfer governs in the enclosure. At lower values of NCR, radiative heat transfer is more powerful in the enclosure. So, the distribution of heat fluxes on the heater surface is similar to the case of pure radiation in the enclosure. As it can be observed, in these situations, the heaters near the edges have higher power due to the closeness to the reflecting walls. This distribution of heater power is similar to the results of Erturk et al. (2002). As they explained, the design elements have higher values of sensitivity with respect to the heater elements near the edges. As NCR increases, conduction heat transfer dominates in the enclosure. As a result, heat flux distribution approaches the result of case of pure conduction. Because of the symmetry in the governing equations and boundary conditions for each heater element, the heater power distribution is almost the same for all heater elements for the case of pure conductive heat transfer. This is slightly different from what is reported by Yang and Chen (1997) and Huang and Wang (1999) due to difference in boundary conditions of the 3D inverse conductive problem. Also calculated heat flux on the design surface q*c is shown in Fig. 4b for different steps of inverse cycle for the case of NCR = 0 (pure radiation). It is interesting to note that initial guess for q*h produces results which are far from the desired values. But as inverse algorithm marches, q*h updates itself and reaches the value which produces the desired temperature and heat flux distribution on the design surface. The maximum relative error between the calculated and the desired value of the last inverse iteration is 0.32% and occurs at design elements of numbers 1, 2, 9 and 10. The reduction in the difference between the calculated and design heat fluxes is a common feature of an inverse algorithm which is observed and reported by many authors such as Pourshaghaghy et al. (2006), Sarvari et al. (2003) and Daun et al. (2006).

Effect of extinction coefficient: Extinction coefficient appears in equations of RTE and energy (in the definition of NCR). Because the effect of NCR studied above, for this case NCR is fixed at 0.01 and does not change with the extinction coefficient.

 Fig. 4(a-b): (a) Effect of conduction radiation parameter on heater power distribution and (b) Convergence history of the inverse problem for the case of pure radiation

 Fig. 5(a-b): (a) Effect of extinction coefficient on heater power distribution and (b) Effects of single scattering albedo on heater power distribution

Figure 5a shows the distribution of heat flux on heater surfaces for various values of extinction coefficient. When β grows, higher amount of radiative energy is absorbed by the medium. Heater elements near the edges transfer the radiative energy surface-to-surface but central heater elements need to transfer the energy through the absorbing medium. That’s why heater elements near the edges (i.e., far from the centre of enclosure) have stronger intensity than those located in the center. When β decreases, the medium absorbs lower amount of radiative energy and central elements produce more thermal energy.

Effect of single scattering albedo: The inverse problem is studied at four different values of the single scattering albedo ω and the results are illustrated in Fig. 5b. The effect of radiative scattering increases when ω enhances. According to Eq. (4), it’s obvious that the gradient of radiative flux diminishes as ω goes up. So, when the value of ω rises, the radiative effects become of minor importance and the same trend for the pure conductive problem is more likely to occur. Therefore, heat fluxes of heater surfaces come close to a constant value. It should be noted that in all these cases, the scattering phase function is assumed linear; i.e., Φ = 1+0.5 (s.s’). Results are the same for isotropic scattering (Φ = 1) but are not included for the sake of conciseness. The distribution of the heater power at lower values of ω has a U-shape form and is almost the same as the diagrams obtained by Pourshaghaghy et al. (2006) who solved a very similar inverse problem in 2D geometry in which radiation was the only mode of heat transfer.

Effect of aspect ratio of the enclosure: With the intention of study the effect of aspect ratio of enclosure, we shift our attention to a more realistic case. The same geometry that is previously investigated by Daun et al. (2006) is reconsidered here with almost identical configuration of heater and design elements. As shown in Fig. 6a the enclosure is a 3D rectangular enclosure with black surfaces. Design surface 6 is located at the south wall of the enclosure and is 0.1 m away from the edges on all sides. The inverse problem is analyzed in two different height-to-width ratios of H/W = 0.5 (as shown) and H/W = 0.3. Surfaces 1, 2 and 3 are heater surfaces while surfaces 4 and 5 as well as the margin of design surface is adiabatic. Heater surfaces are divided into 55 heater elements for the case of H/W = 0.5 (as shown in Fig. 6b) and 46 heater elements for the case of H/W = 0.3. Also as shown in Fig. 6c for both cases design surface is divided into 24 design elements where two simultaneous boundary conditions (T* = 0.5 and q*tot = -0.5) are imposed on all the design elements. Also in both cases NCR = 0.01, β = 1 and radiative scattering of the medium is assumed negligible. The distribution of heat fluxes on the heater surfaces is displayed in Fig. 7 and 8. Decreasing H/W leads to nearness of heaters on the top wall and the design surface. Consequently, it is expected that these heaters have more thermal power compared with the case of higher H/W. For both values of H/W, heater elements on the centre of top surface have more power compared with any other heater element. In general, the results of this study are similar to those reported by Daun et al. (2006) who considered the same geometry in which the radiative heat transfer was the only mode of heat transfer.

 Fig. 6(a-c): (a) Physical domain, (b) Discretized design surface and (c) Discretized heater surfaces

 Fig. 7(a-c): Heater power distribution for H/W = 0.3 on (a) Surface 2, (b) Surface 1 and (c) Surface 3

But because in this study, the convection as well as the radiation contributes in total heat transfer, the difference between the highest and the lowest heat fluxes is slightly less.

In Fig. 9 convergence history of the solution is demonstrated for the case of H/W = 0.5. It is seen that at the first iteration the calculated results are not acceptable at any of design elements because of the inaccuracy of the first guess. But as the solution proceeds, only the design elements at the corners of the design surface do not have the desired heat flux. At final iteration, as it is seen, the difference between the calculated and desired values is insignificant for all design elements. It can be concluded that the inverse algorithm is fully capable to produce the desired heat flux on the design surface.

Regularization of inverse problem by CGM introduces some errors in the inverse calculation.

 Fig. 8(a-c): Heater power distribution for H/W = 0.5 on (a) Surface 2, (b) Surface 1 and (c) Surface 3

 Fig. 9(a-d): Convergence history for the case of H/W = 0.5

To measure it quantitatively, root mean square error (abbreviated as Erms) is defined as follows:

 (18)

Lower value for Erms implies that calculated and desired results are more alike. Table 1 contains the values of Erms for all cases. It is important to notice that for all different problems, these values are less than 2% which is satisfactory for all engineering purposes. The values of Erms are larger than 1% for the scattering cases when the non-linearity of the problem enhances. The highest value of Erms is 1.92% when H/W = 0.5 and the distance between the heater elements and the design surface is relatively large. Additionally, the lowest values of Erms happen at the extreme cases when either pure conduction or pure radiation governs the heat transfer inside the enclosure.

 Table 1: Values of Erms for different problems 2Conduction radiation parameter, 3Extinction coefficient, 4Scattering albedo, 5Root mean square error

The precise prediction of heater power is due to robustness and accuracy of both direct and inverse algorithms.

CONCLUSIONS

In this study 3D inverse boundary design problem of combined conduction-radiation in participating media is numerically investigated. The coupling of FVM and CGM is found suitable and easily programmable for such a problem. In addition, the problem is solved for the extreme cases of pure radiation and pure conduction and it is found that all the heaters must produce the same value of heat flux for the case of pure conduction. The effects of some heat transfer parameters are also studied. It is observed that as NCR increases, the governing equations and results approach the pure conductive case. Moreover, the divergence of radiative heat flux diminishes as ω rises and the conductive heat transfer prevails. An inverse problem with many known and unknown elements is solved precisely by the presented algorithm. In all of the cases, the calculated and desired heat fluxes on the design surface have almost equal values. The results presented here can be useful in heat engineering applications such as furnace design, melting/cooling down process of the glass and ceramics, etc.

AKNOWLEDGMENTS

The authors wish to thank the Azad Islamic University, Abadan Branch, Abadan, Iran for their support throughout this study.

REFERENCES
1:  Abdelkader, S., R. Mebrouk, B. Abdellah and B. Khadidja, 2007. Natural convection in a horizontal wavy enclosure. J. Applied Sci., 7: 334-341.

2:  Esmaeilzadeh, E., A. Alamgholilou and H. Mirzaie, 2008. Numerical investigation on heat transfer enhancement of traverse ribs in 3D turbulent duct flow. Asian J. Applied Sci., 1: 286-303.

3:  Kwaire, A.A., H.M. Falih and H.A. Abdul Bari, 2009. Modelling of three phase system with Non-newtonian liquid using computational fluid dynamics model. Asian J. Ind. Eng., 1: 12-20.

4:  Xingwei, Z. and Z. Chaoying, 2011. Numerical investigation on the aerodynamic characteristics of a forward flight flapping airfoil with nonsymmetrical plunging motion. Inform. Technol. J., 10: 748-758.

5:  Benazza, A., E. Blanco and M. Abidat, 2007. 2D detached-eddy simulation around elliptic airfoil at high reynolds number. J. Applied Sci., 7: 547-552.

6:  Mechighel, F. and M. Kadja, 2007. External horizontally uniform magnetic field applied to steel solidification. J. Applied Sci., 7: 903-912.

7:  Sabeur-Bendehina, A., O. Imine and L. Adjlout, 2005. Optimization of the heat transfer rate in undulated enclosures with multiple partitions. J. Applied Sci., 5: 1237-1243.

8:  Moraveji, M.K., S.A. Kazemi and R. Davarnejad, 2011. CFD modeling of heat and mass transfer in the fluidized bed dryer. Trends Applied Sci. Res., 6: 595-605.

9:  Franca, F.H.R., O.A. Ezekoye and J.R. Howell, 2001. Inverse boundary design combining radiation and convection heat transfer. J. Heat Transf., 123: 884-892.
CrossRef  |

10:  Ozisik, M.N. and H.R.B. Orlande, 2000. Inverse Heat Transfer: Fundamentals and Applications. Taylor & Francis, USA., ISBN: 9781560328384, Pages: 330.

11:  Das, R., S.C. Mishra and R. Uppaluri, 2010. Inverse analysis applied to retrieval of parameters and reconstruction of temperature field in a transient conduction-radiation heat transfer problem involving mixed boundary conditions. Int. Commun. Heat Mass Transf., 37: 52-57.
CrossRef  |

12:  Pourshaghaghy, A., K. Pooladvand, F. Kowsary and K. Karimi-Zand, 2006. An inverse radiation boundary design problem for an enclosure filled with an emitting, absorbing and scattering media. Int. Commun. Heat Mass Transf., 33: 381-390.
CrossRef  |

13:  Sarvari, S.M.H., J.R. Howell and S.H. Mansouri, 2003. Inverse boundary design conduction-radiation problem in irregular two-dimensional domains. Numer. Heat Transf. Part B: Fundamentals, 44: 209-224.
CrossRef  |

14:  Kim, K.W. and S.W. Baek, 2007. Inverse radiation-conduction design problem in a participating concentric cylindrical medium. Int. J. Heat Mass Transf., 50: 2828-2837.

15:  Erturk, H., O.A. Ezekoye and J.R. Howell, 2002. Comparison of three regularized solution techniques in a three-dimensional inverse radiation problem. J. Quant. Spectrosc. Radiat. Transf., 73: 307-316.
CrossRef  |

16:  Abbas, K.A. and O.I. Mohamed, 2007. An inverse estimation method for surface film conductance during cooling of fish packages. Am. J. Food Technol., 2: 281-287.

17:  Chai, J.C., H.S. Lee and S.V. Patankar, 1994. Finite volume method for radiation heat transfer. J. Thermophys. Heat Transf., 8: 419-425.

18:  Talukdar, P., M. Steven, F.V. Issendorff and D. Trimis, 2005. Finite volume method in 3-D curvilinear coordinates with multiblocking procedure for radiative transport problems. Int. J. Heat Mass Transf., 48: 4657-4666.
CrossRef  |

19:  Spalding, D.B., 1972. A novel finite-difference formulation for differential expressions involving both first and second derivatives. Int. J. Numer. Method Eng., 4: 551-559.
CrossRef  |

20:  Payan, S., S.M.H. Sarvari and H. Ajam, 2009. Inverse boundary design of square enclosures with natural convection. Int. J. Therm. Sci., 48: 682-690.
CrossRef  |

21:  Talukdar, P., F.V. Issendorff, D. Trimis and C.J. Simonson, 2008. Conduction-radiation interaction in 3D irregular enclosures using the finite volume method. Heat Mass Transf., 44: 695-704.
CrossRef  |

22:  Daun, K., F. Franca, M. Larsen, G. Leduc and J.R. Howell, 2006. Comparison of methods for inverse design of radiant enclosures. J. Heat Transf., 128: 269-282.