L Baffico

C Grandmont

B Maury

Multiscale modelling of the respiratory tract *

Keywords: Navier-Stokes equations, local existence, coupling of models, ventilation process, Finite Element Method

We propose here a decomposition of the respiratory tree into three stages which correspond to different mechanical models. The resulting system is described by the Navier-Stokes equation coupled with an ODE (a simple spring model) representing the motion of the thoracic cage. We prove that this problem has at least one solution locally in time for any data and, in the special case where the spring stiffness is equal to zero, we obtain an existence result globally in time provided that the data are small enough. The behaviour of the global model is illustrated by three-dimensional simulations.

Introduction, modelling aspects

Breathing involves gas transport through the respiratory tract 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 of the connective tissue framework of the lung (we will refer as the parenchyma in all that follows).

At the time being, the complex fractal geometry of the airway tree makes the air flow simulation on the whole tree unreachable. Besides, the distal airways from generation 7 cannot be visualized 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] Ben-Tal | Simplified models for gas exchange in the human lungs[END_REF][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 is unsuitable to 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, by taking into account the fact that this flow depends on the distal part and is driven by the motion of the diaphragm and the parenchyma. One solution is to find physiologically relevant 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. In this spirit we propose 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 incompressible Navier-Stokes equations hold to describe the fluid flow. The flow incompressibility is valid since the Mach number in the trachea is less than 0.3, even in forced inspiration or expiration,

• 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: The present paper continues and generalizes the multi-model strategies proposed, for air flow modelling, in [START_REF] Grandmont | A multiscale / multimodel approach of the respiration tree[END_REF] or [START_REF] Grandmont | Multiscale modelling of the respiratory tract: a theoretical framework[END_REF]. The same type of multiscale strategies has been developped and investigated for different applications such as blood flow simulations and several groups have successfully coupled three-dimensional models to either resistances or more sophisticated zero-dimensional models (lumped models) or one-dimensional models (see [START_REF] Formaggia | Multiscale modelling of the circulatory system: A preliminary analysis[END_REF][START_REF] Olufsen | Structured tree outflow condition for blood flow in larger systemic arteries[END_REF][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] Quarteroni | Coupling between lumped and distributed models for blood flow problems[END_REF][START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modelling of blood and pressure in arteries[END_REF]).

The paper is organized as follows: in the first part we present the coupled system and its variational formulation from which we derive, at least formally, an energy balance. To obtain energy estimates from this energy balance, the main difficulty is to estimate the nonlinear terms that is to say the flux, at the artificial boundaries Γ i , of the system kinetic energy. For our model the energy balance leads to energy estimates only in dimension two, locally in time and for small data. Nevertheless, in the second part we prove, thanks to a Galerkin approach using a well chosen special basis, that there exists a unique smooth solution, locally in time. Moreover, in the special case where the spring stiffness is zero, we obtain that the solution exists globally in time provided the data are small enough. The proofs are based on the same kind of arguments used in [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF] and rely on uniform estimates in order to pass to the limit in the discrete system. In the third part, we detail the time discretization strategy used and finally we give 3D numerical evidence that this model can reproduce some aspects of normal or pathological breathing.

Problem setting

In the upper part, denoted by Ω (see Fig. 1), we assume that the Navier-Stokes equations hold:

                     ρ ∂u ∂t + ρ(u • ∇)u -µ△u + ∇p = 0, in Ω , ∇ • u = 0, in Ω , u = 0, on Γ ℓ , µ∇u • n -pn = -P 0 n on Γ 0 , µ∇u • n -pn = -Π i n on Γ i i = 1, . . . , N, (1)
where u and p are respectivelly the fluid velocity and the fluid pressure. On the lateral boundary Γ ℓ we impose no-slip boundary conditions on the velocity, whereas on the artificial boundary Γ i , 0 ≤ i ≤ N we consider 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 subtrees should be a dyadic tree in which we assume that the flow is laminar. Thus, by analogy with an electric network, we can consider that the flow is characterized by a unique equivalent resistance (referred to as lumped model) that depends on each resistance of the local branches (see for instance [START_REF] Olufsen | Structured tree outflow condition for blood flow in larger systemic arteries[END_REF], [START_REF] Maury | Outlet Dissipative conditions for air flow in the bronchial tree[END_REF] or [START_REF] Grandmont | A viscoelastic model with non-local damping application to the human lungs[END_REF]). Thus, each of the subtrees is replaced by a cylindrical domain, where the flow satisfies Poiseuille's 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 stands for the alveolar pressure (which is supposed to be uniform). Note that R i depends on the geometric properties (length and diameter) of all the branches of the i-th subtree. Thanks to the relation [START_REF] Bahouri | Ingalités de Hardy précisées[END_REF] the boundary conditions at the outlets Γ i writes

µ∇u • n -pn = -P i n -R i Γ i u • n n on Γ i i = 1, . . . , N. (3)
These are non standard nonlocal boundary conditions that link the fluid stress tensor and its flux. Note that they induce dissipation in the system as we will detail it in the next section. We will call them natural dissipative boundary conditions. Note that similar conditions are also used for blood flow modelling [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]. Next, we assume that all the alveola pressures are equal: P i = P a . Finally, we suppose that these alveoli are embedded in a box filled with an incompressible medium that represents the connective tissus framework of the lung. This incompressibility assumption is valid since the lung parenchyma is made mostly of elastin and collagen fibers and of blood vessels. One part of the box is connected to a spring that governs the diaphragm and parenchyma motion. The equation satisfied by the position x of the diaphragm writes:

mẍ = -kx + f ext + f P , ( 4
)
where m is the total mass of the lung, k is the stiffness of the spring (that characterizes 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, we have

f P = P a S. (5)
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. (6)
Note that, thanks to the fluid incompressibility, we have

N i=1 Γ i u • n = - Γ 0 u • n,
and consequenlty

S ẋ = - Γ 0 u • n. (7)
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 ) × Γ ℓ , µ∇u • n -pn = -P 0 n on (0, T ) × Γ 0 , µ∇u • n -pn = -P a n -R i Γ i u • n n , on (0, T ) × Γ i , i = 1 , . . . , N , mẍ + kx = f ext + SP a , S ẋ = N i=1 Γ i u • n = - Γ 0 u • n. (8)
This system of equations have to be completed by suitable initial conditions

(u, x, ẋ)| t=0 = (u 0 , x 0 , x 1 ), with ∇ • u 0 = 0 , u 0 = 0 on Γ ℓ , Sx 1 = - Γ 0 u 0 • n. (9)
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] 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.

Remark 1.1 From a modelling standpoint, it makes sense to prescribe on Γ 0 also a dissipative boundary condition, to account for the resistance of the upper part of the respiratory tract (nose, pharynx, and larynx). Numerical tests (see Section 4) will be based on this assumption. As it does not change the analysis of the system, we consider here the case of a zero resistance on the inlet, to alleviate notations.

Remark 1.2 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, in particular the fact that the motion of the lung parenchyma is a viscolelastic media and that its motion is limited by the chest wall. We refer to [START_REF] Martin | Modelling of respiratory system mechanics involving gas exchange in the human lungs[END_REF] for a more sophisticated spring model.

Remark 1.3 By setting p = p -P a , by using the expression of the alveolar pressure

P a = m S ẍ + k S x - f ext S
and the volume preserving equation S ẋ = -Γ 0 u • n, the coupled system 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 ) × Γ ℓ , µ∇u • n -pn = -P 0 n - f ext S n - m S 2 d dt Γ 0 u • n n - k S 2 t 0 Γ 0 u • n -Sx 0 n , on (0, T ) × Γ 0 , µ∇u • n -pn = -R i Γ i u • n n , on (0, T ) × Γ i .
(10) Consequently the coupled system (8) reduces to the Navier-Stokes equations, with fluid pressure replaced by the difference between fluid pressure and alveolar pressure, and with generalized natural dissipative boundary conditions. Note that the assumption of incompressiblity is essential here.

Variational formulation

Assuming that all the unknowns are regular enough, we multiply the Navier-Stokes equations by a test-field v which vanishes on Γ ℓ , and the spring equation by

-(1/S) Γ 0 v • n. Using x = x 0 - 1 S t 0 Γ 0 u • n, we obtain                              ρ Ω ∂ t u • v + ρ Ω (u • ∇)u • v + µ Ω ∇u : ∇v + N i=1 R i Γ i u • n Γ i v • n + m S 2 Γ 0 ∂ t u • n Γ 0 v • n + k S 2 t 0 Γ 0 u • n Γ 0 v • n - Ω p∇ • v + P a N i=1 Γ i v • n - Γ 0 v • n = -P 0 Γ 0 v • n - f ext S Γ 0 v • n + k S 2 Sx 0 Γ 0 v • n , ∀v. (11)
Next, considering test functions v that are divergence free, we obtain a second variational formulation of the coupled problem

                     ρ Ω ∂ t u • v + ρ Ω (u • ∇)u • v + µ Ω ∇u : ∇v + N i=1 R i Γ i u • n Γ i v • n + m S 2 Γ 0 ∂ t u • n Γ 0 v • n + k S 2 t 0 Γ 0 u • n Γ 0 v • n = -P 0 Γ 0 v • n - f ext S Γ 0 v • n + k S x 0 Γ 0 v • n , ∀v. (12)
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 S ẋ = -

Γ 0 u • n.
Remark 1.4 Notice that alveolar pressure can be interpreted as the Lagrange multiplier associated to the parenchyma incompressibility constraint S ẋ = N i=1 Γ i u • n.

Remark 1.5 Note that in the previous variational formulation we have used the relation

Γ 0 u • n = -N i=1 Γ i u • n.
This simplifies the weak formulation and will simplify the proof. However in the case where the lateral boundary Γ l is moving this equation does not hold anymore. If we do not use this in-flux /out-flux balance the variational formulation writes:

                         ρ Ω ∂ t u • v + ρ Ω (u • ∇)u • v + µ Ω ∇u : ∇v + N i=1 R i Γ i u • n Γ i v • n + m S 2 N i=1 Γ i ∂ t u • n N i=1 Γ i v • n + k S 2 t 0 N i=1 Γ i u • n N i=1 Γ i v • n = -P 0 Γ 0 v • n - f ext S N i=1 Γ i v • n + k S x 0 N i=1 Γ i v • n , ∀v. (13)
Note that the existence of solution as well as the numerical simulations can be performed also in this case with only a few changes which we will underline in what follows (see remarks 2.5 and 3.2).

Energy balance

We start to derive, at least formally, an energy balance for the coupled system. We take v = u as a test function in [START_REF] De Rochefort | Phase contrast velocimetry with hyperpolarized helium-3 for in vitro and in vivo characterization of airflow[END_REF]. By writing the convective terms, thanks to an integration by parts, as a flux of kinetic energy at the inlet and at the outlets, and remembering that S ẋ = -Γ 0 u • n, we obtain

d dt ρ 2 Ω |u| 2 + m 2 | ẋ| 2 + k 2 |x| 2
Total energy

+ µ Ω |∇u| 2 Dissipation within Ω + N i=1 R i Γ i u • n 2 Dissipation in the subtrees = - ρ 2 N i=0 Γ i |u| 2 (u • n)
In/out-come of kinetic energy

+ P 0 S ẋ Power of inlet pressure + f ext ẋ Power of ext. forces . (14)
Identity ( 14) represents the energy balance of the coupled system. The two first terms represent the variation of the total kinetic energy, the third one being the variation of the mechanical energy of the spring. Moreover, 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 u • n 2 .
This energy balance will provide energy estimates in the case where the flux of kinetic energy at the interface can be correctly estimated by the global kinetic energy and the dissipated energy. Unfortunatly we can control the flux of kinetic energy only locally in time and for small data and only in the case where d = 2. Note that in [START_REF] Grandmont | A multiscale / multimodel approach of the respiration tree[END_REF] this difficulty has been overcome by considering other type of boundary conditions at the inlet and outlets. More precisely the boundary conditions do not involve the fluid stress tensor but the total fluid stress tensor, i.e. where the fluid pressure p has been replaced by the total pressure p + ρ |u| 2 2 . Consequently, the flux of kinetic energy on the boundary do not appear in the energy balance. An alternative approach is proposed in [START_REF] Grandmont | Multiscale modelling of the respiratory tract: a theoretical framework[END_REF], where a more suited estimate of the nonlinear convective term is obtained thanks to an additionnal assumption that states that the velocity profiles at the inlet and outlets are given. Note that this estimate relies on the fact that the trace of the fluid velocity on each Γ i is supposed to be described by a single parameter.

2 Well-posedness issues

Mathematical setting, notations

In this section we put the variational formulation ( 12) onto an abstract form, on which we shall base the well-posedness proof.

Let Ω be a bounded domain in R d , d = 2 or 3, Γ 0 the inlet boundary, Γ i the i-th. outlet, and Γ ℓ the lateral boundary (see Fig. 1). We introduce the following functional spaces:

V = {v ∈ H 1 (Ω) d , ∇ • v = 0, v = 0 on Γ ℓ }. H = V L 2 , We denote by (•, •) 0 the scalar product on H × H defined by (v, w) 0 = ρ Ω v • w + m S 2 Γ 0 v • n Γ 0 w • n , (15)
and by • 0 the associated norm (see Remark 2.4 below for a proper definition of the flux

Γ 0 v • n for v ∈ H). Next we set a 1 (v, w) = µ Ω ∇v : ∇w + N i=1 R i Γ i v • n Γ i w • n , and b(u, v, w) = ρ Ω (u • ∇)v • w,
that are the bilinear and trilinear forms on V associated to our system. Moreover we introduce

a 2 (v, w) = k S 2 Γ 0 v • n Γ 0 w • n .
Finally we introduce an operator A whose eigenvectors will constitute the Galerkin basis used to build our sequence of approximate solutions.

Definition 2.1 (Operator A).

The operator A is defined on H as follows:

D(A) = {v ∈ V , |a 1 (v, w)| ≤ C v 0 , ∀w ∈ V }, (Av, w) 0 = a 1 (v, w), ∀(v, w) ∈ D(A) × V. ( 16
)
Remark 2.2 Operator A is similar to Stokes operator, with special kinds of boundary conditions. More precisely, one can check that Au = f if and only if there exists a pressure field p over Ω such that (u, p) is a weak solution to the Stokes problem

                           -µ△u + ∇p = f in Ω ∇ • u = 0 in Ω u = 0 on Γ ℓ , µ∇u • n -pn = m S 2 Γ 0 f • n n on Γ 0 , µ∇u • n -pn = -R i Γ i u • n n on Γ i , i = 1, . . . , N. (17)
We may now write the variational formulation in an abstract setting, on which we base the notion of solution to our problem.

Definition 2.3

We shall say that u is a solution of ( 8)-( 9

) on [0, T ] if      u ∈ L 2 (0, T ; D(A)) ∩ L ∞ (0, T ; V ) ∩ H 1 (0, T ; H), d dt (u, v) 0 + a 1 (u, v) + b(u, u, v) + a 2 t 0 u, v = ℓ(v) ∀ v ∈ V, (18)
with ℓ(v) = -P 0 Γ 0 v • n - f ext S Γ 0 v • n + k S 2 Sx 0 Γ 0 v • n ,
v • n, g H -1/2 (∂Ω),H 1/2 (∂Ω) = Ω v • ∇g,
where g is any function in H 1 (Ω) such that g = 1 on Γ 0 , and which vanishes on Γ i for i = 1, . . . , N . Such a function exists since in-and out-let boundaries Γ i are not in contact. Note that if v ∈ V then the provided definition corresponds to the flux of v through Γ 0 , because, on the lateral boundary

Γ ℓ , v • n vanishes, thus v • n, g H -1/2 (∂Ω),H 1/2 (∂Ω) = Γ 0 v • n.
We can define in the same way the velocity flux

Γ i v • n on Γ i for v ∈ H.
Remark 2.5 If we had chosen to work with the variational formulation (13) then the scalar product on H and the bilinear form associated to the spring energy have to be modified respectivelly as follows:

(v, w) 0 = ρ Ω v • w + m S 2 N i-1 Γ i v • n N i-1 Γ i w • n , a 2 (v, w) = k S 2 N i-1 Γ i v • n N i-1 Γ i w • n .
The study of the operator A is then very similar and the proof of existence can be easily adapted to this case.

Preliminaries

We gather here some properties, mainly pertaining to the operator A which we introduced in the previous section. We begin with an estimate for divergence free fields that will be essential to control boundary integrals:

Lemma 2.6 It holds Γ 0 v • n ≤ C v L 2 (Ω) , ∀v ∈ H,
where the boundary integral is defined according to Remark 2.4. As a consequence, • 0 is equivalent to the L 2 norm on H.

This a direct consequence of the definition of (H -1/2 , H 1/2 )-duality, based on the Green formula:

Γ 0 v • n := v • n, g H -1/2 (∂Ω),H 1/2 (∂Ω) = Ω g∇ • v + Ω v • ∇g,
where g ∈ H 1 (Ω) is 1 on Γ 0 , and vanishes on Γ i for i = 1, . . . , N .

Proposition 2.7

The operator A has the following properties: i) A ∈ L(D(A) , H) is invertible and its inverse is compact on H.

ii) A is selfadjoint.

As a consequence, A admits a family of eigenfunctions (φ j )

Aφ j = λ j φ j with 0 < λ 1 ≤ λ 2 ≤ . . . λ j ----→ j→+∞ +∞,
which is complete and orthogonal in both H and V .

Proof: The proof of this proposition relies on classical arguments (see for instance [START_REF] Brezis | Analyse Fonctionnelle, Théorie et applications[END_REF]).

As in [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF] (or in [START_REF] Quarteroni | Analysis of a Geometrical Multiscale Model Based on the Coupling of ODEs and PDEs for Blood Flow Simulations[END_REF]), the control of the nonlinear term in a priori estimates is based on regularity properties of operator A: Proposition 2.8 Let Ω be a domain in R d . We suppose that Ω is piecewise smooth (in-/out-let and lateral component of the boundary are smooth, see Fig. 1), and that the Γ i 's meet with lateral boundary Γ ℓ at angle π/2. Let A be defined by Def. 2.1. There exists ε > 0 and a constant C such that, for any

f ∈ H, Au = f , it holds u H 3/2+ε (Ω) ≤ C f 0 .
Proof: First of all, as f is in H, it is divergence free, so that

Γ i f • n ≤ C f 0 . (19)
By Remark 2.2, u solves Stokes problem with right-hand side f , homogeneous Dirichlet boundary conditions on the lateral boundaries, dissipative conditions on the Γ ′ i s, and free boundary conditions on the inlet Γ 0 with an oustide pressure which is controlled by the normal flux of f . By taking u as a test function, and using [START_REF] Grandmont | A viscoelastic model with non-local damping application to the human lungs[END_REF], we have

u H 1 (Ω) ≤ C f 0 . Now consider a auxilliary field p ∈ H 1 (Ω) defined as follows: It is harmonic over Ω, p = m S 2 Γ 0 f • n on Γ 0 , p = -R i Γ i u • n on Γ i , i = 1, . . . , N , and ∂ p ∂n = 0 on Γ ℓ .
As normal fluxes of u and f are both controlled by f 0 , so is the H 1 norm of p. Problem ( 17) can therefore by written

                   -µ△u + ∇(p -p) = ρf -∇p in Ω ∇ • u = 0 in Ω u = 0 on Γ ℓ , µ∇u • n -(p -p)n = 0 on Γ 0 , µ∇u • n -(p -p)n = 0 on Γ i , i = 1, . . . , N. (20)
which is a standard Stokes problem with free out/in-let boundary conditions at the ends, and a new right hand-side with a L 2 norm controlled by f 0 . As H 2 regularity is not an issue away from the Γ ′ i s, we end up with a question of regularity for Stokes problem in a single pipe, with free outlet B.C.'s and L 2 right-hand side.

This regularity is a consequence of recent regularity results for the Stokes problem in nonsmooth domains, for mixed boundary conditions. In the two-dimensional setting, the H 3/2+ε regularity of u, and continuous dependance upon the data f , is given in [START_REF] Orlt | Regularity of viscous Navier-Stokes flows in nonsmooth domains[END_REF]. For the three dimensionnal problem, the key point lies in some weighted estimates for the Stokes problem in a dihedron, which we describe below. We shall disregard here any considerations for conditions at infinity, as they do not play any role in our case.

Up to a C ∞ diffeomorphism, the domain is locally similar quarter-space like domain. Therefore we consider a dihedral domain Ω delimited by flat hypersurfaces Γ 1 and Γ 2 , which meet at angle θ, and (w, q) the solution to

-∆w + ∇q = g ∇ • w = 0;
with homogeneous boundary conditions

∇w • n -qn = 0 on Γ 1 , w = 0 on Γ 2 .
Lemma 2.16 in [START_REF] Maz | Point estimates for Green's matrix to boundary value problems for second order elliptic systems in a polyhedral cone ZAMM Z[END_REF] asserts that u is in some weighted Sobolev space W 2 δ (Ω), with

w W 2 δ ≤ C g 0
as soon as none of the real parts of the solutions to cos(λθ) λ 2 sin 2 θ -cos 2 (λθ) = 0 fits into [0, 1 -δ]. Considering θ = π/2, the solution to the previous equation with the smallest positive real part is λ c ≈ 0.59. As a consequence, δ = 1/2 -ε qualifies as soon as ε is sufficiently small (ε < 0.09).

The norm on W 2 δ is defined by

w 2 W 2 δ = |α|≤2 1≤i≤3 Ω r 2δ |∂ α w i | 2 ,
where r is the distance to the edge of Ω. The expected estimate is a consequence of a Hardy type inequality (see Proposition A.1 in the appendix):

|u| Ḣ2-δ ≤ C u 2 W 2 δ ,
which yields the expected estimate for δ = 1/2 -ε, in the case of a quarter-space domain.

Corollary 2.9 Under the assumptions of Proposition 2.8, the following inequality holds

v L ∞ (Ω) ≤ C Av 0 ,
and there exists θ ∈ (0, 1) such that

v L ∞ (Ω) ≤ C ∇v θ L 2 (Ω) Av 1-θ 0 .
Proof: The first estimate is a direct consequence of H 3/2+ε (Ω) ֒→ L ∞ (Ω). As for the second estimate, we use the interpolation estimate (for 0 < ε ′ < ε),

v H 3/2+ε ′ (Ω) ≤ C ∇v θ L 2 (Ω) Av 1-θ H 3/2+ε (Ω) , valid for θ such that 3/2 + ε ′ = θ + (1 -θ)(3/2 + ε).

Existence and uniqueness of a solution

We establish here two types of existence results: one locally in time for large data, and one globally in time for small enough data in the special case where the stiffness of the spring is equal to zero. Let us first make explicit our assumptions:

Ω is piecewise regular in the sense of Proposition 2.8, which is assumed to hold true (21)

T > 0, P 0 ∈ L 2 (0, T ), f ext ∈ L 2 (0, T ) , u 0 ∈ V , R i ≥ 0 , i = 1, . . . , N. (22)
Theorem 2.10 Under assumptions ( 21)- [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF], there exists a time interval [0, t ⋆ ] on which ( 8)-( 9) has a unique solution in the sense of Def. 2.3. Moreover, in the special case where k = 0 and for small data (initial conditions and applied forces small enough in norm L ∞ in time), then the solution can be defined globally in time, and there exists C > 0 such that for all t > 0

∇u 2 L 2 (Ω) + N i=1 R i Γ i u • n 2 ≤ C.
Proof: We follow the same steps as in [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF]. First we build a sequence of approximated solutions thanks to a Galerkin method. Then, we derive uniform bounds that enable us to pass to the limit in the equation leading to the existence of at least one solution of [START_REF] Grandmont | A multiscale / multimodel approach of the respiration tree[END_REF]. Finally we prove that this solution is unique.

Remark 2.11 In [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF] or [START_REF] Quarteroni | Analysis of a Geometrical Multiscale Model Based on the Coupling of ODEs and PDEs for Blood Flow Simulations[END_REF], the existence of a unique smooth solution of the Navier-Stokes system is proved for other types of boundary conditions, namely mean pressure or mean flux boundary conditions ( [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF]) or boundary conditions that write µ∇u [START_REF] Quarteroni | Analysis of a Geometrical Multiscale Model Based on the Coupling of ODEs and PDEs for Blood Flow Simulations[END_REF]). Note that in the former, the existence is obtained under a condition on the size of the R i 's coefficients which we do not need here, thanks to the choice of the operator A (that takes into account the dissipative part coming from the subtrees resistances). Moreover in [START_REF] Quarteroni | Analysis of a Geometrical Multiscale Model Based on the Coupling of ODEs and PDEs for Blood Flow Simulations[END_REF] the existence of a solution of the Navier-Stokes system coupled with lumped models is proved thanks to a fixed point argument. Note that we can not use here the same kind of argument due to the nature of the coupling. That is the reason why we have considered a non standard scalar product on H taking into account the spring mass in order to treat the system globally.

• n -pn + R i (u • n)n = g ([
Galerkin approximation. For all n ∈ N, we define the following discrete problem: Find u n (t) ∈ span(φ 1 , . . . , φ n ) such that

       ∂u n ∂t , v 0 + a 1 (u n , v) + b(u n , u n , v) + a 2 t 0 u n , v = ℓ(v), ∀v ∈ span(φ 1 , . . . , φ n ) , (u n (0) -u 0 , v) 0 = 0, ∀v ∈ span(φ 1 , . . . , φ n ). ( 23
)
The first step consists in proving that this differential system has a unique solution on a time interval [0, t ⋆ ], where t ⋆ does not depend on n. Let us denote by t n the life-time of the maximal solution u n over [0, T ].

A priori estimates.

In what follows, C stands for a generic constant whose value may change, but it does not depend on n. We write [START_REF] Lin | Characteristics of the turbulent laryngeal jet and its effect on airflow in the human intra-thoracic airways[END_REF] for the test function v = Au n , which is admissible thanks to the choice of the Galerkin basis. Indeed Au n ∈ span(φ 1 , . . . , φ n ). Taking into account the two following equalities

∂u n ∂t , Au n 0 = µ 1 2 d dt ∇u n 2 L 2 (Ω) + 1 2 N i=1 R i d dt Γ i u n • n 2 ,
and

a 1 (u n , Au n ) = Au n 2 0 , we obtain 1 2 µ d dt ∇u n 2 L 2 (Ω) + 1 2 N i=1 R i d dt Γ i u n • n 2 + Au n 2 0 = -b(u n , u n , Au n ) -a 2 t 0 u n , Au n + ℓ(Au n ).
Now thanks to Corollary 2.9 and Lemma 2.6, we can control the nonlinear term:

|b(u n , u n , Au n )| ≤ C u n L ∞ (Ω) ∇u n L 2 (Ω) Au n L 2 (Ω) ≤ C ∇u n 1+θ L 2 (Ω) Au n 2-θ 0
with θ ∈ (0, 1)

≤ C ∇u n 2(1+θ)/θ L 2 (Ω) + 1 4 Au n 2 0 ,
by Young's inequality. Moreover, using Lemma 2.6, the equivalence of the L 2 norm and • 0 on H and Poincaré inequality, we obtain

a 2 t 0 u n , Au n = k t 0 Γ 0 u n • n Γ 0 Au n • n ≤ Ck t 0 ∇u n 2 L 2 (Ω) + 1 4 Au n 2 0 , and
|ℓ(Au n )| ≤ C |P 0 | 2 + |f ext | 2 S 2 + k 2 |x 0 | 2 S 2 + 1 4 Au n 2 0 . Consequently 1 2 µ d dt ∇u n 2 L 2 (Ω) + 1 2 N i=1 R i d dt Γ i u n • n 2 + 1 4 Au n 2 0 ≤ C ∇u n 2(1+θ)/θ L 2 (Ω) + C |P 0 | 2 + |f ext | 2 S 2 + k 2 |x 0 | 2 S 2 + kC t 0 ∇u n 2 L 2 (Ω) . (24)
Next we set

Φ n (t) = µ ∇u n 2 L 2 (Ω) + N i=1 R i Γ i u n • n 2 + k t 0 ∇u n 2 L 2 (Ω) , Ψ n (t) = 1 2 Au n 2 0 , and
f (t) = C |P 0 | 2 + |f ext | 2 S 2 + k 2 |x 0 | 2 S 2 .
Then, by adding

k 2 ∇u n 2
L 2 (Ω) on both sides, [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] writes

d dt Φ n (t) + Ψ n (t) ≤ C Φ (1+θ)/θ n (t) + Φ n (t) + f (t), (25)
Besides, the choice of the initial condition

(u n (0) -u 0 , v) 0 = 0, ∀v ∈ span(φ 1 , • • • , φ n ), implies, choosing v = Au n (0), (u n (0) , Au n (0)) 0 = (u 0 , Au n (0)) 0 , which yields a 1 (u n (0) , u n (0)) = a 1 (u n (0) , u 0 ).
Thus, Cauchy-Schwarz inequality implies

Φ n (0) = µ ∇u n (0) 2 L 2 (Ω) + N i=1 R i Γ i u n (0) • n 2 ≤ µ ∇u 0 2 L 2 (Ω) + N i=1 R i Γ i u 0 • n 2 . (26)
Let us denote by t ⋆ ≤ T a positive time such that the solution F to

F ′ = C F (1+θ)/θ + F + f , F (0) = µ ∇u 0 2 L 2 (Ω) + N i=1 R i Γ i u 0 • n 2 is well-defined over [0, t ⋆ ]
. By a lemma of comparison for differential inequalities (see e.g. [START_REF] Heywood | The Navier-Stokes equations : on the Existence, Regularity and decay of Solution[END_REF]), we deduce from (25) that u n is well-defined as the unique solution to (23) over [0, t ⋆ ] (in particular, t ⋆ < t n ), and that

Φ n (t) = µ ∇u n 2 L 2 (Ω) + N i=1 R i Γ i u n • n 2 + k t 0 ∇u n 2 L 2 (Ω) ≤ F (t), ∀t ∈ [0, t ⋆ ]. (27)
Besides, by integrating (25) over (0, t ⋆ ), we get

t ⋆ 0 Au n 2 0 ≤ C < ∞. (28)
We thus deduce that u n is uniformly bounded in L ∞ (0, t * ; V ) ∩ L 2 (0, t * ; D(A)). Thanks to these estimates we can easily deduce a uniform bound for

t 0 ∂u n ∂τ 2 L 2 (Ω)
dτ , by choosing v = ∂u n ∂t in [START_REF] Grandmont | A multiscale / multimodel approach of the respiration tree[END_REF].

Indeed ∂u n ∂t 2 0 = -ρ Ω (u n • ∇u n ) ∂u n ∂t -(Au n , ∂u n ∂t ) -k t 0 Γ 0 u n • n Γ 0 ∂u n ∂t • n + ℓ ∂u n ∂t ,
and thus, using Lemma 2.6, the equivalence of the L 2 norm and • 0 on H, estimates ( 27) and ( 28), and Corollary 2.9,

∂u n ∂t 2 L 2 (Ω) ≤ C T u n L ∞ (Ω) ∇u n L 2 (Ω) + Au n L 2 (Ω) + |P 0 | + |f ext | S + k|x 0 | S ∂u n ∂t L 2 (Ω) , ∂u n ∂t L 2 (Ω) ≤ C T ∇u n 1+θ L 2 (Ω) Au n 1-θ L 2 (Ω) + Au n L 2 (Ω) + |P 0 | + |f ext | S + k|x 0 | S .
The previous bounds, together with [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF], thus imply that ∂u n ∂t is uniformly bounded in L 2 (0, t * ; L 2 (Ω)). Passage to the limit Consequently there exists a subsequence of (u n ), also denoted (u n ), such that

         u n ⇀ u weakly * in L ∞ (0, t * ; V ) and ∂u n ∂t ⇀ ∂u ∂t weakly in L 2 (0, t * ; L 2 (Ω)). (29)
Futhermore thanks to Aubin's compactess lemma (see [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], p. 57) we have, up to a subsequence (still denoted (u n )),

u n -----→ n→+∞ u in L 2 (0, t * ; L 2 (Ω)).
It enables us to pass to the limit in the discrete problem and to obtain the existence, locally in time, of at least one solution u of (18).

Global existence. Here we consider the special case where the spring stiffness k is equal to zero. Moreover we assume that the data P 0 , f ext are in L ∞ (0, T ). Once again, by taking Au n as a test function in the discrete problem yields

µ 2 d dt ∇u n 2 L 2 (Ω) + 1 2 N i=1 R i d dt Γ i u n • n 2 + Au n 2 0 ≤ ρ u n L ∞ (Ω) ∇u n L 2 (Ω) Au n L 2 (Ω) + C |P 0 | + |f ext | S Au n L 2 (Ω) .
By Corollary 2.9 we have

u n L ∞ (Ω) ≤ C Au n 0 , (30)
thus, using Young inequality and the equivalence of the L 2 norm and

• 0 on H, µ 2
d dt ∇u n 2 L 2 (Ω) + 1 2 N i=1 R i d dt Γ i u n • n 2 + Au n 2 0 ≤ C Au n 2 0 ∇u n L 2 (Ω) + C |P 0 | 2 + |f ext | 2 S 2 + 1 2 Au n 2 0 ,
and consequently

µ d dt ∇u n 2 L 2 (Ω) + N i=1 R i d dt Γ i u n • n 2 + 1 -C ∇u n L 2 (Ω) Au n 2 0 ≤ C |P 0 | 2 + |f ext | 2 S 2 .
Thanks to the definition of A we have ∇u n L 2 (Ω) ≤ C Au n 0 and since

Γ i v • n ≤ C ∇v L 2 (Ω)
for any v ∈ H 1 (Ω), we then obtain an inequality of the type

d dt Φ n + B 1 (1 -B 2 √ Φ n )Φ n ≤ C |P 0 | 2 + |f ext | 2 S 2 ,
which is valid provided 1 -B 2 √ Φ n > 0 and where

Φ n = µ ∇u n 2 L 2 (Ω) + N i=1 R i Γ i u n • n 2 ,
and B 1 , B 2 are positive constants. Consequently, assuming

C |P 0 | 2 + |f ext | 2 S 2 ≤ B 1 8B 2 2 and µ ∇u 0 2 L 2 (Ω) + N i=1 R i Γ i u 0 • n 2 ≤ 1 2B 2
the latter inequality implying that Φ n (0) ≤ 1 2B 2 (see ( 26)), we can check that Φ n is defined over [0, +∞), with

Φ n (t) ≤ 1 2B 2 ∀t.
Under these assumptions, the solution is defined globally in time.

Remark 2.12 Here we assumed that k = 0 because we are not able to control the displacement of the spring. Nevertheless, if we assume a priori that the spring displacement is bounded by a given (small enough) constant then one could obtain the same global existence result in the case of a non zero spring stiffness. On the other hand, global existence of weak solutions can be established without any assumptions on the data (including the stiffness k) if neglected the non linear convective terms of the Navier-Stokes equations. The difficulty here lies in the fact that the nonlinear system we consider (Navier-Stokes + spring-mass) is capable of storing energy in a form which is not counterbalanced by any dissipation phenomenon (potential energy of the spring).

Uniqueness. Let us consider two solutions u 1 , u 2 of (18):

∂u i ∂t , v 0 + a 1 (u i , v) + a 2 t 0 u i , v + b(u i , u i , v) = ℓ(v) ∀v ∈ V , i = 1 , 2,
associated to the same initial data u 0 . Setting w = u 1 -u 2 , substracting the two previous equations and taking v = w as a test function, we obtain 1 2

d dt w 2 0 + k 2S 2 d dt Γ 0 w • n 2 + µ ∇w 2 L 2 (Ω) + N i=1 R i Γ i w • n 2 = b(u 1 , u 1 , w) -b(u 2 , u 2 , w).
We have

b(u 1 , u 1 , w) -b(u 2 , u 2 , w) = b(u 2 , w , w) + b(w , u 2 , w) + b(w , w , w) , and
|b(u 2 , w , w)| ≤ u 2 L ∞ (Ω) ∇w L 2 (Ω) w L 2 (Ω) .
Besides, thanks to Hölder inequality,

|b(w , u 2 , w)| ≤ w L 6 (Ω) ∇u 2 L 2 (Ω) w L 3 (Ω) and |b(w , w , w)| ≤ w L 6 (Ω) ∇w L 2 (Ω) w L 3 (Ω) . Now Sobolev injections for d = 2 , 3 lead to v L 6 (Ω) ≤ C ∇v L 2 (Ω) and v L 3 (Ω) ≤ C v 1 2 L 2 (Ω) ∇v 1 2
L 2 (Ω) . These estimates, along with Young and Poincaré's inequalities, imply 1 2

d dt w 2 0 + k 2S 2 d dt Γ 0 w • n 2 + µ ∇w 2 L 2 (Ω) ≤ C u 2 L ∞ (Ω) w L 2 (Ω) 2 + µ 4 ∇w 2 L 2 (Ω) + C ∇u 2 4 L 2 (Ω) w 2 L 2 (Ω) + µ 4 ∇w 2 L 2 (Ω) + C w 1 2 L 2 (Ω) ∇w 5 2 L 2 (Ω) . Finally ρ 2 d dt w 2 0 + k S 2 d dt Γ 0 w • n 2 + µ -C w 1 2 L 2 (Ω) ∇w 1 2 L 2 (Ω) ∇w 2 L 2 (Ω) (31)
≤ C u 2 2 L ∞ (Ω) + ∇u 2 4 L 2 (Ω) w 2 L 2 (Ω)
and thus on a time interval on which w

1 2 L 2 (Ω) ∇w 1 2 L 2 (Ω) ≤ µ C
, we have, since the L 2 norm and • 0 are equivalent on H,

d dt w 2 0 + k S 2 d dt Γ 0 w • n 2 ≤ C u 2 2 L ∞ (Ω) + ∇u 2 4 L 2 (Ω) w 2 0 . ( 32
)
Taking into account that u 2 is bounded in L ∞ (0, T ; V ) ∩ L 2 (0, T ; D(A)) and the fact that D(A) ⊂ L ∞ (Ω) (see Corollary 2.9), Gronwall lemma implies that w ≡ 0 since w(0) = 0. Note that the same type of estimates could be used to prove that the solution depends continuously on the data. In particular, if u 2 2

L ∞ (Ω) + ∇u 2 4 L 2 (Ω)
is small enough, using [START_REF] Orlt | Regularity of viscous Navier-Stokes flows in nonsmooth domains[END_REF] then any small perturbation w of u 2 decreases exponentially.

Remark 2.13

The operator A which we use to establish the existence result is not the classical Stokes operator: it takes into account the dissipative terms associated to the resistances R i . If we had chosen not to take them into account, the existence result, even locally in time, would have required some control of these resistances (see [START_REF] Quarteroni | Analysis of a Geometrical Multiscale Model Based on the Coupling of ODEs and PDEs for Blood Flow Simulations[END_REF] where a similar proof is proposed).

Numerical method

We present the numerical strategy we adopted, which leads to an implicit implementation of the dissipative boundary conditions (see [START_REF] Maury | Outlet Dissipative conditions for air flow in the bronchial tree[END_REF]).

Note that most of the numerical studies in realistic airways geometries we are aware of (see [START_REF] Croce | In vitro experiments and numerical simulations of airflow in realistic nasal airway geometry[END_REF], [START_REF] De Backer | Computational fluid dynamics can detect changes in airway resistance in asthmatics after acute bronchodilation[END_REF], [START_REF] De Rochefort | Phase contrast velocimetry with hyperpolarized helium-3 for in vitro and in vivo characterization of airflow[END_REF], [START_REF] Fetita | An image based computational model of oscillatory flow in the proximal part of tracheobronchial trees[END_REF]) use commercial softwares, with Neumann or Dirichlet boundary conditions at the inlet and outlets, to perform their flow simulations. Nervertheless most "general purpose" softwares are not designed to take into account easily natural dissipative boundary conditions, nor the coupling with the spring model.

Time discretization

Theoretical results have been obtained in the space of divergence free functions. For the sake of efficiency, the numerical approach relies on a mixed formulation (non divergence free test functions are used, and the divergence free constraint is expressed in a weak form). Consequently, the formulation is obtained by multiplying [START_REF] Dauge | Elliptic boundary value problems on corner domains, Smoothness and asymptotics of solutions[END_REF] with general test functions v ∈ H 1 (Ω) that vanish on Γ ℓ . Thus, we obtain the variational formulation [START_REF] De Backer | Computational fluid dynamics can detect changes in airway resistance in asthmatics after acute bronchodilation[END_REF], with the unknown pressure p = p -P a , along wih the constraint equation Ω q∇ • u = 0, for all q ∈ L 2 (Ω), that is

                                     ρ Ω ∂ t u • v + ρ Ω (u • ∇)u • v + µ Ω ∇u : ∇v + N i=1 R i Γ i u • n Γ i v • n + m S 2 Γ 0 ∂ t u • n Γ 0 v • n + k S 2 t 0 Γ 0 u • n Γ 0 v • n - Ω p∇ • v = -P 0 Γ 0 v • n - f ext S Γ 0 v • n - k S 2 Sx 0 Γ 0 v • n ∀ v ∈ H 1 (Ω), v = 0 on Γ ℓ Ω q∇ • u = 0 ∀ q ∈ L 2 (Ω). (33)
Let δt > 0 be the time step, and t n = nδt, n ∈ N. We denote by u n the approximated solution at time t n . Using a semi-implicit scheme for the non-linear term in the variational formulation, and the following approximation for the time-integral term

t n 0 Γ 0 u • n in (33), t n 0 Γ 0 u • n ≈ n j=0 δt Γ 0 u j • n , (34)
the time discretization reads as follows

                                         ρ δt Ω u n • v + ρ Ω (u n • ∇)u n-1 • v + µ Ω ∇u n : ∇v + N i=1 R i Γ i u n • n Γ i v • n + m S 2 δt + kδt S 2 Γ 0 u n • n Γ 0 v • n - Ω p n ∇ • v =   -P 0 - f n ext S + k S x 0 - k S 2   n-1 j=0 δt Γ 0 u j • n     Γ 0 v • n + m S 2 δt Γ 0 u n-1 • n Γ 0 v • n + ρ δt Ω u n-1 • v ∀v ∈ H 1 (Ω), v = 0 on Γ ℓ Ω q∇ • u n = 0 ∀ q ∈ L 2 (Ω). (35)
             1 δt (u n , v) 0 + a 1 (u n , v) + δta 2 (u n , v) +b(u n , u n-1 , v) + c(v, p n ) = ℓ n (v) ∀ v ∈ H 1 (Ω), v = 0 on Γ ℓ c(u n , q) = 0 ∀q ∈ L 2 (Ω). (36)
where the right-hand side ℓ n (v) is given by

ℓ n (v) =   -P 0 - f n ext S + k S x 0 - k S 2   n-1 j=0 δt Γ 0 u j • n     Γ 0 v • n + 1 δt (u n-1 , v) 0 .
Remark 3.1 The discretization (34) of the integral term t n 0 Γ 0 u • n can be associated to an implicit discretization of the constraint equation S ẋ = -Γ 0 u • n (see Eq. ( 7)). Indeed, using backward difference, we have

x n = x n-1 - δt S Γ 0 u n • n, (37)
so, by downward recursion to initial position x 0 ,

x n = x 0 - δt S n j=0 Γ 0 u j • n, (38)
so that the mass position at time t n is implicitly calculated and can be determined using [START_REF] Vigneron | Espaces de Sobolev d'ordre variable : traces, éclatement, inégalité de Hardy[END_REF] or [START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modelling of blood and pressure in arteries[END_REF] after the resolution of the modified Navier-Stokes system [START_REF] Similowski | Two-compartment modelling of respiratory system mechanics at low frequencies: gas redistribution or tissue rheology?[END_REF].

In the same way, part of the right-hand side of the mixed variational formulation (36) corresponds to the following time discretization of the mass-spring differential equation ( 4)-( 5), coupled with the constraint equation ( 6) or [START_REF] Crooke | A general two-compartment model for mechanical ventilation[END_REF], using an implicit central difference scheme

P n a = 1 S m x n -2x n-1 + x n-2 δt 2 + kx n -f n ext ( 39
)
and using ( 37) and ( 38),

P n a = 1 S   - m δt Γ 0 u n • n - Γ 0 u n-1 • n + k   x 0 - δt S n j=0 Γ 0 u j • n   -f n ext   .
Remark 3.2 Note that another -equivalent-discretization strategy consists in consider the variational formulation obtained by using (8) as the starting system, instead of the system [START_REF] Dauge | Elliptic boundary value problems on corner domains, Smoothness and asymptotics of solutions[END_REF]. In this approach, the surface integral terms on Γ 0 would be written in terms of the sum of surface integrals on Γ i , that is in the semi-discretized variational formulation we would deal with a quadratic term of the form

m S 2 δt + kδt S 2 N i=1 Γ i u n • n N i=1 Γ i v • n .
Hence, in a finite element framework, this would imply that all the degrees of freedom belonging to the outlet boundaries Γ i , i = 1, . . . , N , would be coupled. The two formulations are equivalent as far as discrete solutions are concerned. Yet, the resulting matrices differ in the way degrees of freedom are coupled, and numerical tests suggest that the system resulting from formulation ( 35) is better conditionned.

Space discretization

The solver uses a P 1/P 1 stabilized finite element method (FEM) for the space discretization of Navier-Stokes systems. The linear systems obtained are then solved using a preconditioned GMRES iterative routine.

As seen before, the dissipative conditions coming from the resistance and the mass-spring timediscretization modify the bilinear forms in the variational formulation and, therefore, they would modify the FEM matrix associated to the velocity degrees of freedom. However, from the variational formulation we see that these modifications can be seen as an additive perturbation of the original FEM matrix,

A = A Ω + A res . ( 40
)
where A Ω stands for the standard FEM matrix over Ω, and A res stands for the matrix associated to resistances.

The matrix A res is not assembled : matrix-vector products involving A res are performed as tensor operations based on a vector whose entries depend on Γ j φ i • n (where φ i is a FEM velocity basis function), the resistance values R j , j = 1, . . . , N , and the modified resistance value R0 = R 0 + S -2 (mδt -1 + kδt) (for the entries associated to nodes on Γ 0 ). Remark 3.3 When considering only the Navier-Stokes system with natural dissipative boundary conditions (3) then one could impose these conditions explicitly. In this case the FEM matrice of the system is A Ω . Nervetheless for large resistance value R i this may be unstable, as suggested by the sufficient existence condition that one will obtain by omitting the resistive part in the definition of the bilinear form of the system (see Remark 2.13).

Remark 3.4 If the second variational formulation were used (see Remark 3.2), then the matrix product A res x would be more complicated, since all the outlets' degrees of freedom would be coupled. Moreover each entry associated to an outlet's degree of freedom would depend on the coefficient S -2 (mδt -1 + kδt) coming from the spring-mass system time discretization.

Note on the preconditioner

As the preconditioner for the modified linear system we use the incomplete LU (ILU) preconditioner associated to the original FEM matrices which demonstrates to be sufficient for the physiological numerical tests we performed. Nevertheless for large values of the subtrees resistance this preconditioner may not be appropriated.

Numerical results

In this section we present 3D numerical results of the coupled Navier-Stokes/mass-spring system, performed with an academic code developed at INRIA-Rocquencourt which has implemented the natural dissipative boundary conditions in an implicit way (see [START_REF] Maury | Outlet Dissipative conditions for air flow in the bronchial tree[END_REF]). We 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. We will also test the effect of a modification of the spring stiffness. The goal of this section is to illustrate the feasibility from a numerical point of view of the model using a realistic geometry and realistic data. The precise physiological exploitation of the multiscale system will be the object of a forthcoming paper.

Methodology, data

We have performed numerical calculations of the coupled Navier-Stokes/mass-spring system in a 3D reconstructed bronchial tree (until 6th -7th generation). The tetrahedral meshes (see Fig. 2) were constructed using INRIA GHS3D mesh generator [START_REF]GHS3D, tetrahedal mesh generator[END_REF]. The original 3D surface geometry is the same used in [START_REF] Fetita | An image based computational model of oscillatory flow in the proximal part of tracheobronchial trees[END_REF][START_REF] De Rochefort | In vitro validation of computational fluid dynamic simulation in human proximal airways with hyperpolarized 3 He magnetic resonance phase-contrast velocimetry[END_REF], and it was reconstructed, from CRT medical images, using a marching cube-base adaptive approach (see [START_REF] Fetita | An image based computational model of oscillatory flow in the proximal part of tracheobronchial trees[END_REF] for more details) .

Mesh size

We use a tetrahedral mesh composed of N nodes = 31288 nodes, N tetra = 131562 tetrahedra and N triang = 33950 surface triangles.

Resistance values

The resistance values were estimated using anatomical data from [START_REF] Weibel | Morphometry of the human lung[END_REF], assuming that each outlet Γ i is connected to a dyadic subtree which is assumed to follow a geometric decrease of the pipe sizes, and therefore a geometric increase of individual pipe resistances. As the outlets correspond to different generations (j 0 =4, 5, 6, 7, and 8), the corresponding resistances are computed as partial sums of the type

R = 17 j=j 0 r j 2 j-j 0 +1 , r j = 128µ π l j d 4 j , (41)
where l j and d j stand for the bronchus' length and diameter, respectively. Note that a dissipative boundary condition is prescribed at inlet Γ 0 , to account for the resistance of the upper part of the respiratory tract (nose, pharynx, and larynx). The corresponding resistance value is taken from [START_REF] Ben-Tal | Simplified models for gas exchange in the human lungs[END_REF][START_REF] Martin | Modelling of respiratory system mechanics involving gas exchange in the human lungs[END_REF].

The table below presents the resistance values (in g • cm -4 • s -1 ) for the inlet and for the 58 outlets (the value depends on the generation j 0 only). generation (j 0 ) R 0 1.12E+00 4 2.24E+00 5 3.89E+00 6 6.59E+00 7

1.10E+01 8

1.70E+01

Fluid and spring parameters

The fluid is supposed to be air saturated in water vapor, and its dynamic viscosity and density are given by µ

= 1.98 × 10 -4 g • cm -1 • s -1 , ρ = 1.11 × 10 -3 g • cm -3
so that the kinematic viscosity is given by ν = µ/ρ = 0.18 cm 2 • s -1 . For the mass-spring system, we have used the same values as in [START_REF] Martin | Modelling of respiratory system mechanics involving gas exchange in the human lungs[END_REF]:

m = 300 g, S = 100 cm 2 , k = 3.63 × 10 3 dyn • cm -1 ,

Results

We present the numerical results obtained in the case of normal respiration simulations and in the case of forced maneuvers simulations. In the case of normal respiration simulations we apply a periodic smoothed pulse force (of period T = 4 s) with f max = 1.5 × 10 5 dyn and f min = 0. Figure 3 shows the mass-spring displacement x (in cm) versus time t (in s) and air flow Q (in cm 3 • s -1 ) through boundary Γ 0 versus time t. Note that x max ≈ 4 cm, therefore the tidal volume V T is approximately equal to 400cm 3 , since S = 100 cm 2 , which is the average value for an healthy adult. Figures 4 and5 show the isovalues of the velocity vector field norm, interpolated on two cut planes immediately after the first bifucation (see Fig. 2), at peak inspiration (t = 0.4s), transition (t = 1.6s) and peak expiration (t = 1.9s).

We can observe that, in our experiments, the velocity profiles in these sections differ from inspiration and expiration and that during transition, as expected, more complex velocity profiles are formed in the upper airways. At the peak of inspiration we obtain a typical M-shape after the first bifurcation (see [START_REF] De Rochefort | In vitro validation of computational fluid dynamic simulation in human proximal airways with hyperpolarized 3 He magnetic resonance phase-contrast velocimetry[END_REF] where steady state simulations and experiments are performed on the same geometry).

Note that our numerical simulations are not able to capture the possible turbulent patterns of the flow since it will require a much more finer mesh and a parallel solver to perform DNS (Direct Navier-Stokes) simulations. Note moreover that in [START_REF] Croce | In vitro experiments and numerical simulations of airflow in realistic nasal airway geometry[END_REF], [START_REF] De Backer | Computational fluid dynamics can detect changes in airway resistance in asthmatics after acute bronchodilation[END_REF], [START_REF] Fetita | An image based computational model of oscillatory flow in the proximal part of tracheobronchial trees[END_REF] the authors assume that the flow is laminar, at least for the respiration at rest. It seems to be due to the oscilatory character of the respiration which prevents the development of the turbulence patterns. Nevertheless, the geometries used to perform the simulations in the previous papers (as in the present work) go from the trachea up to few generations of the bronchial tree. It is worth noticing that a recent study [START_REF] Lin | Characteristics of the turbulent laryngeal jet and its effect on airflow in the human intra-thoracic airways[END_REF], where steady state DNS simulations, with prescribed boundary conditions at the outlets and at the inlet, are performed, shows that for a geometry going from the mouth up to the to the 6th generation turbulence patterns appear. Whereas considering a geometry going from the trachea to the 6th generation, for a steady flow, the turbulence is negligible. The model we propose could provide a tool to determine whether or not turbulent patterns appear in a pulsatile regime, by performing DNS simulation in a full geometry including the mouth-larynx system, with physically relevant driving forces, and free inlet boundary conditions.

Parametric studies

In Figure 6 we compare the pressure field of the previous simulation at peak expiration with the case of a partially obstructed bronchial tree at the same time. The obstruction is obtained by increasing some resistance values (in this experiment, the resistances associated to the encircled outlets in Fig. 6 right). It shows that the qualitative behaviour of the flow is modified just by increasing few outlet resistances. Note that during an asthma crisis this is roughly what happens: the resistances of small airways in the distal part increase due to the inflammation process. Next we provide some experiments for forced respiration. In order to simulate the clinical experience of a normal subject performing forced maneuvers, we use a nonlinear constant k = k(x) as in [START_REF] Martin | Modelling of respiratory system mechanics involving gas exchange in the human lungs[END_REF]. Indeed, in this case one needs to take into account the fact that the motion of the lung parenchyma is limited by the chest wall in order to capture the limitation of the inhaled volume. More precisely:

k(x) = k 0 + (f min /x min -k 0 )x/x min , if x < 0 (f max /x max -k 0 )x/x max , if x ≥ 0.
with f min = -11 × 10 5 dyn and f max = 13 × 10 5 dyn and k 0 = 3.63 × 10 3 dyn • cm -1 . The external force belongs to the following intervals: The phase portrait (instantaneous flux vs. inhaled volume) is the curve obtained by spirometry, which provides pneumologists with information on the pathologies a patient may suffer. Figure 7 presents such a curve obtained by measurements, and Figure 8 the numerical results obtained after an appropriate tuning of the forcing term f ext .

f ext ∈ [-
The model allows to investigate the influence of pertubations on the phase diagram, and thus to quantify the observability of different parameters with respect to this phase diagram. In this spirit, we can compute and compare phase diagrams obtained for different values of k (see Fig. 9), and for different values of the distal resistance (Fig. 10). By these simulations we would like to emphasize that the coupled model we propose may be useful to reproduce some pathological features. A precise physiological exploration, in collaboration with lung specialists, still needs to be done.

Conclusions

The proposed well-posed multiscale model enables to describe the air flow behavior in the proximal part of the bronchial tree taking into account the parenchyma and diaphram motions. We have studied the numerical feasibility of the model on a 3D real geometry with realistic data. Note that in this 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 can expect to reproduce some aspects of pathological behaviors such as asthma (increase of the resistances) or emphysema (decrease of k) by fitting to phase portrait and look, then, at the modifications induced for the 3D flow. Nevertheless, we may have to modify the spring model to capture some of the non linear effects of the phenomenon and the physiological exploitation of this multiscale system will be the object of a forthcoming paper.

A Hardy inequality in a quarter space

We establish here a special Hardy inequality which is needed in the proof of Prop. 2.8. We consider the quarter space Ω = x = (x 1 , x 2 , x 3 ) = (x ′ , x 3 ) ∈ R 3 , x 1 ≥ 0 , x 2 ≥ 0 . , where D 2 u is the Hessian matrix of u, and |x ′ | = (x 2 1 + x 2 2 ) 1/2 is the distance to the x 3 -axis. The standard Hardy inequality (which involves u and its gradient) is also valid in the case fractional derivatives:

u |x ′ | δ L 2 (R 2 ) ≤= |D| δ u L 2 (R 2 ) , (42)
for any function u defined over R 2 , and δ ∈ (0, 1). The proof is given in [START_REF] Vigneron | Espaces de Sobolev d'ordre variable : traces, éclatement, inégalité de Hardy[END_REF][START_REF] Bahouri | Ingalités de Hardy précisées[END_REF]. b) Hardy inequality in R 3 , cylindrical setting We prove here that the previous inequality extends to R 3 , with a weight equal to the distance |x ′ | to the x 3 -axis. Let u denote a regular function defined over R 3 . With obvious notations, one has

R 3 |u| 2 |x ′ | 2δ = R R 2 |u| 2 |x ′ | 2δ dx ′ dx 3 R |D| δ u 2 L 2 ([x 3 =s]) ds = |D| δ u 2 L 2 (R 3 ) (43)
c) Hardy inequality with a weight in the right-hand side Let u and v be regular functions over R 3 . One has

|D| 2-δ u, v L 2 (R 3 ) = x ′ δ |D| 2 u, 1 |x ′ | δ |D| -δ v ≤ x ′ δ |D| 2 u L 2 (R 3 ) v L 2 (R 3 )

Figure 1 :

 1 Figure 1: Multiscale model Outlets Γ i , 1 ≤ i ≤ N of the upper part are coupled with Poiseuille flows, which are themselves coupled with the spring motion.The present paper continues and generalizes the multi-model strategies proposed, for air flow modelling, in[START_REF] Grandmont | A multiscale / multimodel approach of the respiration tree[END_REF] or[START_REF] Grandmont | Multiscale modelling of the respiratory tract: a theoretical framework[END_REF]. The same type of multiscale strategies has been developped and investigated for different applications such as blood flow simulations and several groups have successfully coupled three-dimensional models to either resistances or more sophisticated zero-dimensional models (lumped models) or one-dimensional models (see[START_REF] Formaggia | Multiscale modelling of the circulatory system: A preliminary analysis[END_REF][START_REF] Olufsen | Structured tree outflow condition for blood flow in larger systemic arteries[END_REF][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] Quarteroni | Coupling between lumped and distributed models for blood flow problems[END_REF][START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modelling of blood and pressure in arteries[END_REF]).

Remark 2 . 4

 24 and initial condition u(0) = u 0 . The flux term Γ 0 v • n in (15) can be defined for v in H by mean of the standard duality as

Using the bilinear forms a 1 (

 1 •, •) and a 2 (•, •) and the trilinear form b(•, •, •) introduced in Section 2.1, and the bilinear form c(•, •) defined by c(u, p) = -Ω p∇ • u the semi-discretized variational formulation can be written as follows

Figure 2 :

 2 Figure 2: Reconstructed bronchial tree. Cut planes (yellow) after first bifurcation.

Figure 3 :

 3 Figure 3: Mass-spring displacement x and air flow Q versus time.

Figure 4 :Figure 5 :

 45 Figure 4: Isovalues of the velocity vector field norm on cut plane on left first bifurcation, at peak inspiration (t = 0.4s, left), transition (t = 1.6s, middle), peak expiration (t = 1.9s, right)

Figure 6 :

 6 Figure 6: Pressure field in the standard (left) and obstructed (right) situations

1 .0 × 10 5 , 2 . 5 ×

 1525 10 5 ] in dyn (normal), f ext ∈ [-11 × 10 5 , 13 × 10 5 ] in dyn (forced).

Figure 7 :

 7 Figure 7: Measured Displacement (in cm)-Flow (in cm 3 • s -1 ) diagram.

Figure 8 :

 8 Figure 8: Displacement (in cm)-Flow (in cm 3 • s -1 ) diagram. For any u in L 2 (R 3 ), we denote by û its Fourier transform, we set |u| 2 s = R 3 |ξ| 2s |û(ξ)| 2 dξ and we denote by Ḣs (R 3 ) the set of functions for which this quantity is bounded. The homogeneous Sobolev space Ḣs over Ω is defined as Ḣs (Ω) = u |Ω , u ∈ Ḣs (R 3 ) , and it is endowed with the following semi-norm |u| Ḣs (Ω) = inf ũ|Ω =u |ũ| s . Finally, for v defined over R d , sufficiently regular, we define |D| η v by (|D| η v)(x) = (2π) -d

Figure 9 :

 9 Figure 9: Phase diagram for different values of k 0

Acknowledgements The authors would like to acknowledge the support of ACINIM project LeP-oumonVousDisJe, and to thank M. Dauge for fruitful discussions, M. Grasseau for his technical help with the mesh, and P. Gérard for his help in proving Proposition A.1.

supported by the ACINIM project Lepoumonvousdisje.

d) Hardy inequality in the quarter-space Let u be a regular function defined over the quarter space Ω = [x 1 ≥ 0, x 2 ≥ 0]. We assume that u vanishes on [x 2 = 0]. We extend u over R 3 in two steps: firstly, we define u 1 over [x 1 ≥ 0] as follows:

(up to a factor 2). We extend now u 1 over R by following a Babitch process:

This extension preserves again the L 2 norm of |x ′ | δ D 2 u 1 (up to a multiplicative constant). Inequality (44) can be applied to u 2 , which ends the proof.