Preprint
Article

The Basic k-ϵ Model and a New Fundamentally Based Model of Anisotropic Inhomogeneous Turbulence Compared with DNS of Channel Flow at High Reynolds Number

Altmetrics

Downloads

68

Views

18

Comments

0

A peer-reviewed article of this preprint also exists.

This version is not peer-reviewed

Submitted:

02 March 2024

Posted:

04 March 2024

You are already at the latest version

Alerts
Abstract
Predictions of mean values of statistical variables of large scale turbulent flow by the widely used basic k-ϵ model and a new fundamentally based model are verified against published results of Direct Numerical Simulations DNS of the Navier-Stokes equations. The verification concerns turbulent channel flow at shear Reynolds numbers of 950, 2000 and 104. The basic k-ϵ model is largely based on empirical formulations accompanied by calibration constants. This contrasts with the new model where descriptions of leading statistical quantities are based on the general principles of statistical turbulence at large Reynolds number. Predicted values of major output variables such as turbulent viscosity, diffusivity of passive admixture, temperature, and fluid velocities compare well with DNS in case of the new model. Significant differences are seen in case of the basic k-ϵ model.
Keywords: 
Subject: Computer Science and Mathematics  -   Applied Mathematics

1. Introduction

The description of turbulence has been an issue right from the beginning. A solid starting point for analysis are the Navier-Stokes equations which describe the flow. These can be time averaged resulting in equations for mean values of velocity, pressure and temperature. Problem are the average values of the products of the fluctuating quantities in the non-linear convection terms in the equations. They describe the effect of fluctuations on mean flow quantities. Proposed representations of average nonlinear convective fluctuations have been of a hypothetical nature by drawing analogies with molecular chaos. Boussinesq [1] was the first to follow this line of thinking by introducing the gradient hypothesis: the averaged nonlinear fluxes are proportional to the derivative of the mean flow quantity proceeded by a constant termed turbulent viscosity or turbulent diffusion coefficient. Several versions and extensions based on the same idea have since then been put forward by the pioneers of turbulence theory: Taylor [2], Prandtl [3], and Von Karman [4] among others. They form also the basis of many of today’s computer models used in engineering and environmental analysis: Hanjalic and Launder [5], Bernard and Wallace [6].
An offspring of the afore mentioned concept is the basic k- ϵ model widely used in engineering and environmental analysis [5,6]. The average value of momentum fluxes is described by the gradient of the mean velocity. Fluctuations are isotropic and the coefficient proceeding the gradient equals c μ k 2 ϵ where k is mean kinetic energy of fluctuations, ϵ is mean energy dissipation rate and c μ is a calibration constant. The model is completed by equations for k and ϵ . The gradient hypothesis is also applied to several other flux terms in the equations with different values of the calibration constants.
Problem with the gradient based models is their hypothetical origin and lack of uniqueness. Starting from a sketchy analogy with laminar diffusion almost limitless forms of gradient representations and proceeding functional dependencies can be devised. A new development which surpasses these limitations is the statistical description of anisotropic inhomogeneous turbulence: Brouwers [7,8,9]. Turbulence at high Reynolds number is featured by unstable eddies of whirling irregular fluid velocities whose behaviour is governed by domination of the inertial forces in the momentum balance: Monin and Yaglom, Vol. II, Ch.8 [10]. The eddies start at sizes of the configuration and break down to ever smaller ones until the point where viscous forces become into play at sizes of a few millimeter. Here the theory of small viscous scales of Kolmogorov comes into play [10]. Building on this framework of high Reynolds turbulence and using the methods of stochastic theory of Van Kampen [11] and Stratonovich [12] explicit statistical descriptions of governing variables are obtained. They are the leading terms of asymptotic expansions based on small value of the inverse of the universal Lagrangian Kolmogorov constant. Results represent unique descriptions of turbulent flow statistics up to deviations due to truncated higher order terms.
An opportunity for testing the outcome of the turbulence models is provided by the results of Direct Numerical Simulation DNS of the Navier-Stokes equations. With the evolution of modern computer power it has become possible to generate a wealth of accurate data on turbulent flow at high Reynolds number. Most interestingly are the recent DNS data of turbulent channel flow of Hoyas et all [13,14] and Kuerten et all [15]. Fluctuating channel velocities are strongly fundamentally based and averages of flow quantities vary strongly with distance from the wall. Inhomogeneity and anisotropy are the characteristics of turbulence in practice. The DNS data thus provide a meaningful test case for models. In this paper a detailed comparison is presented of the DNS data with the predictions of the basic k- ϵ model and the new fundamentally based model, also referred to as fundamental model.

2. The Basic k- ϵ Model and the New Fundamental Model

Considered is turbulent flow of an incompressible fluid or an almost incompressible fluid, e.g. a liquid or a gas flowing at speeds where the square of the Mach number is small. The density ρ is taken constant. Turbulent fluctuations measured at a fixed point in space are treated as a statistical process which is stationary or almost stationary in time compared to the time of velocity fluctuations. Statistical averages follow from time averaging over sufficiently long time intervals. The time averaged representation of the Navier-Stokes equations is given by
  • Conservation of mass:
u i x i = 0
Conservation of momentum:
u i t + u j u i x j + < u i u j > x j = p x i
Conservation of energy:
θ t + u j θ x j + < u j θ > x j = 0
where angled brackets represent statistical averaging, t and x are time and space coordinate, u, p, and θ are mean or time-averaged values of fluid velocity, pressure, and temperature and u and θ represent fluctuations of velocity and temperature, that is velocity and temperature minus their mean value. As we are concerned with flow at high Reynolds number the contributions of the viscous forces and heat conductivity present in the Navier-Stokes equations have been dropped. Their effect can be disregarded when considering the main flow governed by instability of inviscid flow outside small boundary layers The average temperature in energy equation (3) can also be used to describe the average distribution of passive or almost passive admixture in the fluid. Restricting the formulation of the conservation equations to an incompressible or almost incompressible fluid implies that the solution of (1) and (2) is not affected by the value of the conservative scalar temperature. The representation holds as long as the changes in temperature or admixture imposed at external boundary conditions are of limited magnitude.

2.1. Turbulent Diffusion in the Basic k- ϵ Model

The appearance of turbulent fluxes in the convection terms of the averaged conservation equations results in an unclosed set of equations for mean flow variables. It is known as the closure problem. To resolve this issue a diffusion hypothesis has been introduced in the basic k- ϵ model. In this hypothesis turbulent fluxes are treated as isotropic and are described by [5,6].
< u j u i > = 2 3 k δ i j ν t ( u j x i + u i x j )
< u j θ > = ν t θ x j
where ν t is a scalar which represents diffusivity or turbulent viscosity and which is defined by
ν t = c μ k 2 ϵ
and where k is average kinetic energy of fluctuations, k = 1 2 < u i 2 > and ϵ = 1 2 ν < ( u i x j + u j x i ) 2 > is average energy dissipation rate with ν being kinematic viscosity; c μ in (4) and (5) is a calibration constant whose value is usually taken as 0.09 [5,6]. In Equation (5) correction factors are sometimes added to the diffusion constant, i.e. a turbulent Prandtl number and a turbulent Schmidt number in case of temperature and admixture respectively. But these numbers are generally close to unity and are omitted here.

2.2. Turbulent Diffusion in the Fundamentally Based Model

The derivation of the fundamental model starts from a Langevin equation for fluid particle velocity [8,9,10]. In this equation the limiting form of Kolmogorov’s theory of small scales is implemented, i.e. the limiting form of this theory when times characteristic for the velocity fluctuations are much larger than the time of the viscous scales. This is the case when the Reynolds number is large: Kolmogorov [16]. A further step is an expansion in terms of the inverse of the Kolmogorov constant C 0 . Matching predictions with data of measurements and DNS reveals values of C 0 around 6 7 ; in the present analysis a value of 7 is adopted. Furthermore the area where the Lagrangian description applies can be reduced to a point in the Eulerian flow description in the limit of small C 0 1 . In this way Lagrangian based descriptions are connected to Eulerian ones. The following descriptions for the flux terms are obtained [8,9]
< u j u i > = 2 3 ( k + D n k u n x k ) δ i j D i k u j x k D j k u i x k
< u j θ > = D i k θ x k
where the fundamentally based diffusion tensor D i j is described by
D i J = 2 C 0 1 ϵ 1 σ i n σ n j + 2 C 0 2 ϵ 2 σ l i σ j k u n σ l k x n 4 C 0 2 ϵ 1 σ k j u n x n ( ϵ 1 σ i m σ m k )
and
σ i n = < u i u n >
is co-variance or Reynolds stress. Relations (8) - (10) are part of the description which holds in the entire flow configuration (except from the thin viscous layers at walls) once they are coupled to conservation equations (1) - (3). In this way change of turbulent flux at each point x as described by (7) and (8) is connected to change in kinetic energy and energy dissipation. Although the presence of ν in the expression for energy dissipation ϵ = 1 2 ν < ( u i x j + u j x i ) 2 > may suggest otherwise, ϵ is a characteristic of the main inviscid flow outside the boundary layers [10]. The magnitude of gradients of velocity is governed by the small viscous scales of turbulence. It scales as 1 / ν and makes the magnitude of ϵ independent of ν . This independency is reflected in the equations for k and ϵ presented in the next section.
The above expressions for diffusion of the fundamental model reveal dependency on mean gradients which are of a more complex structure than those of the basic k- ϵ model. It reflects the anisotropy of the fluctuating velocity field and reveals a more complex dependency of flow statistics.
Application of the fundamentally based diffusion approximation to a scalar as done in Equation (8) is only justified if the scalar is a conservative scalar [8,9]. The value of the conservative scalar is constant when following a fluid particle and fluctuates in value in a fixed coordinate system only due to fluctuations of the fluid particle. This is the case for temperature in an incompressible fluid such as liquids. It is approximately correct if the fluid is almost incompressible as is the case in gases flowing at speeds where the square of the Mach number is small. The scalar representation can also be applied to passive or almost passive admixture in fluids such as aerosols in air [7,8,9]. It leads to errors when applied to non-conservative scalars such as kinetic energy and pressure: see section 5.

3. Equations for k and ϵ

Implementing the expressions for turbulent diffusion of the basic k- ϵ model and the fundamental model in the averaged conservation equations introduces two unknowns: the mean kinetic energy k and the mean energy dissipation rate ϵ . Equations for k and ϵ can be obtained from the Navier-Stokes equations [8,9]. Our aim is to describe the flow away from thin viscous controlled layers near walls being in size a few millimeters only. In line with this approach the contributions in the equations for k and ϵ from laminar viscosity will be disregarded ( just as we did in conservation equations (1), (2) and (3)). To provide boundary conditions at the wall the viscous layer is surpassed by applying the solutions of the log layer at x = 0 : [5,6].
  • The equation for k reads as [5,6]
k t + u i k x i + < u i k > x i + ρ 1 < u i p > x i = P ϵ
where P is mean production of turbulent fluctuations defined by
P = σ i j u i x j
and where k and p are the fluctuating parts of kinetic energy and pressure, respectively, that is the kinetic energy and dissipation rates minus their time-averaged values. There are two turbulent flux terms in Equation (11), i. e. the third and fourth term on the LHS of (11), which need to be modelled. In the basic k- ϵ model both terms are lumped together [5,6] and are described by:
  • Basic k-ϵ model
< u i k > + ρ 1 < u i p > = ν t σ k k x i
where σ k is a calibration constant which is usually taken unity: σ k = 1 [5,6].
  • Fundamental model
The theory underlying the fundamental model provides general expressions for turbulent scalar fluxes which are free from calibration factors. However, these expressions are only valid for conservative scalars and lead to disagreement with DNS results when applied to turbulent fluxes of kinetic energy and pressure [8]. A fall back to empirical construction is needed. It reads as
< u i k > + ρ 1 < u i p > = c k D i j k x j
where it is noted that if k were a conserved scalar its diffusion should be described by D i j k x j . The factor c k represents correction for non-conservative behaviour and includes the relative small contribution [8] of pressure diffusion. The above relations for diffusion of kinetic energy and pressure will be compared and calibrated with DNS results in Section 5.
The equation for ϵ conventionally applied in CFD models [5,6] is largely an empirical construction. Its basic k- ϵ form follows from the Navier-Stokes equations. It contains a number of terms which are governed by the small viscous scales. These terms are generally replaced by expressions which meet the criteria of matching to the results of decaying grid turbulence and the log layer of turbulent channel flow. The equation reads as
  • Basic k-ϵ model
ϵ t + u i ϵ x i = x i ( ν t σ ϵ ϵ x i ) + ( c ϵ 1 P c ϵ 2 ϵ ) ϵ k
Fundamental model
ϵ t + u i ϵ x i = x i ( D i j σ ϵ * ϵ x j ) + ( c ϵ 1 P c ϵ 2 ϵ ) ϵ k
The constant c ϵ ensures matching with the case of grid turbulence. Its value is usually taken to be 1.9 which is somewhat less than the theoretical limit value of 2 for infinite Reynolds number: George [17]. The von Karman constant κ is equal to 0.4 . In case of the basic k- ϵ model the calibration constant σ ϵ is usually taken to be around 1.3 and c ϵ 1 is specified by the equation (eq. (8.41) of [16])
c μ σ ϵ κ 2 ( c ϵ 2 c ϵ 1 ) = 1
The values appropriate for c ϵ 1 and σ ϵ * in case of the fundamental model are determined in section 6.

4. Channel Flow

Objective is to compare the results of the basic k- ϵ model and the fundamental model with those of DNS of channel flow. The channel consists of parallel planes in between which the mean velocity u 1 is unidirectional in the direction x 1 parallel to the planes. Its magnitude and the value of statistical averages related to velocity fluctuations only change in the direction x 2 normal to the planes. For channel flow there is an exact solution for the mean pressure [9] which reads as
p ρ = u τ 2 x 1 H σ 22
where u τ is shear velocity and 2 H is channel height. The shear velocity is determined by the pressure drop in direction x 1 at given u 1 . There are theoretical and experimentally confirmed relations for the relationship determining u τ . Another exact result for channel flow is the description of the co-variance σ 12 or < u 1 u 2 > .
σ 12 = u τ 2 ( 1 x 2 H )
which is valid outside the thin viscous layer at the wall [9].
In the subsequent analysis of channel flow we shall make use of dimensionless formulations: u i is made dimensionless by u τ , σ i j and p ρ by u τ 2 , x 1 and x 2 by H, P and ϵ by u τ 3 H ; the subscript2 of x 2 will be dropped. In this new notation Equations (18) and (19) become
p ρ = x 1 σ 22
which is valid outside the thin viscous layer at the wall.
σ 12 = ( 1 x )
Furthermore Equation (12) becomes
P = ( 1 x ) d u 1 d x
The above solutions follow from the averaged momentum equations adopted to the case of channel flow [9]. A complete specification of all variables follows from the expressions for turbulent fluxes and the equations for k and ϵ . In case of channel flow these reduce to:
  • Basic k-ϵ model
σ 11 = σ 22 = σ 33 = 2 3 k
σ 12 = ν t d u 1 d x
< θ u 2 > = ν t d θ d x
ν t = c μ k 2 ϵ
d d x ( ν t σ k d k d x ) + P ϵ = 0
d d x ( ν t σ ϵ d ϵ d x ) + ( c ϵ 1 P c ϵ 2 ϵ ) ϵ k = 0
Fundamental model:
σ 22 = 2 3 ( k + D 12 d u 1 d x )
σ 33 = σ 22
σ 11 = 2 3 ( k 2 D 12 d u 1 d x )
σ 12 = D 22 d u 1 d x
< θ u 2 > = D 22 d θ d x
D 12 = 2 ϵ C 0 σ 12 ( σ 11 + σ 22 )
D 22 = 2 ϵ C 0 ( σ 12 2 + σ 22 2 )
c k d d x ( D 22 d k d x ) + P ϵ = 0
1 σ ϵ * d d x ( D 22 d ϵ d x ) + ( c ϵ 1 P c ϵ 2 ϵ ) ϵ k = 0
In the subsequent sections we shall compare the above descriptions with results of DNS.

5. Testing the Diffusion Representations by DNS

In the derivation of equations for the mean values of flow quantities use has been made of the diffusion representation of mean fluxes, i.e. the mean value of products of fluctuating quantities. Their correctness and accuracy will be assessed by comparison with DNS results published for friction Reynolds numbers R e τ = u τ H ν of 950 [15], 2000 [13] and 10 4 [14]. The comparison is focussed on the flow in the main region, the region outside the thin boundary layer at the wall. In the main region flow statistics are governed by unstable large eddies governed by inertia forces while the effect of viscosity in this region is negligibly small. The presented basic k- ϵ and fundamental model intend to describe these statistics. The boundary layer is governed by viscous forces and is situated in the region x < 100 R e τ , that is, x < 0.1 , 0.05 and 0.01 for R e τ = 950 , 2000 , and 10 4 respectively. The DNS results of the boundary layer near x = 0 are omitted in this analysis.

5.1. Diffusion of Momentum

In case of channel flow turbulent momentum transport is apparent in the equation
< u 1 u 2 > = v D N S d u 1 d x
where ν D N S is turbulent viscosity according to the DNS results, that is the value calculated from eq. (25) when substituting the values of < u 2 u 1 > and d d x u 1 obtained from DNS data. The value according to the basic k- ϵ model ν t is given by Equation (6) and that according to the fundamental model D 22 by Equation (24g) whereby the RHS’s of these equations are evaluated from the DNS data. The DNS data which were used are those of R e τ = 10 4 [14]. The three turbulent viscosities thus calculated have been shown as a function of x in Figure 1. The ratios ν t ν D N S and D 22 ν D N S versus x have been shown in Figure 2. Disregarding the thin viscous layer at the wall, the fundamental model gives satisfactory agreement over the entire x range without the use of calibration factors. Deviations are less than 10 % from the DNS results. They can be ascribed to truncation of the higher order terms in the expansions in powers of C 0 1 which were used in the theory leading to the presented expressions [8,9]. The same conclusion was arrived at when using the DNS data of R e τ = 2000 [9,13]. The basic k- ϵ model on the other hand disagrees quite a lot from DNS. The disagreement is largest at small x and gradually becomes smaller when approaching the central axis of the channel: x = 1 . The dependency on x is apparently not well captured by representation (23d), in contrast with (24g) which shows satisfactory agreement over the entire range.

5.2. Diffusion of Temperature

Temperature is a conservative quantity in incompressible flow at high Reynolds number where heat conductivity by molecular vibration is negligibly small outside thin layers at the walls. From the theory underlying the fundamental model it follows that the thermal diffusion constant equals D i k : cf. eq.(8). In case of channel flow it becomes D 22 and becomes equal to that of turbulent viscosity: cf. Equations (24d)–(24e). This result is confirmed by Lagrangian based DNS of channel flow at R e τ = 950 [15] where the thermal diffusion coefficient according to DNS, < u 2 θ > d θ d x is compared to that of the fundamentally based model D 22 .The agreement closely resembles that shown in Figure 1 and Figure 2, c.q., < u 2 u 1 > d u 1 d x and D 22 . Deviations of the predictions of the basic k- ϵ model are likely equally large as those shown in Figure 1 and Figure 2.

5.3. Diffusion of Kinetic Energy and Pressure

In the case of the basic k- ϵ model fluxes of kinetic energy and pressure are lumped together: Equation (13) For channel flow these become
< k u 2 > + < p u 2 > = ν t σ k d k d x
From eq.(14) we have for the fundamental model
< k u 2 > + < p u 2 > = c k D 22 d k d x
The description is used in the numerical solution of the equation for k in Section 6.
In Figure 3 we have shown the sum of both fluxes versus x when the RHS’s of (26) and (27) are evaluated by the DNS data. Values for the calibration constants σ k and c k of 1 and 1.3 were taken. Also is shown the value of the sums obtained by direct calculation of their values using DNS of R e τ = 10 4 . The empirical construction of Equation (27) based on the fundamental model gives surprisingly good agreement. The agreement extends over the entire x-range and is obtained with the calibration constant c k = 1.3 .

5.4. Diffusion of Energy Dissipation

DNS data do not provide information of average values of products of velocity fluctuations and dissipation fluctuations. Direct verification of the diffusion representation is thus not possible. Further, the differential equation for mean energy dissipation in both models is largely an empirical construction. What is still possible is to verify the values of quantities obtained from the solutions of the coupled differential equations for mean kinetic energy and mean energy dissipation with the DNS data. This is the subject of the next section.

6. Solutions of k and ϵ Compared with DNS

In the previous section predictions of individual components of the models were verified by DNS. Further testing is executed by comparing numerically obtained solutions of the model equations as a whole with DNS.

6.1. Equations and Boundary Conditions

Model equations are given by coupled differential equations for mean kinetic energy and mean energy dissipation. For mean energy dissipation the variable G is introduced which is defined by
G = κ ϵ x
Near the wall outside the viscous layer solutions should comply with the solutions of the log layer. Inertial subrange asymptotics [18] reveal that production P and dissipation ϵ are equal in this region: P = ϵ . Furthermore, κ d u 1 d x = 1 so that G = 1 . DNS results show values of G in the log layer which are between 0.9 and 1 [13,14]. In the present analysis the theoretical value of 1 is taken. Making use of relations (21), (23b), (23d) and (28) differential equations (23e) and (23f) can be transformed into the following coupled equations for k and G appropriate for the basic k- ϵ model.
  • Basic k-ϵ model equations:
κ 2 σ k x G d d x ( A x G d k d x ) + ( 1 x ) 2 A 1 = 0
κ 2 σ ϵ x 2 G 2 d d x ( A x G d d x ( G x ) ) + 1 k ( c ϵ 1 ( 1 x ) 2 A c ϵ 2 ) = 0
where
A = c μ k 2
where values of the empirical constants will be taken in agreement with [5,6]: κ = 0.4 , c μ = 0.09 , σ k = 1 , σ ϵ = 1.3 , c ϵ 1 = 1.49 , c ϵ 2 = 1.9 . The boundary conditions are
x = 0 ; k = 1 c μ ; G = 1
x = 1 ; d k d x = 0 ; d G d x = G
The boundary condition for k at x = 0 follows from relations (21) (23b) and (23d) noting that ϵ = d u 1 d x in the log layer. The boundary condition for G at x = 1 corresponds to zero slope of ϵ at x = 1 .
The equations appropriate for the fundamental model follow from (24h) and (24i) upon using (28), (21), (24a)–(24g).
  • Fundamental model equations:
c k κ 2 x G d d x ( B x G d k d x ) + ( 1 x ) 2 B 1 = 0
κ 2 x 2 k G 2 σ ϵ * d d x B x G d G x d x + c ϵ 1 ( 1 x ) 2 B c ϵ 2 = 0
and
B = 2 C 0 ( σ 12 2 + σ 22 2 ) = 2 C 0 ( ( 1 x ) 2 + σ 22 2 )
In these equations k and σ 22 are related to each other by an algebraic relation which is obtained by systematic elimination of ϵ 1 d u 1 d x , σ 12 , σ 11 , σ 33 , D 12 , D 22 from Equations (21), (24a) and (24g)
k = σ 22 3 σ 22 2 + ( 1 x ) 2 2 ( σ 22 2 ( 1 x ) 2 )
The relation can be used to eliminate σ 22 from (34)–(35) resulting in two differential equations for k and G. However, analysis of the above relation for k shows that in an area close to x = 0 according to (37) σ 22 can have two values for one value of k. Verification of (37) using the results of DNS confirms the correctness of the equation and the double dependency in this area. The other way around is not the case: for each value of σ 22 there is only one value of k possible. The easiest approach is then to substitute for k relation (37) into the two differential equations (34) and (35) and solve for σ 22 and G. The value of k is subsequently obtained from (37). The boundary conditions are
x = 0 ; G = 1 ; σ 22 = C 0 2 1
x = 1 ; d σ 22 d x = 0 ; d G d x = G
where the boundary condition for σ 22 at x = 0 follows from Equations (21), (24d) and (24g) noting that ϵ = d u 1 d x in the log layer. As in the basic k- ϵ model c ϵ 2 = 1.9 . A value of c k of 1.3 was established in Section 5.3 (Figure 3). A value of σ ϵ * of 0.2 is found to lead to the best agreement with DNS.
The value of c ϵ 1 follows from Equation (35) by letting x approach x = 0 . It yields the relation
c ϵ 1 = c ϵ 2 k 0 κ 2 σ ϵ *
which is analogous to eq.(17) in case of the basic k- ϵ model; k 0 is the value of k at x = 0 obtained from eq.(37). For c ϵ 2 = 1.9 , C 0 = 7 , σ ϵ * = 0.2 and k 0 = 4.48 , c ϵ 1 = 1.68 . This suggests negative production. However, for the first term on the LHS of eq.(35) one can write
x 2 d d x ( B x G d G x d x ) = x 2 d d x ( B G d G d x ) x d B d x + B
where the last two terms have the same character as that of production and make the total of production-like-terms positive.

6.2. Numerical Solution

Numerical instability is encountered when solving Equations (29)-(30) and (34)-(35), a feature not uncommon for k- ϵ equations: Lew et all [19]. Way out is to convert the equations into a diffusion problem by adding the terms d k d t and d G d t on the RHS of (29)-(30), d σ 22 d t and d G d t on the RHS of (34)-(35) and solving starting from a suitably chosen initial solution at t = 0 : Borse [20]. After sufficient time the solution has converged to the desired stationary result.
Another difficulty encountered in the execution of the numerical calculations concerned the boundary condition for G at x = 1 : cf. Equations (33) and (39). These were replaced by ( d / d x ) G = 0 at x = 1 . The effect of this simplification is limited to a small region near x = 1 .
In Figure 4 and Figure 5 are respectively shown the distributions versus x of kinetic energy k and energy dissipation rate G according DNS, according the basic k- ϵ model calculated from Equations (29) - (30) and according to the fundamentally based model calculated from Equations (34) - (35).

6.3. Analytical Solution

Distinction is made between the area near the wall referred to as outer region and the area near the centre of the channel referred to as inner region. The outer region is the area where production and dissipation of energy are dominant. Turbulent diffusion of kinetic energy and pressure and turbulent diffusion of dissipation are negligibly small. It is the area where the log layer description for mean flow applies and where production equals dissipation [21]. In the inner region turbulent diffusion of kinetic energy and pressure and turbulent diffusion of dissipation are important and in balance with energy dissipation while production is negligibly small. For each of these regions analytical solutions can be derived which are subsequently matched to arrive at a complete solution.

6.3.1. Solutions for k and G in the Outer Region

Solutions valid in the outer region are obtained by retaining the second and third term in eqs.(28) and eqs.(33) yielding A = ( 1 x ) 2 and B = ( 1 x ) 2 or in terms of k and σ 22
  • Basic k-ϵ model:
k = 1 c μ ( 1 x )
Fundamental model:
σ 22 = ( C 0 2 1 ) ( 1 x )
Using the relation between k and σ 22 at x = 0 according to the fundamental model, cf. Equation (37), we have
  • Fundamental model:
k = k 0 ( 1 x )
where k 0 is the value of k at x = 0
k 0 = 1 2 ( C 0 2 1 ) 3 C 0 4 C 0 4
which amounts to k 0 = 4.48 for C 0 = 7 . To derive the equations for G applicable in the outer region the second and third term of Equation (41) need to be taken into account (and replacing B by A in Equation (41) in case of the basic k- ϵ model). Noting that in the outer region A = ( 1 x ) 2 and B = ( 1 x ) 2 one obtains for G from Equation (30) using relation (17) and from Equation (35) using relation (40) the result
  • Basic k-ϵ model and Fundamental model:
G = ( 1 x ) ( 1 + x )
The above descriptions for k and G are rather simple, They compare well with the numerical results for values of x up to about 0.5 . At greater distances from the wall they start to deviate and fail to meet the conditions of zero slope at x = 1 . To overcome this deficiency solutions valid for the inner region are developed.

6.3.2. Solutions for k in the Inner Region

Solutions for k are developed by describing k by a series of successive powers of η = 1 x . The linear term in η is omitted in order to satisfy the zero-slope condition at η = 0 . Disregarding terms of O ( η 3 ) and higher we write for k
k = k 1 + a η 2
where k 1 is the value of k at x = 1 ; k 1 and a are to be determined. Substituting expression (47) in Equation (29) and (34) and equating the leading terms of O ( η 0 ) we obtain for a in case of the basic k- ϵ and fundamental models respectively the relations
2 a c μ κ 2 k 1 2 σ k = G 1 2
and
16 a c k κ 2 k 1 2 9 C 0 = G 1 2
where G 1 is the value of G at η = 0 . In deriving result (49) we took for σ 22 at η = 0 the value 2 3 k 1 in accordance with Equation (37). The boundary between inner and outer region is defined by η 0 = 1 x 0 . At this boundary the value of k and its slope should match the values according to those of the outer region. This yields the relations:
k 1 + a η 0 2 = η 0 c μ ; k 1 + a η 0 2 = η 0 k 0
2 a η 0 = 1 c μ ; 2 a η 0 = k 0
for basic k- ϵ and fundamental model respectively. From Equations (47)-(51) the following solutions are obtained:
  • For the basic k- ϵ model:
η 0 = 4 σ k G 1 2 c μ κ 2 ; k 1 = η 0 2 c μ ; a = 1 2 c μ η 0
and for the fundamental model
η 0 = 9 G 1 2 C 0 2 κ 2 k 0 3 c k ; k 1 = k 0 η 0 2 ; a = k 0 2 η 0
The above solution parameters can be determined by implementing the values of the system parameters σ k , c μ , κ , C 0 , k 0 and c k . The value of G 1 is taken from the numerical results: G 1 = 0.25 in case of the basic k- ϵ model and G 1 = 0.44 in case of the fundamental model. An alternative approach to determine G 1 is to couple the above analytical solution of k to that of G presented in Section 6.3.3 below. It yields values for G 1 without recourse to the numerical results. The values are practically equal to those of the numerical results.
A complete analytical description of k for outer and inner region is obtained from Equations (41) and (44) when x x 0 and Equations (47) and (52) when x x 0 where x 0 = 0.47 in case of the basic k- ϵ model and x 0 = 0.33 in case of the fundamental model. The result has been shown in Figure 4 by solid blue and red lines.

6.3.3. Solutions for G in the Inner Region

For G in the inner region the following expansion is used
G = G 1 ( 1 η ) + a η 2 .
where the linear term G 1 η is introduced to satisfy the boundary condition at x = 1 : cf. Equations (33) and (39). The values of G 1 and a are to be determined. Substituting description (54) into Equations (30) and (35) and equating the leading terms of O ( η 0 ) yields
a = α G 1 3
where
α = σ ϵ c ϵ 2 2 κ 2 c μ k 1 3
and
α = 9 C 0 σ ϵ * c ϵ 2 16 κ 2 k 1 3
for basic k- ϵ model and fundamental model respectively. At the boundary of inner and outer region η = η 1 the values of G and its slope according to outer and inner regions should become equal. It yields the relations
η 1 ( 2 η 1 ) = ( 1 η 1 ) G 1 + α G 1 2 η 1 2
( 2 3 2 η 1 ) 2 η 1 = G 1 + 2 α G 1 2 η 1
Eliminating G 1 from the above equations yields the following irreducible equation for η 1
α = 2 η 1 4 ( 1 η 1 2 ) 3 ( 1 η 1 4 ) 3 ( 1 3 η 1 4 + η 1 2 4 )
The value of α follows from Equation (56) using the system values for the various parameters and taking for k 1 the value of the numerical calculations, i.e. k 1 = 0.8 , with the result that α = 134 and α = 18 in case of the basic k- ϵ model and the fundamental model respectively. By iteration one finds from eqs (60) for the basic k- ϵ model η 1 = 0.31 ( x 1 = 0.69 ) and for the fundamental model η 1 = 0.47 ( x 1 = 0.53 ) . A complete analytical description of G for outer and inner region is obtained from Equation (46) when x x 1 and Equations (54)-(56) when x x 1 . The result has been shown in Figure 5 by solid blue and red lines.

6.4. Discussion of Results

Numerical and analytical solutions for kinetic energy k and energy dissipation ϵ , the latter in terms of G = κ x ϵ , have been developed for the basic k- ϵ model and a new fundamentally based model. The solutions were compared with DNS results and have been shown in Figure 4 and Figure 5. Conclusions are:
  • Analytical solutions agree in a satisfactory manner with numerical solutions. The analytical solutions reveal the relative contributions of turbulent diffusion, energy production and energy dissipation in outer and inner region of the channel. They show to what extent the descriptions are empirically or fundamentally based and depend on calibration factors.
  • Solutions for k in the in the outer region ( x 0.33 ) according to the fundamental model do not depend on calibration constants. They have a fundamental basis. The differences with DNS can be ascribed to errors as a result of truncation in the expansion with respect to C 0 1 which underlies the fundamental model.
  • Solutions for k in the inner region ( x 0.33 ) according to the fundamental model depend on the calibration factor c k .
  • Solutions for k in both the outer and inner region according to the basic k- ϵ model depend on the calibration factors c μ and σ k .
  • Solutions for G in the outer region according to the fundamental model and the basic k- ϵ model, ( x 0.53 ) and ( x 0.69 ) do not depend on calibration factors. Differences with DNS are due to some deviation between the value of production and dissipation in this area.
  • Solutions for G in the inner region according to the fundamental model and the basic k- ϵ model, ( x 0.69 ) and ( x 0.53 ) respectively, depend on the calibration factors σ ϵ * and σ ϵ .
  • Using standard calibration constants in the solutions of the basic k- ϵ model results in notable deviations compared to DNS. The deviations can be reduced by recalibrating σ ϵ , c μ and σ k . Deviations between diffusion constants remain significant because of different functional dependencies: Figure 1, Figure 2 and Figure 3.

7. Velocity Distributions

The basic k- ϵ model is based on the assumption of an isotropic turbulence field. It cannot predict anisotropic behaviour of co-variances or turbulent stresses. In case of the fundamental model the distribution of σ 22 follows from the solution of eqs.(34) According to Equation (24b) σ 33 = σ 22 while σ 11 = 1 2 k σ 22 where k is given by Equation (37). In Figure 6 the distribution versus x is shown of the calculated values of the standard deviations of the fluctuations σ 11 , σ 22 and σ 33 according to the fundamental model and the results of DNS. Deviations are similar to those shown in [9] and are ascribed to truncation of the expansions of C 0 1 which underlies the theory of the fundamental model. Deviations of the three variances σ 11 , σ 22 , σ 33 from k which describe anisotropy are next to leading order in the C 0 1 expansion (cf. Equation (24a) - (24c)) and therefore more sensitive to truncation error.
For the gradient of the mean flow the following relations can be derived. For the basic k- ϵ model from eqs (21) and (23b)
d u 1 d x = 1 x ν t
For the fundamental model from eqs (21) and (24d)
d u 1 d x = 1 x D 22
In Figure 7 the distributions versus x are shown of κ x d u 1 d x according to DNS, the basic k- ϵ model and the fundamental model. The deviations between DNS and basic k- ϵ model are similar to those shown for ν D N S and ν t shown in Figure 1 and Figure 2. Similarly the agreement between ν D N S and D 22 shown in Figure 1 and Figure 2. According to (21): P = ( 1 x ) d u 1 d x which indicates that disagreement of production with DNS in case of the basic k- ϵ model is similar to that shown in Figure 7 and that agreement of production with DNS in case of the fundamental model is similar to that shown in Figure 7.
The distribution of mean velocity u 1 predicted by the basic k- ϵ model and the fundamental model are obtained by integration of the RHS’s of Equations (61) and (62). Because the models do not describe the velocity in the viscous layer at the wall the integration starts at some distance from the wall for which the position x = 100 R e τ is taken. The value of u 1 at this point according to the DNS results of R e τ = 10 4 is 17.2 which agrees with the values obtained from measurements: Monin and Yaglom, Vol. I, fig. 25 [21]. The results of the integration are shown in Figure 8, where the distributions of mean velocity of the two models and of DNS are presented. It is seen that for equal shear velocity, that is for equal longitudinal pressure gradient in the channel, the fundamental model somewhat overestimates mean velocity, by 5 % while the basic k- ϵ model underestimates by 18 % .

8. Conclusions

An analysis has been presented of two models for determining mean values of turbulent flow variables: velocities u i , pressure p, turbulent diffusion D i j , kinetic energy k and energy dissipation ϵ . The two models are the basic k- ϵ model widely used in engineering applications and environmental analysis [5,6], and a new fundamentally based model which is derived from recently published statistical descriptions of inhomogeneous anisotropic turbulence [7,8,9]. The analysis was focused on turbulent channel flow which is highly anisotropic and inhomogeneous. It is of direct relevance to applications of duct flows, turbulent boundary layers around bodies and the atmospheric surface layer around the earth. Model predictions were verified against published results of direct numerical simulations DNS for shear Reynolds numbers of 950 [15], 2000 [13] and 10 4 [14]. The comparison was focussed on the regions outside the thin viscous boundary layers of a few millimeters at the wall. The development status of DNS is mature; its results can serve as a reliable source for verifying model predictions.
The basic k- ϵ model is built on an empirical isentropic representation of the turbulence field. It contains a gradient hypothesis for turbulent diffusion which is supplemented by calibration constants. Using conventionally proposed values for the calibration constants [5,6] model predictions of several variables are found to deviate remarkably from the results of DNS. A summary has been given in Table 1.
The fundamental model is derived from a theory which has a fundamental basis [7,8,9]. Expressions for turbulent diffusion of momentum and of conservative scalars such as temperature and passive admixture (smoke, aerosols) are of a general nature. They are free from calibration constants. They are found to agree in a satisfactory manner with DNS results at all distances x from the wall (outside the thin boundary layer at the wall). The theory underlying the fundamental model does not provide generally valid expressions for turbulent diffusion of the non-conservative scalars kinetic energy, pressure, energy dissipation. The non-conservative diffusion terms are only important for the determination of k and ϵ in the inner half of the channel. Here non-conservative behaviour has been modeled by using the diffusion expressions of conservative scalars provided with calibration factors. A summary is given in Table 1.
A general limitation of empirical models is that they lack a fundamental basis. It is uncertain whether the modelled terms are correctly described by their dependencies on variables. The addition of calibration factors tempers inaccuracies of prediction to some extent. At the same time they need to be specified for each new case. Models which are derived from generally valid principles do not suffer from these limitations. The presented fundamental model contains descriptions of several variables which satisfy this criterion. Calibration factors are absent and applicability can be expected to exceed that of channel flow. What is still missing in the fundamental model are fundamental expressions for non-conservative scalars. Their development eliminates the remaining empirical factors.

Funding

This research received no external funding.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

B.G.J. Ruis and H.S. Janssen are acknowledged for performing numerical calculations; G.M. Janssen for preparing the manuscript.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Boussinesq, J. Mem.pres.par.div. savantis a l’acad. sci. de Paris, 1877, 23, 1, 1-680, 1877, 24, 2, 726-736.
  2. Taylor, G.I. Phil.Trans. 1915, A215, 1-26, and Proc.Roy.Soc.(London), 1932, A135, 685-701.
  3. Prandtl, L. Zeits. f. Math.u. Mech., 1925, 5, 136-139.
  4. Karman, Th. von. Nachr. Ges. Wiss. Goettingen, Math. Phys., 1930, Klasse 58-76.
  5. Hanjali´c, K.; Launder, B. Modelling Turbulence in Engineering and the Environment: Second-Moment Routes to Closure; Cambridge University Press: Cambridge, UK, 2011. [Google Scholar]
  6. Bernard, P.S.; Wallace, J.K. Turbulent Flow: Analysis, Measurement and Prediction; Wiley: New Jersey, NJ, USA, 2002. [Google Scholar]
  7. Brouwers, J.J.H. Statistical description of turbulent dispersion. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 2012, 86, 066309. [Google Scholar] [CrossRef] [PubMed]
  8. Brouwers, J.J.H. Statistical Models of Large Scale Turbulent Flow. Flow Turbul. Combust. 2016, 97, 369–399. [Google Scholar] [CrossRef]
  9. Brouwers, J.J.H. Statistical Descriptions of Inhomogeneous fundamentally based Turbulence Mathematics 2022, 10, 4619.
  10. Monin, A.S.; Yaglom, A.M. Statistical Fluid Mechanics, Vol II, Ch 8; Dover: New York, NY, USA, 2007. [Google Scholar]
  11. van Kampen, N.G. Stochastic Processes in Physics and Chemistry, 3rd ed.; Elsevier: New York, NY, USA, 2007. [Google Scholar]
  12. Stratonovich, R.L. Topics in the Theory of Random Noise; Gordon and Breach: New York, NY, USA, 1967; Volume 1. [Google Scholar]
  13. Hoyas, S.; Jiménez, J. Scaling of the velocity fluctuations in turbulent channels up to Re_τ=2003. Phys. Fluids 2006, 18, 011702. [Google Scholar] [CrossRef]
  14. Hoyas, S.; Oberlack, M.; Alcántara-Ávila, F.; Kraheberger, S.V.; Laux, J. Wall turbulence at high friction Reynolds numbers. Phys.Rev. Fluids 2022, 7, 014602. [Google Scholar] [CrossRef]
  15. Kuerten, J.; Brouwers, J.J.H. Lagrangian statistics of turbulent channel flow at Re_τ=950 calculated with direct numerical simulation and Langevin models. Phys. Fluids 2013, 25, 105108. [Google Scholar] [CrossRef]
  16. Kolmogorov, A. The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds’ Numbers. Akad.Nauk. Sssr Dokl. 1941, 30, 301–305. [Google Scholar]
  17. George, W.K. The Decay of Homogeneous Isotropic Turbulence. Phys.Fluids A 1992, 4, 1492–1509. [Google Scholar] [CrossRef]
  18. Brouwers, J.J.H. Dissipation equals Production in the Log Layer of Wall Induced Turbulence. Phys.Fluids 2007, 19(10), 101702. [Google Scholar] [CrossRef]
  19. Lew,A.J.; Gustavo, C.; Carrica,B.; Carrida, P.M. Int. J. Computational Fluid Dynamics. 2001,14:3, 201-209. [CrossRef]
  20. Borse,G.J. Numerical Methods with Matlab, Ch.22, PWS Publishing, Boston, USA, 1997. ISBN 0534938221.
  21. Monin, A.S.; Yaglom, A.M. Statistical Fluid Mechanics, Vol I. Dover: New York, NY, USA, 2007.
Figure 1. Turbulent viscosities according to ν D N S at R e τ = 10 4 , fundamental model D 22 and basic k- ϵ model ν t versus dimensionless distance from the wall x. Results from DNS are represented by the solid line, from the fundamental model by dashdot line and from the basic k- ϵ model by dashed line.
Figure 1. Turbulent viscosities according to ν D N S at R e τ = 10 4 , fundamental model D 22 and basic k- ϵ model ν t versus dimensionless distance from the wall x. Results from DNS are represented by the solid line, from the fundamental model by dashdot line and from the basic k- ϵ model by dashed line.
Preprints 100397 g001
Figure 2. Ratio of turbulent viscosity of the fundamental model to that of DNS at R e τ = 10 4 D 22 ν D N S (lower line) and ration of turbulent viscosity of the basic k- ϵ model to that of DNS at R e τ = 10 4 ν t ν D N S (upper line) versus dimensionless distance from the wall x.
Figure 2. Ratio of turbulent viscosity of the fundamental model to that of DNS at R e τ = 10 4 D 22 ν D N S (lower line) and ration of turbulent viscosity of the basic k- ϵ model to that of DNS at R e τ = 10 4 ν t ν D N S (upper line) versus dimensionless distance from the wall x.
Preprints 100397 g002
Figure 3. Sum of turbulent fluxes of kinetic energy and pressure versus dimensionless distance from the wall x. The solid line is the sum according to DNS at R e τ = 10 4 , the dashdot line according to an empirical construction using the fundamental model and the dashed line the empirical construction of the basic k- ϵ model.
Figure 3. Sum of turbulent fluxes of kinetic energy and pressure versus dimensionless distance from the wall x. The solid line is the sum according to DNS at R e τ = 10 4 , the dashdot line according to an empirical construction using the fundamental model and the dashed line the empirical construction of the basic k- ϵ model.
Preprints 100397 g003
Figure 4. Distribution versus x of kinetic energy k according to DNS (solid line), according to the numerically and analytically assessed basic k- ϵ model ( dashed and full line respectively), and according to the numerically and analytically assessed fundamental model (dashdot and full line respectively).
Figure 4. Distribution versus x of kinetic energy k according to DNS (solid line), according to the numerically and analytically assessed basic k- ϵ model ( dashed and full line respectively), and according to the numerically and analytically assessed fundamental model (dashdot and full line respectively).
Preprints 100397 g004
Figure 5. Distribution versus x of kinetic energy G according to DNS (solid line), according to the numerically and analytically assessed basic k- ϵ model ( dashed and full line respectively), and according to the numerically and analytically assessed fundamental model (dashdot and full line respectively).
Figure 5. Distribution versus x of kinetic energy G according to DNS (solid line), according to the numerically and analytically assessed basic k- ϵ model ( dashed and full line respectively), and according to the numerically and analytically assessed fundamental model (dashdot and full line respectively).
Preprints 100397 g005
Figure 6. Distribution versus x of the standard deviations of fluctuating velocities σ 11 , σ 22 and σ 33 according to DNS (solid line) and the fundamental model (dashdot line).
Figure 6. Distribution versus x of the standard deviations of fluctuating velocities σ 11 , σ 22 and σ 33 according to DNS (solid line) and the fundamental model (dashdot line).
Preprints 100397 g006
Figure 7. Distribution versus x of κ x d u 1 d x according to DNS (solid line), the basic k- ϵ model (dashdot line)and the fundamental model ( dashed line).
Figure 7. Distribution versus x of κ x d u 1 d x according to DNS (solid line), the basic k- ϵ model (dashdot line)and the fundamental model ( dashed line).
Preprints 100397 g007
Figure 8. Distribution of mean velocity divided by shear velocity versus x according to DNS (solid line), the basic k- ϵ model (dashdot line)and the fundamental model ( dashed line).
Figure 8. Distribution of mean velocity divided by shear velocity versus x according to DNS (solid line), the basic k- ϵ model (dashdot line)and the fundamental model ( dashed line).
Preprints 100397 g008
Table 1. Predictions of the basic k- ϵ model and the fundamentally based model compared with DNS of turbulent channel flow
Table 1. Predictions of the basic k- ϵ model and the fundamentally based model compared with DNS of turbulent channel flow
variable basic k- ϵ model fundamentally based model
Turbulent viscosity Empirical Equation (23d) Significant deviations Figure 1 and Figure 2 Fundamentally based Equation (24d) Satisfactory agreement: Figure 1 and Figure 2
Turbulent diffusion of temperatures, smokes, aerosols Empirical Equation (23c) Significant deviations similar to Figure 1 and Figure 2 Fundamentally based Equation (24e) Satisfactory agreement similar to Figure 1 and Figure 2
Tubulent diffusion of kinetic energy and pressure Empirical Equation (26) Deviations Figure 3 Empirical Equation (27) Satisfactory agreement Figure 3
Mean value of kinetic energy Empirical Equation (29) Deviations Figure 4 In the outer half of the channel: fundamentally based Equation (55). In the inner half: empirical Equation (55). Satisfactory agreement Figure 4
Mean value of energy dissipation rate In the outer half of the channel: fundamentally based Equation (46) In the inner half: empirical  Equation (55). Deviations: Figure 5 In the outer half of the channel: fundamentally based Equation (46). In the inner half: empirical Equation (55). Satisfactory agreement: Figure 5
RMS values of fluctuations No prediction Qualitative agreement Figure 6
Mean value of velocity Empirical Equation (61) Deviations Figure 7 and Figure 8 Fundamentally based Equation (62) Satisfactory agreement Figure 7 and Figure 8
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated