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 turbulent flows, there are no exact solutions, except for some cases of canonical flows. Due to the high complexity and, at the same time, its high practical importance, this subject has been the object of great efforts of the scientific community around the world.
Due to the great academic and industrial interest related to this theme, many concepts and methods were developed to model the average characteristic information of turbulent flows, incorporating the effects of the turbulence itself, Boussinesq (1877)  , Reynolds (1884)  , Prandtl (1925)  and Kolmogorov (1941)  . Numerous theoretical and experimental studies have been developed, but none, without many limitations, Schlichting (1968)  , Nikuradse (1966)  . 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  and Wilcox  . Also, some works can be found about turbulent flows in flat channels, for instance, Nikuradse  .
Numerous works are found about turbulence closure modeling, such as the Prandtl mixing length model  . In the present work, the turbulence closure problem was modeled using the mixing length cited by Prandtl  as well as the proposal of Nikuradse  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  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  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 correctly model the problem. In the present work, it was proposed to use an additional damping function proposed by Van Driest (1956)  . 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.
2. Plane Channel: Physical and Mathematical Models
2.1. 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)  . Later, it was modeled, in turbulent regime, by Prandtl  and by Nikuradse  . 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 , the half-width of the channel, R, the average speed profile, and the axis of symmetry is illustrated. A Poiseuille flow is characterized by being developed in the x-axis, having components . The flow is maintained with a pressure gradient . 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.
2.2. 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
Figure 1. Schematic illustration of the Poiseuille flow between parallel plates.
particle. These differential equations can be simplified by the assumption of Poiseuille flows, i.e.: , and . Finally, the average flow is considered to be in steady state. Thus, Equation (3) and Equation (4) are obtained:
It is observed that the Equation (4) presents an additional term that is not known. This term makes the model mathematically open. This fact requires an additional equation for the closure of the model.
2.3. 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:
where is called turbulent kinetic energy, is the Cronecker’s operator, and 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.
2.4. 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  . 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 . Substituting Equation (10) into Equation (9) we obtain Equation (11), which is written in a form suitable to be solved.
where is a constant, the term is called Prandtl’s mixing length and its modeling varies depending on the flow. This method is used for solving mixing layer flows, and the most basic solution is to consider linear. Usually it is assumed as , where b is a constant that depends on the flow, where every value is obtained empirically, and shown by Silveira  .
2.5. Variable Mixing Length without Wall Function
Even though linear 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  , 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:
Applying the boundary condition given by Equation (15) we obtain:
Now, rearranging Equation (16), we have the following second order polynomial equation in :
This second-order equation in can be solved and two roots can be determined. Only one provides a physically consistent result, which is given by Equation (19):
In Equation (19), the presence of in the denominator generates indeterminacy, since it can assume null value. To remove this indetermination, the authors proposed the following algebraic manipulation was performed:
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  . For an illustration of the results that are obtained, the following values of the parameters were chosen for convenience and to facilitate solving:
The velocity distribution, as a function of the y-coordinate, taking into account the conditions given above, is shown in Figure 2. The velocity distribution shows a typical behavior of a turbulent flow, showing a strong velocity gradient near the walls of the channel. This behavior results from the strong diffusion of momentum, from the center of the channel towards the walls. As well-known, turbulence intensify the process of molecular diffusion. The Reynolds number that appears in Figure 2 was defined and calculated as bellow:
Equation (20) was nondiminesionalized, resulting in an equation with a single parameter. This procedure facilitates the comparison of the results obtained in the present study with results of other authors. Let’s introduce the following parameters without dimensions:
Two new terms appeared in Equation (23), the friction velocity ( ) and the friction velocity Reynolds number ( ). These parameters are quite used in the study of turbulent channels and ducts. The friction velocity is defined by Equation (24):
Figure 2. Mean velocity profile of a flow with ReD = 40950.
where 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 . Using the terms shown in Figure 3, and taken into account the forces, we obtain:
Taken in Equation (25) and using Equation (24), we obtain:
Using these parameters, Equation (20) can be transformed into Equation (27).
Multiplying both sides of Equation (27) by 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.
Equation (28) requires only one control parameter, the Reynolds number
Figure 3. Schematic illustration of a control volume inside a channel.
based on the friction velocity, , in order to characterize the physical flow to be simulated and compared with other authors. Once we obtain the velocity profile we can obtain the , using Equation (21) and Equation (22), as shown below:
3. Initial Results
3.1. 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  and ICES Turbulence Database  . Four different values were chosen to compare: 180, 395, 590 and 1000.
Beginning with and , 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 and are shown in Figure 5.
For each case it is possible to calculate using Equation (31), therefore: for , ; , ; , ; and , . By Figure 4 and Figure 5 it is possible to notice that the calculated curve better approaches DNS results for higher Reynolds number values.
Figure 4. Comparison of the profiles for (a) and (b).
Figure 5. Comparison of the profiles for (a) and (b).
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.
3.2. 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)  proposed a wall damping function, which is presented here by Equation (32).
However, Equation (32) is only valid for , thus it is necessary to separate the velocity profile in two parts. Cebeci and Bradshaw (1984)  , proposed an expression that combines Equation (12) and van Driest’s damping term given by Equation (32), resulting in Equation (33), which is valid for the whole domain and is already nondimensionalized:
In Figure 6 the mixing length is plotted with and without the damping term. It can be seen that with the damping model, the mixing length dampens more quickly when compared to the curve obtained without the damping function.
4. Final Results
4.1. Mean Velocity Profile
Applying (32) on ODE (28) and solving it, enables to generate the mean velocity profile for the four different stated before, now calculated with the damping function. Comparison with the DNS results from JHTDB  and ICES Turbulence Database  were performed. Beginning by and , or and respectively, the results are shown in Figure 7. After, or , and or . The results are presented in Figure 8.
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 , the 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.
4.2. Reynolds Stress
Combining Equations (6), (8) and (10) it is possible to calculate the Reynolds stress, with Equation (34), for .
Figure 6. Comparison of the mixing length with and without damping.
Figure 7. Comparison of the profile for (a) and (b) with and without damping function.
Figure 8. Comparison of the profile for (a) and (b) with and without damping function.
Table 1. Euclidean norm ( ) for the velocity profile curves comparison with DNS.
Table 2. Infinity norm ( ) for the velocity profile curves comparison with DNS.
Using the parameters in (23), it is possible to obtain the nondimensional Reynolds stress, using the wall stress , It is given by (35):
But, Equation (28) gives us the derivative exactly, and we can obtain the exact solution for :
Therefore, using Equation (36) it is possible to plot the Reynolds stress, the same four different values of Reynolds number. Starting with (a) and (b), plotted in Figure 9, we see a very good agreement with DNS data of  and  .
In Figure 10 the curves for (a) and (b) are plotted, we see that the agreement is still better.
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 t. 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.
Figure 9. Comparison of the Reynolds Stress profile for (a) and (b).
Figure 10. Comparison of the Reynolds Stress profile for (a) and (b).
Table 3. Euclidean norm () for the velocity profile curves comparison with DNS.
Table 4. Infinity norm () for the velocity profile curves comparison with DNS.
These results show how this theoretical method is reliable, since with only one simple numeric solution by Runge-Kutta for velocity and an exact solution for the shear stress, values and plots can be generated that compares very well to DNS (a method that can take days to run). But it is important to stress that the presented method is limited to this isolated physical problem.
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  model for Prandtl mixing length was used. Furthermore, the proposals of Van Driest  and Cebeci and Bradshaw (1984)  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 viscous shear stress rate was solved exactly. The results obtained for four values of the Reynolds number were compared with the results of Direct Numerical Simulation, and excellent agreement was obtained.
The proposed and solved flow does not have necessarily practical usage, since channel and duct flows are being studied experimentally for several years, with various well-defined tables and results, even for different roughnesses in the walls. However, the main objective of the current study, is for didactical reasons, to demonstrate a resolution of one of the most used models in teaching turbulence courses. The present paper may help students, interested in this scientifical area, to visualize a complete solution of a turbulent flow, for which the laminar regime solution they know since the first contact with fluid dynamics courses.
Another very important aspect is that the method presented in the current paper may be used for validation of numerical and computational methods. The presented model can be easily used with agility as validation reference. To date, the world scientific community has determined the rate of convergence of numerical methods using synthesized solutions for problems of purely mathematical nature. With the model proposed in the present work, one can determine the convergence rate of numerical and computational methods for real turbulent flows in a flat channel.
 Poiseuille, J.L.M. (1846) Recherches experimentales sur Ie mouvement des liquides dans les tubes de tres-petits diametres, In Memoires presentes par divers savants a l’Academie Royale des Sciences de l’Institut de France, IX, 433-544.
 Cheever, E. (2015) Linear Physical Systems Analysis Website, Swarthmore College,
 Johns Hopkins Turbulence Database (JHTDB) (2016) DNS Database of Johns Hopkins University.
 ICES Turbulence Database (2007) Institute for Computational Engineering and Sciences of University of Texas.