Next Article in Journal
Physics-Capturing Discretizations for Spectral Wind-Wave Models
Next Article in Special Issue
Effect of Strain Measurement Layout on Damage Detection and Localization in a Free Falling Compliant Cylinder Impacting a Water Surface
Previous Article in Journal
Spectral Early-Warning Signals for Sudden Changes in Time-Dependent Flow Patterns
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Othman Yakhlef
1 and
Cornel Marius Murea
2,*
1
Laboratoires Analyse, Optimisation et Traitement de l’Information, Faculté des Sciences Exactes et Informatique, Université Mohammed Seddik BenYahia, BP 98 Ouled Aissa, Jijel 18000, Algeria
2
Département de Mathématiques, Institut de Recherche en Informatique, Mathématiques, Automatique et Signal, Université de Haute Alsace, 6, Rue des Frères Lumière, 68093 Mulhouse, France
*
Author to whom correspondence should be addressed.
Submission received: 8 November 2020 / Revised: 15 January 2021 / Accepted: 20 January 2021 / Published: 22 January 2021
(This article belongs to the Special Issue Fluid Structure Interaction: Methods and Applications)

Abstract

:
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.

1. 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.

2. Fluid–Structure Interaction with Contact: The Mathematical System

The undeformed structure domain Ω 0 S R 2 has Lipschitz boundary Ω 0 S = Γ ¯ 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 = ] M N [ , Γ C = N P , Γ N = P M .
The structure displacement will be denoted by u = u 1 , u 2 : Ω ¯ 0 S × [ 0 , T ] R 2 . A point X = ( X 1 , X 2 ) in the undeformed structure domain will occupy the position x = φ t X = X + u X , t in the deformed structure domain Ω t S = φ t Ω 0 S . We set Γ 0 = Ω 0 S Γ ¯ D and Γ t = φ t ( Γ 0 ) = φ t ( Γ N ) φ t ( Γ C ) .
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 R S = φ t Γ C Σ 5 is a connected Lipschitz curve. We have φ t Γ C = N R R S S P t , where P t = φ t ( P ) , see Figure 1. The conclusions are the same if we replace R S by a finite union of Lipschitz arches.
The fluid domain D R 2 has the boundary D = i = 1 5 Σ ¯ 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 Ω t S D and the fluid occupies Ω t F = D Ω ¯ t S . The set Ω t F 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 : Ω ¯ 0 S × [ 0 , T ] R 2 , the fluid velocity v ( · , t ) : Ω ¯ t F R 2 and the fluid pressure p ( · , t ) : Ω ¯ t F R such that:
ρ S 2 u t 2 X · σ S u = f S , in Ω 0 S × ] 0 , T [
u = 0 , on Γ D × ] 0 , T [
ρ F v t · σ F v , p = f F , in Ω t F , t ] 0 , T [
· v = 0 , in Ω t F , t ] 0 , T [
v = g , on Σ 4 × ] 0 , T [
v = 0 , on ( Σ 1 Γ D ) Σ 3 × ] 0 , T [
σ F v , p n F = 0 , on Σ 2 × ] 0 , T [
Γ t = φ t ( Γ 0 ) D ¯ , t ] 0 , T [
v = 0 , on Σ 5 { x Σ 5 ; X Γ C , x = φ t ( X ) } , t ] 0 , T [
u t ( X , t ) = v ( φ t ( X ) , t ) , on Γ N { X Γ C ; φ t ( X ) Σ 5 } , t ] 0 , T [ σ S u ( X , t ) n S ( X ) = σ F v , p n F ( φ t ( X ) , t ) ,
on Γ N { X Γ C ; φ t ( X ) Σ 5 } , t ] 0 , T [ σ S u ( X , t ) n S ( X ) = α ( X , t ) n ( φ t ( X ) ) , α ( X , t ) 0 ,
on { X Γ C ; φ t ( X ) Σ 5 } , t ] 0 , T [
u ( X , 0 ) = u 0 ( X ) , X Ω 0 S
u t ( X , 0 ) = u 1 ( X ) , X Ω 0 S
v ( x , 0 ) = v 0 ( x ) , x Ω 0 F
where f S : Ω ¯ 0 S × [ 0 , T ] R 2 are the applied volume forces and n S is the unit outward vector normal to Ω 0 S . Additionally, we define f F ( · , t ) : Ω ¯ t F R 2 and n F ( · , t ) the unit outward vector normal to Ω t F . 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
n S · σ S n S 0 , t S · σ S n S = 0
where t S is the unit tangential vector to Ω 0 S .
If the elastic structure is not in contact with the obstacle Σ 5 , the Equations (9)–(12) are replaced by
v = 0 , on Σ 5 , × ] 0 , T [ u t ( X , t ) = v ( φ t ( X ) , t ) , on ( Γ N Γ C ) × ] 0 , T [ σ S u ( X , t ) n S ( X ) = σ F v , p n F ( φ t ( X ) , t ) , on ( Γ N Γ C ) × ] 0 , T [ .
We have denoted by σ S u : Ω ¯ 0 S R 4 and σ F v , p : Ω ¯ t F 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.

3. Approximation of the Elastodynamics Frictionless Contact Problem

We analyze, in this section, the linear dynamic elasticity equations with a frictionless contact. In Ω 0 S , volume forces f S are imposed and on Γ N surface loads h S are prescribed. The structure is fixed along Γ D . We recall that Ω 0 S = Γ ¯ 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
g r a p h ( ψ ) = ( X 1 , X 2 ) R 2 , X 2 = ψ ( X 1 )
and its epigraph by
e p i ( ψ ) = ( X 1 , X 2 ) R 2 , X 2 ψ ( X 1 ) .
The non-penetration condition of the elastic structure into the obstacle gives
φ t ( Γ C ) e p i ( ψ ) , t ] 0 , T [ .
A point X Γ C belongs to the coincidence set at time instant t if φ t ( X ) g r a p h ( ψ ) . In the case of linear elasticity, see [24], the condition (16) can be replaced by
ψ ( X 1 ) + ψ ( X 1 ) u 1 ( X 1 , X 2 , t ) X 2 + u 2 ( X 1 , X 2 , t ) , ( X 1 , X 2 ) Γ C , t ] 0 , T [ .
Denoting by H 1 ( Ω 0 S ) the first-order Sobolev space, let us introduce the Hilbert space
W S = w S H 1 Ω 0 S 2 ; w S = 0   on   Γ D ,
the bi-linear form a S : W S × W S R ,
a S u , w S = Ω 0 S λ S X · u X · w S + 2 μ S ϵ X u : ϵ X w S d X
and the closed, convex and non-empty set
K = w S = ( w 1 S , w 2 S ) W S ; ψ ( X 1 ) w 1 S ( X 1 , X 2 ) w 2 S ( X 1 , X 2 ) X 2 ψ ( X 1 ) , ( X 1 , X 2 ) Γ C .
We assume f S L 2 0 , T ; L 2 ( Ω 0 S ) 2 , h S L 2 0 , T ; L 2 ( Γ N ) 2 , u 0 K , u 1 L 2 ( Ω 0 S ) 2 . We can write the linear elastodynamics frictionless contact problem as a variational inequality: find u C [ 0 , T ] ; K C 1 [ 0 , T ] ; L 2 ( Ω 0 S ) 2 such that
ρ S Ω 0 S d 2 u ( t ) d t 2 · w S u ( t ) d X + a S u ( t ) , w S u ( t ) Ω 0 S f S ( t ) · w S u ( t ) d X
+ Γ N h S ( t ) · w S u ( t ) d s , w S K , almost everywhere t ] 0 , T [
u ( 0 ) = u 0
d u d t ( 0 ) = u 1 .
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 visco-elastic 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
ρ S Ω 0 S u n + 1 2 u n + u n 1 Δ t 2 · ( w S u n + 1 ) d X + a S ( u n + 1 , w S u n + 1 ) Ω 0 S f S , n + 1 · w S u n + 1 d X + Γ N h S , n + 1 · w S u n + 1 d s , w S K
with initial conditions u 0 = u 0 and u 1 = u 0 + Δ t u 1 .
With the notation
L S n + 1 ( w S ) = Ω 0 S f S , n + 1 · w S d X + Γ N h S , n + 1 · w S d s ,
the variational inequality (22) is equivalent to the optimization problem
u n + 1 = arg inf w S K ρ S 2 Ω 0 S w S · w S Δ t 2 d X + ρ S Ω 0 S 2 u n + u n 1 Δ t 2 · w S d X + 1 2 a S ( w S , w S ) L S n + 1 ( w S ) .
We follow the notations from [18]. Let T h S be a mesh of Ω 0 S of size h, with n v S vertices and n t S triangles. The shape functions ϕ i : T h S R associated to vertex A i are obtained by using the finite element P 1 . We set the basis ϕ i = ( ϕ 1 i , ϕ 2 i ) : T h S R 2 for i = 1 , , 2 n v S defined by
ϕ i = ( ϕ i , 0 ) , for i = 1 , , n v S and ϕ n v S + i = ( 0 , ϕ i ) , for i = 1 , , n v S .
We define the matrix A ˜ S R 2 n v S × 2 n v S and the vector b S n + 1 R 2 n v S by
A ˜ S = ( a ˜ i j S ) , a ˜ i j S = a S ( ϕ j , ϕ i ) + ρ S Δ t 2 Ω 0 S ϕ j · ϕ i d X , i , j = 1 , , 2 n v S
and
( b S n + 1 ) i = L S n + 1 ( ϕ i ) + ρ S Δ t 2 Ω 0 S ( 2 u n u n 1 ) · ϕ i d X , i = 1 , , 2 n v S .
The constraint w S K will be treated weakly. We set the matrix C S R n g C × 2 n v S , where n g C is the number of vertices A i Γ C and the vector c S l b R n g C by
C S = ( c i j S ) , c i j S = Γ C ( ψ ( X 1 ) ϕ 1 j ( X 1 , X 2 ) + ϕ 2 j ( X 1 , X 2 ) ) ϕ i ( X 1 , X 2 ) d s
for j = 1 , , 2 n v S and A i Γ C and
( c S l b ) i = Γ C ( ψ ( X 1 ) X 2 ) ϕ i ( X 1 , X 2 ) d s
for A i Γ C .
In order to impose (2), we set ξ u b , ξ l b R 2 n v S
ξ i u b = ξ n v S + i u b = ξ i l b = ξ n v S + i l b = 0 , if A i Γ D
otherwise ξ i u b = and ξ i l b = . The discrete form of (23) is u h n + 1 = i = 1 2 n v S ξ i ϕ i where
ξ = arg inf ξ R 2 n v S 1 2 A ˜ S ξ , ξ b S n + 1 , ξ
C S ξ c S l b
ξ u b ξ ξ l b .
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
W h S = w h S = i = 1 2 n v S ξ i ϕ i ; such that ξ verify ξ u b ξ ξ l b ,
K h = w h S = i = 1 2 n v S ξ i ϕ i W h S ; such that ξ verify C S ξ c S l b .

4. Approximation of Fluid Equations by Fictitious Domain Method with Penalization

The fluid domain Ω t F 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 Ω t F . The mesh of D is independent on the displacement of the structure.
Let us set the Hilbert spaces
W = w H 1 D 2 , w = 0 on D Σ 2 Q = L 2 D
and the notations
a F : H 1 D 2 × H 1 D 2 R , a F v , w = D 2 μ F ϵ v : ϵ w d x b F : H 1 D 2 × Q R , b F w , p = D · w p d x .
We assume that f F L 2 0 , T ; L 2 ( D ) 2 , g L 2 0 , T ; H 00 1 / 2 Σ 4 2 and ε > 0 is a penalization parameter. Let χ n + 1 S : D R be the characteristic function
χ n + 1 S ( x ) = 1 , x Ω n + 1 S 0 , x D Ω n + 1 S
where Ω n + 1 S is the image of Ω 0 S by the map φ n + 1 ( X ) = X + u n + 1 ( X ) . We denote Ω n + 1 F = D Ω ¯ n + 1 S and we have D ¯ = Ω ¯ n + 1 S Ω ¯ n + 1 F .
The time-integration scheme for the fluid is: for given u n + 1 , u n , find the velocity v ε n + 1 H 1 ( D ) 2 , v ε n + 1 = g n + 1 on Σ 4 , v ε n + 1 = 0 on Σ 1 Σ 3 Σ 5 and the pressure p ϵ n + 1 Q , such that
ρ F D v ε n + 1 v ε n Δ t · w d x + a F v ε n + 1 , w + b F w , p ε n + 1
+ 1 ε D χ n + 1 S v ε n + 1 ( u n + 1 u n ) φ n + 1 1 Δ t · w d x = D f F , n + 1 · w d x , w W
b F v ε n + 1 , q = 0 , q Q .
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 Ω n + 1 S 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 + b u b b l e (we write just P 1 + b ) for the velocity and P 1 for the pressure. Let T h F be a mesh of D of size h, with n v F vertices and n t F triangles. We introduce
W h = w h C 0 ( D ¯ ) 2 ; T T h F , w h | T ( P 1 + b ) ( T ) 2 , w h = 0 on D Σ 2
Q h = q h C 0 ( D ¯ ) ; T T h F , q h | T P 1 ( T ) .

5. Computing the Forces at the Fluid–Structure Interface Using the Fictitious Domain

We recall the notations: φ n + 1 ( X ) = X + u n + 1 ( X ) , Ω n + 1 S = φ n + 1 ( Ω 0 S ) , Ω n + 1 F = D Ω ¯ n + 1 S . We assume in this section that Ω n + 1 S and Ω n + 1 F 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 Γ n + 1 F S = Ω ¯ n + 1 S Ω ¯ n + 1 F . In the case when no contact occurs, the interface is Γ n + 1 F S = Γ n + 1 , where Γ n + 1 = φ n + 1 ( Γ 0 ) , but when the structure touches the rigid obstacle, the interface is Γ n + 1 F S = φ n + 1 ( Γ 0 ) R S . 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 R S 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 to Ω n + 1 F , respectively v ε n + 1 , S for the restriction of v ε n + 1 to Ω n + 1 S if v ε n + 1 H 1 ( D ) 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 Ω n + 1 F in (30), we get
ρ F v ε n + 1 , F v ε n , F Δ t · σ F v ε n + 1 , F , p ε n + 1 , F = f F , n + 1 , in L 2 ( Ω n + 1 F ) 2 .
If v ε n + 1 , F H 2 ( Ω n + 1 F ) 2 and p ε n + 1 , F H 1 ( Ω n + 1 F ) , by the Green formula we get
Γ n + 1 F S σ F v ε n + 1 , F , p ε n + 1 , F n F · w n + 1 , F d s = ρ F Ω n + 1 F v ε n + 1 , F v ε n , F Δ t · w n + 1 , F d x + Ω n + 1 F 2 μ F ϵ v ε n + 1 , F : ϵ w n + 1 , F d x Ω n + 1 F · w n + 1 , F p ε n + 1 , F d x Ω n + 1 F f F , n + 1 · w n + 1 , F d x
for all w n + 1 , F H 1 ( Ω n + 1 F ) 2 such that w n + 1 , F = 0 on Ω n + 1 F D .
Even if v ε n + 1 , F H 1 ( Ω n + 1 F ) 2 and p ε n + 1 , F L 2 ( Ω n + 1 F ) , we can give a sense of
σ F v ε n + 1 , F , p ε n + 1 , F n F on Γ n + 1 F S by the element j F H 00 1 / 2 Γ n + 1 F S 2 , see [35] (Chapter VII, p. 1241) and [36] (p. 325), defined by
j F , γ Γ n + 1 F S ( w n + 1 , F ) Γ n + 1 F S = ρ F Ω n + 1 F v ε n + 1 , F v ε n , F Δ t · w n + 1 , F d x + Ω n + 1 F 2 μ F ϵ v ε n + 1 , F : ϵ w n + 1 , F d x Ω n + 1 F · w n + 1 , F p ε n + 1 , F d x Ω n + 1 F f F , n + 1 · w n + 1 , F d x
for all w n + 1 , F H 1 ( Ω n + 1 F ) 2 such that w n + 1 , F = 0 on Ω n + 1 F D , where · , · Γ n + 1 F S is the duality between H 00 1 / 2 Γ n + 1 F S 2 and H 00 1 / 2 Γ n + 1 F S 2 and γ Γ n + 1 F S is the trace operator on Γ n + 1 F S . We denote by H 00 1 / 2 ( Γ n + 1 F S ) , the space of the trace on Γ n + 1 F S of functions of { v H 1 ( Ω n + 1 F ) ; v = 0 on Ω n + 1 F Γ n + 1 F S } and we have that H 00 1 / 2 ( Γ n + 1 F S ) H 1 / 2 ( Γ n + 1 F S ) .
Using w with a compact support in Ω n + 1 S in (30), we get
ρ F v ε n + 1 , S v ε n , S Δ t · σ F v ε n + 1 , S , p ε n + 1 , S + 1 ε v ε n + 1 , S ( u n + 1 u n ) φ n + 1 1 Δ t = f F , n + 1 , in L 2 ( Ω n + 1 S ) 2 .
Similary, we define j S H 00 1 / 2 Γ n + 1 F S 2 by
j S , γ Γ n + 1 F S ( w n + 1 , S ) Γ n + 1 F S = ρ F Ω n + 1 S v ε n + 1 , S v ε n , S Δ t · w n + 1 , S d x + Ω n + 1 S 2 μ F ϵ v ε n + 1 , S : ϵ w n + 1 , S d x Ω n + 1 S · w n + 1 , S p ε n + 1 , S d x + 1 ε Ω n + 1 S v ε n + 1 , S ( u n + 1 u n ) φ n + 1 1 Δ t · w n + 1 , S d x Ω n + 1 S f F , n + 1 · w n + 1 , S d x
for all w n + 1 , S H 1 ( Ω n + 1 S ) 2 such that w n + 1 , S = 0 on Ω n + 1 S D .
From (30), we get
j S = j F , in H 00 1 / 2 Γ n + 1 F S 2 .
We have (11) at the fluid–structure interface, or equivalently
σ S n S ( X , t ) = ( σ F n F ) ( φ t ( X ) , t ) .
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
ρ F Ω n + 1 S v ε n + 1 , S v ε n , S Δ t · w n + 1 , S d x + Ω n + 1 S 2 μ F ϵ v ε n + 1 , S : ϵ w n + 1 , S d x Ω n + 1 S · w n + 1 , S p ε n + 1 , S d x + 1 ε Ω n + 1 S v ε n + 1 , S ( u n + 1 u n ) φ n + 1 1 Δ t · w n + 1 , S d x Ω n + 1 S f F , n + 1 · w n + 1 , S d x = ρ F Ω 0 S J n + 1 ( v ε n + 1 , S v ε n , S ) φ n + 1 Δ t · w S d X + Ω 0 S J n + 1 σ F v ε n + 1 , S , p ε n + 1 , S φ n + 1 ( F n + 1 ) T : X w S d X + 1 ε Ω 0 S J n + 1 v ε n + 1 , S φ n + 1 ( u n + 1 u n ) Δ t · w S d X Ω 0 S J n + 1 f F , n + 1 φ n + 1 · w S d X
where F n + 1 ( X ) = I + X u n + 1 ( X ) , J n + 1 ( X ) = det F n + 1 ( X ) , w S = w n + 1 , S φ n + 1 W S .
The right-hand side of the above equality, in the case of linear elasticity, could be approached by
ρ F Ω 0 S v ^ ε n + 1 v ^ ε n Δ t · w S d x + Ω 0 S 2 μ F ϵ X v ^ ε n + 1 : ϵ X w S d X Ω 0 S X · w S p ^ ε n + 1 d X + 1 ε Ω 0 S v ^ ε n + 1 ( u n + 1 u n ) Δ t · w S d X Ω 0 S f ^ F , n + 1 · w S d X
where v ^ ε n + 1 = v ε n + 1 , S φ n + 1 , p ^ ε n + 1 = p ε n + 1 , S φ n + 1 and f ^ F , n + 1 = f F , n + 1 φ 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.

6. 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 h n , u h n 1 K h as well as the fluid velocity v ε , h n W h are given. We look for u h n + 1 K h , v ε , h n + 1 W h and p ε , h n + 1 Q h iteratively, as limits of u h k , n + 1 , v ε , h k , n + 1 and p ε , h k , n + 1 , respectively, for k N .
For the stopping test, we use the parameter t o l > 0 and let ω ( 0 , 1 ) be the relaxation parameter.
  • Step 1. Set the initial displacement of the structure u h 0 , n + 1 = u h n and put k = 0 .
  • Step 2. Find the velocity v ε , h k + 1 , n + 1 g h n + 1 + W h and the pressure p ε , h k + 1 , n + 1 Q h by solving the discrete fluid problem
    ρ F D v ε , h k + 1 , n + 1 v ε , h n Δ t · w h d x + a F v ε , h k + 1 , n + 1 , w h + b F w h , p ε , h k + 1 , n + 1 + 1 ε D χ k , n + 1 S v ε , h k + 1 , n + 1 ( u h k , n + 1 u h n ) φ k , n + 1 1 Δ t · w h d x
    = D f h F , n + 1 · w h d x , w h W h
    b F v ε , h k + 1 , n + 1 , q h = 0 , q h Q h
    where χ k , n + 1 S is the characteristic function of the set Ω k , n + 1 S = φ k , n + 1 ( Ω 0 S ) and φ k , n + 1 ( X ) = X + u h k , n + 1 ( X ) .
  • Step 3. Set v ^ ε , h k + 1 , n + 1 = v ε , h k + 1 , n + 1 φ k , n + 1 , p ^ ε , h k + 1 , n + 1 = p ε , h k + 1 , n + 1 φ k , n + 1 and f ^ h F , n + 1 = f h F , n + 1 φ k , n + 1 , then compute
    ( b ˜ S k + 1 , n + 1 ) i = L ˜ S k + 1 , n + 1 ( ϕ i ) , i = 1 , , 2 n v S
    where
    L ˜ S k + 1 , n + 1 w h S = Ω 0 S f h S , n + 1 · w h S d X + ρ S Δ t 2 Ω 0 S ( 2 u h n u h n 1 ) · w h S d X + ρ F Ω 0 S v ^ ε , h k + 1 , n + 1 v ^ ε , h n Δ t · w h S d x + Ω 0 S 2 μ F ϵ X v ^ ε , h k + 1 , n + 1 : ϵ X w h S d X Ω 0 S X · w h S p ^ ε , h k + 1 , n + 1 d X + 1 ε Ω 0 S v ^ ε , h k + 1 , n + 1 ( u h k , n + 1 u h n ) Δ t · w h S d X Ω 0 S f ^ h F , n + 1 · w h S d X .
The last five integrals represent the surface forces from the fluid acting on the structure as in the end of Section 5.
Step 4. Solve the convex constrained optimization problem
ξ = arg inf ξ R 2 n v S 1 2 A ˜ S ξ , ξ b ˜ S k + 1 , n + 1 , ξ ,
subject to (26) and (27), where A ˜ S is the constant matrix given by (24). Put
u ε , h k + 1 , n + 1 = i = 1 2 n v S ξ i ϕ i .
Step 5. If u ε , h k + 1 , n + 1 u h k , n + 1 0 , Ω 0 S < t o l
  • then u h n + 1 = u ε , h k + 1 , n + 1 , v ε , h n + 1 = v ε , h k + 1 , n + 1 , p ε , h n + 1 = p ε , h k + 1 , n + 1 and STOP,
  • else u h k + 1 , n + 1 = u h k , n + 1 + ω ( u ε , h k + 1 , n + 1 u h k , n + 1 ) , put k : = k + 1 and go to Step 2.
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, add c F ( v ε n + 1 , v ε n , w ) where
    c F ( u , v , w ) = ρ F D [ ( u · ) v ] · w d x ,
  • on the left-hand side of (38), add c F ( v ε , h k + 1 , n + 1 , v ε , h n , w h )
  • on the right-hand side of L ˜ S k + 1 , n + 1 w h S at the Step 3, add ρ F Ω 0 S [ ( v ^ ε , h k + 1 , n + 1 · ) v ^ ε , h n ] · w h S d x .

7. 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++.

7.1. 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: Ω 0 S = ] 0 , 0.3 [ × ] 0 , 0.03 [ , D = ] 0 , 0.41 [ × ] 0.5 , 0.5 [ B ¯ ( 0.2 , 0.07 ; 0.05 ) where B ( 0.2 , 0.07 ; 0.05 ) is the open disk of radius 0.05 m and center ( 0.2 , 0.07 ) . 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 : Ω 0 S × [ 0 , T ] R 2 , f S = ( 0 , 0 ) N / m 3 .
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 : D × [ 0 , T ] R 2 , f F = ( 0 , 0 ) N / m 3 .
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 ) .
We impose at the inflow boundary Σ 4 , the nonhomogenuous Dirichlet condition v = g = ( g 1 , g 2 ) where g 1 = 0 and
g 2 ( x 1 , x 2 , t ) = 4 V m a x x 1 ( x 1 L 1 ) L 1 2 1 c o s ( 2 π t / T 1 ) 2 if t T 1 0 if T 1 < t T
where L 1 = 0.41 m , V m a x = 1 m / s , T 1 = 50 s and T = 65 s . On Σ 1 Σ 3 Σ 5 , a non-slip boundary contition is imposed and σ F n F = 0 on Σ 2 , the outflow boundary. Initially, the structure displasement is zero and the fluid is at rest.
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 t o l = 10 3 × a r e a ( Ω 0 S ) = 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.

7.2. Numerical Simulation of Valve Dynamics

The undeformed structure domain Ω 0 S 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: Γ D = ] 0 , H [ × H 1 ,
Γ C = ( X 1 ( t ) , X 2 ( t ) ) R 2 ; X 1 ( t ) = R + H + R cos ( t ) , X 2 ( t ) = H 1 + R sin ( t ) , t π + π 4 , π + π 2
and Γ N = Ω 0 S ( Γ ¯ D Γ ¯ C ) . Let ψ : R R , ψ ( X 1 ) = 0 be the rigid fondation.
The fluid domain is D = ] L 1 , L 2 [ × ] 0 , H 1 [ . For the numerical tests, we used the values: L 1 = 0.5 cm , L 2 = 2 cm , H 1 = 0.95 cm , see Figure 5.
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 : Ω 0 S × [ 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
g 1 ( X 1 , X 2 , t ) = V m a x ( H 1 2 X 2 2 ) 2 H 1 2 1 c o s ( 2 π t / T 1 ) 2 if t T 1 0 if T 1 < t T
with V m a x = 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: v = g on Σ 1 × [ 0 , T ] , v = 0 on Σ 4 × [ 0 , T ] , σ F n F = 0 on Σ 3 × [ 0 , T ] , σ F n F · e 1 = 0 and v 2 = 0 on Σ 2 × [ 0 , T ] with e 1 = ( 1 , 0 ) T , see 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 t o l = 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 m a x = 10 and 584 s for V m a x = 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 m a x = 10 and after about t = 0.370 for V m a x = 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 m a x = 10 .
In Figure 9, Von Mises stress distribution in the structure is shown for V m a x = 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 m a x = 50 , in Figure 10 and Figure 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 m a x = 50 and a non-linear model for the elasticity equations could be used.

7.3. 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 [ × ] 0 , 1.4 [ (bottom), ] 0.7 , 1 [ × ] 1.4 , 2.6 [ (middle) and ] 0 , 2 [ × ] 2.6 , 4 [ (top). 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 i n ( 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 t o l 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 .

8. 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.

Author Contributions

The authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Quarteroni, A.; Formaggia, L. Mathematical modelling and numerical simulation of the cardiovascular system. In Handbook of Numerical Analysis; Ciarlet, P.G., Ed.; North-Holland: Amsterdam, The Netherlands, 2004; Volume XII, pp. 3–127. [Google Scholar]
  2. Tezduyar, T.E.; Sathe, S. Modelling of fluid-structure interactions with the space-time finite elements: Solution techniques. Int. J. Num. Meth. Fluids 2007, 54, 855–900. [Google Scholar] [CrossRef]
  3. Peskin, C.S. The immersed boundary method. Acta Numer. 2002, 11, 479–517. [Google Scholar] [CrossRef] [Green Version]
  4. Boffi, D.; Gastaldi, L.; Heltai, L.; Peskin, C. On the hyper-elastic formulation of the immersed boundary method. Comput. Methods Appl. Mech. Engrg. 2008, 197, 2210–2231. [Google Scholar] [CrossRef]
  5. Glowinsk, R.; Hesla, T.I.; Joseph, D.; Pan, T.; Périaux, J. A distributed Lagrange multiplier/fictitious domain method for the simulation of flow around moving rigid bodies: Application to particulate flow. Comput. Methods Appl. Mech. Engrg. 2000, 184, 241–267. [Google Scholar] [CrossRef]
  6. Boffi, D.; Gastaldi, L. A fictitious domain approach with Lagrange multiplier for fluid-structure interactions. Numer. Math. 2017, 135, 711–732. [Google Scholar] [CrossRef] [Green Version]
  7. Bost, C.; Cottet, G.-H.; Maitre, E. Convergence analysis of a penalization method for the three-dimensional motion of a rigid body in an incompressible viscous fluid. SIAM J. Numer. Anal. 2010, 48, 1313–1337. [Google Scholar] [CrossRef]
  8. Court, S.; Fournié, M.; Lozinski, A. A fictitious domain approach for fluid-structure interactions based on the extended finite element method. ESAIM Proc. Surv. 2014, 45, 308–317. [Google Scholar] [CrossRef]
  9. Alauzet, F.; Fabrèges, B.; Fernández, M.A.; Landajuela, M. Nitsche-XFEM for the coupling of an incompressible fluid with immersed thin-walled structures. Comput. Methods Appl. Mech. Engrg. 2016, 301, 300–335. [Google Scholar] [CrossRef] [Green Version]
  10. Dos Santos, N.; Gerbeau, J.-F.; Bourgat, J.-F. A partitioned fluid-structure algorithm for elastic thin valves with contact. Comput. Methods Appl. Mech. Engrg. 2008, 197, 1750–1761. [Google Scholar] [CrossRef] [Green Version]
  11. Astorino, M.; Gerbeau, J.-F.; Pantz, O.; Traoré, K.-F. Fluid-structure interaction and multi-body contact: Application to aortic valves. Comput. Methods Appl. Mech. Engrg. 2009, 198, 3603–3612. [Google Scholar] [CrossRef] [Green Version]
  12. Mayer, U.M.; Popp, A.; Gerstenberger, A.; Wall, W. 3D fluid-structure-contact interaction based on a combined XFEM FSI and dual mortar contact approach. Comput. Mech. 2010, 46, 53–67. [Google Scholar] [CrossRef]
  13. Kamensky, D.; Hsu, M.-C.; Schillinger, D.; Evans, J.A.; Aggarwal, A.; Bazilevs, Y.; Sacks, M.S.; Hughes, T.J.R. An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves. Comput. Methods Appl. Mech. Engrg. 2015, 284, 1005–1053. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Hecht, F.; Pironneau, O. An energy stable monolithic Eulerian fluid-structure finite element method. Int. J. Numer. Methods Fluids 2017, 85, 430–446. [Google Scholar] [CrossRef]
  15. Kadapa, C.; Dettmer, W.G.; Peric, D. A stabilised immersed framework on hierarchical b-spline grids for fluid-flexible structure interaction with solid-solid contact. Comput. Methods Appl. Mech. Engrg. 2018, 335, 472–489. [Google Scholar] [CrossRef] [Green Version]
  16. Ager, C.; Schott, B.; Vuong, A.-T.; Popp, A.; Wall, W.A. A consistent approach for fluid-structure-contact interaction based on a porous flow model for rough surface contact. Int. J. Numer. Methods Eng. 2019, 119, 1345–1378. [Google Scholar] [CrossRef] [Green Version]
  17. Burman, E.; Fernández, M.A.; Frei, S. A Nitsche-based formulation for fluid-structure interactions with contact. ESAIM M2AN 2020, 54, 531–564. [Google Scholar] [CrossRef]
  18. Yakhlef, O.; Murea, C.M. Numerical procedure for fluid-structure interaction with the structure displacements limited by a rigid obstacle. Appl. Comput. Mech. 2017, 11, 91–104. [Google Scholar] [CrossRef]
  19. Halanay, A.; Murea, C.M.; Tiba, D. Existence and approximation for a steady fluid-structure interaction problem using fictitious domain approach with penalization. Ann. Acad. Rom. Sci. Ser. Math. Appl. 2013, 5, 120–147. [Google Scholar]
  20. Halanay, A.; Murea, C.M.; Tiba, D. Existence of a steady flow of Stokes fluid past a linear elastic structure using fictitious domain. J. Math. Fluid Mech. 2016, 18, 397–413. [Google Scholar] [CrossRef]
  21. Ciarlet, P.G. Mathematical Elasticity. Volume 1: Three Dimensional Elasticity; Elsevier: Amsterdam, The Netherlands, 2004. [Google Scholar]
  22. Halanay, A.; Murea, C.M.; Tiba, D. Extension theorems related to a fluid-structure interaction problem. Bull. Math. Soc. Sci. Math. Roum. 2018, 61, 417–437. [Google Scholar]
  23. Halanay, A.; Murea, C.M. Uniform Estimation of a Constant Issued from a Fluid-Structure Interaction Problem. In System Modeling and Optimization; IFIP Advances in Information and Communication Technology; Bociu, L., Désidéri, J.-A., Habbal, A., Eds.; Springer: Berlin, Germany, 2016; Volume 494, pp. 292–301. [Google Scholar]
  24. Kikuchi, N.; Oden, J.T. Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods; SIAM Studies in Applied Mathematics; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, USA, 1988; Volume 8. [Google Scholar]
  25. Eck, C.; Jarusek, J.; Krbec, M. Unilateral Contact Problems. Variational Methods and Existence Theorems; Chapman and Hall/CRC: Boca Raton, FL, USA, 2005. [Google Scholar]
  26. Duvaut, G.; Lions, J.-L. Les inÉquations en mÉcanique et en Physique; Travaux et Recherches Mathématiques; Dunod: Paris, France, 1972; Volume 21. [Google Scholar]
  27. Lebeau, G.; Schatzman, M. A wave problem in a half-space with a unilateral constraint at the boundary. J. Differ. Equ. 1984, 53, 309–361. [Google Scholar] [CrossRef] [Green Version]
  28. Chau, O.; Shillor, M.; Sofonea, M. Dynamic frictionless contact with adhesion. Z. Angew. Math. Phys. 2004, 55, 32–47. [Google Scholar] [CrossRef]
  29. Eck, C.; Jarusek, J.; Sofonea, M. A dynamic elastic-visco-plastic unilateral contact problem with normal damped response and Coulomb friction. Eur. J. Appl. Math. 2010, 21, 229–251. [Google Scholar] [CrossRef]
  30. Doyen, D.; Ern, A.; Piperno, S. Time-integration schemes for the finite element dynamic Signorini problem. SIAM J. Sci. Comput. 2011, 33, 223–249. [Google Scholar] [CrossRef] [Green Version]
  31. Krause, R.; Walloth, M. Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme. Appl. Numer. Math. 2012, 62, 1393–1410. [Google Scholar] [CrossRef]
  32. Lewy, H. On the coincidence set in variational inequality. J. Differ. Geom. 1971, 6, 497–501. [Google Scholar] [CrossRef]
  33. Kinderlehrer, D.; Stampacchia, G. An Introduction to Variational Inequalities and Their Applications; Reprint of the 1980 Original. Classics in Applied Mathematics; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, USA, 2000; Volume 31. [Google Scholar]
  34. Rodrigues, J.F. Obstacle Problems in Mathematical Physics; North-Holland Publishing Co.: Amsterdam, The Netherlands, 1987. [Google Scholar]
  35. Dautray, R.; Lions, J.-L. Analyse Mathématique et Calcul Numérique pour les Sciences et les Techniques; Masson: Paris, France, 1988; Volume 4. [Google Scholar]
  36. Boyer, F.; Fabrie, P. Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models; Applied Mathematical Sciences; Springer: New York, NY, USA, 2013; Volume 183. [Google Scholar]
  37. Hecht, F. New development in FreeFem++. J. Numer. Math. 2012, 20, 251–265. [Google Scholar] [CrossRef]
  38. Murea, C.M. Numerical simulation of a pulsatile flow through a flexible channel. ESAIM Math. Model. Numer. Anal. 2006, 40, 1101–1125. [Google Scholar] [CrossRef]
Figure 1. The undeformed structure domain (continuous line) and deformed structure domain (pointed line) in contact with the obstacle.
Figure 1. The undeformed structure domain (continuous line) and deformed structure domain (pointed line) in contact with the obstacle.
Fluids 06 00051 g001
Figure 2. Time history of the vertical displacement of two points at Γ C for X 1 = 0.15 (left) and X 1 = 0.20 (right).
Figure 2. Time history of the vertical displacement of two points at Γ C for X 1 = 0.15 (left) and X 1 = 0.20 (right).
Fluids 06 00051 g002
Figure 3. Computed solution at time instant t = 20 s: pressure (left) and zoom of velocity (right).
Figure 3. Computed solution at time instant t = 20 s: pressure (left) and zoom of velocity (right).
Fluids 06 00051 g003
Figure 4. Undeformed structure domain does not satisfy the non-penetration constraint.
Figure 4. Undeformed structure domain does not satisfy the non-penetration constraint.
Fluids 06 00051 g004
Figure 5. The geometrical configuration of fluid–structure.
Figure 5. The geometrical configuration of fluid–structure.
Fluids 06 00051 g005
Figure 6. Time history of iterations number to obtain the convergence at each time step for V m a x = 10 (left) and V m a x = 50 (right).
Figure 6. Time history of iterations number to obtain the convergence at each time step for V m a x = 10 (left) and V m a x = 50 (right).
Fluids 06 00051 g006
Figure 7. The relatively displacement by rapport of the initial configuration of the bottom corner of the tail segment of the structure: horizontal (left) and vertical (right).
Figure 7. The relatively displacement by rapport of the initial configuration of the bottom corner of the tail segment of the structure: horizontal (left) and vertical (right).
Fluids 06 00051 g007
Figure 8. Computed solution for V m a x = 10 , pressure (left) and velocity (right): when the valve was at the highest position t = 0.125 s (top) and when the valve made contact with the foundation t = 0.230 s (bottom).
Figure 8. Computed solution for V m a x = 10 , pressure (left) and velocity (right): when the valve was at the highest position t = 0.125 s (top) and when the valve made contact with the foundation t = 0.230 s (bottom).
Fluids 06 00051 g008
Figure 9. Von Mises stress distribution in the structure for V m a x = 10 .
Figure 9. Von Mises stress distribution in the structure for V m a x = 10 .
Fluids 06 00051 g009
Figure 10. Computed solution for V m a x = 50 , pressure (left) and velocity (right): when the valve was at the highest position t = 0.120 s (top) and when the valve made contact with the foundation t = 0.370 s (bottom).
Figure 10. Computed solution for V m a x = 50 , pressure (left) and velocity (right): when the valve was at the highest position t = 0.120 s (top) and when the valve made contact with the foundation t = 0.370 s (bottom).
Fluids 06 00051 g010
Figure 11. Von Mises stress distribution in the structure V m a x = 50 .
Figure 11. Von Mises stress distribution in the structure V m a x = 50 .
Fluids 06 00051 g011
Figure 12. Vorticity distribution when the valve starts to be in contact with the foundation for V m a x = 10 (left) and V m a x = 50 (right).
Figure 12. Vorticity distribution when the valve starts to be in contact with the foundation for V m a x = 10 (left) and V m a x = 50 (right).
Fluids 06 00051 g012
Figure 13. Velocity magnitude (left) and pressure (right) of the fluid at t = 25 ms, when the structure is at the highest position.
Figure 13. Velocity magnitude (left) and pressure (right) of the fluid at t = 25 ms, when the structure is at the highest position.
Fluids 06 00051 g013
Figure 14. Vertical displacement of the right bottom corner of the structure.
Figure 14. Vertical displacement of the right bottom corner of the structure.
Fluids 06 00051 g014
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Yakhlef, O.; Murea, C.M. Numerical Simulation of Dynamic Fluid-Structure Interaction with Elastic Structure–Rigid Obstacle Contact. Fluids 2021, 6, 51. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6020051

AMA Style

Yakhlef O, Murea CM. Numerical Simulation of Dynamic Fluid-Structure Interaction with Elastic Structure–Rigid Obstacle Contact. Fluids. 2021; 6(2):51. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6020051

Chicago/Turabian Style

Yakhlef, Othman, and Cornel Marius Murea. 2021. "Numerical Simulation of Dynamic Fluid-Structure Interaction with Elastic Structure–Rigid Obstacle Contact" Fluids 6, no. 2: 51. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6020051

Article Metrics

Back to TopTop