Continuous and semi-discrete stability estimates for 3d/0d coupled systems modelling airflows and blood flows

,


Introduction
In the present work, we focus on the analysis and numerical analysis of geometric multiscale models used either for simulating physiological flows such as airflows in the respiratory tract, see e.g.[43,2,20,31,32,28,35] or blood flows in the arterial network, see e.g.[37,42,15,13,24,36].We aim, in particular, at giving a general theoretical framework to numerical observations with respect to numerical stability.
Simulations in patient-specific geometries may provide valuable informations to physicians to improve diagnosis, pulmonary drug delivery [30] or blood surgery [12].Nevertheless, direct simulations of 3D flows in geometries such as the tracheo-bronchial tree or the arterial network are limited by the following constraints: since the whole respiratory tree and the blood network are very complex, with a lot of bifurcations, and with different scales therein, numerical costs related to a full 3D simulation in the whole domain are prohibitive.Not to mention that the image processing of the complete bronchial tree or blood network is out of reach for the time being.Therefore the whole domain is usually truncated, restricting the computational domain to a smaller part which is considered to be the most significant one in terms of flow description at the global scale: the large bronchi for airflows or the aorta region for blood flows.As a countereffect, the removed part has to be taken into account thanks to suitable reduced models in order to describe the global behaviour of the whole system.
Air and blood are commonly modelled as homogeneous, viscous, Newtonian and incompressible fluids.Thus we consider a system of partial differential equations involving the Navier-Stokes equations, which has to be coupled to reduced models to take into account phenomena in the removed part of the domain.In this work, we focus on so-called 0D or Windkessel models that describe how the fluid flux and average pressure on the artificial boundaries is related to the mechanical properties of the truncated part.Note that 1D models (see for instance [14,16]) can also be considered to describe the reduced part.Here, the stability and numerical analysis of 3D/0D coupled systems is investigated, with special attention brought to applications related to, both, airflows and blood flows modelling, which involve different kinds of 0D models sharing nevertheless a similar formalism.The whole resulting system involves Navier-Stokes equations with nonlocal Neumann-type boundary conditions which depend on the chosen 0D model.
Many authors investigated the difficulties related to this kind of problems.From the theoretical point of view and the numerical point of view, one difficulty comes from the lack of energy estimate when considering the Navier-Stokes system with Neumann boundary conditions and more generally mixed Dirichlet-Neumann boundary conditions.Nevertheless existence of strong solutions (global in time for small data or local in time) has been shown in [27] under the assumption that the out/inlets meet the lateral boundary with a right angle and assuming some strong regularity results for the solution of the Stokes problem with mixed boundary conditions.Additionnally, when coupling the Navier-Stokes system with 0D reduced models, we refer to [39] and [22,2] for the same type of wellposedness results of strong solutions.In particular, in [39], the existence result based on a fixed point strategy, is obtained for the Navier-Stokes system with Robin-type conditions under a smallness assumption on the Robin coefficient modeling the resistive part of the 0D model.In [2] the Navier-Stokes system coupled to the resistive 0D model is considered; the regularity assumption that was previously used in [27] is dropped as well as the assumption on the resistance of the 0D model needed in [39].The proof relies on the regularity results for the solution of the Stokes system in polygonal domains with mixed boundary conditions, that have been derived in [34] and on the introduction of a well chosen Stokes-like operator, similar but simpler than the one we will introduce in this paper, that takes into account the 0D model.Note that all these results are valid for small data or in small time.Note moreover that concerning 3D/0D models (see [22]) and 3D/1D models (see [14,16,38]), coupling conditions based on the total pressure enables to have global energy estimates, allowing to prove global existence of weak solutions as in [22].From the numerical point of view, the lack of energy estimates for the Navier-Stokes system with Neumann boundary conditions or coupled with 0D model is linked to numerical instabilities as soon as, for a given physical setting, the applied pressure drop reach a threshold.To overcome this difficulty many strategies have been proposed.For general flows, we refer to the early work of [8] that introduces a whole set of boundary conditions that have been further extended and analysed in [9] where existence of weak solutions is proven.For hemodynamic flows, stabilization methods -some of them similar as the one proposed in [9] -have been introduced [3,24,1,4].These stabilization techniques lead to the modification of the resolved physical system.We also refer to [27,17,32], for reviews on these questions and to [5] where benchmark tests are performed for various stabilization methods.
Concerning the coupling of the Navier-Stokes equations with 0D models (or other reduced models such as 1D models), many strategies have been considered for the theoretical analysis of such problems as for their numerical implementation.For instance, as already stated, the existence result proved in [39] is based on a fixed point theorem, whereas in [2] a global formulation of the coupled system is considered.Moreover, for the modelling of physiological flows, many coupling strategies have been already implemented: the explicit treatment has been used for instance in hemodynamics in [42,37] or for airflows in [10] for 3D/1D coupling; the implicit coupling with Neumann boundary conditions involving the full traction, achieved thanks to an iterative process has been proposed in [24].Still in the context of hemodynamics the implicit coupling and the implicit treatment of the convective part of the fluid equations have been achieved by a Newton algorithm [13,29]; however, in each Newton sub-time-step, the coupling is explicit.The same strategy is used in [35] for mechanical ventilation in a rat bronchial tree.Moreover we refer for instance to [2] in the context of airflow modelling or [6] in the context of blood flows modelling where implicit monolithic coupling schemes are considered.The efficiency of the numerical methods associated to these problems relies on the analysis of two types of numerical difficulties: on the one hand, the explicit/implicit treatment of the nonlocal boundary conditions which couple the 3D and 0D models, which may lead to numerical instabilities and thus possible restrictions on the time-step, even with an implicit treatment of the coupling (when achieved by an iterative procedure).On the other hand, the more intrinsic difficulties coming from the convective term in the Navier-Stokes system which, as already stated, induces a lack of energy estimates and subsequent observed numerical instabilities.
Thus this paper is concerned with the analysis and numerical analysis of the coupled 3D/0D models arising in blood flows in large arteries as well as airflows in the bronchial tree.Here we derive energy or stability estimates for continuous and semi-discrete Stokes or Navier-Stokes system coupled to typical Windkessel models for explicit and implicit couplings, with a special emphasis to the dependance with respect to the physical parameters.The aim is to quantify, depending on the application field, the stability restrictions on the time step or on the data that are sufficient to ensure stability estimates.Note that this paper improves and generalizes partial results that can be found in [19,18].The outline of this article is as follows: in section 2 we introduce the coupled fluid-Windkessel models under study.We choose standard 0D models used in blood flow or airflow and we embed them in a similar formalism.We also present the semi-discretized in time schemes considered for the coupled fluid-Windkessel models, whose stability analysis will be performed in the next sections.Next, in section 3, we derive energy or stability estimates for the continuous or semi-discretized in time Stokes-Windkessel coupled sytem.Implicit and explicit coupling are considered.We exhibit different type of results depending on the considered 0D model.Finally in section 4 we study the Navier-Stokes-Windkessel coupled system.We prove stability estimates both in the continuous case and two semi-discrete cases with either explicit or implicit coupling.Since no energy estimate can be derived, we prove estimates in stronger norms linked to the domain of a new well-chosen Stokes-like operator adapted to the coupled system.Once again we exhibit different types of behavior depending on the considered 0D model as well as on the coupling strategy.
2 Reduced models in airflows and blood flows

Models
In this section, we describe different types of models associated to physiological flows, such as air through the bronchial tree or blood in the arterial network.The bronchial tree and the blood network have a complex structure which can be described as an assembly of tubes in which the biological fluid (air or blood) flows.For instance, the human respiratory tract is a dyadic tree of about 23 generations.The first generation (the trachea) has a length of about 10 centimetres, while the last one is about 1 millimeter.Until the 15th generation, the flow is convective whereas it is mainly diffusive in the acinar region.Moreover, the medical imaging and image processing techniques allow to obtain a good segmented surface and an associated mesh only up to the 6th or 7th generation.In the same way, the arterial network can be described as tube network.
In this context the complexity of the geometries makes it difficult to address direct simulations over the whole domain which then have to be truncated.Nevertheless, the removed parts corresponding to the smaller scales have to be taken into account in the global modelling: this can be done by defining appropriate reduced models.After truncation of the whole domain, we get a domain Ω ⊂ R 3 involving artificial boundaries which are denoted Γ i , with i ∈ {0, • • • , N }, N + 1 being the number of in/outlets.The lateral walls of the respiratory tree or of the aorta are denoted Γ ℓ .In these 3D domains, we assume that the velocity u and the pressure p of the fluid satisfy the following incompressible Stokes or Navier-Stokes system (corresponding respectively to ε = 0 and ε = 1): with u 0 the initial velocity, n the outward unit vector on every part of the boundary ∂Ω and ρ and µ the density and the viscosity of the fluid respectively.In order to model the whole system, i.e. the whole respiratory tree or the whole blood network, taking into account the fluid flow in the removed part, the 3D model has to be completed with a well-chosen reduced model.For instance, the removed part can be condensed into a 0D model (0D in the sense that it does not depend on a space variable) coupled to the 3D model at each outlet Γ i .Here we choose to consider some classical 0D models (also refered to as Windkessel models), used in blood or air flow modelling, but sharing the same formalism.We refer the reader to [40] for a review on geometric multiscale modeling in the context of cardiovascular systems.The coupling between the 3D and the 0D parts can be written as where p i is a constant in space interface pressure that depends on the considered 0D model.In all the studied cases, the 0D pressure is a function of the 0D fluid flux Q i , namely Moreover the mass conservation at the interface Γ i writes As a consequence the coupled system is Stokes or Navier-Stokes system with a generalized Neumann boundary conditions based upon the modelling of phenomena in the truncated part: Note that, as they involve the velocity flux at the artificial boundary, the boundary conditions are nonlocal.The choice of function F i depends on the application field as it is designed to mimick the behaviour of the truncated subtree.For instance in the context of blood networks models, a so-called RCR or LRCR model are used whereas, in the context of airflows in the bronchial tree, a so-called RC model is used.These models will be discussed thereafter.
Remark 1.In this work, several choice have been made: • We express the Neumann condition by using the non-symmetric tensor σ := µ∇u − p I. This choice can be justified by the fact that this quantity is continuous on cross section boundary in a Poiseuille flow in a cylindrical domain.An other choice could be based on the physical symmetric strain tensor σ symm.:= µ(∇u + ∇u t ) − p I and we refer to [27,17] for numerical comparisons between the two versions.Note nevertheless that the analysis performed hereafter remains unchanged when considering the full fluid strain tensor.
• Lateral walls in the 3D part are assumed to be fully rigid and, consequently, we impose the fluid velocity to be equal to zero on Γ ℓ .This assumption is valid in the context of air flows, at least for normal breathing.In the context of bloodflow modelling, the models should be enriched in order to take into account deformable walls: we refer to [21,40] for more sophisticated models involving a deformable domain.
Let us give some details on the considered reduced models.
• The RC model for air flows consists in reducing the truncated subtree into a resistive contribution and a compliant contribution plugged in series at each outlet of the 3D domain, see Figure 1: therefore we introduce a resistance R related to the resistance of the distal network (bronchial subtrees) and a compliance C describing the elasticity property of the surrounding tissues (lung parenchyma).The outlet pressure p i is associated to the current pressure P and the model reduces to the following ODE: As a straightforward consequence, • The RCR model for bloodflows consists in reducing the truncated subtree into a proximal part which is mainly resistive and a distal part which is resistive and compliant.The two parts are plugged in series at each outlet of the 3D domain, see Figure 2. Therefore we introduce a proximal resistance R p , a distal resistance R d and a compliance C, the latter representing the wall elasticity of the blood vessels.The outlet pressure p i is the current pressure P and the model reduces to the following ODE: As a straightforward consequence, As a straightforward consequence, In all these cases, it can be noticed that the proximal pressure that connects the 3D domain to the truncated 0D model only depends on the flux Q i so that function F i (•) takes the following general form: where coefficients α i ≥ 0, β i ≥ 0, γ i ≥ 0 and characteristic time τ i ∈ (0, +∞) are associated to corresponding models, and t → P i (t) is a source term which also depends on themodels.The coefficients α i model dissipation of the flux, β i represent inertia, γ i represent elastance of the 0D models with an associated relaxation time τ i .In particular, Table 2.1 summerizes the possible choices for these parameters related to the previous described models: The link between the different models can be described as follows • The four-element RCRL model with R d = +∞ leads to a so-called RCL model.
• The four-element RCRL model with L = 0 leads to the RCR model.
• The RCR model with R d = +∞ allows us to get the RC model with R = R p .
• The RC model with C = +∞ allows us to get the R model.
Remark 3. The compliance parameter in the RCR and LRCR represent the wall elasticity whereas we choose to consider rigid wall for the 3D part.We refer to [21,40] for more sophisticated models involving a deformable domain.

Variational formulation of the coupled system
Let us now write the variational formulation associated to the coupled problem.Define the following functional spaces Multiplying the first equation of system (1) by v ∈ V , integrating over the whole domain Ω and using integrations by parts with boundary conditions (5), we get: Notice that ε = 1 includes the full Navier-Stokes system whereas ε = 0 allows us to deal with the linear Stokes system.The variational formulation is essential for the derivation of energy estimates (see the forthcoming subsections) which may be obtained by considering v = u.The estimates are easily derived in the linear case but difficulties emerge in the nonlinear case because of the inertial effects.This difficulty is partially overcome by using several tools: construction of a suitable operator, normed space which is deeply associated to the inertia of the system and bilinear form that takes into account the fluid dissipation along with its 0D counterpart, namely the flux dissipation of 0D model.Let us moreover introduce a useful functional space and related property: In this space the following lemma holds true: Lemma 1.There exists κ > 0 such that, We refer the reader to [2] for the proof of this lemma.Note that this estimate is deeply based on the divergence-free property and on the fact that Γ i ∪ Γ j = ∅ for all i = j.Note that, for all v ∈ H, the flux has to be understood in a weak way: indeed it can be defined by means of duality: where g i is a function in H 1 (Ω) such that g i = 1 on Γ i and g i = 0 on Γ j , j = i.Note that such functions exist as the boundaries Γ i are not in contact.Finally, if v ∈ V, the classical flux formula is recovered.
Let us also introduce the following inequality, that can be deduced from the trace inequality: there exists C Γ > 0 such that In what follows, for the sake of simplicity (and without loss of generality for the mathematical analysis), we will consider two artificial boundaries: • one "inlet" Γ 0 with standard Neumann boundary condition, i.e. α 0 = 0, β 0 = 0, γ 0 = 0, P 0 := p 0 , where t → p 0 (t) is a prescribed pressure.
Moreover note that the non-homogeneous Neumann boundary condition on Γ 0 can be reduced to homogeneous Neumann boundary condition by defining a new unknown pressure p − p 0 (which will be still denoted p) and a new Windkessel source term P − p 0 (which will be still denoted P).Thus we will consider p 0 = 0, since the pressure drop is taken into account in the Windkessel source term P.
Consequently we consider the following problem:

Discretization schemes
We investigate the numerical stability of various coupling strategies between the Stokes or Navier-Stokes system and the 0D models.In particular, we aim at deriving stability estimates on the solution of the discretized-in-time and pay attention to the sensitivity of the stability constants or possible smallness conditions with respect to the physiological and numerical parameters.
In what follows, discretized-in-time systems will be referred as semi-discretized systems.Let ∆t > 0 be the time step and t n = n∆t, n ∈ {0, . . ., N }, with N ∆t = T .We denote by (u n , p n ) the approximated solution at time t n of the continuous velocity and pressure fields t → (u(t, •), p(t, •)).If we discretize in time the strong formulation of system (7), using the first order backward Euler scheme for the time derivative, the approximated velocity and pressure u n+1 and p n+1 satisfy: where the choice I ∈ {0, 1} corresponds to a semi-implicit or an implicit treatment of the convection term.Function F ∆t is a time approximation of F , where F is defined by (5).It depends on the approximations (Q k ) 0≤k≤m of the fluxes (Q(t k )) 0≤k≤m .Note that we may consider either explicit coupling with m = n or implicit coupling with m = n + 1.The approximate function is chosen as: Note that, here, the inertance term will be always treated in an implicit way.Indeed an explicit treatment of this added mass term leads to possibly unconditionally unstable schemes, as it has been observed and analyzed in [11] in the context of bloodflows.In what follows we investigate the case β = 0 and the case β = 0 in order to understand the stabilization effect of the inertance on the scheme.
Remark 4. The definition (9) of the approximate function modelling the 0D model is built upon the approximation of the RC and RCR models which leads us to the following quadrature formula Indeed in the case of the RC model (α > 0, β = 0, γ > 0, τ = +∞), the integral t 0 Q(s) ds is approximated by the classical rectangle rule whereas for the RCR model (α > 0, β = 0, γ > 0, τ < +∞), Equation (9) corresponds to the following discretization of system (2) We aim at studying the stability of the coupling schemes in both Stokes and Navier-Stokes regimes.In particular we focus on the derivation of stability conditions for the implicit and explicit schemes.The goal is to quantify the possible restrictions on the data and numerical parameters with respect to the physiological parameters and investigate the differences that can be encountered for various 0D models corresponding to bloodflows and respiratory flows.Consequently since our aim is to differentiate the behaviour of typical models used either in bloodflows or in air flows, we focus on the case γ > 0 in the forthcoming theorems.Note that the case γ = 0 will be treated in remarks pointing out the possible simplifications.In any cases, the variational formulation related to System (8) writes: where P m = P(t m ).In the above formulation, the choice of ε allows us to discuss the Stokes (ε = 0) and Navier-Stokes (ε = 1) cases; the index I corresponds to the implicit (I = 1) or semi-implicit (I = 0) treatment of the convection term; finally the index m allows us to describe an explicit (m = n) or implicit (m = n + 1) coupling between 3D and 0D models.
3 Study of the Stokes-Windkessel coupled system

Energy estimates for the continuous system
Let us now derive the energy estimates related to Problem (7) with ε = 0.This property is an important issue of the analysis of the numerical strategies which are built upon similar principles at the discrete level and is a key ingredient to prove existence of weak solutions.In particular, all the following calculations can be justified thanks to a Galerkin approximation, leading to a rigorous derivation of the existence of weak solution.
Theorem 1 (Energy estimates for the Stokes system).Let T > 0 and µ > 0. Assume that α ≥ 0, β ≥ 0, γ > 0 and 0 < τ ≤ +∞.Any weak solution u of the Problem (7) with ε = 0 satisfies, for 0 ≤ t ≤ T : Introducing the auxiliary volume we obtain easily that V satisfies the following ODE Note that the previous ODE ( 14) is also valid for τ = +∞ since, in this case, V (t) = t 0 Q(s) ds.From equations ( 13) and ( 14), we get Then we obtain Using a trace inequality and Young's inequality, the term P(t) Γ W u(t, •) • n can be controlled by: Finally using (16) to bound the right hand side of ( 15), we obtain We conclude by using Gronwall Lemma and remembering that V (0) = 0.
Remark 5 (Case γ = 0).In this case, the introduction of the auxiliary volume V is not necessary and the estimate ( 12) is still valid with γ = 0.
Remark 6.Note that other energy estimates can be derived.For instance, assume that α = 0. Then the estimate of the source term can be replaced by As a consequence (17) can be replaced by It leads to Note moreover that different estimates can be obtained for other cases: • If µ = 0 and α = 0, then the source term can be bounded as by Lemma 1 and Young's inequality.
• If β = 0, Nevertheless, with the previous inequalities (20) and (21), the energy estimate obtained by Gronwall lemma involves an exponential growth that behaves as e CT with C proportional to the density of the fluid or of the Windkessel model.To summarize, when γ > 0 and τ < +∞, in the case where the system is dissipative in u (µ > 0) or Q (α > 0) then, for zero applied pressures, the energy of the system is decreasing.When µ = α = 0 then the energy of the system is also bounded but with a bound that behaves as e CT , with C behaving like the inertia of either the fluid or the 0D model.
Remark 7. Note that when γ > 0 and τ < +∞, the auxilary volume V defined by (13) is "dissipated" by the model.In particular, in that respect, the RCR model (τ < +∞) and the RC model (τ = +∞) used respectively in blood flow simulations and air flow simulations behave in a different way, the auxilary volume V being dissipated in the RCR case whereas it is not in the RC case.

Energy estimates for the semi-discretized system
In this subsection, we establish energy estimates for the solution of the semi-discretized Stokes system with implicit coupling (see Theorem 2) or explicit coupling (see Theorem 3).Taking u n+1 as a test function in the variational formulation (11) with ε = 0 provides the following equality: We have and The discrete energy balance thus writes: where Defining the dimensionless parameter note that the discrete volume V n+1,n+1 (resp.V n+1,n ) is obtained by the rectangle rule with top-right (resp.top-left) corner approximation of the volume V (t n+1 ): We now distinguish the implicit and explicit cases in the following (resp.m = n + 1 and m = n).

Implicit coupling
Let us first consider the implicit case, namely m = n + 1.In the case of implicit coupling the analysis is nearly the same as in the continuous framework.More precisely the implicit coupling of the Stokes system with any 0D model leads to unconditionally stable schemes in standard energy norms.Denoting V k imp = V k,k , see (24), we easily verify that which corresponds to the time discretization of equation ( 14) by a backward Euler scheme.
Remark 8 (Case γ = 0).As for the continuous case, the same kind of discrete energy estimate for γ = 0 may be derived without introducing the discrete auxiliary volume.
Remark 9. Note that Theorem 2 can be extended to the fully-discretized system, using a standard Lagrange finite element approximation.Moreover different estimates of the source term as the estimates (18) and ( 21) obtained in Remark 6, can be adapted to the semi-discrete and fully-discrete frameworks with straightforward consequences on the global estimate (25).Furthermore in the semi-discrete framework, the estimate (20) is still true, and, in the fully-discrete framework, it is still valid but under some assumptions on the discrete finite element spaces, ensuring that Lemma 1 holds at the discrete finite element level.In particular, Lemma 1 will be satisfied for any and M h contains an internal approximation space of H 1 (Ω).

Explicit coupling
We recall that the β-term that accounts for the 0D model inertia is treated in an implicit way.Using the definition of the auxiliary volume, and denoting V k exp := V k,k−1 , see (24), we can derive discrete equations relating the flow to the discrete auxiliary volume: where δ ∆t is defined by (23).These are the key ingredients to prove: Theorem 3 (Explicit coupling with the Stokes system).Let µ > 0 and T > 0. Assume that α ≥ 0, β ≥ 0, γ > 0 and 0 < τ ≤ +∞.Let ∆t > 0 be the time step and t n = n∆t, n ∈ {0, . . ., N }, with N ∆t = T .The discrete solution u n+1 of Problem (11) with ε = 0 and m = n satisfies the following stability estimate • Under the condition where E 0 is a constant that only depends on the energy norm of the initial conditions and • Assume furthermore that β > 0.Then, under the condition where E 0 is a constant that only depends on the energy norm of the initial conditions and Proof.Considering the energy equality (22), with m = n, we obtain Let us deal with the explicit terms.The first one can be bounded simply as follows The second one, involving the discrete auxiliary volume V n+1 exp , can be rewritten as follows, using equation ( 27) and we get where we have used an estimate similar to (16) to bound the source term.The stability estimate is built upon the control of the following extra terms: Thanks to (27) we have Consequently we obtain and since δ ∆t < 1, we obtain Now we discuss two different cases: β ≥ 0 (general case), β > 0 (0D inertial case).
• In the general case β ≥ 0, and in particular if β = 0, the terms in the right-hand side of ( 29) cannot be controlled by the inertia of the 0D model.However they can be controlled by the inertial term of the fluid.Indeed by Lemma 1, , and, as a consequence, we deduce from (29) and the above inequality that Defining we introduce the roots of this polynomial with, for the sake of convenience, λ i > 0. Then we get As a consequence, we obtain Using the discrete Gronwall lemma [26] and under the condition 0 < ∆t < λ 1 , this provides the following stability estimate where E 0 is a constant that only depends on the energy norm of the initial conditions and • Assume now that β > 0. In that case the terms in the right-hand side of ( 29) can be controlled by the inertia of the 0D model.The estimate (29) yields Defining P (∆t , we introduce the roots of this polynomial with, for the sake of convenience, λi > 0. Then we get As a consequence, we obtain Under the condition 0 < ∆t < λ1 , the discrete Gronwall lemma [26], implies that where E 0 is a constant that only depends on the energy norm of the initial conditions and Note that alternate estimates can be derived following the continuous case that have been developed in Remark 6. Remark 10 (Case γ = 0).In this case, no auxiliary volume is required to derive discrete energy estimates.The sufficient conditions that guarantee the stability of the explicit scheme become: This condition involves the ratio of the inertance of the 3D or 0D system to the resistance of the 0D model.Moreover the exponential growth constants are modified as follows: Remark 11 (Influence of the inertia).When the inertia parameters of the problem, namely ρ and β, tend to +∞, so do the critical times λ 1 and λ1 which implies that in practice no condition on the time step is required to ensure stability.Moreover, the exponential growth remains bounded.Let us discuss the influence of the inertance parameter β on the critical time λ1 : When the inertance of the 0D model is small, so is the critical time λ1 ; nevertheless in this case it is sufficient to impose ∆t ≤ λ 1 that may be less restrictive to ensure the stability of the explicit scheme.
Remark 12. Let us discuss the influence of the characteristic relaxation time τ on the obtained stability estimates and smallness assumption on the time step: • the sufficient conditions on the time step do not depend on the parameter τ .
• when τ → +∞, the contribution e 2∆tT /τ 2 to the exponential growth goes to 1 and thus is uniformly bounded.
• when τ ≪ T , the exponential bound e 2∆tT /τ 2 hides an effective restriction on the time-step.Indeed in order to obtain a uniform bound of the exponential growth, this requires a severe restriction on the time step by choosing ∆t such that ∆tT /τ 2 = O(1), namely ∆t = O( τ 2 T ).
Remark 13 (Influence of the resistance parameter α).When α becomes large then the critical times λ 1 and λ1 behave as Thus the larger α is (which corresponds to the resistive parameter of the 0D model), the more severe is the constraint on the time step together with the exponential growth.
4 Study of the Navier-Stokes-Windkessel coupled system

Estimates for the continuous system
Let us consider the Navier-Stokes system and underline the standard difficulties met when one is interested in analyzing the energy balance when adding nonlinearities to the problem.Note that the estimates we will derive hereafter can enable us to prove existence of strong solution for small time or for small enough data.Let us first review the difficulties coming from the convection term.To fix the idea we consider the Navier-Stokes system (7) coupled with a R model ( Proceeding as in the linear case, we derive the following energy equality: Here we have used the coupling with the R model, as well as the divergence free property of the fluid velocity.We see a term ρ Γ W |u| 2 2 u • n that represents the flux of kinetic energy at the artificial boundary, whose sign is not known a priori.Consequently unlike for the Stokes system one can not derive easily an energy estimate.To obtain a satisfactory energy estimate and existence theorems, one has to be able to control this kinetic energy flux at the interface where Neumann boundary conditions are prescribed.Note that in dimension three we can prove the following bound (see [27]) which does not allow to obtain an energy estimate.Nevertheless existence of a unique strong solution can be proven.In particular, in [2], the existence of a unique strong solution (locally in time or for small data) is derived, based on the same ideas developped in [27] and on regularity results of the solution of the stationary Stokes system with mixed Dirichlet-Neumann boundary conditions in polyhedral domains [33].
Regarding the existence of solutions for the Navier-Stokes system with mixed Dirichlet-Neumann boundary conditions, we refer to [27]: the authors prove the existence of a unique smooth solution which is local in time; under an additional assumption on the smallness of the data, the smooth solution is proven to be global-in-time.Note that the existence of global weak solutions can be derived by choosing appropriate outflow boundary conditions that control the flux of incoming kinetic energy and thus stabilize the system [9].The case of Robin-type boundary conditions which involve the modelling of a local-in-space resistive contribution is analyzed in [39]: existence of a strong solution is obtained under the assumption that the resistance is small enough.In [2] a RC-like model is studied: the existence of a unique local-in-time strong solution for any data is proven; the particular case of a single R model is also investigated, leading to the existence and uniqueness of a global-in-time smooth solution for small data even if the resistance is large.Finally in [39], existence of a local-in-time strong solution for the RCR model is proven for small data.Proofs of the above results are all based upon Galerkin approximations with special bases.Note moreover that they all require that in/outlet meet the lateral boundary with right angles.This framework will be used in the analysis of the semidiscretized Navier-Stokes systems.We point out that the main difficulty in the above references relies on the estimate of the convective term of the Navier-Stokes system.
Let us now focus on the more general 0D model we study here and introduce key tools for the derivation of suitable estimates of the solution of the coupled system.In particular we introduce a new Stokes-like operator adapted to our coupled Navier-Stokes-Windkessel model.

Definition 1 (Stokes operator).
The space H is endowed with the scalar product and we denote • ρ,β the norm associated to this scalar product.Then we define the bilinear form a µ,α as a µ,α : Finally we introduce the operator A µ,α : D(A µ,α ) → H associated to the bilinear form a µ,α by Proposition 2 (Properties of the Stokes operator).The operator A µ,α has the following properties: ), H) is invertible and its inverse is compact on H; 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 [7] for general arguments and [23] for a direct proof in a similar context.

Lemma 3. The following estimates hold:
1.There exists L > 0 such that The constant L depends on the parameters as L := C P ρ + βκ 2 µ , where C P denotes the Poincaré constant.
2. If furthermore the artificial boundaries Γ 0 and Γ W meet the lateral boundaries Γ ℓ at angle π 2 and each boundary is smooth enough, then there exist ε > 0 and M > 0, such that The constant M depends on the parameters as: where Ω are constants which only depend on the domain.
Remark 14.Both constants L and M are proportional to the ratio of a density to a viscosity.
Proof.Both estimates rely on the properties of the Stokes operator.Using the definition of the scalar product on H, we obtain where we have used Lemma 1 and Poincaré inequality.Besides, by definition of the operator A µ,α , we have , which, by simplification, concludes the proof of (33).The proof of estimate ( 34) is based upon a regularity result for the Stokes problem with homogeneous mixed boundary conditions, see [33], for which we need the geometric angular assumption.The problem A µ,α v = f ∈ H can be rewritten as We consider the auxiliary pressure defined by By standard arguments and using Lemma 1 and equation ( 6), we have where C Ω is a constant which only depends on Ω.Using this auxilary pressure, the problem defined by ( 36) can be rewritten as As a consequence, from regularity results that can be found in [33] in the case of right angles, there exists ε > 0 such that where C Ω is a constant which only depends on Ω.Using estimate (37), we get By (33) and since f = A µ,α v, we obtain Ω C (1) which concludes the proof of (34).Now that we have introduced an appropriate Stokes-like operator adapted to our coupled Navier-Stokes-Windkessel system we can derive estimates in suitable norms for this system.The idea relies on the formal choice of u and A µ,α u as test functions and then a linear combination of the obtained inequalities.We will consider two cases.
• the so-called general case for which we prove an estimate valid for small time; • the so-called dissipative case, with τ < +∞, for which we prove that, if the initial data and applied pressures are small enough, an "energy" decrease can be established.
For any smooth solution u of the Problem (7), there exists T * > 0 and t → G r (t) which depend on the data such that: where r is any positive homogeneity constant.
• Dissipative case: global-in-time "energy" bound for small data.Assume furthermore that τ < +∞.Let δ > 0 and η > 0 such that There exists a dissipation parameter D > 0 (defined by (60)) such that, if the initial data and external forces are small enough, namely and then the solution satisfies a stability estimate: Moreover t → t 0 A µ,α u(s, •) 2 ρ,β ds is also bounded on any time interval [0, T ].Proof.Let us derive the estimates in the general case, then in the so-called dissipative case.Note that all the following formal calculations can be justified by using Galerkin approximation with a special basis associated to the eigenfunctions of the operator A µ,α , see Proposition 2.
Taking u as a test function in the variational formulation (7), we proceed as for the Stokes-Windkessel system and obtain estimate (17) with an additional term ρ Ω (u∇)u u that we bound as follows: where C Ω is a constant related to the continuous embedding of H 1 (Ω) onto L 4 (Ω).The estimate writes To control the last term in (41), we need to control u in L ∞ (0, T ; H 1 (Ω)).Thus we take A µ,α u as a test function in the variational formulation (7).By definition of the operator A µ,α , see Definition 1, we have The convection term can be estimated as follows Thanks to the continuous embedding of .
Then, choosing ε ′ < ε where ε is defined in Lemma 3, by a Hilbert interpolation combined with Lemma 3 there exists θ ∈ (0, 1) such that Consequently, Using Young's inequality, we get where δ > 0 will be chosen later on and C Ω,δ is a constant which depends on δ −1 and Ω only.Let us now deal with terms like γV (t)( Γ W A µ,α u(t, •) • n), for which we need a control of the auxillary volume V defined by (13).This control will be provided by estimate (41).By using Lemma 1, the definition of • ρ,β and Young's inequality, we have γV The linear forcing terms can be treated similarly: Thus, from (42) and thanks the previous estimates ( 43), ( 44), (45), we obtain We next choose δ > 0 such that 1 − 3δ > 0, for instance i.e. δ = 1 6 and we denote C Ω := C Ω, 1

6
. Using (47), we add ( 41) and ( 46) that has been multiplied by a positive homogeneity constant r to obtain: Thus we can apply a nonlinear Gronwall lemma by setting which, by (48) and since ∇u(t, •) 2 L 2 (Ω) ≤ 1 µ a µ,α (u(t, •), u(t, •)), satisfies the following inequality: Consequenlty, we obtain a stability estimate at least for a small time T * (depending on the data of the problem).From this bound on ϕ one can deduce that is also bounded on (0, T * ).
Next we further investigate the case τ < +∞.In this case, as already underlined for the Stokes system, the auxiliary volume V , defined by (13), is dissipated by the system.We take advantage of this to derive a stability estimate for any time but for small enough data.When taking u as a test function in the variational formulation (7), we bound the convective term in a coarser way than we did previously: using inequality (33) in Lemma 3 and the definition of • ρ,β we get where, as in inequality ( 40), constant C Ω is related to the continuous embedding of H 1 (Ω) onto L 4 (Ω).The estimate now writes Next we take A µ,α u as a test function in the variational formulation (7).First the convective term can be bounded as follows where we have used the continuity of the embedding H where η > 0 will be chosen later on.The linear forcing terms are bounded as in the general case, see (45).Thus, by ( 42), (51), ( 52) and (45), we obtain By choosing δ and η sufficiently small such that we obtain We multiply (55) by η and add (50) to obtain: Recalling the definition of a µ,α (•, •) (see Definition 1) we get Next we want to make appear some dissipation of even in the case α = 0 and β > 0. Since |Q| ≤ C Γ ∇u L 2 (Ω) , we have that and thus Then, since u L 2 (Ω) ≤ C P ∇u L 2 (Ω) by Poincaré inequality, we thus obtain Estimate ( 56) can be rewritten as with and a dissipation coefficient defined as The constant D stands for the dissipation of the system.In particular, assuming that τ < +∞ ensures that the 0D model is indeed dissipative with respect to the volume V if the data are small enough.Assuming that A − B ∇u(t, •) L 2 (Ω) ≥ 0, we obtain by Gronwall lemma: Consequently if which implies that ∇u 0 L 2 (Ω) ≤ A 2B , then we obtain the following estimate for all time: To conclude, provided that the initial data and the source term are small enough (see conditions (62)), the solution satisfies a stability estimate in suitable norms, namely Finally, combining this estimate with (57) allows us to conclude that Remark 15.This theorem shows as for the Stokes-Windkessel system that typical 0D models used in blood flows and in airflows modeling behave in a different way.Indeed for the RCR model for instance we obtain global-in-time estimates provided the data are small enough whereas for the RC model we obtain only local-in-time estimates.
Remark 16.In the dissipative case, parameters δ and η may be specified, for instance as follows: As a consequence, updating the constants and source terms as we deduce that, if the estimate now reads: Note that the behaviour of the required bound on the initial velocity, namely A B , as well as the upper bound in the estimate (39), namely can be described more precisely with respect to the parameters: A 2B and E tend to 0 when µ → 0, ρ → +∞, β → +∞, γ → +∞ or τ → +∞.In particular the less the 0D model dissipates energy with respect to the auxilary volume V , the more restrictive is the smallness condition on the data.
Remark 17 (Case : γ = 0).Note that, in that case, an estimate can be derived by choosing only A µ,α u as a test function.Indeed no control on the auxiliary volume V (t) defined by (13), is required.More precisely, we easily derive the following estimate In this case the dissipation comes from the term Au(t, •) 2 ρ,β .Assuming that and since Au(t, Thus the analogue of the dissipation coefficient D is, in this case, Consequently, Proceeding as in the proof of Theorem 4, we obtain that if the initial data are small, namely and if the external pressure satisfies then the solution satisfies the following stability estimate Here we recover the results obtained in [27] for the Navier-Stokes system with Neumann boundary conditions as a particular case of this result by setting α = β = 0. Note also that the behaviour of the required bound on the initial velocity, namely 1 M , as well as the upper bound µ M 2 in estimate (65) can be described more precisely with respect to the parameters: D 0 , 1 M and µ M 2 tend to 0 when µ → 0, ρ → +∞, β → +∞.Remark 18.We could have kept (41) instead of (50) and, following nearly the same lines, obtain a stability estimate for small enough data.But it leads to slightly more tedious calculations we choose not to present here for the sake of simplicity.
Remark 19.In the case β > 0, estimates (45) and (52) can be adapted by using the property In that way similar estimates can be derived by using the 0D inertia instead of the fluid intertia, leading to a possibly less restrictive condition on the smallness assumption on the data.In particular adding inertia in the 0D model may stabilize the whole coupled system.

Implicit coupling
We now focus on the semi-discrete Navier-Stokes system implicitly coupled to the Windkessel model which corresponds to the variational formulation (11) with ε = 1 and m = n + 1.We consider a semi-implicit treatment of the convective term, namely I = 0.
• Dissipative case: global-in-time bound for small data.Assume that τ < +∞.Let us consider the constants δ, η, D and the function H as in Theorem 4. If the initial data and external forces are small enough, namely and then the solution of (11) with ε = 1, m = n + 1 and I = 0 satisfies the following estimate: for all n ∈ {0, ..., N } Moreover ∆t N n=0 A µ,α u n 2 ρ,β is also bounded independently on N .• Case τ = +∞: local-in-time bound for small data.Assume that the time step is such that where r is a positive homogeneity constant.Assume furthermore that the initial data, external forces and final time T = N ∆t satisfy then Proof.All the forthcoming calculations can be rigorously justified by using a Galerkin method with a special basis associated to the Stokes-like operator A µ,α , see Proposition 2. Moreover existence of a solution (for which uniqueness could be also proven) can be also derived thanks to the previous estimates by the same Galerkin approximation.
We consider the system (8) with I = 0 for which the convection term is semi-explicit.We prove estimate (66) by induction.We follow the steps of the continuous case by taking u n+1 and A µ,α u n+1 as test functions in the variational formulation (11).Note that the only difference with the continuous case concerns the estimate of the convection term.The discrete analogue of (49) should be read as The discrete stability estimate can thus be written as follows We take A µ,α u n+1 as a test function in the variational formulation (11) for ε = 1, I = 0 and m = n + 1: ρ By definition of operator A µ,α , we have ρ Using again the definition of A µ,α , we have also Now we have to estimate ρ Ω (u n ∇)u n+1 A µ,α u n+1 .We do not follow exactly the same lines as for the continuous case in particular since the convective term is treated in a semi-implicit way: where C Ω comes from the continuous embeddings H 3 2 +ε (Ω) ֒→ W 1,3 (Ω) and H 1 (Ω) ֒→ L 6 (Ω).Then using (45) for the forcing term, (52) for the term involving the volume and estimate (69) for the convection term, the discrete analogue of estimate (55) reads with (δ, η) satisfying (54).By multiplying estimate (70) by η and adding (68), we obtain Recalling that A is defined by (58) and defining B by since , and u n+1 L 2 (Ω) ≤ C P ∇u n+1 L 2 (Ω) by Poincaré inequality, we proceed as in the continuous case in order to derive the discrete analogue of estimate (57): Recalling the definition of D, see (60), and introducing the approximation of Ψ(t n ) defined by we obtain where H n+1 = H(t n+1 ), H being defined by (59).Then, if , for all k ∈ {0, ..., n}, we get Assuming that we conclude that , and Consequently, by induction, we can prove that the solution stays at each time iteration in the same ball defined by estimate (75).
In this case we cannot take advantage of the dissipitation with respect to the volume V n imp .Taking u n+1 as a test function, we obtain Next, once again we take A µ,α u n+1 as a test function.Here we have to control the term without using the dissipative term γ∆t τ (V n+1 imp ) 2 of (68) that is equal to zero in the case τ = +∞.We have, by using Lemma 1 and the definition of • ρ,β : for δ > 0 that will be chosen latter.Remembering estimate (45) of the linear forcing terms and estimate (69) of the convection term, we obtain Ω M∆t ∇u n L 2 (Ω) A µ,α u n+1 2 ρ,β .
Consequently by choosing 2δ = 1 2 , we have Multiplying (78) by a homogeneity coefficient r and adding (76) yields Thus, if we impose ∆t < ρ 8rκ 2 γ = ∆t r , and assuming that Ω M , ∀k ∈ {0, ..., n}, we obtain, thanks to discrete Gronwall Lemma, the following discrete stability estimate Consequently to satisfy (80) for k = n + 1 and obtain the desired result by induction, the data have to verify The above condition requires that the initial conditions, the forcing term as well as the final time T are small enough.
Note that a standard fixed point argument [41,27] allows us to obtain the same kind of stability bound for the solution of system (11) with I = 1 together with the existence of a strong solution.
Remark 20.Let us comment the dependency on the homogeneity constant r.The upper bounds defined by (80) and (82) are increasing with respect to r, are equal to zero for r = 0 and have a finite limit as r goes to +∞.Moreover, at the same time, the critical time step ∆t r goes to zero as r goes to +∞ and so does the exponential growth.Thus, large values for r induce restrictive smallness conditions on the time step and on the data.For small values for r the condition on the time step is dropped but the upper bound on the data goes to 0. Remark 21.We can already note that for typical Windkessel model in blood flow for which τ < +∞, we can ensure global-in-time stability of the semi-discrete solution for small enough initial data and external forces, whereas for typical Windkessel model in airflow for which τ = +∞ we exhibit restrictive sufficient conditions on the time step and on the data (initial data, external forces, final time).In this case, the restriction on the time step writes ∆t < ρ 8rκ 2 γ .In particular the smaller the Windkessel compliance is, the more restrictive the conditions on the data and the time step are.Thus when ρ goes to zero or when γ goes to infinity the time step goes to zero.
Remark 22.Note that we cannot extend easily the proof of Theorem 5 to the fully-discretized system, unlike for the Stokes system.To do so one should introduce the finite element discrete analogue of A µ,α , see [25] in which this type of analysis is done for the Navier-Stokes system with homogeneous Dirichlet boundary conditions.
Remark 23 (Case γ = 0).As for the continuous case, it is sufficient to choose A µ,α u n+1 as a test function in order to derive a stability estimate, as in Remark 17.In this case the system is also dissipative and no condition on the time step is required for the stability.

Explicit coupling
We now focus on the semi-discretized Navier-Stokes-Windkessel model with explicit coupling.More precisely, as for the Stokes-Windkessel model with explicit coupling, the inertia of the 0D model associated to the parameter β is treated implicitly whereas the terms related to the parameters α and γ are treated explicitly.Moreover we consider a semi-implicit treatment of the convective term.It thus corresponds to the problem (11) with ε = 1, m = n and I = 0. We have all the ingredients to study this case, since it will be a mix of the study done for the Stokes system with explicit coupling and the study of the Navier-Stokes system.Theorem 6.Let µ > 0. Assume that the artificial boundaries Γ 0 and Γ W meet the lateral boundaries Γ ℓ at angle π 2 and that each boundary is smooth enough.Assume that α ≥ 0, β ≥ 0, γ > 0 and 0 < τ ≤ +∞.The discrete solution of (11) with ε = 1, I = 0 and m = n satisfies a discrete stability estimate, under restrictive condition on the data (final time, initial data and external forces) and on the time step.More precisely, let r be a positive homogeneity constant, assuming that ∆t < λ 1 (where λ 1 is defined by (31)) and that (the constant C S ∆t being defined in Theorem 3) and then the following discrete estimate holds true Note that here the discrete estimates are derived on the system with a semi-implicit treatment of the convection term.Nevertheless a fixed point procedure could enable to prove a similar estimate in the case of an implicit treatment of this term, namely in the case I = 1.
Proof.By taking u n+1 as a test function in the variational formulation (11) with ε = 1, I = 0 and m = n following the same lines as in the proof of Theorem 3 in the case β ≥ 0 and using moreover the following estimate of the convection term where δ ∆t have been previously defined in the proof of Theorem 3 (see (23)).As the coupling is explicit the next test function to consider is A µ,0 u n+1 .Note that the constant L appearing in the previous estimate (85) does not depend on the parameter α.Following the same lines as in the proof of Theorem 5 leads to Here the constant M 0 is associated to the operator A µ,0 and thus is defined by (35) with α = 0 (see (84)).We next have to estimate the last two terms of the right hand side of (86) corresponding to the explicit coupling.Let us first consider α∆t Q n Γ W A µ,0 u n+1 • n .By Lemma 1 and Young's inequality, we have where δ > 0 will be chosen later.Next we estimate γ∆t V n+1 exp Γ W A µ,0 u n+1 • n using once again Lemma 1: At this stage we do not distinghish two cases as we did previoulsly.In the general case we can not take advantage of the dissipation of the volume, thus we estimate γ∆t V n+1 exp Γ W A µ,0 u n+1 • n as we did in the implicit coupling, see (77): Thus from (86) multiplied by a homogeneity coefficient r > 0, using (87), choosing δ = 1 4 , and adding (85) we obtain Next remembering the definition of the polynomial function P (see (30)) and its positive root λ 1 (see (31)) and using the lower bound (32), we have Thus assuming that Ω L 2 + rC Ω M 0 , ∀k ∈ {0, ..., n}, and if we moreover impose ∆t < λ 1 , we obtain thanks to the discrete Gronwall Lemma that with C NS ∆t,r defined by (83).Consequently if the data satisfy Let us define B 0 by Ω L 2 + C Ω ηM 0 . (95) Note that B 0 corresponds to the constant B defined by (72) with α = 0. Next remembering the definition of the polynomial function P (see (30)) and its positive root λ 1 (see (31)) and using the the lower bound (32), remembering also the definitions of A defined by (58) and of H defined by (59), we have As a conclusion, in the case of the semi-discrete Navier-Stokes system coupled explicitly to the 0D model, the obtained stability estimates are quite similar in both considered cases.In both proofs, we exhibit a similar sufficient condition on the time step (which is the same as for the Stokes system) and restrictive assumptions on the data imposing smallness of the initial data as well as on the applied forces but also on the global time T .In particular in the so called dissipative case when explicity coupled to the Navier-Stokes equation the system does not dissipate energy anymore.
Remark 25.As in the proof of Theorem 3 we can adapt the previous calculations to the case where β > 0 and take advantage of the inertia of the 0D model.For the sake of simplicity we do not reproduce the calculations here.
Remark 26.For the case where τ < +∞ in order to control the exponential growth e 2∆t τ 2 T one could impose to the time step to satisfy ∆t ≤ 2τ 2 T which could be a rather restrictive condition in particular for large time T or small relaxation time τ .
Remark 27 (Case γ = 0).Once again it is sufficient to choose A µ,α u n+1 as a test function.Nevertheless, since no control on the auxiliary volume is required, no restriction on the time step emerges from the subsequent analysis.Yet the stability estimate involves an exponential growth coming from the explicit treatment of the resistive 0D model.Thus a sufficient smallness condition on the initial data, external forces and final time guarantees the stability of the discrete solution.

Conclusion
In this paper we derived energy or stability estimates for Stokes-Windkessel and Navier-Stokes-Windkessel models, both in the continuous setting and in the semi-discrete one with implicit or explicit coupling.One of the key ingredients in the derivation is the control of an auxiliary volume associated to the 0D model.This property is obtained by taking the velocity field u as a test function in the variational formulation for both Stokes and Navier-Stokes regimes.Nevertheless in the Navier-Stokes regime, it is not sufficient to obtain a stability estimate; in this case we have to consider a combination of the previous estimate with another one in which A µ,α u is used as a test function, where A µ,α is a new Stokes-like operator adapted to our coupled system.This enables us to control both the convective term and the auxiliary volume.
In particular, we have shown that for the standard Windkessel model used in bloodflow modeling, namely the RCR model (for which α > 0, β = 0, γ > 0 and 0 < τ < +∞), "energy" dissipation holds true and can be also derived when considering the Navier-Stokes system under smallness assumptions on the data.Meanwhile the standard Windkessel model used in airflow modeling (for which α > 0, β = 0, γ > 0 and τ = +∞) is not dissipative, leading to restrictive conditions -on the time step, on the data, on the final time -for the stability of the Navier-Stokes system even for the implicit coupling; this echoes the numerical simulations of such coupled problems for which instabilities arise leading to the use of stabilization techniques that make it possible to get away from smallness conditions on the data or, at least, give access to a larger range of data, see [5] where benchmark tests are performed for various stabilization methods.When considering an explicit coupling for the Navier-Stokes system, in both cases (RCR and RC models), "energy" dissipation does not hold anymore; stability estimates are derived but they require the same type of restrictive conditions.
We also paid attention to the dependency on the various physical parameters.Even if some of these constants are not explicit and depend on the geometry (such as the Poincaré constant, for instance) and thus on the considered test case, the derived estimates and their related validity conditions give a good insight into the behaviour of coupled systems according to the underlying fluid involved.For instance increasing the inertance in the 0D model (namely taking greater values for β) which is always treated in an implicit way stabilizes the semi-discretized system for the Stokes regime whereas it leads to a more restrictive smallness condition on the data in the Navier-Stokes regime.
All these theoretical achievements shall be further investigated in a forthcoming paper thanks to numerical experiments in order to quantify how sharp the derived estimates are, which range of data (for instance, applied pressures) enable stable simulations for the Navier-Stokes-Windkessel system without stabilization technique, which method has the best performance for each application field.

Figure 1 :
Figure 1: Reduced 0D model: the RC model for air flows

3 2
+ε (Ω) ֒→ L ∞ (Ω) together with estimate(34) of Lemma 3 and the definition of • ρ,β .Let us now deal with the term γV ( Γ W A µ,α u • n).Using Lemma 1, Young's inequality and the definition of • ρ,β , and taking advantage of τ < +∞, leads to Remark 24.Let us discuss the so called dissipative case, namely τ < +∞ for which we could have tried to take advantage of the volume dissipation.In this case we can reproduce (52) to obtain Thus by (87) and (92), the bound (86) becomes By choosing δ and η satisfying (54), the previous estimate (93) writes Consequently the desired discrete stability estimate can be proven by induction assuming that the data of the problem satisfy n ) 2 + ∆tH n .Thus, if ∇u k L 2 (Ω) ≤ A 2B 0 ,for all k ∈ {0, ..., n}, and assuming that ∆t < λ 1 , we obtain thanks to the