A Conservative Pressure-Correction Method on Collocated Grid for Low Mach Number Flows

A novel extension to SMAC scheme is proposed for variable density flows under low Mach number approximation. The algorithm is based on a predictor—corrector time integration scheme that employs a projection method for the momentum equation. A constant-coefficient Poisson equation is solved for the pressure following both the predictor and corrector steps to satisfy the continuity equation at each time step. The proposed algorithm has second order centrally differenced convective fluxes with upwinding based on Cell Peclet number while diffusive flux are viscous fourth order accurate. Spatial discretization is performed on a collocated grid system that offers computational simplicity and straight forward extension to curvilinear coordinate systems. The algorithm is kinetic energy preserving. Further in this paper robustness and accuracy are demonstrated by performing test on channel flow with non-Boussinesq condition on different temperature ratios.


Introduction
Various flow regimes in industrial devices are of low speed nature.Such flows are called incompressible, since the velocities are much smaller than the speed of sound.In non-reacting incompressible flows without heat transfer, the use of a pressure-correction algorithm has proven to be accurate and efficient (e.g.[1,2]).Since density remains constant, no substantial problems are encountered and the solution is straightforward.The mass conservation equation naturally imposes a constraint on the velocity field.However, if density varies strongly in time and space, e.g.due to temperature variation, the set of equations becomes more coupled and an efficient solution is no longer obvious.Various attempts have been made to create efficient solution methods.A basic difficulty stems from the acoustic waves in the compressible formulation.As acoustic waves act at a substantially smaller time scale than the convective phenomena in low Mach number flows, the acoustic modes do not significantly influence the solution and may be regarded as superfluous.The use of larger time steps, corresponding to the convective scales, can therefore strongly improve computational efficiency without loss of relevant information.
Furthermore conservation of kinetic energy in numerical methods has become an important issue in large eddy simulation (LES) and direct numerical simulation (DNS) of turbulence.Kinetic energy conservation in a finite difference formulation is not a consequence of discrete momentum and discrete mass conservation, so conservation of kinetic energy has to be ensured through careful design of the finite difference operators.It is known that dissipative numerical schemes (e.g.up-winding) often introduce too much artificial damping for use in turbulence simulations, because the energy balance in turbulence is rather delicate.In the case of variable density flows, not conserving the kinetic energy can also lead to erroneous temperature and density fields.Much work has been done in the development of kinetic energy conservation algorithms for incompressible flows (see Vasilyev [3]; Gullbrand [4]; Morinishi et al. [5]), but there has been less work on variable density or compressible flows (see Nicoud [6] and Lessani [7]).In low-speed turbulent channel flow applications, the low Mach-number, variable-density approximation of the Navier-Stokes equations is a good basis for simulation, as it supports large density variations while eliminating acoustic waves.This eliminates the need for extremely small time steps driven by the acoustics.This means that the arising velocities are much smaller than the speed of sound, so that density variations due to pressure variations can be neglected.In those so-called low Mach number flows, an efficient way to solve the set of Navier-Stokes equations describing the flow is to use a segregated solver, relying on a pressure-correction algorithm.Here, the pressure is split in a thermodynamic part P 0 and a second order kinematic pressure P 2 , which only appears in the momentum equations.As a result, the momentum equations together with a constraint on the divergence of the velocity decouple from the equations determining the density field.The velocity field is computed from the momentum equations, and is corrected with a pressure-correction to satisfy the divergence constraint.The correction on the pressure is the result of a Poisson-equation, which is elliptic.This paper is organized as follows.Section 2 shows the equations that govern low Mach-number flows, and in Section 3 some details of the numerical method and its implementation are shown.Finally, Section 4 contains test cases and numerical results.

Governing Equation
The low Mach-number approximation of the Navierstokes equations is obtained as the low Mach-number asymptotic limit of the compressible Navier-Stokes equations in which temperature fluctuations are assumed to be of order 1.In this analysis, the pressure is expanded as: In this expansion, P 0 is the spatially uniform thermodynamic pressure, and P 2 is the hydrodynamic pressure fluctuation.Details of the derivation of these equations can be found in Majda and Sethian [8]; Rehm and Baum [9]; Muller [10] and Paolucci [11].The final results of this process are the following equations (excluding the body forces): Conservation of mass: Momentum equation: Conservation of energy: For an open system, the thermodynamic pressure (P 0 ) does not change in time, but in a closed system (sealed enclosure) the thermodynamic pressure can change in time.Notice that the source terms from the energy equation impact the mass conservation equation through the constraint Equation (2.6).

Large Eddy Simulation (LES) Model
In LES, one computes the motion of large-scale structures, while modeling the non-linear interactions with the small-scales.The governing equations for large eddies are obtained after filtering.The filtering operation can be written in terms of convolution integral: Large Eddy Simulations have become an important tool for the study of turbulent transport in environmental and engineering flows as it requires coarser mesh than DNS.The basis of such a technique is the application of a spatial filter to the governing equations.An f turbulent variable is splitted into an f  large component and f  sub grid component.Note that~ corresponds also to the Favre average operator.The non-dimensional forms of the governing Equations (2.2)-(2.5)are then Favre averaged and filtered using implicit filter to obtain governing equations for filtered scale variables.
The flow in the channel is driven by the constant streamwise driving pressure force F.Here the superscript ~refers to the Favre averaged quantities and ij is the resolved strain rate tensor.Where Favre average is defined as In the present work, we consider several different sub grid-scale models for thermal part and kinetic part.For modeling sub-grid stresses, we have used Wall Adaptive Layer Equation (WALE) model suggested by Nicoud and Ducros [12].For thermal part, we have used a dynamic model.In WALE model subgrid scale viscosity accounts for effect of rotation rate and strain rate of the smallest fluctuations.This model has correct wall behavior for sub grid stesses near walls Read [13,14].WALE model approximates SGS eddy viscosity as where L s is a length scale given by where V is the volume of the cell, however 1 3 t V   , is von Karman's constant, z is the distance nearest to the closest wall, w = 0.325 is the wall constant and ij 0.4187 k  C  is the traceless symmetric part of the square of the velocity gradient tensor defined as.
For thermal part we use dynamic Smagorinsky model, the thermal dynamic coupling is taken into account through a similar procedure in order to estimate SGC turbulent Prandtl number.The heat flux after the filtering procedure corresponds to where SGS is given by Equation (2.16) and varies in space and time.ˆ2 2  where C dyn is the constant for dynamic Smagorinsky model and t  is the sub grid scale diffusivity.

Numerical Scheme for Channel Flow
The proposed numerical method is semi-implicit, pressure correction type scheme on a non-staggered structured grid using finite difference scheme for spatial discretisation.The scheme was described by Hirsch [15] and is conceptually similar to the SMAC algorithm described by Amsden and Harlow [16], guided by the work of Cheng and Armfield [17].It is conceptually similar but an extension to the existing scheme, which is made compatible for low Mach number flows.Here we take full Navier-Stokes equation and then removes acoustic modes from it, as acoustic waves act at a substantially smaller time scale than the convective phenomena in low Mach number flows, the acoustic modes do not significantly influence the solution and may be regarded as superfluous.

Temporal Descritisation
The flow field is marched forward in time using a two step predictor-corrector approach.In the predictor step, the time integration of momentum equation is performed using a first order Euler method to obtain the guessed velocity field at the next time step.In the corrector step the guessed velocity fields close to zero.Finally the scheme is given as: and guessed velocity vector are calculated as , where and superscripts n, * , and n + 1 denote the known values at the time level n, the predicted or guessed fields and the values at the new time level n + 1 respectively.The guessed velocity fields do not necessarily satisfy continuity equation.Here, , , is the sum of convective and diffusive fluxes while is the sum of convective and heat transfer fluxes at time level "n".The predicted value for the temperature T * is calculated from Equation (2.4) based on previous values at time level n,   Guessed thermodynamic pressure and guessed de fie ) nsity ld would be In this step guessed velocity field obtained in the predictor step is corrected in a continuity preserving manner, firstly we define correction flux and correction pressure The velocity and temperature field at tim th 1 , , , , e level n + 1 at satisfies continuity is coupled to the pressure field at anew time level n + 1 through a semi explicit discretisation in time of the momentum equations.This is given by,   , The relationship between velocity and pressure co tion (P') can be obtained by subtracting Equation fr Taking divergence of the above the discrete divergence operator defined as below Equation (3.7) with .
Now using continuity equation and defined correction flux we get 1 2 where below mentioned value of density is used son equation For a closed system, like ours, the th odynamic pressure p o is calculated by using erm 0 0

Spatial Descritisation
The spatial discretization is performed using finite difcated mesh with a Carte--winding scheme employs ference methodology on a collo sian coordinate system.The up two points in upstream and one point on the downstream side of the grid under-consideration.Convective fluxes are second order centrally differenced; the choice between up-winding and central differencing is made on the basis of cell Peclet number while diffusive fluxes are viscous fourth order accurate.The discrete Poisson equation is differenced as Convective flux are second order centrally differenced given below with averaged face density values: y averaging: 2 The node velocities are calculated b Diffusive fluxes are viscous fourth order accurate which is calculated as follows: 1, ,  , ,  , , 1 Re

Results and Discussion
Three numerical experiments were performed in o ts implementation.We for TLES (thermal large rder to test the numerical scheme and i used self made FORTRAN code eddy simulation).A mesh of size 96 × 96 × 96 is used for simulation purpose, uniform meshes are used in the stream-wise and span-wise directions and a non-uniform mesh with hyperbolic tangent distribution is used in the wall-normal direction.For the velocities, no slip boundary condition was used on the top and bottom walls and periodic boundary condition was used in x-y direction.The temperatures on the hot and cold walls were set according to R  .Viscosity was calculated using Sutherlands law with reference at T c .The reference Reynolds number based on the reference velocity and channel half width is kep onstant at 2800 while the Prandtl number is fixed at 0.71.
First, a 1D convection-diffusion problem was set to test spatial and temporal convergence in a variable density case.Second t c , a 2D inviscid solenoidal velocity field w lem for this experiment is a temperature m-main was the interval [0, 1].The initial temperature pro-as used to test kinetic energy conservation.Finally, a 3D turbulent channel flow with large temperature gradient is used.

Spatial and Temporal Convergence in 1D
The test prob profile that will be convected in the x-direction.For si plicity the flow was assumed to be inviscid and the do file is a Gaussian given by equation.
  2 0.5 295 50 exp 0.05 The velocity in the x-direction (u) was set equal to 1; the hydrodynamic pressure was set to 0 Pa and t modynamic pressure was set to 1.The work assumed to be air so that its thermal conductivity (k) co l steady quation (4.2). he thering fluid was uld be estimated from polynomial correlations.Since the thermal diffusivity of air is O (1 × 10E−05), the global Peclet number is very high O (1 × 10E−05).For the spatial convergence the grid was changed from 100, 200, 400 to 800 nodes while the time step was held in 0.00125 so that the CFL number changed from 0.125 to 1.For the temporal convergence the mesh was held in 100 grid points and the time was changed from 0.005 to 0.000625 so that the CFL number changed from 0.5 to 0.0625.Figure 1 shows the profile of velocity at t = 1 for four different time steps.The velocity induced by the diffusion process was of the order 1 × 10 −6 .The shape of the u-velocity profile is in agreement with the theory and what is expected from Equation (4.1).Table 1 summarizes the numerical results from this experiment.

Kinetic Energy Conservation
For this numerical experiment a 2-dimensional rectangular domain [0, 1] × [0, 1] is used, with an initia state solenoidal velocity field given by E , sin 2π cos 2π 2) The temperature field was set as a Gaussian random field with a mean value of 397 K and a root mean square fluctuation of 57 K.The density field can be the equation of state.Using Equation (4.2) and the initial computed from random density field, the initial kinetic energy can be computed (KE 0 = 0.2274 J).A mesh of 24 × 24 points was used, so that Δx = Δy = 4.2e−02.According to Nicoud [6], the integration time for this numerical experiment is given by Equation (4.3).Table 2 and Figure 2 show the results for this experiment.It is evident that the scheme conserves global kinetic energy, so that the divergence-free constraint is recovered in the inviscid limit.
The error in Table 2 was calculated using the ratio between the difference in the initial and final kinetic energy i.e. this section, we present detailed numerical results from the test problems that we considered in order to check the robustness and accuracy of the proposed algorithm.This is the case of large-eddy simulation (LES) of non-isothermal, turbulent channel flow with strong temperature gradients due to the temperature difference between the two walls.
In this n t normal directions, respectively.The d domain are 4 πδ × 4/3 πδ × 2δ, with width of the channel.The walls of the channel are normal to the z direction and are held at constant temperature.The boundaries of the domain normal to the x and y directions are periodic in nature.Therefore, the total mass of the system is conserved, i.e., this is an example of flow in a closed domain.A mesh of 64 3 points is used in such a way that Δx + × + + niform meshes are used in the stream-wise and spanise directions and a non-uniform mesh with hyperbolic  tangent distribution is used in the wall-normal direction.
The Wall-Adapting Local Eddy-viscosity (W-ALE) model (Nicoud [12]) is used for modeling the eddy viscosity and dynamic Smagorinsky method is used to model turbulent heat flux.In the homogeneous directions, the convective terms of the momentum and energy equations are calculated using hybrid type upwind.The multi-grid algorithm (Jameson [18]) is used to solve the constant-coefficient pressure Poisson equation.Let  T T = 1.01) and Sutherland law for the case ( h c T T = 2).The R  = 1.01 results are compared with the DNS of Kim et al. [19] and Lessani [20], while R  = 2 results are compared with DNS of Nicoud [21] as shown in Table 3.
Figure 3 shows that R θ = 1.01 is in good agreement with previous incompressible DNS (Kim et al., [19]) for the mean ocity profile.T e expected (fo the f n Reynolds number consid e-wall u + = 2.5 ln(y + ) + 5.5 is recovered.The viscous sub layer is also well resolved as shown in Figure 3. Figure 4 shows a good agreement with previous incompressible DNS (Kim et al., [19]) for the three velocity fluctuations and the Reynolds shea stess.For any instantaneous variable   5 shows mean velocity profile and Figure 6 shows temperature profile for the entire channel.Both are in good agreement with previous DNS of Nicoud [21].Figure 7 shows Mean velocity profile (u + ) with the classical scaling.For the cold side of the channel, law-of-the-wall u + = 2.5 ln (z + ) + 5.5 is recovered but with a different constant while in the hot side of the logarithmic nature of law-of-the-wall is small.This agrees with the previous DNS.The viscous sub layer is also well resolved as shown in Figure 7. Figure 9 shows good agreement of RMS velocities along the wall norm direction.It can be seen that agreement with the pubold channel than e hot channel.There is resent work and DNS.The discrepancy is alerature fluctuations along the wall normal direction.There is a qualitative agreement while the peak fl tuating intensities don't match.There is deviation in the core of the channel s can be due elocities obtained.Further as the comparis is from the DNS which is far more accurate than LES.al lished literature is much better in the c th a good qualitative agreement between p ways less than 5%.Figure 8 shows temp uc and the hotter channel side.These deviation to different bulk and channel centerline v on

Conclusion
In this article, a new algorithm which extends the exist S. M. YAHYA ET AL. 260    particularly useful for unsteady flows with strong temperature gradients.A constant-coefficient Poisson equation, which is computationally more efficient than the variable-coefficient one, is solved for the pressure.A collocated grid is used for spatial discretization and not a ty and straightforward extension to curvinear coordinate systems.The odd-even decoupling problem is avoided by using smartly a continuity equation in conjunction with correction flux.The robustness and accuracy of the algorithm have been assessed through simulations of three test problems; 1D convection-diffusion problem, 2D steady state solenoidal velocity distribution for conservation of kinetic energy and finally the turbulent channel flow with temperature gradients.The results obtained with the proposed algorithm are in very good agreement with the ones reported in earlier studies.The algorithm has a moderate computational cost for solving the Poisson-equation and is stable for high density ratios (at least up to a factor of 10).Though in the examples shown in this paper, P 0 was always a constan pr )

. 12 )
Copyright © 2012 SciRes.WJMwhere the subscripts R and L indicate extrapolated values at the right and left face of the control volume.F order upwinding, and positive values of the velocity,
consider turbulent flow in a channel whose walls are kept at differe t temperatures.This problem is treated numerically with he help of LES.Let x, y and z, denote the stream-wise, span-wise and imensions of the δ being the half Δy × Δz = 33 × 11 × [05:11].U KE 

Figure 2 .
Figure 2. Error of the KE as a function of CFL.
T h and T c denote the temperatures of the hot and cold walls, respectively.Two different cases, corresponding to diffe ent w Speifically, rall temperature ratios, are considered herein.c R  = 1.01 and R  = 2. Here, R  = h c T initialize the flow field.

Figures 5 - 9
Figures5-9show variations of different quantities in the channel for R θ = 2. Figure5shows mean velocity profile and Figure6shows temperature profile for the entire channel.Both are in good agreement with previous DNS of Nicoud[21].Figure7shows Mean velocity profile (u + ) with the classical scaling.For the cold side of the channel, law-of-the-wall u + = 2.5 ln (z + ) + 5.5 is recovered but with a different constant while in the hot side of the logarithmic nature of law-of-the-wall is small.This agrees with the previous DNS.The viscous sub layer is also well resolved as shown in Figure7.Figure9shows good agreement of RMS velocities along the wall norm direction.It can be seen that agreement with the pubold channel than e hot channel.There is resent work and DNS.The discrepancy is alerature fluctuations along the wall normal direction.There is a qualitative agreement while the peak fl tuating intensities don't match.There is deviation in the core of the channel s can be due elocities obtained.Further as the comparis is from the DNS which is far more accurate than LES.

Figure
Figures5-9show variations of different quantities in the channel for R θ = 2. Figure5shows mean velocity profile and Figure6shows temperature profile for the entire channel.Both are in good agreement with previous DNS of Nicoud[21].Figure7shows Mean velocity profile (u + ) with the classical scaling.For the cold side of the channel, law-of-the-wall u + = 2.5 ln (z + ) + 5.5 is recovered but with a different constant while in the hot side of the logarithmic nature of law-of-the-wall is small.This agrees with the previous DNS.The viscous sub layer is also well resolved as shown in Figure7.Figure9shows good agreement of RMS velocities along the wall norm direction.It can be seen that agreement with the pubold channel than e hot channel.There is resent work and DNS.The discrepancy is alerature fluctuations along the wall normal direction.There is a qualitative agreement while the peak fl tuating intensities don't match.There is deviation in the core of the channel s can be due elocities obtained.Further as the comparis is from the DNS which is far more accurate than LES.

Figure 5 .
Figure 5. Mean Stream-wise velocity profile normalized with maximum velocity for the entire channel.

Figure 6
Figure 6.Mean Temperature for  c T T T Δ

Figure 7 .Figure 8 .
Figure 7. Mean velocity profile in log units along the channel for R θ = 2. ing SMAC scheme for low Mach number, variable density flows has been presented.This algorithm can be applied to both open and closed domains.It is based on two-stage predictor-corrector method.This algorithm is a

Figure 9 .
Figure 9. Variation of RMS velocities along the wall normal direction (z) for R θ = 2.