A monolithic fluid-structure interaction formulation for solid
and liquid membranes including free-surface contact
Roger A. Sauer1 and Tobias Luginsland2

arXiv:1710.02128v1 [physics.comp-ph] 5 Oct 2017

Aachen Institute for Advanced Study in Computational Engineering Science (AICES),
RWTH Aachen University, Templergraben 55, 52056 Aachen, Germany

Abstract
A unified fluid-structure interaction (FSI) formulation is presented for solid, liquid and mixed
membranes. Nonlinear finite elements (FE) and the generalized-α scheme are used for the
spatial and temporal discretization. The membrane discretization is based on curvilinear surface
elements that can describe large deformations and rotations, and also provide a straightforward
description for contact. The fluid is described by the incompressible Navier-Stokes equations,
and its discretization is based on stabilized Petrov-Galerkin FE. The coupling between fluid
and structure uses a conforming sharp interface discretization, and the resulting non-linear
FE equations are solved monolithically within the Newton-Raphson scheme. An arbitrary
Lagrangian-Eulerian formulation is used for the fluid in order to account for the mesh motion
around the structure. The formulation is very general and admits diverse applications that
include contact at free surfaces. This is demonstrated by two analytical and three numerical
examples. They include balloon inflation, droplet rolling and flapping flags. They span a
Reynolds-number range from 0.001 to 2000. One of the examples considers the extension to
rotation-free shells using isogeometric FE.
Keywords: arbitrary Lagrangian-Eulerian formulation, contact mechanics, incompressible
Navier-Stokes equations, isogeometric finite elements, nonlinear membranes, surface tension

1

Introduction

Fluid-structure interaction (FSI) problems are challenging problems due to various reasons.
They combine the computational challenges of (generally non-linear) fluid and structural mechanics, and they introduce new challenges, both physical and numerical, due to the coupling.
If the structure is highly flexible, such as a thin membrane, large deformations can be expected.
Those, in turn, have a large influence on the fluid flow. A comprehensive overview of FSI and
its challenges is given by the monographs of Ohayon (2004), Bazilevs et al. (2013) and Bazilevs
and Takizawa (2016). The classical focus in FSI problems is on solid structures. However, some
structures are not solids but rather fluids or fluid-like objects. Examples are liquid menisci,
soap films and lipid bilayers. Lipid bilayers surround biological cells. They are characterized by
both solid-like (i.e. elastic bending) and fluid-like behavior (i.e. in-plane flow). Further, liquid
(and solid) membranes can come into contact with surrounding objects. A classical example
is a liquid droplet rolling on a substrate. The problem is characterized by fluid flow, surface
tension and contact.
While there are various formulations available in the present literature that capture all these
1
2

corresponding author, email: sauer@aices.rwth-aachen.de
current affiliation: Daimler AG, 71059 Sindelfingen, Germany

1

aspects, there is no formulation that unifies them all into a single framework. This is the objective of the present work. In doing so, we build on our recent computational work on contact,
membranes, shells and fluid dynamics.
The presented formulation is based on finite elements (FE) using an interface tracking technique
based on a sharp interface formulation. There is a large literature body on FE-based work on
membrane-FSI that is surveyed in the following. The computational approaches on interactions
between fluids and membrane-like structures can be sorted into two groups. The first group
deals with solid structures like elastic membranes and flexible shells, while the second group
is concerned with liquid membranes and menisci. The first group can be further sorted into
approaches that use surface formulations (based on shell and membrane theories) and contributions that use bulk formulations. The second group can be further sorted into approaches
that only account for the shape equation in order to characterize the liquid membrane (like the
Young-Laplace equation), and approaches that also account for in-plane equations (such as the
surface Navier-Stokes equations). The latter case is necessary for liquid membranes that are
not surrounded by a fluid, and consequently the FSI problem is due to the interplay of membrane shape and surface flow. If a surrounding medium is considered, and no-slip conditions
are applied on the membrane surface, the flow within the membrane is already captured by the
bulk flow, and so no further equations are needed. The method presented here is based on a
surface formulation that accounts for both shape and in-plane equations.
The following references deal with solid membranes using surface formulations. In Liang et al.
(1997) the authors employ a deformable spatial domain space-time FEM to study the interaction
of an incompressible fluid with an elastic membrane. Bletzinger et al. (2006) compute the flow
around a tent structure using a staggered coupling between a shell code and a CFD code.
Tezduyar and Sathe (2007) review their FSI formulation based on space-time FE and introduce
advancements regarding accuracy, robustness and efficiency. Benchmark examples include the
inflation of a balloon, the flow through a flexible diaphragm in a tube as well as a descending
parachute. Parachutes are also analyzed in Karagiozis et al. (2011) and Takizawa and Tezduyar
(2012) using thin-shell formulations. Le et al. (2009) developed an implicit immersed boundary
method for the incompressible Navier-Stokes equations to simulate membrane-fluid interactions.
Their examples include an oscillating spherical ball immersed in a fluid and the stretching of
a red blood cell in a pressure driven shear flow. van Opstal et al. (2015) present a hybrid
isogeometric finite-element/boundary element method for fluid-structure interaction problems
of inflatable structures such as airbags and balloons. Boundary elements are also used in a
recent isogeometric FSI formulation for Stokes flow around thin shells (Heltai et al., 2016).
The following references deal with solid membranes using bulk formulations. Kloeppel and Wall
(2011) numerically investigate the flow inside red blood cells (RBC) by means of monolithically
coupling an incompressible fluid to a lipid bilayer represented by incompressible solid shell
elements. In Franci et al. (2016) the authors develop a monolithic strategy for the description
of purely Lagrangian FSI problems. For the solid, the FEM is used, while the fluid is discretized
using the so-called Particle FEM (Idelsohn et al., 2004). Yang et al. (2016) introduce a finitediscrete element method for bulk solids and combine the developed numerical model with a finite
element multiphase flow model. Only 2D examples are considered, such as a rigid structure
floating on a liquid-gaseous interface.
Recent reviews on computational FSI methods for solids have been given by Dowell and Hall
(2001), van Loon et al. (2007) and Bazilevs et al. (2013). For an introduction to immersedboundary methods as an alternative to conforming FE discretizations we refer to Peskin (2003).
The following references deal with liquid membranes governed only by a shape equation. Walkley
et al. (2005) present an ALE framework for the solution of free surface flow problems including
a dynamic contact line model and show its capabilities for the case of a sliding droplet. Saksono
2

and Perić (2006) propose a 2D finite element formulation for surface tension and apply it to
oscillating droplets and stretched liquid bridges. Montefuscolo et al. (2014) introduce high-order
ALE FEM schemes for capillary flows. The schemes are demonstrated on oscillating and sliding
droplets accounting for the varying contact angles.
The following references deal with liquid membranes governed by shape and in-plane equations. Barett et al. (2015) present a numerical study of the dynamics of lipid bilayer vesicles.
A parametric finite element formulation is introduced to discretize the surface Navier-Stokes
equations. Rangarajan and Gao (2015) introduce a spline-based finite-element formulation to
compute equilibrium configurations of liquid membranes. Sauer et al. (2017) present a new
isogeometric finite element formulation for liquid membranes that accounts for the in-plane viscosity and incompressibility of the liquid.
A general introduction to fluid membranes and vesicles and their configurations observed in nature is given by Seifert (1997). For a review on the droplet dynamics within flows, see Cristini
and Tan (2004).
There is also earlier work on combining a contact and FSI. It can be grouped into two categories:
Either contact is considered between solids submerged within the fluid (e.g. see Tezduyar et al.
(2006); Mayer et al. (2010)), or contact is considered at free liquid surfaces. For liquid surfaces
the same classical contact algorithms as for solid surfaces can be used (Sauer, 2014). An
alternative treatment of free surface contact appears naturally in the Particle FEM (Idelsohn
et al., 2006). Additionally, the contact behavior between liquids and solids is also governed by
a contact angle and its hysteresis during sliding contact. A general computational algorithm
for contact angle hysteresis is given in Sauer (2016).
Existing work is motivated by specific examples that either focus on solid or liquid membranes.
The aim of this paper therefore is to provide a unified FSI formulation that is suitable to describe
solid membranes, such as sheets, fabrics and tissues, liquid membranes, such as menisci and
soap films, and membranes with both solid- and liquid-like character, like lipid bilayers. The
present work considers a monolithic coupling scheme between fluid and structure, and solves
the resulting non-linear system of equations with the Newton-Raphson method. Finite elements
and the generalized-α scheme are used for spatial and temporal discretization. The formulation
uses a conforming interface discretization and an ALE formulation for the mesh motion.
The following aspects are new:
• A unified monolithic FSI formulation for liquid and solid membranes is presented.
• It includes contact on free liquid surfaces, and
• it easily extends to rotation-free shells.
• Two simple analytical FSI examples are presented.
• The formulation is suitable for a wide range of applications, including free-surface flows,
liquid menisci, flags and flexible wings.
• The examples include a flow and contact analysis of a rolling 3D droplet.
The remainder of this paper is structured as follows. Sec. 2 presents the governing theory of
incompressible fluid flow, nonlinear membranes and their coupling. The theory is used to solve
two simple analytical FSI examples in Sec. 3. The computational treatment is then presented
in Sec. 4 using finite elements for the spatial discretization of fluid and membrane, and the
generalized-α scheme for the temporal discretization of the coupled system. Sec. 5 presents
four numerical examples ranging from very low to quite large Reynolds numbers. The paper
concludes with Sec. 6.
3

2

Governing equations

This section summarizes the governing equations for fluid flow, membrane deformation, membrane contact and their coupling. The symbols F and S are used to denote the fluid domain
and the membrane surface, cf. Fig. 1 in Sec. 3.1 and Fig. 12 in Sec. 5.3.

2.1

Fluid flow

The fluid motion is described by an arbitrary Lagrangian-Eulerian (ALE) formulation. It is
therefore necessary to distinguish between the material motion and the mesh motion. An ALE
formulation contains the special cases of a purely Lagrangian description, for which the material
and mesh motion coincide, and a purely Eulerian description, for which the mesh motion is zero.
2.1.1

Fluid kinematics

The material motion of a fluid particle X within domain F is characterized by the deformation
mapping
x = ϕ(X, t)
(1)
and the corresponding deformation gradient (or Jacobian)
∂ϕ
(2)
.
∂X
The volume change during deformation is captured by the Jacobian determinant J := det F .
The velocity of the material is given by the time derivative of x for fixed X, written as
F :=

∂x
(3)
∂t X
and commonly referred to as the material time derivative. It is also often denoted by the dot
notation v = ẋ. An important object characterizing the fluid flow is the velocity gradient
v :=

∂v
∂x

L := ∇v =

(4)

that can also be written as L = Ḟ F −1 , where Ḟ is the material time derivative of the deforma
tion gradient. The symmetric part of the velocity gradient is denoted by D := L + LT /2.
Likewise to Eq. (3), the material acceleration is given by
a := v̇ =

∂v
∂t

X

.

(5)

It is related to the acceleration for fixed x,
v 0 :=

∂v
∂t

x

,

(6)

according to
v̇ = v 0 + L (v − v m ) ,

(7)

where v m is the mesh velocity (Donea and Huerta, 2003). For a purely Lagangian description
v m = v, while for a purely Eulerian description v m = 0.
Remark 2.1: The gradient operator appearing in Eq. (4) (and likewise in Eq. (2)), is defined
here as ∇v := vi,j ei ⊗ ej .3
3
Following index notation, summation is implied on repeated indices. Latin indices range from 1 to 3 and
refer to Cartesian coordinates. Greek indices range from 1 to 2 and refer to curvilinear surface coordinates.

4

2.1.2

Fluid equilibrium

From the balance of linear momentum within F follows the equilibrium equation
div σ + f̄ = ρ v̇

in F ,

(8)

which governs the fluid flow together with the boundary conditions
v = v̄
σn = t = t̄

on ∂x F ,
on ∂t F .

(9)

Here, σ denotes the stress tensor within F, t denotes the traction vector on the surface characterized by normal vector n, and ρ denotes the density, while f̄ , v̄ and t̄ are prescribed body
forces, surface velocities and surface tractions. ∂x F and ∂t F denote the corresponding Dirichlet
and Neumann boundary regions of the fluid domain F. Boundary ∂x F can be split into the
two parts
∂x F = S ∪ ∂x̂ F ,
(10)
where S is the surface of the membrane, which is considered to impose its velocity onto the
fluid, and ∂x̂ F denotes the remaining Dirichlet boundary of the fluid domain. In order to solve
PDE (8) for v(x, t), the initial condition
v(x, 0) = v 0 (x)

(11)

is needed.
2.1.3

Fluid constitution

We consider an incompressible Newtonian fluid with kinematic viscosity ν and dynamic viscosity
η = νρ. In that case the stress tensor is given by
σ = −p 1 + 2η D ,

(12)

where p is the Lagrange multiplier to the incompressibility constraint
g := J − 1 = 0 ,

(13)

div v = 0 .

(14)

which is equivalent to the condition
A consequence of this condition is that the fluid pressure, defined as −tr σ/3, is equal to the
Lagrange multiplier p. It is an additional unknown that needs to be solved for together with v.
In order to uniquely determine the pressure field within F, p needs to be specified at one point
in the fluid domain.
2.1.4

Fluid weak form

In order to solve the problem with finite elements the strong form equations (8), (9.2) and (14)
are reformulated in weak form. They are therefore multiplied by the test functions w and q,
and integrated over the domain F. Function w is assumed to be zero on the Dirichlet boundary
∂x̂ F, but non-zero on the surface S. Functions w and q are further assumed to possess sufficient

5

regularity for the following integrals to be well defined. In the framework of SUPG4 and PSPG5
stabilization, the weak form takes the form
GF := GF in + GF int + Gsupg − GF s − GF ext = 0
GG := Gg + Gpspg = 0
where

∀ w∈W,
∀ q ∈ Q,

(15)

Z
GF in =

w · ρ v̇ dv

(16)

∇w : σ dv

(17)

w · t da

(18)

F

is the virtual work associated with inertia,
Z
GF int =
F

is internal virtual work,
Z
GF s =
S

is the virtual work of the fluid traction t = σn on boundary S,
Z
Z
GF ext =
w · f̄ dv +
w · t̄ da
F

(19)

∂t F

is the external virtual work6 ,
Z
Gg =

q div v dv

(20)

F

is the virtual work associated with incompressibility constraint (14),
Z
Gsupg =
τv fres · ∇w (v − v m ) dv

(21)

F

is the SUPG term,
Z
Gpspg =
F

τp ∇q · fres dv

(22)

is the PSPG term, and
fres := ρ v̇ − div σ − f̄

(23)

is the residual of Eq. (8). Dimensionally the residual is a force per volume. Since in theory
fres = 0, stabilization terms Gsupg and Gpspg do not affect the physical behavior of the system.
In Cartesian coordinates fres · ∇w (v − v m ) = fires wi,j (vj − vmj ). The scalars τv and τp are
stabilization parameters that are discussed in Sec. 4.

2.2

Deforming membranes

This work focuses on pure membranes that do not resist bending and out-of-plane shear. The
description of those membranes is based on the formulation of Sauer et al. (2014), which admits
both solid and liquid membranes. What follows is a brief summary.
4

Streamline upwind/Petrov-Galerkin (Brooks and Hughes, 1982)
Pressure stabilizing/Petrov-Galerkin (Hughes et al., 1986)
6
In the following examples we consider zero Neumann BC (t̄ = 0) and constant gravity loading with f̄ = ρ g.
5

6

2.2.1

Membrane kinematics

The motion of a membrane surface S is fully described by the mapping
x = x(ξ α , t) ,

(24)

where ξ α , for α = 1, 2, are curvilinear coordinates that can be associated with material points
on the surface. They can be conveniently taken from the parameterization of the finite element
shape functions. Based on mapping (24), the tangent vectors aα p
= ∂x/ξ α to surface S, the
metric tensor components aαβ = aα · aβ , the surface stretch Js = det[aαβ ], and the surface
normal n = a1 × a2 /Js can be determined. From the matrix inverse [aαβ ] = [aαβ ]−1 , the dual
tangent vectors aα = aαβ aβ can be defined such that aα · aβ is equal to the Kronecker delta δβα .
In order to characterize deformation, a stress-free reference configuration S0 is introduced. It
will be considered here as the initial membrane surface, i.e. S0 := S|t=0 . In the reference
configuration the tangent vectors, metric tensor components, inverse components and normal
vector are denoted by capital letters, i.e. Aα , Aαβ , Aαβ and N . The in-plane deformation of
surface S is fully characterized by the difference between Aαβ and aαβ .
Following definitions (3) and (5), the membrane velocity v and acceleration a are obtained from
Eq. (24).
2.2.2

Membrane equilibrium

From the balance of linear momentum within S follows the equilibrium equation
(σs aα );α + fs = ρs v̇

in S ,

(25)

which governs the membrane deformation together with the boundary conditions
for x ∈ ∂x S ,

x = x̄

σs ν = ts = t̄s for x ∈ ∂t S ,

(26)

e.g. see Sauer and Duong (2017). Here, σ s denotes the stress tensor within S, ts denotes the
traction vector on the membrane boundary characterized by normal vector ν, and ρs denotes the
membrane density, while x̄ and t̄s are prescribed boundary velocities and boundary tractions.
The body force fs is considered here to have contributions coming from the flow field, contact
and external sources, i.e.
fs = ff + fc + f̄s .
(27)
In order to solve PDE (25) for x(ξ α , t), the initial conditions
x(ξ α , 0) = X(ξ α ) ,
v(ξ α , 0) = v 0 (ξ α ) ,

(28)

are needed.
2.2.3

Membrane constitution

For pure membranes, the stress tensor only has in-plane components, i.e. it has the format
σs = σ αβ aα ⊗ aβ . Two material models are considered in this work. The first,


µ
1 αβ
αβ
αβ
σ =
A − 2a
,
(29)
J
J
7

is suitable for solid membranes. It can be derived from the 3D incompressible Neo-Hookean
material model (Sauer et al., 2014). The second,
σ αβ = γ aαβ ,

(30)

models isotropic surface tension, and is suitable to describe liquid membranes, e.g. see Sauer
(2014). The parameters µ and γ denote the shear stiffness and the surface tension, respectively.
Both are considered constant here.
2.2.4

Membrane contact

This work also considers that sticking contact can occur on the membrane surface Sc ⊂ S.
During sticking contact no relative motion occurs between the membrane and a neighboring
substrate surface Ssub . Mathematically this corresponds to the constraint
g=0

∀ x ∈ Sc ,

(31)

where
g(x) = x − x0p

(32)

denotes the contact gap between the membrane point x ∈ Sc and its initial projection point on
the substrate surface, x0p ∈ Ssub , i.e. x0p is the location where x initially touched Ssub . Here,
constraint (31) will be enforced by a penalty regularization. For this, the contact traction at
x ∈ S is given by
(
− g if g · nc < 0 ,
fc =
(33)
0
else ,
where nc is the surface normal of Ssub . Further details on large deformation contact theory can
be found in the textbooks of Laursen (2002) and Wriggers (2006).
2.2.5

Membrane weak form

In order to employ finite elements, the strong form equations (25) and (26.2) are reformulated
in weak from. As shown in Sauer and Duong (2017), the weak form for the membrane can be
written as
GS := GSin + GSint + Gc − GSf − GSext = 0 ∀ w ∈ W ,
(34)
with the virtual work contributions
Z
w · ρs v̇ da ,

GSin :=
S

Z

σ αβ w;α · aβ da ,
Z
Gc := − w · fc da ,
S
Z
GSf :=
w · ff da ,
S
Z
Z
GSext :=
w · f̄s da +
w · t̄s ds ,
GSint :=

S

S

(35)

∂t S

due to inertia, internal forces, contact forces, fluid forces and external forces acting on S and
∂t S. Test function w is the same as in (15). Therefore, space W needs to additionally satisfy
8

the requirement that all integrals appearing above are well defined. Further w is assumed to
be zero on ∂x S.
Pure membranes are inherently unstable in the quasi-static case (v = v̇ = 0) and therefore
need to be stabilized (Sauer et al., 2014; Sauer, 2014). Here, no stabilization is required as the
fluid forces f f stabilize the membrane, even when ρs = 0. In the numerical examples following
later, f̄s and t̄s , and consequently GSext , are considered zero.
Remark 2.2: It is straight forward to extend weak form (34) to Kirchhoff-Love shells: GSint and
GSext simply need to be extended by the bending moments acting within S and on ∂S (Duong
et al., 2017). Kirchhoff-Love shells are suitable for thin membrane-like surface structures. Such
a structure is considered in Sec. 5.3 using isogeometric finite elements.

2.3

Coupling conditions

The membrane deformation x moves the fluid such that
v = ẋ on S

(36)

is a Dirichlet BC for the fluid. This choice assumes no tangential slip between membrane and
fluid. In response, the flow exerts a traction on the membrane such that
ff = −t on S

(37)

is a ‘body force’ of the membrane. Eq. (36) is the kinematic coupling condition between the
two domains, while Eq. (37) is the kinetic coupling condition. The combined FSI problem is
then characterized by the two governing equations
GF + GS = 0

∀ w ∈ W ∩ Ws ,

GG

∀q∈Q,

= 0

(38)

which can be solved for the unknown velocity v and pressure p in F. The membrane deformation
can then be obtained from integrating v. Coupling condition (37) simply leads to the cancelation
of terms GF s and GSf in the combined weak form (38). This cancelation will carry over to the
discretized weak form, as long as surface S is discretized conformingly on the fluid and membrane
side.

3

Analytical examples

This section presents the analytical solution of two simple examples. They serve as verification
examples for the computational implementation.

3.1

Solid membrane example: Fluid-inflated cylinder

As a first example we consider the radial inflation of a membrane cylinder due to a radial inflow
as is illustrated in Fig. 1. The example is chosen since it can be fully solved analytically and
thus used for verification of the computational formulation, which is then considered in Sec. 5.1.
Given the inflow velocity vin at the inner boundary Rin , the radial fluid velocity at location r
is given by
vin Rin
(39)
v(r) =
r
9

Figure 1: Fluid-inflated cylinder: Membrane deformation S0 → S and fluid velocity v(r) due
to a radial inflow at Rin .

due to continuity. Since v = ṙ, we obtain
r(R, t) =

p
R2 + 2vin (t) Rin t ,

(40)

as the current position of the fluid particle initially at R. The current membrane position is thus
given by rs = r(Rs , t), where Rs is the initial position of the membrane. In vectorial notation,
the flow field can thus be characterized by the position, velocity and acceleration
x(R, t) = r er ,
v(R, t) = v er ,
v2
a(R, t) = − er ,
r

(41)

where er = cos θ e1 + sin θ e2 is the radial unit vector. From this follows
D = − sin2 θ (e1 ⊗ e1 ) − cos2 θ (e2 ⊗ e2 ) − 2 sin θ cos θ (e1 ⊗ e2 + e2 ⊗ e1 ) ,

(42)

such that div D = 0. The equation of motion thus reduces to −∇p = ρ a, which can be
integrated to give the pressure field

ρ
p(R, t) = ps + vs2 − v 2 ,
(43)
2
where vs = v(rs ) is the current membrane velocity, and ps is the pressure acting on the membrane. This pressure equilibrates the membrane stress


1
σ =µ λ− 3
(44)
λ
caused by the membrane stretch λ = rs /Rs according to Eq. (29); see Appendix A. From
ps = σ/rs follows
"
 4 #
µ
Rs
ps =
1−
.
(45)
Rs
rs

3.2

Liquid membrane example: Spinning droplet

As a second example we consider a spinning droplet. This example is considered for comparison
with the computational example of a rolling droplet in Sec. 5.2. At very small length scales
10

the influence of gravity is negligible, so that a rolling droplet remains approximately spherical.
Considering the axis of rotation to be e2 , the motion of a spinning droplet can be expressed as
x(r, t) = r er ,

(46)

where er = cos θ e1 − sin θ e3 , θ = ωt and ω denotes the angular velocity around e2 . Consequently,
v(r, t) = ω r eθ ,
(47)
a(r, t) = −ω 2 r er ,
where eθ = − sin θ e1 − cos θ e3 . Since we can write x1 = r cos θ and x2 = −r sin θ, we find
∇v = ω(e1 ⊗ e3 − e3 ⊗ e1 ) such that D = 0 and
σ = −p 1 .

(48)


The spin tensor, defined as W := L − LT /2, then becomes W = L = ∇v. The axial vector of
W , denoted by ω, thus is ω = ω e2 . It denotes the orientation and magnitude of the droplet’s
spin, and it is equal to half of the vorticity ∇ × v. Solving Eq. (8) (with f̄ = 0) for p now gives
p(r) = ω 2 ρ

r2
+ p0 .
2

(49)

The constant p0 follows from the boundary condition p(r0 ) = 2γ/r0 , where γ is the surface
tension of the droplet and r0 is the droplet radius. This condition enforces the Young-Lapace
equations, which is contained inside Eq. (25), see Sauer (2014). Applying the boundary condition, we find

2γ ρ ω 2 2
(50)
p(r) =
r0 − r2 .
−
r0
2
If desired, the constant velocity v 0 = ω r0 e1 can be added to v(r, t), such that the resulting
velocity is zero at the contact point (where θ = π/2).

4

Finite element formulation

The coupled fluid-membrane problem of Sec. 2 is solved with the finite element method using the
generalized-α scheme. This section presents the required discretization steps and the resulting
algebraic equations.

4.1

Spatial discretization

The computational domain is discretized into nel finite elements, numbered e = 1, ..., nel . Some
of these elements are 3D fluid elements, others are 2D surface elements or 1D line elements.
Element e contains ne nodes and occupies the domain Ωe in the current configuration. Each
fluid element has four degrees-of-freedom (dofs) per node (three velocity components and a
pressure), while the membrane elements each have three unknown displacements per node.
Each fluid element therefore contributes 4ne force components, while each membrane element
contributes 3ne force components that need to be assembled into the global system. Those
elemental forces are discussed in the following two sections.

11

4.1.1

Fluid flow

4.1.1.1 Basic flow variables
Within a fluid element, the fluid velocity is approximated by the interpolation
v≈

vh

=

ne
X

NI v I ,

(51)

I=1

where NI and v I are the nodal shape function and nodal velocity, respectively. In short, this
can also be written as
v ≈ v h = N ve ,
(52)
where N := [N1 1, N2 1, ..., Nne 1] and ve := [v 1 , v 2 , ..., v ne ]T . The corresponding test function
(or variation) is approximated in the same fashion, i.e.
w ≈ wh = N we .

(53)

The fluid pressure is approximated by the interpolation
p ≈ ph = Ñ pe ,

(54)

where Ñ := [N1 , N2 , ..., Nne ]. Likewise,
q ≈ q h = Ñ qe .

(55)

The structure of (52) is also used to interpolate the mesh motion, i.e.
v m ≈ v hm = N vme .

(56)

In the present work, the vme are not treated as unknowns. Instead they will be defined through
the membrane motion.
4.1.1.2 Derived flow variables
As a consequence of the above expressions, we find the approximation of the acceleration (from
Eq. (7))

(57)
v̇ ≈ v̇ h = N ve0 + LN ve − vme ,
the velocity gradient
h

L≈L =

ne
X

v I ⊗ ∇NI ,

(58)

I=1

the pressure gradient
∇p ≈ ∇ph = G pe ,

(59)

div v ≈ div v h = D ve ,

(60)



NI,1
∇NI = NI,2  ,
NI,3

(61)

and the velocity divergence
where

12

G := [∇N1 , ∇N2 , ..., ∇Nne ] and D := [(∇N1 )T , (∇N2 )T , ..., (∇Nne )T ]. Further, we introduce
the classical B-matrix B := [B1 , B2 , ..., Bne ], with


NI,1
0
0
 0
NI,2
0 


 0
0
NI,3 

,
BI := 
(62)

0
N
N
I,3
I,2


NI,3
0
NI,1 
NI,2 NI,1
0
in order to express the symmetric velocity gradient and its corresponding variation in Voigt
notation (indicated by index ‘v’) as
∇s v v ≈ ∇s v hv = B ve ,
∇s wv ≈ ∇s v hw = B we ,

(63)

i.e. arranged as ∇s v v := [v1,1 , v2,2 , v3,3 , v2,3 + v3,2 , v1,3 + v3,1 , v1,2 + v2,1 ]T . The stress tensor,
arranged as σv := [σ11 , σ22 , σ33 , σ23 , σ13 , σ12 ], can thus be written as
σ v ≈ σ hv = C B ve − 1v Ñ pe ,

(64)

with C := diag(2η1, η1) and 1v = [1, 1, 1, 0, 0, 0]T . Here, 1 is the usual identity tensor in R3 .
Due to the symmetry of the stress and since BT 1v = DT , the integrand of GF int becomes
∇wh : σ h = weT BT C B ve − weT DT Ñ pe

(65)

within element Ωe .
In order to represent the SUPG term, we introduce the arrays Bf := [Bf1 , Bf2 , ..., Bfne ], with
the 3 × 3 blocks
BfI := ∇NI ⊗ fres ,
(66)
and Bv := [Bv1 1, Bv2 1, ..., Bvne 1], with
BvI := ∇NI · (v − v m ) .

(67)

The last term can also be used to rewrite the L(v − v m ) term as
Lh (v h − v hm ) = Bv ve .

(68)

4.1.1.3 Weak form contribution of a fluid element
Given the above expressions, the contributions from element Ωe to the fluid weak form (15) can
be written as
e
GeF + GeG = weT fFe + qT
(69)
eg ,
with the (3ne × 1) FE force vector

e
f e + fFe int + fsupg
− fFe extf¯


 F in
−fFe extt̄
fFe :=



−fFe s

for Ωe ⊂ F h ,
for Ωe ⊂ ∂t F h ,

(70)

for Ωe ⊂ S h ,

and the (ne × 1) FE pseudo force vector
e
ge := gge + gpspg
.

13

(71)

They are composed of the FE forces
e ,
:= me ve0 + fcon
Z
ρ NT Bv ve dv ,
:=

fFe in
e
fcon

Ωe

fFe int
e
fsupg

fFe s
fFe extf¯
fFe extt̄

:= ce ve − de pe ,
Z
Z
T
τv BT
τv Bf (v − v m ) dv =
:=
v fres dv ,
Ωe
Ωe
Z
NT t da ,
:=
e
Ω
Z
NT f̄ dv
:=
e
Ω
Z
NT t̄ da ,
:=

(72)

Ωe

the FE pseudo forces
gge
e
gpspg

:= dT
e ve ,
Z
:=
τp GT fres dv ,

(73)

Ωe

and the elemental mass, damping and pressure-force matrices
Z
me :=
ρ NT N dv ,
e
Ω
Z
ce :=
BT C B dv ,
e
Ω
Z
de :=
DT Ñ dv .

(74)

Ωe

The tangent matrices of fFe and ge , needed for linearization, can be found in Appendix B.1.
e
Remark 4.1: One may simply change the sign of both gge and gpspg
in order to highlight the
e
e
symmetry between the second part of fF int and gg .

4.1.1.4 Stabilization terms
e
e
In order to evaluate the residual fres that appears in the stabilization terms fsupg
and gpspg
, we
note that
h + v h ) e = (G2 + H) v ,
2 div D h = (vj,ij
(75)
i
e
i,jj

where G2 := [G21 , G22 , ..., G2ne ], with


NI,11 NI,12 NI,13
G2I := ∇(∇NI ) = NI,21 NI,22 NI,23 
NI,31 NI,32 NI,33

(76)

and H := [H1 1, H2 1, ..., Hne 1], with
HI := tr G2I = NI,11 + NI,22 + NI,33 .

(77)

div σ h = η F ve − G pe ,

(78)

With this we can write

14

where F = G2 + H. Thus we obtain
h
fres ≈ fres
= ρ N ve0 + ρ Bv ve − η F ve + G pe − f̄ .

(79)

e
e
The stabilization parameters τv and τp appearing inside fsupg
and gpspg
are computed from

"
τv = τp =

2
∆t

2


+

2kvk
me he

2


+

4ν
me h2e

2 #− 21

(80)

(Shakib, 1988; Tezduyar, 1992; Rasool et al., 2016), where ∆t is the time step size, he is the
“element length” in the local flow direction taken from
n

e
1
1X
v
=
∇NI ·
he
2
kvk

(81)

I=1

(Tezduyar, 1992) and me depends on the polynomical order of the shape functions. I.e. for L1
(linear Lagrange) and L2 (quadratic Lagrange) elements we have me = 1/3 and me = 1/12,
respectively.7 According to this, parameters τv and τp are local parameters that change from
quadrature point to quadrature point.
4.1.1.5 Transformation of derivatives
In the above expressions
P ∇NI denotes the gradient w.r.t. the current configuration x, which is
h
discretized by x = I NI xmI , where xmI are the nodal positions of the FE mesh. Since it
is convenient to define the shape functions on a master element in ξ = [ξ, η, ζ]T space, where
∂NI /∂ξ is easily obtained, ∇NI needs to be determined from
∇NI =

∂NI
∂NI
= j −T
,
∂x
∂ξ

where

(82)

n

j=

e
∂xh X
∂NI
=
xmI ⊗
∂ξ
∂ξ

(83)

I=1

denotes the Jacobian of the mapping ξ → xm . Likewise, the second derivative G2I = ∇(∇NI )
is obtained from the formula
" n
#
e 
 ∂2N
2N
X
∂
I
J
G2I =
δIJ − ∇NI · xmJ
= j −T
j −1
(84)
∂x ∂x
∂ξ ∂ξ
J=1

that follows from differentiating (82). Eq. (84) is equavalent to the expression given in Dhatt
and Touzot (1984).
4.1.2

Membrane deformation

Following the notation of Eq. (52), the reference position and the current position within a
membrane element are approximated by the interpolations
X ≈ X h = N Xe ,
x ≈ xh = N xe ,
7

(85)

In Eqs. (80) and (81), v is taken from the previous time step in order to avoid the linearization of τv and τp .

15

where Xe and xe are arranged just like ve . From this follows
Aα ≈ Ahα = N,α Xe ,

(86)

aα ≈ ahα = N,α xe ,
where N,α := [N1,α 1, N2,α 1, ..., Nne ,α 1]. Likewise,
w,α ≈ wh,α = N,α we

(87)

follows from Eq. (53). Given Aα and aα , the metric tensor components Aαβ and aαβ can be
determined and the stress can be evaluated as discussed in Sec. 2.2.
Inserting the discretized expressions for v̇, aα , w and wα into the membrane weak form (34)
yields the elemental weak form contribution
GeS = weT fSe ,
with the (3ne × 1) FE force vector
( e
e
e − fe
fSin + fSint
+ fce − fSf
Sextf¯
e
fS :=
e
−fSextt̄

(88)

for Ωe ⊂ S h ,
for Ωe ⊂ ∂t S h ,

(89)

that is composed of
e
fSin
e
fSint

fce
e
fSf
e
fSext
f¯
e
fSext
t̄

Z

ρs NT N dv v̇e ,

:=
Ωe

Z
:=
Ωe

σ αβ NT
,α N,β da xe ,

Z

:= −
NT fc da ,
e
Ω
Z
:=
NT ff da ,
Ωe
Z
:=
NT f̄s da ,
e
Ω
Z
:=
NT
t t̄s ds .

(90)

Ωe

The discretization of the contact traction f c is straight forward, but an active set strategy
needs to be implemented in order to handle the state changes between contact and no contact (Wriggers, 2006). The tangent matrix of fSe , needed for linearization, can be found in
Appendix B.2.
4.1.3

Coupled system

Combining contributions (69) and (88) yields the coupled weak form
e
Ge = weT f e + qT
e g ,

with the (3ne × 1) FE force vector

e

f e + fFe int + fsupg
− fFe extf¯

 F in


 −f e
F extt̄
e
e
e
f := fF + fS =
e
e
e
e


 fSin + fSint + fc − fSextf¯


 −f e
Sextt̄
16

(91)

for Ωe ⊂ F h ,
for Ωe ⊂ ∂t F h ,
for Ωe ⊂ S h ,
for Ωe ⊂ ∂t S h .

(92)

It can be seen that for a conforming FE discretization of surface S, such as is considered here,
e of a membrane element cancels exactly
coupling condition (37) implies that the force vector fSf
e
e and f e
with fF s of the corresponding fluid boundary element. In the coupled system, both fSf
Fs
therefore do not appear anymore.
4.1.4

Double pressure nodes

Since the membrane is described here as a 2D surface that is discretized by 2D surface finite
elements, the membrane nodes carry a special role. Unless the membrane is located at the
boundary of the fluid, it is surround by fluid on both sides and generally supports pressure
jumps. A finite element node on S h therefore must carry two pressure dofs. One for each side
of the membrane. Otherwise, the formulation does not properly account for pressure jumps.
This is especially important for flexible membranes, where pressure jumps tend to become large.
In practice, each FE node on S h that is not located at boundary ∂S h (where both fluid sides
connect), is assigned two pressure dofs.8 When the elemental connectivity is then set up, care
has to be taken in order to connect the element on each side of S h with the correct dofs.
As long as a no-slip condition is considered on both sides of S, as is done here, the velocity field
is continuous across S and no extra velocity degrees of freedom are needed on S h .

4.2

Temporal discretization

The elemental force vectors f e and ge are assembled into the global vectors
f = fF in + fSin + fF int + fSint + fc + fsupg − fext

(93)

g = gg + gpspg ,

(94)

and
T , f T ]T , where
where fext := fF extf¯ + fF extt̄ + fSextf¯ + fSextt̄ . The former can be written as f = [fbr
r
fbr are the boundary reactions of the nodes on ∂x̂ F and ∂x S, and fr are the residual forces of
all the remaining nodes. Accordingly, the global residual vector
" #
fr
r :=
,
(95)
g

can be defined. The finite element forces are in equilibrium if r = 0. In general, r = 0 is a
coupled system of ordinary differential equations for the unknown nodal positions x := [xI ],
velocities v := [v I ], accelerations a := [v 0I ] (for fixed x) and pressures p := [pI ], for I = 1, ..., nno ,
that are all functions of time. The generalized-α scheme (Chung and Hulbert, 1993; Jansen
et al., 2000; Cottrell et al., 2009) is used to discretize r = 0 in time. Instead of solving for
the functions x(t), v(t), a(t) and p(t), the approximations xn ≈ x(tn ), vn ≈ v(tn ), an ≈ a(tn )
and pn ≈ p(tn ) are determined at discrete time steps tn , n = 0, ..., nt . This is based on the
Newmark update formulas for step tn → tn+1

∆t2
(1 − 2β) an + 2β an+1 ,
2

= vn + ∆t (1 − γ) an + γ an+1 ,

xn+1 = xn + ∆t vn +
vn+1
8

(96)

Tezduyar and Sathe (2007) propose to also use double pressure dofs at the boundary of ∂S h in order to
provide additional numerical stability.

17

where β and γ are non-dimensional parameters.9 According to the generalized-α scheme, r is
then evaluated for pn+1 and
xn+αf = xn + αf (xn+1 − xn ) ,
vn+αf = vn + αf (vn+1 − vn ) ,

(97)

an+αm = an + αm (an+1 − an ) ,
where 0 < αm ≤ 1 and 0 < αf ≤ 1 are chosen parameters.10 The global force vectors thus take
the form



f = fF in an+αm , vn+αf + fSin an+αm ) + fF int vn+αf , pn+1 + fSint xn+αf


+ fc xn+αf + fsupg an+αm , vn+αf , pn+1 − fext ,
(98)


g = gg vn+αf + gpspg an+αm , vn+αf , pn+1 .
The temporal inconsistency that is introduced if αm 6= αf 6= 1 is a deliberate feature of the
generalized-α method. The system r = 0 thus reduces to a system of algebraic equations that
can be solved for xn+1 , vn+1 , an+1 and pn+1 given the previous values xn , vn , an and pn .
One option is to pick u := [v, p] as the primary unknowns, solve r = 0 for un+1 , and then
obtain an+1 and xn+1 (which is really only needed for the membrane nodes) from (96). Since
the system r = 0 is non-linear, the Newton-Raphson method is used. This requires the tangent
matrix k that is assembled from the elemental entries
∂re
ke :=
.
(99)
∂un+1
e
It is given in Appendix C for the considered fluid and membrane elements. In the following
computations, the Newmark parameters are taken as (Chung and Hulbert, 1993)
1
− αf + αm ,
2
2
1
β =
1 − αf + αm
4
γ =

(100)

using the generalized-α parameters11
αf =

2
,
3

αm =

5
.
6

(101)

This choice ensures second order accuracy in time and unconditional stability (for linear problems).

4.3

Normalization

In order to implement the above expressions within a computer code they have to be normalized.
We therefore chose a length scale L0 , time scale T0 and force F0 , and use those to normalize
all lengths, times and forces in the system. Velocities, masses, fluid densities, fluid viscosities,
fluid pressures, membrane densities and membrane stresses are then normalized by the scales
v0 :=

L0
,
T0

m0 :=

F0 T02
,
L0

ρ0 :=

m0
,
L30

η0 :=

9

F0 T0
,
L20

p0 :=

F0
,
L20

ρs0 :=

m0
,
L20

γ0 :=

F0
.
L0
(102)

They should not be confused with the physical parameters β and γ used for the surface inclination and
surface tension in other sections.
10
Note that the α introduced by Chung and Hulbert (1993) corresponds to 1 − α here.
11
They are obtained taking a spectral radius of ρ∞ = 12 for the first order system, see Jansen et al. (2000).

18

System (98) can then be expressed in the normalized form
f̄ (ūn+1 ) = f̄F in + f̄Sin + f̄F int + f̄Sint + f̄c + f̄supg − f̄ext ,
(103)

ḡ(ūn+1 ) = ḡg + ḡpspg ,
where a bar denotes normalization with the corresponding scale from above, e.g.
e ,
f̄Fe in = m̄e āe + f̄con

with

Z
m̄e :=

ρ̄ NT N dv̄ ,

Ω̄e
e
f̄con
:=

Z

(104)

(105)
ρ̄ NT B̄v v̄e dv̄ ,

Ω̄e

and ρ̄ = ρ/ρ0 , dv̄ = dv/L30 , B̄v = Bv T0 , v̄e = ve /v0 and āe = ae T0 /v0 . All the other quantities
appearing in (103) are normalized in the same fashion. Solving (103) then gives the normalized
unknowns v̄ = v/v0 and p̄ = p/p0 , while (96) can be solved for x̄ = x/L0 and ā = a T0 /v0 .

4.4

Mesh motion

Apart from the unknown material velocity v and pressure p, the discrete mesh velocity vm
can also be regarded as an unknown. In that case suitable (differential) equations have to be
formulated for vm . A simpler approach is to determine the mesh velocity from the membrane
velocity using linear interpolation: On the membrane surface the mesh motion is considered
Lagrangian, i.e. vm = v, whereas it is treated Eulerian (vm = 0) beyond a certain distance from
the membrane. In-between, simple linear interpolation is used. Details of this are reported in
the following examples. Linear interpolation, and ALE in general, does not work for some FSI
problems. An example are solids revolving within the fluid. For such cases, other techniques
need to be considered.

5

Numerical examples

This section presents three numerical examples that range form very low to quite large Reynolds
numbers. The first example considers a solid membrane (with no bending resistance), the second
example considers a liquid membrane, and the third example considers a solid shell with low
bending resistance.

5.1

Fluid-inflated cylinder

The first numerical example considers the radial inflation of a cylindrical membrane due to
radial inflow. The numerical solution will be compared to the analytical solution derived in
Sec. 3.1. The initial inner radius of the cylinder Rin , the maximum inflow velocity v0 and the
fluid density ρ are used for normalization, such that L0 = Rin , T0 = Rin /v0 and ρ0 = ρ. The
outer radius of the membrane at initialization time t = 0 is taken as Rs = 2L0 . Computationally,
only a quarter of the cylindrical domain is modelled with a chosen height of H = L0 . Sliding
wall conditions12 are applied to all fluid boundaries except the inflow and the membrane surface,
12

The normal velocity and the tangential traction are set to zero.

19

where coupling conditions apply, and the inflow boundary, where the radial inflow velocity
(

1 − cos(πt/T0 ) /2 for t < T0
vin (t) = v0
(106)
1
else
is prescribed. The Reynolds number, Re = ρ vin L0 /η, is chosen as Re = 100 guaranteeing a
purely laminar flow. For water at room temperature ρ ≈ 1000 kg/m3 and η = 1.00 mNs/m2 this
implies v0 = 10 m/s. The membrane is modelled as an incompressible Neo-Hookean, rubber-like
material according to (29). The membrane’s nondimensional shear stiffness is taken as µ̄ = 0.1.
The fluid domain is discretized by Nf = nr × nθ × 1 quadratic volume elements in er , eθ and
e3 direction (see Fig. 1), while the membrane domain is discretized by Ns = nθ × 1 quadratic
surface elements along eθ and e3 . Tab. 1 shows the considered meshes. The time step is chosen
total elements

fluid elements

membrane elements

nodes

dofs

7

6×1×1

1×1

117

495

42

13 × 3 × 1

3×1

567

2,331

100

24 × 4 × 1

4×1

1,323

5,373

Table 1: Fluid-inflated cylinder: Considered FE meshes based on quadratic Lagrange elements.
as ∆t̄ = 0.0025 for all cases. The radial mesh velocity at time step tn+1 is defined by the linear
interpolation

R − Rin
(107)
vm R, tn+1 =
vs (tn ) ,
Rs − Rin
where vs (tn ) is the cylinder’s radial velocity at the previous time step.
Fig. 2 shows the radial flow field and the membrane displacement due to the cylinder inflation
at different time steps. The solid membrane is stretched by more than a factor of 3. For the

(a) t̄ = 0

(b) t̄ = 1

(c) t̄ = 6

(d) t̄ = 11

(e) t̄ = 21

Figure 2: Fluid-inflated cylinder: Radial flow field v̄ = v/v0 and cylinder expansion at various
time steps. Computationally, only a quarter of the system is modelled.
membrane displacement (Fig. 3) and velocity (Fig. 4) the numerical result is in perfect agreement
with the analytical solution derived in Sec. 3.1; see Eqs. (39) & (40). For the pressure shown
in Fig. 5 we observe deviations from the analytical result (43) during the transient part and
again nearly perfect agreement at the final simulation time. The numerical results improve
for a higher mesh resolution. The finite element discretization and its implementation shows
quadratic convergence behavior as expected, see Fig. 3b.

5.2

Rolling droplet

The second example simulates rolling contact of a liquid droplet on an inclined substrate considering a low Reynolds number and a contact angle of 180◦ . As we expect the motion to
20

7

1

6.5

-1

10

5.5

10-2

5

10-3

4.5

e

-r

6

4

10-4

3.5

10-5

3

10-6

2.5
2
0

5

10

15

20

10-7

25

1

10
n

-t

(a) r(t) at R = Rs

100

(b) Convergence

Figure 3: Fluid-inflated cylinder: (a) Membrane position r̄ = r/L0 vs. time t̄ = t/T0 . (Analytical result: green ×, FE solution: red +). (b) Numerical error (L2 -norm) against total number
of L2 elements (r: red +, v: green ×, a: blue ?, p: magenta ). The dash-dotted line marks
quadratic convergence behavior.

0.9

0.35

0.8

0.3

0.7

0.25

0.6
v-

1

0.4

v-

0.45

0.2

0.5

0.15

0.4

0.1

0.3

0.05

0.2

0
0

5

10

15

20

0.1
25

1

2

3

-t

(a) v(t) at r = rs

4
-r

5

6

7

(b) v(t) at t̄ = 21

0.05
0.045
0.04
0.035
0.03
0.025
0.02
0.015
0.01
0.005
0
-0.005

0.05
0
-0.05
-0.1
-0.15
p-

p-

Figure 4: Fluid-inflated cylinder: (a) Normalized membrane velocity vs. time; (b) Normalized
fluid velocity vs. radial position at t = 21 T0 . (Analytical result: green ×, FE solution: red +)

-0.2
-0.25
-0.3
-0.35
-0.4

0

5

10

15

20

-0.45

25

-t

(a) p(t) at r = rs

1

2

3

4
-r

5

6

7

(b) p(t) at t̄ = 21

Figure 5: Fluid-inflated cylinder: (a) Normalized membrane pressure vs. time; (b) Normalized
fluid pressure vs. radial position at t = 21 T0 . (Analytical result: green ×, FE solution: red +)
come close to the spinning solution of Sec. 3.2, a purely Lagrangian FE description is chosen
(v m = v). This also allows to use a classical contact description between droplet and substrate.
21

There is earlier computational work on rolling droplets (Rasool et al., 2012; Li et al., 2013;
Thampi et al., 2013; Wind-Willassen and Sørensen, 2014). But it is either 2D, or non-FE. So
the present study seems to be the first 3D FE simulation of rolling droplets. Novel is also the
way contact is treated here – by using a computational contact algorithm with an active-set
strategy. Within that, a no-slip (sticking) condition is assumed on the contact surface. If slip
occurs, a stick-slip algorithm is needed for the droplet (Sauer, 2016).
The droplet setup considers similar parameters as in Sauer (2016): An initially spherical droplet
with radius R = L0 and volume V = 4πL30 /3 is considered under gravity loading, such that
ρgL30 = γL0 . For water at room temperature, with ρ = 1000 kg/m3 , g = 9.81 m/s2 and
γ = 72.8 mN/m, this corresponds to a droplet with L0 = 2.72 mm and V = 84.6 µl. For
further normalization we choose g0 = g and γ0 = γ, so that T0 = 16.7 ms, F0 = 0.198 mN and
p0 = 26.7 Pa. A high fluid viscosity is chosen, i.e. η = 11.9 Ns/m2 , such that the Reynolds
number becomes very small. A suitable definition for the Reynolds number of a rolling droplet
is
ρ Lc vmean
Re =
,
(108)
η
where Lc is the diameter of the contact surface and vmean is the mean droplet velocity. The
penalty parameter for sticking according to contact model (33) is taken as c = 250 m2 p0 /L0 ,
where m characterizes the FE resolution according to Tab. 2. Quadratic Lagrange elements are
m

fluid elements

membrane elements

nodes

dofs

2

128

48

1,241

4,964

4

832

192

7,407

29,628

8

6,656

768

56,157

224,628

16

53,248

3,072

437,433

1,749,732

Table 2: Rolling droplet: Considered FE meshes based on quadratic Lagrange elements.
used.
Initially the droplet is at rest. Rolling motion is then induced by inclining the substrate considering the time-varying inclination angle

πt

for 0 ≤ t < t1 ,
1 − cos


t1



2
for t1 ≤ t ≤ t2 ,
β0
β(t) =
(109)
π(t − t2 )
2 

1
+
cos
for
t
≤
t
≤
t
+
t
,

2
1
2

t1


0
for t1 + t2 < t ≤ t3 ,
with t1 = 50 T0 , t2 = 200 T0 , t3 = 350 T0 and the two cases:
1. β0 = 10◦ with ∆t = 8 T0 /m, and
2. β0 = 20◦ with ∆t = 4 T0 /m.
Fig. 6 shows the finite element results for the mean droplet velocity vmean for the two cases.13 As
seen the FE results converge upon mesh refinement. The figure also shows that steady rolling
motion is attained at about t = 150 T0 for β0 = 20◦ , while it is attained almost instanteneously
for β0 = 10◦ (i.e. at t = t1 ). The instantaneous response of vmean on β, for low β0 , can be
also seen from the vmean (β)–plot in Fig. 7. Both branches (for increasing β and decreasing
β, respectively) are almost identical. For β0 = 20◦ on the other hand the two branches are
13

The mean droplet velocity vmean is determined by computing the volume average of the fluid velocity v and
then taking its component parallel to the substrate surface.

22

0.03
0.029
mean velocity [L0/T0]

β0 = 20 deg.

0.025
mean velocity [L0/T0]

m=2
m=4
m=8

0.02
0.015
0.01

β0 = 10 deg.

m=2
m=4
m=8
m = 16

0.0285

0.028

0.0275

0.005
0.027
0
0

50

100

150
200
time t [T0]

250

300

150

160

170
180
time t [T0]

190

200

Figure 6: Rolling droplet: Mean droplet velocity vs. time for β0 = 10◦ and β0 = 20◦ using the
meshes from Tab. 2. The right hand side shows an enlargment for β0 = 20◦ . As seen, the FE
results converge upon mesh refinement.
−3

0.03

x 10

0.025
mean velocity [L0/T0]

mean velocity [L0/T0]

10
8
6
4
2

a.

0.015
0.01
0.005

0
0

0.02

0
2

4
6
inclination β [deg.]

8

10

0

b.

5

10
15
inclination β [deg.]

20

Figure 7: Rolling droplet: Mean droplet velocity vs. β for β0 = 10◦ (a) and β0 = 20◦ (b) using
m = 16. The return branch (for decreasing β) is marked by a dashed line.

different.
For further illustration, Fig. 8 shows the droplet deformation and velocity field kvk during
rolling. The deformation is considerable and should not be neglected, as has been done in earlier work (Rasool et al., 2012, 2013). The figure also shows how the contact surface changes.
Initially the contact surface is circular with a diameter of Lc = 1.36 L0 . During steady rolling
the diameter in rolling direction reduces to Lc = 1.04 L0 . Since vmean = 0.0268 L0 /T0 , the
Reynolds number thus becomes Re = 1.04 · 10−3 according to (108). Fig. 8 clearly shows that
the advancing and receding droplet halves are not symmetric during rolling.
This can also be seen from the pressure distribution shown in Fig. 9. The fluid pressure is largest
at the advancing front of the contact surface. Since the contact surface is flat, the fluid pressure
is equal to the contact pressure. Close inspection shows that the pressure is oscillatory in the
vicinity of the contact line C. Those oscillations do not converge with mesh refinement, as the
velocity field does. So it seems that the pressure stabilization scheme, described in Sec. 2.1.4, is
not sufficient to handle the contact boundary of a rolling droplet, even though the static droplet
(at t = 0 and t = 350 T0 ) poses no problem. The problem may be related to the discontinuity of
the contact pressure: it jumps to zero at the contact boundary. The way the fluid velocity, fluid
23

Figure 8: Rolling droplet: Velocity magnitude kvk/v0 at t = 0, t = 50 T0 , t = 100 T0 , t = 200 T0
and t = 350 T0 (left to right) for β0 = 20◦ and m = 8. Only half of the symmetric droplet is
shown. In the top panel the symmetry surface is removed and instead a selected material plane
is tracked during deformation. A single fluid particle is marked by ‘◦’.

Figure 9: Rolling droplet: pressure field p/p0 at t = 0, t = 50 T0 , t = 100 T0 , t = 200 T0 and
t = 350 T0 (left to right) for β0 = 20◦ and m = 8

pressure and contact pressure are interpolated (quadratic Lagrange interpolation is used here)
seem incompatible. It seems that this problem has not yet been addressed in the literature.
Further study is required on the topic. Perhaps C 1 -continuous interpolation, such as is provided
by NURBS, would help. We note that for β = 10◦ , pressure oscillations also appear, but they
are less pronounced.
To remove the pressure oscillations, Gaussian smoothing can be used for post-processing. Selecting the variance of the Gaussian distribution as σ = 1/m, i.e. on the order of the nodal
distance, gives non-oscillatory pressures; see Fig. 10. The smoothed pressure converges with
mesh refinement. The pressure distribution shows that the advancing contact surface carries
most of the droplet weight (component cos β × ρgV ). Component sin β × ρgV is equilibrated
by a tangential sticking force. The moment caused by these external forces is equilibrated by
the internal moment of the fluid stress.
The last plot shows the vorticity (i.e. spin) component 2ω2 := e2 ·(∇×v) (along the axis of rotation e2 ) and the dissipation D = σ : D during rolling; see Fig. 11. Also here smoothing is used.
According to Sec. 3.2 the vorticity of a spinning sphere is a constant vector with magnitude
2ω. In contrast, the vorticity of a rolling droplet is non-constant: A maximum is attained at
the contact boundary and a minimum occurs on the contact surface. Although, away from the
contact surface, the vorticity approaches a constant. The behavior is similar for the dissipation:
Away from the contact surface, the dissipation is zero and thus agrees with the spinning sphere
solution. Non-zero dissipation, associated with shear flow, occurs in the vicinity of the contact
surface, with a maximum occuring at the advancing contact front. For longer rolling droplets,
or for higher β, the shear flow becomes more pronounced, such that an ALE formulation is
24

Figure 10: Rolling droplet: smoothed pressure field at t = 0, t = 50 T0 , t = 100 T0 , t = 200 T0
and t = 350 T0 (left to right) for β0 = 20◦ and m = 8. See also supplementary movie file
drop roll p.mpg.

a.

b.

Figure 11: Rolling droplet: a. smoothed vorticity component 2ω2 at t = 50 T0 and t = 200 T0 ;
b. smoothed dissipation D = σ : D at t = 50 T0 and t = 200 T0 ; both for β0 = 20◦ and m = 8.
The units of 2ω2 are 1/T0 ; the units of D are p0 /T0 .

needed for the mesh. On the free surface (which is tracked explicitly within the present scheme)
such a formulation needs to be Lagrangian in the normal direction but Eulerian in-plane. The
formulation of such an ALE scheme is outside the present scope.

5.3

Flapping flag

The third example simulates the flapping motion of a flag. The problem setup of this example
is shown in Fig. 12. The flag is modeled as a flexible sheet that is supported on the left hand
side. It is excited by a uniform inflow with velocity vin . The length scale L0 , the fluid density
ρ0 and the time scale T0 are used to normalize the problem. The remaining parameters are
chosen according to Tab. 3. Considering L0 = 0.1m, T0 = 1s and ρ0 = 1.2 kg/m3 , the fluid
parameters become ρ = ρ0 and η = 18.37 µNs/m2 , which correspond to the values of air at sea
level and 20◦ C, while the flag parameters become ρs = 0.12 kg/m2 , µ = 5 N/m and c = 0.12 µNm
according to Sec. 4.3.14 The Reynolds number of the problem is
Re =

ρ Lc vin
,
η

(110)

where Lc is the chord length of the flag. For Lc = 3L0 and the considered ρ and η follows
Re = 1960 v̄in . At this Re and density ratio15 , the flag motion can be expected to be chaotic
14
15

Following Sec. 4.3, the bending stiffness needs to be normalized by c0 = F0 L0 , where F0 = ρ0 L40 /T02 .
The density ratio R1 := ρs /(ρLc ), as defined in Shelley and Zhang (2011), is 1/3 here.

25

Figure 12: Flapping flag: Side, top and front view of the problem setup. The flag is fixed on
the left and its lateral displacement and velocity are monitored at point A.
parameter

normalized value

inflow velocity
density of the fluid
viscosity of the fluid
density of the flag
shear stiffness of the flag
bending stiffness of the flag

v̄in = 1
ρ̄ = 1
η̄ = 1.531 · 10−3
ρ̄s = 1
µ̄ = 4.167 · 103
c̄ = 0.01

Table 3: Flapping flag: Considered inflow and material parameters.
according to the phase diagram of Connell and Yue (2007).
The flapping flag example is a good test case since the flag motion and the surrounding flow field
can become very complex, as the experimental data reported in Shelley and Zhang (2011) show.
There have been recent 3D simulations that study the problem in detail (Hoffman et al., 2011;
Banerjee et al., 2015; Gilmanov et al., 2015; de Tullio and Pascazio, 2016). In some of those
works immersed boundary methods are used instead of ALE. Such methods are advantageous
for very large flag motions that may even involve self-contact. In contrast to earlier work, the
flag is discretized here with C 1 -continuous isogeometric shell elements. Their formulation is
e
the same as the one of Eq. (90) with the only exception that fSint
is extended by the internal
bending moments according to the formulation of Duong et al. (2017) using the Canham bending model. A shell formulation is used in order to regularize the system with bending stiffness.
A low stiffness value is used such that the structure remains very flexible. Below a certain
threshold value of c, the flapping behavior becomes independent of c as is shown later.
The fluid domain is discretized with nF el = 8 × 2 × 4 m3 quadratic 3D NURBS elements, while
the flag is discretized with nSel = 3 × 2 m2 quadratic 2D NURBS elements. The number of
nodes and dofs resulting from this discretization16 is listed in Tab. 4. On the surface of the flag,
double pressure dofs are used to account for pressure jumps as described in Sec. 4.1.4. The time
step is taken as ∆t = 0.16 T0 /m.
Fig. 13 shows the flag deformation at selected time steps. Those are snap-shots of the supplementary movie file flag v.mpg. As expected, the structure performs flag-typical oscillations
along its length. Close inspection shows that the flag motion also varies in vertical direction.
16

The number of nodes is nno = (8m + 4)(2m + 3)(4m + 4); the number of dofs is ndof = 4nno + nSel , due to
the double pressure nodes.

26

m

fluid elements

membrane elements

nodes

dofs

2

512

24

1,680

6,744

4

4096

96

7,920

31,776

8

32,768

384

46,512

186,432

Table 4: Flapping flag: Considered FE meshes based on quadratic NURBS elements.

Figure 13: Flapping flag: Deformation at t = 44.36 s, t = 44.92 s, t = 45.48 s, t = 46.04 s
and t = 46.60 s (left to right) for m = 8; see also supplementary movie file flag v.mpg. The
coloring shows the lateral velocity component in the range {−1, 1}v0 (from blue to red). The
streamlines of the flow are also shown.

The pressure field around the flag is shown in Fig. 14. The figure also shows the mesh motion

Figure 14: Flapping flag: Fluid pressure in the mid-plane at t = 44.36 s, t = 44.92 s, t = 45.48 s,
t = 46.04 s, t = 46.60 s and t = 47.16 s (top left to bottom right) for m = 8. The coloring is in
the range {−.7, 1.2}p0 (from blue to red).

around the flag. It is based on the interpolation scheme given in App. D.
For the chosen parameters, the flapping behavior is still (quite) periodic, as Fig. 15 shows. The
period of the main oscillation is 5.60 s. Apart from the main oscillations, there are also fine
scale oscillations, as Fig. 15b shows. Fig. 15 also shows that the simulation results converge
with mesh refinement. For the first 20 seconds, mesh m = 4 already gives quite good results.
The model parameters of Tab. 3 affect the flapping behavior of the flag. The influence of Re
has been discussed in detail in earlier work, e.g. see Shelley and Zhang (2011), so the following
discussion focuses on the membrane parameters. Apart from their usual influence on structural
vibrations, two aspects are noteworthy: 1. For sufficiently low c, the flapping behavior (for given
Re) remains unchanged, i.e. it becomes independent of c. According to Fig. 16a this occurs
below c̄ ≈ 10−3 . Below that c, the flag is effectively a membrane without bending stiffness, and
c is only helpful for regularizing the numerical solution. 2. Increasing µ leads to increased fine
scale oscillations, as Fig. 16b shows. Since µ controls the in-plane stiffness of the flag, those
oscillations can be associated with longitudinal vibrations of the flag.

27

1
0.8

0.8

0.4

flag velocity [L0/s]

flag displacement [L0]

0.6

1.2
m=2
m=4
m=8

0.2
0
−0.2
−0.4
−0.6

m=2
m=4
m=8

0.4
0
−0.4
−0.8

−0.8
−1
0

10

20
30
time t [s]

a.

40

−1.2
0

50

10

b.

20
30
time t [s]

40

50

Figure 15: Flapping flag: Lateral displacement (a) and velocity (b) at point A for various FE
discretizations. Symbol ‘◦’ marks the configurations shown in Fig. 13.

flag displacement [L0]

0.6
0.4

c = 0.1
c = 0.04
c = 0.01
c = 1E−3
c = 1E−4
c = 1E−6

0.8
0.6

0.2
0
−0.2

0.2
0
−0.2

−0.4

−0.4

−0.6

−0.6

−0.8

−0.8

0

5

10

µ = 20835
µ = 4167
µ = 2084

0.4
flag velocity [L0/s]

0.8

15

a.

20
25
time t [s]

30

35

40

0

b.

10

20
30
time t [s]

40

50

Figure 16: Flapping flag: Influence of membrane parameters c (a) and µ (b). The influence of
c vanishes below a threshold value of c. Increasing µ leads to smaller velocities but increased
fine scale oscillations.

6

Conclusion

A unified FSI formulation is presented that is suitable for solid, liquid and mixed membranes.
At free liquid surfaces, sticking contact can be accounted for. The fluid flow and the structure
are discretized with finite elements using a stabilized fluid formulation and a surface-based membrane formulation. A conforming interface discretization is used between fluid and membrane,
which leads to a simple monolithic coupling formulation. On the membrane surface double
pressure nodes are required. The temporal discretization is based on the generalized-α scheme.
Two analytical and three numerical examples are presented in order to illustrate and verify the
proposed formulation. They consider fluid flow at low and high Reynolds numbers.
The proposed formulation is very general and thus suitable as a basis for further research. In
order to increase efficiency, the formulation can be extended to boundary elements (for low Re)
or turbulence models (for high Re). Under current study is the use of enriched finite element
discretizations (Harmel et al., 2017) that are suitable to efficiently capture boundary layers (Rasool et al., 2016). Another extension of the present formulation is to re-examine the pressure
stabilization scheme at contact boundaries. This would be especially important in the presence
28

of sharp contact angles. Such a formulation would then allow for a detailed flow analysis of
droplets on rough surfaces.

Acknowledgements
The authors are grateful to the German Research Foundation (DFG) for supporting this research
under grants GSC 111 and SA1822/3-2. The authors also wish to thank Maximilian Harmel
and Raheel Rasool for proofreading the manuscript.

A

Uniform membrane stretch

For the analytical example of Sec. 3.1, the initial and the current membrane position are described by
X(θ, z) = Rs er + z e3 ,
(111)
x(θ, z) = rs er + z e3 .
From this follows A1 = Rs eθ , a1 = rs eθ and A2 = a2 = e3 with eθ = − sin θ e1 + cos θ e2 . We
further find
 2 
 2 
r 0
Rs 0
,
,
[aαβ ] = s
[Aαβ ] =
0 1
0 1
(112)
 −2 
 −2 
 αβ 
 αβ 
rs
0
Rs
0
,
, a
=
A
=
0 1
0
1
such that J = rs /Rs =: λ and




1
−4
 αβ  µ
1−λ
0 
σ
=  Rs2
.
λ
0
1 − λ−2

(113)


The stress component along eθ is σ := eθ · σ αβ aα ⊗ aβ eθ , which yields expression (44).

B
B.1

FE tangent matrices for the time-continuous system
Fluid element

In order to evaluate the tangent matrix of the finite element force vector fFe defined in (70), we
require
∂ Bv ve )
(114)
Lv :=
∂ve
which can be written as Lv := [Lv1 , Lv2 , ..., Lvne ] with

Therefore

LvI := NI L + BvI 1 .

(115)

h
∂fres
= ρ Lv − η F .
∂ve

(116)

29

Based on this, we find the tangent matrices of the fluid forces defined in (72)-(74)
meF
mesupg
mepspg
cecon

:=
:=
:=
:=

ceF int :=
cesupg

:=

ceg

:=

cepspg

:=

deF int :=
desupg :=
depspg

:=

∂fFe in
∂ve0
e
∂fsupg
∂ve0
e
∂gpspg
∂ve0
e
∂fcon
∂ve
∂fFe int
∂ve
e
∂fsupg
∂ve
∂gge
∂ve
e
∂gpspg
∂ve
∂fFe int
∂pe
e
∂fsupg
∂pe
e
∂gpspg
∂pe

= me ,
Z
=
Ωe

Z

τv ρ BT
v N dv ,
τp ρ GT N dv ,

=
Ωe

Z

ρ NT Lv dv ,

=
Ωe

= ce ,
Z
=
Ωe


T
T
τv BT
f N + ρ Bv Lv − η Bv F dv ,

(117)

= dT
e ,
Z


τp ρ GT Lv − η GT F dv ,

=
Ωe

= − de ,
Z
=
Ωe

Z

τv BT
v G dv ,
τp GT G dv .

=
Ωe

e
e
As seen, a major source of complexity are the stabilization terms fsupg
and gpspg
.

B.2

Membrane element

Linearizing the membrane forces in (91) w.r.t. v̇e and xe yields the mass matrix
Z
e
∂fSin
e
mS :=
=
ρs NT N dv ,
∂ v̇e
e
Ω

(118)

and the stiffness matrix
keS := keSint + kec .
The first term of

keS

(119)

follows from Sauer et al. (2014) as
e
∂fSint
= kegeo + kemat ,
∂xe

keSint :=
with
kegeo :=
kemat

Z
Ωe0

Z

τ αβ NT
,α N,β dA ,
(121)
αβγδ

:=

c
Ωe0

(120)

NT
,α (aβ

τ αβ := Js σ αβ and
cαβγδ := 2
30

⊗ aγ ) N,δ dA ,

∂τ αβ
.
∂aγδ

(122)

Here,
cαβγδ =


2µ αβ γδ
αγ βδ
αδ βγ
a
a
+
a
a
+
a
a
Js2

(123)

for model (29) and
cαβγδ = Js γ aαβ aγδ − aαγ aβδ − aαδ aβγ



(124)

for model (30), see Sauer et al. (2014) and Sauer and Duong (2017). Inserting these into (120),
yields the simpler expression
Z

2µ T α
e
kmat :=
N,α a ⊗ aβ + aβ ⊗ aα + aαβ i N,β dA
(125)
2
Ωe0 Js
for model (29) and
kemat

Z
:=
Ωe0


α
β
β
α
αβ
Js γ NT
,α a ⊗ a − a ⊗ a − a i N,β dA

(126)

for model (30). Here i := aγ ⊗ aγ is the identity tensor on surface S. With this, keSint can be
further simplified, in particular for model (30), see Sauer (2016).
The second term of keS depends on the contact description. Here, sticking contact is considered
with a rigid substrate using the penalty regularization of Eq. (33). For this case, we have
Z
Z

∂fce
T ∂fc
e
N da −
NT fc ⊗ aα N,α da ,
kc :=
=−
N
(127)
∂xe
∂x
Ωe
Ωe
with
∂fc
=
∂x

(

− 1 if g · nc < 0 ,
0

(128)

else .

The front term of kec follows directly from Eqs. (33) and (85), while the rear term is derived in
Sauer and De Lorenzis (2015).

C
C.1

FE tangent matrices for the time-discrete system
Fluid element

For a fluid element Ωe ⊂ F h , the tangent matrix ke defined by Eq. (99) is given by the 4ne ×4ne
matrix
 ∂f e
∂f e 
 ∂ven+1 ∂pn+1

e
,
ke = 
(129)
e
 ∂ge
∂g 
∂ven+1 ∂pn+1
e
with

∂f e
∂f e
=
m
∂an+α
∂ven+1
e
∂ge
∂ge
n+1 =
m
∂an+α
∂ve
e

m
∂an+α
e
∂an+1
e
n+α
∂ae m
∂an+1
e

∂an+1
∂f e
e
+
∂ven+1 ∂ven+αf
∂an+1
∂ge
e
n+1 +
∂ve
∂ven+αf

31

∂ven+αf
,
∂ven+1
∂ven+αf
.
∂ven+1

(130)

Based on (98), (117), (96) and (97) we obtain
∂f e
∂ven+1
∂f e
∂pn+1
e
∂ge
∂ven+1
∂ge
∂pn+1
e

=



αm
meF + mesupg + αf cecon + ceF int + cesupg ,
γ ∆t

= deF int + desupg ,
(131)


αm
=
mepspg + αf ceg + cepspg ,
γ ∆t
= depspg ,

where the individual building blocks are given in (117).

C.2

Membrane element

For a membrane element Ωe ⊂ S h , the tangent matrix ke defined by Eq. (99) is given by the
3ne × 3ne matrix
ke =

m ∂an+1
f ∂xn+1
∂f e ∂an+α
∂f e ∂xn+α
∂f e
e
e
e
e
=
+
n+αf
n+1
n+1
n+1
n+1 .
m
∂an+α
∂ven+1
∂a
∂v
∂x
∂v
∂x
e
e
e
e
e
e

(132)

Based on (98), (117), (96) and (97) we find
ke =

αm
αf β ∆t e
meS +
kS ,
γ ∆t
γ

(133)

where meS and keS are given in Appendix B.2.

D

Mesh motion for flapping flag example

For the flapping flag example in Sec. 5.3, the mesh velocity v m (with cartesian components vmi )
at FE node (i.e. control point) xm is defined by the linear interpolation

for Xm1 ≤ L0 or Xm1 > 6L0 ,
 0


1
for L0 < Xm1 < 4L0 ,
vm1 xm , t = vs1 X p , t
(134)


3 − Xm1 /(2L0 ) for 4L0 < Xm1 < 6L0 ,
for the inflow direction, and
(
vmj



xm , t = vsj X p , t

0

for d ≥ L0 ,

1 − d/L0 for d < L0 ,

(135)

for the other directions (j = 2, 3). Here Xmi are the components of X m = xm t=0 , d(X m ) is the
distance of X m from the flag surface, and v s (X p , t) is the current flag velocity at the initially
nearest membrane gridpoint X p = X p (X m ). Note that vm1 is smooth at Xm1 = 1 since vs1
approaches 0 smoothly as Xm1 → L0 .

32

References
Banerjee, S., Connell, B. S. H., and Yue, D. K. P. (2015). Three-dimensional effects on flag
flapping dynamics. J. Fluid Mech., 783:103–136.
Barett, J. W., Garcke, H., and Nürnberg, R. (2015). Numerical computations of the dynamics
of fluidic membranes and vesicles. Phys. Rev. E., 92:052704.
Bazilevs, Y. and Takizawa, K., editors (2016). Advances in Computational Fluid-Structure
Interaction and Flow Simulation. Springer, Berlin, Heidelberg.
Bazilevs, Y., Takizawa, K., and Tezduyar, T. E. (2013). Computational Fluid-Structure Interaction: Methods and Applications. John Wiley & Sons.
Bletzinger, K.-U., Wüchner, and Kupzok, A. (2006). Algorithmic treatment of shells and free
form-membranes in FSI. In Bungartz, H. and Schäfer, M., editors, Fluid-Structure Interaction.
Lecture Notes in Computational Science and Engineering, vol. 53, pages 336–355, Berlin,
Heidelberg. Springer.
Brooks, A. N. and Hughes, T. J. R. (1982). Streamline upwind/Petrov-Galerkin formulations
for convection dominated flows with particular emphasis on the incompressible Navier-Stokes
equations. Adv. Appl. Mech., 32:199–259.
Chung, J. H. and Hulbert, G. M. (1993). A time integration algorithm for structural dynamics
with improved numerical dissipation: the generalized-α method. J. Appl. Mech., 60:371–375.
Connell, B. S. H. and Yue, D. K. P. (2007). Flapping dynamics of a flag in a uniform stream.
J. Fluid. Mech., 581:33–67.
Cottrell, J. A., Hughes, T. J. R., and Bazilevs, Y. (2009). Isogeometric Analysis. Wiley.
Cristini, V. and Tan, Y.-C. (2004). Theory and numerical simulation of droplet dynamics in
complex flows–a review. Lab Chip, 4:257–264.
de Tullio, M. D. and Pascazio, C. (2016). A moving-least-squares immersed boundary method
for simulating the fluid-structure interaction of elastic bodies with arbitrary thickness. J.
Comput. Phys., 325:201–225.
Dhatt, G. and Touzot, G. (1984). The Finite Element Method Displayed. Wiley.
Donea, J. and Huerta, A. (2003). Finite Element Methods for Flow Problems. Wiley.
Dowell, E. H. and Hall, K. C. (2001). Modeling of fluid-structure interaction. Ann. Rev. Fluid
Mech., 33:445–490.
Duong, T. X., Roohbakhshan, F., and Sauer, R. A. (2017). A new rotation-free isogeometric thin
shell formulation and a corresponding continuity constraint for patch boundaries. Comput.
Methods Appl. Mech. Engrg., 316:43–83.
Franci, A., Oñate, E., and Carbonell, J. M. (2016). Unified Lagrangian formulation for solid
and fluid mechanics and FSI problems. Comput. Methods Appl. Mech. Engrg., 298:520–547.
Gilmanov, A., Le, T. B., and Sotiropoulos, F. (2015). A numerical approach for simulating
fluid structure interaction of flexible thin shells undergoing arbitrarily large deformations in
complex domains. J. Comput. Phys., 300:814–843.

33

Harmel, M., Sauer, R. A., and Bommes, D. (2017). Volumetric mesh generation from T-spline
surface representations. Comput. Aid. Des., 82:13–28.
Heltai, L., Kiendl, J., DeSimone, A., and Reali, A. (2016). A natural framework for isogeometric fluid-structure interaction based on BEM-shell coupling. Comput. Methods Appl. Mech.
Engrg., DOI: 10.1016/j.cma.2016.08.008.
Hoffman, J., Jansson, J., and Stöckli, M. (2011). Unified continuum modeling of fluid-structure
interaction. Math. Mod. Meth. Appl. Sci., 21(3):491–513.
Hughes, T. J. R., Franca, L. P., and Balestra, M. (1986). A new finite element formulation
for computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: A stable
Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations. Comput. Methods Appl. Mech. Engrg., 59(1):85–99.
Idelsohn, S. R., Oñate, E., and Del Pin, F. (2004). The particle finite element method: a
powerful tool to solve incompressible flows with free-surfaces and breaking waves. Int. J.
Num. Meth. Engrg., 61:964–989.
Idelsohn, S. R., Oñate, E., Del Pin, F., and Calvo, N. (2006). Fluid-structure interaction using
the particle finite element method. Comput. Methods Appl. Mech. Engrg., 195:2100–2123.
Jansen, K. E., Whiting, C. H., and Hulbert, G. M. (2000). A generalized-α method for integrating the filtered Navier-Stokes equations with a stabilized finite element method. Comput.
Methods Appl. Mech. Engrg., 190:305–319.
Karagiozis, K., Kamakoti, R., Cirak, F., and Pantano, C. (2011). A computational study
of supersonic disk-gap-band parachutes using large-eddy simulation coupled to a structural
membrane. J. Fluids Struc., 27:175–192.
Kloeppel, T. and Wall, W. A. (2011). A novel two-layer, coupled finite element approach for
modeling the nonlinear elastic and viscoelastic behavior of human erythrocytes. Biomech.
Model. Mechanobiol., 10(4):445–459.
Laursen, T. A. (2002). Computational Contact and Impact Mechanics: Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis. Springer.
Le, D. V., White, J., Peraire, J., Lim, K. M., and Khoo, B. C. (2009). An implicit immersed
boundary method for three-dimensional fluid–membrane interactions. J. Comput. Phys.,
228:8427–8445.
Li, Z., Hu, G.-H., Wang, Z.-L., Ma, Y.-B., and Zhou, Z.-W. (2013). Three dimensional flow
structures in a moving droplet on substrate: A dissipative particle dynamics study. Phys.
Fluids, 25:072103.
Liang, S. J., Neitzel, G. P., and Aidun, C. K. (1997). Finite element computations for unsteady
fluid and elastic membrane interaction problems. Int. J. Num. Meth. Fluids, 24:1091–1110.
Mayer, U. M., Popp, A., Gerstenberger, A., and Wall, W. A. (2010). 3D fluid-structure-contact
interaction based on a combined XFEM FSI and dual mortar contact approach. Comput.
Mech., 46:53–67.
Montefuscolo, F., Sousa, F. S., and Buscaglia, G. C. (2014). High-order ALE schemes for
incompressible capillary flows. J. Comput. Phys., 278:133–147.

34

Ohayon, R. (2004). Fluid-structure interaction problems. In Stein, E., de Borst, R., and Hughes,
T. J. R., editors, Encyclopedia of Computational Mechanics. Vol. 2: Solids and Structures.
Chapter 21. Wiley.
Peskin, C. S. (2003). The immersed boundary method. Acta Numerica, 11:479–517.
Rangarajan, R. and Gao, H. (2015). A finite element method to compute three-dimensional
equilibrium configurations of fluid membranes: Optimal parameterization, variational formulation and applications. J. Comput. Phys., 297:266–294.
Rasool, R., Corbett, C., and Sauer, R. A. (2016). A strategy to interface isogeometric analysis
with Lagrangian finite elements – Application to incompressible flow problems. Computers
& Fluids, 127:182–193.
Rasool, R., Osman, M., and Sauer, R. A. (2013). Computational modeling of liquid droplets
moving on rough surfaces. Proc. Appl. Math. Mech., 13:233–234.
Rasool, R., Sauer, R. A., and Osman, M. (2012). Internal flow analysis for slow moving small
droplets in contact with hydrophobic surfaces. Proc. Appl. Math. Mech., 12:489–490.
Saksono, P. H. and Perić, D. (2006). On finite element modelling of surface tension: Variational
formulation and applications - Part II: Dynamic problems. Comp. Mech., 38(3):251–263.
Sauer, R. A. (2014). Stabilized finite element formulations for liquid membranes and their
application to droplet contact. Int. J. Numer. Meth. Fluids, 75(7):519–545.
Sauer, R. A. (2016). A frictional sliding algorithm for liquid droplets. Comput. Mech., 58(6):937–
956.
Sauer, R. A. and De Lorenzis, L. (2015). An unbiased computational contact formulation for
3D friction. Int. J. Numer. Meth. Engrg., 101(4):251–280.
Sauer, R. A. and Duong, T. X. (2017). On the theoretical foundations of solid and liquid shells.
Math. Mech. Solids, 22(3):343–371.
Sauer, R. A., Duong, T. X., and Corbett, C. J. (2014). A computational formulation for
constrained solid and liquid membranes considering isogeometric finite elements. Comput.
Methods Appl. Mech. Engrg., 271:48–68.
Sauer, R. A., Duong, T. X., Mandadapu, K. K., and Steigmann, D. J. (2017). A stabilized
finite element formulation for liquid shells and its application to lipid bilayers. J. Comput.
Phys., 330:436–466.
Seifert, U. (1997). Configurations of fluid membranes and vesicles. Advances in Physics, 46:13–
137.
Shakib, F. (1988). Finite element analysis of the incompressible Euler and Navier-Stokes equations. PhD thesis, Stanford University, Stanford, USA.
Shelley, M. J. and Zhang, J. (2011). Flapping and bending bodies interacting with fluid flows.
Annu. Rev. Fluid Mech., 43:449–465.
Takizawa, K. and Tezduyar, T. E. (2012).
Computational Methods for Parachute
Fluid–Structure Interactions. Arch. Comput. Meth. Eng., 19:125–169.
Tezduyar, T. E. (1992). Stabilized finite element formulations for incompressible flow computations. Advances in Applied Mechanics, 28:1–44.
35

Tezduyar, T. E. and Sathe, S. (2007). Modelling of fluid–structure interactions with the
space–time finite elements: Solution techniques. Int. J. Num. Meth. Fluids, 54:855–900.
Tezduyar, T. E., Sathe, S., Stein, K., and Aureli, L. (2006). Modeling of fluid-structure interactions with the space-time techniques. In Bungartz, H. and Schäfer, M., editors, Fluid-Structure
Interaction. Lecture Notes in Computational Science and Engineering, vol. 53, pages 50–81,
Berlin, Heidelberg. Springer.
Thampi, S. P., Adhikari, R., and Govindarajan, R. (2013). Do liquid drops roll or slide on
inclined surfaces? Langmuir, 29(10):3339–3346.
van Loon, R., Anderson, P. D., van de Vosse, F. N., and Sherwin, S. J. (2007). Comparison of
various fluid–structure interaction methods for deformable bodies. Comp. Struct., 85:833–
843.
van Opstal, T. M., van Brummelen, E. H., and van Zwieten, G. J. (2015). A finiteelement/boundary-element method for three-dimensional, large-displacement fluid–structureinteraction. Comput. Methods Appl. Mech. Engrg., 284:637–663.
Walkley, M. A., Gaskell, P. H., Jimack, P. K., Kelmanson, M. A., and Summers, J. L. (2005). Finite element simulation of three-dimensional free-surface flow problems with dynamic contact
lines. Int. J. Num. Meth. Fluids, 47:1353–1359.
Wind-Willassen, Ø. and Sørensen, M. P. (2014). A finite-element method model for droplets
moving down a hydrophobic surface. Eur. Phys. J. E, 37:65.
Wriggers, P. (2006). Computational Contact Mechanics. Springer, 2nd edition.
Yang, P., Xiang, J., Fang, F., Pavlidis, D., Latham, J. P., and Pain, C. C. (2016). Modelling
of fluid–structure interaction with multiphase viscous flows using an immersed-body method.
J. Comput. Phys., 321:571–592.

36

