# The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction

Andreas Zilian, Antoine Legay

► To cite this version:

Andreas Zilian, Antoine Legay. The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction. International Journal for Numerical Methods in Biomedical Engineering, 2008, 75 (3), 10.1002/nme.2258. hal-01371129

HAL Id: hal-01371129

https://hal.science/hal-01371129

Submitted on 24 Sep 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction

A. Zilian¹,*,† and A. Legay²

¹Institut für Statik, Technische Universität Braunschweig, Beethovenstraße 51, 38106 Braunschweig, Germany
²Structural Mechanics and Coupled Systems Laboratory, Conservatoire National des Arts et Métiers, 2 Rue Conté, 75003 Paris, France

## SUMMARY

The paper introduces a weighted residual-based approach for the numerical investigation of the interaction of fluid flow and thin flexible structures. The presented method enables one to treat strongly coupled systems involving large structural motion and deformation of multiple-flow-immersed solid objects.

The fluid flow is described by the incompressible Navier-Stokes equations. The current configuration of the thin structure of linear elastic material with non-linear kinematics is mapped to the flow using the zero iso-contour of an updated level set function. The formulation of fluid, structure and coupling conditions uniformly uses velocities as unknowns. The integration of the weak form is performed on a space-time finite element discretization of the domain. Interfacial constraints of the multi-field problem are ensured by distributed Lagrange multipliers. The proposed formulation and discretization techniques lead to a monolithic algebraic system, well suited for strongly coupled fluid-structure systems. Embedding a thin structure into a flow results in non-smooth fields for the fluid. Based on the concept of the extended finite element method, the space-time approximations of fluid pressure and velocity are properly enriched to capture weakly and strongly discontinuous solutions. This leads to the present enriched space-time (EST) method. Numerical examples of fluid-structure interaction show the eligibility of the developed numerical approach in order to describe the behavior of such coupled systems. The test cases demonstrate the application of the proposed technique to problems where mesh moving strategies often fail.

**KEY WORDS:** fluid-structure interaction; enriched space-time finite elements; EST method; level set method

<footnote>*Correspondence to: A. Zilian, Institut für Statik, Technische Universität Braunschweig, Beethovenstr. 51, 38106 Braunschweig, Germany.</footnote>

<footnote>†E-mail: a.zilian@tu-bs.de</footnote>
# 1. INTRODUCTION

The evaluation of fluid-structure interaction effects and the investigation of governing physical phenomena associated with coupled systems are always challenges in problems arising in the field of aero- and hydro-elasticity, life-sciences or bio-engineering. An important field of boundary-coupled interaction problems involves single or multiple thin and very flexible structures, e.g. membranes and plates. Thin structures situated in a surrounding fluid flow or encapsulating a fluid may influence drastically the pattern of the flow and the characteristics of the physical quantities describing the fluid. The flow field may diversify depending on the presence and characteristics of boundary conditions, shear layers or vortex shedding effects. Embedding such a Lagrangian structure of small thickness into an Eulerian fluid field results in discontinuous state variables of the flow. For a velocity-pressure formulation of the Navier-Stokes equations and no-slip interfacial conditions, this produces a strongly discontinuous pressure solution and a weakly (gradient) discontinuous velocity solution along the moving fluid-structure interface. In the case of interfacial slip conditions, the tangential component of the velocity may exhibit a strong discontinuity. Especially for engineering tasks in fluid-structure interaction, it seems to be reasonable to apply simplified models as a first approximation to complex physical processes. As an alternative to time-consuming resolution of all involved length and time scales at moving boundaries, it can be beneficial to treat existent non-smooth solution characteristics, e.g. discontinuities, by taking them explicitly into account in an advanced numerical strategy.

The arbitrary Lagrange–Euler (ALE) method [1–5] is classically applied to treat flow with immersed objects by using the advantages of an Eulerian description of the fluid domain and a Lagrangian description of the solid domain. Combined with a finite element discretization of both domains, a conservative coupling can be achieved by the compatibility of the two meshes (geometry) and approximation (physics) at the common fluid-structure interface. This condition may introduce severe complications, e.g. fluid mesh distortions (see Figure 1), if mesh-moving strategies are applied. Moreover, an ALE-based approach tends to be limited when multiple immersed structures are subjected to rotation or large relative displacements.

In order to avoid these restrictions and to consider physical effects at complex interfaces, several numerical methodologies have been developed. A moving and evolving interface of arbitrary shape can be represented by an implicit method using a fixed background grid. This has the advantage of avoiding the need for mesh-moving or re-meshing procedures. Figure 2 shows explicit and implicit methods of incorporating a thin structure into the fluid domain. Implicit descriptions are used in the immersed boundary method (IBM) by Peskin and McQueen [6] for investigation of fibers in viscous flow, in the fictitious domain method (FDM) developed by Bertrand et al. [7] and Glowinski et al. [8] for moving rigid bodies in a fluid and in the immersed finite element method (IFEM) introduced by Zhang et al. [9] for the simulation of elastic solids in fluids. Within these methods, the conservative and computationally stable coupling of solid and fluid can be realized by higher-order schemes, increasing the numerical effort.

The computational approach to fluid-structure interaction proposed by this contribution focuses on flow-immersed moving thin flexible structures and is based on an implicit interface description. The numerical scheme employs the extended finite element method (XFEM) presented by Belytschko et al. [10] as well as a space-time discretization.

The XFEM combines a level set representation of interfaces [11] and the partition of unity method (PUM) [12]. It enables the representation of special and *a priori* known solution properties within an enriched finite element approximation. It has been applied successfully by Chessa and
<figure><img src="image_2.png" /><figcaption>Figure 1. Distorted ALE meshes due to structural motion.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Interface description: explicit (left) and implicit (right) on a fixed background grid.</figcaption></figure>

Belytschko [13] to flow problems consisting of two immiscible fluids besides many applications in the field of structural mechanics with material interfaces [14] and crack propagation [15]. It has been applied by Legay *et al.* to a compressible fluid bounded by a flexible structure [16] and has been extended to thin immersed structures with prescribed motion in an incompressible and viscous flow [17] using a fractional time stepping scheme. Further, the extended finite element concept has been applied by Kölke *et al.* in conjunction with an enriched space–time (EST) finite element discretization [18] of two-liquid flow problems including surface tension, where discontinuous pressure solutions at evolving and level set-captured interfaces are represented by an extended approximation in space–time. This approach has been discussed by Chessa and Belytschko [19, 20] to arbitrary discontinuities based on space–time finite elements and was realized for spatially one-dimensional shock-dominated problems.

The space–time finite element method for solving the incompressible Navier–Stokes equa-
tions on moving meshes has been introduced by Tezduyar *et al.* [21]. For applications in elas-
todynamics and analyses of wave propagation problems, Hughes and Hulbert [22] introduced
time-discontinuous, stabilized space–time (SST) finite elements. Tezduyar *et al.* developed the
deforming-spatial-domain/stabilized space–time (DSD/SST) finite element formulation [23] for
fluid–structure interaction problems. The space–time finite element discretization of a velocity-
based monolithic model for strongly coupled fluid–structure interaction problems has been shown
to be eligible for a number of different applications by Hübner *et al.* [24, 25] including free-surface
and two-fluid flows [26–28].

The major conceptual advantage of the space–time finite element method is the straightforward applicability of the enriched finite element technology (XFEM) to describe the evolution of non-smooth solutions within the space–time domain more accurately. Many applications of XFEM focus on stationary problems of structural analysis or combine spatial enrichment with semi-discrete time stepping methods for dynamical problems. Therefore, for instationary problems the whole potential of XFEM is achieved by enriching a finite element approximation in both space and time. The combination of level set method, local partition of unity enrichment and a space–time
finite element discretization thus offers an elegant way of describing evolving non-smooth solution features [29, 30] arising from embedding thin structures in a flow field.

The governing equations of fluid- and elastodynamics and interface conditions are introduced in Section 2. Section 3 gives the weighted residual form of the coupled equations. In Section 4 the EST finite element method is proposed. Algorithmic details are given in Section 5. Section 6 introduces numerical test cases of strongly coupled fluid-structure interaction.

## 2. STRONG FORM GOVERNING EQUATIONS

### 2.1. *Coupled fluid-structure system*

The coupled system consists of a thin flexible structure immersed in a fluid as shown in Figure 3. In comparison with the spatial dimensions of the fluid body, the thickness of the solid is assumed to be very small. In this paper the case of an embedded structure that divides the flow domain into two separate parts is considered. A free end of the structure would introduce flow singularities at the tip that requires a treatment by eligible enrichment functions. While this is generally possible, this case is not considered for simplicity and taken out of the focus of this paper.

The behavior of the fluid body within the time interval $I = [t_a, t_b]$ and occupying the space-time domain $Q = \Omega \times I$ is described in the current configuration $Q_t$ using an Eulerian description. A Lagrangian description and the reference configuration $L_0 = \Lambda_0 \times I$ are used for the structure. The fluid domain is bounded by the outer boundary $P_t = \Gamma_t \times I$ and the evolving interface $R_t = \Sigma_t \times I$.

Although the discretization of the fluid-structure interaction problem with the EST method is demonstrated for spatially two-dimensional problems, the approach is applicable to three-dimensional space problems as well. While within the derivation of the governing equations only one immersed structure is considered, the method is extendable to multiple immersed structures. For simplicity, the methodology and applications are restricted to problems where no contact between structures is expected.

### 2.2. *Fluid*

In this paper the flow is modeled as incompressible, viscous and isothermal. The incompressible Navier-Stokes equations are

$$ \rho \mathbf{v}_{,t} + \rho \mathbf{v} \cdot \nabla \mathbf{v} - \nabla \cdot \mathbf{T} - \rho \mathbf{b} = \mathbf{0} \quad \text{in } Q_t \quad (1) $$

$$ \nabla \cdot \mathbf{v} = 0 \quad \text{in } Q_t \quad (2) $$

where $\mathbf{v}$ is the fluid velocity, $\mathbf{b}$ is the external volume force and $\rho$ is the fluid density. The Cauchy stress tensor $\mathbf{T}$ of the Newtonian fluid is given by the following relation:

$$ \mathbf{T} = -p\mathbf{I} + 2\mu\mathbf{D} \quad \text{in } Q_t \quad (3) $$

where $p$ is the hydrostatic pressure, $\mu$ is the kinematic viscosity and $\mathbf{D}$ is the strain rate tensor defined by

$$ \mathbf{D} = \frac{1}{2}(\nabla \mathbf{v} + (\nabla \mathbf{v})^\mathrm{T}) \quad \text{in } Q_t \quad (4) $$
<figure><img src="image_4.png" /><figcaption>Figure 3. Configuration of the coupled problem of flow and immersed thin-walled structure at $t \in I$.</figcaption></figure>

Dirichlet boundary conditions are defined on $P_t^g$ in terms of imposed velocity $\bar{v}$:

$$ \mathbf{v}-\bar{\mathbf{v}}=\mathbf{0} \quad \text{on } P_t^g \qquad (5) $$

where as Neumann boundary conditions are defined on $P_t^h$ in terms of the imposed boundary traction $\bar{\mathbf{t}}$:

$$ \mathbf{T}\mathbf{n}-\bar{\mathbf{t}}=\mathbf{0} \quad \text{on } P_t^h \qquad (6) $$

The outward normal vector $\mathbf{n}$ to the fluid boundary is defined in the current configuration. Furthermore, Equation (1) requires the definition of an initial velocity field $\tilde{\mathbf{v}}$ within $\Omega(t=t_a)$:

$$ \mathbf{v}(t=t_a)-\tilde{\mathbf{v}}=\mathbf{0} \quad \text{on } \Omega(t=t_a) \qquad (7) $$

with respect to the incompressibility constraint formulated in Equation (2).

## 2.3. Thin-walled structure

The structure is assumed to have a constant and very small thickness $h$. The resulting simplified plate-like structure is able to show bending behavior and membrane properties, where the cross sections remain plane and rotate about a defined neutral local axis, but do not remain orthogonal to the deformed axis (Timoshenko model). Out-of-plane strains are supposed to be zero referring to the $(X, Y)$ plane. Additionally, out-of-plane stresses are set to zero to the $(X, Z)$ plane. Figure 4 shows the structural model in its reference configuration $\Lambda_0$ and its current configuration $\Lambda_t$, as well as the main local axis to which the sectional kinematics, strain and stress measures are related. The definition of motion and kinematics follows Felippa [31].

The one-dimensional structure with local coordinates $(X, Y)$ is described in a total Lagrangian setting by the balance of momentum:

$$ \rho \mathbf{G} \dot{\mathbf{w}} - (\mathbf{F}_s)_X - (\mathbf{H}_s) - \rho \mathbf{G} \mathbf{b} = \mathbf{0} \quad \text{in } L_0 \qquad (8) $$

where

$$ \mathbf{w}=[v_X \ v_Y \ \omega]^T \qquad (9) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. Two-dimensional model of a thin-walled plane structure surrounded by a fluid.</figcaption></figure>

is a vector consisting of local velocities along (X, Y) and angular velocity. External loads are collected in vector **b**=[b<sub>X</sub> b<sub>Y</sub> b<sub>ω</sub>]<sup>T</sup> and ρ is the density of the solid. The resulting internal forces are collected in the vector:

$$ \mathbf{s} = [N \ Q \ M]^T \qquad (10) $$

where N, Q and M are the normal force, shear force and moment, respectively. The geometry of the section is described by the diagonal matrix:

$$ \mathbf{G} = \begin{bmatrix} h & 0 & 0 \\ 0 & h & 0 \\ 0 & 0 & h^3/12 \end{bmatrix} \qquad (11) $$

The deformation gradient **F** is given by

$$ \mathbf{F} = \begin{bmatrix} 1+u_{X,X} & -\sin\theta & \theta_{,X}\cos\theta \\ u_{Y,X} & \cos\theta & \theta_{,X}\sin\theta \\ 0 & 0 & \cos\theta+u_{X,X}\cos\theta+u_{Y,X}\sin\theta \end{bmatrix} \qquad (12) $$

and matrix **H** is defined as

$$ \mathbf{H} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & \cos\theta+u_{X,X}\cos\theta+u_{Y,X}\sin\theta & \theta_{,X}(\sin\theta+u_{X,X}\sin\theta-u_{Y,X}\cos\theta) \end{bmatrix} \qquad (13) $$
involving derivatives of the local displacements $u_X$ and $u_Y$ as well as the rotation $\theta$ of the cross section. Since the elasticity matrix

$$ \mathbf{C} = \mathbf{G} \begin{bmatrix} \frac{E}{1-v^2} & 0 & 0 \\ 0 & \frac{E}{2(1+v)} & 0 \\ 0 & 0 & E \end{bmatrix} \qquad (14) $$

with Young's modulus $E$ and the Poisson ratio $\nu$ is time invariant for linear elastic material, the constitutive relation can be given in rate formulation:

$$ \mathbf{C}^{-1} \dot{\mathbf{s}} - \dot{\mathbf{e}} = \mathbf{0} \quad \text{in } L_0 \qquad (15) $$

The vector of local non-linear strain measures $\mathbf{e}=[e_1 \ e_2 \ e_3]^T$ for a thin structure is defined by

$$ e_1 = u_{X,X} + \frac{1}{2}(u_{X,X}^2 + u_{Y,X}^2) \qquad (16a) $$

$$ e_2 = -(1+u_{X,X}) \sin\theta + u_{Y,X} \cos\theta \qquad (16b) $$

$$ e_3 = \theta,_{X}(\cos\theta + u_{X,X}\cos\theta + u_{Y,X}\sin\theta) \qquad (16c) $$

where $e_1$ is the axial strain, $e_2$ is the shear strain and $e_3$ is the curvature. The rate of the introduced strain measures is summarized in the components of the local strain rate vector $\dot{\mathbf{e}}$:

$$ \dot{e}_1 = v_{X,X} + u_{X,X} v_{X,X} + u_{Y,X} v_{Y,X} \qquad (17a) $$

$$ \dot{e}_2 = -v_{X,X} \sin\theta + v_{Y,X} \cos\theta - \omega(\cos\theta + u_{X,X} \cos\theta + u_{Y,X} \sin\theta) \qquad (17b) $$

$$ \begin{aligned} \dot{e}_3 = & v_{X,X} \theta,_{X} \cos\theta + v_{Y,X} \theta,_{X} \sin\theta \\ & + \omega,_{X}(\cos\theta + u_{X,X} \cos\theta + u_{Y,X} \sin\theta) - \omega \theta,_{X}(\sin\theta + u_{X,X} \sin\theta - u_{Y,X} \cos\theta) \end{aligned} \qquad (17c) $$

using the local structural velocities $v_X$, $v_Y$ and the rotational velocity $\omega$. Structural displacements are collected in the vector

$$ \mathbf{d} = [u_X \ u_Y \ \theta]^T \qquad (18) $$

and are computed by time integration as follows:

$$ \mathbf{d}(t) = \int_{t_a}^{t} \mathbf{w}(\tau) \, d\tau \qquad (19) $$

The complete structural boundary is divided into the outer sectional boundary $\tilde{P}_0 = \tilde{\Pi}_0 \times I$ and the boundary $\bar{P}_0 = \bar{\Pi}_0 \times I$ that coincides with the interface to the surrounding fluid domain. On the Neumann part $\tilde{P}_0^h$ of the sectional boundary, external forces

$$ \mathbf{s} - \bar{\mathbf{s}} = \mathbf{0} \quad \text{on} \quad \tilde{P}_0^h \qquad (20) $$
can be imposed. Prescribed generalized boundary tractions $\bar{\mathbf{t}}=[\bar{t}_X \bar{t}_Y \bar{t}_\omega]^T$ are defined by

$$ \mathbf{G}^{-1}\mathbf{s}-\bar{\mathbf{t}}=\mathbf{0} \quad \text{on } \bar{P}_0^h \qquad (21) $$

Equation (8) requires the definition of an initial velocity field $\mathbf{w}$ within $\Lambda(t=t_a)$

$$ \mathbf{w}(t=t_a)-\tilde{\mathbf{w}}=\mathbf{0} \quad \text{on } \Lambda(t=t_a) \qquad (22) $$

The highly non-linear strain measures result from considering non-linear kinematics and allow large deformations and rotations of the thin-walled structure model.

## 2.4. Coupling conditions

Flow domain and structural domain are coupled along the common space-time boundary $R$ representing the fluid-structure interface. At the interface no-slip conditions are applied, requesting continuity of fluid and structural velocities

$$ \mathbf{v}^{\mathrm{F}}-\mathbf{v}^{\mathrm{S}}=\mathbf{0} \quad \text{on } R_t \qquad (23) $$

and therefore ensure geometrical conservation. The vector $\mathbf{v}^S$ contains the velocities $[v_1^S v_2^S]^T$ of the thin structure in the global coordinate system, obtained by a transformation of the local velocity values $[v_X v_Y]^T$. Further, fluid and solid interfacial tractions have to be of the same magnitude and in opposite direction

$$ \mathbf{t}^{\mathrm{F}}+\frac{\mathrm{d}\Lambda_0}{\mathrm{d}\Lambda_t}\mathbf{t}^{\mathrm{S}}=\mathbf{0} \quad \text{on } R_t \qquad (24) $$

in order to fulfill momentum balance at the interface. The referential solid traction $\mathbf{t}^S$ is projected to the current frame. The coupling conditions in Equations (23) and (24) represent no-slip conditions at the interface. Slip interface conditions or frictional effects can be considered by formulating the interfacial balance equations in normal and tangential directions.

## 3. WEAK FORM AND SPACE-TIME FINITE ELEMENT DISCRETIZATION

The weighted residual method is applied to the strong form equations (Sections 2.2–2.4). The resulting weak form of the whole coupled system and the space-time domain are then discretized using the space-time finite element method [32] and a time-discontinuous Galerkin method for integration in time. The basic idea of a space-time discretization is to include the temporal axis in the finite element discretization. For numerical efficiency, the space-time domain $Q$ is divided into a sequence of $N$ time slabs $Q^n=\Omega^n \times [t_n, t_{n+1}]$, as shown in Figure 5, which are solved successively.

At time instant $t_n$, the energy of the discretized system at the end of the previous time slab $t_n^-$ has to be equal to the energy at the beginning of the next time step $t_n^+$. For time-discontinuous approximations of field unknowns, this leads to additional jump terms in the weak form. Moreover, spatial discretizations from $t_n^-$ and $t_n^+$ do not need to be conforming. For first-order ordinary differential equations, the resulting time integration scheme is A-stable and third-order accurate for linear temporal interpolation.
<figure><img src="image_6.png" /><figcaption>Figure 5. Discretization of the continuous space-time domain using space-time finite elements.</figcaption></figure>

### 3.1. Weak form of fluid

The weighted residual formulation of the strong forms (1)–(6) for the incompressible viscous fluid in a space-time slab $Q_t^n$ using the Galerkin method reads as

$$ \int_{Q_t^n} \delta\mathbf{v} \cdot \rho(\mathbf{v}_{,t} + \mathbf{v} \cdot \nabla \mathbf{v}) dQ_t + \int_{Q_t^n} \mathbf{D}(\delta\mathbf{v}) : 2\mu\mathbf{D} dQ_t - \int_{Q_t^n} \nabla \cdot (\delta\mathbf{v}) p dQ_t - \int_{Q_t^n} \delta\mathbf{v} \cdot \rho\mathbf{b} dQ_t \quad (25a) $$

$$ + \int_{Q_t^n} \delta p \nabla \cdot \mathbf{v} dQ_t \quad (25b) $$

$$ + \int_{\Omega_t^n} \delta \mathbf{v}(t_n^+) \cdot \rho(\mathbf{v}(t_n^+) - \mathbf{v}(t_n^-)) d\Omega_t \quad (25c) $$

$$ + \sum_e \int_{Q_t^{n,e}} (\rho \delta \mathbf{v}_{,t} + \rho \mathbf{v} \cdot \nabla (\delta \mathbf{v}) - \nabla \cdot \mathbf{T}(\delta \mathbf{v}, \delta p) \cdot \tau \frac{1}{\rho} (\rho \mathbf{v}_{,t} + \rho \mathbf{v} \cdot \nabla \mathbf{v} - \nabla \cdot \mathbf{T} - \rho \mathbf{b})) dQ_t \quad (25d) $$

$$ - \int_{P_t^{n,h}} \delta \mathbf{v} \cdot \bar{\mathbf{t}} dP_t = 0 \quad \forall \delta \mathbf{v}, \delta p \quad (25e) $$

Line (25a) represents the weak form of conservation of momentum (1), fulfilling the constitutive relation (3) for the Newtonian fluid and the kinematics (4) exactly. The incompressibility constraint (2) is weighted with the pressure in line (25b). Line (25c) ensures the consistent transfer of kinetic energy from the previous time slab end at $t_n^-$ to the current time slab at $t_n^+$. The weighted residual form is stabilized by a Galerkin/least-squares term [33] of the momentum balance in line (25d). The Galerkin/least-squares stabilization suppresses numerical oscillations in solutions to hyperbolic differential equations by the introduction of additional numerical diffusion and allows the application of equal-order approximations of velocities and pressure for the incompressible flow field. The stabilization parameter $\tau$ is determined for each element $e$ using the definition given by Tezduyar et al. [21]. Neumann boundary conditions can be imposed in a weak sense by (25e), while (5) is treated a priori as an essential boundary condition.
## 3.2. Weak form of thin-walled structure

The weighted residual formulation of the geometrically non-linear, linear elastic thin-walled structure on slab $L_0^n = \Lambda_0 \times [t_n, t_{n+1}]$ using the Galerkin method

$$ \int_{L_0^n} \delta \mathbf{w} \cdot \rho \mathbf{G} \dot{\mathbf{w}} dX + \int_{L_0^n} \dot{\mathbf{e}}(\delta \mathbf{w}, \mathbf{d}) \cdot \mathbf{s} dX - \int_{L_0^n} \delta \mathbf{w} \cdot \rho \mathbf{G} \mathbf{b} dX \quad (26a) $$

$$ + \sum_e \int_{L_0^n} \delta \mathbf{s} \cdot (\mathbf{C}^{-1} \dot{\mathbf{s}} - \dot{\mathbf{e}}(\mathbf{w}, \mathbf{d})) dX \quad (26b) $$

$$ + \int_{L_0^n} \delta \mathbf{w}(t_n^+) \cdot \rho \mathbf{G}(\mathbf{w}(t_n^+) - \mathbf{w}(t_n^-)) dX \quad (26c) $$

$$ + \int_{L_0^n} \delta \mathbf{s}(t_n^+) \cdot \mathbf{C}^{-1} (\mathbf{s}(t_n^+) - \mathbf{s}(t_n^-)) dX \quad (26d) $$

$$ -[\delta \mathbf{w} \cdot \bar{\mathbf{s}}]_{\bar{P}_0^{n,h}} - \int_{\bar{P}_0^{n,h}} \delta \mathbf{w} \cdot \bar{\mathbf{t}} dX = 0 \quad \forall \delta \mathbf{w}, \delta \mathbf{s} \quad (26e) $$

applies the local quantities for kinematics and forces introduced in Section 2.3 and describes the elastodynamics of the thin solid in terms of local velocities ($v_X, v_Y$) and rotational velocity $\omega$. Line (26a) gives the integral form of momentum balance after partial integration utilizing the nonlinear strain rate $\dot{\mathbf{e}}$. Local displacements ($u_X, u_Y$) and the angle of rotation $\theta$ are obtained using Equation (19). The linear elastic material is considered in integral form by line (26b) on element level leading to a mixed formulation. The local forces in $\mathbf{s}$ are condensable on the element resulting in a mixed hybrid method, see Section 3.4.2. Kinetic energy is transferred weakly between time slabs with the jump term in line (26c), and internal mechanical energy is conserved in integral form between time slabs by (26d). Sectional natural boundary conditions can be imposed weakly with the first term in line (26e), while boundary tractions (from the surrounding fluid) along the axis are considered by the second term.

## 3.3. Coupling conditions

The coupling model expressed by the interfacial conditions can be introduced using various methods. Penalty methods add the constraints to the set of discrete equations and weight them by a factor $\beta$, providing a certain amount of accuracy while avoiding additional degrees of freedom. However, the heuristic factor considerably increases the condition number of the algebraic system. Alternatively, Nitsche's method or an approach based on Lagrange multipliers could be used.

The weighted residual formulation of the coupling conditions (23) and (24) utilizes interfacial tractions as additional Lagrange multipliers in order to fulfill geometrical continuity and conservation of momentum between fluid and structure. The structural velocities $\mathbf{v}^S$ are used to formulate Dirichlet boundary conditions for the fluid, while interface tractions $\mathbf{t}^F$ of the fluid act as a surface load in the form of a Neumann condition onto the deformed local axis of the thin solid

$$ + \int_{R_t^n} \delta \mathbf{t}^F \cdot (\mathbf{v}^F - \mathbf{v}^S) dR - \int_{R_t^n} \delta \mathbf{v}^F \cdot \mathbf{t}^F dR + \int_{R_t^n} \delta \mathbf{v}^S \cdot \mathbf{t}^F dR \quad (27) $$
<figure><img src="image_7.png" /><figcaption>Figure 6. Bi-linear isoparametric approximation of geometry and physics for fluid and structure.</figcaption></figure>

where $t^F$ are additional unknowns to the coupled system. The integral form of the coupling conditions is defined in the current configuration $R_t^n$ of the fluid-structure interface within the $n$th time slab. The construction of an eligible approximation for the Lagrange multiplier is discussed in Section 4.3.

## 3.4. Space-time finite element types for flow and thin structure

The above deduced modeling of fluid and immersed thin solid body for two space dimensions leads to a three-dimensional model of the flow and a two-dimensional structural model in space-time.

The approximation of geometry and physics of finite elements for fluid and structure, shown in Figure 6, is chosen linear and spans a bi-linear space-time approximation space. It is formed by standard ansatz functions $N_k$ on the support of discrete nodes $k \in \mathcal{N}$ of the set of nodes $\mathcal{N}$ in space-time slab $Q^n$. Independent degrees of freedom are used at discrete time levels $t_n$ of adjacent time slabs.

### 3.4.1. Fluid element.
The six-node space-time fluid element approximates the spatial domain by a 3-node triangle that is extruded in time direction to obtain the prismatic shape of the space-time element. The same bi-linear shape functions are used for standard approximation of both velocity and pressure. Each node $k$, with $k=1...6$, has two velocity unknowns ($v_1^k, v_2^k$) and one pressure unknown $p^k$.

### 3.4.2. Structural element.
The four-node space-time thin-solid element is a typical rectangular finite element, where the space dimension is discretized by linear polynomials as well as the time dimension. Since the structural finite element formulation describes the solid in the reference configuration, deformations of the space-time element do not have to be considered. The same shape functions are used for both velocity and angular velocity since a Timoshenko beam assumption is made. Each node has two velocity unknowns ($v_1^k, v_2^k$) and one angular velocity $\omega^k$. The normal force is approximated by piecewise linear(space)/linear(time) shape functions and is condensed
statically on element level. To prevent shear locking, constant(space)/linear(time) functions are used for the condensed shear force and the bending moment within the element of the mixed-hybrid formulation.

# 4. EST FINITE ELEMENT APPROXIMATION OF THE FLOW

## 4.1. Localization of moving thin structures in the space-time domain by level sets

Taking advantage of an *a priori* knowledge on the behavior and characteristics of the solution to the problem under consideration facilitates the proper enrichment of a space-time approximation with finite elements. If the special solution is of local character, implicit methods of describing interfaces and boundaries, e.g. one or multiple level set functions, can be applied to localize the necessary enrichment. For evolving local characteristics, the level set method introduced by Osher and Sethian has been shown to be an eligible approach for various physical problems.

Figure 7 shows an instant of a surface of discontinuity within the space-time domain and its representation by the zero level set of a space-time defined level set function $\phi(\mathbf{x}, t)$. For application to a fluid-immersed thin solid, the fluid-structure interface is identical to the space-time surface of the evolving discontinuous Eulerian flow field.

## 4.2. Local enrichment of the approximation in the fluid domain

The space-time formulation of the model equations for fluid and structure allows one to track strongly and weakly discontinuous solutions at evolving interfaces (represented by the zero level set) by using the partition of unity concept. Applying the XFEM [10] and their adaption to space-time finite elements [18, 34], the velocity approximation of the fluid domain is locally enriched by an enrichment function $\psi_j(\mathbf{x}, t)$ that realizes a jump-type enrichment:

$$ \mathbf{v}^{\mathrm{F}}(\mathbf{x}, t) = \sum_{k \in \mathcal{N}} N_k^{\mathrm{F}}(\mathbf{x}, t) \mathbf{v}^k \sum_{j \in \mathcal{M}} N_j^{\mathrm{F}}(\mathbf{x}, t) \psi_j(\mathbf{x}, t) \mathbf{w}^j \quad (28) $$

and introduces additional unknowns $\mathbf{w}^j$ at the set of enriched nodes $\mathcal{M}$. The continuity of the formally fully decoupled velocity field at the fluid-structure interface required by Equation (23) is ensured in integral form by the stabilized distributed Lagrange multiplier approach as described in Section 4.3. The combination of an enrichment for strong velocity discontinuities and Lagrange

<figure><img src="image_8.png" /><figcaption>Figure 7. Evolving non-smooth solutions in the space-time domain.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8. Enriched space-time ansatz functions: example of zero level set and enrichment functions.</figcaption></figure>

multipliers for continuity avoids here the introduction of blending elements [35]. Similar to the velocity field, the approximation of the fluid pressure is locally extended by the enrichment function $\psi_j$, enabling strong discontinuities:

$$p^{\mathrm{F}}(\mathbf{x}, t) = \sum_{k \in \mathcal{N}} N_k^{\mathrm{F}}(\mathbf{x}, t) p^k + \sum_{j \in \mathcal{M}} N_j^{\mathrm{F}}(\mathbf{x}, t) \psi_j(\mathbf{x}, t) q^j \quad (29)$$

and introducing additional unknowns $q^j$ at the set of enriched nodes $\mathcal{M}$. The enrichment function $\psi_j$ incorporates the placement of the enrichment within the space-time finite element by using the level set function $\phi(\mathbf{x}, t)$:

$$\psi_j(\mathbf{x}, t) = \frac{1}{2}(1 - \mathrm{sign} \, \phi(\mathbf{x}, t) \cdot \mathrm{sign} \, \phi(\mathbf{x}^j, t^j)) \quad (30)$$

where $(\mathbf{x}^j, t^j)$ are the coordinates of node $j$ in the space-time domain. This enrichment function ensures the Kronecker-$\delta$ property of the approximation.

Figure 8 exemplarily gives the first steps of constructing the enriching ansatz functions for the prismatic space-time finite element. Here, the nodal values $\hat{\phi} = (-0.3, 0.4, 0.8, -0.5, 0.2, 0.6)$ of the interpolated level set function describe (by the associated iso-contour) the location of both: the fluid-structure interface and the discontinuous flow solution within the element. Here, two different types of nodal enrichment functions $\psi_j$ can be obtained from Equation (30). A representation of the used bi-linear space-time ansatz functions is shown in the upper row of Figure 9. The intensity identifies the specific value of the ansatz function and ranges between zero and unity. The bi-linear basis is applied for constructing additional discontinuous shape functions with the help of the nodal enrichment function $\psi_j$. The six additional enriching shape functions are shown in the bottom row of Figure 9.

### 4.3. Distributed perturbed Lagrange multiplier technique

Interfacial constraints (no-slip interface conditions in this paper) are fulfilled in an integral sense (Lagrange multiplier approach) by the weak formulation of the coupled problem as described in Section 3.3. The construction of a proper Lagrange multiplier space along propagating immersed interfaces with respect to the LBB condition and optimal convergence is not obvious. The imposition of constraints for enriched approximations using XFEM has been discussed in [36, 37] and for general non-mesh-matching interfaces in [38].
<figure><img src="image_10.png" /><figcaption>Figure 9. Enriched space-time ansatz functions: bi-linear basis (top) and bi-linear enriched (bottom).</figcaption></figure>

For level set-captured fluid-structure interfaces $R_t$ that are moving and deforming in the space-time domain, an alternative implicit Lagrange multiplier method has been developed. The space-time-defined interface traction $\mathbf{t}^F(\mathbf{x}, t)$ in Equation (27) is represented with the help of a higher-dimensional function $\mathbf{z}^F(\mathbf{x}, t)$ and localized by the zero level set $\phi(\mathbf{x}, t)=0$, so that

$$ \mathbf{t}^F = \mathbf{z}^F \delta_D(\phi) \qquad (31) $$

where $\delta_D$ is the Dirac delta function. The distributed Lagrangian multiplier $z_i$ is defined in the vicinity $Q_t^{n,e}$ of the interface $R_t$, only at nodes whose support is cut by the interface. This results in the implicit formulation of the interfacial fluid tractions $\mathbf{t}^F$, enabling the method to handle constrained motion of fluid-immersed objects in a computationally elegant way.

Since the thin-immersed structure is coupled to the surrounding fluid domain on its upper (+) and lower (-) faces, the coupling conditions need to be ensured at $R_t^+$ and $R_t^-$. The integral form of the fluid-structure coupling conditions (27) is now reformulated by using the distributed Lagrange multiplier $\mathbf{z}^+$ and $\mathbf{z}^-$ to replace the interfacial tractions $\mathbf{t}^F$:

$$ + \int_{R_t^{n,+}} \delta \mathbf{z}^{F+} \cdot (\mathbf{v}^{F+} - \mathbf{v}^{S}) dR - \int_{R_t^{n,+}} \delta \mathbf{v}^{F+} \cdot \mathbf{z}^{F+} dR + \int_{R_t^{n,+}} \delta \mathbf{v}^{S} \cdot \mathbf{z}^{F+} dR \qquad (32a) $$

$$ + \int_{R_t^{n,-}} \delta \mathbf{z}^{F-} \cdot (\mathbf{v}^{F-} - \mathbf{v}^{S}) dR - \int_{R_t^{n,-}} \delta \mathbf{v}^{F-} \cdot \mathbf{z}^{F-} dR + \int_{R_t^{n,-}} \delta \mathbf{v}^{S} \cdot \mathbf{z}^{F-} dR \qquad (32b) $$

$$ + \sum_{e \in \mathcal{E}} \int_{Q_t^{n,e}} \delta \mathbf{z}^{F+} \cdot \tau_z \mathbf{z}^{F+} + \sum_{e \in \mathcal{E}} \int_{Q_t^{n,e}} \delta \mathbf{z}^{F-} \cdot \tau_z \mathbf{z}^{F-} \qquad (32c) $$

where line (32c) introduces an additional stabilizing term to the Lagrange multiplier formulation, equivalent to the perturbed Lagrange multiplier technique. This stabilization is needed since the distributed Lagrange multiplier field $\mathbf{z}$ is defined to be equal to the interface tractions, but is
<figure><img src="image_11.png" /><figcaption>Figure 10. Distributed Lagrange multiplier in prismatic element: (a) domain of integration for isosurface ($R_t^n$) and enriched elements ($Q_t^{n,e} \in \mathcal{L}$); (b) exemplary distributed multiplier $z$; and (c) exemplary interface tractions $z$ along $R_t^n$.</figcaption></figure>

undefined elsewhere in the domain. Uniqueness in the solution is obtained by evaluating the perturbation within space-time finite elements $Q_t^{n,e}$ of the set $\mathcal{L}$ of elements cut by the fluid-structure interface. The perturbation parameter $\tau_z$ is determined sufficiently small by using a heuristic approach based on element size and properties of the fluid.

Figure 10 emphasizes the application of the distributed Lagrange multiplier approach. Numerical integration of the terms in Equation (32) is performed over the domain of the space-time element and the space-time interface, captured by a zero level set. The issue of a specialized numerical quadrature technique is addressed in Section 5.3. For instance, nodal values $\hat{z}$ for the bi-linear space-time shape functions (see upper row if Figure 9) approximating the field of the distributed Lagrange multiplier $z$ within the element are chosen. With the help of the distributed Lagrange multiplier field $z(\mathbf{x}, t)$, defined within the whole space-time finite element, the value of the local Lagrange multiplier on the interface can be determined by interpolating $z(\mathbf{x}, t)$ at the position of the interface/isosurface. The resulting exemplary distribution of interfacial traction values is shown: the interface traction may vary in space and time on the evolving interface.

## 4.4. General monolithic solution strategy

The weak form of the coupled system consisting of the fluid (Section 3.1) and the thin structure (Section 3.2), combined with the space-time discretization of a specific time slab $Q^n$, leads to the system of non-linear algebraic equations:

$$ \hat{\mathbf{A}}(\hat{\mathbf{x}})\hat{\mathbf{x}}=\hat{\mathbf{b}}(\hat{\mathbf{x}}) \qquad (33) $$

with the coefficient matrix $\hat{\mathbf{A}}$, the vector of unknowns $\hat{\mathbf{x}}$ and the right-hand side $\hat{\mathbf{b}}$. The non-linear equations are linearized applying a fix point Picard iteration scheme that results in the system

$$ \hat{\mathbf{A}}(\hat{\mathbf{x}}^{m-1})\hat{\mathbf{x}}=\hat{\mathbf{b}}(\hat{\mathbf{x}}^{m-1}) \qquad (34) $$
to be solved within the mth iteration step of the actual time slab n. The components of the linearized equations of the coupled fluid-structure system

$$ \begin{bmatrix} \hat{M}^S & \hat{S} & 0 & 0 \\ -\hat{S}^T & \hat{Z} & \hat{F}^T & 0 \\ 0 & -\hat{F} & \hat{M}^F & -\hat{H}^T \\ 0 & 0 & \hat{H} & 0 \end{bmatrix} \begin{bmatrix} \hat{v}^S \\ \hat{z} \\ \hat{v}^F \\ \hat{p} \end{bmatrix} = \begin{bmatrix} \hat{b}^S \\ 0 \\ \hat{b}^F \\ 0 \end{bmatrix} \quad (35) $$

are the unknown vectors $\hat{v}^S$, $\hat{v}^F$, $\hat{p}$ and $\hat{z}$ representing structural velocities, fluid velocities, fluid pressure and the distributed Lagrange multipliers, respectively. These vectors of discrete field variables may also contain additional unknowns resulting from the space-time enriched approximation introduced in Section 4.1. External forces on solid and flow are considered by $\hat{b}^S$ and $\hat{b}^F$.

In the computation of the linear structural mass matrix $\hat{M}^S$, sectional-driven inertia forces are included, while inertia effects originating from rotational velocity of the thin solid section are neglected. The mass matrix $\hat{M}^F$ of the Navier-Stokes fluid is linearized in the convective term, by using the velocity $\hat{v}^{m-1}$ within Equation (25a) for the mth iteration step of the utilized outer Picard iteration scheme. The matrix $\hat{H}$ corresponds to the velocity-pressure term needed for the incompressibility constraint of the fluid. The matrices $\hat{S}$ and $\hat{F}$ are the discretized fluid-structure coupling terms in Equations (32a) and (32b). The symmetric matrix $\hat{Z}$ represents the stabilization term (32c) of the distributed Lagrange multiplier approach introduced in Section 4.3. Sub-matrices from the Galerkin/least-squares stabilization used for the convection dominated incompressible fluid are not explicitly shown in (35), but considered in the computational scheme.

The overall solution strategy for N time slabs $Q^n$ of the coupled non-linear fluid-structure problem is shown in Algorithm 1.

<caption>**Algorithm 1.** Computational scheme of the EST method.</caption>

**Data:** finite element mesh of fluid and thin-walled structure, material properties
**Input:** initial conditions $\tilde{v}$, boundary conditions $\bar{v}$ and $\bar{t}$ for fluid and structure
**for** $n < N$ **do** solve time slab $Q^n$:
| **if** $n=0$ **then** $v^{n-1}=\tilde{v}$; use initial conditions
| **else** $\hat{v}^{n-1}=\hat{v}^{n-1,m}$; use converged last time slab solution
| **repeat** perform Picard iteration step $m$ for nonlinearities:
| | **if** $m=0$ **then** $\hat{v}^{n,m-1}=\hat{v}^{n-1,m}$; use converged last time slab solution
| | | (1) calculate actual structural displacements $\hat{d}^{n,m}$ with Equation (19) $\rightarrow L_t^S$;
| | | (2) map $L_t^S$ to $Q_t^F$ by constructing $\phi^{n,m}$ with help of Equation (36) in Section 5.1;
| | | (3) assign immersed structural elements to hosting fluid elements (see Section 5.2);
| | | (4) determine the need for space-time fluid enrichment element-by-element;
| | | (5) initialize system with additional degrees of freedom from enrichment;
| | | (6) numerical integration of the weak form (see Section 5.3);
| | | (7) apply boundary conditions $\bar{v}$ and $\bar{t}$;
| | | (8) solve $m$-th linearized system of equations, Equation (34) $\rightarrow \hat{v}^{n,m}, \hat{p}^{n,m}, \hat{z}^{n,m}$;
| **until** convergence of physical solution and level set function;
| **end**
# 5. ALGORITHMIC DETAILS

## 5.1. Construction of the structure's level set function

The thin structure is embedded into the flow field by using a level set function as in [16]. According to the current configuration of the structure, the level set function $\phi$ is derived following

$$ \phi(\mathbf{x}, t) = \pm \min_{\mathbf{x} \in \Sigma} (\|\mathbf{x}(t) - \mathbf{x}^{\Sigma}(t)\|_2) \quad (36) $$

where the scalar value $\phi$ determines the smallest signed Euclidean distance from a fluid point to the thin structure (line-shaped interface for two-dimensional problems).

Figure 11 shows a sketch of the fluid discretization and an embedded thin structure. The current position of the structure is approximated by the level set function $\phi$ and therefore mapped onto the Eulerian fluid domain. As a result of evaluating (36), at each node of the fluid mesh the value of the level set function represents the relative position to the Lagrangian structure, located at $\phi=0$.

The space-time level set function $\phi(\mathbf{x}, t)$ is interpolated using bi-linear space-time shape functions $N_k^F$ within the prismatic fluid element. It is updated after each non-linear iteration step $m$ and time slab $n$, so that the current position of the Lagrangian structure related to the Eulerian fluid is known at all times.

Mapping a geometrically complex thin structure to a coarse fluid mesh with the help of a level set function may introduce an approximation error $\varepsilon$ depending on the resolution of the fluid discretization. For immersed thin structures of moderate curvature, a fluid meshing taking into account the necessary length scales of dominating flow characteristics seems to be sufficient.

## 5.2. Connectivity update technique

The extended space-time finite element technique avoids repeated updates of the geometrical meshes since the fluid discretization remains topologically fixed, and the structure is immersed by an extension of the physical approximation based on level sets. Therefore, the computational effort is shifted from mesh-update and remeshing cycles to the update of connectivity between the degrees of freedom of the solid velocities, the locally enriched fluid velocities/pressure and the distributed Lagrange multipliers representing interfacial coupling tractions.

Figure 12 shows gradual changes of connectivity (including enrichment degrees of freedom) in the overall algebraic system of equations for the piston problem discussed in Section 6.1

<figure><img src="image_12.png" /><figcaption>Figure 11. Mapping the structural current configuration to the fluid frame.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12. Connectivity changes during calculation of moving immersed structure.</figcaption></figure>

after utilizing the reverse Cuthill–McKee reordering algorithm to the nodes of the geometrical
mesh. During the harmonic oscillation of the spring-supported piston with cycle period *T* and its
crossing of several space–time fluid elements, the degrees of freedom of structural velocity need
to be reconnected to the fluid velocities and interface tractions of the actual fluid element hosting
the immersed structure. The sparsity pattern of the coefficient matrix emphasizes the effect of
distributed Lagrange multipliers to enforce coupling conditions between the fluid and the evolving
immersed solid.

While more sophisticated methods are available (quadtree technique, etc.) for the two-
dimensional implementation of the presented approach, a brute force search algorithm is applied
to determine the current neighborhood and relative positions of structure and fluid elements within
each non-linear Picard iteration *m* of a time step *n*. After the search algorithm has been performed
globally, each space–time fluid element *Q*<sub>*t*</sub><sup>*n*,*m*,*e*</sup> owns an updated list *L*<sub>S</sub><sup>*n*,*m*,*e*</sup> of structural elements
currently passing through. Fluid elements not touched by one or multiple structural elements
have an empty list of assigned immersed solid elements and do not initiate any enrichment of the
approximation. If a space–time fluid element is hosting a part of the structure, all elemental nodes
*j* are assigned to the set of enriched nodes *M*.

5.3. Numerical integration of the weak form with enrichment

The fluid weak form has to be integrated over the flow domain and therefore affects non-enriched,
fully and partially EST fluid elements. Each non-enriched space-time flow element is integrated
applying a standard Gaussian quadrature rule. If the element is enriched by Equations (28), (29),
the integrands of the weak form are discontinuous inside the element. To perform the integration
correctly, the quadrature algorithm has to be modified in these elements: in a first step the prismatic
space-time element is regularly subdivided into tetrahedra as shown in Figure 13. The second step
subdivides the resulting tetrahedra at the zero iso-contour of the linearly interpolated space-time
level set function $\phi(\mathbf{x}, t)$ into the final set of sub-tetrahedra for numerical integration. Within each
elemental sub-tetrahedra, the integrand of the weak form is then smooth and a Gaussian quadrature
rule can be applied for integration.

The subdivision approach is also used for the numerical integration over the spatial domain
only (triangles with stripe pattern in Figure 13) and the fluid-structure interface is represented by
$\phi(\mathbf{x}, t) = 0$. The interface surface is subdivided into a set of planar triangles that approximates the
<figure><img src="image_14.png" /><figcaption>Figure 13. Subdivision of prismatic space-time finite element.</figcaption></figure>

potentially curved zero iso-contour. Numerical quadrature over the interface is then performed for all sub-triangles using Gaussian integration rules.

## 6. EXAMPLES

### 6.1. Motion of a piston

The piston problem consists of a fluid-filled channel with length $b$ and width $a$ and a matching piston-like object modeled by a thin structure (Figure 14). The piston is supported by linear springs of stiffness $k$ and separates the fluid domain with density $\rho_f$ and viscosity $\mu$. The thickness $h$ of the structure, placed at $c$ from the left inlet and of density $\rho_s$, is very small with respect to other spatial system dimensions. The structural displacements are measured in relation to $c$ using $u_1(t)$.

At the boundary of the channel, Dirichlet and Neumann conditions are defined. The vertical velocity $\bar{v}_2$ at the walls of the channel is set to zero. At the left inlet the horizontal tractions $\bar{t}_1^L$ and at the outlet $\bar{t}_1^R=0$ are prescribed. The system is initially at rest, so that for fluid and structure it holds that $v=0$ that and $u=0$.

The physical behavior of the described system is investigated considering the setups R and F of a rigid and flexible piston structure. See Table I for material parameters of fluid and structure. The geometrical dimensions of the system are given by $a=1$ cm, $b=8$ cm, $c=2$ cm and $h=0.02$ cm. The ratio $\xi=m_f/m_s$ of the total fluid mass $m_f=\rho_f ab$ to the total structural mass $m_s=\rho_s ah$ is $\xi=10^3$ and characterizes the piston problem as a fluid-structure system involving a light-weight structural part.

The analytical solution to the time-dependent pressure field $p(x_1, t)$ within the channel is

$$p(x_1, t) = \begin{cases} \bar{t}_1^L - \rho_f v_{1,t}^A(t) x_1, & x_1 < X^A(t) \\ \bar{t}_1^R + \rho_f v_{1,t}^A(t) (b-x_1), & x_1 > X^A(t) \end{cases} \quad (37)$$

which is strongly discontinuous at the location of the rigid piston

$$X^A(t) = c + u_1^A(t) \quad (38)$$
<figure><img src="image_15.png" /><figcaption>Figure 14. Piston system: channel flow with immersed structure.</figcaption></figure>

<table><caption>Table I. Material parameters for the rigid (R) and flexible (F) piston setup.</caption><thead><tr><th>Setup</th><th>ρ<sub>f</sub> (g/cm<sup>3</sup>)</th><th>μ (g/cms)</th><th>ρ<sub>s</sub> (g/cm<sup>3</sup>)</th><th>E (g/cms<sup>2</sup>)</th><th>ν</th><th>k (g/s<sup>2</sup>)</th></tr></thead><tbody><tr><td>R</td><td>1.25 × 10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>5.0 × 10<sup>-2</sup></td><td>—</td><td>—</td><td>0.5</td></tr><tr><td>F</td><td>1.25 × 10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>5.0 × 10<sup>-2</sup></td><td>2.0 × 10<sup>2</sup></td><td>0.0</td><td>0.5</td></tr></tbody></table>

The pressure gradient depends on fluid density $\rho_f$ and current acceleration $v_{1,t}^A(t)$ of the piston. The pressure gradient remains continuous at the position $X^A(t)$ of the piston. The incompressibility constraint for the fluid leads to a spatially constant velocity field of magnitude $v_1^A(t)$.

**6.1.1. Prescribed motion of rigid piston structure.** The state of the fluid part of the coupled system is computed for given motion of a rigid piston (setup R). The prescribed constant acceleration of the piston $a^A = v_{1,t}^A(t) = 0.08 \text{ cm/s}^{-2}$ results in a quadratic displacement and linearly increasing velocity of the piston body within the time interval $I = [0\text{s}, 10\text{s}]$. Traction (and pressure) at the inlet boundary are set to zero.

The channel domain is discretized by $4 \times 32$ EST fluid elements wherein the piston, described by four space-time structure elements, is immersed. For the simulation of a rigid piston, the degrees of freedom for angular velocity $\omega$ of the structure are set to zero. For comparison of results obtained with the EST method, the coupled problem is additionally solved using a monolithic approach [25] to fluid-structure interaction based on space-time finite elements and a traditional mesh moving scheme. The time step size $\Delta t = 1.0 \text{s}$ is kept constant.

As shown in Figure 15, the EST method captures the pressure field and its evolving discontinuity, which is aligned to the location of the accelerated rigid piston. The numerical results obtained by the EST method are spatially exact compared with the analytical solution given by Equations (37)–(38). The strongly discontinuous pressure solution within single elements is captured by the pressure enrichment. While the spatial solution for a specific time instant $t$ can be reproduced exactly, the displacement of the piston is approximated in time by linear shape functions of the associated space-time structural element. Further, the level set function used for representation of the flow-immersed structure is interpolated by linear functions in space and time. The space-time solution to the discontinuous pressure field therefore converges with second order to the correct solution with decreasing time step size $\Delta t$.

Figure 16 compares numerical results with the pressure state within the flow region resulting from a computation using the EST method and a mesh-moving approach for various time instants.
<figure><img src="image_16.png" /><figcaption>Figure 15. Piston with prescribed motion: discontinuous pressure solution.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16. Piston with prescribed motion: pressure obtained by (a) the EST method and (b) a moving mesh approach.</figcaption></figure>

Both methods are able to capture the non-smooth pressure solution within the flow domain. The significant difference is in the much more elaborate grid setup for the computation using a moving mesh compared with the pre-processing for the EST method, which requires only a discretization of the structure and the fluid domain.

### 6.1.2. Vibration of rigid piston structure.
For the situation of a spring-supported, undamped rigid piston immersed in a fluid channel with tractions $\vec{t}_1^L = 2.0 \text{ g/cm/s}^2$ and $\vec{t}_1^R = 0.0 \text{ g/cm/s}^2$, the natural frequency $\omega_0$ and amplitude $a_0$ of the harmonic vibration are given by

$$ \omega_0 = \sqrt{\frac{2k}{m_f + m_s}} \quad \text{and} \quad a_0 = \frac{a(\vec{t}_1^L - \vec{t}_1^R)}{2k} \qquad (39) $$
<figure><img src="image_18.png" /><figcaption>Figure 17. Vibration of rigid piston: discontinuous pressure solution for one period $T = 2\pi/\omega_0$.</figcaption></figure>

while the horizontal piston displacement, velocity and acceleration read

$$u_1^A(t) = a_0(1 - \cos \omega_0 t) \qquad (40)$$

$$v_1^A(t) = a_0 \omega_0 \sin \omega_0 t \qquad (41)$$

$$v_{1,t}^A(t) = a_0 \omega_0^2 \cos \omega_0 t \qquad (42)$$

Owing to the non-zero fluid traction state $\vec{t}$ at the left inlet of the channel, the piston undergoes accelerated motion. This results in a piecewise linear pressure field along the main axis of the channel. With increasing deflection of the piston from its initial position at $c$, the spring forces pull back the piston and lead to its oscillation. Owing to the absence of dampers and fluid shear stresses, the oscillation of the piston is not damped.

The channel domain is discretized uniformly by $8 \times 64$ enriched prismatic space-time fluid elements. The rigid piston is described by four space-time elements for a thin-walled structure. The constant time step size $\Delta t = 0.1$ s in $I = [0 \text{ s} \dots 20 \text{ s}]$ is sufficient to resolve the harmonic oscillation of the piston.

Figure 17 shows the numerical solution to the displacement and velocity of the piston in point A and the pressure field along the main axis of the channel using the EST approach. The behavior of the fully coupled system is well simulated by the EST approach: the spatial discontinuous pressure solution is found exactly. An approximation error is introduced by the time integration scheme since the used linear velocity within a time slab is not able to capture the analytic solution exactly.

The results of the EST method and the mesh-moving approach are compared in Figure 18. As for the prescribed motion of the rigid piston structure, the flow mesh used in the moving mesh technique deforms considerably but remains valid during the calculation. Mesh deformation due to piston motion introduces here to coarse and fine resolved regions within the flow domain. While this may be sufficient for approximation of very smooth functions, local coarsening as a result of
<figure><img src="image_19.png" /><figcaption>Figure 18. Vibration of rigid piston: pressure obtained by (a) the EST method and (b) a moving mesh approach.</figcaption></figure>

mesh deformation can lead to a significant reduction in accuracy. The EST method based on a fixed fluid mesh circumvents these difficulties. The fluid enrichment resolves the strong pressure discontinuity within the flow field. The solution of the smooth pressure states within the fluid sub-regions coincides with the mesh-moving solution.

6.1.3. Forced vibration of flexible piston structure. The previous example is modified using a flexible piston structure (setup F) and is thereby extended to a spatially two-dimensional problem. The light-weight thin piston is subject to the same fluid loading and shows large structural deformations accompanied by the oscillatory motion observed for the rigid piston. The non-linear model for kinematics of thin-walled structures results in an arc-shape deformed piston and non-zero structural velocities in the vertical direction. Consequently, the fluid velocity is neither constant in horizontal direction nor zero in vertical direction as in the previous example. For more viscous fluids (than that considered here by setup F) this leads to a pressure state of two-dimensional character within the flow domain including a strong discontinuity at the immersed structure. The behavior of the coupled system consisting of Navier-Stokes fluid and an embedded deformable structure is not predictable through a closed solution.

The discretization of the situation is realized by 8 × 64 fluid elements and eight structural elements. Owing to the deforming structural part, the time step of $\Delta t = 0.001$ s is chosen throughout the entire time interval $I = [0\text{s}...20\text{s}]$.

Numerical solutions to horizontal displacement and velocity of points A and B on the piston are shown in Figure 19. Since total mass distribution is not modified and damping effects from the viscous fluid are negligible, the frequency of oscillation is almost unchanged. If the piston reaches its initial configuration during a cycle, inertia forces of the distributed-mass piston lead to secondary oscillations of the thin-walled structure, see velocity plot for points A and B in Figure 19. Both approaches, EST method and mesh-moving technique, give similar results to nodal velocity and displacement.
<figure><img src="image_20.png" /><figcaption>Figure 19. Vibration of flexible piston: EST/moving mesh solution to piston displacement and velocity.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 20. Vibration of flexible piston: pressure obtained by (a) the EST method and (b) a moving mesh approach.</figcaption></figure>

Figure 20 underlines the limitations of the moved mesh techniques for fluid-structure problems with considerable large deformations and rigid body motions. The mesh distortions in the middle of a piston oscillation with maximal displacements come very close to situations involving invalid element geometries. In such a pathological case the mesh-moving technique in the coupled solution procedure fails and requires advanced methods, e.g. local re-meshing of the fluid grid. The EST method avoids these restrictions.

## 6.2. Immersed rotating structures

The presented EST method is able to describe coupled systems involving rotating structural parts. The following example considers two differently shaped thin-walled structures of width $h$ embedded
<figure><img src="image_22.png" /><figcaption>Figure 21. Rotor system: channel flow with immersed flexible and rotating structure (shapes A and B).</figcaption></figure>

in a flow channel of dimensions a=0.2cm and b=2.0cm as shown in Figure 21. In horizontal direction each structure ($\rho_s$=10.0g/cm³) is supported by a linear spring of stiffness k and is coupled to the channel flow ($\rho_f$=1.0g/cm³, $\mu_f$=10⁻³g/cm/s). The enclosed domain inside the solid is filled with fluid ($\rho$=$\rho_f$/2 and $\mu$=$\mu_f$/2). Both structures are free to rotate around their hinges (initial positions given by c=0.2cm and d=0.05cm) and are fixed in the vertical direction. At the left boundary of the channel, the horizontal velocity $\bar{v}_1^L$=1.0cm/s is imposed while on the right the fluid boundary tractions $\bar{t}_1^R$ are zero ($p^R$=0).

The curved shapes A and B are generated using the superformula by Gielis [39]:

$$ \frac{1}{r} = \left( \left| \frac{1}{\alpha} \cos \left( \frac{\gamma}{4} \theta \right) \right|^{\lambda_2} + \left| \frac{1}{\beta} \sin \left( \frac{\gamma}{4} \theta \right) \right|^{\lambda_3} \right)^{1/\lambda_1} \quad (43) $$

where r and θ are polar coordinates, α and β are non-zero real numbers and γ, λ₁, λ₂ and λ₃ are real numbers. Table II lists specific values belonging to shapes A and B of the thin-walled structure as well as material parameters of fluid and structural part (rigid and flexible configurations). The maximum radius of the shapes is normalized to $r_{max}$=0.0475. The origin of each shape is identical to the hinge point that is connected to the thin-walled structure by stiffeners as shown in Figure 21. Each shape is composed of $n_p$=120 nodes and linear thin-walled structure elements. Each stiffener consists of one structural element without fluid-immersing properties and therefore not in interaction with the enclosed fluid.

The narrow flow domain around the immersed shapes is regularly discretized by prismatic space-time finite elements with equidistant nodes ($\Delta h$=0.004cm). See Figure 22 for the used mesh. The numerical analysis of the coupled system is performed using a time slab width of $\Delta t$=0.001s for the time interval $T$=[0.0s, 4.0s].

In the beginning of the simulation, the inflow velocity is linearly increased from zero to the final constant $\bar{v}_1^L$ within the time interval $T_v$=[0.0s, 0.5s]. As a consequence the flow field develops ($Re$≈50) and the fluid accelerates the two spring supported structures. The time-dependent fluid forces act on the structures enclosing the lighter fluid and the asymmetric flow situation causes rotational motion of both objects additional to motion in the horizontal direction. Figure 23 shows the solution to horizontal displacements and angle of rotation of the hinges for the rigid and flexible
<table><caption>Table II. Material and shape parameters (sets A and B) for rigid (R) and flexible (F) rotating structures.</caption><thead><tr><th>Setup</th><th><i>E</i> (g/cm s<sup>2</sup>)</th><th><i>ν</i></th><th><i>k</i> (g/s<sup>2</sup>)</th><th><i>h</i> (cm)</th></tr></thead><tbody><tr><td>R</td><td>—</td><td>—</td><td>1.0</td><td>1 × 10<sup>-2</sup></td></tr><tr><td>F</td><td>10<sup>4</sup></td><td>0.0</td><td>1.0</td><td>2 × 10<sup>-3</sup></td></tr></tbody></table><table><thead><tr><th>Shape</th><th>α</th><th>β</th><th>γ</th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th><th>λ<sub>3</sub></th></tr></thead><tbody><tr><td>A</td><td>1.0</td><td>1.0</td><td>3.0</td><td>0.1</td><td>1.8</td><td>1.8</td></tr><tr><td>B</td><td>1.0</td><td>1.0</td><td>4.0</td><td>0.8</td><td>4.0</td><td>4.0</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>Figure 22. Details of spatial discretization of flow domain and immersed rotor shapes A and B.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 23. Horizontal displacement and angle of rotation of immersed rotor shapes A (left) and B (right).</figcaption></figure>

cases. The enriched solution fields of pressure and velocity are given in Figure 24 for selected time instants of the configuration involving flexible structures.

Figure 25 compares relative displacement, rotation and deformation of the immersed structures for the rigid and flexible setup under the same flow conditions. For the flexible configuration (F), considerable deformations of both shapes are observable.
<figure><img src="image_25.png" /><figcaption>Figure 24. Enriched pressure (top) and velocity (bottom) fields for setup F at times $t = 1.0, 2.0$ and $3.0s$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 25. Motion (setup R, upper rows) and deformation (setup F, lower rows) of the immersed rotor shapes.</figcaption></figure>

This fluid-structure interaction problem can be classified as instability-induced and movement-induced structural excitations. Initial motion and deformation are caused by the developing fluid instabilities (shear flow) in the typical Kármàn vortex street. Once the flexible structures are set in motion, the fluid field changes and a strong coupling of fluid and structure is present. The behavior of the whole system is highly non-linear due to the motion and deformations of the involved structures in the viscous Navier-Stokes fluid. The formulation and numerical algorithm of the EST
method presented in this paper are capable of handling this non-linear fluid-structure interaction
situation with large structural motion/rotation and deformation.

7. CONCLUSIONS

The present paper proposes an elegant way to treat a fluid-structure interaction problem in the case of thin immersed structures. The Lagrangian structure meshes are arbitrarily embedded in the Eulerian fluid. The time axis is also included in the mesh discretization for both fluid and structures leading to a space-time discretization. The structures are tracked in the space-time domain by the zero iso-contour of a level set function. The discontinuities in pressure and velocity fluid fields are represented by enriching the fluid field approximations using the partition of unity enrichment concept and the XFEM. The continuity of velocity through the fluid-structure interfaces is ensured by a perturbed Lagrange multiplier technique.

The resulting proposed approach, the enriched space-time (EST) method, enables large displacement of the embedded structures without any fluid mesh update. Potential contacts between structures are not treated within the paper. Although the method has been developed for a fluid-structure problem, it can be applied to any problem of time-evolving interfaces involving weak or strong discontinuities. Furthermore, while the chosen enrichments are sign functions, any kind of enrichment functions may be introduced.

The first application compares the numerical and the analytical solutions of a piston in a one-
dimensional flow. It proofs the accuracy of the method to represent a strong discontinuous pressure
field as well as the correct transfer of momentum and geometrical continuity between fluid and
structure. For a flexible piston structure, an alternative mesh-moving scheme reaches its limits due
to large mesh distortions while the EST method is still able to correctly represent the solution.

The second application presents two different structures immersed in a channel flow. The viscous fluid causes horizontal and rotational motion of both objects that undergo considerable deformation. The potential of the EST method is emphasized in this application where classical ALE methods would fail.

ACKNOWLEDGEMENTS

This publication has been realized within a cooperation between the Institut für Statik at the Technische
Universität Braunschweig and the Structural Mechanics and Coupled Systems Laboratory of the Conser-
vatoire National des Arts et Métiers (CNAM) in Paris. The collaboration is supported by the PROCOPE
programme between France (PAI-EGIDE) and Germany (DAAD). A further financial support of this work
by the Department of Computational Sciences in Engineering (CSE) at TU Braunschweig is gratefully
acknowledged.

REFERENCES

1. Nomura T. ALE finite element computations of fluid-structure interaction problems. *Computer Methods in Applied Mechanics and Engineering* 1994; **146**:91–126.

2. Bathe K-J, Zhang H, Ji S. Finite element analysis of fluid flows fully coupled with structural interaction. *Computers and Structures* 1999; **72**(1–3):1–16.

3. Wall WA, Mok DP, Ramm E. Partitioned analysis approach for the transient coupled response of viscous fluids and flexible structures. In *Proceedings of the European Conference on Computational Mechanics*, ECCM'99, Munich, Wunderlich W (ed.), 1999.
4. Kuhl E, Hulshoff S, de Borst R. An arbitrary Lagrangian–Eulerian finite-element approach for fluid–structure interaction phenomena. *International Journal for Numerical Methods in Engineering* 2003; **57**:117–142.

5. Le Tallec P, Mouro J. Fluid-structure interaction with large structural displacements. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**(24–25):3039–3067.

6. Peskin CS, McQueen DM. A three-dimensional computational method for blood flow in the heart. I. Immersed elastic fibers in a viscous incompressible fluid. *Journal of Computational Physics* 1989; **81**(2):372–405.

7. Bertrand F, Tanguy PA, Thibault F. A three-dimensional fictitious domain method for incompressible fluid flow problems. *International Journal for Numerical Methods in Fluids* 1997; **25**(6):719–736.

8. Glowinski R, Pan TW, Hesla TI, Joseph DD. A distributed Lagrange multiplier/fictitious domain method for particulate flows. *International Journal of Multiphase Flow* 1999; **25**(5):755–794.

9. Zhang L, Gerstenberger A, Wang X, Liu WK. Immersed finite element method. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**(21–22):2051–2067.

10. Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* 2001; **50**(4):993–1013.

11. Sethian JA. Level Set Methods and Fast Marching Methods. *Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science*. Cambridge University Press: Cambridge, 1999.

12. Melenk JM, Babuška I. The partition of unity finite element method: basic theory and applications. *Computer Methods in Applied Mechanics and Engineering* 1996; **139**(1–4):289–314.

13. Chessa J, Belytschko T. An enriched finite element method and level sets for axisymmetric two-phase flow with surface tension. *International Journal for Numerical Methods in Engineering* 2003; **58**(13):2041–2064.

14. Sukumar N, Chopp D, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite element method. *Computer Methods in Applied Mechanics and Engineering* 2000; **190**(46–47):6183–6200.

15. Moës N, Dolbow JE, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**(1):131–150.

16. Legay A, Chessa J, Belytschko T. An Eulerian–Lagrangian method for fluid–structure interaction based on level sets. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**(17–18):2070–2087.

17. Legay A, Tralli A. An Euler–Lagrange enriched finite element approach for fluid–structure interaction. *European Journal of Computational Mechanics* 2007; **16**(2):145–160.

18. Kölke A, Walhorn E, Hübner B, Dinkler D. Numerical modeling of fluid-structure interaction with free surface flows. *Proceedings in Applied Mathematics and Mechanics* 2003; **3**(1):412–413 (Available from: http://www3.interscience.wiley.com/cgi-bin/fulltext/106571654/PDFSTART).

19. Chessa J, Belytschko T. Arbitrary discontinuities in space-time finite elements by level sets and X-FEM. *International Journal for Numerical Methods in Engineering* 2004; **61**(15):2595–2614.

20. Chessa J, Belytschko T. A local space-time discontinuous finite element method. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**(13–16):1325–1343.

21. Tezduyar TE, Behr M, Liou J. A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial–domain/space–time procedure: I. The concept and the preliminary numerical tests. *Computer Methods in Applied Mechanics and Engineering* 1992; **94**(3):339–351.

22. Hughes TJR, Hulbert GM. Space-time finite element methods for elastodynamics: formulations and error estimates. *Computer Methods in Applied Mechanics and Engineering* 1988; **66**(3):339–363.

23. Tezduyar TE, Behr M, Mittal S. A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial–domain/space–time procedure: II. Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders. *Computer Methods in Applied Mechanics and Engineering* 1992; **94**(3):353–371.

24. Hübner B, Walhorn E, Dinkler D. Numerical investigations to bridge aeroelasticity. In *Proceedings of the Fifth World Congress on Computational Mechanics*, Vienna, Mang HA, Rammerstorfer FG, Eberhardsteiner J (eds), 2002. http://wccm.tuwien.ac.at/publications/Papers/fp81407.pdf.

25. Hübner B, Walhorn E, Dinkler D. A monolithic approach to fluid-structure interaction using space-time finite elements. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**(23–26):2087–2104.

26. Kölke A, Dinkler D. Extended space-time finite elements for two-fluid flows in fluid-structure interaction. In *Proceedings of Sixth World Congress on Computational Mechanics*, Beijing, Yao ZH, Yuan MW, Zhong WX (eds), 2004.

27. Walhorn E, Kölke A, Hübner B, Dinkler D. Fluid-structure coupling within a monolithic model involving free surface flows. *Computers and Structures* 2005; **83**(25–26):2100–2111.
28. Kölke A, Dinkler D. Extended space-time finite elements for boundary-coupled multi-field problems on fixed grids. In *International Conference on Computational Methods for Coupled Problems in Science and Engineering*, Santorini, Papadrakakis M, Oñate E, Schrefler B (eds), 2005.

29. Legay A, Kölke A. An enriched space-time finite element method for fluid-structure interaction—Part I: prescribed structural displacement. *Proceedings of the Third European Conference on Computational Mechanics*, Lisbon, 2006.

30. Kölke A, Legay A. An enriched space-time finite element method for fluid-structure interaction—Part II: thin flexible structures. *Proceedings of the Third European Conference on Computational Mechanics*, Lisbon, 2006.

31. Felippa C. *The TL Timoshenko Plane Beam Element*, Chapter 9, Nonlinear Finite Element Methods (ASEN 5107), 2007 (Available from: http://www.colorado.edu/engineering/CAS/courses.d/NFEM.d).

32. Argyris JH, Scharpf D. Finite elements in space and time. *Journal of the Royal Aeronautical Society* 1969; **73**:1041–1044.

33. Masud A, Hughes TJR. A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems. *Computer Methods in Applied Mechanics and Engineering* 1997; **146**(1–2):91–126.

34. Kölke A. Modellierung und Diskretisierung bewegter Diskontinuitäten in randgekoppelten Mehrfeldsystemen. *Ph.D. Thesis*, Technische Universität Braunschweig, 2005 (Available from: http://www.digibib.tu-bs.de/?docid=00001723).

35. Chessa J, Wang HW, Belytschko T. On the construction of blending elements for local partition of unity enriched finite elements. *International Journal for Numerical Methods in Engineering* 2003; **57**(7):1015–1038.

36. Ji H, Dolbow J. On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method. *International Journal for Numerical Methods in Engineering* 2004; **61**:2508–2535.

37. Moës N, Béchet E, Tourbier M. Imposing Dirichlet boundary conditions in the extended finite element method. *International Journal for Numerical Methods in Engineering* 2006; **67**:1641–1669.

38. Hansbo P, Lovadina C, Perugia I, Sangalli G. A Lagrange multiplier method for the finite element solution of elliptic interface problems using non-matching meshes. *Numerische Mathematik* 2005; **100**:91–115.

39. Gielis J. A generic geometric transformation that unifies a wide range of natural and abstract shapes. *American Journal of Botany* 2003; **90**(3):333–338.