A One-Field Monolithic Fictitious Domain
Method for Fluid-Structure Interactions
Yongxing Wang∗, Peter Jimack, Mark Walkley

arXiv:1612.01586v1 [] 5 Dec 2016

School of Computing, University of Leeds, Leeds, UK, LS2 9JT

Abstract
In this article, we present a one-field monolithic fictitious domain (FD) method
for simulation of general fluid-structure interactions (FSI). “One-field” means
only one velocity field is solved in the whole domain, based upon the use of an
appropriate L2 projection. “Monolithic” means the fluid and solid equations are
solved synchronously (rather than sequentially). We argue that the proposed
method has the same generality and robustness as FD methods with distributed
Lagrange multiplier (DLM) but is significantly more computationally efficient
(because of one-field) whilst being very straightforward to implement. The
method is described in detail, followed by the presentation of multiple computational examples in order to validate it across a wide range of fluid and solid
parameters and interactions.
Keywords: Fluid-Structure interaction, Finite element, Fictitious domain,
Monolithic method, One-field fictitious domain method

1. Introduction

5

10

Numerical simulation of fluid-structure interaction is a computational challenge because of its strong nonlinearity, especially when large deformation is
considered. Based on how to couple the interaction between fluid and solid,
existing numerical methods can be broadly categorized into two approaches:
partitioned/segregated methods and monolithic/fully-coupled methods. Similarly, based on how to handle the mesh, they can also be broadly categorized
into two further approaches: fitted mesh/conforming methods and unfitted/nonconforming mesh methods [1].
A fitted mesh means that the fluid and solid meshes match each other at
the interface, and the nodes on the interface are shared by both the fluid and
the solid, which leads to the fact that each interface node has both a fluid velocity and a solid velocity (or displacement) defined on it. It is apparent that
∗ Corresponding

author
Email address: scywa@leeds.ac.uk (Yongxing Wang)

Preprint submitted to Journal of COMPUT METHOD APPL M

March 23, 2018

15

20

25

30

35

40

45

50

55

the two velocities on each interface node should be consistent. There are typically two methods to handle this: partitioned/segregated methods [2, 3] and
monolithic/fully-coupled methods [4, 5, 6]. The former solve the fluid and solid
equations sequentially and iterate until the velocities become consistent at the
interface. These are more straightforward to implement but can lack robustness
and may fail to converge when there is a significant energy exchange between
the fluid and solid [3]. The latter solve the fluid and solid equations simultaneously and often use a Lagrange Multiplier to weakly enforce the continuity of
velocity on the interface [6]. This has the advantage of achieving accurate and
stable solutions, however the key computational challenge is to efficiently solve
the large systems of nonlinear algebraic equations arising from the fully-coupled
implicit discretization of the fluid and solid equations. Fitted mesh methods can
accurately model wide classes of FSI problems, however maintaining the quality
of the mesh for large solid deformations usually requires a combination of arbitrary Lagrangian-Eulerian (ALE) mesh movement and partial or full remeshing
[7]. These add to the computational expense and, when remeshing occurs, can
lead to loss of conservation properties of the underlying discretization [8].
Unfitted mesh methods use two meshes to represent the fluid and solid separately and these do not generally conform to each other on the interface. In
this case, the definition of the fluid problem may be extended to an augmented
domain which includes the solid domain. Similarly to the fitted case, there are
also two broad approaches to treat the solid domain: partitioned methods and
monolithic methods. On an unfitted mesh, there is no clear boundary for the
solid problem, so it is not easy to enforce the boundary condition and solve the
solid equation. A wide variety of schemes have been proposed to address this
issue, including the Immersed Finite Element Method (IFEM) [9, 10, 11, 12, 13]
and the Fictitious Domain method (FDM) [14, 15, 16, 17, 18]. The IFEM developed from the Immersed Boundary method first introduced by Peskin [19],
and has had great success with applications in bioscience and biomedical fields.
The classical IFEM does not solve solid equations at all. Instead, the solid
equations are arranged on the right-hand side of the fluid equations as an FSI
force, and these modified fluid equations are solved on the augmented domain
(occupied by fluid and solid). There is also the Modified IFEM [13], which
solves the solid equations explicitly and iterates until convergence. The FDM
has a similar spirit to IFEM in that it treats the domain occupied by solid as
a fictitious/artificial fluid whose velocity/displacement is constrained to be the
same as that of the solid. The FDM approach usually uses a distributed Lagrange multiplier (DLM) to enforce the constraint [14, 15, 16, 17] whilst the
IFEM typically uses a pseudo body force which is evaluated from the known
deformation of the solid and introduced into the fluid momentum equation. Reference [14] presents a fractional FD scheme for a rigid body interacting with
the fluid, whilst [15] introduces a fractional step scheme using DLM/FD for
fluid/flexible-body interactions. In the case of monolithic methods, [16] uses
a FD/mortar approach to couple the fluid and structure, but the coupling is
limited to a line (2D) representing the structure. Reference [18] uses a mortar
approach to solve fluid interactions with deformable and rigid bodies, and [17]
2

60

65

70

75

80

85

90

95

also solves a fully-coupled FSI system with hierarchical B-Spline grids. There
are also other monolithic methods based on unfitted meshes [20, 21].
It can be seen that the major methods based on unfitted meshes either avoid
solving the solid equations (IFEM) or solve them with additional variables (two
velocity fields and Lagrange multiplier) in the solid domain. However, physically,
there is only one velocity field in the solid domain. In this article, we follow
the one-field spirit and only solve one velocity variable in the whole/augmented
domain. We shall introduce a one-field FD method which can be categorized as
a monolithic approach using an unfitted mesh.
In the one-field spirit, [22] introduces an Eulerian formulation by remeshing
and [23] presents a 1D model using a one-field FD formulation but does not
discuss how to compute the integrals arising from the two different domains.
There are other similar Eulerian formulations for FSI problems, such as the
eXtended Finite Element Method (XFEM) [24], local modification of elements
[25] and other fully Eulerian formulations [26, 27, 28] that are coupled with either
local adaptivity or ALE methods. However these formulations are not in the
spirit of one-field, usually the velocity of the fluid (including fictitious fluid), the
displacement of the solid and the Lagrange multiplier are solved monolithically,
which are three-field formulations (four fields if the moving mesh is solved for
as well).
The main idea of the method presented here is as follows. (1) One-field
formulation: we first discretize the control equations in time, re-write the solid
equation in the form of a fluid equation (using the velocity as a variable rather
than the displacement) and re-write the solid constitutive equation in the updated coordinate system. (2) L2 projection (isoparametric interpolation): we
then combine the fluid and solid equations and discretize them in an augmented
domain. Finally the multi-physics problem is solved as a single field.
The remaining sections are organized as follows. In section 2, the control
equations and boundary conditions for fluid-structure interactions are introduced. Section 3 presents the weak form of the FSI system based on the augmented fluid domain. Section 4 introduces a splitting scheme after discretization
in time. Section 5 and 6 discusses how to linearize the convection step and diffusion step respectively. In section 7, the overall solution algorithm is presented
after discretization in space, which clarifies one of the main differences of the
proposed numerical scheme. In section 8, numerical examples are described to
validate the proposed method across a wide range of flows and material. Some
remarks and observations are discussed in section 9 and finally a brief summary
is presented in section 10.
2. Governing equations for FSI
In the following context, let
Z
(u, v)ω =

3

uvdω,
ω

(1)

100

105

110

where u and v are functions defined in domain ω. Ωft ⊂ Rd and Ωst ⊂ Rd (with
d = 2 in this article) denote the fluid and solid domain repectively which are time
dependent regions as shown in Figure 1. Ω = Ωft ∪Ωst is a fixed domain and Γt =
∂Ωft ∩Ωst is the moving interface between fluid and solid. All subscripts, such as i,
j, and k, represent spatial dimension. If they are repeated in one term (including
the bracket defined in (1)), it implies summation over the spatial dimension; if
they are not repeated, they take the value from 1 to d. All superscripts are
used to distinguish fluid and solid (f and s respectively), distinguish different
boundaries (ΓD and ΓN ) or represent time step (n). For example, ufi and usi
f
s
denote the velocity components of fluid and solid respectively, σij
and σij
denote
n
s
the stress tensor components of fluid and solid respectively, and (ui ) is a solid
velocity component at time tn . Matrices and vectors are denoted by bold letters.

Figure 1:

Schematic diagram of FSI, Ω = Ωft ∪ Ωst , Γ = ΓD ∪ ΓN .

In our model we assume an incompressible fluid governed by the following
equations in Ωft as shown in Figure 1:
ρf

f
σij
= µf

f
∂σij
Dufi
−
= ρf gi ,
Dt
∂xj

∂ufj
= 0,
∂xj
!
∂ufj
∂ufi
f
+
− pf δij = τij
− pf δij .
∂xj
∂xi

(2)

(3)

(4)

We also assume an incompressible solid that is governed by the following
equations in Ωst as shown in Figure 1:
ρs

s
∂σij
Dusi
−
= ρs gi ,
Dt
∂xj

∂usj
= 0,
∂xj
4

(5)

(6)

s
σij
= µs

115

120



∂xsi ∂xsj
− δij
∂Xk ∂Xk



s
− ps δij = τij
− ps δij .

(7)

f
s
In the above τij
and τij
are the deviatoric stress of the fluid and solid ref
s
spectively, ρ and ρ are the density of the fluid and solid respectively, µf is
the fluid viscosity, and gi is the acceleration due to gravity. Note that (5)-(7)
describe an incompressible neo-Hookean model that is based on [16] and is suitable for large displacements. In this model, µs is the shear modulus and ps is
the pressure of the solid (pf being the fluid pressure in (4)). We denote by xi
the current coordinates ofhthe solid
or fluid, and by Xi the reference coordinates
i
∂xi
D
of the solid, whilst F = ∂Xj is the deformation tensor of the solid and Dt
represents the total derivative of time.
On the interface boundary Γt :

ufi = usi ,

(8)

f s
s s
σij
nj = σij
nj ,

(9)

where nsj denotes the component of outward pointing unit normal, see Figure 1.
Dirichlet and Neumann boundary conditions may be imposed for the fluid:
ufi = ūi

on ΓD ,

f
σij
nj = h̄i

on

(10)

ΓN .

(11)

Finally, initial conditions are typically set as:
ufi

t=0

= usi |t=0 = 0,

(12)

though they may differ from (12).
3. Weak formulation
 f

ui in Ωft
pf
Let ui =
s
s and p =
ui in Ωt
ps
following symbolic operations:

in
in

Ωft
. We then perform the
Ωst

(Eq.(2), vi )Ωf − (Eq.(3), q)Ωf + (Eq.(5), vi )Ωs − (Eq.(6), q)Ωs ,
t

125

t

t

t

for independent test functions vi ∈ H01 (Ω) and q ∈ L2 (Ω).
Integrating the stress terms by parts, using constitutive equations (4) and
(7) and boundary condition (11), gives the following weak form for the FSI
system.

5

Find ui ∈ H 1 (Ω) and p ∈ L20 (Ω) such that








Dui
∂vj
∂uj
f ∂vi
ρf
,
, vi
+ τij
− p,
−
,q
Dt
∂xj Ω
∂xj Ω
∂xj
Ω



Ω

Dui
f ∂vi
s
s
f
+ τij − τij ,
, vi
+ ρ −ρ
Dt
∂xj Ωs
Ωst
t


f
s
f
= h̄i , vi ΓN + ρ (gi , vi )Ω + ρ − ρ (gi , vi )Ωs ,

(13)

t

130

135

f
∀vi ∈ H01 (Ω) and ∀q ∈ L2 (Ω). In the above, ρf and τij
are extended to be

2
2
defined over the whole of Ω, and L0 (Ω) = p : p ∈ L (Ω), p|P0 = 0 , where P0
is a reference point. Note that the integrals on the interface (boundary forces)
are cancelled out using boundary condition (9). This is not surprising because
they are internal forces for the whole FSI system considered here.
f
Remark 1 The fluid deviatoric stress τij
is generally far smaller than the
f
s
solid deviatoric stress τij , so we choose to neglect the fluid deviatoric stress τij
s
in Ω in what follows. Note that the classical IFEM neglects the whole fluid
f
stress σij
when computing the FSI force [9]. An equivalent way of interpreting
f
neglecting τij
in Ωs is to view the solid as being slightly visco-elastic, having
the same viscosity as the fluid.

4. Discretization in time
The integrals in equation (13) are carried out in two different domains as
illustrated in Figure 1. We use an Eulerian mesh to represent Ω and an updated
Lagrangian mesh to represent Ωs , therefore the total time derivatives in these
two different domains have different expressions, i.e:
∂ui
∂ui
Dui
=
+ uj
Dt
∂t
∂xj

in Ω,

and

(14)

Dusi
∂usi
=
in Ωs .
(15)
Dt
∂t
Firstly, based on the above two equations (14) and (15), we discretize (13)
in time using a backward finite difference. Then omitting the superscript n + 1,
showing the solution is at the end of the time step, for convenience, we obtain:








∂ui
∂vj
∂uj
ui − uni
f ∂vi
f
+ uj
, vi
+ τij ,
− p,
−
,q
ρ
∆t
∂xj
∂xj Ω
∂xj Ω
∂xj
Ω
Ω




n

u
−
u
∂v
i
i
i
s
(16)
+ ρs − ρf
, vi
+ τij
,
∆t
∂xj Ωs
Ωsn+1
n+1


= h̄i , vi ΓN + ρf (gi , vi )Ω + ρs − ρf (gi , vi )Ωs .
n+1

6

140

Note that in the above we have replaced Ωstn+1 by Ωsn+1 for notational convenience. Using the splitting method of [29, Chapter 3], equation (16) can be
expressed in the following two steps.
(1) Convection step:

 ∗
∗
ui − uni
∗ ∂ui
f
+ uj
, vi
= 0;
(17)
ρ
∆t
∂xj
Ω
(2) Diffusion step:








ui − u∗i
∂vj
∂uj
f ∂vi
f
ρ
, vi
+ τij ,
− p,
−
,q
∆t
∂xj Ω
∂xj Ω
∂xj
Ω


 Ω n 

∂v
u
−
u
i
i
i
s
, vi
+ τij
,
+ ρs − ρf
∆t
∂xj Ωs
Ωsn+1
n+1


f
s
f
= h̄i , vi ΓN + ρ (gi , vi )Ω + ρ − ρ (gi , vi )Ωs .

(18)

n+1

The treatment of the above two steps is described separately in the following
subsections.
145

150

5. Linearization of the convection step
In this section, two methods are introduced to treat the convection equation:
Least-squares method and Taylor-Galerkin method, both of which can be used in
the framework of the proposed scheme. Some numerical results for comparison
between these two methods are discussed subsequently in section 5. Because
the overall scheme is explicit, all non-linear terms are linearized using the values
from the last time step. Of course, the scheme can be made implicit with the
same linearized form by iterating within each time step starting from the value
at the last time step.
5.1. Least-squares method
Linearization of equation (17) gives,


 


∂un
∂u∗
∂un
u∗i + ∆t u∗j i + unj i , vi
= uni + ∆tunj i , vi
.
∂xj
∂xj
∂xj
Ω
Ω

(19)

For Least-squares method [30], we may choose the test function in the following
form:


∂uni
n ∂wi
vi = L (wi ) = wi + ∆t wj
+ uj
,
(20)
∂xj
∂xj Ω
where wi ∈ H01 (Ω). In such a case, the weak form of (17) is:


∂un
(L (u∗i ) , L (wi ))Ω = uni + ∆tunj i , L (wi )
.
∂xj
Ω
155

(21)

In our method a standard biquadratic finite element space is used to discretize
equation (21) directly.
7

5.2. Taylor-Galerkin method
It is also possible to linearize equation (17) as:
 ∗

ui − uni
1 n ∂
∗
n
+ uj
(u + ui ) , vi
=0
∆t
2 ∂xj i
Ω
or



∂un
u∗i − uni
+ unj i , vi
∆t
∂xj

(22)


= 0.

(23)

Ω

Rewriting (23) as
u∗i = uni − ∆tunj

∂uni
,
∂xj

(24)

substituting (24) into equation (22), and applying integration by parts we obtain:
 ∗



ui − uni
∂un
∆t
∂un
∂vi
+ unj i , vi
,
(25)
=−
unk i , unj
∆t
∂xj
2
∂xk
∂xj Ω
Ω
where the boundary integral is neglected because uni is the solution of the previous diffusion step, which means no convection exists on the boundary after the
diffusion step. Finally the weak form of Taylor-Galerkin method [29, Chapter
2] can be expressed, by rearranging the last equation, as:




∂un
∂un
∆t2
∂vi
(u∗i , vi )Ω = uni − ∆tunj i , vi
unk i , unj
.
(26)
−
∂xj
2
∂xk
∂xj Ω
Ω
6. Linearization of the diffusion step
As mentioned above, the overall scheme is explicit, so all the derivatives are
n
computed on the known coordinate (xsi ) (denoted xni for convenience). One
n+1
could also construct xi
at each time step and take derivatives with respect to
xn+1
, however we do not consider such an approach in this article. According
i
s
to the definition of τij
in equation (7),
!

∂xn+1
∂xn+1
j
s n+1
s
i
τij
=µ
− δij .
(27)
∂Xk ∂Xk
The last equation, using a chain rule, can also be expressed as:
!

∂xn+1
∂xn+1
j
s n+1
s
i
− δij
τij
=µ
∂xnk ∂xnk
 n+1 ,
n+1 
∂xj
∂xnk ∂xnl
s ∂xi
− δkl
+µ
n
∂xk
∂Xm ∂Xm
∂xnl

s n+1
can be expressed by coordinate xni as follows:
and then τij
!

∂xn+1
∂xn+1
∂xn+1
∂xn+1
j
j
s n+1
s
s n
i
i
τij
=µ
−
δ
+
(τ
)
.
ij
kl
∂xnk ∂xnk
∂xnk
∂xnl
8

(28)

(29)

Using xn+1
− xni = un+1
∆t, which is the displacement at the current step, the
i
i
last equation may be expressed as:
!


∂un+1
∂un+1
∂un+1
∂un+1
j
j
s
s n+1
s n
i
i
= µ ∆t
+
+ ∆t
τij
+ τij
∂xnj
∂xni
∂xnk ∂xnk
(30)
n+1

∂ujn+1
∂un+1
∂un+1
j
s n
s n
2 ∂ui
s n
i
(τkl )
+ ∆t
.
+ ∆t
τkj + ∆t (τil )
∂xnk
∂xnl
∂xnk
∂xnl
Finally, after linearization of the last equation, the weak form (18) can be expressed as:




 usi − (usi )n
ui − u∗i
ρf
, vi
+ ρs − ρf
, vi
∆t
∆t
Ω
Ωsn+1






∂vj
∂uj
∂ui
∂uj ∂vi
f
− p,
−
+µ
+
,
,q
∂xj
∂xi ∂xj Ω
∂xj Ω
∂xj
Ω


n
n
∂u
∂u
∂u
∂u
∂u
∂u
∂v
i
j
i
j
i
j
+ µs ∆t
+
+ ∆t
+ ∆t i
,
∂xj
∂xi
∂xk ∂xk
∂xk ∂xk ∂xj Ωs
n+1


n
n
∂u
∂u
∂u
∂u
∂v
i
j
i
j
n
n
i
(31)
+ ∆t2
(τ s )
+
(τ s )
,
∂xk kl ∂xl
∂xk kl ∂xl ∂xj Ωs
n+1


∂u
∂v
∂ui s n
j
i
n
τ
+ (τils )
,
+ ∆t
∂xk kj
∂xl ∂xj Ωs
n+1


f
s
f
= h̄i , vi ΓN + ρ (gi , vi )Ω + ρ − ρ (gi , vi )Ωs
n+1


n
n ∂un
n

∂u
∂u
∂u
j
j
i
i
s
2
2
s n
s n ∂vi
+ µ ∆t
+ ∆t
(τ )
− τij ,
.
∂xk ∂xk
∂xk kl ∂xl
∂xj Ωs
n+1

160

The spatial discretization of the above linearized weak form will be discussed
in the following section, along with the overall solution algorithm.
7. Discretization in space and solution algorithm
7.1. Spatial discretization
We shall use a fixed Eulerian mesh for Ω and an updated Lagrangian mesh
for Ωsn+1 to discretize equation (31). First, we discretize Ω as Ωh using P2 P1 elements (the Taylor-Hood element) with the corresponding finite element spaces
as
V h (Ωh ) = span {ϕ1 , · · · , ϕN u } ⊂ H 1 (Ω)
and
Lh (Ωh ) = span {φ1 , · · · , φN p } ⊂ L2 (Ω) .

9

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

h

u (x) =

N
X

p

h

u(xi )ϕi (x),

p (x) =

i=1

N
X

p(xi )φi (x).

(32)

i=1

s
We further discretize Ωsn+1 as Ωsh
n+1 (actually it is discretized once on Ω0
and then updated from the previous mesh) using P1 elements (bilinear triangle
element) with the corresponding finite element spaces as:

1
s
s
Ωsn+1 ,
V sh (Ωsh
n+1 ) = span {ϕ1 , · · · , ϕN s } ⊂ H

and approximate uh (x)

x∈Ωsh
n+1

as:

s

sh

u

(x) =

N
X

s

h

u

(xsi )ϕsi (x)

=

i=1

165

u

N X
N
X

u(xj )ϕj (xsi )ϕsi (x),

(33)

i=1 j=1

where xsi is the nodal coordinate of the solid mesh. Notice that the above

approximation defines an L2 projection Pn+1 from V h to V sh : Pn+1 uh (x) =
ush (x).
Substituting (32), (33) and similar expressions for the test functions vh , q h
and vsh into equation (31) gives the following matrix form:

   
A B u
b
=
,
(34)
p
0
BT 0
where
A = M/∆t + K + DT (Ms /∆t + Ks ) D,

(35)

b = f + DT f s + Mu∗ /∆t + DT Ms Dun /∆t.

(36)

and
In the above, matrix D is the isoparametric interpolation matrix derived
from equation (33) which can be expressed as
 T

P
0
D=
, Pij = ϕi (xsj ).
0 PT
All the other matrices and vectors arise from standard FEM discretization: M
and Ms are mass matrices from discretization of integrals in Ωh (with shape
function ϕi ) and Ωsh (with shape function ϕsi ) respectively, and similarly
for

∂v

stiffness matrices K and Ks . B is from discretization of integral − p, ∂vjj
170

s

in

(31). The force vectors f and f come from discretization of integrals on the
right-hand side of (31) in Ωh and Ωsh respectively. The specific expressions of
these matrices and vectors can be found in Appendix A.

10

7.2. Overall solution algorithm
175

180

185

190

Having derived a discrete system of equations we now describe the solution
algorithm at each time step.
 f n
u
in Ωf
s n
n
(1) Given the solid configuration (x ) and velocity field u =
s n
(u )
in Ωs
at time step n.
(2) Discretize the convection equation (21) or (26) and solve it to get an intermediate velocity u∗ .
(3) Compute the interpolation matrix and solve equation (34) using u∗ and
n
(us ) as initial values to get velocity field un+1 .
n+1
(4) Compute solid velocity (us )
= Dun+1 and update the solid mesh by
n+1
n
n+1
(xs )
= (xs ) + ∆t (us )
, then go to step (1) for the next time step.
Remark 2 The choice of P1 element for an updated domain Ωs is convenient,
because the form of the bilinear shape functions stays the same when updating
n+1
n
n+1
the nodal coordinates using (xs )
= (xs ) + ∆t (us )
.
Remark 3 When implementing the algorithm, it is unnecessary to perform
the matrix multiplication DT Ks D globally, because the FEM interpolation is
locally based. All the matrix operations can be computed based on the local
element matrices only. Alternatively, if an iterative solver is used, it is actually
T s
unnecessary
 to compute D K D. What an iterative step needs is to compute
T s
D K D u for a given vector u, therefore one can compute Du first, then
Ks (Du), and last DT (Ks Du).
8. Numerical experiments

195

200

In this section, we present some numerical examples that have been selected
to allow us to assess the accuracy and the versatility of our proposed numerical scheme. We demonstrate convergence in time and space, furthermore, we
favorably compare results with those obtained using monolithic approaches and
IFEM, as well as compare against results from laboratory experiment.
In order to improve the computational efficiency, an adaptive spatial mesh
with hanging nodes is used in all the following numerical experiments. Readers
can reference [31, 32, 33, 34] for details of the treatment of hanging nodes. The
Least-squares method (section 5.1) is used to treat the convection step in all
tests unless stated otherwise.

205

8.1. Oscillation of a flexible leaflet oriented across the flow direction

210

This numerical example is used by [15, 16, 17] to validate their methods. We
first use the same parameters as used in the above three publications in order
to compare results and test convergence in time and space. We then use a wide
range of parameters to show the robustness of our method. The computational
domain and boundary conditions are illustrated in Figure 2.

11

Figure 2:

Computational domain and boundary conditions, taken from [16].

Fluid
L = 4.0 m
H = 1.0 m
ρf = 100 kg/ m3
µf = 10 N · s/ m2

Leaflet
w = 0.0212 m
h = 0.8 m
ρs = 100 kg/ m3
µs = 107 N / m2

Table 1: Properties and domain size for test problem 8.1
with a leaflet oriented across the flow direction.

215

220

225

230

The inlet flow is in the x-direction and given by ux = 15.0y (2 − y) sin (2πt).
Gravity is not considered in the first test (i.e. g = 0), and other fluid and solid
properties are presented in Table 1.
The leaflet is approximated with 1200 linear triangles with 794 nodes (medium
mesh size), and the corresponding fluid mesh is adaptive in the vicinity of the
leaflet so that it has a similar size. A stable time step ∆t = 5.0 × 10−4 s is
used in these initial simulations. The configuration of the leaflet is illustrated
at different times in Figure 3.
Previously published numerical results are qualitatively similar to those in
Figure 3 but show some quantitative variations. For example, [16] solved a
fully-coupled system but the coupling is limited to a line, and the solid in their
results (Figure 7 (l)) behaves as if it is slightly harder. Alternatively, [15] used
a fractional step scheme to solve the FSI equations combined with a penalty
method to enforce the incompressibility condition. In their results (Fig. 3 (h))
the leaflet behaves as if it is slightly softer than [16] and harder than [17]. In [17]
a beam formulation is used to describe the solid. The fluid mesh is locally refined
using hierarchical B-Splines, and the FSI equation is solved monolithically. The
leaflet in their results (Fig. 34) behaves as softer than the other two considered
here. Our results in Figure 3 are most similar to those of [17]. This may be
seen more precisely by inspection of the graphs of the oscillatory motion of the
leaflet tip in Figure 4, corresponding to Fig. 32 in [17]. We point out here
that Taylor-Galerkin method has also been used to solve the convection step
for this test, and we gain almost the same accuracy using the same time step

12

(a) t = 0.1s

(b) t = 0.2s

(c) t = 0.6s

(d) t = 0.8s

Figure 3:

235

240

Configuration of leaflet and magnitude of velocity on the adaptive fluid mesh.

∆t = 5.0 × 10−4 s. Having validated our results for this example against the
work of others, we shall use this test case to further explore more details of our
method.
We commence by testing the influence of the ratio of fluid and solid mesh
sizes rm =(local fluid element area)/(solid element area). Fixing the fluid mesh
size, three different solid mesh sizes are chosen: coarse (640 linear triangles with
403 nodes rm ≈ 1.5), medium (1200 linear triangles with 794 nodes rm ≈ 3.0)
and fine (2560 linear triangles with 1445 nodes rm ≈ 5.0), and a stable time
13

Figure 4:

245

Evolution of horizontal and vertical displacement at top right corner of the leaflet.

step ∆t = 5.0 × 10−4 s is used. From these tests we observe that there is a slight
difference in the solid configuration for different meshes, as illustrated at t = 0.6s
in Figure 5. Significantly however, the difference in displacement decreases as
the solid mesh becomes finer. Further, we found that 1.5 ≤ rm ≤ 5.0 ensures
the stability of the proposed approach. Note that we use a 9-node quadrilateral
for the fluid velocity and 3-node triangle for solid velocity, so rm ≈ 3.0 means
the fluid and solid mesh locally have a similar number of nodes for velocity.

(a)coarse

(b)medium

(c)fine

Figure 5:

Configuration of leaflet for different mesh ratio rm ,
and contour plots of displacement magnitude at t = 0.6s.

250

255

We next consider convergence tests undertaken for refinement of both the
fluid and solid meshes with the fixed ratio of mesh sizes rm ≈ 3.0. Four different
levels of meshes are used, the solid meshes are: coarse (584 linear triangles
with 386 nodes), medium (1200 linear triangles with 794 nodes), fine (2560
linear triangles with 1445 nodes), and very fine (3780 linear triangles with 2085
nodes). The fluid meshes have the corresponding sizes with the solid at their
maximum refinement level. As can be seen in Figure 6 and Table 2, the velocity
is converging as the mesh becomes finer.
In addition, we consider tests of convergence in time using a fixed ratio of
fluid and solid mesh sizes rm ≈ 3.0. Using the medium solid mesh size and the
same fluid mesh size as above, results are shown in Figure 7 and Table 3. From
14

(a) Coarse

(b) Medium

(c) Fine

(d) Very fine

Figure 6:

Contour plots of horizontal velocity at t = 0.5s.

Between different mesh sizes
coarse and medium
medium and fine
fine and very fine

Difference of maximum
horizontal velocity at t = 0.5s
0.01497
0.00214
0.00190

Table 2: Comparison of maximum velocity for different meshes.

260

these it can be seen that the velocities are converging as the time step decreases.
Steps sizes compared
∆t = 2.0 × 10−3 and ∆t = 1.0 × 10−3
∆t = 1.0 × 10−3 and ∆t = 5.0 × 10−4
∆t = 5.0 × 10−4 and ∆t = 2.5 × 10−4

Difference of maximum
horizontal velocity at t = 0.5s
0.00854
0.00517
0.00263

Table 3: Comparison of maximum velocity for different time step size.
Finally, in order to assess the robustness of our approach, we vary each of
the physical parameters using three different cases as shown in Figure 8. A
medium mesh size with fixed rm ≈ 3.0 is used to undertake all of theses tests.
The dimensionless parameters shown in Figure 8 are defined as: ρr = ρρf , µ̄s =
265

µs
, Re
ρf U 2

f

= ρ µUfH and F r = gH
U 2 where the average velocity U = 10 in this
example. The period of inlet flow is T = 1.
It can be seen from the results of group (a) that the larger the value of
shear modulus µ̄s the harder the solid behaves, however a smaller time step is
15

(a) ∆t = 2.0 × 10−3 s
(breaks down at t = 0.61s).

(b) ∆t = 1.0 × 10−3 s.

(c) ∆t = 5.0 × 10−4 s.

(d) ∆t = 2.5 × 10−4 s.

Figure 7:

(a) ρr = 1, Re = 100 and F r = 0.

(b) ρr = 1, µ̄s = 103 and F r = 0.

(c) Re = 100, µ̄s = 103 and F r = 0.

(d) Re = 100 and µ̄s = 103 .

Figure 8:

270

Contour plots of horizontal velocity at t = 0.5s.

Parameters sets and results, ∆t = 5.0 × 10−4 s for Group (b)∼(d).

required. For the case of µ̄s = 109 , the solid behaves almost like a rigid body,
as we would expect. From the results of group (b) it is clear that the Reynolds
Number (Re) has a large influence on the behavior of the solid. The density
and gravity have relatively less influence on the behavior of solid in this problem
which can be seen from the results of group (c) and group (d) respectively.
16

8.2. Oscillating disc surrounded by fluid
This example is taken from [35] and used to validate the conservation of
mass and energy using the proposed method. The computational domain is a
square [0, 1] × [0, 1] with a homogeneous Dirichlet boundary conditions imposed
for velocity, whilst the pressure is fixed to be zero at the left-bottom corner of
the square. A soft solid disc is initially located in the middle of the square and
has a radius of 0.2. The initial velocity of the fluid and solid are prescribed by
the following stream function
Ψ = Ψ0 sin(ax)sin(by),
275

−2

where Ψ0 = 5.0 × 10 and a = b = 2π. The whole system then evolves from
this initial condition.
We first use the same parameters as used in [35]: ρf = ρs = 1.0, µf = 10−3
and µs = 1.0, then we swith to ρs = 2.0 and ρs = 10, and undertake the
tests with the other parameters unchanged. Three different initial meshes for
the square (coarse: 20 × 20, medium: 40 × 40 and fine: 80 × 80) and 2-level
adaptive refinement near the interface based on these initial meshes are used.
The solid mesh has a similar number of nodes with the fluid mesh near the
interface. A snapshot of the velocity norm on the medium adpative mesh, and
the corresponding deformation of the disc at t = 0.5, are plotted in Figure 9.
The energy of this FSI system is computed as follows.

(b) Defomation of solid.

(a) Velocity norm.

Figure 9: Velocity norm and solid defomation at t = 0.5
for ρs = 1.0, ∆t = 10−3 (medium mesh).
Kinetic energy in Ω:
Ek (Ω) =

ρf
2

Z

2

|u| .

(37)

Ω

Kinetic energy in Ωst :
Ek (Ωst ) =

ρs − ρf
2
17

Z

2

|u| .
Ωst

(38)

Viscous dissipation in Ω:
Z tZ

∂ui
.
∂xj

(39)

(trFFT − d) .

(40)

f
τij

Ed (Ω) =
0

Ω

Potential energy of solid:
Ep (Ωs0 ) =

µs
2

Z
Ωs0

Analytically, the total energy
E = Ek (Ω) + Ek (Ωst ) + Ed (Ω) + Ep (Ωs0 )

280

should be a constant. This is considered in the plots of Figure 10. When the
fluid and solid have the same density, the maximum variation of total energy is
around 1.6% (t = 0.26), as shown Figure 10 (a), and when their densities are
different (ρs = 2.0), as shown in 10 (b), the maximum variation of total energy
is around 2.2% at t = 0.31. For the case of ρs = 10 we have a similar result,
with the maximum variation of total energy being about 4.9% at t = 0.6 using
the same time step ∆t = 10−3 .

(b) ρf = 1.0, ρs = 2.0.

(a) ρf = ρs = 1.0.

Figure 10:

285

290

(41)

Energy evolution of an oscillating disc, ∆t = 10−3 (medium mesh).

We further verify the convergence of both energy and mass, which is clearly
demonstrated in Figure 11.
8.3. Oscillation of a flexible leaflet oriented along the flow direction
The following test problem is taken from [36], which describes an implementation on an ALE fitted mesh. It has since been used as a benchmark to validate
different numerical schemes [17, 18]. The geometry and boundary conditions are
shown in Figure 12.
For the fluid, the viscosity and density are µf = 1.82 × 10−4 and ρf =
1.18 × 10−3 respectively. For the solid, we use shear modulus µs = 9.2593 × 105
18

(a) Engergy convergence.

Figure 11:

Figure 12:

(b) Mass convergence.

Convergence of energy and mass, ρs = 2.0, ∆t = 5.0 × 10−4 .

Computational domain and boundary condition for oscillation of flexible leaflet.

(a) Leaflet displacement and fluid pressure.

Figure 13:

295

(b) Mesh refinement near the structure.

Contour plots of leaflet displacement and fluid pressure at t = 5.44s.

and density ρs = 0.1. The leaflet is discretized by 1063 3-node linear triangles
with 666 nodes, and the corresponding fluid mesh locally has a similar node
density to the leaflet (rm ≈ 3.0). First the Least-squares method is tested and a

19

Figure 14:

300

305

310

315

320

325

Distribution of pressure across the leaflet on the three lines in Figure 13 (b).

stable time step ∆t = 1.0 × 10−3 s is used. Snapshots of the leaflet deformation
and fluid pressure at t = 5.44s are illustrated in Figure 13. In Figure 14, the
distributions of pressure across the leaflet corresponding to the three lines (AB,
CD and EF) in Figure 13 (b) are plotted, from which we can observe that the
sharp jumps of pressure across the leaflet are captured.
The evolution of the vertical displacement of the leaflet tip with respect to
time is plotted in Figure 15(a). Both the magnitude (1.34) and the frequency
(2.94) have a good agreement with the result of [36], using a fitted ALE mesh
and of [17], using a monolithic unfitted mesh approach. Taylor-Galerkin method
is also tested using ∆t = 2.0 × 10−4 s as a stable time step, and a corresponding
result is shown in 15(b). This shows a similar magnitude (1.24) and frequency
(2.86). These results are all within the range of values in [17, Table 4]. Note
that since the initial condition before oscillation for these simulations is an
unstable equilibrium, the first perturbation from this regime is due to numerical
disturbances. Consequently, the initial transient regimes observed for the two
methods (Least-squares and Taylor-Galerkin methods) are quite different. It
is possible that an explicit method causes these numerical perturbations more
easily, therefore makes the leaflet start to oscillate at an earlier stage than when
using Least-squares approach.
8.4. Solid disc in a cavity flow
This numerical example is used to compare our method with the IFEM,
which is described in [11, 35]. In order to compare in detail, we also implement
the IFEM, but we implemented it on an adaptive mesh with hanging nodes, and
we use the isoparametric FEM interpolation function rather than the discretized
delta function or RKPM function of [9, 10].
The fluid and solid properties are chosen to be the same as in [35]: ρf =
s
ρ = 1.0, µf = 0.01 and µs = 0.1. The horizontal velocity on the top boundary
of the cavity is prescribed as 1 and the vertical velocity is fixed to be 0 as shown
in Figure 16. The velocities on the other three boundaries are all fixed to be 0,
and pressure at the bottom-left point is fixed to be 0 as a reference point.

20

(a) Least-squares method.

Figure 15:

(b) Taylor-Galerkin method.

Displacement of leaflet tip as a function of time.

Figure 16:

Figure 17:

Computational domain for
cavity flow, taken from [35].

330

335

340

Adaptive mesh for
cavity flow.

In order to compare our method with IFEM we use the same meshes for fluid
and solid: the solid mesh has 2381 nodes and the fluid mesh locally has a similar
number of nodes (adaptive, see Figure 17). First the Least-squares method is
used to solve the convection step, and the time step is ∆t = 1.0 × 10−3 . Figure
18 shows the configuration of the deformed disc at different stages, from which
we do not observe significant differences of the velocity norm even for a long
run as shown in Figure 18 (b). Then Taylor-Galerkin method is tested, and we
achieve almost the same accuracy by using the same time step (not shown in
the Figure).
We also test different densities, and the cases of µs = 1.0 and µs = 100. For
our proposed method we can use µs = 100 or larger in order to make the solid
behave like a rigid body without changing time step (again, not shown here due
to lack of space). This is not possible for the IFEM for which the simulation
always breaks down for µs = 100, however small the time step, due to the huge
FSI force on the right-hand side of the IFEM system.

21

(a) t = 4.5s.

(b) t = 25.5s.

Figure 18:

345

Velocity norm for a soft solid (µs = 0.1) in a driven cavity flow
using our method (left) and IFEM (right).

8.5. Falling disc in a channel with gravity
The final test that we present in this paper is that of a falling disc in a
channel, as cited by [10, 18] for example, in order to further validate against the
IFEM and a monolithic method respectively. The computational domain and
parameters are illustrated in Figure 19 (a) and Table 4 respectively. The fluid
velocity is fixed to be 0 on all boundaries except the top one.

(a) Computational domain

Figure 19:

(b) Contour of vertical velocity at t = 1s (fine mesh)
Falling disc in a channel with gravity.

22

Fluid
W = 2.0 cm
H = 4.0 cm
ρf = 1.0 g/ cm3
µf = 1.0 dyne · s/ cm2
g = 980 cm/ s2

Disc
d = 0.125 cm
h = 0.5 cm
ρs = 1.2 g/ cm3
µs = 108 dyne/ cm2
g = 980 cm/ s2

Table 4: Fluid and material properties of a falling disc.
There is also an empirical solution of a rigid ball falling in a viscous fluid
[18], for which the terminal velocity, ut , under gravity is given by

  
 r 4 
 r 2
ρs − ρf gr2
L
−
1.7302
, (42)
ln
−
0.9157
+
1.7244
ut =
4µf
r
L
L

350

355

360

where ρs and ρf are the density of solid and fluid respectively, µf is viscosity
of the fluid, g = 980 cm/ s2 is acceleration due to gravity, L = W / 2 and r is
the radius of the falling ball. We choose µs = 108 dyne/ cm2 to simulate a
rigid body here, and µs = 1012 dyne/ cm2 is also applied, which gives virtually
identical results.
Three different meshes are used: the disc boundary is represented with 28
nodes (coarse), 48 nodes (medium), or 80 nodes (fine). The fluid mesh near the
solid boundary has the same mesh size as that of the disc, and a stable time step
t = 0.005s is used for all three cases. A local snapshot of the vertical velocity
with the adaptive mesh is shown in Figure 19 (b). From the fluid velocity
pattern around the disc we can observe that the disc behaves like a rigid body
as expected. In addition, the evolution of the velocity of the mid-point of the
disc is shown in Figure 20, from which it can be seen that the numerical solution
converges from below to the empirical solution.
9. Discussion
In this section, some further remarks and notes concerning the proposed
method are discussed.
9.1. Treatment of the convection equation

365

370

Both the Least-squares method and Taylor-Galerkin method add artificial
diffusive terms in their fomulations to stablize the numerical scheme. Like all
such stabilization approaches this necessarily has an influence on the accuracy,
especially for large Reynolds numbers. In such cases a balance is required between minimizing the artificial dissipation and maintaining a stable time step
size that is acceptable. In our applications, the Reynolds number is around
100 ∼ 500, except for two extreme test cases in section 8.1 whose Reynolds
numbers are 1000 and 5000 respectively (Figure 8 (b)). Even then, in these
cases a minimal amount of diffusion is observed provided we use a small time
23

Figure 20: Evolution of velocity at the center of a falling disc.
(The blue solid line represents the empirical solution from formula (42),
for interpretation of the reference to color, the reader is referred to the web version.)

375

380

385

390

395

step (5 × 10−4 ). Alternatively, an upwind scheme or a discontinuous Galerkin
method could be a better choice. However, we have not yet implemented such
methods on the adaptive mesh with hanging nodes.
9.2. The Lagrangian update of the solid
Updating the solid based upon its velocity could lead to distorted elements,
either in its interior or at its boundary. Should this occur there are advanced
mesh update techniques to improve the quality of solid mesh [37] or discrete
remeshing may be used [7]. However all of the tests undertaken in this article
have been performed based upon published benchmarks using incompressible
solids and a small time step, and we have not encountered the problem of
significantly distorted elements. In other applications our simple Lagrangian
approach may not be adequate and so ALE techniques, possibly including mesh
quality improvement, may also be required.
9.3. Contact between solids and boundaries
In many applications moving solids may run into boundaries (either external
or of other moving bodies). In this article, we have only considered standard
benchmark problems for which contact does not arise. Hence, through the use of
a small time step and an adaptive algorithm to refine the mesh when the solids
are near each other or near the boundaries, we have not needed to implement a
contact test or a contact model. In the future, we do intend to consider adding
a contact model in order to further generalize our method.
9.4. Conditioning of the linear system
f
If ρs ≥ ρf , and negleting τij
, the discretized linear equation system is guaranteed to be well-conditioned. However, this restriction is too stringent to be a
necessary condition. For example, we have implemented and tested a number of
cases for which ρs < ρf and the solid rises in a stable manner due to buoyancy.
24

400

405

9.5. Approximation for pressure
It is well known that the pressure jumps across the interface between the fluid
and solid, and that a high resolution is therefore needed near the interface in
order to capture this jump. In this article, we use an adaptive mesh refinement
near the interface to reduce the error caused by our continuous approximation
(P2 P1 element) for this discontinuous pressure. An alternative or additional
choice is to use P2 (P1 + C) elements (the shape function of pressue is enriched
by a constant) in order to capture an element-based jump of pressure. We intend
to test this element in the future.
10. Conclusion

410

415

420

425

430

In this article we introduce a one-field FD method for fluid-structure interaction, which can be applied to a wide range of problems, from small deformation
to very large deformation and from very soft solids through to very rigid solids.
Several numerical examples, which are widely used in the literature of IFEM and
FD methods with DLM (DLM/FD), are implemented to validate the proposed
method.
The one-field FD method combines features from both the IFEM and DLM/FD.
Nevertheless, it differs from each of them in the following aspects. Firstly, our
one-field FD method solves the solid and fluid equations together while the classical IFEM does not solve the solid equations. Although the implicit form of
IFEM can iteratively solve the solid equations, this is different from our onefield FD method which couples the fluid and solid equations monolithically via
a direct matrix addition as shown in formulas (35) and (36). Secondly, while
both our one-field FD method and DLM/FD solve solid equations, the former
solves for just one velocity field in the solid domain using FEM interpolation,
while the latter solves one velocity field and one displacement field in the solid
domain using Lagrange multipliers. In summary therefore we believe that the
one-field FD method has the potential to offer the robustness and range of operation of DLM/FD, but at a computational cost that is much closer to that of the
IFEM approaches. Expressed another way, we contend that our approach has
all of the advantages of IFEM techniques but the additional robustness usually
associated with more complex monolithic solvers.
Appendix A. Expressions of M, Ms , K, Ks , B, f and f s

435

In this appendix, the specific expressions for the mass matrices M and Ms ,
stiffness matrices K and Ks , matrix B and the force vectors f and f s in equations
(34), (35) and (36) are presented.
(1) M: (k, m = 1, 2, · · · N u )

M11
M = ρf


M22

, (M11 )km = (M22 )km = (ϕk , ϕm )Ωh .

25

(2) Ms : (k, m = 1, 2, · · · N s )

 Ms11
s
s
f
M = ρ −ρ


Ms22

, (Ms11 )km = (Ms22 )km = (ϕsk , ϕsm )Ωsh .

(3) K: (k, m = 1, 2, · · · N u )
K = µf
where


(K11 )km = 2




K12
,
K22

K11
K21

∂ϕk ∂ϕm
,
∂x1 ∂x1




+
Ωh

∂ϕk ∂ϕm
,
∂x2 ∂x2


,
Ωh






∂ϕk ∂ϕm
∂ϕk ∂ϕm
,
+
,
,
∂x2 ∂x2 Ωh
∂x1 ∂x1 Ωh




∂ϕk ∂ϕm
∂ϕk ∂ϕm
=
, (K21 )km = (K12 )mk =
.
,
,
∂x1 ∂x2 Ωh
∂x2 ∂x1 Ωh
(K22 )km = 2

(K12 )km

(4) Ks : (b, m = 1, 2, · · · N s )
Ks =

 s
K11
Ks21


Ks12
,
Ks22

(A.1)

where
 s


∂ϕb ∂ϕsm
∂ϕsb ∂ϕsm
s
+ µ ∆t
,
,
= µ ∆t2
∂x1 ∂x1 Ωsh
∂x2 ∂x2 Ωsh
 s n



∂ϕb ∂u1 ∂ϕsm
∂ϕsb ∂un2 ∂ϕsm
s
2
s
2
+ 2µ ∆t
,
+ µ ∆t
,
∂xk ∂xk ∂x1 Ωsh
∂xk ∂xk ∂x2 Ωsh
 s


 s
∂ϕb s n ∂un2 ∂ϕsm
∂ϕb s n ∂un1 ∂ϕsm
+ ∆t2
(τkl )
,
(τkl )
,
+ 2∆t2
∂xk
∂xl ∂x1 Ωsh
∂xk
∂xl ∂x2 Ωsh
 s
 s


∂ϕb s n ∂ϕsm
∂ϕb s n ∂ϕsm
+ 2∆t
(τ ) ,
+ ∆t
(τ ) ,
.
∂xk k1
∂x1 Ωsh
∂xk k2
∂x2 Ωsh

(Ks11 )bm

s



Ks22 can be expressed by changing the subscript 1 to 2 and 2 to 1 in the
formula of Ks11 .
 s

 n

∂ϕb ∂ϕsm
∂u1 ∂ϕsb ∂ϕsm
(Ks12 )bm = µs ∆t
,
+ µs ∆t2
,
∂x1 ∂x2 Ωsh
∂xk ∂xk ∂x2 Ωsh
 n



s
s
s
s
∂u1 s n ∂ϕb ∂ϕm
2
s n ∂ϕb ∂ϕm
+ ∆t
(τ )
,
+ ∆t (τ1k )
,
,
∂xk kl ∂xl ∂x2 Ωsh
∂xk ∂x2 Ωsh
and (Ks21 )bm = (Ks12 )mb .
(5) B: (k = 1, 2, · · · N p and m = 1, 2, · · · N u )
 


∂ϕm
B1
B=
, (Bi )mk = − φk ,
, (i = 1, 2).
B2
∂xi Ωh
26

(6) f : (m = 1, 2, · · · N u )
 

f
f = 1 , (fi )m = ρf (gi , ϕm )Ωh + h̄i , ϕm ΓN h , (i = 1, 2)
f2
(7) f s : (m = 1, 2, · · · N s )
 s

f
f s = 1s , (fis )m = ρs − ρf (gi , ϕsm )Ωsh
f2


s

∂un ∂unj
∂un s n ∂unj
s n ∂ϕm
+ µs ∆t2 i
+ ∆t2 i (τkl
)
− τij
,
, (i = 1, 2).
∂xk ∂xk
∂xk
∂xl
∂xj Ωsh
References

440

[1] G. Hou, J. Wang, A. Layton, Numerical methods for fluid-structure interactiona review, Commun. Comput. Phys 12 (2) (2012) 337–377. doi:
10.4208/cicp.291210.290411s.
[2] U. Küttler, W. A. Wall, Fixed-point fluid–structure interaction solvers with
dynamic relaxation, Computational Mechanics 43 (1) (2008) 61–72. doi:
10.1007/s00466-008-0255-5.

445

450

[3] J. Degroote, K.-J. Bathe, J. Vierendeels, Performance of a new partitioned procedure versus a monolithic procedure in fluid–structure interaction, Computers & Structures 87 (11-12) (2009) 793–801. doi:10.1016/
j.compstruc.2008.11.013.
[4] M. Heil, An efficient solver for the fully coupled solution of largedisplacement fluid–structure interaction problems, Computer Methods in
Applied Mechanics and Engineering 193 (1-2) (2004) 1–23. doi:10.1016/
j.cma.2003.09.006.
[5] M. Heil, A. L. Hazel, J. Boyle, Solvers for large-displacement fluid–structure
interaction problems: segregated versus monolithic approaches, Computational Mechanics 43 (1) (2008) 91–101. doi:10.1007/s00466-008-0270-6.

455

460

[6] R. L. Muddle, M. Mihajlović, M. Heil, An efficient preconditioner for
monolithically-coupled large-displacement fluid–structure interaction problems with pseudo-solid mesh updates, Journal of Computational Physics
231 (21) (2012) 7315–7334. doi:10.1016/j.jcp.2012.07.001.
[7] R. C. Peterson, P. K. Jimack, M. A. Kelmanson, The solution of
two-dimensional free-surface problems using automatic mesh generation, International journal for numerical methods in fluids 31 (6)
(1999) 937–960. doi:10.1002/(SICI)1097-0363(19991130)31:6<937::
AID-FLD906>3.0.CO;2-p.

27

465

470

[8] M. A. Walkley, P. H. Gaskell, P. K. Jimack, M. A. Kelmanson, J. L.
Summers, Finite element simulation of three-dimensional free-surface
flow problems, J Sci Comput 24 (2) (2005) 147–162. doi:10.1007/
s10915-004-4611-0.
[9] L. Zhang, A. Gerstenberger, X. Wang, W. K. Liu, Immersed finite element method, Computer Methods in Applied Mechanics and Engineering
193 (21) (2004) 2051–2067. doi:doi:10.1016/j.cma.2003.12.044.
[10] L. Zhang, M. Gay, Immersed finite element method for fluid-structure
interactions, Journal of Fluids and Structures 23 (6) (2007) 839–857.
doi:10.1016/j.jfluidstructs.2007.01.001.

475

[11] X. Wang, L. T. Zhang, Interpolation functions in the immersed boundary
and finite element methods, Computational Mechanics 45 (4) (2009) 321–
334. doi:10.1007/s00466-009-0449-5.
[12] X. Wang, C. Wang, L. T. Zhang, Semi-implicit formulation of the immersed
finite element method, Computational Mechanics 49 (4) (2011) 421–430.
doi:10.1007/s00466-011-0652-z.

480

485

490

[13] X. Wang, L. T. Zhang, Modified immersed finite element method for fullycoupled fluid–structure interactions, Computer Methods in Applied Mechanics and Engineering 267 (2013) 150–169. doi:10.1016/j.cma.2013.
07.019.
[14] R. Glowinski, T. Pan, T. Hesla, D. Joseph, J. Périaux, A fictitious domain
approach to the direct numerical simulation of incompressible viscous flow
past moving rigid bodies: Application to particulate flow, Journal of Computational Physics 169 (2) (2001) 363–426. doi:10.1006/jcph.2000.6542.
[15] Z. Yu, A DLM/FD method for fluid/flexible-body interactions, Journal of
Computational Physics 207 (1) (2005) 1–27. doi:10.1016/j.jcp.2004.
12.026.
[16] F. P. Baaijens, A fictitious domain/mortar element method for fluidstructure interaction, International Journal for Numerical Methods in
Fluids 35 (7) (2001) 743–761. doi:10.1002/1097-0363(20010415)35:
7<743::AID-FLD109>3.0.CO;2-A.

495

500

[17] C. Kadapa, W. Dettmer, D. Perić, A fictitious domain/distributed lagrange
multiplier based fluid–structure interaction scheme with hierarchical bspline grids, Computer Methods in Applied Mechanics and Engineering
301 (2016) 1–27. doi:10.1016/j.cma.2015.12.023.
[18] C. Hesch, A. Gil, A. A. Carreño, J. Bonet, P. Betsch, A mortar approach for
fluid–structure interaction problems: Immersed strategies for deformable
and rigid bodies, Computer Methods in Applied Mechanics and Engineering
278 (2014) 853–882. doi:10.1016/j.cma.2014.06.004.
28

[19] C. S. Peskin, The immersed boundary method, Acta numerica 11 (2002)
479–517. doi:10.1016/j.cma.2015.12.023.
505

510

[20] A. Robinson-Mosher, C. Schroeder, R. Fedkiw, A symmetric positive definite formulation for monolithic fluid structure interaction, Journal of Computational Physics 230 (4) (2011) 1547–1566. doi:10.1016/j.jcp.2010.
11.021.
[21] E. Hachem, S. Feghali, R. Codina, T. Coupez, Anisotropic adaptive meshing and monolithic variational multiscale method for fluid–structure interaction, Computers & Structures 122 (2013) 88–100. doi:10.1016/j.
compstruc.2012.12.004.
[22] O. Pironneau, An energy preserving monolithic Eulerian fluid-structure
numerical scheme, arXiv:1607.08083v1 [] 27 Jul 2016.

515

520

[23] F. Auricchio, D. Boffi, L. Gastaldi, A. Lefieux, A. Reali, A study on unfitted
1d finite element methods, Computers & Mathematics with Applications
68 (12) (2014) 2080–2102. doi:10.1016/j.camwa.2014.08.018.
[24] A. Gerstenberger, W. A. Wall, An extended finite element
method/Lagrange multiplier based approach for fluid-structure interaction, Computer Methods in Applied Mechanics and Engineering
197 (19-20) (2008) 1699–1714. doi:10.1016/j.cma.2007.07.002.
[25] S. Frei, Eulerian finite element methods for interface problems and fluidstructure interactions, Ph.D. thesis, Universitt Heidelberg (2016).

525

[26] T. Richter, A fully Eulerian formulation for fluid–structure-interaction
problems, Journal of Computational Physics 233 (2013) 227–240. doi:
10.1016/j.jcp.2012.08.047.
[27] T. Wick, Coupling of fully Eulerian and arbitrary Lagrangian–Eulerian
methods for fluid-structure interaction computations, Computational Mechanics 52 (5) (2013) 1113–1124. doi:10.1007/s00466-013-0866-3.

530

535

[28] T. Dunne, R. Rannacher, Adaptive finite element approximation of fluidstructure interaction based on an Eulerian variational formulation, in: Lecture Notes in Computational Science and Engineering, Springer Science
Business Media, 2006, pp. 110–145. doi:10.1007/3-540-34596-5_6.
[29] O. Zienkiewic, The finite element method for fluid dynamics, 6th Edition,
Elsevier BV, 2005.
[30] P. B. Bochev, M. D. Gunzburger, Least-squares finite element methods,
Vol. 166, Springer Science & Business Media, 2009.

540

[31] A. K. Gupta, A finite element for transition from a fine to a coarse grid,
International Journal for Numerical Methods in Engineering 12 (1) (1978)
35–45. doi:10.1002/nme.1620120104.
29

[32] T.-P. Fries, A. Byfut, A. Alizada, K. W. Cheng, A. Schrder, Hanging nodes
and XFEM, International Journal for Numerical Methods in Engineering
86 (4-5) (2010) 404–430. doi:10.1002/nme.3024.

545

550

[33] W. Bangerth, O. Kayser-Herold, Data structures and requirements for hp finite element software, ACM Transactions on Mathematical Software 36 (1)
(2009) 1–31. doi:10.1145/1486525.1486529.
[34] N. Zander, T. Bog, S. Kollmannsberger, D. Schillinger, E. Rank, Multilevel hp-adaptivity: high-order mesh adaptivity without the difficulties of
constraining hanging nodes, Computational Mechanics 55 (3) (2015) 499–
517. doi:10.1007/s00466-014-1118-x.
[35] H. Zhao, J. B. Freund, R. D. Moser, A fixed-mesh method for incompressible flow–structure systems with finite solid deformations, Journal of
Computational Physics 227 (6) (2008) 3114–3140. doi:10.1016/j.jcp.
2007.11.019.

555

[36] W. A. Wall, Fluid-struktur-interaktion mit stabilisierten finiten elementen,
Ph.D. thesis, Universitt Stuttgart (1999). doi:10.18419/OPUS-127.
[37] Y. Bazilevs, K. Takizawa, T. E. Tezduyar, Computational Fluid-Structure
Interaction: Methods and Applications, Wiley-Blackwell, 2013. doi:10.
1002/9781118483565.

30

