Lattice Boltzmann method for fluid flow around bodies using volume penalization

This paper deals with the implementation of a volume penalization technique in a lattice Boltzmann model, in order to compute flows around obstacles. The penalization term was introduced into the lattice Boltzmann equation via a forcing term. This approach was applied to the one dimensional Burgers equation for motionless and moving obstacles (forced motion, and coupling between the fluid force calculated with the penalized Burgers equation and the motion of the obstacle), and to the two dimensional Navier-Stokes equation, for motionless obstacles (flows over a square obstacle, and past a circular cylinder). A good agreement with numerical results obtained with other techniques, and with results found in literature was obtained.


INTRODUCTION
Flows around moving obstacles are encountered in many applications such as particle laden flows, flows around ships, sailing boats, wind turbines, aircrafts, and so on.Due to the increasing performance of computers, numerical simulation becomes an attractive tool for predicting such flows.The lattice Boltzmann method (LBM), which is easy to implement, involves local computations, and is thus well suited for parallel computing, is a good candidate for solving the equations that describe flows around moving objects.Various techniques have been developed in the past to model particles moving in fluids, in the lattice Boltzmann framework.
In his pioneering work, Ladd [1,2], used the bounce back rule to model the no-slip boundary condition between the particle and fluid, and developed the momentum exchange method to calculate the fluid force exerted on the solid particle.
Noble and Torczynski [3], proposed to modify the collision term in the lattice Boltzmann equation, to account for the interaction between a solid obstacle and the lattice nodes.With this method, the conventional lattice Boltzmann equation was recovered for fluid regions without obstacle, and rigid body motion (i.e. the Zou and He bounce back rule [4]) was obtained in regions occupied only by the solid obstacle.This method was applied successfully to the computation of flow around a periodic line of oscillatory moving cylinders.
The immersed boundary method, which was introduced in the 1970s by Peskin [5] to model blood flow in the heart, and which employs a fixed cartesian mesh for the fluid and lagrangian nodes attached to the moving boundary, and adds a forcing term in the Navier-Stokes equations to model the influence of the boundary on the fluid flow, was adapted to the LBM for flows around rigid particles by Feng and Michaelides [6], and by Ten Cate et al. [7].The forcing term was modeled by Feng and Michaelides with a restauration force acting on the particle boundary.Ten Cate et al. used for this term a variant developed by Goldstein et al. [8], where the restauration force on the boundary can be due, from a physical point of view, to a damped oscillator.The drawback of this method was that results depended on one or two empirical parameters.To avoid this problem, Feng and Michaelides [9], and Dupuis et al. [10], introduced into the IB-LBM (Immersed Boundary -lattice Boltzmann method), the direct forcing approach proposed by Fadlun et al. [11].In this approach, the forcing point was located in the fluid region in the immediate vicinity of the solid boundary, and an interpolation between the solid points and the fluid points was used to evaluate the forcing term.With this approach, the no-slip conditon is not exactly satisfied at the solid boundary, and some streamlines may penetrate in the solid obstacle.To enforce the no-slip condition at the solid boundary, Wu and Shu [12], by using the formulation of the forcing term proposed by Guo et al. [13], developed an implicit procedure, where a velocity correction was calculated for all lagrangian points at the boundary.Owing to the fact that the fluid domain is not remeshed when applying the immersed boundary method, this technique is now frequently used in the lattice Boltzmann method.However, other approaches have been implemented in the lattice Boltzmann framework.
Among these, we can cite the Distributed Lagrange Multiplier/Fictitious Domain approach (DLM/FD) [14].Shi and Phan-Thien [15] applied to the lattice Boltzmann method the DLM/FD formulation developed by Yu [16].In this approach, the fluid equation was solved on a fictitious domain including the fluid domain and the solid domain, and a force (the Lagrange multiplier) was introduced to force the velocity of the fictitious fluid inside the solid domain to be the same as the solid velocity.The solid domain was followed with a Lagrangian description.With this method, the remeshing procedure, which is computationally time consuming, was avoided.
Very recently, Meldi et al. [17] applied the Arbitrary Lagrangian Eulerian formulation to the LBM.In the Arbitrary Lagrangian Eulerian description, the grid is divided into two parts: a fixed one, and a moving one which follows for example the flow around a moving obstacle.The nodes of the moving grid are displaced arbitrarily using a Lagrangian description.In this region, the Navier-Stokes equations are solved on the moving grid [18].Using the Chapman-Enskog procedure, Meldi et al. developed the lattice Boltzmann formulation for a moving grid.Non-deforming grids were chosen.They applied this technique successfully to the computation of flows around non-deforming obstacles.
In this work, we wish to implement the volume penalization technique, to compute flows around obstacles, with the lattice Boltzmann method.The volume penalization technique consists in introducing a volume penalization term into the equation that needs to be solved, in order to take into account the influence of the obstacle on the fluid domain [19,20].Since this equation is solved on both fluid and solid domains, the boundary conditions at the fluidsolid interface are naturally applied.Hence this method appears to be easy to implement, and to parallelize in a lattice Boltzmann framework.In this paper, we applied, in a first step, the volume penalization technique to the one dimensional Burgers equation, which was solved with the lattice Boltzmann method.Several cases were investigated: non moving obstacle, imposed motion of the obstacle, coupling between Burgers equation and the motion of the obstacle.In a second step, the volume penalization approach was implemented in the Navier-Stokes equations, which were solved with the lattice Boltzmann method.Two cases were considered: flow past a square obstacle, and flow past a circular obstacle.In the last section of the paper, conclusions and future work are pointed out.

ONE DIMENSIONAL BURGERS EQUATION
Let us consider a domain Ω, composed of a fluid domain Ω f (t) (wheren t denotes time), and a solid domain Ω s (t), where the fluid solid interface X s (t) may move according to time (see Figure 1).
The volume penalization technique was applied in this paragraph to the one dimensional Burgers equation.The one dimensional Burgers equation, with a volume penalization term, reads: (1 Where u is the velocity according to x direction, v is the kinematic viscosity of the fluid.In the penalization term, χ is a mask function, which is equal to 0 in the fluid domain Ω f (t), and to 1 in the solid domain Ω s (t), η is a permeability coefficient that is very high in the fluid domain and very small in the solid domain, and u s is the obstacle velocity.In the fluid domain, the ratio χ/η is null, hence the penalization term vanishes.In the solid domain, the ratio χ/η is very high, the velocity is thus equal to the solid velocity.
In order to solve this equation with the lattice Boltzmann method, the lattice Boltzmann model for a non linear convection diffusion equation developed by Shi and Guo [21] was selected.The lattice Boltzmann equation, with the single relaxation time Bhatnagar-Gross-Krook (BGK) collision operator, is: Where is the probability of finding a fluid particle with velocity at position and time t, Δt is the time step, τ is the non dimensional relaxation time, is the equilibrium distribution function, F α is the forcing term that takes into account, in our study, the penalization term.For this one dimensional case, the D1Q3 model was used.In this model, the lattice velocities are , where c = Δx/Δt, and x is the lattice spacing, the values of the weight coefficients are ω 0 = 2/3, ω 1 = ω 2 = 1/6, the sound speed is .Shi and Guo showed that the equilibrium distribution function is: (3) C c uI e e c I c , .: and the relationship between the diffusion coefficient in Burgers equation, and the non dimensional relaxation time is: Shi and Guo obtained the following expression for the forcing term: ( With the Chapman-Enskog analysis, they showed that λ is obtained from: When applying the penalization technique, instability occurs if the penalization term is treated explicitly.To circumvent this difficulty, an implicit treatment of the penalization term was adopted here.Using the Chapman-Enskog procedure, this led to the following expression for λ in the forcing term: To test the suitability of the volume penalization technique to handle a moving obstacle with the lattice Boltzmann approach, the one dimensional Burgers equation was solved for three cases: non moving obstacle, imposed motion of the obstacle, coupling between Burgers equation and the motion of the obstacle. For each case, the LBM solution written in dimensional units was compared with a dimensional solution (an analytical solution or a numerical solution computed with the finite difference method).
For this case, we compared the exact solution in the fluid domain: (6) Where: with the numerical solution computed with the penalized Burgers equation and the lattice Boltzmann method.Figure 2 shows the exact solution, and the numerical one, at different times, for v = 0.05 m 2 /s, Δx = 0.003333 m (i.e.301 lattice nodes), and Δt = 1.85 × 10 −5 s.The penalization coefficient was η = 10 −9 .A satisfactory agreement can be noticed.To compare the results given by the lattice Boltzmann method applied to the penalized Burgers equation, we also solved the penalized Burgers equation with the finite difference method (where an implicit scheme was employed for the temporal integration, and 301 nodes were also used for the finite difference computations).In Figure 3, where the solution at different times is presented, we can see that the interface moves at a constant velocity according to time.We can also notice a good agreement between the numerical results computed with the penalized Burgers equation and the lattice Boltzmann method, and those obtained with the finite difference method.

COUPLING BETWEEN THE FLUID FORCE CALCULATED WITH THE PENALIZED BURGERS EQUATION AND THE MOTION OF A SOLID OBSTACLE
In this paragraph, we considered that Burgers equation was coupled with the motion of a spring, governed by the following equation: where x s (t) is the position of the fluid solid interface, m and k are the solid mass and spring constant respectively, X 0 is the spring natural length.F(x s (t)) is the force exerted by the fluid on the obstacle, it was calculated with: The fluid velocity was calculated with the penalized Burgers equation in the lattice Boltzmann framework, in the one dimensional domain Ω = [0; 1.2].The initial condition was: u(x, 0) = sin(πx), x ∈ Ω f (0) = [0; 1], and u(x, 0) = 0, x ∈ Ω s (0) = [1; 1.2].The boundary The results obtained with the lattice Boltzmann method are compared with those computed with the finite difference method (with 301 nodes), where the volume penalization method was also applied.In Figure 4, where the position of the interface, according to time, is shown, we can see that, for the case k = 39.5 N/m, the interface displacement was governed by the spring, but it was damped by the fluid forces.For the case k = 0.395 N/m, the force exerted by the spring was small, and the interface motion was mainly due to the fluid forces.In Figure 5, the dimensional velocity profile obtained at time t = 1 s is also presented.In Figures 4 and 5, we can notice that a good adequation was found between the results computed with the volume penalization method applied to the lattice Boltzmann equation, and those obtained with the finite difference method.

TWO DIMENSIONAL NAVIER-STOKES EQUATIONS
Next, we applied the volume penalization technique to the Navier-Stokes equations, in the lattice Boltzmann framework.The Navier-Stokes equations where the volume penalization technique was applied, read [19]: where μ and ρ are the dynamic viscosity and the density of the fluid considered, is the fluid velocity, and p is the pressure.In this paper where the first step of our work is presented, the Navier-Stokes equations were solved for motionless obstacles The continuity equation was satisfied: (10) The lattice Boltzmann method ( [22], [23], [24]), which can be easily implemented and is naturally suited for parallel computing, was chosen for solving the flow equations.The lattice Boltzmann equation, with the BGK collision operator, is given in equation 2. Two dimensional flows were computed, and the D2Q9 lattice model was selected.In the D2Q9 model, the lattice velocities are for α = 1, 2, 3, 4, for α = 5, 6, 7, 8, where c = Δx/Δt.The equilibrium distribution function is: (11) where ω a are the weighting coefficients (for the D2Q9 model: ω 0 = 4/9, ω α = 1/9 for α = 1, 2, 3, 4, ω α = 1/36 for α = 5, 6, 7, 8, c s is the speed of sound for the D2Q9 model).In the lattice Boltzmann equation (equation 2), the penalization term was modelled using the forcing term proposed by Guo et al. ( [13]): (12) where the fluid velocity is defined according to: uu e e c I c , 1 : As shown by Guo et al., this model avoids spurious spatial and temporal terms in the continuity and momentum equations.In equation 13, the left hand side and the right hand side are taken at the same time instant.Since the right hand side also contains the fluid velocity (via the penalization term), the fluid velocity was calculated according to the following expression (deduced from equation 13): (14) This procedure was applied by Guo and Zhao [25] who built a lattice Boltzmann model in order to compute incompressible flows in porous media, but they did not use the volume penalization technique and treat flows around moving obstacles.
In order to validate the present methodology, two cases were considered: flow past a square obstacle in a channel, and flow past a circular cylinder.For these cases, quantitative results were obtained.The drag force F d and lift force F l were computed with the momentum exchange method [1], and the drag and lift coefficients C d and C l were deduced: (15) (where D is a characteristic dimension of the obstacle, and u max is the maximum velocity in a region where the flow is not influenced by the obstacle).The Strouhal number (where f is the frequency of the vortex shedding) was also computed.

FLOW PAST A SQUARE OBSTACLE
This configuration is depicted in Figure 6.A parabolic velocity profile was applied at the inlet, a zero velocity gradient was applied at the outlet, and no-slip velocity conditions were imposed at the walls.The blockage ratio of this configuration was B = H/D = 8, and the length of the channel was L/D = 50.In the lattice Boltzmann simulations, the Zou and He ( [4]) boundary conditions (i.e.bounce back of the non equilibrium populations) were used.Our results were compared with those obtained by Breuer et al. ( [26]) who performed lattice Boltzmann simulations and employed the simple bounce back boundary conditions to impose the no-slip boundary conditions at the walls (walls and square obstacle).Our simulations were conducted with regular grids, and grid refinement tests were done.In this paper, we show the results obtained for ( where u max is the maximum velocity at the inlet), with the finest grid of 3000 × 480 cells.For this configuration, we carried out simulations using on the one hand the penalization method (with a penalization coefficient of 10 −9 ), and on the other hand the Zou and He bounce back boundary condition at the square obstacle.In Figure 7, we plotted the streamlines computed with the penalization technique, for Re = 40, and Re = 100.For low Reynolds numbers, the flow is steady, and two counterrotating vortices appear symmetrically about the flow axis behind the square obstacle (see In Figures 8 and 9, we compare the profiles of streamwise velocity along the flow axis, and the profiles of streamwise velocity in planes perpendicular to the flow axis, computed with the penalization technique, and with the Zou and He bounce back boundary condition, for the case Re = 100, at the same time instant (i.e. at 666,67s).We can see that these two techniques yielded similar results.
The non dimensional recirculation length, the drag and lift coefficients, and the Strouhal number were gathered in Table 1, as well as Breuer et al.'s results.We can notice a satisfactory agreement between the results computed with the penalization technique, and those computed with the Zou and He bounce back condition.There is a small discrepancy between our results and Breuer et al.'s results that were obtained on a coarser grid (2000 × 320 cells).

FLOW PAST A CIRCULAR CYLINDER
The second configuration is shown in Figure 10.The boundary conditions at the horizontal upper and lower planes were symmetry boundary conditions, and a flat velocity profile was applied at the inlet.The boundary condition at the outlet was the same as in the previous case.Three cases were considered: Re = 20, 40 and 100, where the Reynolds number was built using the cylinder diameter, and the constant velocity at the inlet.The domain dimensions were: L1 = 15 D, L2 = 30 D, H = 60 D. Due to the curved geometry of the cylinder, this configuration was more complex than the previous one.A large number of grids were  necessary to compute this flow.A grid independence test was performed for this case; the results presented in this paragraph were computed with the finest grid (1845 × 2460).For this configuration, the flow around the cylinder was computed with the penalization method on the one hand (with a penalization coefficient of 10 −9 ), and with the bounce back boundary condition proposed by Bouzidi et al. [27] for curved boundaries on the other hand.
The streamlines around the cylinder, computed with the penalization method, are shown in Figure 11.For low Reynolds numbers, symmetric recirculation regions behind the cylinder can be noticed.For higher Reynolds numbers, vortex shedding occurs.With the penalization method, we found that the critical Reynolds number below which the flow is steady was 45.This is in agreement with the critical value of 45 ∼ 49 given in literature.
In Figures 12 and 13 we plotted the profiles of streamwise velocity along the flow axis, and in planes perpendicular to the flow axis, obtained at the same time instant (i.e. at t = 428,32s).A good agreement between the results computed with the penalized Navier-Stokes equations solved with the lattice Boltzmann method, and those obtained with the bounce back boundary condition applied on the fluid/solid interface of the cylinder, can be highlighted.The non dimensional recirculation length, the drag and lift coefficients, and the Strouhal number that we computed with the penalization approach, and the bounce back method, are reported in Table 2.In this table, the simulation results of Zhou et al. [28], He and Doolen [29], Wu and Shu [12], and the experimental results of Tritton [30], and Williamson [31] are included for comparison.The simulation results, found in literature, included in Table 2, were obtained by the lattice Boltzmann method, and the non-equilibrium extrapolation method developed by Guo et al [32] (cf Zhou et al.), an interpolation approach applied to a curvilinear grid (cf He and Doolen), an immersed boundary technique (cf Wu and Shu).Although the results computed with the penalization approach were in agreement with the results in literature, we noticed a slightly higher discrepancy for this more complex geometry than for the square obstacle, which may be due to the high number of grids that is necessary when the penalization technique is applied.

CONCLUDING REMARKS AND FUTURE WORK
In this work, a volume penalization technique was employed to simulate flow past obstacles, using the lattice Boltzmann approach.This technique was applied in a first step to the one dimensional Burgers equation.We compared our results for cases of a motionless obstacle, of the forced motion of an obstacle, and of the coupling between the fluid force calculated with the penalized Burgers equation, and the motion of a solid obstacle, with analytical results, and results obtained with the finite difference method, and we noticed a satisfactory accuracy of these results.The second step of our study focused on the numerical prediction of incompressible flows around motionless obstacles (the Navier-Stokes equations).A satisfactory agreement with results found in literature (lattice Boltzmann results and experimental results) was observed.Regular grids were used, and many cells were necessary to obtain grid independent results.A local grid refinement technique is now introduced into our lattice Boltzmann code, in order to decrease the number of cells, and thus the computing time.This will then enable us to treat more easily flows around moving obstacles.

Figure 2 :
Figure 2: Exact solution (−) and numerical solution (+) (Lattice Boltzmann simulation of the penalized Burgers equation) for a motionless obstacle, at different times

Figure 3 :
Figure 3: Numerical solutions of Burgers equation, for a moving obstacle with a prescribed motion, at different times.Symbols: -Lattice Boltzmann simulation of the penalized Burgers equation, + solution obtained with a finite difference computation of the penalized Burgers equation

Figure 4 :Figure 5 :
Figure 4: Position of the interface, according to time, for the coupling of Burgers equation with the motion of a spring.Symbols: -Lattice Boltzmann simulation of the penalized Burgers equation, + solution obtained with a finite difference computation of the penalized Burgers equation.a) k = 39.5 N/m, b) k = 0.395 N/m

Figure 10 :
Figure 10: Configuration studied (flow past a circular cylinder)

Figure 11 :
Figure 11: Streamlines superposed on the velocity magnitude isocontours, flow past a circular cylinder computed with the penalization method: a) Re = 40, b) Re = 100

Table 1 :
Comparison of recirculation length, drag and lift coefficients, Strouhal number for flow past a square obstacle at Re = 20, 40, 100 (C d − av is the timeaveraged drag coefficient, (Cl max − Cl min ) av /2 is the average of the difference between the maximum and minimum values of the lift coefficient, divided by 2)

Table 2 :
Comparison of recirculation length, drag and lift coefficients, Strouhal number for flow past a circular cylinder at Re = 20, 40, 100 (C d − av is the timeaveraged drag coefficient, (Cl max − Cl min ) av /2 is the average of the difference between the maximum and minimum values of the lift coefficient, divided by 2)