An energy stable one-field monolithic arbitrary Lagrangian-Eulerian formulation for fluid-structure interaction

by   Yongxing Wang, et al.
University of Leeds

In this article we present a one-field monolithic finite element method in the Arbitrary Lagrangian-Eulerian (ALE) formulation for Fluid-Structure Interaction (FSI) problems. The method only solves for one velocity field in the whole FSI domain, and it solves in a monolithic manner so that the fluid solid interface conditions are satisfied automatically. We prove that the proposed scheme is unconditionally stable, through energy analysis, by utilising a conservative formulation and an exact quadrature rule. We implement the algorithm using both F-scheme and d-scheme, and demonstrate that the former has the same formulation in two and three dimensions. Finally several numerical examples are presented to validate this methodology, including combination with remesh techniques to handle the case of very large solid displacement.


page 16

page 20

page 23


An Energy Stable One-Field Fictitious Domain Method for Fluid-Structure Interactions

In this article, the energy stability of a one-field fictitious domain m...

Three-Field Fluid-Structure Interaction by Means of the Variational Multiscale Method

Three-field Fluid-Structure Interaction (FSI) formulations for fluid and...

A monolithic one-velocity-field optimal control formulation for fluid-structure interaction problems with large solid deformation

In this article, we formulate a monolithic optimal control method for ge...

A One-Field Energy-conserving Monolithic Fictitious Domain Method for Fluid-Structure Interactions

In this article, we analyze and numerically assess a new fictitious doma...

A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids

A variational modeling framework for hydraulically induced fracturing of...

Fluid-structure interaction simulations with a LES filtering approach in solids4Foam

We consider the interaction of an incompressible fluid described by a Le...

Numerical Study of Cosserat Fluid-Structure Interaction in a Monolithic Eulerian Framework

We propose a monolithic Eulerian variational formulation in non-classica...

1 Introduction

Numerical methods for Fluid-Structure Interaction (FSI) have been widely studied during the past decades, and a variety of methodologies have been developed in order to address different aspects of the FSI problem. However stability analyses of the existing numerical methods are rare especially when large solid deformation is involved. This paper is dedicated to establishing a robust stability analysis of a one-field monolithic FSI scheme in the Arbitrary Lagrangian-Eulerian (ALE) framework.

Monolithic methods have been regarded as the most robust FSI algorithms in the literature (Heil_2004; Heil_2008; Muddle_2012; Hecht_2017; Wang_2017; Wang_2019; hubner2004monolithic), which solve for the fluid and solid variables simultaneously in one equation system. Among these methodologies for FSI problems, the one-field approaches (Hecht_2017; Wang_2017) express the solid equation in terms of velocity, thus only solve for one velocity in the whole FSI domain. In this case the whole system can be solved similarly to a modified fluid problem, and the coupling conditions at fluid and solid interface are automatically satisfied.

The stability analysis when using the ALE framework is challenging, even for the pure fluid problem, due to the arbitrary moving frame (nobile1999stability; formaggia2004stability; bonito2013time). (Boffi_2016; Boffi_2015) present an energy stable Fictitious Domain Method with Distributed Lagrangian Multiplier (FDM/DLM), and (Hecht_2017; Pironneau_2016) present an energy stable Eulerian formulation by remeshing. In a previous study (Wang_2019) we analysed the energy stability for a one-field FDM method. In this article we extend this one-field idea to the ALE formulation, and the stability result is achieved by expressing the fluid and solid equations in a conservative formulation. In this sense, the formulation is similar to the one introduced in (Hecht_2017). However it differs from (Hecht_2017) in the following perspectives: (1) we formulate the solid in the reference domain and analyse in an ALE frame of reference, in which case the formulation and analysis are exactly the same for two and three dimensional cases, whereas (Hecht_2017) formulates and analyses everything in the current domain, for which the three dimensional case is significantly more complicated (chiang2017numerical)

; (2) we update the solid deformation tensor (the

-scheme) while (Hecht_2017) updates the solid displacement (the -scheme); (3) we implement the scheme by solving an additional solid-like equation at each time step in order to move the mesh, whilst (Hecht_2017) implements their scheme by remeshing which is expensive in the three dimensional case.

The paper is organized as follows. In Section 2 the control equations for the FSI problem are introduced in an ALE framework. In Section 3 the finite element weak formulation is introduced, followed by spatial and time discretisations in Section 4. The main results of energy stability are presented in Section 5. Implementation details are considered in Section 6 and numerical examples are given in Section 7, with some conclusions in Section 8.

2 The arbitrary Lagrangian-Eulerian description for the FSI problem

Let and be the fluid and solid domain respectively (which are time dependent regions), is the moving interface between the fluid and solid, and has an outer boundary , which can be fixed or moving as shown in Figure 1. The Eulerian description is convenient when we observe a fluid from a fixed frame, while the Lagrangian description is convenient when we observe a solid from a frame moving with it. An ALE frame of reference can be adopted when a fluid and solid share an interface and interact with each other as shown in Figure 1, in which case the frame moves arbitrarily from a reference configuration , chosen to be the same as the initial configuration at , to a current configuration . Let us define a family of mappings :


with being the dimensions. We assume that is one-to-one and invertible with continuous inverse . Hence a point has a unique image at time , i.e.


and a point at time has a unique inverse image


We call the Eulerian coordinate, and call its inverse image , via the above mapping , the ALE coordinate. We assume that is differentiable with respect to for all , and define the velocity of the ALE frame as


Given an Eulerian coordinate , its corresponding ALE coordinate should be distinguished from its material (or Lagrangian) coordinate as shown in Figure 1. In fact (not necessarily the same as ) maps to via the Lagrangian mapping, i.e., the trajectory of a material particle at :


and the velocity of the material particle at is defined by

Figure 1: ALE mapping from to . Also shows the comparison between ALE mapping and Lagrangian mapping with Eulerian coordinate , ALE coordinate and material (Lagrangian) coordinate . and , .
Remark 1.

Although the Lagrangian configuration and the ALE configuration are not generally the same, both are chosen to have the initial configuration in this article. We shall also construct the ALE mapping such that coincides with at all boundaries including the fluid-solid interface: and .

Remark 2.

The ALE mapping is the mapping that is actually used to move the domain in this article, and the purpose of introducing the Lagrangian mapping is to discuss its related variables, such as particle velocity and solid deformation tensor , which will be defined in the following context.

Formulated in the current configuration, the conservation of momentum takes the same form in the fluid and solid:


with , , and being the density, gravity acceleration, velocity and Cauchy stress tensor respectively. Here we use the notation , with the superscript and denote fluid and solid respectively, and similar notations are also applied to and . In the above, is the total derivative computed along the trajectory of a material particle at , i.e. via the Lagrangian mapping:


Replacing the above partial time derivative by the total derivative of


leads to the ALE formulation of (7)


We consider here both an incompressible flow and incompressible solid:


with being the deviatoric part of the stress tensor. For a Newtonian fluid in ,


and for a hyperelastic solid (belytschko2013nonlinear) in ,




being the deformation tensor of the solid, being the determinant of F, and being the energy function of the hyperelastic solid material. Combining with the continuity equation


the FSI system is completed with continuity of the velocity and normal stress conditions on the interface :


and (for simplicity of this exposition) homogeneous Dirichlet and Neumann boundaries on and respectively:


with as shown in Figure 1.

3 Finite element weak formulation

Let be the square integrable functions in domain , endowed with norm . Let with the norm denoted by . We also denote by the subspace of whose functions have zero value on the Dirichlet boundary of .

According to equation (2) we construct from , so a function is one-to-one corresponding to a function via


Choosing a test function , the weak formulation may be obtained by multiplying on both sides of equation (10), and integrating the stress term by parts in domain and separately:


We used in the above deduction. Using the boundary conditions (16) and (17), we have the following equation by adding up (19) and (20).


Using Jacobis formula (magnus2019matrix), we have


with Then we can take the time derivative outside the moving domain (conservative formulation (nobile1999stability)),


Substituting (23) into (21), using


and combining the weak form of continuity equation (15), leads to the weak formulation of the FSI problem:

Problem 1.

Given , , and an ALE mapping (consequently given by (4)), : find and , such that , and , , the following equations hold:




with and being the initial interface and outer boundary respectively, as shown in Figure 1, and being the Lagrangian mapping as defined in (5).

4 Discretisation in space and time

Define a stable finite element space, such as the Taylor-Hood elements, for the velocity-pressure pair in :


with and being the number of nodal variables for each velocity component and pressure respectively. Then


Using the backward Euler scheme, equation (25) and (26) can be discretised respectively as follows:




In the above


and is a quadrature formula used to compute . In order to have an unconditionally stable scheme, which will be proved in Section 5, the mid-point integration is adopted for


in the two dimensional case, and the Simpson formula is adopted in the three dimensional case:


Due to the definition of the deformation tensor (14) and ALE velocity (4), we have




Therefore and in (28) can be updated as follows:




Up to now we have not stated how to construct (or ), because very often we only need to construct the ALE mapping at a discrete time level, that is to say computing for at each time step. This will be explained in the rest of this section.

We solve the following static linear elastic equation in in order to compute , and take for . Given the following boundary data:




find such that , the following equation holds:


with and being the Lam constants used here as pseudo-solid parameters. It is well known that the above elliptic problem (37) to (39) has a unique solution (brenner2007mathematical). As a result, we are able to construct a mapping for ,


and further


From the computational point of view, knowing the ALE velocity at the discrete level is sufficient.

Putting all the above together, the discrete ALE-FSI problem reads:

Problem 2.

Given and , find , , and (consequently an ALE mapping by (41)), such that , , , and , the following equation system holds:


with quadrature formula (31) in 2D or (32) in 3D, updating by (35) and updating by (36). In addition, the above FSI system equations are completed with the Dirichlet and Neumann boundary conditions (17) for the momentum and continuity equations (28) and (29), and with the boundary conditions (37) and (38) for the mesh equation (39).

Problem 2 is a highly non-linear system, so we solve it iteratively as described in the following Algorithm 1.

0:  , and a tolerance tol
0:  , and
     1. solve the mesh equation (39) for using boundary conditions (37) and (38)
     2. update using (36)
     3. solve the FSI system (28) and (29) for and
     4. ,
Algorithm 1 Solve Problem 2 for (or ), and

5 Stability analysis

We shall deduce an energy stability result at the end of this section. In preparation for this we first prove the following lemmas.

Lemma 1.

If is the solution of Problem 2, then satisfies the following at .


Noticing that


and integrating by parts:


In the above , thanks to the enclosed flow (17). Using the Sobolev imbedding theorem (mitrovic1997fundamentals, Theorem 6 in Chapter 5), we have in the two dimensional case and in the three dimensional case. Either or is included in because has finite measure. Therefore , and thanks to (29). ∎

Lemma 2.

If is the solution of Problem 2 then, for any , satisfies the following at .


Integrating by parts we get


The boundary integral in (47) is zero due to the enclosed flow condition (17). The second term on the right-hand side of (47) can be expressed as:


we then have (46) by substituting (48) into (47). ∎

Lemma 3.

If is the solution of Problem 2, then






where is the cofactor matrix of . According to the way we construct (40), we know is a polynomial in time of degree (nobile1999stability), with being the space dimension. Also is a constant for , so is linear in time when and quadratic when , and a mid-point integration () or Simpson formula () would exactly compute . This is to say


Noticing that for ,


and using (52), we finally have (49). ∎

Lemma 4.

Define potential energy of the solid:


If is the solution of Problem 2 and is convex on the set of second order tensors (Boffi_2016), then






Due to the convexity assumption of , we have


This gives:


Using (35) we have


which finally leads to (55) by integrating (60) in