Discrete models for fluid-structure interactions: the Finite Element Immersed Boundary Method

The aim of this paper is to provide a survey of the state of the art in the finite element approach to the Immersed Boundary Method (FE-IBM) which has been investigated by the authors during the last decade. In a unified setting, we present the different formulation proposed in our research and highlight the advantages of the one based on a distributed Lagrange multiplier (DLM-IBM) over the original FE-IBM.


Introduction
In this paper we present in a unified setting some results which have been the object of our research on the finite element discretization of the Immersed Boundary Method (ibm) during the last years.This is a survey paper which aims at providing a comprehensive and self-contained presentation of this subject.
The ibm has been introduced by Peskin in the seventies [32,30] for the numerical approximation of biological phenomena involving fluids and solids (typically, the blood flow in the cardiac muscle).Since then, it has been successfully adopted in many application areas; the interested reader can refer to [31,27,13,35,18,17,29,26,23,25,20] and the references therein.One of the main features of the ibm is that the solid is thought as a part of the fluid and its effect on the dynamics of the system is modeled through a Dirac delta functions that has the role of linking Eulerian and Lagrangian variables.The original discretization of the ibm, which makes use of finite differences for the underlying fluid equations, needs a suitable approximation of the delta function.This procedure is a key point of the numerical strategy: tuning the approximation of the delta function is a crucial aspect which has a great influence on how the discontinuities of the solution are well captured.
In this framework it is natural to consider a finite element version of the ibm, where the delta function does not need any approximation, since it can be handled directly by the variational formulation.The first steps in this direction have been presented in [9].We refer to this formulation as fe-ibm [10,11,21,8,6].We recall that our formulation can model systems where the solid and the fluid have the same dimension (i.e., codimension zero case), or where the solid has codimension one.There are also other finite element approaches to the ibm: some of them handle the delta function variationally as we do [22,14,24], other ones use different techniques based on its approximation [33,34,28].
During our studies concerning this topic, we have discussed two important issues: the mass conservation of the overall procedure and the stability of the fully discrete scheme.Concerning the mass conservation, we have observed that, as expected, discontinuous pressure schemes can enforce the mass conservation (expressed by the fact that the fluid velocity is divergence free) much more locally than continuous ones.In this respect, we have analyzed a modification of continuous pressure schemes which improves their performances [3,4].Several numerical tests confirm our theoretical investigations [2].In our scheme, we use a semi-implicit time advancing scheme.For its stability, we have observed that the semi-implicit fe-ibm is subjected to a cfl condition which depends on the parameters involved with the model; in particular, it has been shown that the method is cfl-stable even when the densities of the fluid and the solid are comparable.In this respect, the fe-ibm is superior to the more popular Arbitrary Lagrangian Eulerian (ale) method which has been proved to be unconditionally unstable unless a fully implicit approach is considered [12].
Recently, we have studied a modification of the fe-ibm consisting in the introduction of a Lagrange multiplier associated to the equation coupling solid and fluid velocities.Since the resulting formulation resembles the fictitious domain method with distributed Lagrange multiplier, we refer to the new scheme as dlmibm (see [19,15,16]).First promising numerical results have been presented in [5], showing that the (semi-implicit) scheme is unconditionally stable and (surprisingly) enjoys better conservation properties than the previous one.In [7] we performed a stability analysis for the dlm-ibm, showing rigorously that the semi-implicit scheme is unconditionally stable, and reporting on some additional numerical tests.
The structure of the paper is as follows: in Section 2 we recall the fundamental properties of the original ibm; in Section 3 we describe the fe-ibm and the variational treatment of the Dirac delta function; in Section 4 we show how to introduce the Lagrange multiplier which has the effect of stabilizing the time advancing scheme; in Section 5 we discuss the stability of the proposed scheme and in Section 6 we report some numerical experiments.

The Immersed Boundary Method
In this section we review the formulation of the ibm in its original version introduced by Peskin in [31] for immersed materials modeled as collections of fibers and further extended to cover the case of thick bodies modeled as a hyperelastic material in [8].One of the main difficulties to face when treating fluid-structure interaction problems consists in the fact that the fluid is naturally modeled using Eulerian variables, while for the solid the Lagrangian framework is more appropriate.The ibm is a way to mix Eulerian and Lagrangian variables thanks to the use of the Dirac delta function.The main idea in ibm is to consider the structure as a part of the fluid where additional mass and forces are concentrated.
Let Ω ⊆ R d , d = 2, 3 be a region containing both a viscous incompressible fluid and an immersed elastic structure.We introduce the Navier-Stokes equations describing the dynamics of the fluid with respect to the Eulerian variable denoted by x: (1) Here ρ and ν are positive constants denoting the density and the viscosity of the fluid.The unknowns u and p represent the velocity and the pressure of the fluid, respectively.The right hand side F stands for the forces acting on the fluid, and, in absence of external volume forces, it takes into account the so called fluid-structure interaction forces, that is the forces exerted by the elastic structure on the fluid.The immersed structure is considered as an elastic incompressible material filling at time t a region B t ⊆ Ω of codimension one or zero.Using Lagrangian variables, B t can be represented as the image of a mapping X from a reference domain B ⊂ R m , with m = d or m = d − 1.We denote by s the Lagrangian coordinates in the reference domain B, then X(s, t) represents the position of a point in the current solid domain B t which is labeled s in the reference domain, that is: Since u(x, t) represents the velocity of a material point at position x at time t, we have the following condition (2) u(x, t) = ∂X ∂t (s, t) for x = X(s, t).
The source term in the first equation of (1) can now be written in absence of external volume forces as follows: (3) where f (s, t) is the force density that the immersed material exerts to the fluid and that can be modeled in different ways depending on the application field.
In order to give an idea, we report here the most simple example of an elastic structure represented by a massless closed curve Γ t immersed in the fluid occupying a two dimensional domain Ω.The curve is given in parametric form as X(s, t) for 0 ≤ s ≤ L, with X(0, t) = X(L, t).The local density force applied by the curve to the fluid is given by f = ∂(T τ ) ∂s, where T is the tension and τ is the unit tangent to the curve.Assuming that T depends linearly on |∂X/∂s| we obtain To summarize, the ibm formulation for fluid-structure interaction problems has the following form.
in Ω (10) The above formulation has been derived in [31] by using the principle of least action and it is well suited to the case of a structure described by a system of elastic fibers.In order to treat more general elasticity models for thick structure, it has been observed in [8] that an additional transmission condition along the interface between the immersed body and the surrounding fluid is needed.We give here the formulation of ibm obtained in [6] extending the formulation of [8] to the case of fluid and solid with different densities.
The strong form of the equation of motion can be written as follows (12) where σ stands for the Cauchy stress tensor.Let us consider first the case m = d and let us introduce some assumption describing the characteristics of fluid and solid materials.In the fluid, σ is modeled by means of the Navier-Stokes stress tensor σ f = −pI + ν(∇ u + (∇ u) ).We assume that the solid is composed by a viscous hyperelastic material, therefore σ can be decomposed as the sum of the viscous part σ f and an elastic part σ s , which takes into account the elastic behavior of the material.Hence, the Cauchy stress tensor can be written as follows Since in the description of the deformation of an elastic body the Lagrangian setting is more convenient, we express σ s in Lagrangian variables: this can be done by introducing the first Piola-Kirchhoff stress tensor P defined in such a way that, for any arbitrary smooth portion P of B evolving as where N is the outer normal to the region P in the Lagrangian coordinates.The first Piola-Kirchhoff stress tensor gives the elastic force per unit reference volume (d = 3) or area (d = 2), expressed in the reference space, and its pointwise expression is given by Moreover, the densities of the fluid and the solid could be different so that we set The case of the immersed body occupying a region of codimension one represents a mathematical simplification of a thin body with thickness t s very small with respect to the other space dimensions, so that one can assume that the physical quantities depend only on the variables along the middle section of the body represented by B t and are constant in the orthogonal direction.In this case the thickness t s appears as a multiplicative factor in the expressions of the Piola-Kirchhoff stress tensor and of the density of the solid (see [6] for the details).In order to unify the formulation of the problem we set Using the above definitions ( 14), ( 16) and ( 17) in ( 12), the principle of virtual work provides with some computations the following strong form of the problem.

The finite element Immersed Boundary Method
In this section we review the history of our approach to the finite element Immersed Boundary Method fe-ibm.
The starting point of our analysis has been introduced in [9] (see also [10]).The main idea is that the source term which represents the effects of the structure on the fluid and which involves the presence of a Dirac delta function, can be naturally written in variational form.This leads to a variational formulation of the Immersed Boundary Method which can be used efficiently for the finite element discretization.In order to describe this formulation, let us consider the previously introduced source term defined in terms of the Dirac delta function (see (3)): where the function f is related to elastic properties of the solid expressed in the Lagrangian variable s.The following Lemma shows that F can actually be interpreted as an element of H −1 (Ω).
Lemma 3.1 (see [9,10]).Let B t be Lipschitz continuous for all t ∈ [0, T ] and suppose that f belongs to L 2 (B) for all t ∈]0, T [).Then F (see (19)) is a distribution belonging to H −1 (Ω) defined as This observation made it possible to introduce the following variational formulation for the first fe-ibm approach.We start describing this model in the case of a constant density ρ and when the elastic force can be modeled as κ being the elastic constant of the solid (see ( 4)).
As described in Section 2, the first model has been later modified to include the case of a general hyperelastic material and different densities in the fluid and the solid [8,6].The variational formulation of Problem 2 reads In the last equation we have used the previously introduced definitions for the scaled density difference and the scaled Piola-Kirchhoff tensor P being the standard Piola-Kirchhoff tensor (see (17)).
We are now ready to describe a finite element formulation associated with Problem 4. We start with the space semidiscretization which is a more or less immediate consequence of our variational formulation.
Let T h be a standard triangulation of Ω.It is important to remark that this triangulation will never change during the computation.This is one of the main differences with respect to other strategies, such as the Arbitrary Lagrangian Eulerian (ALE) approach, where a moving mesh is considered.
Then we consider two Stokes-stable finite element spaces . We need a second mesh for the solid: this will be constructed on the reference configuration B and will be denoted by T B h .In general, we use simplicial meshes for the solid and the corresponding finite element space is defined as where T k denotes an element of T B h and M e is the total number of such elements.It turns out that the space semidiscretization of Problem 4 reads Problem 5. Given u h0 ∈ V h and X h0 ∈ S h , for almost every t ∈ ]0, T [, find where M denotes the number of degrees of freedom in the space S h , s i is the i-th vertex, X hi = X h (s i ), and F h = ∇ s X h .
Remark 1.The evaluation of F h in Problem 5 depends on the properties of the solid body and on the approximating spaces.In particular, since we assumed that X h is piecewise linear, than the deformation gradient F h is piecewise constant and so is the approximate Piola-Kirchhoff tensor P(F h ).Hence the source term F h can be computed as follows: (22) where [[P]], the jump of P across the interelement edge e, is defined as: and N + and N − are the normals to the interface e pointing outwards from the "+" or "−" element respectively.For more details, the interested reader is referred to [6].
A fully discretized scheme for the approximation of Problem 4 has been introduced in [6].If we denote by ∆t the time step size and by N the number of time steps, then the approximation of the term d can be performed as follows: where, as usual, the superscript n refers to discrete quantities evaluated after n time steps.
For n = 0, . . ., N − 1 the solution strategy can then be summarized as follows.
Step 1. Compute the source term Step 2. Solve the Navier-Stokes equations: find (u Step 3. Advance the position of the structure (pointwise): We observe that in Step 2 the unknown function u n+1 h (X n h (s)) appears in the integral on the right hand side.By interpolating the basis functions along the structure at time t n , this gives an additional linear contribution to the resulting algebraic system.

The finite element Immersed Boundary Method with distributed Lagrange multiplier
In this section we describe the modification of the fe-ibm introduced in [7].The main idea behind the new formulation consists in a different treatment of the interaction between the fluid and solid velocities.Namely, in Problem 4 the motion of the structure is designed by (2), which for each s ∈ B provides an ordinary differential equation.As a consequence in the fully discretized scheme we used (26) to update the position of each point of the structure.This gives rise to some restrictions on the choice of the discretization parameters as it will be shown in the next section.In this section we write (2) in variational form, so that we can have more flexibility in the choice of the finite elements to be used.
Let us introduce three functional spaces Λ, H 1 and H 2 and two bilinear forms Then equation ( 2) can be written and can be interpreted as a constraint on our system.Therefore we introduce a Lagrange multiplier associated to such constraint and split the first equation in Problem 4 into two separate equations, thus leading to the following dlm-ibm version of the problem.
, and λ(t) ∈ Λ, such that for almost every t ∈]0, T [ it holds (27) The finite element discretization of Problem 6 is straightforward.In addition to the finite element spaces introduced in Sect.3, we consider a finite element space Λ h ⊆ Λ, so that we have the following semidiscrete problem.
and λ h (t) ∈ Λ h , such that for almost every t ∈]0, T [ it holds (28) Remark 2. We observe that in this new formulation we do not need to evaluate the terms d and F, but the third equation represents the elasticity equation with respect to the position of the body.
Let us introduce now the time discretization based on the Euler scheme.As in the previous section, when computing along the structure terms involving functions in V h , we use the value of X at the previous time step, so that we have the following scheme.
Problem 8. Given u h0 ∈ V h and X h0 ∈ S h , suppose X 1 h ∈ S h has been assigned (it can be computed formally by assuming X −1 h = 0 in the scheme we are going to present).Find Problem 8 can be interpreted as a monolithic discretization of the fluid-structure problem and, in the case of a linear model for the Piola-Kirchhoff tensor P(F) = κF = κ ∇ s X, has the following matrix structure: where, denoting by ϕ, ψ, χ and ζ the basis functions respectively in V h , Q h , S h and Λ h , we have used the following notation: We observe that this system can be solved with a block iterative procedure which allows for the use of Navier-Stokes and elasticity solvers.The block structure of the matrix highlights the fact that at each time step we have to solve a saddle point problem whose analysis will be the object of a forthcoming paper.

Stability analysis
In this section we report the stability estimates for the ibm formulations of the fluid-structure interaction problem that we have introduced in the previous sections.The main results concerns the stability in time of the fully discrete schemes.We shall see that the dlm-ibm method is superior to the fe-ibm from this point of view since it is unconditionally stable, while the fe-ibm scheme requires that the discretization parameters are chosen in a appropriate way.
Since the solid is composed by a hyperelastic material, it is characterized by a positive energy density W (F) which depends only on the deformation gradient and the first Piola-Kirchhoff stress tensor can be obtained by derivation with respect to deformation gradient P(F(s, t)) = ∂W ∂F (F(s, t)).We assume that W is a C 1 convex function over the set of second order tensors.Moreover, the elastic potential energy of the body is given by: ( 30) It is not difficult to show the following energy estimate for the continuous versions of fe-ibm and dlm-ibm (see [6,7]).
) d be solution either of Problem 4 or Problem 6, then the following estimate holds true where • 0 and • 0,B denote the norms in L 2 (Ω) and L 2 (B), respectively.
Proof.Let us consider first Problem 4. We take v = u(t) in the first equation of ( 21), obtaining Then the definitions of the source terms d and F give: which concludes the proof.
For the stability of the solution of Problem 6 the proof is even simpler.It is enough to take as test functions v = u(t), q = p(t), Y = ∂X(t)/∂t, and µ = λ(t) in the variational equations in (27) and the result is achieved by summing up the equations and using the same computation as before to deal with the term containing the Piola-Kirchhoff stress tensor.space dim.solid dim.cfl condition x /h s Table 1.cfl condition according to fluid and structure dimensions.
The stability property of the semidiscrete problems can be obtained with the same arguments as in Proposition 1.When we consider the fully discrete schemes, namely ( 24)-( 26) for fe-ibm and Problem 8, the situation changes and we have different type of results.In the case of fe-ibm we state in the following proposition that the energy estimate holds true provided a cfl condition is satisfied.Proposition 3. Assume that the energy density W is a C 1 convex function.Given u 0h ∈ V h and X h0 ∈ S h , for n = 1, . . ., N let u n h ∈ V h , p n h ∈ Q h and X n h ∈ S h satisfy (24)- (26).Then the following energy estimate holds true where ν a is given by x ∆tL n C n e , L n is the maximum distance between any two consecutive vertices of the Lagrangian mesh and C n e stands for the maximum number of Lagrangian elements that touch the same Eulerian element at the given time step.
We summarize the cfl condition for different values of fluid and solid dimension in Table 1.
The fully discrete dlm-ibm scheme is unconditionally stable.
Proposition 4. Assume that the energy density W is a C 1 convex function.Let u n ∈ (H 1 0 (Ω)) d and X n ∈ (H 1 (B)) d for n = 0, . . ., N satisfy Problem 8 with X n ∈ W 1,∞ (B) d , then the following estimate holds true for all n = 0, . . ., N − 1 (35) ρ In this section we report on some numerical tests that we have performed during our research in the finite element approach to ibm.The aim of these tests is to show that the fe-ibm and dlm-ibm methods can be efficiently implemented and used for the approximation of fluid-structure interaction problems.In particular, the presented results have been already published in previous papers or have never been published before even if they have been obtained as the results of previous research.This is the case, for instance, of the snapshots of the presented animations.
We start by showing some snapshots taken from a three dimensional simulation involving the interaction of a codimension one closed solid surface and a fluid contained in a cubic box.The initial configuration is reported in Figure 1 and correspond to an ellipsoid stretched in one of the horizontal directions.The initial fluid is at rest, so that the system is driven only by the elastic force of the solid which is tending to its spherical equilibrium configuration.The evolution of the system is reported in Figure 2: as expected the solid tends to a sphere.
Our second simulation discusses the situation when more than one solid is present.More precisely, the initial configuration is reported in Figure 3 where two circular structures of codimension one are immersed in a fluid confined in a rectangular box.In the top half of the box the fluid is moving from right to left, while in the bottom part the fluid moves in the opposite direction.Figure 4 shows the evolution of the system: it can be appreciated that the two structures change their shape as a consequence of their interaction and the incompressibility of the fluid.
The next set of numerical tests discusses the mass conservation properties of the fe-ibm in the spirit of [3,4].In this case it has been shown that the mass conservation depends on the ability of the Stokes solver to provide a good discretization of the divergence free condition.It is clear that discontinuous pressure finite elements provide a better approximation of the divergence free constraint.In Figure 5 we compare continuous (solid) and discontinuous (dashed) pressure approximation for the fe-ibm.The test is the two dimensional analogue of the one reported in   The mass conservation property is expressed by the conservation of the area inside the structure.We use Hood-Taylor and P 1 -iso-P 2 − P 0 Stokes element for the continuous pressure simulations (see [1]) and the corresponding enhanced elements presented in [3] for the discontinuous case.The slope of the lines corresponds to the rate of mass loss with respect to time.
Figure 6 shows that in this respect computations obtained with dlm-ibm are better than the original fe-ibm.The mesh of the structure (an ellipse tending to a circle) is very coarse in order to emphasize the phenomenon and the used Stokes element is the enriched higher-order Hood-Taylor P 3 − (P c 2 + P 1 ).This aspect will be investigated in our future research.
The last set of numerical experiments corresponds to the stability of the fully discrete scheme.It is well known that non-implicit schemes usually result in severe restrictions to the time step.For the ALE formulation of fluid-structure interactions, it has been shown that this can produce a unconditionally unstable method [12].In Section 5 we recalled that, on the other hand, the semi-implicit formulation of our fe-ibm has been proved to be stable if a suitable cfl condition is satisfied.In Section 5 it is also recalled that the stability properties of the dlm-ibm are even better: in this case unconditional stability occurs, as it has been shown in [7].In Figures 7 and 8 we report the ratio of the total energy at time n over the initial energy.When the stability is violated, the ratio blows up; on the other hand the energy remains bounded if the method is stable.

Figure 1 .
Figure 1.fe-ibm simulation of an elastic ellipsoid: initial configuration Figures 1 and 2: an elliptic elastic string tends to a circular equilibrium configuration.

Figure 2 .Figure 3 .
Figure 2. fe-ibm simulation of an elastic ellipsoid: snapshots of the time evolution (left-to-right and top-to-bottom)

Figure 4 .
Figure 4. fe-ibm simulation of two circular structures moving in opposite directions: snapshots of the time evolution (left-to-right and top-to-bottom)

Figure 7 .
Figure 7. Energy ratio for codimension one structure.The structure elastic constant κ = 5, h x = 1/32, the fluid viscosity ν = 1, δρ = 0.3.The solid line correspond to the dlm-ibm scheme, while the dashed line marks the energy for the fe-ibm scheme.

Figure 8 .
Figure 8. Energy ratio for codimension zero structure.The structure elastic constant κ = 1, h s = 1/8, the fluid viscosity ν = 0.05, δρ = 0.3.The solid line correspond to the dlm-ibm scheme, while the dashed line marks the energy for the fe-ibm scheme.
X h (t) ∈ S h be solution of Problem 5 or Problem 7, then the following estimate holds true