Cédric Galusinski
email: galusins@univ-tln.fr

Paul Vigneaux
email: vigneaux@math.u-bordeaux1.fr

Level-Set method and stability condition for curvature-driven flows

We consider models for the simulation of curvature-driven incompressible bifluid flows, where the surface tension term is discretized explicitly. From this formulation a numerical stability condition arises for which we present a new theoretical estimation for low and medium Reynolds numbers. We illustrate our analysis with numerical simulations of microfluidic flows using Level Set method. Finally, we propose a method to reduce computational cost induced by this stability condition for low flow velocities.

Version française abrégée

Durant la dernière décennie, le champ de la microfluidique s'est drastiquement diversifié : il permet le transport de trés faibles volumes de fluide (nanolitres) à travers des canaux du diamètre d'un cheveu. Cette miniaturisation est de plus en plus utilisée pour construire un vaste panel d'outils dont, par exemple, l'utilisation de microgouttes comme des réacteurs chimiques. La simulation de ces écoulements soulève de nouvelles difficultés en mécanique des fluides numérique.

Dans cette Note, on s'intéresse en particulier à la simulation d'écoulements bifluides incompressibles, pilotés par la tension de surface. En suivant une approche classique pour la modélisation de tels écoulements, on utilise une formulation de suivi d'interface découplant le calcul de la vitesse de l'écoulement et le déplacement de l'interface [START_REF] Demello | DNA amplification moves on[END_REF] (3) (4). Ainsi, la discrétisation du terme de tension de surface est explicite. De cette formulation découle une condition de stabilité dont nous proposons ici une nouvelle estimation théorique, pour des écoulements dont le nombre de Reynolds est faible ou modéré : Proposition 0.1 Supposons que (2) (3) (4) soit discrétisé en temps à l'aide d'une discrétisation explicite du terme de tension de surface. Alors, pour une telle discrétisation en temps et quelque soit la discrétisation en espace, le schéma numérique associé est stable sous la condition

∆t ≤ min (∆t c , ∆t σ ) , avec ∆t c = c 1 u -1 L ∞ (Ω) ∆x et ∆t σ = c 2 η σ ∆x (1)
où ∆t est le pas de discrétisation en temps, ∆x est le pas de discrétisation en espace et c 1 , c 2 ne dépendent pas des données du problème.

Ensuite, à l'aide d'une formulation Level Set pour le suivi d'interface [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF], nous retrouvons numériquement cette condition de stabilité. Nous montrons de plus que son non-respect conduit à la création de courants parasites associés à des oscillations de l'interface.

Pour des écoulements « lents », cette condition de stabilité se révèle être trés restrictive. Nous proposons donc une méthode pour diminuer les temps de calculs. Ces écoulements induisent des formes d'interface stabilisées que nous exploitons pour introduire un référentiel (13) (14) où l'écoulement devient stationnaire. Les coûts de calculs sont alors drastiquement diminués car il ne reste plus qu'à translater les états stationnaires obtenus dans ce référentiel.

A model for bifluid flows

The field of microfluidics has dramatically improved over the last decade : it enables the transport of minuscule amounts of fluid (nanolitres) through channels the diameter of a single human hair. This miniaturization is increasingly employed to build a wide range of technological devices for example, genome sequencing [START_REF] Demello | DNA amplification moves on[END_REF] or using microdroplets as chemical reactors [START_REF] Joanicot | Droplet Control for Microfluidics[END_REF]. But a better understanding of the evolution of droplets as they move through microchannels is necessary in order to achieve the flow control needed in practical applications. Our development of a numerical model for simulating bifluid flows in microchannels leads us to the result of this Note : the derivation of a stability condition induced by surface tension when this term is discretized explicitly. Furthermore, it appears that in microchannels, due to low flow velocities, this stability condition is heavily restrictive ; we thus propose a splitting approach to relax the numerical cost.

We consider the following model based upon the Level Set method [5] [6] proposed by Sussman et al. [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] for incompressible two-phase flows :

ρ(u t + u.∇u) -div(2ηDu) + ∇p = σκδ(φ)n ∀(t, x) ∈ R + × Ω, (2)
div(u) = 0 ∀(t, x) ∈ R + × Ω, (3)
φ t + u.∇φ = 0 ∀(t, x) ∈ R + × Ω, ( 4
)
where Ω is the 2D (or 3D) bounded fluid domain, ∇ denotes the gradient operator with respect to the variable x, Du = (∇u + ∇ T u)/2, the density ρ is assumed constant, η is the variable viscosity such that : η = η(φ) = η i , i = 1, 2 in fluid i, σ is the surface tension coefficient, n is the normal to the interface, κ is the curvature of the interface and δ(φ) is the Dirac delta function localized on the interface. Different boundary conditions on the velocity u close the system.

Surface tension induced stability condition

Because the stability condition induced by surface tension, proposed by Brackbill et al. [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] for Navier-Stokes equation, involves the density but not the viscosity : this condition is only valid in certain flow regimes (high Reynolds numbers). It is therefore necessary to derive another stability condition, for both low and medium Reynolds numbers. In the following, we propose a new stability condition induced by surface tension for Navier-Stokes equation which is also the stability condition for Stokes equation. This stability condition is induced by the explicit discretization of the surface tension term and avoids oscillatory behaviour of the interface around an asymptotic shape of interface. These oscillations are related to what is known in the literature as parasitic currents [START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF].

Stability analysis

Proposition 2.1 Assume that (2) (3) (4) is discretized in time by an explicit discretization of the surface tension term. Then, a numerical scheme, associated to such a time discretization and all space discretizations, is stable under the condition

∆t ≤ min (∆t c , ∆t σ ) , with ∆t c = c 1 u -1 L ∞ (Ω) ∆x and ∆t σ = c 2 η σ ∆x ( 5
)
where ∆t is the time step, ∆x is space step of the discretization, and c 1 , c 2 do not depend on the data of the problem.

Proof : The first constraint ∆t ≤ c 1 u -1 L ∞ (Ω) ∆x is the classical CFL condition where c 1 depends on the choice of the scheme to discretize the transport equation ( 4). The second condition avoids oscillations phenomena of the interface. We perform the analysis on the continuous problem instead of the discrete problem. For that, we assume that the numerical scheme approximates consistently the continuous problem. We perform a stability analysis by an interface perturbation approach. The estimate of the velocity induced by the perturbation is the key point of the analysis.

We consider a perturbation of amplitude δ and wavelength L of the interface Γ 0 . Let us recall that in the Level Set framework, a regularization of Dirac delta function (∇H(φ)) is obtained by considering H ε instead of H, where H ε is a smooth increasing function and approximates the Heaviside function as ε goes to zero,

suppH ′ ε ⊂ (-ε; ε), H ′ ε ∞ ≤ 2 ε .
We then denote w the difference between the velocity induced by the unperturbed interface Γ 0 and the velocity induced by the perturbed interface, Γ. The perturbation of the surface tension term, noted f in the following, can be estimated, in L ∞ (Ω) norm, by :

f L ∞ (Ω) ≤ cσ δ L 3 + ∂ τ κ 0 L ∞ (Γ0) δ ε (6)
where c does not depend on the whole parameters, κ 0 is the curvature of the unperturbed problem and ∂ τ denotes the tangential derivative along Γ 0 . Assuming velocity perturbation responds, in main part, with the same frequency L (which is numerically verified), we have :

∇w L 2 (Ω) ∼ 1 L w L 2 (Ω) (7)
Therefore, plugging [START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF] in a classical energy estimate, one shows there exists a constant C such that :

ρ 2 d dt w 2 L 2 (Ω) + Cη L 2 -ρ ∇u L ∞ (Ω) w L 2 (Ω) ≤ L 2 Cη f 2 L 2 (Ω) (8)
and, if ρ ∇u L ∞ (Ω) ≤ Cη 2L 2 , which is true for smooth flow and sufficiently small wavelength L, we get :

w L 2 (Ω) ≤ L 2 Cη f L 2 (Ω) (9)
Invoking again the profile of the solution (as for assumption [START_REF] Shin | Accurate representation of surface tension using the level contour reconstruction method[END_REF]) and the surface tension term, we have :

w L ∞ (Ω) ≤ L 2 Cη f L ∞ (Ω) (10)
An explicit time discretization of the source term obliges to consider sufficiently small time step so that the displacement (during the time step) of the interface is smaller than the size δ of the perturbation of the interface. If the displacement is bigger than 2δ, perturbations are amplified and oscillates. Finally, we can write the stability condition by saying that the oscillations are removed if :

∆t σ = δ w L ∞ (Ω) ≤ c ′ η σ L 1 + ∂ τ κ 0 L ∞ (Γ0) L 3 ε , ( 11
)
where c ′ does not depend on the physical parameters. As the wavelength L is upper bounded and ε is chosen of size ∆x, this condition is restrictive for the smaller wavelength admissible in the numerical process. We are then concerned with L ∼ ∆x, and it reads :

∆t σ ∼ η σ ∆x. ( 12
)
Remark 1 If the flow velocity is very low then the previous condition on the time step is very restrictive since ∆t c ≫ ∆t σ ( e.g. in microfluidics). Remark 2 There are two ways to relax the stability constraint. First, as mentioned in the reference work of Brackbill et al., this constraint would be removed by an implicit treatment of surface tension term. A step in this direction has been done recently by Hysing [START_REF] Hysing | A new implicit surface tension implementation for interfacial flows[END_REF]. Second, smoothing of the curvature is another way to relax the constraint. Indeed, following the proof, a smoothing allows the artificial increase of the capillary time step ∆t σ .

In order to reduce the numerical cost, we propose another method, which relies on the use of the steady behaviour of the interface in some zones of channels. We will detail this method in the following section.

Numerical results

We focus on numerical simulation of microfluidic flows with very low Reynolds number.

Numerical stability condition

In practical applications, we need to determine the constant c 2 . Running numerical tests with a Stokes code and taking a constant η, we find that c 2 = 6 is the critical value for which the computation is always stable when we refine the mesh. These numerical experiments show that c 2 = 6 is independent of the space step ∆x, and c 2 = 8 is stable except for too fine grids. Figure 1 shows numerical simulations of an unconfined droplet in a microchannel with two values of c 2 .

Figure 1. on the left (velocity field in blue, interface in black and some streamlines in arrowed-line), c 2 = 8 and the computation is stable for a medium grid, whereas on the right (velocity field in red, interface in black), c 2 = 9 and the instability appears on the same grid ; velocity fields are those in the drop frame of reference

Mixer analysis: The drop frame of reference

We want to perform numerical simulations of microdroplets to study mixing inside them by observing the streamlines in the drop frame of reference.

In straight microchannels the shape of the droplets is stabilized and we have to exploit this stationary phenomena in order to reduce the numerical cost. To this end, we work in the drop frame of reference.

Find the scalar droplet's velocity

In order to work in the drop frame of reference, we have to define the scalar droplet's velocity. This velocity has a meaning as soon as the shape of the droplet is stabilized and then there exists a scalar u d such that :

u = u d U + v, v • n = 0 on Γ 0 , (13)
where U is a "unitary" flow parallel to the wall in a regular channel. Then the droplet moves with the global velocity u d U and the velocity in the drop frame of reference is v. In the general case, with stabilized or destabilized droplet shape, we define the local droplet's velocity on the fluid interface where U • n = 0 as

u loc d = u • n U • n on Γ 0 . (14)
We then define the droplet's velocity as the mean value of local droplet velocities.

Splitting of the flow

In straight microchannels, the velocity of a stabilized droplet verifies (13). We reduce the computational cost by using the algorithm which consists in translating the droplet while v • n remains close to zero. This is obtained with a standard CFL condition on the time step with velocity u d U . To sum up, the algorithm consists in two steps : an iterative step which allows to compute the shape modification of the interface due to surface tension (i.e. by transporting φ with v using the restricting time step ∆t σ ) and then, a step which consists in transporting φ with u d U using ∆t c . Note that this second step can be deleted if we only want to get the asymptotic shape of the droplet and not its position in the channel. Comparison of CPU time between a standard method and our splitting method depends on the length of the channel. We decrease heavily CPU time when channel's length is many times droplet's length and for low flows. More precisely, there is no theoretical saving when droplet had not reach its asymptotical shape but when droplet's shape becomes stationary and the droplet is only translating in the channel, time saving is of a factor ∆t σ /∆t c (which, in our microfluidic applications, can be of order 100). Figures 1 and2 show that confinement and injection velocity (among others, which will be presented in a subsequent publication) allow to induce various mixing regimes inside droplets.

Figure 2 .

 2 Figure2. Various regimes for a confined droplet in a microchannel (wall in grey, velocity field in blue, interface in black and some streamlines in arrowed-line) : on the left, injection velocity is lower than on the right.

Acknowledgements

This work was partially supported by ACI-NIM "Calculs de micro-fluides"(CNRS).


Céline Grandmont

Bertrand Maury

Assia Soualah

Multiscale modelling of the respiratory track : a theoretical framework

come

Résumé. Afin de décrire l'écoulement de l'air dans les voies aériennes supérieures et proximales, nous proposons un modéle multiéchelles basé sur la décomposition en trois régions de l'arbre trachéobronchique. Dans chacune de ces régions les comportements mécaniques sont différents. Les équations de Navier-Stokes décrivent l'écoulement de l'air dans la partie supérieure et elles sont couplées à une équation différentielle ordinaire modélisant le mouvement du diaphragme et du parenchyme pulmonaire. Pour ce systeme couplé nous démontrons l'existence de solutions faibles en temps petit pour données quelconques, ainsi que l'existence globale en temps dans le cas, très particulier, où la donnée initiale est suffisamment petite et où aucune force extérieure n'est appliquée au système. De plus, si la raideur du ressort modélisant l'élasticité du poumon est nulle, nous obtenons un résultat d'existence globale pour des données petites (conditions initiales et forces appliquées. Nous illustrons le comportement du modèle par des simulations numériques bi-dimensionnelles.

Introduction, modelling aspects

Breathing involves gas transport through the respiratory track with its visible ends, nose and mouth. Air then streams from the pharynx down to the trachea. The trachea extends from the neck into the thorax, where it divides into right and left main bronchi, which enter the corresponding lungs. The inhaled air is then convected in the bronchus tree which ends in the alveoli embedded in a viscoelastic tissue, made in particular of blood capillaries, and where gaseous exchange occurs. Each lung is enclosed in a space bounded below by the diaphragm and laterally by the chest wall. The air movement is achieved by the displacement of the diaphragm and parenchyma tissue. But the lung may fail to maintain an adequate supply of air. Accidental inhalation of liquid or solid, asthma crisis, pathologies changing the elastic behaviour of the parenchyma (like emphysema or fibrosis), obstructed nose, cancer may occur. The lung is then a place where many exchanges and interactions take place, and a numerical lung would certainly be very helpful in the understanding of some of the deseases and a way to guide the intuition for curative gestures. At the time being, the complex fractal geometry of the airway tree makes the air flow simulation on the whole tree unreachable. Note, moreover, that the distal airways from generation 7 cannot be visualised by common medical imaging techniques. Consequently, it is necessary to find new efficient strategies, including simple but realistic models. One possible choice is to try to describe the evolution of the air flux by a simple ODE as it is done in [START_REF] Martin | Modelling of respiratory system mechanics involving gas exchange in the human lungs[END_REF]. But even if the model can give valuable hints to understand the respiration mechanisms it can not provide precise informations on the full 3D flow. Our aim is to obtain a model that describes accurately the air flow in the proximal part. To achieve this, one can not forget that this flow is dependent of the distal part and driven by the motion of the diaphragm and parenchyma. One solution will be to find physiological boundary conditions. Yet no such "in vivo" pressure or velocity measurement is available. Thus our aim is to obtain a simplified description of the distal part and what we propose is a decomposition of the respiratory tree into three stages where different models will be exploited and in which the mechanical behaviour is quite different:

• the upper part (up to the 6th generation), where the Navier-Stokes equations hold to describe the fluid flow, • the distal part (from the 7th to the 17th generation), where one can assume that the Poiseuille law is satisfied in each bronchiole, • the acini, where the oxygen diffusion takes place and which are embedded in an elastic medium, the parenchyma. We will assume that the pressure is uniform in the acini part and that they are embedded in a box representing the parenchyma. The motion of the diaphragm and the parenchyma is described by a simple spring model. The decomposition can be schematized by the following figure:

Ω Γ l Γ 0 Γ i Ω i Figure 1. Multiscale model
The"outlet" Γ i , 1 ≤ i ≤ N of the upper part are coupled with Poiseuille flows themselves coupled with the spring motion. We will assume that the velocity profile on the "outlets" Γ i , 1 ≤ i ≤ N is given and, for instance, is a parabolic profile. The assumption is reasonnable since we assume that the Poiseuille law holds true from the 7th to the 17th generation of the bronchial tree. Note that the same type of multiscale modelling has been investigated for different applications such as blood flow simulations ( [START_REF] Quarteroni | Analysis of a Geometrical Multiscale Model Based on the Coupling of ODEs and PDEs for Blood Flow Simulations[END_REF], [START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modelling of blood and pressure in arteries[END_REF]) or air flow simulation ( [START_REF] Grandmont | A multiscale / multimodel approach of the respiration tree[END_REF]).

The paper is organized as follows: in a first part we present the coupled system and its variational formulation then we derive, at least formally, a priori error estimates. The main difficulty is to estimate the nonlinear terms that is to say the flux at the artificial boundaries Γ i of system kinetic energy. In a second part we prove, thanks to a Galerkin approach, that there exists at least one solution, locally in time. The proof is based on standard arguments: energy estimates, obtention of additional bounds in order to pass to the limit in the nonlinear convection terms. Finally we give 2D numerical evidence that this model can reproduce some aspects of normal or pathological breathing.

Problem setting

In the upper part, denoted by Ω we assume that the Navier-Stokes equations hold:

                                   ρ ∂u ∂t + ρ(u • ∇)u -µ△u + ∇p = 0, in Ω , ∇ • u = 0, in Ω , u = 0, on Γ l , u = λ i U i on Γ i , i = 0, . . . , N , Γ0 (µ∇u • n -pn) • U 0 = -P 0 Γ0 U 0 • n on Γ 0 , Γi (µ∇u • n -pn) • U i = -Π i Γi U i • n on Γ i i = 1, . . . , N, (1)
where u and p are respectivelly the fluid velocity and the fluid pressure. On the lateral trachea boundary Γ l we impose non-slip boundary conditions on the velocity, whereas on the artificial boundary Γ i , 0 ≤ i ≤ N we consider two types of boundary conditions. First we assume that, on each outlet, the trace of the fluid velocity is proportional to a given velocity profile U i . The proportionality coefficient λ i mesures the velocity flux at the interface and is an unknown of the problem. The second type of boundary condition corresponds to a pressure force exerted on the boundary. The pressure P 0 is given whereas the pressures Π i are unknown that depend on the dowstream parts. Each of the Ωi should be a dyadic tree in which we assume that the flow is laminar. Thus, proceeding as in [START_REF] Maury | Outlet Dissipative conditions for air flow in the bronchial tree[END_REF] and [START_REF] Grandmont | A viscoelastic model with non-local damping application to the human lungs[END_REF] and by analogy with an electric network, we can consider that the flow is caracterized by a unique equivalent resistance that depends on each resistance of the local branches. Thus, each of the Ωi is considered as a cylindrical domain, where the flow satisfies a Poiseuille law:

Π i -P i = R i Γi u • n, R i ≥ 0, (2)
where R i denotes the equivalent resistance of the distal tree and P i is an alveola pressure. Note that R i depends on the geometric properties (length and diameter) of all the branches of the i-th subtree. Moreover, we assume that all the alveola pressures are equal: P i = P . Finally, we suppose that these alveoli are embedded in a box filled with an incompressible medium that represents the parenchyma. One part of the box is connected to a spring that governs the diaphragm and parenechyma motion. The equation satisfied by the position x of the diaphragm writes:

mẍ = -kx + f ext + f P , (3)
where m is the total mass of the lung, k is the stiffness of the spring (that characterises the elastic behavior of the lung) and f ext is the force developped by the diaphragm during inspiration and forced expiration. In order to couple this simple ODE to the upper part of the model, we have to define f P that stands for the pressure force applied by the flow on the elastic medium. If we denote by S the surface of the moving boundary box (diaphragm surface), we have

f P = P S. (4)
Moreover, since the flow is incompressible and since we assume that the parenchyma is made of an incompressible medium, the flow volume variation is equal to the volume variation of the parenchyma box, thus we have

S ẋ = N i=1 Γi u • n = - Γ0 u • n = -λ 0 Γ0 U 0 • n. (5)
Thus the coupled problem can be written as follows:

                                             ρ ∂u ∂t + ρ(u • ∇)u -µ△u + ∇p = 0 , in (0, T ) × Ω , ∇ • u = 0 , in (0, T ) × Ω , u = 0 , on (0, T ) × Γ l , u = λ i U i , on (0, T ) × Γ i , i = 0 , . . . , N , Γ0 (µ∇u • n -pn) • U 0 = -P 0 Γ0 U 0 • n , on (0, T ) × Γ 0 , Γi (µ∇u • n -pn) • U i = -P Γi U i • n -R i Γi U i • n , on (0, T ) × Γ i , i = 1 , . . . , N , mẍ + kx = f ext + SP , S ẋ = -λ 0 Γ0 U 0 • n. (6)
This system of equations have to be completed by initial conditions (u, x, ẋ)| t=0 = (u 0 , x 0 , x 1 ), with Sx 1 = -

Γ0 u 0 • n. (7)
One particularity of this system is that all the outlets Γ i , 1 ≤ i ≤ N are coupled. This is not the case, for instance, in [START_REF] Quarteroni | Analysis of a Geometrical Multiscale Model Based on the Coupling of ODEs and PDEs for Blood Flow Simulations[END_REF] or [START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modelling of blood and pressure in arteries[END_REF] where the same type of multiscale modelling is performed but for blood flow simulations. Note that the elastic behavior of the lung is described by only one degree of freedom. Moreover, in the whole coupled model we have only few parameters to fit: m, k, S, f ext and the resistances R i . In particular by modifying k and R i one could obtain pathological behaviors such as asthma (increase of the resistances) or emphysema (decrease of k). Nevertheless the considered spring model is a very simple one and some aspects of the respiratory cycle can not be reproduced by such a simple model. We refer to [START_REF] Martin | Modelling of respiratory system mechanics involving gas exchange in the human lungs[END_REF] for a more complete spring model.

A priori estimates

We start to derive, at least formally, a priori estimates for the coupled system. We multiply the Navier-Stokes equation by u and integrate over Ω, by taking into account the profile constraint u = λ i U i on Γ i :

ρ 2 d dt Ω |u| 2 + ρ Ω (u • ∇u)u + µ Ω |∇u| 2 + N i=1 R i (λ i ) 2 Γi U i • n 2 + P 0 λ 0 Γ0 U 0 • n +P N i=1 λ i Γi U i • n = 0,
which can be written, by using the fluid incompressibility

ρ 2 d dt Ω |u| 2 + ρ Ω (u • ∇u)u + µ Ω |∇u| 2 + N i=1 R i (λ i ) 2 Γi U i • n 2 + P 0 λ 0 Γ0 U 0 • n -P λ 0 Γ0 U 0 • n = 0.
Moreover we can write the convective terms, thanks to an integration by parts, as a flux of kinetic energy at the inlet and at the outlets:

ρ 2 d dt Ω |u| 2 + ρ 2 N i=0 Γi |u| 2 (u • n) + µ Ω |∇u| 2 + N i=1 R i (λ i ) 2 Γi U i • n 2 + P 0 λ 0 Γ0 U 0 • n -P λ 0 Γ0 U 0 • n = 0. ( 8
)
Then we multiply the spring equation by ẋ

m 2 d dt | ẋ| 2 + k 2 d dt |x| 2 = f ext ẋ + P S ẋ. (9)
Now using the fact that S ẋ = -λ 0 Γ0 U 0 • n and adding ( 8), ( 9) we obtain

ρ 2 d dt Ω |u| 2 + ρ 2 N i=0 Γi |u| 2 (u • n) + µ Ω |∇u| 2 + N i=1 R i (λ i ) 2 Γi U i • n 2 + m 2 d dt | ẋ| 2 + k 2 d dt |x| 2 = P 0 S ẋ + f ext ẋ. ( 10
)
Identity [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] represents the energy balance of the coupled system. In particular, the energy is dissipated: as for the standard Navier-Stokes equations the fact that the flow is viscous contribute to the dissipation of the energy but here there is a second contribution to the dissipation that comes from the resistive part of the bronchial tree, namely

N i=1 R i (λ i ) 2 Γi U i • n 2 = N i=1 R i Γi u • n 2 .
In order to obtain energy estimates we have to estimate this flux of kinetic energy at the interface, that writes, taking into account the profile constraint at the interface

(u = λ i U i on Γ i ) N i=0 Γi |u| 2 (u • n) = N i=0 (λ i ) 3 Γi |U i | 2 (U i • n).
First of all we introduce the following functional spaces:

H = {v ∈ L 2 (Ω) , ∇ • v = 0, v • n = 0 on Γ l }, V = {v ∈ H 1 (Ω) , ∇ • v = 0, v = 0 on Γ l }, W = {v ∈ V , ∃ Λ = ( λi ) 0≤i≤N s. t. v = λi U i on Γ i , i = 0, . . . , N }.
To bound the nonlinear convection term we will use the following lemma:

Lemma 1.1. If U i • n = 0 and if U i ∈ (H 1 2 00 (Γ i )) ′ there exists a constant C i > 0 such that, for all v ∈ W , we have | λi | ≤ C i v L 2 (Ω) , i = 0 , . . . , N.
Proof. It is an easy task to verify, using standard arguments (see for instance [START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF]) that

v • n (H 1/2 00 (Γi)) ′ ≤ C v L 2 (Ω) + ∇ • v L 2 (Ω) , ∀v ∈ H div (Ω) = {v ∈ L 2 (Ω); ∇ • v ∈ L 2 (Ω)}.
Here the space H 1/2 00 (Γ i ) denotes the space of trace function that, if extended by zero over ∂Ω, belongs to H 1/2 (∂Ω). But for all v ∈ W we have

| λi | = v • n (H 1/2 00 (Γi)) ′ U i • n (H 1/2 00 (Γi)) ′ .
Then taking into account the fact that ∇ • v = 0, we have

| λi | ≤ C i v L 2 (Ω) where C i = C U i • n (H 1/2 00 (Γi)) ′ . Lemma 1.1 enables us to estimate the nonlinear convective term N i=0 Γi |u| 2 (u • n) = N i=0 (λ i ) 3 Γi |U i | 2 (U i • n) ≤ ΓC 1 u 3 L 2 (Ω) , (11) with Γ
= N i=0 Γi |U i | 2 |U i • n and C 1 = max 1≤i≤N |C i |.
Note that this estimate relies on the fact that the trace of the fluid velocity on each Γ i is supposed to be described by one parameter. This estimate is, of course, not true in the general case. We will assume, to simplify, that the velocity profils U i and the resistances R i do not depend on time. Nevertheless the energy estimates derived below hold also if

U i ∈ L ∞ (0, T ; H 1 2 00 (Γ i ))) and R i ∈ L ∞ (0, T ).
Using the last estimate and the energy balance (10) yields

ρ 2 d dt Ω |u| 2 + m 2 d dt | ẋ| 2 + k 2 d dt |x| 2 +µ Ω |∇u| 2 + N i=1 R i (λ i ) 2 Γi U i • n 2 ≤ (|SP 0 | + |f ext |)| ẋ| + ρ 2 ΓC 1 u 3 L 2 (Ω) ≤ C(| ẋ| 2 + u 3 L 2 (Ω) + |SP 0 | 2 + |f ext | 2 ). (12)
Thus if P 0 ∈ L 2 (0, T ) and f ext ∈ L 2 (0, T ), we easily verify that there exist a time τ 0 , 0 < τ 0 < T sufficiently small depending on the data and a constant C > 0 also depending on the data such that for any regular solution (u, x):

ρ 2 u 2 L 2 (Ω) (t) + m 2 | ẋ| 2 (t) + k 2 |x| 2 (t) +µ t 0 Ω |∇u| 2 + N i=1 R i t 0 (λ i ) 2 Γi U i • n 2 ≤ C, t ≤ τ 0 . (13)
Thus any regular solution u is bounded by the data in L ∞ (0, τ 0 ; L 2 (Ω)) ∩ L 2 (0, τ 0 ; H 1 (Ω)) that implies in particular that the velocity fluxes and ẋ are bounded in L ∞ (0, τ 0 ).

Note that, due to the bound obtained on the nonlinear convective term, this a priori estimate only holds on a small time interval and will lead to an existence result for large data locally in time.

Nevertheless, if the initial data are small enough and in the special case where no external forces are applied we can prove that the velocity remains in a given ball for any time. Indeed, using Poincaré inequality,

u L 2 (Ω) ≤ C 2 ∇u L 2 (Ω) , (14)
the nonlinear terms can be bounded, as follows:

N i=0 Γi |u| 2 (u • n) ≤ ΓC 1 C 2 2 u L 2 (Ω) ∇u 2 L 2 (Ω) . (15)
In what follows we set K = ΓC 1 C 2 2 . and consequently ρ 2

d dt u 2 L 2 (Ω) + m 2 d dt | ẋ| 2 + k 2 d dt |x| 2 + µ 2 - Kρ 2 u L 2 (Ω) ∇u 2 L 2 (Ω) + N i=1 R i (λ i ) 2 Γi U i • n 2 ≤ 0 (16)
Thus, assuming that u 0 L 2 (Ω) < µ Kρ , we obtain that on any time interval during which

µ 2 - Kρ 2 u L 2 (Ω) > 0 ρ 2 d dt u 2 L 2 (Ω) + m 2 d dt | ẋ| 2 + k 2 d dt |x| 2 ≤ 0. Moreover, if ρ 2 u 0 2 L 2 (Ω) + m 2 |x 1 | 2 + k 2 |x 0 | 2 ≤ µ 2 8ρK 2 , the solution verifies ρ 2 u 2 L 2 (Ω) + m 2 | ẋ| 2 + k 2 |x| 2 ≤ µ 2 8ρK
2 for all time. Remark 1.2. In the special case P 0 = 0, f ext = 0 and under the hypothesis u 0 L 2 (Ω) ≤ µ/ρK, the energy of the coupled system is decreasing. Note that this property cannot be expected in general, because an initial velocity with a large L 3 norm on the inlet is likely to induce, at least transitorily, an increase of global energy by kinetic energy entering the domain. In a similar manner, if the profiles are not prescribed, a highly singular initial velocity on the inlet may produce a similar increase of energy, as small at it may be in L 2 volumic norm.

One can also obtain estimate globally in time in the special case where k = 0. Indeed, in this case we have

ρ 2 d dt u 2 L 2 (Ω) + m 2 d dt | ẋ| 2 +µ ∇u 2 L 2 (Ω) + N i=1 R i (λ i ) 2 Γi U i • n 2 ≤ (|SP 0 | + |f ext |)| ẋ| + ρ 2 ΓC 1 u 3 L 2 (Ω) .
Using Poincaré inequality and since

| ẋ| = λ 0 Γ0 u 0 • n S ≤ C|λ 0 | ≤ C 3 u L 2 (Ω) ≤ C 3 C 2 ∇u L 2 (Ω) ,
we obtain

ρ 2 d dt u 2 L 2 (Ω) + m 2 d dt | ẋ| 2 + µ 2 - Kρ 2 u L 2 (Ω) ∇u 2 L 2 (Ω) + N i=1 R i (λ i ) 2 Γi U i • n 2 ≤ C 2 2 C 2 3 2µ (|SP 0 | 2 + |f ext | 2 ),
and using, once again, Poincaré inequality

ρ 2 d dt u 2 L 2 (Ω) + m 2 d dt | ẋ| 2 + 1 C 2 2 µ 2 - Kρ 2 u L 2 (Ω) u 2 L 2 (Ω) + N i=1 R i (λ i ) 2 Γi U i • n 2 ≤ C 2 2 C 2 3 2µ (|SP 0 | 2 + |f ext | 2 ). Moreover, since | ẋ| 2 ≤ C u 2 L 2 (Ω)
it is easy to obtain an inequality of the type

d dt φ + A( µ 2 -B φ)φ ≤ C 2µ (|SP 0 | 2 + |f ext | 2 ), with φ = | ẋ| 2 + u 2 L 2 (Ω) . Consequently, assuming that φ(0) ≤ µ 4B and that |SP 0 | 2 + |f ext | 2 ≤ Aµ 4 32B 2 C
, it is easy to verify that φ(t) ≤ µ 4B for all t.

Remark 1.3. Note that the case k = 0 can not be treated in the same way because of the diaphragm displacement that can not be bounded easily in L ∞ . Observe, however that from the physical point of view this displacement is bounded (the volume of the chest is given) and one could impose it by considering a more general spring model (see [START_REF] Martin | Modelling of respiratory system mechanics involving gas exchange in the human lungs[END_REF]).

Remark 1.4. The previous assumptions on the data involve µ. Note however that one can obtain conditions involving also the resistances R i (provided that the R i are strictly positive and the flux of the velocity profil Γi U i • n are not equal to zero) by estimating the nonlinear convection term as follows:

N i=0 Γi |u| 2 (u • n) ≤ C 3 u L 2 (Ω) N i=0 |λ i | 2 ,
and by writting λ 0 (or ẋ) as a linear combinaison of the (λ i ) 1≤i≤N . Yet the resulting conditions are of the same type than the one we just obtained, since the resistances R i are proportional to µ.

Variational formulation

In this subsection we give the variational formulation of the coupled system. But, first we precise the assumptions made on the data

(u 0 , x 0 , x 1 ) in H × R 2 with Sx 1 = - Γ0 u 0 • n, P 0 ∈ L 2 (0, T ), f ext ∈ L 2 (0, T ), R i ≥ 0 , i = 1 , . . . , N, U i ∈ H 1 2 00 (Γ i ) , i = 0 , . . . , N. (17)
Let us consider v ∈ {w ∈ H 1 (Ω) , w = 0 on Γ l and ∃(β i ) 0≤i≤N s. t. w| Γi = β i U i } and φ ∈ D([0, t)). In order to obtain a weak formulation of the coupled system we multiply the Navier-Stokes equations by vφ and integrate over Ω × (0, T ). The spring equation is multiplied by -

                                           ρ Ω u(t) • vφ(t) -ρ t 0 Ω u • v φ + ρ t 0 Ω (u • ∇u) • v φ + µ t 0 Ω ∇u : ∇v φ + N i=1 R i t 0 Γi u • n Γi v • n φ - t 0 Ω p∇ • u φ - t 0 P ∂Ω v • n φ - m S 2 t 0 Γ0 u • n Γ0 v • n φ + k S 2 t 0 s 0 Γ0 u • n Γ0 v • n φ = - t 0 P 0 Γ0 v • n φ - t 0 f ext S Γ0 v • n φ - k S x 0 Γ0 v • n t 0 φ -ρ Ω u 0 • v φ(0) - m S x 1 Γ0 v • n φ(0), ∀v ∈ H 1 (Ω) , ∇ • v = 0 , v = 0 on Γ l and ∃( λi ) 0≤i≤N s. t. v| Γi = λi U i . (18)
Next considering test functions that are divergence free, we obtain a second weak formulation of the coupled problem, and we shall say that u is a weak solution of ( 6

) on [0, T ) if -u ∈ L 2 (0, T ; V ) ∩ L ∞ (0, T ; H), -∃λ i ∈ L ∞ (0, T ) such that u i | Γi = λ i U i , for a. e. t, -∀v ∈ H 1 (Ω) , ∇ • v = 0 , v = 0 on Γ l and ∃( λi ) 0≤i≤N s. t. v| Γi = λi U i , ∀φ ∈ D([0, T ))                                        ρ Ω u(t) • vφ(t) -ρ t 0 Ω u • v φ + ρ t 0 Ω (u • ∇u) • v φ + µ t 0 Ω ∇u : ∇v φ + N i=1 R i t 0 Γi u • n Γi v • n φ + m S 2 Γ0 u(t) • n Γ0 v • n φ(t) - m S 2 t 0 Γ0 u • n Γ0 v • n φ + k S 2 t 0 s 0 Γ0 u • n Γ0 v • n φ = - t 0 P 0 Γ0 v • n φ - t 0 f ext S Γ0 v • n φ - k S x 0 Γ0 v • n t 0 φ -ρ Ω u 0 • v φ(0) - m S x 1 Γ0 v • n φ(0). (19)
Note that here we have expressed all the quantities with the help of the fluid velocity. The velocity of the spring can be simply recovered thanks to the equality S ẋ = -

Γ0 u • n = -λ 0 Γ0 U 0 • n.

Main Results

One aim of the paper is to prove the existence of weak solutions as defined in the previous section. We have two types of existence results: one locally in time for large data, and one globally in time for small enough initial data in the special case where no external forces are applied to the coupled system. Note moreover that in the case where the stifness of the spring is equal to zero, one could easily derive an existence result, globally in time, provided that the data (initial data and applied forces) are small enough. More precisely Theorem 1.5. Under the assumptions [START_REF] Similowski | Two-compartment modelling of respiratory system mechanics at low frequencies: gas redistribution or tissue rheology?[END_REF] there exist a time τ 0 , 0 < τ 0 < T sufficiently small and a weak solution u of ( 19) on [0, τ 0 ]. Moreover, there exists C such that the solution satisfies the following estimate:

ρ u 2 L ∞ (0,τ0;L 2 (Ω)) + m S 2 Γ0 u • n 2 L ∞ (0,τ0) + k S 2 t 0 Γ0 u • n 2 L ∞ (0,τ0) +µ u 2 L 2 (0,τ0;H 1 (Ω)) + N i=1 R i τ0 0 Γi u • n 2 ≤ C. (20)
Furthermore, in the special case where the external forces are equal to zero and if we assume that

ρ 2 u 0 2 L 2 (Ω) + m 2 |x 1 | 2 + k 2 |x 0 | 2 ≤ µ 2 8ρK 2 , ( 21
)
where K has been defined in Section 1.1 then there exists at least one weak solution of ( 6) on [0, ∞). Furthermore this solution satisfies:

ρ 2 u 2 L ∞ (0,∞;L 2 (Ω) + m 2S 2 Γ0 u • n 2 L ∞ (0,∞) + k 2 x 0 + 1 S t 0 Γ0 u • n 2 L ∞ (0,∞) ≤ µ 2 8ρK 2 , ( 22
)
and for any

T ρ u 2 L ∞ (0,T ;L 2 (Ω)) + m S 2 Γ0 u • n 2 L ∞ (0,T ) + k S 2 t 0 Γ0 u • n 2 L ∞ (0,T ) +µ u 2 L 2 (0,T ;H 1 (Ω)) + N i=1 R i T 0 Γi u • n 2 ≤ C T . (23)
Remark 1.6. Note that in the case where k = 0 and if the data are small enough one an obtain an existence result globally in time.

Proof of Theorem

The proof of the existence results is standard. In a first step we build a sequence of approximated solutions thanks to a Galerkin method. Then, since the energy bounds are not sufficient to pass to the limit in the weak formulations, we derive additional bounds. Applying a compactness result we then are able to pass to the limit in the system satisfied by the approximated solution leading to the existence of at least one weak solution.

Galerkin Method

The first step is to build a Galerkin basis of the space

W = {v ∈ V , ∃ Λ = ( λi ) 0≤i≤N s. t. v = λi U i on Γ i , i = 0, . . . , N }.
Let w belong to W . We remark first that, taking into account the incompressibility constrain, it implies that

N i=0 α i Γi U i • n = 0 with w| Γi = α i U i . Consequently if Γi U i • n = 0 for all 0 ≤ i ≤ N , the (α i ) belong to a
vector space of dimension N + 1 and if there exists i 0 such that Γi0 U i0 • n = 0 the (α i ) belong to a vectorial space of dimension N . Then, w = z + i∈J α i w i , where

z ∈ {v ∈ H 1 0 (Ω), ∇ • v = 0},
and where J and the w i ∈ V are defined as follows:

-

If Γi U i • n = 0 for all 0 ≤ i ≤ N , then J = {i, 0 ≤ i ≤ N } and            -µ△w i + ∇p i = 0 in Ω ∇ • w i = 0 in Ω w i = 0 on Γ l w i = 0 on Γ j j = i w i = U i onΓ i
Note that this case is improbable in the context of respiration.

-If there exists i 0 such that Γi0 U i0 • n = 0, then J = {i, 0 ≤ i ≤ N, i = i 0 } and we define the functions w i , i = i 0 as follows:

-

for i = i 0                -µ△w i + ∇p i = 0 in Ω ∇ • w i = 0 in Ω w i = 0 on Γ l w i = 0 on Γ j j = i, i 0 w i = β i U i0 onΓ i0 w i = U i on Γ i ,
where

β i = -Γi U i • n Γi 0 U i0 • n . Let (z n ) n∈N be a Galerkin basis of {v ∈ H 1 0 (Ω), ∇ • v = 0}. A Galerkin basis of W is then (z n ) n∈N ∪ (w i ) 0≤i≤N . For all m ∈ N, we define u m (t) ∈ W m = span(z n ) 1≤n≤m ⊕ span(w i ) 1≤i≤N solution of the discrete problem (P m )                                    ρ Ω ∂u m ∂t • v + ρ Ω (u m • ∇u m ) • v + µ Ω ∇u m : ∇v + N i=1 R i Γi u m • n Γi v • n + m S 2 d dt Γ0 u m • n Γ0 v • n + k S 2 t 0 Γ0 u m • n Γ0 v • n = -P 0 Γ0 v • n - f ext S Γ0 v • n - k S 2 x 0 Γ0 v • n, ∀v ∈ W m , u m (0) = u 0m , (24)
where u 0m ∈ W m and is such that u 0m converges to u 0 strongly in H as m goes to infinity.

Lemma 2.1. The discrete problem (P m ) has a unique solution u m ∈ W 1,∞ (0, T ; W m ).

Moreover the following energy estimates are satisfied:

∃τ 0 independent of m, s. t. u m L ∞ (0,τ0;L 2 (Ω)) + u m L 2 (0,τ0;H 1 (Ω)) ≤ C,
where C is a constant that does depends on m. Furthermore, if the initial data satisfy [START_REF] West | Respiratory Physiology -The Essentials[END_REF], and if P 0 = 0 and f ext = 0 then

ρ 2 u m 2 L ∞ (0,∞;L 2 (Ω)) + m 2S 2 Γ0 u m • n 2 L ∞ (0,∞) + k 2 x 0 + 1 S t 0 Γ0 u m • n 2 L ∞ (0,∞) ≤ µ 2 8ρK 2 ,
and for any T

ρ u m 2 L ∞ (0,T ;L 2 (Ω)) + m S 2 Γ0 u m • n 2 L ∞ (0,T ) + k S 2 t 0 Γ0 u m • n 2 L ∞ (0,T ) +µ u m 2 L 2 (0,T ;H 1 (Ω)) + N i=1 R i T 0 Γi u m • n 2 ≤ C T .
Proof. It is clear that there exists a time τ m such that the discrete problem (P m ) has a unique solution

u m ∈ W 1,∞ (0, τ m ; W m ).
In particular the mass matrix defined by

M ij = ρ Ω φ i • φ j + m S 2 Γ0 φ i • n Γ0 φ j • n with φ i , φ j ∈ W m ,
is invertible. Moreover this solution satisfies the estimates (13) derived at Section 1.1 and consequently τ m = τ 0 . Finally, by applying the same argument as the one we used at Section 1.1 we obtain the other desired bounds.

At least for a subsequence of u m (still denoted by u m ), these estimates yield the following weak convergences, as m goes to ∞:

   u m ⇀ u weakly in L 2 (0, T ; V ) u m ⇀ u weakly ⋆ in L ∞ (0, T ; H), (25)
where T = τ 0 in the general case and T can be any time in the case of small enough data. Moreover the limit u satisfies the same energy estimate as u m thus inequality ( 20) is satisfied as well as ( 22), (23) in the case where no external forces are applied and ( 21) is verified. Nevertheless the previous convergences are not sufficient to pass to the limit in the weak formulation (P m ) and in particular in the convection term. Consequently we need to derive additional bounds in order to obtain compactness of the sequence of approximated solution (u m ) m .

Additional bounds

We will follow exactly the same lines as in [START_REF] Girault | Finite Element Approximation of the Navier-Stokes Equations[END_REF][START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]. First we define

H γ (R; B 0 , B 1 ) = {v, v ∈ L 2 (R, B 0 ), D γ t v ∈ L 2 (R, B 1 )},
with D γ t v is the fractionnal derivative of order γ in time of v, defined as the inverse Fourier transform of (2iπτ ) γ û(τ ). The space H γ (R; B 0 , B 1 ) is endowed with the norm

u H γ (R;B0,B1) = ( u 2 L 2 (R;B0) + |τ | γ û 2 L 2 (R;B1) ) 1 2
.

We define also H γ (0, T ; B 0 , B 1 ), as the space of function obtained by restriction to [0, T ] of functions of H γ (R; B 0 , B 1 ). We will apply the following lemma that can be found also in [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]:

Lemma 2.2. Let B 0 , B, and B 1 be three Hilbert spaces such that B 0 ⊂ B ⊂ B 1 and such that B 0 is compactly embedded in B. Then for all γ > 0, the injection H γ (0, T ; B 0 , B 1 ) → L 2 (0, T ; B) is compact.

We are going to apply this lemma with

B 0 = V, B 1 = H, γ = 1 4 -ε and B = H
Indeed the additional estimate holds true:

Lemma 2.3. The sequence u m is bounded in H γ (0, T ; V, H) for 0 ≤ γ ≤ 1 4 -ε.
Proof. We denote by ūm the extention of u m by zero 0 for t < 0 and t > T , and ûm the Fourier transform in time of ūm . Since we already know that u m is uniformly bounded in L 2 (0, T, V ), it remains to prove that

+∞ -∞ |τ | 2γ ûm 2 L 2 (Ω) ≤ C.
Then applying Lemma 2.2 we have that u m is compact in L 2 (0, T ; H).

We have that ūm satisties

ρ Ω ∂ ūm ∂t • v + ρ Ω (ū m • ∇ū m ) • v + µ Ω ∇ū m : ∇v + N i=1 R i Γi ūm • n Γi v • n + m S 2 d dt Γ0 ūm • n Γ0 v • n - k S xm Γ0 v • n = -P0 Γ0 v • n - fext S Γ0 v • n -ρ Ω ūm (T ) • vδ T + ρ Ω ūm (0) • vδ 0 - m S 2 Γ0 u m (T ) • n Γ0 v • n δ T + m S 2 Γ0 u m0 • n Γ0 v • n δ 0 ∀v ∈ W m ,
where

x m = x 0 - 1 S t 0 Γ0
u m • n. We now apply the Fourier transform to the previous equation and take ûm as a test function, it yields

ρ2iπτ Ω |û m (τ )| 2 + ρ Ω Ĝ(τ )û m (τ ) + µ Ω ∇û m (τ ) : ∇û m (τ ) + N i=1 R i Γi ûm (τ ) • n 2 + m S 2 2iπτ Γ0 ûm (τ ) • n 2 + k2iπτ (x m (τ )) 2 = -P0 (τ ) + fext S (τ ) Γ0 ûm (τ ) • n -ρ Ω ūm (T ) • ûm (τ )e -2iπτ T + ρ Ω ūm (0) • ûm (τ ) + m S 2 Γ0 ûm (τ ) • n Γ0 ūm (T ) • n e -2iπτ T - m S 2 Γ0 ûm (τ ) • n Γ0 ūm (0) • n
where Ĝ is the Fourier transform with respect to the time of the nonlinear convection term ūm • ∇ū m .

Taking the imaginary part of the previous equality and using the inequality

Γ0 v • n ≤ C v L 2 (Ω) , for every v ∈ W (see Lemma 1.1), we have |τ | ûm (τ ) 2 L 2 (Ω) ≤ C Ĝm (τ ) V ′ ûm (τ ) V + C ūm (T ) L 2 (Ω) + ūm (0) L 2 (Ω) ûm (τ ) L 2 (Ω) +C P0 (τ ) + fext (τ ) S ûm (τ ) L 2 (Ω) + C Γ0 |ū m (T ) • n| + Γ0 |ū m (0) • n| ûm (τ ) L 2 (Ω) .
Moreover we have

G m V ′ ≤ u m 2 H 1 (Ω) ,
and then, thanks to the energy estimates satisfied by u m ,

+∞ -∞ G m V ′ dt ≤ C T 0 u m (t) 2 H 1 (Ω) dt ≤ C, from which we deduce Ĝm (τ ) V ′ ≤ C.
Futhermore, thanks to the energy estimates satisfied by u m (see Lemma 2.1), we have that u m (T ) L 2 (Ω) ≤ C and since u m (0) L 2 (Ω) ≤ C, we obtain

|τ | ûm (τ ) 2 L 2 (Ω) ≤ C ûm (τ ) H 1 (Ω) + ûm (τ ) L 2 (Ω) C ûm (τ ) H 1 (Ω) .
But, for 0 < γ < 1 4 , we see that

|τ | 2γ ≤ C(γ) 1 + |τ | 1 + |τ | 1-2γ , ∀τ ∈ R, then +∞ -∞ |τ | 2γ ûm (τ ) 2 L 2 (Ω) ≤ C(γ) +∞ -∞ ûm (τ ) H 1 (Ω) 1 + |τ | 1-2γ + C(γ) +∞ -∞ |τ |û m (τ ) H 1 (Ω) 1 + |τ | 1-2γ .
The first integral of right hand side can be estimated

+∞ -∞ ûm (τ ) 2 H 1 (Ω) 1 + |τ | 1-2γ ≤ +∞ -∞ 1 1 + |τ | 1-2γ 2 1 2 +∞ -∞ ûm (τ ) 2 H 1 (Ω) 1 2
The first integral is convergent for any 0 < γ < 1 4 . On the other hand, from Parseval egality

+∞ -∞ ûm (τ ) 2 H 1 (Ω) dτ = T 0 u m (t) 2 H 1 (Ω) dt ≤ C.
With the same kind of arguments the second integral is convergent.

Passage to the limit

The previous compactness result imply the following strong convergence (at least for a subsequence of u m still denoted u m )

u m → u strongly in L 2 (0, T ; L 2 (Ω)).
This convergence result together with (25) enable us to pass to the limit in the following weak formulation (obtained from (P m ) by multiplication by φ ∈ D([0, T )) and integration by parts in time)

                                     ρ Ω u m (t) • vφ(t) -ρ t 0 Ω u m • v φ + ρ t 0 Ω (u m • ∇u m ) • v φ + µ t 0 Ω ∇u m : ∇v φ + N i=1 R i t 0 Γi u m • n Γi v • n φ + m S 2 Γ0 u m (t) • n Γ0 v • n φ(t) - m S 2 t 0 Γ0 u m • n Γ0 v • n φ + k S 2 t 0 s 0 Γ0 u m • n Γ0 v • n φ = - t 0 P 0 Γ0 v • n φ - t 0 f ext S Γ0 v • n φ - k S x 0 Γ0 v • n t 0 φ -ρ Ω u m0 • v φ(0) - m S x m 1 Γ0 v • n φ(0), ∀v ∈ W m , where x m 1 is equal to - 1 S Γ0 u 0m •n and converges towards x 1 = - 1 S Γ0 u 0
•n as m goes to infinity. Consequently we obtain the existence of at least one weak solution locally in time in the general case or globally in time for small enough data.

Numerical illustrations

In this section we present 2D numerical results performed with Freefem++ (see [6]). Our main purpose is to illustrate how this approach makes it possible to investigate the effect of a modification of some resistances in the condensed (distal) part onto the overall flow in the upper (proximal) part of the tree. The numerical approach shall be described and justified in a forthcoming paper. We simply present in the next section the principles it lies on.

Numerical scheme

We propose a time discretization of (6) based on the method of caracteristics for the convection part. Let t 0 = 0 < t 1 < t 2 < • • • < t N = T , with t n+1 -t n = δt = T /N , denote the time steps. We denote by (u n , p n ) approximations of the velocity and pressure at time t n , built recursively according to the following scheme :

ρ 1 δt u n+1 -µ∆u n+1 + ∇p n+1 = ρ 1 δt u n • X n in Ω ∇ • u n+1 = 0 in Ω u n+1 = 0 in Γ l u n+1 = λ n+1 U i sur Γ i , i = 0...N µ∇u n+1 • n -p n+1 n = -P 0 n on Γ 0 µ∇u n+1 -p n+1 n = -P n+1 n -R i Γi u n+1 • n n sur Γ i , i = 1...N P n+1 = 1 S m x n+1 -2x n + x n-1 δt 2 + kx n+1 -F n+1 ext x n+1 -x n δt = - Q n+1 0 S , where Q n+1 0 = Γ0 u n+1
• n, and X n (x) stands for Y n (x, t n ), backward characteristic defined as

   ∂Y ∂t (x , t) = u n (Y ) Y (x , t n+1 ) = x.

Figure 2. Mesh near the first bifurcation

We may now write the variational formulation for one time step: suppose that x n , u n are known, the couple

(u n+1 , p n+1 ) ∈ {v ∈ H 1 (Ω) , v = 0 on Γ l , ∃ Λ = ( λi ) 0≤i≤N s. t. v = λi U i on Γ i , i = 0 , . . . , N } × L 2 (Ω) is computed as the solution to a(u n+1 , v) + b(v, p n+1 ) = l n+1 (v) ∀v ∈ H 1 (Ω) , v = 0 on Γ l , ∃ Λ = ( λi ) 0≤i≤N s. t. v = λi U i on Γ i , i = 0 , . . . , N b(u n+1 , q) = 0 ∀q ∈ L 2 (Ω)
(u n+1 , v) = ρ δt Ω u n+1 • v + µ Ω ∇u n+1 : ∇v + N i=1 R i Γi u n+1 • n Γi v • n + 1 S 2 m δt + kδt N i=1 Γi u n+1 • n N i=1 Γi v • n , l n+1 (v) = ρ δt Ω (u n • X n ) • v -P 0 Γ0 v • n - m δtS 2 Q n 0 + k S x n - f n+1 ext S N i=1 Γi v • n , and b(u n+1 , q) = - Ω ∇ • u n+1 q.
Once u n+1 is known, x n+1 is updated as

x n+1 = x n - δt S Q n+1 0 with Q n+1 0 = Γ0 u n+1 • n.

Numerical tests

Three dimensional computations shall be presented in a forthcoming paper. We simply present here some tests where Navier-Stokes simulations are performed on a 2D domain. From a modelling point of view, the suitability of this approach is highly questionable, as our 2D upper tree cannot be expected to reproduce all features of the real, three-dimensional one. Yet, we tried here to as respectful as possible of the mechanical properties of the human respiratory tree.

The main idea is to construct a two-dimensional tree equivalent to the three-dimensional one such that two conditions are fulfilled: the two trees have the same global resistance and the two-dimensional tree reproduces the resistance repartition of the three-dimensional tree. Consider the situation of geometrical trees, both 2D and 3D, and introduce the ratios of homothety h 2D and h 3D respectively. Knowing that the 2D and 3D resistances are proportional to 1 For the generation zero, the trachea, we keep the real dimensions for the diameter and the length (width = 2.10 -2 , lenght = 12.10 -2 ). As we know the global 3D resistance, we fix the fluid viscosity value such that both 2d and 3d global resistances are equal. Finally, the density is tuned according to the Reynolds number one wants to obtain (about 500 in the present case).

We consider the following set of data: m = 0.4 , S = 0.011 , k = 2 , µ = 4 10 -3 , ρ = 50.

All computations presented here are based on a non structured mesh, a part of which is represented in Fig 2.

Initial conditions are x = ẋ = 0, and the forcing term f ext is taken periodic (period of 5s), piecewise constant, with a 2 second active, inpiratory phase (f ext = 0.1), followed by a passive, expiratory phase (f ext = 0.0). Two situations are considered: in the first one (I), the 8 subtrees are identical, and each of them is symmetric (all resistances are the same in each generation). The second one is obtained as a perturbation of the first one: one of the resistances of subtree 5 (see the arrow on Fig. 4) is multiplied by 10 4 (the corresponding is almost obstructed). Fig. 5 represents the corresponding t → x(t) curves. As expected, the increase in resistance leads to a smaller ventilation amplitude. We plotted in Figure 6 the air fluxes through the 8 exits, at some time during the inspiration phase. The horizontal straight line corresponds to situation I: all fluxes are the same. The step function corresponds to situation II: one can check that the global flux is slightly reduced in the perturbed situation, but all fluxes are increased but through outlet 5, where it is reduced drastically. Note that the perturbation affects significantly outlet 6, and much less outlet 4. It is due to the fact that 5 and 6 are quite closed as ends of the tree, wheras 5 and 4 are not.

1 Sv

 1 Γ0 • n φ and added to the previous contribution. By integrating by parts (in space and in time) the Navier-Stokes equations and by taking into account the constrain S ẋ = -Γ0 u • n we obtain a first weak formulation of the coupled system (6):

Figure 3 .

 3 Figure 3. Pressure field in situation II where a(u n+1 , v) = ρ δt Ω u n+1 • v + µ

Figure 4 .

 4 Figure 4. Resistance perturbation

Figure 5 .Figure 6 .

 56 Figure 5. Simulation of three respiratory cycles, cases I and II

This work has been supported by the ACINIM project Lepoumonvousdisje.