International Journal of Astronautics and Aeronautical Engineering
(ISSN: 2631-5009)
Volume 4, Issue 1
Research Article
DOI: 10.35840/2631-5009/7525
Smooth-Wall Boundary Conditions for Energy-Dissipation Turbulence Models#
DF Hunsaker*, WF Phillips and RE Spall
Table of Content
Figures
Figure 1: Solutions to Eq. (52) with...
Solutions to Eq. (52) with randomly selected wall boundary conditions.
Figure 2: Solutions to Eq. (52) with no slip...
Solutions to Eq. (52) with no slip and no dissipation gradient at the wall.
Figure 3: Grid resolution for the mean...
Grid resolution for the mean velocity predicted from the Launder-Sharma k-ε model.
Figure 4: Grid resolution for the turbulent...
Grid resolution for the turbulent energy predicted from the Launder-Sharma k-ε model.
Figure 5: Grid resolution for the turbulent...
Grid resolution for the turbulent dissipation predicted from the Launder-Sharma k-ε model.
Figure 6: Effects of wall boundary conditions...
Effects of wall boundary conditions on turbulent energy predicted from the Lam-Bremhorst model.
Figure 7: Effects of wall boundary conditions...
Effects of wall boundary conditions on turbulent energy predicted from the Launder-Sharma model.
Figure 8: Effects of wall boundary conditions...
Effects of wall boundary conditions on near-wall dissipation for the Launder-Sharma model.
Figure 9: Effects of wall boundary conditions on the...
Effects of wall boundary conditions on the mean velocity predicted from the Launder-Sharma model.
Figure 10: Effects of wall boundary conditions...
Effects of wall boundary conditions on turbulent energy predicted from the Wilcox 1998 k-ω model.
Figure 11: Effects of wall boundary conditions...
Effects of wall boundary conditions on the turbulent dissipation frequency predicted from the Wilcox 1998 k-ω model.
References
- J Boussinesq (1877) Théore de l'Écoulement Tourbillant, Mém. Présentés par Divers Savants Acad. Sic Inst Fr 23: 46-50.
- WP Jones, BE Launder (1972) The prediction of laminarization with a two-equation model of turbulence. International Journal of Heat and Mass Transfer 15: 301-314.
- CKG Lam, K Bremhorst (1981) A modified Form of the k-ε Model for Predicting Wall Turbulence. Journal of Fluids Eng 103: 456-460.
- BE Launder, BI Sharma (1974) Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc. Letters in Heat and Mass Transfer 1: 131-138.
- VC Patel, W Rodi, G Scheuerer (1985) Turbulence models for near-wall and low reynolds number flows: A reveiw. AIAA Journal 23: 1308-1319.
- PA Durbin (1991) Near-Wall turbulence closure modeling without damping functions. Theoretical and Computational Fluid Dynamics 3: 1-13.
- PA Durbin (1993) Application of a near-wall turbulence model to boundary layers and heat transfer. International Journal of Heat and Fluid Flow 14: 316-323.
- (2006) Standard k-ω Model. FLUENT 6.3 User's Guide. Fluent Inc, Lebanon, 12-31.
- SJ Kline, BJ Cantwell, GM Lilley (1981) Proceedings of the 1980-1981 AFOSR-HTTM-Stanford Conference on Complex Turbulent Flows: Comparison of Computation and Experiments 1-3 (1981) Thermoscience Division, Dept. of Mechanical Engineering, Stanford University.
- HK Myong, N Kasagi, M Hirata (1989) Numerical prediction of turbulent pipe flow heat transfer for various prandtl number fluids with the improved k-ε turbulence model. JSME International Journal Series 32: 613-622.
- HK Myong, N Kasagi (1990) A new approach to the improvement of k-ε Turbulence model for wall-bounded shear flows. JSME International Journal 33: 63-72.
- K Abe, T Kondoh, Y Nagano (1994) A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows - I. flow field calculations. International Journal of Heat and Mass Transfer 37: 139-151.
- DC Wilcox (2008) Formulation of the k-ω turbulence model revisited. AIAA Journal 46: 2823-2838.
- AN Kolmogorov (1942) Equations of turbulent motion of an incompressible fluid. Investia Academy of Sciences, USSR, Physics 6: 56-58.
- DC Wilcox (1998) One-equation and two-equation models, turbulence modeling for CFD. (2nd edn), DCW Industries, Inc, La Canada CA, 103-226.
- DC Wilcox (1988) Numerical accuracy near boundaries, Turbulence modeling for CFD. (2nd edn), DCW Industries, Inc, La Canada CA, 341-344.
- PG Saffman (1970) A model for inhomogeneous turbulent flow. Proceedings of the Royal Society of London A 317: 417-433.
- DC Wilcox (2006) One-equation and two-equation models, turbulence modeling for CFD. (3rd edn), DCW Industries, Inc, La Canada, 107-237.
Author Details
DF Hunsaker*, WF Phillips and RE Spall
Department of Mechanical and Aerospace Engineering, Utah State University, USA
Corresponding author
DF Hunsaker, Assistant Professor, Department of Mechanical and Aerospace Engineering, Utah State University, 4130 Old Main Hill, Logan, Utah 84322-4130, USA, Tel: 435-797-8404.
Accepted: February 26, 2019 | Published Online: February 28, 2019
Citation: Hunsaker DF, Phillips WF, Spall RE (2019) Smooth-Wall Boundary Conditions for Energy-Dissipation Turbulence Models#. Int J Astronaut Aeronautical Eng 4:025.
Copyright: © 2019 Hunsaker DF, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Abstract
It is shown that the smooth-wall boundary conditions specified for commonly used dissipation-based turbulence models are mathematically incorrect. It is demonstrated that when these traditional wall boundary conditions are used, the resulting formulations allow either an infinite number of solutions or no solution. Furthermore, these solutions do not enforce energy conservation and they do not properly enforce the no-slip condition at a smooth surface. This is true for all dissipation-based turbulence models, including the k-ε, k-ω, and k-ζ models. Physically correct wall boundary conditions must force both k and its gradient to zero at a smooth wall. Enforcing these two boundary conditions on k is sufficient to determine a unique solution to the coupled system of differential transport equations. There is no need to impose any wall boundary condition on ε, ω, or ζ at a smooth surface and it is incorrect to do so. The behavior of ε, ω, or ζ approaching a smooth surface is that required to satisfy the differential equations and force both k and its gradient to zero at the wall.
Nomenclature
C1-C5: Arbitrary Constants of Integration, Eq. (43); , : Turbulence Model Closure Coefficients, Eqs. (6) and (7); , : Turbulence Model Closure Coefficients, Eqs. (73) and (74); : Turbulence Model Closure Coefficient, Eqs. (2) and (7); f1, f2: Wall Damping Functions, Eq. (9) or (77); : Wall Damping Function, Eq. (76); : Wall Damping Function, Eq. (8) or (75); E: Wall Damping Function, Eq. (9); E+: Wall-Scaled Dimensionless Wall Damping Function, ; : Jacobian Tensor for a Vector Field; k: Turbulent Kinetic Energy Per Unit Mass, Eq. (1); k+: Wall-Scaled Dimensionless Turbulent Kinetic Energy, Eq. (32); : Arbitrary Dependent Variable, Eqs. (42) and (88); l: Channel Half Width; l+: Wall-Scaled Dimensionless Channel Half Width, ; : Pseudo Mean Pressure, , where is the Mean Pressure, go is the Standard Acceleration of Gravity at Sea Level, , H is Geometric Altitude, and RE is the Radius of the Earth; p+: Wall-Scaled Dimensionless Pseudo Mean Pressure Gradient, Eq. (32); q+: Change of Variables, Eq. (51); Rt: Turbulent Dissipation Reynolds Number, Eq. (35) or (81); Ry: Turbulent Wall Reynolds Number, Eq. (35); : Strain-Rate Tensor for a Vector Field; u+: Wall-Scaled Dimensionless x-Velocity Component, Eq. (32); : Arbitrary Dependent Variable, Eqs. (42) and (88); uT: Friction Velocity, Eq. (28); : Mean Velocity Vector; : Fluctuating Velocity Vector; : Magnitude of the Fluctuating Velocity Vector; : x Component of Mean Velocity Vector; : y Component of Mean Velocity Vector; : x Component of Fluctuating Velocity Vector; x: Axial Coordinate; y: Normal Coordinate Measured Outward from a Wall; y+: Wall-Scaled dimensionless y coordinate, Eq. (32); ε: Turbulent energy-dissipation parameter, Eq. (1); : Turbulent Energy-Dissipation Parameter, ; : Wall-Scaled Dimensionless Turbulent Dissipation Parameter, Eq. (32); : Arbitrary Dependent Variable, Eq. (42); : Wall damping function, Eq. (10); : Wall-Scaled Dimensionless Wall Damping Function, Eq. (32); : Turbulent Energy-Dissipation Parameter, ; : Change of Variables, Eq. (51); : Kinematic Molecular Viscosity; vt: Kinematic Eddy Viscosity; v+: Ratio of the Turbulent Eddy Viscosity to the Molecular Viscosity, Eq. (33); : Fluid Density; : Turbulence Model Closure Coefficient, Eqs. (5) and (7); : Turbulence model closure coefficient, Eqs. (6) and (7); : Wall Shear Stress, Eq. (28); : Turbulent Energy-Dissipation Frequency, ; : Wall-Scaled Dimensionless Turbulent Dissipation Frequency, Eq. (79); : Arbitrary Dependent Variable, Eq. (88).
Introduction
Many of the turbulence models that are now commonly used for computational fluid dynamics (CFD) are based on the analogy between molecular and turbulent transport that was first proposed by Boussinesq [1]. The majority of these turbulence models are usually classified as either k-ε, k-ω, or k-ζ models. Conventional k-ε, k-ω, and k-ζ turbulence models are often thought of as being fundamentally different. Yet, in a larger sense, these three model classifications could all be thought of as energy-dissipation models. This is because all such models are based on the hypothesis that Boussinesq's eddy viscosity is proportional to the product of the root mean square fluctuating velocity, or k1/2 and the dissipation length scale . The parameters k and ε are defined in terms of the fluctuating velocity as
where is the fluctuating velocity vector, is its Jacobian tensor, and the over score denotes an ensemble mean. The eddy-viscosity model that is the foundation for all commonly used k-ε, k-ω, and k-ζ turbulence models is
where is a dimensionless closure coefficient that is nearly universally accepted as being equal to 0.09. The k-ε turbulence models use Eq. (2) directly. The k-ω turbulence models use the change of variables to transform Eq. (2) to the equivalent relation given by . Similarly, the k-ζ turbulence models use the change of variables to transform Eq. (2) to its k- ζ equivalent, . The commonly used k-ε, k-ω, and k-ζ turbulence models are all based on the hypothesis that the characteristic length scale for turbulent transport is proportional to the characteristic length scale for turbulent energy dissipation.
The k-Ɛ turbulence model that is the foundation for most modern Boussinesq-based turbulence models is that of Jones and Launder [2]. In addition to the algebraic equation for the kinematic eddy viscosity that is given by Eq. (2), the Jones-Launder turbulence model comprises the following equations for steady incompressible flow. The continuity equation,
the Boussinesq-based Reynolds-averaged-Navier-Stokes (RANS) equations,
the Boussinesq-based turbulent-energy-transport equation
and a turbulent-dissipation-transport equation obtained by analogy with Eq. (5)
The commonly used closure coefficients for this model are
In this form, the Jones-Launder k-ε turbulence model does not exhibit the proper behavior near a solid wall. Near a no-slip boundary the turbulent velocity fluctuations and turbulent transport are suppressed by the proximity of the solid surface. Accurately modeling this suppression is critical to obtaining accurate predictions for the wall shear stress and heat transfer.
In the attempt to provide realistic results near a wall, the Jones-Launder k-ε model is often implemented with the incorporation of what are called wall damping functions. In a general form, these wall damping functions are added to Eq. (2), Eq. (6), and the definition of ε,
A variety of k-ε turbulence models have been proposed, which differ only in the form of the wall damping functions , f1, f2, E, and ε0. To complete any k-ε model of this form, the wall damping functions are specified as prescribed functions of v, k, and the normal coordinate y, measured outward from the wall. These wall damping functions are simply empirical corrections that are added to force the model to agree more closely with experimental data.
For steady, incompressible, 2-D flow in Cartesian coordinates, the k-ε turbulence model with wall damping functions is specified by
If y is the normal coordinate measured outward from a smooth wall, then the traditional no-slip wall boundary conditions for and are
Likewise, the obvious no-slip wall boundary condition for k at a smooth wall is
The remaining wall boundary condition is assumed to be model dependent and it is the topic of the present paper. The mean velocity and turbulent fluctuations vanish at all points on a smooth wall. Hence, near a smooth wall, changes in the mean velocity and turbulence variables with respect to x are small compared to changes with respect to y, and the near-wall formulation reduces to
Integrating Eq. (19) and applying Eq. (17) yields
Using Eq. (24) in Eq. (21) yields
or
Using Eqs. (24) and (25) in Eq. (20) gives
Because the right-hand side of Eq. (26) is only a function of x, integrating Eq. (26) from the wall to some point y that is still near the wall results in
Because the turbulent eddy viscosity is zero at a smooth wall, the left-hand side of Eq. (27) evaluated at the wall can be written in terms of either the wall shear stress or the friction velocity,
After using Eq. (28) in Eq. (27), the near-wall approximation for the Boussinesq-RANS equations becomes
Similarly, after applying Eq. (24) to Eqs. (22) and (23), the near-wall approximations for the k- and ε-transport equations become
The near-wall formulation given by Eqs. (29)-(31) is commonly called the parallel-flow approximation. This simplification was obtained using only the approximation that changes in the mean velocity and turbulence variables with respect to x are negligible compared to changes with respect to y. This is approximately true for any flow in the region very close to a solid surface. Furthermore, the simplifications used in obtaining Eqs. (29)-(31) hold exactly for fully developed flow in channels.
For attached flows, it is convenient to nondimensionalize the differential equations in Eqs. (29)-(31) using the traditional wall-scaled dimensionless variables as a similarity transformation
Although it is not standard convention, here we shall denote the ratio of the turbulent eddy viscosity to the molecular viscosity as v+. Thus, from Eq. (11)
Applying Eqs. (32) and (33) to Eqs. (25) and (29)-(31) the near-wall formulation is
where . To complete the formulation, the damping functions , f1, f2, E+, , and six boundary conditions must be specified for this coupled sixth-order system.
Although several variations for the k-ε wall damping functions have been proposed, none have been completely successful. In the present paper we shall examine two of the commonly used models. The first k-ε turbulence model to be considered is that developed by Lam and Bremhorst [3]. This is typical of models that use . The second k-ε model to be considered here is that developed by Launder and Sharma [4]. This model is typical of those that use prescribed functions for E+ and , which are nonzero.
The Lam-Bremhorst k-ε Model
For incompressible flow, the Lam-Bremhorst k-ε turbulence model [3] uses the wall damping functions defined by
Hence, we see that for this k-ε model the wall damping functions can be written in terms of only the traditional wall-scaled dimensionless variables.
There has not been universal agreement regarding appropriate wall boundary conditions for the Lam-Bremhorst k-ε model. As a specific example of how such boundary conditions affect the flow, we shall consider the case of fully developed flow in a 2-D channel of half width l. The parallel-flow simplifications used in obtaining Eq. (34) hold exactly for this fully developed flow. Because y was defined to be the normal coordinate measured outward from the solid wall, the symmetry boundary conditions at the channel centerline require
where . The first of the four differential equations in Eq. (34) can be integrated analytically and the first of the three boundary conditions in Eq. (36) can be used with the second differential equation in Eq. (34) to evaluate p+, which gives
Two of the three boundary conditions required at the wall are the obvious no-slip conditions in Eqs. (17) and (18). Hence, using these two boundary conditions and the two remaining symmetry boundary conditions from Eq. (36), Eqs. (34)-(37) reduce to the incomplete one-dimensional fifth-order formulation.
where the ' indicates a derivative with respect to y+.
One additional boundary condition is needed to complete the formulation given by Eq. (38). In their original publication, Lam and Bremhorst [3] give the remaining boundary condition for Eq. (38) as
or (39)
where the " indicates a second derivative with respect to y+. Although this boundary condition is commonly accepted in the literature as being the appropriate smooth-wall boundary condition for Eq. (38), it is in fact mathe-matically incorrect. There are an infinite number of solutions to Eq. (38), which also satisfy Eq. (39).
To see why Eq. (39) is not a viable boundary condition for Eq. (38), consider Eq. (38) in the limit as y+ approaches zero. In this limit, both Rt and Ry go to zero and the wall damping functions and eddy viscosity near a smooth wall reduce to
Using these limiting relations in the differential equations from Eq. (38) produces the near-wall system of equations, which applies in the limit as y+ approaches zero,
From the development of Eq. (41) it can be seen that Eq. (38) satisfies Eq. (39), independent of the fifth boundary condition. Hence, Eq. (39) does not provide the additional information required to obtain a unique solution to the indeterminate system in Eq. (38).
As a further demonstration of why Eq. (39) is not a viable boundary condition for completing the indeterminate system in Eq. (38), consider the similar system
This indeterminate fifth-order system of differential equations with only four boundary conditions is mathematically similar to Eq. (38), yet it is simple enough to permit obtaining a closed-form solution by direct integration. The general solution to Eq. (42) is
After applying the four boundary conditions given in Eq. (42) to eliminate four of the five arbitrary constants, the solution in Eq. (43) yields
As should be expected, there are an infinite number of solutions to any indeterminate fifth-order system of differential equations with only four boundary conditions, such as that specified in either Eq. (38) or Eq. (42). However, if the mathematical logic presented by Lam and Bremhorst [3] is correct, then we should be able to reduce Eq. (44) to a single unique solution by simply applying a boundary condition obtained from the second differential equation in Eq. (42) evaluated at y=0. If we accept this logic, then our final boundary condition for Eq. (42) is
However, the reader may not be surprised to learn that applying Eq. (45) to either Eq. (43) or Eq. (44) yields only the trivial result C4 = C4. Hence, there are an infinite number of solutions to Eq. (42) that also satisfy Eq. (45).
From the discussion presented here, it should be clear that, in developing Eq. (39) as a boundary condition for Eq. (38), Lam and Bremhorst [3] used mathematical logic that is seriously flawed. Equation (39) is certainly a valid near-wall asymptote for the k-transport equation in Eq. (38). Thus, Eq. (39) can be used as an alternative to the k-transport equation for y+ approaching zero, provided that it is combined with five appropriate boundary conditions. However, Eq. (39) cannot be used in lieu of one of the five boundary conditions. If valid boundary conditions could be obtained directly from the differential equations to which they apply, separate boundary conditions would not be needed to isolate a unique solution from a general solution.
If Eq. (39) is not a mathematically viable boundary condition for Eq. (38), it may be fair for a reader to ask, "How is it possible that numerical solutions to Eq. (38) have been obtained using Eq. (39) as a boundary condition?" To answer this question, we must remember that traditional CFD algorithms provide no information regarding the uniqueness of any solutions found. If multiple solutions exist, the numerical algorithm may converge on one of these solutions, but we have no assurance that the solution satisfies the physically correct wall boundary conditions, unless those boundary conditions have all been numerically enforced. It is always the user's responsibility to specify the boundary conditions appropriately, so that any solution found will be unique and physically correct.
In a review of early turbulence models, Patel, Rodi, and Scheuerer [5] point out that using Eq. (39) as a boundary condition for Eq. (38) "is not very convenient since it involves parts of the solution of the system of coupled differential equations." Although Patel, Rodi, and Scheuerer [5] did not state that the boundary condition proposed by Lam and Bremhorst [3] is incorrect, they did suggest a "more convenient boundary condition,"
or (46)
which is also incorrect.
Durbin [6] was the first to point out that, as boundary conditions for Eq. (38), both Eq. (39) and Eq. (46) are incorrect. Durbin [6] presented the correct smooth-wall boundary conditions for Eq. (38), which are
or (47)
With reference to the wall boundary conditions specified by Eqs. (39) and (46), Durbin [6] states, "These conditions must violate the energy balance, and do not ensure satisfaction of conditions (47)." With reference to the k conditions specified in Eq. (47), in a later publication Durbin [7] states, "These two conditions on k suffice to determine the solution for the coupled system of equations; there is no need to impose conditions of ε at the wall - indeed, it would be incorrect to do so."
To demonstrate why the smooth-wall boundary conditions given by Durbin [6] are correct, consider the implications of the no-slip boundary condition on the turbulent velocity fluctuations. At a smooth surface, both the mean and fluctuating velocity components must vanish. The definitions of k and ε are given by Eq. (1). Clearly, the no-slip boundary condition applied to Eq. (1) requires Eq. (18). However, because the no-slip boundary condition places no restriction on the derivative of with respect to y, and ε depends only on the Jacobian of , physics imposes no boundary condition on ε. The second wall boundary condition required for the coupled fourth-order system of k-ε transport equations is obtained by taking the gradient of k, which from the definition in Eq. (1) gives
up>which results in and . On the other hand, the boundary condition for Eq. (42) that is analogous to Eq. (46) is . Applying this boundary condition to the solution of Eq. (42) that is given in Eq. (44) yields the somewhat more troubling result , which may inspire some concern with regard to using Eq. (46) as a boundary condition for Eq. (38).
Examination of the incomplete fifth-order system given by Eq. (42) has revealed that using as the fifth boundary condition results in an infinite number of solutions. On the other hand, Eq. (42) has no solution if is used as the fifth boundary condition. It can be shown that the incomplete fifth-order system in Eq. (38) exhibits very similar behavior. However, solutions to Eq. (38) must be obtained numerically.
Because fully developed flow is one dimensional, a solution to Eq. (38) combined with Eq. (49) can be obtained by direct numerical integration. This permits the use of efficient high-order numerical methods such as the fourth-order Runge-Kutta algorithm. Because such solutions can be obtained quickly on very fine grids, fully developed channel flow provides an excellent benchmark for testing more computationally intensive CFD algorithms.
To facilitate direct numerical integration, the two second-order equations in Eq. (38) can be converted to four first-order equations by using the change of variables
Combining Eq. (38) with Eq. (49), applying the change of variables given in Eq. (51), and eliminating by direct substitution provides the complete one-dimensional fifth-order formulation
It should be noted that the new variable q+ is a dimensionless form of the total diffusive flux of turbulent kinetic energy k, which includes both molecular and turbulent diffusion. Similarly, ?+ is a dimensionless diffusive flux for ε. This brings to light another physical interpretation of the boundary condition given in Eq. (49), which led directly to the equivalent boundary condition in Eq. (52), i.e., . With this interpretation, Eq. (49) can be viewed as a mathematical statement of the simple fact that turbulent kinetic energy cannot be diffused through a solid wall. The formulation for fully developed flow that is given by Eq. (52) requires that vanish at the wall and at the centerline. Thus, all of the turbulent kinetic energy that is generated within this steady flow must also be dissipated within the flow. If a boundary condition obtained from either Eq. (39) or Eq. (46) is used in place of that obtained from Eq. (49), this energy balance is not enforced. This is the origin of Durbin's statement that, "These conditions must violate the energy balance," [6].
A numerical solution to the five first-order differential equations given in Eq. (52) can be obtained using fourth-order Runge-Kutta integration combined with an appropriate numerical root-finding method. Because only three of the five boundary conditions are given at y+ = 0, the solution for and must be obtained from the differential equations. The process is started with initial estimates for and . From these initial estimates, fourth-order Runge-Kutta integration can be used to obtain and . The initial estimates are then refined using an appropriate numerical method until the solution is found, which corresponds to the correct centerline values and .
A few words of caution may be in order here. Some of the terms in Eq. (52) are numerically indeterminate if K+=0 and/or . Notice that a division by zero occurs in the definition of Rt for . Thus, depending on the compiler, conditional relations may be required to enforce and f2=1 for . For most compilers, Eq. (52) is numerically indeterminate for k+=0. In this limit, both Rt and Ry go to zero and the eddy viscosity and wall damping functions reduce to
Hence, for the limit , the formulation given in Eq. (52) should be conditionally replaced with its near-wall asymptote
To demonstrate that Eq. (39) is not a valid boundary condition for completing the formulation given in Eq. (38), the results shown in Figure 1 were obtained from Eq. (52) using randomly selected wall boundary conditions. The no-slip wall boundary conditions were used for both u+ and k+, but the wall boundary conditions for q+, and as well as the wall-scaled dimensionless half width l+ were generated as listed in Figure 1 using the "rand" function, which generates a random number between 0.0 and 1.0. From the results presented in Figure 1, it can be concluded that Eq. (39) is enforced directly by Eq. (38), completely independent of the boundary conditions.
To demonstrate that Eq. (46) is not a valid boundary condition for completing the formulation given in Eq. (38), the results shown in Figure 2 were obtained from Eq. (52) using the no-slip wall boundary conditions for and with the wall boundary conditions for obtained from Eq. (46). For several values of l+ the computed value for q+ at the centerline is plotted as a function of the remaining wall boundary condition . Valid solutions to Eq. (38) could only correspond to those points where these curves intersect the axis . From the results presented in Figure 2, it can be seen that there is only one solution to Eq. (38) that satisfies Eq. (46) and the no-slip wall boundary conditions. That is the trivial laminar solution
There is no turbulent flow solution to Eq. (38) that satisfies Eq. (46) and the no-slip wall boundary conditions.
Examination of the numerical results shown in Figure 1 and Figure 2 reveals that Eq. (38) exhibits behavior very similar to that demonstrated analytically for the hypothetical fifth-order system given by Eq. (42). Using as the fifth boundary condition for Eq. (38) results in an infinite number of solutions. On the other hand, Eq. (38) has no valid turbulent flow solution if is used as the fifth boundary condition. This underscores the critical importance of always using the correct no-slip boundary conditions .
The Launder-Sharma k-ε Model
For incompressible flow, the wall damping functions that are used in Eqs. (8)-(10) are defined for the Launder-Sharma k-ε turbulence model [4] as
It may be worth noting that Launder and Sharma [4] originally defined the wall damping function ε0 in a slightly different but equivalent form,
Following what was done with the Lam-Bremhorst model, the Launder-Sharma wall damping functions can be written in terms of the wall-scaled dimensionless variables that are defined in Eqs. (32) and (33). The result applied to Eq. (34) yields the near-wall formulation for the Launder-Sharma k-ε model
From the first two differential equations in Eq. (57) and the specified relation for v+, it is easily shown that the last term on the right-hand side of the ε-transport equation can be evaluated from
The wall boundary conditions for this model as specified originally by Launder and Sharma [4] are
or (59)
To examine the near-wall behavior of the Launder-Sharma k-ε model, consider the Taylor-series expansions
If we impose the wall boundary conditions and leave as an unknown constant, the near-wall expansions for the damping functions and eddy viscosity are
Using these expansions in Eq. (57), the near-wall expansion for the k-transport equation yields
which is singular at the wall. Notice that enforcing does not enforce . On the other hand, if we apply the wall boundary conditions and treat as an unknown, we obtain
and the near-wall expansion for the k-transport equation becomes
which requires and . Hence, we see that referring to the relation as a boundary condition for the Launder-Sharma k-ε model is a misnomer. If the correct no-slip boundary conditions specified in Eq. (47) are applied with the Launder-Sharma k-ε model, then follows directly from the k-transport equation. In any case, the correct no-slip boundary condition should be enforced with the Launder-Sharma k-ε model, because the k-transport equation is singular at the wall for nonzero values of . With the correct no-slip boundary conditions enforced, the near-wall expansions for the Launder-Sharma damping functions and eddy viscosity are
From the development of Eq. (38), fully developed flow in a 2-D channel of half width l requires and the Launder-Sharma model combined with the boundary conditions given in Eqs. (36) and (47) yields
Using the change of variables defined in Eq. (51) combined with Eq. (58) yields the complete one-dimensional fifth-order formulation
For the limit , the numerically indeterminate terms in Eq. (66) should be conditionally replaced with their near-wall asymptotes
Note that in the limit , the differential equation for q+ reduces to the algebraic equation . Moreover, this differential equation is satisfied for any value of . Thus, is an unknown constant that can be varied along with to enforce the two centerline boundary conditions. If the no-slip boundary condition is not enforced, the differential formulation is mathematically indeterminate.
It should be emphasized that, because physics imposes two wall boundary conditions on k and none on ε, the value of ε throughout the flow field, including its value at the wall, must be determined exclusively from the differential transport equations while enforcing the two wall boundary conditions on k. Because the relation follows directly from the differential equations as shown in Eq. (63), this is certainly a valid relation that can be used to replace a differential transport equation at the wall, provided that it is combined with a complete set of appropriate boundary conditions. However, it cannot be used in lieu of one of the required boundary conditions, as was proposed originally by Launder and Sharma [4] in their presentation of this classical turbulence model.
Numerical Results from CFD Algorithms
Because the zero-gradient boundary condition for k in Eq. (47) is not explicitly enforced in many commonly implemented k-ε turbulence models, solutions obtained from these models are not unique. To demonstrate this fact, the RANS formulations for fully developed channel flow presented in Eqs. (38) and (65) were solved numerically using a second-order finite difference algorithm with successive underrelaxation. Solutions were obtained on the domain extending from the wall to the channel centerline, and grid points were clustered near the wall using logarithmic clustering. To ensure that all results were fully converged, the successive underrelaxation was allowed to continue until the observed changes were reduced to within the double-precision machine accuracy.
To ensure that all results were grid resolved, the grids were uniformly refined until no significant changes were observed with additional grid refinement. For a given axial pressure gradient, the Launder-Sharma model required a somewhat finer grid than was required for the Lam-Bremhorst model. Results of an example grid-resolution study for the Launder-Sharma model are shown in Figure 3, Figure 4 and Figure 5. All results shown in these figures were obtained using the fixed axial pressure gradient, which yields a value of at the centerline equal to 300. For the grid refinements shown in Figure 3, Figure 4 and Figure 5, the four grids produced channel Reynolds numbers (based on the channel width and mean velocity) that were equal to 10,009, 10,653, 10,832, and 10,878, respectively. An additional refinement of the grid to 401 nodes, which is not shown in Figure 3, produced a channel Reynolds number of 10,889. From these and other similar results, it was concluded that for Reynolds numbers on the order of 10,000, the 201-node grid used for Figure 3, Figure 4 and Figure 5 produced adequate grid resolution with both the Lam-Bremhorst and Launder-Sharma turbulence models.
To demonstrate that solutions obtained from commonly implemented k-ε turbulence models are not unique, the second-order successive underrelaxation algorithm was implemented using a slight variant of the wall boundary conditions specified in Eq. (47), which allows the user to specify an arbitrary value for the gradient of k at the wall. Figure 6 and Figure 7 show converged and grid-resolved results obtained from this algorithm using three different gradient boundary conditions for k at the wall: , and .
To demonstrate that traditional implementations of these turbulence models do not necessarily converge to the solution that yields , Figure 6 and Figure 7 also include results obtained from the same algorithm, turbulence models, and grid, but with a traditional implementation of the wall boundary conditions, which uses only and together with an asymptotic relation obtained from the differential equations, i.e., for the Lam-Bremhorst model and for the Launder-Sharma model. For additional comparison, Figure 6 and Figure 7 also show results from the general-purpose finite-volume CFD solver FLUENT [8], which were obtained using the same turbulence models and grid with only the traditional wall boundary conditions implemented.
The results shown in Figure 6 and Figure 7 clearly demonstrate that when the natural boundary condition is not enforced, solutions obtained from commonly implemented k-ε turbulence models are not unique. When the boundary condition is omitted, solutions obtained from the resulting indeterminate formulation are implementation dependent. Notice from Figure 6 that for the Lam-Bremhorst model with traditional implementa-tion of the wall boundary conditions, the finite-difference algorithm converged to a different solution from that obtained using the finite-volume algorithm with the same indeterminate boundary conditions. Neither of these solutions agrees with that obtained from the finite-difference algorithm with the boundary condition enforced. Similarly, we see from Figure 7 that these finite-difference and finite-volume implementations of the Launder-Sharma model converge to different solutions with traditional implementations of the wall boundary conditions. However, for the particular implementation used to obtain the results shown in Figure 7, the indeterminate finite-difference algorithm converged to a solution that is very close to that obtained when the complete set of smooth-wall boundary conditions was enforced. This should not be viewed as an endorsement for implementing the Launder-Sharma turbulence model with mathematically incomplete boundary conditions.
From the near-wall expansions of the Launder-Sharma model given in Eqs. (62) and (63), it was shown that enforcing requires , whereas enforcing does not require . This can also be demonstrated numerically by examining the near-wall behavior of ε obtained from numerical solutions using different gradient boundary conditions for k at the wall. Such results are shown in Figure 8, which were obtained from converged and grid-resolved solutions for the same channel flow that was used to obtain the results shown in Figure 7. Notice that, although does affect the near-wall behavior of all of these solutions satisfy . As can be seen from Eqs. (62), (63), and the discussion here, only the solution corresponding to also satisfies the physically correct no-slip condition at the smooth wall.
Because the wall damping functions for the Lam-Bremhorst model decay rapidly with increasing y+, the wall boundary condition has little impact on the velocity profiles predicted from this turbulence model. On the other hand, the wall damping functions for the Launder-Sharma model decay slower and have a more significant effect on the predicted mean velocity farther from the wall. This can be seen in Figure 9, which displays the velocity profiles for the same solutions that were used to obtain the turbulent energy profiles displayed in Figure 7. It may be worth reiterating that all of the solutions shown in Figure 9 satisfy the traditional wall boundary conditions and together with the asymptotic relation obtained from the differential equations, .
Anyone who has taken time to compare results obtained from different CFD algorithms and different k-ε turbulence models will likely have noticed that there is often a greater difference between the results obtained from two different implementations of the same turbulence model than there is between the results obtained from the same implementation of two different turbulence models [9]. The results shown in Figure 9 may shed some light on the reason for this phenomenon. We should not be too surprised to learn that results obtained from commonly used k-ε turbulence models are implementation dependent, if we recognize that these models are short one boundary condition, and thus are mathematically indeterminate
Because the CFD community has not traditionally implemented two wall boundary conditions on k and none on ε, implementation of the correct smooth-wall boundary conditions first proposed by Durbin [6] has been less than enthusiastically embraced. The actual implementation of these boundary conditions is dependent on the numerical method being used to solve the system of differential equations. However, this implementation should not be difficult using well-known methods in either finite-difference or finite-volume algorithms. For example, results presented in this section were obtained from a finite-difference algorithm. To implement the no-slip wall boundary conditions, the k-transport equation at any wall node was replaced with the boundary condition . At the first node off the wall, the k-transport equation was replaced with a second-order finite-difference approximation for the boundary condition . Because there is no wall boundary condition on ε, the ε-transport equation at any wall node was replaced with the asymptotic relation obtained from the differential equations, i.e., for the Lam-Bremhorst model and for the Launder-Sharma model. The implementation of the ε-transport equation is identical to that of the traditional formulation. The error in the traditional formulation is not that the transport equation for ε is incorrectly implemented. Rather, the error is in assuming that the near-wall asymptote obtained from the differential equations can be used to replace the final boundary condition required at the wall.
As scientists and engineers, we do not have the luxury of choosing boundary conditions for ease of numerical implementation. Boundary conditions are dictated by physics. It is our obligation to understand and implement them correctly if we hope to achieve mathematical formulations that correctly model physics. The fundamental mathematical error of deriving a so called boundary condition directly from the differential equations is not unique to the classical k-ε turbulence models that have been considered here. It is also an important concern for many other turbulence models developed more recently [10-13].
Application to k-ω Models
The k-ω turbulence models that are commonly used for CFD are built on exactly the same dissipation-based eddy-viscosity model that is given in Eq. (2), where k and ε are defined by Eq. (1). These commonly used k-ω turbulence models are based on applying a simple change of variables to Eq. (2), i.e.,
This change of variables applied to Eq. (2) yields an algebraic equation for the kinematic eddy viscosity in terms of only the turbulent kinetic energy per unit mass, k, and the turbulent energy-dissipation frequency, ω,
In addition to this algebraic equation for the kinematic eddy viscosity, the k-ω turbulence model originally proposed by Kolmogorov [14] has been refined to comprise the following equations for steady incompressible flow. The continuity equation combined with the Boussinesq-RANS equations,
the Boussinesq-based turbulent-energy-transport equation obtained by applying the change of variables defined in Eq. (68) to Eq. (5),
and a dissipation-frequency-transport equation obtained by analogy with Eq. (72),
The closure coefficients differ slightly from one version of the model to another and have changed as the model has evolved over the past six decades. In the original k-ω model, Kolmogorov [14] assumed = 0 and he did not include the molecular diffusion term. The closure coefficients often used for the k-ω model [8,15] are
It should be noted that the turbulence variable ω, which is defined by Eq. (68) and referred to here as the turbulent energy-dissipation frequency, is often referred to as the specific dissipation rate
As is the case for the k-ε model, the standard k-ω model does not exhibit the proper behavior near a solid wall. By direct analogy with what has been done with the k-ε model, the k-ω model could also be implemented with the incorporation of wall damping functions. Although this terminology is not commonly used with the k-ω model, to emphasize similarities between the low-Reynolds-number corrections used for the k-ω model and those used for the k-ε model, here we will use exactly the same notation and terminology for both models. Adding wall damping functions to Eqs. (69), (72), and (73) yields
To complete any k-ω turbulence model in this form, the wall damping functions , , f1, and f2, could be specified as prescribed functions of v, k, and ω. As is the case for the k-ε model, these wall damping functions are simply empirical corrections, which are employed to force the model to agree more closely with near-wall experimental data.
Following the development of Eq. (34) for steady, incompressible, 2-D flow in Cartesian coordinates, the near-wall approximation for the k-ω turbulence model with wall damping functions can be written as
Continuing to follow what was done with the k-ε model, the differential equations in this formulation can be nondimensionalized using the wall-scaled dimensionless variables defined in Eqs. (32) and (33) together with the traditional definition for ω+
The result yields the dimensionless near-wall k-ω formulation for steady 2-D incompressible flow
As an example of a k-ω turbulence model that includes wall damping functions, consider what is commonly called the Wilcox 1998 k-ω model [15], which is implemented in FLUENT [8]. Although Wilcox [15] uses a different notation, his formulation is easily rearranged to the format of Eq. (78). For 2-D incompressible flow the resulting wall damping functions are
Notice that the turbulent dissipation Reynolds number Rt, as it is defined in the Wilcox 1998 k-ω model, differs from that defined for the k-ε models, i.e.,
Following the development of Eq. (38), these wall damping functions can be written in terms of the wall-scaled dimensionless variables that are defined in Eqs. (32) and (79). Hence, for fully developed flow in a 2-D channel of half width l, the Wilcox 1998 k-ω model combined with the boundary conditions given in Eqs. (17), (18), and (36) yields
As in the case of Eq. (38), one additional boundary condition is needed to complete the fifth-order formulation expressed in Eq. (82). In the presentation of his 1998 k-ω model Wilcox [15] states, "The final condition follows from examination of the differential equations for k and ω approaching the surface". For a smooth wall in the limit , the boundary condition requires and . Thus, the differential equation for u+ and the ω-transport equation given in Eq. (82) reduce to
Let the leading-order term in the solution for be written as
where A and a are as yet unknown constants. Using Eq. (84) in the near-wall approximation for the ω-transport equation given by Eq. (83) yields
Equating the exponents and coefficients of y+ in the leading-order terms, this near-wall relation requires
Hence, after using Eq. (86) in Eq. (84), the leading-order solution for yields
To minimize numerical truncation error associated with the singularity, Wilcox [16] suggests that Eq. (87) should be used in place of the ω-transport equation "for the first 7 to 10 grid points above the surface." Wilcox also points out that the grid must be fine enough so that "these grid points must lie below ..." In practice, Eq. (87) is often used as the final boundary condition by applying this relation only at the first grid point off the wall [8].
Because the leading-order solution given by Eq. (87) follows exclusively from the ω-transport equation with application of only the single boundary condition k+(0) = 0, all solutions to Eq. (82) will exhibit this asymptotic behavior, completely independent of the fifth boundary condition that is required to obtain a unique solution to this system of equations. Equation (87) is certainly a valid asymptote for the ω-transport equation in Eq. (82) near a smooth wall. Thus, Eq. (87) can be used as an alternative to the ω-transport equation for y+ approaching zero, provided that it is combined with five appropriate boundary conditions. However, Eq. (87) cannot be used as a substitute for one of the five required boundary conditions.
To show that Eq. (87) is not a valid boundary condition for completing the indeterminate system in Eq. (82), consider the similar system
This indeterminate fifth-order system of differential equations with only four boundary conditions is mathematically similar to Eq. (82), yet it is simple enough to yield a closed-form solution. The general solution to Eq. (88) is
After applying the four boundary conditions given in Eq. (88) to eliminate four of the five arbitrary constants, the solution in Eq. (89) yields
Again as should be expected, there are an infinite number of solutions to any indeterminate fifth-order system of differential equations with only four boundary conditions, such as that specified in either Eq. (82) or Eq. (88). The remaining constant of integration C4 can be evaluated only by applying a mathematically appropriate boundary condition. No amount of analysis applied to the differential equations, no matter how sophisticated, will ever yield a result from which the remaining arbitrary constant in Eq. (90) can be determined.
Notice from Eq. (89) that, analogous to the result obtained from Eq. (82), the general solution for approaches y=0 in proportion to y-2. From examination of either Eq. (89) or Eq. (90), it should be clear that none of the integration constants could ever be obtained by using the asymptotic behavior of for as the fifth boundary condition for Eq. (88). In fact, because the behavior of for depends only on the differential equations in Eq. (88), no boundary condition for can be applied to Eq. (88) at y=0. Likewise, because the near-wall behavior of depends only on the differential equations in Eq. (82), no wall boundary condition for can be applied to Eq. (82). The remaining boundary condition for Eq. (88) at y=0 must be applied to . Similarly, the remaining wall boundary condition for Eq. (82) must be applied to k+.
As presented in Eq. (47), at a smooth wall the correct no-slip boundary condition for completing the fifth-order formulation presented in Eq. (82) is . The analogous boundary condition for Eq. (88) is . It is easily shown that applying this wall boundary condition to the solution of Eq. (88) that is given in Eq. (90) yields C4 = -337/420, and the complete unique solution is
Hence, we see that imposing two wall boundary conditions on and none on is sufficient to determine a unique solution to the coupled fifth-order system of differential equations given in Eq. (88). There is no need to impose a wall boundary condition on , and it is incorrect to do so.
It can be shown numerically that the Wilcox 1998 k-ω formulation given in Eq. (82) exhibits behavior similar to that shown analytically for Eq. (88). For example, Figure 10 and Figure 11 show k+ and ω+ for five solutions, which all satisfy both Eqs. (82) and (87). These converged and grid-resolved solutions were obtained using the same second-order successive underrelaxation algorithms that were used to obtain the k-ε solutions shown in Figure 6 and Figure 7. These results demonstrate that it is mathematically incorrect to use Eq. (87) as the sole substitute for the remaining boundary condition, which is required to obtain a unique solution to Eq. (82). Neither Eq. (87) nor any other relation obtained solely from the differential equations can be used to obtain a unique solution from the in determinate k-ω formulation given in Eq. (82). In addition to using Eq. (87) for numerical implementation, all three of the no-slip boundary conditions that are given in Eq. (47) should be explicitly enforced at a smooth wall.
Like the Lam-Bremhorst k-ε model, the wall damping functions for the Wilcox 1998 k-ω model decay rapidly with increasing y+, so the wall boundary condition has little impact on the predicted velocity profiles. Furthermore, as recommended by Wilcox [15], the use of smooth-wall boundary conditions can be avoided com-pletely with the k-ω model by using the rough-wall boundary conditions first suggested by Saffman [17]. For a smooth wall, Wilcox recommends using the rough-wall boundary conditions with a finite wall-scaled dimensionless roughness height less than 5. As Wilcox points out, the ability to implement rough-wall boundary conditions is a key advantage of the k-ω formulation over the k-ε formulation. For the most recent advancements in the k-ω model, including wall boundary conditions for rough and hydraulically smooth surfaces, see Wilcox [13,18].
Conclusions
Despite the comments by Durbin [6], many of the k-ε turbulence models in common use today are based on incorrect wall boundary conditions, e.g., Eq. (39) or Eq. (46). The correct smooth-wall boundary conditions presented by Durbin [6] and given here in Eq. (47) are seldom used with k-ε turbulence models. Furthermore, the smooth-wall boundary conditions given in Eq. (47) should be used with k-ω and k- ζ turbulence models as well. The eddy-viscosity models used with conventional k-ω and k-ζ turbulence models are obtained directly from the k-ε eddy-viscosity model by using simple changes of variables. Hence, applying a wall boundary condition to ω or ζ is no more correct than applying the equivalent wall boundary condition to ε. The correct smooth-wall boundary conditions for all k-ε, k-ω, and k-ζ turbulence models are those given by Eq. (47). Physics imposes no smooth-wall boundary condition directly on ε, ω, or ζ. As has always been the case, boundary conditions must be obtained from physics. They cannot be developed solely from the differential equations, and they cannot be selected arbitrarily for convenience of numerical implementation.
Because the zero-gradient boundary condition for k in Eq. (47) has been omitted from many commonly used turbulence models, these models are short one boundary condition and thus are mathematically indeterminate. Because any solution obtained from these models is not unique, numerical solutions obtained from these models may be highly implementation dependent. Furthermore, the authors have observed that solutions obtained from some implementations of these mathematically indeterminate turbulence models can be very sensitive to the grids and initial estimates used to obtain the solutions. Although one numerical implementation may converge to a solution that agrees closely with the unenforced boundary condition, another implementation of the same turbulence model could converge to a different solution. The fact that some particular numerical implementation converges to a solution that agrees closely with the unenforced boundary condition cannot be used to justify the implementation of any turbulence model that is mathematically indeterminate. Turbulence models must always be implemented with a complete set of physically correct boundary conditions. None of these boundary conditions can ever be replaced with a mathematical relation that has been developed solely from the differential equations.
Many low-Reynolds-number turbulence models, such as the Lam-Bremhorst k-ε model and the Wilcox 1998 k-ω model, have wall damping functions that decay rapidly with increasing y+, so the smooth-wall boundary condition has little effect on the predicted velocity profiles. However, all turbulence models should be implemented in such a way that is mathematically determinate.