Guoneng Li*,Youqu Zheng,Huawen Yang,Wenwen Guo,Yousheng Xu

Department of Energy and Environment System Engineering,Zhejiang University of Science and Technology,Hangzhou 310023,China

1.Introduction

A jetin cross flow(JICF)is one of the classical flow problems.It is encountered in many engineering applications including cooling jets in the turbine blades,fuel jets in combustion systems,industrial mixing processes,and pollutant transportations.Extensive previous experimental studies,theoretical investigations and numerical simulations are reported,exploring different aspects of the fluid flow characteristics,such as the central trajectory[1–4],mixing feature[5–8],vortex evolution[9–12],jet stability[13–16],etc.Detail information and discussions can be found in review papers[17–19].For the conciseness,only selected journal papers regarding to numerical simulations of non-reactive and single phase JICF with low Mach number in the recent five years are reviewed in this work.

In general,convectional Computational Fluid Dynamics(CFD)simulations are basically divided into three categories,i.e.Reynolds Averaged Navier Stokes(RANS)simulation,Large Eddy Simulation(LES),and Direct Numerical Simulation(DNS).Recent CFD studies continue to employ above methods to study JICF.Ivanova et al.[20]found that the unsteady RANS is more accurate than steady RANS,while Javadi and El-Askary concluded that LES gives better predictions than RANS[21].Other investigations showed that LES results in the immediate jet-tocross- flow exhaustion zone still have discrepancy from experimental results[22],but LES is capable of reproducing the coherent vortex structures[23].Recent LES studies found that the well-known transverse jet vertical vortex is close related to the destabilization and evolution of the leading-edge shear-layer rollup[24],and that the mixing between the transverse jet flow and the cross flow can be enhanced considerably as flow pulsations are added to the transverse jet with proper pulsed signal and duty cycle[4].DNS plays an important role in understanding the stability of the transverse jet.Studies show that the primary instability is originated from the destabilization of the leeside shear-layer vortex[16,25],and this kind of instability can be self-sustained under certain jet-to-cross- flow velocity ratio(r),which is close related to the oscillation of upright wake vortex[13].Apulsating transverse jet was also studied in detail with DNS method[3,26],indicating that both the central trajectory and the spread width of the transverse jet can be optimized in order to enhance the mixing process.

Apart from the convectional CFD to study the flow mechanism of JICF,other numerical methods also could be explored to provide engineers with various tools for the prediction of the flow mechanism.One of these alternative methods is the Lattice Boltzmann Method(LBM),which was proven to be a successful computational method in many kinds of physical and/or chemical problems,such as turbulence flow[27],nano- fluid heat transfer[28],mass convection[29–32],and droplet wetting[33].However,very limited literatures on LBM simulation of a JICF system were found.Therefore,an initial and preliminary attempt to develop a LBM model to simulate the fluid flow and the mixing process of a laminar jet in cross flow is presented in this work.On the other hand,this work is a first step to model the combustion and heat transfer process involved a JICF problem using LBM,such as those of our previous works[34,35].Besides,it should be noted that using LBM to simulate some complicated problems that cannot be readily handled by the convectional CFD method is another advantage of LBM,such as modeling the nano-particle behavior in the nano- fluid heat transfer[28].

2.The Lattice Boltzmann Method

2.1.Computational zone

Fig.1 shows the computational zone and the dimensions.A square jet(d=3 mm)with a length of 3d is placed in the bottom wall of the cross flow channel.A moderate Reynolds number,Re=400,based on the cross flow velocity,Uc,and the cross flow hydraulic diameter,D=10d,is selected.The-jet-to-cross- flow velocity ratio,r=Vj/Uc,is set to be 5.0 and 7.5,where Vjis the jet velocity.In the present work,uniform lattice grids with a resolution of Lr=0.1875 mm,i.e.400×160×192(12288000 grids in total),is used.Air is used as the fluid.The kinematic viscosity is v=15.06×10-6m2·s-1,and the Prandtl number is Pr=0.707.The Mach number is low enough to neglect the compressibility of the fluid.

Fig.1.Computational zone of the JICF system.

2.2.LBM model

The basic form of the Lattice Boltzmann Equation(LBE)with an external force for incompressible flows by introducing Bhatnagar-Gross-Krook(BGK)approximations can be written as follows[36]:

where δt is the lattice time step size.τfis the dimensionless collisionrelaxation time for the fluid flow.Fαis the applied force,such as the gravitational force,magnetic force etcis the direction vector,which is shown in Fig.2 for the D3Q19 model shown in Fig.2,defined as

where c is the lattice velocity.fαin Eq.(1)is the Particle Distribution Function(PDF).is the Equilibrium Distribution Functions(EDF)for fα,which is given as

Fig.2.Discrete velocity set of three-dimensional nineteen-velocity(D3Q19)model.

where ρ andare macroscopic density and velocity vector,respectively.where u,v and w are the x,y and z direction velocity component,respectively.wαis a constant,i.e.w0=6/36,wα=2/36 forα=1,2,…,6 and wα=1/36 forα=7,8,…,18.The macroscopic density and velocity vector are calculated by

The Chapman-Enskog expansion for the PDF can recover the continuity equation and the Navier-Stokes equation.The detailed derivation of these equations is given by Hou et al.[37],and will not show here.The lattice viscosity vlbis determined by

where csis the lattice sound speed.Once the macroscopic velocity is chosen,the lattice viscosity can be determined based on the same Reynolds number,which is connected to the physical velocity and dimensionless viscosity.Then the flow collision-relaxation time can be determined.The choice for the viscosity from Eq.(5)makes the above scheme formally a second order method for solving incompressible flows.The positivity of the viscosity requires that τf＞0.5.Eq.(5)can be solved in the following two steps:

wheredenotes the post-collision state of the PDFs.It is noted that Eq.(6)is explicit and easy to implement.

2.3.Boundary conditions

2.3.1.Velocity inlet boundary

The type of Neumann boundary condition(Velocity inlet)constrains a fixed flux at the boundary of a D3Q19 model.The unknown PDFs for the left side of the cross flow inlet at particle direction of 1,7,8,9 and 10 are calculated following the method of Zou and He’s algorithm[38],which is an algorithm for D2Q9.The extending method from a D2Q9 algorithm to a D3Q19 one is briefly introduced as follows[36],

Eq.(7-1)to Eq.(7-4)are totally 4 equations,whereas the there are 6 unknowns,i.e.ρ,f1,f7,f8,f9,and f10.Therefore,we need two more equations.To keep the symmetry of the problem,the bounce back of the nonequilibrium part for all populations are employed,and two new transverse momentum corrections,i.e.are introduced to maintain the conservation of momentum.This idea gives the following equations,

The unknown PDFs at particle direction of 7,8,9 and 10 can be derived easily with Eq.(7-5)to Eq.(7-8)and the Eq.(8-3)to Eq.(8-4).

2.3.2.Outlet boundary

A fully developed outlet boundary condition was adopted,in which the unknown PDFs and EPDFs at particle directions of 2,11,12,13 and 14 are forced to equal to the PDFs at the nearest node in the–x direction,as follows,

where the subscripts are corresponding to each other between the left side and the right side of Eq.(9),e.g.f2(L,y,z,t)is corresponding to f2(L?Lr,y,z,t).

2.3.3.Wall boundaries

The bounce back algorithm is used in the walls to obtain the unknown PDFs[39,40].In this work,the bounce back algorithm is slightly modified.For example,the bounce back algorithm for the top wall boundary is as follows,

The bounce back algorithm is simple,and second-order accuracy can be achieved when the boundaries are located right in the middle of two neighboring lattices[39].Therefore,special efforts was taken to calculate the physical velocities after the calculation of equation(4-1)and(4-2).In the following section,the accuracy of present bounce back algorithm is compared with the standard bounce back method[40].The standard bounce back algorithm for the top boundary is as follows,

where the subscripts has the same rule as Eq.(9).In other words,the standard bounce back algorithm adopted the PDFs in opposite directions and in the previous time step to predict the unknown PDFs in the present time step.

2.3.4.Edge treatment

There is a problem while employing the algorithm of Eq.(10)or Eq.(11),i.e.the modified bounce back algorithm and standard bounce back algorithm.The problem is shown in Fig.3.It seems that the edges between two walls are opened for the fluid to flow outside.The reason is that the required unknown PDFs are not predicted correctly in these nodes,i.e.only the PDFs pointing into the fluid zone are predicted in the edges,as a result,the velocity in these edge nodes is not zero.In this work,an enhanced treatment for these bounce back algorithm is proposed.Take the edge(0,0,a)-(0,L,a)for example,

where the subscripts for Eq.(12-2)has the same rule as Eq.(9).The effect before and after employing this enhanced treatment is shown in Fig.3.

2.4.Convergence criterion

The convergence criterions for the fluid flow is set to be

where n and n-1 denotes the present time step and the previous time step.

2.5.Parallel the calculation using MPI programming

In order to accelerate the calculation,the Message Passing Interface(MPI)is employed to the above code.No other errors were introduced through using the MPI method.The developed MPI code can be used in any number of CPUs.

Fig.3.The Counter rotating Vortex Pair(CVP)before and after the enhanced treatment of the bounce back algorithm.

2.6.Implementation of the computational zone in LBM

In order to implement the computational zone shown in Fig.1,a function,Tag(void),is used to identify the fluid grids.Tag(void)=0.0 when the grids are inside the computational zone,otherwise the Tag(void)=1.0.For those grids with Tag(void)＞0.5,the collision and streaming processes,i.e.Eq.(6-1)and Eq.(6-2)will not be performed.

3.Validation of the Code

3.1.Validation of the code in Poiseuille flow

In order to validate the present code,the fluid flow inside a square channel is simulated.A uniform velocity i.e.u(y,z)=uav=0.1 m·s-1(v=w=0.0 m·s-1),is used in the velocity inlet.The fluid flow will develop into a fully developed velocity distribution after a length larger than the entrance length.The entrance length Lecan be estimated as[41]

Leis estimated to be 3.77aLr.In the present work,the velocity distribution in the crossing section of 4.0 aLrwas used to compare with the exact solution.This validation method was used in a previous work[39].The fully developed Poiseuille velocity inlet in the square channel is as follows[41],

Fig.4 shows the comparison between the numerical results and the exact solutions calculated for Eq.(15-1)to Eq.(15-3).The grid resolution is Lr=0.204 mm,resulting in 250×50×50 grids.which corresponding to physical time step of 5.0×10-5s.The maximum axial velocity(u)is located at the center of the channel.The exact solution is 0.20963 m·s-1,while the modified bounce back numerical result is 0.20834 m·s-1,whereas the standard bounce back numerical result is 0.20169 m·s-1.The errors are 0.62%and 3.79%,respectively.Therefore,the modified bounce back algorithm increases the numerical accuracy.On the other hand,the numerical velocity near the wall from the standard bounce back algorithm is not zero,and this was improved by the modified bounce back algorithm.The errors of mass conservation of the modified bounce back algorithm and the standard bounce back algorithm are 0.39%and 3.71%,respectively.It should be noted that the exact results in Fig.4 has an error of 10-6,which is set artificially in the program,and it should be small enough since the convergence criterions of the LBM method is also 10-6.

The effect of grid resolution Lron the predicted results in the D3Q19 model with modified bounce back algorithm is shown in Fig.5(a).The predicted results are getting close to the exact solutions when increasing the grid resolution.A grid resolution of Lr=0.204 mm(250×50×50)is fine enough to predict the velocity field for the present work,and further fine grids is not necessary.The errors of mass conservation for Lr=0.345 mm(150×30×30)and Lr=0.256 mm(200×40×40)are 4.48%and 1.47%,respectively.The selection of dimensionless collision-relaxation time is also important in predicting fluid flow using a D3Q19 method with BGK model.As shown in Fig.5(b),decreasing the collision-relaxation time has a slightly effect on the predicted results,especially in the peak velocity.The predicted peak velocities are 0.20552 m·s-1and 0.2054 m·s-1for,respectively.These result in errors of 1.96%and 2.02%,respectively.The converging process becomes more and more difficult when the collision relaxation time continues to be lower values.An obvious conclusion can be obtained by comparing Fig.4 and Fig.5.The modified bounce back algorithm is better than the standard bounce algorithm,especially in the near wall region.The modified bounce back algorithm gives much more reasonable results in these regions even with lower grid resolutions.

3.2.Validation of the code in JICF

Fig.4.Comparison between numerical results and exact solutions(τf=0.70331).

The present code is also validated using previous correlations in a JICF system[18,19].The velocity contours and the stream traces in the central plane(y=5d)for the cases r=5.0 and r=7.5 are shown in Fig.6.The central trajectories are taken as the lines with the largest velocity magnitude in the central plane,which are shown in Fig.7.The fitting curves are also presented in Fig.7,which are used to validate the present developed code.For the JICF under r=5.0,the central trajectory is y*=1.35x*0.34,where y*=(z-5d)/(rd),x*=(x-5d)/(rd).The correlation coefficient R equals to 0.994,indicating very good correlation between the fitting curve and the experimental data.For the JICF under r=7.5,the central trajectory is y*=1.28x*0.34 with R=0.995.The fitting equations are consistent with previous correlations,which is obtained from many experimental works by two review papers[18–19],i.e.y*=A x*B,where 1.2＜A＜2.6,0.28＜B＜0.34.Therefore,the developed code can be used to study the flow mechanism of JICF.It should be noticed that the effect of the dimensions of the computational domain(30d×10d×10d)on the central trajectory was not explored.Furthermore,no turbulent models were employed in the present work since the Reynolds number is 400,which indicates the flow is in laminar flow regime.Basically,the developed LBM model is a second order accurate method since Eq.(1)is used to represent the Boltzmann equation.For turbulent JICF,the present LBM code is ready to be extended to include a turbulent model through adding a turbulent viscosity in Eq.(5)[27].

4.Vortex System and Mixing Characteristics

The vortex system is one of the special features for a JICF problem[9–12],which dominants the whole near wake flow region and contributes to the blending of the transverse jet.The vortex system is difficult to be captured unless DNS method[16,25]or LES method with very fine grid resolutions[22,23]are employed.The upright wake vortex is shown in Fig.6,and three dimensional steam traces are shown in Fig.8.As shown in Fig.6 and Fig.8,the upright wake vortex is very import ant in the spreading the transverse jet in the leeside,i.e.this upright fluid flow penetrates into the transverse jet core and pushes the jet to spread out in the±y direction.Meanwhile this process weakens the stiffness of the transverse jet and contributes to the blending the jet.The upright wake vortex originates from the flow entrainment process.The blocking effect of the transverse jet results in a low pressure zone located in the leeside of the jet,where the incoming flow in the cross flow supplies the fluid to this low pressure zone,and flows in the–x direction,then flows up along with the transverse jet.

Fig.5.Effects of the grid resolution L r and effects of the collision-relaxation time.

Fig.6.Contours of the velocity magnitude in the plane of y=5d and the wake vortex.

Fig.7.Central trajectories in the plane of y=5d.

Fig.8.The entrainment effect for the incoming flows.

Fig.9.Counter rotating Vortex Pair in the plane of x=0.02 m.

Another famous coherent vortex structure is the CVP in a JICF problem,which is shown in Fig.9 for the cases with r=5 and r=7.5 in the pane of y=0.02 m(y*=1.67).As shown in Fig.9,the CVPs are clearly identified,and different positions and spreading widths of the CVP can be observed under different jet-to-cross- flow velocity ratios.The central point of the CVP for the case of r=7.5 is higher than the other one since the transverse jet under this case has stronger flow stiffness.On the other hand,the spread width at r=7.5 is smaller than the other one since the blending of transverse jet at this case is weaker than the other case.The CVP plays an important role in the mixing between the transverse jet and the cross flow.A detail data exploration can be seen in Fig.10.

As shown in Fig.10,stream trances at different layers in the height of the incoming cross flows are created,and stream trances in a surface which is perpendicular to the above layers are injected form the transverse jet orifice.It is clearly shown that the incoming fluid and the transverse jet fluid are mixing well in the far field.As shown in Fig.10,it seems that the flow injected fromthis layer with r=7.5 spreads wider than that with r=5.0.However,It's clearly shown in Fig.9 the whole spread width of the transverse jet with r=5.0 is larger than that with r=7.5.Therefore,the flow from different layers have different spread width.Further investigations are needed to find out the relationship between the spread width of the transverse jet and the jet-to-cross- flow velocity ratio,since very limited cases are studied in the present work.An interesting phenomenon can be observed in the near field in the aspect of mixing between the transverse jet fluid and the fluid from different layers of the incoming cross flows.A layer-organized entrainment pattern was found indicating that the incoming fluid at lower position is firstly entrained into the leeside of the jet,and followed by the incoming fluid at the upper position.This phenomenon should be taken into consideration when flow visualization experiments are designed.

5.Conclusions

An incompressible D3Q19 LBM model based on the single relaxation time algorithm was developed to simulate the fluid flow of a JICF problem.The code was validated by the mathematic solution of a Poiseuille flow and by the previous well studied empirical correlations for the JICF problem.Several conclusions are drawn as follows,

(1)The developed LBM modelis able to capture the dominant vortex in the JICF system,i.e.the Counter-rotating Vortex Pair(CVP)and the upright wake vortex.Grid independence test and relaxation time independence test showed that the developed code is numerical stable,and the present code is paralleled showing high computation efficiencies.

Fig.10.Flow mixing between the jet and the cross flow.

(2)Results show that the incoming fluid in the cross flow channel is entrained into the leeside of the jet fluid,which contributes to the blending of the jet,and that the spread width of the jet flow decreases with the velocity ratio under the present studied parameter range.

(3)Alayer-organized entrainment pattern was found indicating that the incoming fluid at lower position is firstly entrained into the leeside of the jet,and followed by the incoming fluid at the upper position.The layer-organized entrainment pattern is important when the mixing is under consideration.

Acknowledgements

We are grateful for all kind suggestions and discussions on the code development by Professor Yildiz Bayazitoglu in RICE University.

[1]S.Muppidi,K.Mahesh,Study of trajectories of jets in cross flow using direct numerical simulations,J.Fluid Mech.530(5)(2005)81–100.

[2]R.Sau,K.Mahesh,Optimization of pulsed jets in cross flow,J.Fluid Mech.653(6)(2010)365–390.

[3]F.Muldoon,S.Acharya,Direct numerical simulation of pulsed jets-in-cross flow,Comput.Fluids 39(10)(2010)1745–1773.

[4]A.Coussement,O.Gicquel,G.Degrez,Large eddy simulation of a pulsed jet in crossflow,J.Fluid Mech.695(1)(2012)1–34.

[5]S.H.Smith,M.G.Mungal,Mixing,structure and scaling of the jet in cross flow,J.Fluid Mech.357(2)(1998)83–122.

[6]J.Y.Fan,S.L.Xu,D.Z.Wang,PDA measurements of two-phase flow structure and particle dispersion for a particle-laden jet in cross flow,J.Hydrodyn.22(1)(2010)9–18.

[7]F.Yao,T.Wang,C.Gu,Experimental investigation on dilute gas-solid multiphase jet in cross flow,J.Mech.Sci.Technol.25(6)(2011)1483–1493.

[8]Z.W.Li,W.X.Huai,Z.D.Qian,Study on the flow field and concentration characteristics of the multiple tandem jets in cross flow,Sci.China Tech.Sci.55(10)(2012)2778–2788.

[9]A.Sau,T.W.H.Sheu,R.R.Hwang,W.C.Yang,Three-dimensional simulation of square jets in cross- flow,Phys.Rev.E 69(6)(2004)066302.

[10]Y.M.Marzouk,A.F.Ghoniem,Vorticity structure and evolution in a transverse jet,J.Fluid Mech.575(3)(2007)267–305.

[11]S.Muppidi,K.Mahesh,Direct numerical simulation of passive scalar transport in transverse jets,J.Fluid Mech.598(3)(2008)335–360.

[12]F.Schlegel,D.Wee,Y.M.Marzouk,A.F.Ghoniem,Contributions of the wall boundary layer to the formation of the counter-rotating vortex pair in transverse jets,J.Fluid Mech.676(6)(2011)461–490.

[13]Bagheri S,Schlatter P,Schmid P J,Henningson D S.Global stability of a jet in cross flow.J Fluid Mech,2009,624(1):33-11

[14]J.Davitian,D.Getsinger,C.Hendrickson,A.R.Karagozian,Transition to global instability in transverse-jet shear layers,J.Fluid Mech.661(10)(2010)294–315.

[15]M.Ilak,P.Schlatter,S.Bagheri,M.Chevalier,D.S.Henningson,Stability of a jet in cross flow,Phys.Fluids 23(9)(2011)091113.

[16]M.Ilak,P.Schlatter,S.Bagheri,D.S.Henningson,Bifurcation and stability analysis of a jet in cross- flow:onset of global instability at a low velocity ratio,J.Fluid Mech.696(4)(2012)94–121.

[17]R.J.Margason,"Fifty years of jet in cross- flow research",In:Computational and experimental assessment of jets in cross- flow,Winchester,UK,1992 34–75.

[18]A.R.Karagozian,"Background on and applications of jets in cross- flow",In:Manipulation and control of jets in cross- flow,Udine,Italy,2001 3–13.

[19]S.Y.No,Empirical correlations for breakup length of liquid jet in uniform cross flow-A review,Korea Inst.Liq.Atomization 18(1)(2013)35–43.

[20]E.M.Ivanova,B.E.Noll,M.Aigner,Computational modeling of turbulent mixing of a transverse jet,J.Eng.Gas Turbines Power 133(2)(2011)021505.

[21]A.Javadi,W.A.El-Askary,Numerical prediction of turbulent flow structure generated by synthetic cross-jet into a turbulent boundary layer,Int.J.Numer.Methods Fluids 69(7)(2012)1219–1236.

[22]D.Cavar,K.E.Meyer,LES of turbulentjet in cross- flow:Part1-A numerical validation study,Int.J.Heat Fluid Flow 36(1)(2012)18–34.

[23]D.Cavar,K.E.Meyer,LES of turbulent jet in cross- flow:Part 2-POD analysis and identification of coherent structures,Int.J.Heat Fluid Flow 36(1)(2012)35–46.

[24]G.Bidan,D.E.Nikitopoulos,On steady and pulsed low-blowing-ratio transverse jets,J.Fluid Mech.714(1)(2013)393–433.

[25]P.Schlatter,S.Bagheri,D.S.Henningson,Self-sustained global oscillations in a jet in cross flow,Theor.comput.Fluid Dyn.25(1-4)(2011)129–146.

[26]R.Sau,K.Mahesh,Optimization of pulsed jets in cross flow,J.Fluid Mech.653(6)(2010)365–390.

[27]L.Jahanshaloo,E.Pouryazdanpanah,N.A.C.Sidik,A review on the application of the lattice Boltzmann method for turbulent flow simulation,Numer.Heat Transf.A 64(11)(2013)938–953.

[28]G.N.Li,M.Aktas,Y.Bayazitoglu,A review on the discrete Boltzmann model for nanofluid heat transfer in enclosures and channels,Numer.Heat Transf.B 67(6)(2015)463–488.

[29]A.Atia,K.Mohammedi,Pore-scale study based on lattice Boltzmann method of density driven natural convection during CO2injection project,Chin.J.Chem.Eng.23(10)(2015)1593–1602.

[30]M.Nazari,L.Louhghalam,M.H.Kayhani,Lattice Boltzmann simulation of double diffusive natural convection in a square cavity with a hot square obstacle,Chin.J.Chem.Eng.23(1)(2015)22–30.

[31]Chen W.,Chen S.Y,Yuan X.G.,Zhang H.,Yu K.,"Three-dimensional simulation of interfacial convection in CO2-ethanol system by hybrid lattice Boltzmann method with experimental validation",Chin.J.Chem.Eng.,23(2):356-365(2015).

[32]B.Fu,X.G.Yuan,B.T.LIU,S.Y.CHEN,H.S.Zhang,A.W.Zeng,G.C.Yu,Characterization of Rayleigh Convection in Interfacial Mass Transfer by Lattice Boltzmann Simulation and Experimental Verification,Chin.J.Chem.Eng.19(5)(2011)845–854.

[33]Y.M.Yong,S.Li,C.Yang,Transport of Wetting and Nonwetting Liquid Plugs in a T-shaped Microchannel,Chin.J.Chem.Eng.21(5)(2013)463–472.

[34]H.Zhou,Q.Tang,T.Ren,G.N.Li,K.F.Cen,Control of thermoacoustic instability by CO2and N2jet in cross- flow,Appl.Therm.Eng.36(1)(2012)353–359.

[35]G.N.Li,Y.Q.Zheng,G.L.Hu,Z.G.Zhang,Experiments on convection heat transfer enhancement of a rectangular flat plate by an impinging jet in cross flow,Chin.J.Chem.Eng.22(5)(2014)489–495.

[36]M.Hecht,J.Harting,Implementation of on-site velocity boundary conditions for D3Q19 lattice Boltzmann simulations,J.Stat.Mech.Theor.Exp.1(1)(2008)01018.

[37]S.Hou,Q.Zou,S.Chen,G.Doolen,A.Cogley,Simulation of cavity flow by the lattice Boltzmann method,J.Comput.phys.118(2)(1995)329–347.

[38]Q.Zou,X.He,On pressure and velocity boundary conditions for the lattice Boltzmann BGK model,Phys.Fluids 9(6)(1997)1591–1598.

[39]R.Mei,W.Shyy,D.Yu,L.S.Luo,Lattice Boltzmann method for 3-D flows with curved boundary,J.Comput.phys.161(2)(2010)680–699.

[40]M.C.Sukop,D.T.Thorne Jr.,Lattice Boltzmann modeling:An introduction for geoscientists and engineers,Spring-Verlag Berlin Heidelberg press,42-44,2006.

[41]White F.M.,Viscous fluid flow(3rd Ed.),Mc Graw Hill press,112-113(2006).