3.2.4.17. Simplified Spherical Harmonics Equations

Radiative transport, which is another mode of heat transfer, describes the spatial variation of radiative intensity corresponding to a given direction and at a given wavelength within a radiatively participating medium. This is governed by the Boltzmann Equation or Radiative Transport Equation (RTE) which represents a balance between absorption, emission and scattering. For the purpose of thermal analysis, we make certain assumptions that simplify the RTE. These are

Gray Media

The optical properties of the material medium are not a function of wavelength i.e., no spectral banding.

Steady State

Energy transport is propagated by photons which travel at the speed of light. This means that there are no transient effects relative to the time scale of other physics.

Isotropic Scattering

Since the scattering phase function is rarely known, we simply assume isotropic scattering where the probability of intensity being propagated in a given direction has equal probability of being scattered in another direction.

Given the above assumptions, we may write the RTE as

(3.76) $\Omega \cdot \grad I\left(\Omega\right) + \left(\sigma_A + \sigma_S \right) I\left(\Omega\right) = \sigma_A I_b + \frac{\sigma_S}{4\pi}G,$

where $\sigma_A$ is the absorption coefficient, $\sigma_S$ is the scattering coefficient, $I(\Omega)$ is the intensity along the direction $\Omega$ , $T$ is the temperature and the angle-integrated intensity is $G=\int_{4\pi}{I\left(\Omega\right)d\Omega}$ . The black body radiation, $I_b$ is defined by $\frac{\sigma T^4}{\pi}$ . Note that for situations in which the scattering coefficient is zero, the RTE reduces to a set of linear, decoupled equations for each intensity to be solved.

Discrete Ordinate (DO) methods such as $S_N$ have been developed to help solve this set of coupled equations by selecting a set of prescribed directions $\Omega_k$ , solving for the discrete intensities $I_k$ and combining these intensities to represent the radiative field. These methods suffer from high computational cost due to the fact that a large number of ordinate directions is required to adequately resolve the intensity field and avoid ray effects. For a given quadrature order N, the number of ordinate directions scales as $N^2$ .

Beta Capability

$SP_N$ is not well tested and should not be considered a production capability.

The simplified spherical harmonics ( $SP_N$ ) approximation to the radiative transport equation (RTE) was first proposed by Gelbard [14] for reactor analysis in the early 1960s. The initial derivation involved replacing the spatial derivatives in the 1D spherical harmonics ( $P_N$ ) approximation with their 3D analogs. It was later shown that a more rigorous theoretical basis was possible and that the $SP_N$ equations may be derived as either an asymptotic correction to the diffusion limit or from the use of certain trial functions in the self-adjoint variational characterization of the even-parity form of the RTE. There are a number of equivalent forms for the $SP_N$ equations in the literature. We choose to use the canonical form which is derived from the 1D even-parity discrete ordinates equations.

(3.77) $-\grad \bcdot \left( \frac{\mu_n^2}{\sigma_T} \grad I_n \right) + \sigma_T I_n = 4\pi\sigma_S \sum_{m=1}^{\frac{N+1}{2}} w_m I_m + \sigma_A I_b$

where $\sigma_T$ is the extinction coefficient, $\sigma_T=\sigma_A+\sigma_S$ , $\mu_n$ is the nth quadrature point in a (N+1)-point Gauss set on [-1, 1], $I_n$ is the angular intensity at quadrature point $n$ , and $w_n$ is the $n^{th}$ quadrature weight.

The appropriate boundary conditions for the canonical SPn equations are derived from the boundary conditions for the 1D even-parity discrete ordinates equations. We consider a Mark boundary condition of the 1st order intensity BC

$I_n = \epsilon I_b + \left(1-\epsilon\right)\sum_{\vec{n} \bcdot \vec{\Omega}_j < 0} I_j w_j \left| \vec{n} \bcdot \vec{\Omega}_j \right|$

such that in canonical form, it yields

(3.78) $-\frac{\mu_n}{\sigma_T} \grad I_n \bcdot \vec{n} = \frac{\epsilon}{2-\epsilon}\left( I_n - I_b \right) + \frac{1-\epsilon}{2-\epsilon} \left[ \frac{ \sum_k \left( I_k - \frac{\mu_k}{\sigma_T} \grad I_k \bcdot \vec{n} \right) \mu_k w_k }{\sum \mu_k w_k} - I_n + \frac{\mu_n}{\sigma_T} \grad I_n \bcdot \vec{n} \right]$

where $\epsilon$ is the emissivity and $\vec{n}$ is the surface normal.

The $SP_N$ equations are solved for the unknown intensities to provide an approximate solution to the RTE. The angle-integrated intensity is then found and used in the source term for the material energy equation. It is to be noted from the form of the $SP_N$ equations that this cannot be used in vacuum or near-vacuum situations where $\sigma_T ~ 0$ . Also analysis by Zheng et al. shows that the existence of a unique solution may be proved when absorption effects are non-negligible and the geometry is small.