Comput Mech(2008)43:61–72
DOI10.1007/s00466-008-0255-5
ORIGINAL PAPER
Fixed-pointfluid–structure interaction solvers with dynamic relaxation Ulrich Küttler·Wolfgang A.Wall
Received:31October2007/Accepted:20January2008/Published online:22February2008
©Springer-Verlag2008
Abstract Afixed-pointfluid–structure interaction(FSI) solver with dynamic relaxation is revisited.New develop-ments and insights gained in recent years motivated us to prent an FSI solver with simplicity and robustness in a wide range of applications.Particular emphasis is placed on the calculation of the relaxation parameter by both Aitken’s∆2 method and the method of steepest descent.The methods have shown to be crucial ingredients for efficient FSI simu-lations.
Keywords Fluid–structure interaction·Fixed-point solver·Dirichlet–Neumman partitioning·Strong coupling 1Introduction
The development of numerical solvers forfluid–structure interaction(FSI)problems has been an active area of rearch in the last decade.Reliable FSI solvers are demanded in areas as diver as aeroelasticity[8,21],civil engineering[41]or hemodynamics[2,16].And as wide as the possiblefield of application are the requirements the solvers are confronted with:aerodynamics applications can couple a light compress-iblefluid to a stiff aircraft wing)or a light incompressiblefluid to a very light parachute or sail),whereas hemodynamics simulations couple incom-pressiblefluids andflexible structures with very comparable density.Thus there cannot be one FSI solver thatfits all needs. Instead a variety of solution procedures is needed.
A particular interesting class of FSI problems,that is the appropriate model in many areas,is the interaction of U.Küttler·W.A.Wall(B)
Chair of Computational Mechanics,TU Munich,
Boltzmannstr.15,85747Garching,Germany
e-mail:wall@lnm.mw.tum.de incompressiblefluids with structures that undergo large deformations.In that ca bothfields,thefluid and the struc-ture,prent computational challenges on their own.And still the coupling is a nontrivial task,as the large structural deformations have a huge impact on theflu
idfield’s size and the coupled solution.Oftentimes there are sophisticatedfield solvers available that should be reud.This prents a fur-ther constraint to a possible coupling scheme.There are,of cour,monolithic solvers that treat the nonlinear coupled problem in one go[2,18,19,36],however the approaches require access to thefield solver internals and cannot be pur-sued with black box solvers.The same goes for partitioned approaches that mimic the behavior of monolithic solvers[7].
One possible partitioning strategy that enables the reu of existingfield solvers is the Dirichlet–Neumann partition-ing,the predominant partition approach for FSI solvers.The most popular coupling methods arefixed-point methods[27, 30,40]and interface Newton Krylov methods[10,15,16]. Other solvers suggested include a block-Newton solver with finite differenced off-diagonal blocks[25]and solvers bad on vector extrapolation methods[26,37].See also[33,34] for a FSI solver framework that includes partitioned block-iterative and quasi-direct coupling solvers as well as mono-lithic solvers.
The most basic and yet highly efficient approach among this variety of methods is thefixed-point method with dynamic relaxation as suggested in[27,40].It is extremely easy to implement and surprisingly efficient and robust. Hence,in many cas this is the method of choice and it is an especially preferred method to u if a new attempt at FSI or other coupled problems is pursued.Unf
ortunately,however, the method has never been published in a journal but only in two hardly available conference proceedings[27,40]so far. This paperfinally provides the detailed treatment of both the Aitken relaxation method and the relaxation via the method
123
of steepest descent in the context offixed-point iteration FSI solvers.We attempt to give as easy an introduction as possi-ble.Yet this is not meant to be a mere reiteration of the orig-inal work[27].The continuous work on FSI problems and the ever growing number of FSI solvers stimulated a subtle change of perception that is reflected in the prent contri-bution.Nowadays there is much more emphasize on the for-mulation of the nonlinear interface problem that is amenable to many nonlinear solution approaches.The nonlinearfield solvers are placed behind the opaque interface equation to keep them parated from the coupling algorithm.
The remainder of this paper is organized as follows:In Sect.2thefield equations and coupling conditions of the FSI problem are prented.The introduction of thefixed-point solver with Dirichlet–Neumann partitioning follows in Sect.3and details on the available relaxation methods are prented in Sect.4.Finally two examples are shown in Sect.5.
2Field equations
The FSI problem domain consists of the nonoverlappingfluid and structural domainsΩF andΩS.Both share a common interfaceΓ.In ca of an incompressiblefluid interacting with an elastic body the velocity u and the pressure p are chon for the unknowns of thefluidfield,whereas the struc-turalfield unknown is the displacement d.
2.1Coupling conditions
Esntial for the interaction of bothfields is the coupling of thefield variables at the interface.Both kinematic and dynamic continuity need to be fulfilled at all times.In the usual ca of non-slip conditions at the interface this amounts to
uΓ=d dΓ
d t
andσSΓ·n=σFΓ·n(1)
where of cour the interface displacement dΓdoes change the interface position xΓ=x0,Γ+dΓrelative
to the start-ing position x0,Γ.Therefore kinematic continuity states that along with the interface velocity uΓof thefluidfield the wholefluid domainΩF changes in time.The dynamic con-tinuity states that the stress equal at the deformed interface where n constitutes the time dependent interface normal.
2.2Structural domain
The structural displacements d are governed by the geomet-rically nonlinear elastodynamics equations
ρS d 2d
d t2=∇·(F·S)+ρS b S inΩS×(0,T),(2)
whereρS and b S reprent the structural density and specific
body force,respectively.The cond Piola-Kirchhoff stress
tensor S is related to the Green-Lagrangian strains via
S=C:E with E=
1
2
F T·F−I
,(3)
where C denotes the material tensor and F=∇d reprents
the deformation gradient.The time dependent problem(2)is
subject to the initial and boundary conditions
d=d0and
d d
d t
=d d0
d t
inΩS at t=0
d=¯d onΓS D,S·n=¯h S onΓS N,
whereΓS D andΓS N denote the Dirichlet and Neumann parti-
tion of the structural boundary,respectively.
鸡蛋饼怎么做好吃2.3Fluid domain
Fluid velocity u and pressure p are governed by the incom-
pressible Navier-Stokes equations
∂u
∂t
+u·∇u−1
ρF
∇·σF=b F inΩF×(0,T),(4)
∇·u=0inΩF×(0,T).(5)
The vectorfield b F denotes the specific body force andρF
the density of thefluid.In ca of FSI simulations thefluid
domainΩF varies in time due to the moving interfaceΓ.
One way to account for the domain changes is to consider
the wholefluid domain to deform continuously,starting with
interface displacement dΓ.That is to prescribe a unique map-
ping
武术基本动作x=ϕ(dΓ,x0,t)(6)
of thefluid domain which matches the interface displace-
ments.This mandates an arbitrary Lagrangian-Eulerian
(ALE)formulation[9,12]of the Navier-Stokes equations
and(4)changes to
∂u
∂t
x0
+c·∇u−1
ρF
∇·σF=b F inΩF×(0,T).(7)
The ALE-convective velocity is given as c=u−u G,with
the domain velocity u G=∂ϕ/∂t.The stress tensor of a
邪恶秀Newtonianfluid is given by
σF=−p I+2µ (u)where (u)=1
2
∇u+∇u T
(8)
denotes the strain rate tensor andµthe viscosity.The kine-
matic viscosity is given byν=µ/ρF.
The partial differential equation(7)is subject to the initial
and boundary conditions
u=u0inΩF at t=0
u=¯u onΓF D,σF·n=¯h F onΓF N.(9)
123
2.4Coupled FSI system
A proper discretization in space and time of equations (2),(5)and (7)and a particular version of the fluid domain map-ping (6)together with discrete versions of the coupling con-ditions (1)lead the following nonlinear system of algebraic equations
F
u n +1,p n +1
,d
G ,n +1 =f F (10)S d n +1 =f S (11)
M
d
G ,n +1
=0(12)C SF d n +1Γ,u n +1Γ,p n +1
=0(13)
C SM d n +1
Γ,d G ,n +1Γ
=0(14)where d G ,n +1constitutes the fluid mesh movement and the coupling C SF between structure and fluid and the coupling C SM between structure and mesh describe conditions to be met at the interface.
It is understood that all operators depend on the current time step n +1.The spatial discretization of all fields is done using finite elements.In detail the operators are:–
Fluid F
u n +1,p n +1,d
G ,n +1 =f F :With one-step-θtime discretization,for simplicity of pre-ntation,the discrete version of Navier-Stokes is
1∆t M
F +θ N F (u n +1)+K F θG
θG
T 0 u n +1p n +1
=
b F +(1−θ)M F ˙u n +1∆t M F u n 0
(15)
where the integration has to be performed on a time vary-ing domain ΩF (t ),that is all terms depend on the mesh deformation d G ,n +1.
Remark 1Additional stabilization terms are required in equation (15).For brevity the terms are not
considered here.See [14]for a throughout discussion of stabilized finite elements for flow calculations.
–
Structure S d
n +1
=f S :Discretized in time using the generalized-αmethod [4]the structural equations read
1−αm β∆t
2M S d n +1
+(1−αf )f int (d n +1)=(1−αf )f n +1ext +αf f n
ext
+M S 1−αm d n +1−αm ˙d n +1−αm −2β2β
¨d n −αf f int (d n ).(16)
–
Mesh M (d G ,n +1)=0:
A simple linear mesh equation might look like this K M d
G ,n +1
=f
M
d n +1Γ
(17)
where f M (d n +1
Γ
)stems only from the prescribed inter-face displacement d n +1
Γand K M is an arbitrary operator that extends the interface displacement to the interior of the fluid domain.
The coupling Eqs.(13)and (14)are couplings that apply at the FSI interface only.There are sophisticated mortar methods available to couple various sorts of meshes,[28].The simplest coupling possible,however,is the coupling of matching meshes where all coupling operators C simplify to diagonal matrices at the FSI interface.For simplicity only matching meshes are considered here.–
Structure to fluid coupling C SF (d n +1Γ,u n +1
Γ)=0:u n +1Γ
=
d n +1Γ−d n Γ∆t
(18)
–
Structure to mesh coupling C SM (d n +1
Γ,d G ,n +1Γ)=0:d G ,n +1Γ=d n +1
Γ
(19)
Remark 2The mesh deformation equation (12)is coupled to the structural equation (11),but not vice versa.So the structural deformation is (of cour)not retained by the mesh equation,but the structure drags the mesh along.
3Dirichlet–Neumann coupling scheme
The system (10)–(14)constitutes a coupled t of nonlinear
algebraic equations to be solved once for each time step.A Dirichlet–Neumann partitioning can be applied to solve this system,where the fluid field becomes the Dirichlet partition
with prescribed interface velocities u n +1
Γ
and the structural field becomes the Neumann partition loaded with interface
forces f F ,n +1
Γ
.This way the field solvers remain independent of each other and black box solvers can be ud.
The particular advantage of a Dirichlet–Neumann cou-pling for FSI problems comes from the deformation of the fluid domain ΩF .The unknown interface displacements d n +1
Γ
determine shape and size of the fluid domain ΩF ,thus the Navier-Stokes equation (15)has to be solved on a domain the size of which depends on the unknown solution.This dif-ficulty of a free boundary value problem aris independent of the particular mesh moving scheme (17),even independent of how the changing fluid domain is described numerically,
123
solvers are unable to solve such free boundary value prob-lems.The Dirichlet–Neumann partitioning,however,cures this problem by prescribing definite interface displacements d n +1
Γ
to the fluid solver.So the Navier-Stokes equation (15)needs to be solved on a domain with prescribed motion.The remaining fluid field unknowns are velocity u n +1and pressure p n +1.3.1Solver coupling
A detailed discussion of the solver coupling in a Dirichlet–Neumann partitioning requires to distinguish between inte-rior degrees of freedom I and degrees of freedom at the coupling interface Γ.With this distinction the structural equation S
d n +1 =f S reads S I I S I ΓS ΓI S ΓΓ
d I d Γ
= f S I
f S Γ
.
(20)
The same applies to the mesh and the fluid equation as well.
However,in ca of the fluid equations all pressure degrees of freedom belong to the domain’s interior and are subsumed in the variable u I in the further notation.Furthermore due to the Dirichlet–Neumann coupling the fluid mesh displace-ment d G ,n +1is known during the fluid solver execution.So the fluid operator turns into
F I I F I Γ
F ΓI F ΓΓ u I u Γ = f F I f F Γ .(21)
With the definitions details of the solver coupling can be given:1.Start with a suitably predicted interface displacement d n +1Γ.
2.
Solve mesh equation
M I I d G ,n +1I
=−M I Γd G ,n +1
Γ(22)
with interface condition d G ,n +1Γ=d n +1
Γ
and calculate resulting grid and interface velocity (e [11,12]).
u
G ,n +1
=
d G ,n +1−d G ,n ∆t and u n +1Γ=d n +1
Γ
−d n Γ∆t
(23)
3.
Solve fluid equation with prescribed interface veloc-ity u n +1
Γ
F I I u n +1I
=f F I −F I Γu n +1Γ(24)
for the interior velocity and pressure values u n +1
I
and calculate coupling forces f F ,n +1Γ=F ΓI u n +1I
+F ΓΓu n +1
Γ(25)
4.
Solve structural equation loaded with coupling forces f F ,n +1Γ S I I S I ΓS ΓI S ΓΓ
d n +1I ˜d n +1Γ
=
f S I
f S
Γ−(1−αf )f F ,n +1Γ−αf f n Γ
(26)
to obtain the structural displacements in the structural
由心而生
fields interior d n +1I
as well as on the interface ˜d
n +1Γ
.To simplify the discussion further the above solution steps
are abbreviated with interface operators that map a given
钢套箱interface displacement d n +1
Γ
to interface forces as follows f F ,n +1Γ=F Γ d n +1Γ and f S ,n +1Γ=S Γ d n +1Γ
(27)where the solution demands equilibrium at the interface.F Γ d n +1Γ =S Γ d n +1
Γ (28)
The fluid interface operator f n +1Γ=F Γ(d n +1Γ)abbreviates steps 2and 3of the above algorithm,whereas the inver
structural interface operator d n +1Γ=S −1Γ(f n +1
Γ)denotes step 4.This way the solver coupling reduces to a single line
˜d n +1Γ=S −1Γ(F Γ(d n +1Γ
)).(29)
It is understood that this line contains the execution of all three field solvers,one after the other.Each field solver carries its own state and depends on its own initial and boundary conditions in addition to the condition given at the coupling interface.
Remark 3With the interface operators (27)defined,the cou-pling can also be stated bad on interface forces ˜f n +1Γ=F Γ(S −1Γ(f n +1Γ
大葱能壮阳吗)).(30)
The resulting algorithm is equivalent to the one bad on interface displacements and will not be discusd further.But e [24]for an application of force bad couplings.3.2Fixed-point coupling algorithm
The Dirichlet–Neumann coupling in some n assumes the structural part to dominate the interaction.If the structure rembles a rigid wall in comparison to the fluid (that is no interaction occurs),the Dirichlet–Neumann assumption is
123
obvious.If,on the other hand,the structure is very light and flexible,fluid forces will have a huge impact on the structure
and the predictor’s guess at the interface displacement d n +1
Γwill most probably not match the result.˜d n +1Γ=d n +1Γ
(31)
In the cas an iterative correction of the interface displace-ment is required.In contrast to a weak coupling scheme that solves the Eq.(29)just once for each time step,iterative or strong coupling schemes iterate (29)until the coupled FSI problem is solved.However,as has been shown in [3,13],a weak coupling Dirichlet–Neumann scheme that couples a flexible structure to an incompressible fluid without subitera-tion will never be unconditionally stable.The pressure has to be implicitly coupled at
least.Thus a strong coupling scheme with iterative correction is mandatory.In the prent contri-bution only fully coupled schemes are considered,however e also [1,29]for iterative coupling methods bad on pres-sure gregation.
The interface operators (29)can be ud to define one FSI cycle inside the coupling iteration ˜d n +1Γ,i +1
=
S −1Γ(F Γ(d n +1
Γ,i ))
(32)
where i indicates the iteration counter.
Remark 4The FSI cycle (32)constitutes one pass in a block
Gauß-Seidel solver for the interface displacements.The
structural solver S −1
Γ
works on the result of the fluid solver F Γ.
In order to define an iterative solver a stopping criteria is required.To this end the interface residual is introduced
r n +1Γ,i +1=˜d n +1Γ,i +1−d n +1Γ,i
(33)
and as convergence criteria the length scaled square norm of the residual is ud.
1√
n eq r n +1Γ,i +1
< (34)For further discussions of the applicable interface norms in FSI calculations e [5].
In order to ensure and accelerate convergence of the iter-ation a relaxation step is needed after each FSI cycle (32)
d n +1Γ,i +1=d n +1Γ,i +ωi r n +1Γ,i +1
(35)
=ωi ˜d n +1Γ,i +1+(1−ωi )d n +1Γ,i
with a variable relaxation parameter ωi .The fixed-point algo-rithm to solve FSI problems consists of the relaxed FSI
cycle (35)with appropriate relaxation parameter and conver-gence criteria (34).Calculation methods for the relaxation parameter ωi will be prented in the next ction.
Remark 5Equation (35)can be reformulated
d n +1Γ,i +1=ωi S −1Γ(F Γ(d n +1Γ,i ))+(1−ωi )d n +1
Γ,i
(36)
=d n +1Γ,i +ωi S −1Γ F Γ(d n +1Γ,i )−S Γ(d n +1Γ,i )
The formulation (36)constitutes a nonstationary Richard-son iteration with the operator F Γ−S Γand the precondi-tioner S −1Γ.Thus it would be possible to approximate S −1
Γ(e [6]),however this mandates to evaluate F Γ−S Γexactly.In particular the evaluation of the structural solver S Γis still required.
Remark 6With the interface residual (33)it is possible to define the interface Jacobian
J Γ=
∂r Γ
∂d Γ(37)and apply Newton’s method [22]to solve the interface sys-tem (28).However,as the interface Jacobian (37)is not easily available,a matrix free Krylov subspace solver [23]has to be applied inside Newton’s method.In that ca the relaxation step (35)is replaced by the solution of
J Γ∆d n +1Γ,i +1=−r n +1
Γ,i +1
(38)
and the update step
d n +1Γ,i +1=d n +1Γ,i +∆d n +1
Γ,i +1.
(39)
For further discussion of Newton Krylov methods in the con-text of FSI problems e for instance [10,16].
Remark 7Multiple evaluations of the FSI cycle (32)within one time step require field solvers that can be ret to calcu-late the same time step veral times.Black box solvers that cannot be ret are in general not suited for strong coupling FSI calculations.
4Relaxation methods
Relaxation of the interface displacements (35)is nothing but
the line arch step of a nonlinear solver [17,22].Con-quently the known solver techniques can be applied here as well.
4.1Fixed relaxation parameter
The simplest and most ineffective method is to choo a fixed parameter ωfor all time steps.The relaxation parameter has to be small enough to keep the iteration from diverging,but as large as possible in order to u as much of the new solu-tion as possible and to avoid unnecessary FSI iterations.The optimal ωvalue is problem specific and not known a priori.Furthermore even the optimal fixed value will lead to more iterations than a suitable dynamic relaxation parameter.
123
4.2Aitken relaxation
FSI problems have already been solved in [27,40]using a
Dirichlet–Neumann partitioned approach combined with a fixed-point solver bad on Aitken’s ∆2method as given by [20].This method has proven to be astonishingly simple and efficient.
The central idea of Aitken’s ∆2method is to u values from two previous iterations to improve the current solu-tion.In the FSI ca two pairs of interface displacements
˜d
Γ,i +1,d Γ,i and ˜d Γ,i +2,d Γ,i +1 are required that fulfill (29).In a scalar ca the improved solution could immedi-ately be given d Γ,i +2=
d Γ,i ˜d
诡辩意思Γ,i +2−˜d Γ,i +1d Γ,i +1d Γ,i −˜d Γ,i +1−d Γ,i +1+˜d Γ,i +2
.
(40)
This is nothing but one step of the cant method.And with d Γ,i +2=d Γ,i +1+ωi +1 ˜d
Γ,i +2−d Γ,i +1
the relaxation factor becomes
ωi +1=d Γ,i −d Γ,i +1
d Γ,i −˜d
Γ,i +1−d Γ,i +1+˜d Γ,i +2(41)
=ωi
d Γ,i −˜d
Γ,i +1d Γ,i −˜d Γ,i +1−d Γ,i +1+˜d Γ,i +2(42)=−ωi
r Γ,i +1
r Γ,i +2−r Γ,i +1(43)
For the vector ca the division by r Γ,i +2−r Γ,i +1is of cour
impossible.Instead [20]suggest to multiply by the vector
inver r Γ,i +2−r Γ,i +1 /
r Γ,i +2−r Γ,i +1 2.
ωi +1=−ωi r Γ,i +1 T r Γ,i +2−r Γ,i +1
r Γ,i +2−r Γ,i +1
2(44)This amounts to a projection of the participating vectors in r Γ,i +2−r Γ,i +1direction and to do the scalar extrapolation with the projected values.
Remark 8The relaxation parameter calculation (44)is valid in every iteration step.This version has already been ud for FSI problems in [27],however a lot of confusion has been caud by citations of [27]that misd the recursion on ωi in (44).The Aitken formula without recursion requires two FSI cycles before one relaxation step (44)is possible.This results in a coupling scheme that needs approximately twice the number of iterations to converge.
Remark 9The relaxation parameter (44)is exact in linear scalar cas.For vector cas things are not that clear,how-ever,the Aitken relaxation parameter works very well for many FSI problems and is extremely cheap to calculate and to implement.
Remark 10Two previous steps are required in (44),thus there is no way to calculate the relaxation pa
rameter after the
first FSI cycle ωn +1
1
.In [20]it is suggested to u the last ωn of the previous time step,however this can sometimes be too large a step already.A better choice is to start with a (problem specific)fixed or at least constrained parameter.ωn +1
1
=max (ωn ,ωmax )(45)
Remark 11It is possible to u more than two history values
to improve the current solution.This idea leads to the known vector extrapolation schemes and will be discusd in a future contribution.
4.3Steepest descent relaxation
The best relaxation parameter possible in (35)is the one that
finds the optimal step length in r n +1
Γ,i +1direction.To find it the existence of a merit function φ(d Γ)is assumed,that is
minimal at the solution d n +1
Γ
and sufficiently smooth,such that the relaxation parameter ωi is given by
ωi =arg min ωi
φ(d n +1Γ,i +ωi r n +1
Γ,i +1).
(46)
And the condition to determine ωi follows
d φd ωi =∂φ(d n +1Γ,i +ωi r n +1
Γ,i +1)Γ
·r n +1Γ,i +1!=0(47)
which leads to
∂φ(d n +1Γ,i +ωi r n +1Γ,i +1)
∂d Γ
=φ (d n +1Γ,i +ωi r n +1
Γ,i +1)!
=0.
(48)
From a Taylor ries expansion of φ(d Γ)
φ(d n +1Γ,i +ωi r n +1Γ,i +1)≈φ(d n +1Γ,i )+ωi φ (d n +1
Γ,i ) T
r n +1Γ,i +1
+ω2i 2
r n +1Γ,i +1 T φ (d n +1Γ,i )r n +1
Γ,i +1(49)
the expression for the optimal ωi is obtained.
ωi =− φ (d n +1Γ,i ) T
r n +1Γ,i +1
r n +1Γ,i +1 T
φ (d n +1Γ,i )r n +1Γ,i +1
(50)
At this point a connection between the merit function φ(d Γ)
and the interface residual (33)is assumed
φ (d n +1Γ,i )=r n +1
Γ,i +1,
(51)
an assumption that constrains the admissible interface resid-uals r n +1Γ,i +1to gradients of a scalar function.A rather vere constraint,that leads to a symmetric interface Jacobian J Γ=
泡发海参的最佳方法∂r Γ(d n +1
Γ,i )
∂d Γ
=φ (d n +1
Γ,i )
(52)
123
