Large Eddy Simulation of Laminar and Turbulent Flow on Collocated and Staggered Grids

Abstract

Momentum and Continuity are the basic equations for fluid flow modeling. The momentum equations in their final form are known as Navier-Stokes equations and can be solved using different numerical methods. There are several approaches such as SIMPLE, PISO, and Fractional Step for solving these equations. In solution procedure, it needs to decide where to store the velocity components. Staggered and Collocated grids can be used to evaluate this problem. On Staggered grids, the velocity components are stored at the cell face, and the scalar variables such as pressure are stored at the central nodes. However, on Collocated grids, all parameters are defined at the same location at the central nodes. The Staggered grids method gives more accurate pressure gradient estimation. However, Collocated grids method is simpler for solving the equations. In this paper, for solving Navier-Stokes equations, Collocated and Staggered grids are employed. Comparison of horizontal and vertical velocities and stream lines at various Reynolds numbers was performed. The results were validated using standard tests such as lid-driven cavity, channel and backward facing step. Discussion is made on accuracy of these methods to estimate horizontal and vertical velocity profiles.

1. Introduction

Thefundamental equations of motion known as Navier-Stokes (N-S) equations are basically valid for laminar flow. To use them for turbulent flow, time, and spatial averaged form of these equations are used [1]. Different approaches such as Large Eddy Simulation (LES) and Direct Numerical Solution (DNS) models can be used to model turbulent flow. The classical models use the Reynolds Average Navier Stokes (RANS) equations. Complication of computational procedure for unsteady turbulent flow simulations is due to the very large number of grid points to be resolved for a wide range of length scales. In DNS, all scales are resolved both spatially and temporally. In LES approach, the equations are filtered to separate large eddies from small ones. In these models, it is assumed that the largest eddies interact strongly with the mean flow and contain most of the energy. Therefore, the main effects of turbulence can be modeled by this approach accurately [2].

There are several methods such as SIMPLE, PISO, and N-S equations. The SIMPLE algorithm was originally developed by Patankar and Spalding (1972) on the Staggered grids arrangement and involves one predictor step and one corrector step. While PISO method can be regarded as an extension of SIMPLE method with further correction steps [1]. The Fractional step method of Kim and Moin (1985) provides an approach that does not need to predict pressure in the predictor step. In this method, the N-S equations split into three equations containing convective, diffusive and pressure terms [3].

In this paper, a turbulent model was developed using Staggered and Collocated grids. To validate the model some standard tests such as lid-driven cavity, Channel flow and Backward-facing step were used.

Several researchers studied these tests. Prasad and Koseff [4] performed a series of experiments in lid-driven cavity for Reynolds number between 3200 and 10000 and spanwise aspect ratio (SAR) between 0.25 : 1 and 1 : 1. They showed that the flow is strongly affected by the SAR and Reynolds number so that this effect is different at low and high Reynolds numbers. Gokarn et al. [5] solved incompressible N-S equation using Large Eddy Simulation and partial-staggered variable arrangement. They validated the numerical scheme by simulating lid-driven cavity flow and flow over a Backward facing step. The results were compared with experimental data and demonstrated good agreement. Kim et al. [6] solved the unsteady N-S equation and DNS method at a Reynolds number of 3300 to simulate the fully developed channel flow. They compared the numerical results with the existing experimental data and discussion was made on origin of errors. Nie and Armaly [7] estimated the boundaries of the reverse flow regions in a three-dimensional (3D) backward-facing step using Laser-Doppler velocity measurement as a function of the Reynolds numbe in a range of 100 to 8000 to cover the laminar, transitional and turbulent flow regimes. Their results well agreed with measurements of Alamry et al. [8] in the laminar and turbulent flow regimes.

2. Governing Equations

The governing equations of incompressible fluids which are Continuity and N-S equations can be expressed as𝜕𝑢𝜕(𝑢)𝜕𝑥+𝜕𝜕𝑣𝑦=0,𝜕𝑡+∇⋅(𝑢𝑉)=−𝜌1𝜕𝜕𝑃𝑥+∇⋅(𝜈∇𝑢)+𝑆^𝑥𝜕(𝑣),𝜕𝑡+∇⋅(𝑣𝑉)=−𝜌1𝜕𝜕𝑃𝑦+∇⋅(𝜈∇𝑣)+𝑆^𝑌,(1) where 𝑡 is time, 𝑉 is the velocity vector, 𝑢 and 𝑣 are velocity components in 𝑋 and 𝑌 direction, respectively, 𝑃 is the static pressure, 𝜈 is the kinematic viscosity, 𝜌 is the fluid density, 𝑆_𝜙 is the source term including gravity acceleration in the 𝜙 direction [1].

In this research, LES model is used to simulate turbulent flow. It relies on the definition of large and small scales in which the former are solved and the later are modeled. Figure 1 shows the schematic view of this procedure.

Figure 1 

Schematic view of the simplest scale separation operator in LES model.

To separate the large scales from the small ones, filtering process is used. A filtered variable, denoted by an over bar, is defined as𝑓(𝑥)=𝑓^(𝑥)𝐺𝑥,𝑥^,Δ𝑑𝑥^,(2) where 𝐷 is the entire domain, 𝐺 is the filter function, and Δ is the filter width, which is the wavelength of the smallest scale. The filter function determines the size and structure of the small scales. If the filtering operation is applied to the governing equation, the filtered equation of motion is obtained. For incompressible flow of a Newtonian fluid the conservation form can be expressed as𝜕𝑢_𝑖𝜕𝑥_𝑖𝜕=0,𝑢_𝑖𝜕𝑡+𝜕𝜕𝑥_𝑗𝑢_𝑖𝑢_𝑗=−1𝜌𝜕𝑝𝜕𝑥_𝑖−𝜕𝜏^𝑖𝑗𝜕𝑥_𝑗+𝜈𝜕²𝑢_𝑖𝜕𝑥_𝑗𝜕𝑥_𝑗.(3) These equations just govern the evolution of the large scales. Therefore the effect of the small scales must be modeled through a subgrid-scale (SGS) stress term as𝜏_𝑖𝑗=𝑢_𝑖𝑢_𝑗−𝑢_𝑖𝑢_𝑗.(4) The Smagorinsky model is one of the various SGS models. This model can be expressed as [9]:𝑣_𝑡(𝑖,𝑗)=𝐶_𝑆Δ²2|||𝑆_𝑖𝑗|||²^1/2,𝑆_𝑖𝑗=12𝜕𝑢_𝑗𝜕𝑥_𝑖+𝜕𝑢_𝑖Δ=Δ_𝑥𝜕𝑥_𝑗,Δ_𝑦^1/2,(5) where 𝑣_𝑡 is turbulence eddy viscosity, 𝐶_𝑠 is a constant assumed to lie in the range of 0.065 to 0.25 [10], Δ is the length of filter, 𝑆_𝑖,𝑗 is the strain rate tensor, 𝑢_𝑖, 𝑢_𝑗 are the filtered velocity in 𝑋 and 𝑌 direction and Δ_𝑥, Δ_𝑦, and Δ_𝑧 are the mesh size in 𝑋, 𝑌 and 𝑍 direction, respectively.

Finally, the total kinematic viscosity is derived as 𝑣(𝑖,𝑗)=𝜇𝜌+𝑣^𝑡(𝑖,𝑗).(6)

3. Staggered and Collocated Grids

The solution algorithm for the governing equations is consisting of discretization of flow domain and storing the velocity vectors. Two methods denoted as Collocated Grids (C.G) and Staggered Grids (S.G) are used in this procedure. On C.G, all velocity components and scalar variables such as pressure are stored in the same locations (see Figure 2).

Figure 2 

Control volume in C.G.

Although application of C.G is simple, it cannot distinguish between a non-uniform pressures field and a uniform one. In S.G method, velocities are stored at the faces of control volume (see Figure 3).

Figure 3 

Control volume in a S.G.

One of the advantages of the S.G arrangement is that velocities are generated at locations where needed for the convection computations. While on C.G, stored velocities at the center of control volumes cannot be used directly. In this situation, for calculation of convection term, Rie-Chow interpolation is used.

4. Solution Algorithm

The solution algorithm employed in this study is a finite volume method based on the PISO predictor corrector algorithm coupled with a pressure correction to discretize the governing equations on structure grids.

In this method, the filtered N-S equations, are discretised using the finite volume method where the domain 𝐷 is divided into cells 𝛿𝑉_𝑖 so that ∑𝛿𝑉_𝑖=𝐷 and ∩_𝑖(𝛿𝑉_𝑖)=𝜙. Integration of the dependent variables over each cell together with application of Gauss’ theorem generates a set of discretised equations. The two-dimensional N-S equations are discritized using Finite volume method as 𝑎_𝑃𝑈𝑢_𝑃=𝑎_𝑁𝐵𝑢_𝑁𝐵−𝑃(𝑖+1,𝑗)−𝑃(𝑖,𝑗)𝑑𝑥+𝑆^𝑢𝑎_𝑃𝑉,𝑣_𝑃=𝑎_𝑁𝐵𝑣_𝑁𝐵−𝑃(𝑖,𝑗+1)−𝑃(𝑖,𝑗)𝑑𝑦+𝑆^𝑣,(7) where 𝑢_𝑁𝐵 and 𝑣_𝑁𝐵 are velocity components in X and Y direction and in neighboring cells (𝑁𝐵=𝐸,𝑊,𝑁,𝑆) and 𝑎_𝑁𝐵 are the related coefficients needs to be estimated. There are several schemes such as upwind, central differencing and power law to estimate these coefficients.

In this research, to estimate the coefficients power law scheme is used as𝑎_𝑖=𝐷_𝑖∗Max0,1−0.1||𝑃𝑒_𝑖||⁵+Max𝐹^𝑖,0,(8) where 𝑃_𝑒𝑖, 𝐷_𝑖, and 𝐹_𝑖 are the Peclet number, diffusion and convection terms on neighboring cells, respectively. The peclet number is the ratio of convection to diffusion terms. Theses terms are estimated on S.G as𝐷_𝐸=𝐷_𝑊=𝜇(𝑖,𝑗)𝑑𝑦𝐷𝑑𝑥_𝑁,=𝐷_𝑆=𝜇(𝑖,𝑗)𝑑𝑥𝐹𝑑𝑦_𝐸,𝐹=𝐹_𝑊=𝜌𝑈𝑑𝑦,_𝑁=𝐹_𝑆=𝜌𝑉𝑑𝑥.(9) Diffusion terms are estimated on C.G similar to S.G. However, convection terms are estimated using interpolation techniques.

The discritized N-S equations must be solved for all control volumes to estimate the velocity components. But in this step, the unknown pressures in the control volumes should be predicted.

In next step, Momentum equation is solved using predicted pressure field and Three Diagonal Matrix Algorithm (TDMA) as𝑎_𝑃𝑈𝑢_𝑃=𝑎_𝑁𝐵𝑢_𝑁𝐵−𝑃^∗(𝑖+1,𝑗)−𝑃^∗(𝑖,𝑗)𝑑𝑥+𝑆^𝑢𝑎_𝑃𝑉,𝑣_𝑃=𝑎_𝑁𝐵𝑣_𝑁𝐵−𝑃^∗(𝑖,𝑗+1)−𝑃∗(𝑖,𝑗)𝑑𝑦+𝑆^𝑣𝑃=𝑃_∗,+𝑃^,(10) where 𝑃^∗, 𝑃, and 𝑃^ are the predicted, corrected and correction pressure, respectively.

In next step, the pressure correction equation is derived as𝑎_𝑃𝑃𝑃^𝑃=𝑎_𝐸𝑃𝑃𝑎^𝐸_𝐸𝑃+𝑎_𝑊𝑃𝑃^𝑊+𝑎_𝑆𝑃𝑃^𝑆+𝑎_𝑁𝑃𝑃^𝑁,𝑎=𝑎_𝑊𝑃=𝜌𝑑𝑦,_𝑆𝑃𝑎=𝑎_𝑁𝑃=𝜌𝑑𝑥,_𝑃𝑃=𝑎_𝐸𝑃+𝑎_𝑊𝑃+𝑎_𝑆𝑃+𝑎_𝑁𝑃.(11) Then, the pressure and velocity components are corrected as𝑢(𝑖,𝑗)=𝑢(𝑖,𝑗)+𝑃𝑑=𝑦𝑃^∗+𝑃^,𝑎_𝑃𝑈(𝑖,𝑗)𝑃^𝑣(𝑖,𝑗)=𝑣(𝑖,𝑗)+𝑑𝑥(𝑖,𝑗)−𝑃^(𝑖+1,𝑗),𝑎_𝑃𝑉(𝑖,𝑗)𝑃^(𝑖,𝑗)−𝑃^(𝑖,𝑗+1).(12) The procedure is continued by solving the second pressure correction equation, rectifying the pressure and velocity with pressure correction and repeating the stages until convergence.

5. Model Validation

In this research, the N-S equations are discretized using a finite volume method of second-order accuracy in space based on Eulerian description at each time step. The N-S solver is developed based on S.G and C.G with a LES model as turbulent.

To validate the model, some standard tests such as lid-driven cavity, Channel flow and Backward-facing step are used.

5.1. Lid-Driven Cavity

The fluid flow is in a rectangular container which moves tangentially to itself and parallel to one of the side walls. Due to the simplicity of the cavity geometry, applying a numerical method on this flow problem in terms of coding is quite easy and straight forward. Despite its simple geometry, the driven cavity flow retains a rich fluid flow physics manifested by multiple counter rotating recirculation regions on the corners of the cavity depending on the Reynolds number.

The geometry of the problem consist of boundary condition of velocity and the successive of eddies which is shown in Figure 4.

(a)

(b)

Figure 4 

Lid-driven cavity test, (a) boundary condition (b) arrangement of eddies.

A sensitivity analysis is performed on mesh size (see Figures 5 and 6). As seen in Figures 5 and 6 by increasing the mesh size, the results become independence of it. For example, for Re=100 and 3200 and for 𝑁𝑋=𝑁𝑌=80 and 150, the results on C.G and S.G become independence of mesh size.

(a)

(b)

Figure 5 

Mesh independency on S.G for (a) and (b) .

(a)

(b)

Figure 6 

Mesh independency on C.G for (a) and (b) .

The model is then performed on S.G and C.G for various Reynolds numbers at a range of 100 to 10000 and results were compared with those of Ghia et al. [11] that is regarded as the true solution of this problem [12]. Figures 7, 8, 9, and 10 show the results for horizontal and vertical velocity components.

(a)

(b)

Figure 7 

Comparison of horizontal velocity at various Reynolds numbers (a) C.G and (b) S.G.

(a)

(b)

Figure 8 

Comparison of vertical velocity at various Reynolds numbers (a) C.G and (b) S.G.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9 

Comparison of horizontal velocity on C.G and S.G and Ghia et al. results.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10 

Comparison of vertical velocity on C.G and S.G and Ghia et al. results.

The streamlines of fluid flow in cavity on C.G and S.G are shown on Figure 11.

(a)

(b)

Figure 11 

Comparison of stream lines for various Reynolds numbers (a) S.G and (b) C.G.

For Re=100 the horizontal velocity at the bottom and top of the domain are zero and 1 m/s, respectively. The maximum negative velocity is −0.2 m/s occurring at the middle point of the domain. The maximum negative velocity is increased due to an increase in Reynolds number. The location of this maximum velocity is moved toward downstream. For example for Re=10000, this velocity is increased from −0.2 m/s to −0.42 m/s and it is location is moved from middle point to 0.07 m above of bottom. So by increasing the Reynolds number, the horizontal velocity gradient near the bottom and top is increased (see Figure 9). The reason for increasing the velocity gradient at two sides of the domain is that the thickness of boundary layer decreased as the Reynolds number increased leading to an increase in flow turbulence.

For Re=100 the vertical velocity at both sides of the domain is zero and maximum negative and positive velocities occurs at 0.25 m and 0.8 m above the bottom equal to 0.17 m/s and −0.24 m/s, respectively. By increasing the Reynolds number, the velocity gradient is increased at both sides of domain. For example at Re=10000, the maximum positive and negative velocities are increased from 0.17 m/s to 0.4 m/s and from −0.24 m/s to −0.62 m/s, respectively. So by increasing the Reynolds number, the velocity gradient of vertical velocity at both sides is increased. The reason for increasing the vertical velocity gradient is the same as that of horizontal velocity gradient.

Figures 7 to 10 show that the horizontal and vertical velocities on low Reynolds number on C.G and S.G are predicted accurately (e.g., for Re=100).

However increasing the Reynolds number increases the inconsistency between results so that the results of S.G are more close to Ghia results. These differences in the upstream and downstream of the domain are greater. While in the middle point, the results of C.G and S.G at various Reynolds numbers are predicted accurately.

Figure 11 illustrates the successive of eddies that are similar on C.G and S.G.

For low Reynolds number the stream lines on C.G and S.G are similar. However as the Reynolds number increases, the shape and the number of eddies change. For example at Re=10000 there is an eddy at the bottom corners of domain on S.G while there are two eddies at the same locations on C.G.

5.2. Channel Flow

In this section, the solution of the channel flow for a long channel with a fully developed laminar flow is presented. The length and height of the channel are selected as 𝐿=40 m and ℎ=1 m and the velocity in the entrance of the channel is 𝑢=1 m/s. The velocity at the end of the channel is calculated on C.G and S.G at various Reynolds numbers in a range of 100 to 1000. Figure 12 compares the results of C.G and S.G. It is evident that a very accurate result in comparison to analytical solution which is 𝑢=6(𝑌−𝑌²) is obtained.

(a)

(b)

(c)

(d)

Figure 12 

Channel test, comparison of fully developed profile on C.G and S.G and Analytical solution.

The extension of boundary layer on C.G and S.G at 𝑋=(0.2−0.8)∗𝐿 in Figure 13 shows that the results have nearly the same accuracy.

(a)

(b)

(c)

(d)

Figure 13 

Channel test, development of Boundary layer on C.G and S.G.

5.3. Flow over a Backward-Facing Step

To assess the accuracy of numerical method, the flow over a backward-facing step in a channel is a good test case. A dissipative scheme cannot predict the correct reattachment length of the recirculation zone downstream of the step [13]. The flow expands downstream of the step and reattaches on 𝑋_𝑟 coordinate. A two-dimensional schematic of the backward-facing step is shown in Figure 14. In this figure, 𝑑 is the height of the channel; ℎ is the height of the step and 𝑋_𝑟 is the reattachment length. In current work 𝑑 and ℎ selected as 1 m and 0.5 m, respectively.

Figure 14 

Geometry of step test.

The regime of fluid flow on Backward facing step, for Re<400, 400<Re<3400, and Re>3400 is laminar, transient and turbulent, respectively [7]. For Re<1000 there is an eddy in the downstream. But as the Reynolds number is increased, secondary or roof eddy is created in the upstream of the step and therefore the reattachment length is decreased [14].

In this section, the governing equations are solved and reattachment length at various Reynolds numbers on C.G and S.G is calculated. Finally the numerical and Armaly et al. results [8] are compared. Figure 15. shows this comparison. Armaly et al. results are regarded as the true solution of this problem.

Figure 15 

Comparison of Reattachment point on C.G and S.G and Armaly result.

It is evident from this figure that the results of S.G are better than C.G. But this result is not accurate for all the Reynolds numbers. For example for Re=5000, the result of S.G is very accurate but for Re=4000 the result of S.G and C.G are not consistent with Almary et al. results.

The reattachment length of the recirculation zone downstream of the step is decreased for Reynolds numbers more than 1000. This is because of the existence of second eddy at the top of step (see Figure 16).

Figure 16 

Step test, existence of second eddy at the top of step for .

6. Discussion and Conclusions

In this paper, three test cases as lid-driven cavity, channel flow and backward-facing step are considered to compare the accuracy of S.G and C.G methods. In the light of different results for N-S solver, following conclusion can be drawn:

In lid-driven cavity test, for low Reynolds numbers the results of S.G and C.G methods are nearly the same. Nevertheless for high Reynolds numbers, the results of S.G demonstrate better agreement. In lid-driven cavity test, for low Reynolds number the stream lines on C.G and S.G are similar. However as the Reynolds number increases, the shape and the number of eddies change. In the channel test, two grids present very good results for fully development and the extension of boundary layer. In backward-facing step, the reattachment length of the recirculation zone downstream of the step is decreased for Reynolds number more than 1000. This is because of the existence of second eddy at the top of step. In backward-facing step, the results of S.G are better than C.G.

It can be concluded from these test cases that, for problems such as lid-driven cavity and backward-facing step associated with eddy on the fluid flow, application of S.G leads to more accurate results. However for other problems such as channel flow, both methods nearly have the same accuracy.

Nomenclature

𝐶𝑠:	A constant in the range [0.025,0.65]
𝑃:	Static pressure
𝑆𝑖,𝑗:	Strain rate tensor
𝑆𝜙:	Source term in the 𝜙 direction
𝑡:	Time
𝑢𝑖	Filtered velocity in 𝑥 direction
𝑢𝑗:	Filtered velocity in 𝑦 direction
𝑢:	Velocity components in 𝑥 direction
𝑣:	Velocity components in 𝑦 direction
𝑉:	Velocity vector.

Greek Symbol

Δ:	Length of filter
Δ𝑥:	Mesh size in 𝑥 direction
Δ𝑦:	Mesh size in 𝑦 direction
Δ𝑧:	Mesh size in 𝑧 direction
𝑣𝑡:	Turbulence eddy viscosity
𝜇:	Dynamic fluid viscosity
𝜌:	Fluid density.