Smooth-Wall Boundary Conditions for Energy-Dissipation Turbulence Models

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-{\epsilon}, k-{\omega}, and k-{\zeta} 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 {\epsilon}, {\omega}, or {\zeta} at a smooth surface and it is incorrect to do so. The behavior of {\epsilon}, {\omega}, or {\zeta} approaching a smooth surface is that required to satisfy the differential equations and force both k and its gradient to zero at the wall.


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 energydissipation 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 , 2 1 k and the dissipation length scale ε 2 3 k . The parameters k and ɛ are defined in terms of the fluctuating velocity as 1 1 where Ṽ is the fluctuating velocity vector, ) (V J   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 µ C 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 ) ( k C µ ε ω ≡ 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 µ f , 1 f , 2 f , E , and ɛ 0 . To complete any k-ɛ model of this form, the wall damping functions are specified as prescribed functions of v, , 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 x V and y V 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 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, 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 + ν . Thus, from Eq. (11) Applying Eqs.
To complete the formulation, the damping functions µ f , 1 f , 2 f , + E , + o ε , 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 + o ε , 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  (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 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 R t and R y 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, 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) 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  [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," 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 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 impli cations 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 noslip boundary condition places no restriction on the derivative of x Ṽ 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 Because the no-slip boundary condition requires 0= V at the wall, Eq. (48) requires Hence, we see that a no-slip wall imposes two boundary conditions on k and none on ɛ. This is sufficient to determine a unique solution to the coupled fourth-order system of k-ɛ transport equations. There is no need to impose a wall boundary condition on ɛ, and it is incorrect to do so. The value of ɛ at a smooth wall is that required to satisfy both Eqs. (18) and (49), as was originally pointed out by Durbin [6].
At this point it may be useful to return to our consideration of Eq. (42), which has a closed-form solution and is mathematically similar to the fifth-order system in Eq. (38). The boundary condition for Eq. (42) that is analogous to Eq. (49) is which results in . Applying this boundary condition to the solution of Eq. (42) that is given in Eq. (44) yields the somewhat more troubling result 23 210 0 − = , 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 fifthorder 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 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 + q 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 ) 0 ( + ε and ) 0 ( + θ must be obtained from the differential equations. The process is started with initial estimates for ) 0 ( + ε and ) 0 ( + θ . From these initial estimates, fourth-order Runge-Kutta integration can be used to obtain ) ( + + l q and ) ( + + l θ . The initial estimates are then refined using an appropriate numerical method until the solution is found, which corresponds to the correct centerline values 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 ,  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 2 From the first two differential equations in Eq. (57) and the specified relation for + ν , 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 To examine the near-wall behavior of the Launder-Sharma k-ɛ model, consider the Taylor-series expansions Using these expansions in Eq. (57), the near-wall expansion for the k-transport equation yields      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.  and (65) were solved numerically using a secondorder 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 + y 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: To demonstrate that traditional implementations of these turbulence models do not necessarily converge  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  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 Similarly, we see from Figure 7 that these finitedifference and finite-volume implementations of the Launder-Sharma model converge to different solutions with traditional implementa tions of the wall boundary conditions. However, for the particular implementation used to obtain the results shown in Figure 7, the indeterminate finitedifference 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.  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 Because the wall damping functions for the Lam-Bremhorst model decay rapidly with increasing y + , the wall boundary condition k 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 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 smoothwall 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 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][11][12][13].

Application to k-ω Models
The k-ω turbulence models that are commonly used for CFD are built on exactly the same dissipationbased eddy-viscosity model that is given in Eq. 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 1 ω C = 0 and he did not include the molecular diffusion term. The closure coefficients often used for the k-ω model [8,15]  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 µ f , k f , 1 f , and 2 f , could be specified as prescribed functions of v, , 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.
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) 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].
Notice that the turbulent dissipation Reynolds number R t , 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 wallscaled dimensionless variables that are defined in Eqs.  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 0 → + y , the boundary condition . Thus, the differential equation for u + and the ω-transport equation given in Eq. (82) reduce to 1 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 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 5 . 2 = + y …" 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  Hence, we see that imposing two wall boundary conditions on k 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 0 ) 0 ( = +' k 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 k has little effect on the predicted velocity profiles. However, all turbulence models should be implemented in such a way that is mathematically determinate.