Discrete Ordinate Method Numerical Solution

The discrete ordinates method solves field equations for radiation intensity that is associated with a fixed direction $s$ , represented by a discrete solid angle.

The following diagram illustrates a discrete solid angle as part of a hemisphere:

A solid angle is the three-dimensional analog of an ordinary angle and is measured in steradians [ $s r$ ]. One steradian is the solid angle that is subtended at the center of a sphere by an area on the surface of the sphere equal to the radius squared. As such, the solid angle of a complete sphere is $4 π$ steradians.

The discrete ordinates calculation requires specification of the number of solid angles or “ordinates” into which the sphere surrounding any particular point is discretized. Increasing the number of ordinates leads to greater accuracy. A detailed description of the method can be found in Modest [398] and Siegel and Howell [419]. As a result of Eqn. (1721), the form of these ordinate equations (for each wavelength band) is:

Figure 1. EQUATION_DISPLAY

s_{i} \cdot \nabla I_{i Δ λ} = - β_{Δ λ} I_{i Δ λ} + k_{a Δ λ} I_{b Δ λ} + \frac{k_{s Δ λ}}{4 π} \sum_{j = 1}^{n} w_{j} I_{j Δ λ} + {\bar{k}}_{p α Δ λ} I_{p b Δ λ} + \frac{{\bar{k}}_{p s Δ λ}}{4 π} \sum_{j = 1}^{n} w_{j} I_{j Δ λ}

(1741)

In the above equation, $Δ λ$ represents a wavelength band from $λ_{m}$ to $λ_{n}$ .

The black body emission in this band is:

Figure 2. EQUATION_DISPLAY

I_{b Δ λ} = \int_{0}^{λ_{n}} I_{b λ} d λ - \int_{0}^{λ_{m}} I_{b λ} d λ

(1742)

Boundary Conditions

The transport equation for each ordinate direction is discretized and solved independently, using the same discretization and iterative solution practices as for the other transport equations. For this reason, global iterations are necessary to include the isotropic in-scattering terms in the RTE and to compute wall boundary conditions. The angularly discretized boundary conditions take the form:

Figure 3. EQUATION_DISPLAY

I_{i Δ λ, w} = ε_{w Δ λ} I_{b Δ λ} + \frac{ρ_{w Δ λ}^{d}}{π_{eff}} \sum_{(n \cdot s^{'} < 0)}^{} I_{j Δ λ} {| n \cdot {s^{'}}_{j} |}_{w} + f ρ_{w Δ λ}^{s} I_{Δ λ} (s^{s}) + τ_{w Δ λ} I_{i Δ λ}

(1743)

where:

$π_{e f f}$ is the half moment for the chosen ordinate set, which is defined as:
Figure 4. EQUATION_DISPLAY

$π_{eff} = \sum_{n \cdot s^{'} > 0}^{} w_{j} {(n \cdot {s^{'}}_{j})}_{w}$

(1744)
$f$ is the correction factor that is used when the direction of specular reflection is not aligned with any discrete ordinate direction.

The diffuse reflection term represents summation over incoming (incident) ordinate directions.

Convergence is measured using the infinity norm of the normalized difference in the incident energy from successive space-angle sweeps.

Radiation Source Term

The radiation solution provides a source term to the fluid dynamic energy equation as given by Eqn. (1727). The discretized form of this source term at any cell is:

Figure 5. EQUATION_DISPLAY

- \nabla \cdot q_{r} = \sum_{λ}^{} k_{a Δ λ} (\sum_{j = 1}^{n} w_{j} I_{j Δ λ} - 4 π I_{b Δ λ})

(1745)

Absorption by Particles

The net absorption rate of radiative energy by all particles (Largrangian parcels) within a given cell is:

Figure 6. EQUATION_DISPLAY

- \nabla \cdot q_{p r} = \sum_{λ}^{} {\bar{k}}_{p a Δ λ} (\sum_{j = 1}^{n} w_{j} I_{j Δ λ} - 4 π I_{p b Δ λ})

(1746)

For Lagrangian particle transport, at convergence this quantity exactly balances the net radiative energy that all particles in the cell absorb.

The equivalent particle absorption and scattering coefficients are given in [398] as:

Figure 7. EQUATION_DISPLAY

{\bar{k}}_{p a Δ λ} = \sum_{i}^{} Q_{a, i} N_{i} \frac{π d_{i}^{2}}{4}

(1747)

Figure 8. EQUATION_DISPLAY

{\bar{k}}_{p s Δ λ} = \sum_{i}^{} Q_{s, i} N_{i} \frac{π d_{i}^{2}}{4}

(1748)

where:

$i$ is the parcel index number.
$Q_{a, i}$ is the absorption coefficient for the $i$ th parcel.
$Q_{s, i}$ is the scattering coefficient for the $i$ th parcel.
$N_{i}$ is the number of particles per unit volume for the $i$ th parcel.
$d_{i}$ is the particle diameter for the $i$ th parcel.

Discrete Ordinate Quadrature

The ordinate set that is required to solve Eqn. (1741) with boundary conditions (Eqn. (1743)) is obtained using the discrete ordinate quadrature scheme. General principles for DOM quadrature schemes [409] are:

All discrete ordinates, $Ω_{m} = (ζ_{m}, η_{m}, μ_{m})$ , where $(ζ_{m}, η_{m}, μ_{m})$ are the direction cosines, must be placed on the unit sphere, that is:
Figure 9. EQUATION_DISPLAY

$(ζ_{m}^{2} + η_{m}^{2} + μ_{m}^{2}) = 1$
(1749)
The weights that are associated with each ordinate direction must be all positive. This requirement ensures that the error of the quadrature scheme is minimized.
The number of photons must be preserved which is equivalent to satisfy:
Figure 10. EQUATION_DISPLAY

$\sum_{m = 1}^{M} w_{m} (n \overline{Ω_{m}}) = 0$
(1750)

where $w_{m}$ is the weight that is associated with the ordinates and $n$ is an arbitrary unit vector.
The quadrature must be invariant to any rotation of the arrangement of the discrete directions around the center of the unit sphere.

Although the first three criteria can be easily fulfilled, the fourth is impossible to be satisfied for any arbitrary rotation by a finite number of ordinates. Hence all ordinate sets must satisfy at least zeroth, first and second moments. For accuracy in radiation calculations, the ordinate set with weights must also satisfy the first moment over a half range along the three principal directions. This expression is:

Figure 11. EQUATION_DISPLAY

\sum_{n \overline{Ω_{m}} > 0}^{} w_{m} (n \overline{Ω_{m}}) = π

(1751)

This additional requirement comes from the following:

Radiative intensity can have directional discontinuity at a wall.
Important radiative heat fluxes at the walls are evaluated through a first moment of intensity over a half range of $2 π$ .

Level Symmetric (Sn) Quadrature

Level symmetric, also known as Sn, quadrature is a quadrature scheme where the discrete ordinates are arranged on the latitudes on the unit sphere. The corresponding latitudes have the same distance in the x-direction, y-direction, and z-direction and therefore this quadrature scheme is called Level Symmetric quadrature. A typical arrangement is illustrated in the figure (symmetric point arrangement for n = 6) [410].

The rotational invariance to the principal octahedron is satisfied implicitly because of the arrangement on latitudes. Due to the latitude-based arrangement, Sn quadrature can be used in solution of discrete ordinate equations both in three-dimensional and axisymmetric two-dimensional problems. Currently only Sn quadratures are available in Simcenter STAR-CCM+. Sn quadratures can be formulated to satisfy either odd or even moments. To find the Sn formulation that satisfies a particular moment, see the following table. Sn(O) is the odd formulation and Sn(E) is the even formulation.


	Type
Order	Sn(O)	Sn(E)
2	0, 1	0, 2
4	0, 1, 2	0, 2, 4
6	0, 1, 2, 3	0, 2, 4, 6
8	0, 1, 2, 3, 4	0, 2, 4, 6, 8
12	0, 1, 2, 3, 4, 5, 6	0, 2, 4, 6, 8, 10, 12
16	0, 1, 2, 3, 4, 5, 6, 7, 8	0, 2, 4, 6, 8, 10, 12, 14, 16

All Sn quadratures ( $N > 2$ ) with odd formulation satisfy the primary requirements of satisfying zeroth, first and second moments whereas the even formulation does not. However odd formulation of Sn incurs negative weights for ( $N \geq 10$ ) which are, in general, not suitable due to inaccuracy. The even formulation of Sn does not have that problem except for N = 10. Hence, in Simcenter STAR-CCM+ odd formulation is available for N = 2, 4, 6, and 8, and even formulation is available for N = 12 and 16. All of the above available Sn quadratures in Simcenter STAR-CCM+ satisfy the half moment criteria along principal directions.