A monolithic fluid-structure interaction formulation for solid and liquid membranes including free-surface contact

by   Roger A. Sauer, et al.
RWTH Aachen University

A unified fluid-structure interaction (FSI) formulation is presented for solid, liquid and mixed membranes. Nonlinear finite elements (FE) and the generalized-alpha 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.



There are no comments yet.


page 24

page 25

page 27


An Immersed Lagrangian-Eulerian Method for Fluid-Structure Interaction

This paper introduces a sharp interface method to simulate fluid-structu...

A consistent and comprehensive computational approach for general Fluid-Structure-Contact Interaction problems

We present a consistent approach that allows to solve challenging genera...

Unified discrete multisymplectic Lagrangian formulation for hyperelastic solids and barotropic fluids

We present a geometric variational discretization of nonlinear elasticit...

A Multigrid Preconditioner for Spatially Adaptive High-order Meshless Method on Fluid-solid Interaction Problems

We present a monolithic geometric multigrid preconditioner for solving f...

On the Application of Total Traction Equilibrium in Topology Optimization of Fluid-Structure Interactions

This work investigates the different techniques of enforcing traction eq...

Coupling of IGA and Peridynamics for Air-Blast Fluid-Structure Interaction Using an Immersed Approach

We present a novel formulation based on an immersed coupling of Isogeome...

Optimization with nonstationary, nonlinear monolithic fluid-structure interaction

Within this work, we consider optimization settings for nonlinear, nonst...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 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 finite-discrete 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 immersed-boundary 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 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.

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 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 within domain is characterized by the deformation mapping


and the corresponding deformation gradient (or Jacobian)


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


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


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


It is related to the acceleration for fixed ,


according to


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 .333Following 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


which governs the fluid flow together with the boundary conditions



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


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


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


where is the Lagrange multiplier to the incompressibility constraint


which is equivalent to the condition


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 non-zero on the surface . Functions and are further assumed to possess sufficient regularity for the following integrals to be well defined. In the framework of SUPG444Streamline upwind/Petrov-Galerkin (Brooks and Hughes, 1982) and PSPG555Pressure stabilizing/Petrov-Galerkin (Hughes et al., 1986) stabilization, the weak form takes the form




is the virtual work associated with inertia,


is internal virtual work,


is the virtual work of the fluid traction on boundary ,


is the external virtual work666In the following examples we consider zero Neumann BC () and constant gravity loading with .,


is the virtual work associated with incompressibility constraint (14),


is the SUPG term,


is the PSPG term, and


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 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.

2.2.1 Membrane kinematics

The motion of a membrane surface is fully described by the mapping


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 stress-free 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 in-plane 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


which governs the membrane deformation together with the boundary conditions


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.


In order to solve PDE (25) for , the initial conditions


are needed.

2.2.3 Membrane constitution

For pure membranes, the stress tensor only has in-plane components, i.e. it has the format . Two material models are considered in this work. The first,


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


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




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


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


with the virtual work contributions


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 quasi-static 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 Kirchhoff-Love shells: and

simply need to be extended by the bending moments acting within

and on (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 moves the fluid such that


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


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


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: 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.

Figure 1: Fluid-inflated cylinder: Membrane deformation and fluid velocity due to a radial inflow at .

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


due to continuity. Since , we obtain


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


where is the radial unit vector. From this follows


such that . The equation of motion thus reduces to , which can be integrated to give the pressure field


where is the current membrane velocity, and is the pressure acting on the membrane. This pressure equilibrates the membrane stress


caused by the membrane stretch according to Eq. (29); see Appendix A. From follows


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


where , and denotes the angular velocity around . Consequently,


where . Since we can write and , we find such that and


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


The constant follows from the boundary condition , where is the surface tension of the droplet and 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


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 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 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 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

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 Basic flow variables

Within a fluid element, the fluid velocity is approximated by the interpolation


where and are the nodal shape function and nodal velocity, respectively. In short, this can also be written as


where and . The corresponding test function (or variation) is approximated in the same fashion, i.e.


The fluid pressure is approximated by the interpolation


where . Likewise,


The structure of (52) is also used to interpolate the mesh motion, i.e.


In the present work, the are not treated as unknowns. Instead they will be defined through the membrane motion. Derived flow variables

As a consequence of the above expressions, we find the approximation of the acceleration (from Eq. (7))


the velocity gradient


the pressure gradient


and the velocity divergence




and . Further, we introduce the classical B-matrix , with


in order to express the symmetric velocity gradient and its corresponding variation in Voigt notation (indicated by index ‘v’) as


i.e. arranged as . The stress tensor, arranged as , can thus be written as


with and . Here, is the usual identity tensor in . Due to the symmetry of the stress and since , the integrand of becomes


within element .
In order to represent the SUPG term, we introduce the arrays , with the blocks


and , with


The last term can also be used to rewrite the term as

(68) 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


with the () FE force vector


and the () FE pseudo force vector


They are composed of the FE forces


the FE pseudo forces


and the elemental mass, damping and pressure-force matrices


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 . Stabilization terms

In order to evaluate the residual that appears in the stabilization terms and , we note that


where , with


and , with


With this we can write


where . Thus we obtain


The stabilization parameters and appearing inside and are computed from


(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


(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.777In 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. 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




denotes the Jacobian of the mapping . Likewise, the second derivative is obtained from the formula


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


where and are arranged just like . From this follows


where . Likewise,


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


with the () FE force vector