On Boundary-Value Problems for RANS Equations and Two-Equation Turbulence Models

Currently, in engineering computations for high Reynolds number turbulent flows, turbulence modeling continues to be the most frequently used approach to represent the effects of turbulence. Such models generally rely on solving either one or two transport equations along with the Reynolds-Averaged Navier–Stokes (RANS) equations. The solution of the boundary-value problem of any system of partial differential equations requires the complete delineation of the equations and the boundary conditions, including any special restrictions and conditions. In the literature, such a description is often incomplete, neglecting important details related to the boundary conditions and possible restrictive conditions, such as how to ensure satisfying prescribed values of the dependent variables of the transport equations in the far field of a finite domain. In this article, we discuss the possible influence of boundary values, as well as near-field and far-field behavior, on the solution of the RANS equations coupled with transport equations for turbulence modeling. In so doing, we defne the concept of a well-defined boundary-value problem. Additionally, a three-dimensional, rather than a simpler one-dimensional analysis is performed to analyze the near-wall and far-field behavior of the turbulence model variables. This allows an assessment of the decay rate of these variables required to realize the boundary conditions in the far field. This paper also addresses the impact of various transformations of two-equation models (e.g., the model of Wilcox) to remove the singular behavior of the dissipation rate (ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document}) at the surface boundary. Finally, the issue of well-posedness regarding the governing equations is considered. A compelling argument (although not a proof) for ill-posedness is made for both direct and inverse problems.


Introduction
In this paper we consider the Reynolds-averaged Navier-Stokes (RANS) equations. The RANS equations are based upon Reynolds expansions for the dependent variables and time (Reynolds) averaging. In practice, this definition is extended to include applying a combination of Reynolds expansions (for the pressure and density) and mass weighted (Favre) expansions for the other dependent flow variables (see for example the textbooks of Wilcox and Pope [1,2]). In so doing, the time-averaged equations can retain the same form as their laminar flow counterparts. The resulting mean flow momentum equations are simply augmented by an additional stress term (a Reynolds stress), and the mean flow energy equation includes an additional heat flux term (Reynolds heat flux) and work term due to the Reynolds stresses. To close the system of equations, some form of modeling is required to determine the Reynolds stresses and Reynolds heat fluxes.
Currently, in engineering applications, partial differential or integral equations of a transport type are generally used in modeling the effects of turbulence (see [3][4][5][6][7][8][9]). One approach for modeling the additional terms appearing in the RANS equations, such as Reynolds stresses, is to invoke the eddy viscosity hypothesis (the Boussinesq approximation) in which the Reynolds stress tensor is proportional to the mean strain rate tensor. This constitutive relation requires the determination of a viscosity coefficient that is called an eddy viscosity. The eddy viscosity depends on two scales: a velocity scale and a length scale. These scales can be determined using transport equations. Usually, turbulence modeling formulations use one or two equations to determine the scales. For one-equation models, the velocity scale comes from solving for either the turbulence kinetic energy k or a quantity that can be related to k. The length scale for these models comes from an algebraic method that depends on physical location and whether the flow is a wall-bounded shear layer or a free shear layer. Two-equation models are often considered since the two scales can be computed from the two transport equations, one for k and the other for a quantity related to the length scale, such as the dissipation rate of k. Two advantages in using one-equation or two-equation models are that they avoid the necessity of solving directly for the Reynolds stresses and simplify considerably the modeling and solution requirements (e.g., typical Reynolds stress models require solving six equations for the Reynolds stresses and one equation for the length scale). The focus in this paper concerns two-equation models.
In the development, transformation and improvement of turbulence models the main focus is the variation of the equations themselves. This approach is interesting from the point of view that the main drivers for solutions are the boundary values. Roughly speaking, the predictions using the model cannot become better than the boundary values enforced. To understand possible shortcomings and behavior of the RANS equations in combination with turbulence models, the impact and influence of the formulation of the boundary conditions cannot be overstated. Both the complete boundary-value problem and the influence of the boundary values on the solution need to be understood. For example, the boundary values need to be chosen in such a way that it is at least possible to compute a solution. When one does not know in which way to choose the boundary values correctly, that is the user needs to make a guess, one cannot in general expect to find a solution of the corresponding boundary-value problem.
Against the background of this approach the following questions arise • In what ways do modifications to the turbulence model cause changes in boundary conditions? • What assumptions are made to derive boundary conditions for the turbulence model.
• What impact do these assumptions have on possible solutions? • Can the prediction behavior of a turbulence model be satisfactory if the formulation of boundary conditions is unsatisfactory?
The principal objective of this paper is the presentation and discussion of the complete boundary-value problem for the RANS equations coupled with two transport equations for turbulence modeling. In this paper, the meaning of a well-defined boundary-value problem is one in which the equations and boundary conditions are defined, as well as the nearfield and far-field behavior of the solution variables. Moreover, the consistency requirements between the specified values of the dependent variables and their decay behavior in the interior of the flow field are addressed in detail. As part of this overall objective, it is a goal to draw attention to which way boundary conditions already establish properties of solutions to the equations. An additional subject pertaining to this system of PDEs is the effect of transforming the transport equations in an effort to remove the singular behavior of the dissipation (or dissipation rate) at a solid surface boundary. In solving any system of PDEs it is critical to establish some understanding of well-posedness. This subject is also considered. An example is given that gives a strong indication of the ill-posedness of the problem governed by the transport equations of the turbulence model.
The question of solvability for the complete boundary-value problem yields directly to important mathematical questions. Namely, it is essential to return to the fundamental issue of constructing a well-posed problem. This issue is existential in solving both the Navier-Stokes and turbulence modeling equations. As pointed out previously, we cannot in general establish a well-posed problem for even the Navier-Stokes equations, in the sense of Hadamard [10] and supplemented with the condition of regularity. This mathematical issue already sets the stage for the issue of an ill-posed problem for the system of transport equations. Such a consequence exposes the probability of anomalous results or multiple solutions. At this point, there has been very little documented concerning the existence of multiple solutions when considering various turbulence models. Moreover, there is currently no systematic way of demonstrating multiple discrete solutions, as well as how many such solutions may exist, when modeling turbulence.
It is important to recognize that there are examples of multiple solutions in fluid mechanics, which have occurred when solving the transonic potential flow equation [11] and the Euler equations [12,13]. There is also the discussion of Temam [14] concerning existence, uniqueness and regularity for the Navier-Stokes equations. In fact, Temam shows an example of non-uniqueness of solutions for the Navier-Stokes equations (i.e., regarding laminar flows) in which he provides proof based on topological arguments.
More recently, in 2014, Kamenetskiy et al. [15] presented numerical evidence of multiple solutions to the RANS equations for turbulent flows. They consider turbulent flow over an extruded two-dimensional (2-D) airfoil geometry and a trap wing configuration from the first AIAA High Lift Workshop [16]. Computations were performed for both the one-equation model of Spalart-Allmaras and the two-equation model of Wilcox. Three primary approaches were considered to find multiple solutions, which are as follows: Varying initial conditions, applying different techniques for enforcement of solid surface boundary conditions, using variations of implicit residual smoothing. For the trap wing case, the number of multiple solutions discovered with the methods depended on the mesh density. The meshes used (designated coarse, medium and fine) were from the AIAA High Lift Workshop. It should be pointed out that in the work of Kamenetskiy et al. numerical results are only classified as solutions when the residuals of the equations have been reduced to machine zero. As discussed previously, reduction of the residual to machine zero is essential. Otherwise, results in which the residual has only been reduced four or five orders of magnitude are affected by solver integration errors that can play a significant role in the behavior of the result, possibly making it quite different from the actual solution. Kamenetskiy et al. call such results pseudosolutions. It is interesting that on a coarse mesh eight solutions were found, while only two solutions were discovered on the fine mesh. This suggests a connection to spatial resolution. Thus, the question immediately arises as to how many multiple solutions, if any, will be obtained when all the relevant scales of the turbulent flow are resolved. With a well resolved flow field (on a sufficiently fine mesh), it is distinctly possible that the number of multiple solutions is much greater than two solutions, which is a consequence of only small variations in the solutions. A final point is that Kamenetskiy et al. observed that the appearance of multiple solutions seems to be closely related to smooth body flow separation that occurs for flows over high-lift configurations. Here, we note that such occurrences of non-unique solutions may only be one member of the set of possible turbulent flows which have multiple solutions.
The authors have encountered both non-unique solutions as well as anomalous (nonphysical) results for the particular case of circulation control airfoil flows [17,18]. Being aware of the fact that the forward problem, that is to find a solution of the boundary-value problem, is possibly ill-posed, we may ask the question vice versa: Is turbulence modeling its self ill-posed? And if so, what kinds of consequences follow from such a conclusion.
The article is organized as follows. Section 2 presents the governing equations of fluid flow, and Sect. 3 considers the additional equations required to model the effects of turbulence. Boundary conditions together with the corresponding boundary-value problems are discussed in Sect. 4. Proposed methods in the literature to deal with problems arising from the boundary conditions of the turbulent flow equations are introduced in Sect. 5. An interpretation of turbulence modeling as an inverse problem is given in Sect. 6. Numerical examples supporting the assertions made in this article are presented in Sect. 7.

Governing Equations of Fluid Flow
To describe flow behavior we consider for the domain D ⊂ R m , m = 2, 3, i.e., an open and connected set, and an interval [0, T ) ⊂ R, T > 0, the RANS equations in conservative form expressed in integral form by where the integral operators V D and R ∂ D are given by and denotes the vector field of conserved variables and n is the unit outward normal on ∂ D. The terms f c and f v describe the convective and viscous contributions Here e i is the ith unit vector. The quantities ρ, u = (u 1 , . . . , u m ) T , E and are the density, the velocity, the specific total energy, and the enthalpy of the fluid. The equation of state defines the pressure p, and γ is the gas dependent ratio of specific heats, which is given by 1.4 for air. Assuming that an effective viscosity is given and, using Stoke's hypothesis, that the bulk viscosity satisfies λ = −2/3μ eff , the viscous stress tensor τ = τ (W ) = τ (W (x, t)) is given by and S denotes the strain rate, which is given by the symmetric part of the total derivative of flow velocity vector u, Throughout the paper Id(x) = x denotes the identity operator. The missing viscous flux term for the energy equation is given by The effective viscosity μ eff and effective conductivity κ eff are computed by and the laminar viscosity is given by Sutherland's law κ l (W ) := c p μ l (W ) Pr l and whereby ρ ∞ > 0 and u ∞ > 0 denote some constant reference density and velocity, L > 0 is some constant reference length scale, Re > 0 is the corresponding Reynolds number, is Sutherland's constant, is the universal gas constant and the laminar Prandtl number is given by Pr l := 0.72. In this article, we restrict ourselves to linear turbulence models represented by differential or integral equations. The solutions of these equations reveal additional quantities in the considered fluid. These occurring variables extend the degrees of freedom given by the conservative variables W by a further unknown function Here N t ∈ N depends on the turbulence model. In this report we have for the Spalart-Allmaras model, for the Wilcox k-ω model and the SST model.
Since we mainly deal with two-equation turbulence models, one can assume N t = 2 throughout this report. The Spalart-Allmaras model is mentioned here for completeness since results used for comparisons are presented in Sect. 7. Note, when for example using an algebraic model N t = 0, no additional unknowns need to be considered. The additional variables required for the turbulence model are then used to to determine the eddy viscosity, which is required for (8). Given the eddy viscosity μ t the turbulent thermal conductivity is described by the algebraic relation Pr t := 0.92.
Formula (12) is also required for (8). The determination of μ t closes the system of the Reynolds-averaged Navier-Stokes equations (1a).

Two Equation Turbulence Models of k-! Type
In this section, we introduce the two-equation k-ω models of interest, which are often used for aerodynamic applications. These are the k-ω model of Wilcox introduced 1988 and the Shear Stress Transport (SST) model by Menter [19]. Both models include differential equations for the unknown, positive functions where k is turbulence kinetic energy per unit mass and ω the dissipation rate of turbulence kinetic energy. Throughout the literature two-equation turbulence models are typically presented in their differential form, though in finite-element and finite-volume codes an integral form is used. Before we state the two transport equations for k and ω, we define t = t i j 1≤i, j≤m and the Reynolds stress tensor τ RS = τ RS i j 1≤i, j≤m , using the strain rate tensors S and S, given in (6) and (5), and the mean-molecular-stress sensor, t = t i j 1≤i, j≤m . Then, Additionally, according to S we define the vorticity as skew-symmetric part of the total derivative of flow velocity u,

Wilcox k-! Model
The k-ω model of Wilcox [5] has the form The eddy viscosity, also called turbulent viscosity, is computed by The source terms for the k-equation and ω-equation are given by Here and throughout the rest of this article the symbol : The constants of the model are There is a 2006 version of the k-ω model of Wilcox given in [6]. For the purposes of this article, it is sufficient to consider the original version of the model.

Menter's SST Model
Menter proposed the Shear-Stress-Transport model [7,19,20]. This model modifies the formula for eddy viscosity μ t , diffusive term and the source term Q when compared with the k-ω model. An update of the SST model was suggested in [8], but this article restricts to the presentation of the original model. The model reads as and the eddy viscosity is defined via Here, F 2 is a blending function defined by and d is the distance to the closest no-slip wall. According to [7], the constants are a 1 = 0.31 and C 1 = 500. The source terms include a diffusion term, where Pr k,SST := τ RS : Pr ω, Di ω,SST : The SST-model involves a blending of a k-ω and a k-ε model. This blending is controlled by a function = (x; ε 1 , ε 2 ). This function is designed to detect the edge of the boundary layer, such that the SST-model behaves inside the boundary layer like a k-ω model and outside like a k-ε model, exploiting the convex combination : To realize smooth blending, the function F 1 : [0, ∞) → [0, 1] is modeled using the hyperbolic tangent, where F 1 is determined by Using the function , functions k , ω , γ , and β required for the viscous flux in (19) are given by and the constants are The last constant v k is well known as the von Kármán constant in the literature.

Boundary-Value Problems
So far, we have only stated the equations of interest. That is, the mean flow equations (1) together with a system of equations to determine the required eddy viscosity μ t , for example (16) or (19). Naturally, for a closed representation, we need to formulate a corresponding boundary-value problem.

Boundary Conditions for Mean Flow Equations
Due to a lack of theoretical understanding of both mean flow and turbulent flow equations, the definition of boundary values and conditions at infinity are not straightforward. For example, for a complete and closed formulation the decay behavior at infinity for ρ, u, p, and additionally even for or k and ω, is required. Since this is in general unknown, we prescribe these values formally. For the representation of the exterior boundary-value problems of interest, we introduce the formal setting, Here 0 ≤ ϕ ≤ π and 0 ≤ ψ < 2π are the Angle of Attack (AoA) and a possible side slip angle. Furthermore, in the sequel, let D ⊂ R m , D = ∅ be a bounded domain, and for the sake of simplicity, we assume that the boundary of ∂ D is connected and that ∂ D is an orientable submanifold of R m of dimension m − 1. Then, the classical (adiabatic) no-slip wall boundary conditions are in the sense of a trace operator. In the following, we also write ∂ D = ∂ D no-slip to clarify the meaning of ∂ D.
Though we have stated the RANS equations in their unsteady form, we are only interested in approximating a steady-state solution. Hence, we formulate the boundary-value problems only for the steady-state. Therefore, no initial condition must be considered as well as timedependent boundary values.

No-Slip Wall Boundary Condition
To realize boundary conditions for k and ω is not straightforward. The only simple boundary condition corresponding to vanishing velocity u ∂ D no-slip = 0 is To derive a boundary condition on no-slip wall for ω we follow the presentation in [1,5].
We assume that near a no-slip wall a solution to Navier-Stokes equations is incompressible and pressure is constant, and as a consequence, convective terms are negligible. Then the equations for k and ω simplify to In a next step, it is assumed that only the velocity gradient in the direction normal to the no-slip wall is dominant. If this direction is identified with x 2 -coordinate, we obtain and as a consequence we have Finally, it is assumed that the production terms are negligible compared to the remaining terms. Then, using (31) to conclude that ν t vanishes on the no-slip wall, all what is left of the ω-equation is the ordinary differential equation which has the solution For generalization of this procedure in the normal direction n, we formally derive under the assumptions mentioned above for ω the boundary condition To investigate the behavior of k in a neighborhood of a no-slip wall, as a direct consequence, we can apply the same assumptions to the equation for k to obtain Searching for a solution of the form k(x) = x α 2 , a comparison of left and right hand side gives and therefore Assuming that the second solution α 2 is non-physical, otherwise k would have singular behavior at a no-slip wall, a contradiction to (31), we obtain which corresponds exactly to the behavior determined in [5,21]. Here C k,no-slip > 0 denotes some constant, which cannot be determined in general. In contrast to ω, the behavior of function k is identified in a neighborhood of a no-slip wall only up to some scaling C k,no-slip . On the other hand, under the assumptions formulated above, we do not only have the boundary condition k ∂ D no-slip = 0, but also, in a neighborhood of a no-slip wall, we have the quantitative behavior lim Remark 1 Under the assumptions mentioned above and imposing boundary condition (32) for ω in a neighborhood of a no-slip wall, k satisfies condition (33), and the value α 1 is determined by the relation of β * and β.
this observation also determines the behavior of the fluctuating part of velocity u in a neighborhood of a no-slip wall. In a neighborhood of a no-slip wall, velocity u satisfies approximately the quantitative behavior Remark 2 Under the assumptions mentioned above and imposing boundary condition (32) for ω, in a neighborhood of a no-slip wall, both the behavior of k and the behavior of u are determined and need to satisfy conditions (33) and (34). To say it directly, the specification of no-slip wall boundary condition and near wall behavior for ω control the near wall behavior of k and u . Hence, physical situations which do not follow conditions (34), (33) and (32) can in general not be simulated using the above boundary condition for ω.

Farfield Boundary Condition
The presented derivation of far-field boundary conditions is similar to the one in [22]. The authors in [22] restrict themselves to a one-dimensional analysis. Here, we generalize the approach for the three-dimensional space.
To determine far-field boundary conditions, (k ∞ , ω ∞ ) we assume for x 2 → ∞ constant velocity u ∞ , and constant density ρ ∞ > 0. Mathematically, this translates into lim where in spherical coordinates ⎛ Then, for sufficiently large x 2 , we get from (16a) and (16b) and assumption (35) Additionally, it seems reasonable to postulate that in the free-stream, far away from the obstacle, the diffusion terms in (16) become small; that is, we assume lim Then, for x 2 → ∞, system (16) simplifies to Using the normalized normal vector Eqs. (37) and (38) can be reformulated by A solution of (40) is given by where a 1 , a 2 , a 3 ∈ R are coefficients. Straightforward differentiation can be used to show that (41) is a solution of (40). For example, this solution corresponds to the one given in formula (9) in [22]. The direct correpondence is derived below and given in (50).
Here, C ω is some constant, which cannot be determined in general. But C ω is exploited to determine reasonable free-stream values for k and ω at the inflow boundary. This is the topic of Sect. 4.3 given below. In this section, we are only interested in the asymptotic decay behavior of ω for x 2 → ∞. Hence, here it is not of importance to quantify C ω .
Using the representation (41) of ω at infinity, we get from (39) the differential equation and with some constant C k > 0. Again, straightforward differentiation can be used to verify that (42) is a solution of (39). To obtain the 1D analysis, for example considered in [22], we may specify n = (1, 0, 0) T . Then, (38) and (39) simplify to their 1D counterparts shows that the derived behavior of k and ω at infinity presented here fits into the 1D analysis and represents therefore a generalization.
From (42) and (41), it is obvious that k and ω have singular behavior on the hyperplane and the singularity is characterized as hyperbolic. According to E, we define To reach the goal of formulating a closed exterior boundary-value problem, two different ways seem to be appropriate.
a) Without loss of generality, we can assume that the domain of interest satisfies D ⊂ V + . Then we need to prescribe for k and ω a certain behavior according to solutions (42) and (41) when approaching E from the interior of V + . b) In a neighborhood of E, assumptions (35) and (36) do not hold, and (39) and (40) are not valid simplifications of either (16) or (19). Therefore, it can be reasonable to assume that k and ω as solutions of either (16) or (19) do not have singular behavior when approaching E, but instead k and ω decay at infinity. Following this argumentation, we can postulate that ω and k as solutions of either (16) or (19) decay at infinity as lim For the k-ω model of Wilcox we have and for the SST model we approximately get lim since lim As a direct result, we obtain for the eddy viscosity Finally, we are in a position to formulate the complete boundary-value problem for the RANS equations for steady-state problems: Considering possibility b) mentioned above, we may formulate: Exterior turbulent flow problem, formulation 1:

Find a function W † that satisfies the steady RANS equations in
and satisfies the (adiabatic) no-slip wall boundary conditions (30), and uniformly for all directions. Additionally, find a function W t that satisfies the k-ω turbulence model in R m \ D, and satisfies the boundary conditions in the sense of (32) and (33) and uniformly for all directions in the sense of (44) for k and (43) for ω. Considering possibility a) mentioned above, we may formulate: Exterior turbulent flow problem, formulation 2: Find a function W † that satisfies the steady RANS equations in V + \ D, that is and satisfies the (adiabatic) no-slip wall boundary conditions (30), and Additionally, find a function W t that satisfies the k-ω turbulence model in V + \ D, and satisfies the boundary conditions in the sense of (32) and (33) and Without knowledge about the behavior of k and ω as solutions of either (16) or (19), it is impossible to say which of the formulations given above is preferred. The second formulation is typically closer to the application.

Practical Considerations for Far-Field Boundary Conditions
To get an estimate of the constant C ω in (41), we fix a point x fs = x fs 1 , x fs 2 , x fs 3 in the free-stream (fs), that is, we assume that x fs is so far away from ∂ D that assumptions (35) and (36) are satisfied. Then denotes the free-stream value of ω at position x fs and determines the constant Inserting this into representation (41), we get Again, considering direction n = (1, 0, 0) T and assuming that x fs = (0, 0, 0), the explicit formula for ω simplifies to which is the three-dimensional counterpart to formula (9) in [22]. To determine the decay behavior of ω along rays through the domain, we define a ray in the inflow direction through the point x fs by and compute ω (y (t)) = ω fs ω fs β In practice, we have a finite domain, and we need to prescribe values for k ∞ and ω ∞ . Then, for example, k fs = k ∞ and ω fs = ω ∞ hold along the outer boundary of the finite domain correspondig to inflow, where we have assumed that the outer boundary is located so far away from ∂ D such that (35) and (36) If ω fs is chosen too small, we have ω (y (t)) ≈ ω fs on the whole computational domain, and no decay behavior can be expected and observed. Thus, we could not expect that the equations for turbulence will have impact on the RANS equations. Therefore, since t corresponds to a length scale which corresponds to the Reynolds number, we obtain from (51) Based on this argumentation, such a choice of ω f s should guarantee a correlation (51) to observe on a given, finite domain the desired decay behavior (43) for ω. Using (48) and (42) we have which gives inserting (49), To determine a reasonable value for k fs we estimate using (52), Simplifying, β * /β ≈ 1, we obtain from (54) Finally, to understand the role of the constant C k , we estimate the turbulence intensity in the free stream by In this sense, C k determines the turbulence intensity in the free-stream.
Relations (55) and (52) that is C k = 0.009. As a consequence, we get This fixed value can in general only hold and be prescribed at the inflow boundary. Assuming the derived decay behavior of k and ω at infinity is correct, the fixed choice (56) is in general wrong along the outflow boundary, and boundary-value problems formulated using (56) along the outflow boundary can in general not be solved.
To demonstrate that the decay behavior (43) for k and (44) for ω is indeed observed in applications, we apply (53) and (50).
To show that in the free stream a solution for k and ω is not significantly influenced by the airfoil and the decay of k and ω behaves as predicted, solutions for k and ω along the line y(t) = tan(AoA)t + const. are shown in Fig. 1. Here the symbol AoA denotes the Angle of Attack, which is used to determine the values for W ∞ and W t,∞ (see Sect. 4.1). Considering the choice (56) and the  [23]. As seen in Fig. 1, these values are approximately observed in the computation. Also, as predicted, the decay rate of (46) for k in the SST model is slower than that of (45) for the k-ω model. A second situation which needs to be considered is a line which follows a path in a neighborhood of the airfoil through the wake. A plot of k and ω for such a situation is given in Fig. 2. From the free-stream value of (56), both k and ω decay as expected. Then, in a neighborhood of the airfoil, both increase suddenly. Inside the wake there is decay as expected, but on a significantly higher level when compared with Fig. 1. Figure 3 gives a detailed view on the decay of k and ω in the wake. These values are not observed in Fig. 3. This is expected, since inside the wake, the assumptions (35) and (36) are not satisfied because the outflow boundary is located too close to the aifoil surface.
As a consequence, in integral form the diffusion terms in (16) and (19) vanish. Moreover, denoting the last inner point at the outflow boundary by x in , we can represent where we have used (57). Hence, the convective part of (16) and (19) can be implemented without prescribing unknown, outflow values. So, using boundary conditions at the outflow (57), we avoid the prescription of unknown values. This procedure is used for both types of outflow boundaries mentioned above. It is in general not straightforward to separate the outer boundary into a section corresponding to inflow and into another one corresponding to outflow. Such separation is usually only possible if the solution of the given problem is known. For the examples considered in this article, we used C-type meshes, where at least a large amount of the outflow boundary can be determined geometrically. A more general approach to determine the outflow boundary is to consider lines of the form (49). Given the angle of attack, for each control volume along the outer boundary the entry intersection and exit intersection can be determined. All entry intersections are marked as inflow and all exit intersections as outflow. The realization of such a procedure and the investigation of the influence on solutions and solution algorithms is part of future work. Finally, we want to close this section with the two following remarks:

Remark 3
One often finds in implementations the possibility to define a certain relation of eddy viscosity to laminar viscosity (see for example [24]), i.e., it is possible to prescribe Input parameter at inflow = μ t,∞ μ l,∞ .
Such an input parameter defines only one condition for determining either k fs or ω fs ; and hence, a second condition is required. This second condition needs to be chosen such that (52) and (55) hold. Though possible, the advantage of prescribing μ t,∞ μ l,∞ is not obvious. Hence, here we directly define k ∞ = k fs and ω ∞ = ω fs along the inflow, and to conclude from this the relation μ t,∞ /μ l,∞ .

Remark 4
The decay behavior for the eddy viscosity in (47) was not formulated in the assumptions of this Sect. 4.2.2, but it is a consequence. The constants β * and β determine the decay behavior at infinity.

Practical Considerations for No-Slip Wall Boundary Conditions
To implement in the discrete problem the quadratic singular behavior (32) for ω is not straightforward. Instead, we follow the idea given in [7,19]. Note, when approaching a smooth no-slip wall, the asymptotic behavior is determined by (32), which can be reformulated as Since h represents the distance to the closest no-slip wall in the normal direction, it was suggested to represent this expression numerically by where p i,bdr y denotes the point on the no-slip wall and p i,n the closest, next discrete point in direction −n p i,bdr y . Then, d i is the distance to the closest no-slip wall of the first discrete point in the interior of the flow field. To take care of the quadratic singular behavior, it was suggested to multiply this value with an additional order of magnitude. Hence, the no-slip boundary value for ω is realized by Now, to obtain the value for ω required to evaluate the boundary flux, a linear extrapolation is done. The increase (i.e., the gradient into normal direction of ω for the boundary edge e i,bdry ) is approximated by Then, in a small neighborhood of p i,bdry , the function ω is approximated by The right-hand-side of (59) is used to evaluate the boundary flux for a no-slip wall. The construction for k is significantly simpler. Given a value for k on the no-slip wall (k i,bdry ), we define as the state for the no-slip wall k no-slip p i,bdry := −k i,bdry .
To evaluate the boundary flux for a no-slip wall these values are averaged; and hence, k = 0 is enforced in the flux, which realizes the boundary condition for k.

Discussion About Existence and Uniqueness of Solutions
Obviously, the boundary conditions for k and ω on a no-slip wall and the decay behavior at infinity influence a solution of the Navier-Stokes equations (1a) together with either (16) or (19). But how restrictive are the constraints on possible solutions? That is, can we expect that under these conditions the exterior turbulent flow problem is well-posed in the sense of Hadamard [10]. Considering the problem not being well-posed it could mean that a) no solution exists or b) the solution is not unique or c) the solution does not depend continuously on the given data.
Even if only one of these points is fulfilled, data produced with a computer code is in general practically useless. And just because the suggested and implemented boundary conditions allow for approximate solutions for a few preselected test cases, it does not allow for the conclusion that they work in general. Moreover, the analysis considered above shows that the decay behavior (44) of k at infinity is different for the original k-ω model and the SST-model. Hence, in general, changes in a model have impact on the behavior of variables at the boundary. This observation makes it obvious that it is generally insufficient to consider changes in the equations of the model without analyzing the impact of these changes to the formulation of boundary values.
Besides these theoretical considerations, in practice values for k and ω and / or their derivatives need to be prescribed at all boundaries and consistently implemented. On a given finite mesh the knowledge of a certain decay behavior is a good indicator, but it is not sufficient to find suitable values. In this sense, a certain implemented choice is only a guess, which has no claim to generality. At the no-slip wall the situation is not much better. The quadratic singular behavior of ω needs to be realized in practice. Any numerical implementation can simulate such behavior. All methods suggested and implemented so far try to deal with the deficiencies of this behavior at the boundary. Furthermore, as shown by the analysis, the behavior of ω determines the decay behavior of k up to a constant. And it is an open question if k can exhibit such behavior in general. Without a good understanding for the choice of boundary values, it is highly possible that problems are ill-posed.
Furthermore, we would like to draw attention to the following fact. The behavior of function k in a neighborhood of the no-slip wall (33) is a consequence of the behavior of ω given in (32). Also, in the free stream, the behavior of k given by (44) is also a consequence of ω given in (43). Hence, both at the no-slip wall and in the free-stream it is the behavior of ω which determines significantly the solution of the turbulence model. Since boundary values are important for determination of the solution to the corresponding boundary value problem, the equation for ω has particularly strong influence.

Variable Substitution
When solving two-equation models such as the k-ω type or k-ε type, one typically encounters the problem of ensuring that turbulence variables remain positive during the iterations (see (13)). In general, negative values of k and / or ω directly yield a breakdown of the iteration resulting in "Not a number". To avoid such problems, it was suggested in [25] to substitute variables k and ε or ω by In particular, the sole substitution for ω found its way into several implementations, most often in the background of discontinuous Galerkin methods. Here, the idea has been picked up originally by Bassi et al. [26] and reused in [27]. Recently, a further application of such substitution was used in [28]. Naturally, substitutions (60) only make sense for nondimensional k, ε and ω. This is not a severe restriction, but depending on the implementation care must be taken. For example, typically implementations want to support flexible restart options. Then, knowledge about the output variable such as its kind of nondimensionalization and scaling is necessary to convert the variable.
Further substitutions discussed in the literature are (see [29][30][31]) These substitutions were in particular introduced to deal with the singular behavior of ω indicated in (32).

ln(k)-and ln (!)-Formulation
To better understand the logarithmic reformulation and its consequences we present a possible derivation. Starting with equation (16b), we substitute ω = e and apply the chain rule to obtain where we have exploited conservation of mass For the diffusive term, one computes in a similar manner div (μ l + σ ω μ t ) grad e = e div ((μ l + σ ω μ t ) grad ) where we have used div grad e = e grad Therefore, formally we obtain from (16b) and using the substitution for ω, the equation for given by Division by e yields the equation Using (18b) the source terms are explicitly given by Using the substitution k = e K in the same way, the differential equation for K can be derived and the source terms are explicitly given by

Equivalence to Original k-! Models
Now, the following questions need to be answered. Consider supplementing the system of equations (16) or ( Is the obtained exterior turbulent flow problem equivalent to the original one; that is, we assume that , or K or ( , K) is a solution of the boundary-value problem with the modified system of equations. Further, does it follow that ω, or k or (k, ω) is a solution of the boundaryvalue problem with non-modified system of equations (16) or (19)?
The answer for equivalence of k-equation is straightforward. Since k represents turbulence kinetic energy, the no-slip wall boundary condition is (33). But the consequence of substituting e K for k is k > 0 rather than 0 at the no-slip wall boundary. Therefore, equation (66) does not allow for solutions of k, which need to satisfy the no-slip wall boundary condition. Hence, for boundary-value problems which rely on such boundary conditions, substituting e K for k cannot be realized equivalently.

Remark 5
Substitution of k by e K cannot be equivalently implemented considering a boundary condition k ∂ D = 0. Such substitution might be considered when wall functions are used. However, due to the discussion above, we do not support that substitution of k is a general and successful way to deal with problems involving k-ω models.
Considerations about substitution of ω = e are more interesting. To realize the boundary condition for ω given by (32), one needs to realize a quadratic singularity, i.e.
Using substitution (60) for ω one obtains for Hence, the quadratic singularity for ω has been converted into a logarithmic singularity for . It can be assumed that a logarithmic singularity is numerically better to realize than a quadratic one. And hence, at first glance the substitution for ω is promising. In detail, at least two questions need to be answered: a) Are the ansatz functions for ω better suited to approximate a function with logarithmic singularity than a function with quadratic singularity? b) How does the numerical error behave when comparing a quadrature rule for a logarithmic singularity to a quadratic singularity?
On the other hand, due to the substitution and transformation, an additional source term appears in (65) which needs to be discretized, namely Because of (67) near a no slip wall we have Therefore, using the substitution ω = e converts the quadratic singularity for ω into a logarithmic for . On the other hand, the singular quadratic behavior was simply shifted to another term. Finally, we shortly summarize the arguments of this paragraph: a) Due to a solid wall boundary condition, the k substitution of K cannot be realized equivalently. b) Quadratic singular behavior at a no-slip wall for ω is converted to logarithmic singular behavior for . Then, the grad quadratic singular behavior needs to be realized.
To show that there is also no advantage of substitutions (61) we observe that due to (43) lim Obviously, the singular behavior has just been shifted from the no-slip wall to an unbounded increase for τ and g in the far field. To deal with increasing values in the far field is just as big a challenge, if not an even bigger one. On a finite mesh, comparably to (56), an arbitrary, generally a guess, is required. Unfortunately, the authors of [29][30][31] leave this question open. Concluding, there is no obvious advantage in using variable substitution. Original problems are only shifted. The equivalence of the original system of equations and the system of equations obtained after variable substitution is at least questionable in the sense that so far there exists no mathematical proof to even analytically ensure the equivalence.

Turbulence Modeling: An Inverse View
To evaluate the formulation and the accuracy of a turbulence model is not an easy task. From our perspective the following four points are a minimum standard one has to consider: a) The full differential or integral formulation of the model, b) its exact implementation, c) a solution algorithm which is able to compute for a given number of degrees of freedom a machine accurate solution and d) mesh converged results.
As soon as one of these criteria is not satisfied, certain doubts about the assertions made concerning a turbulence model arise. Unfortunately, throughout the literature about computational fluid dynamics, the explicit details of the computations are often hidden, and convincing arguments about convergence are also often missing. In particular, information about the actual formulation of boundary conditions as well as the possible impact of certain limitations of variables to stabilize the solution process is often missing. Even when all these criteria are satisfied and there is full evidence about the implementation, a strict conclusion about the accuracy of computed results is difficult to obtain. Typical validation measures are the comparison with measurements, which also come from a process which is inaccurate. This observation about this complex topic motivates one to possibly reconsider turbulence modeling as a parameter identification problem, typically an inverse problem, which is ill-posed.
To look at this point of view, we do not need to consider derivation of the RANS equations. Instead, we consider the system of equations (1a) such that it includes an unknown parameter μ t ≥ 0, which needs to be determined. For the simplest choice, which assumes μ t = 0, it is in general impossible to find steady-state solutions of (1a) for high Reynolds numbers. High Reynolds number computations using equations (1a) and μ t = 0 are called Direct Numerical Simulations (DNS), and the integral curves W lam = W lam (x, t) satisfying (1a) are in general unsteady.
There are many reasons why it is often not possible to numerically approximate W lam directly. Besides the mathematically open problem to prescribe meaningful physical initial values in general, another severe issue is the numerical effort to compute a DNS solution for high Reynolds number turbulent flows. The number of mesh points required (proportional to Re to the 9/4 power) and the corresponding degrees of freedom are extremely large. Thus, in general, such a computation cannot be realized due to computational reqirements and effort. Hence, the question arising is as follows: Does there exist an appropriate modification of the Navier-Stokes equations such that the modified equations exhibit three conditions: a) A steady-state solution exists and it is unique (also for high Reynolds numbers). b) The steady-state solution represents important features of W lam . c) It is possible to compute numerically an approximation to the steady-state solution within an adequate time interval.
Considering that the system of equations (1a) with μ t = 0 allows for steady-state solution when the Reynolds number is sufficiently small, for example 0 < Re ≤ 10 4 , we follow the idea to decrease (at least locally) the Reynolds number. This goes along with the idea that a proper weighting of the diffusion terms needs to be incorporated to enforce steady-state solutions. From the perspective of (1a), this means that artificially the laminar viscosity μ l needs to be increased. This can be realized by adding an additional function, for example denoted by μ t , yielding (8) and (12). Hence, it is the goal to find a a function (parameter) μ t such that the above mentioned criteria a)-c) are satisfied, and additionally, certain specific flow characteristica are satisfied. In applications these are often for example a) pressure distributions, b) skin-friction distributions, c) velocity profiles, d) the location of separations.
We summarize these characteristics as given data Y . In general, these data are measurements; and therefore, perturbed by noise. Hence, we assume that instead of Y , there is only noisy data Y δ available, which satisfies In general, the noise level δ ≥ 0 is unknown.
In this sense, it is the goal of turbulence modeling to prescribe the unknown eddy viscosity μ t ≥ 0 such that characteristics of interest are represented by a solution to the exterior turbulent flow problem. To give the determination of the unknown function μ t a more general view, we make the following assumption. We assume that μ t (x, t) ≥ 0 and μ t is constant outside of the ball B d := {x ∈ R m : x 2 ≤ d}. The characteristics of interest are often data defined on ∂ D, that is, Y ∈ L 2 (∂ D). Within this notation we may formulate turbulence modeling as an inverse problem. The forward problem might be stated as follows. Assume that an eddy viscosity (or at least a method to determine μ t ) is given, and then find a solution of the exterior turbulent flow problem. Note, in this general formulation, equations implicitly defining μ t are not required any more. One can simply assume the existence and knowledge of μ t . Therefore, the finding of a function W t would be eliminated in our formulation for the exterior turbulent flow problem.
The inverse problem is to reconstruct the eddy viscosity μ t from given data Y ∈ L 2 (∂ D). To formulate the inverse problem mathematically, consider the operator which maps the eddy viscosity to the corresponding data, for example a pressure distribution, such that the exterior turbulent flow problem is satisfied. Here, we assume that the unknown function μ t ∈ D(F), where D(F) is the domain of the operator F. Naturally, when we assume that given data are velocity profiles, the image space of F needs to be replaced by a data-adapted one. Though a given turbulence model is an a-priori known methodology to determine μ t , calibration of model parameters and validation is performed by comparison with noisy measurement data for a small number of test cases where data is available such that computed data is in good agreement with these measurements. A given turbulence model is a realization of the operator F −1 . And the problem is well-posed if F is bijective and F −1 is continuous, otherwise ill-posed.
The interpretation of turbulence modeling as an inverse problem gives rise to several questions, the most severe might be the one of well-posedness. Already, the number of turbulence models in the literature yielding similar reconstructions for data such as the pressure distributions and even skin-friction coefficients suggest non-uniqueness of the problem. Also, a small perturbation in the given pressure distribution may produce a significantly different solution of the exterior turbulent flow problem, and therefore, also a significant difference in the eddy viscosity. This suggests additionally that the condition of continuous dependency on the data is violated. These observations indicate that turbulence modeling by its mathematical nature is ill-posed. Therefore, the finding of a general ansatz to represent characteristics of turbulent flows for a whole variety of problems using the RANS equations in combination with a given methodology determining μ t seems to be challenging, and may even be impossible.
In practice, measurements Y δ are available for almost none of the test cases of interest. So, at the moment a practical way for reconstructing μ t requires the knowledge of data independent methods, such as an additional set of equations. Taking into account the possibly ill-posed mathematically nature of the problem, one should be aware of the fact that in general good agreement between noisy measurements and computed data cannot be expected.
Due to the increase in computational power and the ability to realize a scale resolving simulation, we close this section with the following considerations for future work. Assume that a method describing the eddy viscosity μ t ≥ 0 is given, for example an established turbulence model. Let us denote the steady-state solution of the corresponding exterior turbulent flow problem by W † turb . Furthermore, the associated possible time-dependent solution of the flow problem with μ t = 0 is given by W † lam . We introduce the error induced by the turbulence model within the time interval [T 0 , T 1 ] by where . denotes a suitable norm. Such an error measure does not only include surface data but also it includes the entire solution field. We can only expect that err turb (T 0 , T 1 ) can be small if there exists a continuous relationship between the solutions W † lam of the laminar and W † turb of the Reynolds-averaged Navier-Stokes equations. That is, the change of the solution due to the modification of the set of equations by the eddy viscosity should lead to small differences in the solutions, mathematically speaking The evaluation of such errors might be helpful to improve the characterization as well as potentials and shortcomings of established models and those developed from scratch. The argumentation of this paragraph carries over to full Reynolds stress models. Here, instead of reconstructing one parameter μ t , six terms defining the Reynolds stress tensor need to be reconstructed.

Numerical Example
To confirm that different eddy viscosities can produce almost exactly surface data, we consider a numerical example. For the discretization strategy as well as a solution methodology, we refer to [23,[32][33][34][35]. The software used to generate the results is explained in detail in the literature mentioned. As proposed in [19], a production limiter is incorporated into the k-ω turbulence models; that is, the production term of the k-equation is replaced bỹ Pr k,SST := min Pr k,SST , 20De k,SST .
The limiter is constructed such that the value of production does not exceed a multiple of destruction. To demonstrate an effect of these limiters, we consider turbulent flow around the RAE 2822 airfoil with respect to the conditions given in Table 1. The mesh is of size 320×64. This particular airfoil is quite important since it is frequently considered to evaluate turbulence models. Moreover, there is a substantial amount of detailed data from the experimental investigation of Cook et al. [36]. This data provides the necessary basis for discussing additional issues in this paper. Figure 4 shows the computed eddy viscosity for Case 9. Note, both results correspond to fully converged solutions. The corresponding convergence histories are given in Fig. 7. The solution without production limiter shows an increase in eddy viscosity in a neighborhood where the shock impinges on the boundary layer. Such an effect is suppressed when using the production limiter, as shown in Fig. 4 (right). The result is reproducible using another implementation [32]. A structured grid code is applied, which uses a different discretization strategy. The results for the computed eddy viscosity are shown in Figs. 5 and 6. Here a different scaling of the plotted values is used.
The same observation is true for Case 10. Given two fully converged solutions (see Fig. 7 (right)), the one with production limiter shows an increase in eddy viscosity in the region where the shock impinges on the surface boundary layer, which vanishes when the production limiter is used.
These numerical results indicate that both test cases have (approximately) no impact on computed C p and C f distributions. It should be pointed out that essentially the same C p variation implies that the boundary-layer displacement thickness is also about the same. Further, since the surface skin-friction variations are essentially the same and the shear stress plays a key role in determining the law-of-the-wall behavior, this suggest that there is a good correlation with the streamwise velocity profiles. A plot of the C p and C f distributions is given in Fig. 8 for Case 9 and Fig. 9 for Case 10. Obviously, though the eddy viscosity with and without production limiter (68a) is significantly different, its impact on computed C p and C f distributions is not discernible.
To summarize, this short investigation shows that the introduction of a production limiter (68a) can have significant impact on the computed eddy viscosity. On the other hand, this significant change in the computed eddy viscosity, did not yield significant changes in C p and C f distributions. Such an observation can be interpreted in two ways.
a) The influence of computed eddy viscosity is so weak, that if a fully converged solution is possible, no noticeable change in the final solution of mean flow equations can be found. b) Vice versa, it is interesting to notice that significant changes in the computed eddy viscosity do not yield significant differences in C p and C f distributions. Such an observation gives rise to the question: In which way does the eddy viscosity need to be changed such that it influences the solution of the mean flow equations.
In addition, from the viewpoint of an inverse problem, we need to emphasize that this example reveals that the same data can produce multiple solutions for the eddy viscosity.
Returning to the subject of finding a suitable eddy viscosity μ t , we notice that two different functions for μ t generate the same C p and C f surface distributions. That is, the classical way of turbulence modeling to find a methodology to determine an eddy viscosity such that noisy data is matched cannot verify a unique determination of the eddy viscosity. For example, we could make a numerical experiment where we consider the computed C p and C f distributions as given and ask for the reconstruction of a corresponding eddy viscosity. Obviously, both computed eddy viscosities, that is the one with and the one without production limiter are solutions of this question. Hence, this example indicates that the inverse problem for these examples does not have a unique answer.
Moreover, we mention that in formulae (68a) and (68b) a factor of 20 was introduced for limitation. However, throughout the literature one finds other factors such as 10 or 5 for limitation of production compared to destruction. The results obtained were similar.

Concluding Remarks
This article has considered boundary-value problems for the RANS equations in combination with two-equation turbulence models. Classical boundary conditions for the unknowns k and ω have been discussed, including the assumptions formulated under which these boundary conditions were derived. In so doing, we have assigned a meaning to what one can call a welldefined boundary-value problem. Such a boundary-value problem reveals a characterization of the near-field and far-field solutions of the governing equations. This can be extremely important for providing insight into what to expect when solving the corresponding discrete problem.
For example, it has been shown and discussed that determination of boundary values of ω for both the no-slip wall and the far field plays a significant role in determining the behavior at the boundaries for k. Consequently, the behavior of the eddy viscosity μ t in the boundary layer is determined to a large extent by the surface boundary condition for ω. With respect to this knowledge, it can be assumed that reformulation and manipulation of the equations for the turbulence models will only have limited influence.
In this sense, an improvement of turbulence models might only be possible when significantly better understanding of the complete boundary-value problem is available. Any model including an equation for dissipation rate ω is restricted to the shortcomings which are inherent to this equation and its boundary values. In particular, it cannot be excluded that the choice of k and ω at the boundaries yields an ill-posed problem.
Furthermore, it has been demonstrated that proposed variable substitutions do not help in general to remove the singular behavior of the equations. The singular behavior is typically only shifted to other terms. Inclusion of a variable substitution gives generally a rise to many more questions than it answers. For example, there is the important question about equivalence of the original boundary-value problem to the one with a substituted variable.
Finally, it has been discussed in which way classical turbulence modeling can be interpreted as an inverse problem, that is, to reconstruct a parameter in a boundary-value problem such that given noisy measurements are approximated. Based on this approach, it has been discussed that the original mathematical nature of turbulence modeling is possibly ill-posed. A numerical example is used to confirm this assertion.