# Multimode heat transfer with radiation

## Contents

Figure 1: Infinitely long parallel conducting plates with thickness δ and thermal conductivity k.

When radiation is combined with conduction and/or convection, the problems are called conjugate or multi-mode heat transfer. If temperature is one of the unknowns in the problem being studied, then such problems are highly non-linear, in that the fundamental radiation terms involve temperature to the fourth power, while conduction and convection depend on the first power of temperature or its derivatives. Property variations with temperature may add further nonlinearities.

The mathematics of highly nonlinear equations is not well-developed, and almost all combined mode problems are solved numerically; closed-form analytical solutions seldom exist. As an example of a combined-mode radiation-conduction problem, consider the geometry and conditions shown in Fig. 1. Radiation transfer occurs between the parallel plates and from the plates to the environment. Because the plates are not isothermal, conduction occurs within the plates, and will modify the temperature distributions, which are the unknowns in the problem.

For this case, the usual net radiation equations are as derived as following (Faghri et al., 2010):

$\frac{{{{q''}_{1,rad}}({x_1})}}{\varepsilon } = \sigma T_1^4 - \sigma T_2^4{F_{d1 - 2}} + \frac{{1 - \varepsilon }}{\varepsilon }\int_{{A_2}} {{{q''}_{2,rad}}\left( {{x_2}} \right)d{F_{_{d1 - d2}}}\left( {{x_1},{x_2}} \right)} \qquad \qquad(1)$
$\frac{{{{q''}_{2,rad}}\left( {{x_2}} \right)}}{\varepsilon } = \sigma T_2^4({x_2}) - \sigma T_1^4({x_1}){F_{d2 - 1}} + \frac{{1 - \varepsilon }}{\varepsilon }\int_{{A_1}} {{{q''}_{1,rad}}\left( {{x_1}} \right)d{F_{_{d2 - d1}}}\left( {{x_2},{x_1}} \right)} \qquad \qquad(2)$

Note that the q"k,rad is the heat flux required to balance the radiation at the surface. It can be supplied by conduction, convection, internal generation, or a combination of these. In this case, only conduction is present. If conduction is assumed to occur only in the x-direction in the thin plates, then two additional energy equations are written for the conducting plates

${q''_{1,rad}} = {q''_e} + k\delta \frac{{{d^2}{T_1}}}{{dx_1^2}}\qquad \qquad(3)$
${q''_{2,rad}} = k\delta \frac{{{d^2}{T_2}}}{{dx_2^2}}\qquad \qquad(4)$

Equations (3) and (4) contain second-order derivatives, so two boundary conditions are required for each. The appropriate ones are

$\frac{{d{T_1}}}{{d{x_1}}} = \frac{{d{T_2}}}{{d{x_2}}} = 0{\rm{ at }}{x_1},{x_2}{\rm{ }} = {\rm{ }}W/2\qquad \qquad(5)$
$k\frac{{d{T_1}}}{{d{x_1}}} = \varepsilon \sigma T_1^4({x_1} = 0);{\rm{ }}k\frac{{d{T_2}}}{{d{x_2}}} = \varepsilon \sigma T_2^4({x_2} = 0)\qquad \qquad(6)$

Equation (5) is a symmetry condition, and eq. (6) equates the conduction into the plate ends with the radiative loss from the plate end.

Equations (3) and (4) can be substituted into eqs. (1) and (2), resulting in two simultaneous integral-differential equations with the unknowns T1(x1) and T2(x2). The resulting Equations (7) are highly nonlinear, involving temperatures to the fourth power as well as second derivatives in temperature.

$\begin{array}{l} \frac{1}{\varepsilon }\left( {{q_e} + k\delta \frac{{{d^2}{T_1}}}{{dx_1^2}}} \right) = \sigma {T_1}^4\left( {{x_1}} \right) - \int_0^W {\left[ {\sigma T_2^4\left( {{x_2}} \right) - \left( {\frac{{1 - \varepsilon }}{\varepsilon }} \right)k\delta \frac{{{d^2}{T_2}}}{{dx_2^2}}} \right]d{F_{d1 - d2}}} \\ \frac{1}{\varepsilon }\left( {k\delta \frac{{{d^2}{T_2}}}{{dx_2^2}}} \right) = \sigma {T_2}^4\left( {{x_2}} \right) - \int_0^W {\left[ {\sigma T_1^4({x_1}) - \left( {\frac{{1 - \varepsilon }}{\varepsilon }} \right)\left( {{q_e} + k\delta \frac{{{d^2}{T_1}}}{{dx_1^2}}} \right)} \right]d{F_{d1 - d2}}} \\ \end{array}\qquad \qquad(7)$

To simplify eqs. (7), define

${T_{ref}} = {\left( {{q_e}/\sigma } \right)^{1/4}};{N_{CR}} = \frac{{k\delta }}{{\sigma T_{ref}^3{W^2}}},\theta = T/{T_{ref}};{X_1} = {x_1}/W;{X_2} = {x_2}/W\qquad \qquad(8)$

and eqs. (7) become

$\begin{array}{l} \left( {1 + {N_{CR}}\frac{{{d^2}{\theta _1}}}{{dX_1^2}}} \right) = \varepsilon {\theta _1}^4\left( {{X_1}} \right) - \int_0^1 {\left[ {\varepsilon \theta _2^4\left( {{X_2}} \right) - \left( {1 - \varepsilon } \right){N_{CR}}\frac{{{d^2}{\theta _2}}}{{dX_2^2}}} \right]d{F_{d1 - d2}}} \\ \left( {{N_{CR}}\frac{{{d^2}{\theta _2}}}{{dX_2^2}}} \right) = \varepsilon {\theta _2}^4\left( {{X_2}} \right) - \int_0^1 {\left[ {\varepsilon \theta _1^4({X_1}) - \left( {1 - \varepsilon } \right)\left( {1 + {N_{CR}}\frac{{{d^2}{\theta _1}}}{{dX_1^2}}} \right)} \right]d{F_{d1 - d2}}} \\ \end{array}\qquad \qquad(9)$

with boundary conditions [from eqs. (5) and (6)]

$\begin{array}{l} {\left. {\frac{{d{\theta _1}}}{{d{X_1}}}} \right|_{{X_1} = 1/2}} = {\left. {\frac{{d{\theta _2}}}{{d{X_2}}}} \right|_{{X_2} = 1/2}} = 0 \\ {N_{CR}}{\left. {\frac{{d{\theta _1}}}{{d{X_1}}}} \right|_{{X_1} = 0}} = \varepsilon \theta _1^4\left( {{X_1} = 0} \right)\left( {\frac{\delta }{W}} \right){\rm{ }} \\ {N_{CR}}{\left. {\frac{{d{\theta _2}}}{{d{X_2}}}} \right|_{{X_2} = 0}} = \varepsilon \theta _2^4\left( {{X_2} = 0} \right)\left( {\frac{\delta }{W}} \right) \\ \end{array}\qquad \qquad(10)$

In terms of the nondimensional variables, the required configuration factor is (using h = H / W),

$d{F_{d{X_1} - d{X_2}}} = \frac{{{h^2}d{X_2}}}{{2{{\left[ {{{\left( {{X_2} - {X_1}} \right)}^2} + {h^2}} \right]}^{3/2}}}}\qquad \qquad(11)$

The parameters that appear in the problem are the emissivity ε, the conduction/radiation parameter NCR, and the two geometric parameters δ / W and h. The NCR or modifications of it often appear in radiation-conduction problems. Its meaning can be appreciated by multiplying the numerator and denominator by Tref to give

${N_{CR}} = \frac{{k\delta {T_{ref}}}}{{\sigma T_{ref}^4{W^2}}} = \frac{{k\delta \left( {{T_{ref}}/W} \right)}}{{\sigma T_{ref}^4W}}\qquad \qquad(12)$

The last form in eq. (12) is the ratio of a conduction heat transfer to a radiative emission. The conduction is through the plate thickness σ with a temperature gradient of (Tref0) / W. The radiation is emission from a black surface of width W and temperature Tref. The NCR is thus a measure of the importance of conduction to radiation in affecting the temperature distributions on the surfaces. For example, if NCR = 0, eqs. (9) reduce to the radiation-only solutions. If NCR is very large, then the problem becomes a pure conduction problem. (This is shown be reverting to eqs. (3) and (4) and placing them in dimensionless form). It is unlikely that eqs. (9), subject to boundary conditions eqs. (10), will have a closed-form analytical solution, and a numerical approach must be employed.

## Numerical Method

To solve the radiation-conduction problem defined in eqs. (9-11), the integrals in the equations are replaced by summations over elements of width Δx, i.e.,

$\int_0^{N \cdot \Delta {x_n}} {f\left( {{x_m},{x_n}} \right)} d{x_n} \approx \sum\limits_{n = 1}^N {{f_{m,n}}} \Delta {x_n}\qquad \qquad(13)$

Substituting into eq. (9) and using the finite difference form for the second derivatives results in

$\begin{array}{l} \left( {1 + {N_{CR}}\frac{{{\theta _{1,m - 1}} - 2{\theta _{1,m}} + {\theta _{1,m + 1}}}}{{{{\left( {\Delta {X_1}} \right)}^2}}}} \right) = \varepsilon {\theta _{1,m}}^4 - \sum\limits_{n = 1}^N {{\Phi _{1,m,n}}} \Delta {X_2} \\ \left( {{N_{CR}}\frac{{{\theta _{2,n - 1}} - 2{\theta _{2,n}} + {\theta _{2,n + 1}}}}{{{{\left( {\Delta {X_2}} \right)}^2}}}} \right) = \varepsilon {\theta _{2,n}}^4 - \sum\limits_{m = 1}^M {{\Phi _{2,m,n}}} \Delta {X_1} \\ \end{array}\qquad \qquad(14)$

where, including the configuration factor relations, the Φs are

$\begin{array}{l} {\Phi _{1,m,n}} = \left[ {\varepsilon \theta _{2,n}^4 - \left( {1 - \varepsilon } \right){N_{CR}}\frac{{{\theta _{2,n - 1}} - 2{\theta _{2,n}} + {\theta _{2,n + 1}}}}{{{{\left( {\Delta {X_2}} \right)}^2}}}} \right] \\ {\rm{ }} \times \frac{{{h^2}}}{{2{{\left[ {{{\left( {n \cdot \Delta {X_2} - m \cdot \Delta {X_1}} \right)}^2} + {h^2}} \right]}^{3/2}}}} \\ {\Phi _{2,m,n}} = \left[ {\varepsilon \theta _{1,m}^4 - \left( {1 - \varepsilon } \right)\left( {1 + {N_{CR}}\frac{{{\theta _{1,m - 1}} - 2{\theta _{1,m}} + {\theta _{1,m + 1}}}}{{{{\left( {\Delta {X_1}} \right)}^2}}}} \right)} \right] \\ {\rm{ }} \times \frac{{{h^2}}}{{2{{\left[ {{{\left( {m \cdot \Delta {X_1} - n \cdot \Delta {X_2}} \right)}^2} + {h^2}} \right]}^{3/2}}}} \\ \end{array}\qquad \qquad(15)$

The boundary conditions become

$\begin{array}{l} {\theta _{1,(\frac{N}{2} - 1)}} = {\theta _{1,(\frac{N}{2} + 1)}};{\theta _{2,(\frac{M}{2} - 1)}} = {\theta _{2,(\frac{M}{2} + 1)}};M{\rm{ and }}N{\rm{ even}} \\ {N_{CR}}\frac{{{\theta _{1,2}} - {\theta _{1,1}}}}{{{{\left( {\Delta {X_1}} \right)}^2}}} = \varepsilon \theta _{1,1}^4\left[ {\left( {\frac{{\delta /W}}{{\Delta {X_1}}}} \right) + 1} \right] - 1 - \sum\limits_{n = 1}^N {{\Phi _{1,1,n}}} \Delta {X_2} \\ {N_{CR}}\frac{{{\theta _{2,2}} - {\theta _{2,1}}}}{{{{\left( {\Delta {X_2}} \right)}^{\rm{2}}}}} = \varepsilon \theta _{2,1}^4\left[ {\left( {\frac{{\delta /W}}{{\Delta {X_2}}}} \right) + 1} \right] - \sum\limits_{m = 1}^M {{\Phi _{2,m,1}}} \Delta {X_1} \\ \end{array}\qquad \qquad(16)$

Because the boundary elements at n = 1 and m = 1 are finite with width = ΔX, the radiative exchange to the ΔX surfaces must be included in addition to the radiation-conduction balance at the end surfaces. The ΔX approaches zero in the analytical boundary conditions of eq. (10) so these terms do not appear. The result is (M + N) nonlinear algebraic equations with N values of θ1,n and M values of θ2,m. This set can, in principle, be solved by solvers available for treating a set of non-linear algebraic equations. However, experience indicates that such solvers often fail unless some further steps are taken to rearrange the form of the equation.

Let ΔX1 = ΔX2 (equal sized increments on each plate, so that eqs. (14) - (16) become

$\begin{array}{l} \left( {1 + {N_{CR}}\frac{{{\theta _{1,m - 1}} - 2{\theta _{1,m}} + {\theta _{1,m + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right) = \varepsilon {\theta _{1,m}}^4 - \sum\limits_{n = 1}^N {{\Phi _{1,m,n}}} \Delta X \\ \left( {{N_{CR}}\frac{{{\theta _{2,n - 1}} - 2{\theta _{2,n}} + {\theta _{2,n + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right) = \varepsilon {\theta _{2,n}}^4 - \sum\limits_{m = 1}^M {{\Phi _{2,m,n}}} \Delta X \\ \end{array}\qquad \qquad(17)$
$\begin{array}{l} {\Phi _{1,m,n}} = \left[ {\varepsilon \theta _{2,n}^4 - \left( {1 - \varepsilon } \right){N_{CR}}\frac{{{\theta _{2,n - 1}} - 2{\theta _{2,n}} + {\theta _{2,n + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right] \\ {\rm{ }} \times \frac{{{h^2}}}{{2{{\left[ {{{\left( {n - m} \right)}^2}\Delta {X^2} + {h^2}} \right]}^{3/2}}}} \\ {\Phi _{2,m,n}} = \left[ {\varepsilon \theta _{1,m}^4 - \left( {1 - \varepsilon } \right)\left( {1 + {N_{CR}}\frac{{{\theta _{1,m - 1}} - 2{\theta _{1,m}} + {\theta _{1,m + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right)} \right] \\ {\rm{ }} \times \frac{{{h^2}}}{{2{{\left[ {{{\left( {m - n} \right)}^2}\Delta {X^2} + {h^2}} \right]}^{3/2}}}} \\ \end{array}\qquad \qquad(18)$
$\begin{array}{l} {\theta _{1,(\frac{N}{2} - 1)}} = {\theta _{1,(\frac{N}{2} + 1)}};{\theta _{2,(\frac{M}{2} - 1)}} = {\theta _{2,(\frac{M}{2} + 1)}};M{\rm{ and }}N{\rm{ even}} \\ {N_{CR}}\frac{{{\theta _{1,2}} - {\theta _{1,1}}}}{{{{\left( {\Delta X} \right)}^2}}} = \varepsilon \theta _{1,1}^4\left[ {\left( {\frac{{\delta /W}}{{\Delta X}}} \right) + 1} \right] - 1 - \sum\limits_{n = 1}^N {{\Phi _{1,1,n}}} \Delta X \\ {N_{CR}}\frac{{{\theta _{2,2}} - {\theta _{2,1}}}}{{{{\left( {\Delta X} \right)}^{\rm{2}}}}} = \varepsilon \theta _{2,1}^4\left[ {\left( {\frac{{\delta /W}}{{\Delta X}}} \right) + 1} \right] - \sum\limits_{m = 1}^M {{\Phi _{2,m,1}}} \Delta X \\ \end{array}\qquad \qquad(19)$

The parameter NCR / (ΔX)2 now appears in eqs. (12) - (19), and is a guideline for how to approach a numerical solution.

## Conduction Dominated Problems

If conduction dominates (i.e., NCR/(ΔX)2>>1), then eq. (14) can be rearranged to reflect a conduction problem that is perturbed by the presence of radiation:

$\begin{array}{l} {\theta _{1,m - 1}} - 2{\theta _{1,m}} + {\theta _{1,m + 1}} = \frac{{{{\left( {\Delta X} \right)}^2}}}{{{N_{CR}}}}\left[ {\varepsilon {\theta _{1,m}}^4 - 1 - \sum\limits_{n = 1}^N {{\Phi _{1,m,n}}} \Delta X} \right] \\ {\theta _{2,n - 1}} - 2{\theta _{2,n}} + {\theta _{2,n + 1}} = \frac{{{{\left( {\Delta X} \right)}^2}}}{{{N_{CR}}}}\left[ {\varepsilon {\theta _{2,n}}^4 - 1 - \sum\limits_{m = 1}^M {{\Phi _{2,m,n}}} \Delta X} \right] \\ \end{array}\qquad \qquad(20)$

and for this case the right-hand-side of eq. (20) will be small and the solution will be close to the conduction-only solution. Equation (20) (including the insertion of the boundary conditions) is a set of equations of the matrix form

${\mathbf{A\theta (x) = C(x)}}\qquad \qquad(21)$

where A is a tridiagonal matrix of coefficients that need only be inverted once. The solution is then

${\mathbf{\theta (x) = }}{{\mathbf{A}}^{{\mathbf{ - 1}}}}{\mathbf{C(x)}}\qquad \qquad(22)$

Equation (22) is solved by assuming the distributions of θ1,n and θ2,m, using these to evaluate Φ1,m,n and Φ2,m,n, which in turn are used to evaluate C(x). The matrix multiplication indicated in eq. (22) gives a new set of θ values; the process is repeated until convergence. Effectively, this method is solving a set of linear equations by interation in the unknown values of θ, and amounts to a series of matrix multiplications.

## Radiation dominated problems

If radiation dominates so that NCR / (ΔX)2<<1, then eq. (14) can be rearranged to

$\begin{array}{l} {\theta _{1,m}}^4 = \frac{1}{\varepsilon }\left[ {\sum\limits_{n = 1}^N {{\Phi _{1,m,n}}} \Delta X - \left( {1 + {N_{CR}}\frac{{{\theta _{1,m - 1}} - 2{\theta _{1,m}} + {\theta _{1,m + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right)} \right] \\ {\theta _{2,n}}^4 = \frac{1}{\varepsilon }\left[ {\sum\limits_{m = 1}^M {{\Phi _{2,m,n}}} \Delta X - \left( {{N_{CR}}\frac{{{\theta _{2,n - 1}} - 2{\theta _{2,n}} + {\theta _{2,n + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right)} \right] \\ \end{array}\qquad \qquad(23)$

Again, the set of equations can be arranged as a matrix equation of the form

${{\mathbf{A}}_{\mathbf{1}}} \cdot {{\mathbf{\theta }}^{\mathbf{4}}}{\mathbf{ = }}{{\mathbf{C}}_{\mathbf{1}}}{\mathbf{(\theta )}}\qquad \qquad(24)$

Matrix inversion of A1 gives

${{\mathbf{\theta }}^{\mathbf{4}}}{\mathbf{ = }}{\left( {{{\mathbf{A}}_{\mathbf{1}}}} \right)^{{\mathbf{ - 1}}}}{{\mathbf{C}}_{\mathbf{1}}}{\mathbf{(\theta )}}\qquad \qquad(25)$

Again, the A1 matrix need be inverted only once. Equation (25) is solved by assuming the distributions of θ1,n and θ2,m, using these to evaluate Φ1,m,n and Φ2,m,n and the second derivative terms, which in turn are used to evaluate C1(x). The matrix multiplication indicated in eq. (25) gives a new set of θ values; the process is repeated until convergence. Effectively, this method is iteravely solving a set of linear equations in the unknown values of θ4.

## Problems with Both Modes Significant

When NCR is not near either the large or small limit, then the problem is truly nonlinear, and the solution methods described for small or large NCR / (ΔX)2 values often fail. In that case, eq. (14) can be arranged as

$\begin{array}{l} \left( {\varepsilon {\theta _{1,m}}^4 - {N_{CR}}\frac{{{\theta _{1,m - 1}} - 2{\theta _{1,m}} + {\theta _{1,m + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right) = 1 + \sum\limits_{n = 1}^N {{\Phi _{1,m,n}}} \Delta X \\ \left( {\varepsilon {\theta _{2,n}}^4 - {N_{CR}}\frac{{{\theta _{2,n - 1}} - 2{\theta _{2,n}} + {\theta _{2,n + 1}}}}{{{{\left( {\Delta X} \right)}^2}}}} \right) = \sum\limits_{m = 1}^M {{\Phi _{2,m,n}}} \Delta X \\ \end{array}\qquad \qquad(26)$

This equation in matrix form becomes

${{\mathbf{A}}_{\mathbf{2}}}{{\mathbf{\theta }}^{\mathbf{4}}}{\mathbf{ + B\theta = }}{{\mathbf{C}}_{\mathbf{2}}}\left( {\mathbf{\theta }} \right)\qquad \qquad(27)$

We could define ${{\mathbf{A}}_{\mathbf{3}}}\left( {\mathbf{\theta }} \right){\mathbf{ = }}{{\mathbf{A}}_{\mathbf{2}}} \cdot {{\mathbf{\theta }}^{\mathbf{3}}}{\mathbf{ + B}}$ so that eq. (27) becomes

${{\mathbf{A}}_{\mathbf{3}}}{\mathbf{(\theta ) \times \theta = }}{{\mathbf{C}}_{\mathbf{2}}}{\mathbf{(\theta )}}\qquad \qquad(28)$

If this approach is followed, an initial guess for the distributions of θ1,n and θ2,m is used to evaluate A3(θ) and C2(θ), A3(θ) is inverted, and a new set of θ values is found. This method has two problems: First, the A3(θ) matrix must be inverted at each iteration, and A3(θ) is probably a full matrix that can require significant time to invert. Second, this method is often quite unstable, and requires a small relaxation factor to be imposed between iterations to avoid divergence.

To avoid this, a modified Newton-Raphson iteration method can be used. In this case, eq. (27) is rewritten as

${g_j} = \sum\limits_{k = 1}^{M + N} {[{A_{jk}}} \theta _k^{4\left( 0 \right)} + {B_{jk}}\theta _k^{(0)}] - {C_j} \equiv {\rm{ residual}}\qquad \qquad(29)$

The residual is a measure of convergence of the solution, and will approach gj = 0 when solution is complete

Next, the function gjk is found

${g_{jk}} = 4{A_{jk}}{\left[ {\theta _k^{(0)}} \right]^3} + {B_{jk}}{\rm{ }}\left[ { \approx \frac{{\partial (LHS{\rm{ of Eq}}{\rm{.(10}}{\rm{.96}})}}{{\partial {\theta _k}}}} \right]\qquad \qquad(30)$

and gjk is seen to be the gradient in the residual. Now solve the auxiliary equation

$\left[ {{g_{jk}}} \right]\left[ {{\lambda _k}} \right] + \left[ {{g_j}} \right] = 0;{\rm{ }}\left[ {{\lambda _k}} \right] = - {\left[ {{g_{jk}}} \right]^{ - 1}}\left[ {{g_j}} \right]\qquad \qquad(31)$

The next iterative value for θk is found from

$\theta _k^{\left( p \right)} = \theta _k^{\left( {p - 1} \right)} + {\lambda _k}\qquad \qquad(32)$

The values of θkp are now used to find a new residual, and the process is repeated until convergence. The λk values are seen to be those that would make the residual equal to zero on the next iteration if the problem were linear. This approach has been found to be quite stable, although sometimes a relaxation factor is needed on θk between iterations.

Some of the complexities of multi-mode heat transfer when radiation is present have been outlined. A detailed discussion of more advanced numerical techniques is in Hogan and Gartling (2007). There are commercial software packages that can be used for these problems, and they can be particularly useful when convective transfer must be linked to radiation.

## References

Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.

Hogan, R.E. and Gartling, D.K., 2007, “Solution Strategies for Coupled Conduction/Radiation Problems,” Commun. Numer. Meth. Engng,, Vol. 23, pp. 523-542.