Numerical simulation of the laminar boundary layer flow over a flat plate

Understanding the dynamic and thermal characteristics of the boundary layer (BL) flows is an essential step toward the best design and operational conditions for many engineering devices. Flow characteristics within the boundary layer are governed by two forces that are in a mutual race to dominate the flow, the viscous and inertial forces. The state of the flow is determined by the relative domination of these two forces inside the boundary layer zone. Theoretical solutions for many BL flow types existed since the beginning of the 20th century. Theoretical solutions are strictly possible for simple boundary layer flow configurations within the laminar range and before flow separation occurs. For practical situations where boundary layer flow problems are prescribed with complex flow configurations, CFD techniques represent the most suitable approach to tackle such types of flows. In this study, the laminar boundary layer flow over a flat plate has been solved numerically leading to a detailed analysis of the boundary layer flow with accuracy comparable to that of the theoretical solution.


INTRODUCTION
Since its discovery by Ludwig Prandtl in 1904, the boundary layer constitutes one of the main physical phenomenons that retained the greatest attention of the scientific community for its importance, and hence became in its own a major research field in fluid mechanics.The boundary layer phenomenon emerges as a result of viscous effects that tend to suppress fluid motion near the surface of a solid wall.This kind of interaction between the fluid layers adjacent to the wall and the solid surface modifies the dynamic and thermal characteristics of the flow field.For large Reynolds number viscous flows, as viscous diffusion propagates faster across the flow than the advection speed of the fluid particles, viscous diffusion time becomes much shorter than the residence time and hence a thin boundary layer arises.In this region where viscous effects dominate fluid motion the velocity field takes an altered shape constrained by the presence of the wall where the velocity vanishes and by the competing interaction between viscous and inertial effects.Consequently, the flow field can be classified into two distinct regions, the region away from the wall in which viscous effects are negligible and a region immediately adjacent to the walls in which viscous and inertial effects are significant .The mathematical formulation of the boundary layer equations was first introduced by Prandtl in his paper of 1904 which have been studied and solved by Blasius in 1908 for the case of viscous laminar flow over a flat plate [1].For a while, Von Karman (1921) proposed his approach to solve the boundary layer problem using the momentum integral equation [2].Burgers (1925) revealed by experimental observations the velocity distribution across the boundary layer on a flat plate, bringing to light the simultaneous presence of laminar and turbulent regions [3] .Tollmien (1929) considered the Blasius velocity profile near a flat plate and obtained the critical Reynolds number above which the flow becomes unstable to a traveling-wave type of disturbances in a certain frequency range [4].One can also mention some of the early boundary layer works such as the wedge flow boundary layer studied by Flakner -Skan (1931), the linearly retarded flow studied by Hartee (1937) and Stewartson (1954), the wake flow past a circular cylinder studied by Goldstein (1933) and the boundary layer approximations for a jet flow studied by Schlichiting [4].During the era after WW2 (World War 2), special attention was paid to boundary layer problems linked to the aviation and aeroplan industry [5].The stability of compressible boundary layers was first considered by Lees & Lin (1946) and followed by Lees (1947), Dunn & Lin(1955), Lees & Reshotko (1962), and Mack (1965) [6].The subject has increasingly attracts many research works that try to reveal the detailed physics of the boundary layer flows and their instability inducing mechanisms which are full of complicity and rich of physical phenomenon associated to this type of flow.The induced drag force over the surface of an aeroplane wings, ships and missiles, the power generation through the compressor and turbine stages in jet engines, wind turbine capacity to produce clean energy, the effectiveness of air intakes for ram and turbojets and so on are all linked to the concept of the boundary layer and its associated flow physics.Some industrial applications need to take into account the effects of boundary layer flow induced by a stretching surface as such the case of polymer extrusion processes, melt spinning processes, aerodynamic extrusion of plastic sheets, glass fiber production, the cooling and drying of paper and textiles [7].These practical flow problems are generally associated with complex flow configurations and need to be solved in the real three dimensional space.Therefore, classical theoretical techniques can't be used to solve the governing equations of the flow.Here comes the necessity for using the Computational Fluid Dynamics (CFD) techniques which can handle any kind of flow and describe the flow behavior with an accuracy which is very close to that of theoretical solutions.There are three classical numerical techniques used to solve the governing equations of fluid flow problems which are the Finite difference, the finite volume and the finite element methods.For each of these methods, the partial differential equations that govern the flow are transformed explicitly or implicitly into a system of linear equations that can be efficiently solved using algorithms designed for this purpose.The current study presents numerical resolution of a 2D steady boundary layer flow problem over a flat plat using finite difference method.Likwise, in order to validate simulation results the problem has been solved using the classical theoretical approach (Blasius solution).Profiles of all flow variables resulted from the two solution techniques have been compared and analyzed.

GOVERNING EQUATIONS
The governing equations of the boundary layer flow can be derived from the Navier-Stokes (N.S) equations using some realistic physical assumptions that reduce the governing equations into the following parabolic form [8]: With the continuity equation written as Since there is no pressure gradient following x direction the governing equations simplify further to Subjected to the boundary conditions (B.C):

BLASIUS SIMILARITY SOLUTION
The first attempt to solve boundary layer equations were done by Blasius relying on the concept of local similarity of the non-dimensional velocity profile with respect to an appropriate dimensionless similarity variable ( ) Where Using the following definition of a stream function The boundary layer equation can be transformed after some mathematical manipulation to a third order ordinary differential equation Subjected to the following transformed BC's: This third order nonlinear ordinary differential equation (ODE) can be solved numerically by decomposing it into a system of three first-order ordinary differential equations [9].The resulting coupled system of initial value problems can then be integrated easily using second or fourth order Runge-Kutta method.The resulting values of the non-dimensional stream function

( )
f  and its derivatives can then be used to compute all flow variables within the BL (velocity profile, the BL thickness, shear stresses ... etc).

NUMERICAL SOLUTION
We consider a flat plate of length L of infinite width and negligible thickness, that lies in the  −  plane, and whose two edges correspond to  = 0 and  =  as shown in the Figure (1.).The plate is supposed to be subjected to a uniform flow with a velocity field given by   .The flow over the plate is surrounded by appropriate boundary conditions arising from the noslip condition where the velocity components at the surface of the wall vanish ( =  = 0) and the formation of a thin boundary layer, of thickness () , outside which the boundary layers remains effectively inviscid.It follows that the flow external to the boundary layer is unaffected by the presence of the plate.Hence, the streamwise velocity at the outer edge of the boundary layer is U(x)=  .The governing equations (eq:3, eq:2) have been discretized using uniform grid in the streamwise direction and a non-uniform grid that varies as a hyperbolic tangent profile in the normal direction to the wall ( y e ).
The use of non-uniform hyperbolic tangent grid profile in  direction clusters grid points in the flow regions near the wall where the flow field undergoes big variations and changes [10].Grid resolution has been kept the same throughout all simulations realized in this work.
The resolution procedure starts with the discretization of the governing equations using appropriate finite difference formulas to approximate the different derivative terms in these equations.Since flow variables following the normal direction to the boundary exhibit higher variations than that in the streamwise direction, derivatives in the y e direction have been discretized using the central finite difference scheme which is second order accurate in y.Other terms have been discretized using the forward finite difference scheme which is first order accurate.The solution procedure admitted the use of marching technique following downstream direction   , which can be either explicit or implicit.Explicit marching scheme is simple but suffers from the problem of numerical instability unless a very small marching step size ( x  ) is used.Contrary to the explicit marching scheme, implicit marching scheme is unconditionally stable and no restrictions are made to the choice of the marching step size.Therefore, the implicit marching scheme has been adopted to solve the boundary layer equations and model the flow.The equations have been discretized on the rectangular grid as shown by figure (2).The mesh point (,) indicates the actual location where velocity field is completely known or already calculated which would be used to march downstream and predict the next station values at ( + 1,).The discretized governing equation on the non-uniform grid using an implicit marching scheme along   is as follows: i j i j i j u u u u y y i j j i j i j i j j uv i j i j yy x y y j j j j This can be converted into a discrete tridiaogonal form: This corresponds to the following tridiaogonal matrix system: The implicitly discretized momentum equation at grid points along the normal direction performs a tridiagonal matrix system for each downstream location.The constructed tridiagonal matrix system can be solved by any iterative means or using matrix inversion to get the velocity component (  ,   ) at every mesh point.There is no instability constraints regarding the value of the marching step ( x  ).The downstream velocity component can then be used to predict the values of the normal velocity component (  ,   ) at each grid site using the discretized continuity equation: To improve the accuracy of the finite difference modeling of the continuity equation, the first terms u  / x  which represents the first derivative of  with respect to x have to be evaluated at the grid level ( − 1 2 ) as the second term of the continuity equation v  / y  .This approximation is done by calculating the average of u  / x  using its values at the grid level  and the grid level  − 1 as follows: The next downstream value of the normal velocity component can then be obtained: The solution has been conducted using an algorithm written in FORTRAN language and compiled on a Linux system.The advantage of using a Linux system lies on the availability of different FORTRAN compilers and different post-treatment tools free of charge.The resolution algorithm starts by defining the different variables and parameters needed, and then the grid has to be generated followed by setting the boundary conditions.Since the implicit scheme used is based on using the downstream direction as a marching direction, the discretized boundary layer equations at all the nodes in the normal direction constitute a tri-diagonal system and that is for each downstream station.These matrix systems are solved consequently starting from the leading edge of the plate to the trailing edge using Thomas algorithm known as TDMA (Tridiagonal matrix algorithm).The predicted values of downstream velocity component at each downstream location is used to update the normal velocity component using the discretize continuity equation.

RESULTS AND DISCUSSION
The governing equations of the steady laminar boundary layer flow are solved using flow and simulation parameters described by the table (1.)The dynamic viscosity could be determined according to the suggested values of Reynolds number ( Re /  ) which is defined based on the free stream velocity U and the plate length

Boundary layer velocity profile
The main flow feature that acquires much of interest is the characteristics of the velocity field.Many physical flow phenomenon (such as shear stress, convection flux, viscous diffusion, vorticity .. etc) can be described mathematically using flow velocity field.The downstream velocity profile can be regarded as an indicative tool to evaluate the accuracy of flow simulation results when compared to experimental measurements or analytical solutions.Downstream velocity component at 4 Re 5 10 = is compared in Figure (5) to the corresponding values obtained experimentally by (Liepmann (1943)) [11] and also to downstream velocity profile issue of Blasius analytical solution.Velocity profile obtained by numerical simulation superposes exactly to the corresponding profile given by Blasius analytical solution.There are some little discrepancies between both simulation and analytical profiles from experimental values; this may be owed to measurement errors since these measurements are conducted using old measuring techniques such as hot wire probes or Pitot tube probe.Figure ( 6) compares downstream velocity profile of both numerical and Blasius analytical solution at different Reynolds number.The overall behavior of the two solution approaches matches well except at higher values of Reynolds number where there are some discrepancies between analytical and simulation results very close to the wall that may be caused by inappropriate mesh resolution that need to be increased to allow grid resolution captures all features of the flow in this region.The normal velocity component profiles at 4 Re 5 10 = for both numerical and analytical solution are presented by Figure (7).This velocity component as expected acquires low values since the main flow stream is in the downstream direction.Also at this Reynolds number the two velocity profiles take exactly the same values at each non-dimensional location (, ).

8 )
The flow chart presented by Figure (3) resumes the resolution procedure used to solve the considered boundary layer flow problem.Issue 4 (December 2023)

xL
. The parameters , xy NN denote the number of grid points following   and y e directions respectively.Boundary layer flow simulations have been conducted using 10 values   (kg/ 3 m ) ( / ) ms U   (m)   (m)

Figure 6 :
Figure 6 : Streamwise velocity profile for different values of Reynolds numbers.

Figure 7 :
Figure 7: Profile of the normal velocity component