A Numerical study of fluid flow in peristaltic pumps and polymeric elastic pipes

(1)

LAPPEENRANTA UNIVERSITY OF TECHNOLOGY School of Engineering Science

Double Degree Program in Computational Engineering and Mathematics

Gaetano Formato

A NUMERICAL STUDY OF FLUID FLOW IN PERISTALTIC PUMPS AND POLYMERIC ELASTIC PIPES

Examiners: Prof. J. Sorvari Dr. P. Causin Prof. T. Koiranen Prof. T. Kauranne

(2)

ABSTRACT

Lappeenranta University of Technology School of Engineering Science

Double Degree Program in Mathematics and Computational Engineering Gaetano Formato

A NUMERICAL STUDY OF FLUID FLOW IN PERISTALTIC PUMPS AND POLYMERIC ELASTIC PIPES

Master's thesis 2016

105 pages, 85 gures, 4 tables

Examiners: Prof. J. Sorvari, Lappeenranta University of Technology Dr. P. Causin, University of Milan

Prof. T. Koiranen, Lappeenranta University of Technology Prof. T. Kauranne, Lappeenranta University of Technology

Keywords: Computational Fluid Dynamics (CFD), Turbulence Models, Pulsatile ow, Womersley ow, Nonlinear elastic model, Fluid Structure Interaction (FSI), Defective Boundary conditions, Wall Functions, Arbitrary Lagrangian Eulerian approach (ALE), Finite Element Analysis (FEA), ABAQUS, COMSOL, MATLAB, C++

(3)

3 This thesis work deals with a mathematical description of ow in polymeric pipe and in a specic peristaltic pump. This study involves uid-structure interaction analysis in presence of complex-turbulent ows treated in an arbitrary Lagrangian-Eulerian (ALE) framework. The ow simulations are performed in COMSOL 4.4, as 2D axial symmetric model, and ABAQUS 6.14.1, as 3D model with symmetric boundary conditions. In COMSOL, the uid and structure problems are coupled by monolithic algorithm, while ABAQUS code links ABAQUS CFD and ABAQUS Standard solvers with single block- iterative partitioned algorithm. For the turbulent features of the ow, the uid model in both codes is described by RNG k −ϵ. The structural model is described, on the basis of the pipe material, by Elastic models or Hyperelastic Neo-Hookean models with Rayleigh damping properties.

In order to describe the pulsatile uid ow after the pumping process, the available data are often defective for the uid problem. Engineering measurements are normally able to provide average pressure or velocity at a cross-section. This problem has been analyzed by McDonald's and Womersley's work for average pressure at xed cross section by Fourier analysis since '50, while nowadays sophisticated techniques including Finite Elements and Finite Volumes exist to study the ow.

Finally, we set up peristaltic pipe simulations in ABAQUS code, by using the same model previously tested for the uid and the structure.

(4)

Acknowledgments

Hello readers (hopefully scientic readers),

This is the last step in order to ultimate my Master Degree for new paths.

Everything can not be possible without helps and suggestions of people I have met in my personal and academic growing along this Double Degree Program.

I express my gratitude to prof. Joonas Sorvari, to drive me into Industrial Mathematics applications during my thesis work and suggestions on the use of commercial codes, and prof. Paola Causin, to give me the opportunity for Double Degree and to help me to nalize the thesis in the summer period. Furthermore I want to thank you prof.

Tuomas Koiranen to get me more near on engineering approach than pure theoretical approach for the suggested problems.

Another important academic gure is prof. Tuomo Kauranne. I thank him for a help in nalizing thesis, for the opportunity of double degree program and for the Italian

"coee break" tradition trade from Italy to Lappeenranta.

I thank you ECMI for the possibility to attend an international and industrial Math- ematics program. I thank you Lappeenranta people Sebastian, Ramona and Roberta for playing card nights during Aurora Borealis, Dipal for addiction to coee, CEID guys for the discussion in the oce of engineering topics and not, Masudul and the rst semester Italian and international people.

Among Milan people I want to thank you "I Forestieri", for the help in the everyday situation during the rst year out of Naples. In particular I want to thank you Adriana, for introducing in Milan lifestyle, and Maristella, for teaching of Latex. Furthermore, I thank you Corsico people because they make me feel as home.

Last but not least, I want to thank you Naples people. In particular, I want to thank you my family that has trusted in my work since the beginning.

I have to thank you in this academic work people that see my academic beginning: prof.

Gaetano Festa for the discussion and indications for my academic path, Il Loggione for the continuous academic suggestions. I am thankful of all people I grow up with:

Peppe, TR family, Roberta, Lavinia, Calbi R., Cirillo family, Aldo, Gioia, Francesco.

I think the author of this thesis work I am not but all these people. I am only a mean.

Lappeenranta, Month 09, 2016.

Gaetano Formato

(5)

Chapter 1 Introduction

The need of a more accurate control of uid ow for applications in Medicine and Engineering has increased the interest in performing numerical simulations. The control of the uid ow derives from the uid-induced deformation or the boundary moving- induced deformation, that can be caused by the movement of submerged mechanisms or the movement of external boundaries. The second approach is nowadays widely used in chemical and medical processes, because the uid does not directly interact with the pumping mechanisms.

A positive displacement pump makes a uid move by trapping a xed amount of uid and forcing (displacing) that trapped volume into the discharge pipe.

Peristaltic pumps are a positive displacement pumps based on the peristaltic eect.

The uid is pushed forward through the hose pipe with compression and relaxation phases caused by:

1. a drive shaft with circular revolution motion, as in Figure 1.1.

2. a drive shaft with movable rotation axis with epicycloid motion, as in Figure 1.2.

Figure 1.1: A peristaltic pump with revolution motion, by Williamson Manufacturing.

(13)

Figure 1.2: A specic peristaltic pump with epicycloid motion: LPPT32 Flowrox Oy.

By changing kinematic parameters of the rotor motion, the pressure on the hose pipe varies. This variation induces a deformation of the hose pipe. In this framework, the uid ow is controlled by the rotor kinematic.

1.1 Objectives

The study of uid structure interaction involves 2 sub-problems:

• the uid ow subproblem, describing the uid behaviour;

• the structure subproblem, describing the uid boundary deformation.

This thesis work focuses on:

• the analysis of the ow mass rate along an elastic pipes with a pulsatile inlet and comparison between dierent commercial Software and ways to couple the uid and solid physics;

• the simulation of the ow eld for a specic peristaltic pump with comparison between experimental data.

2

(14)

1.2 Research methodology

Elastic pipe ow problems are set up on COMSOL 4.4 and ABAQUS 6.14.1. These commercial Software adopt Finite Element Analysis. The uid part is described by single-phase Reynolds Averaging Navier Stokes model. The structural part is described by the equilibrium law for deformable bodies, including Rayleigh damping features.

Because of complexity of a specic peristaltic pump models that present both structural interactions and uid-structure interaction, the analysis is implemented in ABAQUS 6.14.1 and it is run on the CSC Taito supercluster.

1.3 Organization of the study

This thesis work is divided into 5 chapters. Chapter 2 deals with the mathematical preliminaries for the discussion of the problems. Chapter 3 gives an introduction of a generic uid structure interaction (FSI) problem, by analysing dierent coupling schemes. Chapter 4 describes in mathematical ways specic test cases that involve FSI problems. In Chapter 5, there is a discussion of results and conclusions about the simulated models.

(15)

Chapter 2 Mathematical Preliminaries

In this chapter, we introduce some useful mathematical denitions and theorems.

Denition 2.1 (The Euclidean norm in Rⁿ) The Euclidean norm in Rⁿ is dened as:

∥·∥Rn :X ∈Rⁿ−→

v u u t

n

∑

i=1

X_i² ∈R

Denition 2.2 (The Lebesgue Space) Let Ω⊆Rⁿ be a domain of dimension m ≤ n. We dene the Lebesgue space on Ω of positive integer exponent 1 ≤ p <∞

L^p(Ω) =







v : Ω−→R :





∫

Ω

|v|^p dΠ_m





1/p

<+∞







(2.1) where Πm is the m-dimensional Hausdor measure on Rⁿ.

If p= +∞,

L^∞(Ω) ={v : Ω−→R : |v|<+∞ a. e. in Ω} (2.2) In particular v ∈L^p_loc(Ω) i ∀K ⊂⊂Ω, ∫

K

|v|^p dΠm <+∞

It is possible to equipL^p(Ω) with the norm:

4

(16)

• 1≤p < ∞:

|v|_Lp(Ω):=





∫

Ω

∥v∥^p dΠm



 1

p (2.3)

• p= +∞,

|v|_L^∞_(Ω) := inf{k ∈R:|v(x)| ≤k a.e. inΩ} (2.4) L^p(Ω) is a Banach space respect to the norm ∥·∥_Lp(Ω) for all p. For p= 2, the space is a Hilbert space with scalar product dened as:

v, w∈L²(Ω), (v, w)_L²_(Ω) :=

∫

Ω

v(x)w(x) dΠm (2.5)

Denition 2.3 (The Weak Derivative) Let Ω⊆Rⁿ be a m-dimensional domain, m≤n. Given v ∈L¹_loc(Ω) and α ∈Nⁿ, the α-th weak derivative of v is a function g ∈L¹_loc(Ω) such that ∀ϕ ∈C₀^∞(Ω):

∫

Ω

v ∂^|α|ϕ

∂x^α₁¹...x^α_nⁿ dΠm =−

∫

Ω

gϕ dΠm (2.6)

where |α|=

n

∑

i=1

αi

The space C₀^∞(Ω) is the subset of C^∞(Ω) formed by functions with compact support.

This space is called the test function space.

Denition 2.4 (The Sobolev Space) Let Ω⊆Rⁿ be a m dimensional domain, m ≤ n. We dene the Sobolev space on Ω of positive integer exponents 1 ≤ p≤ ∞ and 1≤k <∞

W^k,p(Ω) = {

v ∈L^p(Ω) : ∂^|α|v

∂x^α₁¹...x^α_nⁿ ∈L^p(Ω), ∀ |α| ≤k }

(2.7) where the symbol of derivative is in weak sense.

The space W^k,p(Ω) is equipped by the norm

• 1≤p < ∞:

∥v∥_W^k,p_(Ω) :=





∑ ∥D^αv∥_L^p_(Ω)





1/p

(2.8)

(17)

• p=∞:

∥v∥_Wk,∞(Ω):= max

0≤|α|≤k∥D^αv∥_Lp(Ω) (2.9)

We dene the space H^k(Ω) = W^k,2, ∀k ∈ N. It is a Hilbert space respect the dened norm for allk.

Theorem 2.1 (The Trace theorem) Let Ω be a m-dimensional domain, m ≤ n, with Lipschitz boundary, Γ, and p ≥ 1. Then there exists a linear and continuous function T:W^1,p(Ω)−→L²(Γ) such that:

if v ∈W^1,p(Ω)∩C(Ω), then Tv =v|Γ;

∃Cf >0 such that ∥v∥_Lp(Γ)≤Cf∥v∥_W1,p(Ω), ∀v ∈W^1,p(Ω) We dene the space W₀^1,p(Ω) :=C₀^∞(Ω)^W

1,p. The space can be characterized:

v ∈W₀^1,p(Ω)⇐⇒v ∈W^1,p(Ω) and Tv = 0 a.e. on Γ (2.10) We dene H⁻¹(Ω) the dual space of H₀¹(Ω) := W₀^1,2(Ω). In [10], it is observed that L²(Ω)⊂H⁻¹(Ω).

The other inclusion is not veried. Some functionals inH⁻¹(Ω)are not functions inL²(Ω). It was proven by J.Lions [7]:

Theorem 2.2 (The Lions Lemma) Let Ω be a domain in Rⁿ and let v ∈ H⁻¹(Ω). If ∂v

∂xi

∈H⁻¹(Ω), then v ∈L²(Ω)

The Lions Lemma is used for the proof of the following useful inequality.

Theorem 2.3 (The Korn inequality) Let Ωbe a domain in R³ and Γ0 ⊂∂Ω with positive area. Then there exists a constant C such that

∥e(v)∥_L²_(Ω)≥C∥v∥_H¹_(Ω) (2.11) for all v∈H¹(Ω)Γ0 ={v ∈H¹(Ω) : v =0 on Γ0 in sense of trace}.

The quantity e(v) := 1

2(∇v+∇v^T).

Theorem 2.4 (The Lax-Milgram Lemma) Let H be a Hilbert space and ˜a be a continuous and coercive bilinear form. Then for each l ∈H^∗, dual space of H, it exists an unique element u∈H such that:

˜

a(u, v) = l(v), ∀v ∈H (2.12) 6

(18)

The Korn inequality and Lax Milgram theorem allow to get the existence and the uniqueness results for solutions of structural problems [7] .

We need to amply the denition of derivative, also for vector spaces of any di- mensions.

Denition 2.5 (The Gateaux derivative) LetU and V be vector spaces, f : V −→U be a function andh̸= 0 andxbe elements ofV. The Gateaux derivative in x along h, ∂f

∂h(x), is dened:

∂f

∂h(x) := lim

ϵ→0

f(x+ϵh)−f(x)

ϵ (2.13)

(19)

Chapter 3 Mathematical modeling of FSI problems

In this chapter, we separately analyze the two physics, uid and solid, by underlying the complexity of the corresponding problems.

We introduce a suitable framework for the problem dierent from classical La- grangian or Eulerian frameworks. Furthermore we discuss how the physics can be coupled and advantages-disadvantages of dierent coupling methods.

3.1 Fluid subproblem

Let us consider a xed uid domain Ω^f ⊂ R³. We use an Eulerian approach in order to model the uid velocity eld, u^f, and the pressure, p^f: the model does not focus on single uid particles but on specic locations in the space where the uid ows. We assume that:

the uid is incompressible: the uid density, ρ^f , does not change in the time and the space;

the thermal eects can be neglected.

Under these assumptions, we derive the following balance equation [16, 25, 28]:

ρ^f∂u^f

∂t +ρ^f(u^f · ∇)u^f =−∇p^f +divσ^f +f^f in Ω^f (3.1a)

∇ ·u^f = 0 inΩ^f (3.1b)

8

(20)

whereσ^f is the uid stress tensor and f^f is the external volumetric force.

The rst equation is called momentum conservation equation and the second one is called mass conservation equation.

The stress tensor should depend on the pressure of the uid and on the sdevia- toric part of the stress τ.

σ^f =−p^fI+τ

If the working uid is Newtonian, the shear tensor may be approximated by the following constitutive law:

τ =µ^f(

∇u^f +(

∇u^f)T)

− 2

3∇ ·u^fI whereµ^f is the dynamic viscosity of the uid.

Thanks to this assumption, we can express the stress tensor as:

σ^f =−p^fI+µ^f(

∇u^f +(

∇u^f)T)

− 2

3∇ ·u^fI We obtain the following form of the Navier Stokes equations, NS:

ρ^f∂u^f

∂t +ρ^f(u^f · ∇)u^f =−∇p^f +µ∆u^f +f^f in Ω^f (3.2a)

∇ ·u^f = 0 inΩ^f (3.2b)

In the following analysis, we suppose that external volumetric forces have no eects on the uid ow, f^f = 0.

3.1.1 Generalities on Turbulence models

For industrial and biological applications, the Navier-Stokes model (NS) is not optimal because it has a heavy computational cost. In these applications the ow is characterized by the presence of small turbulent structures, eddies, that inuence the main ow. In order to use the NS model to get accurate numerical results, we should enormously decrease the mesh size in Finite Element or Finite Volume analysis, [18].

Turbulence models are engineering approximations in order to reduce the complexity in uid ow simulations. The turbulence behaviour in the uid motion

(21)

(22)

If the velocity eld cannot be derived from a potential, we get from equations (3.2) :

ρ^fdω^f

dt +ρ^fu^f ·(∇ω^f) = µ^f∇²ω^f + (ω^f · ∇)u^f (3.4) The term (ω^f · ∇)u^f is called the vortex stretching. One can consider a pipe ow with symmetry for ⃗x axis. If the pipe stretches along ⃗x axis , the vorticity increases for conservation of angular momentum. The eect of vortex stretching is to produce other scales more and more little where the vorticity increases.

In Lagrangian framework, Equation (3.4) gets:

ρ^fDω^f

Dt =µ^f∇²ω^f + (ω^f · ∇)u^f (3.5) If the problem is dened in the ⃗x−⃗y plane, the vorticity is parallel to z direction. In this situation the vortex stretching term reads:

ω^f_z (∂u^f_x

∂z ⃗x+∂u^f_y

∂z ⃗y )

=0

because u^f does not depend on the z coordinate. If the uid is inviscid or has low viscosity, from Equation 3.5 uid particles have constant vorticity,

Dω^f

Dt = 0. The conservation of vorticity for the uid particles is not ad- missible for turbulent phenomena because it implies a regularity of vorticity eld.

This implies that turbulence models cannot be reduced by a 2D models:

turbulence phenomena are intrinsically 3D.

Dierent scale of analysis

In a turbulent ow there are large and small motion scales. These two motion scales are independent within some limits. For a bounded uid ow, the scales can be detected near the boundaries, where big eddies transfer energy to small eddies (energy cascade). This energy is consequentially dissipated into heat by small eddies thanks to molecular viscosity.

3.1.2 Mathematical analysis for k − ϵ model

In order to introduce a mathematical model for turbulent ows, we assume that the ow has two scales so that there is an underlying principle to separate the mean from the oscillations.

(23)

For an incompressible Newtonian uid, the uid ow may be divided into turbulent component and principal ow component. We lter velocity eld, u^f, and pressure, p^f, by time average, called Reynolds average.

Denition 3.1 Given a functionf : Ω^f×I −→Rⁿ, associated to the uid ow, we dene Reynolds Average of f

f(x) = (

T−→+∞lim 1 T

∫ +∞

0

fi(x, t) dt )

i=1,...n

(3.6)

We dene U^f := u^f and P^f := p^f Reynolds averaging velocity and Reynolds averaging pressure. The turbulent component for the uid pressure and velocity are respectively dened as u^′ :=u^f −U^f and p^′ :=p^f −P^f.

In order to get a model for these quantities, we lter the NS equation by time average. The ltered model is called in literature Reynolds Averaging Navier Stokes model (RANS) [29].

The model is not closed because we provide a model for the turbulent component of velocity, u^′. We dene turbulence kinetic energy, k, and turbulent dissipation rate, ϵ, as:

k = 1 2

3

∑

i=1

u^′2i and ϵ = µ^f ρ^f

3

∑

i,k=1

(∂u^′i

∂xk

)2

In the test cases we consider a RANS-based model with k and ϵ. The uid ows are modeled by the renormalization group (RNG)k−ϵ model. It reads:

ρ^f∂U^f

∂t +ρ^f(U^f · ∇)U^f =−∇P^f +(

µ^f +µt

)∆U^f − 2

3∇k (3.7a)

∇ ·U^f = 0 (3.7b)

ρ^f∂k

∂t +ρ^fU^f · ∇k = (µ^f + µt

σk

)∆k−ρ^fϵ+P_k (3.7c)

ρ^f∂ϵ

∂t+ρ^fU^f · ∇ϵ= (µ^f +µt

σϵ

)∆ϵ−C2ϵρ^fϵ²

k +ρ^fϵ

kC1ϵPk−Rϵ (3.7d) where:

µt is the turbulent dynamic viscosity and it depends on k and ϵ;

Pk:= µt

2

[∇U^f +∇(

U^f)T]2

;

Rϵ is a dissipation factor expressed as function ofk,ϵand velocity gradient;

σk, σϵ,C1ϵ, C2ϵ are closure coecients, that are empirically computed:

σ_k =σ_ϵ = 0.7178, 1.3,C_1ϵ = 1.42, C_2ϵ= 1.68 12

(24)

The RNGk−ϵmodel produces good results in complex geometry for the presence of dissipation term Rϵ. The RNG k−ϵ approach is the most suitable for our simulations.

In bounded uid ows, the dissipative viscous phenomena near the wall determine viscous forces are dominant. A dierent treatment near the wall is thus necessary in order for the strong velocity gradient. An articial internal wall is dened,Γh, located at distance y˜ from wall boundary, Γwall. Under the assumption of a tangential velocity behaviour, it is possible to write onΓh:

u^f

uτ

= 1

κlogy⁺+βc, k = u²_τ

√Cµ

,ϵ= u³_τ

κy (3.8)

where

1. uτ is the turbulent velocity that depends on shear wall stress and density;

2. κ= 0.41is the von Karman constant;

3. βc = 5.1;

4. y⁺:= ρy

µuτ, wherey is the height from boundary wall.

In our simulations, we use commercial software ABAQUS and COMSOL, [23, 29].

The wall-normal distance is frequently calculated and updated, because of the moving boundary wall, [1, 8]. This adds more complexity in the numerical model.

3.2 Structural problem: small and nite deforma- tions

For the structural part of FSI problems, we are interested in computing the displacement vector eld.

LetΩ^s⊂R³ be the geometry of a body in reference conguration. Applied forces cause the reference domain to deform . The deformation can be described by the map:

ψ:I×Ω^s −→R³

The deformation function is injective on Ω^s, i.e. we avoid self-contact in the deformed domain along the motion, and without loss of generality orientation

(25)

preserving, i.e. det(∇ψ)>0.

From Linear Algebra, there is a unique decomposition of F:=∇ψ as:

F=RU=VR

where R is an orthogonal matrix, that represents rigid rotation, U and V are the symmetric right and left stretch tensors, respectively. In the following analysis we need to dene the tensors C := F^TF = U^TU = U² called Right Cauchy- Green Deformation Tensor and D := FF^T = VV^T = V² called Left Cauchy- Green Deformation Tensor . From the symmetry and strictly denite positiveness properties, C has 3 positive eigenvalues (λ²₁, ²2, λ²₃) that determines 3 invariants:

I1 =tr(C) = λ²₁+λ²₂+λ²₃ I2 = 1

2

[tr(C)²−tr(C²)]

=λ²₁λ²₂ +λ²₁λ²₃+λ²₂λ²₃ I3 =detC=λ²₁λ²₂λ²₃

The displacement eld is dened asη^s :=ψ−Id, where Id is the identity function onΩ^s.

From Euler structural theory, [7], we get:

Denition 3.2 Given volumetric, or densities, applied force f^s continuous for spacial variable and external applied pressure g_ext on S ⊂ ∂Ω^s, ∃T^s : Ω^s −→

M₃(R) symmetric matrix eld that satises:

ρ^sD²η^s

Dt² − ∇ ·T˜^s=f^s on Ω^s T˜^sn =g_ext on S

whereρ^sis the body density, n is the local normal surface eld andT˜^s= (det∇ψ)⁻¹T^s∇ψ^T .

The tensor T^s is called the rst Piola-Kircho stress tensor eld.

We dene the second Piola-Kircho stress tensor eld Σ : Ω^s −→ M₃(R) as:

Σ:=∇ψ⁻¹T

14

(26)

Another relevant quantity for the analysis is Σ˜ := (det∇xψ)⁻¹Σ∇xψ^T.

It is possible to determine a general form ofΣfor an isotropic and homogeneous elastic material body when the reference conguration is the natural state:

Σ=λ(trE)I+2µE+o(∇η^s) and E:= 1

2(∇ η^s∇η^sT−I) = 1

2(∇η^s+∇η^sT)+ox(∥∇η^s∥) The coecients λ and µ are called Lamè parameters and they are functions of Poisson ration ν and Young modulusE. These parameters are specic for each material.

For small deformations, the innitesimal term is neglectful: the linear model gives a good approximation of the problem.

Without the inertial term and under specic assumptions on the volumetric force, it is possible to give existence and uniqueness results, obtained via the application Korn theorem and Lax Milgramm Lemma.

For nite, or large deformations [3, 22], the contribution of innitesimal term ox(∥∇η^s∥)is not negligible.

Denition 3.3 An elastic material is hyper-elastic if there exists a stored energy function W :I×Ω^s×M⁺

3 −→R such that Σ = ∂W

∂E (3.11)

The derivative is called the Gateaux derivative. The denition includes also elastic material [7].

In the simulations, we consider a model where the stored energy expression is expressed by Neo-Hookean model. It reads:

W =C₁₀(I₁−3) +D₁(J−1) (3.12) where

C10 and D1 are parameters that depend on material: D1 weights the eect of compressibility along the motion, while C10 weights the deformation of the body.;

J :=detF I1 =J²³I1.

The dissipative damping eects are taken on count in the model by considering classical Rayleigh damping approach [4]. The equilibrium law reads:

(27)

(28)

D^ϕf

Dt (x, t) := ∂f(Y, t)

∂t

Y=ϕ⁻¹_t (x)

= ∂f(x, t)

∂t x

+w(x, t)·∇ϕ(x, t) , x ∈Ω^f_t , Y∈Ω^f₀

and

w(x, t) := ∂ϕ(Y, t)

∂t

Y=ϕ⁻¹(x,t)

= ∂x

∂t , x ∈Ω^f_t, Y∈Ω^f₀. (3.13) It is convenient to write in this framework the RNG k−ϵ model, Equation 4.1.

It reads:

ρ^fD^ϕU^f

Dt +ρ^f[(U^f −w)· ∇]U^f =−∇P^f +(

µ^f +µt

)∆U^f − 2

3∇k (3.14a)

∇ ·U^f = 0 (3.14b)

ρ^fD^ϕk

Dt +ρ^f(U^f −w)· ∇k = (µ^f + µt

σk

)∆k−ρ^fϵ+Pk (3.14c) ρ^fD^ϕϵ

Dt +ρ^f(U^f −w)· ∇ϵ= (µ^f +µt

σϵ

)∆ϵ−C2ϵρ^fϵ² k + ϵ

kC1ϵPk−Rϵ

(3.14d)

We consider the smooth function phi obtained as the solution of the harmonic problem , i.e. ∆ϕ = 0. The boundary conditions will be specied for each problem in the next chapter.

(29)

(30)

Figure 3.4: Numerical approach for FSI [17]

The physics can be loosely- or strongly-coupled. In the weak approaches, the equations for dierent elds of two physics are solved separately. For every time step, the increment of uid problem is computed, the uid pressure forces recorded and, eventually, the structural problem is solved for the next time step, Figure 3.5.

(31)

Figure 3.5: Weak coupling from [27]

However,the problem is not unconditionally stable. The stability properties do depend on densities of physics and geometry [5, 11]. For example, body motion in uid shifts some portion of volume of the uid. As a consequence, the uid adds a certain mass to the system. The added mass eect is very relevant when the uid and structure densities are comparable. From numerical point of view, we observe this eect on a limitation of the time step.

On the other hand, a strong coupling scheme is presented when the physics are treated as elements of a single dynamical system, and all of the governing equations are integrated simultaneously, and interactively in the time-domain. A sketch of algorithm is given in the Figure 3.6.

20

(32)

Figure 3.6: Strong coupling from [27]

(33)

Chapter 4 Mathematical modeling of FSI industrial applications

In this chapter we focus on the industrial FSI applications.

First, we mathematically present pipe ow problems and we give the weak formulation of coupled FSI pipe ows. For these applications, we analyse the energy of the whole coupled system.

In the second part, we describe a mathematical model for a specic peristaltic pump and we write the weak formulation of the problem.

4.1 Deformable pipe problems

Let us consider a 3D domain, as Figure 4.1.

Figure 4.1: 3D cylindrical model for pipe ow

The pipe has an initial steel non-movable part, Ω, that links the pipe with the 22

(34)

(35)

Coupling :











w= ∂η^s

∂t on Γ^{f si}_t ∪Γ U^fn=wn on Γ^{f si}_t ∪Γ Jϕ

[(

−P^f − 2 3k

)

I+ (µ^f +µt)(∇U^f + (∇U^f)^T) ]

∇ϕ⁻¹n= ˜T^sn on Γ^{f si}_t ∪Γ

BCs:











U^f = ˜g^f , k = ˜k , ϵ= ˜ϵ on Γin

−P^f − 2

3k+ 2(µ^f +µt)∂U^fz

∂z = 0 , U^fx =U^fy = 0 , ∂k

∂z = 0 , ∂ϵ

∂z = 0 on Γout

U^f ·n= 0 , ∇U^fn =0 , ∇k·n= 0 ,∇ϵ·n = 0 onΓ¹_sym∪Γ²_sym η^s =0 on Γ¹∪Γ²

∂k

∂x = ∂k

∂y = 0, ϵ= u³_τ

κ˜y onΓ^{f si}∪Γ (µ^f +µτ)(∇U^f +∇U^{f T})n=−ρuτ

˜

y [(U^f −w)−((U^f −w)·n)n] onΓ^{f si}_t ∪Γ Tn˜ = 0 on Γ^est

whereg˜^f,k˜ and ˜ϵ are dened functions, we will discuss in the next section. The ALE map presents the following constraints:

ϕn=0 on Γ_in∪Γ_out∪Γ¹_sym∪Γ²_sym We set as initial condition at timet= 0s:

1. U^f =0,P^f = 0,k = 0,ϵ= 0 inΩ^f_t; 2. η^s=0 inΩ^s∪Ω

4.1.1 Inlet conditions: Womersley ow

The problem resolution depends on assigning inlet boundary conditions for˜g^f,k˜ and ˜ϵ. The codes dene a built in functions fork˜ and ˜ϵ that depend ong˜^f:

˜k=cbc

g˜^f

2 ˜ϵ=Cµ

˜k^3/2

lmix onΓin

wherecbc and Cµ are empirical constants and lmix is the mixing length of eddies, dened by Prandtl as "the diameter of the masses of uid moving as a whole in each individual case" [24].

24

(36)

(37)

In Ω˜, no uid structure interaction problem is considered because the domain does not deform. The uid ow is modeled by Eqs.(3.1). For the symmetry of the model, we pass from Cartesian coordinates to cylindrical ones (r, θ, z) and we assume u^f_θ=0 and the components u^f_z and u^f_r do not depend on θ. The NS model reads:

1 r

∂u^fr

∂r + ∂u^fz

∂z = 0 ρ^f

(∂u^fr

∂t +u^fr

∂u^fr

∂r +uz

∂u^fr

∂z )

=−∂p^f

∂r +µ^f∆u^fr

ρ^f (∂u^f_z

∂t +u^f_r∂u^f_z

∂r +u^f_z∂u^f_z

∂z )

=−∂p^f

∂z +µ^f∆u^f_z

On boundaryΓin, dened byz = 0, the ow has direction normal to the surface, u^f_r(r,0, t)=0, and this holds ∂u^f_r

dt (r,0, t)=0 and ∂u^f_r

∂r (r,0, t)=0.

From the last condition and mass conservation equation, we have that:

∂u^fz

∂z (r,0, t) = 0 In particular

ρ^f∂u^fz

∂t (r,0, t) =−∂p^f

∂z (r,0, t) +µ^f1 r

∂

∂r (

r∂u^fz

∂r )

(r,0, t) +µ^f∂²u^fz

∂z² (r,0, t) By experimental observations in [2, 20, 30] 1

r

∂

∂r (

r∂u^f_z

∂r )

(r,0, t)>> ∂²u^f_z

∂z² (r,0, t), that leads to set ∂²u^fz

∂z² (r,0, t) = 0. In order to compute u^f_z(r,0, t), we need to compute ∂p^f

∂z (r,0, t). We assume that ∂p^f

∂z (r,0, t) =−1 c

∂pd

∂t (r,0, t) +∂p0

∂z (r,0, t), wherep0 :=P(0), pd(t) :=P(t)−p0 and cis a constant computed [2]

c= v u u t

K ρ(1 + DK

eE ) where

K is the bulk modulus of uid;

D pipe diameter;

26

(38)

e thickness of pipe;

E is pipe wall material Young modulus.

∂p0

∂z (0) is computed by Darcy-Weisbach equation: ∂p0

∂z (0) = fDρ^fV_∞²/(2D), wherefD is the friction factor, computed by Moody diagram in Figure 4.5,V is the mean velocity magnitude.

Figure 4.5: Moody Diagram. By increasing Re, the friction factor decays in a linear way for laminar pipe ows, while for turbulent ows the friction factor exponentially decays.

In each test casesP(t)∈C(I), whereI is the time interval, and it is possible to consider a Fourier series expansion for pd.

pd(t) =

+∞

∑

n=−∞;n̸=0

pnexp(inωt) whereω = 2π

T and T is the period of the pressure signal.

We can suppose that velocity atΓindepends on radial coordinates and time. The inlet velocity can be represented in the following way :

˜

g^f(r, t) =

+∞

∑

n=−∞

Cn(r)exp(inωt)

whereCn(r)are Fourier coecients space dependent on the ˜g^f. By substituting terms in the last equation, one can get:

+∞

∑

n=−∞

inρ^fωCn(r)exp(inωt) =

+∞

∑ ipn

c ωexp(inωt)− ∂p0

∂z +µ^f

+∞

∑ 1 r

d dr

( r∂Cn

∂r )

exp(inωt)

(39)

Becauseexp(inωt)is a collection of independent functions, we divide the equation into:

n = 0, static pressure equation:

0 = −∂p0

∂z +µ^f1 r

∂

∂r (

r∂C0

∂r )

n ∈Z− {0}, dynamic pressure equation:

inρωCnexp(inωt) = ipn

c ωexp(inωt) +µ1 r

∂

∂r (

r∂Cn

∂r )

By assuming the no slip and axial-symmetric conditions, we compute the solution for the ODEs.

C0(r) = fDρV²

8µD (R²−r²) C_n(r) = Pn

ρc (

1− J0(rαni^3/2) J0(Rαni^3/2)

)

where α_n =√

nωρ^f/µ^f, ∀n∈Z− {0}

The coecientαnis the Womersley number, that represents the ratio between the transient inertial forces and the viscous forces. The function J0 is the Bessel function:

J0(x) =

∞

∑

m=0

(−1)^m m!(m+ 1)!

(x 2

)2m

(4.2)

4.1.2 Pressure defective data

We have used in the simulations the Womersley approach in order to compute the velocity rate when the average transient pressure is known, as in [30]. The Womersley ow is based on empirical assumptions of the velocity eld.

In the following analysis, the assumptions are removed and we get, in a more rigorous mathematical framework, a new formulation of the problem. The problem is dened in the junction part: no deformation domain problems are solved.

Furthermore we suppose that in the domainΩ, Figure 4.6, the uid behaviour is˜ laminar.

28

(40)

Figure 4.6: Domain denition for the Defective inlet pressure problem in pipe ow with radius R. The pipe wall, Γwall, is not movable.

The uid ow presents the following boundary conditions, as [13]:

Defective condition on Γin :







P(t) := 1

|Γ_in|

∫

Γin

p^f(x, t)dS (4.3)

Condition on Γout : {

−p^fn+µ^f∇u^fn=0 (4.4)

Condition on Γwall : {

u^f =0 (4.5)

Condition on Γ¹_sym∪Γ²_sym : {

u^f ·n= 0 , ∇u^fn=0 (4.6)

The conditions are defective because neither Inlet velocity nor the Inlet pressure functions are dened. The conditions describe the mean surface pressure forΓin. In order to compute the uid ow, we make the assumption, as in [14] :

(41)

Assumption 4.1 The normal stress on Γin is aligned with the normal direction and it is constant over the section, that is:

−p^fn+µ^f∇u^fn=cn on Γin (4.7) with c∈R

The uid ow is described inΩ˜ by the generalized Stokes problem [15]

αSTu^f −µ^f∆u^f +∇p^f =0 (4.8a)

∇ ·u^f = 0 (4.8b)

with αST ≥ 0 is a given parameter, related to the choice of the time advancing scheme.

We denote by (u^f(c);p^f(c)) a solution of 4.8 with conditions 4.4 - 4.7. In order to satisfy 4.3 , we need to dene a constrained problem to determine the parameter c.

We test the Eq.(4.12) a on the domain Ω˜ with u∈Vdiv :={

v ∈(H¹( ˜Ω))³ : ∇ ·v = 0, v =0 on Γ_wall,v·n = 0 onΓ¹_sym∪Γ²_sym} and δp∈H¹( ˜Ω). The weak formulation of the problem reads:

αST

∫

Ω˜

u^f·δu dV+µ^f

∫

Ω˜

∇u^f :∇δu dV−

∫

Γin

cδu·n dS−

∫

Ω˜

δp∇·u^f dV= 0 (4.9)

At a xed time stept˜∈I, we dene the following function :

J_p(s) := 1 2



 1

|Γin|

∫

Γin

s dS−P_ist^f





2

(4.10)

wheres ∈L²( ˜Ω) and P_ist^f :=P(˜t).

We look for the minimum of the Jp constrained by Equation 4.9 in order to solve the defective boundary problem. From these considerations we dene the following Lagrangian functional, see [15] :

L(w, s;λ_w, λs;η0) :=Jp(s) +αST

∫

Ω˜

w·λ_w dV+µ^f

∫

Ω˜

∇w:∇λm dV−

∫

Γin

η₀δλ_w·n dS−

∫

Ω˜

λ_s∇ ·w dV

30

(42)

The stationary points of L are described by:

GivenP_ist^f ∈R, nd c∈R, u(c)∈Vdiv, p(c)∈H¹( ˜Ω), λ_u ∈Vdiv andλp ∈H¹( ˜Ω), such that, for all v ∈Vdiv, q ∈H¹( ˜Ω) and ν ∈R,

Derivative with respect toλ_w











∂L

∂v = αST

∫

Ω˜

u·v dV+µ^f

∫

Ω˜

∇u:∇v dV−

∫

Γin

cv·n dS= 0

Derivative with respect toλs

{∂L

∂q =−

∫

Ω˜

q∇ ·u dV= 0

Derivative with respect to w {∂L

∂v =α_ST

∫

Ω˜

v·λ_u dV+µ^f

∫

Ω˜

∇v:∇λ_u dV−

∫

Ω˜

λ_p∇ ·v dV

Derivative with respect tos







∂L

∂q =



 1

|Γin|

∫

Γin

p dS−P_ist^f



 1

|Γin|

∫

Γin

q dS−

∫

Ω˜

q∇ ·λ_u dV=

Derivative with respect toη0

{∂L

∂ν =−

∫

Ω˜

νλu·n dS= 0

The stationary point provides a solution of the defective problem. The problem is well posed for the choise of the functional spaces, we have made, and it admits a unique solution, see [15].

4.1.3 Weak formulation of the turbulent uid and structure models

Before presenting the numerical results, let us write the weak formulation of the model. In order to simplify the discussion, it is convenient to introduce the variable P˜^f :=P^f − 2

3k. Equation for k−ϵ:

Let us consider δk ∈ Vk := {

v ∈H¹(I ×Ω^f₀)

v = 0 on Γin and test the

(43)

equation for k on Ω^f_t. We get:

∫

Ω^f_t

ρ^fD^ϕk

Dt δk dV+

∫

Ω^f_t

ρ^f(U^f −w)∇kδk dV=

−(µ^f + µt

σk

)

∫

Ω^f_t

∇k· ∇δk dV−

∫

Ω^f_t

ρ^fϵδk dV+

∫

Ω^f_t

Pkδk dV (4.11) For the equation for ϵ, we dene test function

δϵ∈Vϵ :={

v ∈H¹(I ×Ω^f₀)

v = 0 onΓin∪Γ^{f si}_t ∪Γ} . The weak formulation for the problem reads:

∫

Ω^f_t

ρ^fD^ϕϵ

Dt δϵdV+

∫

Ω^f_t

ρ^f(U^f −w)∇ϵδϵ dV=

−(µ^f +µt

σϵ

)

∫

Ω^f_t

∇ϵ· ∇δϵdV−

∫

Ω^f_t

C_2ϵϵ²

kδϵ dV+

∫

Ω^f_t

ϵ

kC_1ϵP_kδϵ dV−

∫

Ω^f_t

R_ϵδϵdV (4.12)

RANS model

We consider the test function δu ∈ {

v ∈(H¹(I ×Ω^f₀))³

v = 0 onΓ_in, v·n = 0 onΓ¹_sym∪Γ²_sym, v−(v · n)n = 0 onΓout} and q∈L²(Ω^f_t). We obtain:

∫

Ω^f_t

ρ^fD^ϕU^f

Dt δu dV+

∫

Ω^f_t

ρ^f(U^f −w)∇U^fδu dV=

∫

Ω^f_t

P˜^f∇ ·δu dV

−(µ^f +µt)

∫

Ω^f_t

[∇U^f +(

∇U^{f T})]

:∇δu dV+

∫

Γ^{f si}_t

[ ˜P^fI+ (µ^f +µt)(∇U^f +∇U^{f T})]δu·n dS

(4.13)

∫

Ω^f_t

∇ ·UqdV= 0 (4.14)

Structural part

We consider the model with damping (the model without any damping properties is given by settingαandβequal to 0). We consider test functions δv ∈ Vv := {v ∈(H¹(I×Ω^s))³| v = 0 onΓ¹ ∪Γ²}. We test structural equations and get:

32

(44)

∫

Ω^s

ρ^sd²η^s

dt² δv dV+

∫

Ω^s

ρ^sαdη^s

dt δv dV+

∫

Ω^s

( ˜T^s+βdT˜^s

dt ) :∇δv dV−

∫

Γ^{f si}0

T˜^sn·δv dS+

∫

Γ^{f si}0 ∪Γ^est

βdT˜^s

dt )n·δv dS= 0 (4.15) By coupling the uid and structure problems and considering that δv = δu on Γ^{f si}_t , we get the following results:

∫

Ω^f_t

ρ^fD^ϕU^f

Dt δu dV+

∫

Ω^f_t

ρ^f(U^f −w)∇U^fδudV =

−(µ^f +µt)

∫

Ω^f_t

[∇U^f +(

∇U^f)T]

:∇δu dV+

∫

Ω^f_t

P˜^f∇ ·δu dV+

∫

Ω^s

ρ^sd²η^s

dt² δv dV+

∫

Ω^s

ρ^sαdη^s

dt δv dV+

∫

Ω^s

( ˜T^s+βdT˜^s

dt ) :∇δv dV−

∫

Γ^{f si}0 ∪Γ^est

βdT˜^s

dt n·δv dS

(4.16)

The formulation of ALE problem is non conservative, it means the time derivative of test function along the ALE mapping is not null. With a suitable choice of test function is possible to give conservative formulation of the problem. For more details see [13]. The commercial codes for our simulations use non conservative ALE approach: test functions are not constant along ALE mapping.

4.1.4 Energy estimates

We derive the energy balance equation for the FSI problems with a structural part that presents damping properties. We recall the results from [12] by including damping properties of the structural part. The energy equations should include the uid viscosity and the damping structural eects.

Before introducing the energy equation, we state the following theorem:

Theorem 4.1 (Reynolds transport theorem in ALE description) Let ϕ be the

A Numerical study of fluid flow in peristaltic pumps and polymeric elastic pipes

Acknowledgments

Contents

Chapter 1 Introduction

1.1 Objectives

1.2 Research methodology

1.3 Organization of the study

Chapter 2

Mathematical Preliminaries

Chapter 3

Mathematical modeling of FSI problems

3.1 Fluid subproblem

3.1.1 Generalities on Turbulence models

3.1.2 Mathematical analysis for k − ϵ model

3.2 Structural problem: small and nite deforma- tions

Chapter 4

Mathematical modeling of FSI industrial applications

4.1 Deformable pipe problems

4.1.1 Inlet conditions: Womersley ow

4.1.2 Pressure defective data

4.1.3 Weak formulation of the turbulent uid and structure models

4.1.4 Energy estimates