Mathematical and Numerical Modelling of the Respiratory System. - - PowerPoint PPT Presentation

Mathematical and Numerical Modelling of the Respiratory System. Céline Grandmont

• Enables gaz exchanges
• Pathologies:
• Emphysema
• Fibrosis
• Dust particles / pollution
• Asthma crisis

Weibel

}

Lung Tissue modelling Aerosol modelling

The Respiratory System

Healthy%

A%

Elastase;Treated%

B%

}

Coupling Lung tissue modelling Homogenization Aerosol modelling Air flow modelling

Methodology

• Elaborate representative models

= ⇒ Multiscale/Multiphysics modelling

• Give a mathematical framework

= ⇒ Wellposedness, long time behavior...

• Design efficient numerical schemes

= ⇒ Stability, convergence analysis. . .

• Perform relevant numerical simulations

= ⇒ Experimental validation

3D incompressible viscous flow Laminar flow Simplified elastic model

• L. Baffico, C. G., Y. Maday, B. Maury, A. Soualah

Related works:

[Quarteroni, Tuveri, Veneziani, 00], [Vignon-Cl´ ementel, Figueroa, Jansen, Taylor, 06]

Air Flow Modelling

Γ0 Γi

                     ρ∂u ∂t + ρ(u · ⌅)u µ⇤u + ⌅p = 0, in Ω , ⌅ · u = 0, in Ω , u = 0,

• n Γl ,

µ⌅u · n pn = P0n

• n Γ0 ,

µ⌅u · n pn = Πin

• n Γi

i = 1, . . . , N,

Air flow in the proximal part

In each subtree, we assume that the flow is laminar and thus satisfies a Poiseuille law: Πi Pi = Ri ⇤

Γi

u · n, Ri ⇤ 0, where Ri denotes the equivalent resistance of the distal tree and Pi is an alveola

• pressure. The boundary conditions at the outlets Γi writes

µ⌅u · n pn = Pin Ri ⇤

Γi

u · n ⇥ n on Γi i = 1, . . . , N.

Air flow in the distal part

All the alveola pressures are equal: Pi = Pa. The diaphragm and parenchyma motion are governed by: m¨ x = −kx + fext + SPa, Moreover, since the flow and the parenchyma are incompressible S ˙ x =

N

• i=1

⇥

Γi

u · n = − ⇥

Γ0

u · n,

Lung parenchyma - Diaphragm muscle

Ω

Γ

l

Γ Γi

Ω i

~

Coupled System

⇤ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⇧ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌃ ⌅ ρ∂u ∂t + ρ(u · ⌅)u µ⇤u + ⌅p = 0 , in Ω ⌅ · u = 0 , in Ω , u = 0 ,

• n Γl ,

µ⌅u · n pn = P0n

• n Γ0 ,

µ⌅u · n pn = Pan Ri

• Γi

u · n ⇥ n ,

• n Γi ,

i = 1 , . . . , N , m¨ x + kx = fext + SPa , S ˙ x =

N

⌥

i=1

• Γi

u · n =

• Γ0

u · n.

The Coupled System

Wellposedness

Energy Equality

d dt ✓ρ 2 Z

Ω

|u|2 + m 2 | ˙ x|2 + k 2|x|2 ◆ | {z } Total energy + µ Z

Ω

|ru|2 | {z } Dissipation within Ω +

N

X

i=1

Ri ✓Z

Γi

u · n ◆2 | {z } Dissipation in the subtrees =

• ρ

2

N

X

i=0

Z

Γi

|u|2(u · n) | {z } In/out-come of kinetic energy + P0S ˙ x | {z } Power of inlet pressure + fext ˙ x | {z } Power of ext. forces . and we have

• Z

Ω

(u · ru)u

• 

8 < : CkukL2(Ω)kruk2

L2(Ω), if d = 2,

Ckuk1/2

L2(Ω)kruk5/2 L2(Ω) if d = 3,

• p is replaced by the total pressure p+ρ|u|2/2, existence of weak solutions,

[C. G., B. Maury, Y. Maday, 05]

• u = λiUi on Γi, existence of weak solutions,

[C. G., B. Maury, A. Soualah, Esaim Proc, 08]

– locally in time for any data – for all time for small enough data

• existence of “strong” solutions,

[L. Baffico, C. G., B. Maury, m3as]

– locally in time for any data – for all time for small enough data Related works:

[Heywood, Rannacher, Tureck, 96], [Quarteroni, Veneziani, 03]

Wellposedness

Wellposedness

• u = λiUi on Γi,

and Z

Γi

σf(u, p) · nUi = Pa Z

Γi

Ui · n Ri ✓Z

Γi

Ui · n ◆2 then X

i

Z

Γi

|u|2(u · n) = X

i

(λi)3 Z

Γi

|Ui|2(Ui · n)  C P

i kuk3 X(Γi)

By taking X(Γi) = (H1/2(Γi))0 we obtain X

i

Z

Γi

|u|2(u · n)  Ckuk3

L2(Ω)

Wellposedness

General case

• Consider a special basis (wi)i for Galerkin approximation such that the

(wi)i are the eigenfunctions of the operator A defined by (Av, w)0 = µ Z

Ω

rv : rw + X

i

Ri ✓Z

Γi

v · n ◆ ✓Z

Γi

w · n ◆ , where (v, w)0 = ρ Z

Ω

vw + m S2 ✓Z

Γ0

v · n ◆ ✓Z

Γ0

w · n ◆ , for any v, w 2 {v 2 H1

0,Γ0(Ω), div v = 0}

• Take Au as a test function.

Wellposedness

ρ Z

Ω

∂tu·Au+ m S2 ✓Z

Γ0

∂tu · n ◆ ✓Z

Γ0

Au · n ◆ = µ 2 d dt Z

Ω

|ru|2+ X

i

Ri 2 d dt ✓Z

Γi

u · n ◆2 , µ Z

Ω

ru : rAu + X

i

Ri ✓Z

Γi

u · n ◆ ✓Z

Γi

Au · n ◆ = kAuk2

0,

Nonlinear term: Z

Ω

(u · r)uAu  kukL∞krukL2kAuk0, Question: Is kukL∞ controlled by kAuk0? Answer: Yes, since D(A) ⇢ H3/2+ε for some ε > 0 and there exists θ > 0 such that kukL∞  krukθ

L2kAuk1−θ

.

Wellposedness

• True under the assumption of right angles at the outlets ;
• Relies on regularity result for the solution of the Stokes problem in polyg-
• nal domains [Maz’ya, Rossman, 07] ;
• Nonlinear Gronwall lemma enables to conclud ;
• Additional estimate by taking ∂tu as a test function.

Numerical methods

[Ph.D thesis of J. Incaux-Fouchet]

• Numerical treatment on the artificial boundaries: What are the ”best”

boundary conditions?

• Numerical stability:

– How to deal with the numerical instabilities coming from the convec- tion term? – How to couple the 3D part with the 0D part?

Artificial boundary

Prescribe: pn + ru · n

Prescribe: pn + (ru + (ru)T ) · n

Prescribe: (p + |u|2

2 )n + ru · n

Numerical Stability

Instabilities observed for the Navier-Stokes equations with Neumann boundary conditions

= ⇒ Use of stabilization methods;

8 > > > > > > > > > > > > < > > > > > > > > > > > > : ⇢un+1 un ∆t + ⇢(un+1 · r)un+1 ⌘∆un+1 + rpn+1 = 0 in Ω, r · un+1 = 0 in Ω, un+1 = 0

• n Γ`,

⌘run+1 · n pn+1n = pn+1

in n

• n Γin,

⌘run+1 · n pn+1n = ✓ R Z

Γout

un+1 · n ◆ n

• n Γout,

u0 = u0 in Ω. (2.3.30)

Existence of strong solution for small enough data (initial data, applied force)

Semidiscretize scheme

8 > > > > > > > > < > > > > > > > > : ρ∂tu + ρ(u · r)u η∆u + rp = 0 in Ω, r · u = 0 in Ω, u = 0

• n Γ`,

σ · n = ηru · n pn = pinn

• n Γin,

σ · n = ηru · n pn = pi

• utn
• n Γi
• ut, i = 1, . . . , N,

u(0, ·) = u0(·) in Ω,

pi

• ut(t) = Fi(Qi(t)) = Fi

Z

Γi

• ut

u(s, ·) · n ! , with 0  s  t,

Fi(Qi) = RiQi + 1 Ci Z t Qi = Ri Z

Γi

• ut

u · n + 1 Ci Z t Z

Γi

• ut

u · n,

j and i j if the 3D/0D-interface i

is localized at the genera

Numerical Treatment of 3D/0D Coupling

with

pout R C Q p = 0 pout = RQ + 1 C Z t Q Figure 2.2: The RC reduced model.

Lung Blood

pout Rp C Q Rd Pd Pp

Fi(Qi(t)) = Ri

pQi(t) + P i d,0 e−t/τi + 1

Ci Z t Qi(s)e(s−t)/τids R

• Stokes + implicit treatment =

⇒ unconditional stability

• Stokes + explicit treatment =

⇒ conditional stability – R model: is unconditionally stable but with an upper bound that behaves as exp( R2T

ν ). Moreover if

R > Kν,

• r ∆t ≤ K ρ

R we obtain a better upper bound. – RC model: ∆t ≤ K ρC RC + T – RCR model: same type of sufficient conditions than for the R model.

Parameters identification

Question Are we able to recover the resistances, stiffness parameters thanks to measurements of exhaled volume and flow rate at the mouth ? Partial answer [ Ph.D. thesis of A-C. Egloffe]

• Consider simplified problem ;
• Try to prove stability estimates of the unkown parameters with respect to

the data measurements.

Parameters identification

Steady State Stokes problem:        ν∆u + ru = f, in Ω divu = in Ω ν ∂u

∂n pn

= g,

• n Γe

ν ∂u

∂n pn + ku

= 0,

• n Γ0

We would like to estimate kk1 k2kL2(Γ0) with respect to the difference u1 u2 and p1 p2 on some part of the boundary. We are going to estimate kk1 k2kL2(K), where K is a compact set of Γ0 such that u1 is not equal to zero on K. We have that kk1 k2kL2(K)  C(ku1 u2k1/2

H1(Ω) + kp1 p2k1/2 H1(Ω))

Related work: [Phung, COCV]

Parameters identification

Let 0 < ⌫  1

2 and ! a nonempty open set included in Ω. Then, for all

2

• 0, 1

2 + ⌫

• , there exists c > 0, such that for all ✏ > 0, we have

kukH1(Ω) + kpkH1(Ω)  e

c ✏

kukL2(Γ) + kpkL2(Γ) +

• @u

@n

• L2(Γ)

+

• @p

@n

• L2(Γ)

! + ✏β(kukH

3 2 +⌫(Ω) + kpkH 3 2 +⌫(Ω)),

and kukH1(Ω)+kpkH1(Ω)  e

c ✏

kukH1(ω) + kpkH1(ω)

• +✏β(kukH

3 2 +⌫(Ω)+kpkH 3 2 +⌫(Ω)),

for all couple (u, p) 2 H

3 2 +ν(Ω) ⇥ H 3 2 +ν(Ω) solution of the Stokes problem.

Proof

• Based on local Carleman inequalities applied to u such that ⌫∆u = rp

and to p such that ∆p = 0;

• Idea: estimate the H1 norm in any subset included in Ω;

Parameters identification

It is equivalent to: Let 0 < ν  1

• 2. Let Γ be a nonempty part of the boundary of Ω and ω be a

nonempty open set included in Ω. Then, there exists d0 > 0 such that for all β 2

• 0, 1

2 + ν

• , for all d > d0, there exists c > 0, such that we have

kukH1(Ω) + kpkH1(Ω)  c kukH

3 2 +ν(Ω) + kpkH 3 2 +ν(Ω)

✓ ln ✓ d

kuk

H 3 2 +ν (Ω)

+kpk

H 3 2 +ν (Ω)

kukL2(Γ)+kpkL2(Γ)+k ∂u

∂nkL2(Γ)+k ∂p ∂nkL2(Γ)

◆◆β , (1) and kukH1(Ω) + kpkH1(Ω)  c kukH

3 2 +ν(Ω) + kpkH 3 2 +ν(Ω)

✓ ln ✓ d

kuk

H 3 2 +ν (Ω)

+kpk

H 3 2 +ν (Ω)

kukH1(ω)+kpkH1(ω)

◆◆β , for all couple (u, p) 2 H

3 2 +ν(Ω) ⇥ H 3 2 +ν(Ω) solution the Stokes problem.

Pressure fields

Numerical simulations

Streamlines

Numerical simulations

• 25
• 20
• 15
• 10
• 5

5 10 15 20 1 2 3 4 5 6 7 8 9 10 ’bronche_ns_res0res3_pg0_kx_T20_piston’ using 1:2

Comparison with experimental data

Cardiac' Apical' Le-' Diaphragma2c' Intermediate'

CL' RL' CC' RC' C

R

CI' RI' CA' RA'

P(t)'

Cardiac' Apical' Le-' Diaphragma2c' Intermediate'

C

L'

R

L'

C

C '

R

C '

CD' RD' CI' RI' CA' RA'

P(t)'

Figure'1:'3D,0D'airflow' model'with'outlets'directed' towards'the'5'rat'lobes.'

R ! V + V C = P ! P

peep

• V"="breathing"volume"
• P"="experimental"pressure"
• Ppeep"="1"cm"H2O"

0D Model

0.1 0.2 0.3 0.4 0.5 0.6 0.7 2 4 6 8 10

time, sec pressure, cm H2O

Healthy Emphysema

! !"# !"$ !"% ! & #

'()*+,-*. /012)*+,)3

! !"# !"$ !"% !4! !#! !&! ! &!

'()*+,-*. 5106,78'*+, )39-*.

, , ! !": & &": # !4! !#! !&! ! &!

/012)*+,)3 5106,78'*+, )39-*.

! !": & &": # ! : &!

/012)*+,)3 ;7*--27*+, .),<#= <*81'>? @);>?-*)8

CG, A. Marsden, J. Oakes, I. Vignon-Clementel

Experimental Curves 3D Model

D´ epˆ

• t de particules

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Left Apical Intermediate Diaphragmatic Cardiac

Normalized Particle Deposition

• r Delivery / Alpha

Healthy

Experimental Data Simulated

Homogeneous emphysema/ heterogeneous emphysema

Numerical simulations

3D structure of the flow

3D Viscoelastic model.

[Cazeaux, C. G.] 8 > > > > < > > > > : ∂ttdε − div(σ(dε)) = fS, in Ωε, σ(dε)nε + X

j

Aε

kj

✓Z

Γε,j

∂tdε · nε dγ ◆ nε = 0,

• n Γε,k, ∀k ∈ ZN

ε ,

σ(dε, qε)n = 0,

• n ΓN,

dε = 0,

• n ΓD,

Additional Term:

• Ω
• Ω

Kε(x, x) div∂tdε(x) divφ(x)dxdx

1

F

Y

S

n ! Y

3D Model

r0 r1 r1 r2 r2 r2 r2 rN 1 X0 X1 X2 X3 XN−1 XN Figure 1. Dyadic tree

Hypothesis

• Diadic tree of pipes
• In each tube the Poiseuille law is satisfied =

⇒ the flow is caracterized by a resistance rn at each generation

R0 R0 R1 R1 R1 R1 R2 R2 R2 R2 R2 R2 R2 R2 R3 R3 r0 r1 r1 r2 r2 r2 r2 + + + + = . . . Figure 1. Kernel K(x, y)

1D example

pi =

2N−1

X

j=0

qjRN−νij, with Rn =

n

X

j=0

rj,

Limit problem

dε(x) = d0(x) + εd1(x, x/ε) + ε2d2(x, x/ε) + ... Cell problem - compressible case

8 > > > > > > < > > > > > > : −div yσy(d1) = 0 in Ω × YS σy(d1)n + Z

Ω

K(x, x) „Z

∂YF

∂td1(t, x, y) · n dy « dxn = |YF | Z

Ω

K(x, x)div∂td0(t, x)dxn − σ(d0)n, on Ω × ∂YF d1 Y -periodic.

Cell problem - incompressible case

8 > > > < > > > : −div yσy(d1, q0) = 0 in Ω × YS div yd1 = −div xd0 σy(d1, q0)n = |Y | Z

Ω

K(x, x)div∂td0(t, x)dxn − σ(d0)n, on Ω × ∂YF d1 Y -periodic.

3D Model

d1(x, y, t) = X

16k,l6d

ex(d0)kl(x, t)χkl(y) + λ(x, t)u(y) where χkl are the standard correctors (for a perforated elastic media) and u is corrector related to the pressure effect of the fluid on the structure. The vector λ depends on d0 and (assuming d0 is known and x is a parameter) satisfies a first order ODE. = ⇒ apparition of time delay terms. = ⇒ Limit model: nonlocal in time and space. Remark: In the case of an incompressible elastic medium, there is no time memory effect at the limit.

Comparison with experimental data

• 1000
• 800
• 600
• 400
• 200

200 400 600 800 1000

• 180
• 160
• 140
• 120
• 100
• 80
• 60
• 40
• 20

20 reference shear modulus 100% increase on one side

• 1
• 8
• 6
• 4
• 2

2 4 6 8 1

• 1

8

• 1

6

• 1

4

• 1

2

• 1
• 8
• 6
• 4
• 2

2 r e f e r e n c e r e s i s t a n c e s 1 x i n c r e a s e i n d i s t a l r e s i s t a n c e s

• n
• n

e s i d e

∂tu + (u · ∇)u − µu + ∇p = 0 µ∇u · n − pn = −Πin µ∇u · n − pn = 0 Πi − Pi = Ri∇ · ∂tη R

Γi u · n =

R

Ωi ∇ · ∂tη

∂ttη − ∇ · σ(η) + ∇P = 0 Pi u= 0 Ωi σ(η) · n = −Pmusn

Full Model