Cristóbal Bertoglio
email: cristobal.bertoglio@inria.fr

Alfonso Caiazzo
email: caiazzo@wias-berlin.de

Miguel A Fernández
email: miguel.fernandez@inria.fr

Fractional-step schemes for the coupling of distributed and lumped models in hemodynamics

Keywords: Blood flows, Chorin-Temam projection method, fluid-structure interaction, lumped parameter models, 3D-0D coupling

In three-dimensional (3D) blood flow simulations, lumped parameter models (0D) are often used to model the neglected parts of the downstream circulatory system. We analyze two 3D-0D coupling approaches in which a fractional-step projection scheme is used in the fluid. Our analysis shows that explicit approaches might yield numerical instabilities, particularly in the case of realistic geometries with multiple outlets. We introduce and analyze an implicitly 3D-0D coupled formulation with enhanced stability properties and which requires a negligible additional computational cost. Furthermore, we also address the extension of these methods to fluid-structure interaction problems. The theoretical stability results are confirmed by meaningful numerical experiments in patient specific geometries coming from medical imaging.

Schémas à pas fractionnaire pour le couplage de modèles distribués et localisées en hémodynamique

Résumé : Les modèles à paramètres localisés (0D) sont souvent utilisés dans les simulations tri-dimensionnelles (3D) d'écoulements sanguins pour modéliser l'effet des parties de l'appareil circulatoire négligées en aval. Dans cet article nous analysons deux approches pour le couplage de modèles 3D-0D, dans lesquelles un schéma à pas fractionnaire est utilisé dans le fluide. Notre analyse met en évidence que les schémas de couplage 3D-0D explicites peuvent donner des instabilités numériques dans le cas de sorties multiples. Nous introduisons et analysons une formulation 3D-0D implicite avec des propriétés de stabilité améliorées et un coût de calcul supplémentaire négligeable par rapport au cas explicite. Nous abordons également l'extension de ces méthodes à des problèmes d'interaction fluide-structure. Les résultats théoriques de stabilité sont confirmés par des expériences numériques dans des géométries et données réalistes.

Mots-clés : Écoulements sanguins, schéma de projection de Chorin-Temam, interaction fluidestructure, modèles à paramètres localisés, couplage 3D-0D.

Introduction

In 3D distributed models of blood flow (e.g., Navier-Stokes equations, fluid-structure interaction), downstream pressure boundary conditions are often used to represent the effect of the neglected portion of the vessels. In fact, since pressure measurements are invasive and not always available, the downstream circulation is usually modeled trough lumped parameter (or 0D) models. This results in a set of algebraic-differential equations, relating fluxes and pressures at each outlet boundary of the fluid domain (see, e.g., [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF][START_REF] Quarteroni | Coupling between lumped and distributed models for blood flow problems[END_REF][START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF][START_REF] Vignon-Clementel | Outflow boundary conditions for 3d simulations of non-periodic blood flow and pressure fields in deformable arteries[END_REF] and the references therein). A widely used lumped parameter model is the so-called three-element Windkessel system (see, e.g., [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF][START_REF] Frank | Die Grundform des arteriellen Pulses[END_REF])

The 3D-0D coupling between distributed and lumped models is operated by interface conditions that guarantee the continuity of the fluxes and pressures on the outlet boundaries (3D-0D interfaces). It is well-known that, when this coupling is treated explicitly in time, the whole system might suffer from numerical instability, regardless of the solution scheme in the 3D compartment. This enforces restrictions on the time-step length, that can considerably increase computational cost in realistic applications (see, e.g., the discussion in [START_REF] Quarteroni | Coupling between lumped and distributed models for blood flow problems[END_REF]). On the other hand, implicit 3D-0D coupling schemes overcome this instability issue at the expense of solving a fully coupled 3D-0D system at each time step. The so-called partitioned approaches solve this system by sub-iterating between the 3D and 0D models, which might be inefficient in practice (see, e.g., [START_REF] Malossi | Algorithms for the partitioned solution of weakly coupled fluid models for cardiovascular flows[END_REF]). Alternatively, monolithic procedures solve both models simultaneously and yield system matrices with a modified sparsity pattern, which might lead to preconditioning issues (see, e.g., [START_REF] Moghadam | A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations[END_REF]).

This work is devoted to the formulation and analysis of 3D-0D coupling schemes based on a fractional-step projection time-marching of the fluid (see, e.g., [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Temam | Sur l'approximation de la solution des equations de Navier-Stokes par la méthode des pas fractionaires I[END_REF]). First, an advectiondiffusion problem is solved to recover an intermediate approximation of the velocity field and, then, a suitable pressure field is recovered by solving a Poisson problem. We also consider the case of a 3D fluid-structure interaction modeling, time-discretized by the projection semiimplicit coupling scheme reported in [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]. In this fractional-step framework, the explicit and implicit treatments of the 3D-0D coupling lead to different formulations of the pressure-Poisson projection step. It can be discretized in an explicit fashion by time-advancing the 0D model from the previously computed intermediate velocity flux, which provides an explicit Dirichlet boundary data for the 3D pressure-Poisson problem. For a purely fluid 3D distributed modeling with two or more 3D-0D interfaces, our analysis shows that this approach might compromise stability. We also show that, in the case of fluid-structure interaction, numerical instabilities might appear even with a sole downstream boundary. We propose to overcome these instability issues through the introduction of a implicit 3D-0D formulation with enhanced energy balance across the 3D-0D interface, both for a fluid and a fluid-structure interaction modeling of blood flow. A salient feature of the proposed schemes is that they preserve the two-step velocity/pressure splitting of the original fractional-step scheme, characterized by a 3D-0D coupled problem of reduced size. Moreover, we show that this coupling strategy can be efficiently implemented by considering a single unknown per 3D-0D interface. This yields a computational complexity comparable to an explicit scheme and, hence, overcomes the above mentioned stability/complexity issues.

The rest of the paper is organized as follows. In Section 2 we introduce the 3D fluid equations, its corresponding time discretization (via a fractional-step projection scheme) and the considered Windkessel model. Section 3 is devoted to the formulation and analysis of explicit and implicit 3D-0D coupling schemes with a 3D distributed model based on the Navier-Stokes equations, while Section 4 consider the fluid-structure interaction case. In section 5 we present and discuss the performed numerical experiments, that is, the Navier-Stokes flow in a realistic aortic geometry and with measured clinical data, and the fluid-structure interaction in an ideal abdominal aortic aneurism. At last, Section 6 draws some concluding remarks and lines of future work.

The 3D-0D model problem

We first summarize the main ingredients of a 3D-0D model of blood flow coupling the incompressible Navier-Stokes equation and a lumped parameter description of the downstream boundaries. We then present the time-marching schemes considered in each sub-system and introduce some notation for the spatial discretization of the fluid equations.

3D fluid equations

We consider a domain Ω f ⊂ R 3 with the following partition of its boundary

∂Ω f def = Γ in ∪ Σ ∪ Γ out .
In the context of blood flow simulations, Ω f will represent the lumen of the vessel (see Figure 1), with Γ in , Σ and Γ out denoting, respectively, the inlet, vessel wall and outlet boundaries. We

Ω f Σ Γ in Γ 1 Γ 2 π1 Rd,1 C1 Rp,1 C2 Rp,2
Rd,2 π2

Figure 1: Sketch of the fluid domain Ω f with two outlet boundaries

Γ out = Γ 1 ∪ Γ 2 , (n 0D = 2).
now consider the incompressible Navier-Stokes equations for the velocity u :

Ω f × R + → R 3 and the pressure p : Ω f × R + → R:              ρ f ∂u ∂t + ρ f u • ∇u -∇ • σ(u, p) = 0 in Ω f , ∇ • u = 0 in Ω f , u = u in on Γ in , u = 0 on Σ, (1)
where ρ f stands for the density of the fluid and the fluid Cauchy-stress tensor is given by

σ(u, p) def = -pI + 2µǫ(u) and ǫ(u) def = 1 2
∇u + ∇u T , µ being the dynamic viscosity of the fluid and u in being a given inlet velocity field. The additional boundary conditions on Γ out will be considered in the next subsection.

Inria

0D Windkessel model

In what follows, the oulet boundary Γ out is assumed to be made of n 0D components

Γ out = n0D l=1 Γ l ,
such that Γ i ∩Γ j = ∅ for i, j = 1, . . . , n 0D with i = j (see, e.g., Figure 1). We will consider a threeelement Windkessel model (see, e.g., [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF]Chapter 10] and [START_REF] Frank | Die Grundform des arteriellen Pulses[END_REF]), where the pressure P l : R + → R and the flux Q l : R + → R on the outlet Γ l are related through the following algebraic-differential equations:

   C d,l dπ l dt + π l R d,l = Q l , P l = R p,l Q l + π l , (2)
for l = 1, . . . , n 0D . Here, R p,l and R d,l model the resistance to the flow of the vasculature proximal and distal to each outlet Γ l , respectively, and the capacity C d,l , take into account the deformability of the downstream vessels. The values P l and π l are also called proximal and distal pressures, respectively.

Remark 2.1 We will mainly focus on a three-element Windkessel model, that is one of the most popular choices in computational hemodynamics. However, this choice is due purely illustrative and the methods and the analysis presented in the following sections to can be extended to more complex lumped parameter models (see Section 3.2.4 for details).

3D-0D coupling conditions

The 3D-0D coupling between (1) and ( 2) is defined through the following relations on each

Γ l    Q l = Γ l u • n f , σ(u, p)n f = -P l n f on Γ l , (3)
for l = 1, . . . , n 0D and where n f denotes the exterior unit-vector normal of Ω f .

Energy balance Let the quantity

E(t) def = ρ f 2 u 2 0,Ω f + n0D l=1 C d,l 2 π 2 l
denote the total (kinetic + potential) energy of the 3D-0D coupled system given by ( 1)-(3), while

D(t) def = 2µ t 0 ǫ(u(s)) 2 0,Ω f ds + n0D l=1 t 0 π 2 l (s) R d,l + R p,l Q 2 l (s) ds > 0 (4)
represents the dissipative effects. Assuming that u in = 0 (free system) and using a standard energy argument, we get the following identity

E(t) + D(t) + t 0 Γ out ρ f 2 |u(s)| 2 u(s) • n f ds = E(0). (5)

RR n°7937

Remark 2.2 Since the last term of the left-hand side of Equation (5) can be negative, this expression does not guarantee a correct energy balance across the 3D-0D interface Γ out . This issue is well-known in computational hemodynamics. The interested reader is referred to [4] for a stabilization technique, and to [START_REF] Formaggia | On the stability of the coupling of 3D and 1D fluid-structure interaction models for blood flow simulations[END_REF] for a different 3D-1D coupling. The methods introduced in this paper can be easily adapted to these alternative formulations.

Time semi-discretization

We consider a fractional-step time-marching of the fluid equations (1) and a backward Euler scheme for the lumped parameter model (2). In what follows, the parameter τ denotes the time-step size, we set t n def = nτ for n ∈ N and ∂ τ x n def = (x nx n-1 )/τ stands for the first-order backward difference.

Fractional-step fluid time-marching

Several variants of the original Chorin-Temam projection scheme [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des equations de Navier-Stokes par la méthode des pas fractionaires I[END_REF] have been proposed in the literature (see, e.g., [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for a recent review). The methods presented and analyzed in section 3 below do not a priori depend on the specific formulation considered for the projection scheme. To fix the ideas and without generality loss, we consider the non-incremental pressurecorrection version (see, e.g., [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]Section 4]). Hence, the time semi-discrete approximation of 1 is performed as follows. We set u 0 = u 0 = u 0 and, for n ≥ 1, we compute (u n , p n , u n ) by solving: 1. Viscous step:

         ρ f u n -u n-1 τ + ρ f u n-1 • ∇ u n -2µ∇ • ǫ( u n ) = 0 in Ω f , u n = u in (t n ) on Γ in , u n = 0 on Σ. (6)
where the 2. Projection step:

             ρ f u n -u n τ + ∇p n = 0 in Ω f , ∇ • u n = 0 in Ω f , u n • n f = u in (t n ) • n f on Γ in , u n = 0 on Σ. (7)
From the implementation point of view, the projection step ( 7) is usually reformulated by as the pressure-Poisson problem

     - τ ρ f ∆p n = -∇ • u n in Ω f , τ ρ f ∂p n ∂n f = 0 on Γ in ∪ Σ, (8)
which requires further regularity on the pressure (H 1 (Ω f ) instead of L 2 (Ω f ) in practice). Then, the divergence free (or end-of-step) velocity u n can be eliminated in (6) using the following relation (from (7) 1 )

u n = u n - τ ρ f ∇p n . (9)
It should be noted that the boundary conditions on Γ out have been omitted deliberately in (6)-( 8), since they depend on the type of 3D-0D coupling scheme considered. In particular, they lead to different formulations of the pressure-Poisson equation [START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF] (see Section 3.1).

Inria

Backward-Euler Windkessel time-marching

Without loss of generality, we consider a backward Euler time-discretization of (2), which yields

   C d,l ∂ τ π n l + π n l R d,l = Q n l , P n l = R p,l Q n l + π n l (10)
or, equivalently,

π n l = α l π n-1 l + β l Q n l , P n l = γ l Q n l + α l π n-1 l , (11)
with the notation

α l def = R d,l C d,l R d,l C d,l + τ , β l def = R d,l (1 -α l ) and γ l def = R p,l + β l .

Spatial discretization

In what follows, we will consider the usual Sobolev space H 1 (Ω), for a given domain Ω ⊂ R 3 . Then, for X ⊂ ∂Ω (with meas(X) > 0), we define H 1 X (Ω) the subspace of H 1 (Ω) with vanishing trace on X. The scalar product in L 2 (Ω) is denoted by (•, •) Ω and its associated norm by • 0,Ω .

We consider a family of triangulations {T f,h } 0<h≤1 of the domain Ω f satisfying the usual requirements of finite element approximations (see, e.g., [7]). The subscript h ∈ (0, 1] refers to the level of refinement of the triangulations. In order to ease the presentation, we assume that the family of triangulations is quasi-uniform. For the discretization in space of (1), we consider continuous Lagrange finite element approximations V h and R h of [H 1 (Ω f )] 3 and H 1 (Ω f ), respectively. Other choices of approximation spaces are possible for the projection method (see [START_REF] Guermond | On stability and convergence of projection methods based on pressure Poisson equation[END_REF] for a discussion). For a given X ⊂ ∂Ω f (with meas(X) > 0), we set

V X,h def = V h ∩ [H 1 X (Ω f )] 3 , R X,h def = R h ∩ H 1 X (Ω f ).
3 Fractional-step time-marching and 3D-0D coupling schemes

In this section, we describe two coupling schemes (explicit and implicit) resulting from appropriate time discretizations of the coupling conditions (3).

Explicit 3D-0D coupling scheme

In this case the 3D-0D coupling conditions (3) are time discretized as follows

         Q n l = Γ l u n • n f , p n = P n l on Γ l , 2µǫ( u n )n f = 0 on Γ l , (12)
for l = 1, . . . , n 0D . Note that the continuity of fluxes (3) 1 is treated explicitly by using the flux of the latest computed viscous velocity. For the relation (3) 2 we consider a Dirichlet boundary condition for the pressure, while the viscous part of the fluid stresses is set to zero. This is a standard procedure to decouple the projection and viscous steps in the framework of projection schemes with natural boundary conditions (see, e.g., [START_REF] Guermond | Error analysis of pressure-correction schemes for the time-dependent Stokes equations with open boundary conditions[END_REF] and [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]Section 10]).

Remark 3.1 It is well-known that the artificial Dirichlet boundary condition (12) 2 could lead to sub-optimal approximations since, in general, the relation (12) 2 is not consistent with the solution of continuous problem. For further details on this issue, we refer to the analysis reported [START_REF] Guermond | Error analysis of pressure-correction schemes for the time-dependent Stokes equations with open boundary conditions[END_REF] which suggests the use of the rotational version of the pressure-correction scheme.

The resulting fully discrete time-marching procedure is reported in Algorithm 1. In the viscous-step [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF] we have considered the standard Temam's consistent term, ρ f 2 ((∇• u n-1 ) u n , v) Ω f , which stabilizes the semi-implicit treatment of the convective term. The 3D-0D explicit coupling given by ( 12) allows a fully uncoupled computation of the Windkessel state, fluid pressure and velocity. This is particularly appealing from the implementation and computational efficiency point of view. Nevertheless, as suggested in Section 3.1.1 below (and then confirmed by numerical experiments in Section 5), Algorithm 1 may suffer from stability issues.

Algorithm 1 (Explicit 3D-0D coupling scheme)

Let u 0 def = u 0 , u 0 ∈ V h and π 0 1 , . . . , π 0 n0D ∈ R be given initial data. For n ≥ 1 perform: 1. Viscous step: Find u n ∈ V Σ,h such that          u n | Γ in = u in (t n ), ρ f τ ( u n , v) Ω f + ρ f ( u n-1 • ∇ u n , v) Ω f + ρ f 2 ((∇ • u n-1 ) u n , v) Ω f + 2µ (ǫ( u n ), ǫ(v)) Ω f = ρ f τ u n-1 , v Ω f (13) for all v ∈ V Σ∪Γ in ,h . Thereafter set Q n l def = Γ l u n • n f . 2. Windkessel step: For l = 1, . . . , n 0D , compute (Q n l , π n l , P n l ) ∈ R 3 from Q n l = Q n l , C d,l ∂ τ π n l + π n l R d,l = Q n l , P n l = R p,l Q n l + π n l . (14)
3. Projection step: Find

p n ∈ R h such that    p n | Γ l = P n l , l = 1, . . . , n 0D , τ ρ f (∇p n , ∇q) Ω f = -(∇ • u n , q) Ω f (15)
for all q ∈ R Γ out ,h . Thereafter set

u n def = u n -τ ρ f ∇p n ∈ [L 2 (Ω f )] d .

Stability analysis

Let the quantities

E n def = ρ f 2 u n 2 0,Ω f + n0D l=1 C d,l 2 |π n-1 l | 2 , D n def = 2µ n m=1 τ ǫ( ũm ) 2 0,Ω f + n-1 m=1 n0D l=1 τ |π m l | 2 R d,l + R p,l |Q m l | 2 ,
for n ≥ 1, denote the energy and physical dissipation of the discrete system. Let us also set

E 0 def = ρ f 2 u 0 2 0,Ω f + n0D l=1 C d,l 2 |π 0 l | 2 ,

Inria

We then have the following energy based result (whose proof can be found in [5,Chapter 2]).

Theorem 3.2 Let ( u n , p n ) n≥1 and (Q n l , π n l , P n l ) 1≤l≤n0D n≥1 be the solution given by Algorithm 1 and assume that u in = 0 (free system). The following inequality holds for n ≥ 1

E n + D n + n m=1 ρ f 2 τ ( u m-1 • n f , | u m | 2 ) Γ out ≤ E 0 - n-1 m=1 τ 2 2ρ f ∇p m 2 0,Ω f + n-1 m=1 τ (∇ • u m , φ m ) Ω f + n-1 m=1 τ 2 ρ f (∇p m , ∇φ m ) Ω f , ( 16
)
where φ n ∈ R h is an arbitrary discrete lifting of the (unknown) proximal pressures, namely,

φ n = P n l on Γ l , l = 1, . . . , n 0D . ( 17
)
The left-hand side of estimate ( 16) corresponds to the discrete counterpart of ( 5). Nevertheless, the artificial power introduced by the last two terms of ( 16),

n-1 m=1 τ (∇ • u m , φ m ) Ω f + n-1 m=1 τ 2 ρ f (∇p m , ∇φ m ) Ω f , (18)
cannot be controlled, so that this estimate does not guarantee the energy stability of the approximations provided by Algorithm 1. Two remarks are now in order.

Remark 3.3 It is worth mentioning that the term (18) corresponds to the residual of the projection step (15) (note that φ n / ∈ R Γ out ,h ). In fact, for the space continuous counterpart of (15) we have

(∇ • u n , φ n ) Ω f + τ ρ f (∇p n , ∇φ n ) Ω f = τ ρ f Γ out ∂p n ∂n f φ n = - n0D l=1 Γ l u n • n f - Γ l u n • n f P n l .
Hence, the uncontrolled artificial power involved in the energy estimate [START_REF] Guermond | On stability and convergence of projection methods based on pressure Poisson equation[END_REF] is due to the time-lag in the flux (

Γ l u n • n f instead of Γ l u n • n f )
introduced by the explicit treatment of the continuity of fluxes (12) 1 on the 3D-0D interfaces Γ l . It should be noted that a 3D-0D explicit coupling with a monolithic time-marching scheme in the fluid yields a similar time lag, but now with

Γ l u n-1 • n instead of Γ l u n • n.
Remark 3.4 In the case of a single outlet (i.e., n 0D = 1), we can take φ n = P n in Ω f . From the proof of Proposition 3.2 (see [5,Chapter 2]) one then recovers the following energy estimate for the fluid

ρ f 2 u n 2 0,Ω f + 2µ n m=1 τ ǫ( ũm ) 2 0,Ω f + n m=1 ρ f 2 τ ( u m-1 • n f , | u m | 2 ) Γ out ≤ ρ f 2 u 0 2 0,Ω f - n-1 m=1 τ 2 2ρ f ∇p m 2 0,Ω f .
As a result, the energy stability of the fluid does not depend on the imposed outlet pressure.

RR n°7937

Implicit 3D-0D coupling scheme

In this case the 3D-0D coupling conditions (3) are time-discretized as follows

         Q n l = Γ l u n • n f , p n = P n l on Γ l , 2µǫ( u n )n f = 0 on Γ l , (19)
for l = 1, . . . , n 0D . Note that, in comparison with [START_REF] Frank | Die Grundform des arteriellen Pulses[END_REF], the above coupling scheme treats implicitly the continuity of fluxes on the outlet boundaries. This feature enhances stability, as we will show in Section 3.2.1. However, from the computational point of view, the relations (19) 1,2 and (10) apparently couple the evaluation of ( 8) and [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]. Fortunately, this difficulty can be circumvented via an appropriate reformulation of the pressure boundary condition [START_REF] Moghadam | A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations[END_REF] 2 for the projection step [START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF]. Indeed, by inserting ( 9) into ( 19) 1 , we get

Q n l = Q n l - τ ρ f Γ l ∂p n ∂n f , (20)
which with ( 11) and ( 19) 2 yields the following (implicit) boundary condition for the outlet pressures:

p n | Γ l = γ l Q n l - γ l τ ρ f Γ l ∂p n ∂n f + α l π n-1 l (21)
for l = 1, . . . , n 0D . Note that this expression still enforces p n to be constant on each Γ l . Multiplying (8) 1 by q ∈ R h , integrating by parts, using (8) 1 and the fact that

q| Γ l is constant, we get τ ρ f (∇p n , ∇q) Ω f - τ ρ f n0D l=1 Γ l ∂p n ∂n f q| Γ l = -(∇ • u n , q) Ω f
for all q ∈ R h . We can eliminate the normal derivative of the pressure using [START_REF] Quarteroni | Coupling between lumped and distributed models for blood flow problems[END_REF], which yields the following modified variational formulation for the projection step: Find

p n ∈ R h such that τ ρ f (∇p n , ∇q) Ω f + n0D l=1 (p n | Γ l )(q| Γ l ) γ l = n0D l=1 Q n l + α l π n-1 l γ l q| Γ l -(∇ • u n , q) Ω f (22)
for all q ∈ R h . We can then set P n l = p n | Γ l and retreive (Q n l , π n l ) from [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF], for l = 1, . . . n 0D .

Remark 3.5 The well-posedness of the pressure-Poisson problem [START_REF] Shi | Review of zero-d and 1-d models of blood flow in the cardiovascular system[END_REF] follows from a generalized Poincaré's inequality, which guarantees the coercivity of the left-hand side of (22) in R h .

Remark 3.6 Testing (22) with q = 1, and since

P n l = p n | Γ l , we have n0D l=1 P n l -α l π n-1 l γ l = n0D l=1 Q n l + Ω f ∇ • u n .
Hence, integrating by parts in the last term, using (11) 2 and owing to (6) 2,3 we get the following mass conservation for the Windkessel fluxes:

n0D l=1 Q n l = - Γ in u in (t n ).

Inria

The complete time-marching procedures is reported in Algorithm 2.

Algorithm 2 (3D-0D implicit coupling scheme)

Let u 0 = u 0 , u 0 ∈ V h , and π 0 1 , . . . , π 0 n0D ∈ R be given initial data. For n ≥ 1 perform: 1. Viscous step: Find

u n ∈ V Σ,h such that          u n | Γ in = u in (t n ), ρ f τ ( u n , v) Ω f + ρ f ( u n-1 • ∇ u n , v) Ω f + 2µ (ǫ( u n ), ǫ(v)) Ω f + = ρ f τ u n-1 , v Ω f (23) for all v ∈ V Σ∪Γ in ,h . Thereafter set Q n l def = Γ l u n • n f . 2. Projection-Windkessel step: Find p n ∈ R h and such that τ ρ f (∇p n , ∇q) Ω f + n0D l=1 (p n | Γ l )(q| Γ l ) γ l = n0D l=1 Q n l + α l π n-1 l γ l q| Γ l -(∇ • u n , q) Ω f (24)
for all q ∈ R h . Thereafter, set

P n l = p n | Γ l and compute (Q n l , π n l ) ∈ R 2 from the relations Q n l = P n l -α l π n-1 l γ l , π n l = α l π n-1 l + β l Q n l , l = 1, . . . n 0D (25)
and set

u n def = u n - τ ρ f ∇p n ∈ [L 2 (Ω f )] d .

Stability analysis

The focus of this section is to present the stability result of the formulation ( 23)-(25) summarized in the following proposition. Theorem 3.7 Let ( u n , p n ) n≥1 and (Q n l , π n l , P n l ) 1≤l≤n0D n≥1 be the approximated solution given by Algorithm 2 and assume that u in = 0 (free system). The following energy inequality holds

E n + D n + n m=1 ρ f 2 ( u m-1 • n f , | u m | 2 ) Γ out ≤ E 0 - n-1 m=1 τ 2 2ρ f ∇p m 2 0,Ω f . (26)
Proof. We first test [START_REF] Temam | Sur l'approximation de la solution des equations de Navier-Stokes par la méthode des pas fractionaires I[END_REF] with v = u n and integrate by parts the convective term. This yields the identity

ρ f 2 ∂ τ u n 2 0,Ω f + ρ f 2τ u n -u n-1 2 0,Ω f + 2µ ǫ( u n ) 2 0,Ω f + ∇p n-1 , u n Ω f + ρ f 2 ( u n-1 • n f , | u n | 2 ) Γ out = 0 (27)
for n ≥ 2 and, for n = 1, we get

ρ f 2τ u 1 2 0,Ω f -u 0 2 0,Ω f + ρ f 2τ u 1 -u 0 2 0,Ω f + 2µ ǫ( u 1 ) 2 0,Ω f + ρ f 2 ( u 0 • n f , | u 1 | 2 ) Γ out = 0. ( 28
) RR n°7937
Thereafter, taking (24) at time step n -1, testing with q = p n-1 and integrating by parts in its right-hand side it yields

τ ρ f ∇p n-1 2 0,Ω f + n0D l=1 P n-1 l -α l π n-2 l γ l P n-1 l = u n-1 , ∇p n-1 Ω f .
for n ≥ 2. Hence, the addition and subtraction of suitable terms and the application of the Cauchy-Schwarz and arithmetic-geometric inequalities yields

τ 2ρ f ∇p n-1 2 0,Ω f - ρ f 2τ u n -u n-1 2 0,Ω f -u n , ∇p n-1 Ω f + n0D l=1 P n-1 l -α l π n-2 l γ l P n-1 l ≤ 0. ( 29
)
As a result, the summation of ( 27) and (29) gives

ρ f 2 ∂ τ u n 2 0,Ω f + 2µ ǫ( u n ) 2 0,Ω f + ρ f 2 ( u n-1 • n f , | u n | 2 ) Γ out + n0D l=1 P n-1 l -α l π n-2 l γ l P n-1 l ≤ - τ 2ρ f ∇p n-1 2 0,Ω f (30)
for n ≥ 2. At last, from (25) and its equivalence to (11), we have

P n l -α l π n-1 l γ l P n l = Q n l P n l = R p,l |Q n l | 2 + Q n l π n l ≥ C d,l 2 ∂ τ |π n l | 2 + 1 R d,l |π n l | 2 + R p,l |Q n l | 2 (31)
for n ≥ 1. Hence, by inserting the last inequality of (31) into (30), multiplying by τ and summing over m = 2, ..., n we get the estimate

E n + D n + n m=2 ρ f 2 ( u m-1 • n f , | u m | 2 ) Γ out ≤ E 1 - n-1 m=1 τ 2 2ρ f ∇p m 2 0,Ω f .
for n ≥ 2. We recover the estimate (26) by simply adding to this inequality the expression (28) multiplied by τ , which completes the proof. The estimate (26) corresponds to the discrete counterpart of (5). Note that the right-hand side of (26) is a pure numerical dissipation term (the natural pressure stabilization of the projection scheme). Therefore, the 3D-0D coupling reported in Algorithm 2 does not introduce any uncontrolled artificial power and, hence, a guaranty of numerical stability. This feature will be illustrated in Section 5 via numerical experiments.

Implementation details

In this section we discuss the implementation of the pressure problem [START_REF] Vignon-Clementel | Outflow boundary conditions for 3d simulations of non-periodic blood flow and pressure fields in deformable arteries[END_REF] in a finite element framework. For the sake of simplicity, and without loss of generality, we limit the discussion to the case of a single outlet. We define the arrays P, V ∈ R N corresponding to the degrees of freedom (d.o.f.) of the pressure p n ∈ Q h and of a general test function q ∈ Q h , respectively. The bilinear form (∇p n , ∇q) Ω f , without imposing Dirichlet boundary conditions to p n , can be written in matrix form as

(∇p n , ∇q) Ω f = V T A P = V T I V T O A II A IO A OI A OO P I P O ,

Inria

where the subindexes O and I indicate the elements of the array corresponding to the d.o.f. on Γ out and Ω\Γ out , respectively. Hence, the pressure-Poisson projection step with explicit Dirichlet data can be formulated as

A II P I = F I -A IO 1 O p n |Γ out , (32)
where the notation in right-hand side is such that

V T I V T O F I F O = ρ f τ (∇ • u n , q) Ω f , with 1 O ∈ R N O
denoting a vector of ones, with N O being the number of pressure d.o.f. on Γ out .

The linear system (32) is usually solved by means of a preconditioned conjugate gradient methods (PCG), with the preconditioning operator A -1 II given, for example, by an incomplete Cholesky factorization of A II (alternative preconditioners could be used). With the notations introduced above, the stiffness matrix of the implicit formulation ( 24) can be derived straightforwardly.

Indeed, since V O = 1 O q |Γ out and P O = 1 O p n |Γ out , we obtain A II a a T b P I p n |Γ out =   F I 1 T O F I + ρ f τ Q n + απ n-1 γ   , (33)
with

a def = A IO 1 O and b def = 1 T O A OO 1 O + ρ f τ γ
. It should be noted that the matrix in (33) has a non-standard sparsity pattern. However, since we use a Krylov linear solver this matrix is never assembled in practice (only matrix-vector products are evaluated). In the numerical experiments, we have considered the block-preconditioner given by

A -1 II 0 0 1 b
, which yielded practically the same number of PCG-iterations than in the solution of (32). In the general case of a domain Ω f with multiple outlets, the aforementioned considerations can be extended by considering one additional equation for each 3D-0D interface. Concerning the computational cost, in our numerical simulations we did observe any relevant difference between the implicit and the explicit coupling.

Extension to higher order time-splitting schemes

Although widely used in practice, the original Chorin-Temam projection scheme might suffer of a limited accuracy in time and of spurious boundary layers, due to the unphysical homogenous Neumann boundary condition [START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF] 3 . Among the several variants available (see, e.g., [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for an overview), in this section we describe a possible extension of the implicit 3D-0D coupling (Algorithm 2) in the context of an incremental pressure projection scheme with a second order time discretization.

Following [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], we decompose now the time iteration in a BDF2 time discretization for the viscous step

             ρ f 3 u n -4u n-1 + u n-2 2τ + ρ f u n-1 • ∇ u n -2µ∇ • ǫ( u n ) + ∇p n-1 = 0 in Ω f , u n = u in (t n ) on Γ in , u n = 0 on Σ, 2µǫ( u n )n f = 0 on Γ out , RR n°7937
a projection step for the increment of pressure

           - τ ρ f ∆δp n = - 3 2 ∇ • u n in Ω f , τ ρ f ∂δp n ∂n f = 0 on Γ in ∪ Σ, δp n = P n l -P n-1 l on Γ l , l = 1, . . . , n 0D (34)
completed by the end-of-step updates: p n = p n-1 + δp n and u n = u n -2τ 3ρ f ∇δp n . In order to have a second-order time-discretization for the whole 3D-0D problem, we also discretize the 0D model with a BDF2 scheme, namely

π n l = α l π n-1 l - α l 4 π n-2 l + β l Q n l , P n l = γ l Q n l + α l π n-1 l - α l 4 π n-2 l , (35)
with

α l def = 2R d,l C d,l (3/2)R d,l C d,l + τ , β l def = R d,l τ (3/2)R d,l C d,l +τ , γ l def = R p,l + β l .
Hence, using equation (35) 2 to define the implicit coupling with the projection step, we obtain the following time-stepping method: 1. Viscous step: Find

u n ∈ V Σ,h such that          u n | Γ in = u in (t n ), ρ f 2τ (3 u n , v) Ω f + ρ f ( u n-1 • ∇ u n , v) Ω f + 2µ (ǫ( u n ), ǫ(v)) Ω f + ρ f 2 ((∇ • u n-1 ) u n , v) Ω f = -(∇p n , v) + ρ f 2τ 4u n-1 -u n-2 , v Ω f for all v ∈ V Σ∪Γ in ,h . Thereafter set Q n l def = Γ l u n • n f . 2. Increment-Windkessel step: Find δp n ∈ R h and such that τ ρ f (∇δp n , ∇q) Ω f + n0D l=1 (δp n | Γ l )(q| Γ l ) γ l = n0D l=1 Q n l + α l π n-1 l -ψ l π n-2 l -P n-1 l γ l q| Γ l -(∇ • u n , q) Ω f for all q ∈ R h . 3. End-of-step: Set p n def = p n-1 + δp n , u n def = u n - 2τ 3ρ f ∇δp n ∈ [L 2 (Ω f )] d , P n l = p n | Γ l and compute (Q n l , π n l ) ∈ R 2 from the relations (25).
Remark 3.8 It is well known that the stability of the pressure-incremental projection scheme requires that the finite element spaces for velocity and pressure satisfy an inf-sup condition (see [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for an extended discussion).

Extension to more complex lumped parameter models

The algorithm and analysis presented above can be straightforwardly extended to more complex networks of lumped parameter models, made of resistances, capacitances, and inductances. These networks can be obtained, for instance, by connecting several Windkessel elements (see, e.g., [START_REF] Korakianitis | A concentrated parameter model for the human cardiovascular system including heart valve dynamics and atrioventricular interaction[END_REF][START_REF] Moghadam | A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations[END_REF][START_REF] Shi | Review of zero-d and 1-d models of blood flow in the cardiovascular system[END_REF] and the references therein).

Inria

Let us consider the following general lumped parameter model

     C dP dt + RP = Q + HΨ, L dΨ dt = -H ⊺ P, (36)
with

P ⊺ def = P 1 . . . P n0D π 1 . . . π nπ ,
the collection of outlet P l , l = 1, . . . , n 0D and distal pressures π j , j = 1, . . . , n π , and Ψ are the fluxes through the inductances, both representing the dynamical state of the lumped parameter model. Moreover, we assume that Q(t) ∈ R np has the form

Q ⊺ def = Q 1 . . . Q n0D 0 . . . 0 ,
with Q l the input flux at the outlet Γ l . In (36), C, R, L denote symmetric capacitances, resistances and inductances matrices. This ensures the correct energy balance of the system in the 3D-0D time-space continuos formulation, namely d dt

1 2 P ⊺ CP + 1 2 Ψ ⊺ LΨ = -P ⊺ RP + P ⊺ Q,
since that the last term of the right-hand-side cancels out when coupling (36) with (1). As an example, the three-element Windkessel presented above, with an additional inductance L p parallel to R p , can be rewritten in this format obtaining the following expressions for the system matrices

C = 0 0 0 C d , R = 1/R p -1/R p -1/R p 1/R p + 1/R d , H = -1 1 , L = L p .
More complex models can be represented in this format analagously. Discretizing (36) in time (e.g., using a backward Euler scheme), we obtain

P n = AP n-1 + BQ n + B HΨ n-1 , (37)
with

B def = (C/τ + R) -1 , A def = B(C/τ + τ HL -1 H ⊺ ).
In order to derive the 3D-0D implicit coupling scheme (Algorithm 2), we need explicit expressions for Q n l , l = 1, . . . , n 0D . These can be computed algebraically from (37) by first isolating the degrees-of-freedom on the outlets via

B    Q n 1 . . . Q n n0D    =    P n 1 . . . P n n0D    -ĀP n-1 -B HΨ n-1 ,
with () denoting the first n 0D rows and () the first n 0D rows and columns, respectively. Hence, we obtain the following expression

   Q n 1 . . . Q n n0D    = G    P n 1 . . . P n n0D    -G ĀP n-1 -G B HΨ n-1 ,
with G = ( B) -1 . Finally, by combining the latter with [START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF], we obtain the following formulation of the pressure projection-Windkessel step (compare with ( 24)): Find

p n ∈ R h such that τ ρ f (∇p n , ∇q) Ω f + n0D k,l=1 g k,l (p n | Γ k )(q| Γ l ) = n0D l=1 ( Q n l + b n-1 l q| Γ l )(q| Γ l ) -(∇ • u n , q) Ω f (38)
for all q ∈ R h and with g k,l

def = [G] k,l , b n l def = [G ĀP n + G B HΨ n ] l .
Note that, in this general formulation, the pressure at the outlets can be coupled through the lumped parameter model since generally g k,l = 0 for k = l. Remark 3.9 Using the same arguments than in Section 3.2.1, the unconditional stability results of Theorem 3.7 can be extended to a fractional-step scheme with the generalized pressure projection formulation (38).

Incompressible fluid-structure interaction

Fractional-step time-marching schemes have been a valuable tool for the design of efficient solution methods for incompressible fluid-structure interaction (FSI) problems, yielding the projection semi-implicit coupling schemes [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: Stability analysis and numerics[END_REF][START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]. This coupling approach is based on the following three basic ideas:

• treat explicitly the geometrical non-linearities and the viscous-structure coupling, which reduces computational complexity;

• treat implicitly the presssure-structure coupling, which avoids numerical instability;

• perform this explicit-implicit splitting through a projection scheme in the fluid.

So far, the stability of this method has been analyzed within a simplified framework which enforces null pressure on the outlet boundaries (see [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: Stability analysis and numerics[END_REF][START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]). In this section the analysis is extended to the case of a lumped parameter modeling of the outlet boundaries, with a pressure-Poisson formulation of the projection step, based on the methods of Section 3.

Model problem

For the sake of the analysis (see [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]), we consider as model problem a coupled linear system in which the fluid is described by the Stokes equations, in the fixed domain Ω f , and the structure either by the classical linear elastodynamics equations or by equations based on linear thin-solid models (e.g., plate, shell, etc.). The reference domain of the solid is denoted by Ω s . It will be either a domain or a 2-manifold of R 3 (in this later case the elastic domain is identified to its mid-surface). We denote by Σ def = ∂Ω s ∩ ∂Ω f the fluid-structure interface. In the case the structure is described by thin-solid model we have Ω s = Σ (see Figure 2). The resulting coupled system, describing the fluid velocity u :

Ω f × R + → R 3 , fluid pressure p : Ω f × R + → R and solid Inria Ω f Σ Γ in Γ 1 Γ 2 Ω s Ω s Ω s Figure 2: Examples of geometric configuration Ω s = Σ (Ω s = Σ, see Figure 1). displacement y : Ω s × R + → R 3 , is given by          ρ f ∂ t u -∇ • σ(u, p) = 0 in Ω f , ∇ • u = 0 in Ω f , u = u in on Γ in , u = ẏ in Σ, (39)
ρ s (∂ t ẏ, v s ) Ωs + a s (y, v s ) = -(σ(u, p)n f , v s ) Σ ∀v s ∈ W , ẏ = ∂ t y in Ω s , (40)
completed with the lumped-parameter modeling ( 2)-( 3) on the outlet boundary Γ out (see Figure 2). Here, ρ s denotes the solid density, the abstract bilinear form a s : W × W → R describes the elastic behavior of the structure and W stands for its space of admissible displacements.

It should be noted that the solid problem (40) has been written in weak form, which allows to treat in the same formulation the case of thin and thick solid models. The relations (39) 4 and (40) 1 enforce the kinematic and kinetic interface coupling conditions, respectively. Note that the latter represents also the variational formulation of the structure subproblem. Though simplified, problem (39)-(40) features some of the main numerical issues that appear in complex nonlinear fluid-structure interaction problems involving an incompressible fluid (see, e.g., [START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF]).

Let the quantity E(t)

def = ρ f 2 u 2 0,Ω f + ρs 2 ẏ 2 0,Ωs + 1 2 y 2 s + n0D l=1 C d,l
2 π 2 l denote the total (kinetic + potential) energy of the FSI-0D coupled system given by ( 39)-( 40) and ( 2)-( 3). Here,

• s stands for the elastic energy norm of the solid, defined as y 2 s def = a s (y, y). Assuming that u in = 0 (free system) and using a standard energy argument, we get the following identity

E(t) + D(t) = E(0), (41)
with the dissipative term D(t) > 0 given as in (4).

Note that in this model problem we do not consider the convective term in the fluid, hence we have to solve the fluid in a fixed domain in order to get the right energy balance across the fluid-structure interface Σ.

Spatial discretization

Let {T s,h } 0<h≤1 be a quasi-uniform family of triangulations of the domain Ω s . In order to ease the presentation, we assume that the fluid and solid triangulations T f,h and T s,h match at the interface Σ. For the discretization in space of the solid problem (40), we consider continuous Lagrange finite element approximations, W h of W , which match the fluid velocity discretizations at the interface. Hence,

{v| Σ | v ∈ V h } = {v s | Σ | v s ∈ W h } .
At last, we introduce the standard fluid-sided discrete lifting operator

L h : W h → V h , such that, the nodal values of L h v s vanish out of Σ and (L h v s )| Σ = v s | Σ , for all v s ∈ W h .
In what follows we shall make use of the following continuity estimates (from [9, Lemma 1]) for the discrete lifting operator L h :

L h v s 0,Ω f ≤ C L h 1-α 2 v s 0,Ωs , ∇L h v s 0,Ω f ≤ C L h -1+α 2 v s 0,Ωs (42)
for all v s ∈ W h and with the notation

α def = 0 if Ω s = Σ, 1 if Ω s = Σ. (43) 4.1

.2 Semi-implicit FSI scheme with implicit 3D-0D coupling

In this case we consider numerical approximations of the coupled problem FSI-0D system (39)-( 40) and ( 2)-( 3) by combining the projection semi-implicit coupling scheme reported in [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] with the 3D-0D implicit coupling of Section 3.2. The resulting fully discrete time-marching procedure is reported in Algorithm 3. In the solid subproblem (47), the fluid residual terms are given by

R µ ( u n ), v def = ρ f τ ( u n , v) Ω f + 2µ (ǫ( u n ), ǫ(v)) Ω f - ρ f τ u n-1 , v Ω f , R p (p n , u n ), v def = ρ f τ (u n , v) Ω f - ρ f τ ( u n , v) Ω f -(p n , ∇ • v) Ω f .
Algorithm 3 (Semi-implicit FSI scheme with implicit 3D-0D coupling)

Let u 0 = u 0 , π 0 1 , . . . , π 0 n0D ∈ R and y 0 , ẏ0 ∈ W h be given initial data. For n ≥ 1 perform:

1. Viscous step: Find u n ∈ V h such that        u n | Σ = ẏn-1 | Σ , u n | Γ in = u in (t n ), ρ f τ ( u n , v) Ω f + 2µ (ǫ( u n ), ǫ(v)) Ω f = ρ f τ u n-1 , v Ω f (44) for all v ∈ V Σ∪Γ in ,h . Thereafter set Q n l def = Γ l u n • n f . 2. Implicit step (projection-Windkessel-solid step): Find (p n , u n , y n ) ∈ R h × V h × W h with Inria ẏn = ∂ τ y n and such that      τ ρ f (∇p n , ∇q) Ω f + n0D l=1 (p n | Γ l )(q| Γ l ) γ l = n0D l=1 Q n l + α l π n-1 l γ l q| Γ l -(∇ • u n , q) Ω f -(( ẏn -u n ) • n f , q) Σ , (45)
       u n | Γ in = u in (t n ), u n | Σ = ẏn | Σ , ρ f τ (u n , v) Ω f = ρ f τ ( u n , v) Ω f -(∇p n , v) Ω f , (46)
ρ s (∂ τ ẏn , v s ) Ωs + a s (y n , v s ) = -R µ ( u n ), L h v s -R p (u n , p n ), L h v s (47) for all (q, v, v s ) ∈ R h × V Σ∪Γ out ,h × W h .
Thereafter, set

P n l = p n | Γ l and compute (Q n l , π n l ) ∈ R 2 from the relations Q n l = P n l -α l π n-1 l γ l , π n l = α l π n-1 l + β l Q n l , l = 1, . . . n 0D . (48)
Let the quantities

E n def = ρ f 2 u n 2 0,Ω f + ρ s 2 ẏn 2 0,Ωs + 1 2 y n 2 s + n0D l=1 C d,l 2 |π n l | 2 , D n def = 2µ n m=1 τ ǫ( u m ) 2 0,Ω f + n-1 m=1 n0D l=1 τ |π m l | 2 R d,l + R p,l |Q m l | 2
for n ≥ 1, denote the energy and physical dissipation of the discrete FSI-0D system. We then have the following energy based result.

Theorem 4.1 Let ( u n , p n , y n , ẏn ) n≥1 and (Q n l , π n l , P n l ) 1≤l≤n0D n≥1 be the approximated solution given by Algorithm 3, and assume that u in = 0 (free system) and that the following condition holds

ρ s ≥ 3C L ρ f h 1-α + µτ h 1+α , (49)
with α given by (43). Then, following energy inequality holds

E n + D n E 0 - n-1 m=1 τ 2 ρ f Π ⊥ h (∇p n ) 2 0,Ω f , (50)
with

Π h : L 2 (Ω f ) → V Σ∪Γ in ,h stands for the L 2 -projection operator into V Σ∪Γ in ,h , and Π ⊥ h def = I -Π h for the corresponding orthogonal projection.
Proof. First, testing the viscous step (44

) with v = u n -L h ẏn-1 yields ρ f 2τ u n 2 0,Ω f -u n-1 2 0,Ω f + ρ f 2τ u n -u n-1 2 0,Ω f + 2µ ǫ( u n ) 2 0,Ω f -R µ ( u n ), L h ẏn-1 = 0. ( 51
)

RR n°7937

On the other hand, testing (46

) with v = u n -L h ẏn yields ρ f 2τ u n 2 0,Ω f -u n 2 0,Ω f + ρ f 2τ u n -u n 2 0,Ω f + (∇p n , u n ) Ω f -(p n , ẏn • n f ) Σ -R p (p n , u n ), L h ẏn = 0 (52)
and taking v s = ẏn in (47) yields

ρ s 2 ∂ τ ẏn 2 0,Ωs + ρ s 2τ ẏn -ẏn-1 2 0,Ωs + 1 2 ∂ τ y n 2 s + 1 2τ y n -y n-1 2 s = -R µ ( u n ), L h ẏn -R p (p n , u n ), L h ẏn . ( 53
)
As a result, by adding the equalities (51)-( 53) we get

ρ f 2 ∂ τ u n 2 0,Ω f + 2µ ǫ( u n ) 2 0,Ω f + ρ s 2 ∂ τ ẏn 2 0,Ωs + 1 2 ∂ τ y n 2 s + ρ s 2τ ẏn -ẏn-1 2 0,Ωs + (∇p n , u n ) Ω f -(p n , ẏn • n f ) Σ T1 -R µ ( u n ), L( ẏn -ẏn-1 ) T2 ≤ 0. ( 54
)
Following the argument used in [1, Appendix A], from (45) we infer that

u n = u n + Π ⊥ h ( u n -L h ẏn ) + τ ρ f Π h (∇p n ) . (55)
Thereafter, taking q = p n in (45), integrating by parts in its right-hand side, and since

P n l = p n | Γ l , we have τ ρ f ∇p n 2 0,Ω f + n0D l=1 P n l -α l π n-1 l γ l P n l -( u n , ∇p n ) Ω f + (p n , ẏn • n f ) Σ = 0.
Now, by inserting (55) into this expression and using (31), we get

T 1 = τ ρ f Π ⊥ h (∇p n ) 2 0,Ω f + n0D l=1 P n l -α l π n-1 l γ l P n l -Π ⊥ h ( u n -L h ẏn ) , ∇p n Ω f ≥ τ ρ f Π ⊥ h (∇p n ) 2 0,Ω f + C d,l 2 ∂ τ |π n l | 2 + 1 R d,l |π n l | 2 + R p,l |Q n l | 2 -Π ⊥ h ( u n -L h ẏn ) , ∇p n Ω f T3 .
Therefore, by applying this lower bound to (54) we get

ρ f 2 ∂ τ u n 2 0,Ω f + ρ f 2τ u n -u n-1 2 0,Ω f + 2µ ǫ( u n ) 2 0,Ω f + ρ s 2 ∂ τ ẏn 2 0,Ωs + 1 2 ∂ τ y n 2 s + ρ s 2τ ẏn -ẏn-1 2 0,Ωs + τ ρ f Π ⊥ h (∇p n ) 2 0,Ω f + C d,l 2 ∂ τ |π n l | 2 + 1 R d,l |π n l | 2 + R p,l |Q n l | 2 ≤ T 2 + T 3 . ( 56
)

Inria

Term T 2 can be bounded as in [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF], using (42), which yields

T 2 ≤ ρ f τ u n -u n-1 0,Ω f L h ( ẏn -ẏn-1 ) 0,Ω f + 2µ ǫ( u n ) 0,Ω f ǫ(L h ( ẏn -ẏn-1 )) 0,Ω f ≤ε 1 ρ f 2τ u n -u n-1 2 0,Ω f + ε 2 µ ǫ( u n ) 2 0,Ω f + C L ρ f 2τ ε 1 h 1-α + µ ε 2 h -1-α ẏn -ẏn-1 2 0,Ωs . (57)
Term T 3 can be bounded following the argument used in [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: Stability analysis and numerics[END_REF], which yields

T 3 = Π ⊥ h L h ( ẏn-1 -ẏn ) , Π ⊥ h (∇p n ) Ω f ≤ε 3 τ 2ρ f Π ⊥ h (∇p n ) 2 0,Ω f + ρ f 2τ ε 3 C L h 1-α ẏn -ẏn-1 2 0,Ωs . (58)
Hence, by inserting (57)-( 58) into (56) we get the energy estimate

ρ f 2 ∂ τ u n 2 0,Ω f + ρ f 2τ (1 -ε 1 ) u n -u n-1 2 0,Ω f + µ(2 -ε 2 ) ǫ( u n ) 2 0,Ω f + ρ s 2 ∂ τ ẏn 2 0,Ωs + 1 2 ∂ τ y n 2 s + ρ s 2τ -C L ρ f 2τ h 1-α 1 ε 1 + 1 ε 3 -C L µ ε 2 h -1-α ẏn -ẏn-1 2 0,Ωs + τ ρ f 1 - ε 3 2 Π ⊥ h (∇p n ) 2 0,Ω f + n0D l=1 C d,l 2 ∂ τ |π n l | 2 + 1 R d,l |π n l | 2 + R p,l |Q n l | 2 ≤ 0.
At last, the energy estimate (50) follows by taking in the latter ε 1 = 1 2 , ε 2 = ε 3 = 1, summing over n and using (49), which completes the proof. Proposition 4.1 guarantees the conditional stability of Algorithm 3. Note that the stability condition is similar to the one obtained in [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] with a Darcy-like formulation of the projection step. The estimate (50) corresponds to the discrete counterpart of (41), the right-hand side of (50) is a dissipative numerical term related to the natural pressure stabilization of the projection scheme. Moreover, extensive numerical simulations with the semi-implicit FSI coupling scheme in physiological regimes (see [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF][START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF][START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF]) suggest that condition (49) is not necessary for stability.

Semi-implicit FSI scheme with explicit 3D-0D coupling

We now consider numerical approximations of the coupled problem FSI-0D system (39)-( 40) and ( 2)-( 3) by combining the projection semi-implicit coupling scheme reported in [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] with the 3D-0D explicit coupling of Section 3.1. The resulting fully discrete time-marching procedure is reported in Algorithm 4.

Algorithm 4 (Semi-implicit FSI scheme with explicit 3D-0D coupling) Let u 0 = u 0 , π 0 1 , . . . , π 0 n0D ∈ R and y 0 , ẏ0 ∈ W h be given initial data. For n ≥ 1 perform: 1. Viscous step: Find

u n ∈ V h such that        u n | Σ = ẏn-1 | Σ , u n | Γ in = u in (t n ), ρ f τ ( u n , v) Ω f + 2µ (ǫ( u n ), ǫ(v)) Ω f = ρ f τ u n-1 , v Ω f (59) for all v ∈ V Σ∪Γ in ,h . Thereafter set Q n l def = Γ l u n • n f . RR n°7937 2. Windkessel step: For l = 1, . . . , n 0D , compute (Q n l , π n l , P n l ) ∈ R 3 from          Q n l = Q n l , C d,l ∂ τ π n l + π n l R d,l = Q n l , P n l = R p,l Q n l + π n l . (60)
3. Implicit projection-solid step:

Find (p n , u n , y n ) ∈ R h × V h × W h with ẏn = ∂ τ y n and such that    p n | Γ l = P n l , l = 1, . . . , n 0D , τ ρ f (∇p n , ∇q) Ω f = -(∇ • u n , q) Ω f -(( ẏn -u n ) • n f , q) Σ (61)        u n | Γ in = u in (t n ), u n | Σ = ẏn | Σ , ρ f τ (u n , v) Ω f = ρ f τ ( u n , v) Ω f -(∇p n , v) Ω f , (62)
ρ s (∂ τ ẏn , v s ) Ωs + a s (y n , v s ) = -R µ ( u n ), L h v s -R p (u n , p n ), L h v s (63) for all (q, v, v s ) ∈ R Γ out ,h × V Σ∪Γ in ,h × W h .
The following proposition provides an energy estimate for the approximations provided by Algorithm 4. Theorem 4.2 Let ( u n , p n , y n , ẏn ) n≥1 and (Q n l , π n l , P n l ) 1≤l≤n0D n≥1 be the approximated solution given by Algorithm 4 and assume that u in = 0 (free system). Then, under the condition (49), the following energy inequality holds

E n + D n E 0 - n-1 m=1 τ 2 2ρ f Π ⊥ h (∇p m ) 2 0,Ω f + n-1 m=1 τ (∇ • u m , φ m ) Ω f + n-1 m=1 τ 2 ρ f (∇p m , ∇φ m ) Ω f + n-1 m=1 τ ( ẏm • n f , φ m ) Σ . (64)
Proof. The result follows by combining the arguments involved in the proofs of Propositions 3.2 and 4.1. Hence, only partial details are given. The main difference lies on the estimation of term T 1 in (54). Since q ∈ R Γ out ,h , integration by parts in (61) gives

(∇p n , ∇q) Ω f = ρ f τ ( u n , ∇q) Ω f -( ẏn • n f , q) Σ , so that by testing with q = (p n -φ n ) ∈ R Γ out ,h we get τ ρ f ∇p n 2 0,Ω f -( u n , ∇p n ) Ω f +( ẏn • n f , p n ) Σ + u n - τ ρ f ∇p n , ∇φ n Ω f -( ẏn • n f , φ n ) Σ = 0. ( 65
)

Inria

As in the proof of Proposition 4.1, from (62) we get (55). Hence, inserting this expression into (65) we get

-(∇p n , u n ) Ω f + (p n , ẏn • n f ) Σ + τ ρ f Π ⊥ h (∇p n ) 2 0,Ω f -Π ⊥ h ( u n -L h ẏn ) , ∇p n Ω f = -u n - τ ρ f ∇p n , ∇φ n Ω f + ( ẏn • n f , φ n ) Σ - n0D l=1 Q n l P n l + (∇ • u n , φ n ) Ω f + τ ρ f (∇p n , ∇φ n ) Ω f + ( ẏn • n f , φ n ) Σ (66)
and the estimate (64) follows using the same arguments than in the proof of Proposition 4.1.

Remark 4.3 A comparison of the energy estimates ( 16) and (64) suggests that the fluid-solid interaction introduces and additional destabilizing effect in the explicit splitting of the 3D-0D coupling [START_REF] Frank | Die Grundform des arteriellen Pulses[END_REF], due to the presence of the artificial interface term

n-1 m=1 τ ( ẏm • n f , φ m ) Σ .
In particular, it is worth noting that the observation made in Remark 3.4 for the case of a single outlet is not valid in the FSI framework, since the above term does not vanish for φ m = P m in Ω f . This point will be illustrated trough numerical experiments in Section 5.2.

Numerical experiments

In this section we present two numerical experiments that confirm the analysis of the previous sections.

Blood flow in a patient-specific aorta

Our first example is a pure Navier-Stokes flow within a patient-specific aorta with repaired coarctation (see Figures 3-Left-Center). The geometry comes from the euHeart database segment growing registration algorithm (see [2,3]) was used for the segmentation of the geometry from the medical image. The resulting surface was pre-processed with 3-matic (Materialise) and the final mesh was generated with Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The inflow curve used as boundary condition (Figure 4-Right, black-line) was obtained from the same patient with Phase Contrast MRI. The initial constant pressure was set to 47 mmHg, and the Windkessel parameters (reported in Table 1) where calibrated in order to approximate the measured pressure at the coarctation and the measured flow on each outlet. The physical parameters of the fluid are µ = 0.035 Po and ρ f = 1 gr/cm 3 . For the numerical simulation, we use P 1 finite elements for both pressure and velocity fields (SUPG stabilization in the viscous-convective step, where we assume that the first perform the Chorin-Temam splitting and then we stabilize) and a time step τ = 10 -3 s. A snapshot of the simulation results is shown in Figure 3-Right. Figures 4 and5 show the pressure and flow results for both schemes (Algorithms 1 and 2, respectively). The spurious oscillations in the approximation provided by the explicit 3D-0D coupling scheme are clearly visible for all the outlet pressures and (less pronounced) for some outlet fluxes, while the implicit formulation guarantees stability within the whole cardiac cycle. Hence, these results are in agreement with the stability estimates provided by Theorems 3.2 and 3.7. As further validation of the numerical scheme, Figure 6 shows the pressures and flows curves computed using a monolithic solver, implicitely coupled with the 0D-model, showing a satisfactory agreement of the results. for inlet (black) and outlets 1-4 (blue, green, red, cyan) obtained with a monolithic method (circles) and with the projection method (continuos) using an implicit 3D-0D coupling.

Fluid-structure interaction in an idealized AAA

We consider an idealized abdominal aortic aneurysm (AAA) of length 22.95 cm, minimal diameter 1.7 cm (tubular part), maximal diameter 4.98 cm (aneurysm) and wall thickness 0.2 cm (see Figure 7 (left)). The structure is described by the linear elastodynamics equations and the fluid by the incompressible Navier-Stokes equations in a moving domain (ALE formalism). We considered the non-linear counterpart of Algorithms 3 and 4 using the projection semi-implicit coupling scheme reported in [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]. The solid has Young's modulus 1 M P a, Poisson ratio is 0.46 and density 1.2 g/cm 3 . The fluid viscosity is µ = 0.035 P o and its density 1 g/cm 3 . The inlet and the outlet parts of the solid are clamped. In the fluid, a parabolic velocity profile is enforced at the inlet, with a realistic inflow (peak velocity ≈ 96 cm/s). The Windkessel parameters are R p = 700 dyn • s • cm the pressure is constant and equal to 80 mmHg, whereas all the other state variables are zero. During the whole simulation, the stress received by the structure is corrected by the initial one, and the solid only responds to the difference with the diastolic phase. In this way, the load applied to the structure is kept reasonable small so that the linearity assumptions holds.

The results are summarized in Figure 7 (right), showing the Windkessel pressures in time for a time step P for τ = 0.001 s. Note that the semi-implicit algorithm with explicit 0D-3D coupling (Algorithm 4) is unstable whereas with the implicit 3D-0D treatment (Algorithm 3) the numerical solution remains stable. In fact, from the results one can infer that the interface term outlined in Remark 4.3, namely ( ẏn , P n n f ) Σ , injects a positive artificial power into the system (an increased pressure P n > 0 leads to ẏn • n f > 0).

Conclusions

In this paper we formulated and analyzed the treatment of a 3D-0D coupling between the 3D distributed (fluid, fluid-structure interaction) models and a set of Windkessel 0D models. The key ingredient in the proposed schemes is the use of a fractional-step time-marching in the 3D compartment. For purely fluid problems with multiple outlets, our energy based stability analysis showed that numerical instabilities might appear when using an explicit 3D-0D coupling. Interestingly, this result also holds with a single outlet in the case of fluid-structure interaction. We proposed to overcome this issues via an implicit treatment of the 3D-0D coupling, which involves a negligible additional computational cost with respect to the explicit strategy. These theoretical expectations were confirmed by numerical experiments in realistic geometries and physiological data.

Figure 3 :

 3 Figure 3: Patient-specific aorta. Left: geometry. Center: surface mesh. Right: simulation results for the aorta (vector velocity field and pressure distribution).

Figure 4 :

 4 Figure4: Mean pressures at inlet (black) and outlets 1 to 4 (blue, green, red, cyan). Explicit (left) and implicit (right) schemes.

Inria

Figure 5 :

 5 Figure 5: Fluxes for inlet (black) and outlets 1 to 4 (blue, green, red, cyan). Explicit (left) and implicit (right) schemes.

Figure 6 :

 6 Figure6: Pressures (left) and fluxes (right) for inlet (black) and outlets 1-4 (blue, green, red, cyan) obtained with a monolithic method (circles) and with the projection method (continuos) using an implicit 3D-0D coupling.

- 5 ,

 5 R d = 5 • 10 3 dyn • s • cm -5 and C = 2 • 10 -4 cm 5 • dyn -1 . At t = 0,

Figure 7 :

 7 Figure 7: Left: snapshot of the velocity field (left) inside the idealized AAA, the curves show outlet pressure and displacements of the aneurysm wall. Right: Windkessel pressure for the AAA, comparison between implicit (solid line) and explicit (dashed line) fluid-Windkessel coupling.

Table 1 :

 1 Parameters for the Windkessel's model. The outlets are ordered in direction of the flow.

	1 . A

RR n°7937

www.euheart.eu RR n°7937

RR n°7937RESEARCH CENTRE PARIS -ROCQUENCOURTDomaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

The authors are thankful to Jean-Frédéric Gerbeau and Romain Guibert (INRIA) for valuable discussions, Philipp Beerbaum and Israel Valverde (King's College London), for the acquisition of the patient data, and David Barber, Rod Hose and Cristina Staicu (University of Sheffield) for the segmented geometry. This work was partially founded by the Seventh Framework Program (FP7/2007-2013) of the European Community, under grant agreement 224495 (euHeart project).

Inria