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 ,
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
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
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
2.2.5 Membrane weak form
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 withinand 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.
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
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 contributesforce 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
184.108.40.206 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.
220.127.116.11 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
18.104.22.168 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 .
22.214.171.124 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
(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.
126.96.36.199 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
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,
Inserting the discretized expressions for , , and into the membrane weak form (34) yields the elemental weak form contribution
with the () FE force vector