NOVEL IMMERSED BOUNDARY-LATTICE BOLTZMANN METHOD BASED ON FEEDBACK LAW

Li Xiujuan,Zhao Rongguo,Zhong Chen gwen

(1.College of Civil Engineering and Mechanics,Xiang tan University,Xiangtan,411105,P.R.China;2.National Key Laboratory of Science and Technology on Aerodynamic Design and Research,Northwestern Polytechnical University,Xi’an,710072,P.R.China)

INTRODUCTION

The lattice Boltzmann method(LBM)has been developed to bean alternative and promising numerical scheme for simulating fluid flows and modeling physics in fluids since proposed in 1992.It provides many advantages,including clear physical pictures, easy implementation of boundary conditions, and fully parallel algorithms[1].This method is especially useful for low speed flows, multiphase flows, porous flows, micro channel flows, and particle suspension flows.

The immersed boundary method(IBM)was first proposed by Pesk in in the 1970s when he studied the blood flow in the human heart[2-3].The blood flow is driven by heart valves,and valves immersed in the blood fluid flow can be considered as a deformable boundary.In IBM,the flow field is represented by a set of Eulerian points,which are in fact the fixed Cartesian mesh points,and the boundary of immersed object is represented by a set of Lagrangian points.And this scheme is suitable to simulate fluid structure interaction. The basic idea of IBM is that the effect of boundary on the surrounding fluids can be replaced by the restoring force which is added into governing equations.On the other hand,as an alternative numerical scheme to Navier-Stokes equation solver, LBM has attained wide popularity recently in solving various incompressible flow problems.It is applied to the Cartesian mesh instead of body-fitted mesh which is computationally expensive,so,it can simplify the boundary treatment and improve the computational accuracy.Since the birth of IBM,numerous modifications and a number of variants of this approach have been proposed and widely used in lots of fields,such as bio mechanics,fluid-structure-interaction and Multiphase flow.

LBM is a regular lattice-based scheme,which is a little difficult to deal with the complex geometry.In 2002,IBM is first used in the LBM framework by Feng Zhigang, et al[4], this approach is called immersed boundary-lattice Boltzmann method(IB-LBM).Usually,IBM is coupling with N-S equations.While LBM is a mesoscopic method,its governing equation is the Boltzmann equation. These two methods are completely different. IB-LBM inherits the advantages of LBM:clear physical pictures,easy implementation of boundary conditions,and fully parallel algorithms.

The conventional IB-LBM where the force density is computed explicitly by Hooke’s law and the non-slip boundary condition is approximately satisfied. As a result,obvious flow penetration to the immersed boundary can be observed at the boundary points. Flow penetration implies mass exchange across the boundary.As we know,this can lead into the error.Wu and Shu,et al[5]proposed a scheme called velocity correction to calculate the force term at the boundary points,which guaranteed the non-slip condition in 2009.The capability of that scheme for moving object is well validated in their paper. However, it seems a little complicated to interpolate between two types of coordinates.In this paper,the force density term is obtained by the feedback law.The method not only reduces the computation load for interpolation, but also satisfies the non-slip boundary condition.And,it can be found that the method can be used to simulate incompressible flows around complex geometry object,moving objects and so on.

1 NUMERICAL METHOD FOR IB-LBM

1.1 Lattice Boltzmann method

As show n in Fig.1,the area of viscous flow field is K,which is represented by a set of

Fig.1 Schematics of immersed boundary method

Eulerian points.The elastic boundary isΓ,which is represented by a set of Lagrangian points X(s,t).

In 2000,Hofler and Sch war zer proposed a new IB method to simulate the action of solid particles,which added a body force to the governing equations[6].The governing equations of viscous incompressible flows can be written as

where u and f are the fluid velocity and the force density acting on the fluid,d is the fluid density,p the fluid pressure,and_the dynamic viscosity.

In the lattice Boltzmann context,Eqs.(1,2)can be replaced by the lattice Boltzmann equation[7-8],shown as

where f T is the density distribution function,its corresponding equilibrium state,f the single relaxation time,F the boundary force density of Lagrangian coordinates,F T the force term along the T discrete direction of F,eαthe lattice velocity,and k T the coefficient depending on the selected lattice model.

Eq.(3)is called LBGK modeling.Starting from an initial state,the configuration of particles at each time step evolves in two sequential steps[9]:

(1)Collision,which occurs when particles arriving at a node interact and possibly changes their velocity directions according to scattering rules.It can be written as

(2)Streaming,where each particle moves to the nearest node in the direction of its velocity.It can be written as

The diagram for the particle velocity of D2Q9 is sketched in Fig.2.In every lattice there are sets of particles, every set has different velocity[8].

Fig.2 D2Q9 lattice model

where c=W x/W t,W x and W t are the lattice spacing and time step,respectively.The corresponding equilibrium distribution function is

where csis the sound speed of this model,and it is defined as c s=c/.The weight coefficient k T is written as

The relaxation time is related to the viscosity by

1.2 Velocity correction in IB-LBM

In LBM,the macroscopic variables satisfied N-S equations are obtained by Chapman-Enskog expansion[8],which are defined in terms of the particle distribution function by[4,10]

where u*is the intermediate velocity,and W u the velocity correction.

Then Eq.(4)can be written as

In the conventional IB-LBM,the force density is computed explicitly by the Hook’s law or the direct forcing method[11-14],taking Hook’s law for example,it can be written as

where W(x- X(s,t))is called smoothly approximate function.It can distribute force term f at Eulerian points from the force term F at the boundary(Lagrangian)points.Suppose the force density f is unknown,which is determined by the velocity at the boundary point interpolated from the corrected velocity field satisfying the non-slip boundary condition, otherwise some flows penetrate the boundary.

where u*can be obtained from Eq.(10).Set W Ul

B be the velocity correction vector at every Lagrangian point,so the velocity W u at the Eulerian points can be obtained as same as Eqs.(13,14),shown as

where the Dirac delta function interpolation can be smoothly approximated by continuous kernel distribution,shown as

where W(x-X B(s,t))is proposed by Pesk in[15]in 2002 as

From Eqs.(16-17),the velocity correction at Eulerian points can be expressed as

By parity of reasoning,we have

Substituting Eq.(19)into Eq.(15),we have

輔助資源方面，希望豐富相關(guān)英語視頻和書籍，以豐富學(xué)生的英語輸入類型和數(shù)量，開拓學(xué)生的思維和眼界，激發(fā)學(xué)生學(xué)習(xí)英語的興趣和動(dòng)力。

Substituting Eq.(21)into Eq.(20),we have

Eq.(22)can be also expressed as

where Dij=Dij(x ij-XlB(s,t)),Y=…,B={Δu1,Δu2,… ,Δu l,… ,Δu m}andΔu l=(x ij,t)DijΔxΔy(l=1,2,…,m).

Eq.(23)can be simplified as AY=B.

In this paper,the density and pressure are computed by

The force term f at Eulerian points can be simply calculated as

1.3 Feedback law in IB-LBM

In order to solve two problems caused by velocity correction:the complex interpolation of velocity correction method, and partly interruption of intrinsic parallel nature in solving linear equations of corrected velocity. A new method which introduces feedback law is proposed in this paper.In terms of analysis of mechanics,F denotes the Lagrangian forcing exerted on the boundary by the surrounding fluid. The cylinder used in this paper is in extensible and massive.The governing equation for the cylinder is written in a Lagrangian form.The motion equation is

where X is the position of the cylinder.According to the two hypotheses,Eq.(26)can be presented as

where T is the tension force along the cylinder,V the bending rigidity,d 1 the additional boundary density of cylinder which is different between the structure and the fluid density.

This paper introduces the following characteristic scale:diameter of the cylinder D for the length, free stream velocity U0for the velocity,for the pressure,and d 1for the force term exerted on fluid and boundary,respectively,_U0for the stretching and shearing coefficients, and _U0D2for the bending and twisting coefficients.Thus,Eq.(27)can be written as

And the in extensibility condition[13,16-18]is expressed by

Substituting Eq.(29)into Eq.(28),we can obtain

Note that the tension force of cylinder can be ignored.The tension force T could be omitted from Eq.(30).

The bending force can be expressed as

Bending force F b can be also written as

Set X*denote a pre-position of the cylinder,which is used to substitute Xnin Eq.(32).In order to satisfy the in extensibility condition,the predicted position can be obtained by X*=2Xn-Xn-1. In practice,the use of X*instead of Xnreduces the error further.

Considering Eq.(17),Eq.(33)is similar to Eq.(14),namely

The interaction force can be calculated explicitly by feedback law[19],which can satisfy the non-slip condition.The definition of the law is expressed as

The velocity of the boundary can be computed by Uij,velocity of the points U B nearby are approximate to Uji. As the cylinder is stationary,

The basic solution process of the method can be summarized as follows:

(1)Set initial values.

(2) Use Eq.(3) to get the density distribution function at time level(setting initially),and calculate the macroscopic variable,such as density,pressure,using Eqs.(4,24).

(3)Obtain the force density using Eqs.(32,34),including bending force and Lagrangian force term.

(4)Interpolate the Lagrangian force term into Eulerian force which is exerted on fluid by the boundary.

(5)Compute the equilibrium distribution function using Eq.(7).

(6)Repeat(2)-(5)until convergence is reached.

2 RESULTSAND DISCUSSIONS

2.1 Steady flow around stationary cylinder

The capability of the method for solving fluid-solid boundary problem is well demonstrated through its application to simulate flows around a single circular stationary cylinder. The simulations at Reynolds numbers of 20 and 40 are carried out.The data in the test are shown as follow:

(1)The computational domain is 500×241,and the fine-mesh block covers the area of cylinder and its wake.

(2)The diameter of the cylinder is D=20 and its centre is located at(125,120).

(3)Initial states are U0=0.1,V0=0.0,d0=1.0.

(4)Setting of LBM isΔx=Δy=Δt=1.

Here,Reynolds number and drag coefficient are defined as

where U 0 is the free stream velocity,F d the drag force and it can be calculated by

where f x is the x-component of the force density at the boundary points.

The streamlines of Re=20,40 are presented in Fig.3. From Fig.3,it is clear that the separation flow is stable and symmetrical,and there area coupleof symmetrical vortexes behind thecylin-der.The results show that the length of the vortex is elongated with Reynolds number.

Fig.3 Streamlines around cylinder

The drag coefficient and length of recirculation bubbles(based on the cylinder diameter D)are compared with those of previous literatures in Table 1.

Table 1 Comparison of flow around circular cylinder

The good agreement of the result can be further confirmed by Fig.4,which compares the pressure profile on the surface of cylinder at Re=40.

Fig.4 Pressure distribution on surface of cylinder at Re=40

First of all,LBM is a mesoscopic numerical method. The number of meshes is related to Reynolds number. Then, the number of boundary points has effect on the results,shown in Table 2.Finally,the diameter of the cylinder is related to the accuracy of the results. To demonstrate the effect,numerical simulations for diameter with different mesh sizes are carried out.The results are shown in Table 3.It is shown that the distance of Lagrangian points can affect the accu-racy.If the relationship between arclength and space step satisfies the equationΔs ＜ 0.5Δx,the results will be more accurate.

Table 2 Influence of number of Lagrangian points on drag coefficient at Re=20

Table 3 Influence of cylinder diameter on drag coef ficient at Re=40

2.2 Flow around two circular cylinders in side by side arrangement

For the case presented above,the stationary circular cylinder is single.In order to investigate the ability of the method for complex geometry flows, the simulation of flow around two stationary circular cylinders in a side by side arrangement is carried out at Re=40.

The computational domain is set by 40D×22D with a mesh size of 801×441.The two circular cylinders are set normal to the free stream,and located with the centers at(201,201)and(201,241),respectively.That is to say,the non-dimensional gap spacing g*=g/D=1.The free stream velocity is taken as U0=0.1,V0=0.0 and the fluid density is d 0=1.0.The computation starts with the given free stream velocity.The density distribution function is set in its equilibrium state at the far field boundary and the in flow.

The streamlines of Re=40,g* =1 are presented in Fig.5.From Fig.5,it is shown that the flow is steady and symmetric relative to the centerline,and the reis no vortex shedding behind the cylinders[22].

Fig.5 Streamlines for two side by side cylinders at Re=40,g*=1

3 CONCLUSION

This paper proposes a novel IB-LBM method based on the feedback law. The method is sufficient to capture the important characteristics of flow fields,and the results are within the range of values reported by previous studies.Compared with the conventional IB-LBM or the velocity correction method[23-26],the force density is obtained by feedback law in the method,which not only satisfies the non-slip boundary condition,but also improves the accuracy of simulation for complex geometry by LBM.Furthermore,the method can be easily computed and reserve the advantages of LBM.

[1] Chen S Y,Doolen G D.Lattice Boltzmann method for fluid flows[J]. Annual Review of Fluid Mechanics,1998,30(1):329-364.

[2] Peskin C S.Flow patterns around heart valves:A numerical method[J]. Journal of Computational Physics,1972,10(2):252-271.

[3] Peskin C S.Numerical analysis of blood flow in the heart[J].Journal of Computational Physics,1977,25(3):220-252.

[4] Feng Z G,Michaelides E E.Hydrodynamic force on spheres in cylindrical and prismatic enclosures[J].Journal of Multiphase Flow,2002,28(3):479-496.

[5] Wu J, Shu C, Zhang Y H. Simulation of incompressible viscous flows around moving objects by a variant of immersed boundary-lattice Boltzmann method[J].International Journal for Numerical Method in Fluids,2010,62(3):327-354.

[6] Ho¨fler K,Schwarzer S.Navier-Stokes simulation with constraint forces: Finite-difference method for particle-laden flows and complex geometries[J].Physical Review E,2000,61(6):7146-7160.

[7] Guo Z,Zheng C,Shi B.Discrete lattice effects on the forcing term in the lattice Boltzmann method[J].Physical Review E,2002,65(4):046308.

[8] McNamara G R,Zanetti G.Use of the Boltzmann equation to simulate lattice gas automata[J].Physical Review Letters,1998,61(20):2332-2335.

[9] Qian Y H,d’Humieres D,Lallemand P.Lattice BGK model for Navier-Stokes equation [J].Europhysics Letters,1992,17(6):479-484.

[10]Shu C,Liu N Y,Chew Y T.A novel immersed boundary velocity correction—Lattice Boltzmann method and application to simulate flow past a circular cylinder[J]. Journal of Computational Physics,2007,226(2):1607-1622.

[11]Feng Z G,Michaelides E E.Proteus: A direct forcing method in the simulations of particulateflows[J].Journal of Computational Physics,2005,202(1):20-51.

[12]Fadlun E A, Verzicco R, Orlandi P, et al.Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations[J].Journal of Computational Physics,2000,161(1):35-60.

[13]Niu X D,Shu C,Chew Y T,et al.A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows[J].Physics Letters A,2006,354(3):173-182.

[14]Peskin C S.The immersed boundary method[J].Acta Numerica,2002,11:479-517.

[15]Tornberg A-K, Shelley M J. Simulating the dynamics and interactions of flexible fibers in stokes flows[J].Journal of Computational Physics,2004,196(1):8-40.

[16]Thess A,Zikanov O,Nepomnyashchy A.Finitetime singularity in the vortex dynamics of a string[J].Physical Review E,1999,59(3):3637-3640.

[17]Belmonte A,Shelley M J, Eldakar S T,et al.Dynamic patterns and self-Knotting of a driven hanging chain[J].Physical Review Letters,2001,87(11):114301.

[18]Huang W X,Sung H J. Simulation of flexible filaments in a uniform flow by the immersed boundary method[J]. Journal of Computational Physics,2007,226(2):2206-2228.

[19]Goldstein D,Handler R,Sirovich L.Modeling a noslip flow boundary with an external force field[J].Journal of Computational Physics,1993,105(2):354-366.

[20]Shukla R K,Tatineni M,Zhong X.V ery high-order compact finite difference schemes on non-uniform grids for incompressible Navier-Stokes equations[J].Journal of Computational Physics,2007,224(2):1064-1094.

[21]Nieuw stadt F, Keller H B.Viscous flow past circular cylinders[J].Computers&Fluids,1973,1(1):59-71.

[22]Kang S.Characteristics of flow over two circular cylinders in a side-by-side arrangement at low Reynolds numbers[J].Physics of Fluids,2003,15(9):2486-2498.

[23]Wu J, Shu C, Zhang Y H. Simulation of incompressible viscous flows around moving objects by avari ant of immersed boundary-lattice Boltzmann method[J]. International Journal of Numerical Methods for Heat and Fluid Flow,2010,62(3):327-354.

[24]Wu J,Shu C.An improved immersed boundary lattice Boltzmann method for simulating three dimensional incompressible flows[J]. Journal of Computational Physics,2010,229(13):5022-5042.

[25]Choi J I,Oberoi R C,Edwards J R,et al.An immersed boundary method for complex incompressible flows[J].Journal of Computational Physics,2007,224(2):757-784.

[26]Lu J,Han H,Shi B,et al.Immersed boundary lattice Boltzmann model based on multiple relaxation times[J].Physical Review E,2012,85(1):016711.