An immersed interface-lattice Boltzmann method for fluid-structure interaction
Jianhua Qin Ebrahim M. Kolahdouz Boyce E. Griffith
physics.comp-ph cs.NA math.NA physics.flu-dyn



1,2,*]Jianhua Qin
2]Ebrahim M. Kolahdouz
3,**]Boyce E. Griffith


[1]Key Laboratory of Transit Physics, 
Nanjing University of Science  and Technology, Nanjing, China
[2]Department of Mathematics, University of North Carolina, Chapel Hill, NC, USA
[3]Departments of Mathematics, Applied Physical Sciences, and 
Biomedical Engineering, 
University of North Carolina, Chapel Hill, NC, USA
[*]Corresponding author at: Key Laboratory of Transit Physics, Nanjing University of Science and Technology, Nanjing, China.
[**]Corresponding author at: Departments of Mathematics, Applied Physical Sciences, and Biomedical Engineering,
University of North Carolina, Chapel Hill, NC, USA.
[•]Email: njustqin@gmail.com (J. Qin), boyceg@email.unc.edu (B. E. Griffith)



An immersed interface-lattice Boltzmann method for fluid-structure interaction
    [
    

==============================================================================




An immersed interface-lattice Boltzmann method (II-LBM) is developed for modelling fluid-structure systems. The key element of this approach is the determination of the jump conditions that are satisfied by the distribution functions within the framework of the lattice Boltzmann method when forces are imposed along a surface immersed in an incompressible fluid.  
In this initial II-LBM, the discontinuity related to the normal portion of the interfacial force is sharply resolved by imposing the relevant jump conditions using an approach that is 
analogous to imposing the corresponding pressure jump condition in the incompressible Navier-Stokes equations.
We show that the jump conditions for the distribution functions are the same in both single-relaxation-time and multi-relaxation-time LBM formulations. 
Tangential forces are treated using the immersed boundary-lattice Boltzmann method (IB-LBM). The performance of the II-LBM method is compared to both the direct forcing IB-LBM for rigid-body fluid-structure interaction, and the classical IB-LBM for elastic interfaces.
Higher order accuracy is observed with the II-LBM as compared to the IB-LBM for selected benchmark problems. Because the jump conditions of the distribution function also satisfy the continuity of the velocity field across the interface, the error in the velocity field is much smaller for the II-LBM than the IB-LBM.
The II-LBM is also demonstrated to provide superior volume conservation when simulating flexible boundaries. 




§ INTRODUCTION

The immersed methods are popular in modelling fluid-structure interaction (FSI) because they provide a computationally simple approach to deal with complex interfacial geometries and large deformations <cit.>. For FSI problems involving thin interfaces, immersed methods typically use an Eulerian description of the fluid and a Lagrangian description of the interface. The immersed boundary (IB) method is one of the earliest of these types of methods. It was introduced by Peskin and uses integral transforms with Dirac delta function kernels to couple Eulerian and Lagrangian variables <cit.>. The immersed interface method (IIM) was developed by LeVeque and Li to sharply impose jump conditions in elliptic interface problems <cit.>, and they extended the IIM to solve incompressible Stokes flow problems with flexible boundaries <cit.>. Peskin and Printz <cit.> and Lai and Li <cit.> established the jump conditions for the velocity and pressure in the incompressible Navier-Stokes equations in two and three spatial dimensions, respectively, and Li and Lai extended the IIM to treat the incompressible Navier-Stokes equations <cit.>. Lee and LeVeque developed an alternative approach that combined the IB method and the IIM to model a flexible membrane <cit.>. To facilitate higher-order implementations of the IIM, Xu and Wang systematically derived both spatial and temporal jump conditions for the incompressible Navier-Stokes equations <cit.>, and applied the IIM to model rigid-body FSI <cit.>. Le et al. used the IIM to simulate multiple rigid and flexible structures <cit.>. Meyer et al. combined the IIM and large eddy simulation to capture flow dynamics near fluid-structure interfaces under high Reynolds numbers and turbulent flows using a finite volume method <cit.>. Kolahdouz et al. introduced a version of the IIM for general geometries described by C^0-continuous surface representations <cit.>.

Prior to this work, the immersed boundary-lattice Boltzmann method (IB-LBM) has been adopted to simulate fluid-structure interaction, including rigid body dynamics <cit.>, elastic filaments <cit.>, flexible sheets <cit.> and large deformation of flexible beams <cit.>. The IB-LBM maintains the simplicity of the IB method, and the lattice Boltzmann method (LBM) <cit.> uses the linear lattice Boltzmann equation (LBE) to provide a nearly-incompressible fluid model that is simpler to parallelize than the incompressible Navier-Stokes equations because it permits a fully explicit time stepping scheme that does not require the solution of any global systems of equations. These previous IB-LBM are mostly based on so-called diffusive interface IB methods, however, which effectively regularize discontinuities that generally are present at the immersed boundary. Kang and Hassan proposed a sharp interface IB-LBM by imposing the force density on exterior solid nodes nearest to the boundary <cit.>. However, the physical problems studied in their work are only related to static geometries, and it is not straightforward to generalize their approach to model flexible boundaries. This paper introduces an alternative approach to integrating sharp interface IB approaches with LBM-based IB formulations, resulting in a methodology that is applicable to a wider range of problems. 

In this work, we first derive the jump conditions of the distribution function in the lattice Boltzmann equation that are related to applying a normal force along an immersed boundary. The jump conditions are then used to formulate an immersed interface-lattice Boltzmann method (II-LBM). The jump conditions arising from the normal component of the force are associated with pressure jump conditions in the conventional IIM for the incompressible Navier-Stokes equations. In this initial IIM-LBM formulation, the tangential forces are spread to the background grid using a standard diffuse interface IB approach. In addition, the velocity of the interface is determined using a standard IB interpolation scheme. The resulting methodology is similar to that of Lee and LeVeque <cit.>, except that here we use the LBM to model the fluid flow instead of the incompressible Navier-Stokes equations. Empirical results show that the II-LBM can achieve higher order of accuracy than the IB-LBM. 
Moreover, the comparison between the II-LBM and the IB-LBM using the simple external force term <cit.> and Guo's external force term <cit.> shows that II-LBM has substantially superior volume conservation than the IB-LBM.



§ SINGLE- AND MULTI-RELAXATION-TIME IB-LBM 

This section introduces formulations of single-relaxation-time (SRT) and multi-relaxation-time (MRT) lattice Boltzmann methods and describes the integration of LBM formulations with conventional IB methods. Although the numerical tests considered in this paper are based on the more stable MRT model, the SRT model is included for completeness. Fig. <ref> provides a schematic of a membrane immersed in a background fluid. Throughout this paper, Ω is the whole fluid domain, Γ_t represents the immersed boundary at time t (shown in Fig. <ref>(a)), and the exterior and interior of the membrane at time t are defined as Ω_t^+ and Ω_t^-, respectively, so that Ω = Ω^+_t ∪Ω^-_t (see Fig. <ref>(b)).



 §.§ The immersed boundary-lattice Boltzmann method (IB-LBM)



  §.§.§ Governing equations of the LBM

Before going through the IB-LBM, the LBE governing the fluid flow is introduced. Here, the LBM based on both the SRT operator and the MRT operator are presented. 

The equations for the distribution function f_i of the SRT operator lattice Boltzman equation is

    ∂ f_i/∂ t+e_i·∇ f_i=
    -1/λ(f_i-f_i^(eq)), i=0,1,2,...,Q,

with Q=18 in the three-dimensional model (D3Q19), or by Q=8 in the two-dimensional model (D2Q9). Here, e_i are velocity vectors 
defined by (0,0,0), (±1,0,0), (0,±1,0), (0,0,±1), (±1,±1,±1) for the D3Q19 model and (0,0), (±1,0), (0,±1), (±1,±1) for the D2Q9 model. 
Fig. <ref> shows the e_i vectors associated with the identified point P_1 for the D2Q9 model. λ is the relaxation parameter, and f_i^(eq) is the equilibrium distribution function.

[figure]style=plain,subcapbesideposition=top



By using the MRT operator in the lattice Boltzmann equation <cit.>, the evolution of the distribution functions becomes

    ∂ f_i/∂ t+e_i·∇ f_i=-M^-1S_d(m_i(x,t)-m_i^(eq)(x,t)),

in which m_i represents the moments of the distribution function f_i. 
The matrix operator M relates the distribution functions to the moments of the distribution functions via m=Mf, in which m={m_i} and f={f_i}. 
Let m^(eq)={m_i^(eq)} and f^(eq)={f_i^(eq)} be the moments of equilibrium distribution functions and distribution functions, respectively. They are also related by using M as m^(eq)=Mf^(eq). S_d=diag(s_0,s_1,s_2,...,s_Q)
is a diagonal matrix, the values of the elements in S_d will be given later. Although the SRT based LBM is extremely simple, the MRT based LBM is generally more robust <cit.>.



  §.§.§ Governing equations of the IB-LBM

The equations for the immersed boundary-lattice Boltzmann method (IB-LBM) <cit.> based on the simple external force term <cit.> are

    ∂ f_i/∂ t+e_i·∇ f_i=RHS+
    1/c_s^2ω_ie_i·g(x,t),


    g(x,t)=
    ∫_ΓG(r,s,t) δ(x-X)
     drds,


    ∂X(r,s,t)/∂ t=
    U(r,s,t)=
    ∫_Ωu(x,t)δ (x-X(r,s,t)) dx,


    G(r,s,t)=𝒢(X,U,r,s,t).

Here, Eq. (<ref>) accounts for the momentum of the fluid. The form of RHS is determined by whether the SRT or MRT operator is used. For the SRT operator, RHS = -1/λ(f_i-f_i^(eq)), as in Eq. (<ref>), and for the MRT operator, RHS = -M^-1S_d(m_i(x,t)-m_i^(eq)(x,t)), as in Eq. (<ref>). The physical configuration of the membrane Ω_t is taken as X(r,s,t), in which (r,s) are curvilinear coordinates that parameterize the interface. The Lagrangian force density G(r,s,t) and the Eulerian force density g(x,t) are related to each other via Eq. (<ref>). δ(x-X) is the Dirac delta function. In the numerical tests, we use Peskin's four-point regularized delta function <cit.>. 𝒢 is the functional that determines the interfacial force from the deformations and/or velocities of the immersed interface.



 §.§ Numerical approaches for the IB-LBM

In the discrete equations, Δ x and Δ t are the grid spacing of the fluid and the time step size, respectively. We always use lattice units, so that Δ x=Δ t=1. The grid spacing between Lagrangian points in lattice units is defined as Δ X.


  §.§.§ Discrete LBE

The LBM is used to discretize the LBE for modelling the incompressible fluid flow. The relaxation parameter λ is related to the kinematic viscosity ν by

    ν=c_s^2(λ/Δ t-1/2)Δ t,

in which c_s=1/√(3) is the speed of sound in lattice units, and ω_i is the weighting function defined by

    ω_i=
    1/3,   i=0,
    1/18,   1≤ i ≤6,
    1/36,   otherwise,

for the D3Q19 model, or by

    ω_i=
    4/9,   i=0,
    1/9,   1≤ i ≤4,
    1/36,   otherwise,

for the D2Q9 model.

For an incompressible fluid, the modified equilibrium distribution function <cit.> is:

    f_i^(eq)=ω_i[1/c_s^2p+e_iu/c_s^2
    +(e_iu)^2/2c_s^4-u^2/2c_s^2].

The macroscopic pressure p and velocity u are obtained by

    p=c_s^2∑_i=1^Qf_i

and

    u=∑_i=1^Qe_i f_i.

The gradient of the distribution function (∇ f_i) is discretized via Grad(f_i) and is defined by

    Grad(f_i)=f_i(x,t)-f_i(x-Δx_i,t)/Δx_i,

in which Δx_i=e_iΔ t is the grid spacing in the direction of e_i. The discretized lattice Boltzmann equation with SRT becomes

    f_i(x,t+Δ t)-f_i(x,t)/Δ t+
    e_i·f_i(x,t)-
    f_i(x-Δx_i,t)/Δx_i=
    -1/λ(f_i-f_i^(eq)).

For the D3Q19 MRT model <cit.>, s_0=s_3=s_5=s_7=0,  s_4=s_6=s_8,  s_9=s_11,  s_10=s_12,  s_13=s_14=s_15,  s_16=s_17=s_18. The nonzero elements in the diagonal matrices are s_1=1.19 Δ t,  s_2=s_10=1.4 Δ t,  s_4=1.2 Δ t,  s_9=s_13=Δ t/λ, and s_16=1.98 Δ t. In the D2Q9 MRT model <cit.>, s_0=s_3=s_5=0, s_1=1.63 Δ t, s_2=1.14 Δ t, s_4=s_6=1.92 Δ t, and s_7=s_8=Δ t/λ. Then the discretized lattice Boltzmann equation with MRT becomes

    f_i(x,t+Δ t)-f_i(x,t)/Δ t+
    e_i·f_i(x,t)-
    f_i(x-Δx_i,t)/Δx_i=
    -M^-1S_d(m_i(x,t)-m_i^(eq)(x,t)).




  §.§.§ Discrete IB-LBM

The lattice Boltzmann equation with the simple external force term <cit.> (i.e. Eq. (<ref>)) is discretized by

    f_i(x,t+Δ t)-f_i(x,t)/Δ t+
    e_i·f_i(x,t)-
    f_i(x-Δx_i,t)/Δx_i=
    RHS+1/c_s^2ω_ie_i·g(x,t),

in which the Eulerian force is determined via

    g(x,t)=∑_kG_k δ_h(x-X_k) Δ X.

Here, δ_h(x-X) is a regularized delta function. 

In the direct forcing IB method <cit.>, the Lagrangian force in Eq. (<ref>) is determined by

    G_k(t)=
    η(U_k^D(t)-U_k^*(t)),

in which G_k(t), U_k^D, and U_k^*=
∑_xu^* δ_h(x-X_k) Δ x^2 represent the force density, the desired Lagrangian velocity, and the interpolated Lagrangian velocity for the k^thth (k = 1,2,⋯,N) Lagrangian point X_k, respectively, and u^* is an intermediate Eulerian velocity without the effect the immersed boundary forcing. The penalty parameter η is chosen as η=1/Δ t.

In the classical IB method for a flexible membrane <cit.>, the force on the k^th Lagrangian node is calculated by

    G_k(t)=
    κX_k+1(t)-2X_k(t)
    +X_k-1(t)
    /Δ X^2,

in which κ is a spring constant, and the dynamics of X_k are obtained by discretizing Eq. (<ref>) via

    X_k(t+Δ t)-X_k(t)/Δ t=U_k(t)=
    ∑_xu(x,t) δ_h(x-X) Δ x^2.


The IB-LBM based on Guo's external force terms <cit.> is used for some comparisons with the II-LBM introduced in this paper. Guo's external force terms have been demonstrated to recover the compressible Navier-Stokes equations more accurately. In their formulation, the last term in Eq. (<ref>) is replaced by

    G_i^S=(1-Δ t/2λ)ω_i(e_i-u/c_s^2+e_i·u/c_s^4·e_i)·gΔ t
 
for SRT-LBM and by

    G^M=M^-1(I-1/2S_d)MGΔ t
 
for MRT-LBM, in which G^M=(G_0, G_1,..., G_q)^T and G=(G_0, G_1,..., G_q)^T, with G_i=ω_i(e_i-u/c_s^2+e_i·u/c_s^4·e_i)·g. Here, I is the identity matrix. In this case, the equilibrium function and velocity are

    f_i^(eq)=ω_ip/c_s^2[1+e_iu/c_s^2
    +(e_iu)^2/2c_s^4-u^2/2c_s^2]

and

    p/c_s^2u=∑_i=1^Qe_i f_i+Δ t/2g,

in which p is obtained by Eq. (<ref>).



 §.§ Rigid body dynamics

For the numerical tests involving both the translation and rotation of rigid bodies in two spatial dimensions, the total force and torque acted on the structure are calculated by

    F^ext(t)=-∑_k=1^N G_k(t)Δ X
    +ρ_f V_sU_c(r,s,t)-U_c(r,s,t-Δ t)/Δ t,


    T^ext(t)=   -∑_k=1^N  (X_k(r,s,t)-X_c(r,s,t))
    ×G_k(t)Δ X
       +ρ_f/ρ_s I_sW(r,s,t)-W(r,s,t-Δ t)/Δ t,

in which ρ_f is the density of the fluid, ρ_s is density of the structure, V_s is the volume of the structure, I_s is the moment of inertia, X_c(r,s,t), U_c(r,s,t), and W(r,s,t) are the center of mass, velocity of the center of mass, and angular velocity of the structure at time t, respectively. The last terms in Eqs. (<ref>-<ref>) are to eliminate the effect of inertial forces from within the immersed interface <cit.>.

After obtaining F^ext and T^ext, the acceleration and angular acceleration of the structure are determined by

    ma(t+Δ t)=F^ext(t),
    
    I_s w(t+Δ t)=T^ext(t),

in which m, a, I_s, and w are the mass, acceleration, moment of inertia, and angular acceleration of the rigid body, respectively. Then the velocity and angular velocity of the rigid body are updated via

    U_c(t+Δ t)=U_c(t)+Δ t a(t+Δ t),
    
    W(t+Δ t)=W(t)+Δ t  w(t+Δ t).

Ultimately, the new position and orientation angle of the structure are determined by

    X_c(t+Δ t)   =X_c(t)+U_c(t)Δ t+1/2(Δ t)^2(a(t)+a(t+Δ t)),
    θ(t+Δ t)   =θ(t)+W(t)Δ t+1/2(Δ t)^2(w(t)+w(t+Δ t)),

in which θ(t) is the orientation angle of the rigid body at time t.



§ THE IMMERSED INTERFACE-LATTICE BOLTZMANN METHOD (II-LBM)

The II-LBM depends on the jump condition across the fluid-solid interface. This section formulates the jump conditions in the framework of the LBE. The resulting II-LBM can be viewed as a correction of the IB-LBM considering the formulated jump conditions. 


 §.§ The derivation of the jump conditions

For a scalar field ψ, the jump is defined as

    ψ(X,t)
    =   lim_ϵ ↓ 0ψ(X+ϵn(X,t),t)-
    lim_ϵ ↓ 0ψ(X-ϵn(X,t),t)
    
    
    =   ψ^+(x,t)-ψ^-(x,t),

in which ψ(X,t) is a position along the immersed interface, · is the operator for the jump, and ψ^+(x,t) and  ψ^-(x,t) are limiting values of ψ approaching interface position X from exterior region Ω_t^+ and interior region Ω_t^-, respectively. The unit normal vector at the interface is

    n=∂X/∂ r×∂X/∂ s/|∂X/∂ r×∂X/∂ s|.

For rigid body simulations, |∂X/∂ r×∂X/∂ s| = 1.

The jump conditions of the distribution function based on the simple external force term <cit.> are

    f_i=
    ω_i/c_s^2G(r,s,t)·n/|∂X/∂ r×∂X/∂ s|.

Following Lai and Li <cit.>, we can derive the above jump conditions in three spatial dimensions. Let us choose a banded domain Ω_ϵ,t enclosing the immersed membrane Γ_t with outer and inner subregions Ω_ϵ,t^+ and Ω_ϵ,t^-. Here we denote ϵ as a small distance parameter. Suppose ϕ(x) is any smooth function with compact support on Ω_ϵ,t. Multiplying ϕ(x) on both sides of Eq. (<ref>) and integrating over Ω_ϵ,t, we obtain 

    ∫_Ω_ϵ,t∂ f_i/∂ t ϕ dx+
     ∫_Ω_ϵ,t(e_i·∇ f_i)
    ϕ dx
    
    =∫_Ω_ϵ,tRHS ϕ dx
    +    1/c_s^2ω_ie_i·∫_Ω_ϵ,t(∫_0^L_r∫_0^L_sG(r,s,t)δ (x-X(r,s,t))
     drds)ϕ dx,

in which L_r and L_s are the limits of the parameterized coordinate systems. The first term and third term in the above equation go to zero as ϵ→ 0 because f_i and f_i^(eq) are continuous and bounded. Using the fact that Ω_ϵ,t encloses Γ_t,  the last term in the above equation can be simplified, and the Eq. (<ref>) becomes

    ∫_Ω_ϵ,t(e_i·∇ f_i)
    ϕ dx=
    1/c_s^2ω_ie_i·∫_0^L_r∫_0^L_sG(r,s,t) ϕ drds.

By the divergence theorem, the first term in Eq. (<ref>) is

    e_i∫_Ω_ϵ,t(∇ f_i)ϕ dx=
    e_i
    (
    ∫_Ω_ϵ,t^+ f_inϕ da+
    ∫_Ω_ϵ,t^- f_i(-n)ϕ da-
    ∫_Ω_ϵ,tf_i ∇ϕ dx).

in which da=|∂X/∂ r×∂X/∂ s|drds is the area of the surface element. The last term in Eq. (<ref>) converges to zero as ϵ→ 0 because f_i is bounded.
Then Eq. (<ref>) becomes

    e_i∫_Ω_ϵ,t(∇ f_i)ϕ dx=
    e_i∫_0^L_r∫_0^L_s f_inϕ|∂X/∂ r×∂X/∂ s|
    drds.

Combining Eq. (<ref>) and Eq. (<ref>) and taking advantage of the fact that ϕ is arbitrary, we obtain

    f_in|∂X/∂ r×∂X/∂ s|=ω_i/c_s^2G(r,s,t),    i≠ 0.

Since n is the unit normal vector,

    f_i|∂X/∂ r×∂X/∂ s|=ω_i/c_s^2n·G(r,s,t),    i≠ 0.

Then the only unknown jump condition for the distribution function is
f_0. Apply the jump operator to both sides of Eq. (<ref>):

    p = c_s^2∑_if_i=c_s^2∑_i f_i
           =
    G(r,s,t)·n/|∂X/∂ r×∂X/∂ s|∑_i≠ 0ω_i+c_s^2 f_0
           =
    G(r,s,t)·n/|∂X/∂ r×∂X/∂ s|-
    G(r,s,t)·n/|∂X/∂ r×∂X/∂ s| ω_0+c_s^2 f_0.

Because Eq. (<ref>) only deals with the normal part of the boundary forces, the equations for the jump conditions of the distribution function should correspond to the jump condition of the pressure <cit.>,

    p=G(r,s,t)·n/|∂X/∂ r×∂X/∂ s|.

From the above two equations, the jump condition f_0 is

    f_0=
    ω_0/c_s^2G(r,s,t)·n/|∂X/∂ r×∂X/∂ s|.

Combining Eq. (<ref>) and Eq. (<ref>) gives the jump condition for the distribution functions in Eq. (<ref>). The velocity jump can be determined via

    u=∑_if_ie_i=∑_ie_i f_i=
    ∑_ie_i
    ω_i G(r,s,t)·n/c_s^2|∂X/∂ r×∂X/∂ s|.

After rearrangement, we have

    u=
    G(r,s,t)·n/c_s^2|∂X/∂ r×∂X/∂ s|∑_ie_iω_i.

Because the weighting function ω_i is the same for opposite directions of e_i, we have Σ_ie_iω_i=0. This allows us to reduce Eq. (<ref>) to u=0, which is the continuity of the velocity field across the membrane.



 §.§ Correction of the IB-LBM considering the jump conditions

Jump conditions are imposed in the discretized SRT-LBM by including appropriate correction terms in Eq. (<ref>). Because we consider only the normal component of the force in determining the jump conditions, this is similar to the application of the pressure jump condition on the gradient of the pressure when adopting the finite difference method to solve the incompressible Navier-Stokes equations <cit.>. The discretized gradient of the distribution function (Grad(f_i) in Eq. (<ref>)) can be written considering the jump conditions as

    Grad(f_i)=f_i(x,t)-f_i(x-Δx_i,t)/Δx_i-
     f_i(x,t)/Δx_i.

Therefore, the discretized lattice Boltzmann equation with the simple external force term <cit.> (Eq. (<ref>)) should be corrected by

    f_i(x,t+Δ t)-f_i(x,t)/Δ t+
    e_i·(
    f_i(x,t)-
    f_i(x-Δx_i,t)/Δx_i-
     f_i(x,t)/Δx_i)
    =
         RHS+1/c_s^2ω_ie_i·g^||(x,t),

in which g^||(x,t) is force of the immersed boundary node in the tangential direction. Because the jump conditions determined here for the distribution function only impose the pressure jump condition associated with the normal component of the interfacial force. In the present method, the tangential force is treated via a discretized integral transform with a regularized delta function kernel, as in the conventional IB and IB-LB methods. The external force density on the fluid node is obtained by spreading the Lagrangian force on the IB node

    g^||(x,t)=
    ∑_kG^||_k(t)δ
    (x-X)Δ X,

in which G^||_k=G_k-
G_k·n is the tangential part of the Lagrangian force.
Eq. (<ref>) can be rewritten as

    f_i(x,t+Δ t)=
    f_i(x-Δx_i,t)+RHSΔ t+
     f_i(x,t+Δ t)
    +1/c_s^2ω_ie_i·g^||(x,t+Δ t)Δ t.


In those cases where Guo's external force terms <cit.> are used, the Eulerian force g used to calculate G_i^S in Eq. (<ref>) and G̃_i in Eq. (<ref>) should be replaced by g^||. For the point P_1 in Fig. <ref>, only the distribution functions f_3 and f_7 are updated by Eq. (<ref>) with f_i(x,t)≠ 0. For the point P_2 in Fig. <ref>, only the distributions functions f_1, f_5 and f_8 are updated considering the jump conditions.

To implement the jump conditions of the distribution functions on a certain fluid node in Eq. (<ref>), the corresponding Lagrangian force G(r,s,t) on the structure surface, and the normal direction of the fluid node to the structure surface have to be determined. In this paper, the level set method <cit.> is used when simulating rigid body simulations, and a cubic spline representation is adopted in simulating flexible boundaries <cit.>.









§ RESULTS

This section provides results from several benchmark problems. For the problems studied in this paper, Guo's external force terms are only used in Sec. 4.2 and the multi-relaxation-time operator LBM is used for all the cases.


 §.§ Uniform normal force along a circle

We first consider a circular cylinder immersed in the fluid with evenly distributed force along the normal direction of the surface. This case is the same as the one considered by Li and Lai <cit.>. The schematic of the present calculation is shown in Fig. <ref>. The diameter of the cylinder is D. The fluid domain is [-D,D]× [-D,D] with the cylinder located at its center. The boundaries of the outer domain are considered as outflow and handled by the non-equilibrium extrapolation strategy <cit.>. The Reynolds number and the viscosity are chosen as Re=100 and ν=0.05. The characteristic velocity is U=ν Re/D, and the magnitude of the normal interfacial forces is chosen to be F^ext/ρ_f U^2=1. By varying U and D, different cases are obtained. The normalized L_2-norm and L_∞-norm are defined by

    L_2-norm=√((∑_x 
    (u_x-u_x^e)^2+(u_y-u_y^e)^2)/N_f)/U,
    
    L_∞-norm=max_x max(|u_x-u_x^e|, |u_y-u_y^e|)
    /U,

in which N_f is the number of fluid points of the domain Ω and u_x and u_y are fluid velocities in the x and y directions, respectively. The superscript “e" indicates the exact velocity in the fluid domain. In this case, the exact velocity is u_x=u_y=0.



Errors are assessed at time t^end=D/U. The results of the error and the order of accuracy calculated by the IIM and the direct forcing IB method in this case are shown in Table <ref>. The pressure and velocity fields calculated by using the IB method and the IIM at t^end=D/U are presented in Fig. <ref> and <ref>, respectively. The pressure fields clearly show that the IIM gives piecewise constant pressure whereas the IB method regularizes the pressure discontinuity. The velocity magnitude obtained by the IIM in Fig. <ref>(b) is much smaller than that in Fig. <ref>(a) obtained by the IB method. This is because the jump conditions of the distribution function balance the normal forces while satisfying the continuity of the velocity field across the interface in Eq. (<ref>). On the other hand, the IB method corrects the velocity fields to balance to boundary forces.
  




 





 §.§ Dynamic motion of a thin elastic interface

For the incompressible Navier-Stokes equations, using a staggered-grid spatial discretization <cit.> or specialized differencing operators <cit.> can greatly reduce spurious volume (area) change resulting from the non-divergence free velocity field of the Lagrangian points by the IB method. However, the lattice Boltzmann equation is generally solved with an explicit time integration scheme using a collocated grid. To our knowledge, a staggered-grid LBM has not been developed. Lee and LeVeque <cit.> found that the immersed interface method with the pressure jump condition has better volume conservation than the IB method. To validate that the II-LBM proposed in this paper also preserves the volume for flexible bodies, the dynamic motion of a thin elastic interface is studied.

The schematic of the computational domain is shown in Fig. <ref>, where the initial shape is an ellipse interface and the equilibrium shape is expected to be a circular interface. A square box of side length L is chosen as the computational domain. The boundaries of the outer domain are set to no-slip walls by using the non-equilibrium extrapolation strategy <cit.>. The ellipse interface is initially positioned at (0.5D,0.5D) with semi-major and semi-minor axes of a = 0.35L and b = 0.25^2/0.35L, respectively. The nondimensional spring constant and viscosity are set to κ=κ/U^2L=0.1 and ν=0.02, respectively. The resulting spring constant κ will be used in Eq. (<ref>) to calculate the spring force. The Reynolds number is set to Re=100, and the characteristic velocity is determined as U=ν Re/L. The distance between two adjacent Lagrangian points are taken as Δ X=1/3Δ x for the IB-LBM and Δ X=2Δ x for the II-LBM with respect to the circular interface in its equilibrium state of radius 0.25L. Three different meshes, corresponding to L=80Δ x, 160Δ x, and 240Δ x, are used. The first period of oscillation ends at approximately tU/L=0.553.0625^log_10 (κ)/κ, which is the same as that indicated by Griffith <cit.>. As can be seen from the pressure distribution at tU/D=30 at y/L=0.5 calculated by L=240Δ x in Fig. <ref>, the II-LBM gives a sharp jump at the interface of x=0.25L whereas the other two methods produce a regularized pressure discontinuity. Fig. <ref> shows the changes of the length of the semi-major and semi-minor axes against three different approaches. As expected, both the IB-LBM with multi-relaxation-time using a simple external force term <cit.> and Guo's external force term <cit.> result in substantial volume loss. However, the volume of the ellipse barely changes in time by using the II-LBM, and the final result of L=240Δ x gives a 0.01% error of the semi axis as compared to the analytical result. As seen in Fig. <ref>, although the II-LBM shows some oscillations for the coarsest grid, the oscillations converge rapidly and become unnoticeable for finer grids. To quantify the three different approaches adopted in the performance of preserving volume, Table <ref> is presented to show the area loss per dimensionless time Δ V_loss and total area loss V^total_loss at t^end=30L/U

    Δ V_loss=Δ V_loss/LUΔ t,

and

    V^total_loss=V^total_loss/L^2,

in which Δ V_loss and V^total_loss are dimensional volume of loss per unit time and total volume loss at time t.




  

 









 §.§ Uniform flow around a circular cylinder

Flow around a circular cylinder is frequently used as a benchmark to check the accuracy of a numerical algorithm. The left and right boundaries are taken as inflow and outflow, respectively. The inflow velocity is U. Freestream boundary conditions are used on the side walls. The computational domain is [0,40D]×[0,20D] with the circular cylinder of diameter D located at (10D,10D). The Strouhal number is defined as

    St=fD/U,

in which f is the frequency of the vortex shedding. The drag and lift coefficients are defined as

    C_D=2F^ext_x/ρ_f U^2D,


    C_L=2F^ext_y/ρ_f U^2D,

in which ρ_f is the density of the fluid. F^ext_x and F^ext_y are external forces in the x and y directions, respectively.

Table. <ref> lists the drag coefficient and recirculation length (R_l) and angle of separation (Λ) for the steady flow cases of Re = 20 and 40. As the Reynolds number increases, the flow changes from steady flow (Re ≤ 47) to unsteady flow and the Von Kármán vortex street appears. To quantify the accuracy of the II-LBM in identifying the flow characteristics of unsteady flow, a comparison of average drag coefficient (C_D), lift coefficient (C_L) and Strouhal number (St) is presented in Table. <ref>. The pressure fields in Fig. <ref> show that the II-LBM can capture the discontinuities in the pressure field at the interface. Fig. <ref> present the time history of the drag and lift coefficients at Re=100 and 200. These figures indicate that the periodic vortex shedding is successfully simulated by the II-LBM.


  





  
[figure]style=plain,subcapbesideposition=top


  


 §.§ Prescribed translation of a circular cylinder in a resting background fluid

The study of forced oscillations of the structure is essential for its significance in predicting the mechanism of vortex-induced vibrations. Dütsch et al. <cit.> first studied the prescribed translation of circular cylinder by experiments and numerical simulations. In that work, a rigid circular cylinder with diameter D undergoes a prescribed harmonic motion in the x direction. The motion of the circular cylinder is specified with the following displacement and velocity in the x direction:

    D_c(t)=-Asin(U_maxt/A),
    
    U_c(t)=-U_maxcos(U_maxt/A),

in which D_c(t) and U_c(t) are the displacement and velocity of the cylinder center in the x direction, A is the maximum displacement, and U_max is the maximum translation velocity. The Keulegan-Carpenter number and Reynolds number are KC = 2π A/D and Re=U_maxD/ν, respectively. KC=5 and Re=100 are chosen in our present study to facilitate comparisons to Dütsch et al. <cit.>. The computational domain is the rectangular region [0,55D]×[0,35D], and the circular cylinder is initialized so that its center is located at the origin. Fig. <ref> compares the drag coefficients obtained by the numerical simulation of the II-LBM with the experiment of Dütsch et al. <cit.>. The drag coefficient is defined by Eq. (<ref>). Results obtained using the II-LBM are in good agreement with the prior experimental results of Eutsch et al.<cit.>.
[figure]style=plain,subcapbesideposition=top





 §.§ Vortex-induced vibrations of a circular cylinder in the longitudinal direction

In this section, a circular cylinder with diameter D immersed in a uniform flow is considered to move freely in the y direction. This case was proposed by Ahn and Kallinderis <cit.>, and later studied by Borazjani et al. <cit.> and Bao et al. <cit.>. As the cylinder is placed in the fluid, fluid traction will act on the cylinder due to the no-slip boundary condition of the cylinder surface. On the other hand, the motion of the cylinder will affect the surrounding fluid. Therefore, “two-way" interplay between the fluid and the structure is involved. The computational domain is the rectangular region [0,40D]×[0,20D], and the cylinder is initialized so that its center is located at (10D, 10D). The equation to model vortex-induced vibrations of a cylinder in the longitudinal direction is

    a_y+2ζ(2π/U_red)U_y+(2π/U_red)^2D_y=2C_L/πρ^*,

in which a_y, U_y, and D_y are dimensionless acceleration, velocity, and displacement of the cylinder and are normalized by U^2/D, U, and D, respectively. The reduced velocity is defined as U_red=U/f_ND, where f_N=√(k_sm)/2π is the natural frequency of the spring. k_s and m are the spring constant and mass of the cylinder, respectively. ζ=c/2√(k_sm) is the damping ratio. The lift coefficient C_L is determined in Eq. (<ref>). The density ratio ρ^* is the ratio between the solid density ρ_s and the fluid density ρ_f,

    ρ^*=ρ_s/ρ_f.

Simulation parameters are the same as those used by Ahn and Kallinderis <cit.>. Specifically, the density ratio and damping ratio are set to be 8/π and 0, respectively. The reduced velocity varies in the range of 3≤ U_red≤ 8. The maximum oscillation amplitude of the circular cylinder for the reduced velocity from 3 to 8 are 0.075D, 0.551D, 0.520D, 0.459D, 0.370D, and 0.084D, respectively. The comparison between the present results and previous published results is shown in Fig. <ref>.





 §.§ Sedimentation of an elliptical particle

In this section, the sedimentation of an elliptical particle in a confined channel proposed by Xia et al. <cit.> is considered. A schematic of the model is shown in Fig. <ref>. When sedimenting, the ellipse particle will translate in the x and y directions because of gravity, pressure, and viscous forces. The particle will also rotate. Therefore, this problem involves three degrees of freedom. The computational domain is a confined box of [-4a, 4a]×[-30a, 5a], in which a is the semi major axis of the ellipse particle, and the particle is initialized at (0,0) with an initial orientation angle of θ(t=0)=π/4. The semi-major and semi-minor axes are taken as a=0.025× 10^-2m and b=1/2a, respectively. The outer boundaries are zero-velocity boundary conditions, and handled by the non-equilibrium extrapolation strategy <cit.>. As in the previous section, the density ratio between the ellipse and the fluid is taken as ρ^*. The kinematic viscosity of the fluid is ν_p=0.01 m^2/s and gravitational acceleration is g_e=(0, -9.80) m/s^2. 



We first consider ρ^*=1.1. The grid number along the major axis is 2a=25Δ x, and the distance between adjacent Lagrangian points is the same as the lattice spacing. Fig. <ref> shows the position of the ellipse in the x direction and orientation angle of the ellipse against the position of the ellipse in the y direction. Good agreement with the finite element computations of by Xia et al. <cit.> is achieved. As can be seen in Fig. <ref>, the amplitudes of oscillation in the displacement and orientation angle become smaller as the particle moves. When the density ratio grows to ρ^*=1.5, the motion of the ellipse particle in Fig. <ref> calculated by using a finer grid of 2a=40Δ x indicates that the ellipse approaches to a fluttering motion instead of a steady descent motion.








§ CONCLUSIONS

This paper determines the jump conditions for the distribution functions in the framework of the lattice Boltzmann equation and uses these conditions to develop an immersed interface-lattice Boltzmann method (II-LBM). The derived jump conditions of the distribution function can recover the pressure discontinuity induced by an immersed interface in a fluid described by the Navier-Stokes equations. Comparisons between the direct forcing IB-LBM and the II-LBM indicate that the II-LBM has higher order of accuracy for some problems and much lower errors in the velocity. Moreover, the motion of an elastic interface shows that the II-LBM gives substantially improved volume conservation than the standard IB-LBM with both the simple external force <cit.> and Guo's external force term <cit.>. The II-LBM is also applied to solve flow around a cylinder with fixed position, a circular cylinder undergoing prescribed motion in a resting fluid, vortex-induced vibrations of a circular in one degree of freedom and sedimentation of an ellipse particle. Overall, good agreement with published results is achieved for these validation cases.

As the derivation of the jump conditions of the distribution function is given in three dimensions, it would be straightforward to apply the present method to three dimensional problems.




§ ACKNOWLEDGEMENTS

J. Q. would like to thank Nanjing University of Science and Technology for supporting the living expenses during his stay in the USA. B. E. G. acknowledges research support from NSF Awards DMS 1664645, CBET 175193, OAC 1450327, OAC 1652541, and OAC 1931516. The authors thank Jae Ho Lee for his suggestions and comments.




