A Space-Time Finite Element Approach to the Numerical Simulation of Vascular Fluid-Solid Interaction

(1)

Delft Aerospace Computational Science 0000 0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111

000000000

111111111

ENGINEERING MECHANICS July 2008

A Space-Time Finite Element Approach to

the Numerical Simulation of Vascular

Fluid-Solid Interaction

E.J. Vlijm, E.H. van Brummelen

e.h.vanBrummelen@tudelft.nl

(2)

in whole or in part without the written permission of the author.

TUD/LR/EM

Kluyverweg 1, 2629HS Delft, The Netherlands

P.O. Box 5058, 2600 GB Delft, The Netherlands

Phone: +31 15 278 5460 Fax: +31 15 261 1465

ISSN 1574-6992

(3)

Summary

Numerical studies of cardiovascular diseases like arteriosclerosis have gained increasing attention the last decade. The modeling of blood, blood vessel and their coupling, shows to be a challenging problem. In this thesis a two-dimensional model has been constructed and its behaviour has been investigated. We model the blood vessel by a linear elastic continuum and blood by an incompressible Newtonian fluid. The coupled system is solved with a space-time discontinuous Galerkin method. The solution is allowed to be discontinuous over the discrete time level, but continuity in space is enforced. We calculate the movement of the fluid domain by the equations of linear elasticity, which results in a so-called ’pseudo solid’ mapping. The coupled system is solved with a strongly-coupled solution method. The subiteration procedure occurring within each time step is accelerated by relaxation. Relaxation shows to be necessary when the density of the structure approaches the density of the fluid. The number of iteration steps needed within a time step shows to be quite large and the need for better acceleration becomes apparent. Although the relaxation technique offers some improvement of the subiteration procedure, better techniques are necessary for three-dimensional simulations.

(4)

Introduction

Arterial vessels perform the vital task of supplying blood to every part of the human body. Every second the vessels need to deal with a pressure wave originating from the contraction (systole) and the relaxation (diastole) of the heart. This results in a dynamic motion of the blood as well as the vessel. The vessel’s material must withstand these conditions and is therefore made out of very elastic tissue layers. It can happen that the blood vessel hardens under certain conditions, which is known as atherosclerosis. Atherosclerosis is one of the major cardiovascular diseases caus-ing death in the western world today. It is therefore an important field of study. Hemodynamic (blood dynamic) studies have been frequently performed to increase insight in the origin of diseases like atherosclerosis and to find effective treatments. The modeling of blood and coupling with an interacting vessel, however, is a field that has been studied only since the last couple of years.

This fluid-structure-interaction (FSI) problem requires not only a proper blood flow and blood vessel model, but also adequate coupling of the two. The coupling and the models both influence the behaviour of the coupled system. In this thesis the blood vessel is modeled by a linear elastic solid continuum. Though this simplification omits the visco-elastic behaviour of the vessel, it gives a good first indication of its inter-action properties with blood flow. We model blood by an incompressible Newtonian fluid. A Newtonian fluid assumes a linear dependence between stress and strain rate. The coupling of these models is accommodated in two distinct conditions, viz., the kinematic and dynamic conditions. The kinematic conditions enforce displacements and velocities of fluid and structure to be equal at the interface. The dynamic condi-tion enforces the traccondi-tion of fluid and structure to be equal at the interface.

The resulting equations for fluid and structure are given in Eulerian (spatial) and Lagrangian (material) formulations, respectively. To account for the deformation of the fluid domain, the Arbitrary-Lagrangian-Eulerian (ALE) formulation for the fluid equations is often introduced (see e.g. [4], [10] and [15]). In this thesis we employ a space-time analogy, however, which embeds the ALE-formulation.

To solve the fluid and structure equations we employ the Space-Time Discontinuous Galerkin Finite Element Method (DGFEM). Where the approximate solution is al-lowed to be discontinuous over the discrete time level, but is continuous in space. The

(7)

space-time approach will result in implicit satisfaction of the Geometric Conservation Law, which must be taken into account separately when semi-discrete approaches are employed.

In this thesis we constructed a two-dimensional model employing the existing MPF code (see appendix A). To be able to capture geometric irregularities, we implemented the option to employ unstructured meshes. The available implementation restricted us to employ equal order interpolation functions for both pressure and velocity spaces. Hence, a stabilized formulation for the incompressible Navier-Stokes equations had to be employed and needed to be implemented, including the option to use a deformable domain. The motion of the fluid domain has been modeled by the equations of linear elasticity. The structure model was not present in the existing code and had to be implemented as well. The models have been tested separately.

The coupled system is solved using subiteration. To accelerate the iteration procedure we employ relaxation. We observe that relaxation is necessary when blood flowing through a compliant vessel is simulated. The so-called added mass effect appears to play a considerable role in the FSI-problem. That is, if the fluid density is of the same magnitude as the structure density, the behaviour of the structure is influenced by the fluid considerably. We have analyzed the model for two different test cases. In the first test case we enforce steady inflow condition and in the second test case we enforce a pulsatile inflow condition.

(8)

Chapter 2

Hemodynamics

Modeling blood flow including its interaction with a vessel proves to be one of the most challenging fluid-structure-interaction problems. Not only the coupling of the fluid and structure models introduces complexity, the separate fluid and structure models are very complicated. The mechanical behaviour of blood and blood vessel is very complex and accurate models contain for instance dependencies on age and lifestyle. In this thesis we restrict ourselves to baseline models to give us a good first impression of real-life blood and blood vessel behaviour.

2.1 Blood vessel model

The blood vessel is a thin elastic material which is capable of distending during the systole and diastole. Although various types of blood vessels exist, they all have the same basic structure. One can distinguish different layers of material inside the ves-sel. Basically three layers occur: the intima, media and adventitia (see figure 2.1). The intima is a thin layer which is closest to the lumen (the fluid domain). The layer

Figure 2.1: Basic structure of a blood vessel

consists of the one-cell layer thick endothelial and the elastic interna. The media layer consists of muscle cells and fibers. The adventitia is the outer most layer and consists

(9)

of connective tissue. All three layers behave differently under various circumstances, like stress or a high proportion of cholesterol in the blood. One can imagine that it is hard to integrate the mechanical properties of all three layers into a single model. Though recent studies have been performed on modeling the vessel as an anisotropic composite material including fibers and their dispersities (see e.g. [7]), the complexity of these models demand a large expertise and are therefore not often used in model-ing fluid-structure-interaction problems. A simplification of the vessel which is often applied, is modeling the vessel as a thin walled structure, like the linear Koiter shell model (see e.g.[4]). These models are effective if the vessel is thin. However when considering stenosis, shell type models clearly fail to include the geometric irregu-larities like local thickening of the vessel. To capture these geometric irreguirregu-larities, we need to model the vessel as a solid continuum. To limit the complexity of the model we restrict ourselves to a basic linear elastic solid continuum which is a good first approximation for a real vessel. The material properties of blood vessels depend on the location in the vascular system. Table 2.1 gives the material properties and reasonable values which occur in the human body.

density 0.8 - 2.0 g/cm3

stiffness modulus 1E6 − 1E7 dyn/cm2

Poisson ratio 0.4 - 0.5

Table 2.1: Properties of the blood vessel (dyn = gcm/s2₎

2.2 Blood flow model

An average human body contains 5-6 liters of blood, which can be subdivided into several components. Roughly fifty percent of the blood’s volume consists of blood cells, of which ninety-five percent is occupied by red blood cells and five percent by white blood cells and platelets (responsible for clodding). The other fifty percent volume is occupied by plasma, which is made out of water for ninety percent and ten percent other substances. Again these numbers change under certain circumstances. During an infection, for instance, the number of white blood cells increases dramati-cally. To model these features would be almost impossible and therefore a simplified model will be employed. We model blood as a Newtonian fluid, which states a linear dependence between stress and strain rate. Although blood in reality is an anisotropic visco-elastic fluid, the Newtonian fluid model gives a good impression of its mechan-ical behaviour. Reasonable values for the density and viscosity of blood are given in table 2.2, where P is poise (g/cms). We observe that the difference between the density of the blood and the density of the vessel is small. This will result in the added mass effect. The responsive behaviour of the vessel is considerably influenced by the fluid density. Hence the coupled system shows very different behaviour from cases where the fluid density is much smaller than the structure density, such as air flowing over an airplane wing.

(10)

density 0.8 - 1.2 g/cm3

viscosity 0.01 - 0.05 P Table 2.2: Properties of blood

(11)

Chapter 3

Governing equations

In this chapter the governing equations for the fluid and the structure are presented. Before doing so, the two basic ways of describing a continuum are explained. The structure equations are generally given in Lagrangian description, while for the fluid equations the Eulerian description is preferred. When considering a deforming fluid domain however, the Eulerian description has its shortcomings. The classical ap-proach to overcome these shortcomings is to introduce a description which can be arbitrary Lagrangian or Eulerian (ALE-formulation). In this thesis we will employ a space-time analogy, however, in which the ALE-formulation is embedded.

3.1 Kinematic description of a continuum

Ω Ω x x_t _t x Lagrangian Eulerian χ ALE Ω χ

Figure 3.1: Different continuum formulations

There are two distinct ways to describe a continuum; the Eulerian and Lagrangian description. To illustrate this, let Ω ⊂ ℜ3 _{be a reference configuration of a certain}

(12)

body. This body is allowed to move and deform in time and at time t the body is denoted by Ωt ⊂ ℜ3. We also define an initial stress-free reference configuration

Ω0⊂ ℜ3and let it coincide with the reference configuration Ω. A one-to-one mapping,

sufficiently smooth, from the reference configuration Ω to Ωt can be defined:

χ : Ω × [0, T ] → Ωt (3.1)

If we denote the material points in the reference configuration by x, the position of a material point at time t w.r.t. the reference configuration is given by:

xt= χ(x, t) (3.2)

With these definitions we can define the deformation, the velocity field, deformation gradient and its determinant as:

u(x, t) = χ(x, t) − x, v = ∂χ ∂t, F =

∂χ

∂x, J = det F (3.3) We can define any field quantity φ with values in a certain space V w.r.t. the config-uration Ωt by:

φ = ˜φ(xt, t) : Ωt× [0, T ] → V

This is called the Eulerian description (sometimes called the spatial description). The derivatives in this configuration can be defined as:

∂φ ∂t := ∂ ˜φ ∂t, ∇φ =˜ ∂φ ∂xt := ∂ ˜φ ∂xt (3.4)

Note that this formulation is preferable when the body undergoes very large dis-placements and the field quantities are specified in terms of velocities rather than displacements, which is the case in fluids.

On the other hand we can define the field quantities of φ on the reference configuration Ω by:

φ = ˆφ(x, t) : Ω × [0, T ] → V

This is known as the Lagrangian description. The derivatives in the reference config-uration are defined as:

dφ dt := ∂ ˆφ ∂t, ∇φ =ˆ ∂φ ∂x := ∂ ˆφ ∂x (3.5)

This description is preferred when formulating conservation laws for a solid contin-uum. The Lagrangian and Eulerian description are related to each other through the following relation:

ˆφ(x, t) = ˜φ(χ(x, t), t) (3.6) The time derivative of φ in the Lagrangian and Eulerian formulations are related as follows: d ˆφ dt = ∂ ˜φ ∂χ dχ dt + ∂ ˜φ ∂t = ˜∇ ˜φ · v + ∂ ˜φ ∂t (3.7)

Observe the convective term appearing in the Eulerian formulation.

In the case of a deforming fluid domain the Eulerian description gives problems since it describes material points from a spatial perspective. To overcome these problems, often another referential domain is introduced which can be arbitrary Lagrangian or Eulerian (see figure 3.1). This is the so-called ALE- (Arbitrary Lagrangian- Eulerian) formulation and includes the domain movement.

(13)

3.2 Structure equations

We consider the fluid and structure equations on certain domains Ωf and Ωs, both

subsets of ℜ3 _{with boundaries dΩ}

s and dΩf, the interface Γ and Ωf ∩ Ωs = Ø.

See figure 3.2 for details. We employ the subscript f in the fluid equations and the subscript s in the structure equations. The equations for the structure are given in

Ωs Ωf Ωs ∂Ωs ∂Ωf ∂Ωs Γ Γ

Figure 3.2: Fluid and structure domains Lagrangian formulation and the momentum equation reads:

ρ0 s∂

2_u s

∂t2 − ∇ · S = b on Ωs (3.8)

with the associate boundary conditions:

us= g on ∂Ωs (3.9)

n · S = h on Γ (3.10)

Where S is the first Piola-Kirchhoff stress tensor, b is the body force and ρ0 s is the

density and us is the displacement vector w.r.t. the referential domain. The first

Piola-Kirchhoff stress tensor is defined on the initial reference configuration and is related to the Cauchy stress tensor τ by the transformation:

S = (det F )τF−T _(3.11)

The conservation of mass written in Lagrangian formulation states constant density. The linear elastic constitutive model employed in this thesis expresses interdepen-dence between the stress and strain:

S = Cε (3.12)

where C denotes the fourth order elasticity tensor and the strain tensor ε is defined as:

ε = 1₂(∇us+ ∇uTs) (3.13)

We assume the material to be isentropic and can therefore rewrite the stress-tensor:

S = 2µε + λ tr(ε)I (3.14)

where λ and µ (also known as the two constants of Lam´e) are defined as:

λ = _{(1 + ν)(1 − 2ν)}νE (3.15) µ = E

2(1 + ν) (3.16)

where E is the stiffness modulus and ν Poisson’s ratio. We assume the body to be homogenous which results in constant values for ρ0

(14)

3.3 Fluid equations

The incompressible Navier-Stokes equations are given in Eulerian formulation and read:

ρf∂u_∂tf + ρfuf· ∇uf− ∇ · σf− f = 0 on Ωf (3.17)

∇ · uf = 0 on Ωf (3.18)

with the associated boundary conditions:

uf= g on ∂Ωf (3.19)

uf= h on Γ (3.20)

Where ρfis the density and f is a body force which is assumed zero in this thesis. The

energy equation is discarded since we are only interested in the mechanical aspects of the fluid. Note that uf in these equations is the velocity, whereas in the structure

equations us comprises the displacement. For Newtonian fluids we can write the

stress tensor σf as a function of pressure and strain rate in the following manner:

σf = µ(∇uf+ ∇uTf) − pI (3.21)

Where µ is the dynamic viscosity of the fluid. As discussed in section 3.1, the Eulerian formulation fails to include movements of the domain. To solve this problem the ALE-formulation can be employed. The derivatives of the Eulerian formulation can be rewritten w.r.t. an arbitrary reference domain (see figure 3.2). One obtains a form of equations 3.17 like:

ρfJ˜∂u_∂tf + ρfJu˜ f· ∇(uf− w) − ˜J∇ · σf− ˜Jf = 0 on Ωf (3.22)

∇ · (uf− w) = 0 on Ωf (3.23)

where w is the velocity of the domain and ˜J the Jacobian of the transformation ¯χ from the domain Ωtto the arbitrary reference domain ¯Ω, see figure 3.1.

The ALE-formulation is the standard formulation when employing a semi-discrete approach to solve the fluid equations. However in this thesis we use a space-time analogy which inherently accounts for the deformation of the domain. The ALE-formulation is ”imbedded” in the space-time ALE-formulation and we can therefore use the standard Eulerian formulation (3.17). Musad and Hughes (see [11]) show that if the time domain is mapped onto a specific computational domain the space-time formulation reduces to the standard ALE-formulation. This is exactly the case in this thesis and this aspect will be explained further in section 4.4.

3.4 Interface conditions

The fluid and structure are coupled through two distinct conditions, viz., the kine-matic and dynamic conditions. These conditions ensure continuity and conservation of momentum over the interface.

(15)

3.4.1 Kinematic conditions

In a fluid-structure-interaction problem one can distinguish two kinematic conditions. The first condition ensures the fluid particles to have the same velocities as the struc-ture particles at the interface, which is often referred to as the no-slip condition. The second condition ensures the fluid and structure domain to follow the interface. The two conditions can be formulated as:

˙us= uf on Γ (3.24)

us= ud on Γ (3.25)

3.4.2 Dynamic condition

To ensure conservation of momentum at the interface, the dynamic condition must be satisfied. The dynamic condition enforces the traction of the fluid to be equal to the traction of the structure at the interface. The condition reads:

n · τ = n · σf on Γ (3.26)

We must note that the stress tensors of both structure and fluid in this condition are the Cauchy stress tensors. The structure formulation is defined on the initial reference domain however, hence the tensor must be transformed using relation 3.11.

3.5 Non-dimensionalized equations

We scale the fluid and structure equations using the following expressions: ˆuf =u_Vf, ˆus=u_Ls, ˆx = x_L, ˆt= tV_L, ˆρf =ρ_ρf₀ f, ˆρs= ρs ρ0 s, (3.27) ˆp = p ρ0 fV2, ˆS = S E, ˆσf = σf ρ0 fV2, (3.28) Where V , L, E, ρ0

s and ρ0f denote the characteristic velocity, characteristic length,

characteristic stiffness and characteristic densities. Note that the characteristic ve-locity in the structure equations is assumed to be equal to the characteristic veve-locity of the fluid.

Structure

The structure equation rewritten with the expressions reads: ρ0 sV2 L ˆρs ∂2_ˆu_s ∂ˆt2 − E L∇ · ˆS(ˆuˆ s) = 0 on Ωs (3.29) We non-dimensionalize the equation by dividing it by ρ0sV2

L and note that ˆρs= 1:

∂2_ˆu s ∂ˆt2 − E ρ0 sV2 ˆ ∇ · ˆS(ˆus) = 0 on Ωs (3.30)

(16)

ˆu = _Lg on ∂Ωs (3.31)

n · ˆS = Eh on Γ (3.32)

Observe that there is no difference in form between the non-dimensional and dimen-sional equations, except for the scaling of the boundary conditions and the scaling of the tensor S. Equation 3.30 is characterized by one parameter if we assume homoge-neous material, namely α = E

ρ0_s_V2.

Fluid

The equations for the fluid rewritten become: ρ0

fV2

L

ˆρ∂ˆuf

∂ˆt + ˆρûf· ˆ∇ûf− ˆ∇ · ˆσf(ûf, ˆp) = 0 on Ωf (3.33) V L ˆ ∇ · ûf = 0 on Ωf (3.34)

We non-dimensionalize the equations by dividing the momentum equation by ρ0fV2

L

and note that ˆρf = 1. The continuity equation is divided by V_L and we rewrite ˆσf:

∂ûf ∂ˆt + ûf· ˆ∇ûf− ˆ∇ · 1 Re( ˆ∇ûf+ ˆ∇ûTf) − ˆ∇ˆpf = 0 on Ωf (3.35) ˆ ∇ · ûf = 0 on Ωf (3.36)

ˆuf = _Vg on ∂Ωf (3.37)

ˆuf = _Vh on Γ (3.38)

where

Re = ρf0V L

µ (3.39)

characterizes the equation.

Kinematic and dynamic conditions

The non-dimensionalized form of the dynamic condition includes a ratio of densities and α and reads:

n · ˆτ = M_αn · ˆσf on Γ (3.40)

where

M = ρ_ρf

(17)

Since we employ the same characteristic length and characteristic velocity for both structure and fluid, the non-dimensionalized form of the kinematic conditions is equal to the dimensionalized form.

For the coupled system we can identify three parameters which characterize the sys-tem, namely Re, α and M.

These non-dimensional equations will be used for the fluid-structure simulation in chapter 7.

(18)

Chapter 4

Finite-element formulations

In this chapter the finite-element formulations of the structure and fluid equations are presented. To account for the movement of the fluid domain, the domain itself is modeled by the equations of linear elasticity. The finite-element formulation of these equations is presented as well. To avoid notational confusion we employ the subscript s for the structure unknowns, f for the fluid unknowns and d for the fluid-domain unknowns.

4.1 Space-time discontinuous Galerkin

discretiza-tion

To solve the conservation equations presented in the previous chapter, we introduce a space-time discontinuous Galerkin discretization. The discretization is discontinuous in time and continuous in space. Consider a certain time interval I = [0, T ] on which we seek a solution. We subdivide the time interval into n subintervals, introducing a time slab In = [tn, tn+1] and ∆t = tn+1− tn. The spatial domain at time tn

is denoted by Ωn and can be subdivided into the corresponding fluid and structure

domains Ωf

n and Ωsn. With this notation the complete space-time domain can be

denoted by Q = Ω × I and a space-time slab by Qn= Ω × In. We should note that

the space-time domain deforms in time. The boundary of the spatial domain at time tnis denoted by ∂Ωnand this boundary, extended by a time step, by ∂Qn= ∂Ω × In.

The spatial domain can be subdivided into intervals denoted by Ωe

n, which extruded

in time gives us the space-time elements Qe

n. For notational convenience, we drop the

subscript n indicating the time slab and we subdivide the space-time elements into fluid and structure elements Qe

f and Qes. The interface between the fluid domain and

the structure domain is denoted by Γ = Qf∩ Qs. See figure 4.1 for details.

We will make frequent use of the divergence theorem to rewrite the diffusive term of the momentum equation:

Z Qn div[σ] · w dQn = Z Qn −σ : ∇w + div[σ · w] dQn

(19)

Γ Qe f Qe s In Ω t x y

Figure 4.1: Sketch of a space-time slab = − Z Qn σ : ∇w dQn+ I ∂Qn [σ · w] · n d∂Qn = − Z Qn σ : ∇w dQn+ I ∂Qn σn · w d∂Qn

4.2 Structure finite-element formulation

In this section we present the finite-element formulation for the structure as discussed by Hulbert and Hughes (see [12]). Two different formulations have been presented by Hulbert and Hughes; the single-field formulation and the two-field formulation. The first formulation approximates only the displacement field with interpolation func-tions. The velocities are defined by the temporal derivative of the displacement. This is the straightforward approach. In the case of employing linear interpolation func-tions however, the second derivatives will be zero and the inertial term will vanish. To overcome this problem we will make use of the two-field formulation which ap-proximates the displacement field and velocity field by separate function spaces. This formulation enforces the continuity of the displacement fields between two time slabs weakly via the strain-energy inner product. The velocity fields are enforced weakly between two time slabs through the standard inner product. The velocity and dis-placement fields are related to each other by enforcing the temporal derivative of the displacement to be equal to the velocity through the strain-energy inner product. Before the finite-element formulation is presented, let us define the following function spaces:

(20)

Vs= {vs∈ C0(Qs, ℜd) : vs|Qe s ∈ Pk(Qes, ℜd) ∀ Qes⊂ Qs} (4.1) Ws= {ws∈ C0(Qs, ℜd) : ws|Qe s ∈ Pk(Qes, ℜd) ∀ Qes⊂ Qs} (4.2) Us= {us∈ C0(Qs, ℜd) : us|Qe s ∈ P k_(Qe s, ℜd) ∀ Qes⊂ Qs} (4.3) Ns= {ηs∈ C0(Qs, ℜd) : ηs|Q_se ∈ Pk(Qe_s, ℜd) ∀ Qe_s⊂ Q_s} (4.4)

where d is the dimension of the space-time domain and k the order of the interpolation and test functions. In this thesis, we set k to 1 and we restrict ourselves to problems in two space dimensions and, accordingly, d=3. We must note that the structure formulation is given in the Lagrangian form. The domain Qs on which the

finite-element formulation is defined is the initial domain extruded in time and does not deform.

The finite-element formulation of the structure equation can be written as: find ue s∈

Us and ves∈ Vssuch that ∀ wse∈ Wsand ηes∈ Ns,

S(us, vs, ηs, ws) = s(ws, ηs) (4.5) where S(us, vs, ηs, ws) = X e Z Qe s ρ0 swse· ∂v e s ∂t dQes+ X e Z Qe f ε(we s) : S(ues) dQes +X e Z Qe f ε(ηe s) : S(∂u e s ∂t − vse) dQes +X e Z Ωe s ρ0 swes(t+n) · ves(t+n) dΩes+ X e Z Ωe s ε(ηe s) : S(ues(t+n)) dΩes (4.6) and s(ws, ηs) = X e Z ∂Qe s ∂ws ∂t · t d∂Qes+ X e Z Ωe s ρ0 swes(t+n) · ves(t−n) dΩes +X e Z Ωe s ε(ηe s) : S(ues(t−n)) dΩes (4.7) We define ue s(t+n) and ues(t−n) as: ue s(t+n) = lim_ǫ→0+u e s(t + ǫ) and ue s(t−n) = lim_ǫ→0−u e s(t − ǫ)

Observe that the displacements and velocities are allowed to be discontinuous at the discrete time levels. Their continuity is enforced through the strain-energy inner product and the standard inner product.

4.3 Fluid finite-element formulation

The fluid finite-element formulation employed in this thesis, is the mixed formulation with added stabilizing terms. We will first introduce the mixed formulation, after which the stabilized formulation and the stabilization parameters will be discussed.

(21)

4.3.1 Mixed formulation

To present the finite-element formulation of the incompressible Navier-Stokes equa-tions we introduce the following function spaces:

Pf = {pf ∈ C0(Qf, ℜd) : pf|Qe f ∈ P k_(Qe f, ℜd) ∀ Qef ⊂ Qf} (4.8) Qf = {qf ∈ C0(Qf, ℜd) : qf|Qe f ∈ P k_(Qe f, ℜd) ∀ Qef ⊂ Qf} (4.9) Uf = {uf ∈ C0(Qf, ℜd) : uf|Qe f ∈ P k_(Qe f, ℜd) ∀ Qef ⊂ Qf} (4.10) Wf = {wf ∈ C0(Qf, ℜd) : wf|Qe f ∈ P k_(Qe f, ℜd) ∀ Qef⊂ Qf} (4.11)

The mixed finite-element formulation of the incompressible Navier-Stokes equations can be written as: find ue

f ∈ Uf and pef ∈ Pf such that ∀ wfe∈ Wfe, qef ∈ Qf,

F(uf, pf, wf, qf) = f(wf) (4.12) where F(uf, pf, wf, qf) = X e Z Qe f we f· ρf(∂u e f ∂t + uef· ∇uef)dQef +X e Z Qe f ε(we f) : σf(pef, uef)dQef +X e Z Qe f qe f∇ · uefdQef+ X e Z Ωe f ρfuef(t+n) · wef(t+n) dΩef (4.13) and f(wf) = X e Z Γf we f· tdΓf+ Z Ωe f X e ρfuef(t−n) · wef(t+n)dΩef (4.14)

where strain tensor ε(ue

f) is defined as 12(∇uef+ ∇uefT).

To guarantee a stable and unique solution to eq. (4.12) the velocity and pressure spaces must satisfy the so-called inf-sup condition. This condition restricts the pres-sure space and states that the problem has a unique solution if there exists a constant β such that: inf qf∈Qf_u_fsup_∈U_f b(uf, qf) ||uf|| ≥ β||q|| (4.15) Where b(uf, qf) =R_Q_fqf∇ · ufdQf.

In this thesis we employ linear function spaces for both pressure and velocity and one can find that using equal order function spaces for velocity and pressure will leave this condition unsatisfied. We ”circumvent” this condition by introducing stabilization.

4.3.2 Stabilized formulation

The stabilized formulation allows us to employ equal order interpolation function spaces. The most frequently used stabilization terms include the stream-line upwind/Petrov-Galerkin formulation (SUPG), the upwind/Petrov-Galerkin/least-squares (GLS) formulation on

(22)

in-compressibility condition and the pressure-stabilizing/Petrov-Galerkin (PSPG) for-mulation. The stabilized finite-element formulation including these three terms be-comes: find ue

f∈ Uf and pef ∈ Pf such that ∀ wf ∈ Wfe, qfe∈ Qf,

Fstab(uf, pf, wf, qf) = fstab(wf) (4.16) where Fstab(uf, pf, wf, qf) = X e Z Qe f we f· ρf(∂u e f ∂t + uef · ∇uef) dQef +X e Z Qe f ε(we f) : σf(pef, uef) dQef+ X e Z Qe f qe f∇ · uefdQef +X e Z Qe f 1 ρf[τSUP G(ρf ∂we f ∂t + ρfuef· ∇wef) + τP SP G∇qef] ·ρf ∂u e f ∂t + uef· ∇uef − ∇ · σf(pef, uef) dQe f +X e Z Ωe f ρfuef(t+n) · wef(t+n) dΩef + nel X e=1 Z Qe f τLSIC∇ · wfeρf∇ · uefdQef (4.17) fstab(wf) = X e Z Γe f we f· t dΓef+ X e Z Ωe f ρfuef(t−n) · wef(t+n) dΩef (4.18)

The extra terms improve stability without compromising consistency, i.e. if we replace ue

f and pef by their exact counterparts, equation 4.18 is satisfied. We note that since

we employ linear interpolation functions in this thesis, the second order derivatives in the momentum stabilization term vanish.

The stabilization terms change the space against which we test the fluid equation. Hence eq. (4.18) can be seen as a Petrov-Galerkin formulation. On˜ate [9] showed that this stabilized mixed formulation can be seen as an extension of the infinites-imal theory. The balance equations presented in the previous chapter follow from infinitesimal theory. The elements have a finite dimension, however. If we represent the balance equations over an element using a Taylor series expansion of one order higher than those used in the infinitesimal theory, we can obtain variational forms which are similar to equation (4.18).

4.3.3 Stabilization parameters

The τ parameters appearing in the stabilization terms affect the convergence of the system to a great extent. Numerous definitions of these parameters have been intro-duced (see e.g. [9], [10]). Most of these definitions include the characteristic length of the element in their formulation which changes when the domain deforms in time. Tezduyar et al.[8] defined these parameters without an explicit dependence on the characteristic element length. The characteristic length in their formulation is auto-matically included through the gradient of the test function Na. The deformation

(23)

in the integral, as explained in the next section. Let us first define the following element-wise quantity: hUGN = 2||uef|| nXno a=1 |ue f· ∇Na| −1 (4.19) where nno is the number of nodes per element and hUGN is defined assuming a single

point integration. Since we use more integration points we take the average over all integration points. The stabilization parameters for an element e can be written as follows: τs1=_2||uhUGN_e f|| (4.20) τs2=∆t₂ (4.21) τs3=h 2 UGN 4ν (4.22) τSUP G= 1 τ2 s1 + 1 τ2 s2 + 1 τ2 s3 _−1/2 (4.23) τP SP G= τSUP G (4.24)

τLSIC= hUGN₂ ||uef||z (4.25)

where z is given as:

z = (ReUGN 3 ) if ReUGN≤ 3 1 if ReUGN> 3 with ReUGN = ||u e f||hUGN 2ν

From the definition one can observe that the parameter τP SP G(= τSUP G) is divided

into three components. The convective, temporal or diffusive term dominates the stabilization depending on the flow properties, time step size and mesh width.

4.4 Fluid domain deformation

4.4.1 Domain motion

The deflection of the structure causes the fluid domain to deform in time. To compute the fluid equations at a certain space-time slab, the domain motion of that space-time slab must be take into account. Therefore we need a model which describes the fluid domain motion at that space-time slab. Several methods have been introduced to model the motion of a deformable domain. The deformation method described by Guojun and Yue [5], for instance, calculates the nodal velocities by solving a div-curl equation. Another often applied technique is to model the nodal displacements by solving the Laplace equation. In this thesis we employ the so-called pseudo-solid

(24)

mapping technique which describes the motion of the fluid domain by the equations of linear elasticity:

∇ · S(u) = 0 on Qd (4.26)

u = g on Γ (4.27)

u = 0 on ∂Qd (4.28)

Where S is the first Piola-Kirchhoff stress tensor defined in section 3.2. The function spaces in which we seek a solution can be defined as:

Wd= {wd ∈ C0(Qd, ℜd) : wd|Qe d∈ P k_(Qe f, ℜd) ∀ Qed⊂ Qd} (4.29) Ud= {ud∈ C0(Qd, ℜd) : ud|Qe d∈ P k_(Qe f, ℜd) ∀ Qed⊂ Qd} (4.30)

The finite-element formulation of equation 4.26 becomes: find ue

d ∈ Ud such that ∀ we d∈ Wd , D(ud, wd) = 0 (4.31) where D(ud, wd) = X e Z Qe f ε(we d) : S(ued) dQef (4.32)

For a certain time slab the finite-element formulation (4.32) can be solved, resulting in the deformation of the fluid domain. This deformation can be seen as a mapping χ which maps the computational domain onto the physical domain. The computational domain is defined as the initial spatial domain extruded in time with a time step ∆t (see figure 4.2). Note that χ includes the displacements PN_i=1uiφi, where φi

physical domain computational domain

Figure 4.2: Computational domain versus physical domain

are the interpolation functions which approximate the deformation field. Since the computational domain does not deform in time, it would be convenient to solve the fluid equation on this domain, rather than on the physical domain. Using χ we can rewrite the fluid equations from the physical domain onto the computational domain. Defining the basis of the computational domain as ξi and the basis of the physical

domain as xi, the derivatives in components ξi yields:

∂ ∂ξi = ∂xj ∂ξi ∂ ∂xj (4.33)

(25)

where ∂xj

∂ξi can be determined from χ and the derivatives include a time component.

To write the fluid equations onto the computational domain we need the inverse of the above relation: _

∂xj ∂ξi  −1 ∂ ∂ξi = ∂ ∂xj (4.34)

Using this relation and the Jacobian of the transformation, we can rewrite the fluid equations from the physical domain onto the computational. The left and right-hand sides of the fluid variational formulation on the computational domain become:

˜Fstab(uf, pf, wf, qf) = Z Qe f we f· ρf(∂u e f ∂˜t + uef· ˜∇uef) ˜JdQef + Z Qe f ǫ(we f) : σf(pef, uef) ˜JdQef+ Z Qe f qe f∇ · u˜ efJdQ˜ ef +X e Z Qe f 1 ρf[τSUP G(ρf ∂we f ∂˜t + ρfuef· ˜∇wef) + τP SP G∇q˜ ef] ·ρf ∂u e f ∂˜t + uef· ˜∇uef − ˜∇ · σf(pef, uef) ˜JdQef + Z Ωe f ρfuef(t+n) · wef(t+n) ˜JdΩef + nel X e=1 Z Qe f τLSIC∇ · w˜ feρf∇ · u˜ efJdQ˜ ef (4.35) ˜fstab(wf) = Z Γe f we f· t ˜JdΓef+ Z Ωe f ρfuef(t−n) · wef(t+n) ˜JdΩef (4.36) where ˜∇ and ∂

∂˜tare the derivatives in ξicomponents. The Jacobian of transformation

˜

J is included in every integral. The stabilization parameters defined in the previous section are evaluated using ˜∇ as well.

One note must be made concerning the definition of the function spaces. Since these spaces are defined on the physical domain, they need to be transformed on to the computational as well. The velocity function space for example, in which we look for a solution on the computational domain can be defined as χ ◦ U with U defined as the function space on the physical domain. To ensure we may employ linear interpolation and test functions on the computational domain like on the physical domain, the mapping χ must be a affine mapping. Since χ is defined per element by linear functions, the mapping is affine on the element domain. Hence we can employ the same function spaces defined in 4.8 for the computational domain.

The space-time approach employed in this thesis, results in implicit satisfaction of the so-called Geometric Conservation Law. This law considers the influence of the movement of the domain on the fluid solution. Simply stated, it ensures the volume swept by each face of an element to be equal to the volume change of the element. When employing a semi-discrete approach, this law must be considered explicitly. The space-time approach, however, satisfies this law automatically.

(26)

To show the equivalence of the ALE-formulation with the space-time formulation writ-ten on the computational domain, Musad and Hughes (see [11]) write out relation 4.34 per component. The same ’adjusted’ convective velocity as in the ALE-formulation becomes apparent.

4.4.2 Jacobian-based stiffening

The deformation of the fluid domain is modeled by the equations of linear elasticity as explained in the previous section. Since the deformation varies over the domain, certain elements will deform more than others. It can happen that certain elements deform to such an extent, that the solution quality degrades. To counteract this, one can introduce Jacobian based stiffening which is extensively used in [6]. This method stiffens the smaller elements w.r.t. the larger elements and in this way it reduces the deformation of the smaller elements. The elemental calculations performed to solve the equations of linear elasticity are multiplied by a term:

Z Ω[...]dQf = Z Ω[...]  J0 J1  β_dQ_f _(4.37)

Where J1_{is the Jacobian of a mapping which maps a standard wedge element onto the}

wedge element in question. J0_{is a scaling parameter. The parameter β determines the}

”stiffness” of the smaller elements. Since the wedge element is a triangular element extruded in time, we can take the mapping of the spatial triangular element and obtain the same ratio. J1 _{then becomes a ratio of areas and we rewrite the fraction}

term: _Z Ω[...]  J0 J1  β_{dΩ =}Z Ω[...]  A0 A1  β_dΩ _(4.38)

where A0 _{is a scaling area and A}1 _{is the area of the element in question. Depending}

on the quality of the spatial mesh, the Jacobian based stiffening technique shows to increase the quality of the map χ to a certain extent. The quality increase on a rect-angular domain remains limited, however, if the deformation is not to ’exotic’. Since the problem setup in this thesis didn’t include extreme displacements and difficult geometries, Jacobian based stiffening was not included in the FSI-problem. The dif-ficult geometries do appear in a complete blood-vessel model however, and therefore Jacobian based stiffening has been tested. Appendix C shows the result of some test calculations. When the deformation and the mesh density would vary more, the effect would have been larger as can be seen in [6].

4.5 Solvers

For the structure and fluid-domain equations the finite-element formulations result in a linear algebraic system Au = b. This system has been solved with a direct LU decomposition solver. The tolerance was set to 1 · 10−4_{, unless stated otherwise. For}

the fluid equations the resulting system is a non-linear set of equations P(u) = b. Hence Newton’s method has been employed. The linearized set of equation becomes P(u0) + ∂P(u_∂u0)u = b + ∂P(u_u0)u0. This system has been solved using a GMRES

iterative solver. The tolerance was set to 1 · 10−4_{. The number of Newton iteration}

(27)

Chapter 5

Fluid-structure coupling

The fluid and structure are coupled through the kinematic and dynamic conditions. These coupling conditions ensure continuity and conservation of momentum over the interface. The fluid and structure equations, together with the kinematic and dy-namic conditions, result in a coupled system which can be solved. In general, we can distinguish three methods for solving this coupled system. The weak- and strong-coupled solution methods solve the system in a staggered approach. The monolithic solution method solves the coupled system in a direct approach.

5.1 Kinematic and dynamic conditions

The kinematic and dynamic conditions ensure continuity and conservation of momen-tum over the interface. The first kinematic condition enforces the velocity of the fluid to be equal to the velocity of the solid at the interface (see section 3.24). This con-dition is satisfied by inserting the velocity of the structure as a boundary concon-dition in the fluid equations. The second kinematic condition enforces the displacement of the fluid domain to be equal to the displacement of the structure. This condition is satisfied by inserting the displacement of the structure as a boundary condition in the fluid-domain equations. The dynamic condition is satisfied by insertion of the fluid traction at the interface in one of the right-hand side terms of the structure equation.

5.2 Solution techniques

In general one can distinguish three solution techniques concerning fluid-structure interaction. The loosely- and strongly-coupled partitioned solution methods solve the system in a staggered approach. That is, for a specific time step the fluid equations are solved first, after which the structure equations are solved. These equations are often referred to as the segregated equations. The third solution technique is the so-called monolithic solution method and solves the fluid and structure equations simultaneously, i.e. the aggregated equations.

(28)

5.2.1 Loosely-coupled solution methods

Loosely-coupled solution methods solve the segregated equations once every time step and have the severe drawback that they only conserve energy in an asymptotic sense, i.e. when the mesh width and time step size reduces to zero. Although these methods have been used extensively, this disadvantage limits its application to small deforma-tion and small time steps. As discussed by Michler (see [1]), the energy discrepancy between the fluid and the structure reduce stability of the method. For large defor-mation and large time steps the method is inherently unstable. The advantage of this method is that it requires the least computation power of the three techniques. Another advantage of the staggered approach is the flexibility to use different solvers for fluid and structure. This advantage can play an important factor when one em-ploys different software packages to solve the fluid and the structure equations. To preserve this last advantage and eliminate the disadvantage of energy discrepancy, strongly-coupled solution methods have been introduced.

5.2.2 Strongly-coupled solution methods

Strongly-coupled solution methods solve the same segregated equations as the loosely-coupled solution methods, but solve these equations repeatively. This iterative pro-cedure continuous within each time step, until a certain criteria is satisfied. Due to the iterative procedure the coupling is made fully implicit and the energy discrepancy between the structure and fluid vanishes up to a certain threshold. Since this method has been employed in this thesis we further elaborate the procedure.

To start the iterative procedure we need an initial deflection at the interface uΓ. With

this initial guess uΓ the iteration procedure takes the following steps:

(0) Solve fluid domain: ud∈ Ud: D(ud, wd) = 0 ∀ wd∈ Wd

(1) Solve fluid: uf ∈ Uf, p ∈ P : F(uf, p, wf, q) = f(wf) ∀ wf ∈ Wf, p ∈ P

(2) Solve structure: us∈ Us: S(us, vs, ηs, ws) = s(ws, ηs) ∀ ws∈ Ws, ηs∈ Ns

Since the dynamic and kinematic conditions are enforced strongly they are not pre-sented as separate steps, but are included through boundary conditions and a right-hand side term.

Convergence of the iterative procedure is not guaranteed and if convergence is achieved the number of iterations needed to obtain a converged solution can be very large, es-pecially when the initial guess is not close enough to the actual solution. To accelerate the convergence process we introduce relaxation. This technique ”corrects” the solu-tion at the interface of the previous iterasolu-tion step in the following manner:

¯un+1

Γ = ω˜un+1Γ + (1 − ω)unΓ

where ¯un+1

Γ is the initial deflection for the next iteration step and ˜un+1Γ is the deflection

calculated from the previous iteration step. The parameter ω determines the amount of relaxation. This relaxation parameter can be taken constant over the simulation or can depend on the time step or iteration step. In this thesis ω is taken constant and is hand-fitted.

(29)

5.2.3 Monolithic solution methods

Monolithic solution methods solve the aggregated equations which are in general non-linear. Hence the method typically needs a Newton process to solve the system. This will result in a large linear system of equations and, moreover, it needs derivatives of the functionals D(ud, wd), F(uf, p, wf, q) and S(ws, us). These so-called shape

derivatives give rise to computational costly procedures. Due to different length and time scales between the fluid and the structure the Jacobian matrix may be severely ill-conditioned. Since the functionals must be differentiated w.r.t. the struc-ture unknowns at the interface, one is restricted to simultaneous calculation of these derivatives. Hence the employment of different solvers for the structure and fluid equations is impossible. This is a large disadvantage and it often doesn’t weigh up to the advantage of inherent energy conservation.

(30)

Chapter 6

Convergence analysis

To verify the numerical fluid and structure models, we analyze the convergence be-haviour of the error. To identify the error made by the discretization of the equations, we need an exact reference solution. For the structure and undeformed fluid domain problems typical test cases are available for which we know the exact solution. For the deforming fluid domain problem, however, we do not have an analytical solution. We therefore took the solution calculated on a very fine mesh and used it as the reference solution. Though the reference solution in this case includes an error itself, the convergence behaviour should be the same as if we would have an exact analytical reference solution.

6.1 Error definitions

To measure the error in the displacement and velocity field of the structure we define two norms, the L2-norm and the energy norm. The L2-norm is defined as:

||e||L2(Qs)= Z Qs |e|2_dQ s 1 2 (6.1) with the error defined as:

e = ur s− ucs

where ur

s is the displacement field of the structure calculated on the reference grid

(finest grid) or is the analytical solution and uc

sis the displacement field of the

struc-ture on a coarse grid. A natural way to define the energy norm is the total energy norm: ||e||E(Qs)= Z Qs ε(w) : σ(∇u) dQs+ Z Qs ρ|∂u_∂t|2_dQ s 1 2 (6.2) where the first term on the right-hand side inside the square root equals the strain-energy norm squared and the second term the kinetic-strain-energy norm squared. These norms are defined on the space-time domain. In the same way we define the norms on the spatial domain at the end of a space-time slab n:

||e||E(Ωs)= Z Ωs ε(w) : σ(∇u) dΩs+ Z Ωs ρ|∂u_∂t|2_dΩ s 1 2 (6.3)

(31)

||e||L2(Ωs)= Z Ωs |e|2_dΩ s 1 2 (6.4) For the fluid we define the L2-norms in the same way as for the structure:

||e||L2(Ωf)= Z Ωf |e|2_dΩ f 1 2 (6.5) ||e||L2(Qf)= Z Qf |e|2_dQ f 1 2 (6.6) Since the function space in which we seek a solution is H1 _{in space and L}₂ _{in time,}

we define the norm of the derivatives separately: ||e||T (Qf)= Z Qf ||∂e_∂t||2_dQ f 1 2 ||e||H1_(Q_f₎= Z Qf ||∇e||2_dQ f 1 2 (6.7) where ∇ indicates the gradient operator.

The integration is performed on element level. For all calculations we have employed second order integration accuracy to eliminate integration errors. For extruded un-structured triangular meshes, the difference between fine and course grids also intro-duce an integration error. We therefore employed structured extruded quadrilaterals in the convergence analysis.

6.2 Convergence study

In this section the convergence behaviour of the error w.r.t the mesh width is pre-sented. This is also known as a h-convergence study. For spatial convergence we expect the error to depend on the mesh width in the following manner:

||e||L2 ≤ Chk+1 (6.8)

||e||E ≤ Chk (6.9)

where h is the mesh width and k the polynomial order. Since we employ linear interpolation functions, the polynomial order is one and we expect the error in the L2-norm to converge quadratically and the error in the energy norm to converge

linearly. This convergence behaviour can be influenced by the stabilization however, as discussed by Dettmer and Peri´c (see [14]).

6.2.1 Structure model convergence

To evaluate the convergence rate of the structure model w.r.t the mesh width, we introduce the unit square problem. The unit square problem is very convenient for convergence analysis since we know the exact solution. The problem is defined on a unit square with us= 0 on the boundary. The body force is chosen to be:

f =

(λ + µ)((1 − 2y)(1 − 2x)

ρ0ω2yx(1 − y)(1 − x) − 2µy(1 − y) − 2(λ + 2µ)x(1 − x)

sin(ωt)

(32)

This will result in the exact displacement: us= 0 −xy(1 − x)(1 − y) sin(ωt)

With the displacement field known, we can calculate the exact velocity field and the stress tensor S(us). Note that the displacement field is independent of the stiffness

modulus and Poisson’s ratio. For a very small time step (∆t=0.005) the error in the L2-norm and strain-energy norm has been calculated for various mesh widths.

The error is calculated at an arbitrary time step t=1.25s. The density and stiffness modulus were set to 1.0 and Poisson’s ratio to 0.3 to validate the results with those obtained by Li and Wiberg [16]. Figures 6.1, 6.2 and 6.3 show the order of convergence. The dotted lines plotted in the figures give exact first and second order convergence for reference. We observe that the error in the L2-norm converges quadratically and

in the strain-energy norm linearly. These results are confirmed by the results obtained by Li and Wiberg. To analyze the temporal convergence, we calculate the error for different time-step sizes on a very fine mesh (h=0.005). Figure 6.4 shows the error convergence in the total-energy norm. We observe a third order convergence of the error. 101 10−4 10−3 h−1 ||e|| L2 ( Ωs ) Velocity field max ||e|| L 2(Ωs) O(h) O(h2)

Figure 6.1: Spatial error convergence in L2-norm of the displacement field for the unit

square problem. Second order convergence is observed (∆t = 0.005).

After the convergence behaviour has been verified with the unit square problem, the clamped beam problem has been analyzed. Both ends of the beam were fixed, the body force was set to zero and the initial displacement field attained a maximum at around 0.3L. The initial displacement field is obtained by calculating the steady response on a very fine grid to an overpressure on the lower side, see figure 6.5 for details of the setup. This displacement field is employed as initial condition for calcu-lations with different mesh densities. Figure 6.6 shows the displacement in y-direction at the centre of the beam against time for three different mesh widths. We observe that the displacement experiences a small phase shift when the mesh density changes. This behaviour can be expected, since the eigenfrequency of the discretized model depends on the stiffness matrix. The stiffness matrix changes when the mesh density

(33)

101 10−5 10−4 ||e|| L2 ( Ωs ) h−1 Displacement field max ||e|| L 2(Ωs) O(h) O(h2)

Figure 6.2: Spatial error convergence in L2-norm of the velocity field for the unit

square problem. Second order convergence is observed (∆t = 0.005).

101 10−3 10−2 ||e|| E( Ωs ) h−1 Strain−energy norm max ||e|| E(Ω s) O(h) O(h2)

Figure 6.3: Spatial error convergence in strain-energy norm for the unit square prob-lem. First order convergence is observed (∆t = 0.005).

is changed. This phase shift will cause the error of the solution for different mesh widths w.r.t. the fine grid solution to grow to almost the same value, namely when the displacement field calculated on the fine grid is in anti-phase. This is illustrated by figure 6.7. Hence the error convergence for this setup is somewhat difficult to de-termine and it will not be further analyzed, since the model has already been verified with the unit square problem.

(34)

101 10−5 10−4 10−3 ||e|| E( Ωs ) dt−1

Temporal convergence error in total energy norm

max ||e|| E(Ω s) O(∆t3) O(∆t2) O(∆t)

Figure 6.4: Temporal error convergence in total-energy norm for the unit square problem. Third order convergence is observed (h=0.005).

3L

L/4 0.3L

Figure 6.5: Clamped beam problem set up

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 time (s) Displacement in y−direction coarse mesh (h=0.05) fine mesh (h=0.02) finest mesh (h=0.01)

Figure 6.6: Displacement in y-direction in centre of beam versus time

6.2.2 Fluid model convergence

To investigate the influence of the mesh width and time step size on the solution of the discretized fluid model, calculations have been performed on both non-deforming and

(35)

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 time ||e|| L2 (Q s )

Phase shift in L₂ errors for Solid testcase I coarse1

coarse2

Figure 6.7: Typical behaviour of the error in L2-norm of the displacement field versus

time for two different mesh widths. Mesh density coarse1 < coarse2. The phase shift of the coarser mesh is larger than that of the finer mesh.

deforming domains. For the deforming domain we calculate the error w.r.t. a very fine mesh. Although the convergence behaviour of the stabilized fluid equations is compli-cated, Dettmer and Periˆc (see [14]) give an indication by applying a Fourier analysis on the convection-diffusion equation. They show that the accuracy of the stabilized convection-diffusion equation depends on ∆x, ∆t and τ, including mixed order terms. For linear time-discontinuous elements the lowest orders of the error in the amplitude are O(∆x2_{, ∆t}3_{, τ}2_{, ∆x}2_{τ, ∆tτ). Observe that depending on the time step size and}

mesh width, but also on τ, the order of convergence is determined. Hence, the choice of definition of the stabilization parameters influences the convergence behaviour. Undeformed fluid domain

To study the convergence behaviour of the fluid equations without deformation and especially the h-dependence of the error, we consider Poiseuille flow in a channel. The analytical solution is known and states a quadratic velocity profile: u = A(y − y2₎

and v = 0, where the amplitude A depends on the linear pressure gradient and the viscosity: A = 1

2µ(dxdp). The length of the channel was taken 3L and the diameter L.

The errors are defined as discussed in section 6.1, where the fine-grid solution is re-placed by the analytical solution. Quadratic boundary conditions for the velocity are enforced on the left- and right side of the domain. No-slip boundary conditions are enforced on the bottom and top of the domain. The initial condition was set to the quadratic velocity profile. The following parameter set has been employed: L=1.0, A=4.0 and µ = 0.0357, so Re = 28. The error in the L2-norm of the velocity field

cal-culated on the spatial domain is given in figure 6.8. For a small time step size (0.005) we see the error in the L2-norm of the error in the velocity field converging with a

or-der slightly more than two. The L2-norm of the velocity gradient however (see figure

6.9), converges exactly with O(h). The theoretical convergence behaviour discussed in section 6.2 holds quite well. The stabilization parameters do not contribute much to the convergence behaviour for the steady problem. This can be explained by the fact that the time step was chosen very small. This will cause the stabilization

(36)

pa-rameters to depend only on the time step size and hence the stabilization papa-rameters are constant for all mesh widths.

101 10−4 10−3 ||e|| L2 ( Ωf ) h−1 Spatial error convergence L

2−norm velocity ||e|| L 2(Ωf) O(h) O(h2)

Figure 6.8: Spatial error convergence of velocity in L2-norm. The order of convergence

is even somewhat larger than two.

101 10−2 10−1 ||e|| H 1 ( Ωf ) h−1 Spatial error convergence L

2−norm velocity gradient

||e||

H1

(Ωf)

O(h)

O(h2)

Figure 6.9: Spatial error convergence of velocity gradient in L2-norm. The order of

(37)

Deformed fluid domain

The deformation of the fluid domain is modeled through a ”pseudo solid” as discussed in section 4.4. Hence the error in the fluid solution on a moving domain does not only depend on the discretization of the flow field but also on the discretization of its domain displacement field. When analyzing the convergence rates, we must take into account that both errors can influence the convergence rate. For the convergence analysis of the structure we observed quadratic convergence in the L2-norm for spatial

convergence. Since we employ the same model for the ”pseudo solid”, only without the inertia term, we expect second order convergence for the domain displacement field. The convergence rates of the error in the fluid solution on a moving domain show to be far more complex. We must note that the errors for the fluid solution as defined in section 6.1, have been calculated on the computational domain. The test case setup is similar to the undeformed domain case, except the boundary on the lower side of the domain has a forced motion f(x, t) = 0.1L(1

2cos2πxL −12) sin ωt,

where ω was taken 2π. For details of the setup see figure 6.10. On the left hand side quadratic inflow boundary conditions are enforced, on the right hand side the traction is enforced. This traction is derived from a quadratic outflow velocity profile and zero pressure. The parameter set was taken equal to the undeformed fluid domain test case.

0.1L

L L L

L u=f(y)_v=0 τ

n=g

Figure 6.10: Deformed fluid domain test case

The errors on the spatial domain have been calculated at the end of each time step. Figure 6.11 shows an example of how the error in the L2-norm of the velocity behaves

in time. We observe an oscillating behaviour in the value of the error. This is to be expected, since the lower side of the domain is deformed by an oscillating function. The lower boundary is first deflated after which it is inflated. One can observe that two peaks appear during one oscillation, of which the first is somewhat larger than the second. Although this behaviour can be expected, since the deflation causes the elements to ”grow” which yields a worse approximation of the flow field, more research is needed to draw conclusions on this behaviour.

To study the error convergence w.r.t the time step size the maximum errors are evaluated, that is the error at the peak of the deflation. Figures 6.12-6.15 show the temporal convergence in the different norms. We observe the convergence of the error of the velocity field tends from a linear rate at larger time steps to even something more than quadratic rate at smaller time steps. The error of the velocity gradient shows similar behaviour. The error in the pressure field tends to go to linear convergence for smaller time steps. The error in the temporal derivative of the velocity field shows a behaviour which is one order lower than the velocity itself.

(38)

0 1 2 3 4 5 6 0 0.5 1 1.5 2 2.5 3 3.5x 10 −4 L

2 error of velocity in time

||e||

L2

Ωf

time (s)

Figure 6.11: Typical behaviour of error in L2-norm of velocity field versus time.

102 10−2 dt−1 ||e|| L2 ( Ωf ) velocity field max ||e|| L2(Ωf) O(∆t) O(∆t2)

Figure 6.12: Temporal convergence of error in L2-norm of velocity field. The

conver-gence changes from first order to more than second order for small time steps. The same analysis has been performed for spatial convergence. Figures 6.16-6.19 show the results. The error of the velocity tends toward quadratic convergence for smaller mesh width, and the error of the velocity gradient displays similar behaviour only one order lower. The convergence of the temporal derivative of the velocity shows the same behaviour as the velocity. The convergence rate of the error of the pressure is in between first and second order.

Although the convergence behaviour shows to be varying with mesh width and time step size, we see that almost all convergence orders are between first and second order. The observations indicate that the errors converge in a acceptable way, although more research is needed to draw detailed conclusions.

(39)

102 100 ||e|| L2 ( Ωf ) dt−1 pressure field max ||e|| L 2(Ωf) O(∆t)

Figure 6.13: Temporal convergence of error in L2-norm of pressure field. For smaller

time steps the order of convergence is almost one.

102 10−2 10−1 ||e|| H 1 ( Ωf ) dt−1 gradient velocity field

max ||e||_H1₍_Ω f)

O(∆t)

O(∆t2₎

Figure 6.14: Temporal convergence of error in L2-norm of velocity gradient. The

convergence changes from first order for large time steps to second for smaller time steps

(40)

102 10−2 10−1 ||e|| T( Ωf ) dt−1 temporal derivative velocity field

max ||e||

T(Ω

f)

O(∆t)

O(∆t2₎

Figure 6.15: Temporal convergence of error in L2-norm of temporal derivative of

velocity. The convergence behaviour is one order lower than that of the velocity field.

101 10−2 10−1 h−1 ||e|| L2 ( Ωf ) Velocity field max ||e|| L 2(Ωf) O(h) O(h2)

Figure 6.16: Spatial convergence of error in L2-norm of velocity field. The convergence

(41)

101 100 101 pressure field ||e|| L2 ( Ωf ) h−1 max ||e|| L 2(Ωf) O(h) O(h2₎

Figure 6.17: Spatial convergence of error in L2-norm of pressure field. The

conver-gence changes from first order for large h to second order for small h.

101 10−1

100

gradient velocity field

||e|| H 1 ( Ωf ) h−1 max ||e|| H1 (Ωf) O(h) O(h2)

Figure 6.18: Spatial convergence of error in L2-norm of velocity gradient. We observe

(42)

101 100

101

temporal derivative velocity field

||e|| T( Ωf ) h−1 max ||e|| T(Ω f) O(h) O(h2)

Figure 6.19: Spatial convergence of error in L2-norm of temporal derivative of velocity.

(43)

Chapter 7

Fluid-structure-interaction

simulation

In this chapter the numerical results of the coupled fluid-structure model will be pre-sented. The coupled fluid-structure system is solved using subiteration. To advance the subiteration procedure relaxation has been employed. We will first discuss the convergence behaviour of the subiteration procedure, after which two different test cases are analyzed.

7.1 Model setup

The coupled system has been investigated with two sets of boundary conditions. We start to study the behaviour of the coupled system with the enforcement of a steady quadratic inflow condition. The initial condition for the fluid is chosen to be the same quadratic velocity profile as the inflow condition. At the outflow, the pressure is enforced to be zero. The structure is modeled by a beam with boundary conditions as displayed in figure 7.1.

To simulate the blood flow pulses originating from the heart we also consider unsteady flow boundary conditions. The inflow boundary condition for the fluid is chosen such, that it models a pulsating flow field. The inflow velocity varies in time according to:

uf =A₂ ¯u 0 (1 − cos(ωt))

where ¯u equals the quadratic inflow profile y2_{− y. The outflow boundary condition}

is a zero pressure condition. The blood vessel is modeled by a beam with boundary conditions displayed in figure 7.1.

We define a reference parameter set given in table 7.1. We discretized the upper and lower structure with 80 linear quadrilaterals and the fluid domain with 200 linear quadrilaterals. For the computations with a domain length of 10, the number of elements were doubled.

(44)

p=0 h

L u(t)

xL

Figure 7.1: FSI-model setup

ρs E νs µf A x h L ∆t

85 1.0 · 106 _{0.49 0.035 4.0 5 0.25 1.0 0.02}

Table 7.1: Reference parameter set

α M Re

1.176·104 _1/85 ₂₈

Table 7.2: Characteristic reference parameters

The subiteration procedure solves the segregated equations in an iterative approach (see subsection 5.2.2). The iterative procedure continues within a time step until a certain stopping criteria is satisfied. We adopt the following stopping criteria for every time step:

maxei u, eiv ≤ 10−4 _(7.1) where ei u= ||ui s− ui+1s ||ΓL2 ||u1 s− u0s||ΓL2 (7.2) ei v= ||vi s− vi+1s ||ΓL2 ||v1 s− v0s||ΓL2 (7.3) where i denotes the iteration number within a time step. Since the velocity and displacement field are approximated by separate interpolation functions, we force both fields to reduce the error at the interface in the L2-norm w.r.t the first iteration

step by a factor 104_{. The maximum number of iterations per time step was taken 100}

for most of the calculations.

7.2 Results

7.2.1 Subiteration convergence

To study the convergence of the subiteration method we enforce the steady quadratic inflow condition and employ the reference parameter set. The reference parameter set

A Space-Time Finite Element Approach to the Numerical Simulation of Vascular Fluid-Solid Interaction

000000000

000000000

000000000

000000000

000000000

111111111

111111111

111111111

111111111

111111111

A Space-Time Finite Element Approach to

the Numerical Simulation of Vascular

Fluid-Solid Interaction

TUD/LR/EM

Kluyverweg 1, 2629HS Delft, The Netherlands

P.O. Box 5058, 2600 GB Delft, The Netherlands

Phone: +31 15 278 5460 Fax: +31 15 261 1465

ISSN 1574-6992

Summary

Contents

Chapter 1

Introduction

Chapter 2

Hemodynamics

2.1 Blood vessel model

2.2 Blood flow model

Chapter 3

Governing equations

3.1 Kinematic description of a continuum

3.2 Structure equations

3.3 Fluid equations

3.4 Interface conditions

3.4.1 Kinematic conditions

3.4.2 Dynamic condition

3.5 Non-dimensionalized equations

Chapter 4

Finite-element formulations

4.1 Space-time discontinuous Galerkin

discretiza-tion

4.2 Structure finite-element formulation

4.3 Fluid finite-element formulation

4.3.1 Mixed formulation

4.3.2 Stabilized formulation

4.3.3 Stabilization parameters

4.4 Fluid domain deformation

4.4.1 Domain motion

4.4.2 Jacobian-based stiffening

4.5 Solvers

Chapter 5

Fluid-structure coupling

5.1 Kinematic and dynamic conditions

5.2 Solution techniques

5.2.1 Loosely-coupled solution methods

5.2.2 Strongly-coupled solution methods

5.2.3 Monolithic solution methods

Chapter 6

Convergence analysis

6.1 Error definitions

6.2 Convergence study

6.2.1 Structure model convergence

6.2.2 Fluid model convergence

Chapter 7

Fluid-structure-interaction

simulation

7.1 Model setup

7.2 Results

7.2.1 Subiteration convergence