1 Introduction
Fluidstructure interaction (FSI) problems are challenging problems due to various reasons.
They combine the computational challenges of (generally nonlinear) 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 fluidlike objects.
Examples are liquid menisci, soap films and lipid bilayers.
Lipid bilayers surround biological cells.
They are characterized by both solidlike (i.e. elastic bending) and fluidlike behavior (i.e. inplane 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 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 FEbased work on membraneFSI that is surveyed in the following. The computational approaches on interactions between fluids and membranelike 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 YoungLaplace equation), and approaches that also account for inplane equations (such as the surface NavierStokes 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 noslip 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 inplane equations.
The following references deal with solid membranes using surface formulations.
In Liang et al. (1997) the authors employ a deformable spatial domain spacetime 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 spacetime 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 thinshell formulations.
Le et al. (2009) developed an implicit immersed boundary method for the incompressible NavierStokes equations to simulate membranefluid 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 finiteelement/boundary element method for fluidstructure 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 socalled 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 liquidgaseous 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 and Perić (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 highorder 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 inplane 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 NavierStokes equations. Rangarajan and Gao (2015) introduce a splinebased finiteelement 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 inplane 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 liquidlike character, like lipid bilayers. The present work considers a monolithic coupling scheme between fluid and structure, and solves the resulting nonlinear system of equations with the NewtonRaphson 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 rotationfree shells.

Two simple analytical FSI examples are presented.

The formulation is suitable for a wide range of applications, including freesurface 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.
2 Governing equations
This section summarizes the governing equations for fluid flow, membrane deformation, membrane contact and their coupling. The symbols and 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 LagrangianEulerian (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 within domain is characterized by the deformation mapping
(1) 
and the corresponding deformation gradient (or Jacobian)
(2) 
The volume change during deformation is captured by the Jacobian determinant . The velocity of the material is given by the time derivative of for fixed , written as
(3) 
and commonly referred to as the material time derivative. It is also often denoted by the dot notation . An important object characterizing the fluid flow is the velocity gradient
(4) 
that can also be written as , where is the material time derivative of the deformation gradient.
The symmetric part of the velocity gradient is denoted by .
Likewise to Eq. (3), the material acceleration is given by
(5) 
It is related to the acceleration for fixed ,
(6) 
according to
(7) 
where is the mesh velocity (Donea and Huerta, 2003). For a purely Lagangian description , while for a purely Eulerian description .
Remark 2.1: The gradient operator appearing in Eq. (4) (and likewise in Eq. (2)), is defined here as .^{3}^{3}3Following 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.
2.1.2 Fluid equilibrium
From the balance of linear momentum within follows the equilibrium equation
(8) 
which governs the fluid flow together with the boundary conditions
(9) 
Here,
denotes the stress tensor within
,denotes the traction vector on the surface characterized by normal vector
, and denotes the density, while , and are prescribed body forces, surface velocities and surface tractions. and denote the corresponding Dirichlet and Neumann boundary regions of the fluid domain . Boundary can be split into the two parts(10) 
where is the surface of the membrane, which is considered to impose its velocity onto the fluid, and denotes the remaining Dirichlet boundary of the fluid domain. In order to solve PDE (8) for , the initial condition
(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
(12) 
where is the Lagrange multiplier to the incompressibility constraint
(13) 
which is equivalent to the condition
(14) 
A consequence of this condition is that the fluid pressure, defined as , is equal to the Lagrange multiplier . It is an additional unknown that needs to be solved for together with . In order to uniquely determine the pressure field within , 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 and , and integrated over the domain . Function is assumed to be zero on the Dirichlet boundary , but nonzero on the surface . Functions and are further assumed to possess sufficient regularity for the following integrals to be well defined. In the framework of SUPG^{4}^{4}4Streamline upwind/PetrovGalerkin (Brooks and Hughes, 1982) and PSPG^{5}^{5}5Pressure stabilizing/PetrovGalerkin (Hughes et al., 1986) stabilization, the weak form takes the form
(15) 
where
(16) 
is the virtual work associated with inertia,
(17) 
is internal virtual work,
(18) 
is the virtual work of the fluid traction on boundary ,
(19) 
is the external virtual work^{6}^{6}6In the following examples we consider zero Neumann BC () and constant gravity loading with .,
(20) 
is the virtual work associated with incompressibility constraint (14),
(21) 
is the SUPG term,
(22) 
is the PSPG term, and
(23) 
is the residual of Eq. (8). Dimensionally the residual is a force per volume. Since in theory , stabilization terms and do not affect the physical behavior of the system. In Cartesian coordinates . The scalars and 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 outofplane 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.
2.2.1 Membrane kinematics
The motion of a membrane surface is fully described by the mapping
(24) 
where , for , 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 to surface , the metric tensor components , the surface stretch , and the surface normal can be determined.
From the matrix inverse , the dual tangent vectors can be defined such that is equal to the Kronecker delta .
In order to characterize deformation, a stressfree reference configuration is introduced.
It will be considered here as the initial membrane surface, i.e. .
In the reference configuration the tangent vectors, metric tensor components, inverse components and normal vector are denoted by capital letters, i.e. , , and .
The inplane deformation of surface is fully characterized by the difference between and .
Following definitions (3) and (5), the membrane velocity and acceleration are obtained from Eq. (24).
2.2.2 Membrane equilibrium
From the balance of linear momentum within follows the equilibrium equation
(25) 
which governs the membrane deformation together with the boundary conditions
(26) 
e.g. see Sauer and Duong (2017). Here, denotes the stress tensor within , denotes the traction vector on the membrane boundary characterized by normal vector , and denotes the membrane density, while and are prescribed boundary velocities and boundary tractions. The body force is considered here to have contributions coming from the flow field, contact and external sources, i.e.
(27) 
In order to solve PDE (25) for , the initial conditions
(28) 
are needed.
2.2.3 Membrane constitution
For pure membranes, the stress tensor only has inplane components, i.e. it has the format . Two material models are considered in this work. The first,
(29) 
is suitable for solid membranes. It can be derived from the 3D incompressible NeoHookean material model (Sauer et al., 2014). The second,
(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 . During sticking contact no relative motion occurs between the membrane and a neighboring substrate surface . Mathematically this corresponds to the constraint
(31) 
where
(32) 
denotes the contact gap between the membrane point and its initial projection point on the substrate surface, , i.e. is the location where initially touched . Here, constraint (31) will be enforced by a penalty regularization. For this, the contact traction at is given by
(33) 
where is the surface normal of . 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
(34) 
with the virtual work contributions
(35) 
due to inertia, internal forces, contact forces, fluid forces and external forces acting on and .
Test function is the same as in (15).
Therefore, space needs to additionally satisfy the requirement that all integrals appearing above are well defined.
Further is assumed to be zero on .
Pure membranes are inherently unstable in the quasistatic case () and therefore need to be stabilized (Sauer et al., 2014; Sauer, 2014).
Here, no stabilization is required as the fluid forces stabilize the membrane, even when .
In the numerical examples following later, and , and consequently , are considered zero.
Remark 2.2: It is straight forward to extend weak form (34) to KirchhoffLove shells: and
simply need to be extended by the bending moments acting within
and on (Duong et al., 2017). KirchhoffLove shells are suitable for thin membranelike surface structures. Such a structure is considered in Sec. 5.3 using isogeometric finite elements.2.3 Coupling conditions
The membrane deformation moves the fluid such that
(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
(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
(38) 
which can be solved for the unknown velocity and pressure in . The membrane deformation can then be obtained from integrating . Coupling condition (37) simply leads to the cancelation of terms and in the combined weak form (38). This cancelation will carry over to the discretized weak form, as long as surface 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: Fluidinflated 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 at the inner boundary , the radial fluid velocity at location is given by
(39) 
due to continuity. Since , we obtain
(40) 
as the current position of the fluid particle initially at . The current membrane position is thus given by , where is the initial position of the membrane. In vectorial notation, the flow field can thus be characterized by the position, velocity and acceleration
(41) 
where is the radial unit vector. From this follows
(42) 
such that . The equation of motion thus reduces to , which can be integrated to give the pressure field
(43) 
where is the current membrane velocity, and is the pressure acting on the membrane. This pressure equilibrates the membrane stress
(44) 
caused by the membrane stretch according to Eq. (29); see Appendix A. From follows
(45) 
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 the influence of gravity is negligible, so that a rolling droplet remains approximately spherical. Considering the axis of rotation to be , the motion of a spinning droplet can be expressed as
(46) 
where , and denotes the angular velocity around . Consequently,
(47) 
where . Since we can write and , we find such that and
(48) 
The spin tensor, defined as , then becomes . The axial vector of , denoted by , thus is . It denotes the orientation and magnitude of the droplet’s spin, and it is equal to half of the vorticity . Solving Eq. (8) (with ) for now gives
(49) 
The constant follows from the boundary condition , where is the surface tension of the droplet and is the droplet radius. This condition enforces the YoungLapace equations, which is contained inside Eq. (25), see Sauer (2014). Applying the boundary condition, we find
(50) 
If desired, the constant velocity can be added to , such that the resulting velocity is zero at the contact point (where ).
4 Finite element formulation
The coupled fluidmembrane 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 finite elements, numbered . Some of these elements are 3D fluid elements, others are 2D surface elements or 1D line elements. Element contains nodes and occupies the domain
in the current configuration. Each fluid element has four degreesoffreedom (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
force components, while each membrane element contributes force components that need to be assembled into the global system. Those elemental forces are discussed in the following two sections.4.1.1 Fluid flow
4.1.1.1 Basic flow variables
Within a fluid element, the fluid velocity is approximated by the interpolation
(51) 
where and are the nodal shape function and nodal velocity, respectively. In short, this can also be written as
(52) 
where and . The corresponding test function (or variation) is approximated in the same fashion, i.e.
(53) 
The fluid pressure is approximated by the interpolation
(54) 
where . Likewise,
(55) 
The structure of (52) is also used to interpolate the mesh motion, i.e.
(56) 
In the present work, the 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) 
the velocity gradient
(58) 
the pressure gradient
(59) 
and the velocity divergence
(60) 
where
(61) 
and . Further, we introduce the classical Bmatrix , with
(62) 
in order to express the symmetric velocity gradient and its corresponding variation in Voigt notation (indicated by index ‘v’) as
(63) 
i.e. arranged as . The stress tensor, arranged as , can thus be written as
(64) 
with and . Here, is the usual identity tensor in . Due to the symmetry of the stress and since , the integrand of becomes
(65) 
within element .
In order to represent the SUPG term, we introduce the arrays , with the blocks
(66) 
and , with
(67) 
The last term can also be used to rewrite the term as
(68) 
4.1.1.3 Weak form contribution of a fluid element
Given the above expressions, the contributions from element to the fluid weak form (15) can be written as
(69) 
with the () FE force vector
(70) 
and the () FE pseudo force vector
(71) 
They are composed of the FE forces
(72) 
the FE pseudo forces
(73) 
and the elemental mass, damping and pressureforce matrices
(74) 
The tangent matrices of and , needed for linearization, can be found in Appendix B.1.
Remark 4.1: One may simply change the sign of both and in order to highlight the symmetry between the second part of and .
4.1.1.4 Stabilization terms
In order to evaluate the residual that appears in the stabilization terms and , we note that
(75) 
where , with
(76) 
and , with
(77) 
With this we can write
(78) 
where . Thus we obtain
(79) 
The stabilization parameters and appearing inside and are computed from
(80) 
(Shakib, 1988; Tezduyar, 1992; Rasool et al., 2016), where is the time step size, is the “element length” in the local flow direction taken from
(81) 
(Tezduyar, 1992) and depends on the polynomical order of the shape functions. I.e. for L1 (linear Lagrange) and L2 (quadratic Lagrange) elements we have and , respectively.^{7}^{7}7In Eqs. (80) and (81), is taken from the previous time step in order to avoid the linearization of and . According to this, parameters and are local parameters that change from quadrature point to quadrature point.
4.1.1.5 Transformation of derivatives
In the above expressions denotes the gradient w.r.t. the current configuration , which is discretized by , where are the nodal positions of the FE mesh. Since it is convenient to define the shape functions on a master element in space, where is easily obtained, needs to be determined from
(82) 
where
(83) 
denotes the Jacobian of the mapping . Likewise, the second derivative is obtained from the formula
(84) 
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
(85) 
where and are arranged just like . From this follows
(86) 
where . Likewise,
(87) 
follows from Eq. (53). Given and , the metric tensor components and can be determined and the stress can be evaluated as discussed in Sec. 2.2.
Inserting the discretized expressions for , , and into the membrane weak form (34) yields the elemental weak form contribution
(88) 
with the () FE force vector
Comments
There are no comments yet.