Adaptive Finite Element Methods for Fluid Structure Interaction Problems with Applications to Human Phonation

Interaction Problems with Applications to Human Phonation

NİYAZİ CEM DEĞİRMENCİ

Doctoral Thesis

Stockholm, Sweden 2018

TRITA EECS-AVL-2018:38 ISBN 978-91-7729-764-2

SE-100 44 Stockholm SWEDEN Akademisk avhandling som med tillstånd av Kungl Tekniska högskolan framlägges till offentlig granskning för avläggande av teknologie doktorsexamen i datalogi ons- dagen den 23 mai 2018 klockan 10.30 i Sing-Sing, våningsplan 2, Kungl Tekniska högskolan, Lindstedtsvägen 26, Stockholm.

© NİYAZİ CEM DEĞİRMENCİ, May 2018

Tryck: Universitetsservice US AB

Abstract

This work presents a unified framework for numerical solution of Fluid Structure Interaction (FSI) and acoustics problems with focus on human phonation.

The Finite Element Method is employed for numerical investigation of par- tial differential equations that model conservation of momentum and mass.

Since the resulting system of equations is very large, an efficient open source

high performance implementation is constructed and provided. In order to

gain accuracy for the numerical solutions, an adaptive mesh refinement strat-

egy is employed which reduces the computational cost in comparison to a

uniform refinement. Adaptive refinement of the mesh relies on computable

error indicators which appear as a combination of a computable residual and

the solution of a so-called dual problem acting as weights on computed resid-

uals. The first main achievement of this thesis is to apply this strategy to

numerical simulations of a benchmark problem for FSI. This FSI model is

further extended for contact handling and applied to a realistic vocal folds

geometry where the glottic wave formation was captured in the numerical

simulations. This is the second achievement in the presented work. The FSI

model is further coupled to an acoustics model through an acoustic analogy,

for vocal folds with flow induced oscillations for a domain constructed to cre-

ate the vowel /i/. The comparisons of the obtained pressure signal at specified

points with respect to results from literature for the same vowel is reported,

which is the final main result presented.

Sammanfattning

Detta arbete presenterar en enhetlig ram för numerisk lösning av fluid- strukturinteraktion (FSI) och akustikproblem med fokus på det mänskliga talet. En finita elementmetod används för numerisk lösning av de partiella dif- ferentialekvationer som beskriver konserveringslagar för moment och massa.

Eftersom det resulterande systemet av ekvationer är mycket stort, konstrueras en öppen källkod med hög prestanda. För att få hög noggrannhet i de nu- meriska lösningarna används en adaptiv nätförfiningsstrategi vilken minskar beräkningskostnaden jämfört med en uniform förfining.

Adaptiv förfining av nätet bygger på beräknade felindikatorer som byg- ger på en kombination av en beräkningsbar residual och lösningen av ett så kallat dualt problem. Den första huvudresultatet av denna avhandling är att utveckla en och validera denna strategi för en FSI-modell i ett benchmark problem.

Denna FSI-modell utvidgas vidare för att hantera kontaktmekanik, och används sedan för en realistisk modell av stämbandsstrukturerna där den glottiska vågformationen fångas i de numeriska simuleringarna. Detta är det andra huvudresultatet i det presenterade arbetet.

FSI-modellen kopplas också till en akustikmodell genom en akustisk analogi,

för modell konstruerad för att skapa vokalen / i /. Den erhållna trycksignalen

i ett antal punkter jämförs med resultat från litteraturen, vilket är det slutliga

huvudresultatet som presenteras.

To Nurettin, Süreyya and Ceren Değirmenci.

I would like to thank my supervisors Johan Jansson and Johan Hoffman for providing a liberal research atmosphere and fruitful discussions.

Doghonay Arjmand, Jeannette Hiromi Spühler, Love Lindholm, Henrik Holst and Kaspar Müller thank you for the inspiration at KTH.

Laura Saavedra, Margarida Moragues Ginard, Massimiliano Leoni thank you for the inspiration at BCAM.

Berat Beran Çamak desteğin için teşekkürler. Melih Temel, Arman Toplu,

Türker Küçük ve Emre Atsan iyi ki varsınız.

This thesis consists of an introduction and six papers. The included papers are:

Paper I

Johan Hoffman, Johan Jansson, Rodrigo Vilela De Abreu, Niyazi Cem Degir- menci, Niclas Jansson, Kaspar Müller, Murtazo Nazarov, Jeannette Hiromi Spühler,

"Unicorn: Parallel adaptive finite element simulation of turbulent flow and fluid- structure interaction for deforming domains and complex geometry", Computers &

Fluids, Vol. 80, pp. 310–419, 2013.

The author of this thesis developed the topological mesh adaptation implementa- tion.

Paper II

Johan Jansson, Niyazi Cem Degirmenci, Johan Hoffman, "Framework for adaptive fluid-structure interaction with industrial applications", Int. J. Materials Engineer- ing Innovation, Vol. 4, No. 2, 2013.

The author of this thesis has done the implementation and performed the numerical simulations, and conributed to the writing of the manuscript.

Paper III

Johan Jansson, Niyazi Cem Degirmenci, Johan Hoffman, "Adaptive Unified Con- tinuum FEM Modeling of a 3D FSI Benchmark Problem". Int. J. Numer. Meth.

Biomed. Engng., doi: 10.1002/cnm.2851, 2016.

The author of this thesis has done the implementation of the primal solver, the integration of an existing dual solver into the adaptive algorithm, performed the numerical simulations and authored parts of the manuscript.

vii

Paper IV

Jeannette Hiromi Spühler, Niyazi Cem Degirmenci, Johan Jansson, Johan Hoffman,

"A 3D full-friction contact model for fluid-structure interaction problem", Submit- ted, 2017.

The first and second author performed the implementations and the numerical ex- periments in close cooperation. The modeling and authoring of the manuscript was done in collaboration between all authors.

Paper V

Johan Hoffman, Johan Jansson, Niyazi Cem Degirmenci, Jeannette Hiromi Spühler, Rodrigo Vilela De Abreu, Niclas Jansson, Aurélien Larcher, "FEniCS-HPC coupled multi-physics in computational fluid dynamics", High-Performance Scientific Com- puting: First JARA-HPC Symposium, JHPCS 2016, Aachen, Germany, October 4–5, 2016, Revised Selected Papers, pp. 58–69, 2017.

The author of this thesis has done the implementation and run the numerical sim- ulations for the setting of the vocal folds.

Paper VI

Niyazi Cem Degirmenci, Johan Jansson, Johan Hoffman, Marc Arnela, Patricia Sánchez-Martín, Oriol Guasch, Sten Ternström, "A unified numerical simulation of vowel production that comprises phonation and the emitted sound". Proceedings of INTERSPEECH 2017 (peer-reviewed), Stockholm, Sweden, pp. 3492–3496, Au- gust 20–24, 2017

The author of this thesis has done the implementation and run the numerical sim-

ulations as well as writing parts of the manuscript.

Contents ix

1 Introduction 1

1.1 Research Question . . . . 1

1.2 Background . . . . 2

1.2.1 Voice Production Physics . . . . 2

1.2.2 The Mathematical Model . . . . 4

Mathematical Modeling of the Incompressible Continuum . . 4

Mathematical Modeling of the Acoustics . . . . 5

Mathematical Modeling of the Contact Detection . . . . 6

1.2.3 A First Introduction to Numerical Solution of PDE’s by the Finite Element Method . . . . 6

The Model Problem . . . . 6

Weak Formulation of The Model Problem . . . . 7

Ritz-Galerkin Approximation . . . . 7

1.2.4 A Second Introduction to Solving PDE’s with the Finite El- ement Method . . . . 8

Approximation of Weak Solutions . . . . 8

The Weak Derivative . . . . 8

Sobolev Spaces . . . . 9

The Galerkin Method . . . . 9

The Finite Element Method . . . . 10

1.2.5 Assessing the Scalability of the Implementations . . . . 10

Strong Scaling . . . . 11

Amdahl’s Law . . . . 11

Weak Scaling . . . . 11

Gustafson’s Law . . . . 11

1.3 State of the Art . . . . 12

1.3.1 Existence, uniqueness, regularity results for the solutions . . 12

1.3.2 Numerical results for the problem . . . . 13

2 Methods 15

ix

2.1 Fluid-Structure Interaction Coupling . . . . 15

2.2 Coordinate Systems . . . . 16

2.2.1 Eulerian description of the motion . . . . 17

2.2.2 Lagrangian description of the motion . . . . 17

2.2.3 Arbitrary Lagrangian Eulerian description of the motion . . . 17

Mesh Smoothing . . . . 18

2.3 Stabilized Galerkin Formulation . . . . 19

3 Contributions and Results 21 3.1 A Posteriori Error Estimation . . . . 21

3.2 Contact Modeling . . . . 23

3.3 Acoustic Coupling . . . . 26

3.4 Implementation Details . . . . 27

4 Discussion 29 4.1 Adaptive FEM for FSI problems . . . . 29

4.2 Contact Modeling . . . . 30

4.3 Acoustic Coupling . . . . 31

Bibliography 33

Introduction

1.1 Research Question

The focus of the presented work is to develop an efficient adaptive finite element method (FEM) for fluid-structure interaction (FSI) problems that can be applied to simulate human phonation. The main objectives can be listed as:

• Adaptive FEM for FSI problems:

The objective is to reduce the error in the quantity interest in an efficient way for the numerical approximation, where a particular challenge is the formulation of an estimate for a time dependent nonlinear problem and its implementation. Papers II and III focus on this objective and an introduction is provided in Section 3.1

• Contact Modeling:

The objective is to be able to handle problems when the domain boundaries come in close proximity, with challenges to keep the condition number of the resulting linear system low as well as preventing inverted elements. Paper IV focus on this objective and an introduction is provided in Section 3.2

• Acoustic Coupling:

The objective is to be able to use incompressible FSI-contact computation results as an input to a wave propagation solver for the simulation of human phonation. The challenges involve testing different possible acoustic analo- gies, testing different boundary conditions to prevent numerical artefacts in the form of wave reflections from the domain boundaries, implementation of solvers in the same framework, with particular challenges in mesh partitioning and data transfer for two solvers operating on two different domains and dis- cretizations. Paper VI focus on this objective and an introduction is provided n Section 3.3

1

• HPC implementation:

An implementation with good scalability properties is necessary to be able to approximate solution of 3D problems numerically with high accuracy. Papers I and V focus on this objective and some details of the open source finite element framework are provided in Section 3.4

1.2 Background

1.2.1 Voice Production Physics

Larynx is an organ which is considered to have a key role during the migration of the first vertebrates from an aquatic to a terrestrial environment nearly 350 million years ago in the early Carboniferous period of the Paleozoic era [52]. Similar to that of bichir Polypterus, which possesses both external gills and paired ventral lungs, the ancient "larynx" is thought to be composed of a simple muscle sphincter to protect the lung from the water [60].

The larynx has evolved and changed during time and besides airway protection, and for humans it began to fulfill the additional role of being an acoustic source for sound and speech production with the help of the vocal fold structures contained inside. Figure 1.1 depicts the organ and inner structures as well as the important muscles used for controlling it.

The physics of voice production is rather intricate. In a nutshell, the air em- anating from the lungs impinges on the vocal folds (VF) and pushes them apart until their elasticity takes over and the VF close again. These vibrations result in a pulsating glottal jet with maximum velocity ∼ 14 − 45ms

and a Reynolds number in the range O(10

) − O(10

) [47]. Humans can adduct the VF as well as manipulate the tension and change elastic properties actively by controlling the thyroarytenoid, posterior cricoarytenoid, cricothyroid, and lateral cricoarytenoid muscles. The pressure at the glottis decreases due to the jet flow in the vocal tract (VT) by Bernoulli’s law and a self sustained oscillation is established. See e.g., [59]

for a comprehensive introduction to the topic.

The flow dynamics and VF vibrations result in acoustic sources of mono polar, dipolar and quadrupolar character [66]. The Fundamental frequency is between 80 − 220 Hz. with lower frequencies for male voice due to longer length of the VF.

Acoustic waves are generated, then propagate through the VT and finally become radiated outwards. In the case of vowels, the VT resonances (formants) get excited, which actually allows one to distinguish one vowel sound from another.

The process is described as “chaotic” in [32], referring to a pseudo random

behavior where the precise output is unpredictable and extremely sensitive to slight

differences in initial conditions. The modeling and the numerical simulation of such

complex phenomena poses a big challenge.

(a) (b)

(c) (d)

Figure 1.1: The Larynx (a,b), and the muscles involved in adducting as well as

manipulating the tension (c,d) [54] (Images under public license in their country of

origin.)

1.2.2 The Mathematical Model

Mathematical Modeling of the Incompressible Continuum

In order to describe the motion of the incompressible continuum consisting of a fluid and solid structures, conservation of mass [49], conservation of momentum [48] and constitutive laws derived from the fundamental stress principle of Cauchy- Euler [62] are used as components to formulate a mathematical model, which is elucidated in [27].

The conservation laws can be represented as partial differential equations (PDE) of the form:

∂

∂t u

(x, t) + ∂

∂x

f

(u(x, t)) = s

. (1.1) where u is the vector valued function for the conserved quantity, with u

represent- ing its i

spatial component, f is the flux function and s is the source term. Here the Einstein convention is used where repeated indices signify summation.

The corresponding Cauchy problem for unbounded domains is

∂

∂t u

(x, t) + ∂

∂x

f

(u(x, t)) = s

, (x, t) ∈ R

× T, (1.2) u(x, t

) = u

(x), x ∈ R

,

where T = [t

, ∞), or an initial boundary value problem (IBVP) in a finite spatial domain Ω ∈ R

with a boundary Γ:

∂

∂t u

(x, t) + ∂

∂x

f

(u(x, t)) = s

, (x, t) ∈ Ω × T, (1.3) u(x, t

) = u

(x), x ∈ Ω,

αu(x, t) + β ∂u(x, t)

∂x

n

= u

(t) (x, t), ∈ Γ × T,

where n represents the outward normal for the boundary Γ. Appropriate initial and boundary conditions must be considered for the well-posedness of the resulting problem. The system of equations for the conservation of the momentum and the mass for incompressible Newtonian fluid flow with constant density and with Dirichlet boundary conditions are given as:

ρ(∂

u + (u · ∇)u) − µ∇

u + ∇p − s = 0, (x, t) ∈ Ω × T, (1.4)

∇ · u = 0, (x, t) ∈ Ω × T, u(·, t

) = u

, x ∈ Ω, u(x, t) = u

(x, t) (x, t), ∈ Γ × T,

where u is the velocity, p is the pressure, ρ is the density, µ dynamic viscosity and s

external forces, inside the spatial domain Ω with boundary Γ, and the time domain

T .

The constitutive equations for the solid structure in the context of fluid-structure interaction model corresponds to the relation between internal stresses and the de- formations of the structure. A diligent and step by step derivation of the con- stitutive equations can be found in [21]. We are in particular interested in the Neo-Hookean material model [61], which for the incompressible case reduces to:

σ

= −pI + C

F F

, (1.5)

where σ

is the Cauchy deviatoric, p is the pressure, C

is a material constant (shear modulus) and F is the deformation gradient.

Mathematical Modeling of the Acoustics

We refer to [19] for the modeling of the emanating sound which is done by observing these three principles during the motion of a sound wave:

• The gas moves and changes the density.

• The change in density corresponds to a change in pressure.

• Pressure inequalities generate gas motion.

The resulting equation for the sound assuming that the velocity is independent of the wavelength is given as:

p

= c

p

(1.6)

where c is the speed of the sound and p is the acoustic pressure. The problem can be posed for example as an IVP in 1D such as:

p

= c

p

−∞ < x < ∞, 0 < t < ∞ (1.7)

p(x, 0) = f (x) −∞ < x < ∞

p

(x, 0) = g(x) −∞ < x < ∞

or can also be posed for example as an IBVP for a finite domain Ω together with initial and boundary conditions as in:

p

= c

p

x ∈ Ω, 0 < t < ∞ (1.8)

p(x, 0) = f (x) x ∈ Ω

p

(x, 0) = g(x) x ∈ Ω

p(x, t) = h(x, t) x ∈ ∂Ω, 0 < t < ∞

In this work we are interested in the system of first order PDEs that describe sound propagation in a fluid moving with mean velocity ¯ U :

1

ρ

c

D

p + ∇.u = Q (1.9)

ρ

D

u + ∇p = f (1.10)

Where D

:= ∂

+ ¯ U .∇ is the material derivative, p acoustic pressure, u is acoustic velocity, Q a volume source distribution and f external body force per unit volume, ρ

mean air density and c

speed of the sound, together with appropriate boundary and initial conditions [23].

Mathematical Modeling of the Contact Detection

A problematic phenomenon to model is the contact, for which we defer address- ing the reason for the difficulty as well as the methodology followed to resolve it numerically to later sections.

Here however we present the Eikonal Equation, in order to model the distance from a point in the domain to the domain boundaries. In our application it is used to check the distance of the vocal fold boundaries from each other.

The Eikonal equation used in our work for the domain Ω is of the form:

|∇u(x)| = 1, x ∈ Ω (1.11)

u(x) = 0, x ∈ ∂Ω

together with its boundary conditions. The derivation of the Eikonal equation from the wave equation can be found in many textbooks including [8].

The meaning of the derivatives and solution spaces is the topic of the subsection 1.2.4. Another point that should be clarified is that for the equations (1.4) and (1.5) the constitutive laws for the structure are presented for a material point x, while for the fluid flow part x is a fixed spatial coordinate. The coordinate systems are discussed in the section 2.2

1.2.3 A First Introduction to Numerical Solution of PDE’s by the Finite Element Method

In this section we are going to examine a model PDE and apply the classical steps to obtain a numerical solution using the FEM without the intention of being complete. Some of the missing foundation for the ideas are explained in the second introduction part.

The Model Problem

Consider the boundary value problem in the interval I = [0, 1]

− ∂

u

∂x

= f x ∈ (0, 1) (1.12)

u(0) = u(1) = 0 (1.13)

which is in strong form and will be used for our demonstration purpose.

Weak Formulation of The Model Problem

If we multiply the strong form with any sufficiently smooth function v with v(0) = v(1) = 0 and perform integration by parts we get:

− Z

I

∂

u

∂x

v dx = Z

I

f v dx (1.14)

− ∂u

∂x v|

+ Z

I

∂u

∂x

∂v

∂x dx = Z

I

f v dx (1.15)

Z

I

∂u

∂x

∂v

∂x dx = Z

I

f v dx (1.16)

where the boundary terms disappear due to v(0) = v(1) = 0.

By using the space V = {v ∈ L

(0, 1) : v(0) = v(1) = 0, R

I

( ∂v

∂x )

dx < ∞} and the definitions:

a(u, v) :=

Z

I

∂u

∂x

∂v

∂x dx (1.17)

L(v) :=

Z

I

f v dx (1.18)

we can formulate the problem as find u ∈ V such that:

a(u, v) = L(v) ∀v ∈ V (1.19)

This is the weak formulation of the same problem and the equivalence of both formulations under regularity assumptions on f and u can be shown [6].

Ritz-Galerkin Approximation

The space V defined in the previous subsection has infinite dimensions. Consider a finite dimensional subspace S ⊂ V and the Ritz-Galerkin approximation of the weak solution in this space is formulated as follows:

Find u

∈ S s.t:

a(u

, v

) = L(v

) ∀v

∈ S (1.20)

Using the information that S has a finite basis {φ

: 1 ≤ i ≤ n}, we can construct a finite linear system of equations by testing 1.20 against each φ

u

:=

n

X

j=1

U

φ

(1.21)

A

:= a(φ

, φ

) (1.22)

b

:= L(φ

) (1.23)

and finally getting the system

AU = b (1.24)

It can be shown that this system has a unique solution [6].

In this work we use the Ritz-Galerkin method to obtain a numerical approxi- mation of the weak solutions. This method is also a very useful tool to show the existence of a weak solution for 1.19 by showing the existence of solutions for a series of problems on finite dimensional spaces converging to the original problem in infinite dimensions as shown in section 1.2.4.

1.2.4 A Second Introduction to Solving PDE’s with the Finite Element Method

Here we provide some more foundation to the ideas presented previously and apply them for the partial differential equations resulting from the conservation laws.

Approximation of Weak Solutions

Since finding solutions to PDEs in the classical sense is not always possible, a way to proceed is to search for a solution in the distributional sense, by testing the equations against smooth functions with compact support (C

) and thus relaxing the requirements for the point-wise regularity.

Assuming φ ∈ C

(R

× T ), the weak formulation of the problem (1.2) is stated as:

Z

Rⁿ×T

∂

∂t φ(x, t)u

(x, t) + ∂

∂x

φ(x, t)f

(u(x, t)) + φ(x, t)s

dxdt+ (1.25) Z

Rⁿ

(u

)φ(x, t

)dx = 0

which is constructed by multiplying the original PDE by φ and integrating in R

×T and then using integration by parts. A solution u is called a weak solution if it satisfies the equation (1.25) for every φ ∈ C

(R

× T )

However searching for solutions in special complete normed vector spaces (Ba- nach spaces) has the advantages of being able to use powerful functional analysis tools.

The Weak Derivative

Assume U ⊂ R

is an open and bounded domain, u, v ∈ L

(U ) are locally inte- grable functions on every compact subset of U , then v is called the α

weak partial derivative of u if:

Z

U

uD

φdU = (−1)

Z

U

vφdU

for all φ ∈ C

(U ). Here for the multi index α = (α

, ...α

) D

is defined as :

D

φ = ∂

∂x

...∂x

We also write v = D

u.

Sobolev Spaces

The Sobolev space W

(U ) is defined as

W

(U ) := {u : u ∈ L

(U ), D

u ∈ L

(U ) ∀|α| ≤ k}

||u||

:=

( ( P

|α|≤k

R

U

|D

u|

dU )

(1 ≤ p < ∞) P

|α|≤k

ess sup

|D

u| (p = ∞)

The Sobolev Spaces are complete [15] and the closure of C

(U ) in W

(U ) is denoted by W

(U ).

Now taking T

= (t

, t

) with t

< t

< ∞, the problem (1.4) can be stated as find u

, p ∈ W

(Ω × T

) for i = 1, .., n such that:

Z

Ω×T⁰ n

X

i=1

ρ(∂

u

v

+ u · ∇u

v

) + µ ∂u

∂x

∂v

∂x

+ ∂p

∂x

v

− s

v

+ ∂u

∂x

q dx dt = 0 (1.26) is satisfied ∀v

, q ∈ W

(Ω × T

).

The Galerkin Method

One way to construct a weak solution to the PDE is by finding solutions on a series of finite dimensional subspaces. When these subspaces converge to the original Sobolev space, it is possible to pass to limits for the finite dimensional solutions under certain conditions.

The summary of the method in [15] after adapting the notation to the problem (1.3) is as follows:

Let f be a first order partial differential operator in (1.3) and {w

}

be an orthogonal basis of W

(Ω) and orthonormal basis of L

(Ω). For each m, we seek a function u

: [t

, t

] → W

(Ω) of the form:

u

(t) :=

m

X

k=1

d

(t)w

(1.27)

where the coefficients d

(t)(t

≤ t

, k = 1, · · · , m) satisfy d

(0) =

Z

Ω

u

w

dx (1.28)

Z

Ω

u

w

− f

(u) ∂w

∂x

− s

w

dx = 0 (i = 1 · · · n) (1.29) .

If it is possible to show the existence of each u

and prove necessary energy estimates that depend on f and s; one can show existence and uniqueness of weak solutions for (1.3).

The Finite Element Method

Using piecewise polynomials as the basis {w

} is a natural option to build finite dimensional subspaces of the previously mentioned Sobolev spaces in a systematic manner.

The following definition is taken from [6] which is following Ciarlet’s definition of a finite element [9] .

Let:

• K ⊂ R

be a bounded closed set with nonempty interior and piecewise smooth boundary (the element domain)

• P

be a finite-dimensional space of functions on K (the space of shape func- tions)

• N

= {N

, N

· · · N

} be a basis for the dual space of P

Then (K, P

, N

) is called a finite element.

A set t

= {K

}, where ¯ U = ∪

K

, for a finite number of element domains K

with the requirement that i 6= j ⇒ ˚ K

∩ ˚ K

= ∅, is called a mesh. In the simple case of U ⊂ R

, the element domains may be triangles, and if U ⊂ R

the element domains may be tetrahedrons. If the mesh further satisfies the geometrical restriction that 2 neighbor elements share only one facet (no hanging nodes), then it is called a geometrically admissible mesh.

Let t

be a geometrically admissible mesh and V

= {v

: v

|K ∈ P

}. If for every K ∈ t

, P

⊂ W

(K) and V

⊂ C

( ¯ Ω) then V

∈ W

(Ω). If in addition x ∈ ∂Ω ⇒ v

(x) = 0, ∀v

∈ V

, then V

⊂ W

(Ω) [10].

1.2.5 Assessing the Scalability of the Implementations

An important property of the parallel algorithms is how well they scale when the

number of processes increase. In order to quantify this, mostly strong and weak

scaling concepts are used which are defined below. Paper I and V contains these

evaluations.

Strong Scaling

Strong scaling evaluates the algorithms by relating the total execution time to the number of processes for a constant problem size. The strong scaling of the ideal algorithm is the linear function S

(N ) = N . It uses the Speedup in Latency concept S(N ):

S(N ) := t

t

(1.30)

where N is the total number of processes, t

is the time needed for one process to accomplish the task, and t

is the time needed for N processes to accomplish the same task.

Amdahl’s Law

Amdahl has used 1.30 as an argument against parallel architectures in [2]. With 0 < α < 1 representing the proportion of the total time spent in the portion of the algorithm which is not parallelizable, he rewrites the speedup in latency as:

S(N ) := t

t

= t

αt

+ (1 − α)t

(1.31)

= t

αt

+ (1 − α)t

N

= 1

α + (1 − α) N

(1.32)

lim

N →∞

S(N ) = 1

α (1.33)

.

The Speedup in latency is therefore bounded by 1/α. The speedup graph for various α with respect to N is given in Figure 1.2.

Weak Scaling

A more optimistic perspective on the evaluation is proposed by Gustafson in [24]

where the author argues that the problem size will likely grow with the number of processes for many industrial applications instead of being constant as in the previous perspective.

Weak scaling evaluates the algorithms by relating the total execution time to the number of processes for a constant problem size per process. The weak scaling of the ideal algorithm is the constant function S

(N ) = 1.

Gustafson’s Law

Rather than using the Speedup in Latency, Gustafson uses the definition of scaled

speedup SS(N ) which is the time a single process would need to crunch the same

amount of data that N processes can process in unit time.

10⁰ 10¹ 10² 10³ 10⁴ 10⁵ N

0 2 4 6 8 10 12 14 16 18 20

S(N)

α=0.5 α=0.2 α=0.05

Figure 1.2: Speedup in latency versus number of processes for different levels of α.

SS(N ) := αt

+ N (1 − α)t

t

= N + (1 − N )α (1.34)

N →∞

lim SS(N ) = ∞ (1.35) .

Assuming unbounded problem sizes available, the scaled speedup is not bounded.

This perspective is applicable to most of the industrial applications where FEM is applied, as bigger solution subspaces are necessary to reduce the error as much as possible.

1.3 State of the Art

1.3.1 Existence, uniqueness, regularity results for the solutions The global existence, uniqueness and regularity of solutions of the incompressible Navier-Stokes equations for dimensions of three is an open problem [1].

Leray introduced a global weak solution [44] that is shown to be unique and regular for the case of dimension 2 [57, 41]. Strong solutions with more regularity have been shown either locally or when the initial data and force term is small, a review of important advances on this problem can be found in [20].

For the Fluid-Structure Interaction (FSI) problem an existence result of weak

solutions for a 2 dimensional fluid domain interacting with a 1 dimensional structure

without contact is shown in [5], as well as an overview of recent advances focusing

on existence, uniqueness and regularity results.

Existence of strong solutions for a similar setting is shown in [42, 43]. As stated in these studies, the contact phenomenon for the incompressible continuum is a source of difficulty for theoretical as well numerical studies, which is relevant for the mathematical modeling of vocal folds.

Existence, uniqueness and regularity results for wave equations can be found in [15] among other textbooks.

1.3.2 Numerical results for the problem

As for numerical test results for the FSI problem, a full-scale wind turbine applica- tion is presented in [4], which is one of the most advanced industrial FSI applications with turbulence as a key phenomenon.

Three dimensional numerical modeling of human vocal folds is studied in [39, 53], with prescribed vocal folds motion.

Mittal and co-workers in [65, 67] carry out Direct Numerical Simulation (DNS), with the Reynolds number artificially reduced by one order of magnitude, and with a partitioned FSI approach where the structure is represented as an immersed boundary in a fluid simulation. Contact is enforced by a kinematic constraint, which is not described in detail.

Jiang and co-workers in [34] simulate full coupling with acoustics for the problem

using a compressible continuum model as well as a multi layered structure model for

the self oscillating vocal folds. For these simulations however the Reynolds number

is again reduced and contact is prevented by always ensuring some gap between the

vocal folds by limiting the movement of the vocal folds.

Methods

2.1 Fluid-Structure Interaction Coupling

In order to model the coupling of the fluid and the structure, employing a weak coupling strategy [16, 30] is the first possibility. This strategy is characterized by independent discretization of the fluid and the solid domains together with an explicit time stepping scheme and exchanging information over the fluid-structure interface. This strategy makes it possible to reuse existing solvers focused on solving each problem separately. This however comes with the cost of introducing extra instability sources, and these methods are not suitable for the numerical solution of every problem.

In order to provide stability, implicit time stepping methods (also known as strong coupling) [56], or semi-implicit time stepping methods [18] using a combina- tion of implicit and explicit time stepping for different terms is possible.

Another direction is employing a monolithic strategy [63] where the fluid-structure continuum is discretized as one, allowing straightforward description of error con- trol techniques and providing stability at the expense of solving a bigger system of equations.

A monolithic description based on a fully Eulerian coordinate formulation of a fluid-structure model is described in [14]. Here the continuum model consists of the conservation equations with a general stress variable, with an additional structure displacement variable and an additional equation for interface capturing. This last equation is necessary to prevent interface smearing on a fixed computational mesh.

This way an initial position˝ set(IP) is explicitly tracked. An adaptive method is presented with 2D numerical test results.

The theoretical foundation for a goal-oriented adaptive method in an FSI-ALE setting on a fixed reference domain is presented in [50]. The method is verified and validated for the 2D FSI-I Hron and Turek benchmark in [64] where the deformation is very small, and also for a 3D test problem, again with very small deformations.

15

(a) (b)

(c) (d)

Figure 2.1: Two choices for the description of the motion of a deforming solid in a stationary fluid flow: initial configuration (a), a bent structure in an Eulerian coordinate system together with a mesh (b), a bent structure on a Lagrangian coordinate system together with a mesh (c), a bent structure on an ALE coordinate system with a a smoothed mesh (d)

2.2 Coordinate Systems

It is possible to present the same PDE in different coordinate systems assuming

existence of a one to one mapping between the two coordinate systems at each time

instant. In the discretized problem we let the mesh track the deformation of the

coordinate system so that the velocity of the coordinate system β is related to the

mesh velocity m.

2.2.1 Eulerian description of the motion

The Eulerian description of motion is defined with respect to a fixed coordinate system, which is often preferred in fluid mechanics.

The main problem with this setting in the case of describing motion of a solid body is that since the spatial resolution is finite, smearing of the interface of the body occurs. This is of course not a problem in the case of a homogeneous fluid flow problem, where the exact location of a certain particle is uninteresting. A problematic case in fluid-structure interaction however is illustrated in Figure 2.1b, where the motion of a bending bar is examined and the shape in bent state is often interesting.

2.2.2 Lagrangian description of the motion

The Lagrangian description presents the equations for a particle and the coordinate system follows the particles in motion. This setting is preferred in solid mechanics and often the convective terms –that appear in the Eulerian description of motion–

disappear leaving a simpler system of equations to solve. Since the mesh tracks the individual particles, the tracking of interfaces becomes trivial.

The main problem with this setting in the case of large deformations is that the quality of the mesh decreases, often leading to a bad condition number for the system of the resulting linear equations. This problem is illustrated in Figure 2.1c. In the extreme cases inverted elements may appear where the condition i 6= j ⇒ ˚ K

∩ ˚ K

= ∅ is not satisfied. In such cases an expensive re-meshing operation becomes necessary, and the projection of the variables from the previous time step introduces an additional source of error.

2.2.3 Arbitrary Lagrangian Eulerian description of the motion The Arbitrary Lagrangian Eulerian (ALE) description is designed to overcome the shortcomings of the previously described perspectives. The mesh can be moved with an arbitrary velocity such that some phenomenon is captured in some region (like the interface of a bending bar) and in other parts the velocity may be chosen to optimize the quality of the mesh. In the case of conservation laws, the ALE description adds convective terms to the partial differential equations with the velocity of the coordinate system [13]. In Figure 2.1d the coordinate velocity is equal to the structure velocity in the solid part to capture the interface, and in the fluid part the coordinate velocity is decided using a mesh smoothing algorithm.

With Ω ⊂ R

the spatial domain, and Q = Ω × I is a space-time domain with

I = (t

, t

] a time interval and ∂Ω = Γ

∪ Γ

, ˚ Γ

∩ ˚ Γ

= ∅, the FSI equations

for the unified continuum model with unknowns u, p, σ

, θ can be given in an ALE

coordinate system with β velocity as:

ρ(u

+ ((u − β) · ∇)u) − ∇ · σ − ρg = 0 (x, t) ∈ Q (2.1)

∇ · u = 0 (x, t) ∈ Q u(x, t

) = u

(x) x ∈ Ω

u(x, t) = u

(t) (x, t) ∈ Γ

× [t

, t

]

σ = ¯ σ − pI (x, t) ∈ Q (2.2)

¯

σ = θ¯ σ

+ (1 − θ)¯ σ

(x, t) ∈ Q

¯

σ

= µ

(∇u + ∇u

) (x, t) ∈ Q

p(x, t) = p

(t) (x, t) ∈ Γ

× [t

, t

]

¯

σ

= 2µ

 + ∇u¯ σ

+ ¯ σ

∇u

(x, t) ∈ Q (2.3) θ

+ ((u − β) · ∇)θ = 0 (x, t) ∈ Q

¯

σ

(x, t

) = σ

(x) x ∈ Ω

θ(x, t

) = θ

(x) x ∈ Ω

where µ

is the shear modulus (solid) and µ

is the dynamic viscosity (fluid) and θ marks the continuum either as fluid if θ = 1 or solid if θ = 0.

Mesh Smoothing

In order to capture the structure interface, β is defined to be equal to the continuum velocity u when θ = 0. For the fluid part, the velocity of the mesh is obtained by solving two additional equations, one Laplace type equation and one nonlinear equation to improve the mesh quality as explained in Paper IV.

Linear Smoother The Linear Smoother solves a linear elastic equation for the mesh velocity m where the vertices are diffusively relocated over the domain,

˙

m = ∇ · τ

(2.4)

τ

= 2µ

(m). (2.5)

This is a simple and fast method similar to the FSI-GST in [58]. We measure that

adjusting the mesh quality by advancing 1 time step with the Non-Linear Smoother

typically takes 10 times longer than 1 time step with the Linear Smoother. How-

ever, there is no certainty that the mesh quality is enhanced since no information

with respect to shape of the cells is taken into account.

Non-Linear Smoother In the Non-Linear Smoother the deformation of the mesh is formulated as a time-dependent non-linear elasticity problem where we apply the incompressible Neo-Hookean constitutive model to describe the elastic property of the mesh. The stiffness of the model is weighted by the quality Q(K) of the cell element K ∈ T

:

Q(K) := ||F ||

d · det(F )

, (2.6)

where d stands for the dimension of the mesh. F is the deformation gradient between a tetrahedral element K ∈ T

and an ideal element ˆ K with optimal quality and ||F ||

indicates its Frobenius norm.

To calculate the deformation velocity m of the elastic model, the following partial differential equations are solved:

F

= deformation gradient between K ∈ T

and ˆ K.

F ˙

= −F

∇m,

 = 1

2 (1 − FF

), strain tensor τ

= 2µ

 + λtr()I, stress tensor

˙

m = ∇ · (qτ

), (2.7)

where q(x) = Q(K) for x ∈ K, I is the identity matrix and λ the Lamé’s first parameter. The quality of the mesh is improved towards its goal of optimal shape as the elasticity problem is approaching a stationary solution.

2.3 Stabilized Galerkin Formulation

It is known that instabilities may arise when seeking numerical solutions to (1.27,1.28) for finite dimensional discretizations. A typical example for hyperbolic systems is the oscillations caused by convection dominated solutions [36]. Another example for mixed problems is violation of the Ladyshenskaya, Babuska, Brezzi (LBB) con- dition [6] if a non suitable pair of finite elements are used. Stabilization terms of streamline diffusion type are developed e.g. in [7], [35], [36], [37], [38].

The stabilized finite element formulation of the Unified Continuum Model taken from Paper IV is stated in the rest of this section.

Let Ω

⊂ R

be a polygonal domain. We introduce a sequence of discrete time steps 0 := t

< t

< · · · < t

:= ˆ t where I

:= [t

, t

) is a time interval of length k

:= t

− t

.

T

= {K} specifies the spatial discretization of Ω

and h

identifies the maxi-

mal diameter of the cell elements K ∈ T

. We also introduce the space-time slab

S

:= Ω

× I

.

The Sobolev-space H

(Ω

) and W

the finite dimensional function space of polynomials of continuous piecewise linear functions on S

are defined as follows:

H

(Ω

) : = {v ∈ L

(Ω

)| ∂v

∂x

∈ L

(Ω

) k = 1, 2, 3} (2.8) W

: = {v ∈ H

(S

)| v ∈ C

(S

), v ∈ P

(K ) × P

(I

), ∀K ∈ T

} (2.9)

W

: = {v ∈ W

| v = 0 on ∂Ω} (2.10)

W

: = [W

]

. (2.11)

The space of piecewise constant in space functions D

is defined as

D

: = {v ∈ L

(S

)| v ∈ P

(K ) × P

(I

), ∀K ∈ T

} (2.12)

D

: = [D

]

. (2.13)

We identify the discrete solution for velocity, pressure and solid stress as ˆ U = (U, P, τ

), the stress for both fluid and solid as T = −P I + θ2µ

(U) + (1 − θ)τ

, the discrete mesh velocity as M and the test function as ˆ v = (v, q, s).

We now formulate the spatially and temporally discretized variational formula- tion of the continuum model (2.1), (2.2), (2.3) based on these definition by using the midpoint quadrature rule in time, we obtain a Crank-Nicholson time stepping scheme: for each t

, find (U

, P

, τ

) := (U(t

), P (t

), τ

(t

)) with U

∈ W

, P

∈ W

and τ

∈ D

and T

= T (t

), such that:

(ρk

(U

− U

) + (ρ( ¯ U

− M

) · ∇) ¯ U

, v) + (T

: ∇v) + (∇ · ¯ U , q) (2.14) + SD

( ¯ U

, M

, P

, v, q, ρ) = 0,

(k

(τ

− τ

: s)) = (2µ

(u

) + ∇u

τ

+ τ

(∇u

)

: s)(2.15) for ∀(v, q, s) ∈ W

× W

× D

where ¯ U

=

(U

+ U

) and

(v, w) = X

K∈Tⁿ

Z

K

v · w dx. (2.16)

We use a simplified Galerkin/least-square method, to stabilize the convection dom- inated problem, where the time derivative term is dropped since the test functions are piecewise constant in time:

SD

( ¯ U

, M

, P

, v, q, ρ) = (δ

∇ · ¯ U

, ∇ · v)+ (2.17) (δ

ρ((( ¯ U

− M

) · ∇) ¯ U

+ ∇P

), ρ(( ¯ U

− M

) · ∇)v + ∇q) The stabilization parameters are chosen as δ

= κ

ρh|U

| and δ

= κ

ρ

(k

+

|U

− M

|

h

)

, where κ

, κ

are problem independent positive constants

of order O(1).

Contributions and Results

3.1 A Posteriori Error Estimation

As the dimension of the finite element subspace increases, the system of equations that should be solved numerically also grows. Combined with the restrictions on the time step size of the type of a Courant, Friedrichs, Lewy (CFL) condition [12], it becomes essential to employ a strategic method to optimize the mesh for controlling the error of a quantity of interest in particular.

For simplicity let A be a linear operator i the problem A(u) = f , for u the exact solution in a Hilbert space V . We are interested in the error in the form of the bounded functional M (u − U

), for the computed solution U

on a discretization with cell diameter h such that:

M (u − U

) = (u − U

, ζ)

A dual problem A

(φ) = ζ with initial and boundary conditions is constructed such that the bilinear equality (Au, v) = (u, A

v) is satisfied ∀u, v ∈ V . Then one can get:

|M (u − U

)| = |(u − U

, ζ)|

= |(u − U

, A

(φ))|

= |(A(u − U

), φ)|

= |(R(U

), φ)|

= |h(R(U

), h

(φ − π

φ))|

≤ kh(R(U

)kkCDφk

where π

φ is the linear interpolant of the dual solution in the finite dimensional solution subspace V

⊂ V and R(U

) = f − A(U

) is the residual. Due to the bilinear equality the dual equation has to be solved backwards in time if the primal

21

equation is time dependent. In the case of a non-linear operator, the adjoint (dual) operator is linearized at the u and U

.

The algorithm for adaptive mesh refinement based on a posteriori estimates becomes:

Algorithm 1 Adaptive Mesh Refinement Algorithm 1. solve the primal problem until final time.

2. solve the dual problem backwards in time.

3. calculate the error indicators using computed residuals and dual problem gradient.

4. if the sum of the indicators is below a given tolerance terminate.

5. refine a specified percentage of the cells with highest error indicator value.

6. return to step 1.

Paper II focuses on challenges in the derivation of the dual problem for the time dependent non-linear FSI problem, and solving it backwards in time with exactly the same time stepping and data of the primal problem, obtaining successful results for a serial computation.

Paper III shows an application of the method for a parallel implementation

obtaining good results, even though with simplifications in the dual problem and

data of the primal problem.

3.2 Contact Modeling

θ(t

) θ(t

)

U (t

) U (t

)

Figure 3.1: Visualizations of the vocal folds and quantities computed for the contact method, at two time points: t

- opened folds and t

- closed folds. Phase function θ (first row) and velocity function U at the mid plane (second row)

One difficulty for the simulation of human vocal folds is to sustain mesh quality when vocal folds get close to each other. A way to solve this problem is using re-meshing. In the Unified Continuum model the approach is solving an additional Eikonal equation to decide if the solids are close enough and to change the type of continuum from fluid to solid to sustain mesh quality. This approach is detailed in Paper IV. In the same paper it is also shown that the contact model is in particular suitable for describing the motion of the continuum for the self oscillating vocal folds problem as it qualitatively forms the glottal waves as shown in Figure 3.3.

The contact model is presented here as Algorithm 2 and its application to vocal

folds simulations is given in Figures 3.1, 3.2.

Algorithm 2 Unified Continuum Contact Algorithm 1. Mark all cells K as non-contact.

2. Solve the Eikonal equation |∇D

| = 1 for the discrete distance D

= D

(x), using an artificial viscosity stabilized cG(1) method with D

= 0 on the boundary of the fluid sub domain Ω

.

3. Compute |∇D

|, and define the medial axis M as: M = n x 

 |∇D

(x)| ≤ γ o , with the threshold parameter γ < 1.

4. Define the contact medial axis: ˆ M = n x 

 x ∈ M, x 6∈ Ω

, D

(x) < αˆ h o , with ˆ h the minimum cell size in the mesh.

5. Solve the Eikonal equation |∇D

| = 1 for the discrete distance D

, using an artificial viscosity stabilized cG(1) method with D

= 0 for the distance from ˆ M .

6. Mark all fluid cells as contact which fulfill: C = n x

 D

≤ βˆ h o

.

θ(t

) θ(t

)

D(t

) D(t

)

|∇D(t

)| |∇D(t

)|

D

(t

) D

(t

)

Figure 3.2: Visualizations of the vocal folds and quantities computed for the contact

method, at two time points: t

- opened folds and t

- closed folds.

a b c d

e f g h

a b c d

e f g h

a b c d

e f g h

Figure 3.3: Comparison of the glottal wave between a schematic figure (top row) of the expected contact pattern in the oscillatory cycle and simulations from [55], with a slice through the center of the vocal folds (second row), a clip through the solid phase (third row) and a clip through the solid phase together with volume rendering of the magnitude of the velocity (bottom row).

3.3 Acoustic Coupling

The computed flow field u and pressure field p from the FSI solver is used to con-

struct the acoustic source terms for a wave operator, following an acoustic analogy

approach. In an ALE setting, the standard wave equation becomes inappropriate for that purpose [22] and one has to resort to a mixed formulation [11] for the acoustic pressure, p

, and acoustic particle velocity, u

, to account for the moving boundaries (in our case the VF). Let Q

= Ω

× I stand for the space-time acoustic domain with β denoting again the ALE velocity. The equations to be solved are given by [22],

1

ρ

c

∂

p

− 1

ρ

c

β · ∇p

+ ∇ · u

+ α

p

= Q

in Q

, (3.1) ρ

∂

u

− ρ

β · ∇u

+ ∇p

+ α

u

= f

in Q

. (3.2) Q

(x, t) represents a volume source distribution and f

(x, t) an external body force per unit volume. To perform the one-way coupling between FSI and the acoustics equations, the terms Q

= (ρ

c

)

∂

p and f

= 0 have been chosen.

That corresponds to the mixed form of the analogy in [51], which as said, can be viewed as a simplification of the acoustic perturbation equations in [31, 23].

Equation (18) is to be supplemented with the following boundary and initial conditions

u

(x, t) · n = γ

p

x ∈ Γ

, u

(x, t) · n = 0 x ∈ Γ

,

u

(x, 0) = 0, (3.3)

corresponding to absorbing boundary conditions on walls Γ

and homogeneous Dirichlet boundary conditions on Γ

where we have taken γ

= 0.05/c

ρ

, c

= 350, ρ

= 1.225 for the numerical example in Paper VI.

The boundary conditions are imposed partly strongly and partly weakly as explained in variational form III of [3] The coupling methodology was successful for correctly reproducing the first 3 formants of vowel /i/ as presented in Paper VI

3.4 Implementation Details

The implementation is part of FEniCS-HPC, an open source framework for auto- mated solution of PDE on massively parallel architectures, providing automated evaluation of variational forms given a high-level description in mathematical no- tation, duality-based adaptive error control, implicit turbulence modeling by use of stabilized FEM and strong linear scaling up to thousands of cores as shown in Paper I, Paper V and [33, 40, 46, 28, 29]. FEniCS-HPC is a branch of the FEniCS [45, 17]

framework focusing on high performance on massively parallel architectures.

The framework is based on components with clearly defined responsibilities.

The main components are the following, with their dependencies shown in the dependency diagram in figure 3.4:

• The Finite Element Automated Tabulator (FIAT): Automated generation of

finite elements/basis functions