Testing a one-closure equation turbulence model in neutral boundary layers

Highlights

• New universal wall function for 1-closure RANS turbulence model for the channel flow.

• The modeling of the new boundary condition for $k$, as a Neumann Boundary condition.

• Boundary condition on $k$ yields a robust mathematical structure that is analyzed.

• The derivation of a universal formula for the mixing length l from DNS profiles.

• Simulations in a boundary layer to test the model reproduces standard profiles.

Abstract

We aim to test the performances of an incompressible turbulence Reynolds-Averaged Navier–Stokes one-closure equation model in a boundary layer (NSTKE model). We model a new boundary condition for the turbulent kinetic energy $k$ (TKE), and we achieve the mathematical analysis of the resulting NSTKE model. A series of direct numerical simulation are performed, with flat and non trivial topographies, to obtain by interpolation a generic formula for the Prandtl mixing length $\ell =\ell (R{e}_{\star },z)$, $R{e}_{\star }$ being the frictional Reynolds number, and $z$ the distance to the wall. This allows us to carry out numerical simulations at high Reynolds numbers with this turbulence model, in order to discuss its ability to properly reproduce the standard profiles observed in neutral boundary layers, and to assess its advantages, its disadvantages and its limits.

Introduction

The simulation of a turbulent flow by a direct numerical simulation (DNS) using the Navier–Stokes Equations (NSE) remains today – and likely for a very long time – out of reach for a high Reynolds number $Re$. Indeed, the Kolmogorov’s laws imply that $O\left(R{e}^{9∕4}\right)$ degrees of freedom are necessary to do so, which is too large in term of computing power for realistic turbulent flows, such as geophysical flows, the Reynolds number of which is larger than $10^8$. This is why turbulence models are inescapable. Among all turbulence models, two main classes can be distinguished: the Large Eddy Simulation models (LES), such as Smagorinsky’s model, and the Reynolds-Averaged Navier–Stokes (RANS) models, such as the $k-\epsilon$ model (see Chacon–Lewandowski [1], Schlichting–Gersten [2], Pope [3] and Sagaut [4]).

This paper combines modeling, mathematical analysis and numerical simulations. The aim is to investigate the ability a basic incompressible RANS model to faithfully reproduce a neutral boundary layer based on series of numerical simulations, after having modeled suitable boundary conditions and carried out the mathematical analysis of the resulting PDE system. The starting point is a by-product of the $k-\epsilon$ model with only one closure equation, specified by the following equations (Chacon–Lewandowski [1] chapter 6, Wilcox [5]).:

$$\left\{\begin{array}{cc}(\overline{\mathbf{v}}\cdot \nabla )\overline{\mathbf{v}}-\nabla \cdot \left[(2\nu +C_v\ell \sqrt{k})D\overline{\mathbf{v}}\right]+\nabla \overline{p}=\mathbf{f},\hfill & \hfill \left(i\right)\\ \nabla \cdot \overline{\mathbf{v}}=0,\hfill & \hfill \left(ii\right)\\ \overline{\mathbf{v}}\cdot \nabla k-\nabla \cdot ((\mu +C_k\ell \sqrt{k})\nabla k)=C_k\ell \sqrt{k}{\left|D\overline{\mathbf{v}}\right|}^2-{\ell }^{-1}k\sqrt{\left|k\right|},\hfill & \hfill \left(iii\right)\end{array}\right.$$

where “$\nabla \cdot$” is the divergence operator and

• $\overline{\mathbf{v}}=(\overline{u},\overline{v},\overline{w})$ is the long time average of the flow as considered in Berselli–Lewandowski [6] and Lewandowski [7] (or any stationary statistical mean, which does not make any difference thanks to the ergodic assumption about turbulent flows, see for instance in Frisch [8]), $\overline{p}$ the mean pressure, $k$ the turbulent kinetic energy (TKE), $D\overline{\mathbf{v}}=(1∕2)(\nabla \overline{\mathbf{v}}+\nabla {\overline{\mathbf{v}}}^T)$ the deformation tensor,

• $\nu >0$ is the kinematic viscosity of the flow, $\mu >0$ a diffusion coefficient, $\mathbf{f}$ a source term expressing possible external forces,

• In system (1.1), the functions

  $${\nu }_t={\nu }_t\left(k\right)=C_v\ell \sqrt{k},{\mu }_t={\mu }_t\left(k\right)=C_k\ell \sqrt{k},$$

  are the eddy viscosity and the eddy diffusion, $\ell$ is the Prandtl mixing length, $C_v>0$ and $C_k>0$ are dimensionless constants. These eddy functions ${\nu }_t$ and ${\mu }_t$ may be truncated for large values of $k$ for reasons related to mathematical analysis, which will be discussed in the rest of this article.

• the term $C_k\ell \sqrt{k}{\left|D\overline{\mathbf{v}}\right|}^2$ in Eq. ((1.1), (iii)) is the dissipation of the mean flow, generating turbulent kinetic energy, whereas $\epsilon ={\ell }^{-1}k\sqrt{\left|k\right|}$ is the mean dissipation of the fluctuations, damping the TKE.

This type of one-closure equation model can be a good alternative to the full two-closure equations $k-\epsilon$ model, which is expensive and very hard to implement numerically, although very accurate and effective. Evolutionary versions of (1.1) have been used for large scale oceanic simulations (See Blanke–Delecluse [9], Lewandowski [10]), and also in marine engineering to simulate a 2D flow around a fishing net, which has been studied by Lewandowski–Pichot [11].

The model is known as the Prandtl one-equation model (see in Wilcox [5]). Others one equation closures exist in literature such as the Baldwin–Barth model [12], Spalart–Allmaras model [13]. All such models generate correct results for a flat plane bounded flow. An extensive review of RANS model is available in the publication of Argyropoulos and Markatos [14]. The mathematical foundation of the models are available in Wilcox’s book [5]. Recently Fares and Schröder [15] developed another one equation models derived from the $k-\omega$ two equations model with good performance.

To fix ideas, let us consider a 3D flow over a plate, placed at $z=0$. The flow domain is the half space $\{z\ge 0\}$, divided in two regions:

• the boundary layer $\{0\le z\le z_0\}$, $z_0$ being the height of the boundary layer,

• $\{z_0\le z\}$, which is the computational domain and where a turbulent model is implemented.

This raises the questions of the boundary conditions, in particular at the bottom of the computational domain $z=z_0$. In this place, the boundary condition satisfied by the mean velocity is usually a wall law, one of the most popular being the Glaucker–Manning law (see in Chacon–Lewandowski [1, Chapter 5])³ :

$$\overline{\mathbf{v}}\cdot \mathbf{n}=0,-{\left[(2\nu +{\nu }_t)D\overline{\mathbf{v}}\cdot \mathbf{n}\right]}_{\tau }={\alpha }_v\left|\overline{\mathbf{v}}\right|\overline{\mathbf{v}},\text{at }z=z_0.$$

A more general set of wall functions can be found in Kalitzin et al. [16].

When the friction Reynolds number is small, the Dirichlet boundary condition can be applied for the velocity (Wilcox [5] Parente–Gorlé [17]).

The natural boundary condition for the TKE is,

$$k={\left|\overline{\mathbf{v}}\right|}^2\text{at }z=z_0,$$

which connects in a non linear way the TKE to the velocity at the bottom of the computational domain.

All this, equations and boundary conditions, yields a turbulence model. If the main question is to assess the ability of such a model to simulate turbulent flow from a practical point of view, we must begin to question its mathematical structure, which leads to purely theoretical considerations.

Mathematical motivation and new boundary condition. We know that the resulting boundary value problem given by the system (1.1) and the boundary conditions (1.3) (wall law for $\overline{\mathbf{v}}$) and (1.4) (boundary condition that links $\overline{\mathbf{v}}$ and $k$), for a given smooth function $\ell =\ell \left(\mathbf{x}\right)$, yields serious mathematical and numerical complications mainly because of the relation (1.4) between $\overline{\mathbf{v}}$ and $k$ at the bottom. In particular, we have never been able to demonstrate the existence of a solution to this problem in the evolutionary case, and although we have a proof in the steady-state case as considered here (see in Chacon–Lewandowski [1]), the proof is particularly technical and long, which is not satisfactory. This is the reason why we looked for an alternative condition to (1.4), which decorrelates the velocity of the TKE, while preserving the physical characteristics of the flow and leading to a better treatment from the point of view of mathematical analysis.

It is important to note that this mathematical analysis that we undertake in this paper, namely the theoretical question of the existence of solutions to the system of partial differential equations used to simulate the flow, has a considerable impact on the question of stability and consistency of the numerical schemes used for the numerical simulations, which amounts to a theoretical validation of these simulations, regardless of DNS results and/or any in situ data, which unfortunately we do not have here.

By adapting the reasonings which allow the obtaining of the wall law for the velocity such as (1.3), we get in this paper the following nonlinear law for the TKE at the bottom of the numerical domain:

$${\mu }_t\left(k\right)\frac{\partial k}{\partial n}={\alpha }_{k}k\sqrt{k},\text{at }z=z_0,$$

where ${\mu }_t\left(k\right)$ is the eddy viscosity diffusion, such as given by Eq. (1.2), ${\alpha }_k$ a dimensionless coefficient that we have calibrated according to the DNS. Considering this boundary condition instead of (1.4), we get a more affordable mathematical structure than that provided by the relation (1.4) between $k$ and $\overline{\mathbf{v}}$. In particular, we are able to prove in this paper, by a proof not too complicated, the existence of a weak solution to the system of Eqs. (1.1), with the boundary conditions (1.3) (wall law for the velocity) and (1.5) (new boundary condition for the TKE). This is Theorem 3.1 below, the proof of which is developed in Section 3, which is the main theoretical result of this paper.

Remark 1.1

The resulting model, given by Eqs. (1.1), (1.3), (1.5), is appointed by the acronym NSTKE.

Remark 1.2

Homogeneous Neumann boundary condition may also be prescribed for the TKE at $z=z_0$, that is

$$\frac{\partial k}{\partial n}=0,$$

see in Cindori-Juretié [18], Liu [19], which comes within the scope of our theoretical study, in a simplified way.

Numerical framework and mixing length. Once the structure of the new model is understood, the next step is to evaluate its numerical performances inside a boundary layer, which means taking $z_0$ at least of order of the height of the viscous sublayer. We will do so in the case of a flat bottom, which is the simple case with which to start, then in the case of a non trivial topography as displayed in Fig. 1, called the rough case. The aim is to test topology that are small compared to the turbulent boundary layer height such as grass, sand (or wave for oceanic boundary layer). Problems related to buildings or hills are not treated in this paper.

Another motivation is that the model is easy to implement. The mixing length $\ell$ is computed by the product of two functions that depend only on the distance, determined by DNS and an extrapolation method, which yields the need to determine four dimensionless coefficients, that we assume to only depend on the friction Reynolds number. For comparison, Spalart–Allmaras [13] model is closed with eight closure coefficients and three damping functions. Baldwin–Barth model seven closure coefficients and three damping functions, Wilcox [5]. Furthermore, the proposed model is universal. The implementation of the wall function is also easy in solvers like OpenFoam.

It is important here to explain how the mixing length $\ell$ is determined, which is one of the main contributions of this study. In the full $k-\epsilon$ model, it is deduced from $k$ and $\epsilon$ by the standard formula

$$\ell =\frac{k^{3∕2}}{\epsilon },$$

From DNS data of small and high Reynolds turbulent channel flow, $k$ and $\epsilon$ are obtained assuming that $k$, $\epsilon$ and $\ell$ are homogeneous in the $x-y$ axes, we evaluate $\ell =\ell (z,R{e}_{\star })$ at the grid points, properly averaging the data raised from the DNS. We then interpolate the collected sets of numerical values to get a general formula in both flat and rough cases (see formula (4.17), (4.18), complemented by (4.19), (4.21), (4.22) for the calculation of the different dimensionless coefficients).⁴ Our DNS are compared to the DNS of Moser et al. [20], [21], which serves as the benchmark for our results.

To return to performance, comparing to other models, note that eddy viscosity computed by both Spalart-Allmaras and Baldwin-Barth models do not decay at the end of the turbulent boundary layer (see in Wilcox [5]). In our model, given the expression for the turbulent mixing length, the computed TKE and the turbulent mixing length decrease (and remain positive) at the edge of the boundary layer leading to a decaying eddy viscosity. This fundamental property avoids the extra dissipation at the summit of the boundary layer as seen in algebraic models (0 equation closure) Smith–Cebeci [22] and Baldwin–Lomax models [23].

Numerical results. With numerical formula (4.17), (4.18) for $\ell$, several numerical simulations with the NSTKE model have been performed up to $R{e}_{\star }=10000$, in both flat and rough cases, after having evaluated the roughness coefficients contained in the boundary conditions. There is no universal turbulence model, and this one is neither better nor worse than another. However, compared to others, the existence of weak solutions provides us some guarantees. Note also that the numerical simulation is as important as the model itself, and consistency analysis of the scheme used should be also performed. This an open problem left for future studies. In this work only a numerical assessment of the convergence has been performed.

The NSTKE model behaves properly in the flat case, which validates our approach. However, the results are less good in the rough case. This does not mean that the model and our approach have reached their limits, and the present study opens several questions. We clearly do not have enough DNS in the rough case for the determination of $\ell$. Moreover, in this case, the assumption $\ell =\ell (z,R{e}_{\star })$ must be called into question, in favor of $\ell =\ell (x,y,z,R{e}_{\star })$. Last but not least, this first series of results shows that the topography should be strongly taken into account in the calculation of the roughness coefficients, and it remains an open question to know how to find a universal and simple way for doing this, which is one of the main challenge in the field (see the “bulk algorithm” in Pelletier [24] for instance).

Organization. This paper is organized as follows. We carry out in Section 2 the modeling of the new TKE boundary condition (1.5), based on a standard assumption about the eddy viscosity in the subviscous layer and a Taylor expansion, following the usual outline for deriving wall laws.

We prove in Section 3 an existence result to the resulting NSTKE system. We briefly review the state of the art about NSTKE systems. Then we carefully set the functional spaces and the norm we use (see Remark 3.1), that allows us to define the notion of weak solutions for this system (see (3.63), (3.64), (3.65)), following the same pattern as in Chapter 7 in Chacon–Lewandowski [1]. One of the main characteristics of the problem is the equation for k which is an elliptic equation with a right hand side (r.h.s) in $L^1$ whose boundary condition is a non linear Neumann Boundary condition. Such elliptic equations with r.h.s in $L^1$ and nonlinear Neumann BC have not been considered before so far we know. This is why we prove in Lemma 3.2 an estimate à la Boccardo–Gallouët [25] for such equations, which is the main novelty in this analysis part.

Section 4 is devoted to display several DNS, that yields to numerical formula for the Prandtl mixing length $\ell$ (see (4.17), (4.18)), as well as the setting of the constants.

We show in Section 5 the results of the simulations from the NSTKE model, using the recursive algorithm given in Section 5.1, both in the flat and in the rough cases. We also carry out a numerical convergence analysis of the algorithm.