The k-ε Turbulence Model
The k-ε model is one of the most used turbulence models for industrial applications. This module includes the standard k-ε model (Ref. 1). The model introduces two additional transport equations and two dependent variables: the turbulent kinetic energy, k, and the turbulent dissipation rate, ε. The turbulent viscosity is modeled as
(3-76)
where Cμ is a model constant.
The transport equation for k reads:
(3-77)
where the production term is
(3-78)
The transport equation for ε reads:
(3-79)
The model constants in Equation 3-76, Equation 3-77, and Equation 3-79 are determined from experimental data (Ref. 1) and the values are listed in Table 3-4.
Table 3-4: Model Constants
Cμ
Cε1
Cε2
σk
σε
Mixing Length Limit
Equation 3-77 and Equation 3-79 cannot be implemented directly as written. There is, for example, nothing that prevents division by zero. The equations are instead implemented as suggested in Ref. 10. The implementation includes an upper limit on the mixing length, :
(3-80)
The mixing length is used to calculated the turbulent viscosity. should not be active in a converged solution but is merely a tool to obtain convergence.
Realizability Constraints
The eddy-viscosity model of the Reynolds stress tensor can be written
where δij is the Kronecker delta and Sij is the strain-rate tensor. The diagonal elements of the Reynolds stress tensor must be nonnegative, but calculating μT from Equation 3-76 does not guarantee this. To assert that
the turbulent viscosity is subjected to a realizability constraint. The constraint for 2D and 2D axisymmetry without swirl is:
(3-81)
and for 3D and 2D axisymmetry with swirl flow it reads:
(3-82)
Combining equation Equation 3-81 with Equation 3-76 and the definition of the mixing length gives a limit on the mixing length scale:
(3-83)
Equivalently, combining Equation 3-82 with Equation 3-76 and Equation 3-80 gives:
(3-84)
This means there are two limitations on lmix: the realizability constraint and the imposed limit via Equation 3-80.
The effect of not applying a realizability constraint is typically excessive turbulence production. The effect is most clearly visible at stagnation points. To avoid such artifacts, the realizability constraint is always applied for the RANS models. More details can be found in Ref. 4, Ref. 5, and Ref. 6.
Model Limitations
The k-ε turbulence model relies on several assumptions, the most important of which is that the Reynolds number is high enough. It is also important that the turbulence is in equilibrium in boundary layers, which means that production equals dissipation. These assumptions limit the accuracy of the model because they are not always true. It does not, for example, respond correctly to flows with adverse pressure gradients and can result in under-prediction of the spatial extent of recirculation zones (Ref. 1). Furthermore, in simulations of rotating flows, the model often shows poor agreement with experimental data (Ref. 2). In most cases, the limited accuracy is a fair tradeoff for the amount of computational resources saved compared to using more complicated turbulence models.
Wall Functions
The flow close to a solid wall is for a turbulent flow very different from the free stream. This means that the assumptions used to derive the k-ε model are not valid close to walls. While it is possible to modify the k-ε model so that it describes the flow in wall regions (see The Low Reynolds Number k-ε Turbulence Model), this is not always desirable because of the very high resolution requirements that follow. Instead, analytical expressions are used to describe the flow near the walls. These expressions are known as wall functions.
When using wall functions in COMSOL Multiphysics, a theoretical lift-off from the physical wall is assumed as shown in Figure 3-7. Expressed in viscous units, the wall lift-off is defined as
The first argument is derived from the law of the wall. The second argument is the distance from the wall, in viscous units, where the logarithmic layer meets the viscous sublayer (or rather would meet it if there were no buffer layer in between). This lower limit ensures that the wall functions remain non-singular for all Reynolds numbers.
The wall lift-off, δw, is defined as
where the friction velocity, uτ, is defined by
where in turn, κv, is the von Kárman constant (default value 0.41) and B is a constant that by default is set to 5.2. The two arguments for uτ are under some addition assumptions theoretically identical in the logarithmic layer (Ref. 1), but deviate in stagnation points and when the local Reynolds number becomes low. The definition of δw is such that it becomes h/2 when , but it can become larger when the lower limit for , 11.06, takes effect.
Wall functions give reasonable predictions as long as is lower than some upper limit that depends on the turbulent Reynolds number (Ref. 28). The upper limit is hardly ever lower than 50, and in many practical applications as high as a few hundred. Highest accuracy is obtained if is also everywhere larger than 25, which approximately corresponds to the beginning of the logarithmic layer. It can also be worthwhile to check that δw is small compared to the dimensions of the geometry, especially if for a significant fraction of the wall area.
The boundary conditions for the velocity is a no-penetration condition u n = 0 and a shear stress condition
where
is the viscous stress tensor.
The turbulent kinetic energy is subject to a homogeneous Neumann condition n ⋅ ∇= 0 and the boundary condition for ε reads:
See Ref. 10 and Ref. 11 for further details.
Wall Functions for Rough Walls
The physics interfaces Wall and Interior Wall have an option to apply wall roughness by modifying the wall functions. Cebeci (Ref. 21) suggested a model which adjusts the friction velocity for surface roughness,
(3-85)
where
is the roughness height in viscous units,
The roughness height, ks, is the peak-to-peak value of the surface variations and the wall is relocated to their mean level.
Figure 3-8: Definitions of the roughness height and the modified wall location.
Hence, when Equation 3-85 is used the lift-off is modified according to,
where h+ is the height of the boundary mesh cell in viscous units. Cs is a parameter that depends on the shape and distribution of the roughness elements. When the turbulence parameters κν and B have the values 0.41 and 5.2, respectively, and Cs = 0.26, ks corresponds to the equivalent sand roughness height, kseq, as introduced by Nikuradse (Ref. 22). A few characteristic values of the equivalent sand roughness height are given in Table 3-5 below,
50 μm
Use other values of the roughness parameter Cs and roughness height ks to specify generic surface roughnesses.
Initial Values
The default initial values for a stationary simulation are (Ref. 10),
where is the mixing length limit. For time-dependent simulations, the initial value for k is instead
Scaling for Time-Dependent Simulations
The k-ε equations are derived under the assumption that the flow has a high enough Reynolds number. If this assumption is not fulfilled, both k and ε have very small magnitudes and behave chaotically in the manner that the relative values of k and ε can change by large amounts due to small changes in the flow field.
A time-dependent simulation of a turbulent flow can include a period when the flow is not fully turbulent. A typical example is the startup phase when for example an inlet velocity or a pressure difference is gradually increased. To sort out numerical fluctuations in k and ε during such periods, the default time-dependent solver for the k-ε model employs unscaled absolute tolerances for k and ε. The tolerances are set to
(3-86)
where Uscale and Lfact are input parameters available in the Advanced Settings section of the physics interface node. Their default values are 1 m/s and 0.035 respectively. lbb,min is the shortest side of the geometry bounding box. Equation 3-86 is closely related to the expressions for k and ε on inlet boundaries (see Equation 3-149).
The practical implication of Equation 3-86 is that variations in k and ε smaller than kscale and εscale respectively, are regarded as numerical noise.