Partitioned FSI strategy for simulations of a thin elastic valve

c

TU Delft, The Netherlands, 2006

PARTITIONED FSI STRATEGY FOR SIMULATIONS OF A

THIN ELASTIC VALVE

N. Diniz dos Santos∗, J.-F. Gerbeau†, J.-F. Bourgat

Institut National de Recherche en Informatique et en Automatique Rocquencourt, B.P. 105, 78153 Le Chesnay, France

∗_{e-mail: nuno-miguel.diniz-dos-santos@inria.fr} web page: http://www-rocq.inria.fr/REO/

† _{e-mail: jean-frederic.gerbeau@inria.fr}

Key words: partitioned strategy, Fictitious Domain, Arbitrary Lagrangian Eulerian, Lagrange Multiplier, Valve movement, FSI.

Abstract. We present a Fictitious Domain (FD)/Lagrange multiplier method [10, 15, 16] to approximate a thin valve movement immersed in an incompressible fluid. We developed a partitioned FSI algorithm that is able to keep the fluid and structure codes independent and which involve only minor modifications to each one [3]. We propose a quantitative comparison between FD and ALE [6], when the displacement is moderated. We show that the two approaches give results in good agreement. We also illustrate FD in a simulation involving very large displacements.

1 INTRODUCTION

further information see [11]). To enforce the continuity of the fluid-structure velocity fields an augmented problem is implemented by means of a Lagrange multiplier. We refer, for example, to [9, 10, 15] for a discussion of Lagrange Multiplier/Fictitious Domain methods in the case of moving particles in a fluid and to [1, 5, 4, 13] for a version of this approach applied to the aortic valve simulation. A similar method applied to biological valves has been developed in [1, 4, 13, 5]. Our approach differs from their work in two points: our structure is indeed a line (in 2D) whereas they consider a 2D structure whose thickness is neglected when the coupling with the fluid is performed; from a computational viewpoint they use a monolithic ad hoc solver whereas we adopt a coupling strategy between two independent solvers.

2 FLUID-STRUCTURE MODELLING

The fluid is considered to be homogeneous, incompressible and viscous. Its domain is defined as ΩF(t) = Ω \Σ(t), as can be seen in Fig. 1

Σ

Γ

Γ

Γ

Γ

Ω

Ω

Figure 1: Bi-dimensional configuration.

The fluid is governed by the incompressible Navier-Stokes equations:    ρ ∂u ∂t + u · ∇u  − η∆u + ∇p = 0 in ΩF(t) × (0, T ), div u = 0 in ΩF(t) × (0, T ), (1) where u is the fluid velocity, p the pressure, ρ the fluid density and η the dynamic viscosity. The Cauchy stress tensor is denoted by

σ = −p I + 2η D(u),

where I is the identity tensor and D(u) is the strain rate (∇u + ∇uT_{)/2. We define}

n as the outward normal on ∂ΩF(t). The total stress is prescribed at the inlet Γin and

the outlet Γout. No-slip boundary conditions are imposed on the wall Γ0 and symmetry

boundary conditions on the central axis Γc (see Fig. 1).

 σ · n(·, t) = −pin(·, t) n on Γin× (0, T ),

where pin is a given function. The system is completed with suitable initial conditions.

The valve Σ(t) defines a natural partition of ΩF(t) into two subdomains Ω−(t) and

Ω+_{(t) (see Fig. 1). On Σ(t), we define n}+ _{(resp. n}−_{) as the outgoing normal on ∂Ω}+

(resp. ∂Ω−_{), and:}

f_Σ = −(σ+· n++ σ−· n−), (2)

where σ+_{(x) (resp. σ}−_{(x)), for x ∈ Σ(t), is the limit of σ(x − εn}+_{) as ε goes to 0}+_(resp.

0−_).

We consider a 1D inextensible structure parametrized by the curvilinear abscissa s. The elastic deformation energy is given by

W (x) = 1 2 Z L 0 EI(s)     ∂2_x ∂s2     2 ds (3)

where EI denotes the flexural stiffness, the Young Modulus E multiplied by the inertia I, x(s) is the position vector of a point on the structure’s axis. We suppose that the external forces exerted on the elastic structure depend solely on the curvilinear s-coordinate and the position of its axis.

The solutions will be searched within the following set K(t) = {x ∈ [H2_{(0, L)]}2_;∂2x ∂s2 ∈ [L 2_{(0, L)]}2_;     ∂x ∂s     2 = 1; x(0, t) = 0} (4) which is the set of the allowed structure configurations at time t that satisfy the boundary and the inextensibility conditions.

Thus, the configurations of the immersed structure are the solutions of the following variational formulation:                 

Find x(s, t) ∈ K(t) such that for all t Z L 0 m∂ 2_x ∂t2 ξds + Z L 0 EI∂ 2_x ∂s2 ∂2_ξ ∂s2ds = Z L 0 f_Σξds, ∀ξ ∈ K(t) x(s, 0) = x0 ∂x ∂t(s, 0) = ˙x 0 x(0, t) = 0 (5)

where m denotes the linear mass of the structure and f_Σ in defined in (2).

The fluid-structure coupling takes place on the fluid-structure interface. In the con-figuration at hand, the structure being immersed and having one dimension less than the fluid, the interface coincides with the structure domain Σ(t). We introduce the trace operator on Σ

TrΣ : [H1(Ω)]2 −→ [H1/2(Σ)]2

v 7−→ TrΣv = v|Σ.

∂x

∂t (s, t) = u(x(s, t), t). (6)

In the sequel we will denote by uΣ(xΣ, t) the structure velocity on a point xΣ ∈ Σ(t).

3 DISCRETIZATION

3.1 Fluid problem

3.1.1 Time discretization of the augmented formulation

Problem (1) is discretized in time by a semi-implicit Euler scheme. The velocity con-dition (6) is imposed via a Lagrange multiplier. We introduce the following functional spaces:

X = {v ∈ (H1(Ω))2, v = 0 on Γ0, v · n = 0 on Γc},

V = {v ∈ X, TrΣ(v) = 0, div v = 0},

V (uΣ) = {v ∈ X, TrΣ(v) = uΣ, div v = 0}.

We suppose that Σn+1 _{and u}n+1

Σn+1 are given. A first variational formulation of the fluid problem is: Find un+1_{∈ V (u}n+1

Σn+1) such that, for all v ∈ V , Z Ω ρu n+1_{− u}n δt · v + Z Ω ρ un· ∇ un+1· v + Z Ω 2η D(un+1) : D(v) = − Z Γin pinn· v. (7)

In view of the discretization, it is much more convenient to work with space X. Thus, following [9], we introduce the Lagrange multiplier spaces

M = L2(Ω),

Ln+1 = (H−1/2(Σn+1))2,

corresponding to the two constraints which define the space V (un+1

Σn+1). We then consider the following variational formulation: find (un+1_{, p}n+1_{, λ}n+1_{) ∈ X × M × L}n+1 _{such that,}

where h·, ·i denotes the duality pairing on (H−1/2_(Σn+1₎₎2 _{× (H}1/2_(Σn+1₎₎2_{. If Green’s}

formula is applied in the integrals over Ω− _{and Ω}+ _{(see Fig. 1) and, using the fact that}

(un+1_{, p}n+1_{) solves problem (8), we obtain:}

λn+1 = fn+1_Σ . (11)

In other words, the Lagrange multiplier λn+1 corresponding to the constraint un+1 ₌

un+1_Σn+1 represents the jump of the hydrodynamic stress through the valve. 3.1.2 Space discretization

The spaces X and M are approximated by finite element spaces Xh and Mh. The

sim-ulations presented in this paper have been performed withP1/P1 stabilized finite elements.

The 1D structure seen from the fluid is defined by NΣdiscretization points (xn+1i )i=1,...,NΣ on Σn+1_{. The Lagrange multipliers space L}n+1 _{is approximated by}

Ln+1_h = {µ_h measure on Σn+1, µ_h =

NΣ X

i=1

µ_iδ(xn+1_i ), µ_i ∈R2},

where δ(xn+1_i ) denotes the Dirac measure at point xn+1_i . Note that since Xh ⊂ (C0(Ω))2

the quantity < µh, T rΣn+1(v_h) >= NΣ X i=1 µ_ivh(xi), for µh ∈ Lh and vh ∈ Xh (12) is well defined. 3.2 Structure problem

3.2.1 Time discretization and solution strategy

We use a Houbolt scheme for time discretization. The solution is found solving the following problem:

 



Find xn _{∈ K}n _{such that}

Z L 0 m2x n_{− 5x}n−1_{+ 4x}n−2_{− x}n−3 δt2 ξds + Z L 0 EI∂ 2_xn ∂s2 ∂2_ξ ∂s2ds = Z L 0 fn_Σξds (13) where Kn _{is the allowed configurations set at time t = nδt. The numerical solution of}

(13) is equivalent to finding the minimum of the following functional J(x) = 2 δt2 Z L 0 m |x|2ds + 1 2 Z L 0 EI|x′′|2ds − Z L 0 fn_Σxds − 1 δt2 Z L 0 m(−5xn−1+ 4xn−2− xn−3)xds (14)

3.3 Solution strategy

To enforce the fluid-structure coupling, a fixed-point method is used at each time iter-ation. Problems (9),(11) and (13) are solved successively. The iterations are accelerated by the Aitken formula (see [6, 8] for more details).

4 NUMERICAL TESTS

4.1 ALE comparison

We wished to compare the displacement and the jump of stres through the immersed valve given by ALE and FD formulations. We consider a stiff valve (Fig. 2) that ex-periences moderate displacement and thus would not force the ALE method to resort to remeshing due to elements distortion. The data are: ρ = 1.0 g cm−3_{, η = 0.1 poise,}

EI = 0.7 ×10−2_{g cm}3_s−2 _{(flexural stiffness), m = 2.5 g and a sinusoidal pressure function}

with amplitude A = 40 g cm−1_s−2 _{is imposed at the inlet.}

hs Method error L/9 ALE 0.00441579 FD 0.00502732 L/18 ALE 0.00345255 FD 0.0045414 L/27 ALE 0.00216656 FD 0.00333814 Table 1: The relative error of the displacement on the L∞

(0, T ; L∞

(Σ)) norm, taking the ALE with hs=

L

45 as the reference. The results are shown for different structural space steps hs

The length of the valve is L = 0.45 cm, the tube is 6cm long and 1cm high. We took as benchmark the ALE result with hs = L/45 = 0.01. Then we compared the position of

the apex along the x-axis (see Fig. 3 and Fig. 32).

All the results, even those obtained with 9 space discretization steps on the valve, are quite close to the actual results(see Tab. 1). For the loads on the structure, using the L∞

norm we found a 6% error between both methods. 4.2 Two valves with large displacement

We intended on verifying the robustness of the algorithm in presence of very large displacements. To this extend we introduced two valves clamped on the outer wall. In the inlet Γin we defined a periodic pressure function which represents approximately the

Figure 2: Comparison of the iso-values of the velocity obtained with ALE and FD methods.

5 CONCLUSIONS

We have presented a Fictitious Domain and a partitioned strategy to solve the cou-pling between an immersed thin structure and an incompressible fluid. Numerical tests have demonstrated a good agreement with ALE simulations in the case of moderate dis-placements and a good robustness in the case of very large disdis-placements. Moreover, the present approach is well-suited to manage contact. This will be presented in a forthcoming work [6].

6 Acknowledgments

2.85 2.9 2.95 3 3.05 3.1 3.15 3.2 0 1 2 3 4 5 6 ALE 45 FD 9 FD 18 FD 27 2.85 2.9 2.95 3 3.05 3.1 3.15 4.7 4.8 4.9 5 5.1 5.2 5.3 ALE 45 FD 9 FD 18 FD 27

Figure 4: Two valves immersed in a 6 by 1cm tube where a pressure difference is imposed. The pressure and velocity fields are shown respectively by colour and arrows.

REFERENCES

[1] F.P.T. Baaijens. A dlm/fd method for fluid/flexible-body interactions. A fictitious domain/mortar element method for fluid-structure interaction, 35:743–761, 2001. [2] J.-F. Bourgat, J.M. Dumay, and R. Glowinski. Large displacement calculations of

flexible pipelines by finite element and nonlinear programming methods. Journal on Scientific Computing, 1:34–81, 1980.

[3] P. Causin, N. Diniz dos Santos, J.-F. Gerbeau, C. Guivier, and P. Metier. An embedded surface method for valve simulation. application to stenotic aortic valve estimation. ESAIM Proceedings, 14, 2005.

[5] J. De Hart, G.W.M. Peters, P.J.G. Schreurs, and F.P.T. Baaijens. A three-dimensional computational analysis of fluid-structure interaction in the aortic valve. Journal of Biomechanics, 36:103–112, 2003.

[6] N. Diniz dos Santos, J.-F.. Gerbeau, and J.-F. Bourgat. A partitioned fictitious domain approach for the 2d simulation of thin elastic valves. in preparation, 2006. [7] D. Garcia, P. Pibarot, J. Dumesnil, F. Sakr, and L.-G. Durand. Assessment of aortic

valve stenosis severity, a new index based on the energy loss concept. Circulation, 101:765–771, 2000.

[8] J.-F. Gerbeau and M. Vidrascu. A quasi-newton algorithm based on a reduced model for fluid-structure interactions problems in blood flows. Math. Model. Num. Anal., 37(4):631–648, 2003.

[9] R. Glowinski, T.-W. Pan, T.I. Hesla, and D.D. Joseph. A distributed lagrange multiplier/fictitious domain method for particulate flows. International Journal of Multiphase Flow, 25:755–794, 1999.

[10] R. Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph, and J. P´eriaux. A fictitious domain approach for the direct simulation of incompressible fluid flow past moving rigid bodies: Application to particulate flow. Journal of Computational Physics, 169:363–426, 2001.

[11] R. Glowinski, T.-W. Pan, T.I. Hesla, and J. P´eriaux. A fictitious domain method for dirichlet problem and applications. Computational Methods in Applied Mechanics and Engeneering, 111:283–303, 1994.

[12] Y.G. Lai, B.C. Krishnan, and J. Lemmon. Numerical simulation of mechanical heart valve closure fluid dynamics. Journal of Biomechanics, 35:881–892, 2002.

[13] R. Loon, P.D. Anderson, J. De Hart, and F.P.T. Baaijens. A combined finctitious domain/adaptative meshing method for fluid-structure interaction in heart valves. International Journal for Numerical Methods in Fluids, 46:533–544, 2004.

[14] B. Maury. Direct simulations of 2d fluid-particle flows in biperiodic domains. Journal of Computational Physics, 156:325–351, 1999.

[15] N.A. Patankar, P Singh, D.D. Joseph, R. Glowinski, and T.-W. Pan. A new formu-lation of the distributed lagrange multiplier/fictitious domain method for particulate flows. International Journal of Multiphase Flow, 26:1509–1524, 2000.