Turbulence modeling

Introduction

The MOOSE Navier-Stokes module includes experimental turbulence modeling capabilities. Presently, the models are insufficient for stand-alone predictive simulation, and we recommend that the user tunes the model parameters for their problem of interest using experimental data or a higher-fidelity code for reference solutions.

Reynolds averaging

One common approach in turbulent flow modeling is to use the Reynolds averaging procedure. In this procedure, all of the relevant variables are decomposed into an average value and a fluctuating value. For example, velocity will be decomposed as where is the average value and is the fluctuating component. It is difficult to rigorously define what is meant by "average" and what is meant by "fluctuating", but these concepts can be roughly defined as: the average component is the component that we intend to explicitly resolve with the numerical solver and the fluctuating component is the remainder which varies on a scale that is smaller than the solver mesh size and/or faster than the solver time step. For a more detailed discussion, refer to the books by Bailly and Comte-Bellot (Bailly and Comte-Bellot, 2015) or Moukalled et al. (Moukalled et al., 2016).

The averaging procedure can also be understood by the relations that it obeys which include,

Equations for turbulent flow can be derived by applying the Reynolds averaging procedure to the relevant equations. As an example, consider the conservation of linear momentum in the -direction which can be written as, where is the fluid density, is the velocity, and is the -component of the force felt by the fluid.

The first step in the Reynolds averaging process is to apply the decomposition. Here the analysis will be restricted to incompressible flows which means so . Applying the decomposition to the other variables gives,

Because the solver will only resolve the average behavior, the next step is to apply the averaging operator to each term, (1)

The quantity, , seen in Eq. (1) is known as the Reynolds stress. It models a change in the resolved momentum due to a covariance in the unresolved velocity field.

A model is needed to relate the Reynolds stress to known quantities in the system (like the average velocity). One model frequently used in engineering analysis is the Boussinesq hypothesis which treats the Reynolds stress with a functional form similar to the one used for viscous stress.

Boussinesq hypothesis and turbulent momentum transfer

Forces due to viscous friction can be analyzed in terms of a viscous stress tensor, ,

For a Newtonian fluid, the stress tensor is given by, where is the dynamic viscosity, is the bulk viscosity, and is the identity matrix.

The Boussinesq hypothesis models the Reynolds stress in the analogous form, (Bailly and Comte-Bellot, 2015; Moukalled et al., 2016) where is the eddy viscosity and is the turbulent kinetic energy which is itself defined by,

The Boussinesq hypothesis can also be stated in terms of , the eddy diffusivity for momentum as, (2)

Note that is analogous to the kinematic viscosity, , so is also referred to as the eddy viscosity in some sources.

As a practical matter, the second term on the right-hand-side of Eq. (2) is frequently neglected and implicitly subsumed into the pressure term of the momentum equation. The new pressure variable can be referred to as the turbulent pressure instead of the thermodynamic pressure (Moukalled et al., 2016). Neglecting this term results in the model, (3)

This stress model is implemented in the INSFVMixingLengthReynoldsStress kernel.

Turbulent energy transfer

A model for the effects of turbulence on the conservation of energy can be derived by following the same process that was used for the conservation of momentum. First, the Reynolds-averaging procedure is applied to the conservation equation. Then, the flux due to the fluctuating covariance term is replaced with a model analogous to the one used for molecular processes.

The conservation of energy can be expressed as where is the specific heat capacity, is the heat flux due to molecular conduction, and captures any sources and losses of heat. Note that this equation only conserves energy for the case of constant and .

The conduction term can be modeled with Fourier's law as where is the thermal conductivity.

Substituting in Fourier's law gives the form,

The Reynolds averaging procedure can now be applied to find,

As before with the Boussinesq hypothesis for turbulent momentum transfer, we can model turbulent heat transfer with an analogy to the molecular mechanism, where is the eddy diffusivity for heat (Lienhard IV and Lienhard V, 2020).

The diffusivity for heat is related to the diffusivity for momentum by, where is the turbulent Prandtl number (Lienhard IV and Lienhard V, 2020).

Turbulent transfer of trace species

Some problems also require modeling the transport of trace chemicals (e.g. a dye) in the fluid. Here, the analysis is limited to very dilute substances that have no significant impact on bulk fluid properties like and .

The conservation of these passive species can be expressed as where is the concentration of the species in number of molecules (or moles) per unit volume, is the molecular diffusion coefficient (Fick's law has already been applied in this equation), and is an unspecified source and/or loss term.

Applying the Reynolds averaging procedure leads to the equation,

The molecular analogy is then, where is the turbulent diffusivity for the passive species which is in turn modeled with the relationship, where is the turbulent Schmidt number (Tominaga and Stathopoulos, 2007).

Mixing-length models

There are many models used in CFD codes to compute the eddy viscosity, . Some of the most popular are the isotropic 2-equation models, - and -. However, the MOOSE Navier-Stokes module only supports simple algebraic mixing-length models at present.

Derivations of mixing length models generally consider 2D geometry of the near-wall region. Let be a direction tangential to the wall and be the direction perpendicular to the wall. With this notation, the eddy diffusivity can be modeled as, (Bailly and Comte-Bellot, 2015; Todreas and Kazimi, 2011) where is the mixing length, a characteristic length for turbulent eddies in the fluid. The mixing length will be discussed further below.

This model is useful for analysis of simple geometries, but it is inconvenient to apply for general meshes where the walls may take an arbitrary shape. This is because there is no clear way to define directions parallel and perpendicular to the wall in the general case (e.g. consider a mesh cell with walls on two sides). Consequently, we use Smagorinsky's velocity scale:

where

This momentum diffusivity model is implemented in the INSFVMixingLengthReynoldsStress kernel. The corresponding model for diffusivity of passive scalars (like energy) is implemented in the INSFVMixingLengthScalarDiffusion kernel.

A model is then needed for the mixing length itself. One popular choice is to assume the mixing length is proportional to the distance from the nearest wall, (Bailly and Comte-Bellot, 2015; Todreas and Kazimi, 2011) where is the wall-distance and is known as the von Kármán constant. A von Kármán of is often used for the near-wall region (Todreas and Kazimi, 2011).

A modification to this model was done by (Escudier, 1966), who claims that the mixing length grows linearly in the boundary layer region and then it takes a constant value. This prevents an excessive growth of Prandtl's original mixing length model. The equations for the mixing length are then

where is the Von Karman constant, as in Escudier's model and has length units and represents the thickness of the velocity boundary layer. Note that for large values of , the mixing length in Prandtl's original model is obtained.

This mixing length is implemented in the WallDistanceMixingLengthAux auxiliary kernel. Note that the wall-distance calculation can be expensive for large meshes. If the mesh is constant in time, this cost can be amortized by setting the execute_on parameter to initial so that the wall distance is computed only once at the beginning of the simulation.

Tuning the mixing length

In practice, the mixing length distribution is highly problem-dependent. There are also known deficiencies with the mixing length model itself. For example, this model predicts no turbulent mixing where the velocity gradient is zero.

Consequently, we recommend that users compare results generated with this model against reference solutions (either from experiment or from software with higher fidelity models). These comparisons can also be used to tune an appropriate mixing length model to the system of interest.

The examples directory includes a simple circular pipe problem and shows how can be tuned so that the simulated pressure drop matches a correlation. Also note that other MOOSE auxkernels can be used to implement different mixing length models. For example, if the near-wall region and the core region of the pipe are assigned to different mesh blocks, then a WallDistanceMixingLengthAux can be combined with a ConstantAux to model a mixing length that is proportional to the wall-distance in the near-wall region and constant elsewhere.

Wall Function Boundary Condition for Velocity

In wall-bounded flows, the velocity gradients present in the near wall region are quite steep. Therefore, the mesh resolution needed to capture these gradients can become expensive. To economize computer time and storage, wall functions have been proposed, which are equations empirically derived to satisfy the physics in the near wall region under certain assumptions.

Experimental and dimensional analysis shows that for high values () the wall shear stress is related to the mean velocity parallel to the wall through the so-called logarithmic law of the wall: (Moukalled et al., 2016) (Launder and Spalding, 1983)

where:

• is the dimensionless wall-tangential velocity component.

• is the wall-tangential velocity component at the centroid of the adjacent cell to the wall.

• is the friction velocity.

• is the dimensionless wall distance

• is the von Karman constant

• is the log law offset

The log law offset is set for smooth walls to a value of 9.0 (Launder and Spalding, 1983).Theoretically, in the viscous sublayer for low values (), the functional form of the velocity is linear . The log law offset is obtained by matching the viscous sublayer profile with the logarithmic law of the wall at . (Launder and Spalding, 1983). Rough walls are not currently supported. Nonetheless, this could be done through a modification of the log law offset.

The logarithmic law of the wall provides an implicit equation for the magnitude of the wall shear stress. The orientation of the wall shear stress is dictated by the orientation of the tangential velocity, which is parallel to the wall. Its sense is opposite to the velocity's sense. Its implementation is present in INSFVWallFunctionBC.

This shear stress is used in solving the momentum equation by invoking its value as a flux source term. The second condition needed is the non-penetrability of the velocity through the wall, such that the mass flux through that face is zero (Segal et al., 1993).

The location of the first cell centroid must be such that . Auxiliary Kernel WallFunctionYPlusAux can be used to obtain the value. Another Auxiliary Kernel allows to obtain the wall shear stress, WallFunctionWallShearStressAux.

This standard wall function was originally calibrated for fully developed turbulent flows in steady state, in absence of an external force and an adverse pressure gradient. It may not provide accurate results under separation, recirculation or low Reynolds number flows.

References