A consistent and versatile computational approach for general fluid‐structure‐contact interaction problems

We present a consistent approach that allows to solve challenging general nonlinear fluid‐structure‐contact interaction (FSCI) problems. The underlying formulation includes both “no‐slip” fluid‐structure interaction as well as frictionless contact between multiple elastic bodies. The respective interface conditions in normal and tangential orientation and especially the role of the fluid stress within the region of closed contact are discussed for the general problem of FSCI. A continuous transition of tangential constraints from no‐slip to frictionless contact is enabled by using the general Navier condition with varying slip length. Moreover, the fluid stress in the contact zone is obtained by an extension approach as it plays a crucial role for the lift‐off behavior of contacting bodies. With the given continuity of the formulation, continuity of the discrete system of equations for any variation of the coupled system state (which is essential for the convergence of Newton's method) is reached naturally. As topological changes of the fluid domain are an inherent challenge in FSCI configurations, a noninterface fitted cut finite element method (CutFEM) is applied to discretize the fluid domain. All interface conditions, that is the “no‐slip” FSI, the general Navier condition, and frictionless contact are incorporated using Nitsche based methods, thus retaining the consistency of the model. To account for the strong interaction between the fluid and solid discretization, the overall coupled discrete system is solved monolithically. Numerical examples of varying complexity are presented to corroborate the developments. In a first example, the fundamental properties of the presented formulation such as the contacting and lift‐off behavior, the mass conservation, and the influence of the slip length for the general Navier interface condition are analyzed. Beyond that, two more general examples demonstrate challenging aspects such as topological changes of the fluid domain, large contacting areas, and underline the general applicability of the presented method.

will still remain at this interface but would only be visible on a smaller length scale. The no-adhesive force condition of classical contact mechanics can thereby be reformulated based on the difference between contact and fluid traction.
As a result of this setup increasing the fluid pressure can lift-off two contacting bodies not only in the physical world but also in the numerical model. This effect demonstrates the importance of a defined fluid stress state at every point in the contact zone. Depending on different criteria, fluid models of variable complexity and accuracy on the reduced dimensional manifold of the contact zone can be considered. These criteria include the microstructure of the contacting surfaces (smooth vs. rough), the quantities/processes of interest (e.g., leakage of a sealing), and the macroscopic problem configuration (closed vs. open fluid domain). As this work intends to provide a general FSCI framework, a simple extension based approach to model the fluid in the contact zone is applied herein. The general framework then allows to include also various types of alternative and more sophisticated fluid models in the contact zone in the future. For a model with greater physical depth, the reader is, for example, referred to Reference 3, where a homogenized poroelastic layer formulation is applied to model the fluid flow through contacting rough surfaces.
When looking at the available literature, one can note that a large portion of literature on FSCI formulations is either interested in the analysis of heart valves 2,4-12 and only a smaller set in solving a more general problem setup of FSCI. 3,[13][14][15][16][17] However, while most of those formulations work well for certain problem setups, they suffer from some restrictions preventing their application to more general complex problem classes. In the following, we try to give a brief overview on existing formulations but will not describe the different approaches in detail but rather point to certain features and especially to assumptions or restrictions.
In References 2 and 9, contact in surrounding fluid does not need to be considered due to the chosen problem setup with geometrically separated contact and fluid-structure interfaces. A penetration of the solid bodies is accepted in Reference 4 since contact is not treated explicitly. In Reference 8, contact is included, but in the presented computations only the valve opening phase without significant influence of the contact formulation is analyzed. In Reference 15, no explicit contact formulation is considered and a minimal distance of one mesh cell still remains between two flaps. Using reduced modeling with included contact of the heart valve, Reference 11 avoids the requirement for a general FSCI formulation. In Reference 3, a very general FSCI formulation for contact of rough surfaces is presented, where the interface is modeled via a homogenized poroelastic layer. Such a formulation is very powerful and also well motivated by the involved physical phenomena but it is also more complex and not always needed. In such frequent cases, the approach presented in this article will be a better alternative.
Explicit treatment of the contact is considered in References 3,5-7, and 14 by Lagrange multiplier based contact methods, in References 8,10,13, and 16 by methods based on penalty contact contributions, and in Reference 12 by an approach based on enforcement of equal structural velocity. Interface fitted computational meshes for discretization of the fluid domain are enabled by approaches that require to enforce a nonvanishing fluid gap between approaching bodies and therefore avoid topological changes in the fluid domain preventing degenerated elements. 13,16 Approaches enabling the use of a noninterface fitted discretization, which allow the consideration of "real" contact scenarios and the resulting topological changes of the fluid domain directly, are applied in References 3-8, 10-12,14, and 15. The majority of these formulations consider dimensionally reduced structural models (i.e., membranes and shells), [4][5][6][7][8][10][11][12] whereas bulky structures (i.e., structures of significant thickness as compared to the spatial resolution of the computational discretization in the fluid domain) are only considered in References 3,14,15,and 17. The restriction to slender bodies of the noninterface fitted approaches is often related to issues concerning system conditioning and mass conservation errors close to the fluid-structure interface. This is due to the fact that the discontinuity of the fluid stress between two sides of a submersed solid typically is not represented by the discrete formulation (see, e.g., References 10 and 11), which prevents the analysis of configurations including large pressure jumps. This issue is not a fundamental limitation for noninterface fitted FSI as shown, for example, in References 18 and 19 (without contact), but increases the complexity of such a formulation including the underlying algorithm.
It is surprising that, except for References 3 and 17, none of the referenced works include a substantial discussion concerning the requirement for the fluid state in the contact zone as elaborated earlier. Most works include contact just as a constraint additional to the incorporation of FSI conditions, which is still enforced in closed contact. If such an approach is carried out properly, it can result in a continuous FSCI formulation. Still, the strategy to recover the fluid state in the contact zone, which is required to enforce the FSI conditions, remains an open question. For formulations which circumvent topological changes of the fluid domain 13,16 this issue does not arise directly, as there is always a numerically motivated fluid domain in between the contacting bodies. For all other formulations that support the actual contact of surfaces, the different approaches for the numerical solution of the fluid problem provide a nonphysical fluid state outside the fluid domain in the contact zone by default, which builds the basis to incorporate the FSI interface conditions. Depending on the underlying numerical method and the respective FSCI problem, this can provide, but not necessarily is, a sufficiently good approximation of the fluid state in the contact zone. We would like to point out that this fluid state in the contact zone has an essential influence, for example, to the detaching behavior of contacting solid bodies, for example, a high/low fluid pressure supports/prohibits the separation of two bodies. Another important aspect that should be considered when applying such a strategy of incorporating contact on top of the FSI conditions is that the FSI traction includes a tangential component. In contrast to the normal component of the FSI traction, which serves merely as an offset to normal contact traction, the tangential component directly acts on the contacting surfaces potentially deteriorating the solution accuracy.
In this paragraph, we would like to comment briefly on the so called "collision paradox," which states that for incompressible, viscous fluid with "no-slip" condition on smooth boundaries, contact between submersed bodies cannot occur in finite time (see, e.g., References 20 and 21). This is in contrast to the observation of contacting solid bodies which is also observed when bodies are submersed in fluid. A physical explanation for this paradox can be found in the missing consideration of nonsmooth surfaces (analyzed in, e.g., References 22 and 23) arising from the rough microstructure, or other effects on the micro scale, not considered in the macroscopic fluid model of the "no-slip" boundary condition. 24 Nevertheless, even for the computational analysis of physical models where these effects are not considered, contact still has to be considered. This is a result of the numerical solution approaches which are always accompanied by approximations of the underlying physical model when considering general configurations. If there is no explicit contact treatment within the FSI formulation only fluid forces in the gap keep the bodies apart. But when an artificial collision of solid bodies occurs, for example, during an iterative nonlinear solution procedure, there is no separation force acting as there is no remaining fluid between these bodies. This is shown in References 4 and 7, where a penetration of surfaces can be observed as no contact formulation is considered.
In this contribution, we present a general FSCI formulation considering flows based on the incompressible Navier-Stokes equations interacting with nonlinear elastic solids that are not restricted to slender structural bodies. Therein, contact is not considered as a mere additional constraint on the FSI problem, but focus is rather put on the mutual exclusiveness of fluid-solid and solid-solid coupling. Thus, the application of fluid forces on the interface in the zone of active contact, where typically no good representation of the fluid solution is available, is automatically eliminated. The approach is applied to noninterface fitted discretizations for the fluid domain by the cut finite element method (CutFEM), due to its ability to represent sharp discontinuities of the solution at the interface. This is of essential importance for the discrete representation of the prevalent discontinuity of the fluid stress between opposite sides of (potentially thin) structural bodies. Hence, nonphysically high gradients arising from a continuous fluid solution representation can be avoided (see remark in Reference 6, p 1753). Crucial for the discrete form is the continous transition from fluid-solid to solid-solid interaction, which is achieved by the use of Nitsche-based methods for both constraints.
The CutFEM in general enables the use of noninterface fitted, fixed Eulerian meshes for the fluid discretization in complex and deforming domains. This method is perfectly suited for handling of large interface motion and topological changes of the computational domain, typically occurring for FSCI problems, and therefore is applied here. To enable a determined, continuous transition from the "no-slip" condition to frictionless contact, a relaxation of the tangential constraint is proposed, while retaining the mass-balance in normal direction. This is enabled by a flexible formulation capable of handling the "no-slip" and the "full-slip" limit on the fluid-solid interface. CutFEM has seen great progress in recent years and meanwhile enjoys a solid mathematical base. Initial analysis was performed for the Poisson equation, 25 extended to the Stokes equation, 26,27 and finally, including advection, on the Oseen equation. 28,29 Therein, a so-called "ghost penalty" stabilization 30 guarantees a well-conditioned formulation for arbitrary interface positions. Successful applications of the CutFEM on two-phase flow and FSI are presented in References 31-33 and 18,34,35, respectively. Therein, the "no-slip" interface condition is applied weakly by a Nitsche-based method. The basis for the general Navier interface condition applied in this work was presented for the Poisson equation in Reference 36, extended to the general Navier boundary condition for the Oseen equation in Reference 37, and applied to enforce the tangential coupling condition on the interface of an poroelastic solid and a viscous fluid in References 3 and 38.
To obtain a continuous transition of the discrete formulation from fluid-structure to contact interaction, both FSI and the treatment of contact are enforced via Nitsche's method. A first application of Nitsche's method to contact problems was presented in Reference 39. More recently, the development of Nitsche-type methods for contact problems gained more attention due to the mathematical analysis of symmetric and skew-symmetric Nitsche methods provided by References 40-42 for small deformation frictionless and frictional contact problems. In addition, Reference 43 analyzed a penalty-free variant for the Signiorini-problem. Based on these works, Reference 44 extended Nitsche's method to nonlinear elasticity at finite deformation and Reference 45 to nonlinear thermomechanical problems. Most classical contact formulations F I G U R E 2 Fluid-structure-contact interaction (FSCI) problem setup for two contacting bodies "1" and "2": fluid domain Ω F , solid domain Ω S = Ω S 1 ∪ Ω S 2 , fluid-structure interface Γ FS = Γ FS 1 ∪ Γ FS 2 , the active (closed) contact interface Γ S,c = Γ S 1 ,c ∪ Γ S 2 ,c , overall coupling interface Γ = Γ 1 ∪ Γ 2 , and outer boundaries Γ F,D , Γ F,N , Γ S,D , Γ S,N employ a so-called slave-master concept introducing an inherent bias to the formulation by a (user-defined) choice of the slave and master surface. In the context of Nitsche methods, References 44 and 46 introduced an unbiased variant. The method proposed in Reference 45 is based on a harmonic weighting of the contact stress resulting in an almost unbiased approach as the only bias is introduced by the applied integration rule. In this work, we will extend the method of Reference 45 to a completely unbiased form by integrating not only on the slave but also on the master surface similar to so-called two-half-pass algorithms. 47 Finally, the transition between active and inactive contact has to be balanced carefully with the ambient fluid traction.
The resulting formulation is discretized in time by the one-step-scheme. Finally, the nonlinear system of equations is solved for all unknowns, that is, nodal structural displacements, fluid velocities and pressures, by a Newton-Raphson based procedure. Due to the strong interaction of all involved physical domains for the FSCI problem this is done simultaneously, that is, a monolithic procedure is applied (see, e.g., Reference 48).
Recently, Reference 17 presented a Nitsche-based formulation for FSCI similar to the one derived in this article. Therein, linear Stokes flow and linear elastic solids based on a fully Eulerian description in combination with contact to a rigid, straight obstacle at the fluid boundary is analyzed and stability results for this formulation are shown. In contrast to the extrapolation based strategy proposed in this contribution, which allows for complete topological changes, two strategies based on a thin remaining fluid film are presented in Reference 17 to obtain the fluid stress in the contacting zone.
The article is organized as follows. In Section 2, the governing equations, comprised of the structural and fluid mechanics model as well as conditions on the interface in normal an tangential direction of the FSCI problem, are given. This is followed by a presentation of the discrete formulation, including all volume and interface contributions, and the solution strategy in Section 3. Different numerical examples, capable of analyzing different aspects of the formulation, are presented in Section 4. Finally, in Section 5, a short summary and an outlook are given.

GOVERNING EQUATIONS
In this section, we discuss the governing equations and conditions for all physical domains and interfaces of the FSCI problem. A typical configuration for such a problem is shown in Figure 2. The domain Ω of the overall FSCI problem includes the fluid domain Ω F and the solid domain Ω S . The overall coupling interface Γ consists of the fluid-structure interface Γ FS and the active (closed) contact interface Γ S,c . The different boundaries on the outer boundary Ω are denoted by Γ F,D , Γ F,N , Γ S,D , and Γ S,N . In the following, all quantities * , * with additional "zero"-index * 0 , * 0 are described in the undeformed reference/material configuration, whereas a missing index indicates the current configuration (see Reference 49 for details). An additional "hat"-symbol * , * indicates time-dependent prescribed quantities at the boundaries and in the domains. Prescribed quantities at the initial point in time t 0 are indicated by the"ring"-symbol * , * . The outer boundary of a domain Ω * is specified by Ω * .

Structural domain S
The displacements of every point in the hyperelastic structural domain are governed by the transient balance of linear momentum: Therein, the displacement vector u = x S − X S describes the motion of a material point (with position X S at initial time t = t 0 ), due to deformation of the elastic body, to the current position x S . The structural density in the undeformed configuration is denoted by S 0 , the material divergence operator by 0 ⋅ * , the deformation gradient by F, the second Piola-Kirchhoff stress tensor by S S , and the body force per unit mass byb S 0 . A hyperelastic strain energy function characterizes the nonlinear material behavior and hence provides the stress-strain relation. Therein, the strain is quantified by the Green-Lagrange strain tensor E. The Cauchy stress can be expressed by S = 1 J F ⋅ S S ⋅ (F) T , with J being the determinant of the deformation gradient F. This representation of solid stress S in the current configuration will be required for coupling of the solid domain and the fluid domain on their common interface. Additional initial conditions for the displacement fieldů and velocity fieldv S are required: Finally, to complete the description of the initial boundary value problem for nonlinear elastodynamics, adequate boundary conditions on the outer boundary Ω 0 ∩ Ω S 0 have to be specified with the predefined displacementû on Dirichlet boundaries Γ S,D 0 and the given tractionĥ S,N 0 on Neumann boundaries Γ S,N 0 : The outward-pointing reference unit normal vector on the boundary Ω S 0 is specified by n S 0 . Conditions on the remaining subset of the structural boundary Γ FS 0 ∪ Γ S,c 0 = Ω S 0 ⧵ (Γ S,D 0 ∪ Γ S,N 0 ), where the structural domain is coupled to the fluid domain or contact occurs will be discussed in Sections 2.3 and 2.4. This remaining subset is not part of the outer boundary of the FSCI problem Ω 0 ∩ (Γ FS 0 ∪ Γ S,c 0 ) = ∅.

Fluid domain F
In the fluid domain transient, incompressible, viscous flow is considered. Therefore, the governing equations are the incompressible Navier-Stokes equations which include the balance of mass and linear momentum: Therein, the velocity and the pressure of the fluid continuum at a specific point in space are denoted by v and p, respectively. The constant fluid density is denoted by F , the constant dynamic viscosity by , and the prescribed body force per unit mass byb F . Further, the symmetric strain-rate tensor is defined by Due to the present derivative of the velocity in time, the initial velocity fieldv has to be prescribed: By prescribing adequate boundary conditions on the outer boundary Ω ∩ Ω F , the description of the fluid problem is completed. Thereby the fluid velocityv on Dirichlet boundaries Γ F,D , or the fluid tractionĥ F,N on Neumann boundaries Herein, the Cauchy stress F = −pI + 2 (v) and the outward unit normal n F of the fluid domain is utilized. Again, conditions on the remaining subset of the fluid boundary Γ FS = Ω F ⧵ (Γ F,D ∪ Γ F,N ), which equals the common interface of fluid and structural domain, will be discussed in Sections 2.3 and 2.4. This remaining subset is not part of the outer boundary of the FSCI problem Ω ∩ Γ FS = ∅.

The fluid extension operator
In order to formulate the interface conditions at any point in space x on the overall coupling interface, an extension operator  x ∶ Γ FS → Γ from the fluid-structure interface Γ FS to the overall interface Γ is required. This extension is applied for all quantities solely defined in the fluid domain Ω F and thus for all quantities on the fluid-structure interface Γ FS which are required for the formulation of the interface constraints on Γ. In the following, the extension of any quantity * is denoted by an additional index *  . Exemplary, the extension of the normal fluid stress F nn to a position x on Γ is defined as follows: where the extension origin position x  is properly chosen on Γ FS . The last line in (9) represents the continuity of the extension operator. The applied extension operator for all presented numerical examples is discussed in Section 3.4.4. Alternative approaches to obtain fluid quantities on the overall interface Γ are briefly discussed in the Remarks 11 and 12.

Conditions on the overall coupling interface in normal direction
For the formulation of the interface constraints, which are split in the interface normal direction and in the tangential plane, the solid outward unit normal n = n S will be considered. The normal component of the respective Cauchy stress is denoted as: S nn = S :P n and F nn = F :P n , with the normal projection operator being specified as P n ∶= n ⊗ n. The conditions in the normal direction for purely nonadhesive structural contact configurations are given by the classical Hertz-Signiorini-Moreau (HSM) conditions: which ensure the nonpenetration, the absence of adhesive contact forces, and the complementarity between the contact pressure and normal gap g n . To obtain the normal gap g n , the pointx(x) is obtained as the projection of x along its normal n onto the opposite solid surface; in the case that no such projection exists, we assume g n → ∞. All quantities * evaluated at this projection point will be denoted by a check * .
In the case contacting bodies are surrounded by fluid, the fluid flow in the contacting zone has to be considered properly as discussed in the introduction. Applying the classical HSM conditions (10)- (12) directly would result in the implicit assumption that fluid does not fill the contact zone. For such a configuration, an instantaneous change from zero traction to the traction arising from the ambient fluid in the contact opening zone on the solid boundary would occur. Considering, on the contrary, the presence of (physically reasonable) fluid in the contact zone (on a smaller length scale and not resolved but just modeled at the current macroscopic scale) leads to modified HSM conditions, where a lifting of both bodies occurs for vanishing relative traction of contact (solid) traction and ambient fluid traction. These conditions result in a continuous transition from the FSI traction to the contact traction on the common interface of fluid and solid and vice versa. Then, the conditions on the interface Γ formulated for a specific point x on Γ are: Condition (13) enforces a positive or vanishing gap g n between two solid bodies. In condition (14), a negative or vanishing relative traction has to be guaranteed, at least in the case without adhesive forces that is considered here. Finally, in Equation (15), either a vanishing gap in the contact case of solid-solid interaction or a vanishing relative traction in the case of FSI is enforced. Additionally, the dynamic equilibrium between two contacting bodies has to be formulated: In the contact case, due to the vanishing gap g n , the normal fluid traction equals its projection F nn, =F nn, and therefore the classical dynamic equilibrium between both contacting bodies is recovered. For the FSI case, due to the vanishing relative traction S nn = F nn, , both sides of the equilibrium vanish and as a result Equation (16) is automatically fulfilled. Finally, the mass balance for the motion of solid bodies connected to a fluid domain is given as: Herein, a vanishing normal relative velocity v rel n is enforced solely on the interface Γ FS , which is part of the fluid outer boundary Ω F . Applying an extension to the normal relative velocity v rel n, , this condition is automatically fulfilled on the remaining subset of the interface Γ S,c and hence on the entire Γ.
Remark 1 (Influence of the fluid extension operator). It should be highlighted, that conditions (14) to (16) are expressed by an extension of the fluid stress from the fluid-structure interface Γ FS to the contact interface Γ S,c . The fluid stress extension has an essential influence only close to the condition changing point/curve (Γ FS ∩ Γ S,c ). This point is contained in the origin from which the extension is constructed, namely the fluid domain. Thus, even the application of a simple continuous extension strategy of the fluid stress, which is by definition more accurate close to the fluid domain, provides a sufficiently accurate fluid stress representation for a wide range of problem configurations. Still, we would like to emphasize that the continuous extension operator is considered in this work especially to enable a clear presentation due to its simplicity. In the case that a more accurate physical fluid solution is required in the contact zone, alternative extension based strategies can be considered or appropriate equations to describe the fluid flow in this zone can be solved.

Conditions on the overall coupling interface in tangential direction
In the tangential direction, frictionless solid-solid contact in combination with the general Navier boundary condition as a kinematic constraint between solid bodies and the fluid domain is considered for simplicity of presentation. Then, the following conditions have to be fulfilled on the interface Γ: Herein, the tangential projection operator is specified by P t ∶= I − n ⊗ n. While condition (18) represents the vanishing tangential traction component on the contact interface Γ S,c , condition (19) enforces the dynamic equilibrium between solid and fluid on interface Γ FS . As these two conditions can coincide at the common point Γ S,c ∩ Γ FS only in the case of a vanishing tangential fluid traction ( F ⋅ n ⋅ P t = 0), the general Navier boundary condition (20) with a varying slip length is applied. This condition includes the no-slip boundary condition for a vanishing slip length = 0, which is the common interface condition, successfully applied for macroscopic problem setups. Nevertheless, on smaller scales, due to characteristics such as surface roughness or wettability, an interfacial velocity slip can be observed in a large number of experiments. 24 In this contribution, the main emphasis of applying the general Navier boundary condition is to guarantee continuity for transitions between FSI and frictionless contact solid-solid interaction and to enable a relaxation of the tangential constraint close to the contacting zone. To obtain these properties, an infinite slip length = ∞ is specified close to the common point Γ S,c ∩ Γ FS , whereas a vanishing slip length still allows the consideration of the no-slip condition for the majority of the fluid-structure interface Γ FS representing the macroscopic modeling point of view. Further details on the specification of the slip length for the presented formulation are given in Section 3.3.
Remark 2 (Continuity of the formulation considering frictional contact). It should be pointed out that also for the case when frictional contact is considered, specific treatment of the tangential constraints will be required to result in a continuous problem formulation. This issue arises due to the fact that the fluid wall shear stress on a fluid-structure interface is not automatically equal to the tangential stress resulting from sliding friction of two contacting structures on a macroscopic view. To specify one potential solution to this, in the case of a friction model with vanishing tangential interface traction at the condition changing point/curve Γ FS ∩ Γ S,c , applying the presented strategy directly results in a continuous problem also for frictional contact. The presented general Navier conditions yields a zero tangential fluid traction at the condition changing point Γ FS ∩ Γ S,c . Hence, to ensure continuity, a solid contact friction model has to provide a vanishing tangential traction at this point as well. For instance, this can be achieved using a Coulomb friction law (friction coefficient ) based on the relative normal stress with the total friction bound ⋅ ( S nn − F nn, ).

DISCRETE FORMULATION
In this section, the discrete formulation applied to the numerical solution of the FSCI problem is presented. The spatial discretization of the continuous problem, presented in the previous section, is based on the FEM and temporal discretization by the one-step-scheme is applied. First, the semi-discrete weak forms directly derived from the governing equations, including additional fluid stabilization operators, are given. To account for topological changes in the fluid domain, an elementary feature occurring for the FSCI problems, the CutFEM is applied to the discretization of the fluid equations and is thus discussed in the following. Therein, details on the determination of a consistent discrete set of fluid domain and fluid-structure interface for the contact case are given. The interface conditions, which are split in normal and tangential direction, are incorporated by Nitsche-based approaches. For the normal direction, a single interface traction representation is proposed, automatically incorporating the fluid-structure and contact conditions. A detailed explanation of the resulting contributions by this normal interface traction is given by analyzing the different cases. Further, a Nitsche-based formulation to incorporate the tangential fluid-structure interface condition including potential slip is presented. The specification of the slip length parameter on the interface to enable a continuous transition from fluid-structure coupling to frictionless contact is discussed. Finally, all contributions are treated in a single global system of equations and solved monolithically. Additional details on the solution procedure of the FSCI problem are given at the end of this section. To shorten the presentation only some aspects that help understanding the approach are discussed here, while many more details can be found in the referenced literature for the particular building block methods.
In the following, all quantities, including the primary unknowns, the test functions in the weak form, the domains and interfaces, as well as derived quantities, are discretized in space. Still, no additional index h is added to these discrete quantities for the sake of brevity of presentation. The expressions ( * , * ) Ω and ⟨ * , * ⟩ Ω denote the  2 -inner products integrated in the domain Ω and on the boundary/interface Ω, respectively.

Weak forms for the domains S , F
The weak forms for the structural domain  S , the fluid domain  F , and the overall coupled problem  FS can be derived from Equations (1) and (5)- (6), respectively.
Herein, ( u, v, p) are the corresponding test functions of the primary unknowns (u, v, p). The discrete solution space is created by a spatial discretization consisting of elements containing piece-wise polynomials in an element reference coordinate system which are continuous on the inter-element boundaries. For the pressure p and for each component of the vector-valued quantities fluid velocity v and solid displacement u the discrete approximation space is directly constructed by these functions. All test functions are discretized by the same space as their corresponding primal unknowns. Modifications of these spaces for the incorporation of strong Dirichlet boundary conditions on Γ S,D 0 and Γ F,D are performed in the usual way. For the structural displacements, an interface fitted discretization is applied, meaning that the elements fill the entire domain Ω S . Details on the discretization of the fluid domain, which is noninterface fitted, are given in Section 3.1.2.
As a central part of the formulation, the unique interface traction n introduces the Nitsche-based methods to incorporate the FSI and contact interface conditions. Its specification will be discussed in Sections 3.2 and 3.3. Inserting n in (23) results basically in a consistent stress term and an additional penalty term integrated on the interface Γ. This allows to include the respective dynamic equilibrium in normal direction (14)- (16), as well as in tangential direction (18), and (19) directly in the weak form. As the interface conditions (13)-(20) require a separate treatment of normal and tangential constraints, the normal component nn and the tangential component n ⋅ P t of the interface traction n = nn ⋅ n + n ⋅ P t are treated separately in Sections 3.2 and 3.3. It should also be noted that the interface integrals on the overall discrete interface Γ in (23) are always evaluated on the entire interface specified by the solid domain boundaries. Thus, even for geometrically overlapping contact interfaces, the interface integrals are evaluated twice (with an opposite normal vector).
To extend the interface contribution on Γ FS arising from partial integration of the viscous and pressure contributions in domain Ω F to the overall interface Γ, an additional definition of the fluid test functions on the whole interface Γ is consulted. For the additional interface contributions in (23) vanishing fluid test functions outside of the fluid domain Ω F are considered:

Stabilization of the discrete fluid formulation
In addition to the naturally arising terms of the fluid weak form (22), discrete stabilization operators have to be added to control convective instabilities, to ensure discrete mass conservation, and to guarantee inf-sup stability for equal order interpolation of velocity and pressure: which are connected to one element in  Γ FS , the "ghost penalty" stabilization is applied

The CutFEM utilized for discretization of the fluid domain F
As discussed in the introduction, the CutFEM is applied for the discretization of the fluid domain allowing for a fixed Eulerian computational mesh. Herein, the boundaries and interfaces of the fluid domain are not required to match the boundary of the computational discretization. This beneficial feature of the CutFEM allows the direct handling of large motion or deformation of the solid domain Ω S and even topological changes of the fluid domain Ω F as it is typically occurring for FSCI problems. The discretization concept is visualized for an exemplary contacting configuration in Figure 3.
The typical small penetration of contacting solid bodies in the discrete solution is visualized exaggerated in this figure. This aspect is left aside here and is discussed in detail in Section 3.1.3. All solid domains Ω S 1 and Ω S 2 are discretized boundary and interface matching. The fluid discretization is specified to cover the entire fluid domain Ω F and is not matching to the interface Γ FS . As shown by the exemplary configuration in Figure 3, the outer boundaries of the fluid domain often match the discretization boundary, which does not necessarily have to be. Then, the physical fluid domain Ω F results from "cutting out" the nonfluid domain which is specified by the boundary of the solid domain Ω S and potential nonmatching outer boundaries.
In the following, a brief overview on the most important aspects for application of the CutFEM to the FSCI problem is given. The treatment of all interface conditions is not included here, but presented in Sections 3.2 and 3.3. A general overview of this method is given in Reference 51 including references for further details.
The integration of the  2 -inner products in the fluid weak form (22) has to be performed solely in the physical fluid domain. This domain is described by the outer fluid boundaries Γ F,D and Γ F,N as well as the deforming position of the interface Γ FS including its solid outward unit normal vector n. By separation of the fluid discretization, which is constant in time, in different sets of elements, the numerical integration of (22) can be realized. The computational fluid mesh consists of the sets of elements which are not intersected by the interface Γ FS and affiliated to the fluid domain  F or affiliated to the nonfluid (solid) domain  0 as well as the set of all remaining elements in  Γ FS which are intersected by the interface Γ FS . This domain  Γ FS is split into the physical fluid part  F Γ FS and the nonfluid (solid) part  0 Γ FS , which can be identified by the unit solid outward solid normal vector n. For the "nonintersected" elements in  F Γ FS standard Gaussian quadrature is applied, whereas no integration has to be performed on elements in  0 Γ FS . For the numerical integration of the physical fluid sub-domain  F Γ FS of the intersected elements, the method described in Reference 52, where the divergence theorem is utilized repeatedly, is applied. No integration has to be performed on the remaining sub-domain  0 Γ FS . Due to the arbitrary relative position of the deformed interface Γ FS and the fixed computational fluid mesh, any geometric intersection configuration has to be treated properly. In fact, intersections leading to very small contributions of single discrete degrees of freedom to the weak form (22) are critical. If not handled appropriately, these configurations can lead to an ill-conditioned resulting system of equations or a loss of discrete stability arising from the weak incorporation of interface conditions presented in Sections 3.2 and 3.3. These issues can be tackled by additional weakly consistent stabilization operators added to the weak form (22). Therein, in principle, any nonsmoothness of the discrete extension of the solution into the nonfluid domain  0 Γ FS is penalized. Single degrees of freedom with vanishing contribution in the weak form (22) are then still defined by the smooth extension of the solution, even if there is no physical relevance left. This kind of stabilization is called "ghost penalty" stabilization and was first presented in Reference 30 for the Poisson's problem. The method which is applied here for the stabilization of the fluid equations is analyzed in Reference 29. The operators (26) added to the fluid weak form therein penalizes jumps of normal derivatives of the velocity v and the pressure p: These operators are integrated on a selected set of inner element faces  Γ FS marked in Figure 3 by red lines.
Remark 3 (Existence of the discrete fluid test functions in the ghost domain). It should be highlighted that the discrete test functions ( v, p), in contrary to definition (24), do not vanish in the ghost domain  0 Γ FS , and are evaluated on the inter-element faces for the face-oriented stabilization and "ghost penalty" stabilization (25) and (26) outside of the physical fluid domain Ω F .

Consistent fluid domain F and fluid-structure interface FS representation for the contacting case
The weak form (22) is solely integrated in the physical domain Ω F . This domain is characterized by the nonmoving outer boundaries Γ F,D and Γ F,N as well as the moving fluid-structure interface Γ FS . The discrete motion of the interface Γ FS is given by the general interface Γ and hence by the motion of the solid domain Ω S . It is essential to evaluate the overall weak form on a consistent pair of domain Ω F and interface Γ FS . This aspect is straight-forward as long as no contact between solid bodies occurs, but should be discussed in detail for the case of contacting discrete bodies. The contacting scenario illustrated in Figure 4 results in partial overlap of both solid domains due to the discrete approximation. Therefore, in a first step all parts of the interface Γ which are overlapping-identified by the solid unit outward solid normal vector n-are removed from the "intersection" interface. The corresponding fluid domain Ω F ∪ Ω F * potentially includes small fluid fractions occurring from the discrete contact formulation. To avoid these "islands," the purely numerically occurring segments on the current "intersection" interface are removed additionally, leading to the consulted interface Γ FS . For sufficiently spatially resolved computational meshes, the identification can be simply performed by a predefined maximal ratio of the element size compared to the actual size of the bounding box containing a single fluid fraction. Finally, the intersection of the computational fluid mesh is performed with this interface Γ FS , resulting in a physical fluid domain Ω F which does not include the domain Ω F * . The discrete contact interface is then defined by: Γ S,c = Γ ⧵ Γ FS .

Nitsche-based method on the overall coupling interface in normal direction
All interface conditions in normal direction (13)- (17) are incorporated by Nitsche-based methods. To demonstrate that the formulation enables a continuous transition between the FSI and contact conditions, one(!) representative interface traction nn = n ⋅ n in normal direction is specified which fulfills the following criteria. Insertion of nn in the weak form (23) results always in a Nitsche-based method including a consistent interface stress term and a consistent penalty term which comply with the interface conditions (13)- (17). The value of nn changes continuously for any change of the system state (u, v, p). It is unique in the sense that it depends only on the system state (u, v, p) at the spatial position (which potentially has to be evaluated by an extension) and the interface normal direction and, for example, not on the type of interface on which it is evaluated. Defining the normal interface traction to: with two sufficiently large parameters F  > 0 and S > 0 and a representative solid stress S nn , allows the fulfillment of these requirements. The first term of the minimum in (27) corresponds to enforcing the classical FSI conditions as it is "active" for an existing fluid domain g n > 0 in combination with a vanishing relative traction S nn − F nn, ≈ 0 and relative velocity v rel n, ≈ 0. The second term of the minimum is "active" for a vanishing gap g n ≈ 0 and a negative relative traction S nn − F nn, < 0 and enforces the classical conditions of closed contact. The transition of these two types of conditions, which are combined by the minimum function, is continuous. If no feasible projection exists or the distance exceeds a certain predefined limit and as a consequence the projection is not evaluated, we assume g n → ∞ and as a result the FSI condition is enforced. In the following, a detailed discussion for both sides of the minimum function is given.
In the case that the second term of the minimum in (27) is "active," the conditions of closed contact are incorporated. For this case, the balance of linear momentum across the closed contact interface, in which condition (16) reduces to S nn =̌S nn , is accommodated for by using the same representative solid stress S nn on both sides of the potential contact surfaces. In the most simple case, one of the two potentially contacting solids, for example, Ω S 1 is designated as a so-called slave side and the representative solid stress is chosen as the discrete stress representation of that side  45,53 Additionally, the interface contribution resulting from S g n is required to ensure a stable numerical formulation and to enforce the no-penetration constraint required by (13) in combination with (15). At this point, we would like to highlight that the integration of the (structural) interface integral is performed on the entire interface Γ with the displacement test function at this respective point. For bodies in contact, this entails that one actually integrates twice on the (closed) contact interface, that is, once for each of the contacting surfaces with only the respective test function taken into account. This is in contrast to classical applications of Nitsche's method for two body contact problems 44,46 which integrate over one or both of the contacting surfaces with the jump of the test functions. Our method relates more closely to the so-called two-half-pass algorithms used in Gauss-point-to-segment (GPTS) penalty contact formulations 47 with the difference being the inclusion of the consistent stress term. The approach presented here includes two advantages: first, multi-body and self contact are naturally included. Second, and more importantly, it allows to perform decisions of "active" and "inactive" contact separately on the two surfaces in contact. This is crucial to enable a continuous transition to the FSI condition since in general the effective fluid traction ( F nn, + F  v rel n, ) is slightly different on the two surfaces for the discrete formulation. Appendix A discusses the relation of the presented contact treatment in comparison to existing Nitsche contact algorithms.
In the case that the FSI condition is enforced, the normal interface traction is represented uniquely by the normal fluid traction F nn . Thus, the essential dynamic equilibrium (14) in the case equal to zero and equilibrium (16) due to vanishing contributions on both sides separately are fulfilled . For this choice, a properly scaled, consistent penalty contribution F v rel n is added to guarantee discrete stability of the formulation and to enforce the constraint (17). In addition to the resulting traction and penalty contribution, a skew-symmetric adjoint consistency term is added to the weak form (23): This term allows the direct balance of the contribution of the fluid pressure in addition to the viscous contribution, when introducing F nn in (23). Due to the inherent constraint (17), this additional contribution does not alter the consistency of the formulation. When enforcing the FSI conditions, also a representation of the interface traction by the corresponding solid stress would be possible, but is not considered in the following. The resulting Nitsche-based formulation given by (27), (28), and (23) is similar to the approach in Reference 17, which was applied to a problem of reduced complexity.
Remark 4 (Activation of the contact conditions). Basically, the use of the representative solid stress S nn could lead to a unintended activation of the second term of the minimum (contact formulation) in (27). For two interfaces Γ 1 and Γ 2 with a significant distance between them (g n > 0) and a different value of the fluid stress on both sides (e.g., nn + S g n holds true. As a consequence the contact conditions would be enforced on Γ 1 . Nevertheless, this aspect is more of theoretical interest and due to several reasons this behavior was never observed in any numerical computation. Since a sufficiently large Nitsche penalty parameter S is required anyway, a quite significant difference of the fluid stress F 1 nn and F 2 nn is required to trigger this behavior for an essential gap. Assuming a reasonable computational accuracy, the difference of the fluid stress F 1 nn and F 2 nn will be small for a small gap g n and increase with increasing gap g n . On the other hand, for configurations which include a large distance between both interfaces in combination with almost independent fluid solutions on both interfaces these considerations have practically no relevance since the distance between the interface exceeds a predefined limit and, thus, g n → ∞ is assumed.

A demonstration of the different resulting interface contributions
This paragraph does not add new components to the numerical formulation, but it is meant to provide a boarder discussion of the previously introduced contributions and its implications. To demonstrate the arising interface contributions from incorporation of the normal interface traction (27) into the weak form (23), the boundary integral on the interface Γ is split into the solid-solid contact "+" and the FSI "−" parts: Remark 5 (Relation between the interfaces Γ S,c , Γ FS and Γ, +, Γ, −). For the continuous problem presented in Section 2, integration on the interface subsets Γ, + and Γ, − coincidences with an integration on the contact interface Γ S,c and the fluid-structure interface Γ FS , respectively. Due to the discrete error, this relation does not hold for the discrete formulation, where in general a deviation between these interfaces will occur.
In definition (29), the sign of the scalar  indicates, which side of the min[] function in (27) represents the normal interface traction. In addition to this split of interface Γ in the "+" and "−" parts, a purely geometric split into interfaces Γ FS and Γ S,c was described in Section 3.1.3. As the interface Γ FS is part of the outer fluid boundary Ω F , the fluid state (v, p) and the corresponding test functions ( v, p) are directly defined on this interface without any extension required. Combining these two different subdivisions when performing the integration of the normal traction (23) on the interface Γ leads to four cases (I − IV) which finally needs to be dealt with: A visualization of these four cases for a specific discrete contact configuration is given in Figure 4. In the following, the resulting contributions, which have to be evaluated, are depicted. Vanishing contributions are not included and the extension operator is just applied in the case no direct representation of the corresponding quantity is available on the relevant segment. Further, the skew-symmetric adjoint consistency term introduced in (28) Herein, contribution I equals the classical Nitsche-based method for the imposition of the mass conservation on an fluid-structure interface as applied in References 18 and 35. This method includes an interface traction representation by the fluid stress, a penalty term which is consistent due to the includes mass conservation (17), and the skew symmetric viscous and pressure adjoint consistency term which also includes (17).
The evaluated terms in summand II coincide in principle with Nitsche-based methods for classical contact problems, for example, applied in References 44-46. Here, the interface traction is represented by a one-sided or two-sided weighted solid stress of both contacting bodies with an additional penalty term including the no penetration condition included in (13) and (15). No adjoint consistency terms are applied. Due to the vanishing fluid test functions ( v ∅ , p ∅ ), no contribution to the fluid weak form occurs.
Finally, contributions III and IV arise solely close to the condition changing point/curve  = 0 and the common point/curve of both interface Γ S,c ∩ Γ FS . The impact of these summands compared to contributions I and II is relatively small and so making use of a simple extension of the fluid quantities in (34) is acceptable. Still, both contribution have to be applied to ensure continuity of the discrete problem for a transition from FSI to contact and to guarantee geometrically fitting interface conditions applied onto the fluid domain.
By comparison of contributions (33) and (35), the coincidence of both formulations at the condition-changing point  = 0 and for fulfilled mass balance (17) can be directly seen. As the impact of contribution III is generally small and arises solely close to  = 0, this modification does not have a significant influence onto the performance of the presented formulation.
Remark 7 (Determination of the solid penalty parameter S ). For Nitsche's method, the parameter S = S 0 S with a sufficiently large, positive constant S 0 is required to establish discrete stability of the formulation. Therein, material-and mesh-dependencies of S are considered in S by a local generalized eigenvalue problem as presented, for example, in Reference 45. Larger values of S 0 improve the constraint enforcement (g n = 0), while smaller values typically reduce the error of the consistent stress representation ( S nn ). For the FSCI problem, additional aspects have to be considered. The influence of case IV (see (34)) should be minimized, as it incorporates the extended fluid solution and switching between the cases II and IV during the nonlinear solution procedure should be reduced. A small penalty parameter S supports this behavior as it turns out to reduce the ratio of  ≤ 0 on the interface Γ S,c . As a result, a small but still numerical stable constant S 0 is beneficial for solving the FSCI problem. This aspect is not critical as the constant S 0 is not problem dependent for a properly defined scaling S and the same value can be kept for all computation ( S 0 = 1.0 for all presented numerical examples).
Remark 8 (Determination of the fluid penalty parameter F ). The penalty term in (27) with the parameter F = F 0 F n h −1 Γ balances viscous, convective, and temporal contributions according to Reference 29 and so enables a discrete stable formulation. Therein, F 0 is a sufficiently large positive constant, h Γ an appropriate element volume to interface area ratio, and F v a scaling taking into account the different flow regimes. For the determination of the constant F 0 , constraint enforcement as well as the resulting interface stress error is important. Additionally, for the computed numerical examples, it was observed, that a small penalty parameter F 0 is beneficial for the FSCI problem as it incorporates an inherent relaxation of the kinematic constraints especially close to the point of changing conditions ( = 0) and hence improves the performance of the nonlinear solution procedure. The relevance of this aspect depends highly on the complexity of the considered problem configuration and increases for a reduced accuracy of the applied numerical integration procedure on the interface Γ.
Remark 9 (Applied numerical integration procedure on the interface Γ). For the numerical integration of the contributions (31)-(35) on the interface, nonsmooth and noncontinuous functions on single solid boundary elements have to be integrated. These kinks and jumps potentially occur due to element boundaries of the contact partner or on the intersection of the interface with fluid element boundaries. To enable an accurate numerical integration, each solid boundary element has to be split by all other element boundaries involved and a numerical integration rule has to be specified, for example, by triangulation, on these segments (see, e.g., Reference 54). For the numerical examples presented in the following, this most accurate approach was not applied. Instead, the numerical integration points on the interface Γ are constructed to account for the intersection of the interface with fluid element boundaries solely. To account for the integration of the discontinuous solid stress in the contacting case, an increased number of integration points is applied.

Nitsche-based method on the overall coupling interface in tangential direction
The tangential component of the interface traction n has to fulfill the traction free condition (18) due to the consideration of frictionless contact on the contact interface Γ S,c : Further, the dynamic equilibrium (19) and the Navier slip boundary condition (20) have to be fulfilled on the fluid-structure interface Γ FS . A representation of the unique tangential interface traction by: complies with these condition. For the limit cases no-slip (slip length = 0) and free-slip (slip length = ∞), the tangential interface traction reduces to n ⋅ P t = ⋅ P t and n ⋅ P t = 0, respectively. Incorporating of the tangential interface traction in the weak form (23) and adding an additional consistent skew-symmetric adjoint term results in the contributions: It can be directly seen that this formulation is consistent, as the term in the first line includes the naturally arising fluid stress applied on fluid and solid boundary due to the balance (19) and the additional terms include directly the constraint (20). Theses additional terms are present to guarantee a discrete stable formulation and to enforce the kinematic constraint. They balance the destabilizing effects of the viscous boundary integral occurring in line one and the term of similar structure in line two. The penalty parameter in tangential direction F t,0 needs to be a positive and sufficiently large constant. This Nitsche-based contribution for the general Navier interface condition is based on the formulation presented and analyzed in Reference 36 for the Poisson problem and Reference 37 for the linearized fluid problem. It was successfully applied to enforce the coupling conditions between a poroelastic structure and fluid flow in References 3 and 38.

Definition of the slip length
As motivated already in Section 2.4 for the overall problem, the no-slip interface condition = 0 on Γ FS has to be applied. Solely close to the contacting zone, a relaxation of this constraint is designated. A continuous transition between the tangential FSI condition (19)- (20) and the tangential frictionless contact condition (18) can be guaranteed for an infinite slip length = ∞ on Γ FS ∩ Γ S,c . The applied interpolation between these limiting points is given by Herein, the minimal value of the gap g n between two solid interfaces to apply the no-slip interface condition is specified by the fluid element size h. The interpolation function can be specified by the constant reference slip length 0 . It should be pointed out that for a reduction of the fluid element size h, also the range of influence for this modification compared to a pure no-slip condition is reduced. For small scales, an alternative formulation for the slip length in relation (40) due to the underlying physical slip can improve the accuracy of the interface condition. A second advantage of allowing a certain amount of slip on the interface is to avoid blockage of single fluid elements between approaching surfaces due to the insufficient discrete solution space. This aspect is less essential for a weak enforcement of the interface condition by Nitsche's method with an appropriately chosen penalty parameter F t,0 than for a strong enforcement of interface conditions (see Reference 55 for a comparison of strong enforcement and weak imposition of fluid boundary conditions). To give an example, in Reference 5, this issue is resolved by a modification of the fluid-structure interface constraint close to contact.

Overall formulation for the coupled FSCI problem
Finally, by making use of the corresponding interface traction representation in normal (27) and in tangential (36)-(37) direction in the weak form (23) of the overall coupled problem, and summing up all additional contributions including interface adjoint consistency terms and discrete stabilization contributions, the following semi-discrete weak form of the FSCI problem has to be solved. Find (u, v, p) such that ∀( u, v, p): For discretization of the weak form (41) in time, the one-step-scheme is applied to the occurring time derivatives in solid (21) and fluid weak form (22) with an equal time integration factor . This procedure leads to a nonlinear system of equations of form  n = 0 for each discrete instance in time (index n) in the interval [t 0 , t E ]. An iterative solution scheme based on a Newton-Raphson procedure is applied to solve this nonlinear system: Herein, all equations arising from the overall weak form (41)  For a sufficiently small value of a residual norm || i n || < , the current iteration state approximates the solution state for this timestep i+1 n = n . Based on the previously computed solution state, the solution at the next discrete instance in time is computed by another Newton-Raphson iteration procedure. It should be highlighted that the discrete solution space of the fluid pressure and velocity potentially changes between different iteration steps or time steps. As a consequence, a procedure to transfer the previously evaluated solution state to the current solution space is required for the nonlinear solution procedure and the discrete time integration. Details on this procedure and the overall solution algorithm for the CutFEM-based FSI, which builds the basis for the algorithm applied herein, are presented in Reference 35. From an algorithmic point of view, solely the evaluation of different contributions on the interface varies from the presentation therein.
In the following, we present algorithmic details for the monolithic solution algorithm, applied for the computation of the subsequently presented numerical examples. The configurations of the presented examples (except for the first validation example) are chosen to challenge the presented FSCI formulation. Thus, it is ensured that the fluid-structure coupling, the solid-solid contact, as well as the change between these conditions have an essential impact to the overall problem. To guarantee a strong interaction between the structures and the fluid, equal initial densities are considered within all domains. Specific strategies to enable a robust solution of the resulting highly nonlinear problems are discussed in the following.

3.4.1
Damping strategy for the update procedure

Update strategy for geometric intersection
A simple procedure to avoid deterioration of the convergence behavior in the Newton-Raphson procedure due to "algorithmic"-discontinuities arising from geometric tolerances in the algorithm intersecting the interface Γ FS and the computational fluid domain Ω F is applied. Herein, the geometric intersection (includes the creation of numerical integration points in the physical fluid domain Ω F and on the fluid-structure interface Γ FS ) is just updated as long as the maximal displacement increment ||Δ i n || ∞ > geom in an iterations step exceeds a specified valued. For the remaining iteration steps, the intersection information of the previous iteration step is applied.

Solution space update strategy
As explained in detail in Reference 35, the solution space is updated dynamically within the iterative solution procedure for solving the system of nonlinear equations. For classical FSI computations without structural contact, applying this procedure typically results in a constant solution space after few iterations. Nevertheless, including contact increases the sensitivity of the formulation with respect to changing solution spaces. This aspect can result in periodically repeating changes of the solution space within the iterative solution procedure for specific geometric configurations and so prohibits the convergence of the scheme. The reason for this behavior is a discontinuity in the discrete formulation, which arises due to the change in the considered set of faces in the weakly consistent "ghost penalty" stabilization (26) and the face-oriented fluid stabilization (25). The influence of this effect onto the convergence of the Newton-Raphson based procedure is especially relevant in the case when two physical fluid domains are merged or separated. Thus, when exceeding a maximum number of iterations in the nonlinear solution procedure, no reduction in the computational nodes carrying fluid degrees of freedom is performed anymore. The fluid solution space is then just enlarged within the actual timestep.
To retain a solvable system of equations, the "ghost penalty" stabilization and the face-oriented fluid stabilization has to include the faces connecting all additional degrees of freedom to the physical domain. This strategy leads to a constant set of faces considered for the stabilization during the Newton-Raphson based procedure and as a result avoids the occurring discontinuity in the discrete formulation. With this modification, the consistency of the formulation is not touched. Only some additional fluid degrees of freedom, which represent an extension of the solution in the nonphysical domain, are appended to the system. To ease the use of this strategy, the discrete fluid solution space is constructed by a maximum of one set of fluid unknowns on each node in the presented computations in the following. As long as no slender solid bodies are considered, this restriction still results in an appropriate discrete fluid solution space (for more details on multiple sets of fluid unknowns on single computational nodes see Reference 28).
Remark 10 (Discontinuity of the discrete interface contribution). To give a comprehensive overview it should be mentioned that also the viscous interface contributions in (23) can be discontinuous in the FSI case when the interface is parallel to the underlying fluid element boundary which is crossed. For this specific case, an instantaneous change of the interface area in the fluid element occurs also for small interface motion. Since the viscous traction is based on the fluid velocity gradient, which is discontinuous between fluid elements, a discontinuity of the fluid residual occurs.

Applied extension operator
In Sections 2.3 and 3.2, an extension operator  x is required to extend the scalar fluid quantities, normal stress, relative normal velocity, and fluid Nitsche penalty parameter from the fluid-structure interface to the contact interface Γ FS → Γ S,c .
In the numerical examples presented, a very simple operator is applied. Herein, starting from the coordinate x on interface Γ S,c , the closest point x  to Γ FS ∩ Γ c is computed. At this point, the scalar quantity is evaluated. Then a constant extension is applied and as a result the computed value of the scalar quantity equals the extension. This approach, which is based on a single point on Γ FS ∩ Γ c for each extension operation, potentially results in discontinuities of the extended scalar field. Nevertheless, since the extended quantities are mostly only essential close to the fluid domain this aspect is often not critical. The constant extension from the closest point is only one possibility to specify the extension operator. As long as the continuity criteria in (9) is fulfilled, alternative specifications can be utilized without requiring any modification to the presented overall formulation.

NUMERICAL EXAMPLES
In

Falling, contacting, and lifting of a rounded stamp
The first presented numerical example, a simple configuration, including the falling, contacting, and lifting of a rounded stamp, is considered to analyze basic properties of the presented formulation. Due to the symmetry of this configuration, just the half rounded stamp Ω S 2 and fluid domain Ω F are considered.

Problem description
The geometric setup and basic boundary conditions are visualized in Figure 5. The solid domain Ω S 1 is rigid and fixed in space by a Dirichlet boundary condition on the overall domain. In the initial phase, the stamp is exposed solely to a prescribed constant-in-time Neumann load on the boundary Γ S,N in negative y-direction (see Figure 5 (right)), which describes the material behavior of the solid domain Ω S 2 , with Young's modulus E = 100 and Poisson's ratio = 0.0. To analyze the presented formulation, two different spatial resolutions are applied. For the "coarse" variant, the fluid mesh consists of 16 × 24 = 384 elements (in −1.5 ≤ x ≤ 0 and 0 ≤ y ≤ 1) and the solid mesh of domain Ω S 2 is created by 400 elements. In the "fine" variant, 64 × 96 = 6144 fluid elements and 6400 elastic solid elements are used. The weighting of the solid contact stress is purley based on the domain Ω S 2 due to the rigid domain Ω S 1 . The reference slip length is set to

Numerical results and discussion
In Figure 6, the vertical displacement of the spatial point with initial position X S = (0, a) of the solid domain Ω S 2 is depicted. Comparing the "fine" and "coarse" discretizations shows a good agreement down to a gap of approximately two coarse fluid elements (vertical displacement u y (0, a) = −0.375), where both variants start to deviate significantly. Solely the "slip" variant leads to a fundamentally different impact behavior, which is clear due to the nonphysical boundary condition applied to the viscous fluid. All variants lead to contact in finite time, even though this phenomenon is not expected for the no-slip variants theoretically (see References 20 and 21). The explanation for this (realistic) behavior lies in the inherent constraint relaxation arising from the weak imposition by Nitsche's method. As soon as the solution can no longer be resolved sufficiently, a tangential slip occurs numerically also for the no-slip model. To substantiate this explanation, a variant with increased tangential penalty parameter F t,0 = 1000 F,std t,0 (with the parameter of the standard configuration given by F,std t,0 ) is computed, which reduces the numerical slip and thus results as expected in a slower approach velocity.
In the following, the difference between the no-slip condition and the Navier slip condition of the computed solution is discussed. As expected, the Navier slip variant results in an increased velocity, starting from fluid gaps smaller than one fluid element (see definition of the slip length in Section 3.3). Still, the difference between both approaches is not substantial (compared to the error between "coarse" and "fine" resolution). While this simple configuration allows to solve the FSCI problem for both interface conditions, applying the Navier slip condition seems to be beneficial for general configurations in two aspects. First, independent of the approach applied for the imposition of the interface condition, a controlled way of relaxation of the tangential constraint can be incorporated. Second, this type of condition is required to allow for a continuous problem formulation on the interface.

F I G U R E 6
Vertical displacement of the spatial point with initial position X S = (0, a) of the solid domain Ω S 2 for different computed variants over time (left overview, right detail): "noslip fine" (slip length on the interface Γ FS specified to = 0, computed with the "fine" discretization variant), "navslip fine" (slip length on the interface Γ FS as defined in Section 3.3 ( 0 = 0.1), computed with the "fine" discretization variant), "noslip" (slip length on the interface Γ FS specified to = 0, computed with the "coarse" discretization variant), "navslip" (slip length on the interface Γ FS as defined in Section 3.3 ( 0 = 0.1), computed with the "coarse" discretization variant), "noslip incpen" (configuration as "noslip" with increased tangential Nitsche penalty constant F t,0 by a factor of 1000), "navslip incpen" (configuration as "navslip" with increased tangential Nitsche penalty constant F t,0 by a factor of 1000), "slip" (slip length on the interface Γ FS specified to = ∞, computed with the "coarse" discretization variant). The horizontal black dash-dotted lines (thick lines for coarse mesh) indicate the fluid element boundaries F I G U R E 8 Comparison of computed flow rates and flow rate errors for the "coarse" and "fine" mesh resolution for the no-slip interface condition. Herein, Φ is the flow rate on boundary Γ F,N , Φ 1 err the overall flow rate error, and Φ 2 err the flow rate error on the interface Γ FS Herein, Φ 1 err indicates errors in the overall mass balance, and Φ 2 err characterizes the mass balance errors due to the weak imposition of the interface condition by the Nitsche method. When analyzing the overall mass balance Φ 1 err , an unexpectedly small error for this mesh resolution can be observed. An explanation to this effect is given in the following. The discrete fluid mass balance is comprised of the divergence term in (22), the weakly consistent face-oriented stabilization operators (25) and "ghost-penalty" stabilization operators (26), and the skew-symmetric adjoint consistency term on the interface (28). Partial integration of the divergence term in (22) for the fluid balance of mass is performed and the resulting terms are combined with adjoint consistency term the (28) in (46).
Discrete fluid balance of mass: It can be observed that the fluid velocity in the boundary integral in the second line is replaced by the solid velocity on the interface Γ FS . The skew-symmetric adjoint consistency term (28) acts therefore as a compensation term for the violation of the balance of mass on the fluid-structure interface. Hence, the error Φ 1 err is not influenced by the accuracy of the FSI constraint enforcement but is solely attributed to the stabilization terms from the CIP and the GP stabilization. In addition, the finite convergence tolerance of the nonlinear solution procedure yields perturbations in the error level depending on the remaining residual. Finally, the interface error Φ 2 err is observed to be significantly larger than the overall error Φ 1 err . Comparing the "coarse" and the "fine" mesh resolution allows the analysis of the spatial convergence of this error. For the time range with similar flow rates (Φ coarse ≈ Φ fine), a reduction in the error, approximately of second order with respect to the fluid mesh element size h, can be observed.
To give a comprehensive view of the balance of mass for this FSCI formulation, the results for the application of the no-slip condition on the entire interface are also given in Figure 8. No significant difference between both results can be observed. Due to the logarithmic axis scaling, a deviation for the small flow rates (600 ≤ t ≤ 1000) after contact established can be observed. As this difference does not essentially influence the lift-off procedure afterward, the principal discussion done for the Navier slip condition holds also for the no-slip condition.
In the following, the spatial convergence behavior of the formulation with a focus on contact is analyzed. For studies that do not include contact the reader is referred, for example, to References 18 and 35. Therefore, various computational meshes for the fluid and structural discretization are created based on the discretization variants "coarse" and "fine" by merging (four into one) or splitting (one into four) all bilinear elements. The respective element size for the considered meshes is approximately h ≈ 1/8, 1/16, 1/32, 1/64. In Figure 9, the vertical displacement of the initial point X S = (0, a) of the solid domain Ω S 2 is shown. While in the initial phase no difference between the variants is observable, for a remaining gap smaller than 2h, a significant deceleration due to the insufficient fluid solution space occurs. The finer the spatial resolution, the later contact establishes. Still, this configuration also demonstrates that there is no lower limit of the element size h such that contact will not occur in finite time. Due to the inherently included contact treatment in the FSCI formulation, this does not lead to a significant deviation from the exact solution where no penetration of solid bodies F I G U R E 9 Vertical displacement of the spatial point with initial position X S = (0, a) of the solid domain Ω S 2 for computational meshes with different spatial resolution over time. Computed with the Navier slip interface condition F I G U R E 10 Computed average overall interface traction, FSI traction, and Contact traction for computational meshes with different spatial resolution over time. Computed with the Navier slip interface condition. The overall traction includes all contributing of cases I − IV, the FSI traction includes case I, and the contact traction includes cases II − IV (specified in the interface contributions (31)- (35)). The average traction is reconstructed from the nodal interface contributions of (31)- (35)) to the overall weak form on the solid mesh divided by the interface area on the circular part of Γ FS 2 F I G U R E 11 Error of the average overall interface traction of the coarser computational meshes based on a comparison to the values computed with the finest discretization h ≈ 1/64 occurs. For the computational mesh with h ≈ 1/8, a small penetration of the solid bodies in the contact phase can be observed. A comparison of the different computations for the lift-off process reveals similar behavior as for the contacting process.
In Figure 10, the computed average interface traction ∫ Γ FS 2 circ || n ||dΓ∕∫ Γ FS 2 circ 1dΓ on the circular part of the interface Γ FS 2 circ ⊂ Γ FS 2 for all considered computational meshes is shown. Due to the different points in time where contact or lift-off is initiated, the FSI traction and the contact traction differ significantly for all performed computations. Nevertheless, the overall traction, which is effectively acting on the interface, is almost equal for all variants. This aspect demonstrates that including a contact formulation allows to retain physical meaningful behavior of the formulation also for the case of prevalent numerical errors. The difference of the computed average overall interface traction of the coarser discretizations to the finest computation mesh is shown in Figure 11. In general, a convergence behavior for increasing spatial resolution can be observed. Nevertheless, due to the temporal shift of the essential processes especially in the initial and final phase, this is not true for every point in time.
To analyze the influence of the Nitsche penalty constant S 0 , computations with different parameters S 0 = 0.1, 1, 10,100 are performed. For computations with a smaller value of the penalty constant S 0 = 0.01, the nonlinear solution procedure did not converge, which is in agreement with the requirement for a sufficiently large contact S . In Figure 12, the computed local gap g n and the effective interface traction S nn + S g n for contact, corresponding to case II in Section 3.2, F I G U R E 12 Computed local gap g n (left) and effective interface traction S nn + S g n (right) at t = 1000 for different Nitsche penalty constants S 0 = 0.1, 1, 10,100. The dotted line corresponds to the interface Γ FS 1 and the solid line to interface Γ FS 2 are shown. For the smallest value S 0 = 0.1, a large penetration of the contacting bodies occurs and additionally in a relatively large interface area with positive gap contact is enforced. With increasing penalty constant S 0 , the amount of penetration of the contacting bodies reduces drastically. Nevertheless, as discussed in Remark 7, to reduce the influence of the extended fluid state in the contact zone and to allow for a certain amount of constraint relaxation, which is beneficial for the performance of the nonlinear solution procedure, a small but numerically robust parameter is preferred. For the effective interface traction S nn + S g n , no significant difference of the variant with S 0 ≥ 1 can be observed. For these variants, the effective interface traction is almost equal on both contacting interfaces (comparing the dotted and continuous lines) and, thus, the balance of linear momentum in the contact zone is fulfilled. In contrast to that, for S 0 = 0.1 and x ≤ −0.25, a deviation of the effective interface traction on both sides of the interface due to the deviating interface normal vectors is observable. As a consequence of these findings, the constant is set to S 0 = 1 for all computations.
In Figure 13, a detailed view of the contacting zone for different points in time is given. Three different types of traction are visualized by arrows, namely, the overall traction, the FSI traction, and the contact traction. At t = 100 (first row in Figure 13), the body Ω S 2 approaches Ω S 1 and as a result a high pressure peak occurs in the smallest constriction. This peak is almost equal to the FSI traction concluding that viscous traction is not significant. At t = 340 (second row in Figure 13), the majority of the external load is carried by the contact traction. For the overall traction, the continuous transition of FSI traction and contact traction can be seen. An essential part of the external load at t = 2020 (third row in Figure 13) is carried by the FSI traction, but due to the fluid inertia there is still contact at the area around x = 0. Finally, at t = 2420, the structural body Ω S 2 completely lifted again and so the lowest pressure and FSI traction can be identified in the smallest constriction.

Elastic pump
In the following example, an elastic fluid pump powered by an external load is analyzed. This configuration includes large deformation of the solid domain and a periodically changing topological connection of the fluid domain. Large fluid pressure discontinuities when crossing the valves occur which need to be represented properly by the fluid solution space.

Problem description
The geometric setup and basic boundary conditions are depicted in Figure 14.

Numerical results and discussion
In Figure  as the periodic response to the periodic load with A = 1 and was already discussed in detail previously. Now, analyzing the subsequent interval t ∈ [0.15, 0.3] with A = 1.5, after a transition phase in the fourth load cycle the pump exhibits again an almost periodic behavior for the last two load cycles. It can be seen that oscillations with higher frequencies occur than for the smaller load amplitude, which is tackled by a reduced time step size in the time integration scheme.
To make a statement on the performance of the pump, the volume transported through the pump is presented in Figure 16 (right). It can be seen that in each cycle the transported volume through Γ in at first is smaller than through Γ out mainly due to the volume change in the fluid chamber. The difference in the transported volume is smaller for A = 1 than for A = 1.5 as larger deformation occurs. Analyzing the transported volume per cycle, it can be seen that, for the smaller amplitude, each cycle transports approximately 0.074, whereas the higher load amplitude leads to a slight transport opposite to the design flow direction.
To generate understanding for this phenomenon, four exemplary points in time with load amplitude A = 1.5 are selected and shown in Figure 17. Compared to the load with amplitude A = 1.0, higher fluid velocities occur leading to higher pressures and finally an increase of the interface traction, at t = 0.167. This fluid state leads to a nonsymmetric deformation of the flaps in the right valve. As it can be seen at t = 0.174, finally the lower flap snaps through and as a result the right valve does not prevent flow properly anymore. For t ∈ [0.1791, 0.1822], contact between the upper and lower part of the fluid chamber occurs, prohibiting the flow in the chamber. At t = 0.19, the lower flap of the right valve starts to snap back, whereas the left valve is exposed to large nonsymmetric deformation. Finally at t = 0.194, the left valve has snap-through, allowing for flow opposite to the design flow direction. In short, a load amplitude of A = 1.5 is beyond the maximal load resulting in a proper operation of the elastic pump.

F I G U R E 18 Geometry and boundary
conditions for the flow-driven squeezed elastic structure (left). Visualization of the discretization for the solid domains Ω S 1 and Ω S 2 with black lines indicating element boundaries (right) Nevertheless, from a computational point of view, it is noteworthy that the presented formulation demonstrates to be applicable also for these rather complex scenarios and hence promises to be a rather general tool. Processes beyond the intended design can be computed without requiring changes to the problem setup. In this example, unexpected deformation and topological changes to the fluid domain were handled without any modifications to the problem setup.

Flow-driven squeezed elastic structure
In the following, a configuration is considered where an initially cylindrical elastic body Ω S 2 is squeezed through an elastic structure Ω S 1 by the load of the surrounding fluid flow. This configuration is designed to test the formulation's capability to handle frequent changes between the fluid-structure interface and the contact interface including large contacting areas and essential topological changes.

Problem description
The problem setup of the this example, including all dimensions and basic boundary conditions, is visualized in Figure 18 (left). All dimensions of solid body Ω S 1 , which are not explicitly indicated in this figure, are defined by symmetry and replication of the given dimensions (e.g., all unspecified radii are equal to r 2 ). On the inflow boundary Γ in , a time-dependent, parabolic velocity profilev = [0, −100(1 − x 2 )4000t] T for t ∈ [0.0, 0.00025] andv = [0, −100(1 − x 2 )] T for t ∈ [0.00025, 0.016] is prescribed as Dirichlet boundary condition. On the outflow boundary Γ out , a zero traction Neumann boundary condition is prescribed.
The material properties of the incompressible fluid are specified by the density F = 10 −6 and the dynamic viscosity = 10 −5 . The initial density in both solid domains equals the fluid density S 0 = F = 10 −6 . Similar to the numerical examples presented previously, a Neo-Hookean material model with strain energy function (43) is considered for both solids. The parameters of the material model in the squeezed domain Ω S 2 are given by E S 2 = 100 and S 2 = 0.3, whereas the outer domain Ω S 1 has an increased stiffness by E S 1 = 200 and S 1 = 0.3.
The structured computational mesh applied for the discretization of the fluid domain consists of 120 × 300 = 36000 bilinear elements(in −1 ≤ x ≤ 1 and 0 ≤ y ≤ 5). The solid domain is discretized fitted to the interface Γ FS by 4890 elements in domain Ω S 1 and by 1562 elements in domain Ω S 2 . The solid discretization is given in Figure 18 (right). Compared to the examples presented previously, the penalty parameters constants F 0 and F t,0 are divided by a factor of 7, in order to relax the kinematic constraints and thus support the nonlinear solution procedure (see Remark 8). With this modification, the penalty parameters are still large enough to provide discrete stability of the formulation. A contact stress based on  Figure 20 harmonic weighting between the stress representation of both solid domains, as discussed in Section 3.2, is applied. The reference slip length is set to 0 = 0.1. The temporal discretization is preformed with = 1 and a time step size of Δt = 0.00002 for t ∈ [0, 0.0056] and Δt = 0.000005 for t ∈[0.0056, 0.016] to account for the varying dynamic of the coupled system.

Numerical results and discussion
In Figure 19, the computed fluid velocity and the computed deformation of solid domains are presented for different points in time. Following the different snapshots, the motion of the solid domain Ω S 2 can be observed. In the initial phase (0 < t < 0.00324), a vertical motion of Ω S 2 is induced by the fluid flow. At t = 0.00324 contact between Ω S 2 and the right part of Ω S 1 occurs. Starting from t = 0.00386, additional contact between Ω S 2 and the left part of Ω S 1 establishes. Therefore, the topology of the fluid domain changes from one connected domain, to two separated fluid domains. In the subsequent phase (0.00386 < t < 0.006), the pressure in the upper fluid domain increases, which leads to a squeezing process of Ω S 2 and a deformation of Ω S 1 and thus a storage of elastic energy. For t > 0.0065, an acceleration in vertical direction of Ω S 2 can be observed by the transfer of the elastic energy via contact forces. Finally at t = 0.00668, contact between both solid bodies is released and a single connected fluid domain reoccurs. Reestablishing contact at t = 0.00713 of Ω S 2 and the right part of Ω S 1 , this principal process repeats for two additional cycles. Nevertheless, due to the varying geometric setup around the three smallest constrictions, the physical process is not repeated exactly and thus the robustness of the algorithm is tested for this challenging configuration. Finally, at t = 0.015155, both solid domains separate for the last time. In the remaining period, the fluid traction is exclusively acting on the interface Ω S 2 .
To give a more detailed view of the computed process, the fluid solution as well as the interface traction for four exemplary points in time are shown in Figure 20. First, the point in time just before contact occurs t = 0.003 is discussed. Due to the small cross-section of the connection between the upper and lower part of the fluid domain, the pressure in the upper part is already increased. Therefore, an essential fluid flow can be observed between Ω S 2 and the left part of Ω S 1 . The distance in the smallest constriction for the right part leads to an increased fluid pressure compared to the ambient pressure and thus an FSI traction separating the two bodies occurs. At t = 0.005, contact between both solid domains is established in two positions. Due to the inflow on Γ in , the pressure in the upper part of the fluid domain is increased, which leads to an increased FSI traction on the affected part of the interface. Although the maximal contact traction is significantly higher than the FSI traction, there is a continuous transition along the interface. The y-components of the resulting FSI force and contact force are almost in balance, and as a result only a very slow motion of the system (see fluid velocity) is observed, continuously adapting to the increasing pressure difference. At t = 0.0065, this state changed fundamentally. Due to the deformation based change of the contact interface orientation, the resulting contact F I G U R E 20 Visualization of the computed fluid velocity and pressure, the computed deformation of solid domains, and the interface traction. The color code represents the fluid pressure and the black arrows in domain Ω F indicate the fluid velocity. The blue arrows on Γ represent the FSI traction (case I) and the orange arrows on Γ visualize the contact traction (case II − IV). The visualization of the traction is reconstructed from the nodal interface contributions of (31)- (35) to the overall weak form on the solid mesh. Four points in time are represented from top-left to bottom-right t = 0.003, t = 0.005, t = 0.0065, and t = 0.007. The position of each detailed view in the overall problem is marked in Figure 19 by a red frame force accelerates the solid body in Ω S 2 , and with it the surrounding fluid, in negative y-direction. The fluid pressure in the upper part of the flow domain drops, whereas the pressure in the lower part increases resulting in an almost constant FSI traction acting on Ω S 2 . Finally, at t = 0.007, contact is released and the structural body in Ω S 2 approaches the second barrier. This process leads again to an increased local fluid pressure and thus a growth of the related FSI traction. Due to the structural motion the pressure in the left chamber is raised, which results in a fluid flow out of the fluid chamber. This description of the computed physical process highlights the capabilities of the presented formulation to predict the physical processes in FSCI without requiring a specific treatment whenever topological changes occur.
In Figure 21 (left), the computed flow rates at the inflow boundary and outflow boundary are presented. While the prescribed flow rate at the inflow is constant in time after the initial ramp up phase, the flow rate at the outflow boundary is massively influenced by the overall system. Three phases can be observed where a lower outflow rate (than the inflow rate) is followed by a peak of the flow rate. These can be identified as the phases where the solid domains are compressed due to increasing pressure as Ω S 2 blocks the flow. These phases are always followed by the highly dynamic process of squeezing through. To analyze the overall balance of mass, the transported volume through the inflow-and outflowboundary is given in Figure 21 (right). The difference between the transported volume of outflow and inflow results from the compression or expansion of the solid domains. As no systematical increase of this difference in time can be recognized, no relevant loss in mass occurs. This behavior is expected as discussed in the first presented numerical example in Section 4.1.

CONCLUSION
In this contribution, we presented a consistent formulation to solve general FSCI problems numerically. Topological changes of the fluid domain are enabled by the CutFEM with noninterface fitted discretization. A weak incorporation of F I G U R E 21 Computed flow rates at the inflow boundary Γ in and the outflow boundary Γ out . The normal vector therein is oriented in negative y-direction, which is the main flow direction of the overall configuration (left). Transported volume through the inflow boundary Γ in and the outflow boundary Γ out computed in a post-processing step where an integration in time of the flow rates is performed (right) the governing interface conditions by approaches based on Nitsche's method allows the formulation of a continuous discrete problem even for changing interface conditions. To specify the fluid stress in the region of closed contact, we propose and apply an extension approach. The continuous transition between the "no-slip" and frictionless contact condition in tangential interface orientation is enabled by a general Navier interface condition with a specific law for the slip length. In a first numerical example, the fundamental process in FSCI problems, the contacting and lifting of a convex elastic structure in fluid is analyzed. Therein, the suitability of applying the general Navier interface condition in comparison to a "no-slip" interface condition is evaluated. The positive effect of a skew-symmetric fluid adjoint consistency interface term to the fluid mass conservation is observed and discussed. In two more general examples, the treatment of challenging aspects by the formulation is demonstrated. This includes the representation of large discontinuities of the fluid stress between opposite sides of the structure. Further, large deformation and essential topological changes of the fluid domain as well as large contacting areas are considered.
Still, some aspects for solving general FSCI problems are left to future work. This includes strategies for improving the spatial resolution close to the interface, for examining physically more sophisticated fluid stress representations in the contact zone compared to the simple extension based approach used so far, and for extending the formulation to frictional contact.

ACKNOWLEDGMENT
Open access funding enabled and organized by Projekt DEAL.

APPENDIX A. DISCUSSION OF THE TWO-HALF-PASS NITSCHE CONTACT ALGORITHM
To illustrate the proposed two-half-pass Nitsche contact algorithm, let us have a closer look on the resulting formulation. Consider two elastic bodies with the domains Ω S 1 and Ω S 2 with the potential contact boundaries Γ 1 and Γ 2 in the current configuration. At this point, it is not necessary to include FSI, multi-body, or self-contact. Instead, the discussion is based on a purely structural, two-body contact problem. In that case, (41) reduces to Therein,  S [ u, u] includes the inertia, internal forces, and Neumann boundary conditions of both bodies, see (21). Moreover, the inner product on Γ = Γ 1 ∪ Γ 2 (see Figure 4) includes integration over both potential contact surfaces. With (27) In a continuous problem, in closed contact (i.e., the involved min is unequal to zero), we have that the contacting part of both interfaces Γ 1 and Γ 2 coincide and the two normal vectors are exactly opposing one another: n 1 = −n 2 . In the discrete setting, this may no longer hold exactly. Still, for smooth surfaces in contact, we can still assume n 1 ≈ −n 2 and geometrically overlaying interfaces Γ 1 and Γ 2 , at least in the vicinity of closed contact. Therefore, (A2) can be approximated as wherein integration of the contact integral is performed over Γ 1 only. When choosing S nn = S 1 nn , this is exactly the biased two-body Nitsche contact method discussed in Reference 44 and for a contact stress S nn based on harmonic weights one obtains the method discussed in Reference 45.
The transition from (A3) to (A2) allows for independent decisions of "active" and "inactive" contact on the two contacting surfaces which is crucial in the case of FSCI (see discussion in Section 3.2). This, however, comes at the expense of n 1 ≈ −n 2 holding only approximately. As a consequence, the balance of linear momentum is no longer guaranteed exactly across the contact interface, which is a common drawback of two-half-pass algorithms. 47