Radiation Modeling with Participating Media

In AcuSolve enclosure or surface to surface radiation models, effect of media between surfaces is ignored. This assumption is acceptable when you are dealing with lower temperature fluids. Nonetheless, while you have semi-transparent media like glass or high temperature gases like in flames, effect of media should be considered in heat transfer analysis. Two models are available in AcuSolve: A simple one equation P₁ model and more detailed but expensive discrete ordinates (DO) model.

Radiative Transfer Equation (RTE)

Radiative energy balance in a participating media is governed by the following integro-differential equation, known as the radiative transfer equation (RTE):

\underset{rate ​ ​ of change}{\underset{︸}{Ω \cdot \nabla I (r, Ω)}} + \underset{Absorption}{\underset{︸}{κ I (r, Ω)}} + \underset{Scattering loss}{\underset{︸}{σ_{s} I (r, Ω)}} = \underset{Emission}{\underset{︸}{f (r)}} + \underset{Scattering addition}{\underset{︸}{\frac{σ_{s}}{4 π} \int_{4 π} I (r, Ω^{'}) ϕ (Ω^{'}, Ω) d Ω^{'}}}

where

I

is the radiation intensity,

r

is the spatial vector,

Ω

is the unit directional vector,

Ω^{'}

is the scattering directional vector,

κ

is the absorption coefficient, and

σ_{s}

is the scattering coefficient,

f (r)

is an emission, n the refractive index, and

σ

is the Stefan-Boltzman constant (5.67 × 10^-8 W m^-2 K^-4). T is the temperature (K).

P1 radiation model
Discrete Ordinates (DO) model

P1 Radiation Model

The P₁ model is the lowest order P_N (spherical harmonics) type radiation model. The method reduces the five independent variables of the Radiative Transfer Equation (RTE) into a PDE that is relatively simple in comparison. The model is the most computationally efficient of the radiation models in AcuSolve, but it can lose accuracy, under certain conditions, for optically thin media. It performs best in scenarios where the radiative intensity is near isotropic.

Governing Equation

P₁ approximation and assumptions

The P₁ model is derived from the general P_N formulation (a spherical harmonic series expansion of the radiative intensity for the angular variable) by assuming that the series is limited to four terms and integrating over all solid angles. From the first harmonic in the series approximation, the divergence of the radiative flux ( $q$ ) can be derived by integrating the RTE over all solid angles as

\nabla \cdot q = κ (4 n^{2} σ T^{4} - G)

where $G$ is the incident radiation $\frac{σ_{s}}{4 π} \int_{4 π} I (r, Ω^{'}) ϕ (Ω^{'}, Ω) d Ω^{'}$ , $κ$ is the absorption coefficient, $n$ the refractive index, and $σ$ is the Stefan-Boltzman constant (5.67 × 10^-8 W m^-2 K^-4), $q$ is the radiative flux.

Additionally, a second vector equation can be derived from the other three harmonic terms for the radiative flux

q = Γ ​ \nabla G

where $Γ$ is the diffusion coefficient.

By taking the divergence of (2), and substituting into this the right hand side of equation (1), leads to elimination of the heat flux. The final diffusion reaction equation describing the transport of incident radiation is given by

Γ \nabla^{2} G - κ (4 n^{2} σ T^{4} - G) = 0

where $Γ$ is a diffusion coefficient given by

Γ = \frac{1}{3 (κ + σ_{s}) - A_{1} σ_{s}}

where $κ$ is the absorption coefficient, $σ_{s}$ is the scattering coefficient, and $A_{1}$ is the linear-anisotropic phase function.

Anisotropic scattering

The implementation in AcuSolve includes the ability to model linear anisotropic scattering

ϕ (Ω^{'} \cdot Ω) = 1 + A_{1} Ω^{'} \cdot Ω

where

Ω

is the unit directional vector in the scattering direction,

Ω^{'}

is the unit directional vector in the incident radiation direction. This term is included in the scattering term of the RTE along with the four term spherical harmonic expansion of the radiative intensity field (the first term being isotropic and the other three anisotropic). The final form of the P₁ approximation in equation (3) includes this term. The values of the phase function A₁ have the following meaning:

A₁ = 1: More radiation is scattered in the forward direction
A₁= -1: More radiation is scattered in the backward direction
A₁= 0: Isotropic scattering

Coupling to energy equation

The source term in equation (1) can be substituted into the energy equation as a negative source since a local increase in radiative heat flux is due to a local decrease in thermal energy.

In AcuSolve the default stagger sequence when the P₁ radiative heat transfer solver is enabled is:

Solve the energy equation with source term ( κ(4n²σT⁴-G) ), where G is zero for the time step
Pass temperature to radiation solver and solve equation (3)
Repeat until converged

Boundary Conditions

Marshak's boundary condition

Based on the assumption that the walls are diffused gray surfaces, that is, independent of wavelength, the appropriate wall boundary condition is

n \cdot \nabla G = \frac{ε}{2 (2 - ε)} (4 n^{2} σ T_{w}^{4} - G_{w})

where $ε$ is the surface emissivity, $T_{w}^{}$ is the wall temperature, and $G_{w}$ the wall incident radiation. The boundary radiative heat flux can be calculated from the incident radiation and the temperature at the wall:

q_{w} = \frac{ε}{2 (2 - ε)} (G_{w} - 4 n^{2} σ T_{w}^{4})

Discrete Ordinates (DO) Model

Governing equation

The governing equation is the radiative transfer equation limited to a finite number of directions (or ordinates)

Ω \cdot \nabla I (r, Ω) + κ I (r, Ω) + σ_{s} I (r, Ω) = f (r) + \frac{σ_{s}}{4 π} \int_{4 π} I (r, Ω^{'}) ϕ (Ω^{'}, Ω) d Ω^{'}

f (r) = κ I_{b} (r) = κ \frac{σ T^{4}}{π}

Scattering term (source term)

The integral over all the directions in equation (1) is replaced by a numerical quadrature for $N$ different ordinate directions ( $Ω_{i}$ )

S = \frac{σ_{s}}{4 π} \sum_{j = 1}^{N} w_{j} I (r, Ω_{j}) ϕ (Ω_{j}, Ω_{i})

The phase function, $ϕ (Ω_{j}, Ω_{i})$ , is given by

ϕ (Ω^{'} \cdot Ω) = 1 + A_{1} (Ω^{'} \cdot Ω)

Boundary conditions (RTE)

Diffused surface

If a surface emits and reflects diffusely, the exiting intensity is directionally independent and is given by

I (r_{w}, Ω_{i}) = ε (r_{w}) I_{b} (r_{w}) + \frac{1 - ε (r_{w})}{π} \sum_{n \cdot Ω_{j} < 0} w_{j} I (r_{w}, Ω_{j}) |n \cdot Ω_{j}|, n \cdot Ω_{j} > 0

Specular and diffuse surface

The diffused fraction defines the proportion of reflected radiation intensity at a surface which is diffused, that is, the reflection may also have a specular component. If the radiation intensity reflection coefficient at the surface is defined by

ρ = ρ^{S} + ρ^{D} = 1 - ε

then the diffused reflection coefficient, $ρ^{D}$ , is defined in terms of the diffused fraction $α$ and the emissivity of the surface by

ρ^{D} = α (1 - ε)

and the specular reflection coefficient

ρ^{S} = (1 - α) (1 - ε)

If $α$ =1, then the reflection at the surface is completely diffused. If $α$ =0 then the reflection is specular. The outgoing intensity, I, at the surface in terms of the above two reflection coefficients is given by

I (r_{w}, Ω_{i}) = ε (r_{w}) I_{b} (r_{w}) + \frac{ρ^{D} (r_{w})}{π} \sum_{n \cdot Ω_{j} < 0}^{N} w_{j} I (r_{w}, Ω_{j}) |n \cdot Ω_{j}| + ρ^{S} (r_{w}) I (r_{w}, Ω_{j}), n \cdot Ω_{j} > 0

where the first terms represent emission from the surface, the second term the diffused component incoming radiation heat flux and the third the specular component. The diffused component represents a sum over all radiation intensities along ordinates that are incident to the surface (that is, a hemisphere of incoming radiation to the surface); $n$ is the normal into the domain and $Ω_{j}$ the jth ordinate direction. The ordinate direction ( $Ω$ ), the total number of ordinate directions ( $N$ ) and the weights ( $w$ ) are automatically defined by the order of the radiation_quadrature (S2, S4, S6, S8 & S10). The specular ordinate direction ( $Ω_{S}$ ) is the direction that the radiation intensity must strike the surface to reflect in a specular fashion along the outgoing ordinate direction, $Ω_{i}$ , and is given by

Ω_{s} = Ω_{i} - 2 (Ω_{i} \cdot n) n

which means the angle that incident radiation intensity strikes the surface equals the angle of reflection.

Boundary conditions (Energy equation)

Interface and outflow/inflow boundary conditions

At an opaque interface between participating and non-participating media or outflows/inflows a radiative heat flux must be added to the boundaries in the energy equation. This flux is given by

Q_{r a d} = ε (\sum_{Ω_{j} \cdot n > 0} w_{j} I (r_{w}) |Ω_{j} \cdot n| - n^{2} σ T_{w}^{4})

For an opening, that is, outflow or inflow, the black body intensity used in the calculation of the outgoing intensity at the surface is based on the opening temperature of the surrounding:

I_{b} = \frac{κ n^{2} σ T_{o p e n}^{4}}{π}

while for an interface it is based on the current temperature solution.

Output metrics

Two directionally integrated output metrics can be derived from the radiative intensities: incident radiation and radiative heat flux.

Incident radiation

Incident radiation is the total intensity impinging on a point from all directions and is given by

G = \sum_{i = 1}^{N} w_{i} I (r)

where $I$ is the intensity in direction i, $N$ the number of ordinates, $w$ the weights.

Interface Between two Semitransparent Media

At the interface between two semitransparent media (referred to as medium 1 and medium 2 below), radiative intensity is both transmitted and reflected. The proportion of transmitted and reflected intensity at the interface depends on: the refractive indices (n₁, n₂) of the two media; the incident angle that radiative intensity strikes the surface; and the diffuse fraction of the surface.

Reflection and transmission for specular interfaces

Reflection at an interface is governed by the angle of incidence of a radiative intensity and the refractive indices of the two media. The image below shows the different refracted and reflected rays between two media.

Figure 1. Reflected and refracted directions at the interface between two participating media of different refractive indices (n₁ < n₂). For the medium of higher refractive index if the incoming direction is greater than the critical angle total internal reflection occurs (no transmission occurs across the interface). This is represented by the gray dashed lines in the image.

The cosine of the incident angle for the incoming ordinate is given by

\cos θ_{1} = Ω_{I}^{1} \cdot n

where $n$ is the outward facing normal direction at the interface (towards the second medium) and $Ω_{I}^{1}_{1}$ is the unit direction of incoming radiative intensity to the surface, given by

Ω_{I}^{1}_{1} = Ω_{R}^{1} - 2 (Ω_{R}^{1} \cdot n) n

where $Ω_{I}^{1}_{1}$ is the unit reflected ordinate direction vector and also represents the current ordinate direction being solved. The equivalent calculation can also be performed for medium two.

Radiative intensity that is transmitted into a second medium undergoes refraction governed by Snell's law,

n_{1} \sin θ_{1} = n_{2} \sin θ_{2}

where n₁ and n₂ are the refractive indices of mediums. θ₁ and θ₂ are the angles of incidence and refraction of radiative intensity relative to the interface normal, respectively. This can also be represented in vector form by

n_{1} (n \times Ω_{I}^{1}) = n_{2} (n \times Ω_{R}^{2})

The incoming direction vector in medium two for a ray refracted from medium two to one is given by

Ω_{I}^{2} = \frac{n_{1}}{n_{2}} Ω_{R}^{1} + (\frac{n_{1}}{n_{2}} \cos θ_{1} - \sqrt{1 - {(\frac{n_{1}}{n_{2}})}^{2} (1 - \cos θ_{1})}) n

The above expression is valid providing the expression under the radicand is greater than zero; otherwise total internal reflection occurs.

The actual reflected and refracted directions differ slightly from the calculated direction since these directions will unlikely coincide with a discrete ordinate direction. Since the number of directions is governed by the order of radiation quadrature, higher quadrature orders are more accurate for interface problems.

Figure 2. Octant of angular quadrature with transmission direction. The calculated transmission direction (Ω_T) is shifted to match the nearest quadrature point.

Depending on the refractive indices of the two media and the angle of incidence, θ₁, the proportion of radiation intensity that is reflected or transmitted will vary. If $n_{1} < n_{2}$ , then the radiative intensity in medium one will be partially reflected and partially transmitted into a cone defined by the critical angle, θ_c, which is given by:

θ_{c} = \sin^{- 1} \frac{n_{1}}{n_{2}} when n_{1} < n_{2}

and defines a cone in 3D (see the images below). The critical angle is defined by a ray that grazes the surface on the side of medium 2 and is transmitted exactly at the critical angle in medium 1. As an example, if

n_{1} > n_{2}

and

θ_{1} > θ_{c}

, the direction of incoming radiation is greater than the critical angle and total internal reflection will occur. This means the incoming ray is reflected at the same angle of incidence and no transmission occurs.

Figure 3. Total internal reflection of intensity rays at an interface ( $n_{1} > n_{2}$ ). The critical angle cone defines the angle outside which total internal reflection occurs, that is, $θ_{1} > θ_{c}$ .

However, if

θ_{1} < θ_{c}

, the outgoing intensity will include transmission from medium 2 to medium 1, as shown in the image below.

Figure 4. Reflection and transmission of intensity rays at an interface ( $n_{1} > n_{2}$ ) for $θ_{1} < θ_{c}$

At the interface, the intensity is partially reflected and transmitted into the other medium if $θ_{1} < θ_{c}$ or the outgoing direction is in medium 2. The reflected proportion from $Ω_{I}^{1} \to Ω_{R}^{1}$ , or reflectance, is given by

r_{12} = \frac{1}{2} {(\frac{n_{1} \cos θ_{1} - n_{2} \cos θ_{2}}{n_{1} \cos θ_{1} + n_{2} \cos θ_{2}})}^{2} + \frac{1}{2} {(\frac{n_{2} \cos θ_{1} - n_{1} \cos θ_{2}}{n_{2} \cos θ_{1} + n_{1} \cos θ_{2}})}^{2}

and the transmitted proportion from the $Ω_{I}^{2} \to Ω_{R}^{1}$ , or transmittance, is given by

τ_{21} = 1 - r_{12}

In the second medium, for the current scenario where $n_{2} > n_{1}$ , if $θ_{2} < θ_{c}$ the radiative intensity is, as for medium one, partially reflected and partially transmitted. The reflection coefficient is as described above since $r_{21} = r_{12}$ . If $θ_{2} > θ_{c}$ , then total internal reflection occurs and $r_{21} = 1.0$ and $τ_{12} = 0.0$ , meaning no transmission of radiative intensity into the second medium or from the first medium. This is shown in the image above with the gray dashed lines.

From the above, the outgoing radiative intensity on the side one of the interface for the current ordinate direction, $Ω_{R}^{1}$ is given by

I (Ω_{R}^{1}) = r_{12} I (Ω_{I}^{1}) + τ_{21} {(\frac{n_{1}}{n_{2}})}^{2} I (Ω_{I}^{2})

where for medium one, the first term on the right-hand side represents the reflected intensity in medium one and the second term represents the transmitted intensity from medium two to one. For medium two, if the current ordinate direction is $Ω_{R}^{2}$ then the intensity outgoing radiative intensity is given by

I (Ω_{R}^{2}) = r_{21} I (Ω_{I}^{1}) + τ_{12} {(\frac{n_{2}}{n_{1}})}^{2} I (Ω_{I}^{1})

For $n_{2} < n_{1}$ , the subscripts of the above analysis must be exchanged, and total internal reflection will now occur in medium one.

Reflection and transmission for diffuse interfaces

If the interface is diffused, for example, diffused_fraction = 1.0, the reflectivity of the interface is given by the hemispherically averaged reflectance:

r_{D, 12} = \frac{1}{2} + \frac{(n - 1) (3 n + 1)}{6 {(n + 1)}^{2}} - \frac{2 n^{3} (n^{2} + 2 n + 1)}{{(n^{4} - 1)}^{2} (n^{2} + 1)} + \frac{8 n^{4} (n^{4} + 1)}{{(n^{4} - 1)}^{2} (n^{2} + 1)} \ln n + \frac{n^{2} {(n^{2} - 1)}^{2}}{{(n^{2} + 1)}^{3}} \ln (\frac{n - 1}{n + 1})

where

n = n_{1} / n_{2}

is the ratio of refractive indices.

Note:

n_{1}

always represents the medium with higher refractive index and

n_{2}

the medium of lower refractive index.

The transmission from medium one to two is given by

τ_{D, 12} = 1 - r_{D, 12}

For the reverse direction the reflectance and transmittance are given by:

τ_{D, 21} = 1 - \frac{1}{n^{2}} (1 - r_{D, 12})

and

τ_{D, 21} = \frac{1}{n^{2}} τ_{D, 12}

respectively.

The incoming radiative intensity to the interface is given by the hemispherically averaged intensity for medium one and two:

Q_{1} = \sum_{j = 1}^{N} w_{j} I_{j} |n \cdot Ω_{j}| (n \cdot Ω_{j}) when n \cdot Ω_{j} > 0

Q_{2} = \sum_{j = 1}^{N} w_{j} I_{j} |n \cdot Ω_{j}| (n \cdot Ω_{j}) when n \cdot Ω_{j} \leq 0

where $n$ is the outward facing normal. From these fluxes, the outgoing radiative intensity at the wall for the current ordinate direction, $Ω$ , is given by

I (Ω) = r_{D, 12} \frac{Q_{1}}{π} + τ_{D, 21} \frac{Q_{2}}{π}

Reflection and Transmission for partially specular and partially diffuse interfaces

For partially specular and partially diffuse interfaces 0.0 < diffused fraction < 1.0.

Interfaces between semi-transparent media are typically not 100 percent diffused or specular and the diffuse fraction lies somewhere between zero and one. In this range the outgoing radiative intensity is treated as a linear combination of the specular and diffuse components, for example:

I (Ω) = (1 - α) I^{S} (Ω) + α I^{D} (Ω)

where $α$ is the diffuse fraction, $I^{S} (Ω)$ is the outgoing specular component of radiative intensity, and $I^{D} (Ω)$ is the outgoing diffuse component of radiative intensity. For example, in medium one in the image above the components would be:

I^{S} (Ω) = I (Ω_{R}^{1}) = r_{12} I (Ω_{I}^{1}) + τ_{21} {(\frac{n_{1}}{n_{2}})}^{2} I (Ω_{I}^{2})

and

I^{D} (Ω) = r_{D, 12} \frac{Q_{1}}{π} + τ_{D, 21} \frac{Q_{2}}{π}

Specular and diffuse interfaces

For the interface equations to be applied weakly, I_w must be applied in both mediums. If the current ordinate direction, $Ω$ , is outgoing from the interface in medium 1 then I_w is equal to the proportion of radiative intensity reflected and transmitted. From the analysis in the previous section, this would be:

I_{w} = r_{12} I (Ω_{I}^{1}) + τ_{21} {(\frac{n_{1}}{n_{2}})}^{2} I (Ω_{I}^{2})

In medium 2 since the current ordinate direction, $Ω$ , is incoming to the surface no boundary flux is added to the equation, that is, $η {(I_{w}, v)}_{Γ} = 0$ .

Reflection and Transmission for diffuse interfaces of Type External

Exchange of radiative intensity occurs for external surfaces when the medium inside the computational domain is semitransparent. That is the medium surrounding, which is not modeled using a computational mesh, participates in radiative transfer. For this case a mathematical model of external radiation is used. The model assumes that the surrounding medium has uniform radiative intensity in all directions, that is, the radiative flux is isotropic. The isotropic radiative intensity is given by the following blackbody source:

I_{ext} = \frac{ε_{ext} σ T_{ext}^{4}}{π}

where $ε_{ext}$ is the external emissivity and is set to one, $σ$ is the Stefan-Boltzmann constant, $T_{ext}$ is the temperature of the surrounding medium. At the external interface $I_{ext}$ is transferred into the medium. This condition can only be applied to boundaries as the interface is only modeled mathematically.