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 NavierStokes 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 NavierStokes 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 semiimplicitly, 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 NavierStokes flow past a Saint VenantKirchhoff elastic bar are performed. In addition to the stability, we also confirm the first order temporal accuracy of our explicit interface advancing scheme.
luid Structure Interaction; Arbitrary Lagrangian Eulerian; NavierStokes Equations; Saint VenantKirchhoff;
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:
(1.1)  
(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,
(1.3) 
If the solid is elastic, then there is a strain energy with density :
(1.4) 
so that the stress is determined by . Here is the deformation gradient. As a result, satisfies
(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 VenantKirchhoff material with . Simple calculation shows
(1.6) 
If the deformation of the solid is small, one popular choice of is
(1.7) 
where is the displacement vector. We can define an associated with and check . We can also verify the convexity of :
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)])
(1.8)  
(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
(1.10) 
Here , , and are prescribed fluid velocity, fluid traction, solid position and solid traction on the ’s.
The wellposedness 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 fluidstructure coupling to pressurestructure coupling after a temporal discretization; [17] proposes a beautiful stable splitting scheme for fluidmembrane 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 NavierStokes 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 codimension1, 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 wellknown one field formulation of multifluid 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.

(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 NavierStokes 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 :
(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 NavierStokes 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.
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 nonconservative 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 ,
(2.12) 
where is defined on and satisfies .
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 :
(2.13) 
(2.14) 
With the same idea, (1.9) can be rewritten as
(2.15) 
As the above equation is true for any , the integrand must be zero. Then, since , we get
(2.16) 
Now, putting (2.13),(2.14),(2.16) together, we obtain (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:
(2.19)  
(2.20)  
(2.21)  
(2.22)  
(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
(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.
The system (2.19)–(2.23) has the following weak form: We are looking for with , , such that for any ,
(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),
(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
(2.27) 
Recall the following Reynolds transport theorem ([1, page 488, (15.23)]): for any ,
(2.28) 
So,
(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
(2.30) 
We have used in the last step. Plugging (2.30) into (2.29), we get (2.27).
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
(2.32) 
With this fixed, we set in (2.25) and choose test function satisfing
(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
(2.34) 
As is exactly the velocity of domain , we can apply the Reynolds transport theorem (2.28) with to get
(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
(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 spacetime domain of the fluid, and satisfying and , for any defined in and satisfying and