Theoretical Study of Fully Developed Turbulent Flow in a Channel, Using Prandtl’s Mixing Length Model

This research used the common decomposition of the velocity and pressure in an average part and a fluctuating part, for high Reynolds number, of the Navier-Stokes equation, which leads to the classic problem of turbulent closure. The Prandtl’s mixing length model, based on the Boussinesq hypothesis and traditionally used for free shear flows, was chosen and adapted for internal flows to solve the closure problem. For channel flows, Johann Nikuradse proposed a model for the Prandtl mixing length. In the present paper, which has an academic character, the authors made a return to the model of the mixing length of Prandtl and the model of Nikuradse, to solve turbulent flows inside a plane channel. It was possible to develop an ordinary differential model for the velocity in the direction of the flow whose solution occurs computationally in a simple but extremely accurate way when compared with Direct Numerical Simulation databases. For the viscous stress on the wall, it was possible to determine the exact mathematical solution of the ordinary differential equation. It is a model of great academic value and even to be used as reference for verification of computational codes destined to the solution of complete numerical and computational models.


Introduction
Turbulent flows are more common in nature and in industrial applications when compared to laminar flows. It is possible to determine exact and numerical solutions for laminar flows. However, for transition to turbulence and for fully tur-but none, without many limitations, Schlichting (1968) [5], Nikuradse (1966) [6]. The studies by material experimentation are limited by the constructive difficulties of the experimental apparatus, as well as by the difficulties of instrumentation for the development of the measures. Computational experiments are limited by the methodologies used to solve differential models, as well as by the processing and storage power of computers. Finally, it is important to comment on the methods of exact solution, which are feasible only for some canonical flows, due to the complexity of the differential mathematical models associated to the turbulent flows.
Studies involving exact solutions of canonical turbulent flows, such as, free shear flows, can be found in the literature, Schlichting [5] and Wilcox [7]. Also, some works can be found about turbulent flows in flat channels, for instance, Nikuradse [6].
Numerous works are found about turbulence closure modeling, such as the Prandtl mixing length model [3]. In the present work, the turbulence closure problem was modeled using the mixing length cited by Prandtl [3] as well as the proposal of Nikuradse [6] for modeling the characteristic length of turbulence for flows in flat channels and in circular ducts. Thus, the main objective of the present work is related to the semi-analytical study of turbulent flows in flat channels. We present the differential model for the mean velocity field as well as the closure model based on the Prandtl and Nikuradse proposals. The model of the mixing length of Prandtl [3] was, in the present work, modified, using damping functions, to consider the presence of walls. With these functions, it was ensured that the effects of turbulence near walls were damped to zero, as expected from the physical point of view. The solution of the differential model for the mean velocity was obtained by numerical and computational methods.
The solution to the viscous stress on the wall was obtained exactly.
It is important to emphasize that for internal flows, in the near regions of walls, the viscous molecular effects grow dramatically as it approaches the walls, which dampens the effects of the turbulence. The model of the mixing length of Prandtl, by itself, does not correctly model this physical characteristic of the internal flows. Nikuradse [6] proposed a polynomial function for evaluation of the mixing length along the whole channel, including the wall regions. However, as will be shown in the results chapter, this polynomial function is not enough to Journal of Applied Mathematics and Physics correctly model the problem. In the present work, it was proposed to use an additional damping function proposed by Van Driest (1956) [8]. It will be shown that the use of this additional damping function of the mixing length of Prandtl, promotes an additional damping, apparently small, but that presents very important consequences in the prediction of average flow over the whole channel.

Physical Model
Poiseuille laminar and turbulent flows in flat and circular channels are among the first flows that were solved by the scientific community. Initially this type of flow was solved in laminar regime by Poiseuille (1846) [9]. Later, it was modeled, in turbulent regime, by Prandtl [3] and by Nikuradse [6]. The physical model for the flat Poiseuille flow is shown in Figure 1, where the walls, the entrance and the exit, the coordinate axes system ( ) , , x y z , the half-width of the channel, R, the average speed profile, ( ) u y and the axis of symmetry is illustrated. A Poiseuille flow is characterized by being developed in the x-axis, having The flow is maintained with a pressure gradient p x ∂ ∂ .
On average, the flow is considered in permanent regime. The fluid is Newtonian, and the flow is incompressible. This is the physical basis that will be used to establish the differential mathematical model, presented in the following item.

Differential Mathematical Model
Laminar and turbulent, isothermal flows of Newtonian fluids can be modeled by the continuity and Navier-Stokes equations. These equations, for incompressible flows, when submitted to the mean temporal operator, are written below, using indicial notation: Equation (1) represents the mean mass balance over a fluid particle, while Equation (2) represents the linear momentum balance over the same fluid It is observed that the Equation (4)

Boussinesq Hypothesis
The decomposition of velocities and pressure leads to the appearance of the Reynolds stress tensor, which, for the three-dimensional formulation, is given by Equation (5): .
Boussinesq proposed to model this tensor using his hypothesis, which consists of making an analogy to the Stokes model for the molecular viscous stress tensor.
The model is represented by Equation (6), for one-dimensional flows, as is the case of the flow analyzed in the present work: is called turbulent kinetic energy, xy δ is the Cronecker's operator, and t ν was defined as the kinematic turbulent viscosity. While molecular viscosity is a physical property of fluids, turbulent viscosity is a property of the flow. Substituting this model into Equation (4) and simplifying terms, results in Equation (5): The turbulent viscosity modeling is presented in the following section.

Prandtl Mixing Length Model
The turbulent viscosity that appears in Equation (7) remains to be modeled and calculated. One of the first proposals for this was made by Prandtl [3]. This proposal was made for parallel flows, which is given by Equation (8) Substituting Equation (8) into Equation (7) gives equation (9): Analyzing Figure 1 we conclude that: Equation (10) is valid for any value of the coordinate 0 y R ≤ ≤ . Substituting Equation (10) into Equation (9) we obtain Equation (11), which is written in a form suitable to be solved.

Variable Mixing Length without Wall Function
Even though linear m l gives a satisfactory result for free-shear flows, like mixing layer, for internal flows, inconsistent results are generated. Since the molecular viscous effects are more relevant near the wall than the turbulent effects, Nikuradse proposed a polynomial empirical equation in [6], Rewritten by Equation (12). This model gives a mixing length model that varies with y coordinate, zeroing near the walls.
Now, Equation (11) can be rewritten using Equation (12): Considering that Equation (13) is a second order differentiation equation, its solution requires two boundary conditions. On the walls the non-slip condition must be applied by Equation (14). In the center of the channel, a second specie boundary condition of Neumann, that is, a condition of symmetry was used, given by Equation (15).
Integrating Equation (13) we have: Journal of Applied Mathematics and Physics ( ) Applying the boundary condition given by Equation (15) we obtain: Now, rearranging Equation (16), we have the following second order poly- The Equation (20) is a first order, homogeneous, non-linear ordinary differential equation, with no exact solution. Therefore, a numerical tool is required, and the software MATLAB was used to resolve this equation, using the fourth order Runge-Kutta method already implemented in the software as a function, and well explained in Cheever [11]. For an illustration of the results that are obtained, the following values of the parameters were chosen for convenience and to facilitate solving:  where w τ is the shear stress on the wall, which is illustrated in Figure 3. It is possible to balance the forces around the control volume, to find the value of w τ . Using the terms shown in Figure 3, and taken into account the forces, we obtain: Taken lim 0 x δ → in Equation (25) and using Equation (24), we obtain: Using these parameters, Equation (20) Multiplying both sides of Equation (27) by R u τ we obtain the final nondimensional form: The new nondimensional boundary conditions are given by Equation (29) and Equation (30). Reminding that equation (30) was already applied.

Results without Damping Function
To validate the method presented in the present work, DNS data already existent in the literature, was used. The results to be compared are from DNS simulations of a turbulent channel done by JHTDB [12] and ICES Turbulence Database [13].
Beginning with 180 Re τ = and 395 Re τ = , the mean velocity profiles are plotted in Figure 4. The results are compared with DNS results. We see that the results obtained in the present work do not match with the reference results.
Results for 590 Re τ = and 1000 Re τ = are shown in Figure 5.  That happens because the method used in the present paper is for completely turbulent flows. It is not valid for laminar and transitional flows. Thus, the lower is the Reynolds values, worst is the precision of the method. However, even with high Reynolds number values the results do not agree with DNS data base, since there is still a big difference between the presented curves and DNS curves.

Results with Damping Function
As previously stated, the Prandtl mixing length should tend to zero near the walls, since the viscous effects are more relevant than the turbulent effects. The mixing length function given by Equation (12), effectively zeroes on the walls, but not as fast as it should. This requires the application of an additional damping function. Van Driest (1956) [8] proposed a wall damping function, which is presented here by Equation (32).
The results are presented in Figure 8.
The euclidean norm ( 2 L ) and infinity norm ( L ∞ ) was also calculated for a quantitative comparison between the no damping and damping curve, to the DNS. The results are presented in Table 1 and Table 2.
With the damping equation the profiles and norm practically coincide with the DNS. As presented, for the mixing length without wall function, the higher the Reynolds value, the higher the precision of the method. In the case 180 Re τ = , the D Re is only 6000, which characterize a low Reynolds number turbulence, being out of limit for the use of this method. With the correction of the damping function, the results fit better with the turbulent flow characteristics.
In Figure 10 the curves for 590 Re τ = (a) and 1000 Re τ = (b) are plotted, we see that the agreement is still better.
For these case the euclidean norm is also calculated, as well as the infinity norm ( L ∞ ), obtaining Table 3 for 2 L and Table 4 for L ∞ .
The same analysis performed for the mean velocity profile can be performed here, since it is possible to see that increasing the Reynolds number, the precision of the method is also increased. The closer to full turbulent flow, the closer is the results of the present work as compared with DNS data.
For every Reynolds number, it was possible to notice that the Reynolds stress had the expected behavior. Increasing the Reynolds, more turbulent it became, and consequently higher turbulent stress are obtained. Therefore, as expected, the higher the Reynolds number, the higher the peak of τ. Also, as the flow gets more turbulent, the velocity profile flattens, due to the momentum exchange in the y direction. Thus, the viscous region near the wall decreases, and as shown on the plots, the turbulent Reynolds stress peak, moves closer to the wall, as expected.

Conclusions
In the present work, one of the first proposals of turbulence closure modeling was used: the mixing length model of Prandtl to solve turbulent flows in a plane channel. The Nikuradse [6] model for Prandtl mixing length was used. Furthermore, the proposals of Van Driest [8] and Cebeci and Bradshaw (1984) [14] were also used to model the required damping near the walls.
A very simple model, composed of an ordinary differential equation, was proposed for the mean velocity of the flow as well as for the viscous stress on the channel walls. The differential equation for the mean velocity was solved using the Runge-Kutta's integration method. However, the differential equation for