Most of nowadays CFD simulations are conducted with the Reynolds Averaging approach. Reynolds Averaged Navier-Stokes (RANS) simulation, the “working horse” of industrial CFD is based on the Reynolds decomposition according to which a flow variable is decomposed into mean and fluctuating quantities. When the decomposition is applied to Navier-Stokes equation an extra term known as the Reynolds Stress Tensor arises and a modelling methodology is needed to close the equations. The “closure problem” is apparent as higher and higher moments of the set of equations may be taken, more unknown terms arise and the number of equations never suffices. This is of course an obvious consequence to the fact that taking these higher moments is simply a mathematical endeavor and has no physical contribution what so ever.
THEME FIGURE (NSE IN INTEGRAL FORM SCRIBBLED ON A NOTE) TAKEN FROM LINKEDIN AND CREATED BY THE CREATIVE AARON PERRY (presumably by rumor while paying his restaurant bill).
The credit for taking the picture goes for James Shaw while on vacation a year or so ago , for realizing a cool opportunity.
Levels of modeling are related to the number of differential equations added to Reynolds Averaged Navier-Stokes equations in order to “close” them.
0-equation (algebraic) models are the simplest form of turbulence models, a turbulence length scale is specified in advance through experimenting. 0-equations models are very limited in applications as they fail to take into account history effects, assuming turbulence is dissipated where it’s generated, a direct consequence of their algebraic nature.
1-equation and 2-equations models, incorporate a differential transport equation for the turbulent velocity scale (or the related the turbulent kinetic energy) and in the case of 2-equation models another transport equation for the length scale (or time scale), subsequently invoking the “Boussinesq Hypothesis” relating an eddy-viscosity analog to its kinetic gasses theory derived counterpart (albeit flow dependent and not a flow property) and relating it to the Reynolds stress through the mean strain.
In this sense 2-equation models can be viewed as “closed” because unlike 0-equation and 1-equation models (with exception of 1-equations transport for the eddy viscosity itself as described in a former post: Understanding The Spalart-Allmaras Turbulence Model) these models possess sufficient equations for constructing the eddy viscosity with no direct use for experimental results.
2-equations models do however contain many assumptions along the way for achieving the final form of the transport equations and as such are calibrated to work well only according to well-known features of the applications they are designed to solve. Nonetheless although their inherent limitations, today industry need for rapid answers dictates CFD simulations to be mainly conducted by 2-equations models whose strength has proven itself for wall bounded attached flows at high Reynolds number (thin boundary layers) due to calibration according to the law-of-the-wall.
The Law of the Wall
In what follows I shall present a fairly heuristic description of the most widely-used concepts from the classical theory of turbulence
The velocity profile for a fully-developed (time mean) turbulent flow in a duct is very different from the laminar profile:
Shear stress at the wall is calculated as:
It can be seen that the velocity gradient at the wall, and hence also the wall shear stress, corresponding to a fully-developed Poiseuille flow in a duct is much larger for the turbulent case than for the laminar one. It can also be understood just from looking at the turbulent flow velocity profile that there must be at least two length scales associated with this flow.
The first corresponds to the rather thin region adjacent to the walls in which the velocity profile is nonuniform, with large gradients (keep in mind the non-slip boundary condition), and exhibits a near linear velocity profile. This layer, dominated by viscosity effects is termed the viscous sublayer. In light of the dominance of the viscous effects it is often associated with the Kolmogorov scale (These scales were predicted on the basis of dimensional analysis as part of the K41 theory, corresponding to small length scales large wavenumbers where kinetic energy of fluid motion is converted to thermal energy – mind that such association is quite subtle). The outer region shows a near constant velocity with distance from the wall.
I will note though, that following some more rigor motivating the more heuristic treatment), relying on perturbation analysis (a powerful mathematical tool intended on matching solutions in the form of an asymptotic expansions in our case) it could be shown that the inner profile (nearly linear) which satisfies the no-slip boundary condition at the wall does not correctly asymptote the outer (nearly uniform) solution which at the same time cannot satisfy the no-slip boundary condition, hence, a third solution is to match them (derived by means of asymptotic expansion). For the heuristic treatment, somewhat following H. Tennekes and J. L. Lumley. A First Course in Turbulence the above formalism is not required.
and we denote a velocity scale for the inner region as:
The first is associated with a large advective length scale (as the hight of the duct):
and the latter with small viscous length scale which by pure dimensional analysis we can express through:
In order for an intermediate scale to be realizable, the ratio of the two length scales must be large:
Now we identify a range distance y from the wall for which:
within this range the heuristic argument claims that ν/uτ is to small enough and h is large enough to not be able control the flow dynamics for the first nor result in direct interaction for the latter. In other words we claim that y is the only length scale corresponding for this region but there are two velocity scales:
Basing only on dimensional analysis we relate the two according to:
with the constant determined from experimental data. We define dimensionless quantities:
so the above relation between the two velocities becomes:
Now we integrate to achieve the following:
This is the well-know “log-law” matching the inner and outer layer (the second constant added is an integration constant also to be evaluated experimentally).
The friction velocity chosen earlier is defined as:
Below the law of the wall is finally displayed clearly and we can see four different length scales:
We can also distinguish what we actually mean by the “inner portion of the turbulent boundary-layer”, of which we shall focus most of our attention upon trying to obtain acceptable CFD-based predictions of certain flow applications:
- Viscous sublayer – the portion of the boundary layer of which y+<5, and is often also called the laminar sublayer. While dominated by viscous effects, this does not mean that the flow is laminar here, because that would imply that there are no fluctuations whatsoever. In the viscous sublayer velocity fluctuations do occur, but they are induced by the turbulence above the viscous sublayer.
In the viscous sublayer we get much more dissipation than production of turbulent kinetic energy.
- buffer layer – smoothly connects the viscous sublayer to the inertial sublayer and might be thought of as a region of which the dissipation and inertial effects are somewhat balanced.
- Inertial sublayer (log-layer) – where the log profile holds.
If we plot the different tems in the turbulent kinetic energy transport equation which is part of the most popular 2-equation closure models (such as k-ω SST and k-ε ), we find that the turbulent kinetic energy production and dissipation are nearly equal in this region, and that the diffusion is almost zero.
- Defect layer – sometimes simply called the “outer layer,” is the region of a turbulent boundary layer beyond the log layer (see below), and beginning approximately one-tenth the boundary layer thickness from the wall.
The law of the wall may serve during CFD analysis in prior estimation of grid requirement as explained below:
Near Wall Modeling
Near wall modeling is arguably the most problematic area in turbulence modeling. Dealing with near wall modeling means focusing on the turbulent boundary-layer.
Walls are the main source of vorticity in turbulence. We wouldn’t have any turbulence, if it wouldn’t have a wall generating it. The presence of walls gives rise to turbulent momentum boundary layers, which have particularly very steep gradients in the inner-most portion of the turbulent boundary-layer.
Successful prediction of viscous drag in aerodynamic flows, pressure drop due to separation, and heat transfer are merely a few examples critically influenced by how we decide to model this inner-most portion of the turbulent boundary-layer.
However, from a practitioner standpoint, it’s not just the modeling choice which is important, but also what kind of mesh resolution is needed in the boundary layer to achieve acceptable predictions. The model of choice and the grid resolution are tightly coupled in the sense that a very good model for a specific application will achieve acceptable predictions only if the mesh resolution is such that it matches the specific model.
Shortcomings In Wall Treatment for the Standard k-ε Turbulence Model
Although being perhaps the most popular and known turbulence model, the standard k-ε turbulence model carries along some harsh shortcomings which are important to acknowledge.
First, it is important to note that the model is essentially a high Reynolds model, meaning the law of the wall must be employed and provide velocity “boundary conditions” away from solid boundaries (what is termed “wall-functions”). From a mathematical standpoint, even if one could impose Dirichlet conditions for ε on solid boundary, after meshing it would still be difficult to numerically approach the problem due to what is termed in numerical analysis as stiffness of the numerical problem, partially related to the high gradients.
In order to integrate the equations through the viscous/laminar sublayer a “Low Reynolds” approach must be employed. This is achieved as additional highly non-linear damping functions are needed to be added to low-Reynolds formulations (low as in entering the viscous/laminar sublayer) to be able to integrate through the laminar sublayer (y+<10). This again produces numerical stiffness and in case is problematic to handle in view of linear numerical algorithms.
Furthermore the low Reynolds methodology should not be confused with transitional Reynolds modeling, as sometimes many practitioners rely on low Reynolds models to achieve a measure of transition from laminar to turbulence prediction. As the low Reynolds methodology is devised to handle the near wall viscous/laminar sublayer, there is no reason to expect it also satisfactory transition predictions, especially noting that there are many mechanisms for transition onset. Perhaps the only predictions which shall be close to satisfactory by low Reynolds model (maybe “pseudo-transition” behavior) are the ones related to bypass transition of which high levels of turbulence in the free stream occur and transition is dominated by diffusion effects (a brief description of the mechanism in my former post: A Forest of Hairpins – on the quest for turbulence coherent structures )
Another major drawback is the model lack of sensitivity to adverse pressure-gradient. It was observed that under such conditions it overestimates the shear stress and by that delays separation. Menter’s k-ω SST alleviates this drawback through the Shear Stress Transport (SST) concept. This drawback is evident in almost all eddy-viscosity models, relating the Reynolds stress to the mean flow strain and in fact is the major difference between such a modeling approach and a full Reynolds-stress model (RSM). The RSM approach accounts for the important effect of the transport of the principal turbulent shear-stress. The ingenious idea of Menter to include it in the revised k-ω model (termed the Baseline (BSL) model) is related observed success in implementing what is termed as the Bradshaw’s assumption, that the shear-stress in the boundary layer is proportional to the turbulent kinetic energy (a brief evaluation of the model in former post: Understanding The k-ω SST Model ).
Remedies in ANSYS Fluent for Near-Wall Treatment
Upon approaching a wall bounded turbulence flow RANS simulation a decision for the level of fidelity by which the turbulent boundary layer needs to be resolved needs to be made. The apparent question at hand is to whether or not resolving the viscous sublayer (y+<5-10) is mandatory to capture the physical phenomena at hand.
Commonly encountered situations where the flow demands integrating through the viscous sublayer include: heat transfer, calculation of aerodynamic drag, flows with adverse pressure gradients, etc’…
Integrating through the viscous sublayer would generally impose resolving the mesh to y+<1.
Typical situations for which the flow does not demand integrating the solution through the viscous sublayer are cases where wall-bounded effects are secondary, or where the flow separates to a dominant shear flow due to geometry such as the case of bluff body surrounding flow. For such cases the first cell away from the wall shall be placed at a y+ of 30-300.
As mentioned above unlike the standard ε-equation, the ω-equation can be integrated through the viscous sublayer without the need for a two-layer approach. This feature can be utilized for a y+-insensitive wall treatment by blending the viscous sublayer formulation and the logarithmic layer formulation based on y+. This formulation is the default for all ω. Such is also the case for the one equation eddy-viscosity transport Spalart-Allmaras turbulence model. While the k-ω model does that by taking an “elliptic” near wall behavior (partial differential wise), meaning that it has an inherent nature of being able to “communicate” with the wall and actually has Dirchlet (as in no-slip in this case) boundary conditions, The Spalart-Allmaras does that by indirectly incorporating the distance from the wall. See “Understanding The k-ω SST Model” and “Understanding The Spalart-Allmaras Turbulence Model” for in-depth review.
As the above also suggests, this is not the case for the k-ε model and traditionally, there are two approaches to modeling the near-wall region. In the first, the viscous sublayer and buffer layer (viscosity-affected inner region) are not resolved. Instead, semi-empirical formulas called “wall functions” are used to bridge the viscosity-affected region between the wall and the fully-turbulent region. The use of wall functions obviates the need to modify the turbulence models to account for the presence of the wall.
In the second approach, the turbulence models are modified to enable the viscosity-affected region to be resolved all the way to the wall, this is termed “Low-Reynolds Modeling”.
The main shortcoming of standard wall functions treatment is that the numerical results deteriorate under refinement of the grid in the wall normal direction. y+ values less than O(10) will gradually result in unbounded errors in wall shear stress and wall heat transfer.
On the other hand, the traditional low Reynolds approach demands for the first grid point away from the wall to be located at y+<1. For cases of which y+>1, the low-Re formulation results in inaccurate wall values for the wall shear stress (and heat transfer). This has is a major drawback on simulations demanding an integration through the viscous sublayer by such a restriction which depicts itself in highly unwanted computational demand that stems from the fact that for most industrial uses, complex geometries are involved.
Moreover, even when a y+<1 is achieved, lets say for a geometry which allows such control of the mesh, the methodology by which it is achieved induces numerical outcomes which are numerically stiff as the ε-equation can’t be integrated through the viscous sublayer to the wall and additional highly non-linear damping functions are needed to construct low-Reynolds formulations (low as in entering the viscous/laminar sublayer) that are able to bridge that gap.
To alleviate the above restriction ANSYS Fluent proposes a few remedies that circumvent the above presented shortcomings.
In what follows, it should be noted that in ANSYS Fluent, the law-of-the-wall for mean velocity and temperature are based on the wall unit, y*, rather than:
An alternative scaling for near wall variables – Why Star?
There is one problem with the scaling we’ve used for the boundary layer. The friction velocity is defined as:
Then when we go to separation, characterized by the wall shear going to zero, the friction velocity is also zero:
and all the scaling breaks down.
Therefore, an alternative scaling is proposed in ANSYS Fluent, by which is based on the turbulent kinetic energy (note that u* and the friction velocity are equal for equilibrium flows):
where we can also build the logarithmic profile based on u*, thus avoiding the singularity in the log profile and other quantities, and compute the wall shear stress:
This means that as the law-of-the-wall for mean velocity yields:
where U* is dimensionless velocity defined as:
The dimensionless distance from the wall, y*, may be defined as:
It is further important to note that in ANSYS Fluent the log-law is employed when y*>11.225. For values below 11.225 a laminar stress-strain relationship is employed (i.e. y*=U*).;
ANSYS Fluent proposes the following near wall treatment approaches:
- Standard wall functions:
The standard wall functions in ANSYS Fluent are based on the Launder and Spalding route.
Wall functions are applicable for a suitable range of y* and by that to the flow’s Reynolds. As the purpose is to allow not integrating through the viscous sublayer, the lower limit of y* is 11 hence standard wall functions should not be used below that limit as the solution’s accuracy might deteriorate in an uncontrolled manner.
The upper limit, as the range of the integral layer in the law-of-the-wall depends on the Reynolds number on a much more physical level. This conceptual understanding of the law-of -the wall is important and should be confronted upon approaching a wall-bounded turbulence problem. For very high Reynolds numbers (such as submarines, ships, airplanes, etc’…) the logarithmic layer can extend to values as high as several thousand, while for low Reynolds applications (such as UAVs, turbine blades, etc’…) the upper limit for the integral layer may be as low as 100. Moreover, for the later, it is much less profound, meaning that the level by which the law-of the-wall and the integral layer could be recognized is unclear.
The last notion implies that the wall functions approach (especially in the standard form) is most suitable for high Reynolds flows of which integrating through the viscous sublayer has only a slight impact on the desired outcome.In k-ε models, the turbulence kinetic energy is solved in the entire domain under the boundary condition zero normal turbulence kinetic energy flux.
The source terms for the kinetic turbulence are its production Gk ,and its dissipation rate, ε. Those are computed by the assumption that they are equal in the near wall control volume – this is termed the local equilibrium hypothesis. Following the above explanation, the production of the turbulence kinetic energy is based on the log-law and may be computed as:and the dissipation equation not solved at the wall-adjacent cells, but instead is computed using:
The standard wall functions work reasonably well for a broad range of wall-bounded flows – the flows inside the range of applicability of constant-shear and local equilibrium assumptions that do not incorporate physical phenomena of interest such as flows involving separation, reattachment, and impingement where the mean flow and turbulence are subjected to pressure gradients and rapid changes.
Scalable wall functions:
Industrial problem often contain complex geometries making the mission of defining consistently the first grid point away from the wall at y*>11 quite a conundrum. To avoid deterioration of the solution by standard wall functions in situations where it’s unavoidable for the first grid point to be located at y*<11, ANSYS Fluent proposes wall functions that produce consistent results for arbitrary grid refinement by forcing the usage of the log law in conjunction with the standard wall functions approach:
- Non-equilibrium wall functions:
As mentioned above standard wall functions will not work well outside their range of applicability (i.e. constant-shear and local equilibrium assumptions).
For flows flows just before separation, just after reattachment, impinging flows, flows of which the mean flow is subjected to adverse pressure gradients and rapid changes, improvements can be obtained, particularly in the prediction of wall shear and heat transfer characteristics. This is achieved in ANSYS Fluent by means of sensitizing the log-law for mean velocity to pressure gradient effects and by the use of the two-layer-based concept to compute the reciprocal relations between turbulence kinetic energy production, Gk, and its dissipation ε, by non-equilibrium means. The log-law mean velocity sensitized mean velocity to pressure gradients is achieved as follows:
the physical viscous sublayer thickness is computed from:
The sensitized mean velocity may be drawn from:
The wall-neighboring cells are assumed to consist of a viscous sublayer and a fully turbulent layer in order to employ the two-layer concept in computing the reciprocal relations between turbulence kinetic energy production, Gk, and its dissipation ε, by non-equilibrium means to solve the turbulence kinetic energy equation at the wall-neighboring cells.
To be able to that the following profile assumptions for turbulence quantities are made:
What is done next is the computation of cell-averaged production of turbulence kinetic energy and cell-averaged dissipation of turbulence kinetic energy from the volume average of Gk and ε of the wall-adjacent cells. These proportions between viscous sublayer and the fully turbulent layer which varies widely from cell to cell in highly non-equilibrium flows actively determine
The turbulence kinetic energy budget for the wall-neighboring cells.
Therefore to some extent the non-equilibrium formulation for the wall functions takes into account the effect of pressure gradients on the distortion of the velocity profiles and by that account for some non-equilibrium effects that are not accounted for in the standard formulation due to the local equilibrium assumption.
Enhanced Wall Treatment ε-Equation (EWT-ε):
If the near-wall mesh is fine enough to be able to resolve the viscous sublayer (typically with the first near-wall node placed at y+<1), then the enhanced wall treatment will be identical to the traditional two-layer zonal model.
However, as noted above the restriction that the near-wall mesh must be sufficiently fine everywhere be limiting from a computational resources standpoint.
The best scenario by concept is to have a near wall treatment suitable for viscous sublayer integration when the mesh is at y+<1 and to use wall functions formulation when the first grid away from the wall falls in the log-layer and the viscous sublayer solution is of minor consequence. Moreover, the scenario should prevent excessive error to occur for the intermediate meshes where the first near-wall node is placed neither in the fully turbulent region, where the wall functions are suitable, nor in the direct vicinity of the wall at y+<1, where the low-Reynold approach is adequate.
This is achieved in ANSYS Fluent by incorporating the two-layer model with enhanced wall functions.
The separation of the two regions is conducted via a wall-distance-based, turbulent Reynolds number:
Where y is the distance to the nearest wall and defined uniquely for complex shaped domains and independent on mesh topology by applying the following on the union of all boundaries involved:
For regions where the above defined Reynolds number is above 200 the original standard k-ε is employed.
By the same token if the above Reynolds number is below 200, a one equation for the transport of turbulence kinetic energy is employed (Wolfstein’s k-equation), while the turbulent viscosity is computed from:
The length scale in the above:
The above intention is two blend the enhanced wall treatment eddy viscosity with that of high Reynolds at the outer region by:
the blending functions is such that it is 1 away from the wall and 0 in their viscinity:
In the viscosity affected region the 1-equation for the turbulence kinetic energy is completed by algebraically computing the dissipation, ε, as:
The length scale in the above (Chen & Patel):
The constants in the above are achieved by calibration to be:
In order for the method’s applicability to be extended to the entire near wall region Fluent also blends the different near wall regions to a single wall-law formulation. One outcomes is that this approach allows the fully turbulent law to be easily modified and extended to take into account other effects such as pressure gradients and variable properties due to compressibility or heat transfer.
Menter-Lechner ε-Equation (ML-ε):
In the above methodology it was shown how the two layer model may serve as a basis for a y+ insensitive model.
We also recall a definition of a wall-distance-based, turbulent Reynolds number to clearly separate two regions:
The use of such a definition hides some shortcomings within. One of its major ones is that regions with small turbulence kinetic energy might also have a turbulent Reynolds number below 200 and as such treated with a near-wall formulation even though they are actually away from the wall.
The Menter-Lechner ε-Equation (ML-ε) approach offers a y+ insensitive which is a not based on a two layer model assuming a sufficient resolution of the boundary layer, predicts in a y+ independent manner the wall shear stress and wall heat flux. Such a formulation switches gradually from wall functions to a low-Reynolds formulation when the mesh is refined.
The Menter-Lechner ε-Equation (ML-ε) concept is achieved by adding a a source term to the transport equation of the turbulence kinetic energy that accounts for near-wall effects:
This source term (marked in red) is active in the viscous sublayer and becomes zero in the logarithmic layer.
This closes the near wall treatment for what generally be referred to as ε-equation models.
To conclude, RANS shall take part in near wall modeling for long time to come. Large-Eddy Simulation (LES) allows high fidelity simulations for many applications but has severe limitations in the near wall regions, as the computational effort required to reliably model the innermost portion of the boundary layer (sometimes constituting more than 90% of the mesh) where turbulence length scale becomes very small is far from the resources available to the industry. Anecdotally, best estimates speculate that a full LES simulation for a complete airborne vehicle at a reasonably high Reynolds number will not be possible until approximately 2050…
Nevertheless even with the enormous reduction in nowadays computing resources limitations encountered with RANS it should be noted that the practitioner should always strive to gain more knowledge on the physical phenomena at hand before clearly preferring one methodology over the other.