2. Supported Equation Set

This section provides an overview of the currently supported equation sets. Equations will be described 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 omitting a turbulence model specification.

2.1. Conservation of Mass

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

$$\int \frac{\partial \overline{\rho }}{\partial t}dV+\int \overline{\rho }{\tilde{u}}_i{n}_{i}dS=0$$

Since Nalu-Wind 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{array}{r}\int \frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial t}\mathrm{d}V+\int \overline{\rho }{\tilde{u}}_i{\tilde{u}}_j{n}_j\mathrm{d}S=\int {\tilde{\sigma }}_{ij}{n}_j\mathrm{d}S-\int {\tau }_{ij}^{sgs}{n}_j\mathrm{d}S\\ +\int (\overline{\rho }-{\rho }_{\circ })g_i\mathrm{d}V+\int {\mathrm{f}}_i\mathrm{d}V,\end{array}$$

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

(2.2)$${\tau }_{ij}^{sgs}\equiv \overline{\rho }(\widetilde{u_i{u}_j}-{\tilde{u}}_i{\tilde{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 {\tilde{S}}_{ij}^{\ast }-\overline{P}{\delta }_{ij}$$

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

$$\begin{array}{r}{\tilde{S}}_{ij}^{\ast }={\tilde{S}}_{ij}-\frac{1}{3}{\delta }_{ij}{\tilde{S}}_{kk}\\ ={\tilde{S}}_{ij}-\frac{1}{3}\frac{\partial {\tilde{u}}_k}{\partial x_k}{\delta }_{ij}.\end{array}$$

In a low Mach flow, as described in the low Mach theory section, the above pressure, $\overline{P}$ is the perturbation 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 modified to accommodate acoustics, this pressure is interpreted at the thermodynamic pressure itself.

For LES, ${\tau }_{ij}^{sgs}$ 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 }_{ij}^{sgs}={\tau }_{ij}^{sgs}-\frac{1}{3}{\delta }_{ij}{\tau }_{kk}^{sgs}$$

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

$$\overline{\rho }k=\frac{1}{2}\overline{\rho }(\widetilde{u_k{u}_k}-{\tilde{u}}_k{\tilde{u}}_k).$$

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

$${\tau }_{kk}^{sgs}=2C_I\overline{\rho }{\mathrm{\Delta }}^2|\tilde{S}{|}^2.$$

Above, $|\tilde{S}|=\sqrt{2{\tilde{S}}_{ij}{\tilde{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{array}{rl}& \int \frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial t}\mathrm{d}V+\int \overline{\rho }{\tilde{u}}_i{\tilde{u}}_j{n}_j\mathrm{d}S+\int (\overline{P}+\frac{2}{3}\overline{\rho }k)n_i\mathrm{d}S=\\ & \int 2(\mu +{\mu }_t)({\tilde{S}}_{ij}-\frac{1}{3}{\tilde{S}}_{kk}{\delta }_{ij})n_j\mathrm{d}S+\int (\overline{\rho }-{\rho }_{\circ })g_i\mathrm{d}V,\end{array}$$

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\overline{\rho }{ϵ}_{ijk}{\mathrm{\Omega }}_j{u}_k.$$

Here, $\mathrm{\Omega }$ is the Earth’s angular velocity vector, and ${ϵ}_{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 $\varphi$ (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{array}{r}2\overline{\rho }\omega \left[\begin{array}{c}{u}_n\sin \varphi -u_u\cos \varphi \\ -u_e\sin \varphi \\ u_e\cos \varphi \end{array}\right],\end{array}$$

where $\omega \equiv ||\mathrm{\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 \varphi$ 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-Wind. 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 (\overline{\rho }-{\rho }_{\circ })g_i\mathrm{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\mathrm{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 \overline{\rho }\tilde{Z}}{\partial t}\mathrm{d}V+\int \overline{\rho }{\tilde{u}}_j\tilde{Z}{n}_j\mathrm{d}S=-\int {\tau }_{Z,j}^{sgs}{n}_j\mathrm{d}S+\int \overline{\rho }D\frac{\partial \tilde{Z}}{\partial x_j}{n}_j\mathrm{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 }_{Z,j}^{sgs}\equiv \overline{\rho }(\widetilde{Z{u}_j}-\tilde{Z}{\tilde{u}}_j).$$

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

$${\tau }_{Z,j}^{sgs}=-\overline{\rho }{D}_t\frac{\partial Z}{\partial x_j},$$

where $D_t$ is the turbulent mass diffusivity, modeled as $\overline{\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 $\overline{\rho }D=\mu /\mathrm{Sc}$ so that the final modeled form of the filtered mixture fraction transport equation is

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

In integral form the mixture fraction transport equation is

$$\int \frac{\partial \overline{\rho }\tilde{Z}}{\partial t}dV+\int \overline{\rho }{\tilde{u}}_j\tilde{Z}{n}_{j}dS=\int (\frac{\mu }{\mathrm{Sc}}+\frac{{\mu }_t}{{\mathrm{Sc}}_t})\frac{\partial \tilde{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{array}{rl}\int \frac{\partial \overline{\rho }\tilde{h}}{\partial t}\mathrm{d}V+\int \overline{\rho }\tilde{h}{\tilde{u}}_j{n}_j\mathrm{d}S& =-\int {\overline{q}}_j{n}_j\mathrm{d}S-\int {\tau }_{h,j}^{sgs}{n}_j\mathrm{d}S-\int \frac{\partial {\overline{q}}_i^r}{\partial x_i}\mathrm{d}V\\ & +\int (\frac{\partial \overline{P}}{\partial t}+{\tilde{u}}_j\frac{\partial \overline{P}}{\partial x_j})\mathrm{d}V+\int \overline{{\tau }_{ij}\frac{\partial u_i}{\partial x_j}}\mathrm{d}V+\int S_{\theta }\mathrm{d}V.\end{array}$$

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 }$).

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

$${\overline{q}}_j=-\overline{[\frac{\kappa }{C_p}\frac{\partial h}{\partial x_j}-\frac{\mu }{\mathrm{P}\mathrm{r}}\sum _{k=1}^K{h}_k\frac{\partial Y_k}{\partial x_j}]}-\overline{\frac{\mu }{\mathrm{S}\mathrm{c}}\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 ${\overline{q}}_j=-\frac{\mu }{\mathrm{Pr}}\frac{\partial \tilde{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{array}{rl}\overline{{\tau }_{ij}\frac{\partial u_i}{\partial x_j}}& =((\mu +{\mu }_t)(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i})-\frac{2}{3}(\overline{\rho }\tilde{k}+{\mu }_t\frac{\partial {\tilde{u}}_k}{\partial x_k}){\delta }_{ij})\frac{\partial {\tilde{u}}_i}{\partial x_j}\\ & =[2\mu {\tilde{S}}_{ij}+2{\mu }_t({\tilde{S}}_{ij}-\frac{1}{3}{\tilde{S}}_{kk}{\delta }_{ij})-\frac{2}{3}\overline{\rho }\tilde{k}{\delta }_{ij}]\frac{\partial {\tilde{u}}_i}{\partial x_j}.\end{array}$$

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

$${\tau }_{h{u}_j}\equiv \overline{\rho }(\widetilde{h{u}_j}-\tilde{h}{\tilde{u}}_j).$$

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

$${\tau }_{h,j}^{sgs}=-\frac{{\mu }_t}{{\mathrm{Pr}}_t}\frac{\partial \tilde{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{array}{rl}\int \frac{\partial \overline{\rho }\tilde{h}}{\partial t}\mathrm{d}V+\int \overline{\rho }\tilde{h}{\tilde{u}}_j{n}_j\mathrm{d}S& =\int (\frac{\mu }{\mathrm{P}\mathrm{r}}+\frac{{\mu }_t}{{\mathrm{P}\mathrm{r}}_t})\frac{\partial \tilde{h}}{\partial x_j}{n}_j\mathrm{d}S-\int \frac{\partial {\overline{q}}_i^r}{\partial x_i}\mathrm{d}V\\ & +\int (\frac{\partial \overline{P}}{\partial t}+{\tilde{u}}_j\frac{\partial \overline{P}}{\partial x_j})\mathrm{d}V+\int \overline{{\tau }_{ij}\frac{\partial u_j}{\partial x_j}}\mathrm{d}V.\end{array}$$

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)$$\mathrm{L}\mathrm{e}=\frac{\mathrm{S}\mathrm{c}}{\mathrm{P}\mathrm{r}}=\frac{\alpha }{D}.$$

If the diffusion rates of energy and mass are equal,

(2.15)$$\mathrm{S}\mathrm{c}=\mathrm{P}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{L}\mathrm{e}=1.$$

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

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

(2.17)$$\mathrm{P}\mathrm{r}=\frac{c_p\mu }{\kappa }=\frac{\mu }{\rho \alpha },$$

and

(2.18)$$\mathrm{S}\mathrm{c}=\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}\mathrm{d}V+\int q_j{n}_j\mathrm{d}S=\int S\mathrm{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-Wind, 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=⟨{\mathrm{u}}_i⟩+{\mathrm{u}}_i^{\prime }$ and $\overline{\rho }=⟨\overline{\rho }⟩+{\overline{\rho }}^{\prime }$. A decomposition of the plane averaged momentum at a given instance in time is then

$$⟨\overline{\rho }{\mathrm{u}}_i⟩=⟨\overline{\rho }⟩⟨{\mathrm{u}}_i⟩+⟨{\overline{\rho }}^{\prime }{\mathrm{u}}_i^{\prime }⟩.$$

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

$$⟨\overline{\rho }{\mathrm{u}}_i^{\ast }⟩=⟨\overline{\rho }⟩⟨{\mathrm{u}}_i^{\ast }⟩+⟨{\overline{\rho }}^{\prime }{{\mathrm{u}}_i^{\ast }}^{\prime }⟩,$$

where ${\mathrm{u}}_i^{\ast }$ is an unknown reference velocity field whose volume average is the desired velocity $⟨{\mathrm{u}}_i^{\ast }⟩={\mathrm{U}}_i$. Since the correlation $⟨{\overline{\rho }}^{\prime }{{\mathrm{u}}^{\ast }}_i^{\prime }⟩$ is unknown, we assume that

$$⟨{\overline{\rho }}^{\prime }{{\mathrm{u}}^{\ast }}_i^{\prime }⟩=⟨{\overline{\rho }}^{\prime }{\mathrm{u}}_i^{\prime }⟩$$

such that the momentum source can now be defined as:

(2.20)$${\mathrm{f}}_i={\alpha }_u\left(\frac{⟨\overline{\rho }⟩{\mathrm{U}}_{\mathrm{i}}-⟨\overline{\rho }⟩⟨{\mathrm{u}}_i⟩}{\mathrm{\Delta }t}\right)$$

where $⟨⟩$ denotes volume averaging at a certain time $t$, ${\mathrm{U}}_i$ is the desired spatial averaged velocity, and $\mathrm{\Delta }t$ is the time-scale between when the source term is computed (time $t$) and when it is applied (time $t+\mathrm{\Delta }t$). This is typically chosen to be the simulation time-step. In the case of an ABL simulation with flat terrain, the volume 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-Wind 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 }_{\mathrm{r}\mathrm{e}\mathrm{f}}-⟨\theta ⟩}{\mathrm{\Delta }t}\right)$$

where ${\theta }_{\mathrm{r}\mathrm{e}\mathrm{f}}$ is the desired spatial averaged temperature, $⟨\theta ⟩$ is the spatial averaged temperature, $C_p$ is the heat capacity, ${\alpha }_{\theta }$ is the relaxation factor, and $\mathrm{\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 \overline{\rho }{\tilde{Y}}_k}{\partial t}\mathrm{d}V+\int \overline{\rho }{\tilde{Y}}_k{\tilde{u}}_j{n}_j\mathrm{d}S=-\int {\tau }_{Y_k,j}^{sgs}{n}_j\mathrm{d}S-\int \overline{\rho Y_k{\hat{u}}_{j,k}}{n}_j\mathrm{d}S+\int \overline{{\dot{\omega }}_k}\mathrm{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 }_{Y_{k}j}^{sgs}$ is defined as

$${\tau }_{Y_k,j}^{sgs}\equiv \overline{\rho }(\widetilde{Y_k{u}_j}-\widetilde{Y_k}{\tilde{u}}_j).$$

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

$${\tau }_{Y_k,j}^{sgs}=-\frac{{\mu }_t}{{\mathrm{Sc}}_t}\frac{\partial {\tilde{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 \overline{\rho }{\tilde{Y}}_k}{\partial t}\mathrm{d}V+\int \overline{\rho }{\tilde{Y}}_k{\tilde{u}}_j{n}_j\mathrm{d}S=\int (\frac{\mu }{\mathrm{S}\mathrm{c}}+\frac{{\mu }_t}{{\mathrm{S}\mathrm{c}}_t})\frac{\partial {\tilde{Y}}_k}{\partial x_j}{n}_j\mathrm{d}S+\int {\bar{\dot{\omega }}}_k\mathrm{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

$${\tilde{Y}}_n=1-\sum _{j\ne n}^n{\tilde{Y}}_j.$$

Finally, the reaction rate source term is expected to be added based on an operator split approach whereby 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 $\mathrm{\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}_k^{\ast }-\rho Y_k}{\mathrm{\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 \overline{\rho }{k}^{\mathrm{sgs}}}{\partial t}\mathrm{d}V+\int \overline{\rho }{k}^{\mathrm{sgs}}{\tilde{u}}_j{n}_j\mathrm{d}S=\int \frac{{\mu }_t}{{\sigma }_k}\frac{\partial k^{\mathrm{sgs}}}{\partial x_j}{n}_j\mathrm{d}S+\int (P_k^{\mathrm{sgs}}-D_k^{\mathrm{sgs}})\mathrm{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 {\tilde{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_{ϵ}\frac{{k^{\mathrm{sgs}}}^{\frac{3}{2}}}{\mathrm{\Delta }},$$

where the grid filter length, $\mathrm{\Delta }$, is given in terms of the grid cell volume by

$$\mathrm{\Delta }=V^{\frac{1}{3}}.$$

The subgrid turbulent eddy viscosity is then provided by

$${\mu }_t=C_{{\mu }_{ϵ}}\mathrm{\Delta }{k^{\mathrm{sgs}}}^{\frac{1}{2}},$$

where the values of $C_{ϵ}$ and $C_{{\mu }_{ϵ}}$ 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. RANS Model Suite

Although Nalu-Wind is primarily expected to be a LES simulation tool, RANS modeling is supported through the activation of different two-equation RANS models: the Chien $k-ϵ$ model [Chi82], the Wilcox 1998 $k-\omega$ model [W+98], and the SST model. For the first two models, the reader is referred to the reference papers and the NASA Turbulence Modeling Resource for the Chien and Wilcox models. The SST model is explained in more details below.

2.10.1. Shear Stress Transport (SST) Formulation

It has been observed that standard 1998 $k-\omega$ models display a strong sensitivity to the free stream value of $\omega$ (see Menter, [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-ϵ$ away from the wall. Because of the relationship between $\omega$ and $ϵ$, 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 Menter 2003 model are then

$$\int \frac{\partial \overline{\rho }k}{\partial t}\text{d}V+\int \overline{\rho }k{\tilde{u}}_j{n}_j\text{d}S=\int (\mu +{\hat{\sigma }}_k{\mu }_t)\frac{\partial k}{\partial x_j}{n}_j+\int (P_k^{\omega }-{\beta }^{\ast }\overline{\rho }k\omega )\text{d}V,$$

$$\begin{array}{r}\int \frac{\partial \overline{\rho }\omega }{\partial t}\text{d}V+\int \overline{\rho }\omega {\tilde{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{\overline{\rho }{\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}\text{d}V\\ +\int (\frac{\hat{\gamma }}{{\nu }_t}{P}_k^{\omega }-\hat{\beta }\overline{\rho }{\omega }^2)\text{d}V.\end{array}$$

where the value of ${\beta }^{\ast }$ is 0.09.

The model coefficients, ${\hat{\sigma }}_k$, ${\hat{\sigma }}_{\omega }$, $\hat{\gamma }$ and $\hat{\beta }$ must also be blended, which is represented by

$$\hat{\varphi }=F{\varphi }_1+(1-F){\varphi }_2.$$

where ${\sigma }_{k1}=0.85$, ${\sigma }_{k2}=1.0$, ${\sigma }_{\omega 1}=0.5$, ${\sigma }_{\omega 2}=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 (ar{g}_1^4),$$

where

$$ar{g}_1=min(max(\frac{\sqrt{k}}{{\beta }^{\ast }\omega y},\frac{500\mu }{\overline{\rho }{y}^2\omega }),\frac{4\overline{\rho }{\sigma }_{\omega 2}k}{C{D}_{k\omega }{y}^2}).$$

The final parameter is

$$C{D}_{k\omega }=max(2\overline{\rho }{\sigma }_{\omega 2}\frac{1}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j},{10}^{-10}).$$

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

$${\mu }_t=\frac{a_1\overline{\rho }k}{max(a_1\omega ,S{F}_2)},$$

where $F_2$ is another blending function given by

$$F_2=\tanh (ar{g}_2^2).$$

The final parameter is

$$ar{g}_2=max(\frac{2\sqrt{k}}{{\beta }^{\ast }\omega y},\frac{500\mu }{\overline{\rho }\omega y^2}).$$

The Menter SST Two-Equation Model with Controlled Decay (SST-SUST) is also supported, [SR07]. Two new constants are added that are incorporated into additional source terms for the transport equations:

$$+\int \left({\beta }^{\ast }\overline{\rho }{k}_{amb}{\omega }_{amb}\right)\text{d}V,$$

$$+\int \left(\hat{\beta }\overline{\rho }{\omega }_{amb}^2\right)\text{d}V.$$

where the constants are $k_{amb}$ and ${\omega }_{amb}$. Typically these are set to $k_{amb}={10}^{-6}{U}_{\mathrm{\infty }}^2$ and ${\omega }_{amb}=\frac{5U_{\mathrm{\infty }}}{L}$, where $L$ is a defining length scale for the particular problem, and $U_{\mathrm{\infty }}$ is the freestream velocity. The value chosen for these constants should match the values for $\omega$ and $k$ at the inflow BC.

2.10.2. SST Mixing Length Limiter

When using SST to model the Atmospheric Boundary Layer with the Coriolis effect, a mixing length limiter should be included. The limiter included here is based on the limiter for the k-epsilon model in [Kob13]. An analogous limiter was derived for the SST model. The SST limiter was presented in [AdFM+21] and will be written up in a future publication.

The mixing length limiter replaces the SST model parameter, $\gamma$, in the $\omega$ equation with ${\gamma }^{\ast }$. ${\gamma }^{\ast }$ is a blend of ${\gamma }_1^{\ast }$ and ${\gamma }_2^{\ast }$ using the SST blending function, $F$

$${\gamma }^{\ast }=F{\gamma }_1^{\ast }+(1-F){\gamma }_2^{\ast }.$$

${\gamma }_i^{\ast }$ for $i=1,2$ is calculated from $C_{\epsilon 1,i}^{\ast }$ as

$${\gamma }_i^{\ast }=C_{\epsilon 1,i}^{\ast }-1.$$

$C_{\epsilon 1,i}^{\ast }$ is calculated by applying the mixing length limiter to $C_{\epsilon 1,i}$ as

$$C_{\epsilon 1,i}^{\ast }=C_{\epsilon 1,i}+(C_{\epsilon 2,i}-C_{\epsilon 1,i})\frac{l_t}{l_e}.$$

$C_{\epsilon 1,i}$ is calculated from the SST model constant $gamm{a}_1$ as

$$C_{\epsilon 1,i}={\gamma }_i+1.$$

$C_{\epsilon 2,i}$ is calculated from the SST model constants ${\beta }_i$ and ${\beta }^{\ast }$ as

$$C_{\epsilon 1,i}=\frac{{\beta }_i}{{\beta }^{\ast }}+1.$$

The maximum mixing length, $l_e$ is calculated as

$$l_e=.00027G/f_c,$$

where $G$ is the geostrophic (freestream) velocity and $f_c$ is the Coriolis force. The Coriolis force is calculated as

$$f_c=2\mathrm{\Omega }sin\lambda ,$$

where $\mathrm{\Omega }$ is the earth’s angular velocity and $\lambda$ is the latitude.

2.11. Laminar-Turbulent Transition Model Formulation

To account for the effects of the laminar-turbulent boundary layer transition, Menter’s one-equation $\gamma$ transition model [MSLA15] is coupled with the SST model. The model consists of single transport equation for intermittency

$$\frac{D(\rho \gamma )}{Dt}=P_{\gamma }-D_{\gamma }+\frac{\partial }{\partial x_j}[(\mu +\frac{{\mu }_t}{{\sigma }_{\gamma }})\frac{\partial \gamma }{\partial x_j}]$$

The production term, $P_{\gamma }$, and destruction term, $D_{\gamma }$, are as:

$$P_{\gamma }=F_{length}\rho S\gamma (1-\gamma )F_{onset},D_{\gamma }=C_{a2}\rho \mathrm{\Omega }\gamma F_{turb}(C_{e2}\gamma -1)$$

The model constants are:

$$F_{length}=100,c_{e2}=50,c_{a2}=0.06,{\sigma }_{\gamma }=1.0$$

The transition onset criteria of the model are defined as:

$$F_{onset1}=\frac{R{e}_v}{2.2R{e}_{\theta c}},F_{onset2}=(F_{onset1},2.0)$$

$$F_{onset3}=max(1-{\left(\frac{R_T}{3.5}\right)}^3,0),F_{onset}=max(F_{onset2}-F_{onset3},0)$$

$$F_{turb}=e^{-{\left(\frac{R_T}{2}\right)}^4},R_T=\frac{\rho k}{\mu \omega },R{e}_v=\frac{\rho d_w^{2}S}{\mu }$$

The transition onset occurs once the scaled vorticity Reynolds number, $R{e}_v/2.2$, exceeds the critical momentum thickness Reynolds number, $R{e}_{\theta c}$, from the empirical correlations.

The output intermittency from the transition model is applied to both the production and destruction terms of the turbulent kinetic energy transport equation.

2.12. Detached 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 described 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 }^{\ast }\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.

2.13. Active Model Split (AMS) Formulation

The AMS approach is a recently developed hybrid RANS/LES framework by Haering, Oliver and Moser [HOM20]. In this approach a triple decomposition is used, breaking the instantaneous velocity field into an average component $<u_i>$, a resolved fluctuation $u_i^{>}$ and an unresolved fluctuation $u_i^{<}$. The subgrid stress is split into two terms, ${\tau }_{ij}={\tau }_{ij}^{SGRS}+{\tau }_{ij}^{SGET}$, with ${\tau }_{ij}^{SGRS}$ modeling the mean subgrid stress, taking on the form of a standard RANS subgrid stress model and ${\tau }_{ij}^{SGET}$ representing the energy transfer from the resolved to the modeled scales. In addition, a forcing stress $F_i$ is added to the momentum equations to induce the transfer of energy from the modeled to the resolved scales. Thus the AMS approach solves the following momentum equation

(2.26)$$\begin{array}{rl}& \int \frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial t}\mathrm{d}V+\int \overline{\rho }{\tilde{u}}_i{\tilde{u}}_j{n}_j\mathrm{d}S+\int (\overline{P}+\frac{2}{3}\overline{\rho }k)n_i\mathrm{d}S=\\ & \int 2\mu ({\tilde{S}}_{ij}-\frac{1}{3}{\tilde{S}}_{kk}{\delta }_{ij})n_j\mathrm{d}S+\int (\overline{\rho }-{\rho }_{\circ })g_i\mathrm{d}V+\\ & \int 2{\mu }_t(<S_{ij}>-\frac{1}{3}<S_{kk}>{\delta }_{ij})n_j\mathrm{d}S+\int ({\mu }_{ik}^t{\tilde{u}}_j+{\mu }_{jk}^t{\tilde{u}}_i)n_k\mathrm{d}S+\int f_i\mathrm{d}V.\end{array}$$

2.13.1. Split subgrid model stress

In a typical Hybrid RANS/LES approach, the observation that the LES and RANS equations take on the same mathematical form is leveraged, relying simply on a modified turbulent viscosity coefficient that takes into account the ability to resolve some turbulent content. Due to deficiencies in this approach as discussed in Haering et al. [HOM20], an alternative approach where the modeled term is split into two contributions, one representing the impact on the mean flow and one the impact on the resolved fluctuations, from the unresolved content, is used in the Active Model Split (AMS) formulation.

Starting by substitution of a triple decomposition of the flow variables, $\varphi =⟨\varphi ⟩+{\varphi }^{>}+{\varphi }^{<}$, with $⟨\cdot ⟩$ representing a mean quantity, ${\varphi }^{>}$ a resolved fluctuation and ${\varphi }^{<}$ an unresolved fluctuation and dropping all terms that have an unresolved fluctuation in them (since by definition, these terms cannot be resolved and thus must be modeled) we get the following momentum equation:

$$\frac{\partial \overline{u_i}}{\partial t}+\frac{\partial \overline{u_i}\overline{u_j}}{\partial x_j}=-\frac{1}{\rho }\left(\frac{\partial \bar{P}}{\partial x_i}\right)+\nu \frac{{\partial }^2\overline{u_i}}{\partial x_j\partial x_j}+\frac{\partial {\tau }_{ij}^M}{\partial x_j}+F_i$$

Note that here, $\bar{\varphi }=⟨\varphi ⟩+{\varphi }^{>}$ represents an instantaneous resolved quantity and $F_i$ is a forcing term discussed in Sec. AMS forcing.

The model term in AMS, ${\tau }_{ij}^M$ is split into two pieces, the first representing the impact of the unresolved scales on the mean flow, referred to as ${\tau }_{ij}^{SGRS}$, since this mimics the purpose of RANS models and seeks to model the subgrid Reynolds Stress (SGRS). The second term represents the impact of the unresolved scales on the resolved fluctuations which acts to capture energy transfer from the resolved fluctuations to the unresolved fluctuations, which as Haering et al. points out, is really the primary function of typical LES SGS models. As this term models subgrid energy transfer (SGET), it is referred to as ${\tau }_{ij}^{SGET}$.

${\tau }_{ij}^{SGRS}$ is modeled using a typical RANS model, but since in the hybrid context, some turbulence is resolved, the magnitude of the stress tensor is reduced through a derived scaling with $\alpha ={\beta }^{1.7}$, $\beta =1-k_{res}/k_{tot}$, where $k_{tot}$ is the total kinetic energy, taken from the RANS model and $k_{res}=0.5⟨u_i^{>}{u}_i^{>}⟩$, a measure of the average resolved turbulent kinetic energy.

${\tau }_{ij}^{SGET}$ is modeled using the M43 SGS model discussed in Haering et al. [HOM19]. This uses an anisotropic representation, ${\tau }_{ij}={\nu }_{jk}{\partial }_k{u}_i+{\nu }_{ik}{\partial }_k{u}_j$, of the stress tensor and a tensorial eddy viscosity, ${\nu }_{ij}=C(\mathcal{M})⟨ϵ{⟩}^{1/3}{\mathcal{M}}_{ij}^{4/3}$, with $C$, a coefficient determined as a function of the eigenvalues of $\mathcal{M}$, a metric tensor measure of the grid and $⟨ϵ⟩$ inferred from the RANS model.

So this produces the final form for the AMS model term,

$$\begin{array}{r}\begin{array}{rl}{\tau }_{ij}^M& ={\tau }_{ij}^{SGRS}+{\tau }_{ij}^{SGET}\\ & =2\alpha (2-\alpha ){\nu }_t^{RANS}⟨S_{ij}⟩-\frac{2}{3}\beta k_{tot}{\delta }_{ij}+C(\mathcal{M})⟨ϵ{⟩}^{1/3}({\mathcal{M}}_{jk}^{4/3}\frac{\partial u_i^{>}}{\partial x_k}+{\mathcal{M}}_{ik}^{4/3}\frac{\partial u_j^{>}}{\partial x_k}).\end{array}\end{array}$$

The AMS model terms are implemented for the edge based scheme in MomentumSSTAMSDiffEdgeKernel. The isotropic component, $2\beta k_{tot}{\delta }_{ij}/3$ is brought into the pressure term.

2.13.2. Averaging functions

The averaging algorithms are invoked as part of the AMSAlgDriver and are called from the pre_iter_work function so that they are executed at the beginning of each Picard iteration. The AMSAlgDriver is invoked last, so to ensure that this is also done initially, so that the initial step has correct average quantities, the averaging functions are also called in the initial_work function.

The main averaging algorithm is SSTAMSAveragesAlg. The averaging function is solving a simple causal average equation for the intermediate (or final) quantity:

$$\frac{\partial ⟨{\varphi }^{\ast }⟩}{\partial t}=\frac{1}{T_{RANS}^{\ast }}({\varphi }^{\ast }-⟨{\varphi }^n⟩)$$

Here $⟨\cdot ⟩$ refers to a mean (time-averaged) quantity and $T_{RANS}$ is the timescale of the turbulence determined by the underlying RANS scalars ($1/({\beta }^{\ast }\omega )$ in SST). $({)}^n$ refers to a previous timestep quantity and $({)}^{\ast }$ refers to an intermediate quantity. Note that currently the time scale is stored in a nodal field.

We can discretize the causal average equation explicitly to produce the implemented form:

$$\begin{array}{r}\begin{array}{rl}⟨{\varphi }^{\ast }⟩& =⟨{\varphi }^n⟩+\frac{\mathrm{\Delta }t}{T_{RANS}^{\ast }}({\varphi }^{\ast }-⟨{\varphi }^n⟩)\\ ⟨{\varphi }^{\ast }⟩& =\frac{\mathrm{\Delta }t}{T_{RANS}^{\ast }}{\varphi }^{\ast }+(1-\frac{\mathrm{\Delta }t}{T_{RANS}^{\ast }})⟨{\varphi }^n⟩\end{array}\end{array}$$

The averaging is carried out for velocities ($u_i$), velocity gradients ($\partial u_i/\partial x_j$), pressure ($P$), density ($\rho$), resolved turbulent kinetic energy ($k_{res}=0.5u_i^{>}{u}_i^{>}$) and the kinetic energy production $({\mathcal{P}}_k=⟨S_{ij}⟩({\tau }_{ij}^{SGRS}-u_i^{>}{u}_j^{>}))$. Note that ${}^{>}$ is used to denote a resolved fluctuation, i.e. ${\varphi }^{>}=\varphi -⟨\varphi ⟩$.

2.13.3. Dynamic Hybrid Diagnostics

Typically in a hybrid model, it is necessary to have diagnostics that assess the ability of the grid to resolve turbulent content and to aid in its introduction. In AMS, there are two main diagnostic quantities, $\alpha ={\beta }^{1.7}=(1-k_{res}/k_{tot}{)}^{1.7}$ and a resolution adequacy parameter, $r_k$, which is used to evaluate where in the flow the grid and the amount of resolved turbulent content is inconsistent.

$\beta$ is a straight-forward calculation. Limiters are applied to $\beta$ to realize the RANS and DNS limits. In the RANS limit, $k_{res}=0$ and thus $\beta =1$, so $\beta$ is limited from evaluation above 1. In the DNS limit, ideally, the ratio of the approximate Kolmogorov velocity scale to total TKE would be used as a lower bound,

$$\beta =max(1-\frac{k_{res}}{k_{tot}},\frac{(\nu ϵ{)}^{1/2}}{k_{tot}}),$$

but that has shown some performance issues near the wall when using SST with AMS. Currently an adhoc lower bound of $\beta =0.01$ is used instead. The resolution adequacy parameter is based on the ratio between the anisotropic grid metric tensor, $\mathcal{M}={\mathcal{J}}^{T}J$, where $\mathcal{J}$ is the mapping from a unit cube to the cell, and the length scale associated with the model production. It takes the form,

$$r_k={\left(\frac{3}{2⟨{\bar{v}}^2⟩}\right)}^{3/2}\underset{\lambda }{max}({\mathcal{P}}_{ik}^{SGS}{\mathcal{M}}_{kj}).$$

For the RANS models used in Nalu-Wind, $⟨{\bar{v}}^2⟩\approx 5{\nu }_{RANS}/T_{RANS}$. ${\mathcal{P}}_{ij}^{SGS}=\frac{1}{2}({\tau }_{ik}\partial {\bar{u}}_k/\partial x_j+{\tau }_{jk}\partial {\bar{u}}_k/\partial x_i)$ is the full subgrid production tensor, with ${\tau }_{ij}={\tau }_{ij}^{SGRS}+{\tau }_{ij}^{SGET}+\frac{2}{3}\alpha k_{tot}{\delta }_{ij}$.

2.13.4. Forcing Term

When the grid is capable of resolving some turbulent content, the model will want to reduce the modeled stress and allow for resolved turbulence to contribute the remaining piece of the total stress. As discussed in Haering et al. [HOM20] and the observation of “modeled-stress depletion” in other hybrid approaches, such as DES, a mechanism for inducing resolved turbulent fluctuations at proper energy levels and timescales to match your reduction of the modeled stress is needed. AMS resolves this issue through the use of an active forcing term, designed to introduce turbulent fluctuations into regions of the grid where turbulent content can be supported. The implications of the specific form and method of introduction for this forcing term is still an area of ongoing research, but for now, empirical testing has shown great success with a simple turbulent structure based off of Taylor-Green vortices.

The forcing term $F_i$ is determined by first specifying an auxiliary field based off of a Taylor-Green vortex:

$$\begin{array}{r}\begin{array}{rl}{h}_1& =\frac{1}{3}\cos (a_x{x}_1^{\prime })\sin (a_y{x}_2^{\prime })\sin (a_z{x}_3^{\prime }),\\ h_2& =-\sin (a_x{x}_1^{\prime })\cos (a_y{x}_2^{\prime })\sin (a_z{x}_3^{\prime }),\\ h_3& =\frac{2}{3}\sin (a_x{x}_1^{\prime })\sin (a_y{x}_2^{\prime })\cos (a_z{x}_3^{\prime }),\end{array}\end{array}$$

with ${\mathbf{x}}^{\prime }=\mathbf{x}+⟨\mathbf{u}⟩t$ and $a_i=\pi /{\mathbb{P}}_i$. $\mathbb{P}$ is determined as follows,

$$\begin{array}{r}\begin{array}{rl}l& =\frac{4-(1-max(\beta ,0.8))}{0.4}\frac{(\beta k{)}^{3/2}}{ϵ}\\ l& =min(max(l,70\frac{{\nu }^{3/4}}{{ϵ}^{1/4}}),d)\\ \\ l_i^{\prime }& =min(l,L_{p_i})\\ f_i& =\mathrm{nint}\left(\frac{L_{p_i}}{l_i^{\prime }}\right)\\ {\mathbb{P}}_i& =\frac{L_{p_i}}{f_i},\end{array}\end{array}$$

where $L_{p_i}$ is related to the periodic domain size and is user specified. With the initial TG vortex field, $h_i$, determined, we now determine a scaling factor ($\eta$) for the forcing.

$$\begin{array}{r}\begin{array}{rl}{T}_{\beta }& =max(\beta k/ϵ,6\sqrt{\nu /ϵ})\\ F_{tar}& =8\sqrt{\alpha {\bar{v}}^2}/T_{\beta }\\ \\ {\mathcal{P}}_r& =\mathrm{\Delta }t{F}_{tar}\left(h_i{u}_i^{>}\right)\\ \\ {\beta }_K& =min(\sqrt{\nu ϵ/k},1)\\ \\ \hat{\beta }& =\{\begin{array}{rl}\frac{1-\beta }{1-{\beta }_K}& {\beta }_K<1\\ 10000& \mathrm{else}\end{array}\\ \\ C_f& =-1\tanh (1-\frac{1}{\sqrt{min(⟨r_k⟩},1)})\\ \\ C_f& =C_f(1.0-min(\tanh (10(\hat{\beta }-1))+1,1))\\ \\ \eta & =\{\begin{array}{rl}{F}_{tar}{C}_f& ⟨r_k⟩<1,{\mathcal{P}}_r\ge 0\\ 0& \mathrm{else}\end{array}\end{array}\end{array}$$

Now the final forcing field, $F_i=\eta h_i$. Since this is being added as a source term to the momentum solve, we are not projecting onto a divergence free field and are instead allowing that to pass into the continuity solve, where the intermediate velocity field with the forcing will then be projected onto a divergence free field. This is implemented in the node kernels as MomentumSSTAMSForcingNodeKernel.

2.13.5. AMS with SST Mixing Length Limiter

When using AMS with SST as the mean (RANS) contribution to model the Atmospheric Boundary Layer with the Coriolis effect, SST should include a mixing length limiter. The mixing length limiter is described in SST Mixing Length Limiter. For consistency, when using the limiter the RANS time scale, $T_{RANS}^{\ast }$, should depend on the mixing length rather than $\omega$ to account for the effect of the limiter. The time scale becomes

$$T_{RANS}^{\ast }=\frac{l_t}{\sqrt{k}}.$$

2.14. 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.27)$$\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.28)$${\sigma }_{ij}=\mu (\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})+\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.29)$$\mu =\frac{E}{2(1+\nu )},$$

and

(2.30)$$\lambda =\frac{E\nu }{(1+\nu )(1-2\nu )}.$$

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

(2.31)$$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.32)$$v_i=\frac{\partial u_i}{\partial t}.$$

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

(2.33)$$v_i=\frac{{\gamma }_1{u}_i^{n+1}+{\gamma }_2{u}_i^n+{\gamma }_3{u}_i^{n-1}}{\mathrm{\Delta }t}$$

2.15. 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 $\varphi$ can be written as:

(2.34)$$\int \frac{\rho \varphi }{\partial t}dV+\int \rho \varphi (u_j-v_j)n_{j}dS+\int \rho \varphi \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 transferred 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.16. 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.35)$$s_i\frac{\partial }{\partial x_i}I\left(s\right)+({\mu }_a+{\mu }_s)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 coefficient, $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 linear, 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.36)$$\frac{\partial q_i^r}{\partial x_i}={\mu }_a[4\sigma T^4-G],$$

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_i^r={\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 2.1.

2.1 Ordinate Direction Definition, $\mathbf{s}=\sin {\theta }_{zn}\sin {\theta }_{az}\mathbf{i}+\cos {\theta }_{zn}\mathbf{j}+\sin {\theta }_{zn}\cos {\theta }_{az}\mathbf{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.37)$$I\left(s\right)=\frac{1}{\pi }[\tau \sigma T_{\mathrm{\infty }}^4+ϵ\sigma T_w^4+(1-ϵ-\tau )K],$$

where $ϵ$ is the total normal emissivity of the surface, $\tau$ is the transmissivity of the surface, $T_w$ is the temperature of the boundary, $T_{\mathrm{\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 +ϵ=1.$$

where it is implied that $\alpha =ϵ$.

2.17. Wall Distance Computation

RANS and DES simulations using $k-\omega$ SST or SST-DES equations requires the specification of a wall distance for computing various turbulence parameters. For static mesh simulations this field can be generated using a pre-processing step and provided as an input in the mesh database. However, for moving mesh simulations, e.g., blade resolved wind turbine simulations, this field must be computed throughout the course of the simulation. Nalu-Wind implements a Poisson equation ([G03]) to determine the wall distance $d$ using the gradients of a field $\varphi$.

$$\begin{array}{rl}{\mathrm{\nabla }}^2\varphi & =1\\ d& =\pm \sqrt{\sum _{j=1,3}{\left(\frac{\partial \varphi }{\partial x_j}\right)}^2}+\sqrt{\sum _{j=1,3}{\left(\frac{\partial \varphi }{\partial x_j}\right)}^2+2\varphi }\end{array}$$