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 nonNewtonian 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). SabeurBendehina 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.
Recently, researchers have shown much interest in the use of inverse techniques
for solving various thermal design problems which due to illposedness, cannot
be mathematically formulated by direct approaches. In many industrial processes
two simultaneous thermal conditions (i.e., both temperature and heat flux) are
imposed on some of the boundaries, called the design surface of the system.
In contrast, no boundary condition is imposed on some other boundaries which
are usually called the heater surface. That is the designer’s responsibility
to find the boundary condition on heater surface in order to satisfy both conditions
on the design surface. Having a trialanderror nature, the conventional design,
also named forward design, is based on the direct solution of the problem where
one and just one condition is imposed on each boundary (Franca
et al., 2001). So, designer has to guess the thermal condition on
the unconstrained elements, solve the direct problem and check both conditions
of the design surface (Franca et al., 2001).
If it is not satisfactory, a new guess is made and the calculations are repeated.
As expected, this may lead to a great number of iterations until reaching the
acceptable results (Franca et al., 2001). To
avoid this cumbersome and timeconsuming procedure, inverse algorithms are developed
to help designer find the desired output with a mathematically robust procedure.
Fundamentals and a good review on the available inverse methods for solving
heat transfer problems can be found by Ozisik and Orlande
(2000). There are many works with regards to inverse heat transfer including
purely radiative problems and radiation in combination with other modes of heat
transfer. Das et al. (2010) have considered the
problem of inverse simultaneous estimation of parameters in a rectangular enclosure
where conduction and radiation were the modes of heat transfer. They solved
the inverse problem by genetic algorithm. Pourshaghaghy
et al. (2006) have investigated the problem of inverse boundary design
in emitting, absorbing and scattering media. They used both Conjugate Gradients
Method (CGM) and variable metric method as inverse algorithm and compared the
accuracy and solution cost of both methods. Inverse boundary design of conductionradiation
problem in irregular 2D geometries has been reported in a work by Sarvari
et al. (2003). They defined a least square based function and tried
to minimize it by CGM. Kim and Baek (2007) have studied
the inverse conductionradiation problem in a concentric cylinders filled with
a participating media. As they stated, LevenbergMarquardt algorithm was selected
as inverse solver. 3D inverse boundary design in a purely radiative cubic is
studied by Erturk et al. (2002). They compared
CGM with two different methods to deal with the inverse problem and concluded
that CGM is an efficient method to solve inverse problem. The prediction of
the surface film conductance by transient temperature measurements in fish flesh
samples during air cooling in a chilled duct at a constant temperature of 1°C
is done by Abbas and Mohamed (2007). They utilized an
SFS method to solve the inverse problem. Their predicted temperature profile
was similar to the measured one and they could successfully calculate the surface
film conductance of the fish flesh samples.
While, many works in thermal design are dedicated to 2D problems, the effect
of 3D combined conductionradiation 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 conductionradiation
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.
The steady state energy equation for combined conductionradiation is as follows:
It is assumed that thermal conductivity k is independent of temperature. The
divergence of radiative heat flux is given by:
where, κ is the absorption coefficient, I_{b} is the blackbody intensity and G is the incident radiation. The dimensionless parameters for this problem are introduced as:
So, the dimensionless energy equation is expressed as:
Here, ω is the single scattering albedo and N_{CR} is the conductionradiation 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:
where, the source function can be defined as:

Fig. 1(ac): 
(a) Physical domain, (b) A typical threedimensional 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 incoming (s’) and outgoing (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):
The coefficients and the source term can be found by Talukdar
et al. (2005).
The radiative flux at enclosure walls can be calculated as:
Where:
H_{w} represents the irradiation and by definition is the integration of all incoming directional intensities.
The total heat flux q_{tot} is the summation of both conductive and
radiative heat fluxes. In dimensionless form, total heat flux on any wall can
be formulated as:

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:
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:
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:
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 conductionradiation 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 q_{w}
term) and is solved using the NewtonRaphson method (Deuflhard,
2004). The solution is terminated when the steady state condition is reached
and the following condition is satisfied:
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:
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):
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 Z_{m,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:
where, k represents the step number. The process of obtaining λ^{k},
the search step size and d^{k}, 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 q^{0}_{h} 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 q_{c}^{*k}
is calculated and stopping criteria is checked. If ∇F^{k}
= 2Z^{T} (q^{*}_{d}q_{c}^{*k})5x10^{6}
the solution is terminated; otherwise we proceed to the next step 
• 
The heat flux of the next step q_{h}^{k+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(ac): 
(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 conductionradiation parameter: The problem is solved for
five different values of N_{CR} 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 N_{CR
}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 N_{CR}, 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 N_{CR} 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 N_{CR} = 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 N_{CR}). Because the
effect of N_{CR} studied above, for this case N_{CR} is fixed
at 0.01 and does not change with the extinction coefficient.

Fig. 4(ab): 
(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(ab): 
(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 surfacetosurface 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 Ushape 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 heighttowidth 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 N_{CR} = 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(ac): 
(a) Physical domain, (b) Discretized design surface and (c)
Discretized heater surfaces 

Fig. 7(ac): 
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(ac): 
Heater power distribution for H/W = 0.5 on (a) Surface 2,
(b) Surface 1 and (c) Surface 3 

Fig. 9(ad): 
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:
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 nonlinearity 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 

^{2}Conduction radiation parameter, ^{3}Extinction
coefficient, ^{4}Scattering albedo, ^{5}Root 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 conductionradiation 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 N_{CR} 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.