Combined field formulation and a simple stable explicit interface advancing scheme for fluid structure interaction

# Combined field formulation and a simple stable explicit interface advancing scheme for fluid structure interaction

## Abstract

We develop a combined field formulation for the fluid structure (FS) interaction problem. The unknowns are , being the fluid velocity, fluid pressure and solid velocity. This combined field formulation uses Arbitrary Lagrangian Eulerian (ALE) description for the fluid and Lagrangian description for the solid. It automatically enforces the simultaneous continuities of both velocity and traction along the FS interface. We present a first order in time fully discrete scheme when the flow is incompressible Navier-Stokes and when the solid is elastic. The interface position is determined by first order extrapolation so that the generation of the fluid mesh and the computation of are decoupled. This explicit interface advancing enables us to save half of the unknowns comparing with traditional monolithic schemes. When the solid has convex strain energy (e.g. linear elastic), we prove that the total energy of the fluid and the solid at time is bounded by the total energy at time . Like in the continuous case, the fluid mesh velocity which is used in ALE description does not enter into the stability bound. Surprisingly, the nonlinear convection term in the Navier-Stokes equations plays a crucial role to stabilize the scheme and the stability result does not apply to Stokes flow. As the nonlinear convection term is treated semi-implicitly, in each time step, we only need to solve a linear system (and only once) which involves merely if the solid is linear elastic. Two numerical tests including the benchmark test of Navier-Stokes flow past a Saint Venant-Kirchhoff elastic bar are performed. In addition to the stability, we also confirm the first order temporal accuracy of our explicit interface advancing scheme.

F

luid Structure Interaction; Arbitrary Lagrangian Eulerian; Navier-Stokes Equations; Saint Venant-Kirchhoff;

## 1 Introduction

The interaction of a rigid or deformable solid with its surrounding fluid or the fluid it enclosed gives rise to a very rich variety of phenomena. For example, a rotating fan, a vibrating aircraft wing, a swimming fish. The daily activities of many parts of our human body are also intimately related to this fluid solid or more commonly called fluid structure (FS) interaction. For example, the heart beating, respiration, speaking, hearing, and even snoring. To understand those phenomena, we need to model both the fluid and the solid. In this paper, we assume the fluid is incompressible and the solid is deformable. The governing equations for FS interaction are then as follows:

 ρf(∂t\boldmathu+\boldmathu⋅∇\boldmathu)=∇⋅\boldmathσf+ρf\boldmathgf,∇⋅\boldmathu=0 in % Ωf(t), (1.1) ρs∂tt\boldmathφ=∇⋅% \boldmathσs+ρs\boldmathgs in Ωs. (1.2)

In the above equations, is the fluid domain at time and is the initial configuration of the solid. is the velocity at a spatial point in . is the position at time of a material point in . Constants and are the densities of the fluid and the solid. As one can tell, the fluid is described by the Eulerian (spatial) description while the solid is described by the Lagrangian (material) description with reference configuration . is the stress of the fluid in the Eulerian description and is the stress of the solid in the Lagrangian description. If the fluid is assumed to be viscous Newtonian,

 \boldmathσf=\boldmathσf(\boldmathu,p)=2ρfνf\boldmathϵ(\boldmathu)−pI=ρfνf(∇\boldmathu+∇\boldmathu⊤)−pI. (1.3)

If the solid is elastic, then there is a strain energy with density :

 Is(\boldmathφ(⋅,t))=∫ΩsW(∇% \boldmathφ(\boldmathz,t))d\boldmathz, (1.4)

so that the stress is determined by . Here is the deformation gradient. As a result, satisfies

 ddϵ∣∣∣ϵ=0Is(\boldmathφ+ϵ\boldmathϕ)=⟨\boldmathσs(% \boldmathφ),∇\boldmathϕ⟩Ωs∀\boldmathφand\boldmathϕ. (1.5)

Like in , we write instead of for simplicity. Here and . Let be the strain tensor. If the strain of the solid is small, we can model the solid as Saint Venant-Kirchhoff material with . Simple calculation shows

 \boldmathσs(\boldmathφ)=2μs% \boldmathF\boldmathE+λs(tr\boldmathE)\boldmathF. (1.6)

If the deformation of the solid is small, one popular choice of is

 \boldmathσs(\boldmathφ)=μs(∇% \boldmathη+∇\boldmathη⊤)+λs(∇⋅\boldmathη)\boldmathI, (1.7)

where is the displacement vector. We can define an associated with and check . We can also verify the convexity of :

 Missing or unrecognized delimiter for \left

But we will see later that is not convex (see (3.40)). Indeed, it is not even polyconvex [9, sec 3.9]

The fluid and the solid are coupled by the continuity of velocity and traction across the FS interface. Mathematically, it can be written as (see Antman [1, p.485, , p.489, (15.31)])

 \boldmathu(\boldmathφ(\boldmathz,t),t)=∂t\boldmathφ(\boldmathz,t)∀\boldmathz∈Γ, (1.8) ∫\boldmathφ(\boldmathe,t)% \boldmathσf(\boldmathx,t)\boldmathnfda(% \boldmathx)+∫\boldmathe\boldmathσs(% \boldmathz,t)\boldmathnsda(\boldmathz)=0∀\boldmathe⊂Γ. (1.9)

Here is the FS interface at and is any part of it. is the image of at time under the mapping and it is part of the fluid boundary. and are the outward normals. Besides the interface conditions (1.8) and (1.9), we also need boundary conditions on the rest of the boundaries. These fixed boundaries are called to : , . We require

 \boldmathu|Σ1=\boldmathub,% \boldmathσf\boldmathnf|Σ2=\boldmathσfb,\boldmathφ|Σ3=\boldmath% φb,\boldmathσs\boldmathns|Σ4=\boldmathσsb. (1.10)

Here , , and are prescribed fluid velocity, fluid traction, solid position and solid traction on the ’s.

The well-posedness of FS interaction problems has been studied for various models (see [11, 16, 10] and the references in [10]). When the solid is deformable, the analysis uses Lagrangian description for both the fluid and the solid [16, 10]. But for computational efficiency, we use Arbitrary Lagrangian Eulerian (ALE) description for fluid in this paper. We learn from [14] the idea of using conservative ALE description with proper intermediate mesh or meshes so that the stability bound obtained does not explicitly depend on the mesh movement. But [14] only considers the convection diffusion equation on a domain whose motion is known.

We will show that (1.9) is equivalent to (2.18) while the latter is more popular in the literature. Schemes based on {(1.1),(1.2),(1.8), (2.18), (1.10)} can be classified into two classes: partitioned schemes and monolithic schemes. In a partitioned scheme ([5, 12]), different solvers are used for fluid and elastic equations. For stability reason, one may hope that the two interface conditions (1.8) and (2.18) are satisfied simultaneously. But this cannot be achieved with one single iteration between the fluid and the solid solvers. Various improvements have been proposed: [13, 2] reduce the fluid-structure coupling to pressure-structure coupling after a temporal discretization; [17] proposes a beautiful stable splitting scheme for fluid-membrane interaction problem; [15] discusses how to achieve fast convergence by properly choosing the boundary conditions for different solvers. So far, all the numerical analysis requires that the fluid domain does not move and the flow is Stokes type [13, 17, 15]. These assumptions conceal the important fact that the convection term in the Navier-Stokes equations indeed stabilizes the scheme when the domain moves. In a monolithic scheme ([18, 4]), the governing equations for the fluid and the solid as well as the governing equation for the displacement of the mesh are coupled and solved all together. Hence a large nonlinear system is required to solve in each time step [4, remark below (3.11)] [18, section 3.5]. Since the interface conditions (1.8) and (2.18) are built into the system, they are automatically satisfied once the system is solved. Then it is easy to believe that monolithic schemes will preserve the stability of the associated continuous models. However, as the discrete schemes so far proposed are very complicated, we are not aware of any proof of existence and stability of numerical solutions in the literature.

The celebrated immersed boundary method of Peskin [25] uses delta function to represent the force at the FS interface. When the solid is codimension-1, i.e., a surface in or a curve in , there comes the immersed interface method of LeVeque and Li [19] which is spatially more accurate. In this method, instead of using delta function, the forces are translated into the jump conditions across the interface. These jump conditions are then taken care by changing the discretization of the differential operators at stencils acrossing the interface [19, 3]. To our point of view, there are still some aspects left to be improved for the methods initiated by Peskin, LeVeque and Li: Immersed boundary method in general is only first order in space [25, p.500,p.509], [3, p.4]; Immersed interface method cannot handle solid with a finite volume.

In this paper, we will present a new monolithic scheme. It is based on a new formulation of FS interaction which in some aspect is similar to the well-known one field formulation of multi-fluid flows (see [29, (1)] for example, but you will not see any delta function in our method). In our scheme, the unknown is , being the fluid velocity, fluid pressure and solid velocity. There are many nice features of our scheme:

• (Explicit mesh moving) Our scheme is explicit interface advancing which means that the FS interface at time is constructed explicitly using only information of the solid at time . So, determining the FS interface and computing are decoupled.

• (Smaller system) Let be the fluid velocity. Let and be the solid velocity and position. Let and be the fluid mesh velocity and position. Then, ignoring the pressure for simplicity, [18, (37)]’s unknown is , and [4, (3.10)]’s unknown is . Our unknown is . In our scheme, the mesh related information and the solid position are all treated explicitly.

• (Linear system which is solved only once) Our scheme only asks to solve a linear system once per time step when it is Navier-Stokes flow and linear elastic structure coupling. But [18, 4] ask to solve a nonlinear system with unknowns mentioned before per time step.

• (Easy implementation) As our scheme is Jacobian free and does not require characteristics, it is very easy to implement.

• (Stability) The most important feature of our scheme is its stability: Even though the fluid mesh is constructed explicitly, the total energy of the fluid and the solid at time is bounded by the total energy at time . The stability bound hence obtained does not explicitly depend on the fluid mesh velocity which however is explicitly used in the computation through the ALE formulation (see Theorem 4 which is our main theorem). This feature is shared by the original continuous model (see Theorem 2.3). Surprisingly, the nonlinear convection term in the Navier-Stokes equations plays a crucial role in proving Theorem 4 and the stability result does not apply to Stokes flow.

However, as the energy bound alone at is not enough to prevent the FS interface from colliding with itself or other fixed boundaries at , we have to assume (but see Remark 4.1) that the we take satisfies the following collision free condition so that we are able to construct the fluid mesh at time :

 the time step Δt that we take will not make % the FSinterface collides with itself or other fixed boundaries. (1.11)

As the interface is determined by extrapolation, one may wonder what is the point to use a monolithic scheme — can the unknowns be easily solved for separately? If one try that, the scheme becomes a loosely coupled partitioned scheme. Consequently, there will be time lag in the enforcement of the two continuity conditions. How this time lag affects the stability will not be addressed in this paper and we refer to [5, 12] and the references therein. By paying the price of solving a larger system, monolithic scheme enables the satisfaction both (1.8) and (2.18) at the same time which is the key to get stability in all density ratio regime. How to solve the system for efficiently will not be addressed in this paper, even though we do have a linear system when it is Navier-Stokes and linear elasticity coupling.

As extrapolation is used to determine the interface position, one may wonder whether it will damage the accuracy. We note that the position of any point satisfies where is the solid velocity. We know using the right end point rule (hence implicit) to approximate the integral is not necessarily more accurate than using the left end point rule (hence explicit). Obviously, left end point rule is cheaper. Less obviously, as we have mentioned and will prove later, left end point rule is also stable if it is handled properly. We will present numerical test that verifies the first order temporal accuracy in Section 5 (see Table 2 and Fig. 3). Higher order schemes will be discussed in a forthcoming paper.

The rest of the paper is organized as follows: In Section 2 we introduce the combined field formulation and its weak form for FS interaction, and also discuss the conservative ALE description of it. We present our scheme in Section 3, starting from the discretization of the solid part. Then we discuss the fluid mesh construction, the ALE mapping and its applications. In particular, we can define the backward in time ( ) and the forward in time () extension (see (3.52) and (3.63)). Before we introduce our scheme (3.68), we mention two results (Lemma 3.2.4 and Corollary 3.2.4) which are related to Geometric Conservation Law and present Lemma 3.2.4 which contributes to some crucial cancelation that will be used in the stability proof. After this long preparation, we prove the stability of our scheme in Section 4. This is our main theoretical result (Theorem 4). Two numerical tests will be discussed in Section 5. They verfies the stability as well as the first order temporal accuracy of our scheme.

We emphasize that (1.9) does not contain any Jacobian and it fits perfectly into the framework of (isoparametric) finite element methods. For (1.9), the author thanks Dr. Stuart Antman for the teaching of elasticity around year 2002 which motivates most of the work in this paper.

## 2 Combined field formulation using ALE description

We show that the FS system {(1.1),(1.2),(1.8),(1.9),(1.10)} has a very clean and simple weak formulation. Then we put it into ALE format using the conservative formulation. To save space, we will not discuss the non-conservative formulation but refer the readers to [21].

### 2.1 Traction boundary condition

We first show that we can insert test functions into the traction boundary condition (1.9) ([1, page 489]). {lemma} Let be the FS interface at . If (1.9) is true, then for any and any that is defined on ,

 ∫\boldmathφ(\boldmathe,t)\boldmathϕf(\boldmathx)⋅(\boldmathσf(% \boldmathx,t)\boldmathnf)da(\boldmathx)+∫\boldmathe\boldmathϕs(\boldmathz)⋅(\boldmathσs(\boldmathz,t)\boldmathns)da(\boldmathz)=0, (2.12)

where is defined on and satisfies .

{proof}

Consider . Suppose is the image of the mapping with . Then, when , the mapping maps to and then to . So, we can change all the integrations to :

 ∫\boldmathφ(\boldmathe,t)\boldmath% ϕf(\boldmathx)⋅(\boldmathσf(\boldmathx)\boldmathnf)da(\boldmathx)=∫\boldmathπ\boldmathϕf(% \boldmathx(\boldmaths))⋅(\boldmathσf(\boldmathx(\boldmaths))(∂\boldmathx∂s1×∂% \boldmathx∂s2))ds1ds2, (2.13)
 ∫\boldmathe\boldmathϕs(\boldmathz)⋅(\boldmathσs(\boldmathz)\boldmathns)da(\boldmathz)=−∫\boldmathπ\boldmathϕs(\boldmathz(\boldmaths))⋅(\boldmath% σs(\boldmathz(\boldmaths))(∂\boldmathz∂s1×∂\boldmathz∂s2))ds1ds2. (2.14)

With the same idea, (1.9) can be rewritten as

 ∫\boldmathπ\boldmathσf(\boldmathx(\boldmaths))(∂\boldmathx∂s1×∂\boldmathx∂s2)−\boldmathσs(\boldmathz(\boldmaths))(∂\boldmathz∂s1×∂% \boldmathz∂s2)ds1ds2=0. (2.15)

As the above equation is true for any , the integrand must be zero. Then, since , we get

 Missing or unrecognized delimiter for \right (2.16)

Now, putting (2.13),(2.14),(2.16) together, we obtain (2.12).

Remark: Since and , (2.15) can be written as

 ∫\boldmathπ\boldmathσf(\boldmathx)(det∂\boldmathx∂\boldmath% z)(∂\boldmathx∂\boldmathz)−⊤\boldmathns−\boldmathσs(% \boldmathz)\boldmathnsda(\boldmathz)=0. (2.17)

So, we can use

 (det∂\boldmathx∂\boldmathz)\boldmathσf(\boldmathx(\boldmathz))(∂\boldmathx∂\boldmathz)−⊤\boldmathns−\boldmathσs(% \boldmathz)\boldmathns=0 on Γ (2.18)

as the boundary condition for the continuity of traction (e.g. [1, page 489, (15.34)]). However, as (2.18) contains Jacobian, it is not as friendly to numerics as (2.12).

### 2.2 Combined field formulation

We will use velocity field as the unknown for both fluid and solid because it leads to a simple enforcement of the boundary condition (1.8). The combined field formulation for time dependent FS interaction is as follows:

 ρf(∂t\boldmathu+\boldmathu⋅∇\boldmathu)=∇⋅\boldmathσf(% \boldmathu,p)+ρf\boldmathgf,∇⋅% \boldmathu=0 in Ωf(t), (2.19) ρs∂t\boldmathv=∇⋅% \boldmathσs(\boldmathφ0+∫t0\boldmathv(τ)dτ)+ρs\boldmathgs in% Ωs, (2.20) \boldmathu(\boldmathφ(\boldmathz,t),t)=\boldmathv(\boldmathz,t)∀\boldmathz∈Γ, (2.21) ∫Γ\boldmathϕs(\boldmathz)%\boldmath$σ$s(\boldmathz,t)\boldmathnsda(% \boldmathz)+∫\boldmathφ(Γ,t)\boldmath% ϕf(\boldmathx)\boldmathσf(\boldmathx,t)\boldmathnfda(\boldmathx)=0∀% \boldmathϕs, (2.22) \boldmathu|Σ1=\boldmathub,\boldmathσf\boldmathnf|Σ2=% \boldmathσfb,\boldmathv|Σ3=∂t\boldmathφb,\boldmathσs% \boldmathns|Σ4=\boldmathσsb, (2.23)

where is defined by through for any . is the FS interface at . . and are fixed boundaries. At the initial time , we are given , and being the initial fluid velocity, initial solid velocity and initial solid position.

Now we will derive a weak form of (2.19)–(2.23). Suppose we know a function defined on , introduce

 V\boldmathφ(t)={(\boldmathu;p;\boldmathv): \boldmathu∈H1(Ωf(t)),p∈L2(Ωf(t)),\boldmathv∈W1,∞(Ωs), \boldmathu(\boldmathφ(\boldmathz))=\boldmathv(\boldmathz)∀\boldmathz∈Γ}. (2.24)

Note that and are defined on different domains (see Fig. 1 for an illustration). is given from the very beginning and will never change. is varying with respect to . We require since can be rather nonlinear. But for linear elasticity (1.7), requiring is enough.

{theorem}

The system (2.19)–(2.23) has the following weak form: We are looking for with , , such that for any ,

 ⟨ρf(∂t\boldmathu+% \boldmathu⋅∇\boldmathu),\boldmathϕf⟩Ωf(t)+⟨\boldmathσf(\boldmathu,p),∇\boldmathϕf⟩Ωf(t)−⟨∇⋅\boldmathu,qf⟩Ωf(t) +⟨ρs∂t\boldmathv,\boldmathϕs⟩Ωs+⟨\boldmathσs(% \boldmathφ0+∫t0\boldmathv(τ)dτ),∇\boldmathϕs⟩Ωs =⟨ρf\boldmathgf,\boldmathϕf⟩Ωf(t)+⟨\boldmathσfb,% \boldmathϕf⟩Σ2+⟨ρs\boldmathgs,\boldmathϕs⟩Ωs+⟨\boldmathσsb,\boldmathϕs⟩Σ4 (2.25)

for any satisfying , . The in is defined by . {proof} The boundary condition (2.21) has been built into the function space . We dot (2.19a) with , dot (2.19b) with and dot (2.20) with . After integration by parts for the and terms, we add the resulting equations together. The boundary integrals on and cancel each other because of (2.22) and the definition of .

### 2.3 Stability identity for (2.25)

We want to derive some stability identity for (2.25) when the solid has a strain energy (1.4). Once again, we denote by . With the defined in (1.4),

 ddtIs(\boldmathφ(⋅,t))=∫Ωs∂W(\boldmathF)∂\boldmathF:∇∂t\boldmathφ(\boldmathz,t)d\boldmathz=⟨\boldmathσs(\boldmathφ(⋅,t)),∇% \boldmathv(⋅,t)⟩Ωs. (2.26)

We have the following identity for divergence free velocity field which can be applied to the first term in (2.25) when and : {lemma} Consider a divergence free velocity field defined on a time varying domain . Assume at any point on , either or where is the outward normal and is the velocity that moves. Then

 12ddt∥\boldmathu∥2Ωf(t)=⟨∂t\boldmathu+\boldmathu⋅∇\boldmathu,\boldmathu⟩Ωf(t). (2.27)
{proof}

Recall the following Reynolds transport theorem ([1, page 488, (15.23)]): for any ,

 ddt∫Ωf(t)η(\boldmathx,t)d% \boldmathx=∫Ωf(t)∂tη(\boldmathx,t)+∫∂Ωf(t)η(\boldmathx,t)\boldmathw(\boldmathx,t)⋅\boldmathn. (2.28)

So,

 12ddt∥\boldmathu∥2Ωf(t)= 12ddt∫Ωf(t)|\boldmathu(\boldmathx,t)|2d\boldmathx=∫Ωf(t)∂t\boldmathu⋅\boldmathu+12∫∂Ωf(t)|\boldmathu|2\boldmathw⋅\boldmathn = ∫Ωf(t)∂t\boldmathu⋅\boldmathu+12∫∂Ωf(t)|% \boldmathu|2\boldmathu⋅\boldmathn. (2.29)

In the last step we have used the condition either or on . Using the divergence free condition, the last surface integral can be rewritten as

 12∫∂Ωf(t)|\boldmathu|2\boldmathu⋅\boldmathn=∫Ωf(t)12∇⋅(|\boldmathu|2\boldmathu)d\boldmathx=∫Ωf(t)(\boldmathu⋅∇\boldmathu)⋅\boldmathud\boldmathx. (2.30)

We have used in the last step. Plugging (2.30) into (2.29), we get (2.27).

Combining (2.26) and (2.27), if we let , and in (2.25), we obtain the following result: {theorem} When the solid has strain energy (1.4) and when , , in (2.23), the solution of (2.25) satisfies

 ρf2 ddt∥\boldmathu(⋅,t)∥2Ωf(t)+ρfνf2∥∇\boldmathu+∇\boldmathu⊤∥2Ωf(t)+ρs2ddt∥% \boldmathv(⋅,t)∥2Ωs +ddtIs(\boldmathφ0+∫t0\boldmathv(τ)dτ)=⟨ρf\boldmathgf,\boldmathu⟩Ωf(t)+⟨ρs\boldmathgs,\boldmathv⟩Ωs+⟨\boldmathσsb,\boldmathv⟩Σ4. (2.31)

### 2.4 Conservative formulation of Arbitrary Lagrangian Eulerian (ALE) description

When the time reaches (which can be any number), we choose the as the reference domain and construct a backward in time mapping that maps to for any . See (3.47) for a way of constructing . Note that . Define

 \boldmathwn(\boldmathx,t)=∂tΦn(% \boldmathx,t). (2.32)

With this fixed, we set in (2.25) and choose test function satisfing

 ddt\boldmathϕf(Φn(\boldmathx,t),t)=0∀x∈Ω(tn). (2.33)

In a finite element method, requiring (2.33) for a basis function means the function is always attached to the nodal point it starts from no matter how the mesh moves. See (3.52),(3.53). Evaluating (2.33) at and using , we obtain

 ∂t\boldmathϕf(\boldmathx,tn)=−(∂tΦn(\boldmathx,tn))⋅∇\boldmathϕf(Φn(\boldmathx,tn),tn)=−\boldmathwn(% \boldmathx,tn)⋅∇\boldmathϕf(\boldmathx,tn). (2.34)

As is exactly the velocity of domain , we can apply the Reynolds transport theorem (2.28) with to get

 ddt∣∣∣t=tn⟨\boldmathu(⋅,t),\boldmathϕf(⋅,t)⟩Ωf(t)=∫Ωf(tn)(∂t\boldmathu)⋅% \boldmathϕf+\boldmathu⋅(∂t\boldmathϕf)+∇⋅(\boldmathu⋅\boldmathϕf\boldmathwn), (2.35)

where the integrand on the right hand side is evaluated at . Because of (2.34), the sum of the last two terms in the integrand equals . So, using , (2.35) leads to

 ⟨∂t\boldmathu(⋅,tn)+% \boldmathu(⋅,tn)⋅∇\boldmathu(⋅,tn),% \boldmathϕf(⋅,tn)⟩Ωf(tn) = ddt∣∣∣t=tn⟨\boldmathu(⋅,t),\boldmathϕf(⋅,t)⟩Ωf(t)+⟨∇⋅(\boldmathu⊗(\boldmathu−% \boldmathwn(⋅,tn))),\boldmathϕf(⋅,tn)⟩Ωf(tn) (2.36)

To summarize, if with , is a solution to (2.19)–(2.23), then satisfies the following: For any , for any given backward in time mapping , for any defined in the space-time domain of the fluid, and satisfying and , for any defined in and satisfying and