2. Supported Equation Set

This section provides an overview of the currently supported equation sets. Equations will be decribed in integral form with assumed Favre averaging. However, the laminar counterparts are supported in the code base and are obtain in the user file by ommitting a turbulence model specification.

2.1. Conservation of Mass

The continuity equation is always solved in the variable density form.

\[\int \frac{\partial \bar{\rho}} {\partial t}\, dV + \int \bar{\rho} \widetilde{u}_i n_i\, dS = 0\]

Since Nalu uses equal-order interpolation (variables are collocated) stabilization is required. The stabilization choice will be developed in the pressure stabilization section.

Note that the use of a low speed compressible formulation requires that the fluid density be computed by an equation of state that uses the thermodynamic pressure. This thermodynamic pressure can either be computed based on a global energy/mass balance or allowed to be spatially varying. By modification of the continuity density time derivative to include the \(\frac{\partial \rho}{\partial p}\) sensitivity, an equation that admits acoustic pressure waves is realized.

2.2. Conservation of Momentum

The integral form of the Favre-filtered momentum equations used for turbulent transport are

(2.1)\[\begin{split}\int \frac{\partial \bar{\rho} \widetilde{u}_i}{\partial t} \, {\rm d}V + \int \bar{\rho} \widetilde{u}_i \widetilde{u}_j n_j \, {\rm d}S = \int \widetilde{\sigma}_{ij} n_j \, {\rm d}S -\int \tau^{sgs}_{ij} n_j \, {\rm d}S \\ + \int \left(\bar{\rho} - \rho_{\circ} \right) g_i \, {\rm d}V + \int \mathrm{f}_i \, {\rm d}V,\end{split}\]

where the subgrid scale turbulent stress \(\tau^{sgs}_{ij}\) is defined as

(2.2)\[\tau^{sgs}_{ij} \equiv \bar{\rho} ( \widetilde{u_i u_j} - \widetilde{u}_i \widetilde{u}_j ).\]

The term \(\mathrm{f}_i\) is a body force used to represent additional momentum sources such as wind turbine blades, Coriolis effect, driving forces, etc. The Cauchy stress is provided by,

\[\sigma_{ij} = 2 \mu \widetilde S^*_{ij} - \bar P \delta_{ij}\]

and the traceless rate-of-strain tensor defined as follows:

\[\begin{split}\widetilde S^*_{ij} = \widetilde S_{ij} - \frac{1}{3} \delta_{ij} \widetilde S_{kk} \\ = \widetilde S_{ij} - \frac{1}{3} \frac{\partial \widetilde u_k }{\partial x_k}\delta_{ij}.\end{split}\]

In a low Mach flow, as described in the low Mach theory section, the above pressure, \(\bar P\) is the purturbation about the thermodynamic pressure, \(P^{th}\). In a low speed compressible flow, i.e., flows confined to a closed domain with energy or mass addition in which the continuity equation has been modifed to accomodate accoustics, this pressure is interpreted at the thermodynamic pressure itself.

For LES, \(\tau^{sgs}_{ij}\) that appears in Equation (2.1) and defined in Equation (2.2) represents the subgrid stress tensor that must be closed. The deviatoric part of the subgrid stress tensor is defined as

(2.3)\[\tau^{sgs}_{ij} = \tau^{sgs}_{ij} - \frac{1}{3} \delta_{ij} \tau^{sgs}_{kk}\]

where the subgrid turbulent kinetic energy is defined as \(\tau^{sgs}_{kk} = 2 \bar \rho k\). Note that here, k represents the modeled turbulent kinetic energy as is formally defined as,

\[\bar \rho k = \frac{1}{2} \bar\rho ( \widetilde{u_k u_k} - \widetilde u_k \widetilde u_k).\]

Model closures can use, Yoshizawa’s approach when k is not transported:

\[\tau^{sgs}_{kk} = 2 C_I \bar \rho \Delta^2 | \widetilde S | ^2.\]

Above, \(| \widetilde S | = \sqrt {2 \widetilde S_{ij} \widetilde S_{ij}}\).

For low Mach-number flows, a vast majority of the turbulent kinetic energy is contained at resolved scales. For this reason, the subgrid turbulent kinetic energy is not directly treated and, rather, is included in the pressure as an additional normal stress. The Favre-filtered momentum equations then become

(2.4)\[\begin{split}&\int \frac{\partial \bar{\rho} \widetilde{u}_i}{\partial t} {\rm d}V + \int \bar{\rho} \widetilde{u}_i \widetilde{u}_j n_j {\rm d}S + \int \left( \bar{P} + \frac{2}{3} \bar{\rho} k \right) n_i {\rm d}S = \\ & \int 2 (\mu + \mu_t) \left( \widetilde{S}_{ij} - \frac{1}{3} \widetilde{S}_{kk} \delta_{ij} \right) n_j {\rm d}S + \int \left(\bar{\rho} - \rho_{\circ} \right) g_i {\rm d}V,\end{split}\]

where LES closure models for the subgrid turbulent eddy viscosity \(\mu_t\) are either the constant coefficient Smagorinsky, WALE or the constant coefficient \(k_{sgs}\) model (see the turbulence section).

2.2.1. Earth Coriolis Force

For simulation of large-scale atmospheric flows, the following Coriolis force term can be added to the right-hand-side of the momentum equation ((2.1)):

(2.5)\[\mathrm{f}_i = -2\bar{\rho}\epsilon_{ijk}\Omega_ju_k .\]

Here, \(\Omega\) is the Earth’s angular velocity vector, and \(\epsilon_{ijk}\) is the Levi-Civita symbol denoting the cross product of the Earth’s angular velocity with the local fluid velocity vector. Consider an “East-North-Up” coordinate system on the Earth’s surface, with the domain centered on a latitude angle \(\phi\) (changes in latitude within the computational domain are neglected). In this coordinate system, the integrand of (cor-term), or the Coriolis acceleration vector, is

(2.6)\[\begin{split}2 \bar{\rho} \omega \begin{bmatrix} u_n \sin\phi - u_u \cos\phi \\ -u_e \sin\phi \\ u_e \cos\phi \end{bmatrix},\end{split}\]

where \(\omega \equiv ||\Omega||\). Often, in geophysical flows it is assumed that the vertical component of velocity is small and that the vertical component of the acceleration is small relative to gravity, such that the terms containing \(\cos\phi\) are neglected. However, there is evidence that this so-called traditional approximation is not valid for some mesoscale atmospheric phenomena cite{Gerkema_etal:08}, and so the full Coriolis term is retained in Nalu. The implementation proceeds by first finding the velocity vector in the East-North-Up coordinate system, then calculating the Coriolis acceleration vector ((2.6)), then transforming this vector back to the model \(x-y-z\) coordinate system. The coordinate transformations are made using user-supplied North and East unit vectors given in the model coordinate system.

2.2.2. Boussinesq Buoyancy Model

In atmospheric and other flows, the density differences in the domain can be small enough as to not significantly affect inertia, but nonetheless the buoyancy term,

(2.7)\[\int \left(\bar{\rho} - \rho_{\circ} \right) g_i ~{\rm d}V,\]

may still be important. The Boussinesq model ignores the effect of density on inertia while retaining the buoyancy term in Equation (2.1). For the purpose of evaluating the buoyant force, the density is approximated as

(2.8)\[\frac{\rho}{\rho_{\circ}} \approx 1 - \beta (T-T_{\circ}),\]

This leads to a buoyancy body force term that depends on temperature (\(T\)), a reference density (\(\rho_{\circ}\)), a reference temperature (\(T_{\circ}\)), and a thermal expansion coefficient (\(\beta\)) as

(2.9)\[\int -\rho_{\circ} \beta (T-T_{\circ}) g_i ~{\rm d}V.\]

The flow is otherwise kept as constant density.

2.3. Filtered Mixture Fraction

The optional quantity used to identify the chemical state is the mixture fraction, \(Z\). While there are many different definitions of the mixture fraction that have subtle variations that attempt to capture effects like differential diffusion, they can all be interpreted as a local mass fraction of the chemical elements that originated in the fuel stream. The mixture fraction is a conserved scalar that varies between zero in the secondary stream and unity in the primary stream and is transported in laminar flow by the equation,

(2.10)\[\frac{\partial \rho Z}{\partial t} + \frac{ \partial \rho u_j Z }{ \partial x_j} = \frac{\partial}{\partial x_j} \left( \rho D \frac{\partial Z}{\partial x_j} \right),\]

where \(D\) is an effective molecular mass diffusivity.

Applying either temporal Favre filtering for RANS-based treatments or spatial Favre filtering for LES-based treatments yields

(2.11)\[\int \frac{\partial \bar{\rho} \widetilde{Z}}{\partial t} {\rm d}V + \int \bar{\rho} \widetilde{u}_j \widetilde{Z} n_j {\rm d}S = - \int \tau^{sgs}_{Z,j} n_j {\rm d}S + \int \bar{\rho} D \frac{\partial \widetilde{Z}}{\partial x_j} n_j {\rm d}S,\]

where sub-filter correlations have been neglected in the molecular diffusive flux vector and the turbulent diffusive flux vector is defined as

\[\tau^{sgs}_{Z,j} \equiv \bar{\rho} \left( \widetilde{Z u_j} - \widetilde{Z} \widetilde{u}_j \right).\]

This subgrid scale closure is modeled using the gradient diffusion hypothesis,

\[\tau^{sgs}_{Z,j} = - \bar{\rho} D_t \frac{\partial Z}{\partial x_j},\]

where \(D_t\) is the turbulent mass diffusivity, modeled as \(\bar{\rho} D_t = \mu_t / \mathrm{Sc}_t\) where \(\mu_t\) is the modeled turbulent viscosity from momentum transport and \(\mathrm{Sc}_t\) is the turbulent Schmidt number. The molecular mass diffusivity is then expressed similarly as \(\bar{\rho} D = \mu / \mathrm{Sc}\) so that the final modeled form of the filtered mixture fraction transport equation is

\[\frac{\partial \bar{\rho} \widetilde{Z}}{\partial t} + \frac{ \partial \bar{\rho} \widetilde{u}_j \widetilde{Z} }{ \partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \frac{\mu}{\mathrm{Sc}} + \frac{\mu_t}{\mathrm{Sc}_t} \right) \frac{\partial \widetilde{Z}}{\partial x_j} \right].\]

In integral form the mixture fraction transport equation is

\[\int \frac{\partial \bar{\rho} \widetilde{Z}}{\partial t}\, dV + \int \bar{\rho} \widetilde{u}_j \widetilde{Z} n_j\, dS = \int \left( \frac{\mu}{\mathrm{Sc}} + \frac{\mu_t}{\mathrm{Sc}_t} \right) \frac{\partial \widetilde{Z}}{\partial x_j} n_j\, dS.\]

2.4. Conservation of Energy

The integral form of the Favre-filtered static enthalpy energy equation used for turbulent transport is

(2.12)\[\begin{split} \int \frac{\partial \bar{\rho} \widetilde{h}}{\partial t} {\rm d}V + \int \bar{\rho} \widetilde{h} \widetilde{u}_j n_j {\rm d}S &= - \int \bar{q}_j n_j {\rm d}S - \int \tau^{sgs}_{h,j} n_j {\rm d}S - \int \frac{\partial \bar{q}_i^r}{\partial x_i} {\rm d}V \\ &+ \int \left( \frac{\partial \bar{P}}{\partial t} + \widetilde{u}_j \frac{\partial \bar{P}}{\partial x_j} \right){\rm d}V + \int \overline{\tau_{ij} \frac{\partial u_i}{\partial x_j }} {\rm d}V + \int S_\theta {\rm d}V.\end{split}\]

The above equation is derived by starting with the total internal energy equation, subtracting the mechanical energy equation and enforcing the variable density continuity equation. Note that the above equation includes possible source terms due to thermal radiatitive transport, viscous dissipation, pressure work, and external driving sources (\(S_\theta\)).

Above, \(\frac{\partial \bar{q}_i^r}{\partial x_i}\) represents the divergence of the radiative flux and is active when participating media radiation (PMR) coupling is active. At present, there is no sophisticated absorption coefficient model supported in the gas phase. Therefore, this source term is fully applied to the enthalpy along with any possible linearizations that solely originate from the PMR realm solve and transfer.

The simple Fickian diffusion velocity approximation, Equation (2.22), is assumed, so that the mean diffusive heat flux vector \(\bar{q}_j\) is

\[\bar{q}_j = - \overline{ \left[ \frac{\kappa}{C_p} \frac{\partial h}{\partial x_j} - \frac{\mu}{\rm Pr} \sum_{k=1}^K h_k \frac{\partial Y_k} {\partial x_j} \right] } - \overline{ \frac{\mu}{\rm Sc} \sum_{k=1}^K h_k \frac{\partial Y_k}{\partial x_j} }.\]

If \(Sc = Pr\), i.e., unity Lewis number (\(Le = 1\)), then the diffusive heat flux vector simplifies to \(\bar{q}_j = -\frac{\mu}{\mathrm{Pr}} \frac{\partial \widetilde{h}}{\partial x_j}\). In the code base, the user has the ability to either specify a laminar Prandtl number, which is a constant, or provide a property evaluator for thermal conductivity. Inclusion of a Prandtl number prevails and ensures that the thermal conductivity is computed base on \(\kappa = \frac{C_p \mu}{Pr}\). The viscous dissipation term is closed by

\[\begin{split}\overline{\tau_{ij} \frac{\partial u_i}{\partial x_j }} &= \left(\left(\mu + \mu_t\right) \left( \frac{\partial \widetilde{u}_i}{\partial x_j} + \frac{\partial \widetilde{u}_j}{\partial x_i} \right) - \frac{2}{3} \left( \bar{\rho} \widetilde{k} + \mu_t \frac{\partial \widetilde{u}_k}{\partial x_k} \right) \delta_{ij} \right) \frac{\partial \widetilde{u}_i}{\partial x_j} \\ &= \left[ 2 \mu \widetilde{S}_{ij} + 2 \mu_t \left( \widetilde{S}_{ij} - \frac{1}{3} \widetilde{S}_{kk} \delta_{ij} \right) - \frac{2}{3} \bar{\rho} \widetilde{k} \delta_{ij} \right] \frac{\partial \widetilde{u}_i}{\partial x_j}.\end{split}\]

The subgrid scale turbulent flux vector \(\tau^{sgs}_{h}\) in Equation (2.12) is defined as

\[\tau_{h u_j} \equiv \bar{\rho} \left( \widetilde{h u_j} - \widetilde{h} \widetilde{u}_j \right).\]

As with species transport, the gradient diffusion hypothesis is used to close this subgrid scale model,

\[\tau^{sgs}_{h,j} = - \frac{\mu_t}{\mathrm{Pr}_t} \frac{\partial \widetilde{h}}{\partial x_j},\]

where \(\mathrm{Pr}_t\) is the turbulent Prandtl number and \(\mu_t\) is the modeled turbulent eddy viscosity from momentum closure. The resulting filtered and modeled turbulent energy equation is given by,

(2.13)\[\begin{split}\int \frac{\partial \bar{\rho} \widetilde{h}}{\partial t} {\rm d}V + \int \bar{\rho} \widetilde{h} \widetilde{u}_j n_j {\rm d}S &= \int \left( \frac{\mu}{\rm Pr} + \frac{\mu_t}{{\rm Pr}_t} \right) \frac{\partial \widetilde{h}}{\partial x_j} n_j {\rm d}S - \int \frac{\partial \bar{q}_i^r}{\partial x_i} {\rm d}V \\ &+ \int \left( \frac{\partial \bar{P}}{\partial t} + \widetilde{u}_j \frac{\partial \bar{P}}{\partial x_j}\right){\rm d}V + \int \overline{\tau_{ij} \frac{\partial u_j}{\partial x_j }} {\rm d}V.\end{split}\]

The turbulent Prandtl number must have the same value as the turbulent Schmidt number for species transport to maintain unity Lewis number.

2.5. Review of Prandtl, Schmidt and Unity Lewis Number

For situations where a single diffusion coefficient is applicable (e.g., a binary gas system) the Lewis number is defined as:

(2.14)\[{\rm Le} = \frac{\rm Sc}{\rm Pr} = \frac{\alpha}{D}.\]

If the diffusion rates of energy and mass are equal,

(2.15)\[{\rm Sc = Pr \ and \ Le = 1}.\]

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

(2.16)\[\alpha = \frac{\kappa}{\rho c_p},\]

(2.17)\[{\rm Pr} = \frac{c_p \mu }{\kappa} = \frac{\mu}{\rho \alpha},\]

and

(2.18)\[{\rm Sc} = \frac{\mu }{\rho D},\]

where \(c_p\) is the specific heat, \(\kappa\), is the thermal conductivity and \(\alpha\) is the thermal diffusivity.

2.6. Thermal Heat Conduction

For non-isothermal object response that may occur in conjugate heat transfer applications, a simple single material heat conduction equation is supported.

(2.19)\[\int \rho C_p \frac{\partial T} {\partial t} {\rm d}V + \int q_j n_j {\rm d}S = \int S {\rm d}V.\]

where \(q_j\) is again the energy flux vector, however, now in the following temperature form:

\[q_j = -\kappa \frac{\partial T}{\partial x_j}.\]

2.7. ABL Forcing Source Terms

In LES of wind plant atmospheric flows, it is often necessary to drive the flow to a predetermined vertical velocity and/or temperature profile. In Nalu, this is achieved by adding appropriate source terms \(\mathrm{f}_i\) to the momentum equation (2.1) and \(S_\theta\) to the enthalpy equation (2.12).

First, the momentum source term is discussed. The main objective of this implementation is to force the volume averaged velocity at a certain location to a specified value (\(<\mathrm{u}_i>=\mathrm{U}_i\)). The brackets used here, \(<>\), mean volume averaging over a certain region. In order to achieve this, a source term must be applied to the momentum equation. This source term can be better understood as a proportional controller within the momentum equation.

The velocity and density fields can be decomposed into a volume averaged component and fluctuations about that volume average as \(\mathrm{u}_i = \left< \mathrm{u}_i \right> + \mathrm{u}_i'\) and \(\bar{\rho} = \left< \bar{\rho} \right> + \bar{\rho}'\). A decomposition of the plane averaged momentum at a given instance in time is then

\[\left< \bar{\rho} \mathrm{u}_i \right> = \left< \bar{\rho} \right> \left< \mathrm{u}_i \right> + \left< \bar{\rho}' \mathrm{u}'_i \right>.\]

We now wish to apply a momentum source based on a desired spatial averaged velocity \(\mathrm{U}_i\). This can be expressed as:

\[\left< \bar{\rho} \mathrm{u}_i^* \right> = \left< \bar{\rho} \right> \left< \mathrm{u}^*_i \right> + \left< \bar{\rho}' {\mathrm{u}^*_i}' \right>,\]

where \(\mathrm{u}_i^*\) is an unknown reference velocity field whose volume average is the desired velocity \(\left< \mathrm{u}_i^* \right> = \mathrm{U}_i\). Since the correlation \(\left< \bar{\rho}' \mathrm{u^*}'_i \right>\) is unknown, we assume that

\[\left< \bar{\rho}' \mathrm{u^*}'_i \right> = \left< \bar{\rho}' \mathrm{u}'_i \right>\]

such that the momentum source can now be defined as:

(2.20)\[{\mathrm{f}_i} = \alpha_u \left( \, \frac{\left< \bar{\rho} \right> \mathrm{U_i} - \left< \bar{\rho} \right> \left< \mathrm{u}_i \right>} {\Delta t}\right)\]

where \(\left< \right>\) denotes volume averaging at a certain time \(t\), \(\mathrm{U}_i\) is the desired spatial averaged velocity, and \(\Delta t\) is the time-scale between when the source term is computed (time \(t\)) and when it is applied (time \(t + \Delta t\)). This is typically chosen to be the simulation time-step. In the case of an ABL simulation with flat terrain, the voulme averaging is done over an infinitesimally thin slice in the \(x\) and \(y\) directions, such that the body force is only a function of height \(z\) and time \(t\). The implementation allows the user to prescribe relaxation factors \(\alpha_u\) for the source terms that are applied. Nalu uses a default value of 1.0 for the relaxation factors if no values are defined in the input file during initialization.

The enthalpy source term works similarly to the momentum source term. A temperature difference is computed at every time-step and a forcing term is added to the enthalpy equation:

\[S_\theta = \alpha_\theta C_p \left( \frac{\theta_{\rm ref} - \left< \theta \right>}{\Delta t} \right)\]

where \(\theta_{\rm ref}\) is the desired spatial averaged temperature, \(\left< \theta \right>\) is the spatial averaged temperature, \(C_p\) is the heat capcity, \(\alpha_\theta\) is the relaxation factor, and \(\Delta t\) is the time-scale.

The present implementation can vary the source terms as a function of time and space using either a user-defined table of previously computed source terms (e.g., from a precursor simulation or another model such as WRF), or compute the source term as a function of the transient flow solution.

2.8. Conservation of Species

The integral form of the Favre-filtered species equation used for turbulent transport is

(2.21)\[\int \frac{\partial \bar{\rho} \widetilde{Y}_k}{\partial t} {\rm d}V + \int \bar{\rho} \widetilde{Y}_k \widetilde{u}_j n_j {\rm d}S = - \int \tau^{sgs}_{Y_k,j} n_j {\rm d}S - \int \overline{\rho Y_k \hat{u}_{j,k}} n_j {\rm d}S + \int \overline{\dot{\omega}_k} {\rm d}V,\]

where the form of diffusion velocities (see Equation (2.22)) assumes the Fickian approximation with a constant value of diffusion velocity for consistency with the turbulent form of the energy equation, Equation (2.12). The simplest form is Fickian diffusion with the same value of mass diffusivity for all species,

(2.22)\[\hat{u}_{j,k}= - D \frac{1}{Y_k} \frac{\partial Y_k}{\partial x_j} .\]

The subgrid scale turbulent diffusive flux vector \(\tau^{sgs}_{Y_kj}\) is defined as

\[\tau^{sgs}_{Y_k,j} \equiv \bar{\rho} \left( \widetilde{Y_k u_j} - \widetilde{Y_k} \widetilde{u}_j \right).\]

The closure for this model takes on its usual gradient diffusion hypothesis, i.e.,

\[\tau^{sgs}_{Y_k,j} = - \frac{\mu_t}{\mathrm{Sc}_t} \frac{\partial \widetilde{Y}_k}{\partial x_j},\]

where \(\mathrm{Sc}_t\) is the turbulent Schmidt number for all species and \(\mu_t\) is the modeled turbulent eddy viscosity from momentum closure.

The Favre-filtered and modeled turbulent species transport equation is,

(2.23)\[\int \frac{\partial \bar{\rho} \widetilde{Y}_k}{\partial t} {\rm d}V + \int \bar{\rho} \widetilde{Y}_k \widetilde{u}_j n_j {\rm d}S = \int \left( \frac{\mu}{\rm Sc} + \frac{\mu_t}{{\rm Sc}_t} \right) \frac{\partial \widetilde{Y}_k}{\partial x_j} n_j {\rm d}S + \int \overline{\dot{\omega}}_k {\rm d}V .\]

If transporting both energy and species equations, the laminar Prandtl number must be equal to the laminar Schmidt number and the turbulent Prandtl number must be equal to the turbulent Schmidt number to maintain unity Lewis number. Although there is a species conservation equation for each species in a mixture of \(n\) species, only \(n-1\) species equations need to be solved since the mass fractions sum to unity and

\[\widetilde{Y}_n = 1 - \sum_{j \ne n}^{n} \widetilde{Y}_j .\]

Finally, the reaction rate source term is expected to be added based on an operator split approach wherebye the set of ODEs are integrated over the full time step. The chemical kinetic source terms can be sub-integrated within a time step using a stiff ODE integrator package.

The following system of ODEs are numerically integrated over a time step \(\Delta t\) for a fixed temperature and pressure starting from the initial values of gas phase mass fraction and density,

\[\dot{Y}_k = \frac{\dot{\omega}_k \left( Y_k \right) }{\rho} \ .\]

The sources for the sub-integration are computed with the composition and density at the new time level which are used to approximate a mean production rate for the time step

\[\dot{\omega}_k \approx \frac{\rho^{\ast} Y^{\ast}_k - \rho Y_k}{\Delta t} \ .\]

2.9. Subgrid-Scale Kinetic Energy One-Equation LES Model

The subgrid scale kinetic energy one-equation turbulence model, or \(k^{sgs}\) model, [Dav97], represents a simple LES closure model. The transport equation for subgrid turbulent kinetic energy is given by

(2.24)\[\int \frac{\partial \bar{\rho}{k^\mathrm{sgs}}}{\partial t} {\rm d}V + \int \bar{\rho}{k^\mathrm{sgs}} \widetilde{u}_j n_j {\rm d}S = \int \frac{\mu_t}{\sigma_k} \frac{\partial {k^\mathrm{sgs}}}{\partial x_j} n_j {\rm d}S + \int \left(P_k^\mathrm{sgs} - D_k^\mathrm{sgs}\right) {\rm d}V.\]

The production of subgrid turbulent kinetic energy, \(P_k^\mathrm{sgs}\), is modeled by,

(2.25)\[P_k \equiv -\overline{\rho u_i'' u_j''} \frac{\partial \widetilde{u}_i}{\partial x_j},\]

while the dissipation of turbulent kinetic energy, \(D_k^\mathrm{sgs}\), is given by

\[D_k^\mathrm{sgs} = \rho C_{\epsilon} \frac{{k^\mathrm{sgs}}^{\frac{3}{2}}}{\Delta},\]

where the grid filter length, \(\Delta\), is given in terms of the grid cell volume by

\[\Delta = V^{\frac{1}{3}}.\]

The subgrid turbulent eddy viscosity is then provided by

\[\mu_t = C_{\mu_{\epsilon}} \Delta {k^\mathrm{sgs}}^{\frac{1}{2}},\]

where the values of \(C_{\epsilon}\) and \(C_{\mu_{\epsilon}}\) are 0.845 and 0.0856, respectively.

For simulations in which a buoyancy source term is desired, the code supports the Rodi form,

\[P_b = \beta \frac{\mu^T}{Pr} g_i \frac{\partial T}{\partial x_i}.\]

2.10. Shear Stress Transport (SST) RANS Model Suite

Although Nalu is primarily expected to be a LES simulation tool, RANS modeling is supported through the activation of the SST equation set.

It has been observed that standard 1998 \(k-\omega\) models display a strong sensitivity to the free stream value of \(\omega\) (see Mentor, [MKL03]). To remedy, this, an alternative set of transport equations have been used that are based on smoothly blending the \(k-\omega\) model near a wall with \(k-\epsilon\) away from the wall. Because of the relationship between \(\omega\) and \(\epsilon\), the transport equations for turbulent kinetic energy and dissipation can be transformed into equations involving \(k\) and \(\omega\). Aside from constants, the transport equation for \(k\) is unchanged. However, an additional cross-diffusion term is present in the \(\omega\) equation. Blending is introduced by using smoothing which is a function of the distance from the wall, \(F(y)\). The transport equations for the Mentor 2003 model are then

\[\int \frac{\partial \bar{\rho} k}{\partial t} \text{d}V + \int \bar{\rho} k\widetilde{u}_{j} n_{j} \text{d} S = \int {(\mu + \hat \sigma_k \mu_{t})} \frac{\partial k}{\partial x_{j}} n_{j} + \int \left(P_{k}^{\omega} - \beta^* \bar{\rho} k \omega\right) \text{d} V,\]

\[\begin{split}\int \frac{\partial \bar{\rho} \omega}{\partial t}\text{d} V + \int \bar{\rho} \omega \widetilde{u}_{j} n_{j} \text{d}S = \int {(\mu + \hat\sigma_{\omega} \mu_{t})} \frac{\partial \omega}{\partial x_{j}} n_{j} + \int {2(1-F) \frac{\bar{\rho}\sigma_{\omega2}} {\omega} \frac{\partial k}{\partial x_j} \frac{\partial \omega}{\partial x_j} } \text{d}V \\ + \int \left(\frac{\hat\gamma}{\nu_t} P_{k}^{\omega} - \hat \beta \bar{\rho} \omega^{2}\right) \text{d}V.\end{split}\]

The model coefficients, \(\hat\sigma_k\), \(\hat\sigma_{\omega}\), \(\hat\gamma\) and \(\hat\beta\) must also be blended, which is represented by

\[\hat \phi = F\phi_1+ (1-F)\phi_2.\]

where \(\sigma_{k1} = 0.85\), \(\sigma_{k2} = 1.0\), \(\sigma_{\omega1} = 0.5\), \(\sigma_{\omega2} = 0.856\), \(\gamma_1 = \frac{5}{9}\), \(\gamma_2 = 0.44\), \(\beta_1 = 0.075\) and \(\beta_2 = 0.0828\). The blending function is given by

\[F = \tanh(arg_{1}^{4}),\]

where

\[arg_{1} = \min \left( \max \left( \frac{\sqrt{k}}{\beta^* \omega y}, \frac{500 \mu}{\bar{\rho} y^{2} \omega} \right), \frac{4 \bar{\rho} \sigma_{\omega2} k}{CD_{k\omega} y^{2}} \right).\]

The final parameter is

\[CD_{k\omega} = \max \left( 2 \bar{\rho} \sigma_{\omega2} \frac{1}{\omega} \frac{\partial k}{\partial x_{j}} \frac{\partial \omega}{\partial x_{j}}, 10^{-10} \right).\]

An important component of the SST model is the different expression used for the turbulent viscosity,

\[\mu_{t} = \frac {a_1 \bar{\rho} k} {\max\left( a_1 \omega, S F_2 \right) },\]

where \(F_2\) is another blending function given by

\[F_2 = \tanh(arg_{2}^{2}).\]

The final parameter is

\[arg_{2} = \max\left( \frac{2 \sqrt{k}}{\beta^* \omega y}, \frac{500 \mu}{\bar{\rho} \omega y^{2}} \right).\]

2.11. Direct Eddy Simulation (DES) Formulation

The DES technique is also supported in the code base when the SST model is activated. This model seeks to formally relax the RANS-based approach and allows for a theoretical basis to allow for transient flows. The method follows the method of Temporally Filtered NS formulation as decribed by Tieszen, [TDB05].

The SST DES model simply changes the turbulent kinetic energy equation to include a new minimum scale that manipulates the dissipation term.

\[D_k = \frac{\rho k^{3/2}} {l_{DES}},\]

where \(l_{DES}\) is the min(\(l_{SST}, c_{DES}l_{DES}\)). The constants are given by, \(l_{SST}=\frac{k^{1/2}}{\beta^* \omega}\) and \(c_{DES}\) represents a blended set of DES constants: \(c_{{DES}_1} = 0.78\) and \(c_{{DES}_2} = 0.61\). The length scale, \(l_{DES}\) is the maximum edge length scale touching a given node.

Note that the production term appearing in both the turbulent kinetic energy and specific dissipation rate equation is limited to a user-supplied scaling of the above dissipation term. Currently, all SST variants use a production to dissipation ratio limiting of ten.

2.12. Solid Stress

A fully implicit CVFEM (only) linear elastic equation is supported in the code base. This equation is either used for true solid stress prediction or for smoothing the mesh due to boundary mesh motion (either through fluid structure interaction (FSI) or prescribed mesh motion).

Consider the displacement for component i, \(u_i\) equation set,

(2.26)\[\rho \frac{\partial^2 u_i} {{\partial t}^2} - \frac{\partial \sigma_{ij}}{\partial x_j} = F_i,\]

where the Cauchy stress tensor, \(\sigma_{ij}\) assuming Hooke’s law is given by,

(2.27)\[\sigma_{ij} = \mu \left ( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) + \lambda \frac{\partial u_k}{\partial x_k} \delta_{ij}.\]

Above, the so-called Lame coefficients, Lame’s first parameter, \(\lambda\) (also known as the Lame modulus) and Lame’s second parameter, \(\mu\) (also known as the shear modulus) are provided as functions of the Young’s modulus, \(E\), and Poisson’s ratio, \(\nu\); here shown in the context of a isotropic elastic material,

(2.28)\[\mu = \frac{E}{2\left(1+\nu\right)},\]

and

(2.29)\[\lambda = \frac{E \nu}{\left(1+\nu\right) \left(1-2 \nu \right)}.\]

Note that the above notation of \(u_i\) to represent displacement is with respect to the classic definition of current and model coordinates,

(2.30)\[x_i = X_i + u_i\]

where \(x_i\) is the position, relative to the fixed, or previous position, \(X_i\).

The above equations are solved for mesh displacements, \(u_i\). The supplemental relationship for solid velocity, \(v_i\) is given by,

(2.31)\[v_i = \frac{\partial u_i}{\partial t}.\]

Numerically, the velocity might be obtained by a backward Euler or BDF2 scheme,

(2.32)\[v_i = \frac{\gamma_1 u^{n+1}_i + \gamma_2 u^n_i + \gamma_3 u^{n-1}_i}{\Delta t}\]

2.13. Moving Mesh

The code base supports three notions of moving mesh: 1) linear elastic equation system that computes the stress of a solid 2) solid body rotation mesh motion and 3) mesh deformation via an external source.

The linear elastic equation system is activated via the standard equation system approach. Properties for the solid are specified in the material block. Mesh motion is prescribed by the input file via the mesh_motion block. Here, it is assumed that the mesh motion is solid rotation. For fluid/structure interaction (FSI) a mesh smoothing scheme is used to propagate the surface mesh displacement obtained by the solids solve. Simple mesh smoothing is obtained via a linear elastic solve in which the so-called Lame constants are proportional to the inverse of the dual volume. This allows for boundary layer mesh locations to be stiff while free stream mesh elements to be soft.

Additional mesh motion terms are required for the Eulerian fluid mechanics solve. Using the geometric conservative law the time and advection source term for a general scalar \(\phi\) can be written as:

(2.33)\[\int \frac {\rho \phi } {\partial t}\, dV + \int \rho \phi \left ( u_j - v_j \right) n_j\, dS + \int \rho \phi \frac{\partial v_k}{\partial x_j}\, dV,\]

where \(u_j\) is the fluid velocity and \(v_j\) is the mesh velocity. Mesh velocities and the mesh velocity spatial derivatives are provided by the mesh smoothing solve. Activating the external mesh deformation or mesh motion block will result in the velocity relative to mesh calculation in the advection terms. The line command for source term, “\(gcl\)” must be activated for each equation for the volume integral to be included in the set of PDE solves. Finally, transfers are expected between the physics. For example, the solids solve is to provide mesh displacements to the mesh smoothing realm. The mesh smoothing procedure provides the boundary velocity, mesh velocity and projected nodal gradients of the mesh velocity to the fluids realm. Finally, the fluids solve is to provide the surface force at the desired solids surface. Currently, the pressure is transfered from the fluids realm to the solids realm. The ideal view of FSI is to solve the solids pde at the half time step. As such, in time, the \(P^{n+\frac{1}{2}}\) is expected. The fsi_interface input line command attribute is expected to be set at these unique surfaces. More advanced FSI coupling techniques are under development by a current academic partner.

2.14. Radiative Transport Equation

The spatial variation of the radiative intensity corresponding to a given direction and at a given wavelength within a radiatively participating material, \(I(s)\), is governed by the Boltzmann transport equation. In general, the Boltzmann equation represents a balance between absorption, emission, out-scattering, and in-scattering of radiation at a point. For combustion applications, however, the steady form of the Boltzmann equation is appropriate since the transient term only becomes important on nanosecond time scales which is orders of magnitude shorter than the fastest chemical.

Experimental data shows that the radiative properties for heavily sooting, fuel-rich hydrocarbon diffusion flames (\(10^{-4}\)% to \(10^{-6}\)% soot by volume) are dominated by the soot phase and to a lesser extent by the gas phase. Since soot emits and absorbs radiation in a relatively constant spectrum, it is common to ignore wavelength effects when modeling radiative transport in these environments. Additionally, scattering from soot particles commonly generated by hydrocarbon flames is several orders of magnitude smaller that the absorption effect and may be neglected. Moreover, the phase function is rarely known. However, for situations in which the phase function can be approximated by the iso-tropic scattering assumption, i.e., an intensity for direction \(k\) has equal probability to be scattered in any direction \(l\), the appropriate form of the Botzmann radiative transport is

(2.34)\[s_i \frac{\partial}{\partial x_i} I\left(s\right) + \left(\mu_a + \mu_s \right) I\left(s\right) = \frac{\mu_a \sigma T^4}{\pi} + \frac{\mu_s}{4\pi}G,\]

where \(\mu_a\) is the absorption coeffiecient, \(\mu_s\) is the scattering coefficeint, \(I(s)\) is the intensity along the direction \(s_i\), \(T\) is the temperature and the scalar flux is \(G\). 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 liniear, decoupled equations for each intensity to be solved.

The flux divergence may be written as a difference between the radiative emission and mean incident radiation at a point,

(2.35)\[\frac{\partial q_i^r}{\partial x_i} = \mu_a \left[ 4 \sigma T^4 - G \right] ,\]

where \(G\) is again the scalar flux. The flux divergence term is the same regardless of whether or not scattering is active. The quantity, \(G/4\pi\), is often referred to as the mean incident intensity. Note that when the scattering coefficient is non-zero, the RTE is coupled over all intensity directions by the scalar flux relationship.

The scalar flux and radiative flux vector represent angular moments of the directional radiative intensity at a point,

\[G = \int_{0}^{2\pi}\!\int_{0}^{\pi}\! I\left(s\right) \sin \theta_{zn} d \theta_{zn} d \theta_{az} ,\]

\[q^{r}_{i} = \int_{0}^{2\pi}\!\int_{0}^{\pi}\! I\left(s\right) s_i \sin \theta_{zn} d \theta_{zn} d \theta_{az} ,\]

where \(\theta_{zn}\) and \(\theta_{az}\) are the zenith and azimuthal angles respectively as shown in Figure Fig. 2.1.

Fig. 2.1 Ordinate Direction Definition, \({\bf s} = \sin \theta_{zn} \sin \theta_{az} {\bf i} + \cos \theta_{zn} {\bf j} + \sin \theta_{zn} \cos \theta_{az} {\bf k}\).

The radiation intensity must be defined at all portions of the boundary along which \(s_i n_i < 0\), where \(n_i\) is the outward directed unit normal vector at the surface. The intensity is applied as a weak flux boundary condition which is determined from the surface properties and temperature. The diffuse surface assumption provides reasonable accuracy for many engineering combustion applications. The intensity leaving a diffuse surface in all directions is given by

(2.36)\[I\left(s\right) = \frac{1}{\pi} \left[ \tau \sigma T_\infty^4 + \epsilon \sigma T_w^4 + \left(1 - \epsilon - \tau \right) H \right] ,\]

where \(\epsilon\) is the total normal emissivity of the surface, \(\tau\) is the transmissivity of the surface, \(T_w\) is the temperature of the boundary, \(T_\infty\) is the environmental temperature and \(H\) is the incident radiation, or irradiation (incoming radiative flux). Recall that the relationship given by Kirchoff’s Law that relates emissivity, transmissivity and reflectivity, \(\rho\), is

\[\rho + \tau + \epsilon = 1.\]

where it is implied that \(\alpha = \epsilon\).