Numerical Simulation of Dynamic Fluid-Structure Interaction with Elastic Structure–Rigid Obstacle Contact

An implicit scheme by partitioned procedures is proposed to solve a dynamic fluid–structure interaction problem in the case when the structure displacements are limited by a rigid obstacle. For the fluid equations (Sokes or Navier–Stokes), the fictitious domain method with penalization was used. The equality of the fluid and structure velocities at the interface was obtained using the penalization technique. The surface forces at the fluid–structure interface were computed using the fluid solution in the structure domain. A quadratic optimization problem with linear inequalities constraints was solved to obtain the structure displacements. Numerical results are presented.


Introduction
The Arbitrary Lagrangian Eulerian (ALE) method was successfully employed for solving fluid-structure interaction problems, see [1]. The fluid equations were written in the moving mesh which matches the structure displacement. This method was not adapted for a structure-rigid obstacle contact when the topology of the fluid domain changed from double to simply connected. Some authors have introduced a gap between the elastic structure and the rigid obstacle as in [2]. Some methods exist for solving fluid-structure interaction problems using a fixed mesh for the fluid domain. We recall some of them. The immersed boundary method, see the survey paper [3], was designed originally for a thin structure. It was extended to a thick, viscous hyper-elastic structure in [4], where it was assumed that the fluid and structure densities and viscosities were the same. The fictitious domain with distributed Lagrange multiplier [5] was employed for the simulation of flow around moving rigid bodies. The extension to a visco-elastic structure was proposed in [6]. The densities, respectively, the viscosities of fluid and visco-elastic materials were not the same. For a rigid thick body immersed in an incompressible fluid, the convergence of a penalization method was presented in [7] and the extended finite element method (XFEM) was used in [8]. Nitsche-XFEM was used in [9] for a thin elastic structure immersed in an incompressible fluid.
Concerning the fluid-structure interaction with a structure-rigid obstacle contact, we can cite some works. In [10], for a 1D elastic structure, the Lagrange multiplier was employed to compute the interface forces. The Uzawa algorithm was used to handle the contact. An extension to a 3D nonlinear shell was presented in [11]. An approach using the extended finite element method (XFEM) and a mortar contact formulation was proposed in [12]. In [13], an immersogeometric variational framework for fluid-structure interactions with application to a 3D heart valve was presented. In recent papers, a monolithic Eulerian framework with remeshig was used in [14], a stabilised immersed methodology on hierarchical b-spline grids was employed in [15], the cut finite element method was used in [16], a Nitsche-based formulation with artificial fluid in the gap between structure and obstacle was presented in [17].
In [18], the fictitious domain method with penalization presented in [19,20] was used in order to handle the contact between a linear elastic structure and a rigid obstacle in a fluid-structure interaction problem. The surface forces at the fluid-structure interface were computed using the fluid solution in the structure domain. The equality of the fluid and structure velocities at the interface was obtained using the penalization technique. In the present paper, we present a dynamic fluid-structure interaction problem when the elastic structure is in contact with a rigid obstacle. The fluid was modeled by the Stokes as well as by the Navier-Stokes equations.

Fluid-Structure Interaction with Contact: The Mathematical System
The undeformed structure domain Ω S 0 ⊂ R 2 has Lipschitz boundary ∂Ω S 0 = Γ D ∪ Γ N ∪ Γ C . On Γ D the displacement is zero, on Γ N surface loads are imposed and a subset of Γ C will touch a rigid obstacle, after deformation. In Figure 1, we have Γ D =]MN[, The structure displacement will be denoted by u = (u 1 , A point X = (X 1 , X 2 ) in the undeformed structure domain will occupy the position . The bounded domain of boundary Σ 5 is the rigid obstacle. In the case of an elastic structure-rigid obstacle contact, we have ϕ t (Γ C ) ∩ Σ 5 = ∅. To simplify, we assume that Figure 1. The conclusions are the same if we replace RS by a finite union of Lipschitz arches.
The fluid domain D ⊂ R 2 has the boundary ∂D = ∪ 5 i=1 Σ i and Γ D ⊂ Σ 1 . In Figure 1, D is a rectangle with a hole of boundary Σ 5 . The top side Σ 4 is the inflow, the bottom side Σ 2 is the outflow, the left side Σ 1 and the right side Σ 2 represent the wall.
We suppose that Ω S t ⊂ D and the fluid occupies The set Ω F t is not necessary a Lipschitz domain when the structure touches the rigid obstacle.
The system to be solved is to find the structure displacement u : Ω S 0 × [0, T] → R 2 , the fluid velocity v(·, t) : Ω F t → R 2 and the fluid pressure p(·, t) : Ω F t → R such that: where f S : Ω S 0 × [0, T] → R 2 are the applied volume forces and n S is the unit outward vector normal to ∂Ω S 0 . Additionally, we define f F (·, t) : Ω F t → R 2 and n F (·, t) the unit outward vector normal to ∂Ω F t . The Equations (13)-(15) represent the initial conditions, u 0 , u 1 , v 0 are given.
In (5), g : Σ 4 × [0, T] → R 2 is the imposed velocity. If no-slip boundary conditions are prescribed at ∂D, since the fluid is incompressible, then the volume of the structure is constant, too. In our case, it is not necessary to add the same restriction on the volume of the structure, even when g = 0 because we use (7) at the outflow boundary Σ 2 .
For (12), we followed [21] (Theorem 5.3-1, p. 210), n(ϕ t (X)) is the unit vector normal to Σ 5 at the point ϕ t (X) ∈ Σ 5 , oriented to the exterior of D. We point out that the value of α(X, t) is not given. The meaning of (12) is that the force acting on the structure on the contact zone is parallel to n(ϕ t (X)) and it has opposite direction. In the case of linear elasticity, (12) could be approached by where t S is the unit tangential vector to ∂Ω S 0 . If the elastic structure is not in contact with the obstacle Σ 5 , the Equations (9)-(12) are replaced by : Ω F t → R 4 the Cauchy stress tensors of the structure and fluid, respectively. The structure equations are written using Lagrangian coordinates σ S (u) = λ S (∇ X · u)I + 2µ S X (u), λ S , µ S > 0 are the Lamé coefficients, X (u) = 1 2 ∇ X u + (∇ X u) T and I is the unit matrix. The Eulerian coordinates have been used for fluid equations σ F (v, p) = −p I + 2µ F (v), µ F > 0 is the viscosity and (v) = 1 2 ∇v + (∇v) T . The mass densities are denoted by ρ S > 0 and ρ F > 0. To the best of our knowledge, there are no reported results for fluid-structure interaction with a contact. However, without a structure-obstacle contact, there are reported results for fluid-structure interactions. In general, in the literature, the domain of the fluid is assumed to be Lipschitz. However, when the structure touches the rigid obstacle, the fluid domain is not necessary a Lipschitz domain. Existing results for fluid-structure interactions are presented in [19,20,22,23] and the references are given there. In [19,20,22], the fictitious domain method was used and this technique could be more appropriate to handle the structure-obstacle contact in the fluid-structure interaction framework.

Approximation of the Elastodynamics Frictionless Contact Problem
We analyze, in this section, the linear dynamic elasticity equations with a frictionless contact. In Ω S 0 , volume forces f S are imposed and onΓ N surface loads h S are prescribed. The structure is fixed along Γ D . We recall that ∂Ω S 0 = Γ D ∪ Γ N ∪ Γ C . Let ψ ∈ C 1 (R) be a function describing a part of the top boundary of the obstacle. We set its graph by graph(ψ) = (X 1 , X 2 ) ∈ R 2 , X 2 = ψ(X 1 ) and its epigraph by The non-penetration condition of the elastic structure into the obstacle gives A point X ∈ Γ C belongs to the coincidence set at time instant t if ϕ t (X) ∈ graph(ψ). In the case of linear elasticity, see [24], the condition (16) can be replaced by Denoting by H 1 (Ω S 0 ) the first-order Sobolev space, let us introduce the Hilbert space and the closed, convex and non-empty set We assume f S ∈ L 2 0, T; We can write the linear elastodynamics frictionless contact problem as a variational inequality: The existence for dynamic linear elasticity with frictionless contact in an arbitrary domain is still open, see the monograph ([25] Section 4.1). Existence and uniqueness results for an elastodynamics contact problem were obtained: in [26] for linear and visco-elastic models with Coulomb friction, in [27] for wave problem in a half-space, in [28] for viscoelastic body with frictionless adhesion and in [29] for elastic-visco-plastic equations with Coulomb friction. Several discretization schemes for elastodynamics contact problems have been developed, see the survey papers [30,31]. Let ∆t > 0 be the time step and we note by u n , f S,n , h S,n approximations of u(t n ), f S (t n ), h S (t n ), respectively, for t n = n∆t. In this paper, we use the following implicit time-integration scheme: find u n+1 ∈ K such that with initial conditions u 0 = u 0 and u 1 = u 0 + ∆t u 1 .
With the notation the variational inequality (22) is equivalent to the optimization problem We follow the notations from [18]. Let T S h be a mesh of Ω S 0 of size h, with nvS vertices and ntS triangles. The shape functions φ i : T S h → R associated to vertex A i are obtained by using the finite element P 1 . We set the basis We define the matrixÃ S ∈ R 2 nvS×2 nvS and the vector b n+1 S ∈ R 2 nvS bỹ and The constraint w S ∈ K will be treated weakly. We set the matrix C S ∈ R ngC×2 nvS , where ngC is the number of vertices A i ∈ Γ C and the vector c lb S ∈ R ngC by The optimization problem (25)-(27) has a unique solution, because the cost function is strictly convex and the constraints define a convex set. We set

Approximation of Fluid Equations by Fictitious Domain Method with Penalization
The fluid domain Ω F t can change the topology from double, when no contact occurs, to simply connected, when the structure touches the obstacle. The ALE technique can not by applied in this case. We use the fictitious domain approach and we write the fluid equations in the fixed domain D including Ω F t . The mesh of D is independent on the displacement of the structure.
Let us set the Hilbert spaces and the notations We assume that f F ∈ L 2 0, T; L 2 (D) 2 , g ∈ L 2 0, T; H 1/2 00 (Σ 4 ) 2 and ε > 0 is a penalization parameter. Let χ S n+1 : D → R be the characteristic function The time-integration scheme for the fluid is: for given u n+1 , u n , find the velocity v n+1 The problem has a unique solution, see [19,20] for example in the steady case. As a consequence of the penalization term in (30), we get the weak equality of fluid velocity and structure velocity over the whole structure domain Ω S n+1 and it implies the condition (10) on the fluid-structure interface. In addition, the boundary condition at the outflow (7) is weakly verified.
We employ the stable finite elements P 1 + bubble (we write just P 1 + b) for the velocity and P 1 for the pressure. Let T F h be a mesh of D of size h, with nvF vertices and ntF triangles. We introduce

Computing the Forces at the Fluid-Structure Interface Using the Fictitious Domain
We recall the notations: . We assume in this section that Ω S n+1 and Ω F n+1 are Lipschitz domains. If no contact occurs, this assumtion is true if the displacement is Lipschitz, see [23].
The fluid-structure interface is defined by In the case when no contact occurs, the interface is Γ FS n+1 = Γ n+1 , where Γ n+1 = ϕ n+1 (Γ 0 ), but when the structure touches the rigid obstacle, the interface is Γ FS n+1 = ϕ n+1 (Γ 0 ) \ RS. For a related problem, called a thin obstacle problem, it is proved in [32] that the coincidence set is a finite union of closed disjoint intervals. The conclusions of this section are the same if we replace RS by a finite union of Lipschitz arches. In the case of a thick obstacle problem, we can find results on the topology of the coincidence set in [33,34]. For example, let u : Ω → R be an elastic membrane. If Ω ⊂ R 2 is strictly convex and the obstacle ψ : Ω → R is strictly concave, it is possible to prove that the coincidence set is a simply connected domain (without holes), see [33] (Theorem 6.2, p. 176).
We employ the notations v n+1,F ε for the restriction of v n+1 2 and similarly for p n+1,F ε and p n+1,S ε if p n+1 ε ∈ Q, w n+1,F and w n+1,S if w ∈ W. Using w with compact support in Ω F n+1 in (30), we get , see [35] (Chapter VII, p. 1241) and [36] (p. 325), defined by for all w n+1, Similary, we define for all w n+1,S ∈ H 1 (Ω S n+1 ) 2 such that w n+1,S = 0 on ∂Ω S n+1 ∩ ∂D. From (30), we get We have (11) at the fluid-structure interface, or equivalently Then, the surface forces from the fluid acting on the structure are computed by using the right-hand side of (35). As in [19,20], by using the change of variable formula, we get The right-hand side of the above equality, in the case of linear elasticity, could be approached by where v n+1

Remark 1.
The five terms of (37) represent the surface forces from the fluid acting on the structure. This formula holds in the case of contact as well as if no contact occurs.

Time Integration Scheme by Fixed Point Iterations
The fixed point algorithm is a well known method for solving a fluid-structure interaction without contact. Additionally, it was successfully used in [18] when the structure comes into contact with a rigid obstacle, in the steady case.
Algorithm from Time Instant n to n + 1 We assume that the structure displacements u n h , u n−1 h ∈ K h as well as the fluid velocity v n ε,h ∈ W h are given. We look for u n+1 ε,h and p k,n+1 ε,h , respectively, for k ∈ N. For the stopping test, we use the parameter tol > 0 and let ω ∈ (0, 1) be the relaxation parameter.
Step 1. Set the initial displacement of the structure u 0,n+1 h = u n h and put k = 0.
Step 3. Set v k+1,n+1 The last five integrals represent the surface forces from the fluid acting on the structure as in the end of Section 5.

Remark 2.
All the results can be extended to the Navier-Stokes equation for the fluid: • on the left-hand side of (3), add ρ F (v · ∇)v, • on the left-hand side of (30) and Section 5, • on the left-hand side of (38), add c F (v k+1,n+1

Numerical Results
We use the software FreeFem++ [37] for the numerical tests. The linear system at the Step 2 will be solved by the library "MUltifrontal Massively Parallel sparse direct Solver" (MUMPS) and the constrained optimization problem at the Step 4 will be solved by the library "Interior Point OPTimizer" (IPOPT) implemented in FreeFem++.

Dynamic Fluid-Structure Interaction when the Displacements Are Limited by a Rigid Disk
We use a dynamic version of the example from [18]. The geometrical configuration is: We set for the structure: ρ S = 1100 Kg/m 3 , E S = 3 × 10 5 N/m 2 , ν S = 0.3, the mass density, Young's modulus, Poisson's ratio, respectively and the applied volume forces on the structure f S : For the fluid, we use: ρ F = 1000 Kg/m 3 , µ F = 0.0035 N · s/m 2 , the mass density, dynamic viscosity, respectively, and the applied volume forces on the fluid f F : In the fluid domain, we denote by Σ 1 , Σ 2 , Σ 3 , Σ 4 the left, bottom, right, top boundaries, respectively, and by Σ 5 the boundary of the circular obstacle B(0.2, −0.07; 0.05).
The triangular finite elements P 1 +bubble and P 1 have been used for the approximation of the fluid velocity and pressure and the finite element P 1 was used for the structure problem. We use a fluid mesh of 37,124 triangles and 18,885 vertices, fluid mesh size h F = 0.005, a structure mesh of 768 triangles and 451 vertices, structure mesh size h S = 0.005, the time step ∆t = 0.1 s and the number of time steps N = 650.
The penalization parameter is ε = 10 −5 , the relaxation parameter is set to ω = 0.1 and for the stopping test tol = 10 −3 × area(Ω S 0 ) = 9 × 10 −6 . At Step 2, we solve a linear system of dimension 93,749 and at Step 4, we solve a constrained optimization problem of dimension 902, subject to 59 inequalities corresponding to (26) and 14 components are fixed to zero corresponding to (27).
In Figure 2, it can be seen that the slope of the curves was zero when the structure was in contact with the obstacle. The structure traveled down and made contact with the rigid disk from 9 s to 23 s, then it ascended. At 32 s, the position was similar to the undeformed shape, it continued to ascend up to maximal position at 43 s and then descended. At each time step, we solve the fluid-structure interaction problem iteratively using the fixed point method. At each fixed point iteration, one fluid sub-problem and one structure sub-problem are solved. The number of iterations at each time step gives an indicator of the computational effort. In this simulation, the number of fixed point iterations was less than 5, except for five time steps, where the number was 7, 10 and for three time steps it was 50. Slow convergence of the fixed point method can be explained when the contraction constant of the map was close to 1. Other methods for implicitly solving the fluid-structure interaction problem are Newton or quasi-Newton methods, see [38] for a comparative study.
In Figure 3, we show the computed solution at time instant t = 20 s when the structure starts to touch the obstacle. The fluid velocity was maximal near the right end of the structure. In the zone limited by the left side of the fluid domain, the bottom side of the structure and the obstacle, the fluid velocity was very small.

Numerical Simulation of Valve Dynamics
The undeformed structure domain Ω S 0 is a quarter of a ring with interior radius R = 1 cm, thickness H = 0.03 cm and center (R + H, H 1 ), where H 1 = 0.95 cm, see Figure 4. Its boundaries are as follows:  The mechanical properties of the elastic structure are the following: mass density ρ S = 1.1 g/cm 3 , Young's modulus E S = 4 × 10 5 dyne/cm 2 , Poisson's ratio ν S = 0.4, applied volume forces on the structure f S : Ω S 0 × [0, T] → R 2 , f S = (0, 0) dyne/cm 3 . The parameters for the fluid are: mass density ρ F = 1 g/cm 3 , dynamic viscosity µ F = 0.035 dyne · s/cm 2 , applied volume forces on the fluid f F : D × [0, T] → R 2 , f F = (0, 0) dyne/cm 3 . The inflow velocity profile is g = (g 1 , g 2 ) : Σ 1 × [0, T] → R 2 where g 2 = 0 and with V max = 10 cm/s or 50 cm/s, T 1 = 0.2 s and T = 0.5 s. The boundary conditions for the fluid are as follows: Figure 5. The initial conditions are as follows: u 1 = 0, v 0 = 0. Since, the undeformed structure domain does not satisfy the non-penetration constraint, we solved a steady elasticity contact problem to obtain u 0 . We used a fluid mesh of 13,090 triangles and 6718 vertices, fluid mesh size h F = 0.02, a structure mesh of 40 triangles and 42 vertices, structure mesh size h S = 0.078, the time step ∆t = 0.005 s and the number of time steps N = 100. The penalization parameter was ε = 10 −3 , the relaxation parameter was set to ω = 0.1 and for the stopping test tol = 10 −6 . At Step 2, we solved a linear system of dimension 33,244 and at Step 4, we solved a constrained optimization problem of dimension 84, subject to 9 inequalities corresponding to (26) and 4 components were fixed to zero corresponding to (27).
We observe in Figure 6, the number of fixed point iterations at each time step, which proves the efficiency of the method. The structure velocity was near to zero, when the opening of the valve was maximal or when the valve was closed, then only 1-2 iterations were necessary to achieve the convergence. The total computational time was 323 s for V max = 10 and 584 s for V max = 50. In Figure 7, we plot the horizontal and vertical displacement of the structure tail. We observe in this figure that the valve was closed after about t = 0.230 for V max = 10 and after about t = 0.370 for V max = 50. At the initial time, the valve was closed and prestressed, since the undeformed structure domain did not satisfy the non-penetration constraint. It ascended to its highest opening position, see Figure 8 the top images, then the valve descended and made contact with the rigid foundation, see Figure 8 the bottom images, for V max = 10.
In Figure 9, Von Mises stress distribution in the structure is shown for V max = 10. We observe that Von Mises stress was larger near to the top fixed boundary even when the valve was in contact with the bottom wall. The reason for this is that the structure was prestressed when the valve was closed. Similar behavior we observe for V max = 50, in Figures 10 and 11.
We also observe a vortex downstream of the closed valve, see also Figure 12, where the vorticity is plotted.     The structure deformations are important when V max = 50 and a non-linear model for the elasticity equations could be used.

Check Valve Interacting with Navier-Stokes Flow
As we seen in Remark 2, the extension from Stokes to Navier-Stokes equations for the fluid is straightforward. Now, we present the numerical results for the check valve problem adapted from [15]. In our case, the fluid is governed by the Navier-Stokes equations but the structure is still governed by the linear elasticity equations.
All the units are derived from centimeter (cm), gram (g), millisecond (ms). The fixed computational fluid domain D was the union of three rectangles: ]0, 2 . The mass density of the fluid was ρ F = 0.8 g/cm 3 , the dynamic viscosity was µ F = 0.005 g/(cm · ms) and the boundary conditions were as follows: p in (t) = 10 sin 2πt 100 g/(cm · (ms) 2 ) prescribed surface stress at the bottom boundary, do-nothing at the top boundary and non-slip boundary condition otherwise. There were no applied volume forces in fluid. The undeformed structure domain was a rectangle of thickness H = 0.05 cm and length L = 0.9 cm and the coordinates of the left bottom corner were (0.3, 2.6). The mass density, Young's modulus, Poisson's ratio of the structure were as follows: ρ S = 8 g/cm 3 , E S = 1 × 10 5 g/(cm · (ms) 2 ), ν S = 0.3, respectively. There were no applied volume forces in structure.
In the paper [15], the structure is modeled with Neo-Hookean material of thickness H = 0.025 cm and E S = 1 × 10 4 g/(cm · (ms) 2 ). Our technique to compute the forces at the fluid-structure interface using integral over the structure domain requires a structure that is not too thin. For this reason, we used a structure of thickness H = 0.05 cm. A linear model for the structure with E S = 1 × 10 4 g/(cm · (ms) 2 ) gives important displacements, for this reason we worked with E S = 1 × 10 5 g/(cm · (ms) 2 ).
The fluid mesh had 34,432 triangles and 17,602 vertices, the structure mesh had 156 triangles and 121 vertices, the time step was ∆ t = 0.05 ms. The penalization parameter ε, the relaxation parameter ω, the parameter for the stopping test was tol and the finite elements were the same as in the previous simulations.
The contour plots of the velocity and pressure in Figure 13 are very similar to [15]. The maximal velocity magnitude 4.8 cm/ms is comparable with 42.75 mm/ms = 4.275 cm/ms in [15] and the minimal pressure was obtained in the right top corner of the middle rectangle of the fluid domain: −4.4 g/(cm · (ms) 2 ) in our test and −0.5 g/(mm · (ms) 2 ) = −5 g/(cm · (ms) 2 ) in [15]. The vertical displacement of the structure tail in Figure 14 is correlated to the surface stress at the inflow (41): the valve ascended to its highest position t = 25 ms, then descended and made contact with the rigid obstacle after t = 51 ms. The valve was closed in the time interval t ∈ [51, 100] ms. One difference from the report of [15], was that in our test, the shape of the vertical displacement of the structure tail did not have two humps in the time interval t ∈ [0, 50] ms.

Conclusions
We have presented a numerical method for solving by partitioned procedures the dynamic behavior of an elastic thick structure immersed in an incompressible fluid. This method handles the contact between a linear elastic structure and a rigid obstacle. For the fluid flow, the fictitious domain method with penalization was used. The surface forces at the fluid-structure interface were computed using the fluid solution in the structure domain. Numerical tests were presented.