[Fluent Inc. Logo] return to home search
next up previous contents index

11.3.6 The Discrete Ordinates (DO) Radiation Model

The discrete ordinates (DO) radiation model solves the radiative transfer equation (RTE) for a finite number of discrete solid angles, each associated with a vector direction ${\vec s}$ fixed in the global Cartesian system ( $x,y,z$). The fineness of the angular discretization is controlled by you, analogous to choosing the number of rays for the DTRM. Unlike the DTRM, however, the DO model does not perform ray tracing. Instead, the DO model transforms Equation  11.3-1 into a transport equation for radiation intensity in the spatial coordinates ( $x,y,z$). The DO model solves for as many transport equations as there are directions ${\vec s}$. The solution method is identical to that used for the fluid flow and energy equations.

The implementation in FLUENT uses a conservative variant of the discrete ordinates model called the finite-volume scheme [ 41, 203], and its extension to unstructured meshes [ 179].

The DO Model Equations

The DO model considers the radiative transfer equation (RTE) in the direction ${\vec s}$ as a field equation. Thus, Equation  11.3-1 is written as

\begin{displaymath} \nabla \cdot (I({\vec r},{\vec s}){\vec s}) + (a + \sigma_s)... ...ec s} \; ') \; \Phi({\vec s}\cdot {\vec s} \; ') \; d{\Omega}' \end{displaymath} (11.3-37)

FLUENT also allows the modeling of non-gray radiation using a gray-band model. The RTE for the spectral intensity $I_{\lambda} ({\vec r},{\vec s})$ can be written as

\begin{displaymath} \nabla \cdot (I_{\lambda}({\vec r},{\vec s}){\vec s}) + (a_{... ...ec s} \; ') \; \Phi({\vec s}\cdot {\vec s} \; ') \; d{\Omega}' \end{displaymath} (11.3-38)

Here $\lambda$ is the wavelength, $a_{\lambda}$ is the spectral absorption coefficient, and $I_{b \lambda}$ is the black body intensity given by the Planck function. The scattering coefficient, the scattering phase function, and the refractive index $n$ are assumed independent of wavelength.

The non-gray DO implementation divides the radiation spectrum into $N$ wavelength bands, which need not be contiguous or equal in extent. The wavelength intervals are supplied by you, and correspond to values in vacuum ( $n=1$). The RTE is integrated over each wavelength interval, resulting in transport equations for the quantity $I_{\lambda}\Delta \lambda$, the radiant energy contained in the wavelength band $\Delta \lambda$. The behavior in each band is assumed gray. The black body emission in the wavelength band per unit solid angle is written as

\begin{displaymath}[F(0 \rightarrow n\lambda_2 T) - F(0 \rightarrow n \lambda_1 T)]n^2 \frac{\sigma T^4}{\pi} \end{displaymath} (11.3-39)

where $F(0 \rightarrow n \lambda T)$ is the fraction of radiant energy emitted by a black body [ 174] in the wavelength interval from 0 to $\lambda$ at temperature $T$ in a medium of refractive index $n$. $\lambda_2$ and $\lambda_1$ are the wavelength boundaries of the band.

The total intensity $I({\vec r},{\vec s})$ in each direction ${\vec s}$ at position ${\vec r}$ is computed using

\begin{displaymath} I({\vec r},{\vec s}) = \sum_k I_{\lambda_k}({\vec r}, {\vec s}) \Delta \lambda_k \end{displaymath} (11.3-40)

where the summation is over the wavelength bands.

Boundary conditions for the non-gray DO model are applied on a band basis. The treatment within a band is the same as that for the gray DO model.

Angular Discretization and Pixelation

Each octant of the angular space $4\pi$ at any spatial location is discretized into $N_\theta \times N_{\phi}$ solid angles of extent $\omega_i$, called control angles. The angles $\theta$ and $\phi$ are the polar and azimuthal angles respectively, and are measured with respect to the global Cartesian system $(x,y,z)$ as shown in Figure  11.3.3. The $\theta$ and $\phi$ extents of the control angle, $\Delta \theta$ and $\Delta \phi$, are constant. In two-dimensional calculations, only four octants are solved due to symmetry, making a total of $4N_\theta N_{\phi}$ directions in all. In three-dimensional calculations, a total of $8 N_\theta N_{\phi}$ directions are solved. In the case of the non-gray model, $4N_\theta N_{\phi}$ or $8 N_\theta N_{\phi}$ equations are solved for each band.

Figure 11.3.3: Angular Coordinate System
\begin{figure} \psfig{file=figures/coord.ps,width=3in} \end{figure}

When Cartesian meshes are used, it is possible to align the global angular discretization with the control volume face, as shown in Figure  11.3.4. For generalized unstructured meshes, however, control volume faces do not in general align with the global angular discretization, as shown in Figure  11.3.5, leading to the problem of control angle overhang [ 179].

Figure 11.3.4: Face with No Control Angle Overhang
\begin{figure} \psfig{file=figures/no-overhang.ps,width=4in} \end{figure}

Figure 11.3.5: Face with Control Angle Overhang
\begin{figure} \psfig{file=figures/overhang.ps,width=4in} \end{figure}

Essentially, control angles can straddle the control volume faces, so that they are partially incoming and partially outgoing to the face. Figure  11.3.6 shows a 3D example of a face with control angle overhang.

Figure 11.3.6: Face with Control Angle Overhang (3D)
\begin{figure} \psfig{file=figures/overhang3d.ps,width=4in} \end{figure}

The control volume face cuts the sphere representing the angular space at an arbitrary angle. The line of intersection is a great circle. Control angle overhang may also occur as a result of reflection and refraction. It is important in these cases to correctly account for the overhanging fraction. This is done through the use of pixelation [ 179].

Each overhanging control angle is divided into $N_{\theta_p} \times N_{\phi_p}$ pixels, as shown in Figure  11.3.7.

Figure 11.3.7: Pixelation of Control Angle
\begin{figure} \psfig{file=figures/pixels.ps,height=2.75in} \end{figure}

The energy contained in each pixel is then treated as incoming or outgoing to the face. The influence of overhang can thus be accounted for within the pixel resolution. FLUENT allows you to choose the pixel resolution. For problems involving gray-diffuse radiation, the default pixelation of $1\times1$ is usually sufficient. For problems involving symmetry, periodic, specular, or semi-transparent boundaries, a pixelation of $3\times3$ is recommended. You should be aware, however, that increasing the pixelation adds to the cost of computation.

Anisotropic Scattering

The DO implementation in FLUENT admits a variety of scattering phase functions. You can choose an isotropic phase function, a linear anisotropic phase function, a Delta-Eddington phase function, or a user-defined phase function. The linear anisotropic phase function is described in Equation  11.3-14. The Delta-Eddington function takes the following form:

\begin{displaymath} \Phi({\vec s} \cdot {\vec s} \; ') = 2f\delta({\vec s} \cdot {\vec s} \; ') + (1-f)(1 + C {\vec s} \cdot {\vec s} \; ') \end{displaymath} (11.3-41)

Here, $f$ is the forward-scattering factor and $\delta({\vec s} \cdot {\vec s} \; ')$ is the Dirac delta function. The $f$ term essentially cancels a fraction $f$ of the out-scattering; thus, for $f=1$, the Delta-Eddington phase function will cause the intensity to behave as if there is no scattering at all. $C$ is the asymmetry factor. When the Delta-Eddington phase function is used, you will specify values for $f$ and $C$.

When a user-defined function is used to specify the scattering phase function, FLUENT assumes the phase function to be of the form

\begin{displaymath} \Phi({\vec s} \cdot {\vec s} \; ') = 2f\delta({\vec s} \cdot {\vec s} \; ') + (1-f)\Phi^*({\vec s} \cdot {\vec s} \; ') \end{displaymath} (11.3-42)

The user-defined function will specify $\Phi^*$ and the forward-scattering factor $f$.

The scattering phase functions available for gray radiation can also be used for non-gray radiation. However, the scattered energy is restricted to stay within the band.

Particulate Effects in the DO Model

The DO model allows you to include the effect of a discrete second phase of particulates on radiation. In this case, FLUENT will neglect all other sources of scattering in the gas phase.

The contribution of the particulate phase appears in the RTE as:

\begin{displaymath} \nabla \cdot (I{\vec s}) + (a + a_{p} + \sigma_p) I({\vec r}... ...c s} \; ') \; \Phi({\vec s} \cdot {\vec s} \; ') \; d{\Omega}' \end{displaymath} (11.3-43)

where $a_p$ is the equivalent absorption coefficient due to the presence of particulates, and is given by Equation  11.3-17. The equivalent emission $E_p$ is given by Equation  11.3-16. The equivalent particle scattering factor $\sigma_p$, defined in Equation  11.3-20, is used in the scattering terms.

For non-gray radiation, absorption, emission, and scattering due to the particulate phase are included in each wavelength band for the radiation calculation. Particulate emission and absorption terms are also included in the energy equation.

Boundary Condition Treatment at Gray-Diffuse Walls

For gray radiation, the incident radiative heat flux, $q_{\rm in}$, at the wall is

\begin{displaymath} q_{\rm in} = \int_{{\vec s} \cdot {\vec n} > 0} I_{\rm in} {\vec s} \cdot {\vec n} d \Omega \end{displaymath} (11.3-44)

The net radiative flux leaving the surface is given by

\begin{displaymath} q_{\rm out} = (1 - \epsilon_w)q_{\rm in} + n^2 \epsilon_w \sigma T^{4}_{w} \end{displaymath} (11.3-45)

where $n$ is the refractive index of the medium next to the wall. The boundary intensity for all outgoing directions ${\vec s}$ at the wall is given by

\begin{displaymath} I_0 = \frac{q_{\rm out}}{\pi} \end{displaymath} (11.3-46)

For non-gray radiation, the incident radiative heat flux $q_{{\rm in},\lambda}$ in the band $\Delta \lambda$ at the wall is

\begin{displaymath} q_{{\rm in}, \lambda} = \Delta \lambda \int_{{\vec s} \cdot ... ... n} > 0} I_{{\rm in},\lambda} {\vec s} \cdot {\vec n} d \Omega \end{displaymath} (11.3-47)

The net radiative flux leaving the surface in the band $\Delta \lambda$ is given by

\begin{displaymath} q_{{\rm out}, \lambda} = (1 - \epsilon_{w \lambda})q_{{\rm i... ... T_w) - F(0 \rightarrow n \lambda_1 T_w)] n^2 \sigma T^{4}_{w} \end{displaymath} (11.3-48)

where $\epsilon_{w \lambda}$ is the wall emissivity in the band. The boundary intensity for all outgoing directions ${\vec s}$ in the band $\Delta \lambda$ at the wall is given by

\begin{displaymath} I_{0 \lambda} = \frac{q_{{\rm out},\lambda}}{\pi \Delta \lambda} \end{displaymath} (11.3-49)

The Diffuse Fraction

FLUENT allows you to specify the fraction of incoming radiation that is treated as diffuse at diffuse boundaries. If $q_{\rm in}$ is the amount of radiative energy incident on the wall, then

where $f_d$ is the diffuse fraction and $\epsilon_w$ is the wall emissivity. See below for more information about specular walls.

For non-gray radiation, you can specify the diffuse fraction separately for each band.

Boundary Condition Treatment at Semi-Transparent Walls

FLUENT allows the specification of both diffusely and specularly reflecting semi-transparent walls. You can prescribe the fraction of the incoming radiation at the semi-transparent wall which is to be reflected and transmitted diffusely; the rest is treated specularly.

For non-gray radiation, this treatment is applied on a band basis. The radiant energy within a band $\Delta \lambda$ is transmitted, reflected, and refracted as in the gray case; there is no transmission, reflection, or refraction of radiant energy from one band to another.

Specular Semi-Transparent Walls

Consider a ray traveling from a semi-transparent medium $a$ with refractive index $n_a$ to a semi-transparent medium $b$ with a refractive index $n_b$ in the direction ${\vec s}$, as shown in Figure  11.3.8. Side $a$ of the interface is the side that faces medium $a$; similarly, side $b$ faces medium $b$. The interface normal ${\vec n}$ is assumed to point into side $a$. We distinguish between the intensity $I_{a}({\vec s})$, the intensity in the direction ${\vec s}$ on side $a$ of the interface, and the corresponding quantity on the side $b$, $I_{b}({\vec s})$.

Figure 11.3.8: Reflection and Refraction of Radiation at the Interface Between Two Semi-Transparent Media
\begin{figure} \psfig{file=figures/semi-transparent-interface.ps,width=3in} \end{figure}

A part of the energy incident on the interface is reflected, and the rest is transmitted. The reflection is specular, so that the direction of reflected radiation is given by

\begin{displaymath} {\vec s}_r = {\vec s} - 2\left({\vec s} \cdot {\vec n}\right){\vec n} \end{displaymath} (11.3-50)

The radiation transmitted from medium $a$ to medium $b$ undergoes refraction. The direction of the transmitted energy, ${\vec s}_t$, is given by Snell's law:

\begin{displaymath} \sin \theta_{b} = \frac{n_a}{n_b} \sin \theta_{a} \end{displaymath} (11.3-51)

where $\theta_{a}$ is the angle of incidence and $\theta_{b}$ is the angle of transmission, as shown in Figure  11.3.8. We also define the direction

\begin{displaymath} {\vec s}\; ' = {\vec s}_t - 2\left( {\vec s}_t \cdot {\vec n}\right){\vec n} \end{displaymath} (11.3-52)

shown in Figure  11.3.8.

The interface reflectivity on side $a$ [ 174]

\begin{displaymath} r_a({\vec s}) = \frac{1}{2} \left(\frac{n_a \cos \theta_b - ... ...\cos \theta_b}{n_a \cos \theta_a + n_b \cos \theta_b}\right)^2 \end{displaymath} (11.3-53)

represents the fraction of incident energy transferred from ${\vec s}$ to ${\vec s}_r$.

The boundary intensity $I_{w,a}({\vec s}_r)$ in the outgoing direction ${\vec s}_r$ on side $a$ of the interface is determined from the reflected component of the incoming radiation and the transmission from side $b$. Thus

\begin{displaymath} I_{w,a}({\vec s}_r) = r_a({\vec s}) I_{w,a}({\vec s}) + \tau_b({\vec s} \; ')I_{w,b}({\vec s} \; ') \end{displaymath} (11.3-54)

where $\tau_b({\vec s} \; ')$ is the transmissivity of side $b$ in direction ${\vec s} '$. Similarly, the outgoing intensity in the direction ${\vec s}_t$ on side $b$ of the interface, $I_{w,b}({\vec s}_t)$, is given by

\begin{displaymath} I_{w,b}({\vec s}_t) = r_b({\vec s} \; ') I_{w,b}({\vec s} \; ') + \tau_a({\vec s})I_{w,a}({\vec s}) \end{displaymath} (11.3-55)

For the case $n_a < n_b$, the energy transmitted from medium $a$ to medium $b$ in the incoming solid angle $2\pi$ must be refracted into a cone of apex angle $\theta_c$ (see Figure  11.3.9) where

\begin{displaymath} \theta_c = \sin^{-1} \frac{n_a}{n_b} \end{displaymath} (11.3-56)

Figure 11.3.9: Critical Angle $\theta_c$
\begin{figure} \psfig{file=figures/fig-theta-c.ps,width=3in} \end{figure}

Similarly, the transmitted component of the radiant energy going from medium $b$ to medium $a$ in the cone of apex angle $\theta_c$ is refracted into the outgoing solid angle $2\pi$. For incident angles greater than $\theta_c$, total internal reflection occurs and all the incoming energy is reflected specularly back into medium $b$.

When medium $b$ is external to the domain, $I_{w,b}({\vec s} \; ')$ is given in Equation  11.3-54 as a part of the problem specification. This boundary specification is usually made by providing the incoming radiative flux and the solid angle over which the radiative flux is to be applied. The refractive index of the external medium is assumed to be unity.

Diffuse Semi-Transparent Walls

In many engineering problems, the semi-transparent interface may be a diffuse reflector. For such a case, the interfacial reflectivity $r({\vec s})$ is assumed independent of ${\vec s}$, and equal to the hemispherically averaged value $r_d$. For $n = n_a/n_b > 1$, $r_{d,a}$ and $r_{d,b}$ are given by [ 235]

$\displaystyle r_{d,a}$ $\textstyle =$ $\displaystyle 1 - \frac{(1-r_{d,b})}{n^2}$ (11.3-57)
$\displaystyle r_{d,b}$ $\textstyle =$ $\displaystyle \frac{1}{2} + \frac{( 3n+1)(n-1)}{6(n+1)^2} + \frac{n^2(n^2-1)^2}{(n^2+1)^3} \ln \left(\frac{n-1}{n+1} \right) - \; \;$  
    $\displaystyle \frac{2n^3(n^2 + 2n -1)}{(n^2+1)(n^4-1)} + \frac{8n^4(n^4+1)}{(n^2+1)(n^4-1)^2} \ln(n)$ (11.3-58)

The boundary intensity for all outgoing directions on side $a$ of the interface is given by

\begin{displaymath} I_{w,a} = \frac{r_{d,a} q_{{\rm in},a} + \tau_{d,b} q_{{\rm in},b}}{\pi} \end{displaymath} (11.3-59)

Similarly for side $b$,

\begin{displaymath} I_{w,b} = \frac{r_{d,b} q_{{\rm in},b} + \tau_{d,a} q_{{\rm in},a}}{\pi} \end{displaymath} (11.3-60)


$\displaystyle q_{{\rm in},a}$ $\textstyle =$ $\displaystyle -\int_{4\pi}I_{w,a} {\vec s} \cdot {\vec n} d\Omega, \; \; \; {\vec s}\cdot {\vec n} <0$ (11.3-61)
$\displaystyle q_{{\rm in},b}$ $\textstyle =$ $\displaystyle \int_{4\pi}I_{w,b} {\vec s} \cdot {\vec n} d\Omega, \; \; \; {\vec s}\cdot {\vec n} \ge 0$ (11.3-62)

As before, if medium $b$ is external to the domain, $q_{{\rm in},b}$ is given as a part of the boundary specification.

Beam Irradiation

As mentioned above, FLUENT allows the specification of the irradiation at semi-transparent boundaries. The irradiation is specified in terms of an incident radiant heat flux (W/m $^2$). You can specify the solid angle over which the irradiation is distributed, as well as the vector of the centroid of the solid angle. To indicate whether the irradiation is reflected specularly or diffusely, you can specify the diffuse fraction.

For non-gray radiation, FLUENT allows you to specify the irradiation at semi-transparent boundaries on a band basis. The irradiation is specified as an incident heat flux (W/m $^2$) for each wavelength band. As in the gray case, you can specify the solid angle over which the irradiation is distributed, as well as the vector of the centroid of the solid angle.

The Diffuse Fraction

At semi-transparent boundaries, FLUENT allows you to specify the fraction of the incoming radiation that is treated as diffuse. The diffuse fraction is reflected diffusely, using the treatment described above; the transmitted portion is also treated diffusely. The remainder of the incoming energy is treated in a specular fashion.

For non-gray radiation, you can specify the diffuse fraction separately for each band.

Boundary Condition Treatment at Specular Walls and Symmetry Boundaries

At specular walls and symmetry boundaries, the direction of the reflected ray ${\vec s}_r$ corresponding to the incoming direction ${\vec s}$ is given by Equation  11.3-50. Furthermore,

\begin{displaymath} I_w({\vec s}_r) = I_w({\vec s}) \end{displaymath} (11.3-63)

Boundary Condition Treatment at Periodic Boundaries

When rotationally periodic boundaries are used, it is important to use pixelation in order to ensure that radiant energy is correctly transferred between the periodic and shadow faces. A pixelation between $3\times 3$ and $10\times 10$ is recommended.

Boundary Condition Treatment at Flow Inlets and Exits

The treatment at flow inlets and exits is described in Section  11.3.3.

next up previous contents index Previous: 11.3.5 The Rosseland Radiation
Up: 11.3 Radiative Heat Transfer
Next: 11.3.7 The Surface-to-Surface (S2S)
© Fluent Inc. 2003-01-25