\documentclass[XUPS,XML,SOM,Unicode,francais, NoFloatCountersInSection, NoEqCountersInSection]{cedram}
\usepackage{xups22-01}
\setcounter{tocdepth}{2}
%\XUPScorrections

\hyphenation{modè-le modè-les liné-ai-re liné-ai-res dépen-dre résul-tat résul-tats coef-fi-cient coef-fi-cients notée noté notés appli-qué appli-quée appli-qués appli-quées milieu homo-gé-né-isé patients varia-bi-li-té}

\begin{document}
\frontmatter

\title[Modélisation de la ventilation]{Modélisation mathématique et\nobreakspace numérique de la ventilation}

\author[\initial{C.} \lastname{Grandmont}]{\firstname{Céline} \lastname{Grandmont}}
\address{Inria, 2 Rue Simone Iff, 75012 Paris, France\\
\& LJLL, Sorbonne Université, CNRS,
4 Place Jussieu, 75005 Paris, France}
\email{celine.grandmont@inria.fr}
\urladdr{https://team.inria.fr/commedia/grandmont/}

\thanks{Journées X-UPS 2022. Mathématiques et biologie. Éditions de l'École polytechnique, 2022}

\begin{abstract}
Le poumon est un organe complexe, à la fois multi-physique et multi-échelle, dont les fonctions principales sont le transport de l'air de la bouche aux alvéoles pulmonaires et les échanges gazeux entre les gaz atmosphériques et les gaz du sang. Dans ce texte nous nous intéressons plus particulièrement à la ventilation et nous présentons une hiérarchie de modèles mathématiques permettant de la décrire et de la simuler. Les différents modèles peuvent être utilisés pour comprendre les phénomènes physiologiques en jeu, explorer différents scénarios ou encore développer des outils d'aide au diagnostique.
\end{abstract}

\maketitle
\enlargethispage{5pt}
\noindent
\begin{minipage}{\textwidth}
\tableofcontents
\end{minipage}

\mainmatter

\Changel

Dans ce texte, nous souhaitons illustrer comment la \hbox{modélisation} mathématique, l'étude mathématique et numérique des modèles développés et leurs simulations numériques peuvent apporter une meilleure compréhension de certains phénomènes physiologiques ou biologiques. Outre une meilleure compréhension, les modèles peuvent également aider à explorer différents scénarios sans avoir nécessairement recours à des expériences \emph{in vivo} ou \emph{in vitro}, permettre de prédire des comportements, ou même aider au diagnostic ou à l'optimisation de protocoles thérapeutiques. Ces modèles participent donc à éclairer un contexte physiologique ou médical particulier et représentent un outil supplémentaire qui complète les expériences ou les observations. Ils soulèvent, par ailleurs, bien sûr, des questions mathématiques et numériques propres, qui vont nécessiter des développements méthodologiques spécifiques. Ces développements mathématiques sont indispensables afin de s'assurer, par exemple, de l'existence et de l'unicité des solutions, explorer leur comportement qualitatif, s'assurer des bonnes propriétés des schémas numériques mis en \oe uvre.

\begin{figure}[htb]
\centering
\includegraphics[scale=.8]{poumon_These_Nicolas_Pozin}
\caption{Maillage du parenchyme pulmonaire et arbre bronchique -- Thèse de Nicolas Pozin \cite{Pozin-2017}}
\end{figure}

Ainsi nous présenterons ici une hiérarchie de modèles déterministes constitués d'équations différentielles ordinaires ou d'équations aux dérivées partielles permettant de décrire la respiration et plus spécifiquement la ventilation.
La modélisation dépend bien sûr des besoins ou questionnements qui la motivent et le cas échéant, des données disponibles: particulièrement dans le contexte médical où certains paramètres sont parfois difficilement accessibles ou mesurables, et cette limitation est à prendre en compte dans la démarche de modé\-lisation. Le but ici n'est pas de donner un tableau exhaustif, ni des preuves complètes et rigoureuses mais plutôt d'esquisser la démarche de modélisation, les questions mathématiques et numériques qui en découlent et la possible utilisation des modèles pour éclairer la physiologie.

\section{Contexte physiologique}
La fonction principale de l'appareil respiratoire est d'assurer les échanges gazeux entre l'air ambiant et le sang. L'air inspiré, riche en oxygène, est transporté dans les voies aériennes: pharynx, \hbox{larynx}, trachée, bronches, bronchioles jusqu'aux alvéoles pulmonaires où ont lieu la diffusion des gaz entre l'air ambiant et le sang à travers la membrane alvéolo-capillaire. L'air est donc transporté à travers l'arbre bronchique qui, chez l'humain, est un arbre dyadique à environ 23 générations, dont les échelles sont très différentes puisque le diamètre de la trachée est d'environ deux centimètres tandis que celui des bronchioles terminales est de l'ordre du millimètre. La première partie de l'arbre bronchique, jusque vers la quinzième génération, est purement conductrice et a pour fonction
principale le transport de l'air. Au-delà, des alvéoles viennent se greffer sur les bronchioles et permettent un échange de gaz (oxygène et dioxyde de carbone) avec les capillaires sanguins qui tapissent l'extérieur de quelque 300 millions d'alvéoles. Le mouvement de l'air dans ce réseau de tubes est assuré par la contraction du diaphragme et le déplacement de la cage thoracique, véritables moteurs de la respiration. Cela entraîne, à l'inspiration, l'augmentation du volume du poumon et induit une différence de pression entre la bouche (ou le nez) et les alvéoles pulmonaires, créant alors un flux d'air dans l'arbre. Au repos en respiration normale, l'expiration est quant à elle passive et l'air est expulsé grâce à l'élasticité du tissu pulmonaire appelé \emph{parenchyme pulmonaire}. La~fonction du poumon est donc à la fois de transporter l'air et d'assurer les échanges gazeux. Nous allons ici nous concentrer sur le phénomène de transport en nous attachant à décrire le flux d'air dans les bronches et l'action des tissus élastiques pulmonaires. Nous allons présenter plusieurs modèles décrivant l'irrigation du parenchyme pulmonaire par l'arbre bronchique partant d'un modèle simple décrivant la dynamique du volume d'air inspiré pour aboutir à des modèles tridimensionnels décrivant soit le comportement élastique des tissus pulmonaires et leur interaction avec l'arbre bronchique soit l'écoulement de l'air dans la partie proximale de l'arbre en prenant en compte la dynamique du parenchyme. Il s'agit d'une démarche \emph{bottom-up} où l'on enrichit les modèles afin de prendre en compte des phénomènes physiques négligés, intégrer des non-linéarités, la dimension spatiale... Nous essayerons d'illustrer comment chacun de ces modèles peuvent apporter une meilleure compréhension de la ventilation dans des cas sains ou pathologiques et quelles sont leurs limitations et leurs propriétés mathématiques.

\section{Un modèle d'équation différentielle ordinaire}
Dans cette section, nous allons adopter un point de vue macroscopique qui prend en compte des données accessibles en routine clinique, par exemple lorsque des tests de spirométrie sont effectués afin de diagnostiquer des pathologies telles l'asthme ou en service de réanimation lorsque le patient ou la patiente sont placés sous respirateur mécanique. Dans ce cas, on a accès à la dynamique des volumes expirés et même, dans le cas de la ventilation mécanique, à la pression à la bouche. Ainsi on cherche un modèle qui permet de décrire la dynamique du volume pulmonaire, noté $V(t)$ à l'instant $t$. Le poumon humain peut être vu comme un ballon élastique caractérisé par son élastance, ici notée $E$, et le réseau de tubes constituant l'arbre bronchique peut-être décrit par un unique tube caractérisé par sa résistance, ici notée $R$.

\begin{figure}[ht]
\centering
\includegraphics[height=5cm]{0D2}
\caption{Représentation schématique du poumon}
\end{figure}

Plus précisément, l'air traverse un tube et son débit $Q(t)$
est proportionnel à la différence de pression à l'entrée, supposée être la pression atmosphérique $P_{\atm}$, et à la sortie de celui-ci, pression qui correspond à la pression dans les alvéoles pulmonaires que l'on suppose uniforme et notée $P_{\alv}$. On obtient
\begin{equation}
\label{eq_tube}
P_{\atm}-P_{\alv} = R Q(t) = R V'(t).
\end{equation}

Ce terme de résistance hydraulique vient des frottements du fluide sur les parois et, en particulier, si l'on suppose que l'écoulement est un écoulement laminaire de type Poiseuille la résistance $R$ d'un cylindre de diamètre $d$ et de longueur $l$ s'écrit
\begin{equation}
\label{R_Poiseuille}
R=\frac{128\nu}{\pi}\frac{l}{d^4},
\end{equation}
où $\nu$ est la viscosité du fluide. Nous verrons dans la section suivante comment calculer une résistance équivalente de l'arbre bronchique.

Le comportement élastique du ballon est exprimé par une relation linéaire qui lie le saut de pression et la variation de volume par rapport au volume à vide $V_0$, qui correspond au volume du poumon en fin d'inspiration dans les cas sains, soit
\begin{equation}
\label{eq_ballon}
P_{\alv}-P(t) =E(V(t)-V_0),
\end{equation}
où $P(t)$ est la force appliquée sur le poumon venant de la déformation du diaphragme.
Ainsi on obtient l'équation différentielle linéaire du premier ordre:
\begin{equation}
\label{eq_RC}
R V'(t)+ E(V(t)-V_0)= P_{\atm}-P(t).
\end{equation}

On décrit donc la ventilation pulmonaire par un simple système $RC$, où $C=\sfrac{1}{E}$ est la capacité/compliance du système.

Si on suppose que la donnée initiale est $V_\ri$ et on suppose, pour simplifier, que $P$ est continue en temps, le problème de Cauchy admet une unique solution donnée par:\vspace*{-3pt}\enlargethispage{\baselineskip}
\begin{equation}
\label{sol_RC}
V(t)= V_0 +(V_\ri-V_0) e^{-\sfrac{t}{RC}}+\int_0^t e^{-\sfrac{(t-s)}{RC}}(P_{\atm}-P(s)) \text{d}s.
\end{equation}

En particulier le temps caractéristique de relaxation de ce simple système est le produit $RC$.
L'équation \eqref{eq_RC} est celle utilisée par les médecins et lorsque les volumes, débits, et pressions peuvent être mesurés, les valeurs $R$ de la résistance des voies aériennes et $E$ de l'élastance des tissus sont estimées (en utilisant, par exemple, en routine clinique et pour des patients et des patientes sous ventilateur mécanique, certains instants du cycle respiratoire). Ce sont en effet des indicateurs de l'état du ou de la malade: les asthmatiques auront, entre autre, une résistance des voies aériennes augmentée car l'asthme réduit la lumière bronchique et donc le diamètre $d$ des bronches (voir l'expression de la résistance de Poiseuille \eqref{R_Poiseuille} qui est inversement proportionnelle à $d^4$). Dans ce cas le temps de relaxation est en particulier augmenté et c'est également le cas lorsque $C$ augmente (en cas d'emphysème qui diminue l'élasticité des tissus). Ainsi si l'on admet cette EDO comme premier modèle elle permet d'avoir une indication de l'état du patient et éventuellement de l'évolution de sa maladie.

D'un point de vue mathématique, sur cette simple équation différentielle ordinaire linéaire, on peut déjà se demander s'il existe une solution périodique quand la pression exercée est périodique, et lorsqu'on considère cette solution périodique, quelles sont les conditions sur les paramètres et la période pour que le volume initial associé à cette solution périodique soit proche du volume à vide $V_0$. On peut également tester ces résultats avec des paramètres physiologiquement pertinents ou leurs limites quand $R$ ou $E$ tendent vers l'infini. Par exemple si $P$ est périodique de période $T$, si on considère $P_{\atm}=0$ pour simplifier, il existe une unique solution périodique et celle-ci est associée à la condition initiale\vspace*{-3pt}
\[
V_\ri=V_0 -\frac{1}{R(1-e^{\sfrac{-T}{RC}})}\int_0^T e^{-\sfrac{(T-s)}{RC}} P(s) \text{d}s.
\]
Ainsi le volume à la fin d'un cycle respiratoire (ou en début d'inspiration) n'est pas exactement le volume à vide du poumon. Le volume additionnel 
\[
-\frac{1}{R(1-e^{\sfrac{-T}{RC}})}\int_0^T e^{-\sfrac{(T-s)}{RC}} P(s) \text{d}s
\]
est négligeable pour un cas sain.
En revanche, si la résistance des voies aériennes tend vers $+\infty$ et que les efforts sont maintenus constants, $V_\ri$ tend vers $V_0-\Psfrac{C}{T}\int_0^T P(s) \text{d}s$. Ainsi, afin de vider le poumon, il~est alors nécessaire d'avoir une expiration active.

Cette simple équation différentielle ordinaire est bien sur un modèle très simplifié de la réalité où on a négligé l'inertie du poumon (qu'il peut être important de prendre en compte pour des patients obèses), où on a supposé que le comportement des tissus élastiques était linéaire, ce qui n'est pas le cas. Des modèles plus riches peuvent être envisagés incluant la prise en compte de non-linéarités en faisant dépendre par exemple l'élastance de l'état d'inflation du poumon et donc de $V$ ou encore en prenant en compte les effets inertiels dans l'écoulement dans l'arbre bronchique et en faisant dépendre la résistance du débit, soit $V'$.

Nous ne détaillerons pas les extensions possibles et renvoyons au livre de Bertrand Maury \cite{Livre-Bertrand} qui dresse un panorama très complet de modèles possibles et de questions mathématiques qui en décou\-lent, et propose également des exercices accessibles et variés en lien avec la physiologie. Il est à noter que si la résistance et la compliance sont des paramètres estimés en routine clinique grâce au monitorage de la pression et des volumes inspirés et expirés, certains paramètres auxquels les médecins voudraient pourtant avoir accès comme le volume \og à vide\fg du poumon nécessitent d'autres mesures et protocoles médicaux. Enrichir les modèles en faisant intervenir d'autres paramètres ne fait que renforcer cette problématique.

\section{Un modèle jouet mono-dimensionnel}
Dans cette section nous allons présenter un modèle jouet où le parenchyme pulmonaire est décrit par un système de masses-ressorts, représentant le parenchyme pulmonaire, système lui-même \og irrigué\fg par un arbre dyadique représentant l'arbre bronchique dans lequel l'air circule. La description et l'obtention de ce modèle sont issus de~\cite{Grandmont2006}, dans lequel sont détaillées les preuves des différents résultats énoncés. Le semi-groupe associé et le comportement en temps long des solutions a par ailleurs été également analysé dans la thèse de Christine Vannier \cite{these-christine}.

\subsection{L'arbre bronchique comme réseau de résistances}
Dans le modèle précédent, l'arbre bronchique était décrit par un unique scalaire représentant la résistance globale de l'arbre bronchique. Afin de prendre en compte la structure arborescente de l'arbre bronchique, nous allons maintenant assimiler celui-ci à un réseau de résistances. On considère un arbre dyadique à $N$ générations, chacune des branches de l'arbre étant associée à une résistance. La génération $n=0$ correspond à la trachée qui se sépare en deux bronches et ainsi de suite, jusqu'aux branches terminales. Il y a donc $2^N$ sorties (on~rappelle que l'arbre bronchique humain possède environ 23 générations soit plus de huit millions de sorties).

Dans chaque tube, on suppose que la loi de Poiseuille est vérifiée et donc que la différence de pression est proportionnelle au débit, le coefficient de proportionnalité étant la résistance du tube. Par exemple dans le réseau suivant:

\begin{figure}[ht]
\centering
\includegraphics[scale=.25]{Arbre}
\vspace*{-.7\baselineskip}
\caption{Arbre à une génération}
\vspace*{-.7\baselineskip}
\end{figure}
\noindent on a les relations suivantes:\vspace*{-3pt}
\[
p_0-p= r_0 q_0, \quad p-p_1=r_1 q_1, \quad p-p_2=r_2 q_2,
\]
où les $q_i$ désignent les débits dans chacune des branches. De plus, la~conservation de la masse impose:\vspace*{-3pt}
\[
q_0=q_1+q_2
\]
Si maintenant on suppose que $r_1=r_2$ et que $p_1=p_2$ alors\vspace*{-3pt}
\[
p_0-p_1=(r_0+\sfrac{r_1}{2}) q_0,
\]
et
$r_0+\sfrac{r_1}{2}$ est la résistance équivalente du réseau.
Nous allons supposer que les résistances de l'arbre bronchique que l'on considère sont constantes par génération et la résistance à la génération $n$ est notée~$r_n$,~$0\leq n\leq N$.
Ainsi des calculs analogues permettent de calculer sa résistance équivalente soit
\begin{equation}
\label{res_eq}
R^N_{\eq}=\sum_{n=0}^N \frac{r_n}{2^n}.
\end{equation}
Si $r_n$ est une résistance de Poiseuille, alors
\[
r_n=\frac{128\nu}{\pi}\,\frac{l_n}{d_n^4},
\]
où $l_n$ est la longueur des branches à la génération $n$ et $d_n$ leur diamètre. En supposant que l'arbre est géométrique c'est-à-dire qu'il existe $\lambda$ tel que $l_n=\lambda^n l_0$ et $d_n=\lambda^n d_0$ alors
\begin{equation}
\label{res_eq_geo}
R^N_{\eq}=r_0\sum_{n=0}^N \frac{1}{(2\lambda^3)^n}.
\end{equation}
Cette série converge dès que $\lambda > \sfrac{1}{\sqrt[3]{2}}$ or les mesures physiologiques montrent que pour le poumon humain ce coefficient de réduction est plutôt, au moins pour la zone de conduction, de l'ordre de $0.85$ et ainsi appartient à la zone de convergence de la série géométrique et $ R^\infty_{\eq}<+\infty$.
\begin{remark}
Si maintenant on s'intéresse au volume occupé par l'arbre bronchique celui-ci s'écrit comme la somme de termes du type $(2\lambda^3)^n$ et donc la série associée est, elle, divergente.
\end{remark}
Dans \cite{mauroy} les auteurs se sont intéressés à l'optimalité de l'arbre bronchique. En reprenant l'expression \eqref{res_eq} et en supposant que $d_n=c\, l_n$, on peut se poser la question sous la forme d'un problème de minimisation sous contrainte de volume, qui s'écrit
\begin{equation}
\label{pb_min}
\min_{\substack{l=(l_0, l_1, \dots, l_N)\in \mathbb{R}^{N+1}\\\sum_{n=0}^N {2^n l_n^3}\leq V_{\max}}}\ \sum_{n=0}^N \frac{1}{2^n l_n^3}
\end{equation}
Il s'agit donc de minimiser la résistance (ou de manière équivalente la puissance dissipée dans l'arbre) sous une contrainte de volume.
On montre facilement que le minimiseur, s'il existe, appartient à un compact $[m, M]^{N+1}$ de $\mathbb{R}^{N+1}$, avec $m>0$. On en déduit l'existence du minimiseur. Celui-ci peut-être calculé en écrivant l'équation d'Euler-Lagrange de ce problème de minimisation sous contrainte d'inégalité (contrainte saturée à l'optimum) et on obtient alors que\vspace*{-3pt}
\[
l_0=\sqrt[3]{\spfrac{V_{\max}}{N+1}},\quad l_n=l_0\bigl(\sfrac{1}{\sqrt[3]{2}}\bigr)^n.
\]
De ce point de vue, le poumon humain est donc sous-optimal. Cependant $\sfrac{1}{\sqrt[3]{2}}$ est la valeur critique pour laquelle à la fois la série qui définit la résistance équivalente et celle qui définit le volume des bronches sont divergentes. Ainsi cette sous-optimalité, en terme de résistance, assure en fait une certaine robustesse vis-à-vis de variations des diamètres des bronches par exemple.

Nous allons maintenant exprimer ce réseau de résistances comme un opérateur liant les flux en bout d'arbre (et donc irriguant le parenchyme pulmonaire) et les pressions également en bout d'arbre (et~donc s'appliquant sur le parenchyme pulmonaire). On indexe les sorties par $0, 1, \dots, 2^N-1$. On appelle $(q^n_{k})_{k\in \{0,\dots, 2^n-1\}}$ les débits dans les branches de la génération $n\geq 1$, les pressions sont notées $(p^n_{k})_{k\in \{0,\dots, 2^n-1\}}$. On va exprimer le saut de pression de la racine à la sortie $0$. La pression à l'entrée de la première bifurcation est notée $p_0^0$ et le débit dans la trachée $q_0^0$.
Pour la branche mère, on a\vspace*{-3pt}\enlargethispage{.5\baselineskip}
\[
p_{\atm}-p_0^0= r_0 q_0^0.
\]
Le débit $q_0^0$ est le débit total et donc\vspace*{-3pt}
\[
q_0^0=\sum_{k=0}^{2^N-1} q^N_k.
\]
A l'étape suivante, on a de même\vspace*{-8pt}
\[
p_0^0-p_0^1=r_1 q_0^1, \quad\text{avec}\quad q_0^1=\sum_{k=0}^{2^{N-1}-1} q^N_k
\]
On obtient finalement\vspace*{-8pt}
\[
p_{\atm}-p_0^N= \sum_{n=0}^N r_n \sum_{k=0}^{2^{N-n}-1} q^N_k.
\]
On peut refaire le même raisonnement pour chacune des sorties de l'arbre et en considérant le vecteur des pressions de sorties $p=(p^N_{k})_{k\in \{0,\dots, 2^N-1\}}$, et celui des débits à la dernière génération $q=(q^N_{k})_{k\in \{0,\dots, 2^N-1\}}$, il vient
\[
p_{\atm} e- p= A_N q,
\]
où $e$ le vecteur dont toutes les composantes sont égales à $1$ et $A_N$ la matrice de résistances.
On choisit dans la suite $p_{\atm}=0$.

En introduisant $J_n$ matrice $2^n\times 2^n$ dont tous les éléments sont égaux à $1$, la matrice $A_N$ est donnée par
\begin{multline*}
A_N=r_0 J_N+\begin{pmatrix}r_1 J_{N-1} & 0\cr 0 & r_1 J_{N-1} \end{pmatrix}\\
+\begin{pmatrix}r_2 J_{N-2} & 0 & 0 & 0\cr 0 & r_2 J_{N-2} & 0 & 0\cr 0 & 0 &r_2 J_{N-2} & 0\cr 0 &0 & 0 &r_2 J_{N-2} \end{pmatrix} +\dots + r_n I_{2^N}.
\end{multline*}

Le but ici est de regarder la limite quand le nombre $N$ de géné\-ra\-tions de l'arbre tend vers $+\infty$.
Pour cela, nous allons maintenant écrire la relation entre les pressions en sorties et les débits en sorties à l'aide d'un opérateur sur l'espace $L^2(0,1)$. On considère l'intervalle $I={}]0, 1[$ et on associe à chaque sortie $i$ de l'arbre le sous intervalle \hbox{$]i/2^N, (i+1)/2^N[$}.
A toute fonction $u\in L^2(0, 1)$, on associe le vecteur de flux $q^N\in \mathbb{R}^{2^N}$ défini par
\[
q^N= (q^N_i)_{0\leq i\leq 2^N-1},\quad\text{avec}\quad q^N_i=\int_{i/2^N}^{(i+1)/2^N} u(s) \rd  s.
\]
On définit également les pressions associées à travers l'arbre, soit
\[
p^N= (p^N_i)_{0\leq i\leq 2^N-1}=- A_N q^N,
\]
et la fonction $p\in L^2(0, 1)$ constante par morceaux sur chacun des sous intervalles prenant la valeur $p_i^N$ sur $]i/2^N, (i+1)/2^N[$. Cela définit un opérateur $\mathcal{R}_N$ de $L^2(0, 1)$ dans lui-même, ce qu'on écrit $\mathcal{R}_N\in\mathcal{L}(L^2(0, 1))$.
\begin{propo}
\label{prop-noyau}
L'opérateur $\mathcal{R}_N$ s'exprime comme un opérateur à noyau:\vspace*{-10pt}
\[
p=-\mathcal{R}_N u\ \iff\ p(x)=-\int_0^1 K_N(x, y) u(y) \rd  y,
\]
où le noyau $K_N\in L^1(I\times I)$ est constant par morceaux et de la forme\vspace*{-3pt}
\[
K_N(x, y) = \sum_{n=0}^N r_n \sum_{i=0}^{2^n-1} \chi^n_i (x,y),
\]
où $ \chi^n_i (x,y)$ est la fonction caractéristique du pavé\vspace*{-3pt}
\[
]i/2^n, (i+1)/2^n[{}\times{}]i/2^n, (i+1)/2^n[.
\]
\end{propo}
On remarque que la norme de $K_N$ dans $L^1(I\times I)$ est précisément la résistance équivalente de l'arbre \eqref{res_eq}.
On introduit maintenant le noyau limite\vspace*{-8pt}
\begin{equation}
\label{noyau_limite}
K(x, y) =
\sum_{n=0}^{+\infty} r_n \sum_{i=0}^{2^n-1} \chi^n_i (x,y),
\end{equation}
Cette somme est en fait toujours une somme finie à $(x,y)$ fixé.
On a le théorème suivant:

\begin{thrm} La suite d'opérateurs $(\mathcal{R}_N)$ converge vers $\mathcal{R}$ dans $\mathcal{L}(L^2(0, 1))$ si et seulement si\vspace*{-3pt}
\begin{equation}
\label{rest_cv}
\sum_{n=0}^{+\infty}\frac{r_n}{2^n}<+\infty.
\end{equation}
Quand l'opérateur $\mathcal{R}$ existe, il est associé au noyau $K$ défini par~\eqref{noyau_limite}, qui est alors dans $L^1(I\times I)$.
De plus, si \eqref{rest_cv} est vérifiée, $\mathcal{R}$ est auto\-adjoint, compact et monotone (c'est-à-dire $(\mathcal{R}u, u)_{L^2(0, 1)}\geq 0$ pour $ u\in L^2(0, 1)$).

\end{thrm}
En particulier la compacité de $\mathcal{R}$ vient du fait qu'il est limite d'une suite d'opérateurs de rang fini. De plus, en utilisant \eqref{noyau_limite}, il vient facilement\vspace*{-5pt}
\begin{equation}
\label{Monotone_R}
(\mathcal{R}u, u)_{L^2(0, 1)}=\sum_{n=0}^{+\infty} r_n \sum_{i=0}^{2^n-1}\biggl(\int_{i/2^n}^{(i+1)/2^n} u(x)\rd  x\biggr)^2\geq 0.
\end{equation}
Cette quantité représente l'énergie dissipée dans l'arbre.\medskip

Replaçons nous dans le contexte du poumon humain où, comme nous l'avons vu, la condition \eqref{rest_cv} est satisfaite. Nous retrouvons ici la robustesse du poumon humain vis-à-vis de perturbations et en particulier vis-à-vis du nombre de générations dans l'arbre. Lorsque
$
\sum_{n=0}^{+\infty}\sfrac{r_n}{2^n}=+\infty,
$
alors le noyau $K_N$ se concentre sur la diagonale et pour tout $q\in L^2(0, 1)$,\vspace*{-4pt}
\[
\int_0^1 \frac{K_N(x, \cdot)}{\|K_N\|_{L^1(I\times I)} }\,q(x) \rd  x
\]
converge fortement vers $q(\cdot)$ dans $L^2(0, 1)$.

Sous l'hypothèse \eqref{rest_cv}, la décomposition spectrale de l'opérateur $\mathcal{R}$ est de plus immédiate. En effet, introduisons la base de Haar, qui est une base hilbertienne de $L^2(0, 1)$. Soit $\psi\in L^2(0, 1)$ définie par \hbox{$\psi=1_{]0, 1/2[}- 1_{]1/2, 1[}$}. Pour tout $n\geq 0$ et $0\leq k\leq 2^n-1$, on pose $\psi_{n, k}(x)=2^{n/2}\psi(2^n x-k)$, on note également $\psi_0$ la fonction identiquement égale à un sur $(0, 1)$. La famille $(\psi_0, (\psi_{n, k})_{n\geq 0,\, 0\leq k\leq 2^n-1})$ est alors la base de Haar et un calcul direct montre que
\begin{propo}
La famille $(\psi_0, (\psi_{n, k})_{n\geq 0,\, 0\leq k\leq 2^n-1})$ est la famille de vecteurs propres de $\mathcal{R}$ associés
aux valeurs propres\vspace*{-5pt}
\[
\lambda_0=\sum_{l=0}^{+\infty}\frac{r_l}{2^l}, \quad \lambda_{n, k}=\rho_{n+1}=\sum_{l=n+1}^{+\infty}\frac{r_l}{2^l}, \quad \forall n\geq 0, \ 0\leq k\leq 2^n-1.
\]
\end{propo}

Si l'arbre n'est pas symétrique (\ie si les résistances ne sont pas constantes à chaque génération) alors les vecteurs propres et valeurs propres ne sont plus si faciles à calculer. Nous renvoyons notamment aux travaux de Bertrand Maury, dont un résumé se trouve dans le livre \cite{Livre-Bertrand} mais également à \cite{seb-bertrand-res} pour les questions de sensibilité de l'arbre par rapport aux variations des résistances.\medskip

Dans le cas où l'arbre est géométrique, c'est-à-dire où les résistances s'écrivent $r_n = \sfrac{r_0}{(\lambda^3)^n}$, on peut déduire des propriétés supplémentaires sur $\mathcal{R}$. On pose $\alpha=\lambda^{-3}$ et on se place dans le cas convergent où $\alpha<2$. Les valeurs propres de $\mathcal R$ sont donc $\lambda_0=\spfrac{r_0}{1-\alpha/2}$ et $\rho_{n+1}=\lambda_0 (\sfrac{\alpha}{2})^{n+1}$. On peut alors exprimer le taux de convergence de $\mathcal R_N$ vers $\mathcal R$ en remarquant que, pour tout $u$ dans $L^2(0,1)$\vspace*{-3pt}
\begin{align*}
\|\mathcal R u- \mathcal R_N u\|_{L^2(0,1)}&\leq \|K-K_N\|_{L^1(I\times I)} \|u \|_{L^2(0, 1)}\\
&\leq \rho_{N+1} \|u \|_{L^2(0, 1)}= \lambda_0 \left(\sfrac{\alpha}{2}\right)^{N+1} \|u \|_{L^2(0, 1)}.
\end{align*}
Lorsque $\alpha\in{}]\sqrt{2}, 2[$, on peut également montrer que $\mathcal R$ envoie $L^2(0, 1)$ dans l'espace de Sobolev\footnote{Pour $s\in{}]0,1[$, l'espace de Sobolev $H^s(0,1)$ est, par définition, l'espace des fonctions $u\in L^2(0,1)$ telles que la fonction
\[
(x,y)\mto \dfrac{u(x)-u(y)}{|x-y|^{s+\sfrac{1}{2}}}
\]
appartient à $L^2((0,1)\times (0,1))$.} $H^s(0, 1)$ avec $s=1-\sfrac{\ln(\alpha)}{\ln(2)}$.

Cette propriété vient du fait que $H^s(0, 1)$, lorsque $0\leq s< \sfrac{1}{2}$, est entièrement caractérisé par les coefficients de la décomposition de tout vecteur de $L^2(0, 1)$ dans la base orthonormée de Haar. Plus précisément
\[
H^s(0, 1) =\biggl\{u\in L^2(0, 1)\Bigm| |c_0(u)|^2+\sum_{n=0}^{+\infty}\sum_{k=0}^{2^n-1} 2^{2sn}|c_{nk}(u)|^2<+\infty\biggr\},
\]
avec
\[
c_0(u)=\int_0^1 u(x){\rd }x, \quad c_{nk}(u)=\int_0^1 u(x)\psi_{n, k}(x){\rd }x.
\]
Ainsi pour tout $u\in L^2(0, 1)$
\[
\mathcal R u=\lambda_0 c_0(u)\psi_0 + \lambda_0 \sum_{n=0}^{+\infty}\sum_{k=0}^{2^n-1} \left(\sfrac{\alpha}{2}\right)^{n+1}c_{nk}(u)\psi_{n, k},
\]
qui appartient à $H^s(0, 1)$ pour $s~=~1-\sfrac{\ln(\alpha)}{\ln(2)}$.

\subsection{Un système couplé arbre - masses - ressorts}
Nous allons, dans cette sous-section, connecter cet arbre résistif à un réseau de masses-ressorts représentant le parenchyme pulmonaire.
On considère que l'intervalle $I$ est composé d'une succession de masses et de poches d'air (les alvéoles pulmonaires), reliées les unes aux autres par un arbre dyadique à $N$ générations. On définit les $2^N$ cellules par
\[
C_i^N = ]i/2^N, (i+1)/2^N[, \quad I=\cup_{i=0}^{2^N-1} C_i^N.
\]
De plus, on pose $x_i^N=\sfrac{i}{2^N}$, $h_N=\sfrac{1}{2^N}$, et on définit une partie solide $S_i^N$ et une partie fluide $T_i^N$ en posant
\begin{align*}
S_i^N&={} \bigl]x_i^N-\Psfrac{1}{2}h_N\alpha_S, x_i^N+\Psfrac{1}{2}h_N\alpha_S\bigr[, \quad i\in \{1, \dots 2^N-1\},
\\
S_0^N&= {}\bigl]x_0^N, x_0^N+\Psfrac{1}{2}h_N\alpha_S\bigr[, \quad S_{2^N}^N= {}\bigl]x_{2^N}^N-\Psfrac{1}{2}h_N\alpha_S, x_{2^N}^N\bigr[,
\\
T_i^N&={} \bigl]x_i^N+\Psfrac{1}{2}h_N\alpha_S, x_{i+1}^N-\Psfrac{1}{2}h_N\alpha_S[, \quad i\in \{0, \dots 2^N-1\},
\end{align*}
où $\alpha_S$ est la proportion de solide. Chaque masse est connectée à ses voisines par un ressort de raideur $k_N$. De plus on suppose que la densité totale de solide est constante et indépendante de $N$ ainsi la masse $m_N$ dans chaque solide $S_i^N$ est égale à $\sfrac{\rho\alpha_S}{2^N}$.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=.41]{ressort2}
\end{center}
\caption{Système masse-ressort}
\label{fig:ressort}
\end{figure}

On désigne par $u_i^N$ le déplacement de $x^i_N$ et $f^i_N$ la force extérieure appliquée. Sur chaque masse s'exerce par ailleurs les forces de pressions venant de l'arbre résistif. La force exercée sur la $i$ème masse s'écrit $-p_i^N+p_{i-1}^N$, où $(p_i^N)_{1\leq i\leq 2^n-1}$ est le vecteur des pressions en bout d'arbre. La relation fondamentale de la dynamique s'écrit alors, pour $i\in \{1, \dots, 2^N-1\}$,
\begin{multline}
\label{rfd0}
m_N\ddot{u}_i^N(t)-k_N\bigl({u}_{i+1}^N(t)-2 {u}_i^N(t)+ {u}_{i-1}^N(t)\bigr)\\+ p_i^N(t)-p_{i-1}^N(t)=m_N f^i_N(t),
\end{multline}
où $\dot u$ désigne ici la dérivée par rapport au temps.
On sait d'après la section précédente que $-p^N= A_N q^N$, où $q^N$ est le vecteur des débits en bout d'arbre (compté positivement quand l'air circule de la trachée vers les branches terminales). Ce vecteur peut s'exprimer en fonction des vitesses de chaque masse entourant une alvéole donnée et on a
\[
q_j^N=\dot{u}_{j+1}^N-\dot{u}_{j}^N,\quad i\in \{0, \dots, 2^N-1\}.
\]
On obtient finalement pour $i\in \{1, \dots, 2^N-1\}$
\begin{multline}
\label{rfd}
m_N\ddot{u}_i^N(t)-k_N\bigl({u}_{i+1}^N(t)-2 {u}_i^N(t)+ {u}_{i-1}^N(t)\bigr)\\
-\sum_{j=1}^{2^N} A^N_{i, j-1} (\dot{u}_{j}^N-\dot{u}_{j-1}^N)=m_N f^i_N(t).
\end{multline}
Les deux premiers termes du membre de gauche de \eqref{rfd} représentent respectivement l'accélération de la masse et les forces de rappels des ressorts de part et d'autre de celle-ci. Le dernier terme décrit \hbox{l'action} de l'arbre sur la masse. Ce terme lie toutes les masses entre elles à travers le réseau résistif. Il s'agit d'un terme non local. Ce système est complété par des conditions aux limites et des conditions initiales. Il~s'agit ainsi d'un système couplé d'équations différentielles ordinaires du second ordre pour lequel l'existence d'une unique solution est immédiate.

\subsection{Limite quand \texorpdfstring{$N$ tend vers $+\infty$}{Ninf}}
La question que l'on se pose maintenant est: quel est le modèle limite lorsque le nombre de générations de l'arbre, et donc le nombre de masses, tend vers l'infini?
Afin de passer à la limite quand $N$ tend vers l'infini, un point clé est l'obtention d'estimations \emph{a priori} vérifiées par la solution. La démarche est standard et consiste à multiplier~\eqref{rfd} par $\dot{u}_i^N$ et appliquer le lemme de Gronwall. Cependant nous allons d'abord reformuler le système d'EDO à l'aide d'une équation aux dérivées partielles. Considérons $u$ une fonction continue sur $[0, 1]$, affine dans chaque cellule, pour tout temps, et telle que $u(x_i^N, t)= u_i^N(t)$ où les $u_i^N$ sont solutions de \eqref{rfd}. \medskip

On rappelle que l'espace de Sobolev $H^1(0, 1)$ est l'espace des fonctions $u\in L^2(0,1)$ dont la dérivée au sens des distributions (ou encore dérivée faible) appartient à $ L^2(0,1)$. Muni de la norme
\[
\Vert u\Vert_{H^1(0,1)}=\sqrt{\Vert u\Vert^2_{L^2(0,1)}+\Vert u'\Vert^2_{L^2(0,1)}},
\]
c'est un espace de Hilbert. Il s'avère que tout élément de $H^1(0, 1)$ est une fonction continue sur $[0,1]$ et le sous-espace des $u\in H^1(0,1)$ telles que $u(0)=u(1)=0$ est noté $H^1_0(0, 1)$ (pour plus de détails on pourra se référer à \cite{brezis1999analyse}).

On introduit alors le sous-espace fermé de $H^1_0(0, 1)$ défini par
\begin{multline*}
V_N=\bigl\{v: [0,1]\rightarrow\mathbb R \mid v \in C^0([0, 1]), \\
v \text{ affine sur } C_i^N \text{ pour } i\in \{0, \dots, 2^N-1\} \text{ et } v(0)=v(1)=0 \bigr\},
\end{multline*}
et le sous-espace suivant de $L^2(0, 1)$
\begin{multline*}
M_N=\{\nu: ]0,1[{} \rightarrow \mathbb R \mid
\nu \text{ constante par morceaux: } \\ \exists\nu_i^N \in \mathbb R,\; \nu(x)=\nu_i^N \text{ sur } C_i^N \}.
\end{multline*}
Toutes les dérivées partielles considérées ci-dessous sont des dérivées partielles au sens des distributions.

Soit $v\in V_N$. Pour tout $x\in C_i^N$, on a $\partial_x v(x)=2^N(v(x_{i+1}^N)-v(x_i^N))$. Ainsi $\partial_x v\in M_N$.

De plus, si $\nu\in M_N$, on a
\[
\partial_x\nu=\sum_{i=1}^{2^N-1}(\nu_i^N-\nu_{i-1}^N)\delta_{x_i^N},
\]
où $\delta_{x_i^N}$ désigne la mesure de Dirac en $x_i^N$.

Par conséquence en considérant $u^N$ et $p^N$ telles que, pour $t\in{}]0, T[$,
\begin{itemize}
\item
$u^N(t, \cdot)\in V_N$ avec $u^N(t, x_i^N)=u_i^N(t)$,
\item
et $p^N(t, \cdot)\in M_N$ avec $p^N(t, x)=p_i^N(t)$, pour $x\in C_i^N$,
\end{itemize}
le système d'EDO \eqref{rfd0} s'écrit
\begin{multline*}
m_N \sum_{i=1}^{2^N-1} \partial^2_{t} u(t, x)\delta_{x_i^N}(x) -\frac{k_N}{2^N} \partial^2_{x} u^N(t, x)+\partial_x p^N(t, x)\\[-5pt]
=m_N\sum_{i=1}^{2^N-1} f_i^N(t)\delta_{x_i^N}(x).
\end{multline*}
Mais
\[
q_i^N=\partial_t u^N(x_{i+1}^N)-\partial_t u^N(x_i^N)=\int_{x_i^N}^{x_{i+1}^N}\partial_x \partial_t u^N(t, y)\rd  y,
\]
ce qui entraîne, d'après la Proposition \ref{prop-noyau}, que
\[
p^N(t, x)=-\int_0^1 K_N(x, y) \partial_x \partial_t u^N(t, y)\rd  y
\]
et donc
\begin{multline}
\label{rfd-continu}
m_N \sum_{i=1}^{2^N-1} \partial^2_{t} u^N(t, x_i^N)\delta_{x_i^N}(x) -\frac{k_N}{2^N} \partial^2_{x} u^N(t, x)\\
-\partial_x\left(\int_0^1 K_N(x, y) \partial_x \partial_t u^N(t, y)\rd  y\right)=m_N\sum_{i=1}^{2^N-1} f_i^N(t)\delta_{x_i^N}(x).
\end{multline}
Nous pouvons maintenant démontrer des estimations \emph{a priori}. En supposant que $\sfrac{k_N}{2^N}=k$ et que les données (forces extérieures appliquées et conditions initiales) satisfont des hypothèses de convergence quand $N$ tend vers l'infini, on peut montrer:

\begin{lemma}\label{estim-uN} Il existe une constante $C>0$ telle que pour tout $N$, la~solution $u^N$ vérifie
\[
\|\partial_t u^N(t)\|_{L^2(0, 1)}+\|u^N(t)\|_{H^1(0, 1)}\leq C,
\]
pour presque tout $t\in(0,T)$.
\end{lemma}
La démonstration de ce lemme découle essentiellement de l'égalité d'énergie suivante, obtenue en multipliant \eqref{rfd-continu} par la vitesse $\partial_t u^N$ et en intégrant sur $]0, 1[$, soit
\begin{multline}
\label{egalite-energie-continu}
\underbrace{\frac{m_N}{2}\frac{d}{dt}\sum_{i=1}^{2^N-1} |\partial_t u^N(t, x_i^N)|^2} _{\substack{\text{dérivée en temps}\\\text{de l'énergie cinétique}}}+\underbrace{ \frac{k}{2} \frac{d}{dt}\|\partial_x u^N(t, \cdot)\|^2_{L^2(0, 1)}}_{\substack{\text{dérivée en temps}\\\text{de l'énergie mécanique}}}\\
+\underbrace{(\mathcal R_N \partial_x\partial_t u^N, \partial_x\partial_t u^N)_{L^2(0, 1)}}_{\text{énergie dissipée dans l'arbre}}= m_N\sum_{i=1}^{2^N-1} f_i^N(t)\partial_t u^N(t, x_i^N).
\end{multline}
On remarque que l'énergie dissipée dans l'arbre
\[
(\mathcal R_N \partial_x\partial_t u^N, \partial_x\partial_t u^N)_{L^2(0, 1)}
\]
est positive (\cf \eqref{Monotone_R}). Le membre de droite est majoré en utilisant les propriétés de $u_i^N(t)=u^N(t, x_i^N)$ et l'hypothèse de convergence sur les forces extérieures $f_i^N$, ce qui donne le résultat voulu.

En particulier, l'énergie dissipée dans l'arbre appartient au segment $[0,C]$ pour tout $N$ et presque tout $t$. Ceci permet
d'obtenir un peu de régularité de $ \partial_x\partial_t u^N$ (\ie cette fonction appartient à un certain espace de Sobolev que nous ne préciserons pas ici). De cette régularité, on déduit un résultat de convergence dans le sens suivant.\medskip

Pour $X$ un espace de Hilbert, et $p\in[1,\infty]$, nous introduisons l'espace $L^p(0,T;X)$ des fonctions $w$ de $(0,T)$ dans $X$ telles que $t\mto \Vert w(t)\Vert_X$ appartienne à $L^p(0,T)$. \medskip

Par suite, $L^\infty(0,T;X)$ est le dual de $L^1(0,T;X)$, et on dit qu'une suite de fonctions $(w_n)_n$ de $L^\infty(0,T;X)$ converge $\star$-faiblement vers une fonction $w$ si, pour tout $v\in L^1(0,T;X)$, la suite $(\langle w_n,v\rangle)_n$ converge vers $\langle w,v\rangle$.
\medskip

\noindent Le lemme \ref{estim-uN} donne
\[
\partial_t u^N\in L^\infty(0,T; L^2(0, 1))\quad\text{et}\quad u^N\in L^\infty(0,T; H^1(0,1)).
\]
En supposant que le nombre de générations $N$ tend vers l'infini, que $k=\sfrac{k_N}{2^N}$ et que la résistance globale de l'arbre $\sum_{n=0}^N\sfrac{ r_n}{2^n}$ converge, on déduit alors que la fonction $u^N$, affine par morceaux sur les cellules et valant $u_i^N$ aux n\oe uds $x^N_i$ (et sa dérivée $\partial_t u^N$), converge $\star$-faiblement vers une fonction $u\in H^1_0(0,1)$ (respectivement $\partial_t u$) solution faible (c'est-à-dire, au sens des distributions) de
\begin{equation}
\label{onde-visqueux}
\alpha_S\rho\partial^2_{t} u -k\partial^2_{x} u-\partial_x\mathcal R\partial_x\partial_t u=\rho\alpha_S f,
\end{equation}
où la fonction $f$ est déterminée par l'hypothèse de convergence des forces extérieures $f_i^N$.\medskip

L'équation obtenue est une équation des ondes avec un terme de dissipation non local (en espace). La question de sa stabilisation et du comportement asymptotique des solutions est une question naturelle et cette question pour les équations des ondes avec des termes d'amortissement est très explorée dans la littérature.

La démarche que l'on a suivie ici est: comment, dans un cadre idéalisé, peut-on décrire l'interaction arbre-parenchyme pulmonaire et dans quelle mesure ceci nous apporte néanmoins des éléments sur la physiologie du poumon? Nous avons obtenu un système discret, ainsi que sa limite quand le nombre de générations de l'arbre tend vers l'infini: le modèle obtenu est un modèle visco-élastique, où l'arbre apparaît comme un terme non local qui dissipe l'énergie du système, permettant de montrer qu'en l'absence de forces extérieures, celle-ci tend vers zéro asymptotiquement en temps.

Malgré sa simplicité et son caractère très académique, ce modèle (ou plutôt sa version à un nombre fini de générations) a été utilisé et adapté par Michael Brunengo au cours de sa thèse \cite{These-Brunengo} dans une collaboration industrielle afin d'explorer
et d'optimiser les techniques d'oscillations hautes fréquences de la paroi thoracique.

\section{Généralisation tridimensionnelle}
Nous présentons dans cette section une généralisation du modèle précédent en dimension $3$, les détails des calculs et des démonstrations se trouvent dans \cite{cazeaux-grandmont}. L'écoulement dans l'arbre
bronchique est toujours décrit par un réseau de résistances qui alimente maintenant des alvéoles tridimensionnelles du parenchyme
pulmonaire dont la dynamique est décrite par les équations de l'élastodynamique. On~se propose donc de décrire le parenchyme pulmonaire comme un milieu poreux constitué de cellules perforées réparties périodiquement (les~alvéoles) et connectées par un réseau de résistances dyadiques. Soit~$d$ la dimension, $d=2$ ou $3$.

On fixe une cellule $Y$ (par exemple $]0,1[^d$ ou toute autre cellule permettant le pavage de $\mathbb R^d$), de mesure égale à $1$.

On suppose que $Y$ est constitué de deux ouverts, l'un simplement connexe $Y_F$, tel que $\overline{Y_F} \subset Y$, qui représentera l'alvéole, et $Y_S = Y \setminus \overline{Y_F}$ qui représente la partie élastique. L'interface entre les deux parties est notée $\Gamma = \partial Y_F$.

Pour $\varepsilon>0$ et $k \in \mathbb{Z}^d$, on translate et redimensionne la cellule $Y$ comme suit
\begin{equation}
\begin{aligned}
Y^{\varepsilon,k} &= \varepsilon (Y + k),& Y_F^{\varepsilon,k} &= \varepsilon (Y_F + k),\\
Y^{\varepsilon,k}_S &= \varepsilon (Y_S + k), & \Gamma^{\varepsilon,k} &= \partial Y^{\varepsilon, k}_F.
\end{aligned}
\end{equation}
En utilisant $\mathbb{Z}^d_\varepsilon = \{ k \in \mathbb{Z}^d \mid \overline{Y^{\varepsilon, k}} \subset \Omega \}$ on peut définir les domaines perforés fluide et structure, et leur interface
\begin{equation}
\Omega^\varepsilon_S = \Omega \setminus \bigcup_{k\in \mathbb{Z}^d_\varepsilon} \overline{Y_F^{\varepsilon,k}},\qquad
\Omega^\varepsilon_F =\bigcup_{k\in \mathbb{Z}^d_\varepsilon} Y_F^{\varepsilon,k} ,\qquad
\Gamma^\varepsilon = \bigcup_{k\in \mathbb{Z}^d_\varepsilon} \Gamma^{\varepsilon,k}
\end{equation}

Finalement, on appelle $\chi_F$ et $\chi_S$ les fonctions caractéristiques respectivement de $Y_F$ et $Y_S$, et $\chi_F^\varepsilon$ et $\chi_S^\varepsilon$ les fonctions caractéristiques respectivement de $\Omega^\varepsilon_F$ et $\Omega^\varepsilon_S$ (voir la figure \ref{fig:domaine}).
\begin{figure}[htb]
\begin{center}
\def\JPicScale{.5}
\input{./xups22-01_figures/Domaine.tex}
\caption{\label{fig:domaine}}
\end{center}
\end{figure}

On considère que la partie élastique satisfait aux équations de l'élastodynamique linéarisées (hypothèse de petits déplacements). On~va supposer que
sur une partie de la frontière externe $\Gamma_D$, on applique des conditions de Dirichlet homogènes (c'est-à-dire que l'on suppose que le déplacement est nul sur cette partie du bord), tandis que sur l'autre partie $\Gamma_N $, on applique des conditions de Neumann (c'est-à-dire que l'on suppose que l'on se donne les contraintes normales au bord, ici une force de pression notée $p_{\pleu}$).
\[
% \label{eq:vfdiv0}
\begin{cases}
\rho \partial_t^2 \dd^\varepsilon - \div \ \sigma(\dd^\varepsilon) = \f_S , &\text{dans } \Omega^\varepsilon_S ,\\
\sigma(\dd^\varepsilon) \cdot \nn^\varepsilon = -p_k^\varepsilon\nn^\varepsilon, & \text{sur chaque } \Gamma^{\varepsilon,k},\\
\sigma(\dd^\varepsilon) \cdot \nn =- p_{\pleu}\nn, & \text{sur } \Gamma_N, \\
\dd^\varepsilon = 0, & \text{sur } \Gamma_D,
\end{cases}
\]
où $\nn$ et $\nn^\varepsilon$ désignent les normales unitaires sortantes du domaine et le tenseur des contraintes $\sigma$ est donné par la loi de Hooke
\begin{equation}
\label{hooke}
\sigma(\dd) = \lambda \div (\dd) I + 2 \mu \mathrm{e}(\dd),\quad \mathrm{e}(\dd)=\frac{1}{2}(\nabla \dd+ (\nabla \dd)^T),
\end{equation}
les paramètres $\lambda$ et $\mu$ étant les coefficients de Lamé du matériau caractérisant son comportement élastique. La pression $p_k^\varepsilon$ est la pression du fluide à la sortie de l'arbre correspondant à l'alvéole $Y_F^{\varepsilon,k}$.
On~construit ensuite un arbre et une décomposition du domaine pour lesquelles les sorties de l'arbre sont connectées aux alvéoles. Une telle construction est illustrée par la figure suivante
\begin{figure}[ht]
\centering
\includegraphics[height=3cm]{Domaine_perfore}
\caption{H-Tree}
\label{fig:photo}
\end{figure}

Comme pour le modèle mono-dimensionnel, il existe une matrice qui lie les pressions en sorties et les débits en sorties de l'arbre et on aura
$-p^\varepsilon= A^\varepsilon q^\varepsilon$, et un opérateur à noyau $\mathcal R_\varepsilon$, qui vérifie, comme dans le cas mono-dimensionnel,\vspace*{-3pt}
\[
\mathcal R_\varepsilon u_\varepsilon =\int_\Omega K_\varepsilon(x, y) u_\varepsilon(y) \rd  y,
\]
avec $K_\varepsilon$ constant par morceaux sur les pavés de la décomposition dyadique de $\Omega\times \Omega$ (correspondant aux cellules périodiques).

La composante $k$ du vecteur de flux s'exprime en fonction du flux de la vitesse de la structure élastique et on a\vspace*{-3pt}
\[
q^\varepsilon_k=-\int_{ \Gamma^{\varepsilon,k} }\partial_t \dd^\varepsilon\cdot \nn^\varepsilon.
\]

Par conséquent nous obtenons le système suivant
\[
\begin{cases}
\rho \partial_t^2 \dd^\varepsilon - \div \ \sigma(\dd^\varepsilon) = \f_S, &\text{dans } \Omega^\varepsilon_S,\\
\sigma(\dd^\varepsilon) \cdot \nn^\varepsilon = -\sum_j A^\varepsilon_{k j} \biggl(\dis\int_{\Gamma_j^\varepsilon} \partial_t\dd^\varepsilon \cdot \nn^\varepsilon\biggr) \nn^\varepsilon, & \text{sur chaque } \Gamma^{\varepsilon,k},\\
\sigma(\dd^\varepsilon) \cdot \nn =-p_{\pleu}\nn , & \text{sur } \Gamma_N, \\
\dd^\varepsilon = 0, & \text{sur } \Gamma_D.
\end{cases}
\]
Les conditions aux limites sur le bord des alvéoles sont non locales: elles dépendent d'une part, du flux de la vitesse sur le bord de la cellule fluide considérée, et d'autre part, du flux sur toutes les autres cellules fluides par l'intermédiaire de l'arbre.
Si on suppose que l'on a prolongé $\dd^\varepsilon$ à l'intérieur des alvéoles et donc sur $\Omega^\varepsilon_F$ et que l'on désigne ce prolongement également par $\dd^\varepsilon$, le flux se réécrit
\[
q^\varepsilon_k=\int_\Omega \chi_F^{\varepsilon, k} \div \partial_t\dd^\varepsilon.
\]
En multipliant la première équation du système par une fonction test~$\bb$ (\ie $\bb\in H^1(0,1)$) et en intégrant sur $\Omega^\varepsilon_S$, on obtient alors formellement
\begin{multline*}
\rho \int_{\Omega^\varepsilon_S}\partial_{t}^2\dd^\varepsilon \cdot \bb+ \int_{\Omega^\varepsilon_S} (\sigma(\dd^\varepsilon):\nabla \bb)\\
+\int_\Omega\int_\Omega K_\varepsilon (x, y)(\chi_F^{\varepsilon}\div \partial_t\dd^\varepsilon)(x)(\chi_F^{\varepsilon}\div \bb)(y)\rd  x\rd  y\\
=\int_{\Omega_S^\varepsilon } \f_S\cdot\bb -\int_{\Gamma_N} p_{\pleu} \bb\cdot \nn,
\end{multline*}
où $A:B=\mathrm{tr}(AB^T)$.
Ainsi l'opérateur à noyau agit sur la divergence de la vitesse et l'égalité d'énergie vérifiée par la solution, qui s'obtient en prenant $\bb=\partial_t\dd^\varepsilon$ s'écrit
\begin{multline*}
\underbrace{\frac{\rho}{2} \frac{d}{dt} \int_{\Omega^\varepsilon_S}|\partial_{t}\dd^\varepsilon|^2}_{\substack{\text{dérivée en temps}\\ \text{de l'énergie cinétique}}}+ \underbrace{\frac{1}{2}\frac{d}{dt} \int_{\Omega^\varepsilon_S} (\sigma(\dd^\varepsilon): \mathrm{e}(\dd^\varepsilon))}_{\substack{\text{dérivée en temps}\\ \text{de l'énergie mécanique}}}\\
+\underbrace{\int_\Omega\int_\Omega K_\varepsilon (x, y)(\chi_F^{\varepsilon}\div \partial_t\dd^\varepsilon)(x)(\chi_F^{\varepsilon}\div \partial_t\dd^\varepsilon)(y)\rd  x\rd  y}_{\text{énergie dissipée dans l'arbre}\;\geq 0}\\=\underbrace{\int_{\Omega_S^\varepsilon} \f_S\cdot\partial_t\dd^\varepsilon -\int_{\Gamma_N} p_{\pleu} \partial_t\dd^\varepsilon \cdot \nn.}_{\text{travail des forces extérieures}}
\end{multline*}
Cette égalité est similaire à \eqref{egalite-energie-continu} obtenue dans le cas mono-dimen\-sionnel.

La question est maintenant de savoir quel est le modèle homogénéisé quand $\varepsilon \to 0$.
L'intérêt de tels modèles homogénéisés est qu'ils décrivent le comportement élastique macroscopique tout en prenant en compte la microstructure, mais évitent une résolution du problème microscopique complet.

\subsection{Modèle homogénéisé}
L'idée ici est considérer un développement asymptotique de la solution suivant deux échelles: l'échelle \og lente\fg $x$ et l'échelle \og rapide\fg~$x/\varepsilon$, soit
\[
\dd^\varepsilon(t, x) = \dd^0(t, x) + \eps \dd^1(t, x,x/\eps) + \eps^2 \dd^2(t, x,x/\eps) + \cdots
\]

Les développements asymptotiques peuvent être effectués formellement pour obtenir les problèmes vérifiés par les $\dd^i$ et le passage quand $\varepsilon\to 0$ peut être justifié, par exemple, avec des outils de convergence deux échelles introduit dans \cite{allaire-92a}. Pour une justification de ce passage à la limite on pourra se référer à \cite{cazeaux-grandmont}. Ce passage à la limite se fait sous les mêmes hypothèses sur le réseau de résistances et la convergence de l'opérateur à noyau qui décrit l'action de l'arbre.

On peut montrer que $\dd^1$ est de la forme
\[
\dd^1 (t, x,y)= \sum_{1\leqslant k,l \leqslant d} \mathrm{e}_x(\dd^0)_{k l}(t, x) \boldsymbol{\chi}^{k l}(y) + \partial_t\pi(t, x) \boldsymbol{\chi}_0(y),
\]
où les fonctions $\boldsymbol{\chi}^{k l}$ sont les correcteurs standards obtenus en théorie de l'homogénéisation de matériaux perforés et donc sont ceux que l'on obtiendrait si l'arbre n'était pas présent. Le correcteur $\boldsymbol{\chi}_0$ prend en compte, lui, l'effet de l'arbre sur le comportement mécanique (qui agit comme une pression).

En désignant par $\theta$ la proportion de la structure dans le domaine, $\dd^0$ et $\pi$ satisfont à
\[
\begin{cases}
\theta \rho \partial_{t}^2 \dd^0 - \div \left (\sigma^{\hom} (\dd^0, \pi) \right) = \theta \f_S, & \text{dans } (0,T) \times \Omega, \\
\sigma^{\hom}(\dd^0, \pi) = \mathcal{A}^{\hom} \mathrm{e}(\dd^0) {- \partial_t \pi \mathcal{B}^{\hom}}, & \text{dans } (0,T) \times \Omega, \\
{ \pi+ \tau_{\hom} \partial_t (\mathcal{R} \pi)} = - \mathcal{R} \left (\mathcal{B}^{\hom}: \mathrm{e}(\dd^0)\right) , & \text{dans } (0,T) \times \Omega, \\
\tau_{\hom}= \int_{\partial Y_F}\boldsymbol{\chi}_0\cdot \nn,\\ \mathcal{B}_{kl}^{\hom} =(1-\theta)\delta_{kl} -\int_{\partial Y_F} \boldsymbol{\chi}_{kl}\cdot \nn.
\end{cases}
\]

Les équations macroscopiques obtenues peuvent être vues comme une perturbation du modèle auquel on aboutirait sans arbre. Ici, l’effet de l’arbre se traduit par la présence d’un terme de pression supplémentaire qui contribue à dissiper l’énergie du système, cette pression ayant sa propre dynamique en temps. Cette dynamique dépend du réseau de résistances et induit une relaxation en temps: le modèle est à la fois non local en espace et non local en temps. Comme dans le cas mono-dimensionnel nous obtenons ainsi un modèle visco-élastique.

\begin{remark}[Cas incompressible]
La même démarche peut être effec\-tuée en considérant que le matériau de départ qui constitue le paren\-chyme pulmonaire est incompressible. Cette hypothèse est satis\-faite pour la plupart des tissus vivants qui peuvent être considérés comme quasi-incompressibles. Dans le contexte des petites déformations, cela correspond à ajouter une contrainte sur la divergence du déplacement $\div (\dd^\varepsilon)=0$. Dans ce cas, le matériau homogénéisé structure+arbre est bien sur \og compressible\fg puisque l'air peut s'échapper et cette contrainte n'est plus vérifiée à la limite. On obtient alors un modèle très similaire au modèle mono-dimensionnel:
\[
\begin{cases}
\theta \rho \partial_{t}^2 \dd^0 - \div  \bigl (\sigma_\mathrm{inc}^{\hom} (\dd^0) \bigr) = \theta \f_S, \\[5pt]
\sigma_\mathrm{inc}^{\hom}(\dd^0) = \mathcal{A}_\mathrm{inc}^{\hom} \mathrm{e}(\dd^0) + (1-\theta) \mathcal R(\div  \partial_t\dd^0) I.
\end{cases}
\]
Cette équation est l'équivalent de \eqref{onde-visqueux} dans le cas pluri-dimensionnel.
\end{remark}

\begin{remark}[Approximation numérique]
Pour déterminer des solutions approchées de ce systèmes d'EDP, une difficulté spécifique vient du terme non local en espace qui est non standard et couple tous les points du domaine de calcul par l'intermédiaire de l'arbre. Le~système linéaire associé fait ainsi intervenir des matrices pleines. La~structure arborescente de l'arbre permet néanmoins de définir des algorithmes rapides pour calculer les produits matrices-vecteurs en jeu (voir \cite{cazeaux-grandmont}).
\end{remark}

\subsection{Illustrations numériques: portrait de phase}
Le comportement numérique des solutions montre que lorsque l'on diminue les résistances du réseau, on voit apparaître des effets non locaux de plus en plus marqués. Quand on applique une onde de pression, en plus de se propager dans le milieu élastique, elle se propage à travers l'arbre: c'est ce qui pourrait expliquer pourquoi des impacts violents aux poumons peuvent endommager des tissus loin du point d'impact.

Par ailleurs une question que l'on doit se poser est: est-ce que ce modèle est capable de reproduire des données physiologiques?
On a cherché à reproduire des expériences de spirométrie qui sont effectuées en routine clinique en particulier pour diagnostiquer l'asthme. La patien\-te, ou le patient, respire dans un dispositif médical qui mesure les volumes et débits à la bouche: après une série de respirations normales, il ou elle doit inspirer profondément et lentement, puis enchaîner par une expiration forcée et rapide, afin d'expulser l'air le plus possible des poumons. Les courbes qu'ont l'habitude de lire les praticiens sont des portraits de phase débits--volumes. Nous en avons reproduit ici. On applique une force extérieure au système homogénéisé incompressible qui permet de reproduire une inspiration profonde et une expiration forcée.

\begin{figure}[htb]
\centering
\includegraphics[height=6cm]{Asymmetric}
\caption{Portraits de phase volume-débit - augmentation hétérogène des résistances }
\end{figure}

Contrairement à la convention de signe que nous avons utilisée jusqu'à présent, les débits négatifs (\resp positifs) correspondent à la phase d'inspiration (\resp d'expiration).
Lorsque que l'on augmente les résistances \emph{distales} (\ie des résistances dans les dernières générations de l'arbre) sur une partie de l'arbre, on constate une diminution du débit de pointe à l'expiration et une courbe en fin d'expiration qui s'incurve légèrement se distinguant de la courbe de référence dont la pente est constante. Les volumes inspirés ne sont que peu impactés.

Lorsque que l'on augmente le coefficient de cisaillement sur une partie du domaine élastique, la modification majeure est la quantité d'air inspiré qui diminue fortement. Cette augmentation peut apparaître en cas de fibrose, pathologie qui rigidifie les tissus et donc ne permet plus au patient de ventiler correctement.

\begin{figure}[htb]
\centering
\includegraphics[height=6cm]{Stiffening}
\caption{Portraits de phase volume-débit -- augmentation non uniforme du coefficient de cisaillement ($\mu$)}
\end{figure}

On retrouve ainsi des portraits de phase qui ont l'allure de ceux obtenus en clinique. Cependant pour obtenir des quantités physiologiquement pertinentes, l'intensité de la force appliquée semble être plus petite que celle développée en pratique: en effet une estimation de celle-ci est obtenue grâce à la mesure de la pression \oe sophagienne dont l'ordre de grandeur est 10 fois supérieur à celui considéré ici. Il est donc vraisemblable que certains phénomènes importants aient été négligés et que l'on ne peut pas se contenter de modèles linéaires pour décrire des man\oe uvres d'expiration forcées.

\section{Écoulement de l'air dans les voies aériennes}
Dans les sections précédentes nous avons décrit l'arbre bronchique comme un réseau de résistances en supposant que, dans chaque tube, l'écoulement était de type Poiseuille. Cependant, cette hypothèse peut paraître restrictive puisque l'arbre bronchique peut être schématiquement divisé en trois parties:

\begin{itemize}
\item les voies aériennes supérieures (jusqu'à la septième génération), dans lesquelles on peut supposer que l'air est newtonien visqueux incompressible et donc que les équations de Navier-Stokes sont satis\-faites. L'hypothèse d'incompressibilité est valide car le nombre de Mach est très inférieur à $0.3$, même en respiration forcée;

\item les bronches segmentaires et bronchioles (de la huitième génération à la dix-septième génération) dans lesquelles on peut supposer l'écoulement laminaire;
\item les acinis et sacs alvéolaires (de la dix-septième génération à la vingt-troisième génération) entourés du parenchyme qui sont le lieu de la diffusion des gaz dans le sang.
\end{itemize}

Ainsi dans la partie supérieure de l'arbre, dont la géométrie est accessible grâce à l'imagerie médicale, l'écoulement de l'air peut être décrit par les équations de Navier-Stokes tridimensionnelles. Enrichir le modèle de cette manière permet en particulier de prendre en compte certaines non-linéarités qui n'étaient pas prises en compte précédemment. Toutefois on peut difficilement envisager de mener une simulation directe de ces équations sur l'arbre complet: d'une part à cause des coûts de calcul, d'autre part à cause des limites de l'imagerie médicale classique qui ne permet de segmenter l'arbre bronchique au delà de la dixième génération.
Il s'agit donc de déterminer des conditions aux limites permettant de prendre en compte la partie non simulée ou encore d'obtenir des modèles réduits que l'on couplera ensuite au reste de l'arbre.
Une première solution est de caler des pressions en sortie de l'arbre (typiquement $2^7$ sorties si on considère un arbre à $7$ générations). Malheureusement de telles mesures expérimentales ne sont pas disponibles. Une deuxième solution est de dériver des modèles réduits permettant de décrire la partie distale de l'arbre de manière simplifiée et facilement paramétrable. Ici,~on adopte donc une perspective inverse de la démarche précédente: \hbox{décrire} précisément l'écoulement de l'air dans la partie proximale et prendre en compte la partie distale de manière simplifiée. Le modèle de la partie distale va s'inspirer du modèle d'EDO déjà présenté.\enlargethispage{-2\baselineskip}

\subsection{Modélisation}
Soit $\Omega$ un ouvert de ${\mathbb R}^d$ avec $d=2, 3$, représentant les premières générations de l'arbre bronchique. On suppose que, dans cette partie proximale, l'air est newtonien visqueux incompressible et satisfait aux équations de Navier-Stokes. Dans les arbres distaux, on fait l'hypothèse que la loi de Poiseuille est vérifiée dans chacune des branches. Ainsi l'écoulement dans le $i$-ème sous-arbre est caractérisé par une résistance équivalente $R_i$.
On suppose ensuite que l'arbre bronchique débouche dans une boîte dont un des bords de surface $S$, mobile dans une seule direction, est relié à un ressort de raideur $k$. Ce bord mobile modélise le comportement du diaphragme, moteur de la respiration.
Ainsi, le modèle couplé peut être schématiquement décrit par la figure~\ref{fig:model}.

\begin{figure}[htb]
\centering
\includegraphics[scale = .55]{geogif2}
\caption{Modèle multi-échelles}
\label{fig:model}
\end{figure}

Dans le domaine $\Omega$, la vitesse $\uu$ et la pression $p$ du fluide vérifient les équations de Navier-Stokes avec une condition de Dirichlet homogène sur la paroi latérale $\Gamma_l$ et des conditions de Neumann sur les sorties $\Gamma_i$. Les pressions $p_i$ aux sorties artificielles $\Gamma_i, i\neq 0$, sont inconnues et liées à la pression alvéolaire $P_{\ra}$ par la loi de Poiseuille\enlargethispage{\baselineskip}
\[
p_i-P_{\ra}= R_i\biggl(\int_{\Gamma_i}\uu\cdot\nn\biggr),
\]
où
$\int_{\Gamma_i}\uu\cdot\nn$ est le débit du fluide à la sortie $\Gamma_i$.
Les conditions aux bords s'écrivent donc
\begin{equation}
\label{conditions:0D}
\nu\frac{\partial \uu}{\partial \nn}- p \nn = -P_{\ra} \nn -\biggl(R_i\int_{\Gamma_i}\uu\cdot\nn\biggr)\nn, \quad \text{sur }\Gamma_i,
\end{equation}
où $\sfrac{\partial \uu}{\partial \nn}= \nabla \uu\cdot \nn$ et $\nu$ désigne la viscosité du fluide.

Ces conditions sont des conditions \og dissipatives\fg qui correspondent à un couplage $3D$-$0D$. Elles représentent le modèle réduit le plus simple appelé \emph{lumped model}. Elles sont utilisées également pour la modélisation des écoulements sanguins (voir par exemple \cite{vignon-figuera-jansen-taylor-06}, \cite{veneziani-vergara-05}, \cite{quarteroni-veneziani-03}, \cite{quarteroni-ragni-veneziani-01}). Dans le cas des écoulements sanguins, des modèles réduits hyperboliques sont également utilisés \cite{formaggia-gerbeau-nobile-quarteroni-01}, \cite{formaggia-moura-nobile-07}, correspondant au phénomène de propagation d'une onde de pression. Cependant dans le cas de la ventilation, de~tels phénomènes de propagation d'ondes de pression n'apparaissent pas et l'hypothèse selon laquelle le fluide vérifie la loi Poiseuille dans la partie distale est justifiée. La pression $P_{\ra}$ ainsi qu'une force extérieure $f_{\ext}$, qui représente la force développée par le diaphragme, agissent sur le piston, dont la position est notée~$x$. Enfin on suppose que le tissu élastique qui remplit la boîte est incompressible, ce qui se traduit par la relation
\begin{equation}
\label{incompressibilite:parenchyme}
V'= Sx'=\sum_{i=1}^N\int_{ \Gamma_i} \uu\cdot \nn=-\int_{ \Gamma_0} \uu\cdot \nn,
\end{equation}
où $S$ représente la surface de la boîte.
\begin{rem}
Dans l'égalité précédente, on a utilisé l'incompressibilité du fluide et le fait que la vitesse est nulle sur la paroi latérale considérée comme fixe. Si on avait tenu compte du mouvement de la paroi alors la vitesse de celle-ci serait apparue dans \eqref{incompressibilite:parenchyme}.
\end{rem}
Les équations globales du système s'écrivent donc:\vspace*{-5pt}
\begin{equation}
\label{eq:Pe1}
\begin{cases}
\rho_f\partial_t\uu+\rho_f
(\uu\cdot\nabla)\uu-\nu \Delta \uu + \nabla p ={\boldsymbol 0}, &\!\text{dans } (0,T)\times\Omega,\ \\
\nabla \cdot\uu =0, &\!\text{dans} (0,T)\times\Omega\\
\uu= 0, &\!\text{sur } (0,T)\times\Gamma_l,\\
\nu \nabla \uu \cdot \nn
- p \nn =-P_0\nn &\!\text{sur } (0,T)\times\Gamma_0,\\
\nu \nabla \uu \cdot \nn
- p \nn =-P_{\ra}\nn
- R_i \bigl (\int_{\Gamma_i} \uu\cdot \nn\bigr)\nn,
&\!\biggl\{\begin{array}{l}
\text{sur }(0,T)\times\Gamma_i,\\
 i=1,\ \dots,\ N,
\end{array}\\
mx''+kx=f_{\ext}+SP_{\ra}, \\
Sx'=\dis\sum_{i=1}^{N}\int_{\Gamma_i} \uu\cdot\nn=-\int_{\Gamma_0} \uu\cdot\nn,
\end{cases}
\end{equation}
où $\rho_f$ désigne la densité du fluide et $\nu$ sa viscosité. La première équation traduit la conservation de la quantité de mouvement du fluide, la deuxième équation l'incompressibilité du fluide. Sont ensuite prescrites les conditions aux limites, qui elles-mêmes dépendent de la dyna\-mique du ressort de masse $m$ et de raideur $k$ (cette raideur $k$ est l'analogue de l'élastance $E$ introduite précédemment).

Ce problème admet comme formulation variationnelle, en supposant que $\vv\in H^1(\Omega)$ est une fonction test à divergence nulle, la formulation suivante:\vspace*{-5pt}\enlargethispage{\baselineskip}
\begin{multline*}
\label{eq:vfdiv0}
\int_\OO \rho_f\partial_t \uu \cdot \vv+\int_\OO \rho_f(\uu\cdot\nabla)\uu \cdot\vv+\nu\int_\OO \nabla\uu:\nabla\vv\\[-3pt]
\shoveright{+\sum_{i=1}^N R_i \biggl(\int_{\GG_i}\uu\cdot \nn\biggr)\biggl(\int_{\GG_i}\vv\cdot \nn\biggr)}\\
\shoveright{+\frac{m}{S^2}\biggl(\int_{\GG_0}\partial_t\uu\cdot \nn\biggr)\biggl(\int_{\GG_0}\vv\cdot \nn\biggr)+\frac{k}{S^2}\biggl(\int_0^t\int_{\GG_0}\uu\cdot \nn\biggr)\biggl(\int_{\GG_0}\vv\cdot \nn\biggr)}\\
\shoveright{=- P_0\int_{\GG_0}\vv\cdot \nn-\frac{f_{\ext}}{S}\int_{\GG_0}\vv\cdot \nn +\frac{k}{S}x_0\biggl(\int_{\GG_0}\vv\cdot \nn\biggr),}\\
\forall \vv \; \text{telle que } \nabla\cdot\vv=0, \vv=0 \; \text{sur } \Gamma_l.
\end{multline*}
Tout a été écrit en fonction de la vitesse du fluide en utilisant l'incompressibilité globale du système \eqref{incompressibilite:parenchyme}. Il est à noter que la forme bilinéaire\vspace*{-5pt}
\begin{equation}
\label{forme:bi:R}
a(\uu, \vv)=\nu\int_\OO \nabla\uu:\nabla\vv
+\sum_{i=1}^N R_i \biggl(\int_{\GG_i}\uu\cdot \nn \biggr) \biggl(\int_{\GG_i}\vv\cdot \nn \biggr)
\end{equation}
associée au système n'est pas la forme bilinéaire standard associée
au problème de Stokes. Elle contient des termes supplémentaires associés à la dissipation venant des sous-arbres et faisant intervenir le produit des flux sur les sorties artificielles $\Gamma_i, i\neq 0$. Cette particularité nécessite l'élaboration de stratégies numériques adaptées.

En prenant $\uu$ comme fonction test et en rappelant que\vspace*{-3pt}
\[
V'=Sx' =-\int_{\Gamma_0}\uu\cdot \nn,
\]
on obtient formellement l'égalité d'énergie suivante:\vspace*{-5pt}
\begin{multline*}
\label{eq:bilan_energie}
\frac{d}{dt}\underbrace{\left (\dis\frac{\rho_f}{2}\dis\int_{\OO}\dis|\uu|^2+\frac{m}{2}|x'|^2+\frac{k}{2}|x|^2\right) }_{\text{énergie totale}}
\dis+\underbrace{\nu\dis\int_{\OO}|\nabla\uu|^2}_{\substack{\text{dissipation}\\\text{dans }\OO}}
+\underbrace{\dis\sum_{i=1}^NR_i\biggl(\int_{\GG_i}\uu\cdot\nn\biggr)^2
}_{\substack{\text{dissipation}\\\text{dans les sous-arbres}}}
\\[-5pt]
=-\underbrace{\dis\sum_{i=0}^{N}\int_{\GG_i} \frac{\rho_f}{ 2}|\uu|^2(\uu\cdot\nn)
}_{\text{flux d'énergie cinétique}}
+
\underbrace{P_0S x'}_{\text{puissance à l'entrée}} +\underbrace{f_{\ext}x'}_{\substack{\text{puissance}\\\text{des forces extérieures}}}.
\end{multline*}
Le signe du terme cubique en $\uu$ représentant le flux d'énergie cinétique aux frontières est indéterminé.
Cela peut poser un problème aussi bien d'un point de vue mathématique pour l'obtention d'estimations \emph{a priori}
que d'un point de vue numérique pour la stabilité des schémas \cite{these-irene}.

\subsection{Existence de solutions}
Si l'on se pose la question de l'existence de solutions pour ce problème couplé, une des premières difficultés est l'obtention d'estimations d'énergie et donc le contrôle des termes de flux d'énergie cinétique aux frontières artificielles. On a
\begin{multline*}
\biggl|\sum_{i=0}^{N}\int_{\GG_i}\frac{ \rho_f} {2}\, |\uu|^2(\uu\cdot\nn)\biggr|=\biggl|\int_{\Omega} \rho_f(\uu\cdot \nabla)\uu \uu\biggr|\\
\leq C \begin{cases}
\|\uu\|_{L^2(\Omega)} \|\nabla\uu\|^2_{L^2(\Omega)}& \text{si } d=2,\\[2ex]
\|\uu\|^{\sfrac{1}{2}}_{L^2(\Omega)} \|\nabla\uu\|^{\sfrac{5}{2}}_{L^2(\Omega)}&\text{si } d=3.
\end{cases}
\end{multline*}
Par conséquent une estimation d'énergie pour le système couplé peut être obtenue uniquement dans le cas $d=2$, pour des données petites et un temps petit.
Il est donc nécessaire soit de chercher d'autres estimations soit de modifier le problème afin d'obtenir des estimations d'énergie.
Une première possibilité est de considérer la pression totale $p+\Psfrac{\rho_f}{2}|\uu|^2$ à la place de la pression $p$ dans les contraintes fluides (voir \eqref{incompressibilite:parenchyme}) ou d'ajouter un terme du type $\Psfrac{\rho_f}{2}\uu (\uu\cdot \nn)$ aux contraintes normales fluides. On peut alors obtenir l'existence de solu\-tions faibles.
Une deuxième possibilité (voir \eqref{conditions:0D}) est de contraindre les vitesses du fluide à être proportionnelles à un profil donné:\vspace*{-3pt}
\[
\uu(t, \xx)= \lambda_i (t) {\boldsymbol U}_i (\xx) \quad\text{sur } \Gamma_i,
\]
où $ {\boldsymbol U}_i $ est donné et $\lambda_i$ est une inconnue du problème.

Dans ce cas, les conditions aux limites en contraintes s'écrivent\vspace*{-3pt}\enlargethispage{\baselineskip}
\[
\int_{\Gamma_i}\Bigl(\nu\, \frac{\partial \uu}{\partial \nn}-p \nn\Bigr)\cdot {\boldsymbol U}_i = -P_{\ra} \int_{\Gamma_i}{\boldsymbol U}_i\cdot \nn - R_i \lambda_i\biggl(\int_{\Gamma_i} {{\boldsymbol U}_i\cdot \nn}\biggr)^2.
\]
Sous ces conditions\vspace*{-3pt}
\[
\biggl|\sum_{i=0}^{N}\int_{\GG_i}\frac{ \rho_f} {2}\, |\uu|^2(\uu\cdot\nn)\biggr|= \sum^N_{i=0} \gamma_i \lambda_i^3,
\]
où les $\gamma_i$ sont des constantes dépendant des données. Ainsi le terme de convection s'écrit à l'aide d'un nombre fini de degrés de liberté. Or chacun des coefficients $\lambda_i$ peut être contrôlé par la norme $L^2$ de la vitesse $\uu$. En effet, en supposant que ${\boldsymbol U}_i\cdot\nn \neq 0$ et en tenant compte
de l'incompressibilité du fluide, on peut montrer que
\[
|\lambda_i|\leq C_i \|\uu\|_{L^2(\Omega)}.
\]
Cette estimation permet donc de contrôler le terme convection par $\|\uu\|_{L^2(\Omega)}^3$ et entraîne donc l'obtention d'estimations \emph{a priori} au moins en temps petit. L'existence de solutions faibles peut alors être démontrée par la méthode de Galerkin. La compacité forte des vitesses, nécessaire pour passer à la limite dans le système approché, est obtenue de manière classique \cite{girault-raviart-86} grâce à des estimations supplémentaires sur des dérivées fractionnaires en temps de $\uu$ et par application du lemme d'Aubin.
L'existence pour tout temps est, quant à elle, obte\-nue pour des données (forces et conditions initiales) suffisamment petites.

\begin{rem}
Il est à noter que numériquement une telle contrainte stabilise également le système \cite{taylor}.
\end{rem}

Maintenant, si l'on considère le cas général non contraint, on peut démontrer l'existence d'une unique solution régulière \cite{baffico-grandmont-maury}.

Le cadre fonctionnel dans lequel on se place est le suivant:
\begin{align*}
V&=\{\vv\in H^1(\Omega)^d\mid \nabla\cdot\vv=0,\, \vv=0 \mbox{ sur } \Gamma_l\},
\\
H&=\overline{V}^{L^2}\quad\text{(adhérence de $V$ dans $L^2(\Omega)$).}
\end{align*}
Sur $H$, on ne va pas mettre le produit scalaire usuel de $L^2(\Omega)$ mais un produit scalaire prenant en compte
la partie inertielle venant du ressort. Soit $(\cdot\,, \cdot)_{H}$ le produit scalaire $H\times H $ défini par
\begin{equation}
\label{eq:defsp0}
(\vv, \ww)_{H}=\rho_f\int_\OO \vv\cdot\ww+ \frac m {S^2}\biggl(\int_{\GG_0}\vv\cdot\nn\biggr)\biggl(\int_{\GG_0}\ww\cdot\nn\biggr),
\end{equation}

\begin{rem}
Le flux $\int_{\GG_0}\vv\cdot\nn$ d'une fonction $\vv$ dans $H$ n'est en fait pas défini. Cependant on peut en donner une définition rigoureuse par dualité $ (H^{\sfrac{1}{2}})'-H^{\sfrac{1}{2}}$, définition correspondant au flux dans le cas de fonctions régulières.
\end{rem}
De plus, on introduit un opérateur $A$ adapté au problème dont les modes propres constitueront la base de Galerkin utilisée pour définir un problème approché.
\BE
\label{def-op}
\begin{split}
D(A)&=\{\vv \in V\mid |a(\vv, \ww)|\leq C \|\vv\|_H, \, \forall \ww\in V\},\\
(A\vv, \ww)_{H}&=a(\vv, \ww), \quad \forall (\vv, \ww)\in D(A)\times V,
\end{split}
\EE
où $a$ est la forme bilinéaire définie en \eqref{forme:bi:R}.

Cet opérateur est un opérateur de type Stokes avec des conditions aux limites mixtes de Dirichlet dissipatives définies par \eqref{conditions:0D}. Il est auto-adjoint, inversible de $D(A)$ dans $H$ et d'inverse compact.

On dira alors que $\uu$ est une solution de \eqref{eq:Pe1} si
\BE
\label{FV:NS:piston}
\begin{cases}
\uu\in L^2(0,T;D(A))\cap L^{\infty}(0,T;V),\quad \partial_t\uu \in L^2(0, T; H),\\
\dis \frac{d}{dt}(\uu, \vv)_{H}+a(\uu, \vv)+ b(\uu, \uu, \vv)\\
\hspace*{4.5cm}+c\biggl (\dis\int_0^t \uu, \vv\biggr)= l(\vv)\quad \forall \vv\in V, \\
\end{cases}
\EE
où $b$ est la forme trilinéaire donnée par
\[
b(\uu, \vv, \ww)= \rho_f \int_{\OO} (\uu\cdot\nabla)\vv\cdot\ww,
\]
$c$ est la forme bilinéaire associée à l'énergie mécanique du ressort
\[
c(\vv, \ww)=\frac{k}{S^2} \biggl(\int_{\GG_0}\vv\cdot \nn\biggr)\biggl(\int_{\GG_0}\ww\cdot \nn\biggr),
\]
et
\[
l(\vv)=- P_0\int_{\GG_0}\vv\cdot \nn-\frac{f_{\ext}}{S}\int_{\GG_0}\vv\cdot \nn +\frac{k}{S^2}\,Sx_0\biggl(\int_{\GG_0}\vv\cdot \nn\biggr).
\]

On démontre:
\begin{thrm}
Sous l'hypothèse que $\Gamma_l$ et les $\Gamma_i$ se rencontrent à angle droit, si la vitesse initiale $\uu_0\in V$ alors il existe une unique solution de \eqref{eq:Pe1} au sens défini par \eqref{FV:NS:piston}.
\end{thrm}

La démonstration de l'existence de solutions est basée sur une méthode de Galerkin, où les fonctions de la base sont les vecteurs propres de $A$. En particulier en prenant formellement successivement $A\uu$ et $\partial_t\uu$ comme fonctions tests on obtient des estimations d'énergie. Il est à noter que $A\uu$ est, en particulier, à divergence nulle grâce au choix de $A$.

La première estimation (en prenant $\vv=A\uu$), valable en temps petit, repose sur des propriétés de régularité de la solution de $A\ww=\f$ avec $\f \in H$. En effet, on estime le terme de convection de la façon suivante
\[
\left|b(\uu, \uu, A\uu)\right|\leq \|\uu\|_{L^\infty(\Omega)}\|\nabla\uu\|_{L^2(\Omega)}\|A\uu\|_{L^2(\Omega)}.
\]
De plus, si on se place en dimension 3, si $A\uu\in H$, on peut montrer qu'il existe $0<\varepsilon<\sfrac{1}{2}$ tel que $\uu\in H^{3/2+\varepsilon}(\Omega)$ (pour un ouvert $\Omega$ de $\mathbb R^3$, l'espace $H^{3/2+\varepsilon}(\Omega)$ est l'ensemble des fonctions $\vv\in L^2(\Omega)$, telles que $\nabla\vv\in L^2(\Omega)$ et $(x, y)\mto \psfrac{\nabla \vv(x)-\nabla\vv(y)}{|x-y|^{2+\varepsilon}}$ appartient à l'espace $L^2(\Omega \times \Omega)$). En dimension $3$, $H^{3/2+\varepsilon}(\Omega)$ s'injecte dans $L^\infty(\Omega) $. Alors, par interpolation hilbertienne, il existe $\theta\in (0, 1)$ tel que
\[
\|\uu\|_{L^\infty(\Omega)}\leq C \|\nabla\uu\|^{\theta}_{L^2(\Omega)}\,\|A\uu\|^{1-\theta}_{L^2(\Omega)}.
\]
Par conséquent
\[
\left|b(\uu, \uu, A\uu)\right|\leq C \|\nabla\uu\|^{1+\theta}_{L^2(\Omega)}\,\|A\uu\|^{2-\theta}_{L^2(\Omega)},
\]
ce qui permet d'obtenir, grâce à un lemme de Gronwall non liné\-aire, des estimations en temps petit de $\uu$ dans l'espace $L^2(0, T; D(A))\cap L^\infty(0, T; V)$.
\begin{rem}
L'étude des propriétés de régularité des solutions de $A\ww=\f$ se ramène à l'étude de la régularité des solutions d'un problème de Stokes avec conditions mixtes Dirichlet-Neumann. Dans le cas où les $\Gamma_i$ et la frontière latérale $\Gamma_l$ sont régulières et se rencontrent à angle droit, on montre que $\ww\in H^{3/2+\varepsilon}(\Omega)$ (\cite{orlt-sandig-95}, \cite{rossmann-mazja-07}).

\end{rem}
On montre par ailleurs l'unicité de la solution en raisonnant par l'absurde et en considérant la différence de deux solutions comme fonction test. L'existence globale en temps est, quant à elle, démontrée dans le cas où les conditions initiales et les forces appliquées sont suffisamment petites et si $k=0$. Cette restriction vient du fait que, lorsque $k\neq 0$, le système n'est pas dissipatif.

\subsection{Remarque sur la discrétisation du système}
Comme nous l'avons fait remarquer précédemment, à cause des conditions aux limites non standard sur les frontières artificielles, la~forme bilinéaire \eqref{forme:bi:R} associée au problème fait intervenir le flux des vitesses sur les $\Gamma_i$. Ainsi, si l'on cherche à discrétiser le problème continu à l'aide de la méthode des éléments finis, ces termes couplent tous les degrés de liberté des interfaces $\Gamma_i$.
On peut pour les traiter envisager plusieurs stratégies:
\begin{itemize}
\item Effectuer un traitement explicite. Cette solution permet d'utiliser les fonctionnalités de n'importe quel code résolvant les équations de Navier-Stokes. Cependant dans notre cas particulier, cela implique que l'on ne résout pas les équations du ressort et pose des problèmes de stabilité numérique importants. De façon générale, traiter le terme de flux explicitement dans les conditions aux limites \eqref{conditions:0D} pose des problèmes de stabilité des schémas quand les résistances sont grandes. D'autre part, si l'on n'intègre pas les termes résistifs dans la forme bilinéaire du système, l'existence de solutions peut être obtenue mais sous une condition de petitesse des $R_i$, illustrant la possibilité d'instabilités numériques dans le cas de traitement explicite. Les couplages explicites ont été analysés dans \cite{grandmont-martin}.
\item Effectuer un traitement implicite. Cela implique de modifier la matrice de rigidité du problème qui est la matrice associée à la forme bilinéaire qui prend en compte les résistances des sous-arbres ou, si la résolution du système linéaire repose sur un algorithme itératif, modi\-fier uniquement le produit matrice-vecteur. En effet, cette matrice n'a alors pas une structure éléments finis usuelle et prend en compte la dissipation venant des sous-arbres. Se pose alors la question du conditionnement et de la construction de préconditionneurs adaptés.
\item Élaborer de nouvelles méthodes numériques qui permettent de résoudre le problème en ne résolvant que des problèmes standards. On renvoie à la thèse d'Anne Devys \cite{ADevys} pour un algorithme de ce type.\enlargethispage{-\baselineskip}
\end{itemize}

\subsection{Simulations}
Nous allons maintenant présenter quelques simulations numériques démontrant que le modèle permet d'obtenir des écoulements instationnaires tridimensionnels en géométrie réelle dont le seul moteur est la force extérieure exercée sur le ressort que l'on choisit en créneau en temps (qui modélise la force musculaire appliquée par le diaphragme). La géométrie utilisée est la même que celle employée dans \cite{DeRochefort}, \cite{Fetita-etal}. De plus, les simulations illustrent la capacité de ce modèle, en jouant sur les paramètres, à reproduire certains effets de pathologies comme une crise d'asthme (augmentation des résistances $R_i$ des sous-arbres distaux) ou un emphysème (diminution de la raideur $k$).

\begin{figure}[htb]
\begin{center}
\includegraphics[scale = .23]{P_exp_norm_autocb_view2}\qquad
\includegraphics[scale = .23]{P_exp_obs_autocb_view2_circ}
\end{center}
\caption{Champs de pression au pic d'expiration dans un cas normal (à gauche) et obstrué (à droite). Seules les sorties pour lesquelles la résistance $R_i$ a été augmentée sont entourées.}\label{fig:pres}
\end{figure}

Pour les valeurs des coefficients utilisés et physiologiquement pertinents, on renvoie à \cite{baffico-grandmont-maury} d'où proviennent les illustrations des figures \ref{fig:pres}, \ref{fig:mesh1}, \ref{fig:cut1}. On a en particulier estimé les valeurs des résistances à l'aide des données anatomiques de \cite{Wei63}. On constate une différence quantitative de la pression moyenne dans la trachée, et ce, uniquement en ayant augmenté quatre résistances sur 68 (figure \ref{fig:pres}).

Nous avons d'autre part regardé les isovaleurs du module de la vitesse dans deux plans de coupe après la première bifurcation (figu\-res~\ref{fig:mesh1} et \ref{fig:cut1}).

\begin{figure}[htb]
\centering \includegraphics[scale = .8]{bronche_leftrightplane}
\vspace*{-2mm}
\caption{Géométrie reconstruite. Plans de coupe après la première bifurcation.}
\label{fig:mesh1}
\end{figure}

\begin{figure}[htb]
\begin{center}
\resizebox{0.8\textwidth}{!}{
\includegraphics[height=0.4\hsize,angle=90]{bronche_leftplaneins_cb_new_dxeps}
\includegraphics[height=0.4\hsize,angle=90]{bronche_leftplanetrans_cb_new_dxeps}
\includegraphics[height=0.4\hsize,angle=90]{bronche_leftplaneexp_cb_new_dxeps}}
\end{center}
\vspace{-2mm}
\caption{Isovaleurs de la vitesse dans le plan de coupe droite,
au maximum de l'inspiration ($t=0.4s$, à gauche), à la transition ($t=1.6s$, au milieu), au maximum de l'expiration ($t=1.9s$, à droite)}
\label{fig:cut1}
\end{figure}

\begin{figure}[htb]
\begin{center}
\resizebox{0.8\textwidth}{!}{
\includegraphics[height=0.45\hsize,angle=90]{bronche_rightplaneins_cb_new_dxeps}
\includegraphics[height=0.45\hsize,angle=90]{bronche_rightplanetrans_cb_new_dxeps}
\includegraphics[height=0.45\hsize,angle=90]{bronche_rightplaneexp_cb_new_dxeps}}
\end{center}
\vspace{-2mm}
\caption{Isovaleurs de la vitesse dans le plan de coupe gauche,
au maximum de l'inspiration ($t=0.4s$, à gauche), à la transition ($t=1.6s$, au milieu), au maximum de l'expiration ($t=1.9s$, à droite)}
\label{fig:cut2}
\end{figure}

On constate la présence de la forme en \emph{M} caractéristique de la transition (voir par exemple \cite{DeRochefort2} où des expériences in vitro d'écoulement dans un moulage de bronches et des simulations numériques ont été réalisées sur une même géométrie).

D'autre part, des expériences d'expiration forcée (avec un modèle de ressort non linéaire, voir \cite{MM07}) peuvent également être réalisées. Le but est d'obtenir des portraits de phase volume-débit à la bouche.

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=.8\textwidth]{knorm_kvar_pp}
\end{center}
\caption{ Courbes déplacement--débit (en m et $\mathrm{m}^3\cdot \mathrm{s}^{-1}$)
pour différentes valeurs de $k$: $ k\!=\!40.172 + 10^{s-3}$, $s\!=\!0, \dots, 6$.
La zone des débits positifs correspond à la phase d'expiration et en inspiration les débits sont comptés négativement. Le~début de l'expiration se fait au point de volume minimal. Dans une première phase le patient respire normalement, puis inspire profondément, expire et inspire une nouvelle~fois.}
\label{fig:pp_knorm_kvar}
\end{figure}
Ces courbes illustrent la capacité du modèle à reproduire des expé\-riences de spirométrie. Mais une des questions est: quelle est l'adéquation réalité/simulations en ce qui concerne l'écoulement tridimensionnel? De plus encore une fois, les efforts appliqués pour retrouver les bons ordres de grandeurs des volumes et débits semblent être en dessous des ordres de grandeurs réels bien que des non-linéarités aient été prises en compte: non-linéarités venant du fluide grâce à la réso\-lu\-tion des équations de Navier-Stokes dans la partie proximale et prise en compte d'une raideur non linéaire (non explicitée ici). Cela indique encore une fois que le modèle néglige certains phénomènes cruciaux à prendre en compte pour la description de man\oe uvres de respiration forcée. En particulier lors de l'expiration, l'hypothèse de rigidité des bronches est sûrement questionnable.

\section{Applications des modèles couplés arbre-parenchy\-me}
Ici on présente des résultats obtenus lors de la thèse Cifre de Nico\-las Pozin en collaboration avec Air Liquide \cite{Pozin-2017} où un modèle d'arbre irriguant un domaine $3D$ élastique incorporant des données d'imagerie médicale a été mis en \oe uvre numériquement et simulé afin d'explorer deux questions:
\begin{itemize}
\item Peut-on voir un impact sur la ventilation en fonction des mélanges gazeux inspirés?
\item Peut-on utiliser nos modèles afin de détecter des pathologies qui affectent l'arbre bronchique, tel l'asthme ou les broncho-pneumopathies chroniques obstructives (BPCO)?
\end{itemize}

Dans les deux cas, le modèle considéré est le modèle suivant: on se donne $\Omega\subset \mathbb R^3$ et on décompose $\Omega$ en $N$ sous-domaines $\Omega_i$ et chacun de ses sous-domaines est irrigué par une sortie d'un arbre résistif. Dans ce domaine, on suppose que l'on a un matériau élastique et on cherche à résoudre le problème suivant:\vspace*{-3pt}\enlargethispage{1.5\baselineskip}
\begin{equation}
\label{syst_fini}
\begin{cases}
\rho_s\partial_{t}^2\dd -\div  \sigma(\dd) +\nabla p = \mathbf{0}, \quad\text{dans } \Omega,\\
\dis p_{|_{\Omega_i}}=-\sum_i^N A_{ij}\int_{\Omega_j} \div  \partial_t \dd,\\
\dd=\uu_d, \quad\text{sur } \partial \Omega,
\end{cases}
\end{equation}
où $A_{ij}$ représente l'arbre, et $\sigma$ est le tenseur des contraintes qui satisfait à~\eqref{hooke}.
Le système est discrétisé en temps et résolu par la méthode des éléments finis.
Les hypothèses que l'on fait sur les données sont les suivantes:
\begin{itemize}
\item la géométrie est reconstruite à partir d'images médicales: à la fois l'enveloppe du parenchyme pulmonaire, les premières bronches de l'arbre et, dans la première application, les résistances qui leurs sont associées (basées sur une estimation des rayons et des longueurs);
\item on connaît $\uu_d$ qui est le déplacement de la surface du poumon. En effet contrairement à la pression appliquée au poumon qui n'est pas mesurable, le déplacement de la surface peut être déduit à partir de données d'imagerie. Ici il a été reconstruit à partir de deux images à deux niveaux d'inflation différents;
\item on se donne les paramètres élastiques du parenchyme pulmonaire (ici $\lambda$ et $\mu$).
\end{itemize}

\subsection{Application à l'étude de l'influence des gaz porteurs}
Le but de cette étude était de regarder si, lorsque l'on remplace l'azote par l'hélium dans les mélanges inhalés par les patients et patien\-tes, ceux-ci et celles-ci en tirent un bénéfice en terme d'effort respiratoire. La différence entre l'air et les mélanges hélium-oxygène est en particulier la densité du mélange bien plus faible dans le cas de l'hélium-oxygène. Or la densité du gaz n'apparaît pas dans le modèle précédent. Le modèle a donc été modifié afin de prendre en compte les effets inertiels du fluide en considérant que la résistance dans un tube cylindrique de longueur $l$ et de diamètre $d$ s'écrit\vspace*{-3pt}
\[
R=\max (R_{\pois}, \gamma(\Re\,\sfrac{d}{l})^{1/2}R_{\pois}),
\]
où la résistance de Poiseuille $R_{\pois}$ est défini par \eqref{R_Poiseuille}, et $\Re$ désigne le Reynolds dans la bronche: $\Re=\sfrac{4\rho_f Q}{\pi\nu d}$.
Cette formule est une formule calée à partir d'expériences (voir \cite{pedley}). Pour des débits faibles la résistance est la résistance de Poiseuille et pour des débits élevés on voit alors apparaître les effets inertiels.
Nous avons donc calculé numériquement une valeur approchée de la solution de~\eqref{syst_fini}, sur la même géométrie mais en faisant varier les rayons des branches pour simuler des broncho-constrictions. Le paramètre de sortie est le travail de la force exercée par le patient pour respirer à l'inspiration, noté $W$.
\begin{figure}[!ht]
\begin{center}
\includegraphics[scale=.43,angle=-90]{heliox}
\end{center}
\caption{Moyennes et écarts-types sur un cycle respira\-toire dans différentes configurations: un patient en bonne santé, cinquante configurations d'asthmatiques avec bron\-cho-cons\-trictions,
cinquante patients avec seulement des bron\-cho-cons\-trictions modérées diffuses (plug free) et cinquante patients avec seulement des
broncho-cons\-trictions importantes isolées (plug uniquement). Comparaison des efforts pour de
l'air ou un mélange hélium-oxygène.}
\label{Heliox}
\end{figure}
On constate ainsi une grande variabilité entre \og patients\fg simulés, ce qui est constaté également dans les études cliniques. Le mélange d'hélium-oxygène apporte une amélioration de 24\% en moyenne et jusqu'à 30\%. On note que $W$
est augmenté en moyenne d'un facteur~$6.7$ des cas sains au cas d'asthme en respirant de l'air mais d'un
facteur $5.2$ lors d'inhalation d'hélium-oxygène. La~variabilité est \hbox{réduite} lorsque seules des broncho-constrictions modé\-rées
sont simulées, mais il y a un bénéfice à respirer de l'hélium-oxygène dans tous les cas pathologiques. Noter qu'en moyenne, la respiration est plus difficile avec de nombreuses constrictions modérées qu'avec quelques graves; mais la variabilité est beaucoup plus élevée dans ce dernier~cas.%\enlargethispage{\baselineskip}

\subsection{Application à la détection de broncho-constriction}
Dans cette dernière sous-section, on illustre comment des simulations numériques peuvent être utilisées comme base d'apprentissage: on parle alors de données synthétiques. Elles peuvent contribuer à l'identification de paramètres du modèle adaptés au patient.

On souhaite ici pouvoir identifier les, ou certaines, résistances de l'arbre bronchique. On suppose que pour un patient ou une patiente les variations de volumes dans les sous-régions $\Omega_i$ du parenchyme pulmonaire sont connues grâce à des données d'imagerie. On suppose également que l'on sait mesurer le déplacement de chaque point de $\Omega$ que l'on note $\dd_{\mes}$. Si on avait un moyen de calculer les pressions en chaque n\oe ud de l'arbre, on aurait alors accès à toutes les résistances de l'arbre grâce à une relation de type \eqref{eq_tube}.

Dans un premier temps nous utilisons un algorithme de classification afin de réduire le nombre de bronches qui ont le plus probablement des constrictions, les résistances des autres branches étant supposées connues et fixées à une valeur saine, de façon à ce que l'on soit dans une situation où, connaissant les pressions en bout d'arbre, on puisse remonter aux pressions aux n\oe uds les branches identifiées comme pathologiques.
Cet algorithme d'apprentissage s'appuie sur des données simulées: on utilise le modèle adapté au patient (\ie avec une géométrie du parenchyme et de l'arbre et un dépla\-cement de la surface du poumon provenant de données médicales) et on simu\-le des scénarios de respiration en faisant varier les résistances de l'arbre. Cela constitue la base d'apprentissage. Les données réelles de ventilation sont ensuite comparées à ces données simulées et une classification des branches pathologiques déduite.

La deuxième étape est d'estimer les pressions en bout d'arbre. Pour cela le problème de minimisation suivant est résolu\vspace*{-3pt}\enlargethispage{\baselineskip}%
\[
\min_{\substack{p=\sum_{i=1}^{N}p_i 1_{\Omega_i}\\ p=p_{i_0} \text{ dans } \Omega_{i_0}}} \|\rho_s\partial_{t}^2\dd_{\mes} -\div  \sigma(\dd_{\mes}) +\nabla p \|_{H{^{-1}}(\Omega)},
\]
où le domaine $i_0$ correspond à un domaine irrigué par un chemin sain dans l'arbre et $H^{-1}(\Omega)$ désigne le dual de l'espace $H^1_0(\Omega)$.

\Changelback

Ainsi on connaît une estimation de la pression en bout d'arbre, ce~qui permet de calculer les résistances identifiées comme pathologiques.
Cette méthodologie a été testée, mais malheureusement unique\-ment sur des données synthétiques bruitées, avec un taux de 85.2\% de détection des broncho-constrictions.

\vspace*{-.3\baselineskip}\vskip0pt
\section{Conclusion}
Nous avons ainsi présenté une série de modèles permettant de décrire la ventilation et d'explorer différents scénarios. Ces modèles sont bien sûr imparfaits, cependant leur étude soulève à la fois des questions mathématiques et numériques originales, et ils apportent également de nouveaux outils pour aider à la compréhension de la physiologique, ou même construire de nouvelles approches pour aider au diagnostic.

\vspace*{-.7\baselineskip}\vskip0pt
\backmatter
\bibliographystyle{jepalpha+eid}
\bibliography{xups22-01}
\end{document}
