3.6. Radiation

3.6.1. Surface Radiation

A surface may exchange energy with its surroundings through thermal radiation. Any incident surface radiation will be either transmitted, reflected or absorbed. Letting $\tau$ , $\rho$ and $\alpha$ represent the fractions of the incident flux in each category then

$1 = \tau + \rho + \alpha$

and for no transmission

$1 = \rho + \alpha \text{ or } \rho = 1 - \alpha \;.$

Using the Kirchhoff identity

$\alpha = \epsilon$

then the reflectance is

(3.100) $\rho = 1 - \epsilon$

where $\epsilon$ is the emissivity.

In order to understand the radiative energy balance at a surface one considers the rate at which energy streams away from the surface, the radiosity, defined as

(3.101) $J = \epsilon E_{b} + \rho G$

where $E_{b}$ is the blackbody emissive power and $G$ is the irradiation. Substituting for the reflectance (3.100) then

(3.102) $J = \epsilon E_{b} + (1 - \epsilon) G \;\;.$

The surface flux $q$ is the difference between the energy that radiates away and the incident energy

(3.103) q = J - G

and substitution for the radiosity we find that

$q = \epsilon ( \sigma T^{4} - G )\;\;.$

When $G$ is derived from an external temperature interaction this boundary condition is often called far-field radiation, since it usually models the radiative transfer of energy between a surface and the external environment. However, the boundary condition is found to have more general utility when one considers its role in modeling radiative transfer between two surfaces, $1$ & $2$ , as shown in Fig. 3.3, where $A_{1}$ is analogous to $\partial \Omega_i$ .

For the case in which the temperature $T_{2}$ is known and independent of temperature $T_{1}$ then using the emissive power $\sigma T^{4}$ the normal flux per unit area across $A_{1}$ may be written as

(3.104) $q = \sigma \epsilon F \left( T_{1}^4 - T_{2}^4 \right),$

where $\sigma$ denotes the Stefan-Boltzmann constant, $\epsilon$ is the emissivity of the surface $A_{1}$ and $F$ is the form factor. We remark that (3.104) is a nonlinear boundary condition, since the unknown temperature, $T_{1}$ is raised to the fourth power and furthermore the emissivity may be a function of temperature. It is important to note that the form factor $F$ may differ from the more familiar view factor $F_{12}$ encountered in enclosure radiation problems. The question often asked is how does one determine the appropriate value of form factor $F$ .

The form factor $F$ can be best described using a network analogy of radiative transfer between two surfaces as shown below.

Fig. 3.4 Radiative Transfer Circuit Model

Using the emitted energy $E$ and radiosity $J$ the network the heat flux can be written in terms of a thermal resistance $R$ as

(3.105) $q = \sigma \frac{T_{1}^{4} - T_{2}^{4}}{R}$

and from the network model

$R = \frac{1-\epsilon_{1}}{ \epsilon_{1}} + \frac{1}{F_{12}} + \frac{ 1-\epsilon_{2}}{ \epsilon_{2}} \frac{A_{1}}{A_{2}} = \left( \frac{1}{\epsilon_{1}} - 1 \right) + \frac{1}{F_{12}} + \left( \frac{1}{\epsilon_{2}} - 1 \right) \frac{A_{1}}{A_{2}} \;\;.$

For $A_{2} \gg A_{1}$ the third term of $R$ can be neglected and $F_{12} = 1$ . Comparing expressions (3.104) and (3.105) then

$R = \frac{1}{\epsilon_{1}} \text{ and } F = 1 \;\;.$

so estimation of $F_{12}$ is not required.

For $\epsilon_{2} = 1$ (black receiving surface) but $A_{2}\; \text{ not } \gg A_{1}$ and once again the third term of $R$ can be neglected so that

(3.106) $R = \frac{ 1 - \epsilon_{1} }{\epsilon_{1}} + \frac{1}{F_{12}} = \frac{(1 - \epsilon_{1}) F_{12} + \epsilon_{1} }{\epsilon_{1}F_{12}} \text{ and } F = \frac{F_{12}}{( 1 - \epsilon_{1} ) F_{12} + \epsilon_{1}} \;\;.$

If both surfaces are black $\epsilon_{1} = \epsilon_{2} = 1$ then from the previous expression (3.106) we find that $F = F_{12}$ .

During a simulation the surface heat flux is integrated over the spatial discretization of surface $A_{1}$ . Here we note that defining the flux on a per unit area basis enables us to apply the radiative flux (3.104) consistently with $F$ evaluated for the entire surface $A_{1}$ even when the surface is discretized.

3.6.2. Enclosure Radiation

When energy radiates from one portion of a surface to another, and the intermediate medium is transparent (i.e., it does not absorb any energy), then enclosure radiation may be used to model the heat flux on the surface. Using the net radiation method [18], the normal flux at a particular location on the surface may be written as the difference between the emitted radiative heat flux leaving the surface, and the absorbed incident radiative flux due to the rest of the enclosure, namely

(3.107) $q_n = \sigma \epsilon \vartemp^4 - \alpha G,$

where $\alpha$ denotes the absorptivity of the surface, and $G$ represents the surface irradiation. Under the additional assumption that the emissivity, absorptivity, and reflectivity are independent of direction and wavelength, we may write

(3.108) $\epsilon = \alpha = 1 - \rho,$

where $\rho$ is the reflectivity.

Without loss of generality, we can regard the enclosure, $\Gamma_{\mathcal{E}}$ , as a union of $E$ surfaces,

$\Gamma_{\mathcal{E}} = \Gamma_1 \cup \Gamma_2, \ldots \Gamma_{E-1}\cup\Gamma_E$

This situation is illustrated in Fig. 3.6, where the radiosity for surface $i$ in the enclosure is defined to be

(3.109) $J_i = \sigma \epsilon_i \vartemp^4_i + \rho_i G_i,$

where $\vartemp_i$ is the area averaged constant temperature on facet $\Gamma_i$ computed via

(3.110) $\vartemp_i = \frac{ \sum \limits_{j=1}^N \vartemp_j \int_{\Gamma_i} \psi_j \,d\Gamma} {\int_{\Gamma_i} d\Gamma }$

The surface irradiation for surface $i$ is determined by the radiosity of all the other surfaces in the enclosure through the relation

(3.111) $G_i = \sum\limits_{j = 1}^E {F_{ij} J_j},$

where $F_{ij}$ denotes the geometric viewfactor of surface $i$ with respect to surface $j$ . The viewfactor may be considered as the fraction of energy that leaves surface $j$ and arrives at surface $i$ . For a given pair of surfaces shown in Fig. 3.5, the viewfactor is defined as the following integral

(3.112) $F_{ij} = \frac{1}{A_i} \int_{A_{i}} \int_{A_{j}} \frac{ cos \theta_{i} \cos \theta_{j}}{\pi r^{2}} \;\; \delta_{ij} dA_{i} dA_{j}\;\;.$

Viewfactors are computed using the Chaparral library [19]. Determination of the viewfactors is a compute intensive endeavor. As such, extraneous calculations are eliminated based upon the geometry. One example of this would be excluding this calculation for surfaces which are not visible to each other. Moreover, from a geometric view of enclosure surfaces, Fig. 3.5, one can conclude that legitimate interactions between surfaces are those for which $n_{i}$ and $n_{j}$ are opposed. Thus an important feature of the enclosure model is the notion of inward facing normals. This convention effectively defines the interaction between the enclosure subfacets.

Note

For closed surfaces (watertight enclosures), for each facet $i$ the row-sum over all surface $N$ facets should equal one, i.e., $\sum_{j=1}^N F_{ij}= 1.0$ .

Surface to surface interaction in enclosure radiation — Fig. 3.6 Two arbitrary facets radiating energy to one another in a *radiation enclosure*. The energy exchanged depends on: the shape, orientation, distance, area $A_i$ , $A_j$ , temperatures $\vartemp_i$ , $\vartemp_j$ , and radiative properties of the facets $\epsilon_i$ , $\epsilon_j$ .

Fig. 3.6 Two arbitrary facets radiating energy to one another in a *radiation enclosure*. The energy exchanged depends on: the shape, orientation, distance, area $A_i$ , $A_j$ , temperatures $\vartemp_i$ , $\vartemp_j$ , and radiative properties of the facets $\epsilon_i$ , $\epsilon_j$ .

Upon substitution of equations (3.111) and (3.108) into (3.109), the radiosity may be written as

(3.113) $J_i - (1 - \epsilon_i) \sum\limits_{j = 1}^N {F_{ij} J_j} = \sigma \epsilon_i \vartemp^4_i,$

Finally, the first $J_i$ term in (3.113) may be moved inside the summation to yield

(3.114) $\sum\limits_{j = 1}^N {\left[ \delta_{ij} - (1 - \epsilon_i) F_{ij} \right] J_j} = \sigma \epsilon_i \vartemp^4_i,$

where

$\delta_{ij} = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j \end{cases}$

(3.114) is a system of equations that can be solved for the unknown radiosity at each face. Finally, we may rewrite (3.107) to express the radiative heat flux on surface $i$ as

(3.115) $q_n = \sigma \epsilon_i \vartemp_i^4 - \epsilon_i G_i,$

where $G_i$ is given by (3.111) and depends on the unknown radiosity. This heat flux is applied as a boundary condition to the appropriate energy equations, thereby fully coupling the enclosure radiation problem with the volume energy equations.

Instead of solving the monolithic equation system by forming the jacobian for the combined coupled system of equation systems, a segregated approach is used. At each nonlinear iteration, the radiosity equation is solved separately with the right hand side computed from the average facet temperature (3.110) using the current nonlinear temperature solution. Subsequently, the irradiation is computed from the radiosity in order to compute the heat flux contribution given by (3.115). A Newton step is then performed where the jacobian contribution of the radiative flux is computed as

(3.116) $\dfrac{\partial R_i}{\partial T_j} = \int \phi_i \dfrac{\partial q_n}{\partial T_j} d\Gamma,$

where $R_i$ is the $i^{th}$ nonlinear temperature residual, $\phi_i$ the test function and where

(3.117) $\dfrac{\partial q_n}{\partial T_j} = \dfrac{\partial q_n}{\partial T^e}\dfrac{\partial T^e}{\partial T_j}.$

Here, the facet temperature $T_i$ given by (3.110) has been relabeled to $T^e$ to avoid confusion with the nodal finite element temperature solution. The sensitivity of the radiative heat flux with respect to the facet temperature is approximated as

(3.118) $\dfrac{\partial q_n}{\partial T^e} = 4 \sigma \epsilon^e \left[T^e\right]^3$

while the sensitivity of the facet temperature with respect to the nodal temperature dof is

(3.119) $P^e_j = \dfrac{\partial T^e}{\partial T_j} = \dfrac{\int_{\Gamma^e} \phi_j d\Gamma}{A^e},$

where $A^e$ is the area of a facet. In terms of the element contribution, this can be compactly written as

(3.120) $\dfrac{\partial R_i^e}{\partial T_j} = 4 \sigma \epsilon^e \left[T^e\right]^3 A^e P^e_i P^e_j$

Note

The average facet temperature is $T_e = \sum_j P^e_j T_j$ . For computing the cube $\left[T^e\right]^3$ , experience shows projecting the raised power i.e., $\left[T^e\right]^3 = \sum_j P^e_j \left[T_j\right]^3$ leads to less spurious oscillations.

3.6.3. Banded Wavelength Enclosure Radiation

The enclosure radiation model previously presented considers a net response over all wavelengths assuming that all surfaces behave as grey bodies. Here the surface response with regard to different wavelengths is implicitly included in the model by using temperature dependent emissivities and allowing the surfaces to emit as grey bodies. In this case the surface flux previously mentioned (3.107) can be alternatively expressed as

(3.121) $q_n = \int_{0}^{\infty} q(\lambda) \;\; d \lambda .$

In many applications engineering materials respond differently to different portions of the thermal energy spectrum. Thus modeling the thermal radiation response of these materials using the entire blackbody spectrum results in poor characterization of the enclosure response. For this purpose more specialized approaches for radiation modeling of wavelength dependent surfaces have been developed using continuous representations of emissivity as a function of wavelength and temperature while integrating over the wavelength spectrum. For numerical modeling the simplest approach involves discretization of the wavelength spectrum into a few representative bands $N_{b}$ and integrating over each of the $k$ bands. Using this approach the methods previously described follow directly except that the emissivity and emissive power now have independent representations in a wavelength band, each of which contributes some surface flux $\Delta q_{k}$ thus

(3.122) $q_n = \sum\limits_{k=1}^{N_{b}} \Delta q_{k \lambda} .$

Recalling that flux can also be expressed in terms of radiosity then similarly the radiosities are obtained by solving a system of equations

(3.123) $\sum\limits_{k=1}^{N_{b}} { \sum\limits_{j = 1}^N {\left[ \delta_{(ik)j} - (1 - \epsilon_{(ik)}) F_{ij} \right] \Delta J_{jk} } } \Delta \lambda_{k} = \sum\limits_{k=1}^{N_{b}} { {\mbox{F}}_{\Delta \lambda_{k}} \sigma \epsilon_{(ik)} \left[\vartemp_i \right] ^4 } \Delta \lambda_{k}$

where $\Delta J_{jk}$ is the incremental radiosity and $\mbox{F}_{\Delta \lambda_{k}}$ is the blackbody fraction for band $k$ . Finally we note that given a fixed facet temperature field this solution can be carried out independently for each wavelength band and the radiosities $\Delta J_{jk}$ can be accumulated to obtain the net facet radiosity $J_{j}$ and net flux $q_{n}$ . Thus the banded wavelength model very much resembles the more simplified enclosure radiation model. From a simulation point of view, the major difference being that the modeler must supply information describing the band discretization in addition to emissivity models for each band.

Construction of the banded wavelength model begins by first prescribing the emissivity variation for each surface of the enclosure as a function of wavelength $\lambda$ . The overall band discretization is then obtained by collective consideration of the wavelengths at which the emissivity changes on any of the given surfaces. Thus the number of bands for the enclosure will likely be more than those for any single surface and is demonstrated in Fig. 3.7 below.

Note

Each modeled band must have an emissivity model supplied for each surface, even when its emissivity is not changing.

Fig. 3.7 Emissivity bands for two-surface problem.