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

HAL
open science

# A Nitsche-based formulation for fluid-structure interactions with contact

Erik Burman, Miguel Angel Fernández, Stefan Frei

► To cite this version:

Erik Burman, Miguel Angel Fernández, Stefan Frei. A Nitsche-based formulation for fluid-structure interactions with contact. ESAIM: Mathematical Modelling and Numerical Analysis, In press. hal-01784841v2

HAL Id: hal-01784841

https://inria.hal.science/hal-01784841v2

Submitted on 7 May 2018

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.
informatics mathematics
Inria

A Nitsche-based
formulation for
fluid-structure
interactions with contact

Erik Burman, Miguel A. Fernández, Stefan Frei

RESEARCH
REPORT

Nº 9172

May 2018

Project-Team REO
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

A Nitsche-based formulation for
fluid-structure interactions with contact

Erik Burman*, Miguel A. Fernández†‡, Stefan Frei*

Project-Team REO

Research Report n° 9172 — May 2018 — 38 pages

**Abstract:** We derive a Nitsche-based formulation for fluid-structure interaction (FSI) problems with contact. The approach is based on the work of Chouly and Hild [SIAM Journal on Numerical Analysis. 2013;51(2):1295–1307] for contact problems in solid mechanics. We present two numerical approaches, both of them formulating the FSI interface and the contact conditions simultaneously in equation form on a joint interface-contact surface. The first approach uses a relaxation of the contact conditions to allow for a small mesh-dependent gap between solid and wall. The second alternative introduces an artificial fluid below the contact surface. The resulting systems of equations can be included in a consistent fashion within a monolithic variational formulation, which prevents the so-called “chattering” phenomenon. To deal with the topology changes in the fluid domain at the time of impact, we use a fully Eulerian approach for the FSI problem. We compare the effect of slip and no-slip interface conditions and study the performance of the method by means of numerical examples.

**Key-words:** fluid-structure interaction, contact mechanics, Eulerian formalism, Nitsche's method, slip conditions.

RESEARCH CENTRE
PARIS

2 rue Simone Iff - CS 42112
75589 Paris Cedex 12

<footnote>* University College London, London WC1E 6BT, United Kingdom</footnote>

<footnote>† Inria Paris, 75012 Paris, France</footnote>

<footnote>‡ Sorbonne Université, Lab. Jacques-Louis Lions, 75252 Paris, France</footnote>
Une méthode de Nitsche pour des problèmes d'interaction
fluide-structure avec contact

**Résumé** : Nous dérivons une méthode de type Nitsche pour des problèmes d'interaction fluide-structure (IFS) avec contact. L'approche est basée sur la méthode proposée par Chouly et Hild [*SIAM Journal on Numerical Analysis. 2013; 51 (2): 1295-1307*] pour les problèmes de contact en mécanique du solide. Nous présentons deux méthodes numériques qui traite le couplage IFS et les conditions de contact simultanément. La première approche utilise une relaxation des conditions de contact dépendant du maillage. La deuxième alternative introduit un fluide artificiel en dessous de la paroi de contact. Les systèmes d'équations qui en résultent peuvent être inclus, de façon consistante et sans inconnues supplémentaires, dans une formulation variationnelle IFS monolithique, ce qui empêche le phénomène dit de "chattering". Les changements de topologie dans le domaine des fluides au moment du contact sont traités par une formulation entièrement eulérienne du problème IFS. Nous comparons l'effet des différentes conditions d'interface, avec et sans glissement, et étudions la performance de la méthode dans quelques exemples numériques.

**Mots-clés** : interaction fluide-structure, mécanique du contact, formalisme Eulerien, méthode de Nitsche, conditions de glissement.
# 1 Introduction

Contact problems have to be considered in many physical processes in engineering, medicine and nature. To name only a few consider for example the contact of balls and races in roller bearings, closing and opening heart valves or falling and jumping bouncing balls. While an extensive amount of literature exists for the numerical simulation of contact in a purely mechanical context (see for example Wohlmuth [1] for an overview), i.e. disregarding the gas or liquid that mostly lies between contacting structures, much less works can be found considering full fluid-structure interaction with contact. The flow between contacting surfaces might however be of great importance for the contact dynamics. In the example of heart valves, the pulsating blood flow is even the driving force that enables opening and closure. In the case of ball bearings, fluid forces in the lubricant between ball and bearing may have a significant influence on the performance and wear of the bearing [2, 3].

Contact between different structures is typically formulated by means of variational inequalities [1]. In the context of full fluid-structure interaction (FSI), first results and algorithms can be found using either an artificial penalty force [4] or Lagrange multipliers [5, 6, 7] to obtain a well-posed and computationally feasible variational formulation.

However, these approaches have certain drawbacks: The use of a penalty force prevents real contact. The force is typically an artificial force and involves the choice of penalty parameters. If they are chosen too small, the structures might overlap in a numerical simulation. If they are chosen too large, the contact dynamics might be significantly perturbed [8]. In the case of Lagrange multipliers, additional variables are introduced on the contacting surfaces and an inf-sup condition is needed to ensure the well-posedness of the system. To tackle the variational inequality numerically, an extra loop is usually used in each time step (for example within an active-set strategy), meaning in particular that the system of equations has to be solved several times in each time step.

Recently, a new contact formulation using Nitsche's method [9] was derived by Chouly and co-workers [10, 11, 12], following the ideas of Alart and Curnier [13]. Their approach is based on an equivalent re-formulation of the contact conditions in equality form. In the case of the contact of an elastic body with a wall, the mechanical contact conditions on the contact surface $\Gamma_C(t)$ read

$$d_n \le 0, \quad \sigma_{s,n}(d) \le 0, \quad \sigma_{s,n}(d)d_n = 0,$$

where $d_n$ denotes the solid displacement in normal direction and $\sigma_{s,n}(d)$ is the normal stress component. It can be shown that these conditions are equivalent to the equality

$$\sigma_{s,n}(d) = -\frac{1}{\gamma} \max\{0, d_n - \gamma\sigma_{s,n}(d)\} \quad \text{on } \Gamma_C(t) \quad (1)$$

for arbitrary $\gamma > 0$ (Chouly & Hild [10]). The authors incorporate this equality condition weakly in the variational formulation using Nitsche's method. This approach has the advantage that it is fully consistent and hence the contact dynamics will not be perturbed. Furthermore, no additional variables have to be introduced and no additional loop within each time-step is needed. The authors were able to prove numerical convergence in a series of papers for friction-free and frictional contact [11, 14, 15], disregarding however the fluid that usually lies between the structures.

Following these works, Burman et al. used the re-formulation (1) to derive a Galerkin Least Squares formulation in equality form for the obstacle problem [16] and a Galerkin Least Squares and a Lagrange multiplier formulation for membrane contact [17].
Modelling of contact in an FSI context brings along a further issue: It is unclear, whether the incompressible Navier-Stokes equations are an appropriate model in the fluid part, when it comes to contact. Theoretical studies show, that for a smooth, rigid solid body, no contact with an exterior wall can happen, when no-slip conditions are used on the interface and the outer boundary of the fluid domain, see Hillairet [18] and Hesla [19] in 2 space dimensions and Hillairet et al. [20, 21] in 3 space dimensions. This changes, when slip- or Navier-slip conditions are used on both the interface and the wall [21] or when the boundary of the solid is non-smooth [22, 23]. Gerard-Varet and Hillairet [22] found in a model example that it comes to contact for a solid with a $C^{1,\alpha}$-parametrised boundary for $\alpha < 1/2$, while no contact happens for $\alpha \ge 1/2$. In the context of fluid-structure interactions, the regularity of the solid boundary depends on the solid displacement $d$, for which such a regularity can usually not be guaranteed.

For a full FSI problem with a thin-walled structure, a no-collision result has been shown by Grandmont and Hillairet [24] in the no-slip case. For an overview on further results regarding existence of fluid-structure interaction problems, we refer to Grandmont et al [25]. Recently, Muha & Čanić [26] showed the well-posedness of a fluid-structure interaction system with slip-conditions.

Motivated by these theoretical results, we will study both no-slip and slip conditions on the FSI interface $\Gamma(t)$ in this work. It will turn out that the latter transits naturally into a “no-friction” condition when it comes to contact, while the prior leads to frictional contact. In this work, we will therefore consider friction-free contact, when a slip-condition is used on $\Gamma(t)$ and the specific frictional contact condition that follows from the transition of the interface conditions, when a no-slip condition is used on $\Gamma(t)$. For recent works on the incorporation of different friction laws (in particular Coulomb and Tresca friction), we refer to Chouly and coworkers [14, 15]. Moreover, we will study only contact of a deformable elastic structure with a fixed and straight wall for simplicity. Efficient algorithms to treat contact between more complex structures can be found for example in Puso [27] or Yang et al [28].

The change of topology in the fluid domain causes additional numerical difficulties. Standard numerical approaches as the Arbitrary Lagrangian Eulerian method are not able to deal with topology changes, as the map from the reference domain to the Eulerian domain degenerates necessarily in this situation. The dynamics shortly before the impact can only be handled robustly, when a Eulerian description of the fluid equations is used.

In the last years, several numerical approaches have been developed that are able to deal with topology changes. The methods can be split into two categories, according to the coordinate systems that are used for the solid system: Fully Eulerian approaches, where also the structure equations are formulated in Eulerian coordinates [29, 30, 31, 8, 32]; and Euler-Lagrangian techniques, where Lagrangian coordinates are used for the solid equations [33, 34, 35, 36, 6, 37, 38, 39, 40], such as in the *Immersed Boundary* or *Immersed Finite Element* methods.

Regarding the *Euler-Lagrangian* techniques, one can further distinguish between methods using Lagrange multipliers for the coupling of fluid and structure (Legay et al. [36], Gerstenberger & Wall [37]) and methods based on Nitsche-techniques (Hansbo et al. [41], Burman & Fernández [38], Alauzet et al. [39], Massing et al. [40], Kamensky et al. [42]). For the latter, a theoretical stability and convergence analysis has been derived [38]. The reader is referred to Boilevin-Kayl et al. [43] for a comparative study on the accuracy of some of these approaches.

The FSI approach we use here is based on the monolithic *Fully Eulerian* approach (Dunne & Rannacher [29], Cottet et al. [30], Richter [31], Frei & Richter [8], Hecht & Pironneau [32]). As the complete system of equations is formulated in Eulerian coordinates, the incorporation of contact conditions is straight-forward by means of variational principles. While some of the early works in this context suffered from stability and accuracy issues, see e.g. Dunne [44], accurate and robust discretisation and stabilisation techniques have been developed recently
(Frei & Richter [45, 46, 8], Hecht & Pironneau [32]). We remark, however, that the algorithms we derive to incorporate contact can be combined in a straight-forward way with different FSI coupling techniques, e.g. Fictitious Domain or Immersed Boundary methods.

Concerning discretisation, we allow both for unfitted and fitted finite element approaches. For the unfitted case, so-called “ghost penalty” stabilisations can be used to guarantee the coercivity of the system [47, 38]. In order to simplify the presentation, we concentrate on fitted discretisations in this work and will comment on the unfitted case in a remark. In the numerical examples at the end of this paper, we will use the fitted locally modified finite element method [45].

The remainder of this paper is organised as follows: In Section 2, we first introduce the equations and the contact model. Then, we derive a variational formulation for the problem of an obstacle within the fluid domain in Section 3.1, where we already have to deal with variational inequalities, but a topology change in the fluid domain is avoided. In Section 3.2, we study contact with an exterior wall and discuss in particular the effect of slip- and no-slip interface conditions. Then, we show a stability result in Section 4. We show detailed numerical studies in Section 5, investigating the influence of contact parameters, interface conditions and different contact formulations as well as convergence under mesh refinement. We conclude in Section 6.

## 2 Model

We begin by presenting the models for the fluid part, the solid part and the fluid-structure interaction on one hand in Section 2.1 and the contact model in Section 2.2 on the other hand. For both models, Nitsche-based variational formulations are introduced. Different possibilities to combine the two models will then be presented in Section 3.

### 2.1 Fluid-structure interaction without contact

We consider a fluid-structure interaction problem that is given on an overall domain $\Omega \subset \mathbb{R}^2$ which is split into a (variable) sub-domain $\Omega_f(t)$ occupied by a viscous fluid, a sub-domain $\Omega_s(t)$ occupied by an elastic solid and a lower-dimensional interface $\Gamma(t)$ separating them, such that

$$ \Omega = \Omega_f(t) \cup \Gamma(t) \cup \Omega_s(t). $$

The boundary of the fluid domain is partitioned as follows $\partial\Omega_f(t) = \Gamma_{fsi}(t) \cup \Gamma_f^D \cup \Gamma_f^N$, where $\Gamma_{fsi}(t)$ stands for the fluid-solid interface. As regards the solid boundary, we assume that $\partial\Omega_s(t) = \Gamma(t) \cup \Gamma_s^D \cup \Gamma_s^N$, where the boundary part $\Gamma(t) = \Gamma_{fsi}(t) \cup \Gamma_C(t)$ is decomposed in the terms of $\Gamma_{fsi}(t)$ and the contact zone $\Gamma_C(t)$ (see the left sketch of Figure 1). The restriction to two dimensions is made only to simplify the presentation. A generalisation to three space dimensions is straight-forward.

In this work we will use a Eulerian description for the complete FSI problem. As already mentioned in the introduction this is not necessary for the contact algorithms derived below, but one convenient way to deal with (possible) topology changes in the fluid domain $\Omega_f(t)$. In an Eulerian description, the solid sub-domain and the interface are implicitly defined by the (unknown) solid displacement $d$

$$ \Omega_s(t) = \{x \in \Omega \mid T(x, t) \in \Omega_s(0)\}, \quad \Gamma(t) = \{x \in \Omega \mid T(x, t) \in \Gamma(0)\}, \quad (2) $$

where $T: \Omega(t) \to \Omega$ is a bijective map, that is given by $T(x, t) = x - d(x, t)$ in the solid domain $\Omega_s(t)$ and by an arbitrary (smooth) extension in $\Omega_f(t) = \Omega \setminus (\Omega_s(t) \cup \Gamma(t))$. For the details, we refer to the textbook of Richter [48] or Frei [49].
In the variable fluid sub-domain $\Omega_f(t)$, we consider the linear incompressible Stokes equations

$$ \partial_t u - \operatorname{div} \sigma_f(u, p) = f_f, \quad \operatorname{div} u = 0, $$

where the Cauchy stress tensor $\sigma_f$ is defined by

$$ \sigma_f(u, p) = \nu_f (\nabla u^T + \nabla u) - pI, $$

$u$ denotes the fluid velocity, $p$ stands for pressure and $\nu_f > 0$ is a constant viscosity. In the solid-subdomain $\Omega_s(t)$, we assume a linear elastic material

$$ \partial_t \dot{d} - \operatorname{div} \sigma_s(d) = f_s, \quad \partial_t d = \dot{d}, $$

where the Cauchy stress tensor $\sigma_s$ is given by

$$ \sigma_s(d) = 2\mu_s E(d) + \lambda_s \operatorname{tr}(E(d)), \quad E(d) = \frac{1}{2} (\nabla d + \nabla d^T), $$

$\dot{d} = \partial_t d$ denotes the solid velocity and $\lambda_s, \mu_s > 0$ are positive constants.

For the coupling across the fluid-solid interface $\Gamma_{\text{fsi}}(t)$, the continuity of velocities and normal stresses

$$ \dot{d} = u, \quad \sigma_f(u, p)n = \sigma_s(d)n \quad \text{on} \quad \Gamma_{\text{fsi}}(t) $$

is typically considered for viscous fluids, where $n = n_s$ denotes the outer normal vector of the solid domain. We recall that, since in this section there is no contact in the solid (i.e., $\Gamma_C(t) = \emptyset$), we have $\Gamma(t) = \Gamma_{\text{fsi}}(t)$. When it comes to contact, it is however questionable, whether this condition is still a reasonable approximation of the underlying physics. Theoretical works [22] show, that the Navier-Stokes equations in combination with these “no-slip” conditions do not allow for contact. Therefore, we will study slip-conditions in this work as well (see Section 3.3), where the continuity across $\Gamma_{\text{fsi}}(t)$ is only imposed for the normal velocity

$$ u \cdot n = \dot{d} \cdot n, \quad \tau^T \sigma_f n = 0, \quad \sigma_f n = \sigma_s n \quad \text{on} \quad \Gamma_{\text{fsi}}(t). \qquad (3) $$

In order to close the system of equations, we define exterior boundary conditions for the fluid and solid

$$ u = 0 \quad \text{on} \quad \Gamma_f^D, \quad \sigma_f(u, p)n = 0 \quad \text{on} \quad \Gamma_f^N $$

$$ d = 0 \quad \text{on} \quad \Gamma_s^D, \quad \sigma_s(d)n = 0 \quad \text{on} \quad \Gamma_s^N $$

and the initial conditions

$$ u(x, 0) = u^0(x) \text{ in } \Omega_f(0), \quad d(x, 0) = d^0(x), \quad \dot{d}(x, 0) = \dot{d}^0(x) \text{ in } \Omega_s(0). $$

We introduce the finite element spaces $V_h$, $Q_h$ and $W_h$ on a quasi-uniform family of triangulations $(\mathcal{T}_h)_{h>0}$ and use Nitsche's method to combine both equations and interface conditions into a monolithic variational formulation (see e.g. Hansbo et al. [41], Burman & Fernandez [38]). Denoting by $n = n_s$ the outer normal vector of the solid domain $\Omega_s(t)$, the monolithic system of equations reads in the no-slip case: Find $u(t) \in V_h, p(t) \in Q_h, d(t) \in W_h$, such that $\dot{d} = \partial_t d$ and

$$ \mathcal{A}_{\text{fsi, no-slip}}(u, p, d, \dot{d})(v, q, w) = (f_f, v)_{\Omega_f(t)} + (f_s, w)_{\Omega_s(t)} \quad \forall v, q, w \in V_h \times Q_h \times W_h, \quad (4) $$
where

$$ \mathcal{A}_{\text{fsi, no-slip}}(u, p, d, \dot{d})(v, q, w) := (\partial_t u, v)_{\Omega_f(t)} + (\sigma_f(u, p), \nabla v)_{\Omega_f(t)} + (\text{div } u, q)_{\Omega_f(t)} \\ + S(p, q) + (\partial_t \dot{d}, w)_{\Omega_s(t)} + (\sigma_s(d), \nabla w)_{\Omega_s(t)} - (T_f(u, p, \dot{d}), w - v)_{\Gamma(t)} - (\dot{d} - u, \sigma_f(v, -q)n)_{\Gamma(t)} \quad (5) $$

where the numerical fluid traction on the interface is defined by $T_f(u, p, \dot{d}) := \sigma_f(u, p)n - \gamma_{\text{fsi}}(\dot{d} - u)$. The Nitsche parameter is chosen as $\gamma_{\text{fsi}} := \gamma_{\text{fsi}}^0 \nu_f h^{-1}$. The term $S(p, q)$ stands for a pressure stabilisation term that is non-zero in case that the discrete fluid spaces do not fulfil a discrete inf-sup condition. Note that at the FSI interface, we have used the following relation for the interface terms arising from integration by parts

$$ (\sigma_s n, w)_{\Gamma(t)} - (\sigma_f n, v)_{\Gamma(t)} = (\sigma_f n, w - v)_{\Gamma(t)} + (\llbracket \sigma n \rrbracket, w)_{\Gamma(t)}, \quad (6) $$

where we have dropped the dependencies of $\sigma_f$ and $\sigma_s$ for better readability. In the absence of contact, the jump of the stresses defined by

$$ \llbracket \sigma n \rrbracket := \sigma_s n - \sigma_f n, \quad (7) $$

vanishes everywhere on $\Gamma(t)$. Furthermore, we have added the term

$$ - (\dot{d} - u, \sigma_f(v, -q)n)_{\Gamma(t)} $$

as in Burman & Fernandez [38] for stability reasons.

Using a slip-condition at the FSI interface, the variational formulation reads: Find $u \in V_h, p \in Q_h, d \in W_h$, such that $\dot{d} = \partial_t d$ and

$$ \mathcal{A}_{\text{fsi, slip}}(u, p, d, \dot{d})(v, q, w) = (f_f, v)_{\Omega_f(t)} + (f_s, w)_{\Omega_s(t)} \quad \forall v, q, w \in V_h \times Q_h \times W_h, $$

where

$$ \begin{aligned} \mathcal{A}_{\text{fsi, slip}}(u, p, d, \dot{d})(v, q, w) := & (\partial_t u, v)_{\Omega_f(t)} + (\sigma_f(u, p), \nabla v)_{\Omega_f(t)} + (\text{div } u, q)_{\Omega_f(t)} + S(p, q) \\ & + (\partial_t \dot{d}, w)_{\Omega_s(t)} + (\sigma_s(d), \nabla w)_{\Omega_s(t)} \\ & - (T_f(u, p, \dot{d}) \cdot n, (w-v) \cdot n)_{\Gamma(t)} - ((\dot{d} - u) \cdot n, n^T \sigma_f(v, -q)n)_{\Gamma(t)}. \end{aligned} \quad (8) $$

## 2.2 Contact model without fluid

We assume that the solid is at a positive distance to the boundary at initial time and that contact can only happen with the lower wall (see Figure 1, left sketch, where the situation at contact is shown)

$$ \Gamma_w = \{(x_1, x_2) \in \partial\Omega \mid x_2 = 0\}. $$

We denote the outer normal vector of the fluid domain $\Omega_f(t)$ at $\Gamma_w$ by $n_w = -e_2$. Moreover, let $g_0(x_2) > 0$ be the function describing the initial distance of a point $(x_1, x_2) \in \Gamma(t)$ to the wall $\Gamma_w$.

When contact with $\Gamma_w$ occurs on a part $\Gamma_C(t) \subset \Gamma(t)$, suitable contact conditions are (Alart & Curnier [13], Chouly & Hild [10])

$$ d \cdot n_w \le g_0, \quad \sigma_{s,n} := n_w^T \sigma_s n \le 0, \quad (d \cdot n_w - g_0) \sigma_{s,n} = 0 \quad \text{on } \Gamma(t). \quad (9) $$
<figure><img src="image_4.png" /><figcaption>Figure 1: Left: Body in contact with the wall $\Gamma_w$. Right: Virtual obstacle line $\Gamma_o$ within the fluid domain $\Omega_f$.</figcaption></figure>

The first inequality in (9) ensures that the solid can not pass through $\Gamma_w$, the second inequality describes that the normal stress is zero (in the absence of contact) or negative (during contact) and the third condition is a complementarity condition that guarantees that at least one of the inequalities is “active”.

In the model setting, that is considered here, the normal and tangential vectors $n_w$ and $\tau_w$ are constant. For later purposes we distinguish between the normal and tangential vectors of the wall $n_w$, $\tau_w$ and the normal and tangential vector of the FSI interface $n, \tau$. The vectors are equal on $\Gamma_C(t)$, where the body is in contact with the wall, but different before contact.

For arbitrary $\gamma_C > 0$ the first line in (9) is equivalent to [10, 13]

$$ \sigma_{s,n}(d) = -\gamma_C[d \cdot n_w - g_0 - \frac{1}{\gamma_C}\sigma_{s,n}(d)]_+ := -\gamma_C[P_{\gamma,s}(d)]_+, \quad (10) $$

where $[f]_+ := \max\{f, 0\}$. The equivalence can be shown by simple calculations for each of the cases $P_{\gamma,s}(d) \le 0$ and $P_{\gamma,s}(d) > 0$, see Chouly & Hild [10].

As tangential contact condition, Chouly & Hild [10] used the “no-friction” condition

$$ \tau_w^T \sigma_s n = 0 \quad \text{on } \Gamma_C(t) \quad (11) $$

for the pure solid problem. Choosing $\gamma_C = \gamma_C^0 \mu_s h^{-1}$, the variational formulation reads: Find $d(t) \in W_h$ such that $\dot{d} = \partial_t d$ and

$$ \left( \partial_t \dot{d}, w \right)_{\Omega_s(t)} + (\sigma_s(d), \nabla w)_{\Omega_s(t)} + \gamma_C ([P_{\gamma,s}(d)]_+, w \cdot n_w)_{\Gamma(t)} = (f_s, w)_{\Omega_s(t)} \quad \forall w \in W_h. \quad (12) $$

We will discuss the tangential contact conditions for the case of fluid-structure interaction with contact below.

We close this section by mentioning that Chouly, Hild & Renard [11] proposed a more general contact formulation that makes use of the consistency of the term

$$ (\gamma_C[P_{\gamma,s}(d)]_+ + \sigma_{s,n}(d), \sigma_{s,n}(w))_{\Gamma(t)}. $$

For $\theta \in [-1, 1]$, the contact term can be generalised to

$$ \begin{aligned} \gamma_C ([P_{\gamma,s}(d)]_+, w \cdot n_w)_{\Gamma(t)} & - \theta (\gamma_C[P_{\gamma,s}(d)]_+ + \sigma_{s,n}(d), \sigma_{s,n}(w))_{\Gamma(t)} \\ & = (1 - \theta)\gamma_C ([P_{\gamma,s}(d)]_+, w \cdot n_w)_{\Gamma(t)} + \theta (\gamma_C[P_{\gamma,s}(d)]_+, P_{\gamma,s}(w))_{\Gamma(t)} \quad (13) \\ & \quad - \theta (\sigma_{s,n}(d), \sigma_{s,n}(w))_{\Gamma(t)}. \end{aligned} $$

For $\theta = 0$, we recover the formulation (12). Besides that, the case $\theta = 1$ is of particular interest, as it yields a symmetric formulation, for which a stability result has been shown [11].
# 3 Fluid-structure interaction with contact

We now address the case where both fluid-structure interaction and contact occur, by combining the ideas described in Sections 2.1-2.2. In the next section, we consider first a simplified setting in which the solid can enter into contact with a virtual obstacle invisible to the fluid. This simplification allows us to consider the numerical treatment of the interface conditions without accounting for the issues related to topology changes within the fluid. Then, in Section 3.2, we show how this numerical setting can be extended to model physically realistic contact and fluid-structure interaction.

## 3.1 Virtual obstacle within the fluid domain

In order to simplify the presentation we introduce the combined FSI-contact formulation first for a model problem with a horizontal obstacle line $\Gamma_o$ within the fluid domain (see the right sketch of Figure 1), before we discuss the full FSI-contact problem including topology changes. The obstacle is invisible to the fluid, but an obstacle to the solid. One may consider for example a membrane that is perfectly permeable for the fluid, but not for the solid or a magnetic field below the obstacle that prevents the solid from crossing the line. We assume for simplicity that the obstacle line is parallel to the fluid boundary and denote its distance by $\alpha > 0$. The initial distance of the FSI interface $\Gamma(0)$ to the obstacle line is then given by $g_\alpha(x_2) = g_0(x_2) - \alpha$. We assume that $g_\alpha \ge 0$.

When considering fluid-structure interactions, the body is pre-loaded before it reaches the obstacle by means of the balance of normal forces

$$ \sigma_s n = \sigma_f n \quad \text{on } \Gamma(t). \qquad (14) $$

If the interface $\Gamma(t)$ reaches the obstacle $\Gamma_o$, the additional constraint $d \cdot n_w \le g_\alpha$ has to be fulfilled. This artificial non-penetration condition can be formulated by means of the complementarity conditions

$$ d \cdot n_w \le g_\alpha, \quad \lambda \le 0, \quad (d \cdot n_w - g_\alpha)\lambda = 0. \qquad (15) $$

This gives rise to an additional surface force in terms of the Lagrange multiplier $\lambda$ on $\Gamma(t)$, that acts in the direction $-n_w = e_2$ (normal to $\Gamma_o$). The equality (14) will be perturbed to

$$ \llbracket \sigma n \rrbracket - \lambda n_w = 0 \quad \text{on } \Gamma(t). \qquad (16) $$

Splitting into normal and tangential contributions, (16) reads

$$ n_w^T \llbracket \sigma n \rrbracket - \lambda = 0, \quad \tau_w^T \llbracket \sigma n \rrbracket = 0 \quad \text{on } \Gamma(t). \qquad (17) $$

Using the same trick as in Section 2.2, (15) is equivalent to

$$ \lambda = -\gamma_C \left[ d \cdot n_w - g_\alpha - \frac{1}{\gamma_C} \lambda \right]_+ =: -\gamma_C [P_\gamma(\lambda, d)]_+ \quad \text{on } \Gamma(t). \qquad (18) $$

Eliminating the Lagrange multiplier by using (17), this reads

$$ \lambda = -\gamma_C \left[ d \cdot n_w - g_\alpha - \frac{1}{\gamma_C} \llbracket \sigma_n(u, p, d) \rrbracket \right]_+ =: -\gamma_C [P_\gamma(\llbracket \sigma_n(u, p, d) \rrbracket, d)]_+ \quad \text{on } \Gamma(t), \qquad (19) $$
where we have used the abbreviation $\llbracket \sigma_n \rrbracket := n_w^T \llbracket \sigma n \rrbracket$.

The natural formulation in the discrete setting is to consider the numerical stress $T_f$ in (7) and in (16). Let us derive these conditions first for the no-slip case. Adding the additional surface force to the FSI-Nitsche formulation (4), the discrete variational formulation reads

$$ \mathcal{A}_{\text{FSI, no-slip}}(u, p, d, \dot{d})(v, q, w) - (\lambda n_w, w)_{\Gamma(t)} = (f_f, v)_{\Omega_f(t)} + (f_s, w)_{\Omega_s(t)} \quad (20) \\ \forall v, q, w \in \mathcal{V}_h \times \mathcal{Q}_h \times \mathcal{W}_h. $$

Due to the additional Nitsche interface terms in (16), this formulation includes the interface condition

$$ \llbracket \tilde{\sigma}_n(u, p, d) \rrbracket - \lambda = 0 \quad \text{on } \Gamma(t), \quad (21) $$

where the numerical stress jump in the FSI-Nitsche formulation across the interface is given by

$$ \llbracket \tilde{\sigma}_n(u, p, d) \rrbracket := \sigma_{s,n}(d) - T_f(u, p, \dot{d}) \cdot n_w. \quad (22) $$

Eliminating the Lagrange multiplier by means of (21), the identity (18) reads

$$ \lambda = -\gamma_C \left[ d \cdot n_w - g_\alpha - \frac{1}{\gamma_C} \llbracket \tilde{\sigma}_n \rrbracket \right]_+ =: -\gamma_C [P_\gamma(\llbracket \tilde{\sigma}_n(u, p, d) \rrbracket, d)]_+ \quad \text{on } \Gamma(t). \quad (23) $$

For the definition of numerical stress in the slightly more complicated case of slip-interface conditions, we refer to Section 3.3. We will in the following analyse both contact formulations (19) and (23) and in particular what the different possibilities mean with respect to the weakly imposed interface conditions. In order to avoid too much repetition, we use a general formulation that includes the Lagrange multiplier $\lambda = \lambda(u, p, d)$, keeping in mind that $\lambda(u, p, d)$ will be either chosen as the jump of normal stresses $\llbracket \sigma_n(u, p, d) \rrbracket$ or the jump of normal discrete stresses $\llbracket \tilde{\sigma}_n(u, p, d) \rrbracket$.

Using (18), the variational formulation reads:

**Variational Formulation 1.** Find $u(t) \in \mathcal{V}_h, p(t) \in \mathcal{Q}_h, d(t) \in \mathcal{W}_h$ such that $\dot{d} = \partial_t d$ and

$$ \mathcal{A}_{fsi,*}(u, p, d, \dot{d})(v, q, w) + \gamma_C ([P_\gamma(\lambda, d)]_+, w \cdot n_w)_{\Gamma(t)} = (f_f, v)_{\Omega_f(t)} + (f_s, w)_{\Omega_s(t)} \quad (24) \\ \forall v, q, w \in \mathcal{V}_h \times \mathcal{Q}_h \times \mathcal{W}_h. $$

where the bilinear form is one of the forms $\mathcal{A}_{fsi, \text{no-slip}}$ or $\mathcal{A}_{fsi, \text{slip}}$ defined in (5) and (8), respectively, and the contact parameter is chosen [10] $\gamma_C = \gamma_C^0 \mu_s h^{-1}$.

### 3.1.1 Weakly imposed interface conditions

Let us now analyse which interface conditions on $\Gamma(t)$ are implicitly included in Variational Formulation 1. For simplicity, we analyse the formulation with the bilinear form $\mathcal{A}_{fsi, \text{no-slip}}$ corresponding to no-slip conditions. Therefore, we integrate by parts in (24) and consider only the interface terms by formally neglecting the bulk and inter-element terms. For better readability, we drop all the dependencies of $\sigma_f$ and $\sigma_s$. Collecting all terms with the fluid test function $v$, we obtain

$$ -\sigma_f n + \sigma_f n - \gamma_{\text{fsi}}(\dot{d} - u) = -\gamma_{\text{fsi}}(\dot{d} - u) = 0 \quad \text{on } \Gamma(t), \quad (25) $$

i.e. the kinematic condition $\dot{d} = u$. Next, we collect the interface terms for the solid part $w$ and split into a normal part $(w \cdot n_w)$ and a tangential part $(w \cdot \tau_w)$. We recall that since the
boundary $\Gamma_w$ is flat, the extension to $\Gamma(t)$ of its tangential and normal vectors are trivial. For the tangential part, we obtain as usual for Nitsche-based FSI

$$ \tau_w^T \sigma_s n - \tau_w^T \sigma_f n + \gamma_{\text{fsi}}(\dot{d} - u) \cdot \tau_w = 0 \quad \text{on } \Gamma(t). \qquad (26) $$

For the normal part, we have

$$ \sigma_{s,n} - \sigma_{f,n} + \gamma_{\text{fsi}}(\dot{d} - u) \cdot n_w + \gamma_C[P_\gamma(d, \lambda)]_+ = 0 \quad \text{on } \Gamma(t). $$

Let us first consider the case that the contact force is not active. We obtain, as in the standard FSI-Nitsche formulation

$$ \llbracket \sigma_n \rrbracket + \gamma_{\text{fsi}}(\dot{d} - u) \cdot n_w = 0. $$

If the contact force is active, we get

$$ \begin{aligned} 0 &= \llbracket \sigma_n \rrbracket + \gamma_{\text{fsi}}(\dot{d} - u) \cdot n_w + \gamma_C(d \cdot n_w - g_\alpha) - \lambda \\ &= \begin{cases} \gamma_C(d \cdot n_w - g_\alpha) + \gamma_{\text{fsi}}(\dot{d} - u) \cdot n_w, & \lambda = \llbracket \sigma_n \rrbracket \\ \gamma_C(d \cdot n_w - g_\alpha), & \lambda = \llbracket \tilde{\sigma}_n \rrbracket \end{cases} \end{aligned} \qquad (27) $$

In the first case, this is a combination of the "active" contact condition $d \cdot n_w = g$ and the continuity of velocities. As the continuity of velocities is imposed from the fluid side (25), this is not an issue for the model problem considered here. It will however lead to problems, when we consider contact of the solid with the lower wall $\Gamma_w$ in Section 3.2. There, the second formulation, based on the discrete stresses (21) will be needed. As can be seen in (27), the pure contact condition $d \cdot n_w = g$ is valid from the solid side during contact.

## 3.2 Contact with the boundary of the fluid domain

In this section, we will derive two numerical approaches for the full FSI-contact problem, considering contact with the lower fluid boundary $\Gamma_w$. First, we note that the formulation derived in the previous section has a simple extension to the case of contact with the fluid boundary, if we relax the contact condition by a small $\epsilon > 0$. see Section 3.2.1. Then, we introduce an artificial fluid below the contact line in Section 3.2.2, where we drive the fluid velocity to zero by means of a volume penalty approach. This enables us to include the FSI interface and contact conditions fully implicitly without relaxation.

### 3.2.1 Relaxed contact formulation

The idea of the relaxed formulation is to place the virtual obstacle of Section 3.1 at a distance $\alpha = \epsilon(h)$ from $\Gamma_w$, see Figure 2 on the left. We assume that $\epsilon(h) \to 0$, as the mesh size $h$ tends to zero. Now we can use the contact formulation derived for an obstacle within the fluid domain above using the function $g_{\epsilon(h)}(x_2) = g_0(x_2) - \epsilon(h)$. The corresponding variational formulation is Variational Formulation 1 using $g_{\epsilon(h)}$ in the definition of $P_\gamma$ (19) instead of $g_\alpha$.

Besides its simplicity, the main advantages of this contact formulation are:

* The numerical difficulties related with a topology change of the fluid domain are avoided.

* No-slip conditions can be used on both $\Gamma(t)$ and the lower wall $\Gamma_w$.
<figure><img src="image_5.png" /><figcaption>Figure 2: Illustration of the two possible approaches to include contact. Left: Relaxed contact formulation: Body in contact with an obstacle line $\Gamma_{\epsilon}$ close to $\Gamma_{w}$. Right: Introduction of an artificial fluid domain $\Omega_{f}^{C}$ below $\Gamma_{w}$.</figcaption></figure>

On the other hand, from a modelling point of view the contact conditions (15) make only sense, if we assume that an infinitesimal fluid layer remains between the fluid-structure-interface $\Gamma(t)$ and the lower wall $\Gamma_{w}$. Then the fluid stresses $\sigma_{f}$ appearing in the jump terms $[\sigma n]$ remain well-defined during contact. This assumption might be justified, if no-slip conditions are used on interface $\Gamma(t)$ and boundary $\Gamma_{w}$, see the theoretical works of e.g. Gerart-Varet & Hillairet [22]. In this framework, the choice of $\lambda$ in terms of the physical or numerical stresses appears not to be essential.

### 3.2.2 Contact formulation using an artificial fluid

The idea of this second approach is to add an artificial fluid domain $\Omega_{f}^{C}$ below $\Gamma_{w}$, see Figure 2, right sketch. In $\Omega_{f}^{C}$, we will drive the velocity $u$ to zero by means of a penalty term.

The variational formulation reads:

**Variational Formulation 2.** Find $u \in \mathcal{V}_h, p \in \mathcal{Q}_h, d \in \mathcal{W}_h$ such that $\dot{d} = \partial_t d$ and

$$ 
\begin{aligned}
\mathcal{A}_{FSI,*}^{C}(u, p, d, \dot{d})(v, q, w) + \gamma_C ([P_{\gamma}(\lambda, d)]_+, w \cdot n_w)_{\Gamma(t)} \\
&= (f_f, v)_{\Omega_f(t)} + (f_s, w)_{\Omega_s(t)} \quad \forall v, q, w \in \mathcal{V}_h \times \mathcal{Q}_h \times \mathcal{W}_h,
\end{aligned}
\quad (28) $$

where in the no-slip case

$$ 
\begin{aligned}
\mathcal{A}_{\text{fsi, no-slip}}^{C}(u, p, d, \dot{d})(v, q, w) &:= (\partial_t u, v)_{\Omega_f(t) \cup \Omega_f^C} + (\sigma_f(u, p), \nabla v)_{\Omega_f(t) \cup \Omega_f^C} + (\text{div } u, q)_{\Omega_f(t) \cup \Omega_f^C} \\
&\quad + S(p, q) + \gamma_a(u, v)_{\Omega_f^C} + (\partial_t \dot{d}, w)_{\Omega_s(t)} + (\sigma_s(d), \nabla w)_{\Omega_s(t)} \\
&\quad - (T_f(u, p, \dot{d}), w - v)_{\Gamma(t)} - (\dot{d} - u, \sigma_f(v, -q)n)_{\Gamma(t)}
\end{aligned}
\quad (29) $$

and for slip interface conditions

$$ 
\begin{aligned}
\mathcal{A}_{\text{fsi, slip}}^{C}(u, p, d, \dot{d})(v, q, w) &:= (\partial_t u, v)_{\Omega_f(t) \cup \Omega_f^C} + (\sigma_f(u, p), \nabla v)_{\Omega_f(t) \cup \Omega_f^C} + (\text{div } u, q)_{\Omega_f(t) \cup \Omega_f^C} \\
&\quad + S(p, q) + \gamma_a(u, v)_{\Omega_f^C} \\
&\quad + (\partial_t \dot{d}, w)_{\Omega_s(t)} + (\sigma_s(d), \nabla w)_{\Omega_s(t)} - (T_f(u, p, \dot{d}) \cdot n, (w-v) \cdot n)_{\Gamma(t)} \\
&\quad - ((\dot{d} - u) \cdot n, n^T \sigma_f(v, -q)n)_{\Gamma(t)}.
\end{aligned}
\quad (30) $$

Note the presence of the penalty term $\gamma_a(u, v)_{\Omega_f^C}$ within the artificial fluid, where $\gamma_a := \gamma_a^0 h^{-2}$.
**Discussion of the weakly imposed interface conditions: The no-slip case** We have already derived the weakly imposed interface conditions for the no-slip case in Section 3.1.1. Here, however, we have to consider that the fluid below the contact line is artificial, and we should in particular make sure that there is no feedback from the artificial fluid to the solid. Owing to (27), this naturally motivates the choice of $\lambda$ in terms of the numerical stress.

First, we note that the continuity of velocities (25) is imposed from the fluid side, such that no feedback to the solid is included from this equation. Considering the contact condition for the normal contact (27), we obtain

$$ \gamma_C(d \cdot n_w - g_0) + \gamma_{\text{fsi}}(\dot{d} - u) \cdot n_w = 0, \quad (31) $$

when choosing $\lambda = \llbracket \sigma_n \rrbracket$. Instead of the condition $d \cdot n_w = g$, this induces an influence from the artificial velocity $u$ from $\Omega_f^C$ onto the solid displacement. Moreover, if $u$ is driven to zero in $\Omega_f^C$, $\dot{d} \cdot n_w$ goes to zero as well, which might prevent the body from releasing from contact. On the other hand, using the jump of fluxes $\lambda = \llbracket \tilde{\sigma}_n \rrbracket$, we obtain the "pure" contact condition

$$ d \cdot n_w = g_0 \quad (32) $$

as desired.

The weakly imposed tangential contact condition (26) reads

$$ \tau_w^T \sigma_s n + \gamma_{\text{fsi}} \dot{d} \cdot \tau_w = 0, \quad (33) $$

when considering that the fluid velocity $u$ is driven to zero asymptotically. As $\gamma_{\text{fsi}} \to \infty$ for $h \to 0$, this means (asymptotically) that the solid is not allowed to slide along the line $\Gamma_w$. While this might seem restrictive at first sight, this condition is in fact in some sense inherited from the no-slip condition at $\Gamma(t)$ before contact. This is due to the continuity of velocities on $\Gamma_{\text{fsi}}(t)$ and the fact that the velocity is driven to zero in $\Omega_f^C$ (and hence on $\Gamma_w$). Moreover, the no-slip condition on the fluid part of $\Gamma_w$ implies that the solid can not slide on the end points $x_{C,1}$ and $x_{C,2}$ of the contact interval (see Figure 2, right sketch). Altogether, this shows in agreement with a number of theoretical works (e.g. Gerard-Varet et al [21]) that the no-slip interface conditions are not an appropriate model for the case that it comes to contact with an exterior wall.

## 3.3 Slip conditions

Motivated by the above considerations, we study slip conditions in this section. The arguments that follow can be applied to the formulations proposed in Sections 3.1 and 3.2. As mentioned in the introduction, Gerard-Varet, Hillairet and Wang showed that contact can not happen for a smooth rigid body falling down to the ground by means of gravity, if the no-slip boundary conditions are used for the velocities on both the fluid-structure interface $\Gamma_{\text{fsi}}(t)$ and the contact wall $\Gamma_w$, but that contact happens in this example, if slip boundary conditions are used on both $\Gamma_{\text{fsi}}(t)$ and $\Gamma_w$ [21].

Observe that in the continuous case the discussion of (14)-(16) remains valid in the case of contact with slip conditions (noting that the tangential stresses vanish on both sides of the fluid-solid interfaces). However, at the discrete level, care has to be taken to use only the normal component of the numerical stress. Indeed, in this case the relation (16) becomes

$$ (n^T \llbracket \sigma n \rrbracket) n - \lambda n_w = 0 \quad \text{on } \Gamma(t), \quad (34) $$

which yields $\lambda = (n^T \llbracket \sigma n \rrbracket)(n \cdot n_w) =: \llbracket \sigma_{n, \text{slip}} \rrbracket$.
At the discrete level, using the numerical stress this expression translates to

$$ \lambda = (n^T[\tilde{\sigma}n])(n \cdot n_w) =: [\tilde{\sigma}_{n,slip}], \quad (35) $$

where the jump operator is given by (22).

The resulting discrete formulations are given by Variational Formulation 1 or Variational Formulation 2, respectively, with the respective choice of $\lambda$.

**Weakly imposed interface conditions: The slip case** Let us consider again which interface conditions are implicitly included in the variational formulation, when using the numerical stresses (35) for $\lambda$. Considering the interface terms with fluid test function $v$ yields as usual the interface conditions $\tau^T \sigma_f n = 0$ and $(\dot{d} - u) \cdot n = 0$. Let us therefore concentrate on the terms with the solid test function $w$

$$ (\sigma_s n, w)_{\Gamma(t)} - (T_f \cdot n, w \cdot n)_{\Gamma(t)} + \gamma_C([\tilde{P}_\gamma]_+, w \cdot n_w)_{\Gamma(t)} = 0. $$

In the case without contact, the last term vanishes and hence we retrieve the standard consistency of Nitsche's method for fluid-structure interaction with slip conditions. On the other hand, by developing the solid test function $w$ in the local basis of $\Gamma(t)$ we have

$$ (\tau^T \sigma_s n, w \cdot \tau)_{\Gamma(t)} + ((n^T[\tilde{\sigma}n]), w \cdot n)_{\Gamma(t)} + \gamma_C([P_\gamma(d, \lambda)]_+, w \cdot n_w)_{\Gamma(t)} = 0. $$

Now, we use the identity $n = (n \cdot n_w)n_w + (n \cdot \tau_w)\tau_w$ to get

$$ (\tau^T \sigma_s n, w \cdot \tau)_{\Gamma(t)} + ((n^T[\tilde{\sigma}n]), (w \cdot \tau_w)(n \cdot \tau_w))_{\Gamma(t)} + ((n^T[\tilde{\sigma}n]), (w \cdot n_w)(n \cdot n_w))_{\Gamma(t)} \\ + \gamma_C([P_\gamma(d, \lambda)]_+, w \cdot n_w)_{\Gamma(t)} = 0. $$

Hence, in the case of contact by using the definition of $P_\gamma$, we get

$$ (\tau^T \sigma_s n, w \cdot \tau)_{\Gamma(t)} + ((n^T[\tilde{\sigma}n]), (w \cdot \tau_w)(n \cdot \tau_w))_{\Gamma(t)} + \gamma_C(d \cdot n_w - g_0, w \cdot n_w)_{\Gamma(t)} = 0. \quad (36) $$

When it comes to contact, we have in the asymptotic limit $\tau_w = \tau$ and $n \cdot \tau_w = 0$, so that (36) reduces to the “no-friction” condition $\tau^T \sigma_s n = 0$ and the non-penetration condition $d \cdot n_w = g_0$. This is the desired contact condition, as in the no-slip case, see (32).

We observe that in (36) both the tangential and the normal components (with respect to $\tau_w$ and $n_w$ have asymptotically vanishing perturbations. It is possible to eliminate the perturbation in the normal component, from the solid stress term $(\tau^T \sigma_s n, w \cdot \tau)_{\Gamma(t)}$, by adding the corresponding term to the definition of $\lambda$, i.e.

$$ \lambda = [\tilde{\sigma}_{n,slip}] + (\tau^T \sigma_s n)(\tau \cdot n_w). \quad (37) $$

In this case the relation (36) takes the form

$$ (\tau^T \sigma_s n, w \cdot \tau_w(\tau \cdot \tau_w))_{\Gamma(t)} + ((n^T[\tilde{\sigma}n]), (w \cdot \tau_w)(n \cdot \tau_w))_{\Gamma(t)} + \gamma_C(d \cdot n_w - g_0, w \cdot n_w)_{\Gamma(t)} = 0. $$

Testing with $w = n_w$ we see that here the non-penetration condition $d \cdot n_w = g_0$ is imposed unperturbed for all $h > 0$. By moving the perturbation to $\lambda$ as in (37) it is instead the form $P_\gamma$ that is perturbed, which implies a (weakly consistent) perturbation of the contact zone.
# 4 Stability

In this work, we will use fitted finite elements, i.e. we assume that both the interface $\Gamma(t)$ and the wall $\Gamma_w$ are resolved by mesh lines. The approaches presented, however, can be combined with unfitted finite elements as well, with the only difference that some more stabilisation terms have to be added to the variational formulation. In order to simplify the presentation, we will concentrate on the fitted case first and discuss the extension to unfitted finite elements in a remark afterwards.

We will use equal-order finite elements in combination with a pressure stabilisation term $S_p(p, q)$ for the fluid equations. For the stability analysis, the only requirement on $S_p$ is that it leads to a well-posed discrete fluid problem. Possibilities include the Brezzi-Pitkäranta stabilisation [50], local projections (LPS) [51], the pressure-stabilised Petrov-Galerkin approach (PSPG) [52] or the Continuous Interior Penalty method (CIP) [53]).

In order to present the stability analysis in a general setting, we introduce some further notation splitting the contact force variable $\lambda$ into a fluid part $\lambda_f$ and a solid part $\lambda_s$, such that $\lambda = \lambda_s - \lambda_f$. For the no-slip case, we define $\lambda_s(w) = \sigma_{s,n}(w)$ and

$$ \lambda_f(v, q, w) = \begin{cases} \sigma_{f,n}(v, q), & \text{if } \lambda = \llbracket \sigma_n \rrbracket, \\ T_f(v, q, w) \cdot n, & \text{if } \lambda = \llbracket \tilde{\sigma}_n \rrbracket, \end{cases} \quad (38) $$

see (17) and (21). For the slip case, we define $\lambda_s(w) = \sigma_{s,n}(w)(n \cdot n_w)$ and

$$ \lambda_f(v, q, w) = \begin{cases} \sigma_{f,n}(v, q)(n \cdot n_w), & \text{if } \lambda = \llbracket \sigma_{n, \text{slip}} \rrbracket, \\ (T_f(v, q, w) \cdot n)(n \cdot n_w), & \text{if } \lambda = \llbracket \tilde{\sigma}_{n, \text{slip}} \rrbracket, \end{cases} \quad (39) $$

see (34) and (35).

## 4.1 Generalised contact formulation

Before we conduct the stability analysis, let us introduce a generalised contact formulation, following the ideas of Chouly et al. [11]. We have already briefly discussed their ideas for the case of a pure solid with contact in (13). The generalisation of (13) to the FSI-contact system (Variational Formulation 2) would be to add the terms

$$ -(\gamma_C[P_\gamma(\lambda, d)]_+ + \lambda, \lambda(\partial_t v, \partial_t q, w))_{\Gamma(t)}. \quad (40) $$

The time derivatives on the test functions $v$ and $q$ are motivated by the stability analysis below, where we have to test the variational form with $v = u, q = p$ and $w = \partial_t d$, in order to show stability (see also Burman & Fernández [38]).

On the other hand, the term (40) is not usable within a time-stepping scheme due to the time derivatives on the test functions. A remedy is to shift the time derivatives to the first integrand (ignoring the boundary terms), i.e. adding the consistent terms

$$ -(\gamma_C[P_\gamma(\lambda, d)]_+ + \lambda, \lambda_s(w))_{\Gamma(t)} - (\partial_t(\gamma_C[P_\gamma(\lambda, d)]_+ + \lambda), \lambda_f(v, q, w))_{\Gamma(t)}. $$

This yields the variational formulation:

**Variational Formulation 3.** Find $u(t) \in V_h, p(t) \in Q_h, d(t) \in W_h$, such that $\dot{d} = \partial_t d$ and

$$ \begin{aligned} \mathcal{A}(u, p, d)(v, q, w) &:= \mathcal{A}_{*,FSI}^*(u, p, d)(v, q, w) + \gamma_C([P_\gamma(\lambda, d)]_+, w \cdot n_w)_{\Gamma(t)} \\ & \quad -\theta(\gamma_C[P_\gamma(\lambda, d)]_+ + \lambda, \lambda_s(w))_{\Gamma(t)} \\ & \quad -\theta(\partial_t(\gamma_C[P_\gamma(\lambda, d)]_+ + \lambda), \lambda_f(v, q, w))_{\Gamma(t)} \\ & \quad = (f_f, v)_{\Omega_f(t)} + (f_s, w)_{\Omega_s(t)} \quad \forall v, q, w \in V_h \times Q_h \times W_h, \end{aligned} \quad (41) $$
where $\mathcal{A}_{*,FSI}^*$ is one of the bilinear forms $\mathcal{A}_{\text{no-slip,FSI}}$, $\mathcal{A}_{\text{slip,FSI}}$, $\mathcal{A}_{\text{no-slip,FSI}}^C$ or $\mathcal{A}_{\text{slip,FSI}}^C$.

## 4.2 A stability result

In this section, we will investigate stability of the discrete formulation for different values of $\theta$. In particular, we will show a stability result for the symmetric formulation ($\theta = 1$) and stability up to a specific term in the general case (including $\theta = 0$). These results correspond to the results that have been obtained by Chouly et al. [11] for the pure solid case. For the stability analysis, we will assume infinitesimal displacements, i.e. the sub-domains $\Omega_f$ and $\Omega_s$ as well as the interface $\Gamma$ are fixed. We introduce the notation $\tilde{\Omega}_f$ for the combined fluid and artificial fluid domain in the case of $\mathcal{A}_{*,FSI}^C$ and for $\Omega_f$ else.

**Theorem 4.1.** Let $u, p, d \in V_h$ and $\dot{d} = \partial_t d$. We have the following stability result for the form $\mathcal{A}$ defined in (41), where $\theta \in [0, 1]$ and $\gamma_C^0$ sufficiently large

$$ \begin{aligned} & \|u(T)\|_{\tilde{\Omega}_f}^2 + \|\dot{d}(T)\|_{\Omega_s}^2 + \|d(T)\|_{H^1(\Omega_s)}^2 + \int_0^T \left( \nu_f \|\nabla u\|_{\tilde{\Omega}_f}^2 + S_p(p, p) + \gamma_a \|u\|_{\Omega_f}^2 \right) dt \\ & \qquad + \int_0^T \gamma_{fsi} \|(\dot{d} - u) \cdot n\|_{\Gamma}^2 dt + \theta \left\| \gamma_C^{-1/2} \lambda(T) + \gamma_C^{1/2} [P_{\gamma}(\lambda, d)]_+ (T) \right\|_{\Gamma}^2 \\ & \le C \left( \int_0^T \mathcal{A}(u, p, d; u, p, \dot{d}) - (1-\theta)\gamma_C ([P_{\gamma}(\lambda, d)]_+, \dot{d} \cdot n_w)_\Gamma dt + \|u_0\|_{\tilde{\Omega}_f}^2 \right. \\ & \qquad \left. + \|\dot{d}_0\|_{\Omega_s}^2 + \|d_0\|_{H^1(\Omega_s)}^2 + \theta \left\| \gamma_C^{-1/2} \lambda_0 + \gamma_C^{1/2} [P_{\gamma}(\lambda_0, d_0)]_+ \right\|_{\Gamma}^2 \right), \end{aligned} $$

where we have used the abbreviation $\lambda_0 := \lambda(u_0, p_0, d_0)$. For the no-slip case, the term $\|(\dot{d} - u) \cdot n\|_{\Gamma}$ can be replaced by $\|(\dot{d} - u)\|_{\Gamma}$.

**Remark 4.1.** (Contact terms) The second line gives us control over the satisfaction of the FSI-contact condition for $\theta > 0$. In contrast to the work by Chouly et al. [11] for a pure solid problem, here we obtain discrete stability for $\theta = 1$, for the following positive discrete energy

$$ E(T) := \|u(T)\|_{\tilde{\Omega}_f}^2 + \|\dot{d}(T)\|_{\Omega_s}^2 + \|d(T)\|_{H^1(\Omega_s)}^2 + \left\| \gamma_C^{-1/2} \lambda(T) + \gamma_C^{1/2} [P_{\gamma}(\lambda, d)]_+ (T) \right\|_{\Gamma}^2. $$

For $\theta \neq 1$ on the other hand, the contact term $(1-\theta)\gamma_C ([P_{\gamma}(\lambda, d)]_+, \dot{d} \cdot n_w)_\Gamma$ appears on the right-hand side. The last term on the right-hand side vanishes, if we assume that the contact conditions are fulfilled at initial time, for example if the solid is not in contact with at $t = 0$.

**Proof.** We test (41) with $w = \dot{d} = \partial_t d, v = u$ and $q = p$ and integrate in time. We start by deriving a lower bound for $\mathcal{A}_{*,FSI}^*$. For the fluid part, we use the techniques from Burman & Fernandez [38], to show coercivity of the Stokes part including the coupling terms. For the no-slip case, the authors have shown

$$ \begin{aligned} & (\sigma_f(u, p), \nabla u)_{\tilde{\Omega}_f} + (\operatorname{div} u, p)_{\tilde{\Omega}_f} + S_p(p, p) + \gamma_a(u, u)_{\Omega_f^C} \\ & \qquad - (\sigma_f(u, p)n, \dot{d} - u)_\Gamma - (\dot{d} - u, \sigma_f(u, -p)n)_\Gamma + \gamma_{fsi} \|\dot{d} - u\|_{\Gamma}^2 \\ & \ge c \left( \nu_f \|\nabla u\|_{\tilde{\Omega}_f}^2 + \gamma_{fsi} \|\dot{d} - u\|_{\Gamma}^2 + S_p(p, p) + \gamma_a \|u\|_{\Omega_f^C}^2 \right). \end{aligned} $$
Analogously, one can show in the slip-case that

$$ (\sigma_f(u,p), \nabla u)_{\bar{\Omega}_f} + (\operatorname{div} u, p)_{\bar{\Omega}_f} + S_p(p,p) + \gamma_a(u, u)_{\Omega_f^c} \\ - (n^T \sigma_f(u,p)n, (\dot{d}-u) \cdot n)_\Gamma - ((\dot{d}-u) \cdot n, n^T \sigma_f(u,-p)n)_\Gamma + \gamma_{\text{fsi}} \|(\dot{d}-u) \cdot n\|_{\Gamma}^2 \\ \geq c \left( \nu_f \|\nabla u\|_{\bar{\Omega}_f}^2 + \gamma_{\text{fsi}} \|(\dot{d}-u) \cdot n\|_{\Gamma}^2 + S_p(p,p) + \gamma_a \|u\|_{\Omega_f^c}^2 \right). $$

Using the symmetry of $\sigma_s$, integration in time and a Korn's inequality, we obtain for the solid part

$$ \int_0^T (\sigma_s(d), \nabla \dot{d})_{\Omega_s} dt = \frac{1}{2} \int_0^T \partial_t (\sigma_s(d), \nabla d)_{\Omega_s} dt = \frac{1}{2} ((\sigma_s(d(T)), \nabla d(T))_{\Omega_s} - (\sigma_s(d(0)), \nabla d(0))_{\Omega_s}) \\ \geq c_1 \|\nabla d(T)\|_{\Omega_s}^2 - c_2 \|\nabla d_0\|_{\Omega_s}^2. $$

Moreover, we have

$$ \int_0^T (\partial_t u, u)_{\bar{\Omega}_f} + (\partial_t \dot{d}, \dot{d})_{\Omega_s} dt = \frac{1}{2} \left( \|u(T)\|_{\bar{\Omega}_f}^2 + \|\dot{d}(T)\|_{\Omega_s}^2 - \|u_0\|_{\bar{\Omega}_f}^2 - \|\dot{d}_0\|_{\Omega_s}^2 \right). $$

Together, we have shown that

$$ \|u(T)\|_{\bar{\Omega}_f}^2 + \|\dot{d}(T)\|_{\Omega_s}^2 + \|d(T)\|_{H^1(\Omega_s)}^2 + \int_0^T \nu_f \|\nabla u\|_{\bar{\Omega}_f}^2 + S_p(p,p) + \gamma_a \|u\|_{\Omega_f^c}^2 + \gamma_{\text{fsi}} \|\dot{d}-u\|_{\Gamma}^2 dt \\ \leq C \left( \int_0^T A_{*,FSI}^*(u,p,d;u,p,\dot{d}) dt + \|u_0\|_{\bar{\Omega}_f}^2 + \|\dot{d}_0\|_{\Omega_s}^2 + \|d_0\|_{H^1(\Omega_s)}^2 \right) $$

Let us now estimate the contact terms. We split the principal contact term into

$$ \int_0^T \gamma_C ([P_\gamma(\lambda, d)]_+, \dot{d} \cdot n_w)_\Gamma dt = \int_0^T \theta \gamma_C ([P_\gamma(\lambda, d)]_+, \dot{d} \cdot n_w)_\Gamma + (1-\theta) \gamma_C ([P_\gamma(\lambda, d)]_+, \dot{d} \cdot n_w)_\Gamma dt. $$

We have to estimate the terms

$$ \theta \int_0^T \underbrace{\gamma_C ([P_\gamma(\lambda, d)]_+, \dot{d} \cdot n_w)}_{I_1} \Gamma - \underbrace{\gamma_C [P_\gamma(\lambda, d)]_+ + \lambda, \lambda_s(\dot{d})}_ {I_2} \Gamma \\ - \underbrace{(\partial_t (\gamma_C [P_\gamma(\lambda, d)]_+ + \lambda), \lambda_f(u, p, \dot{d}))}_{I_3} \Gamma dt \quad (42) $$

From the definition of $P_\gamma$ we can write $d \cdot n_w = P_\gamma(\lambda, d) + \gamma_C^{-1}\lambda$. Hence, since the lower wall is assumed to be time independent, we have

$$ \int_0^T I_1 dt = \gamma_C \int_0^T ([P_\gamma(\lambda, d)]_+, \partial_t(P_\gamma(\lambda, d) + \gamma_C^{-1}\lambda))_\Gamma dt \\ = \gamma_C \int_0^T ([P_\gamma(\lambda, d)]_+, \partial_t([P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda))_\Gamma dt. $$
In the second line, we have used that [11]

$$ \frac{1}{2} \partial_t [\phi]_+^2 = [\phi]_+ \partial_t [\phi]_+ = [\phi]_+ H(\phi) \partial_t [\phi]_+ = [\phi]_+ \partial_t \phi $$

where $H$ is the Heaviside function. We insert $\pm\gamma_C^{-1}\lambda$ and integrate by parts

$$ \begin{aligned} \int_0^T I_1 dt &= \gamma_C \int_0^T ([P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda, \partial_t([P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda))_\Gamma dt \\ &\quad - \int_0^T \underbrace{(\lambda, \partial_t([P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda))}_\Gamma dt \\ &= \underbrace{\frac{\gamma_C}{2} \| [P_\gamma(\lambda, d)(T)]_+ + \gamma_C^{-1}\lambda(T) \|_{0,\Gamma}^2}_{D_T} - \underbrace{\frac{\gamma_C}{2} \| [P_\gamma(\lambda_0, d_0)]_+ + \gamma_C^{-1}\lambda_0 \|_{0,\Gamma}^2}_{D_0} - \int_0^T I_4 dt. \end{aligned} $$

Now the idea is to split the contribution from the term $I_4$ into fluid and solid stresses and to apply integrating by parts in time (only) in the solid stress contribution. By definition, we have

$$ - \int_0^T I_4 dt = - \int_0^T (\lambda_s(d), \partial_t([P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda))_\Gamma + (\lambda_f(u, p, \dot{d}), \partial_t([P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda))_\Gamma dt $$

and, by integrating by parts in the first term of the right-hand side, we have

$$ - \int_0^T I_4 dt = \underbrace{\int_0^T (\lambda_s(\dot{d}), [P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda)_\Gamma dt}_{I_5} - \underbrace{\int_0^T (\lambda_s(d(T)), [P_\gamma(\lambda, d)(T)]_+ + \gamma_C^{-1}\lambda(T))_\Gamma dt}_{I_6} \\ + \underbrace{\int_0^T (\lambda_s(d_0), [P_\gamma(\lambda_0, d_0)]_+ + \gamma_C^{-1}\lambda_0)_\Gamma dt}_{I_7} + \underbrace{\int_0^T (\lambda_f(u, p, \dot{d}), \partial_t([P_\gamma(\lambda, d)]_+ + \gamma_C^{-1}\lambda))_\Gamma dt}_{I_8} $$

The terms $I_5$ and $I_8$ cancel with the terms $I_2$ and $I_3$ in (42). The term $I_6$ is treated in a standard fashion using Young's inequality, an inverse inequality and the dissipation provided by $D_T$ and the elastic energy $\frac{1}{2}a^s(d(T), d(T))$ for $\gamma_C^0$ sufficiently large:

$$ \begin{aligned} I_6 &\geq -\gamma_C^{-1} \|\lambda_s(d(T))\|_{0,\Gamma}^2 - \frac{\gamma_C}{4} \|[P_\gamma(\lambda(T), d(T))]_+ + \gamma_C^{-1}\lambda(T)\|_{0,\Gamma}^2 \\ &\geq -\frac{1}{4} \|d(T)\|_{H^1(\Omega_s)}^2 - \frac{1}{2} D_T. \end{aligned} $$

In the same way, we obtain for $I_7$

$$ I_7 \leq \frac{1}{4} \|d_0\|_{H^1(\Omega_s)}^2 + \frac{1}{2} D_0, $$

which completes the proof. $\square$

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

**Remark 4.2.** (Unfitted finite elements) When using unfitted finite elements [54, 38], additional stabilisation terms $S_u$ and $S_d$ are needed, if the interface $\Gamma(t)$ is not resolved by mesh lines. Their purpose is to extend the coercivity of the fluid system from $\Omega_f$ (resp. $\Omega_s$) to the extended domains $\Omega_h^f$ (resp. $\Omega_h^s$) that consists of all element $T \in \mathcal{T}_h$, with a non-empty intersection with the respective sub-domain $(T \cap \Omega_i \neq \emptyset)$. Suitable “ghost penalty” operators have been defined in Burman [47]. The same stability result as in Theorem 4.1 can then be shown with an analogous argumentation.
**Remark 4.3.** (Newton convergence) While the symmetric formulation ($\theta = 1$) seems beneficial from the theoretical point of view, the additional terms in (41) can cause severe difficulties for the non-linear solver. The reason is that the additional contact terms are not only highly non-linear, but also non-smooth, especially due to the time derivative acting on the maximum operator $\partial_t[\cdot]_+$. In our numerical tests, we were not able to obtain numerical convergence for different versions of generalised Newton methods [55]. The investigation of the case $\theta \neq 1$ by means of numerical tests and in particular the construction of a robust non-linear solver are subject to future research.

# 5 Numerical Results

In this section, we show some numerical results to analyse and to compare the different contact formulations. As mentioned in Remark 4.3, we were not able to obtain results for the generalised contact formulation with $\theta \neq 0$, due to divergence of the generalised Newton-type methods we have tried. Therefore, we only show results for $\theta = 0$, where Newton convergence was not an issue, at least when the time step $\delta t$ was chosen reasonably small. In the computations made for this paper, the Newton algorithm needed 1-2 iterations per time step to reduce the initial residual by a factor of $10^{-7}$, if the contact force was not getting active during the iteration, and 1-5 iterations per time step in and around the interval of contact.

We first give some details in Section 5.1 on the fitted, equal-order finite element discretisation and the stabilisations we use. Then, in Section 5.2, we study the problem of a virtual obstacle within the fluid domain introduced in Section 3.1. The purpose of this example is to isolate the effect of the contact terms from issues related to discretisation during contact and the topology change in the fluid domain $\Omega_f(t)$. Then, we study in Section 5.3 a model problem with contact with the boundary of the fluid domain, where we compare among other aspects the two contact formulations introduced in Section 3, the different possibilities to choose the fluxes $\lambda$ and the effect of slip and no-slip boundary and interface conditions.

## 5.1 Details on discretisation and stabilisation

For the numerical results in this paper, we will use a monolithic *Fully Eulerian* approach on a global mesh $\mathcal{T}_h$ covering $\Omega(t)$. In order to resolve the interface $\Gamma(t)$ within the discretisation, we use the *locally modified finite element method* introduced by Frei & Richter [45]. The idea of this approach is to use a fixed coarse triangulation $\mathcal{T}_{2h}$ of the overall domain $\Omega = \tilde{\Omega}_f(t) \cup \Gamma(t) \cup \Omega_s(t)$ that is independent of the position of the interface $\Gamma(t)$. Then, in each time step, this coarse grid is refined once by splitting each so-called “patch” element in either eight triangles or four quadrilaterals to resolve the interface in at least a linear approximation, see Figure 3 for an illustration.

The finite element space $V_h$ is then defined as a combination of piece-wise linear and piece-wise bi-linear finite elements on the patches. It can be guaranteed that a maximum angle condition is fulfilled in each of the sub-cells, leading to optimal-order interpolation and error estimates [45].

For temporal discretisation, we split the time interval $I$ into $m$ equidistant-distant time intervals $I_j = (t_{j-1}, t_{j-1} + k]$ and use a time-stepping scheme that is based on a modified discontinuous Galerkin time discretisation of lowest order (dG(0)), see Frei & Richter [46]. The displacement-velocity relation $\partial_t d = \dot{d}$ is included by means of the $L^2$-projection

$$ (\partial_t d, z)_{\Omega_s(t)} - (\dot{d}, z)_{\Omega_s(t)} = 0 \quad \forall z \in V_h^s, $$

where $V_h^s$ denotes the (modified) finite element space that is spanned by the degrees of freedom of the elements in the solid part $\Omega_h^s$.
<figure><img src="image_7.png" /><figcaption>Figure 3: Left: Fixed triangulation $\mathcal{T}_{2h}$ of the domain $\Omega$. Right: Subdivision of the patches $P \in \mathcal{T}_{2h}$ such that the interface $\Gamma(t)$ is resolved in a linear approximation by the discrete interface $\Gamma_h$.</figcaption></figure>

The domain affiliation of a point $x \in \Omega(t)$ is determined by means of the *Initial Point Set/Backward Characteristics method* [29, 30], that uses the displacement $d(t)$ in the solid domain and an extension to $\Omega_f(t)$ in order to trace back points to their initial position in $\Omega(0)$, following the definition (2).

For pressure stabilisation, we use an anisotropic variant of the Continuous Interior Penalty method, see Frei [56, 49]. In addition, we add the temporal pressure stabilisation term

$$S_{pt}(p, q) = \gamma_{pt} h(p^m - p^{m-1}, q)_{\Gamma(t)}$$

in each time interval $I_m$. This additional stabilisation is needed, as the mesh $\mathcal{T}_h(t_m)$, and hence the finite element spaces, change from time-step to time-step. The solution $u^{m-1}$ from the previous time-step $t_{m-1}$ is therefore not discrete divergence-free with respect to the new mesh $\mathcal{T}_h(t_m)$, which gives rise to pressure oscillations, see for example Besier & Wollner [57].

All the following results have been obtained using the finite element library *Gascoigne 3d* [58].

## 5.2 Virtual obstacle within the fluid domain

In order to isolate the effect of the contact model from issues related to the Navier-Stokes-contact paradox and the topology changes in the fluid domain, let us first study the simplified problem introduced in Section 3.1. We define the initial fluid and solid domains as

$$\Omega_f(0) = (0, 1) \times (0, 0.5), \quad \Omega_s(0) = (0, 1) \times (0.5, 0.6)$$

and a lower-dimensional obstacle $\Gamma_w = (0, 1) \times (0.25, 1)$ within the fluid domain. We consider a moving interface $\Gamma(t)$, which is resolved using the locally modified finite element method. The sub-domains $\Omega_s(t)$ and $\Omega_f(t)$ and the interface $\Gamma(t)$ depend on the solid displacement $d(t)$, see (2).

The constraint for the solid displacement is given by

$$d \cdot n_w \le 0.25 \quad (=: g_\alpha). \qquad (43)$$

We use the elasticity parameters $\lambda_s = \mu_s = 2 \cdot 10^6$ and the fluid viscosity $\nu_f = 1$. The structure is pulled towards the bottom by fluid forces due to a prescribed pressure mean value
at the left and right boundary of the fluid domain

$$ \int_{\Gamma_{f,\text{left}}} p \, ds = \int_{\Gamma_{f,\text{left}}} \bar{P} \, ds, \quad \int_{\Gamma_{f,\text{right}}} p \, ds = \int_{\Gamma_{f,\text{right}}} \bar{P} \, ds $$

where $\bar{P} := 1.3 \cdot 10^5$. We consider the Variational Formulation 1 with no-slip conditions and $\lambda = [[\bar{\sigma}_n]]$, see (21), on a Cartesian mesh that consists of 5120 elements and with a small time step $\delta t = 10^{-5}$. The Nitsche constant at the FSI interface is chosen $\gamma_{\text{fsi}}^0 = 10^3$ and temporal pressure stabilisation with $\gamma_{pt} = 10^{-2}$ is used.

To analyse the results, we define the “minimal distance”

$$ d_{\min} := \min_{x \in \Gamma(t)} x_2 - 0.25 $$

of the interface $\Gamma(t)$ to $\Gamma_w$. To be precise the term “minimal distance” is only correct before contact, as $d_{\min}$ gets negative in case of an overlap. Moreover, we define the following functionals in order to analyse the pressure $p$ and the contact force at the interface $\Gamma(t)$

$$ J_p := \left| \int_{\Gamma(t)} p \, ds \right|, \quad J_{P_\gamma} := \gamma_C \int_{\Gamma(t)} [P_\gamma(\lambda, d)]_+ ds. $$

In Figure 4, we plot these three functionals over time for two contact periods and the contact parameters $\gamma_C^0 = 10, 10^2, 10^3$. In the top left plot, we observe that the solid is pulled down until it reaches $\Gamma_w$ at $t \approx 1.5 \cdot 10^{-3}$. After a short contact period, it is released again due to its elastic properties before it reaches the obstacle for a second time at $t \approx 5 \cdot 10^{-3}$.

The contact condition $d \cdot n_w \le 0.25$ is only significantly violated for the smallest contact parameter, where $d_{\min}$ reaches a minimum value of around $-8 \cdot 10^{-4}$, see the zoom-in around the contact interval on the top left. This value is more than an order of magnitude smaller than the mesh size in vertical direction $h \approx 1.4 \cdot 10^{-2}$. For the larger values of $\gamma_C^0$, the minimal value of $d_{\min}$ is even much closer to zero. On the other hand, we observe that even for the largest value $\gamma_C^0 = 10^3$, the contact condition is slightly relaxed, allowing for very small overlaps of solid and contact line.

In the second row of Figure 4, we observe that the pressure shows a peak at the beginning ($t \approx 1.5 \cdot 10^{-3}$ and $t \approx 5 \cdot 10^{-3}$) of the contact periods, followed by some small oscillations. The peak is caused by the fluid dynamics and will be discussed below. The oscillations get smoother for larger values of $\gamma_C^0$ and are barely visible for $\gamma_C^0 = 10^3$.

Similarly, the contact force $J_{P_\gamma}$ shows oscillations for $\gamma_C^0 = 10$ and a much smoother behaviour for $\gamma_C^0 \ge 10^2$. Note that this does not contradict the stability result in Theorem 4.1, where we have assumed that $\gamma_C^0$ is large enough. Besides these oscillations, the contact force is roughly of the same size for different $\gamma_C^0$.

**Investigation of the pressure peak** The pressure peak at the beginning of the contact interval can be explained as follows. As the fluid does not "see" the obstacle before reaching it, the solid is pulled down towards it without reducing its velocity. At the moment when the obstacle is reached, its vertical velocity $\dot{d} \cdot n_w$ has to decrease to zero in an instant. Due to the continuity of velocities, the same happens for the fluid velocity $u \cdot n_w$ at the interface, and due to the incompressibility constraint the velocity has to change globally in the fluid domain $\Omega_f(t)$. The pressure can be seen as a Lagrange multiplier and more specifically as sensitivity of a (jumping) energy functional with respect to the incompressibility constraint, which explains the peak.
<figure><img src="image_8.png" /><figcaption>Figure 4: Top row: 'Minimal distance' $d_{min}$ to $\Gamma_w$ over time with two contact periods. Left: Total time interval. Right: Zoom-in at the first contact interval. Bottom left: Pressure functional $J_p$ over time. Bottom right: Functional $J_{P_\gamma}$ measuring the contact force around the contact interval over time.</figcaption></figure>

To substantiate this explanation numerically, we add an artificial penalty for the velocity on the sub-domain $\Omega_f^0$ below the contact line

$$S_a(u, v) := \gamma_a^0 h^{-2}(u, v)_{\Omega_f^0}.$$

For $\gamma_a := \gamma_a^0 h^{-2} \to \infty$, the fluid velocity is driven to zero below the obstacle. As this is already the case before contact, no abrupt changes in the fluid velocity are expected at the moment of the impact. Note that the problem with $\gamma_a^0 > 0$ is purely artificial, as the pressure mean values are still applied on the whole fluid boundary, including the boundary of $\Omega_f^0$.

In Figure 5 we compare the minimal distance and pressure functionals for computations without penalty ($\gamma_a^0 = 0$) to results for $\gamma_a^0 = 10$. First, we note that the vertical displacement is significantly influenced by the penalty, which has to be expected as the fluid dynamics are altered. Contact happens later at $t \approx 0.0022$ with the artificial penalty. Moreover, we observe indeed that the initial pressure peak at the time of impact is significantly reduced for $\gamma_a^0 = 10$.
<figure><img src="image_9.png" /><figcaption>Figure 5: Minimal distance $d_{min}$ to the contact line $\Gamma_w$ (left) and mean pressure over the interface $J_p$ (right), plotted over time with an artificial penalty for the velocity below the virtual obstacle.</figcaption></figure>

## 5.3 Contact problem

Next, we study a problem, where it comes to real contact with the wall $\Gamma_w = \{(x, y) \in \Omega, y = 0.25\}$. At time $t = 0$, we define

$$ \Omega_f(0) := (0, 1) \times (0.25, 0.5), \quad \Gamma(0) := (0, 1) \times 0.5, \quad \Omega_s(0) := (0, 1) \times (0.5, 0.6). $$

Below $\Gamma_w$, we define a fixed artificial fluid domain $\Omega_f^C := (0, 1) \times (0, 0.25)$.

We apply again a pressure mean value $\bar{P}$ on the lateral boundaries $\Gamma_{f,\text{left}}$ and $\Gamma_{f,\text{right}}$ of the fluid domain $\Omega_f(t)$. As the size of $\Gamma_{\text{left}}$ and $\Gamma_{\text{right}}$ is smaller and the viscous fluid forces acting against the “real” contact are stronger than in the previous example, we have to set a larger pressure force $\bar{P} = 3 \cdot 10^5$ in order to obtain contact. On the other hand, the contact was never released again in our numerical experiments, when we used this constant boundary force for all times. Therefore, we decrease $\bar{P}$ linearly from $t = 10^{-3}$ on until it reaches zero at $t = 1.2 \cdot 10^{-3}$. In order to avoid the issues related to no-slip conditions and contact, we use slip-interface conditions first, i.e. the Variational Formulation 2 with $A_{\text{slip,FSI}}^C$. Unless stated differently, $\lambda$ is chosen as the jump of numerical stresses $[\tilde{\sigma}_{n,\text{slip}}]$. Moreover, we use again a Cartesian mesh that consists of 5120 elements, a time step $\delta t = 10^{-5}$ and temporal pressure stabilisation with $\gamma_{pt} = 10^{-2}$. Unless explicitly stated, the Nitsche parameters are chosen as $\gamma_{\text{fsi}}^0 = \gamma_{\text{C}}^0 = 10^3$ and the penalty in the artificial fluid as $\gamma_a^0 = 10^2$.

The results on a coarser mesh are illustrated in Figure 6 at four time instants. Contact happens after the pressure on the lateral boundaries is released, as the solid continues moving downwards for some time. During contact, there is a very small overlap of the solid with the artificial fluid $\Omega_f^C$. As the overlap is of order $10^{-5}$, it can barely be seen in the bottom left picture. Notice however the triangular cells in $\Omega_f^C$ that are used only, when a patch is cut by $\Gamma(t)$.

**Comparison of the two contact formulations** First, we compare the two contact strategies derived in Section 3.2, i.e. the relaxed contact formulation introduced in Section 3.2.1 with a small gap of size $\epsilon(h) = h/10$ between the solid and $\Gamma_w$ and the strategy using an artificial fluid derived in Section 3.2.2. A comparison of the results for $d_{min}$, the pressure norm $\|p\|_{L^2(\Omega_{f,\text{mid}}(t))}$, where

$$ \Omega_{f,\text{mid}}(t) := \{x \in \Omega_f(t), 0.4 \le x_1 \le 0.6\} $$
<figure><img src="image_10.png" /><figcaption>Figure 6: Illustration of the contact problem at four time instances on a coarse mesh: $t = 0$ (top left), $t = 1.2 \cdot 10^{-3}$ (top right), $t = 2 \cdot 10^{-3}$ (bottom left) and $t = 2.5 \cdot 10^{-3}$ (bottom right). The grey part corresponds to the structure $\Omega_s(t)$, the white part is the artificial fluid $\Omega_f^C$. In the fluid domain $\Omega_f(t)$, values of the pressure $p$ are visualised.</figcaption></figure>

denotes the central part of the fluid domain, and the contact force $J_{P,\gamma}$ are shown in Figure 7 on two different meshes with 5120 and 20480 elements, respectively.

First, we observe from the plots in the top row that the interface stays at a distance to $\Gamma_w$ of about $\epsilon \approx h/10 \approx 1.4 \cdot 10^{-3}$ on the coarser and $\epsilon \approx h/10 \approx 7 \cdot 10^{-4}$ on the finer mesh for the relaxed formulation. The much smaller overlap with $\Omega_f^C$ in the artificial fluid formulation is not visible, not even in the zoom-in on the right.

While the curves for $d_{\min}$ look similar in the global picture (left), the zoom-in shows significant differences already before the impact. The contact happens earlier for the artificial fluid formulation: on the coarser mesh at time $t_{C,a} = 1.87 \cdot 10^{-3}$ compared to $t_{C,r} = 2.02 \cdot 10^{-3}$ for the relaxed formulation. This deviation is already much smaller on the finer mesh, where $t_{C,t} - t_{C,a} = 2 \cdot 10^{-5}$. The reason for this deviation is that in the artificial fluid formulation the wall $\Gamma_w$ is only asymptotically for $\gamma_a \to \infty$ impermeable for the fluid. Therefore, in practice, the fluid forces acting against the contact, in particular the pressure $p$, are smaller for this formulation. As $\gamma_a = \gamma_a^0 h^{-2} \to \infty$ for $h \to 0$ the difference is significantly reduced on the finer mesh.

To substantiate this explanation, we plot the pressure norm $\|p\|_{L^2(\Omega_{f,\mathrm{mid}}(t))}$ over the central part of the fluid domain $\Omega_{f,\mathrm{mid}}(t)$ and the time period before the impact in the bottom left figure. The maximum value of the norm on the coarser mesh at time $t = 1.1 \cdot 10^{-3}$ is approximately 20.563 for the relaxed formulation and about 19.065 for the artificial fluid version. On the finer mesh, the pressure values are much closer. After that time the functional values decrease because the domain $\Omega_{f,\mathrm{mid}}(t)$ gets smaller.

In the next paragraph, we will study the performance of both contact formulations under
<figure><img src="image_11.png" /><figcaption>Figure 7: Comparison of the relaxed and the artificial fluid contact formulation. Top: Minimal distance $d_{min}$ to $\Gamma_w$. Right: Zoom-in around the contact interval. Bottom: Pressure norm $\|p\|_{L^2(\Omega_f, \text{mid})}$ over the central part of the fluid domain before contact over time.</figcaption></figure>

mesh refinement.

**Convergence under mesh refinement** We solve the same problem on three different meshes with 1.280, 5.120 and 20.480 mesh elements, where the finer meshes are obtained from the coarsest one by global mesh refinement. The plots of the functionals $d_{min}$ and $J_{P\gamma}$ as well as the functionals

$$J_{\text{contact}} := \frac{\|\gamma_C^{1/2}[P_{\gamma}(\lambda, d)]_+ + \gamma_C^{-1/2}\lambda\|_{\Gamma(t)}}{\|\lambda\|_{\Gamma}}, \quad J_{\text{vel,fsi}} := \frac{\|(\dot{d} - u) \cdot n\|_{\Gamma_{\text{fsi}}(t)}}{\|u \cdot n\|_{\Gamma_{\text{fsi}}} + \|\dot{d} \cdot n\|_{\Gamma_{\text{fsi}}}}$$

measuring the fulfilment of the contact condition and the continuity of velocities on the part $\Gamma_{\text{fsi}}(t)$ of $\Gamma(t)$ that is not in contact with $\Gamma_w$

$$\Gamma_{\text{fsi}}(t) := \{x \in \Gamma(t) \mid P_{\gamma}(\lambda, d)(x) \le 0\}.$$

are shown in Figure 8 for the artificial fluid formulation and in Figure 9 for the relaxed contact formulation over time. The quantities $\|\cdot\|_*$ that are used to scale the functionals are temporal averages of the respective norms over the interval $I = [0, 0.004]$, computed on the finest grid.

First, we observe for both formulations in the plots on the top left that the contact happens later, the finer the discretisation is, as the fluid forces which act against the closure of the fluid
channel are better resolved on the finer meshes (see also Figure 7 and the related discussion above).

The curves for the contact force $P_γ$ on the top right of both figures show significant differences between the two formulations. While the functional values seem to converge for the artificial fluid formulation (if we neglect the time shift), the contact force gets larger under mesh refinement for the relaxed formulation. The larger values for the relaxed formulation are due to the presence of the fluid forces $\sigma_{f,n}$ during the whole contact interval, that are not penalised in this formulation. As a Lagrange multiplier for the incompressibility constraint, the continuous pressure $p$ gets singular when it comes to contact. The discrete pressure $p_h$ gets larger and larger under mesh refinement in our computations.

On the other hand, the functional $J_{contact}$, that measures the difference between $-\gamma_C^{1/2}[P_γ]_+$ and $\gamma_C^{-1/2}\lambda$ decreases under mesh refinement for both formulations. Besides the differences in the contact force $J_{P_γ}$, the functional values on each of the mesh levels are actually very similar for the two formulations. The reason must be that the fluid forces $\sigma_{f,n}$ enter in both $\lambda$ and $[P_γ]_+$. We conclude that the increase in the functional $J_{P_γ}$ seems not to be an issue for the contact dynamics.

Both $J_{contact}$ and the functional $J_{vel,fsi}$ are controlled by the stability estimate in Theorem 4.1 for $\theta = 1$. Although the parameter $\theta = 0$ is used here, we observe that both functionals decrease with mesh refinement before and during contact. While the convergence for the contact functional is quite slow, the values of the velocity functional indicate a convergence order $O(h^\alpha)$ with $0.5 \le \alpha \le 1$ for both formulations. Note that in contrast to the term $J_{contact}$, $J_{vel,fsi}$ is controlled in the stability estimate in Theorem 4.1 even with the pre-factor $(\gamma_{fsi}^0\mu_f)^{1/2}h^{-1/2}$.

**Flux formulations** Next, we compare the different choices for $\lambda$. We show results exemplarily for the artificial fluid formulation with slip interface conditions. We will compare results using the jump of stresses $\lambda = [[\sigma_{n,slip}]]$ (34), the jump of the numerical fluxes $\lambda = [[\tilde{\sigma}_{n,slip}]]$ (35) and the extended fluxes $\lambda = [[\tilde{\sigma}_{n,slip}]] + \tau^T\sigma_s n(\tau \cdot n_w)$ (see (37)). As the results for the latter two choices are nearly identical in this example, we show here only plots for the jump of stresses and the (non-extended) numerical fluxes. We use the artificial fluid formulation (Variational Formulation 2) and the previously used mesh with 5120 elements.

In Figure 10, we show the minimal distance $d_{min}$ to $\Gamma_w$, the contact force $J_{P_γ}$ and the integral over the velocity difference across the contact part $\Gamma_C(t)$ of the interface over time

$$J_{vel,C} := \int_{\Gamma_C(t)} (\dot{d} - u) \cdot n \, ds, \quad \Gamma_C(t) := \{x \in \Gamma(t) \mid P_γ(\lambda, d)(x) \le 0\}.$$

The fluid velocity $u$ is here artificial as it comes from $\Omega_f^C$. When choosing $\lambda = [[\tilde{\sigma}_{n,slip}]]$, we ensure that there is no feedback from this artificial velocity to the solid, see (32). For the jump of stresses $\lambda = [[\sigma_{n,slip}]]$, we obtain a mixture of the solid contact condition and the continuity of normal velocities and a feedback might result. This follows analogously to the no-slip case, see (27).

In the left sketch of Figure 10, we see that the minimal distance in the stress-based formulation shows oscillations during the whole contact interval, especially in the second half. The interface jumps back and forth over the contact line many times. The curve corresponding to the formulation using discrete fluxes is much smoother. Similarly, the contact force $J_{P_γ}$ looks smoother, when the flux formulation is used. The reason for this behaviour is the mixture of the interface conditions during contact. On the bottom of Figure 10, we see that the velocity difference $J_{vel,C}$ shows wild oscillations for the stress formulation, while it looks much smoother when using $[[\tilde{\sigma}_{n,slip}]]$. As the artificial velocity in $\Omega_f^C$ has no physical meaning, it is not a drawback
<figure><img src="image_12.png" /><figcaption>Figure 8: Convergence studies under mesh refinement for the artificial fluid formulation by means of the following functionals over time: Top left: Minimal distance $d_{\min}$ of $\Gamma(t)$ to $\Gamma_w$, top right: contact force $J_{P,\gamma}$, Bottom left: Fulfilment of the contact condition $J_{\text{contact}}$. Bottom right: Continuity of velocities $J_{\text{vel,fsi}}$.</figcaption></figure>

that the absolute values of $J_{\text{vel,C}}$ are larger. Due to the feedback of this velocity to the contact conditions, the oscillations appear in the displacement as well.

On the other hand, we should mention that the oscillations are relatively small. Especially those in $d_{\min}$ are almost by a factor $10^3$ smaller than the mesh size $h \approx 1.4 \cdot 10^{-2}$ in vertical direction in this example and are therefore still acceptable.

**Influence of the contact parameter $\gamma_C^0$**

Next, we study the effect of different contact parameters $\gamma_C^0$ for the artificial fluid formulation and $\lambda = \|\tilde{\sigma}_{n,\text{slip}}\|$ on the mesh with 5120 elements. In Figure 11, we show the 'minimal distance' $d_{\min}$ (top) and the contact force $J_{P,\gamma}$ over time for different contact parameters $\gamma_C^0$. The results are similar to the corresponding results for the virtual obstacle problem in Figure 4. For the smallest contact parameter $\gamma_C^0 = 10$, the contact condition is violated throughout the contact interval ($d_{\min} < 0$). The maximum overlap into the artificial fluid domain is again approximately by a factor 30 smaller than the mesh size $h \approx 1.4 \cdot 10^{-2}$. This violation gets smaller, the larger the contact parameter is chosen. The instabilities for the smallest parameter are still much better visible in the contact force $J_{P,\gamma}$. At time $t = 1.89 \cdot 10^{-3}$ the functional shows a huge peak, as the contact condition $d \cdot n_w \le g_0$ is severely violated and it vanishes from $t = 1.94 \cdot 10^{-3}$ to $t = 1.97 \cdot 10^{-3}$, when the contact is in fact shortly released.
<figure><img src="image_13.png" /><figcaption>Figure 9: Convergence studies under mesh refinement for the *relaxed contact formulation* by means of the following functionals over time: *Top left:* Minimal distance $d_{\min}$ of $\Gamma(t)$ to $\Gamma_w$, *top right:* contact force $J_{P_\gamma}$, *Bottom left:* Fulfilment of the contact condition $J_{\text{contact}}$. *Bottom right:* Continuity of velocities $J_{\text{vel,fsi}}$.</figcaption></figure>

For the larger values $\gamma_C^0 \ge 10^2$, the curves are relatively smooth and very similar. Altogether, this shows again that the assumption “$\gamma_C^0$ sufficiently large” in Theorem 4.1 is necessary in order to ensure stability.

**Slip vs no-slip conditions** Next, we compare the effect of slip- and no-slip boundary and interface conditions in Figure 12. Due to the difficulties associated with the artificial fluid formulation and no-slip interface and boundary conditions (see the discussion at the end of Section 3.2.2), we use the relaxed contact formulation on the mesh with 5120 elements in this paragraph.

We show results for

• Slip conditions on the interface $\Gamma(t)$ and the lower wall $\Gamma_w$

• A slip condition on $\Gamma(t)$ and a no-slip condition on $\Gamma_w$

• No-slip conditions on $\Gamma(t)$ and $\Gamma_w$.

Note that the second option is possible, as for the relaxed contact formulation $\Gamma_w \cap \Gamma(t) = \emptyset$.

We observe that the contact condition (or more precisely the relaxed condition $d \cdot n_w \le g_\epsilon$) is earlier active, when using slip-conditions: at $t_C = 1.42 \cdot 10^{-3}$ for slip/slip conditions compared
<figure><img src="image_14.png" /><figcaption>Figure 10: Comparison of the different possibilities to choose the fluxes $\lambda$. Minimal distance $d_{\min}$ to $\Gamma_w$ (top left), contact force $J_{P,\gamma}$ (top right) and velocity difference $J_{vel,C}$ (bottom) integrated over the contact part $\Gamma_C(t)$ of $\Gamma(t)$ over time for computations with $\lambda = [[\sigma_{n,slip}]$ and $\lambda = [[\tilde{\sigma}_{n,slip}]$ for the artificial fluid formulation.</figcaption></figure>

to $t_C = 2.02 \cdot 10^{-3}$ for slip interface and no-slip boundary conditions and at $t_C = 2.23 \cdot 10^{-3}$ for no-slip conditions on interface and boundary. The reason is that the fluid forces, and in particular the pressure, that act against the contact are larger for no-slip conditions, as the fluid can not “slip” out of the contact zone easily. This can be seen in the pressure plot on the right. The pressure is considerably larger from $t \approx 5 \cdot 10^{-4}$ for the no-slip conditions until contact is reached for the slip/slip case at $t_C = 1.42 \cdot 10^{-3}$.

As we are allowing for a small gap between the solid and the ground, these results do not contradict the theoretical results by Gerard-Varet et al [21] discussed in Section 3.3, who showed that (in their configuration with a rigid body) contact can not happen, when no-slip conditions are used on the interface and/or the boundary. As discussed in Section 3.2.1, the basic assumption of the relaxed formulation is that a small or infinitesimal fluid layer remains during contact. On the contrary, the results confirm that contact is more likely to happen for slip-conditions, which is in agreement with the theoretical results.

**Comparison with an explicit ad hoc approach** The probably simplest possibility to combine the FSI model introduced in Section 2.1 and the contact approach described in Section 2.2 is to split $\Gamma$ explicitly in each time-step into a fluid-structure interface $\Gamma_{fsi}(t_{m-1})$ and a contact surface $\Gamma_C(t_{m-1})$ based on the displacement $d(t_{m-1})$ of the previous time-step and to use the
<figure><img src="image_15.png" /><figcaption>Figure 11: Parameter studies for the contact parameter $\gamma_C^0$: Minimal distance of $\Gamma(t)$ to the wall $\Gamma_w$ (Top left: Total time interval, top right: zoom around the contact interval) and contact force $J_{P_\gamma}$ (bottom) over time.</figcaption></figure>

interface condition (14) on $\Gamma_{\text{fsi}}(t_{m-1})$ and the contact condition (10) on $\Gamma_C(t_{m-1})$. A strategy of this type has been used by Hecht & Pironneau [32]. The system of equations reads in the slip case: Find $u \in \mathcal{V}, p \in \mathcal{Q}, d \in \mathcal{W}$ such that $\dot{d} = \partial_t d$ and

$$ (\partial_t u, v)_{\Omega_f(t)} + (\sigma_f(u, p), \nabla v)_{\Omega_f(t)} + (\text{div } u, q)_{\Omega_f(t)} + (\partial_t \dot{d}, w)_{\Omega_s(t)} + (\sigma_s(d), \nabla w)_{\Omega_s(t)} \\ - (n^T \sigma_f(u, p)n, (w-v) \cdot n)_{\Gamma_{\text{fsi}}(t_{m-1})} + \gamma_{\text{fsi}}((\dot{d}-u) \cdot n, (w-v) \cdot n)_{\Gamma_{\text{fsi}}(t_{m-1})} \\ - ((\dot{d}-u) \cdot n, n^T \sigma_f(v, -q)n)_{\Gamma_{\text{fsi}}(t_{m-1})} \\ + \gamma_C(P_{\gamma,s}(d), w \cdot n_w)_{\Gamma_C(t_{m-1})} = (f_f, v)_{\Omega_f(t)} + (f_s, w)_{\Omega_s(t)} \quad \forall v, q, w \in \mathcal{V} \times \mathcal{Q} \times \mathcal{W}, \quad (44) $$

where $P_{\gamma,s}$ is defined in Section 2.2 for the pure solid problem, i.e. without any fluid contributions. We use the same numerical parameters as for the contact formulations presented in this work.

To compare this approach with the artificial fluid formulation we show the minimal distance to the ground $d_{\min}$ and the integral over the normal solid stresses over $\Gamma(t) = \Gamma_{\text{fsi}}(t) \cup \Gamma_C(t)$

$$ J_{\sigma_{s,n}} = \int_{\Gamma(t)} \sigma_{s,n} ds $$
<figure><img src="image_16.png" /><figcaption>Figure 12: Comparison of slip- and no-slip interface/boundary conditions by means of the minimal distance $d_{\min}$ of the interface $\Gamma(t)$ to $\Gamma_w$ around the contact interval (left) and the $L^2$-norm of the pressure over a region $\Omega_{f,\text{mid}}(t)$ around the contact surface before contact (right) over time. Due to the larger pressure before contact, the impact happens later when using no-slip conditions.</figcaption></figure>

on the finer mesh with 20480 elements in Figure 13. While the curves for $d_{\min}$ over the total time interval shown on the top left look similar, a zoom-in on the right shows again that the presence of the artificial fluid leads to an earlier time of impact. Moreover, we observe chattering for the ad-hoc approach at the beginning of the contact interval, i.e. contact is released twice again before the solid stays in contact with $\Gamma_w$. The interface jumps back to the fluid domain, with a (relatively small) minimal distance of approximately $1.6 \cdot 10^{-5}$.

In fact, the functional $d_{\min}$ is not a good indicator to investigate stability for the ad hoc approach, as it is zero, as soon as one point of the interface lies on $\Gamma_w$. Note that this is different for the approaches presented in this work, where the interface can go beyond $\Gamma_w$ (or $\Gamma_{\epsilon}$ for the relaxed approach). In the actual computation, the interface oscillates considerably in each time-step and contact is released and renewed frequently in different points. The functional $J_{\sigma_{s,n}}$ on the bottom left of Figure 13 serves to get a better impression of the instabilities during contact. It oscillates throughout the contact interval including a huge peak at $t = 2.32 \cdot 10^{-3}$. Moreover, we see that the elastic dynamics after the contact are also significantly influenced by these instabilities. Compared to the artificial fluid approach the oscillations in the displacement are significantly larger after contact.

We have also tried to iterate for the splitting into $\Gamma_C(t)$ and $\Gamma_{\text{fsi}}(t)$ within each time-step of the ad-hoc approach, which can be seen as an active-set strategy. This did however not cure the problem, as cycling between different active sets is not prevented.

# 6 Conclusions

We have presented two consistent formulations for fluid-structure interactions with contact, both including a continuous switch between the FSI interface and the contact condition depending on the contact force $P_{\gamma}$. In contrast to certain penalty approaches, the contact force is physically motivated and included in a consistent way in the variational formulations. Our numerical results indicate that the two proposed formulations have better stability properties than the usual ad hoc approaches and no chattering was observed in our computations.

Moreover, we have derived analytically a stability result for a generalised formulation includ-
<figure><img src="image_17.png" /><figcaption>Figure 13: Comparison of an explicit ad hoc approach to include the contact and interface conditions with the approach using an artificial fluid described in Section 3.2.2. We show the minimal distance $d_{\min}$ of the interface $\Gamma(t)$ to $\Gamma_w$ on the top left and a zoom-in on the top right. Chattering, i.e. an alteration between contact and no contact is visible for the ad hoc approach at the beginning of the contact interval. The integral over the normal solid stresses $J_{\sigma_{s,n}}$ shown in the graph on the bottom shows large instabilities for the ad hoc approach.</figcaption></figure>

ing a parameter $\theta \in [0, 1]$. As in the pure solid case (Chouly & Hild [10]), this result implies stability for $\theta = 1$ and stability up to a term including the contact force for $\theta \neq 1$. In our computations, we have however not observed any stability issues for the choice $\theta = 0$ either.

The contact formulations were derived here for the simplified configuration of contact with a fixed and straight wall and using linear models for the fluid and solid sub-problems. The algorithms can be applied to more complex contact configurations by using approaches from the literature to compute the projection and the distances between different surfaces [1, 28, 59]. In particular, the extension to the incompressible Navier-Stokes equations in the fluid and to non-linear elasticity in the solid can be addressed by combining the proposed approach with the arguments recently reported in Mlika et al.[60]. Moreover, Coulumb or Tresca friction can also be incorporated by following Chouly et al. [15, 14].
Acknowledgements

The first author acknowledges support by the EPSRC grant EP/P01576X/1. The third author was supported by the DFG Research Scholarship FR3935/1-1.

References

[1] B. Wohlmuth, Variationally consistent discretization schemes and numerical algorithms for contact problems, Acta Numerica 20 (2011) 569–734.

[2] S. Knauf, S. Frei, T. Richter, R. Rannacher, Towards a complete numerical description of lubricant film dynamics in ball bearings, Computational Mechanics 53 (2) (2014) 239–255.

[3] V. Bruyere, N. Fillot, G. E. Morales-Espejel, P. Vergne, Computational fluid dynamics and full elasticity model for sliding line thermal elastohydrodynamic contacts, Tribology International 46 (1) (2012) 3 – 13.

[4] T. E. Tezduyar, S. Sathe, Modeling of fluid-structure interactions with the space-time finite elements: solution techniques, International Journal for Numerical Methods in Fluids 54 (2007) 855–900.

[5] U. M. Mayer, A. Popp, A. Gerstenberger, W. A. Wall, 3d fluid-structure-contact interaction based on a combined XFEM FSI and dual mortar contact approach, Computational Mechanics 46 (1) (2010) 53–67.

[6] N. D. dos Santos, J.-F. Gerbeau, J.-F. Bourgat, A partitioned fluid-structure algorithm for elastic thin valves with contact, Computer Methods in Applied Mechanics and Engineering 197 (19) (2008) 1750–1761.

[7] M. Astorino, J.-F. Gerbeau, O. Pantz, K.-F. Traoré, Fluid-structure interaction and multi-body contact: Application to aortic valves, Computer Methods in Applied Mechanics and Engineering 198 (45-46) (2009) 3603 – 3612.

[8] S. Frei, T. Richter, An accurate Eulerian approach for fluid-structure interactions, in: S. Frei, B. Holm, T. Richter, T. Wick, H. Yang (Eds.), Fluid-Structure Interaction: Modeling, Adaptive Discretization and Solvers, Radon Series on Computational and Applied Mathematics, Walter de Gruyter, Berlin, 2017.

[9] J. A. Nitsche, Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind, Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg 36 (1970) 9–15.

[10] F. Chouly, P. Hild, A Nitsche-based method for unilateral contact problems: numerical analysis, SIAM Journal on Numerical Analysis 51 (2) (2013) 1295–1307.

[11] F. Chouly, P. Hild, Y. Renard, Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments, Mathematics of Computation 84 (293) (2015) 1089–1112.

[12] F. Chouly, P. Hild, Y. Renard, A Nitsche finite element method for dynamic contact: 1. space semi-discretization and time-marching schemes, ESAIM: Mathematical Modelling and Numerical Analysis 49 (2) (2015) 481–502.
[13] P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Computer Methods in Applied Mechanics and Engineering 92 (3) (1991) 353–375.

[14] F. Chouly, An adaptation of nitsche's method to the tresca friction problem, Journal of Mathematical Analysis and Applications 411 (2014) 329–339.

[15] F. Chouly, P. Hild, V. Lleras, Y. Renard, Nitsche-based finite element method for contact with Coulomb friction, https://hal.archives-ouvertes.fr/hal-01654487/file/authorv3.pdf, (2017).
URL https://hal.archives-ouvertes.fr/hal-01654487

[16] E. Burman, P. Hansbo, M. G. Larson, R. Stenberg, Galerkin least squares finite element method for the obstacle problem, Computer Methods in Applied Mechanics and Engineering 313 (Supplement C) (2017) 362–374.

[17] E. Burman, P. Hansbo, M. G. Larson, Augmented Lagrangian and Galerkin least squares methods for membrane contact, ArXiv e-prints http://adsabs.harvard.edu/abs/2017arXiv171104494B. arXiv:1711.04494.

[18] M. Hillairet, Lack of collision between solid bodies in a 2d incompressible viscous flow, Communications in Partial Differential Equations 32 (9) (2007) 1345–1371.

[19] T. I. Hesla, Collisions of smooth bodies in viscous fluids: A mathematical investigation, Ph.D. thesis, University of Minnesota. (2004).

[20] M. Hillairet, T. Takahashi, Collisions in three-dimensional fluid structure interaction problems, SIAM Journal on Mathematical Analysis 40 (6) (2009) 2451–2477.

[21] D. Gerard-Varet, M. Hillairet, C. Wang, The influence of boundary conditions on the contact problem in a 3d Navier-Stokes flow, Journal de Mathématiques Pures et Appliquées 103 (2015) 1–38.

[22] D. Gérard-Varet, M. Hillairet, Regularity issues in the problem of fluid structure interaction, Archive for Rational Mechanics and Analysis 195 (2) (2010) 375–407.

[23] C. Wang, Strong solutions for the fluid-solid systems in a 2-d domain, Asymptotic Analysis 89 (3-4) (2014) 263–306.

[24] C. Grandmont, M. Hillairet, Existence of global strong solutions to a beam-fluid interaction system, Archive for Rational Mechanics and Analysis 220 (3) (2016) 1283–1333.

[25] C. Grandmont, M. Lukáčová-Medvidová, Š. Nečasová, Mathematical and numerical analysis of some FSI problems, Springer Basel, Basel, 2014, pp. 1–77.

[26] B. Muha, S. Čanić, Existence of a weak solution to a fluid-elastic structure interaction problem with the Navier slip boundary condition, Journal of Differential Equations 260 (12) (2016) 8550–8589.

[27] M. A. Puso, A 3d mortar method for solid mechanics, International Journal for Numerical Methods in Engineering 59 (3) (2004) 315–336.

[28] B. Yang, T. A. Laursen, X. Meng, Two dimensional mortar contact methods for large deformation frictional sliding, International Journal for Numerical Methods in Engineering 62 (9) (2005) 1183–1225.
[29] T. Dunne, R. Rannacher, Adaptive finite element approximation of fluid-structure interaction based on an Eulerian variational formulation, in: H.-J. Bungartz, M. Schäfer (Eds.), *Fluid-Structure Interaction: Modeling, Simulation, Optimization*, Lecture Notes in Computational Science and Engineering, Springer, 2006, pp. 110–145.

[30] G.-H. Cottet, E. Maitre, T. Milcent, Eulerian formulation and level set models for incompressible fluid-structure interaction, ESAIM: Mathematical Modelling and Numerical Analysis 42 (3) (2008) 471–492.

[31] T. Richter, A fully Eulerian formulation for fluid-structure interactions, Journal of Computational Physics 233 (2013) 227–240.

[32] F. Hecht, O. Pironneau, An energy stable monolithic eulerian fluid-structure finite element method, International Journal for Numerical Methods in Fluids 85 (7) (2017) 430–446.

[33] C. S. Peskin, Flow patterns around heart valves: a numerical method, Journal of Computational Physics 10 (2) (1972) 252–271.

[34] D. Boffi, L. Gastaldi, A finite element approach for the immersed boundary method, Computers & Structures 81 (8) (2003) 491–501.

[35] L. Zhang, A. Gerstenberger, X. Wang, W. K. Liu, Immersed finite element method, Computer Methods in Applied Mechanics and Engineering 193 (21) (2004) 2051–2067.

[36] A. Legay, J. Chessa, T. Belytschko, An Eulerian-Lagrangian method for fluid-structure interaction based on level sets, Computer Methods in Applied Mechanics and Engineering 195 (2006) 2070–2087.

[37] A. Gerstenberger, W. A. Wall, An extended finite element method/Lagrange multiplier based approach for fluid-structure interaction, Computer Methods in Applied Mechanics and Engineering 197 (19) (2008) 1699–1714.

[38] E. Burman, M. A. Fernández, An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes, Computer Methods in Applied Mechanics and Engineering 279 (2014) 497–514.

[39] F. Alauzet, B. Fabréges, M. A. Fernández, M. Landajuela, Nitsche-XFEM for the coupling of an incompressible fluid with immersed thin-walled structures, Computer Methods in Applied Mechanics and Engineering 301 (2016) 300–335.

[40] A. Massing, M. Larson, A. Logg, M. Rognes, A nitsche-based cut finite element method for a fluid-structure interaction problem, Communications in Applied Mathematics and Computational Science 10 (2) (2015) 97–120.

[41] P. Hansbo, J. Hermansson, T. Svedberg, Nitsche's method combined with space-time finite elements for ale fluid-structure interaction problems, Computer Methods in Applied Mechanics and Engineering 193 (39-41) (2004) 4195–4206.

[42] D. Kamensky, M. C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, T. J. Hughes, An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves, Computer Methods in Applied Mechanics and Engineering 284 (2015) 1005–1053.
[43] L. Boilevin-Kayl, M. A. Fernández, J.-F. Gerbeau, Numerical methods for immersed FSI with thin-walled structures, Research Report RR-9151, Inria Paris, https://hal.inria.fr/hal-01704575 (2018).

[44] T. Dunne, Adaptive finite element approximation of fluid-structure interaction based on Eulerian and arbitrary lagrangian-Eulerian variational formulations, Ph.D. thesis, Heidelberg University (2007).

[45] S. Frei, T. Richter, A locally modified parametric finite element method for interface problems, SIAM Journal on Numerical Analysis 52 (5) (2014) 2315-2334.

[46] S. Frei, T. Richter, A second order time-stepping scheme for parabolic interface problems with moving interfaces, ESAIM: Mathematical Modelling and Numerical Analysis 51 (4) (2017) 1539-1560.

[47] E. Burman, Ghost penalty, Comptes Rendus Mathematique 348 (21-22) (2010) 1217-1220.

[48] T. Richter, Finite Elements for Fluid-Structure Interactions. Models, Analysis and Finite Elements., Vol. 118 of Lecture Notes in Computational Science and Engineering, Springer, 2017.

[49] S. Frei, Eulerian finite element methods for interface problems and fluid-structure interactions, Ph.D. thesis, Heidelberg University, http://www.ub.uni-heidelberg.de/archiv/21590 (2016).

[50] F. Brezzi, J. Pitkäranta, On the stabilization of finite element approximations of the stokes equations, in: Efficient solutions of elliptic systems, Springer, 1984, pp. 11–19.

[51] R. Becker, M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo 38 (4) (2001) 173-199.

[52] T. J. Hughes, L. P. Franca, M. Balestra, A new finite element formulation for computational fluid dynamics: V. circumventing the Babuška-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Computer Methods in Applied Mechanics and Engineering 59 (1) (1986) 85-99.

[53] E. Burman, P. Hansbo, Edge stabilization for the generalized Stokes problem: a continuous interior penalty method, Computer Methods in Applied Mechanics and Engineering 195 (19) (2006) 2393-2410.

[54] E. Burman, P. Hansbo, Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method, Applied Numerical Mathematics 62 (4) (2012) 328-341.

[55] S. Mandal, A. Ouazzi, S. Turek, Modified newton solver for yield stress fluids, in: B. Karasözen, M. Manguoğlu, M. Tezer-Sezgin, S. Göktepe, Ö. Uğur (Eds.), Numerical Mathematics and Advanced Applications ENUMATH 2015, Springer International Publishing, 2016, pp. 481–490.

[56] S. Frei, An edge-based pressure stabilisation technique for finite elements on arbitrarily anisotropic meshes., http://www.homepages.ucl.ac.uk/~ucahfre/Edgestab_aniso.pdf, submitted ((2017)).

[57] M. Besier, W. Wollner, On the pressure approximation in nonstationary incompressible flow simulations on dynamically varying spatial meshes, International Journal for Numerical Methods in Fluids 69 (6) (2012) 1045-1064.
[58] R. Becker, M. Braack, D. Meidner, T. Richter, B. Vexler, The finite element toolkit Gascoigne3d, http://www.gascoigne.uni-hd.de.

[59] K. Poulios, Y. Renard, An unconstrained integral approximation of large sliding frictional contact between deformable solids, Computers & Structures 153 (2015) 75–90.

[60] R. Mlika, Y. Renard, F. Chouly, An unbiased Nitsche's formulation of large deformation frictional contact and self-contact, Computer Methods in Applied Mechanics and Engineering 325 (2017) 265–288.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>Model</td><td>5</td></tr><tr><td>2.1</td><td>Fluid-structure interaction without contact</td><td>5</td></tr><tr><td>2.2</td><td>Contact model without fluid</td><td>7</td></tr><tr><td>3</td><td>Fluid-structure interaction with contact</td><td>9</td></tr><tr><td>3.1</td><td>Virtual obstacle within the fluid domain</td><td>9</td></tr><tr><td>3.1.1</td><td>Weakly imposed interface conditions</td><td>10</td></tr><tr><td>3.2</td><td>Contact with the boundary of the fluid domain</td><td>11</td></tr><tr><td>3.2.1</td><td>Relaxed contact formulation</td><td>11</td></tr><tr><td>3.2.2</td><td>Contact formulation using an artificial fluid</td><td>12</td></tr><tr><td>3.3</td><td>Slip conditions</td><td>13</td></tr><tr><td>4</td><td>Stability</td><td>15</td></tr><tr><td>4.1</td><td>Generalised contact formulation</td><td>15</td></tr><tr><td>4.2</td><td>A stability result</td><td>16</td></tr><tr><td>5</td><td>Numerical Results</td><td>19</td></tr><tr><td>5.1</td><td>Details on discretisation and stabilisation</td><td>19</td></tr><tr><td>5.2</td><td>Virtual obstacle within the fluid domain</td><td>20</td></tr><tr><td>5.3</td><td>Contact problem</td><td>23</td></tr><tr><td>6</td><td>Conclusions</td><td>31</td></tr></table>
informatics mathematics
Inria

RESEARCH CENTRE
PARIS

2 rue Simone Iff - CS 42112
75589 Paris Cedex 12

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399