An improved immersed-boundary algorithm for fluid-solid interaction in attice-boltzmann simulations

G. Boroni^†, J. Dottori^†, D. Dalponte^‡, P. Rinaldi^† and A. Clausse^†‡

† CONICET and Universidad Nacional del Centro, Tandil, Argentina.
‡ Comisión Nacional de Energía Atómica, Argentina.
gboroni@conicet.gov.ar, jdottori@exa.unicen.edu.ar, diego.dalponte@gmail.com, rinaldipablo@gmail.com and clausse@exa.unicen.edu.ar.

Abstract— An improved algorithm combining the features of the lattice-Boltzmann and the immersed-boundary methods is presented. Following previous formulations, the method represents a fluid constrained by flexible boundaries by means of a force term acting on the cells adjacent to the boundary, which in turn is moved by the fluid. The present algorithm introduces a more efficient iteration procedure to calculate the fluid-boundary interaction, which facilitates the implementation and improves performance. The simulations were validated against experimental and analytical data showing good agreement and demonstrating the performance of the method to simulate different kind of fluid-solid interaction.

Keywords— Immersed Boundary Method; Lattice Boltzmann Method; Fluid-structure Interaction.

INTRODUCTION

The lattice-Boltzmann method (LBM) is an evolution of the lattice-gas automata, a discrete particle kinetics method originated from Boltzmann's kinetic theory of gases. In the method, time and space are discretized to solve the Boltzmann equation for particle velocity distribution functions on a regular grid. The LBM is an alternative to traditional numerical methods to obtain the flow solution. The LBM is rather simpler and more suitable to add mesoscopic kernels to simulate new phenomena as compared to conventional approaches that solve the flow problem directly discretizing the partial differential equations of the flow.

In the past few decades LBM has achieved great success in simulating various incompressible flow problems with simple bounce-back rules to represent solid walls. However, it is not a straightforward task to simulate complex geometries over a regular grid with simple streaming-collision rules.

On the other hand, the simulation of flows in flexible channels is very important in some applications, like hemodynamics. Among various numerical solvers, the development of Cartesian mesh solvers is an active research area due to their simplicity to simulate flows interacting with moving boundaries (Fadlun et al., 2000; Ge and Sotiropoulos, 2007; Gilmanov and Sotiropoulos, 2005; Johansen and Colella, 1998; Mittal and Iaccarino, 2005; Peskin, 1977; Udaykumar et al., 2001). These problems involve complicated interaction between a viscous fluid, deformable body, and free moving boundary, making them difficult to discern. Consequently, analytic solutions are few and almost nonexistent, but a computational approach is a viable possibility. Nowadays there exist a variety of methods developed to handle this kind of problems. Traditionally solving a fluid interacting with complex geometry relies on the unstructured grid method or the overset grid method (Zhang, 2002 and 2008).

For solid wall representation in LBM, the classic idea consists in treating the lattice cells as fluid or solid. When the bounce-back rule is applied to the solid cells, the incoming distribution particles are returned back to the fluid. Some subsequent works present evolutions of this basic concept, Zou and He (1997) generalized a second order bounce-back scheme to simulate rigid straight walls. Other approaches, like Guo et al. (2007), use extrapolation techniques to implement rigid curved walls in LBM.

Recently, Cheng and Zhang (2010) successfully coupled LBM with the Immersed Boundary method (IB) to simulate moving curved walls in a mitral heart valve. IB (Peskin, 2002) was originally proposed to solve the fluid-flexible-structure interaction problem. This method decouples the solutions of governing equations and boundary conditions. The boundary effect into the surrounding fluid is replaced by adding forces to the fluid equations.

The solution proposed by Cheng and Zhang (2010) uses a regular grid contained in the flow field. The fluid-particles velocity field is solved by adding a force density term into the LBM equation. This method preserves the advantages of LBM in tracking particles and offers an alternative to handling the solid-fluid boundary conditions. The method also solves the problems of forces and velocities fluctuations on the particles when boundary conditions are applied.

In the present article, the Cheng and Zhang's algorithm has been redesigned with a more efficient iteration scheme for calculation of forces at the boundary points. Furthermore, different object-fluid interactions were tested including straight and curved immersed walls and flexible curved walls. The simulations were validated against experimental or analytical data showing good agreement and demonstrating the performance of the method to efficiently simulate different kind of fluid-solid interaction.

II. THE LATTICE-BOLTZMANN METHOD

The Lattice-Boltzmann Method (LBM) is a class of algorithms that produce fast running simulations of fluid flows based in a mesoscopic kinetic representation (Sukop and Thorne, 2006). The method operates on a regular grid representing the fluid by a set of particles that populates and moves between the cells of the grid. The most common model for 2D simulations is D2Q9 (Sukop and Thorne, 2006), which uses a square lattice with 9 velocity directions (Fig. 1). The state of every cell is given by a particle density function f_α(x,t) (not directly observable), representing the number of particles in the node x at time t moving with velocity e_α, which is treated as an internal variable. The function f_α(x,t) changes its value according to predetermined rules simulating the mechanisms of transport, collision and sources of the particles. The physical observables are macroscopic variables generated by moments of f_α(x,t) respect to the internal variable e_a, namely

(1)

(2)

Figure 1. Internal discrete velocities of D2Q9.

where ρ is number of particles and u is average velocity.

For problems with external forces, the discrete LBM equation with is (Chen and Doolen, 1998):

(3)

where (x,t) is the equilibrium distribution function, g_α(x,t) is a source term, δ_t is the time increment, and τ is a relaxation time that is related with the viscosity ν by:

(4)

Since the viscosity is always positive, τ should be greater than 1/2.

The equilibrium density function (x,t) plays the role of local velocity redistribution of the fluid particles due to collisions, and is calculated in terms of the macroscopic variables. This function is the key to LBM dynamics: it determines the partial differential equations that the macroscopic quantities will follow. To approximate the Navier-Stokes (N-S) equations the equilibrium distribution function corresponding to the Bhatnagar-Gross-Krook (BGK) collision operator for D2Q9 is defined by (Chen and Doolen, 1998):

(5)

where

(6)

III. THE IMMERSE BOUNDARY METHOD

The Immerse Boundary (IB) Method assumes that the boundary consists of massless particles, so that the force generated by distortions of the boundary can be calculated and transferred to the fluid (Cheng and Zhang, 2010). Figure 2 shows a 2D example with a closed immersed boundary. The boundary and the fluid domain are denoted by _b and Ω_f, respectively.

Figure 2. Closed immerse boundary in a 2D lattice.

The state of the boundary is represented by X(s,t), a Lagrangian vector function of arc length s and time t, which returns the location of the boundary nodes on _b. The boundary influence is represented by a force density F(s,t) at the boundary point X(s,t). Thus, F(s,t) is determined by the configuration of X(s,t) and it is transferred into the force term g in Eq. (3), which in turn determines the flow velocity and pressure throughout domain Ω_f.

For a boundary immersed in a viscous fluid, the IB is given by the following set of equations (Cheng and Zhang, 2010):

	(7)
	(8)
	(9)
	(10)
	(11)

where u is the flow velocity, ρ the fluid density, p the pressure, µ the fluid viscosity, f_b the external force, X the boundary coordinate, s the boundary fiber length, U the boundary speed, x the fluid flow coordinate, S_f the boundary force generation operator, and δ(r) the Dirac delta function. Equations (7) and (8) are the N-S equations with external force f_b, while Eqs. (9) and (10) are the IB equations. Equation (11) and the right part of Eq. (9) represent the interaction of boundary and fluid.

IV. DISCRETE INTERACTION

The interaction between fluid nodes and boundary points is ruled by integration of the continuous Eqs. (9) and (11). The discretized versions of Eqs. (9) and (11) using a regularized discrete delta function δ_h are

	(12)
	(13)

where h=Δx=Δy is the fluid node spacing and Δs_k is the boundary segment length. The delta function δ_h is an approximation to the Dirac delta function δ(r) (Peskin, 2002), that is:

	(14)
	(15)

The boundary force density F is defined by the boundary configuration. For b with tension, bending and fastening forces, F is expressed as

(16)

where k_c is the tension stiffness, k_γ the bending rigidity, k_f the fastening stiffness and Z the target position of the boundary. The discretized version of Eq. (16) can be expressed as

(17)

V. IB-LBM COUPLING

The external forcing term g_α (Eq. 3) has first order of convergence (Feng and Michaelides, 2004). For problems with slow moving boundary or flexible boundary with small pressure gradient, this term does not affect the result. However, a higher order method is needed for fast moving boundaries or flexible boundaries with large pressure gradients. Cheng and Doolen (1998) proposed an alternative second-order convergence scheme given by:

(18)

where

(19)

Equations (18) and (19) has a second order accuracy for spatial resolution of fluid, matching the order of the original LBM model. The price paid is that Eq. (18) is implicit. In the following section, an efficient iterative procedure to solve Eq (18) is proposed.

VI. ALGORITHM OF IB-LBM METHOD

The IB-LBM coupling procedure needs an iteration process because the scheme of the IB method is implicit (Eq. 18). The proposed solution is to combine the LBM method to solve the fluid equations and the force formula calculated with the IB method. Basically, it uses the external force given by Eq. (12) to propagate the boundary force to the lattice, and Eq. (13) to interpolate the boundary speed in the immediate node lattice. The algorithm, which combines the IB coupling cycle interaction with each step of LBM method, is as follows:

For t=0

1.

(Set target position of the boundary Z_k, the mass density r, and the average velocity.)

For t≠0

2.

(Calculate the equilibrium distribution function.)

3.

(Calculate the collision step)

4.

(Initialize step iteration n = 0 and the boundary point in the t time.)

5. Repeat until convergence and calculate the external forcing term.

	a.
	b.
	(Calculate the force density at the boundary points.)
	c. Initialize
	For each ,
	d.
	(Distribute the force density to the fluid external force.)
	e.
	(Calculate the external forcing terms for velocity vectors.)
	f.
	f. with and
	(Updates distribution functions in the t time for the n iteration.)
	g.
	(Calculate the macroscopic variables.)
	h.

(Update the boundary points positions using the fluid velocity.)

6. Until

(Compare the force density in two iterations to check convergence.)

7. (Updates for ƒ_αij time t.)

8.

(Move to next time.)

It can be seen that in the present algorithm the initialization of forces and the fluid state is explicit, which is important in order to ensure consistency on the initial state of the automata. This was not explained in previous works. Also, the number of cells used here interacting with each IB point is much lower than the total number of cells. In effect, the computation time is improved over Cheng and Zhang (2010) by using only the cells related to each IB point. Another important implementation change proposed here is that is pre-computed (step 3) before the internal iteration (over the same memory space). Also the Dirac delta function can be pre-computed before the internal iterations. These optimizations improve significantly the computation time of the algorithm.

VII. RESULTS

Three examples with IB-LB scheme were simulated to test the performance of the procedure. All the cases were carried out considering a uniform grid.

A. Elliptical deformable contour immersed in a fluid

Let us consider a closed immersed boundary with 400 boundary points forming an ellipse A with radius r_Ax = r_Ay=20 located at the center of a 2D grid with Lx = Ly = 100 cells, ρ_0ij = 0.05, τ = 1.5 and ν = 1/3. The parameters used in Eq. (17) are k_f = 0.01, k_c=0.01, k_γ get two values in different experiments (k_γ = 0.0 and k_γ=0.1). The ellipse is immersed in a viscous fluid that moves from left to right forced by an imposed pressure gradient (i.e., ρ_1,i=0.055, ρ_100,i=0.045).

Figure 4 shows the position of the IB at iterations 100, 200, 400 and 800. Figure 5 shows the positions of the main radio of the ellipse in each case. The IB oscillates during a transient period and then remains stable after about 1500 steps.

Figure 4. Temporal evolution of an elliptical deformable immersed boundary (X_k) in a viscous fluid (L_x=L_y=100 cells, ρ_0ij=0.05, τ = 1.5, k_f = 0.01, k_c = 0.01, k_γ = 0.0.)

Figure 5. Temporal evolution of horizontal and vertical radius of the ellipse (L_x=L_y=100 cells, ρ_0ij=0.05, τ = 1.5)

B. Poiseuille flow

The plane Poiseuille flow is a pressure-driven flow in a duct of length L and width H (Fig. 6). We assume p₁ > p₂, and that pressure is uniform in the lateral direction. The analytical solution of the flow is:

(20)

Figure 6. Posseuille flow in a rectangular channel.

where the pressure gradient G is a constant related to the centerline velocity u₀ by

(21)

The Reynolds number is defined as

(22)

In LBM, the pressure is related to the density by

(23)

Two IB of 100 points each are located into a lattice to represent the channel walls as shown in Fig. 7. The grid dimensions are 110 cells long and 61 cells wide. The IB lines are separated 5.5 cells away from the grid limits, delimiting between them a channel of H = 50 cells.

Figure 7. Representation of a channel by means of IB.

The initial density ρ₀ was set at 1.0 and Δρ=0.00001 corresponding to a Δp/Δx = 1/3×10^-7. The relaxation parameter τ was set to 1.5 resulting in a cinematic viscosity ν = 1/3. According to Eqs. (21) and (22) the velocity in the center of the channel is 0.003125 and Re = 46.87. Figure 8 compares the calculated velocity profile with Eq. (20). It shows that the IB is able to represent the rigid walls of the channel correctly. The order of convergence was studied calculating the dependence of the error with the distance between the IB points, Δs (Fig. 9). The error was calculated as the maximum value of the difference between the numerical and the analytical solutions. The dashed line shown in the graphic is a fitting showing that the order of convergence is 2 as expected (Cheng and Zhang, 2010).

Figure 8. Comparison between the analytical and numerical velocity profile of a Poiseuille flow (k_f = 0.01, k_c = 0.0, k_γ = 0.0).

Figure 9. Order of convergence for the Poisseuille channel. The error is the maximum absolute value of the difference between the numerical and the analytical solutions. Δs is the distance between IB points. The line corresponds to a dependence Δs².

C. Vortex shedding

The third case studied is a 2D channel flow with a cylindrical obstacle, downstreaming a street of von Karmán vortexes. The channel size is 256×128 cells. The pressure drop is Δp = 5 10^-4, the average density is ρ_o = 0.05 and τ = 0.65. The obstacle is a circle of diameter d = 22 centered at position (64, 64) inside the channel. Its perimeter is represented by an IB of 200 points, which guarantees non-slip boundary condition with no fluid leakage having a distance between points Δs≤h/2 (Cheng and Zhang, 2010).

The empirical correlation (valid for Re >100) between the frequency of vortex detachments and the Reynolds number is given by the Strouhal number fd/V (Von Karman, 1964):

(24)

where ƒ is the vortex shedding frequency, V is the fluid velocity and Re is the Reynolds number defined as:

(25)

Figure 10 shows the dependence of the Strouhal number with the Reynolds number. It can be seen that the method is able to follow the experimental trend. Fig. 11 shows the temporal evolution of the y-profile of the velocity module (in grey scale) at x =192. The vortex shedding can be recognized after time step 2500. Figure 12 shows the contour map of the velocity module at time step 5000, where the vortex train is evidenced. Finally, the IB spatial resolution influence over the St was studied, in Fig. 13 can be observed the convergence to a specific value.

Figure 10. Dependence of the Strouhal number with the Reynolds number (d=22). Empirical correlation (line) and Lattice-Boltzmann calculation (points).

Figure 11. Temporal evolution of the y-profile of |u| at x =192, represented in grey scale (darker means lower velocity).

Figure 12. Contour map of |u| at time step 5000.

Figure 13. Influence of the IB resolution over the St.

VIII. CONCLUSIONS

An efficient algorithm to implement the immersed boundary method in a Lattice Boltzmann scheme for the Navier-Stokes equations was presented. The IB permits the generalization of the rigid regular grids used required by Lattice Boltzmann to domains with complex boundaries. In addition, it was shown that the present implementation is capable of simulating flows with deformable particles and solid-structure interaction problems, also reproducing flow instabilities.

REFERENCES
1. Chen, S. and G.D. Doolen, "Lattice Boltzmann Method for Fluid Flows," Annual Reviews Fluid Mechanics, 30, 329-364 (1998).         [ Links ]
2. Cheng, F. and H. Zhang, "Immersed boundary method and lattice Boltzmann method coupled FSI simulation of mitral leaflet flow," Computers & Fluid, 39, 871-881 (2010).         [ Links ]
3. Fadlun, E.A., R. Verzicco, P. Orlandi and J. Mohd-Yusof, "Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations," J. Comput. Phys., 161, 35-60 (2000).         [ Links ]
4. Feng, Z.G. and E.E. Michaelides, "The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems," J Comp Phy., 195, 602-628 (2004).         [ Links ]
5. Ge, L. and F. Sotiropoulos, "A numerical method for solving the 3D unsteady incompressible Navier-Stokes equations in curvilinear domains with complex immersed boundaries," J. Comput. Phys., 225, 1782-1809 (2007).         [ Links ]
6. Gilmanov, A. and F. Sotiropoulos, "A hybrid Cartesian/ immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies," J. Comput. Phys., 207, 457-492 (2005).         [ Links ]
7. Guo, Z., C. Zheng and B. Shi, "An extrapolation method for boundary conditions in lattice Boltzmann method," Phys. Fluids, 14, 2007-2010 (2002)         [ Links ]
8. Johansen, H. and P. Colella, "A Cartesian grid embedded boundary method for Poisson's equation on irregular domains," J. Comput. Phys., 147, 60-85 (1998).         [ Links ]
9. Mittal, R. and G. Iaccarino, "Immersed boundary methods," Annu. Rev. Fluid Mech., 37, 239-261 (2005).         [ Links ]
10. Peskin, C.S., "Numerical analysis of blood flow in the heart," J. Comput. Phys., 25, 220-252 (1977).         [ Links ]
11. Peskin, C.S. "The immersed boundary method," Acta Numer, 11, 479-517 (2002)         [ Links ]
12. Sukop, M. and D. Thorne, Lattice Boltzmann Modeling, Springer (2006).         [ Links ]
13. Udaykumar, H.S., R. Mittal, P. Rampunggoon and A. Khanna, "A sharp interface Cartesian grid method for simulating flows with complex moving boundaries," J. Comput. Phys., 174, 345-380 (2001).         [ Links ]
14. Von Karman, T., Aerodynamics, McGraw-Hill (1964).         [ Links ]
15. Zhang, X., D. Schmidt and B. Perot, "Accuracy and conservation properties of a three-dimensional unstructured staggered mesh scheme for fluid dynamics," J. Comput. Phys., 175, 764-791 (2002).         [ Links ]
16. Zhang, X, S.Z. Ni and G.W. He, "A pressure-correction method and its applications on an unstructured Chimera grid," Comput. Fluids, 37, 993-1010 (2008).         [ Links ]
17. Zou, Q. and X. He, "On pressure and velocity boundary conditions for the lattice Boltzmann BGK model," Phys. Fluids, 9, 1591-1598 (1997)         [ Links ]

Received: December 1, 2011
Accepted: October 5, 2012
Recommended by Subject Editor: Adrián Lew