Participating Media Radiation

Simcenter STAR-CCM+ accounts for participating media effects by using the Discrete Ordinate Method (DOM) or Spherical Harmonics.

The DOM is a numerical model for describing radiative heat transfer in participating media. It simulates thermal radiation exchange between diffuse/specular surfaces forming a closed set.

The Spherical Harmonics model simulates thermal radiation propagating anisotropically, the pattern of propagation being expressed in terms of spherical harmonic functions. It does not simulate refraction or specular reflection.

The surface radiative properties are quantified in terms of emissivity, specular and diffuse reflectivity, transmissivity, and radiation temperature. These properties do not depend on direction. For the multiband spectrum model, the surface properties can depend on wavelength.

The medium that fills the space between the surfaces can also absorb, emit, or scatter radiation. Therefore, the amount of radiation that each surface receives and emits depends on this effect, as well as the optical properties of the surface and the thermal boundary conditions that are imposed on it.

As radiation travels through a medium, the intervening material absorbs and increments its radiant intensity $I$ in the $Ω$ direction. The radiative transfer equation (RTE) governs this process. This equation can be written in terms of radiant intensity for a specific wavelength $λ$ :

Figure 1. EQUATION_DISPLAY

\frac{d I_{λ}}{d s} = - β_{λ} I_{λ} + k_{a λ} I_{b λ} + \frac{k_{s λ}}{4 π} \int_{4 π}^{} I_{λ} Ω d Ω + k_{p a λ} I_{p b λ} + \frac{k_{p s λ}}{4 π} \int_{4 π}^{} I_{λ} Ω d Ω

(1721)

where:


$I_{λ}$	radiative intensity at wavelength $λ$ [ $W / (m^{2} s r m^{- 1})$ ]
$I_{b λ}$	black body intensity at wavelength $λ$
$I_{p b λ}$	particle black body intensity at wavelength $λ$ and current particle temperature.
$s$	distance in the $Ω$ direction
$β_{λ}$	extinction coefficient, which is defined as $β_{λ} = k_{a λ} + k_{s λ} + k_{p a λ} + k_{p s λ}$
$k_{a λ}$	absorption coefficient at wavelength $λ$ $(m^{- 1})$
$k_{s λ}$	scattering coefficient at wavelength $λ$ $(m^{- 1})$
$k_{p a λ}$	particle absorption coefficient at wavelength $λ$ $(m^{- 1})$
$k_{p s λ}$	particle scattering coefficient at wavelength $λ$ $(m^{- 1})$
$Ω$	solid angle

The black body intensity is given by:

Figure 2. EQUATION_DISPLAY

I_{b λ} = \frac{2 C_{1}}{λ^{5} (e^{C_{2} / λ T} - 1)}

(1722)

where $C_{1} = 0.595522 \times 10^{- 16}$ $W m^{2} / s$ and $C_{2} = 0.01439$ $m K$ .

The in-scattering component is assumed to be isotropic.

When the medium has several components, the absorption coefficient $k_{a λ}$ is given by:

Figure 3. EQUATION_DISPLAY

k_{a λ} = \sum_{i}^{} k_{a λ i}

(1723)

where $i$ indexes the participating components, such as CO₂, H₂O, or soot.

When the absorption and scattering coefficients of the medium are independent of wavelength, the medium is called gray. In that case, the RTE can be integrated over wavelength (or, equivalently, wavenumber) to produce a wavelength-independent equation.

The particle radiation is valid for both gray and multiband radiation for continuous media. However, for the case of multiband (spectral) radiation, the particle properties are still gray. The surface and continuum properties can be spectral, but the particle properties are gray (the absorption efficiencies are independent of wavelength). Particle scattering is isotropic.

Boundary Conditions: Reflection, Refraction, and Radiosity

The boundary condition that is applied to the RTE for each wavelength $λ$ is:

Figure 4. EQUATION_DISPLAY

I_{λ} (s) = ε_{w λ} I_{b λ} + \frac{E_{w λ}^{d}}{π} + \frac{ρ_{w λ}^{d}}{π_{e f f}} \int_{(n \cdot s^{'} < 0)}^{} I_{λ} (s') | n \cdot s^{'} | d Ω + ρ_{w λ}^{s} I_{λ} (s^{s}) + τ_{w λ} (s^{i}) I_{λ} (s^{i})

(1724)

where:

$ε_{w λ}$ is the diffuse emissivity
$\frac{E_{w λ}^{d}}{π}$ is the diffuse non-Planck emission term
$ρ_{w λ}^{d}$ is the diffuse reflectivity
$π_{e f f} = \sum_{s_{i} \cdot n_{i}} w_{i} (s_{i} \cdot n_{i})$ is the half moment for the ordinate set chosen
$ρ_{w λ}^{s}$ is the specular reflectivity
$s'$ is the unit vector along the distance coordinate leaving the boundary
$s^{s} = s - 2 (s \cdot n) n$ is the incoming ray direction as determined by the laws of reflection
$s^{i}$ is the incoming ray direction as determined by the laws of refraction or Snell's law
$n$ is the surface normal
$τ_{w λ}$ is the transmissivity. Only interfaces can be partially transmissive; all other external boundaries are modeled as opaque

The $w$ subscript indicates the state at the boundary. Surface normal points toward the surface. Radiative properties such as specular reflectivity and transmissivity at the boundary become direction-dependent for refraction.

The boundary radiosity is:

J_{w} = \int_{n \cdot s < 0} ϵ_{w} I_{b w} | n \cdot s | d Ω + \int_{n \cdot s > 0} ρ_{w}^{d} (s^{'}) I (s') | n \cdot s^{'} | d Ω + \int_{n \cdot s^{s} > 0} ρ_{w}^{s} (s^{s}) I (s^{s}) | n \cdot s^{s} | d Ω

(1725)

where:

$I_{b w}$ is the total, spectrally integrated blackbody intensity at the boundary.
$ρ_{w}^{d}$ and $ρ_{w}^{s}$ are the diffuse and specular reflectivities at the boundary.

Boundary radiosity on external boundaries does not include contributions due to transmission from the environment. However, for interface boundaries, radiosity includes contributions due to transmission from adjacent regions. This is done to keep DOM consistent with S2S. For DOM, radiosity is a post-processing parameter.

Radiant Heat Flux

The radiant heat flux in a particular direction $q_{r}$ is calculated using the integration of the radiant intensity over all solid angles and over the wavelength spectrum:

Figure 5. EQUATION_DISPLAY

q_{r} (r) = \int_{0}^{\infty} \int_{4 π} I_{λ} (s) s d Ω d λ

(1726)

The radiation solution is coupled to the fluid dynamic solution through the divergence of the radiative heat flux. This term exchanges energy between the fluid and the radiant energy field. Given the intensity field, the divergence of the heat flux is computed as:

Figure 6. EQUATION_DISPLAY

\nabla \cdot {\dot{q}}_{r} = \int_{0}^{\infty} k_{a λ} (4 π I_{b λ} - \int_{4 π}^{} I_{λ} d Ω) d λ

(1727)

Refractive Index and Blackbody Intensity

Local blackbody intensity $I_{b}$ is augmented by the square of the absolute refractive index $n$ :

I_{b} = \frac{n^{2} σ T^{4}}{π}

(1728)

This equation applies both within the medium and at boundaries. In the medium, blackbody intensity is multiplied by absorption coefficient, while at boundaries, it is multiplied by boundary emissivity.

When radiation is transmitted or refracted from one medium to another, the intensity changes by:

\frac{I_{2}}{n_{2}^{2}} = \frac{I_{1} τ_{1}}{n_{1}^{2}}

where subscript $1$ indicates values on the incoming side of the boundary and subscript $2$ , for the outgoing side. The increase of intensity does not violate the emission maximum given by Planck's Law. Only the local intensity rises, not the global intensity for the whole medium.

Refraction, Reflection, and Critical Angle

Consider an interface boundary between two continua with refractive indices $n_{1}$ and $n_{2}$ , with a light ray passing from the first to the second. There are two possible cases:

Case 1: the ray passes from low to high refractive index, $n_{2} > n_{1}$
Case 2: the ray passes from high to low refractive index, $n_{2} < n_{1}$

$n_{i}$ and $n_{o}$ are the inner and outer unit normal vectors
$s_{i}$ and $s_{t}$ are the unit direction vectors of the incident and transmitted rays
$θ_{i}$ and $θ_{t}$ are the angles of incidence and transmission

The angles are related to the indices of refraction by Snell's law:

n_{1} \sin (θ_{i}) = n_{2} \sin (θ_{t})

In Case 1, the transmitted ray moves closer to the normal on the outer side, so all incident rays are transmitted through the surface.

In Case 2, the transmitted ray moves further from the normal on the outer side. As $θ_{i}$ increases, $θ_{t}$ eventually reaches 90° to the normal and is parallel to the surface. Beyond this point, the ray is reflected back into the first medium. The angle at which this change occurs is called the "critical angle" $θ_{c}$ .

θ_{c} = \sin^{-1} (\frac{n_{2}}{n_{1}})

(1729)

The parameter calculated and stored for critical angle checking during run time is:

\cos^{2} (θ_{c}) = 1 - {(\frac{n_{2}}{n_{1}})}^{2}

For Case 1, because $n_{2} > n_{1}$ , the righthand side of this equation is negative, meaning there is no critical angle and no internal reflection.

For Case 2, because, $n_{2} < n_{1}$ , the righthand side of this equation is between 0 and 1, meaning a critical angle and internal reflection exist.

The vectors of the incident and transmitted rays are related as follows:

Figure 7. EQUATION_DISPLAY

s_{t} = \frac{n_{1}}{n_{2}} (s_{i} - (s_{i} \cdot n_{i}) n_{i}) - \sqrt{1 - \sin^{2} (θ_{t})} n_{o}

(1730)

where:

\sin^{2} (θ_{t}) = {(\frac{n_{1}}{n_{2}})}^{2} (1 - {(s_{i} \cdot n_{i})}^{2})

Correction Factors for Reflection and Refraction

Since the outgoing direction calculated by the laws of reflection and refraction does not always line up with the preset ordinate directions, Simcenter STAR-CCM+ uses a nearest neighbor ordinate direction. This modeling creates slight misalignment in the direction of outgoing energy, which requires a correction factor to conserve total energy over that boundary face integrated over the hemisphere.

The correction factor for specular reflection is:

Figure 8. EQUATION_DISPLAY

f^{r e f l} = \frac{\int_{n \cdot s > 0} ρ_{w}^{s} (s^{i}) I (s^{i}) | n \cdot s | d Ω}{\int_{n \cdot s < 0} ρ_{w}^{s} (s_{q u a d}^{i}) I (s_{q u a d}^{i}) | n \cdot s | d Ω}

(1731)

where:

$ρ_{w}^{s}$ is specular reflectivity. For boundaries without refraction, this is a constant spatially and directionally, and drops out of the integral
$s^{i}$ is the incoming direction
$s_{q u a d}^{i}$ is the incoming direction related to each of the outgoing quadrature directions, found using the laws of reflection and nearest neighbor quadrature

For semi-transparent boundaries with refraction, specular reflectivity is not always constant over the hemisphere. If the critical angle $θ_{c}$ is less than 90°, specular reflectivity outside that cone of incidence is equal to specular reflectivity inside plus transmissivity, $ρ_{w}^{s} + τ_{w}$ .

The correction factor for refraction transmission is:

Figure 9. EQUATION_DISPLAY

f^{refrac} = \frac{\int_{n \cdot s > 0} τ_{w} (s^{r}) I (s^{r}) | n_{1} \cdot s^{r} | d Ω}{\int {}_{n \cdot s < o}{τ_{w} (s_{q u a d}^{r}) I (s_{q u a d}^{r}) | n_{2} \cdot s |} d Ω}

(1732)

where:

$s^{r}$ is the incoming direction of transmission
$s_{q u a d}^{r}$ is the incoming direction of transmission related to each of the outgoing quadrature directions, found using the laws of reflection and nearest neighbor quadrature

To find $s_{q u a d}^{r}$ , first calculate a theoretical incident ray $i$ using:

i = \frac{n_{2}}{n_{1}} (s_{t} - (s_{t} \cdot n_{o}) n_{o}) - \sqrt{1 - \sin^{2} (θ_{t})} n_{i}

with

\sin^{2} (θ_{i}) = {(\frac{n_{2}}{n_{1}})}^{2} (1 - {(s_{t} \cdot n_{o})}^{2})

This theoretical incident ray $i$ is then used to find the nearest quadrature direction within the critical angle $θ_{c}$ .

For the transmissivity $τ_{w}$ , if the critical angle $θ_{c}$ is less than 90°, transmissivity is zero outside the cone of incidence.

The outgoing intensity is then given by Eqn. (1724) with the DOM direction correction included:

Figure 10. EQUATION_DISPLAY

I_{λ} (s) = ε_{w λ} I_{b λ} + \frac{ρ_{w λ}^{d}}{π_{e f f}} \int_{n \cdot s' < 0} I_{λ} (s') | n \cdot s^{'} | d Ω + f^{refl} ρ_{w λ}^{s} (s_{q u a d}^{i}) I_{λ} (s_{q u a d}^{i}) + f^{refrac} τ_{w λ} (s_{q u a d}^{i}) I_{λ} (s_{q u a d}^{r})

(1733)

Full Spectrum k-Distribution Formulation

The k-Distribution Thermal Radiation model is used to model spectral variations of H₂O and CO₂ within participating media. The radiative transfer equations (RTE) are solved on a spectrally reordered basis when the model is selected. See Modest [411]. The Correlated k-Distribution Method describes how the absorption coefficient values are calculated and what they represent. The full-spectrum formulation then applies these coefficients at each spectral location to give the spectrally integrated total radiative heat flux:

Figure 11. EQUATION_DISPLAY

\frac{d I_{k}}{d s} = k f (T, k) I_{b} + k_{p a} I_{b} + (k + k_{p s} + k_{s} + k_{p a}) I_{k} + \frac{k_{s}}{4 π} \int_{4 π}^{} I_{k} (Ω) d Ω + \frac{k_{p s}}{4 π} \int_{4 π}^{} I_{k} (Ω) d Ω

(1734)

where:

$k_{p a}$ is the particle absorption coefficient.
$k_{p s}$ is the particle scattering coefficient.
$k_{s}$ is the scattering coefficient.
$k_{p a}$ , $k_{p s}$ , and $k_{s}$ are not functions of wavelength.
$Ω$ is the direction vector.

with the boundary condition

Figure 12. EQUATION_DISPLAY

I_{k} = ε_{w} {f (T_{w}, k) I}_{b w} + \frac{E_{w}^{d}}{π} + \frac{ρ_{w}^{d}}{π_{e f f}} \int_{(\hat{n} \cdot \hat{s})} I_{k} (s') | n \cdot s' | d Ω + ρ_{w}^{s} I_{k} (s^{i})

(1735)

where:

Figure 13. EQUATION_DISPLAY

I_{k} = \int_{0}^{\infty} I_{λ} δ (k - k_{a}) d λ

(1736)

where the total spectrally integrated intensity is given by:

I = \int_{0}^{\infty} I_{b λ} ⅆ λ = \int_{0}^{\infty} I_{k} d k

To facilitate spectral integration, Eqn. (1734) and Eqn. (1735) are divided by Eqn. (1736) evaluated at a given reference state. Assuming isotropic scattering, this results in:

Figure 14. EQUATION_DISPLAY

\frac{d I_{g}}{d s} = (k^{*} (T_{0}, ϕ, g) + {\bar{κ}}_{p}) [a (T, T_{0}, g) I_{b} (T) - I_{g}] - k_{s} (I_{g} - \frac{1}{4 π} \int_{V̵ Ω} I_{g} d Ω')

(1737)

with the boundary condition

Figure 15. EQUATION_DISPLAY

I_{g} = ε a_{w} I_{b λ w} + \frac{ρ_{w}^{d}}{π_{e f f}} \int_{n \cdot s < 0} I_{g} (s') | n \cdot s | d Ω + ρ_{w}^{s} I_{g} (s^{s}) + τ_{w} (s^{i}) I_{g} (s^{i})

(1738)

where:

I_{g} = I_{k} / f (T_{0}, ϕ_{0}, k)

g (T, ϕ_{0}, k) = \int_{0}^{k} f (T, ϕ_{0}, k) d k = g (T, ϕ_{0}, k^{*})

a (T, T_{0}, g) = \frac{f (T, ϕ_{0}, k)}{f (T_{0}, ϕ_{0}, k)}

a_{w} = \frac{f (T_{w}, ϕ_{0}, k)}{f (T_{0}, ϕ_{0}, k)}

Gaussian Quadrature Scheme

The k-distribution method differs from the WSGG method in that it is formulated assuming a continuous distribution rather than finite steps, allowing for higher order integration techniques. Wang and Modest have proposed the use of Gauss-Chebyshev integration methods for the selection of spectral points (in reordered g-space), and their subsequent integration. The abscissas of a Gaussian quadrature between the limits of -1 and 1 are the zeroes of a particular orthogonal polynomial. For the nth rank Chebyshev polynomial of the second kind, the zeroes are:

Figure 16. EQUATION_DISPLAY

x_{n}^{k} = \cos (θ_{k}), k =, 1, 2, ..., n

with the corresponding Gaussian weight functions as:

w_{n}^{k} = \frac{4 \sin (θ_{k})}{n + 1} \sum_{l = 1}^{(n + 1) / 2} \frac{\sin (2 l - 1) θ_{k}}{2 l - 1}

(1739)

where:

θ_{k} = \frac{k π}{n + 1}

For the k-distribution application, g-space varies from 0 to 1, so the quadrature points are picked from the upper half, and weights are scaled. Depending on the problem it may be advantageous to bias this distribution so that more points are in the spectral range where most heat transfer occurs. For radiative transport problems this generally occurs in the optically intermediate zone (0.5< kL <100), which for k-distribution problem often happens in g -space > 0.9. By introducing a biasing factor, the user will be able to select the quadrature to obtain higher accuracy per number of bands selected. Wang and Modest introduce the transformation:

Figure 17. EQUATION_DISPLAY

{x'}_{n} = 1 - {(1 - x_{n})}^{α}

(1740)

{w'}_{n} = α w_{n} {(1 - x_{n})}^{α - 1}

Wang and Modest recommend a value for the shape factor $α$ in the range $1.5 < α < 2.5$ . The default value is 1.5.

The number of quadrature points and the shape factor are chosen to accurately resolve the variation in the absorption coefficient encountered in a particular problem. In most cases, the default values of eight points and a shape factor of 1.5 is best. In some rare cases, if the absorption coefficient is small for all bands, a lower number of points with a higher shape factor is sufficient to resolve the radiative transfer.