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

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

The Respiratory System • Enables gaz exchanges • Pathologies: - Emphysema } Lung Tissue modelling - Fibrosis - } Dust particles / pollution Aerosol modelling - Weibel Asthma crisis A % B % Healthy% Elastase;Treated%

Homogenization Coupling Lung tissue modelling Air flow modelling Aerosol modelling

Methodology • Elaborate representative models = ⇒ Multiscale/Multiphysics modelling • Give a mathematical framework = ⇒ Wellposedness, long time behavior... • Design e ffi cient numerical schemes = ⇒ Stability, convergence analysis. . . • Perform relevant numerical simulations = ⇒ Experimental validation

Air Flow Modelling L. Baffico, C. G., Y. Maday, B. Maury, A. Soualah Simplified elastic 3D incompressible Laminar flow model viscous flow Related works: [Quarteroni, Tuveri, Veneziani, 00], [Vignon-Cl´ ementel, Figueroa, Jansen, Taylor, 06]

Air flow in the proximal part Γ 0 ρ∂ u  ∂ t + ρ ( u · ⌅ ) u � µ ⇤ u + ⌅ p = 0 , in Ω ,       = 0 , in Ω , ⌅ · u     = 0 , on Γ l , u Γ i    µ ⌅ u · n � p n = � P 0 n on Γ 0 ,       = on Γ i i = 1 , . . . , N, µ ⌅ u · n � p n � Π i n 

Air flow in the distal part In each subtree, we assume that the flow is laminar and thus satisfies a Poiseuille law: ⇤ Π i � P i = R i u · n , R i ⇤ 0 , Γ i where R i denotes the equivalent resistance of the distal tree and P i is an alveola pressure. The boundary conditions at the outlets Γ i writes �⇤ ⇥ µ ⌅ u · n � p n = � P i n � R i n on Γ i i = 1 , . . . , N. u · n Γ i

Lung parenchyma - Diaphragm muscle All the alveola pressures are equal: P i = P a . The diaphragm and parenchyma motion are governed by: m ¨ x = − kx + f ext + SP a , Moreover, since the flow and the parenchyma are incompressible N ⇥ ⇥ � S ˙ x = u · n = − u · n , Γ i Γ 0 i =1

The Coupled System Γ 0 Γ i Γ Ω l ~ Ω i Coupled System ρ∂ u ⇤ ∂ t + ρ ( u · ⌅ ) u � µ ⇤ u + ⌅ p = 0 , in Ω ⌃ ⌃ ⌃ ⌃ = 0 , in Ω , ⌅ · u ⌃ ⌃ ⌃ ⌃ = 0 , on Γ l , ⌃ u ⌃ ⌃ ⌃ µ ⌅ u · n � p n = � P 0 n on Γ 0 , ⌃ ⌃ ⌃ ⌃ �� ⇥ ⇧ = on Γ i , µ ⌅ u · n � p n � P a n � R i n , u · n Γ i ⌃ ⌃ i = 1 , . . . , N , ⌃ ⌃ ⌃ ⌃ m ¨ x + kx = f ext + SP a , ⌃ ⌃ ⌃ ⌃ N ⌃ � � ⌃ ⌃ ⌥ ⌃ S ˙ = u · n = � x u · n . ⌃ ⌃ ⌅ Γ i Γ 0 i =1

Wellposedness Energy Equality N ✓ ρ ◆ ✓Z ◆ 2 Z Z d | u | 2 + m x | 2 + k X 2 | x | 2 | r u | 2 2 | ˙ + + u · n µ R i 2 dt Ω Ω Γ i i =1 | {z } | {z } | {z } Dissipation within Ω Total energy Dissipation in the subtrees N Z X ρ | u | 2 ( u · n ) = � + P 0 S ˙ + f ext ˙ x x . 2 | {z } | {z } Γ i i =0 Power of inlet pressure Power of ext. forces | {z } In/out-come of kinetic energy and we have 8 C k u k L 2 ( Ω ) kr u k 2 L 2 ( Ω ) , if d = 2 , � � Z < � � ( u · r u ) u �  � � C k u k 1 / 2 L 2 ( Ω ) kr u k 5 / 2 � : L 2 ( Ω ) if d = 3 , Ω

Wellposedness • p is replaced by the total pressure p + ρ | u | 2 / 2, existence of weak solutions, [C. G., B. Maury, Y. Maday, 05] • u = λ i U i 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. Ba ffi co, 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 • u = λ i U i on Γ i , ◆ 2 ✓Z Z Z and σ f ( u , p ) · nU i = � P a U i · n � R i U i · n Γ i Γ i Γ i then Z Z X X | u | 2 ( u · n ) ( λ i ) 3 | U i | 2 ( U i · n ) = Γ i Γ i i i i k u k 3  C P X ( Γ i ) By taking X ( Γ i ) = ( H 1 / 2 ( Γ i )) 0 we obtain Z X | u | 2 ( u · n )  C k u k 3 L 2 ( Ω ) Γ i i

Wellposedness General case • Consider a special basis ( w i ) i for Galerkin approximation such that the ( w i ) i are the eigenfunctions of the operator A defined by ✓Z ◆ ✓Z ◆ Z X ( A v , w ) 0 = µ r v : r w + v · n w · n R i , Ω Γ i Γ i i where ✓Z ◆ ✓Z ◆ Z vw + m ( v , w ) 0 = ρ v · n w · n , S 2 Ω Γ 0 Γ 0 for any v , w 2 { v 2 H 1 0 , Γ 0 ( Ω ) , div v = 0 } • Take A u as a test function.

Wellposedness ◆ 2 ✓Z ◆ ✓Z ◆ ✓Z Z ∂ t u · A u + m = µ d Z R i d X | r u | 2 + ∂ t u · n A u · n u · n ρ , S 2 2 dt 2 dt Ω Γ 0 Γ 0 Ω Γ i i ✓Z ◆ ✓Z ◆ Z X = k A u k 2 r u : r A u + u · n A u · n µ R i 0 , Ω Γ i Γ i i Nonlinear term: Z ( u · r ) u A u  k u k L ∞ kr u k L 2 k A u k 0 , Ω Question: Is k u k L ∞ controlled by k A u k 0 ? Answer: Yes, since D ( A ) ⇢ H 3 / 2+ ε for some ε > 0 and there exists θ > 0 such that k u k L ∞  kr u k θ L 2 k A u k 1 − θ . 0

Wellposedness • True under the assumption of right angles at the outlets ; • Relies on regularity result for the solution of the Stokes problem in polyg- onal domains [Maz’ya, Rossman, 07] ; • Nonlinear Gronwall lemma enables to conclud ; • Additional estimate by taking ∂ t u 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 3 D part with the 0 D part?

Artificial boundary Prescribe: ( � p + | u | 2 2 ) n + r u · n Prescribe: � p n + ( r u + ( r u ) T ) · n Prescribe: � p n + r u · n

Numerical Stability Instabilities observed for the Navier-Stokes equations with Neumann boundary conditions = ⇒ Use of stabilization methods;

Semidiscretize scheme ⇢ u n +1 � u n 8 + ⇢ ( u n +1 · r ) u n +1 � ⌘ ∆ u n +1 + r p n +1 = 0 in Ω , > > > ∆ t > > r · u n +1 = 0 > in Ω , > > > > u n +1 = 0 > on Γ ` , > < ⌘ r u n +1 · n � p n +1 n = � p n +1 on Γ in , in n > > > ✓ ◆ Z ⌘ r u n +1 · n � p n +1 n = � u n +1 · n > > R on Γ out , n > > > > Γ out > u 0 = u 0 > > in Ω . : (2.3.30) Existence of strong solution for small enough data (initial data, applied force)

Numerical Treatment of 3D/0D Coupling 8 ρ∂ t u + ρ ( u · r ) u � η ∆ u + r p = 0 in Ω , > > > r · u = 0 in Ω , > > > > > u = 0 on Γ ` , < σ · n = η r u · n � p n = � p in n on Γ in , > > > σ · n = η r u · n � p n = � p i on Γ i out , i = 1 , . . . , N, out n > > > > > u (0 , · ) = u 0 ( · ) in Ω , : with Z ! p i out ( t ) = F i ( Q i ( t )) = F i u ( s, · ) · n , with 0  s  t, Γ i out Blood Lung C R p out Q R p R d Z t P d Q p out p out = RQ + 1 Q C P p p = 0 0 C Figure 2.2: The RC reduced model. Z t Z t F i ( Q i ) = R i Q i + 1 u · n + 1 Z Z Q i = R i u · n , C i C i Γ i Γ i 0 0 Z t out out d , 0 e − t/ τ i + 1 F i ( Q i ( t )) = R i p Q i ( t ) + P i Q i ( s ) e ( s − t ) / τ i ds j and j if the 3D/0D-interface i i is localized at the genera C i 0 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( R 2 T ν ). Moreover if or ∆ t ≤ K ρ R > K ν , R we obtain a better upper bound. – RC model: ρ C ∆ t ≤ K RC + T – RCR model: same type of su ffi cient conditions than for the R model.

Parameters identification Question Are we able to recover the resistances, sti ff ness parameters thanks to measurements of exhaled volume and flow rate at the mouth ? Partial answer [ Ph.D. thesis of A-C. Eglo ff e] • 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 + r u = in Ω f ,   div u = 0 in Ω  ν ∂ u ∂ n � p n = on Γ e g ,   ν ∂ u ∂ n � p n + k u = on Γ 0 0 ,  We would like to estimate k k 1 � k 2 k L 2 ( Γ 0 ) with respect to the di ff erence u 1 � u 2 and p 1 � p 2 on some part of the boundary. We are going to estimate k k 1 � k 2 k L 2 ( K ) , where K is a compact set of Γ 0 such that u 1 is not equal to zero on K . We have that k k 1 � k 2 k L 2 ( K )  C ( k u 1 � u 2 k 1 / 2 H 1 ( Ω ) + k p 1 � p 2 k 1 / 2 H 1 ( Ω ) ) Related work: [Phung, COCV]