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

In this article, the energy stability of a one-field fictitious domain method is proved and validated by numerical tests in two and three dimensions. The distinguishing feature of this method is that it only solves for one velocity field for the whole fluid-structure domain; the interactions remain decoupled until solving the final linear algebraic equations. To achieve this the finite element procedures are carried out separately on two different meshes for the fluid and solid respectively, and the assembly of the final linear system brings the fluid and solid parts together via an isoparametric interpolation matrix between the two meshes. The weak formulations are introduced in the continuous case and after discretization in time. Then the stability is analyzed through an energy estimate. Finally, numerical examples are presented to validate the energy stability properties.

• 7 publications
• 5 publications
• 4 publications
05/10/2017

### 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...
03/08/2020

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

In this article we present a one-field monolithic finite element method ...
10/07/2020

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

Three-field Fluid-Structure Interaction (FSI) formulations for fluid and...
04/28/2021

### Existence, uniqueness, and approximation of a fictitious domain formulation for fluid-structure interactions

In this paper we describe a computational model for the simulation of fl...
03/23/2021

06/02/2021

### 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...
03/29/2022

### A Multigrid Preconditioner for Spatially Adaptive High-order Meshless Method on Fluid-solid Interaction Problems

We present a monolithic geometric multigrid preconditioner for solving f...

## 1 Introduction

Three major questions arise when considering a finite element method for the problem of Fluid-Structure Interactions (FSI): (1) what kind of meshes are used (interface fitted or unfitted); (2) how to couple the fluid-structure interactions (monolithic/fully-coupled or partitioned/segregated); (3) what variables are solved (velocity and/or displacement). Combinations of the answers of these questions lead to different types of numerical method. For example, (1; 2) solve for fluid velocity and solid displacement sequentially (partitioned/segregated) using an Arbitrary Lagrangian-Eulerian (ALE) fitted mesh, whereas (3; 4; 5) use an ALE fitted mesh to solve for fluid velocity and solid displacement simultaneously (monolithic/fully-coupled) with a Lagrange Multiplier to enforce the continuity of velocity/displacement on the interface. The Immersed Finite Element Method (IFEM) (6; 7; 8; 9; 10; 11; 12) and the Fictitious Domain Method (FDM) (13; 14; 15; 16; 17; 18) use two meshes to represent the fluid and solid separately. Although IFEM could be monolithic (6), the classical IFEM only solves for velocity, while the solid information is arranged on the right-hand side of the fluid equation as a prescribed force term. Although the FDM may be partitioned (18), usually the FDM approach solves for both velocity in the whole domain (fluid plus solid) and displacement of the solid simultaneously via a distributed Lagrange multiplier (DLM) to enforce the consistency of velocity/displacement in the overlapped solid domain. In the case of one-field and monolithic numerical methods for FSI problems, (19) introduces a 1D model using a one-field FDM formulation based on two meshes, and (20; 21) introduces an energy stable monolithic method (in 2D) based on one Eulerian mesh and discrete remeshing.

In a previous study (22), we present a one-field monolithic fictitious domain method (subsequently referred to as the one-field FDM) which has the following main features: (1) only one velocity field is solved in the whole domain, based upon the use of an appropriate projection; (2) the fluid and solid equations are solved monolithically. Our motivation for proposing the one-field FDM is based on comparing its features with those of existing numerical schemes. Compared with IFEM the classical IFEM does not solve the solid equation (8; 9; 10; 11; 12). Instead, the solid information is arranged on the right-hand side of the fluid equation as a prescribed force. The one-field FDM solves the solid equation together with the fluid equation in one discretized linear algebraic system. The similarity is that both methods only solve for velocity and pressure fields (no solid displacement). DLM/FDM methods (13; 14; 15; 16; 17; 18) solve the solid equation, but for a displacement field, and couple this displacement with the velocity of the fictitious fluid via a Lagrange multiplier. This leads to a large discretized linear algebra system. The one-field FDM rewrites the solid equation in terms of a velocity variable and couples the fictitious fluid through a finite element interpolation. Monolithic Eulerian methods (20; 21) also express the solid equation in terms of velocity, and the fluid and solid are coupled naturally on an interface-fitted mesh. The one-field FDM uses two meshes to represent the fluid and solid respectively. Consequently, before discretization in space, these two methods have many similarities, the advantage of the one-field FDM being that interface fitting is not required.

The main developments in this paper, following from (22), are as follows. The energy preserving property in the continuous case is proved. The energy nonincreasing property after time discretization is proved, and the same property is also proved after spatial discretization. The implementation in this paper is based on an

-scheme, i.e., the solid deformation tensor

is updated (see section 4), while the previous paper uses a -scheme (see equation (29) in (22)). The advantage of this -scheme is that the integral is expressed in the reference domain so that it becomes linear for the neo-Hookean solid model (see equation (15)), which however is nonlinear if expressed in the current domain for the -scheme (see equation (29) in (22)). The methodology and analysis is demonstrated to extend to the three-dimensional case.

The paper is organized as follows. Control equations and weak formulation are introduced in section 2 and 3 respectively. The time discretized weak form is then presented in section 4. Stability of the proposed scheme is analyzed in section 5. Space discretization is discussed in section 6. Numerical examples are given in section 7, and conclusions are presented in section 8.

## 2 Control equations

In the following context, and with denote the fluid and solid domain respectively which are time dependent regions as shown in Figure 1. is a fixed domain (with outer boundary ) and is the moving interface between fluid and solid. We denote by the reference (material) coordinates of the solid, by the current coordinates of the solid, and by the initial coordinates of the solid.

Let denote the density, viscosity, velocity, pressure and stress tensor respectively. We assume both an incompressible fluid and incompressible solid, then the conservation of momentum and conservation of mass take the same form as follows:

Momentum equation:

 ρdudt=∇⋅σ, (1)

Continuity equation:

 ∇⋅u=0. (2)

An incompressible Newtonian constitutive equation in can be expressed as:

 σ=σf=μfDuf−pfI, (3)

where . We shall use an incompressible neo-Hookean solid in (14; 20), and in common with previous work (10; 12) we also assume the solid has the same viscosity as the fluid. The constitutive equation may be expressed as:

 σ=σs=c1J−1(FFT−I)+μfDus−psI, (4)

where = is the deformation tensor of the solid, and is the determinant of . Finally the system is complemented with the following boundary and initial conditions.

 uf=usonΓt, (5)
 nsσf=nsσsonΓt, (6)
 uf=0onΓ, (7)
 uf∣∣t=0=uf0, (8)
 us|t=0=us0. (9)

Other boundary conditions are possible on but (7) are used here for simplicity.

###### Remark 1.

The corresponding energy function for the hyperelastic stress in (4) is defined by (16):

 Ψ(F)=c12(trFFT−d)−c1ln(J). (10)

## 3 Weak formulation

The finite element weak form discussed in this section is almost the same as that in (22), the only difference is that we integrate the solid stress in the reference domain, because we shall update the deformation tensor (-scheme) rather than the solid stress as done in (22) (-scheme). In the following context, 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 values on the boundary of , and denote by the subspace of whose functions have zero mean value.

Let . Given , we perform the following symbolic operations:

 ∫ΩftEq.(???)⋅vdx+∫ΩstEq.(???)⋅vdx.

Integrating the stress terms by parts, the above operations, using constitutive equation (3) and (4) and boundary condition (6), gives:

 ρf∫Ωdudt⋅vdx+μf2∫ΩDu:Dvdx−∫Ωp∇⋅vdx+ρδ∫Ωstdudt⋅vdx+c1∫ΩstJ−1(FFT−I):∇vdx=0, (11)

where . Note that the integrals on the interface are cancelled out using boundary condition (6). This is not surprising because they are internal forces for the whole FSI system considered here.

Transforming the integral of the last two terms of (11) to the reference coordinate system, combined with the following symbolic operations for ,

 −∫ΩftEq.(???)qdx−∫ΩstEq.(???)qdx,

leads to the weak form of the FSI system as follows.

###### Problem 1.

Given and , find , and , such that for , , the following two equations hold:

 ρf∫Ω∂u∂t⋅vdx+ρf∫Ω(u⋅∇)u⋅vdx+μf2∫ΩDu:Dvdx−∫Ωp∇⋅vdx+ρδ∫ΩsX∂u∂t⋅vdX+c1∫ΩsXF:∇XvdX−c1∫ΩstJ−1∇⋅vdx=0, (12)

and

 −∫Ωq∇⋅udx=0. (13)
###### Remark 2.

Because domain is stationary (the Eulerian description will be used) and is transient which will be updated by its own velocity (the updated Lagrangian description), there is a convection term from the total derivative of time in , but there is no convection term in .

###### Remark 3.

Problem 1 is equivalent to the equation (12) in (22).

## 4 Discretization in time

We may use the backward Euler method to discretize Problem 1, and update coordinates of the solid by . As a result, is updated by , and so,

 ∫ΩsXFn+1:∇Xv=∫ΩsXFn:∇Xv+Δt∫ΩsX∇Xun+1:∇Xv. (14)

Using equation (14), the discretized weak form corresponding to Problem 1 may be expressed as:

###### Problem 2.

Given , and , find , and , such that for , , the following four relations hold:

 (15)
 −∫Ωq∇⋅un+1dx=0, (16)
 Ωsn+1={x:x=xn+Δtun+1,xn∈Ωsn}, (17)

and

 Fn+1=Fn+Δt∇Xun+1. (18)
###### Remark 4.

We shall use a fixed-point iteration at each time step to construct implicitly.

###### Remark 5.

Problem 2 is similar to equation (30) in (22), however here the discretized weak form is expressed as an implicit scheme, and the solid deformation tensor is updated rather than the solid stress in (22).

## 5 Stability by energy estimate

### 5.1 Energy conservation in the continuous case

In this section we shall prove that the weak forms (12) and (13), associated with Problem 1, preserve energy.

###### Lemma 3.

The energy function for the hyperelastic stress satisfies:

 c1∫t0∫ΩsXF:∇XudX−c1∫t0∫ΩstJ−1∇⋅udx=∫ΩsXΨ(F)dX. (19)
###### Proof.

Since and . Using the fact that ( and are arbitrary matrices), we have:

 ddt∫ΩsXΨ(F)dX=∫ΩsX∂Ψ∂F:dFdtdX=c1∫ΩsX(F−F−T):ddt(I+∇Xd)dX=c1∫ΩsXF:∇XudX−c1∫ΩstJ−1∇⋅udx,

where is displacement of the solid at time . ∎

###### Lemma 4.

If is the solution pair of Problem 1, then

 ∫Ω(u⋅∇)u⋅udx=0. (20)
###### Proof.

First,

 ∫Ω(u⋅∇)u⋅udx=∫Ω∇(u⊗u)⋅udx−∫Ω|u|2∇⋅udx. (21)

Integrate by parts:

 ∫Ω∇(u⊗u)⋅udx=∫Γ|u|2u⋅ndΓ−∫Ω(u⋅∇)u⋅udx. (22)

According to a Sobolev imbedding theorem (23, Theorem 6 in Chapter 5) and the inclusion between spaces ( if ), we know (for both 2D and 3D). Therefore , i.e., . That is to say . Then we have from (13). We also have from the boundary condition (7). Substituting these two equations into (21) and (22) gives equation (20). ∎

###### Proposition 5 (Energy Conservation).

Let be the solution pair of Problem 1, then

 ρf2∫Ω|u|2dx+μf2∫t0∫ΩDu:Dudx+ρδ2∫ΩsX|u|2dX+∫ΩsXΨ(F)dX=0. (23)
###### Proof.

We first let in (12) and integrate from time to , then let in (13) and substitute into (12). Finally we can construct the above equation of energy balance due to Lemma 3 and 4. ∎

### 5.2 Stability analysis after time discretization

We next demonstrate a similar energy stability result for Problem 2.

###### Lemma 6.

The trace function satisfies:

 12tr(Fn+1FTn+1)−12tr(FnFTn)=ΔtFn+1:∇Xun+1−Δt22|∇Xun+1|2, (24)

where for an arbitrary matrix .

###### Proof.
 Fn+1FTn+1−FnFTn=Fn+1FTn+1−(Fn+1−Δt∇Xun+1)(Fn+1−Δt∇Xun+1)T=ΔtFn+1∇TXun+1+Δt∇Xun+1FTn+1−Δt2∇Xun+1∇TXun+1.

Lemma 6 holds due to

 12tr(Fn+1FTn+1−FnFTn)=Δt⋅tr(Fn+1∇TXun+1)−Δt22|∇Xun+1|2.

###### Lemma 7.

The log-determinant function satisfies:

 ln(detFn+1)−ln(detFn)≥Δt∇⋅un+1−Δt22∣∣F−1n+1∇Xun+1∣∣2.
###### Proof.

Use the fact that function is concave over the set of positive definite matrices (24, Chapter 3). Let , and , then

According to the property of concave functions, we have , this is to say:

In the above, we use the trace property of cyclic permutations: . ∎

From the above two lemmas, we have:

###### Proposition 8.

The energy function for the hyperelastic stress satisfies:

 ∫ΩsXΨ(Fn+1)dX−∫ΩsXΨ(Fn)dX≤Δtc1∫ΩsXFn+1:∇Xun+1dX−Δtc1∫Ωsn+1J−1n+1∇⋅un+1dx+Rn+1, (25)

where

 Rn+1=c1Δt22∫ΩsX(∣∣F−1n+1∇Xun+1∣∣2−|∇Xun+1|2)dX. (26)

Similarly to Lemma 4, we have:

###### Lemma 9.

If is the solution pair of Problem 2, then

 ∫Ω(un+1⋅∇)un+1⋅un+1dx=0. (27)
###### Proposition 10 (Energy Nonincreasing).

Let be the solution pair of Problem 2. If , then

 ρf2∫Ω|un+1|2dx+ρδ2∫ΩsX|un+1|2dX+∫ΩsXΨ(Fn+1)dX+Δtμf2n+1∑k=1∫ΩDuk:Dukdx≤ρf2∫Ω|un|2dx+ρδ2∫ΩsX|un|2dX+∫ΩsXΨ(Fn)dX+Δtμf2n∑k=1∫ΩDuk:Dukdx+Rn+1, (28)

where is defined in equation (26).

###### Proof.

Let in (15) and multiply on both side of the equation, and then let in (16) and substitute into equation (15), we get:

 ρf∫Ω(un+1−un)⋅un+1dx+Δtμf2∫ΩDun+1:Dun+1dx+ρδ∫ΩsX(un+1−un)⋅un+1dX+c1Δt∫ΩsXFn+1:∇Xun+1dX−c1Δt∫Ωsn+1∇⋅un+1dx=0. (29)

Using the Cauchy-Schwarz inequality and the fact , we have:

where or . Substituting the above relation into (29), we get (28) due to Proposition 8 and Lemma 9. ∎

###### Remark 6.

Relation (28) does not exactly show energy nonincreasing, because we do not know whether is greater or less than . However, is and this will be demonstrated in section 7 by numerical tests. In order to test the energy property, let us use the following notation for the different contributions to the total energy in (28): (1) Kinetic energy of fluid plus fictitious fluid (2) Kinetic energy of solid minus fictitious fluid (3) Viscous dissipation (4) Potential energy of the solid Denote the total energy as and the energy ratio as:

 Eratio=Etotal(tn)Etotal(t0). (30)

We shall numerically demonstrate that is nonincreasing in section 7.

## 6 Discretization in space

We shall use a fixed Eulerian mesh for and an updated Lagrangian mesh for to discretize Problem 2. First, we discretize as with the corresponding finite element spaces as

 Vh(Ωh)=span{φ1,⋯,φNu}⊂H10(Ω)

and

 Lh(Ωh)=span{ϕ1,⋯,ϕNp}⊂L20(Ω).

The approximated solution and can be expressed in terms of these basis functions as

 uh(x)=Nu∑i=1u(xi)φi(x),ph(x)=Np∑i=1p(xi)ϕi(x). (31)

We further discretize as with the corresponding finite element spaces as:

 Vsh(Ωsh0)=span{φs1,⋯,φsNs}⊂H1(Ωs0),

then move the vertices of each element of by their own velocities to get , and approximate as:

 ush(x)=Ns∑i=1uh(xsi)φsi(x)=Ns∑i=1Nu∑j=1u(xj)φj(xsi)φsi(x), (32)

where is the nodal coordinate of the solid mesh. Notice that the above approximation defines an projection from to : ,

We then discretize Problem 2 in space as follows.

###### Problem 11.

Given , and , find , and , such that for , , the following four relations hold:

 ρf∫Ωhuhn+1−uhnΔt⋅vdx+ρf∫Ωh(uhn+1⋅∇)uhn+1⋅vdx+μf2∫ΩhDuhn+1:Dvdx−∫Ωhphn+1∇⋅vdx+ρδ∫ΩshXushn+1−ushnΔt⋅vsdX+c1Δt∫ΩshX∇Xushn+1:∇XvsdX−c1∫Ωshn+1J−1n+1∇⋅vsdx=−c1∫ΩshXFshn:∇XvsdX, (33)
 −∫Ωq∇⋅uhn+1dx=0, (34)
 Ωshn+1={x:x=xn+Δtushn+1,xn∈Ωshn}, (35)

and

 Fshn+1=Fshn+Δt∇Xushn+1, (36)

where and .

###### Remark 7.

The proof of the energy estimate (Proposition 10) for the spatially continuous case can also be applied to the discrete case (see A).

###### Remark 8.

There are two sources of nonlinearity in Problem 11: the convection term and the moving solid domain. We can accommodate these by moving the convection term to the right-hand side of the equation, and using a fixed-point iteration to construct in order to solve the nonlinear system at each time step. For other methods to treat convection, readers may refer to (25; 26). We shall only use this fully implicit implementation to consider low Reynolds number () cases in this paper in order to test the energy stability.

###### Remark 9.

A two-step explicit splitting scheme (-scheme) is discussed in B with corresponding energy analysis. This scheme is similar to that in (22) (-scheme), which may be adapted to problems at large Reynolds number (see (22) for more examples).

## 7 Numerical experiments

In this section, we focus on validation of the energy stability of the proposed numerical method in two and three dimensions. For more two-dimensional numerical examples and validation of the basic algorithm see (22). We shall use linear triangles (2D) or linear tetrahedra (3D) to discretize the solid domain . In domain , the elements will be used, i.e., the standard - element is enriched by a constant for approximation of the pressure. This element has the property of local mass conservation and the constant may better capture the element-based jump of pressure (27; 28). We shall demonstrate the improvement of mass conservation and energy conservation by using the elements compared to the elements. We shall also validate that the total energy is nonincreasing as stated in Proposition 10 and Remark 6.

### 7.1 Oscillating disc driven by an initial kinetic energy (activated disc)

In this test, we consider an enclosed flow () in with a periodic boundary condition. A solid disc is initially located in the middle of the square and has a radius of . The initial velocity of the fluid and solid are prescribed by the following stream function

 Ψ=Ψ0sin(ax)sin(by),

where and . In this test, , , and . In order to visualize the flow a snapshot of the velocity and deformation fields is presented in Figure 2, and the evolution of energy is presented in Figure 3 using a mesh (biquadratic squares for the fluid velocity and bilinear triangles for the solid velocity).

We commence by comparing elements and elements. The evolution of mass variation and energy ratio are demonstrated in Figure 4, from which it can be seen that the enrichment of the pressure field by a constant has an effect of stabilizing the mass and energy evolution. In addition, this enrichment of the pressure field dramatically improves the mass conservation, although the effect for energy conservation is not obvious. Then using element , time convergence of the total energy can be observed from Figure 5 (a), from which we can see a nonincreasing energy and a first order time convergence for both the implicit and explicit scheme (see B for the energy estimate of the explicit scheme). It can be seen from Figure 5 (b) that the residual term defined in (26) is very small and converges rapidly to zero when reducing ().

### 7.2 Oscillating disc driven by an initial potential energy (stretched disc)

In the previous example, the disc oscillates because a kinetic energy is prescribed for the FSI system at the beginning. In this test, we shall stretch the disc and create a potential energy in the solid, then release it causing the disc to oscillate due to this potential energy. The computational domain is a square . One quarter of a solid disc is located in the left-bottom corner of the square, and initially stretched as an ellipse as shown in Figure 6. Notice the equation of an ellipse and its area , hence we ensure that this stretch does not change mass of the solid.

We choose , , and . The fluid adopts a mesh of biquadratic squares, and the solid has similar node density ( linear triangles) as the fluid. A snapshot of pressure on the fluid mesh and corresponding solid deformation with its velocity norm are displayed in Figure 7, and the evolution of energy is presented in Figure 8. The nonincreasing total energy can be observed from Figure 9 (a) for both the implicit and explicit scheme (see B for energy estimate of the explicit scheme). It can be seen from Figure 9 (b) that the residual term defined in (26) is very small and converges rapidly to zero when reducing .

### 7.3 Oscillating ball driven by an initial kinetic energy

In this section, we consider a 3D oscillating ball, which is an extension of the example in section 7.1. The ball is initially located at the center of with a radius of . Using the property of symmetry this computation is carried out on of domain : . The initial velocities of and components are the same as that used in section 7.1 and the component is set to be 0 at the beginning. We adopt the same parameter and mesh size defined in section 7.1 (with the same mesh size in the z direction). A snapshot of the solid ball and the corresponding fluid velocity norm are presented in Figure 10, and the nonincreasing energy property is presented in Figure 11.

## 8 Conclusions

In this article, we first introduce an implicit version of (22) for the one-field fictitious domain method (one-field FDM) based upon updating the solid deformation tensor . Then the energy-preserving property for this one-field FDM is proved on the continuous level, and the energy-nonincreasing property is proved after discretization in time and space. The energy property for an explicit scheme is also analyzed in B. Finally, a selection of numerical tests are presented to demonstrate this theoretical energy estimate in both two and three dimensions. It has therefore been demonstrated that the proposed one-field FDM is a stable, robust and computationally efficient technique for the solution of a wide range of fluid-structure interaction problems.

## Appendix A Stability analysis after space discretization

As with the previous stability estimate (Proposition 10) after time discretization, we have the following estimate after space discretization.

###### Proposition 12.

Let be the solution pair of Problem 11. If , then

 ρf2∫Ωh∣∣uhn+1∣∣2dx+ρδ2∫ΩshX∣∣ushn+1∣∣2dX+∫ΩshXΨ(Fshn+1)dX+Δtμf2n+1∑k=1∫ΩhDuhk:Duhkdx≤ρf2∫Ωh∣∣uhn∣∣2dx+ρδ2∫ΩshX∣∣uhn∣∣2dX+∫ΩshXΨ(Fshn)dX+Δtμf2n∑k=1∫ΩhDuhk:Duhkdx+Rhn+1, (37)

where

 Rhn+1=c1Δt22∫ΩshX(∣∣∣(Fshn+1)−1∇Xushn+1∣∣∣2−∣∣∇Xushn+1∣∣2)dX. (38)
###### Proof.

Let in (33) and multiply on both side of the equation, and then let in (34) and substitute into equation (33), we get:

 ρf∫Ωh(uhn+1−uhn)⋅uhn+1dx+Δtμf2∫ΩhDuhn+1:Duhn+1dx+ρδ∫ΩshX(ushn+1−ushn)⋅ushn+1dX+c1Δt∫ΩshXFshn+1:∇Xushn+1dX−c1Δt