Development of a source code to analyze the effect of Reynolds number on a square lid driven cavity

Numerical simulations of two-dimensional, steady, incompressible lid driven flow in a square cavity were investigated in this work. A commercial finite volume package of ANSYSFLUENT was used to analyze and visualize the nature of the flow inside the cavity at different Reynolds Numbers. In addition, a MATLAB code was developed and validated by comparing the results with the reference values from literature. Staggered grid was employed in the discretization of the cavity to avoid checkerboard pressure while developing the code. The governing equations were discretized in terms of velocity and pressure fields. The artificial compressibility method was used to de-couple the pressure and velocity terms in the governing equations. A 129×129 grid system was used in both cases. The effects of Reynolds number (100 on the flow characteristics were illustrated through an analysis of stream function, velocity vector, pressure co-efficient and velocity contours. The thinning of the wall boundary layers with an increase in the Reynolds number is evident from the u and v velocity profiles along the vertical and horizontal lines at the geometric center, although the rate of this thinning is very slow for Re > 5000.


INTRODUCTION
The lid-driven cavity flow is the motion of a fluid inside a rectangular cavity created by a constant translational velocity of one side while the other sides remain at rest. Fluid flow behaviors inside lid driven cavities have been the subject of extensive computational and experimental studies over the past years. Applications of lid driven cavities are widely seen in material processing, dynamics of lakes, metal casting and galvanizing.
The investigation of a lid-driven flow in a rectangular cavity is motivated by three major factors. First, this flow is an idealized representation of several engineering situations, such as the flow over cutouts and repeated slots on the walls of heat exchangers or on the surface of aircraft bodies. Second, the lid-driven cavity flow offers the opportunity to study a "stationary, captive" primary vortex. Third, the regular geometry and easily posed boundary conditions have made this flow a popular test case for computational schemes, giving rise to a need for an accurate data base against which such schemes may be validated [2]. Research into the liddriven cavity flow structure is an area of continuing interest and was selected for a benchmark study in a major international workshop [3]. Hou et al. [4] have used the lattice Boltzmann method for simulation of the cavity flow. They have used 256 × 256 grid points and presented solutions up to Reynolds numbers of Re = 7500. In recent times, Barragy and Carey [5] have used a p-type finite element scheme on a 257 × 257 strongly graded and refined element mesh. The authors have obtained a high solution for steady cavity flow solution for up to Reynolds numbers of Re = 12,500. Although Barragy and Carey have presented qualitative solution for Re = 16,000, they have concluded that their solution for Re = 16, 000 was under-resolved and needed a far greater mesh size. According to Botella and Peyret [6], highly-accurate solutions for the lid-driven cavity flow are computed with the Chebyshev collocation method. Schreiber and Keller [7] have introduced an efficient numerical technique for steady viscous incompressible flows. The non-linear differential equations are solved by a sequence of Newton and chord iterations. Upon increasing the Reynolds number, the flows in the cavity of interest were found to comprise a transition from a strongly two-dimensional character to a truly three-dimensional flow and, subsequently, a bifurcation from a stationary flow pattern to a periodically oscillatory state. Finally, travelling wave instability further induced longitudinal vortices, which are essentially identical to Taylor-Gortler vortices [8]. Ercan Erturk presented that physically, the flow in a driven cavity is neither two-dimensional nor steady, most probably, even at Re=1,000. At high Reynolds numbers, when the incompressible driven cavity flow is considered as two dimensional and also steady, then the considered flow is a fictitious flow. It would be needless to study the hydrodynamic stability of a fictitious flow, i.e. the 2-D steady incompressible flow in a driven cavity at high Reynolds numbers [9]. Bruneau and Saad [10], have pointed out that driven cavity flows exhibit almost all the phenomena that can patterns, chaotic particle motions, instability, and turbulence. Critical comparison with the former numerical experiments confirms the high-accuracy of the method, and gives extensive results for the flow at Reynolds number Re = 1000. In this paper, flow analysis is done for a square lid driven cavity using a commercial software called ANSYS 16.1. The effect of Reynolds number is investigated and solutions are obtained for configurations with several Reynolds numbers such as 100, 400, 1000, 3200, 5000, 7500 and 10,000 and meshes consisting of as many as 129 × 129 points. Then, a self-written code in MATLAB is developed with proper discretization. Finally, a comparison of the two simulation results is presented.

MATLAB CODE DEVELOPMENT
A self-written MATLAB code was developed to simulate the flow in a square lid driven cavity. The whole procedure can be divided into mainly three parts which are the development of the code, simulation and lastly matching the results with the benchmark for validation. The problem statement is discussed with Figure 1 given below. The square domain is taken where the left wall, the right wall and the bottom wall are completely stationary. Therefore, on those walls, u velocity and v velocity that is the horizontal and the vertical velocity, respectively is 0 while the top wall is moving and because it is moving in the x-direction, only the u component of the velocity has some value which has been shown in Figure 1 that is u=U and the vertical velocity is zero. The first step to solve the equation is to discretize the domain or to create the mesh which means the whole domain is divided into small parts and an analysis will be done on every single small part individually. Figure 2 shows a very simple Cartesian and structured mesh that is created in where the whole square is divided into a number of smaller squares.
The macroscopic variables of interest in conventional numerical methods, such as velocity and pressure are usually obtained by solving the Navier-Stokes equation. Such numerical methods for two dimensional steady incompressible Navier-Stokes equations are often tested for code validation. The one sided lid-driven square cavity flow has been used as a benchmark problem for many numerical methods due to its simple geometry and complicated flow behaviors [11]. Initial attempts to solve the Navier-Stokes equations employed straightforward centered finite differences to the spatial operators on a regular grid, with the pressure and velocity components being unknown at the corners of each cell. Two types of instabilities are soon discovered, associated with this type of spatial discretization. The pressure can be highly oscillatory or even undetermined by the discrete system, although the corresponding velocities may be well approximated. The reason for this phenomenon is that the symmetric difference operator will annihilate checkerboard pressures, i.e., pressures which oscillate between 1 and -1 on each grid line connecting the grid points. In fact, if the vertices are colored in a checkerboard pattern, then the pressure at the black vertices will not be related to the pressure at the white vertices. Hence, the pressure is undeter-mined by the discrete system and wild oscillations or overflow will occur.
The remedy for oscillatory or checkerboard pressure solutions is, in a finite difference context, to introduce a staggered grid in space [12]. This means that the primary unknowns, the pressure and the velocity components, are sought at different points in the grid such as the displacement of velocity components rather being on the nodal points. Figure 3 displays such a grid in 2D. In this way, staggering of u and v velocity is done, which means that instead of defining the u and v velocity at the nodal point, they are defined at the half distance from the node. Finally, 7 by 8 mesh and 8 by 7 mesh will be created for the whole domain for u and v velocity respectively. Pressure points are located at the center of the squares. Here, an 8 by 8 mesh will be created for pressure points. The points of all these three variables are shown in Figure 3. It becomes very complicated when all the things are put together. Figure 3 contains three small volumes colored by blue, green and brown which describes the relationship between u, v and pressure point locations in the whole domain. The partial differential form of the governing equations are converted to system of algebraic equations form for each control volume. The analysis is done on the grids one by one. From every analysis, a new velocity and pressure will be obtained which will be used to find the consecutive values. Now, the conversion of X momentum equation is described through the green Colored area of the   The capital letters N, P, S, W & E denote the u velocities at the center and around the control volume. From the X momentum equation, it is seen that the first term is a time derivative: . (1) This term includes the unsteadiness in local flow. Equation (2) to Equation (6) are discretization of each term in Equation (1): , , .
The next step is to solve the continuity equation. Originally, the continuity equation is the divergence of the velocities is zero: But here, an extra term is added in the equation and for this additional term, this method is called an artificial compressibility method. This term is added on purpose. One of the early techniques proposed for solving incompressible Navier-Stokes equations in a primitive variable form is the artificial compressibility method of Chorin [13]. In this method, the continuity equation is modified to include an artificial compressibility term which vanishes when a steady state solution is obtained. The artificial term is the form of and for this it behaves like compressible flow. This particular scheme is more simple and easier than other methods like pressure correction method. In the equation, denotes the artificial compressibility constant. Usually, this value is taken as unity as described by Chorin [13]. The conversion of continuity equation is shown by the second equation: .
The updated velocities are already known from X and Y momentum equation. From continuity equation, it will be able to know the updated pressure. The artificial term, (artificial compressibility constant) actually makes this possible. It decouples the pressure and velocity terms in incompressible equation very simply for which there is no need to use any pressure correction scheme. Now, after completing discretization and conversion, boundary condition should be employed properly,; otherwise it is not possible to get the correct result. As it is already known that the left, bottom and top boundaries are stationary that means u and v velocities should be made zero. This is the reason behind the extra layer of u velocity. To maintain the velocity zero at the bottom wall, this extra layer of velocity (-u) m/s is used so that the average will be made zero at the points on bottom wall. But the top wall is moving with a velocity for example, the value is taken as 1 m/s. So, to keep the velocity 1 m/s at the points on the top wall, the upper and lower velocities of the top wall should be maintained as velocity 1 m/s so that the average will be 1 m/s. For the v velocity, similar procedure is applied. Across the boundaries, there should be no pressure difference as no flow occurs across the boundaries. To maintain no flow across the boundaries, pressure at interior and exterior points across each boundaries should be equal so that there will be no pressure gradient and no flow will occur.

Defining the variables
Set the value of x grid y grid, height of the geometry, length of the cavity geometry, velocity of the lid, residual:  grid_x = 129,  grid_y = 129,  length = 1; %Lenght of the geometry,  height = .5; %Height of the geometry,  steps = 1; %Number of steps that the calculation run. Initially we assume it as one,  velocity = .00010048; %Velocity of the top wall.
User can change it. It will dynamically change;  residual = .000001; %this is residual. How accurate we are trying to get the results.

Set the value of known property
Density of the water(ρ), dynamic viscosity (μ),Time setps (dt):  row = 998.2; %It is Greek letter row. which represents the density of water;  miu = .001003; % It's a Greek letter miu. which represents the dynamic viscosity of water;  dt = 0.001; %Time steps.

Calculating the required properties
Reynolds number from the following equation: .
Reynolds number is a dimensionless value where L is the cavity length, U is the flow velocity attached to the upper lid, is the density of fluid and is the dynamic viscosity. Unit length of the grid along x direction (dx), Unit height of the gird along y direction (dy):  Re = (row * velocity * alength) / miu; %We are calculating the Reynolds number here;  dx = alength/(grid_x-1); %Lenght of each small grid;  dy = height/(grid_y-1); %Height of each small grid.

Setting the initial u & v velocity
For u velocity the value of all the points will be zero except for the velocity of top two point, which will be the same as the velocity that set in step one. For v velocity, the value of all the points will be zero. As the lid doesn't have any velocity component on vertical direction. And for setting this velocity, no extra code is needed as the array was defined with MATLAB built in function zeros.

Solving the u velocity
Here the new u velocity is calculated by using the algebraic form of x-momentum equation. Which is described above. This new u velocity is stored in another array un.
for j=2:grid_x un(1,j) = 0;%This will set velocity 0 at the left wall un(grid_x,j) = 0;%This will set velocity 0 at the right wall. end Bottom wall is also stationary. So the average of u(i,1) and un(i,2) will be zero: , (10) . (11) The top wall is moving. Its velocity was set on step one. Its velocity needs to be same on all the iteration of the while loop.

Solving the continuity equation
The algebraic form of the continuity equation is solved here. For the proper result, the updated u and v velocity which has been stored in un and vn is used. By solving the continuity equation, the pressure can be determined. Relevant code: for i=2:grid_x for j=2:grid_y pn(i,j) = p(i,j)-dt*delta*( ( un(i,j)-un(i-1,j) )/dx + ( vn(i,j)vn(i,j-1) ) /dy ); end end As no fluid passes through the wall so the pressure will be same between the two grid of the wall. So pn(i,1)=pn(i,2). for i=2:grid_x pn(i,1) = pn(i,2);%This will set the bottom imaginary pressure as the bottom pressure. pn(i,grid_y+1) = pn(i,grid_y);%This will set the top imaginary pressure as the top pressure. end for j=1:grid_y+1 pn(1,j) = pn(2,j);%This will set the left imaginary pressure as the left pressure. pn(grid_x+1,j) = pn(grid_x,j);%This will set the right imaginary pressure as the right pressure. end

Calculating the error
Initially, the error is set to 0 to replace the previous error value in while loop. The error of any point is the divergence of the velocity. The error is calculated for all the points. And all the error value of all the point is summed up. error = 0; for i=2:grid_x for j=2:grid_y m(i,j) = ( ( un(i,j)-un(i-1,j) )/dx + ( vn(i,j)-vn(i,j-1) )/dy ); error = error + abs( m(i,j) ); end end

Displaying error & redefining u,v and p
The error is shown after every 1000 steps.
if rem(steps,1000) == 0 fprintf( 'Error is %f for the step %f\n',error,steps); end The value of u,v and p is set from the value un, vn, pn. So that for next iteration of while loop the updated value of u,v and p is used.

Calculating the grid point value
So far, the staggered velocity and pressure have been calculated. Staggered velocity and pressure is averaged to get the grid point value.

RESULTS VALIDATION AND DISCUSSION
The present work can be divided into three parts. A MATLAB code is developed at first to simulate the driven cavity flow, secondly, code validation is done comparing the results with experimental values from Ghia et al. [1] for Reynolds number 100, 1000 and 5000 and also with the results from commercial software ANSYS FLUENT 16.1. Thirdly, the effect of Reynolds number on the driven cavity flow is investigated. The ANSYS simulation was done first.
Three types of mesh were used that are 129x129, 161x161 and 257x257. The results from each mesh size match very closely with one another. Thatis why 129x129 mesh size is used for every simulation in this study as it becomes easier and less time consuming to simulate the flow with MATLAB code. The reference paper Ghia et al. [1] also represented their results with the mesh size 129x129. The u and v profile results done by MATLAB along the vertical line and horizontal line respectively through the geometric center for Re=100, 1000 and 5000 are compared with Ghia and ANSYS simulation in Figure 6 and 7.
As velocity changes with the change of Re where, , (L is the cavity length) velocity ranges from .0001 to .01 as Re increases from 100 to 10,000 in this investigation. Each graph shows the same manner of velocity change at various Reynolds numbers. The thinning of the wall boundary layers with an increase in Reynolds number is evident from the profiles in Figure 6 and 7, although the rate of this thinning is very slow for Re > 5000. The near linearity of these velocity profiles in the central core of the cavity is indicative of the uniform vorticity region that develops here for large Reynolds number. The high Reynolds number profiles of u exhibit a kink near y = 1, while a similar behavior is observed for the v profiles near x = 1. Such behavior has been reported by some previous investigators, and it is seen to persist in the present ANSYS and MATLAB solutions. This would imply that the velocity distributions near these walls are not extremely sensitive to mesh size. The dimensionless pressure drop coefficient for the cavity is an important parameter in the design of the cavity. The pressure drop coefficient is a strong function of the position and the Reynolds number. From Figure 8 and 9, It can be noted that the pressure drop coefficient is higher for a lower Reynolds number. The pressure drop coefficient attains its maximum value within the above half cavity, whereas for most Re, it gains the minimum value with the bottom half cavity.
The streamline contours for the cavity flow configurations with Re increasing from 100 to 7500 are shown in Figure 10. All the Figures show that a primary vortex is generated. Except for the top wall, other walls are stationary. Thus, fluid cannot move across these boundaries. Boundary layers are formed on the solid walls. Near the solid boundaries, viscous force is dominant with respect to the inertia force of the fluid. But, far from the solid boundaries, inertia forces become dominant and as boundary layers form on the left, right and bottom walls, due to inertia forces of the fluid faces an obstacle towards the flow near the boundary layers and for this reason, a change of momentum of fluid occurs that in turn creates rotation around the center. Eventually, the rotation of center turns into a vortex around which fluid rotates. In case of lower Re, the primary vortex is created near the top wall. As, Re increases and boundary layers near the solid walls are becoming thinner, the location of primary vortex changes downward towards the geometric center of the cavity. In the range of Re = 100 to 1000, the location of primary vortex changes slowly. At Re = 3200, the location of primary vortex is near the geometric center and another small vortex is seen at the right-bottom corner of the cavity. This vortex is named as secondary vortex 1. This secondary vortex becomes larger with the increase of Re. Another secondary vortex is created at Re = 7500 at the left-bottom corner of the cavity. Also, a tendency to generate the third secondary vortex is observed at the top-left corner with the increase of the Re which is due to the fact that this vortex only becomes really apparent after Re=3200.

CONCLUSIONS
The investigation led to several conclusions: 1. The thinning of the wall boundary layers with increase in Reynolds number is evident from u and v velocity profiles along the vertical and horizontal line at the geometric center, although the rate of this thinning is very slow for Re > 5000.