Fluid Structure Interaction modelling using Lattice Boltzmann and Volume Penalization method

TL;DR

This paper introduces a novel fluid-structure interaction modeling approach combining Volume Penalization and Lattice Boltzmann methods, enabling efficient computation of rigid body dynamics without complex interface tracking.

Contribution

The paper presents a new VP-LBM method that simplifies FSI simulations by avoiding interface tracking, using GPU acceleration for efficient computation.

Findings

Accurately computes drag and lift coefficients for various FSI cases

Achieves results in good agreement with existing literature

Demonstrates efficiency with GPU implementation

Abstract

This paper proposes an approach combining the Volume Penalization (VP) and the the Lattice Boltzmann method (LBM) to compute fluid structure interaction involving rigid bodies. The method consists in adding a force term in the LBM formulation, and thus considers the rigid body similar to particular porous media. Using a characteristic function for the solid domain avoids the expensive tracking of the fluid-solid interface employed commonly in LBM to treat FSI problems. The method is applied to three FSI problems and solved using a Graphics Processor Units (GPU) device. The applications focus on the capacity of the method to compute the drag and lift coefficients for various cases : the imposed displacement of a cylinder, the particle sedimentation at a very low Reynolds number and the VIV of a cylinder in a transverse fluid flow. For all cases the VP-LBM approach yields results which are in a good agreement with those of literature.

Figures

Tables

Table 1. Table 1: Average drag coefficient obtained with the VP-LBM and Code _ _ \_ Saturne

•	VP-LBM	Code_Saturne
$\overline{C_D}$	1.653	1.577

Equations

$$\begin{array}[]{ll}\displaystyle\nabla\cdot\mathbf{u}=0\\ \displaystyle\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}=\displaystyle-\frac{1}{\rho}\nabla p+\nu\nabla^{2}\mathbf{u}-\dfrac{\chi_{\Omega_{s}}}{\eta}\left(\mathbf{u}-\mathbf{u_{s}}\right)\\ \end{array}$$

$$\chi_{\Omega_{s}}\left(\boldsymbol{x},t\right)=\left\{\begin{array}[]{ll}1\;\;\textrm{if}\;\;\boldsymbol{x}\in\Omega_{s}\left(t\right)\\ 0\;\;\textrm{otherwise}\end{array}\right.;\hskip 28.45274pt\eta\ll 1\;\;\textrm{penalization factor}$$

$$\frac{\partial f}{\partial t}+\mathbf{c}\cdot {\nabla }_{x}f=\Omega \left(f\right)$$

$$f_{\alpha }\left(\mathbf{x}+{\mathbf{c}}_{\alpha }△t,t+△t\right)-f_{\alpha }\left(\mathbf{x},t\right)={\Omega }_{\alpha }\left(f\right)+△t{F}_{\alpha }$$

$${\Omega }_{\alpha }\left(f\right)=-\frac{1}{\tau }\left(f_{\alpha }\left(\mathbf{x},t\right)-f_{\alpha }^{\text{eq}}\left(\mathbf{x},t\right)\right)$$

$$\nu =c_s^2△t\left(\tau -\frac{1}{2}\right)$$

$${\Omega }_{\alpha }\left(f\right)=-\sum _{\beta }{\left(M^{-1}SM\right)}_{\alpha \beta }\left(f_{\beta }\left(\mathbf{x},t\right)-f_{\beta }^{\text{eq}}\left(\mathbf{x},t\right)\right)$$

$$\begin{gathered} |f_{\alpha }(\mathbf{x}+{\mathbf{c}}_{\alpha }\Delta t,t+\Delta t)⟩-|f_{\alpha }(\mathbf{x},t)⟩= \\ -M^{-1}\left(S\left(|m_{\alpha }\left(\mathbf{x},t\right)⟩-|m_{\alpha }^{eq}\left(\mathbf{x},t\right)⟩\right)-\left(I-\frac{S}{2}\right)M\Delta t|F_{\alpha }\left(\mathbf{x},t\right)⟩\right) \end{gathered}$$

$$|m_{\alpha }(\mathbf{x},t)⟩=M|f_{\alpha }(\mathbf{x},t)⟩\Rightarrow |f_{\alpha }(\mathbf{x},t)⟩=M^{-1}|m_{\alpha }(\mathbf{x},t)⟩$$

$$\mathbf{c_{\alpha}}=\left\{\begin{array}[]{lll}\left(0,0\right)&\alpha=0\\ \left(\textrm{cos}\left(\left(\alpha-1\right)\dfrac{\pi}{2}\right),\textrm{sin}\left(\left(\alpha-1\right)\dfrac{\pi}{2}\right)\right)c&\alpha=1,2,3,4\\ \left(\textrm{cos}\left(\left(2\alpha-9\right)\dfrac{\pi}{4}\right),\textrm{sin}\left(\left(2\alpha-9\right)\dfrac{\pi}{4}\right)\right)\sqrt{2}c&\alpha=5,6,7,8\end{array}\right.$$

$$f_{\alpha }^{\text{eq}}={\omega }_{\alpha }\rho \left(1+\frac{{\mathbf{c}}_{\alpha }\cdot \mathbf{u}}{c_s^2}+\frac{\mathbf{uu}:\left({\mathbf{c}}_{\alpha }{\mathbf{c}}_{\alpha }-c_s^{2}I\right)}{2c_s^4}\right),$$

$$F_{\alpha }=\left(1-\frac{1}{2\tau }\right){\omega }_{\alpha }\left(\frac{{\mathbf{c}}_{\alpha }-\mathbf{u}}{c_s^2}+\frac{{\mathbf{c}}_{\alpha }\cdot \mathbf{u}}{c_s^4}{\mathbf{c}}_{\alpha }\right)\cdot \rho \mathbf{F}$$

$$\omega_{\alpha}=\left\{\begin{array}[]{lll}4/9,&&\alpha=0\\ 1/9,&&\alpha=1,2,3,4\\ 1/36&&\alpha=5,6,7,8\\ \end{array}\right.,$$

$$\begin{array}[]{ll}|m_{\alpha}^{eq}\rangle&=\displaystyle\left(\rho,e^{eq},\varepsilon^{eq},j^{eq}_{x},q^{eq}_{x},j^{eq}_{y},q^{eq}_{y},p^{eq}_{xx},p^{eq}_{xy}\right)\\ &=\displaystyle\rho\left(1,-2+3(u^{2}_{x}+u^{2}_{y}),1-3(u^{2}_{x}+u^{2}_{y}),u_{x},-u_{x},u_{y},-u_{y},u^{2}_{x}-u^{2}_{y},u_{x}u_{y}\right)^{T}\end{array}$$

$$\begin{split}\begin{array}[]{cccccccccccccc}\displaystyle M=\begin{pmatrix}1&\hskip 5.69046pt1\hskip 5.69046pt&1\hskip 5.69046pt&1\hskip 5.69046pt&\hskip 5.69046pt1&\hskip 5.69046pt1&\hskip 5.69046pt1&\hskip 5.69046pt1&\hskip 5.69046pt1\\ -4&-1&-1&-1&-1&2&2&2&2\\ 4&-2&-2&-2&-2&1&1&1&1\\ 0&1&0&-1&0&1&-1&-1&1\\ 0&-2&0&2&0&1&-1&-1&1\\ 0&0&1&0&-1&1&1&-1&-1\\ 0&0&-2&0&2&1&1&-1&-1\\ 0&1&-1&1&-1&0&0&0&0\\ 0&0&0&0&0&1&-1&1&-1\end{pmatrix}\end{array},\end{split}$$

$$S=diag(s_0,s_1,s_2,s_3,s_4,s_5,s_6,s_7,s_8)$$

$$|{\tilde{f}}_{\alpha }⟩=M^{-1}\left(|m_{\alpha }⟩-S\left(|m_{\alpha }⟩-|m_{\alpha }^{eq}⟩\right)-\left(I-\frac{S}{2}\right)M\Delta t|F_{\alpha }⟩\right)$$

$$|f_{\alpha }\left(\mathbf{x}+{\mathbf{c}}_{\alpha }△t,t+△t\right)⟩=|{\tilde{f}}_{\alpha }\left(\mathbf{x},t\right)⟩$$

$$\rho =\sum _{\alpha }{f}_{\alpha }\rho \mathbf{u}=\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }+\frac{△t}{2}\rho \mathbf{F}$$

$$\rho \mathbf{u}=\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }-\frac{△t}{2}\rho \frac{{\chi }_{{\Omega }_s}}{\eta }\left(\mathbf{u}-{\mathbf{u}}_{\mathbf{s}}\right)$$

$$\rho \left(1+\frac{△t}{2}\frac{{\chi }_{{\Omega }_s}}{\eta }\right)\mathbf{u}=\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }+\frac{△t}{2}\rho \frac{{\chi }_{{\Omega }_s}}{\eta }{\mathbf{u}}_{\mathbf{s}}$$

$$\mathbf{u}=\frac{\sum _{\alpha }{\mathbf{c}}_{\alpha }{f}_{\alpha }+\frac{△t}{2}\frac{{\chi }_{{\Omega }_s}}{\eta }\rho {\mathbf{u}}_{\mathbf{s}}}{\rho +\frac{△t}{2}\frac{{\chi }_{{\Omega }_s}}{\eta }\rho }$$

$$m\frac{d^2{\mathbf{x}}_{\mathbf{G}}}{d{t}^2}+c\frac{d{\mathbf{x}}_{\mathbf{G}}}{dt}+k\left({\mathbf{x}}_{\mathbf{G}}-{\mathbf{x}}_0\right)={\mathcal{F}}_{\mathbf{f}}+{\mathcal{F}}_{\mathbf{ext}}$$

$$I\frac{d^2\mathbf{\theta }}{d{t}^2}+c_{\theta }\frac{d\mathbf{\theta }}{dt}+k_{\theta }\left(\mathbf{\theta }-{\mathbf{\theta }}_0\right)={\mathcal{T}}_{\mathbf{f}}+{\mathcal{T}}_{\mathbf{ext}}$$

$$\mathbf{F}\left({\mathbf{x}}_{\mathbf{\Gamma }}\right)=\left({\mathbf{e}}_{\alpha }-{\mathbf{u}}_{\mathbf{\Gamma }}\right){\tilde{f}}_{\alpha }\left({\mathbf{x}}_{\mathbf{f}}\right)-\left({\mathbf{e}}_{\overline{\alpha }}-{\mathbf{u}}_{\mathbf{\Gamma }}\right){\tilde{f}}_{\overline{\alpha }}\left({\mathbf{x}}_{\mathbf{s}}\right)$$

$${\mathcal{F}}_{\mathbf{f}}=\sum \mathbf{F}\left({\mathbf{x}}_{\mathbf{\Gamma }}\right)$$

$${\mathcal{T}}_{\mathbf{f}}=\sum \left({\mathbf{x}}_{\mathbf{\Gamma }}-{\mathbf{x}}_{\mathbf{G}}\right)\times \mathbf{F}\left({\mathbf{x}}_{\mathbf{\Gamma }}\right)$$

$$S=diag\left(\frac{1}{\tau },1.1,1.25,\frac{1}{\tau },1.8,\frac{1}{\tau },1.8,\frac{1}{\tau },\frac{1}{\tau }\right)$$

$$y_G\left(t\right)=A+B\text{cos}\left(\omega t\right)=\left(y_G\left(0\right)-\frac{D}{4}\right)+\frac{D}{4}\text{cos}\left(\omega t\right)$$

$${\omega }^{\star }=\omega \frac{D}{U_0}{\text{B}}^{\star }=\frac{\text{B}}{D}{\mathbf{u}}^{\star }=\frac{\mathbf{u}}{U_0}$$

Full text

Keyword

Lattice Boltzmann Method, Fluid Structure Interaction, Volume Penalization, Vortex Induced Vibration (VIV)

1 Introduction

In the following section, the mathematical background is presented. This part deals with the Lattice Boltzmann Method and more particularly the Multiple Relaxation Time (MRT) approach, the Volume Penalization and how these two methods are combined. The last section presents computational applications computed on a GPU device. For the first tested case, the imposed displacement of a cylinder in a transverse fluid flow at a Reynolds of $100$, is not a real case of FSI. This application enables to validate the capacity of the method to compute drag and lift forces. The second example deals with particle sedimentation under gravity, the particle is free to move under fluid and gravity constraints and complex trajectories can be obtained. The last application focuses on the behavior of the VP-LBM method for FSI problems when a parameter is varied. The stiffness parameter of a free oscillating cylinder in a transverse fluid flow is changed, and the capacity of the VP-LBM to reproduce results from the literature is evaluated.

2 Governing equations

In this section, the numerical models are exposed. The following notations are used : $\rho$ and $\mathbf{u}$ are the macroscopic density and velocity, and bold characters denote vectors.

2.1 Volume penalization

Let us consider a fluid domain $\Omega_{f}$, a solid domain $\Omega_{s}$, $\Gamma$ the fluid-solid interface, and let us note $\Omega=\Omega_{f}\cup\Omega_{s}\cup\Gamma$. The Volume Penalization (VP) method consists in extending the Navier-Stokes equations on the whole domain $\Omega$, and considering the solid domain as a porous medium with a very small permeability. The method was introduced by Angot et al. [11] and already applied to macroscopic equations for moving bodies [14]. The small permeability of the solid domain is modelled using a penalization coefficient, hence the desired boundary conditions at the fluid-solid interface are naturally imposed. With this method, the incompressible Navier-Stokes equations are written as follows :

[TABLE]

where

[TABLE]

$\mathbf{u}$ denotes the velocity field, $p$ is the pressure field, $\rho$ and $\nu$ are the density and the viscosity of the fluid. $\mathbf{F}=\dfrac{\chi_{\Omega_{s}}}{\eta}\left(\mathbf{u}-\mathbf{u_{s}}\right)$ is the penalization term, and $\mathbf{u_{s}}$ is the velocity field in the solid domain.

2.2 Lattice Boltzmann method

Based on the Boltzmann equation (equation (3)) proposed in the context of the Kinetic Gaz Theory by L. Boltzmann in 1870, the Lattice Boltzmann Method has been successfully used to model fluid flow since the 90’s

[TABLE]

This equation models the transport of $f\left(\boldsymbol{x},t,\mathbf{c}\right)$, a probability density function to find a particle at location $\boldsymbol{x}$ and time $t$ with the velocity $\mathbf{c}$, $\Omega\left(f\right)$ being a collision operator. The link between the Boltzmann equation and the Navier-Stokes equations is well-known since the Chapmann-Enskog expansion proposed in 1915.

The Lattice Boltzmann method considers the discretization of equation (3) according to space and velocity and leads to the following discretized equations :

[TABLE]

where $\displaystyle f_{\alpha}\left(\boldsymbol{x},t\right)=\displaystyle f\left(\boldsymbol{x},\mathbf{c_{\alpha}},t\right)$, $\displaystyle F_{\alpha}$ is a forcing term related to the discrete velocity $\mathbf{c_{\alpha}}$. The first model proposed by Bhatnagar et al. [15] is the BGK model which is based on a linear collision operator with a single relaxation time :

[TABLE]

where $f^{\textrm{eq}}$ is the equilibrium function and $\tau$ is the non dimensional relaxation time which is linked to the fluid viscosity as follows (6).

[TABLE]

We propose in this paper to use the Multiple Relaxation Time (MRT) model, introduced by d’Humières [16] for stability reasons. This scheme consists in using a transformation matrix $M$ to work with macroscopic quantities in the moment space. For the MRT model, the collision operator is :

[TABLE]

where $S$ is a diagonal relaxation matrix.

The transformation matrix $M$ enables to express the moments $m_{\alpha}(\boldsymbol{x},t)$ according to the distribution functions $f_{\alpha}(\boldsymbol{x},t)$. The lattice Boltzmann equation with a source term, as given by Lu et al. [17], becomes:

[TABLE]

In this equation, $|\bullet\rangle$ denotes a column vector and $I$ is the identity matrix. The moments $|m_{\alpha}\left(\boldsymbol{x},t\right)\rangle=\left(m_{0},m_{1},\ldots,m_{8}\right)^{T}$ are deduced from:

[TABLE]

The model used most commonly to simulate two-dimensional flows is the nine-velocity square lattice model D2Q9 (figure 1) [18].

[TABLE]

Where $c=\dfrac{\triangle x}{\triangle t}$. Usually $\triangle x=\triangle y=\triangle t=1$ are chosen. In addition, the equilibrium distribution function for the D2Q9 model is :

[TABLE]

and the forcing term $\displaystyle F_{\alpha}$ in equation (4) is [19] :

[TABLE]

where $\omega_{\alpha}$ are the weighting coefficients commonly used for the D2Q9 model :

[TABLE]

$c_{s}$ is the speed of sound, which for the D2Q9 model is $c_{s}=\dfrac{c}{\sqrt{3}}$.

For the D2Q9 model the corresponding equilibria in the moment space are given by :

[TABLE]

where $u_{x}$ (respectively $u_{y}$) is the horizontal (respectively vertical) component of velocity $\mathbf{u}$.

The transformation matrix $M$ is :

[TABLE]

and $S$ is a diagonal matrix that contains the relaxation rates of each moment:

[TABLE]

At each time-step, the LBM algorithm consists in solving first a collision step :

[TABLE]

and next a streaming step :

[TABLE]

Finally, the macroscopic quantities are computed according to the following expressions :

[TABLE]

In the present approach, the volume penalization term is added :

[TABLE]

To avoid instabilities, the term including $\mathbf{u}$ in the penalization force is moved to the left hand side of equation (20)

[TABLE]

This leads to the modified update step to compute the macroscopic velocity field :

[TABLE]

In the fluid domain, where $\chi_{\Omega_{s}}=0$ the classical LBM equation is obtained whereas in the solid domain, where $\chi_{\Omega_{s}}=1$, equation (22) forces the velocity field to approach $\mathbf{u_{s}}$.

2.3 Structure displacement

In this study, only structures which can be modelled as rigid bodies coupled with springs and dampers are considered. The center of gravity moves according to the following equation :

[TABLE]

where $\boldsymbol{x_{G}}$ are the coordinates of the center of gravity of the solid body, $\boldsymbol{x_{0}}$ is the equilibrium position, $m$ is the mass, $c$ and $k$ are the damping and stiffness coefficients, $\boldsymbol{\mathcal{F}_{f}}$ are the fluid forces and $\boldsymbol{\mathcal{F}_{ext}}$ are external forces (gravity for example). The rotation of the solid is solved using the following equation :

[TABLE]

where $\boldsymbol{\theta}$ is the rotation vector, $I$ is the inertia moment of the body, $c_{\theta}$ and $k_{\theta}$ are the damping and stiffness coefficients for rotation, $\boldsymbol{\mathcal{T}_{f}}$ is the torque induced by the fluid and $\boldsymbol{\mathcal{T}_{ext}}$ is an external torque.

The fluid forces are computed using the momentum exchange method (MEM) proposed by Wen at al. [20]. We note $\mathbf{x_{f}}$ a boundary node in the fluid domain and $\mathbf{x_{s}}$ the image of this boundary node through the solid interface by a lattice velocity $\mathbf{e_{\alpha}}$, also called incoming velocity( cf. figure 2). The intersection point between the fluid-solid interface and the link $\mathbf{x_{f}}-\mathbf{x_{s}}$ is $\displaystyle\mathbf{\displaystyle x_{\Gamma}}$, and the outgoing lattice velocity is denoted $\displaystyle\mathbf{\displaystyle e_{\overline{\alpha}}}=-\mathbf{e_{\alpha}}$.

The local force at $\mathbf{\displaystyle x_{\Gamma}}$ is computed using the following expression :

[TABLE]

and the total fluid force acting on the solid domain is :

[TABLE]

The torque is obtained with

[TABLE]

Giovacchini and Ortiz [21] showed that the MEM does not depend on the way the boundary conditions at the solid domain are implemented.

For each time step the fluid-structure problem is solved according to algorithm 1.

In a previous work, Benamour et al. [12] showed that the volume penalisation method combined with the LBM gives good results for computing flows around fixed bodies. In the following section, the VP-LBM is applied to different cases of fluid flows around moving bodies.

3 Numerical applications

The VP-LBM is applied to three cases. The first one considers the imposed transverse displacement of a cylinder in a fluid flow at a Reynolds number of $100$. The lift (and drag) coefficients computed with the VP-LBM are compared with those obtained using a classical CFD code (Code$\_$Saturne, [22]). The second case is the study of particle sedimentation and the last one focuses on vortex-induced vibration (VIV) of a circular cylinder. All computations were run on a NVIDIA QUADRO P500 GPU card, using a CUDA implementation. For all computations, the following relaxation rates were chosen :

[TABLE]

where the relaxation time $\tau$ is related to the fluid viscosity thanks to the equation (6), and a value of penalization factor $\eta=10^{-6}$ was selected.

In the remaining of this paper l.u. will refer to lattice length units and t.s. to lattice time units

3.1 Imposed displacement of a cylinder in a transverse fluid flow

To validate the VP-LBM, the imposed displacement of a cylinder in a transverse flow at a Reynolds number of $100$ was first simulated (figure (3)). For that case, a constant velocity profile was imposed at the inlet using the classical half-way Bounce-Back method, and the outflow boundary condition at outlet was modelled using the convective condition [23]. Symmetry boundary conditions ( $\mathbf{u}\cdot\mathbf{n}$) were imposed at the other boundaries, in order to apply the same boundary conditions as those employed with the finite volume code (Code$\_$Saturne).

The computational procedure is the following : first, a computation was run without solid displacement to obtain a well-established fluid flow. This state is considered as the initial time $\left(t=0\right)$. Next, the motion of the cylinder was imposed according to the following expression :

[TABLE]

where $y_{G}$ is the $\mathbf{y}$ coordinate of the center of gravity, and $D$ the cylinder diameter, and the fluid flow was calculated.

To test the ability of the VP-LBM to calculate flows around moving obstacles, and to validate the computation of the hydrodynamic force (equation (25)) exerted on the moving body, the results are compared with those obtained by a conventional computational code, the finite volume software Code$\_$Saturne [22]. To compare the results obtained with the VP-LBM, and with the finite volume method, the same non dimensional numbers were used :

[TABLE]

In this study: $\omega^{\star}=1.55$, and the LBM parameters (in lattice units) are

[TABLE]

For the computation performed with Code$\_$Saturne, a non uniform grid that was very fine in the vicinity of the cylinder, and coarser in the remaining of the fluid domain, was built (see figure 4(a)). Such a mesh enabled a decrease in computational time. For the LBM computations, a regular grid was used.

Figure 5 shows the norms of the non-dimensional velocity field computed with Code$\_$Saturne software and with the VP-LBM at the same non-dimensional time $t^{\star}=64$, where $t^{\star}=t\dfrac{U_{0}}{D}$. It can bee seen that both computations yield the same velocity field. More particularly, the vortex which appears behind the moving cylinder is located at the same place.

In figures 6, the lift $\left(C_{L}\right)$ and drag $\left(C_{D}\right)$ coefficients are compared :

[TABLE]

In this figure, it can be noticed that spurious oscillations occur when LBM is used. However, when using LBM, a cartesian grid is commonly used, and this behavior has been highlighted for computations of FSI problems carried out with LBM, combined with any method chosen for simulating flows around curved boundaries (Bounce-back or immersed boundary methods) ([20]). A good agreement between both methods was obtained for the lift coefficient, whereas small differences are observed for the drag coefficient. These differences are not discriminating for the VP-LBM approach, because as shown in table 1, the difference based on the average drag is less than $5\%$, and this gap is in the range of what can be expected when two different numerical models are compared.

To conclude with this first example, the VP-LBM predicts accurately the fluid forces exerted on a solid whose motion is imposed. In order to test the validity of the VP-LB, the following examples deal with real cases of fluid structure interaction.

Remark

In order to reduce the size of the LBM problem ($L_{1}\times L_{2}$) and thus the computational time, the value of the parameter $\tau=0.56$ was chosen very close to the stability limit $0.5$. An inlet velocity of $0.048780$ ensured a small Mach number suitable for the LBM approach. For that simulation, $41$ lattice units were used for the cylinder diameter, which is a small value for computing a flow of Reynolds number $100$ with LBM. This can explain the spurious oscillations.

3.2 Sedimentation of a particle under gravity

The next case focuses on the sedimentation of a particle under gravity in an infinite channel (figure 7) for centered and non-centered configurations. These problems have been extensively used for model validations and are very useful to test the capacity of a method to capture complex trajectories [24, 25, 20].

A circular particle of diameter $D$ falls under gravity $\boldsymbol{g}$ in a fluid of density $\rho$ in a vertical channel of width $H$. At the initial state, the particle is at a distance $x_{0}$ from the left wall, a distance $y_{0}$ from the top of the channel and the velocity of the particle is equal to zero. For that case, the particle displacement can be described using equations (23 ) and (24) modified as follows :

[TABLE]

where $\rho$ denotes the fluid density and $\rho_{s}$ the solid density. The last term in equation (31) represents the weight and the buoyancy (Archimedes’ principle) acting on the particle.

For small Reynolds numbers and a large non dimensional width $\tilde{H}=\dfrac{H}{D}$, the particle reaches a steady state, where the drag force can be approximated according to the following expression [26]:

[TABLE]

where $\kappa$ denotes the correction factor which represents the channel confinement effect,

[TABLE]

$\mu$ is the fluid viscosity and $\mathbf{v_{s}}$ is the gravitational settling of the particle :

[TABLE]

3.2.1 Centered particle : $x_{0}=\dfrac{H}{2}$

We consider the same physical properties as those chosen in previous works [24, 27]. The diameter of the particle is $D=0.24\;\textrm{cm}$, the fluid density is $\rho=1\;\textrm{g}\cdot\textrm{cm}^{-3}$, the fluid viscosity is $\mu=0.1\;\textrm{g}\cdot\textrm{cm}^{-1}\cdot\textrm{s}^{-1}$, the gravitational acceleration is $\|\boldsymbol{g}\|=980\;\textrm{cm}\cdot\textrm{s}^{-2}$ and the non dimensional width of the channel is $\tilde{H}=5$.

The problem was modelled with a regular $120\times 1200$ lattice grid, $24$ nodes across the diameter of the cylinder ($D_{LBM}=24\;\textrm{l.u.}$), and a relaxation time $\tau=0.8$. Initially, the particle is located at the center of the channel : $x_{0}=0.5H$ and $y_{0}=0.5L$. No-slip boundary conditions were imposed on the left and right walls. A zero velocity boundary condition was applied at the inlet (top of the channel) and free stream conditions were applied at the outlet (bottom). A large value of $L$ was chosen, so that the inlet and outlet did not influence the behavior of the particle.

Computations were performed for various mass ratios $\rho_{r}=\dfrac{\rho_{s}}{\rho}$ taken from the following list :

[TABLE]

For all cases the particle velocities and the drag coefficients are plotted in figure 8 and compared with the analytical solution (equations (33) and (35)).

For $\rho_{r}=\left\{0.98;0.99;1.01;1.02\right\}$, results match with the analytical solution. For larger ratios ($\rho_{r}=\left\{0.95;1.05\right\}$) differences are observed. As mentioned in previous studies [25, 27], the reason could be that the analytical solutions are available only for small mass numbers $\rho_{r}$ and not for larger ratios such as $1.05$ or $0.95$, which explains the gap between the analytical solution and the results obtained with the VP-LBM approach. As noticed in the previous example, small fluctuations on the drag coefficient are observed when the velocity increases $\left(\rho_{r}=0.95\;\textrm{or}\;1.05\right)$, but the averages fit with the analytical solutions.

3.2.2 Non centered particle : $x_{0}=0.75D$

The following case deals with a particle whose initial position is not at the center of the channel. The fluid properties are $\rho=1\;\textrm{g}\cdot\textrm{cm}^{-3}$, and $\mu=0.1\;\textrm{g}\cdot\textrm{cm}^{-1}\cdot\textrm{s}^{-1}$ and the physical problem concerns a particle of diameter $D=0.1\;\textrm{cm}$, a mass ratio $\rho_{r}=1.03$ and $\|\boldsymbol{g}\|=980\;\textrm{cm}\cdot\textrm{s}^{-2}$. The Reynolds number based on the final velocity of the particle is $Re=8.33$.

For the LBM computations a $125\times 1550$ lattice grid was used, the cylinder diameter was $31\;\textrm{l.u.}$, the non dimensional relaxation time was $\tau=0.6$ and the particle was released at $y_{0}=12.5D$. The boundary conditions were the same as in the previous case.

Results are plotted in figures 9. $\mathbf{v_{g}}$ denotes the vector velocity of the gravity center of the particle, and $\omega_{g}$ its rotational velocity. The time-dependent position, horizontal, vertical and rotational velocities are compared with results from Tao et al. [25] and good agreements are found between them. Small spurious oscillations are observed for the rotational velocity (figure 9(d)), but these oscillations are also observed by Tao et al. [25].

Figures 10 and 11 show the fluid velocity and the vorticity field around the particle at four different times. The dynamics of the flow field and the particle can be analyzed using the velocity magnitude and the vorticity. The particle goes first to the right and rotates in a positive direction. Next a brief oscillation occurs around the central line of the channel and finally the particle stays in the middle of the channel with a steady velocity.

This example shows that the VP-LBM method is able to predict a complex trajectory for a real case of fluid structure interaction at a very low Reynolds number.

3.3 Vortex Induced Vibration of a cylinder

Let us consider the case depicted in figure 3. In this paragraph the displacement of the cylinder is not imposed, but driven by the fluid forces. The displacement is let free according to the $\mathbf{y}$ axis. Initially, the cylinder is at rest, the force exerted by the spring on the cylinder is equal to [math]. Due to the vortex shedding, the fluid applies a force according to the vertical axis and the cylinder begins to oscillate. The rigid displacement of the cylinder is solved using equation (23) projected onto the $\mathbf{y}$ axis, with $c=0$ and $\boldsymbol{\mathcal{F}_{ext}}=0$ :

[TABLE]

where $y_{G}$ is the position of the center of gravity of the cylinder according to the $\mathbf{y}$ axis, $y_{0}$ is the position at rest, and $\mathcal{F}_{y}$ denotes the fluid forces acting on the cylinder in the vertical direction. Equation (24) is not used.

To simplify the analysis, the non dimensional form of equation (37) is considered :

[TABLE]

with the non dimensional numbers :

[TABLE]

To scale the results, the effective stiffness introduced by Shiels et al. [28] is used :

[TABLE]

where $\omega^{\ast}$ is computed from the analysis of the cylinder displacement. The effective stiffness enables to use only one plot to represent the results, because it combines the reduced mass $m^{\ast}$ and the reduced stiffness $k^{\ast}$. For a given Reynolds number, $k^{\ast}_{\textrm{eff}}$ completely determines the system.

The computational parameters are identical to those related to the first test case (section 3.1)

[TABLE]

$29$ values of the stiffness parameter $k^{\ast}$ were considered, from $1$ to $30$. For each $k^{\ast}$, the choice of $\omega^{\ast}$ was based on the the analysis of the drag coefficient, and with these values, $k^{\ast}_{\textrm{eff}}$ was computed.

In figure 12 the maximum of the non dimensional amplitude $y^{\ast}_{max}-y^{\ast}_{0}$ versus $k^{\ast}_{\textrm{eff}}$, obtained obtained by the VP-LBM methods and compared with data from [28] is plotted. It can be seen that the proposed method is able to reproduce the lock-in phenomena when $k^{\ast}$ is close to $m{\omega^{\ast}}^{2}$, and when important displacement of the cylinder occurs (more than $50\%$ of the diameter). Figure 13 depicts the temporal evolution of the non-dimensional amplitude for $3$ values of $k^{\ast}$. The velocity, pressure and vorticity fields are presented in figures 14-16. Different behaviors according to the value of $k^{\ast}$ were obtained. For high values of $k^{\ast}$ the displacement of the cylinder is small, and the fluid flow shows similar features as seen for a fixed cylinder. However, the behavior of the solid is more chaotic, several frequencies can be noticed when examining the displacement of the body (figure 13(c)). For $1\leq k^{\ast}_{eff}\leq 5$ the flow pattern changes drastically. The vortex shedding frequency increases and the velocity decreases in the wake of the cylinder (figure 14(b)).

To conclude, this last application shows the capacity of the VP-LBM method to capture various physical behaviors induced by changes of the parameter $k^{\ast}$.

4 Conclusion

A combined approach coupling the Volume Penalization and the Lattice Boltzmann method for fluid structure interaction was proposed. The method consists in adding a force term, which is similar to Darcy’s law, into the Boltzmann equation. The advantage of this method is that no explicit computation of the fluid structure interface is needed, the use of a characteristic function is sufficient. In addition, the fluid forces exerted on the structure are computed using the classical momentum exchange method. The method was implemented on a GPU architecture device and tested on three cases. The first case which deals with the imposed displacement of a cylinder in a transverse fluid flow at a Reynolds of $100$, validates the capacity of the method to compute drag and lift forces. The second application focuses on the sedimentation of a particle at a very low Reynolds number in a channel. The proposed method succeeds to capture the complex trajectory of the particle, which is composed of translational and rotational components. In the last application, the vortex induced vibration of a cylinder in a transverse fluid flow is considered. The capacity of the method to predict the physics of the fluid flow and the structure behavior for different values of the stiffness parameter of the cylinder is tested. Results are in a good agreement with those of literature. All these applications validate the VP-LBM as an efficient tool to model FSI in case of rigid bodies.

5 Bibliography