Ying He ·Jie Shen

Abstract

Keywords Fluid-structure interaction·Pressure correction·Stability analysis

1 Introduction

Fluid-structure interaction (FSI) plays an important role in many scientific/engineering applications,e.g.,design of engineering systems,blood flow in human arteries,etc.It has been extensively studied in recent years both analytically and computationally (cf.[6,9,11,18] and the references therein).

There are mainly three approaches,monolithic,partitionedand semi-implicit projection,for solving FSI problems numerically.The partitioned approach (cf.,for instance,[2,4,10,24]) solves the fluid and structure dynamics separately with explicit interface conditions.While each subproblem can be solved efficiently by existing algorithms,the explicit treatment of the interface condition may lead to instability in the presence of strong addedmass eff ect [5] and requires very restrictive time step constraint.In contrast,the monolithic approach (cf.,for instance,[19,20,23]) simultaneously solves the fluid and structure dynamics coupled by the implicit interface conditions.This type of schemes usually have good stability properties,but at each time step,a nonlinear coupled system has to be solved and,due to the presence of the pressure in the coupled system,it is usually difficult to design an eff ective iterative scheme to solve the nonlinear coupled system.On the other hand,the semi-implicit projection approach was first proposed in [12].It decouples the computation of fluid velocity from that of the pressure and structure displacement by using a projection method.This method appears to have some computational advantage over the partitioned or monolithic approaches (cf.,for instance,[1,3,12]).

In this paper,we shall construct a diff erent class of semi-implicit projection schemes which decouple the computation of pressure from that of the velocity and structure displacement.Our schemes will be computationally very efficient.More precisely,in the first step of our schemes,we solve a coupled,but elliptic,system for an intermediate fluid velocity and the structure displacement; then in the second step,we solve a Poisson equation for the fluid pressure and obtain the fluid velocity with a simple correction.Furthermore,we shall also prove rigorously that these schemes are unconditionally stable.

For fluid problems,an eff ective approach to decouple the computation of the pressure from that of the velocity is to use a projection-type method,originally proposed by Chorin and Temam in the late 1960s [7,28].A comprehensive review on various projection-type methods can be found in [14].However,a main difficulty in the design of a projection method for the FSI problem is to assign a boundary condition for the pressure at the interface.It is well known that a proper boundary condition for the pressure Poisson equation in a projection-type method,at the Dirichlet part of the boundary,is the homogeneous Neumann boundary condition.Indeed,most existing projection-type schemes for the FSI problem also use,explicitly or implicitly,the Neumann-type boundary condition for the pressure Poisson equation at the Dirichlet part of the boundary as well as at the interface.However,imposing a Neumann-type boundary condition for the pressure at the interface appears to aff ect,to a certain degree,the stability of the scheme,and we are not aware of any proof of unconditional stability for this type of projection schemes,only a conditional stability has been proved in [12] for a linear FSI problem.In a previous paper [17],the authors constructed an unconditionally stable scheme for a linear FSI problem.The aim of this paper is to extend it to a nonlinear FSI problem.

In [13],the authors proposed and analyzed pressure-correction schemes for Navier–Stokes equations with open boundary where the usual stress-free boundary condition is applied.It is shown that the proper boundary condition at the open boundary is of Dirichlet type instead of Neumann type.Two schemes are constructed in [13],one is based on the standard pressure correction which leads to poor accuracy at the open boundary,while the other is based on the rotational pressure correction and with a proper Dirichlet boundary condition at the open boundary.It is shown in [13] that both the standard and rotational pressure-correction projection schemes,when applied to the time-dependent Stokes problem,are unconditionally stable,but the rotational version leads to much better accuracy.Since one of the matching interface conditions for the FSI problem is related to the stress,it makes sense to extend the approach in [13] for problems with open boundary to the FSI problem.

Besides the difficulty associated with the pressure boundary condition on the interface,another major difficulty is to prove the unconditional stability of the rotational pressure-correction scheme for the nonlinear FSI problem.The original stability proof of the rotational pressure-correction scheme in [15] was only valid for Stokes problems.An essential step of the proof was to take the “discrete time derivative” of the scheme.Unfortunately,this proof cannot be extended to the nonlinear case.In [8],the authors constructed an unconditionally stable rotational velocity-correction scheme for the Navier–Stokes equations.However,they only provided a stability proof for the linear Stokes equations,while showing numerically that the scheme was unconditionally stable.In [25],the author proposed a Gauge–Uzawa approach for the rotational pressure-correction scheme of the Navier–Stokes equations,and proved that the scheme was unconditionally stable.We shall extend the approach in [25] for the Gauge–Uzawa scheme of the Navier-Stokes equations to the rotational pressure-correction schemes for the FSI problem.

To fix the idea,we consider in this paper a simple model of the FSI problem where the movement of the interface is assumed infinitesimal so the interface is treated as fixed.This nonlinear FSI problem captures many of the essential difficulties of the more general FSI problems with moving interface,and its well-posedness has been studied in [22].

The rest of the paper is organized as follows.In the next section,we describe the governing equations for our FSI model,formulate its weak form and the energy dissipation law.In Sect.3,we construct a standard and rotational pressure-correction scheme for the FSI problem and prove their unconditional stability.Then,in Sect.4,we describe a generic approach for spatial discretization as well as a Fourier–Legendre method for a special case of a periodic channel.We present some numerical results in Sect.5 to validate our numerical schemes and to demonstrate their temporal accuracy.Some concluding remarks are given in Sect.6.

2 Governing Equations

We consider the following model for interaction of a viscous fluid with an elastic body in a two- or three-dimensional bounded domain Ω ,with the fluid region Ωf,the solid region Ωsand the interface Γc,so we have Ω=Ωf∪Ωs∪Γc.We also denote Γf=Ωf∖Γcand Γs=Ωs∖Γc(cf.Fig.1).

We assume that the interface undergoes infinitesimal displacements,i.e.,Γcis fixed.The more complicated situation with moving interface will be considered in a forthcoming paper.

Fig.1 Geometry discription for fluid-structure problem

In the fluid region Ωf,we have the Navier–Stokes equations: whereudenotes the fluid velocity,pthe fluid pressure,u0is the given initial velocity,f1is the given body force per unit mass,is the strain tensor,andandare the constant fluid density and viscosity.

In the solid region Ωs,we have the wave equation for linear elasticity:

wherewdenotes the displacement of the solid,w0andw1are the given initial data,and(w) is the elastic stress tensor,given by

f2is the given loading force per unit mass,andare the Lamé constants,andis the constant solid density.

and

For simplicity,we take in this paperWe further take=1 and=0 which implies div(w)=Δw,and the interface condition (2.4) reduces to

To derive a weak formulation for (2.1)–(2.2),we need to introduce some notations.Let us denote byHk(Ω) andHk0(Ω) (fork≥0 ) the standard Sobolev spaces,equipped with the standard norm ‖⋅‖k,Ω.In particular,we denoteL2(Ω)=H0(Ω) with the associated norm ‖⋅‖ .We will useto denote the vector-valued Sobolev spaces.We also denote

Then,a weak solution (u,p,w) for (2.1)–(2.2) will satisfy

with the interface conditions (2.3) and (2.5) on Γc.

Using (2.5),we can reformulate the above as

withu=wton the interface Γc.

and summing up the two resultant equations,we obtain

or equivalently

where

is the total energy of the FSI system.

For the well-posedness of the system (2.7),we refer to [22].

3 Time Discretization

For FSI problems,it is very important to design numerical schemes which have good,preferably unconditional,stability property.Usually,this is achieved by fully coupled,implicit schemes which require solving,at each time step,a coupled,nonlinear,saddle-point system.

We construct in this section time discretization schemes based on the standard and rotational pressure-correction approach for (2.7).These schemes are unconditionally stable and lead to,at each time step,a coupled,linear elliptic system in Ω and a pressure Poisson equation in Ωf,which can be efficiently solved by standard numerical methods.The stability analysis for each scheme is carried out in this section.

3.1 Standard Pressure-Correction Scheme

We first construct a first-order scheme for the FSI problem based on the standard pressurecorrection scheme for the Navier–Stokes problem with the open boundary condition [13]:

Step 1Given (un,pn,wn) ,computeandby solving

This is a coupled,linear elliptic system forwith the coupling condition at the interface Γc.Hence,it can be efficiently solved,for example,by a standard domain decomposition approach (cf.,for instance,[26,29]).

Step 2Computeby solving

We observe that a Dirichlet boundary condition is imposed forpn+1on the interface Γc,as opposed to the usual Neumann boundary condition in a pressure-correction formulation.This is due to the interface condition (2.5) which is similar to the open boundary condition considered in [13].

Hence,we only have to solve a Poisson equation at this step.

For the above scheme,we have the following result:

Theorem 3.1The scheme (3.1)–(3.3),withp0|Γc=0 ,is unconditionally stable.More precisely,if we define the discrete energy

then we have,for alln≥0,

ProofTo simplify the notations,we define,for any sequence {uk} ,the discrete time deriva-tives

Taking and taking the inner product of (3.2a) withq=2Δt∇pn,summing up the three relations,we obtain

Rewrite (3.2a) as

Taking the inner product with itself from both sides and integrating by parts,thanks topk|Γc=0 for allk(due top0|Γc=0 ),andwe obtain

Summing up (3.5) and (3.7),we obtain

which implies the desired result.

We recall that due to the artificial Dirichlet boundary condition for the pressure in (3.2c),a higher-order discretization for the velocity will not increase the accuracy.Hence,to obtain a higher-order scheme,one needs to resort to the rotational pressure-correction (cf.[13]).

3.2 Rotational Pressure-Correction Schemes

3.2.1 First-Order Scheme

We start by constructing a first-order scheme.

Step 1Given (un,vn,wn,pn),computeby solving

Step 2Computeby solving

We observe that the main diff erence of the rotational scheme (3.8)–(3.9) with the standard scheme (3.1)–(3.3) is the additional termThis term replaces the artificial Dirichlet B.C.by an improved B.C.On the other hand,the numerical procedures for the two schemes are essentially identical.

The proof of unconditional stability for the rotational scheme is much more difficult.The original stability proof of the rotational pressure-correction scheme in [15] was carried out only for Stokes problems,and an essential step of the proof was to take the “discrete time derivative” of the scheme.Unfortunately,this proof cannot be extended to the nonlinear case.However,we can prove that the above rotational scheme is unconditionally stable using a similar procedure to that in the proof below for the second-order rotational scheme.We omit the details for the sake of brevity.

3.2.2 Second-Order Scheme

We observe that it is not straightforward to construct a second-order version of (3.8)–(3.9) using the usual backward diff erence formula (BDF).Hence,we first introduce an additional variablev=wtand rewrite the FSI equations as

with the boundary condition

and the initial condition

We can now construct a second-order rotational pressure-correction scheme as follows:

Step 1Given (un,wn,vn,pn),computeandby solving

Step 2Compute (un+1,pn+1) by solving

• One observes that all the terms,except the pressure,are discretized with a secondorder BDF or Adams–Bashforth formula.We recall that a first-order treatment of the pressure term,coupled with second-order treatment for other terms,can lead to second-order accuracy for the velocity [14].

• It is clear that,at each time step,the numerical procedure for solving (3.13)–(3.14) is essentially the same as for the first-order scheme (3.1)–(3.3).

• In [25],the author proved the unconditional stability for a Gauge–Uzawa scheme of the Navier–Stokes equations.A useful idea in [25] is to introduce a sequence {qn} defined by

We shall also use this sequence in our stability proof below.

Theorem 3.2The scheme (3.13)–(3.14),withp−1|Γc=p0|Γc=0 ,is unconditionally stable.More precisely,if we define the discrete energy as

then we have

ProofFor any sequencewe have

LetIn1(u) ,In2(u) andIn3(u) be the last three terms in the right-hand side.Using the algebraic identities

and

we find

Using the first equation in (3.14a),we have

where,by (3.20),

Using (3.15),we can rewrite (3.14a) as

Taking the inner product with itself from both sides of the above equation,integrating parts and using (3.15) and the fact that (pk+qk)|Γc=…=(p0+q0)|Γc=0 ,we obtain

Multiplying the above by Δtand adding it to (3.21),we obtain

Thanks to (3.15),we have

where we have used the well-known Korn’s inequalitywithd=2 or 3 being the space dimension.

Finally,using the above inequality,(3.20) and (3.23) in (3.25),we find

which implies the desired result.

Remark 3.3With the stability results established in this section,it is also possible to derive similar error estimates for these schemes as in [13].

4 Galerkin-Type Spatial Discretization and Implementation

We briefly describe a general procedure to implement the time discretization schemes constructed in the last section.Letandbe some finite dimensional approximation spaces,withpreferably satisfying the Babuska–Brezzi inf-sup condition.We also denoteWe note that one can generalize the stability proofs for the semi-discretized schemes in the last section to their full discretized versions using the above discrete settings; we refer to [14] for more detail in this regard.

To fix the idea,we take the scheme (3.8)–(3.9) as an example.The other schemes can be treated by using exactly the same procedure.

4.1 A General Setup

A Galerkin approximation of the scheme (3.8)–(3.9) is as follows:

Step 1LetThen we look forsuch that

and

Define

and

Then,we can rewrite (4.1) as: Findsuch that

Thus,the equation (4.4) can be viewed as a two-domain approximation to a linear elliptic problem with the discontinuous coefficientNote that from (4.1b),is continuous at Γc.Hence,one can efficiently solve the coupled linear system using a standard domain decomposition approach.In particular,in the two-dimension case,one can form the Schurcomplement to solve the unknown at the interface first,and then solve for the velocity in the fluid region and displacement in the solid region separately (cf.,for instance,[26,29] and a simple example in the next subsection).

Step 2Findsuch that

whereQhis anL2-projection operator ontoMh.

We note that (4.5) is just a discrete Poisson equation in Ωfwith the homogeneous Dirichlet boundary condition on Γc,and (4.6) involves only a projection,so they can be effi-ciently solved.

4.2 An Example with a Fourier–Legendre Approximation

As an example,we consider a two-dimensional periodic channel withWe denotebyWe assume that all functions are periodic in thex-direction.

Leth=(M,N),whereMis the number of equally spaced points in thex-direction,andN+1 is the number of Legendre–Gauss–Lobatto points in they-direction of Ωfand Ωs.For simplicity,we use the same number of points in they-direction of Ωfand Ωs.LetPNbe the set of all polynomials of degree less than or equal toN.We set

For the sake of efficiency and to take full advantage of periodicity in thex-direction,we shall treat the nonlinear convective term in (4.1) explicitly.To this end,we modify (4.2) to

With this modification we expand all the functions in discrete Fourier series,e.g.,

The system (4.4) reduces to: Form=−M∕2,…,0,1,…,M∕2 ,findsuch that

where

Next we construct a set of basis functions forX0N.

We define,fori=0,1,...,N−2,

and the basis function at the interface is

Then,

Then,writing

It is clear that (4.5) will reduce to a sequence of one-dimensional problems inwhich can be easily solved by a Legendre-spectral method.

5 Numerical Results

To examine the correctness and accuracy of the proposed numerical schemes,we consider the following non-homogeneous problem:

with the boundary condition

We set the exact solution to be

The functionsf,g,hcan then be computed accordingly.

We employ the Fourier–Legendre method presented in the last section and choose (M,N) large enough so that the errors are dominated by the time discretization.In the following examples,we choose=0.5 ,which is a preselected parameter introduced in (3.9a) and (3.14a).

In Fig.2,we plot theL2-errors for the pressure and for the velocity and displacement with the second-order standard and rotational pressure-correction schemes.We observe that the rotational scheme performs much better than the standard scheme.

In Fig.3,we plot the convergence rate of the second-order rotational scheme.We consider ending timeT=2 and vary the step size from Δt=0.1 to Δt=0.000 1 .We observe that theL2errors for the fluid velocity,the structure displacement and the pressure all converge at a rate close to 3 / 2.Due to the Dirichlet boundary condition used for the pressure at the interface,the second-order rotational scheme does not achieve full second-order accuracy for the velocity.This is consistent with the error estimates derived and convergence rates observed for Stokes equations with open boundary in [13].

Next,we examine the energy stability of our schemes by solving the homogeneous (withf,gandhbeing zero) FSI problem with the same initial conditions as in the last example.

Fig.2 Δt=0.01 ,T=2; second-order scheme; top: standard; bottom: rotational

Fig.3 L2 Error for second-order rotational scheme

We take the second-order rotational scheme as an example and plot in Fig.4 the discrete energy for the cases with Δt=0.01 and Δt=0.05 .We observe that the discrete energy indeed decays monotonically.

Fig.4 Energy decay for time step 0.01 and 0.05

6 Conclusions

We constructed in this paper standard and rotational pressure correction schemes for the FSI problem with a fixed interface and proved rigorously that they are unconditionally energy stable.These schemes are new and fundamentally diff erent from existing schemes for the FSI problem.Besides the unconditional stability,they are also computationally very efficient: at each time step,they lead to (i) a coupled linear elliptic system for the velocity and displacement,with the coupling condition at the interface between the fluid and solid regions,which can be efficiently solved by using a standard domain decomposition (with two domains) approach; and (ii) a discrete Poisson equation in the fluid region.

We validated these schemes by using a Fourier–Legendre spatial discretization for the FSI problem in a periodic channel.In particular,our numerical results indicate that the convergence rates of the second-order rotational scheme for the velocity,pressure and displacement inL2-norm are close to 3/2-order.

Although we only considered the FSI problem with fixed interface,we believe that the essential approaches used here in constructing our numerical schemes can be extended to the FSI problem with moving interface [21],which we plan to address in a future endeavor.

References

1.Badia,S.,Codina,R.: On some fluid-structure iterative algorithms using pressure segregation methods.Application to aeroelasticity.Int.J.Numer.Methods Eng.72(1),46–71 (2007)

2.Badia,S.,Nobile,F.,Vergara,C.: Fluid-structure partitioned procedures based on Robin transmission conditions.J.Comput.Phys.227(14),7027–7051 (2008)

3.Badia,S.,Quaini,A.,Quarteroni,A.: Splitting methods based on algebraic factorization for fluidstructure interaction.SIAM J.Sci.Comput.30(4),1778–1805 (2008)

4.Burman,E.,Fernández,M.A.: Stabilized explicit coupling for fluid-structure interaction using Nitsche’s method.C.R.Math.Acad.Sci.Paris 345(8),467–472 (2007)

5.Causin,P.,Gerbeau,J.F.,Nobile,F.: Added-mass eff ect in the design of partitioned algorithms for fluid-structure problems.Comput.Methods Appl.Mech.Eng.194(42–44),4506–4527 (2005)

6.Chakrabarti,S.K.,Hernandez,S.,Brebbia,C.A.: Fluid structure interaction and moving boundary problems,volume 43 of Advances in Fluid Mechanics.WIT Press,Southampton,2005.Edited papers from the 3rd International Conference on Fluid Structure Interaction and the 8th International Conference on Computational Modelling and Experimental Measurements of Free and Moving Boundary Problems held in La Coruna,September 19–21 (2005)

7.Chorin,A.J.: Numerical solution of the Navier–Stokes equations.Math.Comp.22,745–762 (1968)

8.Dong,S.,Shen,J.: An unconditionally stable rotational velocity-correction scheme for incompressible flows.J.Comput.Phys.229(19),7013–7029 (2010)

9.Farhat,C.,Lesoinne,M.,LeTallec,P.: Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: momentum and energy conservation,optimal discretization and application to aeroelasticity.Comput.Methods Appl.Mech.Eng.157(1–2),95–114 (1998)

10.Felippa,C.A.,Park,K.,Farhat,C.: Partitioned analysis of coupled mechanical systems.Comput.Methods Appl.Mech.Eng.190,3247–3270 (2001)

11.Fernández,M .A.: Coupling schemes for incompressible fluid-structure interaction: implicit,semiimplicit and explicit.SeMA J 55,59–108 (2011)

12.Fernández,M.A.,Gerbeau,J.-F.,Grandmont,C.: A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid.Int.J.Numer.Methods Eng.69(4),794–821 (2007)

13.Guermond,J.L.,Minev,P.,Shen,J.: Error analysis of pressure-correction schemes for the timedependent Stokes equations with open boundary conditions.SIAM J.Numer.Anal.43(1),239–258 (2005).(electronic)

14.Guermond,J.L.,Minev,P.,Shen,J.: An overview of projection methods for incompressible flows.Comput.Methods Appl.Mech.Eng.195(44–47),6011–6045 (2006)

15.Guermond,J.L.,Shen,J.: On the error estimates for the rotational pressure-correction projection methods.Math.Comput.73(248),1719–1737 (2004).(electronic)

16.He,Y.,Nicholls,D.P.,Shen,J.: An efficient and stable spectral method for electromagnetic scattering from a layered periodic structure.J.Comput.Phys.231(8),3007–3022 (2012)

17.He,Y.,Shen,J.: Unconditionally stable pressure-correction schemes for a linear fluid-structure interaction problem.Numer.Math.Theory Methods Appl.7(4),537–554 (2014)

18.Hou,G.,Wang,J.,Layton,A.: Numerical methods for fluid-structure interaction–a review.Commun.Comput.Phys.12(2),337–377 (2012)

19.Hron,J.,Turek,S.: A Monolithic FEM/Multigrid Solver for an ALE Formulation of Fluid-Structure Interaction with Applications in Biomechanics.Springer,Berlin,Heidelberg (2006)

20.Hübner,B.,Walhorn,E.,Dinkler,D.: A monolithic approach to fluid-structure interaction using spacetime finite elements.Comput.Methods Appl.Mech.Eng.193(23),2087–2104 (2004)

21.Ignatova,M.,Kukavica,I.,Lasiecka,I.,Tuff aha,A.: On well-posedness for a free boundary fluidstructure model.J.Math.Phys.53(11),115624 (2012)

22.Kukavica,I.,Tuff aha,A.,Ziane,M.: Strong solutions to a nonlinear fluid structure interaction system.J.Diff er.Equ.247(5),1452–1478 (2009)

23.Küttler,U.,Wall,W.A.: Fixed-point fluid-structure interaction solvers with dynamic relaxation.Comput.Mech.43(1),61–72 (2008)

24.Matthies,H.G.,Steindorf,J.: Partitioned strong coupling algorithms for fluid-structure interaction.Comput.Struct.81,805–812 (2003)

25.Pyo,J.-H.: Error estimates for the second order semi-discrete stabilized Gauge-Uzawa method for the Navier-Stokes equations.Int.J.Numer.Anal.Model.10(1),24–41 (2013)

26.Quarteroni,A.,Valli,A.: Domain decomposition methods for partial diff erential equations.Numerical Mathematics and Scientific Computation.The Clarendon Press Oxford University Press,New York (1999)

27.Shen,J.: Efficient spectral-Galerkin method.I.Direct solvers of second- and fourth-order equations using Legendre polynomials.SIAM J.Sci.Comput.15(6),1489–1505 (1994)

28.Temam,R.: Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires i.Arch.Rat.Mech.Anal.32,135–153 (1969)

29.Toselli,A.,Widlund,O.: Domain decomposition methods—algorithms and theory.Springer Series in Computational Mathematics,vol.34.Springer,Berlin (2005)