3.2.2. Laminar Flow Equations

Laminar transport equations are not used for fire problems, but they are important for other classes of problems such as manufacturing. The low Mach number approximation is assumed (see Low Mach Number Equations).

3.2.2.1. Conservation of Mass

The mass conservation equation of a mixture of gases is given by

(3.25) $\int {{\partial \rho} \over {\partial t}} {\rm d}V + \int \rho u_j n_j {\rm d}S = 0,$

where $u_j$ is the mass average velocity of the mixture [6].

3.2.2.2. Conservation of Momentum

The conservation of momentum equations are given by

(3.26) $\int {{\partial \rho u_i} \over {\partial t}} {\rm d}V + \int \rho u_i u_j n_j {\rm d}S + \int P n_i {\rm d}S = \int \tau_{ij} n_j {\rm d}S + \int \left(\rho - \rho_{\circ} \right) g_i {\rm d}V ,$

where the viscous stress tensor is

(3.27) $\tau_{ij} = \mu \left( {{\partial u_i} \over {\partial x_j}} + {{\partial u_j} \over {\partial x_i}} \right) - {2 \over 3} \mu {{\partial u_k} \over {\partial x_k}} \delta_{ij} .$

The pressure, $P$ , in the momentum equations deserves a special note as this quantity can represent either the dynamic, i.e., the second term in the Mach number expansion in the case of the low Mach number assumption, or the static pressure in the case of formally compressibility. In either case, as shown above the hydrostatic pressure gradient has been removed which gives rise to the far-field density, $\rho_\circ$ , in the buoyancy body force. Optionally, we allow for the following sets of buoyancy models:

Boussinesq buoyancy approximation where the density difference is approximated as

(3.28) $\left( \rho -\rho_\circ \right) \approx - {{\rho_\circ} \over {T_\circ}} \left( T - T_\circ \right),$

Standard buoyant model in which case the pressure above does include the hydrostatic pressure and the buoyancy right-hand-side source term is,

(3.29) $\rho g_i ,$

A Boussinesq approximation for a binary mixture in which case the right-hand-side contribution is:

(3.30) $\rho MW^{ref}\left( {1\over MW_1 } - {1\over MW_2 }\right)\left[Y_1 - Y^{ref}\right] g_i.$

The user is referred to the Fuego user manual for exact line commands for each of these buoyancy options.

Note that zero pressure is almost always a convenient initial condition for a low Mach fluid flow. However, in cases without buoyancy, it can be anything, as the value only defines the additive constant for the pressure solve. However, one must ensure that the value matches for both initial and boundary condition specifications.

For buoyant flow, specifying zero pressure is convenient in tandem with the “differential” buoyancy option. This buoyancy term subtracts off the hydrostatic contribution such that the source term is written as

(3.31) $g \left(\rho - \rho_\circ\right)$

One can see that using this term along with a zero pressure initial condition allows one to avoid specifying initial and boundary conditions as the hydrostatic pressure, i.e., as a function of height.

3.2.2.3. Conservation of Energy

The conservation of energy equation in terms of enthalpy (including a source term due to radiation absorption and emission) is

(3.32) $\int {{\partial \rho h} \over {\partial t}} {\rm d}V + \int \rho h u_j n_j {\rm d}S = \nonumber - \int q_j n_j {\rm d}S - \int {{\partial q_i^r} \over {\partial x_i}} {\rm d}V$

(3.33) $+ \int \left({{\partial P} \over {\partial t}} + u_j {{\partial P} \over {\partial x_j}}\right){\rm d}V + \int \tau_{ij} {{\partial u_i} \over {\partial x_j }} {\rm d}V,$

where the energy diffusion flux vector is given by

(3.34) $q_j = - \kappa {{\partial T} \over {\partial x_j}} + \sum_{k=1}^K \rho h_k Y_k \hat{u}_{j,k} ,$

and $\hat{u}_{j,k}$ is the diffusion velocity of species $k$ in the $j$ direction.

This form of the energy equation is derived by starting with the energy equation and supplemental relationships of internal energy and total enthalpy provided in Asymptotic Expansion. The time term and convection term due to kinetic energy are expanded using the chain rule and simplified by enforcing the continuity equation. The remaining kinetic energy terms and gravitational force term are removed by dotting velocity with the momentum equation (to obtain the mechanical energy equation) and subtracting it from the energy equation. This procedure provides the full material derivative of pressure and the expanded viscous dissipation term.

The last two terms of (3.33) are only active when formal compressibility (in an acoustic sense) are important (see the Fuego user manual for the appropriate command lines to activate the low speed compressible and high speed compressible form in Fuego).

For a low Mach number flow, the time derivative of the pressure appearing above is substituted by the thermodynamic reference pressure, $P_{th}$ , that can only be nonzero in a closed volume with energy addition or subtraction. However, the low Mach number approximation mandates that the thermodynamic pressure is always spatially uniform.

The enthalpy of the mixture, $h$ , is a mass-average of the component enthalpies, $h_k$ , given by

(3.35) $h = \sum_{k=1}^{K} h_k Y_k .$

The energy diffusion flux vector includes a scaled gradient of temperature whereas the independent field to be solved in (3.33) is enthalpy. The form of the gradient of temperature is derived by first taking the gradient of (3.35) and using the chain rule,

(3.36) ${{\partial h} \over {\partial x_j}} = \sum_{k=1}^K Y_k {{\partial h_k} \over {\partial x_j}} + \sum_{k=1}^K h_k {{\partial Y_k} \over {\partial x_j}}.$

Given the thermodynamic definition of specific heat, the above equation is given by,

(3.37) ${{\partial h} \over {\partial x_j}} = \sum_{k=1}^K Y_k {C_p}_k {{\partial T} \over {\partial x_j}} + \sum_{k=1}^K h_k {{\partial Y_k} \over {\partial x_j}}$

(3.38) ${{\partial h} \over {\partial x_j}} = C_p {{\partial T} \over {\partial x_j}} + \sum_{k=1}^K h_k {{\partial Y_k} \over {\partial x_j}}.$

This equation is rearranged,

(3.39) ${{\partial T} \over {\partial x_j}} = {{1} \over {C_p}} \left( {{\partial h} \over {\partial x_j}} - \sum_{k=1}^K h_k {{\partial Y_k} \over {\partial x_j}} \right),$

and substituted into the energy diffusion flux vector to obtain,

(3.40) $q_j = - { {\kappa} \over {C_p}} \left({{\partial h} \over {\partial x_j}} - \sum_{k=1}^K h_k {{\partial Y_k} \over {\partial x_j}} \right) + \sum_{k=1}^K \rho h_k Y_k \hat{u}_{j,k}.$

Commonly, the last two terms in the above equation can be canceled when a simple diffusion model is assumed (see Conservation of Species, (3.46)) in the limit where the ratio of thermal and mass diffusion is equal (unity Lewis number, or equivalently speaking the Prandtl number equals the Schmidt number, i.e.,

(3.41) $Le^{unity} = {{Sc} \over {Pr}} = {{\alpha} \over {D}} = 1.$

For completeness, the thermal diffusivity, Prandtl and Schmidt number are defined by,

(3.42) $\alpha = {{\kappa} \over {\rho C_p}},$

(3.43) $Pr = {{C_p \mu } \over {\kappa}} = {{\mu} \over {\rho \alpha}}.$

and

(3.44) $Sc = {{\mu } \over {\rho D_{ab}}}.$

3.2.2.4. Conservation of Species

The mass conservation equation for species $k$ in a mixture of $K$ gas phase species is

(3.45) $\int {{\partial \rho Y_k} \over {\partial t}} {\rm d}V + \int \rho Y_k u_j n_j {\rm d}S = - \int \rho \hat{u}_{j,k} Y_k n_j {\rm d}S + \int \dot{\omega}_k {\rm d}V,$

where $\dot{\omega}_k$ is the mass generation rate of species $k$ per unit volume by homogeneous chemical reactions. We allow several approximations for the diffusion velocity, $\hat{u}_{j,k}$ .

The simplest form is Fickian diffusion with the same value of mass diffusivity for all species,

(3.46) $\hat{u}_{i,k}= - D {1 \over {Y_k}} {{\partial Y_k} \over {\partial x_i}} .$

This form is used for the Reynolds-averaged form of the equations for turbulent flow.

A more accurate approximation uses a mixture-averaged diffusion coefficient, $\bar{D}_k$ , for each species diffusion velocity,

(3.47) $\hat{u}_{i,k} = - \bar{D}_k {1 \over {X_k}} {{\partial X_k} \over {\partial x_i}} = - \bar{D}_k \left( {1 \over {Y_k}} {{\partial Y_k} \over {\partial x_i}} + {1 \over {W}} {{\partial W} \over {\partial x_i}} \right) .$

However, the form above does not enforce the requirement that the sum over all species diffusional fluxes is equal to zero. To achieve this, we decompose the fluxes ( $j_{i,k}$ ) into a modeled form and a correction to the modeled form in terms of the diffusional velocities

(3.48) $j_{i,k} = \rho Y_k \left(\hat{u}^m_{i,k} + \hat{u}^c_{i} \right)$

where $\hat{u}^m_{i,k}$ is defined by (3.47). To define the correction velocity $\hat{u}^c_{i}$ we sum over all species and apply the zero flux condition to get

(3.49) $\sum\limits_{k=1}^K j_{i,k} = -\sum\limits_{k=1}^K \frac{\rho W_k \bar{D}_k}{W} \frac{\partial X_k}{\partial x_i} + \rho \hat{u}^c_{i} \sum\limits_{k=1}^K Y_k = 0$

which can be simplified to

(3.50) $\hat{u}^c_{i} = \sum\limits_{k=1}^K \frac{\rho W_k \bar{D}_k}{W} \frac{\partial X_k}{\partial x_i}$

3.2.2.5. Conservation of Momentum, Axisymmetric with Swirl

Axisymmetric flows, with or without swirl, are described by two-dimensional equations in cylindrical coordinates. All azimuthal derivatives are zero (i.e., $\partial / \partial \theta = 0$ ). The axial coordinate is $x$ , the radial coordinate is $r$ , and the azimuthal coordinate is $\theta$ . The radius is retained in the equations and the purpose will become more clear in the discussion of the discrete integral form. The axial velocity is $u$ , the radial velocity is $v$ , and the azimuthal velocity is $w$ .

Axial-Momentum

(3.51) ${{\partial \rho u r} \over {\partial t}} + {{\partial} \over {\partial x}} \left( \rho u^2 r \right) + {{\partial} \over {\partial r}} \left( \rho u v r \right) + r {{\partial p} \over {\partial x}} = {{\partial} \over {\partial x}} \left( r \tau_{xx} \right) + {{\partial} \over {\partial r}} \left( r \tau_{xr} \right) + \rho r g_x$

Radial-Momentum

(3.52) ${{\partial \rho v r} \over {\partial t}} + {{\partial} \over {\partial x}} \left( \rho u v r \right) + {{\partial} \over {\partial r}} \left( \rho v^2 r \right) + r {{\partial p} \over {\partial r}} - \rho w^2 = {{\partial} \over {\partial x}} \left( r \tau_{rx} \right) + {{\partial} \over {\partial r}} \left( r \tau_{rr} \right) - \tau_{\theta \theta} + \rho r g_r$

Azimuthal-Momentum

(3.53) ${{\partial \rho w r} \over {\partial t}} + {{\partial} \over {\partial x}} \left( \rho u w r \right) + {{\partial} \over {\partial r}} \left( \rho v w r \right) + \rho v w = {{\partial} \over {\partial x}} \left( r \tau_{\theta x} \right) + {1 \over r} {{\partial} \over {\partial r}} \left( r^2 \tau_{\theta r} \right)$

The viscous stress terms for the cylindrical equations are

(3.54) $\tau_{xx} = \mu \left[ 2 {{\partial u} \over {\partial x}} - {2 \over 3} \left( {{\partial u} \over {\partial x}} + {{\partial v} \over {\partial r}} + {v \over r} \right) \right]$

(3.55) $\tau_{rx} = \mu \left[ {{\partial v} \over {\partial x}} + {{\partial u} \over {\partial r}} \right]$

(3.56) $\tau_{rr} = \mu \left[ 2 {{\partial v} \over {\partial r}} - {2 \over 3} \left( {{\partial u} \over {\partial x}} + {{\partial v} \over {\partial r}} + {v \over r} \right) \right]$

(3.57) $\tau_{\theta\theta} = \mu \left[ 2 {v \over r} - {2 \over 3} \left( {{\partial u} \over {\partial x}} + {{\partial v} \over {\partial r}} + {v \over r} \right) \right]$

(3.58) $\tau_{r\theta} = \mu r {{\partial} \over {\partial r}} \left( {w \over r} \right)$

(3.59) $\tau_{x\theta} = \mu {{\partial w} \over {\partial x}}$

The azimuthal equation can be simplified by relating the swirl velocity to the angular velocity, $w=r\omega$ . The momentum equation, written in terms of the angular velocity, is

(3.60) ${{\partial \rho \omega r} \over {\partial t}} + {{\partial} \over {\partial x}} \left( \rho u \omega r \right) + {{\partial} \over {\partial r}} \left( \rho v \omega r \right) + 2 \rho v \omega = {{\partial} \over {\partial x}} \left( r \mu {{\partial \omega} \over {\partial x}} \right) + {{\partial} \over {\partial r}} \left( r \mu {{\partial \omega} \over {\partial r}} \right) + 2 \mu {{\partial \omega} \over {\partial r}}.$

The production term that is used in the turbulence model is

(3.61) $\Phi = 2 \left[ \left( {{\partial u} \over {\partial x}} \right)^2 + \left( {{\partial v} \over {\partial r}} \right)^2 + \left( {{v} \over {r}} \right)^2 \right] + \left( {{\partial u} \over {\partial r}} + {{\partial v} \over {\partial x}} \right)^2 - {2 \over 3} \left( {{\partial u} \over {\partial x}} + {{\partial v} \over {\partial r}} + {{v} \over {r}} \right)^2 .$

3.2.2.6. Laminar Flow Boundary Conditions

The laminar flow math models require boundary conditions for velocity, pressure, temperature and enthalpy variables, and mixture composition.

3.2.2.6.1. Inflow

There are three types of inflow boundary conditions. For velocity-specified inflow, Dirichlet conditions are applied to velocities in the momentum equations, temperature in the energy equation, and mass fractions in the species equations. The mass flow rate at the boundary is specified for the continuity equation. The pressure floats to a consistent value. Alternatively, a control volume balance is retained at the boundary nodes and the convection fluxes are specified.

For pressure-specified inflow, the outflow boundary condition is applied with the added condition that the flow must enter the domain normal to the mesh boundary. Transport equations are solved for the momentum, energy and species equations.

3.2.2.6.2. Outflow

The pressure is specified at integration points on the outflow boundary. The specified pressure is used in the surface integration procedure for approximation nodal gradients. The pressure gradients are used to construct an interpolation for the mass flow rate at the boundary. Transport equations are solved for the momentum, energy and species equations. Upwind extrapolation is used for the scalars if the flow is leaving the domain. The boundary values of velocity and specified far-field values of scalars are used if the flow is entering the domain.

3.2.2.6.3. Wall

It is assumed that there is no mass flow through the wall. The velocity is specified as a Dirichlet boundary condition in the momentum equations. The temperature is specified as a Dirichlet boundary condition in the energy if the wall is isothermal. We currently do not support heterogeneous chemical reactions at a surface, so there should be no boundary condition applied to the mass fractions.

3.2.2.6.4. Symmetry Plane

There is no mass flow rate through the symmetry plane and there is no transport of scalar variables. The normal stress (pressure and viscous) at the symmetry plane is applied in the momentum equations.