Numerical procedure for fluid-structure interaction with structure displacements limited by a rigid obstacle

A fixed point algorithm is proposed to solve a fluid-structure interaction problem with the supplementary constraint that the structure displacements are limited by a rigid obstacle. Fictitious domain approach with penalization is used for the fluid equations. The surface forces from the fluid acting on the structure are computed using the fluid solution in the structure domain. The continuity of the fluid and structure velocities is imposed through the penalization parameter. The constraint of non-penetration of the elastic structure into the rigid obstacle is treated weakly. A convex constrained optimization problem is solved in order to get the structure displacements. Numerical results are presented. c © 2017 University of West Bohemia. All rights reserved.


Introduction
Fluid-structure interaction is encountered in diverse medical and engineering applications: blood flow in artery [2,6], motion of pharmaceutical liquid capsule enclosed by elastic membrane in flow [5], impact of tsunami wave on coastal protection [10].
We study the behavior of an elastic structure immersed in a viscous incompressible fluid with the supplementary constraint that the structure displacements are limited by a rigid obstacle. We use Stokes equation to model the flow motion and the displacement of the structure will be governed by linear elasticity equations. This choice is in order to simplify the computational complexity. A perspective is the extension towards to the Navier-Stokes equations.
The Arbitrary Lagrangian Eulerian (ALE) formulation consisting in using moving mesh which follow the structure displacement was successfully used for solving fluid-structure interaction problem without contact. But the ALE method is not adapted when the structure is in contact with a rigid obstacle and the topology of the fluid domain changes. In [13], an algorithm for protecting the quality of the fluid mesh in the gap between the elastic structure and the rigid obstacle is employed.
Other methods enter in the category of fixed mesh. In [4], for a 1D elastic structure, surface Lagrange multiplier was used in order to compute the forces from the fluid acting on the structure and to impose the kinematic boundary condition at the interface. The contact is handled by the Uzawa algorithm. An extension to 3D nonlinear shell which manages non-convex contact constraint is proposed in [1]. An approach combining an extended finite element method (XFEM) and a dual finite deformation mortar contact formulation is proposed in [12].
In this paper, we have adapted the fictitious domain method with penalization presented in [7] and [8] in order to handle the contact. The surface forces from the fluid acting on the structure are computed using the fluid solution in the structure domain. The continuity of the fluid and structure velocities is imposed through the penalization parameter. We treat weakly the constraint of non-penetration of the elastic structure into the rigid obstacle. To manage the contact in linear elasticity, we solve a convex constrained optimization problem, not a variational inequality.

Setting for fluid-structure interaction with contact
Let Ω S 0 ⊂ R 2 be the undeformed structure domain and we assume that it is an open, bounded and connected subset. We suppose that its boundary is Lipschitz and admits the decomposition are relatively open subsets, mutually disjoint. The structure is fixed along the portion Γ D , it is subjected to volume forces applied in Ω S 0 and to surface loads on Γ N . After deformation, a portion of Γ C will be in contact with a rigid obstacle. This portion is unknown beforehand. For a better understanding, we consider for the undeformed structure domain the rectangle Ω  Suppose that the structure is elastic and denote by u = (u 1 , u 2 ) : Ω S 0 → R 2 its displacement. A particle of the structure whose initial position was the point X = (X 1 , X 2 ) will occupy the position With the notations ϕ(P ) = P and ϕ(Q) = Q , we have ϕ (Γ N ) = P Q ∪ Q M, where P Q is a curve of ends P and Q . The disk of boundary Σ 5 is a rigid obstacle. When the structure is in contact with the obstacle, then Fig. 1. We point out that the rigid obstacle is not a part of Ω S u . Let D ⊂ R 2 be a bounded open domain with Lipschitz boundary such that ∂D = ∪ 5 i=1 Σ i and Γ D ⊂ Σ 1 . In Fig. 2, D is a rectangle with a circular hole of boundary Σ 5 and Σ 1 , Σ 2 , Σ 3 , Σ 4 are the left, bottom, right, top sides, respectively. We admit that Ω S u ⊂ D for all admissible displacements and the fluid occupies To resume, the contact zone between and the common boundary of the fluid and obstacle is Σ 5 \ RS. We point out that the set Ω F u is not a Lipschitz domain when the structure is in contact with the rigid obstacle. Also it can change the topology, before contact Ω F u is double connected and it could be simply connected after contact, as in Fig. 2.
The fluid equations are described using Eulerian coordinates, while for the structure equations, the Lagrangian coordinates are employed. The gradients with respect to the Eulerian coordinates x ∈ Ω S u of a scalar field q or a vector field w are denoted by ∇q, ∇w. The divergence operators with respect to the Eulerian coordinates of a vector field w and of a tensor σ are denoted by ∇ · w and ∇ · σ. Similarly, when the derivatives are with respect to the Lagrangian coordinates X = ϕ −1 (x) ∈ Ω S 0 , we use the notations: ∇ X u, ∇ X · u, ∇ X · σ. We assume that the structure verifies the linear elasticity equation. The stress tensor of the structure written in the Lagrangian framework is σ S (u) = λ S (∇ X · u) I + 2μ S X (u), where λ S , μ S > 0 are the Lamé coefficients, X (u) = 1 2 ∇ X u + (∇ X u) T and I is the unit matrix. We assume that the fluid is governed by the Stokes equations, then the Cauchy stress tensor is given by

Strong formulation
The problem is to find the structure displacement u : Ω S 0 → R 2 , the fluid velocity v : Ω F u → R 2 and the fluid pressure p : Ω F u → R such that: where f S : Ω S 0 → R 2 are the applied volume forces on the structure and n S is the structure unit outward vector normal to ∂Ω S 0 . Similarly, we define f F : Ω F u → R 2 and n F the fluid unit outward vector normal to ∂Ω F u . In (5), g : Σ 4 → R 2 is a prescribed velocity. For (12), we can consult [3], Theorem 5.3-1, p. 210. In (12), n(ϕ(X)) is the unit vector normal to Σ 5 at the point ϕ(X) ∈ Σ 5 , oriented to the exterior of D.
We have denoted by σ S (u) : the Cauchy stress tensors of the structure and fluid, respectively. We point out that the stress tensor of the structure is defined on the undeformed structure domain Ω S 0 , while the Cauchy stress tensors of the fluid is defined in the deformed domain Ω F u . In view of (9), non-slip boundary condition is imposed to the flow velocity at the fluid-structure interface Γ u \ RS as well as at the common boundary of fluid and obstacle Σ 5 \ RS. If the structure doesn't come into contact with Σ 5 the boundary of rigid obstacle, the equations (9)-(12) are replaced by

Contact problem without friction in linear elasticity
First, we recall some definitions and notations. Let Ω be a bounded, connected, open set of R 2 with Lipschitz boundary ∂Ω and let Γ be an open subset of ∂Ω with positive measure. We denote by L 2 (Ω) the class of all measurable functions v : Ω → R such that Ω v 2 (x) dx < ∞. As usual, we denote the first order Sobolev space by The space of trace of functions of . Now, we pay attention to the elastic structure in contact with a rigid foundation modeled as the graph of a function, as in Fig. 3. The body Ω S 0 is subjected to volume forces f S applied in Ω S 0 and to surface loads h S on Γ N . We recall that ∂Ω S 0 = Γ D ∪ Γ N ∪ Γ C . Let ψ ∈ C 1 (R) be the rigid foundation and we denote its graph by and its epigraph by We assume that the undeformed structure domain is in the interior of epi(ψ). Since the structure is fixed along Γ D , then The foundation is rigid, then the constraint of non-penetration of the elastic structure into the obstacle gives In the case of small displacements, we can approach and the condition (14) will be replaced by We Let us introduce the Hilbert space Following [11], Chapter 5, it is possible to define the order ≤ for the functions in Sobolev space H 1/2 (Γ C ), consequently we can define the set which is closed, convex and non empty. The problem is to find u ∈ K ψ solution of the constrained optimization problem It is known, see for example [11], that the problem has a unique solution and verifies the variational inequality Let T S h be a triangulation of Ω S 0 of size h, with nvS vertices and ntS triangles. We construct the shape functions φ i : T S h → R associated to vertex A i , which are piecewise linear function, globally continuous. For the two-dimension displacements, we introduce the basis We define the matrix A S ∈ R 2nvS×2nvS and the vector b S ∈ R 2nvS by The constraint w S ∈ K ψ will be treated weakly. So that, we introduce the matrix C S ∈ R ngC×2nvS , where ngC is the number of vertex A i ∈ Γ C and the vector c lb S ∈ R ngC by In order to impose (13), we set ξ ub , ξ lb ∈ R 2nvS The optimization problem (20)-(22) will be solved by the library IPOPT implemented in the software freefem++, see [9]. We set

Stokes equations using fictitious domain approach with penalization
The set occupied by the fluid Ω F u is not a Lipschitz domain when the structure is in contact with the rigid obstacle and it can change the topology from double to simply connected. The ALE framework can not by applied. We use the fictitious domain approach consisting to write the fluid equations in the fixed domain D which includes Ω F u . The mesh of D does not depend on the displacement u, it is fixed.
Let us introduce the Hilbert spaces and the notations We assume that f F ∈ (L 2 (D)) 2 , g ∈ H 1/2 00 (Σ 4 ) 2 and ε > 0 is a penalization parameter. Let χ S u : D → R be the characteristic function For given u, the fluid problem is to find: the velocity v ∈ (H 1 (D)) This problem has a unique solution, see [7,8] for example. Let T be a triangle of vertices A, B, C of coordinates ( . We denote by G its barycenter, i.e. the point of coordinates we introduce the barycentric coordinates λ 1 (x), λ 2 (x), λ 3 (x), as the solution of the linear system ⎛ ⎝ We define the bubble function b T : T → R by b T (x) = λ 1 (x)λ 2 (x)λ 3 (x) and the spaces of functions The finite element denoted by P 1 + bubble or P 1 + b is the triplet and the finite element P 1 is the triplet (T, {A, B, C}, P 1 (T )).
Let T F h be a triangulation of D of size h, with nvF vertices and ntF triangles. We introduce

Computing the surface forces from the fluid acting on the structure
Supposing for instant that structure is not in contact with the obstacle Σ 5 . We employ the notations v + ε for the restriction of v ε to Ω F u , v − ε for the restriction of v ε to Ω S u and similarly for p + ε and p − ε . Using w with compact support in Ω F u in (25), we get and using w with compact support in Ω S u , we get Even if v + ε and p + ε are not smooth, we can give a sense of σ F (v + ε , p + ε ) n F on Γ u by the element From (25), we get that which at least formally means But at the fluid-structure interface, we have the boundary conditions (10), (11) or Finally, we can compute the surface forces from the fluid acting on the structure using the right-hand side of (30). As in [7,8], by using the change of variable formula, we get where ϕ(X) = X + u(X), F(X) = I + ∇ X u(X), J(X) = det F(X), w S = w − • ϕ ∈ W S . In the case of small displacement, the right-hand side of the above equality could be approached by where v ε = v ε • ϕ and p ε = p ε • ϕ.

Fixed point iterations
The fixed point algorithm was successful used for fluid-structure interaction without contact. We adapt this method in the case where the structure comes in contact with a rigid obstacle. Let tol > 0 be a parameter for the stopping test, ω ∈ (0, 1) the relaxation parameter and k the iterations counter.

Algorithm
Step 1. Choose an initial displacement of the structure u 0 h ∈ K ψ,h and put k := 0.
Step 2. Find the velocity v k+1 ε,h ∈ g + W h and the pressure p k+1 ε,h ∈ Q h by solving the discrete fluid problem Step 3. where The last four terms in the right-hand side of (33) represent the approximate surface forces from the fluid acting on the structure.

Linear elasticity with frictionless contact
The rigid foundation is given by ψ : The parameters for the fluid are: dynamic viscosity μ F = 0.0035 N · s/m 2 , the applied volume forces on the fluid f F : D → R 2 , f F = (0, −9.81 × 100) N/m 3 . The inflow velocity profile is g = (g 1 , g 2 ) where The considered boundary conditions for the fluid are: Fig. 2 for the notations of the boundaries. For the approximation of the fluid velocity and pressure we have employed the triangular Table 1. Numerical parameters for fluid-structure interaction with obstacle finite elements P 1 +bubble and P 1 . The finite element P 1 was used in order to solve the structure problem. The characteristic function was approached by P 0 , the piecewise constant discontinuous finite element. In Table 1, k is the number of iterations in order to achieve the stopping test with tol = 10 −4 . For ε = 10 −6 , the algorithm does not satisfy the stopping test. We remark that in the deformed structure domain, the fluid velocity is almost zero, see Table 1, column v ε,h 0,Ω S u . For a fluidstructure interaction without obstacle, it is proved that v ε 0,Ω S u /ε is bounded independent of ε, see [8], inequality (5.10). We report the computed values in the last column of Table 1. Fig. 9 shows computed solutions for different values of the parameters: ε > 0 penalization parameter, ω ∈ (0, 1) relaxation parameter, h F size of fluid mesh, h S size of structure mesh.

Conclusions
A partitioned procedure based on fixed point iterations is proposed to solve a fluid-structure interaction problem with contact. Employing fictitious domain approach with penalization, it is possible to compute the surface forces from the fluid acting on the structure by using the fluid solution in the structure domain. We treat weakly the constraint of non-penetration of the elastic structure into the rigid obstacle. Solving the contact elasticity problem consists in finding the solution of a convex constrained optimization problem. In a future work, we extend the actual results to the dynamic case: interaction between unsteady Navier-Stokes equations and linear elastodynamics with contact.