Modelling and simulation of mechano-chemical fluid-structure - - PowerPoint PPT Presentation

Modelling and simulation of mechano-chemical fluid-structure interaction with application to atherosclerotic plaque growth

Stefan Frei University College London Thomas Richter Otto von Guericke University Magdeburg Thomas Wick Leibniz University Hannover Modeling, Simulation and Optimization of the Cardiovascular System Magdeburg, October 24, 2018

Gefördert durch

Plaque growth in blood vessels

Mechano-chemical fluid-structure interaction Solid growth depending on concentration of monocytes/foam cells Migration of monocytes depends on wall shear stress

Transport of Monocytes Transendothelial migration and differentiation Formation of foam cells Plaque Growth Ωs(t) Ωs(t) Ωf (t) Γi(t) Γi(t)

Difficulties Strongly-coupled FSI problem: ρf ≈ ρs Large deformations up to full clogging Different time scales of fluid dynamics (milliseconds-seconds) and plaque growth (days-months)

• S. Frei | Mechano-chemical FSI

2

Overview

1

Model

2

Temporal two-scale approach

3

Fluid-structure interaction

4

Discretisation

5

Numerical results

• S. Frei | Mechano-chemical FSI

3

Overview

1

Model

2

Temporal two-scale approach

3

Fluid-structure interaction

4

Discretisation

5

Numerical results

• S. Frei | Mechano-chemical FSI

4

Equations

Mechano-chemical FSI system based on Yang et al. 2016 (simplified) ρf (∂tvf + vf · ∇vf ) − div σf = 0 div vf = 0

• in F(t)

ρs(∂tvs + vs · ∇vs) − div σs(cs)= 0 ∂tus + vs · ∇us − vs= 0

• in S(t)

σf nf + σs(cs)ns = 0 vf = vs

• n Γi(t)

∂tcs = γ(σWS

f

) Material laws (simplifications) Newtonian fluid: σf = ρf νf (∇v + ∇T v) − pf I Solid material law based on St.Venant-Kirchhoff material σs(cs) Dependence of the foam cell concentration cs on the wall shear stress σWS

f

γ(σWS

f

) = γ0

• 1 + σWS

f

σ

−1

, σWS

f

=

• Γi

|σf nf · e1| do

• S. Frei | Mechano-chemical FSI

5

Coupling FSI-Chemistry

Concentration of foam cells cs depends on the wall shear stress σWS

f

Feedback: Multiplicative decomposition of the deformation gradient Fs = Fg(cs)Fe (Rodriguez et al. 1994)

σWS

f

cs FSI (ODE) Fs = Fg(cs)Fe

Idea ˆ Fs = I + ˆ ∇ˆ us ˆ V ˆ Vg V Fs = I − ∇us ˆ Fe = F −1

e

ˆ Fg = F −1

g

• S. Frei | Mechano-chemical FSI

6

Solid material law with growth in Lagrangian coordinates

Multiplicative decomposition of the deformation gradient ˆ Fe = ˆ Fs ˆ F −1

g

(cs) = (I + ˆ ∇ˆ us)ˆ F −1

g

(cs)

• St. Venant-Kirchhoff in Lagrangian coordinates

ˆ Ee = 1 2 (ˆ F T

e ˆ

Fe − I), ˆ Σe = 2µstr(ˆ Ee) + λs ˆ Ee. Solid system of equations ˆ ρ0

s ∂tˆ

vs − div ˆ Fe(cs)ˆ Σe(cs) = 0 ∂tˆ us − ˆ vs= 0 in ˆ S

Yang et al (2016), S.F., T.Richter, T.Wick, JCP (2016)

• S. Frei | Mechano-chemical FSI

7

Solid material law with growth in Eulerian coordinates

Multiplicative decomposition of the deformation gradient Fe = F −1

g

(cs)Fs = F −1

g

(cs)(I − ∇us)

• St. Venant-Kirchhoff in Eulerian coordinates

Ee = 1 2 (F −T

e

F −1

e

− I), Σe = 2µstr(Ee) + λsEe, σs = JsF −1

e

ΣeF −T

s

, where Js = det Fs. Density: ρ = J−1

g

(cs)Jsρ0

s , where Jg= det Fg.

Solid system of equations J−1

g

(cs)Jsρ0

s (∂tvs + vs · ∇vs) − div

JsF −1

e

(cs)Σe(cs)F −T

s

• = 0

∂tus + vs · ∇us − vs= 0 in S(t)

S.F., T.Richter, T.Wick, JCP (2016)

• S. Frei | Mechano-chemical FSI

8

Overview

1

Model

2

Temporal two-scale approach

3

Fluid-structure interaction

4

Discretisation

5

Numerical results

• S. Frei | Mechano-chemical FSI

9

Time scales

Different time scales for fluid dynamics (tshort = 1s, period of inflow data) and plaque growth (Tlong = 1 month) To resolve the short scale (k ≤ 1/20s) over Tlong would require > 107 time-steps Assumption: cs is approximately constant on the short scale Consider the averaged quantity ¯ cs = t

t−1s cs(r) dr

There holds ∂t¯ cs =

t

t−1s

∂tcs(r) dr =

t

t−1s

γ(σWS

f

(v, p)) dr =: ¯ γ(σWS

f

(v, p)). Problem: No initial values available for u, v, as we can not resolve the short scale We assume the existence of a periodic solution vcs (t − 1s) = vcs (t) for fixed cs For Stokes flow

• v(t) − vcs(t)
• H1(Ω) +

p(t) − pcs(t)

• L2(Ω) = O(α),

α = tshort Tlong ∂t¯ cs ≈ ¯ γ(σWS

f

(vcs , pcs ))

• S. Frei | Mechano-chemical FSI

10

Time scales

Different time scales for fluid dynamics (tshort = 1s, period of inflow data) and plaque growth (Tlong = 1 month) To resolve the short scale (k ≤ 1/20s) over Tlong would require > 107 time-steps Assumption: cs is approximately constant on the short scale Consider the averaged quantity ¯ cs = t

t−1s cs(r) dr

There holds ∂t¯ cs =

t

t−1s

∂tcs(r) dr =

t

t−1s

γ(σWS

f

(v, p)) dr =: ¯ γ(σWS

f

(v, p)). Problem: No initial values available for u, v, as we can not resolve the short scale We assume the existence of a periodic solution vcs (t − 1s) = vcs (t) for fixed cs For Stokes flow

• v(t) − vcs(t)
• H1(Ω) +

p(t) − pcs(t)

• L2(Ω) = O(α),

α = tshort Tlong ∂t¯ cs ≈ ¯ γ(σWS

f

(vcs , pcs ))

• S. Frei | Mechano-chemical FSI

10

Time scales

Different time scales for fluid dynamics (tshort = 1s, period of inflow data) and plaque growth (Tlong = 1 month) To resolve the short scale (k ≤ 1/20s) over Tlong would require > 107 time-steps Assumption: cs is approximately constant on the short scale Consider the averaged quantity ¯ cs = t

t−1s cs(r) dr

There holds ∂t¯ cs =

t

t−1s

∂tcs(r) dr =

t

t−1s

γ(σWS

f

(v, p)) dr =: ¯ γ(σWS

f

(v, p)). Problem: No initial values available for u, v, as we can not resolve the short scale We assume the existence of a periodic solution vcs (t − 1s) = vcs (t) for fixed cs For Stokes flow

• v(t) − vcs(t)
• H1(Ω) +

p(t) − pcs(t)

• L2(Ω) = O(α),

α = tshort Tlong ∂t¯ cs ≈ ¯ γ(σWS

f

(vcs , pcs ))

• S. Frei | Mechano-chemical FSI

10

Averaging of FSI

Averaging for linear FSI under periodicity 0 =

t

t−1s

∂tvcs (r) + A(ucs (r), vcs (r), pcs (r)) dr = vcs (t) − vcs (t − 1s) + A(¯ ucs , ¯ vcs , ¯ pcs ) = A(¯ ucs , ¯ vcs , ¯ pcs ). For non-linear FSI, additional terms would arise on the right, that will be neglected in this work Averaged system of equations (T) A(¯ ucs , ¯ vcs , ¯ pcs ) = 0 ∂t¯ cs = ¯ γ(σWS

f

(vcs , pcs )) Short-scale influence enters the right-hand side of the ODE

• S. Frei | Mechano-chemical FSI

11

Two-scale algorithm

Strategy to design a two-scale algorithm Discretise the total time interval in sub-intervals Im = [Tm, Tm+1], K = Tm − Tm−1 Discretise for each m a short scale interval [Tm, Tm + 1s] into sub-intervals in = [tn, tn+1], k = tn+1 − tn For m = 1, ..., M: Solve a short-scale problem on [Tm, Tm + 1s] (non-stationary FSI) and compute ¯ γ(σWS

f

) Advance cs by solving the ODE and (v, p, u) by solving the averaged stationary FSI problem In practice: Iterate on short-scale until v(tm + 1s) ≈ v(tm) for periodicity Theoretically O(k2 + αK 2 + α) for a Stokes system with growing boundary and second-order time-stepping (S.F., T. Richter, in preparation), open for nonlinear FSI

• S. Frei | Mechano-chemical FSI

12

Two-scale and pure long-scale approach

We will compare our approach to a pure long-scale approach (a) Pure long-scale algorithm: A(¯ u, ¯ v, ¯ p) = 0 ∂t¯ cs = γ(σWS

f

(¯ v, ¯ p))

FSI (ODE) ¯ σWS

f

¯ cs F = Fg(¯ cs)Fe

(b) Temporal two-scale-algorithm: A(¯ u, ¯ v, ¯ p) = 0 ∂t¯ cs = ¯ γ(σWS

f

(vcs , pcs ))

Long ODE Short Stat. Nonstat. FSI FSI ¯ cn

s

γ(σWS

f

) ¯ cn

s , ¯

un−1

s

Fs = Fg(¯ cn

s )Fe

• S. Frei | Mechano-chemical FSI

13

Short-scale and long-scale problem

Short-scale (t) problem in [Tm, Tm + 1s] Periodic initial values (iteration), concentration cs fixed Inflow vin

f ,x(t) = cω(1 + sin(2πt)), initial displacement d(Tm)

ρf (∂tv + v · ∇v) + div σf = 0, div vf = 0 in F(t), JsJ−1

g

(cs)ρ0

s (∂tv + v · ∇v) + div

JsF −1

e

(cs)Σe(cs)F −T

s

• = 0,

∂tu + v · ∇u − v= 0 in S(t), σf nf + JsF −1

e

(cs)Σe(cs)F −T

s

= 0, vf = vs in Γi(t) Compute averaged growth function: γ(σWS

f

) = Tm+1s

Tm

γ(σWS

f

) dt, γ(σWS

f

) = γ0

• 1 +

σWS f σ0

−1

, σWS

f

=

Γi

|σf nf · e1| do Long-scale (T) problem in Im = [Tm, Tm+1], e.g. K = Tm=1 − Tm = 1 day Accumulation of foam cells: ∂t¯ cs = ¯ γ(σWS

f

). Averaged inflow vin

f ,x = cω

ρf ¯ vf · ∇¯ vf + div ¯ σf = 0, div ¯ vf = 0 in F(T) div¯ Js ¯ F −1

e

(¯ cs)¯ Σe(¯ cs)¯ F −T

s

• = 0

in S(T) ¯ σf nf + ¯ Js ¯ F −1

e

(¯ cs)¯ Σe(¯ cs)¯ F −T

s

ns= 0, ¯ vf = ¯ vs in Γi(T)

• S. Frei | Mechano-chemical FSI

14

Overview

1

Model

2

Temporal two-scale approach

3

Fluid-structure interaction

4

Discretisation

5

Numerical results

• S. Frei | Mechano-chemical FSI

15

Coupling approaches for FSI

Arbitrary Lagrangian Eulerian (ALE) method: Solid equations in Lagrangian coordinates, fluid equations an arbitrary (matching) reference frame Immersed approaches (Solid in Lagrangian, fluid in Eulerian coordinates) Fully Eulerian Approach (Dunne 2006): Monolithic Eulerian formulation

• S. Frei | Mechano-chemical FSI

16

Fully Eulerian approach: Variational formulation

Monolithic variational formulation: Find global velocity v ∈ vd + V, V := H1

0(Ω; Γd v )

solid displacement us ∈ ud

s + Ws, Ws := H1 0(Ωs(t); Γd u)

fluid pressure pf ∈ Lf := L2

0(Ωf (t)) Γi(t) Γi(t) Ωf (t) Ωs(t) Ωs(t)

such that (ρ(∂tv + v · ∇v), φ)Ω + (σ, ∇φ)Ω = (ρf , φ)Ω ∀φ ∈ V, (∂tus + v · ∇us − v, ψs)Ωs(t) = 0 ∀ψs ∈ Ws, (div v, ξf )Ωf (t) = 0 ∀ξf ∈ Lf . Coupling conditions included by means of variational principles Global velocity v ∈ vD + V → vf |Γi (t) = vs|Γi (t) Continuity of normal stresses by global test space: φ ∈ V

Dunne & Rannacher 2006, Cottet, Maitre & Milcent 2006

• S. Frei | Mechano-chemical FSI

17

Initial point set/Backward characteristics

How to capture the moving interface? Initial point set function traces back points to their initial position: ΦIPS(x, t) :=

• x − us(x, t)

x ∈ Ωs(t), x − ext(us)(x, t) x ∈ Ωf (t).

Ωs(t) x ˆ Ω = Ωs(0) ˆ x φIPS(x) = x − us(x)

The domain affiliation is given by: ΦIPS(x, t) ∈ Ωs(0) ⇔ x ∈ Ωs(t), ΦIPS(x, t) ∈ Ωf (0) ⇔ x ∈ Ωf (t). Capable to capture sharp interfaces No additional equation has to be solved Good mass conservation properties

• S. Frei | Mechano-chemical FSI

18

Fully Eulerian approach: Advantages and challenges

Advantages: Monolithic formulation: Strongly-coupled problems, Error control Large displacements up to contact are possible Challenges for the discretisation: Solid and interface move over grid cells Global velocity v ∈ H1(Ω), but ∇v, ∂tv discontinuous on Γi(t) Reduced convergence and stability issues for standard FE discretisation (not considering the interface position)

• S. Frei | Mechano-chemical FSI

19

Overview

1

Model

2

Temporal two-scale approach

3

Fluid-structure interaction

4

Discretisation

5

Numerical results

• S. Frei | Mechano-chemical FSI

20

Finite element discretisation: A locally modified finite element method

Fixed, regular patch mesh independent of the interface location Split interface patches into 8 triangles to resolve the interface

Γh Ω1 Ω2 Γ

Piecewise linear FE in interface patches, piecewise bilinear FE else Simple implementation

Only standard basis functions Local number of dofs is always 9 → Sparsity pattern of the system matrix independent of the interface location

• S. Frei | Mechano-chemical FSI

21

Locally modified FEM

4 different types of cut patches s r r s s Anisotropic cells might occur, e.g. r → 0

• Γi

Maximum angle condition: All interior angles bounded by 3π/2 < 2π. → Nodal interpolation of second order O(h2

P)

→ Optimal order error estimates (S.F., T.Richter, SINUM (2015)) Condition number bounded by O(h−2

P ) using a hierarchical FE basis

FSI: Equal-order locally modified finite elements with an anisotropic variant of the interior-penalty pressure stabilisation (S.F., arXiV (2018))

• S. Frei | Mechano-chemical FSI

22

Locally modified FEM

4 different types of cut patches s r r s s Anisotropic cells might occur, e.g. r → 0

• Γi

Maximum angle condition: All interior angles bounded by 3π/2 < 2π. → Nodal interpolation of second order O(h2

P)

→ Optimal order error estimates (S.F., T.Richter, SINUM (2015)) Condition number bounded by O(h−2

P ) using a hierarchical FE basis

FSI: Equal-order locally modified finite elements with an anisotropic variant of the interior-penalty pressure stabilisation (S.F., arXiV (2018))

• S. Frei | Mechano-chemical FSI

22

Time discretisation

Solution v not smooth across Γi(t) (in space and time) Ansatz: Continuous/discontinuous Galerkin approach (cG(r)/dG(r)) in time: X cG,r

k

= v ∈ C(I, V )

• v|Im ∈ Pr(Im, V )

, X dG,r

k

= v ∈ L2(I, V )

• v|Im ∈ Pr(Im, V )

. Find wk ∈ X cG,r+1

k

• X dG,r

k

• s.t. B(wk, φk) = (f , φk)

∀φk ∈ X dG,r

k

Problem: Approach does not consider the reduced regularity across Γi(t) → wk ∈ C∞ across Γi(t), but v only H1(Q).

x t Γi(t) Im+1 Im Im−1

• S. Frei | Mechano-chemical FSI

23

Time discretisation

Solution v not smooth across Γi(t) (in space and time) Ansatz: Continuous/discontinuous Galerkin approach (cG(r)/dG(r)) in time: ˜ X cG,r

k

= v ∈ C(I, V )

(v◦Tm)|Im ∈ Pr(Im, V )

, ˜ X dG,r

k

= v ∈ L2(I, V )

(v◦Tm)|Im ∈ Pr(Im, V )

. Find wk ∈ ˜ X cG,r+1

k

˜

X dG,r

k

• s.t. B(wk, φk) = (f , φk)

∀φk ∈ ˜ X dG,r

k

Idea: Define polynomials on trajectories Tm inside the subdomains Tm(Ωm−1

i

, t) = Ωi(t), Tm(Γm−1

i

, t) = Γ(t).

x t Im+1 Im Im−1 Tm

Second-order convergence for the modified CG(1)-approach (S.F. & T.Richter, ESAIM M2AN, 2017)

• S. Frei | Mechano-chemical FSI

23

Overview

1

Model

2

Temporal two-scale approach

3

Fluid-structure interaction

4

Discretisation

5

Numerical results

• S. Frei | Mechano-chemical FSI

24

Configuration

Γi(t) Ωs(t) Γout

f

Γd

s

Γd

s

Γd

s

Γd

s

Γin

f

Ωf (t) Γi(t) Γd

s

Γd

s

Boundary conditions vf = vin

f

• n Γin

f ,

ρf νf ∂nv − pn = 0 on Γin

f ,

us = 0 on Γd

s

Parameters Fluid parameters: ρf = 1

g cm3 , νf = 0.3 cm2 s

Solid material parameters: ρs = 1

g cm3 , µs = 104 dyn cm2 , λs = 4 · 104 dyn cm2

Growth function: γ(σWS

f

) = γ0

• 1 +

σWS

f

σ0

−1

, σ0 = 50

g cm s2 , γ0 = 5 · 10−7 = O(α)

Growth tensor: ˆ Fg(ˆ x, τ) = 1 + cs(T) exp(−ˆ x2

1 )(2 − |ˆ

x2|)I

• S. Frei | Mechano-chemical FSI

25

Numerical parameters

Fully Eulerian approach Equal-order locally modified FE (P1/Q1) with anisotropic edge-oriented stabilisation Modified dG(0) time discretisation ALE approach for comparison Equal-order Q1 elements with LPS stabilisation Implicit Euler time-stepping Harmonic/Biharmonic extension for the ALE map Temporal two-scale approach Time steps: k = 0.02s, K = 0.1 days = 8640s Periodicity: 3 iterations to get (approximately) periodic solution

• S. Frei | Mechano-chemical FSI

26

Example with full clogging

Pure long-scale computation Parabolic inflow profile with mean value vin

f ,x(T) = 5ω(T) cm/s

depending on the minimal channel width ω(T). (Plaque growth)

• S. Frei | Mechano-chemical FSI

27

Comparison with ALE calculations

Meshes at T ≈ 40 days for ALE (left) and Fully Eulerian (right) Functionals: Channel width and vorticity Jvort =

F(T)(∂yv1 − ∂xv2)2 dx over time ALE Euler Channel width over time in cm 60 50 40 30 20 10 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 ALE Euler Vorticity over time in cm2/s2 60 50 40 30 20 10 70 60 50 40 30 20 10

• S. Frei | Mechano-chemical FSI

28

Convergence

If we do not allow the inflow to vanish vin

f ,x(T) = 0.1 + 5ω(T) cm/s

the channel does not close completely:

ALE, harmonic ALE biharmonic Euler Channel width over time in cm 180 160 140 120 100 80 60 40 20 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 ALE harmonic ALE biharmonic Euler Vorticity over time in cm2/s2 180 160 140 120 100 80 60 40 20 7500 7000 6500 6000 5500 5000 4500 4000 3500 3000

Convergence studies in space at T = 50 days (Fully Eulerian) #patches Wall Stress Width Vorticity Outflow Euler 256 1.033 · 102 1.092 3.408 · 103 9.251 1024 1.050 · 102 1.064 3.457 · 103 9.547 4096 1.060 · 102 1.052 3.472 · 103 9.648 Extrapol. 1.074 · 102 1.047 3.479 · 103 9.700 Conv. 0.77 1.81 1.71 1.55

• S. Frei | Mechano-chemical FSI

29

Pure long-scale vs two-scale algorithm

Growth function γ(σWS

f

) of the short-scale problem over δt = 3s Strong underestimation in pure long-scale simulation Third period approximately periodic

1 × 10−7 2 × 10−7 3 × 10−7 4 × 10−7 5 × 10−7 0.5 1 1.5 2 2.5 3 long- term avg. short- term Averaging of the growth function short-scale Euler short-scale ALE Long-Short, ALE Long, ALE Long-Short, Euler Long, Euler Wall shear stress in g/s2 200 150 100 50 1800 1600 1400 1200 1000 800 600 400 200 Long-Short, ALE Long, ALE Long-Short, Euler Long, Euler Channel width in cm 200 150 100 50 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2

• S. Frei | Mechano-chemical FSI

30

Conclusion and outlook

Conclusion Models for mechano-chemical FSI with solid growth in Eulerian and Lagrangian coordinates Two-scale algorithm in time to deal with different time-scales Fully Eulerian approach for FSI with large displacements and contact Accurate and robust discretisation techniques in space and time for interface problems Outlook More realistic models for fluid, solid and plaque growth Analysis of the temporal two-scale approach Consideration of ruptures

Main references: S.F., T. Richter, T.Wick: Long-term simulation of large deformation, mechano-chemical fluid-structure interactions in ALE and fully Eulerian coordinates, J Comp Phys 321 (2016), p. 874-891 S.F: Eulerian finite element methods for interface problems and fluid-structure interactions, PhD thesis, Heidelberg University (2016)

• S. Frei | Mechano-chemical FSI

31