Theory

This page introduces the governing equations and the numerical discretization at a high level. nekRS includes models for incompressible flow, a partially compressible low-Mach formulation, the Stokes equations, and the \(k\)-\(\tau\) RANS equations.

Incompressible Flow Model

The governing equations for conservation of mass, momentum, and energy for an incompressible fluid are

\[\nabla\cdot\mathbf u=0\]

\[\rho\left(\frac{\partial\mathbf u}{\partial t}+\mathbf u\cdot\nabla\mathbf u\right)=-\nabla P+\nabla\cdot\tau+\rho\ \mathbf f\]

\[\rho C_p\left(\frac{\partial T}{\partial t}+\mathbf u\cdot\nabla T\right)=\nabla\cdot\left(k\nabla T\right)+\dot{q}\]

In these equations, \(t\) is time, \(\mathbf u\) is the velocity vector, \(\rho\) is the density, \(P\) is the pressure, \(\tau\) is the viscous stress tensor, \(\mathbf f\) is a force vector, \(C_p\) is the isobaric specific heat capacity, \(T\) is the temperature, \(k\) is the thermal conductivity, and \(\dot{q}\) is a volumetric heat source. If the viscosity is constant, the viscous stress tensor can be contracted to give

\[\nabla\cdot\tau=\nabla^2\mathbf u\]

This is referred to as the “no-stress” formulation. In the general case for non-constant viscosity, the viscous stress tensor is given by the Navier-Stokes closure as

\[\nabla\cdot\tau=\nabla\cdot\left\lbrack\mu\left(\nabla\mathbf u+\nabla\mathbf u^T\right)\right\rbrack\]

Non-Dimensional Formulation

It is often advantageous to solve these equations in non-dimensional form. Here, we introduce the non-dimensional form in as general a manner as possible, assuming the use of variable material properties for density, viscosity, specific heat capacity, and thermal conductivity that are functions of temperature, \(\rho=\rho(T)\), \(\mu=\mu(T)\), \(C_p=C_p(T)\), and \(k=k(T)\). For simplicity, the functional notation is omitted throughout.

Introduce the non-dimensional variables \(\mathbf x^\dagger=\frac{\mathbf x}{L}\), \(\mathbf u^\dagger=\frac{\mathbf u}{U}\), \(t^\dagger=\frac{tU}{L}\), and \(\mathbf f^\dagger=\frac{\mathbf f L}{U^2}\). For the material properties, we non-dimensionalize based on some reference temperature \(T_0\), such that \(\rho^\dagger=\frac{\rho}{\rho_0}\), \(\mu^\dagger=\frac{\mu}{\mu_0}\), \(k^\dagger=\frac{k}{k_0}\), and \(C_p^\dagger=\frac{C_p}{C_{p,0}}\). Here, a subscript of \(0\) is shorthand notation that indicates that the property is evaluated at \(T_0\), such that \(k_0\equiv k(T_0)\). Finally, for convection-dominated flows, the pressure is non-dimensionalized in terms of the dynamic pressure as \(P^\dagger=\frac{P}{\rho_0 U^2}\).

Inserting these non-dimensional variables into the conservation of mass and momentum equations gives

\[\frac{\partial u_i^\dagger}{\partial x_i^\dagger}=0\]

\[\rho^\dagger\left(\frac{\partial u_i^\dagger}{\partial t^\dagger}+u_j^\dagger\frac{\partial u_i^\dagger}{\partial x_j^\dagger}\right)=-\frac{\partial P^\dagger}{\partial x_i^\dagger}+\frac{1}{Re}\frac{\partial \tau_{ij}^\dagger}{\partial x_j^\dagger}+\rho^\dagger f_i^\dagger\]

In these equations, the \(\nabla\) are expanded to explicitly show that all derivatives are taken with respect to the nondimensional space variable \(\mathbf x^\dagger\). \(Re\) is the Reynolds number

\[Re\equiv\frac{\rho_0 UL}{\mu_0}\]

To non-dimensionalize the energy conservation equation, use the previous non-dimensional variables in addition to a non-dimensional temperature, \(T^\dagger=\frac{T-T_0}{\Delta T}\), where \(\Delta T\) is a reference temperature rise relative to a baseline temperature \(T_0\). The heat source is non-dimensionalized as \(\dot{q}^\dagger=\frac{\dot{q}}{\rho_0 C_{p,0} U\Delta T/L}\), which arises naturally from the simple formulation of bulk energy conservation of \(Q=\dot{m}C_p\Delta T\), where \(Q\) is a heat source (units of Watts) and \(\dot{m}\) is a mass flowrate.

Inserting these non-dimensional variables into the energy conservation equation gives

\[\rho^\dagger C_p^\dagger\left(\frac{\partial T^\dagger}{\partial t^\dagger}+u_i^\dagger\frac{\partial T^\dagger}{\partial x_i^\dagger}\right)=\frac{1}{Pe}\frac{\partial}{\partial x_i^\dagger}\left(k^\dagger\frac{\partial T^\dagger}{\partial x_i^\dagger}\right)+\dot{q}^\dagger\]

where \(Pe\) is the Peclet number,

\[Pe\equiv\frac{LU}{\alpha}\]

and \(\alpha\) is the thermal diffusivity,

\[\alpha\equiv\frac{k_0}{\rho_0 C_{p,0}}\]

Low-Mach Partially-Compressible Model

Stokes Equations

RANS Models

The RANS equations are derived from the conservation of mass, momentum, and energy equations by expressing each term in the equation as the sum of a mean and a fluctuation. Because nekRS is based on the incompressible flow model, all such averages (even for the energy equation) are based on the notion of Reynolds averaging, where each field \(f\) is expressed as the sum of a time mean \(\overline{f}\) and a time fluctuation, \(f'\),

\[f(\mathbf x, t)=\overline{f}(\mathbf x)+f'(\mathbf x,t)\]

where the time averaged is defined as

\[\overline{f}=\lim_{S\rightarrow\infty}\frac{1}{S}\int_{t}^{t+S}f(\mathbf x,t)dt\]

For compressible flows in which energy conservation affects density, the RANS equations are instead derived with Favre averaging, where each field \(f_i\) is expressed as the sum of a density-weighted time average \(\tilde{f}_i\) and a fluctuation \(f_i^{''}\). It is therefore an important distinction here that we only consider Reynolds averaging, which leads to a simpler formulation of the RANS energy conservation equation that the compressible flow case.

Inserting the above “Reynolds decomposition” for \(\mathbf u\), \(P\), and \(T\) into the governing equations and averaging in time then gives the RANS equations. For the incompressible flow equations in Incompressible Flow Model, the RANS mass, momentum, and energy equations are

\[\frac{\partial\overline{u_i}}{\partial x_i} = 0\]

\[\rho\left(\frac{\partial\overline{u_i}}{\partial t}+\overline{u_j}\frac{\partial\overline{u_i}}{\partial x_j}+\frac{\partial}{\partial x_j}\overline{u_i'u_j'}\right)=-\frac{\partial \overline{P}}{\partial x_i}+\frac{\partial}{\partial x_j}\left(2\mu \overline{S_{ij}}\right)+\rho\overline{\mathbf f}\]

\[\rho C_p\left(\frac{\partial\overline{T}}{\partial t}+\overline{u_i}\frac{\partial\overline{T}}{\partial x_i}+\frac{\partial\overline{u_i'T'}}{\partial x_i}\right)=\frac{\partial}{\partial x_i}\left(k\frac{\partial\overline{T}}{\partial x_i}\right)+\overline{\dot{q}}\]

where \(\overline{S_{ij}}\) is the mean strain rate tensor,

\[\overline{S_{ij}}=\frac{1}{2}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial\overline{u_j}}{\partial x_i}\right)\]

The mass, momentum, and energy conservation equations have the same form as the instantaneous flow equations in Incompressible Flow Model except for the addition of another stress tensor to the momentum equation - \(\rho \overline{u_i'u_j'}\), and the addition of another heat flux vector to the energy equation - \(\rho C_p\overline{u_i'T'}\). The stress term in the momentum equation is referred to as the Reynolds stress tensor; \(\rho\ \partial(\overline{u_i'u_j'})/\partial x_j\) represents the time-averaged rate of momentum transfer due to turbulence. The heat flux term in the energy equation is referred to as the turbulent heat flux; \(\rho C_p\partial\overline{u_i'T'}/\partial x_i\) represents the time-averaged rate of energy addition due to turbulence. The objective of RANS models is to provide closures for the Reynolds stress tensor and turbulent heat flux vector in terms of the mean properties such that the time-averaged equations can be solved for the mean flow.

Boussinesq Approximation

The RANS models in nekRS are based on the Boussinesq eddy viscosity approximation, which assumes that the momentum flux that induces the Reynolds stresses shares the same functional form as the momentum flux that induces the molecular stresses. In other words, the Navier-Stokes closure that was used to relate the deviatoric stress tensor \(\tau_{ij}\) to the strain rate tensor,

\[\tau_{ij}=\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\underbrace{\frac{2}{3}\mu\frac{\partial u_i}{\partial x_i}\delta_{ij}}_\text{$=\ 0$ if incompressible}\]

is assumed applicable to the Reynolds stress tensor, but with instantaneous velocities replaced by mean velocities and the molecular viscosity replaced by the turbulent eddy viscosity \(\mu_T\),

\[\rho\overline{u_i'u_j'}=\mu_T\left(\frac{\partial\overline{u_i}}{\partial x_j}+\frac{\partial\overline{u_j}}{\partial x_i}\right)-\underbrace{\frac{2}{3}\mu\frac{\partial \overline{u_i}}{\partial x_i}\delta_{ij}}_\text{$=\ 0$ if incompressible}-\frac{2}{3}\rho k\delta_{ij}\]

Here, \(k\) is the turbulent kinetic energy,

\[k\equiv\frac{1}{2}\left(\overline{u_1'u_1'}+\overline{u_2'u_2'}+\overline{u_3'u_3'}\right)\]

The final term in the Boussinesq approximation for the Reynolds stress tensor simply ensures that the trace of the Reynolds stress tensor equals \(2k\), because otherwise, for incompressible flows, the trace of the Reynolds stress tensor would be zero. Inserting the Boussinesq eddy viscosity model for the Reynolds stress tensor into the incompressible flow mean momentum equation then gives

\[\rho\left(\frac{\partial\overline{u_i}}{\partial t}+\overline{u_j}\frac{\partial\overline{u_i}}{\partial x_j}\right)=-\frac{\partial \overline{P}}{\partial x_i}+\frac{\partial}{\partial x_j}\left\lbrack 2\left(\mu+\mu_T\right) \overline{S_{ij}}-\frac{2}{3}\rho k\delta_{ij}\right\rbrack+\rho\overline{\mathbf f}\]

In nekRS, as well as many other RANS codes, it is commonplace to combine the gradient of \(P+\frac{2}{3}\rho k\) terms together into a single reduced pressure,

\[P_r\equiv\overline{P}+\frac{2}{3}\rho k\]

such that the term proportional to \(\rho k\) can be bundled into a single pressure gradient kernel,

\[-\frac{\partial\overline{P}}{\partial x_i}-\frac{\partial}{\partial x_j}\left(-\frac{2}{3}\rho k\delta_{ij}\right)\rightarrow-\frac{\partial P_r}{\partial x_i}\]

Warning

The pressure solution, available on the nrs->P object and written to output under the name “pressure,” represents this reduced pressure. To obtain \(\overline{P}\), you should subtract \(\frac{2}{3}\rho k\) from nrs->P.

Turbulent Prandtl Number

Closure for the turbulent heat flux is typically motivated from considerations of the analogy between momentum and energy transfer; while the Boussinesq approximation was used to introduce a relationship between the Reynolds stress tensor \(\rho\overline{u_i'u_j'}\) in terms of the mean strain rate, the turbulent heat flux is assumed proportional to the mean temperature gradient via a gradient diffusion approximation,

\[\rho C_p\overline{u_i'T'}=k_T\frac{\partial \overline{T}}{\partial x_i}\]

where \(k_T\) is the turbulent conductivity. \(k_T\) is related to \(\mu_T\), the turbulent momentum diffusivity, by the turbulent Prandtl number \(Pr_T\),

\[Pr_T\equiv\frac{\nu_T}{\alpha_T}\]

where \(\nu_T\equiv\mu_T/rho\) and \(\alpha_T\) is the turbulent thermal diffusivity,

\[\alpha_T\equiv\frac{k_T}{\rho C_p}\]

Inserting this gradient diffusion approximation into the incompressible flow mean energy equation then gives

\[\rho C_p\left(\frac{\partial\overline{T}}{\partial t}+\overline{u_i}\frac{\partial\overline{T}}{\partial x_i}\right)=\frac{\partial}{\partial x_i}\left\lbrack\left(k+\frac{\mu_T}{Pr_T}C_p\right)\frac{\partial\overline{T}}{\partial x_i}\right\rbrack+\overline{\dot{q}}\]

The \(k\)-\(\tau\) Model

nekRS uses the \(k\)-\(\tau\) turbulence model to close the mean flow equations [Kalitzin]. Because the \(k\)-\(\epsilon\), \(k\)-\(\omega\), and \(k\)-\(\omega\) SST models tend to dominate the RANS space, extra discussion is devoted here to motivating the use of this particular model. Because \(Pr_T\) is typically taken as a constant, often 0.90 [Wilcox], the objective of incompressible flow RANS models is to compute the eddy viscosity and \(k\) needed to close the mean momentum and energy equations.

Note

Take care not to confuse the inverse of the specific dissipation rate, \(\tau\), with the deviatoric molecular stress tensor, which is also represented here as \(\tau\) due to convention.

The \(k\)-\(\tau\) model is a modification of the standard \(k\)-\(\omega\) turbulence model that bases the second transport equation on the inverse of the specific dissipation rate \(\omega\),

\[\tau\equiv\frac{1}{\omega}\]

rather than the on \(\omega\). The \(k\)-\(\tau\) model attempts to retain two important features of the \(k\)-\(\omega\) model -

1. Good predictions for flows with adverse pressure gradients and separation, and

2. Reasonable prediction of boundary layers and near-wall behavior without wall functions or special low-\(Re_t\) treatments.

These two aspects contribute to better predictions of complex flows with reduced numerical complexity associated with wall functions or damping functions that can cause stiff behavior [Kok] and inaccurate flow predictions. By introducing the definition of \(\tau\equiv 1/\omega\), the \(k\)-\(\tau\) model attemps to improve upon the \(k\)-\(\omega\) model in two main ways -

1. Simplify wall boundary conditions for the second transport equation, and

2. Bound the source terms in the second transport equation in near-wall regions.

As \(y\rightarrow 0\), \(\omega\rightarrow y^{-2}\), while \(k\rightarrow 0\) [Kok]. Therefore, while \(\omega\) is infinite at walls, \(\tau\) is zero. Traditionally, this singular behavior in \(\omega\) was treated by applying “rough wall” boundary conditions to \(\omega\) with the wall roughness set to a “small enough” value to simulate a hydraulically smooth wall [Kok]. However, this ad hoc approach retains a strong dependence on the near-wall mesh resolution, often requiring prohibitively fine elements to accurately predict boundary layer properties [Kalitzin]. And, such an approach retains near-singular behavior in the first and second derivatives of \(\omega\). Applying a zero boundary condition to \(\tau\) on solid walls is comparatively trivial.

With regards to the second point, the \(\omega\) transport equation contains a source term propotional to \(\omega^2\); because \(\omega\rightarrow y^{-2}\) as \(y\rightarrow 0\), this source term displays singular behavior as \(y\rightarrow 0\). Singular behavior of the source terms can result in large numerical errors and stiffness that negatively affects the convergence of the computational solution. Conversely, all source terms in the \(\tau\) equation are bounded near walls [Kok].

With this motivation, the \(k\) and \(\tau\) equations are described next. A slightly lengthier description is provided for each in order to give greater context to the genesis of this model.

The \(k\) Equation

The \(k\) equation is a model version of the true \(k\) equation. The true \(k\) equation is derived by taking the trace of the Reynolds stress equation, a process that is itself motivated by recognition that the trace of the Reynolds stress tensor is equal to \(2k\),

\[\overline{u_i'u_i'}=2k\]

The true \(k\) equation contains terms that depend on the mean flow velocity, the turbulent kinetic energy, and the dissipation, in addition to more exotic terms such as \(\overline{u_i'u_i'u_j'}\) and \(\overline{P'u_j'}\). These additional fluctuating terms do not bring the true \(k\) equation any closer to a tractable solution, so Prandtl introduced a \(\partial k/\partial x_j\) gradient diffusion approximation for the turbulent transport and pressure diffusion terms (\(\frac{1}{2}\rho\overline{u_i'u_i'u_j'}+\overline{P'u_j'}\)) with a diffusion coefficient of \(\mu_T/\sigma_k\), where \(\sigma_k\) is a constant [Wilcox]. With this gradient diffusion model, the true \(k\) equation is simplified to a tractable model \(k\) equation [Launder],

\[\frac{\partial(\rho k )}{\partial t}+\nabla\cdot\left(\rho k\overline{\mathbf u}\right)=\nabla\cdot\left\lbrack\left(\mu+\frac{\mu_T}{\sigma_k}\right)\nabla k\right\rbrack+\mathscr{P}-\rho\epsilon\]

where \(\mathscr{P}\) is the production of turbulent kinetic energy by velocity shear,

\[\mathscr{P}\equiv\rho\overline{u_i'u_j'}\frac{\partial\overline{u_i}}{\partial x_j}\]

and \(\epsilon\) is the dissipation per unit mass,

\[\epsilon\equiv\nu\overline{\frac{\partial u_i'}{\partial x_j}\frac{\partial u_i'}{\partial x_j}}\]

The production term represents the rate at which energy is transferred from the mean flow to the turbulent flow, while the dissipation term represents the rate at which turbulent kinetic energy is converted to heat. Note that the only difference between this model \(k\) equation and the true \(k\) equation is the introduction of the gradient diffusion approximation for the turbulent transport and pressure diffusion terms.

The \(k\) equation used in the \(k\)-\(\tau\) model is then obtained as a simple transformation of the standard \(k\) equation by the following relationship [Kok],

\[\omega\equiv\frac{\epsilon}{\beta^*k}\]

where \(\beta^*\) is a constant. Inserting \(\omega\beta^*k\) for \(\epsilon\) in the dissipation term \(\rho\epsilon\) gives the \(k\) equation used in nekRS,

\[\frac{\partial(\rho k )}{\partial t}+\nabla\cdot\left(\rho k\overline{\mathbf u}\right)=\nabla\cdot\left\lbrack\left(\mu+\frac{\mu_T}{\sigma_k}\right)\nabla k\right\rbrack+\mathscr{P}-\rho\beta^*\frac{k}{\tau}\]

The \(\tau\) Equation

In two-equation RANS turbulence modeling, the greatest source of uncertainty is the proper choice of the second transport variable. While a true \(k\) equation is often used as the starting point for developing the model \(k\) equation, it is commonplace to start immediately from an ad hoc, “fabricated,” model equation for the second turbulence variable. Of course, “exact” equations can always be derived for the second turbulence variable through various operations on the mean Navier-Stokes equation or the Reynolds stress equation, but the exact equations for \(\epsilon\), \(\omega\), or other turbulence quantities tend to be far more complex than the exact equation for \(k\) shown earlier.

In 1942, Kolmogorov was the first to propose the \(k\)-\(\omega\) model [Wilcox]. His formulation was very heuristic - from the Boussinesq approximation, it is likely that \(\nu_T\propto k\), which requires another variable with dimensions inverse time. Based on the work of Kolmogorov and many subsequent researchers of the \(k\)-\(\omega\) model, inserting \(\tau\equiv 1/\omega\) into the \(\omega\) equation gives the \(\tau\) equation - this approach is very similar to that used to obtained the \(k\) equation. The \(\tau\) equation used in nekRS is [Kok]

\[\frac{\partial(\rho\tau)}{\partial t}+\nabla\cdot\left(\rho\tau\overline{\mathbf u}\right)=\nabla\cdot\left\lbrack\left(\mu+\frac{\mu_T}{\sigma_\tau}\right)\nabla \tau\right\rbrack-\alpha\frac{\tau}{k}\mathscr{P}+\rho\beta-2\frac{\mu}{\tau}\nabla\tau\cdot\nabla\tau\]

where \(\sigma_\tau\), \(\alpha\), and \(\beta\) are constants. The last term on the right-hand side of the \(\tau\) equation is in practice implemented in the form

\[\frac{2}{\tau}\nabla\tau\cdot\nabla\tau\rightarrow 8\nabla\sqrt{\tau}\cdot\nabla\sqrt{\tau}\]

in order to reduce the discretization error associated with the computation of gradients of a term that scales as \(y^2\) as \(y\rightarrow 0\) [Kok].

The Eddy Viscosity

The objective of RANS models is to estimate the eddy viscosity \(\mu_T\) that appears in the Boussinesq approximation. The particular form for \(\mu_T\) can be understood here in terms of the standard \(k\)-\(\epsilon\) model [Launder], for which \(\mu_T\) is given as

\[\mu_T=C_\mu\rho\frac{k^2}{\epsilon}\]

where \(C_\mu\) is a constant. Inserting \(\tau\equiv 1/\omega\) and \(\epsilon=\beta^*\omega k\) gives

\[\mu_T=\rho k\tau\]

which presumes that \(C_\mu\) and \(\beta^*\) are really the same constant, but with different notation developed separately by the \(k\)-\(\epsilon\) researchers and the \(k\)-\(\tau\) researchers [Kok].

Closure Coefficients and Other Details

Table RANS Coefficients shows the values for the various constants used in nekRS’s \(k\)-\(\tau\) model.

Table 1 RANS Coefficients

Coefficient	Value	Source
\(\sigma_k\)	\(5/3\)
\(\sigma_\tau\)	\(2.0\)
\(Pr_T\)	user-selected	—

A limiter is applied to both \(k\) and \(\tau\) to prevent negative values of either \(k\) or \(\tau\),

\[k = \max{\left(k, 0.01|k|\right)}\]

\[\tau = \max{\left(\tau, 0.01|\tau|\right)}\]

Warning

nekRS’s \(k\)-\(\tau\) implementation currently requires that the laminar dynamic viscosity and the density are constant, because the setup routines can only accept constant values. See RANS Plugin for more information.

Note

Even if the molecular viscosity is constant, you must set stressFormulation = true in the input file because the total viscosity (molecular plus turbulent) will not be constant.

Boundary Conditions

On walls, because the asymptotic behavior of \(\omega\) is \(\omega\propto y^{-2}\) as \(y\rightarrow0\), and because the instantaneous velocity \(u_i\equiv \overline{u_i}+u_i'\) must be zero due to the no-slip condition, both \(k\) and \(\tau\) should be set to zero on no-slip boundaries.

On turbulent inlets, however, both \(k\) and \(\tau\) are generally nonzero, and estimates for \(k\) and \(\tau\) must be provided. Turbulent inflow conditions are usually unknown unless the modeler is fortunate enough to have experimental data - therefore, the specification of inlet conditions on \(k\) and \(\tau\) tends to be fairly ad hoc. An inlet condition for \(k\) can be estimated by prescribing the turbulent intensity at the inlet. The turbulent intensity \(I\) is defined as the root mean square fluctuating velocity normalized by the magnitude of the mean velocity, or

\[I\equiv\frac{\sqrt{\frac{1}{3}\left(u_i'u_i'\right)}}{\sqrt{\overline{u_j}\ \overline{u_j}}}\]

which can equivalently be written in terms of the turbulent kinetic energy as

\[I\equiv\frac{\sqrt{\frac{2}{3}k}}{\sqrt{\overline{u_j}\ \overline{u_j}}}\]

Therefore, if an inlet turbulent intensity can be prescribed, the inlet turbulent kinetic energy is

\[k\equiv\frac{3}{2}I^2\overline{u_j}\ \overline{u_j}\]

For instance, it is common to assume a uniform turbulent intensity over an inlet of between 1 and 5% for pipe flows [Russo]. Experiments have also quantified the scaling of turbulent intensity with Reynolds number. From simulations and experiments of both incompressible and compressible flows, the turbulent intensities on the axis of a smooth circular pipe are [Russo]:

\[\begin{split}I_\text{axis}=\begin{cases}0.0853Re^{-0.0727} & \text{incompressible}\\ 0.0550Re^{-0.0407} & \text{compressible}\end{cases}\end{split}\]

With a slightly different definition based on the turbulent intensity averaged over the cross-sectional area of a circular pipe, the turbulent intensity instead scales as [Russo]:

\[\begin{split}I_\text{area}=\begin{cases}0.140Re^{-0.0790} & \text{incompressible}\\ 0.227Re^{-0.1} & \text{compressible}\end{cases}\end{split}\]

Similar correlations have also been developed for rough pipes [Basse]. It is also commonplace to apply conditions on \(k\) in terms of the friction velocity \(u_\tau\), defined as

\[u_\tau\equiv\sqrt{\frac{\tau_w}{\rho}}\]

where \(\tau_w\) is the wall shear stress. The wall shear stress can then be estimated using a friction factor correlation for \(f_D\), the Darcy friction factor, which is defined as

\[\frac{\Delta P}{l}\equiv f_D\frac{1}{D}\frac{\rho \overline{u_j}\overline{u_j}}{2}\]

where \(\Delta P/l\) is a pressure gradient and \(D\) is a hydraulic diameter. For circular pipes, the friction factor is related to \(\tau_w\) as

\[f_D=8\frac{\tau_w}{\rho\overline{u_j}\overline{u_j}}\]

Typically, a simple friction factor correlation such as the Blasius model for pipes, is selected

\[f_D=0.316Re^{-0.25}\]

Combining \(u_\tau\), \(f_D\), and one of the previous estimates for \(I_\text{area}\), such as \(I_\text{area}=0.317Re^{-0.11}\) [Basse] gives the following relationship between \(I_\text{area}\) and \(f_D\),

\[I_\text{area}=0.526f_D^{0.44}\]

Then, inserting the relationship between \(u_\tau\) and \(f_D\) gives, after some manipulation,

\[k=2.56 u_\tau^2\left(\overline{u_j}\ \overline{u_j}\right)^{0.24}\]

While this algebraic exercise hasn’t actually introduced any new information or closures beyond what has already been discussed, it is information to see \(k\) expressed in terms of the friction velocity, because some nekRS inputs set \(k\) on inlets by first computing the \(u_\tau\) from a friction correlation and then using an expression similar to above.

In any case, one of these correlations for turbulent intensity, or simply a fixed turbulent intensity of, say, 5%, can be used to prescribe a uniform value of \(k\) on an inlet. However, including some spatial variation in \(k\) on the inlet may reduce Gibbs phenomena if the inlet turbulent intensity enforces the physical zero wall value. Spatial fits, such as those developed for circular pipes [Russo], may improve the numerical stability and accuracy of your simulation.

Finally, on outlets, “free-stream” boundary conditions are usually applied to \(k\) and \(\tau\).

Non-Dimensional Formulation

Now that the \(k\)-\(\tau\) model has been presented in its full form, the non-dimensional formulation is provided. Because nekRS’s \(k\)-\(\tau\) model is currently limited to constant densities and laminar viscosities, the non-dimensional formulation is slightly simpler than the more general case shown in Non-Dimensional Formulation that derived the non-dimensional form of the instantaneous Navier-Stokes equations.

Introduce the non-dimensional variables \(\mathbf x^\dagger=\frac{\mathbf x}{L}\), \(\overline{\mathbf u}^\dagger=\frac{\overline{\mathbf u}}{U}\), \(t^\dagger=\frac{tU}{L}\), and \(\overline{\mathbf f}^\dagger=\frac{\overline{\mathbf f} L}{U^2}\). For convection-dominated flows, the pressure is non-dimensionalized in terms of the dynamic pressure as \(P^\dagger=\frac{P}{\rho U^2}\). Finally, the turbulent kinetic energy and inverse dissipation rate are non-dimensionalized as \(k^\dagger=\frac{k}{U^2}\) and \(\tau^\dagger=\frac{\tau U}{L}\).

Inserting these non-dimensional variables into the mean flow mass and momentum conservation equations gives

\[\frac{\partial \overline{u_i}^\dagger}{\partial x_i^\dagger}=0\]

\[\frac{\partial \overline{u_i}^\dagger}{\partial t^\dagger}+\overline{u_j}^\dagger\frac{\partial \overline{u_i}^\dagger}{\partial x_j^\dagger}=-\frac{\partial \overline{P}^\dagger}{\partial x_i^\dagger}+\frac{1}{Re}\left(1+\frac{\mu_T}{\mu}\right)\frac{\partial\overline{\tau_{ij}}^\dagger}{\partial x_j^\dagger}-\frac{\partial}{\partial x_j^\dagger}\left(\frac{2}{3}k^\dagger\delta_{ij}\right)+f_i^\dagger\]

To non-dimensionalize the energy conservation equation, use the previous non-dimensional variables in addition to a non-dimensional temperature, \(T^\dagger=\frac{T-T_0}{\Delta T}\), where \(\Delta T\) is a reference temperature rise relative to a baseline temperature \(T_0\). The heat source is non-dimensionalized as \(\dot{q}^\dagger=\frac{\dot{q}}{\rho C_{p} U\Delta T/L}\), which arises naturally from the simple formulation of bulk energy conservation of \(Q=\dot{m}C_p\Delta T\), where \(Q\) is a heat source (units of Watts) and \(\dot{m}\) is a mass flowrate. Inserting these non-dimensional variables into the energy conservation equation gives

\[\frac{\partial \overline{T}^\dagger}{\partial t^\dagger}+\overline{u_i}^\dagger\frac{\partial \overline{T}^\dagger}{\partial x_i^\dagger}=\frac{1}{Pe}\left(1+\frac{\mu_T/Pr_T}{k}C_p\right)\frac{\partial}{\partial x_i^\dagger}\frac{\partial \overline{T}^\dagger}{\partial x_i^\dagger}+\dot{q}^\dagger\]

To non-dimensionalize the \(k\) and \(\tau\) equations, define \(\mathscr{P}^\dagger=\frac{\mathscr{P}}{\rho U^3/L}\) and \(\epsilon^\dagger=\frac{\epsilon}{U^3/L}\). With previous non-dimensional variables already defined, the non-dimensional \(k\) and \(\tau\) equations become

\[\frac{\partial k^\dagger}{\partial t^\dagger}+\overline{u_i}^\dagger\frac{\partial k^\dagger}{\partial x_i^\dagger}=\frac{1}{Re}\left(1+\frac{\mu_T/\sigma_k}{\mu}\right)\frac{\partial}{\partial x_i^\dagger}\frac{\partial k^\dagger}{\partial x_i^\dagger}+\mathscr{P}^\dagger-\beta^*\frac{k^\dagger}{\tau^\dagger}\]

\[\frac{\partial\tau^\dagger}{\partial t^\dagger}+\overline{u_i}^\dagger\frac{\partial\tau^\dagger}{\partial x_i^\dagger}=\frac{1}{Re}\left(1+\frac{\mu_T/\sigma_\tau}{\mu}\right)\frac{\partial}{\partial x_i^\dagger}\frac{\partial\tau^\dagger}{\partial x_i^\dagger}-\alpha\frac{\tau^\dagger}{k^\dagger}\mathscr{P}^\dagger+\beta-\frac{2}{Re}\frac{1}{\tau^\dagger}\frac{\partial\tau^\dagger}{\partial x_i^\dagger}\frac{\partial\tau^\dagger}{\partial x_i^\dagger}\]

Finally, the eddy viscosity is computed as

\[\mu_T=\rho k\tau\]

which after inserting the non-dimensional variables becomes

\[\mu_T=\rho ULk^\dagger\tau^\dagger\]

such that the non-dimensional eddy viscosity can be written as \(\mu_T^\dagger=k^\dagger\tau^\dagger\). Therefore, the overall diffusion coefficients in the mean momentum equation, mean energy equation, \(k\) equation, and \(\tau\) equation are, respectively

\[\frac{1}{Re}\left(1+\frac{\mu_T}{\mu}\right)\rightarrow\frac{1}{Re}+\mu_T^\dagger\]

\[\frac{1}{Pe}\left(1+\frac{\mu_T/Pr_T}{k}C_p\right)\rightarrow\frac{1}{Pe}+\frac{\mu_T^\dagger}{Pr_T}\]

\[\frac{1}{Re}\left(1+\frac{\mu_T/\sigma_k}{\mu}\right)\rightarrow\frac{1}{Re}+\frac{\mu_T^\dagger}{\sigma_k}\]

\[\frac{1}{Re}\left(1+\frac{\mu_T/\sigma_\tau}{\mu}\right)\rightarrow\frac{1}{Re}+\frac{\mu_T^\dagger}{\sigma_\tau}\]