2. Theory

PMR solves the radiative transport equation (RTE) using the method of discrete ordinates with user-selectable quadrature rules. The quadrature rule and order defines a set of discrete directions in which intensity is solved. Scattering is not currently supported, so this results in a decoupled set of linear PDEs in each ordinate direction.

For each ordinate direction ( $\symSK$ ), a weight is provided ( $\symWK$ ) by the quadrature rule selected. For each ordinate direction, the ordinate intensity $\symIK$ is solved for and overall scalar flux ( $\symScalarFlux$ ) and radiative heat flux ( $\vec{q}$ ) can then be defined by summing over ordinate directions as

$\symScalarFlux = \sum_k \symIK \symWK$

and

$\vec{q} = \sum_k \symIK \symSK \symWK$

The continuous form of the PDE evaluated for intensity (excluding scattering terms) in a given ordinate direction is

$\symSK \cdot \nabla \symIK = \symRadSrc - \symAbs \symIK$

where the radiative source term ( $\symRadSrc$ ) and absorption coefficient ( $\symAbs$ ) are calculated by the fluid code and transferred to PMR. The radiative source term is usually similar to the form below, but may include subgrid effects where relevant:

$\symRadSrc \simeq \frac{\symAbs \symSB T^4}{\pi}$

On the domain boundaries, the following condition is enforced for intensity

$\symIB = \frac{\symBoundarySrc + \left(1 - \symEmis - \symTrans\right)\symIrrad}{\pi}$

where $\symIrrad$ is the irradiation (sum of incoming intensities on the boundary) and the boundary source ( $\symBoundarySrc$ ) is sent from the fluids code and usually has a form similar to the following equation:

$\symBoundarySrc = \symSB \left( \symEmis T_{bc}^4 + \symTrans T_{env}^4 \right)$

In the case where $\symEmis + \symTrans < 1$ the boundary is partially reflective and the problem becomes nonlinear and coupled across ordinate directions through $H$ . When this is the case, multiple nonlinear iterations should be used to achieve a converged solution. The irradiation is calculated on all boundaries as

$\symIrrad = \sum_k \begin{array}{ll} \symWK \symIK \vec{n} \cdot \symSK & \mbox{if $\vec{n} \cdot \symSK > 0$} \\ 0 & \mbox{if $\vec{n} \cdot \symSK \le 0$} \end{array}$

where $\vec{n}$ is the outward facing normal vector.

The primary numerical choices when solving this equation are what type of stabilization to use for the advection-like term (upwind or SUCV) and whether to use an edge-based or element-based stencil.

2.1. Quadrature Rules

2.1.1. Thurgood

The Thurgood [1] quadrature rule partitions space into $8 N^2$ ordinate directions, where $N$ is the specified quadrature order. The order must be at least 2 for this quadrature rule.

2.1.2. Level-Symmetric

The Level-Symmetric [2] quadrature rule partitions space into $N(N+2)$ ordinate directions, where $N$ is the specified quadrature order. The order must be one of: 2, 4, 6, 8, 12, or 16. No other orders are valid for this quadrature rule.

2.1.3. PNTN

The Legendre-Chebyshev (PN-TN) [3] quadrature rule partitions space into $N(N+2)$ ordinate directions, where $N$ is the specified quadrature order. The order must be even for this quadrature rule.