W Ubachs
email: w.m.g.ubachs@vu.nl

M L Diouf

F M J Cozijn

B Darquié

E J Salumbides

Lamb-dips and Lamb-peaks in the saturation spectrum of HD

Lamb-dips and Lamb-peaks in the saturation spectrum of HD

Vibrational transitions in the hydrogen deuteride (HD) molecule, associated with the small permanent dipole moment of this hetero-nuclear species, were first detected by Herzberg [1]. Thereafter transition frequencies of a number of lines in the (2-0) overtone band were spectroscopically investigated early on [2,3], and over the decades at increasing accuracy [4]. A vast number of spectroscopic studies have been performed on the HD vibrational spectra, but all were limited by Doppler broadening. In two recent studies laser precision experiments were performed that, for the first time, saturated the absorption spectrum of the very weak R(1) line at 1.38 µm in the (2-0) overtone band in HD and acquired a Doppler-free resonance, resulting in orders-of-magnitude improved accuracies. In one study the method of noiseimmune-cavity-enhanced optical heterodyne molecular spectroscopy (NICE-OHMS) was employed [5], while in a second study a Lamb-dip was observed via cavity ringdown spectroscopy [6]. Where in both experiments the spectroscopy laser was locked to a frequency-comb laser for achieving accuracy at 10 -10 levels, the resulting transition frequencies deviated by 900 kHz, corresponding to 9σ discrepancy for the combined uncertainties. In view of the importance of such precision measurements for testing quantum electrodynamics in hydrogen molecular systems [7,8], for probing fifth forces [9] and extra dimensions [10], we have reinvestigated this R(1) line of HD at improved background suppression and signal-to-noise ratio, and under varying pressure conditions with the aim of finding the cause of this discrepancy.

The frequency-comb locked NICE-OHMS setup described in the previous study [5] is significantly modified on both the signal and frequency acquisition. It now involves a high-speed lock-in amplifier (Zurich Instruments, HF2LI), which allows parallel demodulation to extract all relevant signal components simultaneously as both the down-converted high-frequency modulation (f F SR ) and low-frequency dither (f W M ) all fall inside the lock-in amplifier bandwidth (50 MHz). A schematic layout of the experimental setup is displayed in Fig. 1.

FIG. 1. Schematic layout of the experimental setup used

in the present study. The ECDL laser is modulated twice by the EOM at both the FSR and PDH frequencies, of which the latter is used for short term stabilization onto the high-finesse cavity. Additionally, the laser is wavelengthmodulated through a dither on the cavity piezo for derivative detection. A slow lock with the frequency comb is performed for long term stabilization.

The spectroscopy laser is directly stabilized onto the high-finesse cavity through a Pound-Drever-Hall (PDH) lock for short term stabilization, while long term stability and accuracy (1 kHz) is obtained by a direct comparison with a Cs-stabilized frequency comb. To remove the effect of the periodic frequency dither (f W M = 415 Hz, 90 kHz peak-to-peak amplitude) in the acquisition, the counter operates at a gate time set to an integer value of the dither period. The measured beatnote value is used in a digital feedback loop to feedback the cavity piezo. The DeVoe-Brewer signal, required to stabilize the modulation frequency onto the cavity Free-Spectral-Range f F SR = 305 MHz, is also retrieved within the lock-in amplifier. It occurs at both frequencies of f F SR ± f P DH , which are located at 5 MHz or 45 MHz within the lock-in amplifier due to the down-conversion. To allow retrieval of this error function, it is required that both the internal oscillators of the lock-in amplifier and the PDH frequency generator are all stabilized to an external reference to ensure a fixed phase relationship.

Saturated absorption spectra of the HD R(1) line were measured as a function of pressure in the range 0.25 to 10 Pa, recordings of which are displayed in Fig. 2. With a cavity finesse of ∼ 150 000, the circulating power is estimated to be ∼ 100 W. The spectra represent the average of 8 scans, each taken at 20 minutes recording time (with 50-kHz frequency steps), yielding more than a five-fold improvement in signal-to-noise ratio with respect to the previous study [5]. It should be noted that the NICE-OHMS technique employed, which is essentially a form of frequency modulation spectroscopy in an optical cavity, yields a dispersive spectral line shape [START_REF] Axner | NICE-OHMS-frequency modulation cavityenhanced spectroscopy-principles and performance[END_REF][12][13][14]. However, the recordings are performed with an additional slow wavelength modulation at frequency f W M , which improves the signal-to-noise ratio and allows for detecting the extremely weak HD signal in saturation. Hence, after 1f demodulation at f W M , a derivative of the NICE-OHMS dispersion channel is expected in the form of a symmetric line shape. Indeed, a symmetric saturated absorption line is observed for C 2 H 2 (an R(9) line at 271 377 365 327 (5) kHz), recorded under the same NICE-OHMS modulation conditions and displayed in the right panel of Fig. 2. The C 2 H 2 line reflects the expected line shape, and the characteristic Lamb-dip of saturation spectroscopy. [5] (Amsterdam), [6] (Hefei), and [15] (Caserta), as well as the theoretical value from Ref. [16]. The shaded areas indicate the uncertainty estimates of the reported line positions.

The measurements performed on the HD R(1) line exhibit a variation in line profiles at different pressures. At low pressures (0.25-1 Pa) the line profile is dominated by a feature with a reversed sign when compared to the C 2 H 2 spectrum. These HD-signals at low pressure form Lamb-peaks, that were interpreted in the previous study [5] as the saturated line from which the centre frequency of the R(1) line was derived, at 217 105 181 895 (20) kHz.

At the lowest pressures (0.25 Pa), besides the marked enhanced absorption (peak) at low frequencies some reduced absorption at higher frequencies is observable. The reduced absorption feature, or a Lamb-dip, becomes more pronounced as a more intense and gradually broader signal with increasing pressure. With these two features, the resonance resembles a dispersive line shape. This deviation from an expected symmetric line shape is the central finding of the present experimental study.

While the Lamb-peak signal at low frequencies corresponds to the feature analyzed in Ref. [5], a broader feature, of the same sign as the C 2 H 2 feature is also observed. This feature, having the sign of a Lamb-dip, is located in the range of higher frequencies, where the Lamb-dip of the HD R(1) line was found in the cavity ring-down study, at 217 105 182.79(8) MHz [6]. These observations on pressure-dependent line shapes may form the basis for resolving the discrepancy between the two results previously published [5,6]. The line position reported in a recent study, a Doppler-limited measurement using frequency-comb assisted cavity ring-down spectroscopy [15], is also indicated in Fig. 2. In the following, an analysis of the line shape is presented that should quantitatively support this assumption, although a number of approximations and hypotheses must be made. Crucial for our proposed interpretation of the line shape is the underlying hyperfine structure of the R(1) transition. The nuclear spins for the deuteron I D = 1 and the proton I H = 1/2 give rise to a splitting in 5 hyperfine substates in the J = 1 ground state and into 6 hyperfine substates in the J = 2 excited state. This gives rise to 21 possible hyperfine components within the saturated absorption profile of the R(1) line. The level structure for the v = 0, J = 1 level follows from the analysis of RF-spectra by Ramsey and Lewis [17]. For the v = 2 level hyperfine coupling constants were calculated via ab initio theory, and from these the hyperfine level splitting in v = 2, J = 2 were derived and assumed to be accurate within 5% [18]. Relative line intensities of the 21 components were calculated via angular momentum algebra [18]. The calculated hyperfine structure is shown as a stick spectrum in Fig. 3. This graph also includes the locations of so-called cross-over resonances that couple non-zero velocity classes in a narrow saturation peak. Between the 5 ground state levels and the 6 excited state levels there exist 32 resonances where two ground state levels are connected to a common upper state (Λ-type cross overs), and 36 resonances that couple two excited levels to a common ground state (V-type cross overs). The occurrence of cross-over resonances is well documented in saturation spectroscopy [19][20][21][22]. They give rise to sign reversals, hence Lamb peaks, depending on the V-type or Λ-type connection of hyperfine levels [23,24].

To study the general features of the observed HD R(1) line profile, an approximate model was set up using the density matrix formalism that includes the 5 hyperfine sublevels of the v = 0, J = 1 ground level and the 6 hyperfine sublevels of the v = 2, J = 2 excited level. The coupled Bloch equations involving populations of excited sublevels ρ jj and ground sublevels ρ ii and the coherences ρ ij are:

d dt ρ ii = j ρ jj γ pop,ij - i h j (V ji ρ ji -ρ ji V ij ), (1)
d dt ρ jj = - i ρ jj γ pop,ij - i h i (V ij ρ ij -ρ ij V ji ), (2)
d dt ρ ij = -(i∆ ij + γ coh,ij )ρ ij + i h V ij (ρ ii -ρ jj ) - i h k =j V ik ρ kj , (3)
d dt ρ jk = (i∆ jk -γ coh,jk )ρ jk - i h (V ij ρ ji -ρ ki V ki ), (4)
with the detuning

∆ ij = ω L -(ω ij + k • v)
, where ω L is the laser frequency, ω ij is the transition frequency between states i and j, and k • v is the Doppler shift for the velocity class v. The population relaxation rates γ pop,ij connecting i and j states with allowed dipole transitions and relaxation rates for coherences γ coh,ij effectively describe radiative and non-radiative processes. The coupling V ij = µ ij (E + + E -) + c.c. consists of electric field contributions from forward +k and backward -k propagating laser beams with equal intensities. The electronic transition dipole moment as well as relative intensities of the hyperfine transitions are included in µ ij . In the calculation, the detuning ∆ of the applied laser field is scanned over the vicinity of the resonances, with the Doppler shift accounted for both +k and -k beams and for each of the n v velocity classes around v = 0. We have neglected the effect of the two laser sidebands at f F SR = ± 305 MHz (at field strengths of about 1% of the carrier), to simplify the treatment. The occurrence of a recoil doublet [25], with a splitting of 68 kHz for the R(1) line, is ignored as well.

Another complication in explaining the observed spectra lies in the reduction of transit-time broadening. We note here, as was discussed in [5], that the observed linewidth of the HD R(1) line is observed much narrower than expected from transit time broadening for a roomtemperature velocity distribution. This was attributed to the mechanism of selection of cold molecules under conditions of very weak saturation, which is a known phenomenon [26][27][28]. In the model description no attempt was made to explicitly explain this behavior.

The coupled Bloch equations were solved by numerical integration using a Python code and associated numerical and scientific libraries. The computation was performed in a cluster computer utilizing 32 nodes with 16 cores. For n v ∼ 4400 velocity classes that cover Doppler shifts in the range of [-1.8, 1.8] MHz around each resonance, a calculation of a spectrum typically takes around 18 hours.

A first calculation was performed that was restricted to the v = 0 velocity fraction (v = 0 along the propagating laser beams). Its result shows a spectrum displaying all hyperfine resonances at the expected positions indicated as proper hyperfine resonances in Fig. 3 as Lambdips. Subsequently, more elaborate calculations were performed, integrating the spectrum over extended velocity classes. Depending on the parameters invoked for γ coh and γ pop , this treatment produces cross-over resonances with the same sign as a Lamb-dip, but also cross-over features with a reversed-sign. The simulations reveal the effects of velocity selective optical pumping that are observed at particular velocity classes where the cross-over interaction between the counter-propagating beams occur. From a parameter-space analysis it is found that the sign-reversed cross-over peaks occur only under certain values for the input parameters. Clear anti-crossovers are produced only when the population decay γ pop occurs at a faster rate than the Rabi frequency that was estimated at V /h = 20 kHz in the measurements [5].

This suggests that there must be a mechanism transferring population from the excited v = 2 state to the lower v = 0 state in the molecule, a refilling mechanism. This rate should be much larger than the radiative lifetime of the upper HD level, being on the order of Γ ∼ 1 Hz. Hence, in this specific case of a saturated transition with a long upper state lifetime in a molecule, similar to the case of NH 3 [29], the population decay (γ pop ) cannot be induced by spontaneous radiative decay, as in the case of atoms [20,21,24].

In order to effectively produce anti-crossovers or Lamb peaks an hypothesis must be invoked as to the refilling mechanism for ground state population. In some exam-ples, that bear some similarity, the main contribution to the rate of population decay γ pop was attributed to optical pumping in a multi-level system, with the transittime effect as the absorbers traverse the beam playing a role [23,24]. In a model for saturation in coherent anti-Stokes spectroscopy in H 2 the mechanism of velocitychanging collisions was adopted as a possible transfer mechanism [30]. The decay value in the latter CARS experiment was close to the value of the pressure broadening coefficient of γ P = 35 kHz/Pa for saturated absorption extracted in Ref. [5], while Ref. [15] reports a coefficient of γ P = 10 kHz/Pa for a Doppler-limited study (without saturation). For the pressure P ∼ 1 Pa accessed in the measurements, the collision rate contribution of γ P P to γ coh,ij is ∼ 10 -35 kHz. Without further speculation on the nature of the mechanism, a range of values for the refilling rate was used in solving the density matrix framework. Figure 4 shows numerical results obtained for conditions that mimic the observations of spectra at different pressures. The line profiles qualitatively agree with experimental observations. In particular the Lamb-peak and dip features correspond to the measurements, although the intensities and widths are not well reproduced. Both the experiments and simulations demonstrate that the Lamb-peak persists at the lowest pressures (lowest values of γ coh,ij ), while the Lamb-dip decreases in intensity at the low pressures.

The simulations show that the saturation line profiles of HD R(1) exhibit a composite line shape with a Lambpeak at low frequency and a broader Lamb-dip occurring at higher frequencies. We postulate this as an explanation for the results obtained in previous studies, where a low-frequency Lamb-peak was analyzed [5]. We assume that the high-frequency component was measured in the cavity ring down study [6], presumably at some-what higher pressures. We emphasize that a number of approximations were made in the numerical treatment which could affect the intensities and widths, such as the neglect of the recoil doublet, the inclusion of approximate results on the hyperfine structure from calculations, the lack of a quantitative description of the mechanism selecting cold molecules resulting in larger transit times, absence of cavity standing-wave effects, and in particular the lack of understanding of the refilling mechanism.

A comparison can be made between the observed line shapes and the composite profiles from the numerical treatment, i.e. in a fit. Since the numerical analysis is based on the hyperfine components dispersed at wellknown frequencies, such a fit results in a value for the hyperfineless HD R(1) rovibrational transition. It is interesting to note that a pressure-induced shift of -10 kHz/Pa for the Lamb peak line position in the simulations (Fig. 4) is close to that extracted in the experimental results. By taking into account profiles for various pressures in the simulations, the extrapolated collisionfree Lamb peak line center can be related to the hyperfineless transition frequency. Adopting this treatment to the measurements yields a value of 217 105 181 901 kHz for the R(1) transition frequency. In view of the numerous assumptions made and the pressure-dependent line shape we take a conservative estimate for the uncertainty of the hyperfineless transition frequency of some 50 kHz. The analysis shows also that the zero-pressure extrapolation for the Lamb-peak, as determined in [5], is rather close to the hyperfine center-of-gravity to within 20 kHz. The present result agrees, within its 50 kHz uncertainty, with our previous determination [5].

The presently obtained transition frequency may be compared to the theoretical values obtained from advanced molecular QED calculations. While a calculation of early 2018 yielded a value of 217 105 174.8 (1.8) MHz, included in [6], deviating by 7 or 8 MHz from both experimental values, recently a theoretical value was reported at 217 105 180.2 (0.9) MHz [16]. The latter value, that was based on an improved treatment of relativistic corrections in non-adiabatic perturbation theory, is within 2σ from the presently obtained experimental value. This allows for the conclusion that agreement is obtained between experiment and theory of infrared transitions in hydrogen at the 5 × 10 -9 level, where the present result is of higher accuracy than current state-of-the art calculations. Meanwhile a theoretical framework has been developed that should produce a more accurate value [8], to be confronted with the experimental value. On the experimental side the observation of the R(0) line in the (2-0) band, exhibiting a much simpler hyperfine structure with a near-degenerate ground state [31], could be measured. Unfortunately that line is overlaid by a water resonance prohibiting its measurements in the current setup.

The authors thank J.M.A. Staa (VU) for assisting in the measurements and P. Dupré (Dijon) for making available calculations on the hyperfine structure. The work

FIG. 2 .

 2 FIG.2. Observed 1f -spectra of the HD R(1) transition under conditions of varying pressure are shown in the left panel. For comparison in the right panel, an R(9) acetylene line (f0 = 217 377 365 327 kHz) recorded at 0.1 Pa with identical settings of modulation and demodulation phases for the HD measurements. The acetylene spectrum reflects the expected 1f -derivative of the dispersion profile for a saturated NICE-OHMS Lamb dip. The vertical lines in the left panel indicate the line positions reported in Refs.[5] (Amsterdam),[6] (Hefei), and[15] (Caserta), as well as the theoretical value from Ref.[16]. The shaded areas indicate the uncertainty estimates of the reported line positions.

FIG. 3 .

 3 FIG. 3. A stick spectrum representing the 21 hyperfine components (in blue) and the 68 cross-over resonances (in red and green) in the R(1) line of HD. The zero on the x-axis represents the location of the hyperfineless rovibrational transition. The line strengths of the cross-over resonances represent the root of the product of hyperfine line intensities of the two connecting transitions.

FIG. 4 .

 4 FIG. 4. Comparison between observed and modeled spectral line shapes of the R(1) line of HD; (a) Experimental spectra for 0.25, 0.5, 1.0 and 10 Pa on an absolute frequency scale; (b) Calculated spectra on a relative frequency scale with respect to the hyperfineless point (at 0 MHz) for different values of γ coh and γpop in units of kHz and the Rabi frequency fixed at 20 kHz.

was supported by the Netherlands Organisation for Scientific Research (NWO) via the program "The mysterious size of the proton". SURFsara (www.surfsara.nl) is acknowledged for the support in using the Lisa Compute Cluster. WU acknowledges the European Research Council for an ERC-Advanced Grant under the European Unions Horizon2020 research and innovation programme (Grant agreement No. 670168).


Éric Canon

Frédéric Chardard

Grigory Panasenko

Olga Štikonienė

Numerical solution of the viscous flows in a network of thin tubes: equations on the graph

Keywords: Navier-Stokes Equations, Asymptotic models, weakly singular kernels, multi-scale problems. AMS classification: 78M35, 74G10, 35Q30, 35R09, 34K28. 1

A non-stationary flow in a network of thin tubes is considered. Its one-dimensional approximation was proposed in a paper by G.Panasenko and K.Pileckas, Flows in a tube structure: equation on the graph, JMP (2014). It consists of a set of equations with weakly singular kernels, on a graph, for the macroscopic pressure. A new difference scheme for this problem is proposed. Several variants are discussed. Stability and convergence are carefully investigated, theoretically and numerically. In addition, numerical results are compared to the direct numerical solution of the full dimension Navier-Stokes equations.

Introduction

Modeling and effective computations of viscous flows in networks of thin pipes are important topics, notably in microfluidics and flows in blood vessels. Because of the complexity of the geometry of such 2D or 3D networks, it is important to have simplified, but relevant, models for these flows, and then efficient schemes for solving them numerically. The main aim of the present work is to design and study numerical schemes for such a simplified model which was obtained by Panasenko & Pileckas, and consists in a time dependent, non local in time diffusion equations for the pressure, set on a 1D-graph, with appropriate (Kirchhoff) junction conditions on the vertices of the graph. It is described in [START_REF] Panasenko | Flows in a tube structure: equation on the graph[END_REF] and mathematically justified in [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF], by letting the diameters of the pipes tend to zero, with appropriate scalings on the data. Also, we compare numerically the original multidimensional problem with its simplified version, in order to establish the relevance of this limit model [START_REF] Panasenko | Flows in a tube structure: equation on the graph[END_REF].

Let us go first more into detail about this modeling. Stationary laminar incompressible viscous flows in a thin tube were first studied empirically [START_REF] Louis | Recherches sur les causes du mouvement du sang dans les vaisseaux capillaires[END_REF] and theoretically [START_REF] Boussinesq | Essai sur la théorie des eaux courantes[END_REF] in the nineteenth century. They were described by the stationary Poiseuille profile. However, the time scale of the phenomena may prevent stationary approximation to be valid as noticed by Womersley [START_REF] John R Womersley | Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known[END_REF] for blood flows. In that case, one needs to use a time-dependent Poiseuille type profile: a Womersley profile. Pileckas et al. [START_REF] Pileckas | Existence of a nonstationary Poiseuille solution[END_REF][START_REF] Giovanni | On the unsteady Poiseuille flow in a pipe[END_REF][START_REF] Panasenko | Asymptotic analysis of the nonsteady viscous flow with a given flow rate in a thin pipe[END_REF] studied the behaviour of the incompressible Navier-Stokes equations in an infinite tube when the diameter of the pipe is small, proving the existence and relevance of the time-dependent Poiseuille profile in a quite general setting. Panasenko & Pileckas studied the incompressible Navier-Stokes equation in a network of thin tubes, when the flow is stationary [START_REF] Panasenko | Partial asymptotic decomposition of domain: Navier-Stokes equation in tube structure[END_REF][START_REF] Panasenko | Asymptotic expansion of the solution of Navier-Stokes equation in tube structure and partial asymptotic decomposition of the domain[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure. I. The case without boundary-layer-in-time[END_REF][START_REF] Panasenko | Divergence equation in thin-tube structures[END_REF] or timedependent [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure. I. The case without boundary-layer-in-time[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF]. In that case of non-steady Navier-Stokes equations in the tube-structure, similarly to [START_REF] Pileckas | Existence of a nonstationary Poiseuille solution[END_REF][START_REF] Giovanni | On the unsteady Poiseuille flow in a pipe[END_REF][START_REF] Panasenko | Asymptotic analysis of the nonsteady viscous flow with a given flow rate in a thin pipe[END_REF], they obtain ( [START_REF] Panasenko | Flows in a tube structure: equation on the graph[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF]) a limit, simplified model retaining this Womersley effect. The model consists in a one dimensional, nonlocal in time, diffusion equation on the 1Dskeleton (or graph) of the initial structure for the pressure, coupled with heat equations set in rescaled cross sections of the tubes of this initial structure, with Kirchhoff type conditions at the vertices of the graph and the continuity condition for the pressure.

The main aim of the present paper is to propose and study a numerical scheme, with several variants, for this model. We show convergence and stability without any special restriction on the shape of the cross sections of the tubes. An important point in designing the schemes is that in the coupling mentioned above the heat equations in the cross sections can be solved explicitly with respect of the gradient of the pressure, allowing to rewrite the equations for the pressure on the edges as diffusion equations with convolution in time. In designing the schemes, we have to pay a particular attention on how to discretize this couplings, or convolution products. Different possibilities are investigated : exact kernels for special geometry of the cross-sections, numerical approximation on the kernels for general cases, with order 1 or 2 in time, with of without corrections for small times, using asymptotics for small times obtained in the second part of this work [START_REF] Canon | Numerical solution of the viscous flows in a network of thin tubes: equations on the cross-section[END_REF].

At this point, let us underline that this work is in two parts: properties of the kernels are investigated in detail in a forthcoming paper ([8]). However, as mentioned above, some of these properties are used in the numerics of the present paper. Last, we also check the convergence numerically and compare the asymptotic model with the directly numerically solved Navier-Stokes equations.

Let us mention well-known alternative models for network of pipes, considered for instance in [START_REF] Ashmetkov | A boundary value problem for the linearized haemodynamic equations on a graph[END_REF][START_REF] Borzov | Conservative algorithm of substance transport over a closed graph of cardiovascular system[END_REF][START_REF] Formaggia | One dimensional models for blood flow in arteries[END_REF][START_REF] Gamilov | Virtual fractional flow reserve assessment in patient-specific coronary networks by 1D hemodynamic model[END_REF][START_REF] Robertson | A director theory approach for modeling blood flow in the arterial system: An alternative to classical 1D models[END_REF][START_REF] Spencer | Computational modelling of 1D blood flow with variable mechanical properties and its application to the simulation of wave propagation in the human arterial system[END_REF][START_REF] Vassilevski | Numerical issues of modelling blood flow in networks of vessels with pathologies[END_REF]. These models are hyperbolic on the graph and are designed for a different context (high Reynolds number, elastic walls, as in arteries) than ours (modest Reynolds number, rigid walls as in arterioles). See [START_REF] Panasenko | Flows in a tube structure: equation on the graph[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF] and the just mentioned references for detail. Now, the plan of the paper is as follows. In Section, we recall the original and limit models in [START_REF] Panasenko | Flows in a tube structure: equation on the graph[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF] and give a brief presentation of our scheme. In Section we present a few notations and the weak form of the equations on the graph and prove a few additional regularity results that are needed in this paper. In Section 3 we go more into detail about the discretization: the scheme as presented in Section 1 is derived from the weak form as a Galerkin method; we show coercivity (and continuity) of that approximate bilinear form; we give a few words on the discretization in the cross sections (the details been given in the second part of this work [START_REF] Canon | Numerical solution of the viscous flows in a network of thin tubes: equations on the cross-section[END_REF]); and last some remarks on the implementation of the scheme. Section 4 is devoted to the main results : convergences of the scheme and rates of convergence, in relation with the accuracy on the approximation of the kernels: Theorems 1 and 2 and Corollary 1. Section 5 is devoted to the presentation of numerical results : comparison with exact test cases and comparisons with full 3D Navier-Stokes equations.

Thin tubes structure

We first describe a d-dimensional thin tubes structure, where d ∈ {2, 3}. Let O 1 , . . . O N be different vertices of R d , e 1 , . . . , e M closed segments e j = [O ij O kj ] connecting these vertices. We assume that these edges can only intersect at vertices. The end points of the skeleton of the graph, that is the vertices that belong to a single edge are denoted by O 1 , . . . , O N1 (0 ≤ N 1 ≤ N -1) while O N1+1 , . . . , O N corresponds to junctions between the pipes. We fix a positive orientation for all edges e j = [O ij O kj ] as the direction from O ij to O kj . For simplicity, we use the same notation e j for the vector ----→ O ij O kj as for the edge. Let σ 1 , . . . , σ M be smooth domains of R d-1 representing the cross-sections of the pipes. The variable in the cross-section direction is denoted with a hat: x ∈ R d-1 . Let R (e1) , . . . , R Let ω N1+1 , . . . , ω M be bounded domains of R d containing 0 which represent the junctions of the tubes and define the initial d-dimensional network of tubes as :

Ω ε = N i=N1+1 (O i + εω i )   N j=1 O ij + R (ej ) (εσ j ×]0, |e j |[)   .
We assume that, for ε small enough, Ω ε is connected, smooth except at the end sections

N1 j=1 R (ej ) (ε∂σ j × {0}). Last, B = M j=1 e j
is the graph of the structure.

We now, states the equations of the model we are dealing with. These equations are set on the graph B and were obtained in [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF] by passing to the limit as ε tends to 0 on the following incompressible Navier-Stokes equations with Dirichlet's boundary conditions, in Ω ε :

                           ∂u ε ∂t + (u ε • ∇u ε ) -∆u ε = -∇p ε on Ω ε , ∇ • u ε = 0 on Ω ε , u ε (O i + R (ej ) (x, 0), t) = R (ej ) g i (ε -1 x, ε -2 t) on O i + R (ej ) (εσ j × {0}), i = 1, ..., N 1 , O i ∈ e j , u ε = 0 on ∂Ω ε \ N1 i=1 (O i + R (ej ) (εσ j × {0})| Oi∈ej ) , u ε | t=0 = 0, p ε (O 1 , t) = 0.
(1.1) Here, g i : σ i × R + → R d denotes the velocity profile in a reference frame at the free extremity of the tube i for i ∈ {1, . . . , N 1 }. Let ∂ ∂x (ej ) denote the derivative along the vector |e j | -1 e j . According to [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure. I. The case without boundary-layer-in-time[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF], under some hypotheses on g i , p ε (x, ε 2 t) can be approximated at the leading order by P (x, t) where P is the solution to the following problem on B:

                 - ∂ ∂x (ej ) L (σj ) ∂P ∂x (ej ) (x, t) = F (x, t) for x ∈ e j , Oi∈ej α i,j L (σj ) ∂P ∂x (ej ) (O i , t) = -Ψ i (t) for i = 1, . . . , N, P is continuous on the graph, P (O 1 , t) = 0. (1.2)
In (1.2) 2 : the summation is over all edges e j having O i as an end point; α i,j = 1 if the orientation e j starting from O i is positive, and α i,j = -1 otherwise. The right hand sides in (1.2) 1 are given by :

Ψ i (t) =    σj g i (x, t), ( 0, α i,j ) dx if i ∈ {1, . . . , N 1 }, O i ∈ e j 0 if i ∈ {N 1 + 1, . . . , N },
where •, • denotes the Euclidean product.

The operator L (σj ) are nonlocal in time operators relating the macroscopic equation for the pressure (1.2) to the microscopic structures σ j . For each such cross section σ they are defined:

L (σ) : L 2 (0, +∞) → H 1 0 (0, +∞) by: ∀t > 0, (L (σ) q)(t) = σ U (x, t)dx (1.3)
where U is solution of

       ∂U ∂t -∆ xU = q(t) on σ, U = 0 on ∂σ, U = 0 if t = 0 on σ. (1.4)
Remarks on (1.2).

(i) In usual applications F = 0, but it is useful, at least for building test problems with exact solutions, to consider more general F .

(ii) Equation(1.2) 2 is the Kirchhoff condition which expresses the conservation of mass at the junctions and at the ends of the graph. If the system has a solution, it satisfies the following compatibility condition:

∀t ∈ [0, T ], N =1 Ψ (t) + B F (x, t)dx = 0. (1.5)
where

B
denotes the integral over all edges of the graph (see (2.13)). If this compatibility condition is satisfied then, the equation on Ψ 1 is redundant and can be omitted in the discretization.

(iii) Pressure can be obtained by solving system (1.2) and then, the velocity profile at some distance from the junctions by solving (1.4). This is discussed in Section 3.5.3. To obtain the next order of approximation, the continuity of pressure (1.2) 3 may be replaced by a coupling with Navier-Stokes equations in regions surrounding the junctions (see [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF]). Such couplings are also discussed and numerically tested in [START_REF] Bertoglio | A Stokes-residual backflow stabilization method applied to physiological flows[END_REF][START_REF] Bertoglio | Benchmark problems for numerical treatment of backflow at open boundaries[END_REF][START_REF] Bertoglio | Fractional-step schemes for the coupling of distributed and lumped models in hemodynamics[END_REF].

Numerical scheme and convergence

The main goal of this article is to propose a numerical method to solve problem (1.2). This is done in the present section. First, the operator L (σ) can be rewritten as a convolution operator:

(L (σ) q)(t) = t 0 K (σ) (t -τ )q(τ )dτ, (1.6)
where the kernel K (σ) is defined by

K (σ) (t) = σ V (x, t)dx,
V being the solution of the heat equation:

       ∂V ∂t -∆V = 0 on σ×]0, +∞[, V (., 0) = 1 on σ, V (x, t) = 0 t > 0, x ∈ ∂σ.
(1.7)

The discretization is then as follows. Let (t 0 = 0, t 1 = k, . . . , t Q = kQ = T ) be a subdivision of [0, T ] with a uniform step k. For each edge e, we introduce a subdivision (x

(e) 0 = 0, x (e) 1 = h (e) , . . . , x (e)
S (e) = S (e) h (e) = |e|) of [0, |e|] with a uniform step h (e) . Let X (ei)

s = O i + R (e) ( 0, x (e) s ). Let P (e)
s,q be an approximation of P (X (e) s , t q+1/2 ), where t q+1/2 = t q + t q+1 2 . Consider the following difference scheme:

f (ej ) s,q+1 = -k q q=0 K(σj) q-q P (ej ) s+1,q -2P (ej ) s,q + P (ej ) s-1,q (h (ej ) ) 2 if      1 ≤ j ≤ M, 0 < s < S (ej ) 0 ≤ q < Q , (1.8)
Ψ ,q+1 = 1≤j≤M,O ∈ej   - h (ej ) 2 f (ej ) s,q+1 -k q q=0 K(σj) q-q P (ej )
s+α ,j ,q -P (ej ) s,q

h (ej )   if 2 ≤ ≤ N 0 ≤ q < Q , (1.9)
P (e) s,q = P (ẽ) s,q if e, ẽ ∈ {e 1 , . . . , e M }, s ∈ {0, S (e) }, s ∈ {0, S (ẽ) }, X (e) s = X (ẽ) s , 0 ≤ q < Q (1.10)

P (e) 0,q = 0 if X (e) 0 = O 1 , 0 ≤ q < Q, (1.11)
where:

• K(σ) q is either K (σ) q = 1 k tq+1 tq K (σ) (t)
dt for q ≥ 0 when the kernel can be computed explicitly (see (3.22)-(3.25)) or a suitable approximation of it (see Section 3.4).

• Ψ ,q = Ψ (t q ) for 0 ≤ q ≤ Q and 1 ≤ ≤ N ,

• f (ej ) s,q = F (X (ej ) s , t q ). • in (1.9), s = 0, α ,j = 1 if X (ej ) 0 = O and s = S (ej ) , α ,j = -1 if X (ej )
S (e j ) = O . Note that X (ej ) s+α ,j is the node closest to O among those on the edge e j .

A similar numerical scheme can be found in [START_REF] Mette S Olufsen | Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions[END_REF][START_REF] Qureshi | Numerical simulation of blood flow and pressure drop in the pulmonary arterial and venous circulation[END_REF][START_REF] Qureshi | A computational study of pulmonary hemodynamics in healthy and hypoxic mice[END_REF] when B is a tree, the cross-section a disc, and all the data are periodic in time. However, the convergence is not addressed in [START_REF] Mette S Olufsen | Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions[END_REF][START_REF] Qureshi | Numerical simulation of blood flow and pressure drop in the pulmonary arterial and venous circulation[END_REF][START_REF] Qureshi | A computational study of pulmonary hemodynamics in healthy and hypoxic mice[END_REF].

Let:

θ(k) = max 1≤j≤M |K (σj ) 0 - K(σj) 0 | + Q-1 q=1 |K (σj ) q -K (σj ) q-1 -K(σj) q + K(σj) q-1 |. (1.12) Note that θ(k) = 0 if K(σ) q = K (σ) q .
In Section 3 below, we first reformulate this scheme as a Galerkin method. In Propositions 1 and 2 of Section 3.2, we then give two different sufficient stability conditions on K(σ) q , which are verified when θ(k) = 0. Last in Section 4, we prove our main result: the scheme is convergent if θ(k) → 0 and is of L 2order 1 in time and H 1 -order 1 in space if θ(k) = 0. This is Theorem 1. Furthermore, with additional regularity assumptions, these orders of convergence can be improved. This is Theorem 2 and Corollary 1.

Continuous problem on the graph

In this section, we get results on the continuous problem (1.2) that are used in the proof of Theorem 1.

Vector spaces and notations

To avoid summations, let us denote:

B f (x)dx = M j=1 |ej | 0 f (O i + R j ( 0, x d ))dx d .
(2.13)

The operators L (σj ) and ∂ ∂x (ej ) will be applied to functions which are defined on B. To shorten notations and to be able to use the previous notation, let us denote σ(x) = σ j and e(x) = e j when

x ∈ e j and x is not a vertex, so that ∂f ∂x (e) (x) = ∂f ∂x (ej ) (x) and L (σj ) g(x, t) = L (σ) g(x, t). Define the norms and the vector spaces:

f 2 L 2 (B) = B |f | 2 dx, f 2
H k dc (B) = B |f | 2 + ∂f ∂x (e) 2 + . . . + ∂ k f ∂x (e) k dx, L 2 (B) = {f : B → R | f L 2 (B) < +∞}, H 1 (B) = {f : B → R | f (O 1 ) = 0, f continuous, f H 1 dc (B) < +∞}, ∀k ∈ N * , H k dc (B) = {f : B → R | f 2 H k dc (B) < +∞}.
Notice that a function of H k dc (B) may be discontinuous at the junctions O 1 , O 2 , . . . , O . For any Banach space E and r ∈ N * , let us define:

H r (0, T, E) = f : [0, T ] → E| T 0 |f (t)| 2 E + |f (t)| 2 E + . . . + |f (r) (t)| 2 E dt < +∞ , H r 00 (0, T, E) = f ∈ H r (0, T, E)|f (0) = 0, . . . , f (r-1) (0) = 0 , H r (0, T ) = H r (0, T, R), H r 00 (0, T ) = H r 00 (0, T, R).

Weak form

Let us assume that F ∈ H 1 00 (0, T, L 2 (B)) and that Ψ 1 , . . . , Ψ N ∈ H 1 00 (0, T ) and satisfy the compatibility condition (1.5). The well-posedness in L 2 (0, T, H 1 (B)) of system (1.2) was proved by means of the Lax-Milgram lemma in [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure. I. The case without boundary-layer-in-time[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF]. The weak form of (1.2) is as follows: find P ∈ L 2 (0, T, H 1 (B)) such that:

∀ψ ∈ L 2 (0, T, H 1 (B)), a T (P, ψ) = b T (ψ) (2.14)
where

a T (P, ψ) = T 0 B ∂ 2 (L (σ) P ) ∂τ ∂x (e) ∂ψ ∂x (e) dxdt, (2.15)
b T (ψ) = T 0 B ∂F ∂τ ψdxdt + T 0 N =1 ∂Ψ ∂τ ψ(O , t)dt.
According to [START_REF] Panasenko | Flows in a tube structure: equation on the graph[END_REF], b T is continuous on L 2 (0, T, H 1 (B)), a T is continuous and coercive on L 2 (0, T, H

T,s (ψ) = T 0 B ∂ s+1 F ∂τ s+1 ψdxdt + T 0 N =1 ∂ s+1 Ψ ∂τ s+1 ψ(O , t)dt.
Now let P s be the solution of a T (P s , .) = b T,s (.). Then, using the uniqueness of the solution to Problem (2.14), we get that P s (., t) = t 0 P s+1 (., τ )dτ whenever 0 ≤ s < k.

Indeed, for any f ∈ L 1 (0, T ), letting F (t) = t 0 f (τ )dτ , using integration by parts and F (0) = 0, we get, for any σ ∈ {σ 1 , . . . , σ M }:

L (σ) F (t) = t 0 L (σ) f (τ )dτ.
Since P k ∈ L 2 (0, T, H 2 dc (B)), we conclude that P ∈ H k 00 (0, T, H 2 dc (B)).

Galerkin method

In this section, we reformulate the scheme (3.4.2) as a Galerkin method in order to prove the convergence result.

Time discretization

Let k > 0 be a time step, such that Q = T k is an integer. Let us denote t q = qk for q ∈ {0, 1, . . . , Q},

V k = P 0 k (0, T, H 1 (B)) ⊂ L 2 (0, T, H 1 (B)
) where P 0 k (0, T ) denotes the set of piecewise constant functions over the subdivision (t 0 = 0, . . . , t Q = T ). In the sequel the funtions K (σ) , L (σ) (φ q ), L (σ) (φ q ) are extended by zero on R * -for σ ∈ {σ 1 , . . . , σ M }. A numerical scheme on the graph is then defined as follows: find p ∈ V k such that:

∀ψ ∈ V k , a T,k (p, ψ) = b T (ψ) (3.17)
where a T,k is an approximation of a T to be specified below.

We use separation of space and time variables. Let us denote φ q = 1 [tq,tq+1[ . Then, for any u, v ∈ H 1 (B), q, q ∈ {0...Q}, We get:

a T (uφ q , vφ q ) = tq+1 tq B ∂(L (σ) φ q ) ∂τ ∂u ∂x (e) ∂v ∂x (e) dxdτ = B (L (σ) φ q )(t q+1 ) -(L (σ) φ q )(t q ) ∂u ∂x (e) ∂v ∂x (e) dxdτ = t q +1 t q B K (σ) (t q+1 -τ ) -K (σ) (t q -τ ) ∂u ∂x (e) ∂v ∂x (e) dxdτ
Denoting for any σ ∈ {σ 1 , . . . , σ M }:

K (σ) q = 1 k tq+1 tq K (σ) (t)dt for q ≥ 0 and K (σ) q = 0 for q ≤ 0, then: a T (uφ q , vφ q ) = k B (K (σ) q-q -K (σ) q-q -1 ) ∂u ∂x (e)
∂v ∂x (e) dx. Note that:

• K (σ) q -K (σ) q-1 = k 1 -1 (K (σ) ) (k(q + t))(1 -|t|)dt if q > 1, • K (σ) 0 -K (σ) -1 = K (σ) 0 , • K (σ) q -K (σ)
q-1 = 0 if q < -1. Now, the approximation a T,k of a T is any bilinear form defined on V k by:

a T,k (uφ q , vφ q ) = k B ( K(σ) q-q - K(σ) q-q -1 ) ∂u ∂x (e)
∂v ∂x (e) dx (3.18) such that the K(σ) q are approximations of the K (σ) q satisfying K(σ) q = 0 for q ≤ 0. Note in addition that:

b T (vφ q ) = B (F (x, t q+1 ) -F (x, t q ))v(x)dx + N =1 (Ψ (t q+1 ) -Ψ l (t q ))v(O ).
The time discretized problem (3.17) has a unique solution if a T,k is continuous and coercive.

Continuity and coercivity of the approximate bilinear form

We give two different sufficient conditions for a T,k to be continuous and coercive: one which relies on the smallness of θ(k) defined in (1.12); and another one which relies on properties of ( K(σ) q ) only. These conditions can be seen as stability conditions for the numerical scheme and will be used to prove Theorems 1 and 2.

Lemma 2. Let (A q ) 0≤q≤Q-1 be given real numbers, and a be a bilinear form defined by: a :

       P 0 k (0, T ) × P 0 k (0, T ) → R (b, c) = Q-1 q=0 b q φ q , Q-1 q=0 c q φ q → k Q-1 q=0 b q q q=0
A q c q-q .

Then:

|a(b, c)| ≤ b L 2 c L 2 Q-1 q=0 |A q |.
Proof: Let a, b, c be as in the lemma. Using a change of variable and Cauchy-Schwarz inequality, we get:

|a(b, c)| = Q-1 q=0 A q k Q-1 q=q b q c q-q ≤ k Q-1 q=0 |A q | Q-1 q=0 b 2 q 1/2 Q-1 q=0 c 2 q 1/2 = Q-1 q=0 |A q | b L 2 c L 2 .
Lemma 3. The bilinear form a T,k defined above is continuous and satisfies, for any u, v in P 0 (0, T, H 1 (B)):

|a T,k (u, v)| ≤ max 1≤i≤M Q-1 q =0 | K(σi) q - K(σi) q -1 | u L 2 (0,T,H 1 (B)) v L 2 (0,T,H 1 (B)) .
Proof:

Let u = Q-1 q=0 u q φ q , v = Q-1 q=0 v q φ q ∈ P 0 (0, T, H 1 (B)). Then, |a T,k (u, v)| ≤ M i=1 k Q-1 q=0 q q =0 | K(σi) q - K(σi) q -1 | ei ∂u q ∂x (ei) ∂v q-q ∂x (ei) dx ≤ M i=1 k Q-1 q=0 q q =0 | K(σi) q - K(σi) q -1 | u q H 1 (ei) v q-q H 1 (ei) .
Using Lemma 2 and Cauchy-Schwarz inequality, we get:

|a T,k (u, v)| ≤ M i=1 Q-1 q =0 | K(σi) q - K(σi) q -1 | u L 2 (0,T,H 1 (ei)) v L 2 (0,T,H 1 (ei)) ≤ max 1≤i≤M Q-1 q =0 | K(σi) q - K(σi) q -1 | M i=1 u L 2 (0,T,H 1 (ei)) v L 2 (0,T,H 1 (ei)) ≤ max 1≤i≤M Q-1 q =0 | K(σi) q - K(σi) q -1 | u L 2 (0,T,H 1 (B)) v L 2 (0,T,H 1 (B)) .
As a direct application of the lemma, we get the following estimate on the approximation of a T .

Proposition 1. The following estimate holds :

a T,k -a T ≤ θ(k)
where the norm is the norm for bilinear forms on

L 2 (0, T, H 1 (B)). In addition, if a T is α-coercive with α > θ(k), then a T,k is (α -θ(k))
-coercive, and continuous.

In Proposition 1, α depends on T . Hence, for a given k, a T,k may not be coercive for large T . The following lemma gives a condition under which a T,k is coercive for any k and T . It is a refinement of the Lemma 4.2 of [START_REF] Mclean | Numerical solution of an evolution equation with a positivetype memory term[END_REF].

Lemma 4. Let k, C, E, T m ∈ R + * and (K q ) q∈Z such that:

• k < min{T m , T }, • 0 ≤ K q ≤ C if q ≥ 0, • K q = 0 if q < 0, • K q+1 -2K q + K q-1 ≥ 0 if q ≥ 1, • E ≤ K q+1 -2K q + K q-1 k 2 if T m ≤ qk ≤ 2T m .
Let a be the bilinear form defined as in Lemma 2 with

A q = K q -K q-1 . Then |a(b, c)| ≤ C b L 2 c L 2 , a(b, b) ≥ α T b 2 L 2 ,
where:

α T = 1 2 T 2 m E min T 2 m 16T 2 , 1 40 .
Remark: The assumptions on (K q ) are discrete analogues of the fact that the (K (σ) ) are nonnegative and bounded, (K (σ) ) are nonnegative and the (K (σ) ) is larger than some constant E > 0 on an interval [T m , 2T m ], for any σ ∈ {σ 1 , . . . , σ M }. This fact arises straightforwardly from the generalized Dirichlet series expression in Lemma 2 in [START_REF] Giovanni | On the unsteady Poiseuille flow in a pipe[END_REF] (see also the forthcoming second part of this article).

Proof: (K q+1 -K q ) q≥0 is increasing and bounded by 2C and therefore convergent. Its limit is zero, otherwise (K q ) would not be bounded. Thus, (K q+1 -K q ) q≥0 is negative and so (K q ) q≥0 is decreasing. We have, from Proposition 1 that

|a(b, c)| ≤ b L 2 c L 2 (K 0 + Q-1 q=1 |K q -K q-1 |) ≤ b L 2 c L 2 (K 0 + K 0 -K Q-1 ) ≤ 2C b L 2 c L 2 . Let b = Q-1 q=0 b q φ q , c = Q-1 q=0
c q φ q . Let b, ĉ, G be the discrete Fourier transform of (kb q ) q , (kc q ) q ,

(K q -K q-1 ) q : b(ξ) = k Q-1 q=0 b q e -2πiqkξ , ĉ(ξ) = k Q-1 q=0 c q e -2πiqkξ , G(ξ) = q∈Z (K q -K q-1 )e -2πiqkξ .
Then, according to Parseval equality for Fourier series, we have:

T 0 b(t) 2 dt = k Q-1 q=0 b 2 q = 1 2k -1 2k | b(ξ)| 2 dξ, a(b, c) = 1 2k -1 2k G(ξ) b(ξ)ĉ(ξ)dξ. Let D n (x) = 1 + 2 n q=1
cos(qx) be the Dirichlet kernel. By summing by parts, and using the decrease of (K q ) q≥0 , we obtain:

Re G(ξ) = ∞ q=1 (K q -K q-1 )(cos(2πqkξ) -1) + lim q→∞ (K 0 -K q ) ≥ (K 1 -K 0 )(1 -1) + ∞ q=1 (K q+1 -2K q + K q-1 ) q q=0 (1 -cos(2π qkξ)) ≥ ∞ q=1 (K q+1 -2K q + K q-1 ) q + 1 2 - 1 2 D q (2πkξ) . • For x ∈] - π 4n , π 4n [ \ {0}, D n (x) ≤ - n q=1 q 2 cos π 4 = - √ 2 2 n(n + 1)(2n + 1) 6 ≤ - n 3 6
. By integrating twice, we get that D n (x) ≤ 2n + 1 -

n 3 12 x 2 . • For x ∈ ]-π, π[ \ ] -π 2n , π 2n [, we have | sin x 2 | ≥ | sin π 4n | = π 4n | sinc π 4n | ≥ π 4n | sinc π 4 | ≥ 2 3n ,
and then

|D n (x)| = sin(n + 1 2 )x sin x 2 ≤ 1 | sin x 2 | ≤ 3n 2 (sin, sinc are monotone on [0, π 4 ],
where sinc is defined by: sinc(x) = sin x x ).

• For

x ∈] - π 2n , π 2n [ \ ] -π 4n , π 4n
[, D n (x) has the same sign as -x. Hence

D n (x) ≤ 2n + 1 - n 3 12 ( π 4n ) 2 ≤ 2n + 1 - n 20 , q + 1 2 - 1 2 D q (2πkξ) ≥ min 1 24 q 3 (2πkξ) 2 , q 40 , q + 1 2 - 3q 4 for ξ ∈] -1 2k , 1 2k [. Then Re G(ξ) ≥ Tm≤qk≤2Tm (K q+1 -2K q + K q-1 ) min 1 24 q 3 (2πkξ) 2 , q 40 , Re G(ξ) ≥ T m k k 2 E min 1 24
T m k 3 (2πkξ) 2 , 1 40
T m k , Re G(ξ) ≥ T 2 m E min T 2 m ξ 2 , 1 40 .
Using Cauchy-Schwarz inequality, we also have:

b ∞ ≤ √ T b L 2 .
Hence,

1 4T -1 4T | b(ξ)| 2 dξ ≤ 1 2T b 2 ∞ ≤ 1 2 b 2 L 2 and 1 2 b 2 L 2 ≤ ]-1 2k , 1 2k [ \ ] 1 4T ,-1 4T [ | b(ξ)| 2 dξ
, and therefore:

a(b, b) ≥ ]-1 2k , 1 2k [ \ ]-1 4T , 1 4T [ Re G(ξ)| b(ξ)| 2 dξ ≥ ]-1 2k , 1 2k [ \ ]-1 4T , 1 4T [ 2α T | b(ξ)| 2 dξ ≥ α T b 2 L 2 .
Proposition 2. Assume that ( K(σ) q ) q∈Z satisfies the hypotheses of Lemma 4 for each σ ∈ {σ 1 , . . . , σ M }.

Let D = M i=1 |e i |.
Then a T,k , as defined in (3.18), is continuous and coercive:

∀p, c ∈ L 2 (0, T, H 1 (B)), |a T,k (p, c)| ≤ 2C p L 2 (0,T,H 1 (B)) c L 2 (0,T,H 1 (B)) , ∀p ∈ L 2 (0, T, H 1 (B)), a T,k (p, p) ≥ α T √ 1 + D 2 p 2 L 2 (0,T,H 1 (B))
.

Proof:

Continuity is a direct consequence of Proposition 1 and Lemma 4.

For c ∈ H 1 (B), |c(M ) -c(O 1 )| ≤ B ∂c ∂x (e) dx ≤ B dx B ∂c ∂x (e) 2 dx.
It follows that the Poincaré inequality:

p 2 H 1 (B) ≤ (D 2 + 1) B ∂p ∂x (e)
2 dx holds true.

Then, using Lemma 4 with a = a T,k , for p ∈ L 2 (0, T, H 1 (B)), we have that

a T,k (p, p) ≥ α T B T 0 ∂p ∂x (e) 2 dxdt ≥ α T T 0 1 √ D 2 + 1 p(., t) 2 H 1 (B) dt ≥ α T √ D 2 + 1 p 2 L 2 (0,T,H 1 (B)) .

Total discretization on the graph

Let h be a space step, h > 0.

For each edge e, let us denote h (e) the largest space step such that h (e) ≤ h and

S (e) = |e| h (e) is an integer. Let 0 = x (e) 0 , x (e) 1 = h (e) . . . , x (e)
S (e) = S (e) h (e) = |e| be the subdivision [0, |e|] with step h (e) . For s ∈ {0, . . . , S (e) }, let X (e) s designates the point of e with last local coordinate equal to x (e) s (so that if the local coordinate system has O i as origin then X (e) s = O i + R (e) ( 0, x (e) s )). Let P 1 h (B) be the subset of continuous functions ψ :

B → R such that ψ is affine on [X (e) r , X (e)
r+1 ] for each r ∈ {0, . . . , S (e) -1} (i.e., in local coordinates, ψ(O i + R (e) ( 0, x (e) )) is affine on [x (e) s , x

s+1 ]). Assume that a T,k is continuous coercive. Let us define p h,k ∈ P 0 k (0, T, P 1 h (B)) by:

∀ψ ∈ P 0 k (0, T, P 1 h (B)), a T,k (p h,k , ψ) = b T,h (ψ) (3.19)
where b T,h is an approximation of b T defined below. Let P h be the set of all X (e) s , with e an edge of the graph and s ∈ {0, . . . , S (e) }, with O 1 excluded. For X ∈ P h , let us define ψ X ∈ P 1 h (B) by ψ X (X) = 1 and ψ X (Y ) = 0 if Y ∈ P h and Y = X. We approximate the integral on B in b T by the trapezoidal rule. So, we define b T,h by:

b T,h (ψ X φ q ) = (F (X, t q+1 ) -F (X, t q )) B ψ X (x)dx + N =1 ψ X (O )(Ψ (O , t q+1 ) -Ψ (O , t q )),
for every X ∈ P h and q ∈ {0, . . . , Q -1}. Note that if F = 0, then b T,h = b T , and that for any ψ ∈ P 1 h (B) and f ∈ H r dc (B), r ∈ {1, 2}, we have:

B f ψdx - X∈P h f (X)ψ(X) B ψ X dx ≤ h r B ∂ r (f ψ) (∂x (e) ) r dx.
In the right handside, B is to be understood as a sum of integrals on each step of the space subdivision: over these intervals, ψ is smooth and the r th -derivative of f ψ is in L 1 . Furthermore:

B ∂(f ψ) ∂x (e) dx ≤ h( f H 1 dc (B) ψ L 2 (B) + f L 2 (B) ψ H 1 dc (B) ). B ∂ 2 (f ψ) (∂x (e) ) 2 dx ≤ h 2 ( f H 2 dc (B) ψ L 2 (B) + 2 f H 1 dc (B) ψ H 1 dc (B) + 0) if r = 2.
Hence, if F ∈ H 1 (0, T, H r dc (B)), then, by using Cauchy-Schwarz inequality:

|(b T,h -b T )(ψφ q )| ≤ 3h r √ k F H 1 ([tq,tq+1],H r dc (B)) ψ H 1 (B) ,
and therefore, for ψ ∈ P 0 k (0, T, P 1 h (B)) and F ∈ H 1 (0, T, H r dc (B)), r ∈ {1, 2}, we have:

|(b T,h -b T )(ψ)| ≤ 3h r F H 1 (0,T,H r dc (B)) ψ L 2 (0,T,H 1 (B)) . (3.20) Now, rewriting p h,k = X∈P h ,0≤q≤Q-1 P X,q ψ X φ q ,
and introducing the notations

K q,X,Y = B K(σ) q ∂ψ X ∂x (e) ∂ψ Y ∂x (e) dx, R X,q = b T,h   ψ X q-1 q=0 φ q   = b T,h (ψ X 1 [t0,tq] ),
the totally discretized problem (3.19) is equivalent to solve for the coordinates P q,X of p h,k the following system:

∀r ∈ {0, ..., Q-1}, X ∈ P h : Y ∈P h   K 0,X,Y P Y,r + r-1 q=0 (K r-q,X,Y -K r-q-1,X,Y )P Y,q   = R X,r+1 -R X,r .
(3.21) Since R 0,X = 0, by summation over r this can be reformulated as:

∀q ∈ {0, ..., Q -1}, X ∈ P h : q q=0 Y ∈P h K q-q,X,Y P Y,q = R X,q+1 . Note that R X,q =          Ψ (t q ) + 1 2 1≤j≤M,O ∈ej h (ej ) F (O , t q ) if X = O , h (e) F (X (e) s , t q+1 ) if X = X (e)
s , s ∈ {1, . . . , S (e) -1}, and

K q,X,Y =                              1≤j≤M,O ∈ej K(σj) q h (ej ) if X = Y = O , 2 K(σj) q h (ej ) if X = Y = X (ej ) s , s ∈ {1, . . . , S (ej ) -1}, - K(σj) q h (ej ) if {X, Y } = {X (ej ) s , X (ej )
s+1 }, s ∈ {0, . . . , S (ej ) -1}, 0 otherwise.

Hence, we have the same scheme as described in the introduction, Section 1.2.

Discretization in the cross-section

If we want to use Theorem 1 and 2 of Section 4, we must determine an approximation of K (σ) q such that θ(k) is small. In Section 5, we test the scheme on the graph with the exact value of the kernel when known, and with the three different approximations for the kernel, described below: (3.26),(3.29),(3.30).

Exact kernels

When the cross-section σ is a rectangle, an equilateral triangle or a disc, we express in[8] K (σ) as a computable generalized Dirichlet series:

K (σ) (t) = ∞ k=-∞ 4 (2k + 1) 2 π 2 e -π 2 (2k+1) 2 t if σ =]0, 1[, (3.22)
K (σ) (t) = abK (]0,1[) (a -2 t)K (]0,1[) (b -2 t) if σ =]0, a[×]0, b[, (3.23)
K (σ) (t) = 4π +∞ i=1 1 µ 2 i e -µ 2 i t if σ = x ∈ R 2 ; x 2 < 1 , (3.24)
K (σ) (t) = 3 √ 3 2π 2 +∞ k=1 1 k 2 exp - 16π 2 3 k 2 t if σ = T (3.25)
where µ i denotes the i-th zero of the zero-th Bessel function J 0 and T denotes the triangle with vertices (0, 0), (0, 1), (

√ 1 2 ,
.

If these exact expressions are used, then θ(k) = 0. For the square and the triangle cross section, we obtain series expansions for the kernels that converge rapidly for small times, by using the Poisson summation formula. Furthermore, one can show that (K (σ) q ) satisfies the hypotheses of Proposition 2. Hence, with this choice, the scheme on the graph is unconditionally stable.

Numerical approximation of the kernels

For an arbitrary smooth domain σ, let us consider a quasi-uniform triangulation of σ and denote by H the diameter of its largest triangle. Let S (σ) H ⊂ L 2 (σ) be the space of P r-1 finite elements where r ≥ 2. Let T (σ)

H : L 2 (σ) → S (σ)
H be the operator discretizing the inverse of the Dirichlet Laplacian, using Nitsche method for the boundary conditions. Then, (T

(σ) H -τ I) : S (σ) H → S (σ)
H is invertible whenever τ ≥ 0. Consider the following numerical scheme which use the implicit Euler method as time integrator:

           V 0 H,k = W (σ) , T (σ) H V q+1 H,k -V q H,k k = V q+1 H,k , K (σ) q = σ V q H,k dx, (3.26)
where W (σ) is an approximation of 1. Note that there is a special method of resolution which is presented in Section 3.5.3.

In the present paper, we test the case where W (σ) is the interpolant of (I -k∆) -1 1 and can be computed explicitly.

If W = 1, then we show in [START_REF] Canon | Numerical solution of the viscous flows in a network of thin tubes: equations on the cross-section[END_REF] that ( K(σ) q ) satisfies the hypotheses of Proposition 2. Hence, with this choice, the scheme on the graph is unconditionally stable. Using techniques in Thomée's book [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF] we also prove that when h = k γ , γ > 0, and k small enough, the scheme is convergent with:

θ(k) ≤ C max{H 2 , k} 1 3 .
(3.27)

The numerical experiments in [START_REF] Canon | Numerical solution of the viscous flows in a network of thin tubes: equations on the cross-section[END_REF] indicate that the power 1/3 may be improved and replaced by 1/2 and the use of L-stable higher-order methods do not improve the accuracy, unless an approximation for small times is used as in the next paragraph.

Corrected approximation

When the numerical approximation is used, most of the error is done for small time. Hence, if an asymptotic expansion near t = 0 + is known, one can use it to obtain a more accurate approximation of the kernel. When σ is smooth and simply connected, using the boundary layer theory presented in Gie, Jung & Temam [START_REF] Gung-Min Gie | Recent progresses in boundary layer theory[END_REF], we obtain in [START_REF] Canon | Numerical solution of the viscous flows in a network of thin tubes: equations on the cross-section[END_REF] an asymptotic expansion of K.

Let us denote γ : [0, |∂σ|] → ∂σ, a parametrization by the arclength of ∂σ (i.e. γ = 1), κ(s) = det(γ (s), γ (s)) the curvature. Then:

K (σ) (t) = |σ| - 2 √ π |∂σ|t 1 2 + πt + 1 6 √ π t 3 2 |∂σ| 0 κ(s) 2 ds + O t→0 + (t 2 ). (3.28)
Let us consider a modified version of scheme (3.26), and take the BDF2 time integrator (see [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]) and P 1 -elements in space.

                   V 0 H,k = 1, T (σ) H V 1 H,k -V 0 H,k k = V 1 H,k , T (σ) H 3V q+2 H,k -4V q+1 H,k + V q H,k 3k = V q+2 H,k if q ≥ 0 K(σ) q,num = σ V q H,k + V q+1 H,k 2 dx.
(3.29)

Let us use the asymptotic expansion (3.28) to correct the previous approximation for small times.

From (3.28),

τ 0 K (σ) dt = χ(t) + O t→0 + (t 5
2 ), with

χ(t) = |σ|t - 4 3 √ π |∂σ|t 3 2 + π 2 t 2 .
We take K(σ)

q,cor to be such that:

   K (σ) q,cor = K (σ) q,num if Q 1 3 ≤ q ≤ Q -1, K (σ) q,cor -K (σ) q-1,cor k = χ((q + 1)k) -2χ(qk) + χ((q -1)k) k 2 if 1 ≤ q ≤ Q 1 3 .
(3.30)

In this case, the hypotheses of Proposition 2 are not automatically satisfied, but we can refine (3.27) and get, when

H = k γ , γ ≥ 1 2
and k small enough:

θ(k) ≤ C max{k, H} 2 3 . (3.31)
Numerical experiments suggest that the power 2/3 may be replaced by 1, or more if all the terms from (3.28) are used.

Remarks on the solving of the scheme

Let us describe three methods to solve the scheme (3.21) on the graph. The first two methods can be used for any type of kernel approximation but the last one will only work for (3.26).

Let S be the number of elements of P h .

Using the block triangular structure

The first computational strategy uses the block-triangular structure of the matrix with respect to the time variable and solves the corresponding system of equations by the method of forward substitution. In this case, only the sparse matrix (K 0,X,Y ) X,Y ∈P h needs to be inverted. It is tridiagonal, except for the junctions. Hence, its Choleski factorisation takes O(N S + N 3 ) flops.

The Choleski factor has O(N S + N 2 ) nonzero elements. The computation of (P X,q ) X∈P h at the time step number q, when values (P X,0 ) X∈P h , . . ., (P X,q-1 ) X∈P h are known, costs O(q(N S + N 2 )) flops.

Hence, the total computational cost is O N S + N 3 + Q 2 (N S + N 2 ) .

Using the Fast Fourier Transform

Let us describe the second strategy. Denote (A q ) 0≤q≤Q-1

(B q ) 0≤q≤Q-1 =   q q=0 A q B q-q   0≤q≤Q-1 .
This discrete convolution can be computed by using the Discrete Fourier Transform with only O(Q log(Q)) multiplications. First, let us solve the equation:

(K Z) q,X,Y = δ X,Y if q = 0, 0 if q = 0.
The solution Z can be computed by log 2 Q Newton's iterations:

Z = Φ log 2 Q (A), A q = ((K 0,X,Y ) X,Y ∈P h ) -1 if q = 0 0 if q = 0 , Φ(B) = B (2I -B K B).
Then P X,q is given by: ∀q ∈ {0, ..., Q -1}, X ∈ P h :

P X,q = Y ∈P h q q=0 Z q,X,Y R q+1,Y .
The computational cost is here O(Q log(Q) 2 S 3 ).

Remark: When there is no forcing along the pipe (F = 0), the second method is more efficient for large Q since one can take one space step per pipe without loss of precision, which implies S = N .

Coupling with the cross-section

The following formulation of the scheme is of interest when we want to recover the velocity inside the tubes. Besides, contrary to the two previous methods, its complexity scales linearly with the number of time steps. For x ∈ B, x ∈ σ(x), let U (x, x, t) such that:

       ∂U ∂t (x, x, t) -∆ xU (x, x, t) = - ∂P ∂x e (x), x ∈ σ, t > 0, U (x, t) = 0, x ∈ ∂σ, t > 0, U (x, 0) = 0, x ∈ σ.
Then [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure. I. The case without boundary-layer-in-time[END_REF][START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF], at some distance from the junctions:

u ε (O ij + R (ej ) (εx, x d ), tε 2 ) R (ej ) ( 0, U (x, O ij + R (ej ) ( 0, x d ), t).
In the sequel, U

(ej ) s+ 1 2 ,q designates an approximation of U •, X (ej ) s+ 1 2 , t q+1/2 , with X (ej ) s+ 1 2 = X (ej ) s+1 + X (ej ) s 2 .
Let K(σ) q as defined in (3.26). Then, the finite difference scheme is equivalent to the continuity condition (1.10) and Dirichlet condition for the pressure (1.11) coupled with following equations:

U (ej ) s+ 1 2 ,q ∈ S (σj ) H ⊂ L 2 (σ j ) if 1 ≤ j ≤ M, 0 ≤ s < S (ej ) , 0 ≤ q < Q, U (ej ) s,0 = 0 if 1 ≤ j ≤ M, 0 ≤ s < S (ej ) , U (ej ) s+ 1 2 ,q+1 = (I -k(T (σj ) H ) -1 ) -1 U (ej ) s+ 1 2 ,q+1 -k P (ej ) s+1,q+1 -P (ej ) s,q+1 h (ej ) W (σj ) if      1 ≤ j ≤ M, 0 ≤ s < S (ej ) , 0 ≤ q < Q, f (ej )
s,q+1 = σj U (ej ) s+ 1 2 ,q+1 -U (ej ) s-1 2 ,q+1 h (ej ) dx if 1 ≤ j ≤ M, 1 < s < S (ej ) , 0 ≤ q < Q, Ψ ,q+1 = 1≤j≤M,O ∈ej - h (ej ) 2 f (ej ) s,q+1 + α ,j σj U (ej ) s+ 1 2 α ,j ,q+1 dx if 1 ≤ ≤ N, 0 ≤ q < Q,
where, in the sum, s is such that X (ej ) s = O . This formulation can be solved in the following way. When the U (ej ) s,q are known for a given q ∈ {0, 1 . . . , Q -2}, P (ej ) s,q+1 is the unique solution to (1.10)-(1.11) coupled with the system:

k K(σ) 0 P (ej ) s+1,q+1 -2P (ej ) s,q+1 + P (ej ) s-1,q+1 (h (ej ) ) 2 = -f (ej ) s,q+1 - σj (I -k(T (σj ) H ) -1 ) -1 U (ej ) s+ 1 2 ,q -U (ej ) s-1 2 ,q h (ej ) dx if 1 ≤ j ≤ M, 0 ≤ s < S (ej ) , 0 ≤ q < Q, Ψ ,q+1 = 1≤j≤M,O ∈ej   - h (ej ) 2 f (ej ) s,q+1 -k K(σj) 0 P (ej ) s+α ,j ,q+1 -P (ej ) s,q+1 h (ej ) +α ,j σj (I -k(T (σj ) H ) -1 ) -1 U (ej ) s+ 1 2 α ,j ,q dx , if 1 ≤ ≤ N, 0 ≤ q < Q,
where in the sum, s is such that X (ej ) s = O . Once (P (ej ) s,q+1 ) are known, (U (ej ) s+ 1

Convergence of the scheme

The following theorem is the main result of this paper. It proves that the proposed scheme on the graph is convergent with order one in space and time, provided the discretization in the crosssection is well-chosen: either θ(k) = 0 (exact kernels) or the order is limited by the accuracy of the approximations of the kernels discussed above (3.27), (3.31). With additional assumptions, we improve the order of convergence in Theorem 2.

Theorem 1. Let θ(k) = max 1≤j≤M |K (σj ) 0 - K(σj) 0 | + Q-1 q=1 |K (σj ) q -K (σj ) q-1 -K(σj) q + K(σj) q-1 |.
Let p h,k ∈ L 2 (0, T, H 1 (B)) be the function such that:

• p h,k is affine on [X (ej ) s , X (ej )
s+1 ] × {t} for t ∈ [0, T ], 1 ≤ j ≤ M , s ∈ {0, . . . , S (ej ) -1}.

• p h,k (X (ej ) s , t) = P (ej ) s,q for t ∈]t q , t q+1 [, 1 ≤ j ≤ M , s ∈ {0, . . . , S (ej ) }, q ∈ {0, . . . , Q -1}.

Assume that

θ(k) → 0 as k → 0, F ∈ H 1 00 (0, T, H 1 dc (B)), Ψ 1 , . . . , Ψ N ∈ H 1 00 (0, T ).
Then:

p h,k → P in L 2 (0, T, H 1 (B)) when (h, k) → (0, 0).
Furthermore, there exists C 1 , C 2 such that if:

θ(k) < C 2 , F ∈ H 2 00 (0, T, H 2 dc (B)), Ψ 1 , . . . , Ψ N ∈ H 2 00 (0, T ), then P ∈ H 1 (0, T, H 2 dc (B)
) and:

p h,k -P L 2 (0,T,H 1 (B)) ≤C 1 h P L 2 (0,T,H 2 dc (B)) + h C 2 -θ(k) F H 1 (0,T,H 1 dc (B)) +k P H 1 (0,T,H 1 (B)) + P L 2 (0,T,H 1 (B)) θ(k) C 2 -θ(k)
.

Proof: We use the following lemma.

Lemma 5. There exists C > 0 such that, for every u ∈ H 1 (0, T, H 1 (B) ∩ H 2 dc (B)):

inf ψ∈P 0 k (0,T,P 1 h (B)) ψ -u L 2 (0,T,H 1 (B)) ≤ C h u L 2 (0,T,H 2 dc (B)) + k u H 1 (0,T,H 1 (B))
Proof of the Lemma: If u ∈ H 1 (B)∩H 2 dc (B) and P h is the continuous piecewise linear interpolant in P 1 h (B) at the points of the subdivision, then:

u -P h u H 1 (B) ≤ Ch u H 2 dc (B) , P h u H 1 (B) ≤ C u H 1 (B) .
Now, let v ∈ H 1 (0, T ) and let Q k be the L 2 -projector on P 0 k (0, T ). Then:

v -Q k v L 2 (B) ≤ Ck v H 1 (B)
and

p -P h Q k p L 2 (0,T,H 1 (B)) ≤ p -Q k p L 2 (0,T,H 1 (B)) + Q k p -P h Q k p L 2 (0,T,H 1 (B)) ≤ Ck p H 1 (0,T,H 1 (B)) + Ch Q k p L 2 (0,T,H 2 dc (B)) ≤ Ck p H 1 (0,T,H 1 (B)) + Ch p L 2 (0,T,H 2 dc (B))
. This proves the lemma. Now, we prove the convergence theorem.

Let α be the coercivity constant of a T . From Proposition 1 in Section 3.2, we know that:

a T,k -a T ≤ θ(k).
Thus, for k small enough, a T,k is (α -θ(k))-coercive. Let α(k) be the constant of coercivity of a T,k . Let p h,k ∈ P 0 k (0, T, P 1 h (B)) be the unique solution of:

∀ψ ∈ P 0 k (0, T, P 1 h (B)), a T (p h,k , ψ) = b T (ψ).
According to Céa's lemma, we have:

P -p h,k L 2 (0,T,H 1 (B)) ≤ a T α inf ψ∈P 0 k (0,T,P 1 h (B)) P -ψ L 2 (0,T,H 1 (B)) . ( 4
a T,k (p h,k -p h,k , ψ) = a T (p h,k , ψ) -a T,k (p h,k , ψ) + b T,h (ψ) -b T (ψ).
By taking ψ = p h,k -p h,k and using coercivity on the left-hand side and continuity on the right-hand side, we get:

α(k) p h,k -p h,k L 2 (0,T,H 1 (B)) ≤ θ(k) p h,k + b T,h -b T L 2 (0,T,H 1 (B)) .
Hence, we have:

P -p h,k L 2 (0,T,H 1 (B)) ≤ a T α inf ψ∈P 0 k (0,T,P 1 h (B)) P -ψ L 2 (0,T,H 1 (B)) + θ(k) α(k) p h,k L 2 (0,T,H 1 (B)) + 1 α(k) b T,h -b T L 2 (0,T,H 1 (B)) . (4.33)
Let us now show that lim (h,k)→(0,0)

P -p h,k L 2 (0,T,H 1 (B)) = 0. As H 1 (0, T, H 1 (B) ∩ H 2 dc ((B)) is dense in L 2 (0, T, H 1 ((B)), Lemma 5 implies that lim (h,k)→(0,0) inf ψ∈P 0 k (0,T,P 1 h (B)) P -ψ L 2 (0,T,H 1 (B)) = 0.
From equation (3.20), lim h→0 b T,h -b T L 2 (0,T,H 1 (B)) = 0, and, as α(k) ≥ α -θ(k) implies that lim inf α(k) > 0, we conclude that lim (h,k)→(0,0)

P -p h,k L 2 (0,T,H 1 (B)) = 0.
This is the first part of Theorem 1. Now, if we assume in addition that F ∈ H 2 00 (0, T, H 2 dc (B)), Ψ 1 , . . . , Ψ N ∈ H 2 00 (0, T ), then, from the regularity resulsts of Section 2.3, we have that P ∈ H 1 00 (0, T, H 2 dc (B)). Hence, using Lemma 5:

inf ψ∈P 0 k (0,T,P 1 h (B)) P -ψ L 2 (0,T,H 1 (B)) ≤ Ch P L 2 (0,T,H 2 dc (B)) + Ck P H 1 (0,T,H 1 (B)) . Besides, if p h,k = 0, α p h,k L 2 (0,T,H 1 (B)) ≤ 1 p h,k L 2 (0,T,H 1 (B)) a T (P, p h,k ) ≤ a T P L 2 (0,T,H 1 (B))
and from equation (3.20), we get:

b T,h -b T L 2 (0,T,H 1 (B)) ≤ 3h F H 1 (0,T,H 1 dc (B)) . Also, α(k) ≥ α -θ(k).
Combining all these inequalities into (4.33), we get the second part of Theorem 1.

Let us introduce the following lemma which is used to prove Theorem 2.

Lemma 6. Let p : [0, T ] → R a C 2 function. Let p = Q-1 q=0 p(t q+1/2 )φ q , where φ q = 1 [tq,tq+1[ . Assume that K : [0, T ] → R is a continuous function, continuously differentiable on ]0, T ], such that K ≥ 0, K ≤ 0 and that K (t) = O t→0 + (t -1/2 ). Let Lp(t) = t 0 K(t -s)p(s)ds.
Then, for any φ ∈ P k (0, T ):

T 0 (L(p -p)) (t)φ(t)dt ≤ Ck 2 ( p ∞ log(T /k) + p ∞ ) φ L 2 .
If furthermore p (0) = 0, then:

T 0 (L(p -p)) (t)φ(t)dt ≤ Ck 2 p ∞ φ L 2 .
Proof:

Let µ(t) = t k/2 (p -p)(τ )dτ , so that µ is C 2 on [0, T ]\{t 0 , t 1 , . . . , t Q } and µ (t q+1/2 ) = 0 for q ∈ {0, . . . , Q -1}.
Besides, for any t ∈]0, T -k[\{t 0 , t 1 , . . . , t Q }, there exists q ∈ {0, 1, . . . , Q} such that t q ∈]t, t + k[. Then:

|µ (t + k) -µ (t)| = |p (t + k) -p (t)| ≤ k p ∞ .
Integrating from s = t q-1/2 to s = t, we get:

|µ (t + k) -µ (t)| ≤ k 2 p ∞ . Integrating from s = t q-1 to s = t, using |µ(t q ) -µ(t q -k)| = | tq tq-1 p(t) -p(t q+1/2 )dt| ≤ k 3 p ∞ ,
we get:

|µ(t + k) -µ(t)| ≤ k 3 p ∞ + |µ(t q ) -µ(t q -k)| ≤ 2k 3 p ∞ .
Let us now introduce µ in the integral to be bounded:

T 0 (L(p -p)) (t)φ q (t)dt = tq+1 tq K(0)(p -p)(t) + t 0 K (τ )(p -p)(t -τ )dτ dt = K(0)(µ(t q+1 ) -µ(t q )) + tq+1 tq tq 0 K (τ )(p -p)(t -τ )dτ dt + tq+1 tq t tq K (τ )(p -p)(t -τ )dτ dt = K(0)(µ(t q+1 ) -µ(t q )) + tq 0 K (τ )(µ(t q+1 -τ ) -µ(t q -τ ))dτ + tq+1 tq t tq K (τ )(p -p)(t -τ )dτ dt Since -Ct -1 2 ≤ K (t) ≤ 0 and thus T 0 |K (t)|dt ≤ K(0)
, and since p -p ≤ k p , we get:

T 0 (L(p -p)) (t)φ q (t)dt ≤ K(0)k 3 p ∞ + 2K(0)k 3 p ∞ + k p ∞ tq+1 tq t tq Ct -1 2 dτ dt ≤ 3K(0)k 3 p ∞ + 2Ck 5 2 (q + 1) -1 2 p ∞ .
From this, we conclude that, for φ = Q-1 q=0 a q φ q ∈ P 0 k (0, T ):

T 0 (L(p -p)) (t)φ(t)dt ≤ 2 Q-1 q=0 |a q |C(q + 1) -1 2 k 5 2 p ∞ + 3 Q-1 q=0 |a q |k 3 K(0) p ∞ .
With Cauchy-Schwarz inequality this yields:

T 0 (L(p -p)) (t)φ(t)dt ≤ 2C φ L 2 √ k k 5 2 p ∞ Q-1 q=0 (q + 1) -1 1/2 + 3 φ L 2 √ k T k K(0)k 3 p ∞ , T 0 (L(p -p)) (t)φ(t)dt ≤ 2Ck 2 p ∞ (log(T /k)) 1/2 + 3 √ T K(0)k 2 p ∞ φ L 2 .
If furthermore p (0) = 0, then, for t

∈ [0, k], |(p -p)(t)| ≤ k 2 p ∞ .
We conclude that, in this case:

T 0 (L(p -p)) (t)φ q (t)dt ≤ 3K(0)k 3 p ∞ + k 2 p ∞ tq+1 tq t tq |K (τ )|dτ dt ≤ 3K(0)k 3 p ∞ + K(0)k 3 p ∞ Hence: T 0 (L(p -p)) (t)φ(t)dt ≤ 4K(0) √ T k 2 φ L 2 .
This ends the proof of Lemma 6. Now, we prove Theorem 2 below.

Theorem 2. Let θ(k) and p h,k be defined as in Theorem 1. Let ph,k ∈ L 2 (0, T, H 1 (B)) be the function such that:

• ph,k is affine on [X (ej ) s , X (ej )
s+1 ] × {t} for t ∈ [0, T ], 1 ≤ j ≤ M , s ∈ {0, . . . , S (ej ) -1}. • ph,k (X (ej ) s , t) = P X (ej ) s , t q+1/2 for t ∈]t q , t q+1 [, 1 ≤ j ≤ M , s ∈ {0, . . . , S (ej ) }, q ∈ {0, . . . , Q -1}.
Assume that the kernels K (σ) satisfy the hypotheses of Lemma 6; assume that P is a C 4 function on each edge of the graph; assume that θ(k) → 0 as k → 0.

Let β(k) be defined by

β(k) = k 2 if ∂P ∂t (•, 0) is constant, β(k) = k 2 log(T /k) otherwise.
Then, there exists two positive constants C 1 , C 2 such that for θ(k) < C 2 :

p h,k -ph,k L 2 (0,T,H 1 (B)) ≤ C 1 C 2 -θ(k) (β(k) + h 2 + θ(k))m(P ),
where

m(P ) = ∂ 2 P ∂x∂t ∞ + ∂ 3 P ∂x∂t 2 ∞ + ∂ 4 P ∂x 4 ∞ .
Remark The fact that K (σ) satisfies the assumptions of Lemma 6 for any smooth σ and also for all the examples given in Section 3.4.1 is proven in [START_REF] Canon | Numerical solution of the viscous flows in a network of thin tubes: equations on the cross-section[END_REF].

Proof Let ph (t, x) = X∈P h P (X, t)ψ X . Since from equation (1.2) ∂F ∂t (x, t) = -K (σ) (0) ∂ 2 P (∂x (e) ) 2 (x, t) - t 0 (K (σ) ) (t -τ ) ∂ 2 P (∂x (e) ) 2 (x, τ )dτ,
we have that F ∈ H 1 00 (0, T, H 2 dc (B)) and that F H 1 (0,T,H 2 dc (B)) ≤ Cm(P ). Similarly, Ψ 1 , . . . , Ψ n ∈ H 1 00 (0, T ). Let φ ∈ P 0 k (0, T ) and ψ ∈ P 1 h (B). From the previous lemma, it follows that:

|a T (p h,k -ph , φψ)| ≤ C φ H 1 ψ L 2 β(k)m(P )
where a T is defined in (2.15). Let us compute now a T (P -ph , φψ). As

[X (e) s ,X (e) s+1 ]
∂ 2 L (σ) (P -ph ) ∂t∂x (e) ∂ψ ∂x (e) dx =

∂(L σ P ) ∂t (X (e) s+1 , t) - ∂(L σ P ) ∂t (X (e) s , t) - ∂(L (σ) ph ) ∂t (X (e) s+1 , t) + ∂(L (σ) ph ) ∂t (X (e) s , t) ∂ψ ∂x (e) (X (e) s+ 1 2
)h (e) , = 0 we get that a T (P -ph , φψ) = 0, so that

|a T (P -ph,k , φψ)| ≤ Cβ(k) φ H 1 ψ L 2 m(P ).
Hence, for any ψ ∈ P 0 k (0, T, P 1 h (B)):

|a T (P -ph,k , ψ)| ≤ Cβ(k) ψ L 2 (0,T,H 1 (B)) m(P ).
We have:

a T,k (p h,k -ph,k , ψ) = a T,k (p h,k , ψ) + (a T -a T,k )(p h,k , ψ) + a T (P -ph,k , ψ) -a T (P, ψ) = (b T,h -b T )(ψ) + (a T -a T,k )(p h,k , ψ) + a T (P -ph,k , ψ)
where a T,k is defined in (3.18). Taking ψ = p h,k -ph,k , using the coercivity of a T,k on the left handside, using continuity of a T , a T -a T,k and the previous inequality on the right handside, and dividing by the norm of p h,k -ph,k , we get:

α(k) p h,k -ph,k 2-1 L 2 (0,T,H 1 (B)) ≤ Ch 2 F H 1 (0,T,H 2 dc (B)) + θ(k) ph,k L 2 (0,T,H 1 (B)) + Cβ(k)m(P ). With F H 1 (0,T,H 2 dc (B)) ≤ Cm(P ), p h,k H 1 (0,T,H 2 dc (B)) ≤ Cm(P ) and α(k) ≥ a T -θ(k), the proof is complete. Corollary 1 If P is of class C 4 , then: ph,k -p h,k L ∞ (0,T,H 1 (B)) ≤ Cm(P )(θ(k)k -1 2 + h 2 k -1 2 + k 3 2 log(T /k)).
If P is of class C 5 and ∂ r P ∂t r (•, 0) = 0 for r ∈ {0, . . . , 5}, then:

ph,k -p h,k L ∞ (0,T,H 1 (B)) ≤ Cm( ∂P ∂t )(θ(k) + h 2 + k 2 )
Proof: The first part is obtained by using that, for φ ∈ P 0 (0, T ):

√ k φ L ∞ (0,T ) ≤ φ L 2 (0,T ) .
Hence:

√ k p h,k -ph,k L ∞ (0,T,H 1 (B)) ≤ p h,k -ph,k L 2 (0,T,H 1 (B)) ,
and Theorem 2 yields the result.

For the second part of the corollary, let us extend P , p h,k , ph,k by zero for negative times. Let us

denote ∆ k φ(t) = φ(t) -φ(t -k) k . Then ∆ k P is a solution to (1.2) with F, Ψ 1 , . . . , Ψ N replaced by ∆ k F, ∆ k Ψ 1 , . . . , ∆ k Ψ N .
Applying Theorem 2 to ∆ k P over [0, T ], we get for sufficiently small k:

∆ k p h,k -∆ k ph,k L 2 (0,T,H 1 (B)) ≤ Cm(∆ k P )(θ(k) + h 2 + k 2 ). But m(∆ k P ) ≤ m( ∂P ∂t
), and for φ ∈ P 0 (0, T ) extended by zero,

φ L ∞ (0,T ) ≤ √ T ∆ k φ L 2 (0,T ) . Hence: p h,k -ph,k L ∞ (0,T,H 1 (B)) ≤ Cm( ∂P ∂t )(θ(k) + h 2 + k 2 ).
This completes the proof of the corollary.

22

In Section 5.1, we test the numerical accuracy of the scheme (3.4.2) with the variants in the cross-sections. In Section 5.2, we then compare the asymptotic model (1.2) with the initial 3D-Navier-Stokes equations.

Let us assume that an error E converges to 0 when a parameter ε tends to 0 and that E ∼ ε→0 Cε α for some constants C and α. The parameter α is the order of convergence of E with respect to ε. Let us consider a strictly decreasing sequence (ε n ) which converges to 0 and let E n be the corresponding sequence of errors. Then,

α = lim n→+∞ log( En-1 En ) log( εn-1 εn )
.

Hence, in all the tables below, when we report an error (E n ) n and a varying parameter (ε n ) n , we also report log( En-1 En ) log( εn-1 εn )

on a column on the right to estimate the order of convergence.

Numerical accuracy

To evaluate the numerical accuracy of the scheme (3.4.2), we first build an exact test case solution to the asymptotic system (1.2).

Test case on a single tube

Let d ∈ {2, 3}. We consider the case of a single tube (M = 1) of length 1 with two extremities

O 1 = ( 0, 0), O 2 = ( 0, 1) (N 1 = N = 2). Let the cross-section of the tube be σ = - 1 2 , 1 2 if d = 2, and σ = {x ∈ R 2 ; x 2 < 1, } if d = 3.
Let us take P ((

where β ∈ {0, 1}. When β = 1, P and all its time derivatives are zero when t → 0. Then, the flow at the left extremity O 1 of the pipe is given by:

Ψ 1 (t) = - t 0 K (σ) (s)(1 -(t -s)) exp - β t -s ds.
At the right extremity O 2 of the pipe, it is given by:

Ψ 2 (t) = t 0 K (σ) (s)(1 -(t -s)) exp 1 - β t -s ds.
The force applied along the pipe is:

F (( 0, x (e)
3 ), t) = -

t 0 K (σ) (s)(1 -(t -s)) 2 exp (1 -(t -s))x (e) 3 - β t -s ds.
The theoretical data Ψ 1 , Ψ 2 , F can be numerically computed with a high accuracy by using the change of variable τ = exp(1 -1 s )t to remove the singularity at s = 0. The integrand is then smooth and the Clenshaw-Curtis integration method yields accurate results. For the sake of completeness, we chose to evaluate the numerical accuracy of the scheme with an approximate kernel ( K(σ) q ) as proposed in Section 3.5.3, in the case d = 2, while for the case d = 3 we consider a tube with a disc cross-section and numerical approximation in the three following way: with the exact kernel given in Section 1.3, with a numerical approximation based on the discretization of equation (1.7), and last with this numerical approximation with a correction that uses the asympotic expansion (3.28).

Table 1: β = 1, p(x, t) = e (1-t)x-1/(t+ε) . In the last three lines, the order is computed with respect to H. For d = 2, we test the scheme with the an approximate kernel ( K(σ) q ) defined as proposed in Section 3.5.3.

Let W (x) = 1 - sinh(( 1 2 + x)/ √ k) + sinh(( 1 2 -x)/ √ k) sinh(1/ √ k)
be the solution of:

U -kU = 1 in σ = ]-1/2, 1/2[ , U = 0 on ∂σ = {-1/2, 1/2} .
Let H be a space step for the cross-section such that 1/H is an integer. We consider the following difference scheme in the cross-section:

             V 0,r = W - 1 2 + rH if r ∈ {1, . . . , 1/H -1}, V q,r -V q+1 k = V q,r-1 -2V q,r + V q,r-1 H 2 if q ∈ {1, . . . , Q}, r ∈ {1, . . . , 1/H -1}, V q,0 = V q,1/H = 0 if q ∈ {1, . . . , Q}.
Then we take K(σ)

q = H 1/H-1 r=1 V q,r .
We compute the maximum norm of the error of the numerical solution with respect to the exact solution. The results of the numerical test for different step sizes are listed in Table 1 when β = 1 and in Table 2 when β = 0. We observe that the numerical order for the time step is 1 when β = 0 and 1/2 when β = 1. Corollary 1 can be applied in both cases, but the accuracy is limited by the precision of the approximation of the kernel. When β = 1, we observe that the order of the error with respect to the cross-section space step H is 2. 4 (top: β = 0) and Table 3 (bottom: β = 1). On each graph, only one parameter varies, the two others are set by default to h • p = 0 and v(M, t) colinear to e 1 at the beginning O 1 of e 1 .

= 2 -10 , k = 0.1 • 2 -14 , H = π2 -10 . O 4 (1, 0, 0) e 1 O 1 (0, 0, 0) e 2 O 2 (1 + cos π 4 , sin π 4 , 0) e 3 O 3 (1 + cos π 4 , -sin π 4 , 0) O 1 (0, 0, 0) e 1 O 4 (1, 0, 0) e 4 O 5 (1 + cos π 4 , sin π 4 , 0) e 2 O 2 (1 + cos π 4 , 3 2 , 0) e 6 O 6 (1 + cos π 4 , -sin π 4 , 0) e 3 O 3 (1 + cos π 4 , - 3
•

v(M, t) = - e 2 |e 2 | v 0 sin ε -2 40t 1 -4ε -2 O 2 M 2 at the beginning O 2 of e 2 . The flux through this tube is then 1 8 πv 0 ε 2 sin ε -2 40t ; • v(M, t) = 2 e 3 |e 3 | v 0 sin ε -2 40t 1 -4ε -2 O 3 M 2 at the beginning O 3 of e 3 .
The flux through this tube is then π 4 v 0 ε 2 sin ε -2 40t ;

• v = 0 on the rest of the boundary.

Using the mass conservation law, we conclude that the flux through this first tube is: 3 8 πv 0 ε 2 sin(ε -2 40t).

Results

Simulations with COMSOL software run over the interval [0, 0.0875ε 2 ] for ε ∈ {0.2, 0.1, 0.05, 0.025}. Let q ε be the continuous in space function, affine over each edge, such that q ε (O ) = p ε (O ) for ∈ {1, . . . , N } We compute the L 2 (0, T, H 1 (B))-distance between p ε restricted to the graph, q ε obtained with these simulations and the pressure P predicted by the asymptotic model. The results are reported in Tables 5 and6. We also show in Figures 4,5, 6 and 7 the velocity magnitude and the pressure inside Ω ε . In the figures 6 and 7, the velocity for the asymptotic problem is obtained using the scheme proposed in Section 3.5.3. Notice that that the velocity magnitude is not convex in the cross-section, unlike the stationary Poiseuille profile.

According to [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF], the junction boundary layers have a space width of order of ε. Thus, the order of convergence of p ε -q ε , p ε -P in L 2 (0, T, H 1 (B)) norm is 1/2. The order of convergence of q ε -P in L 2 (0, T, H 1 (B)) norm is 1. This is compatible with what we observe. However, we are limited by the accuracy of the numerical resolution of the Navier-Stokes equations. Anyway, the asymptotic model shows a good agreement with this numerical solution, with only 6% error when ε = 0.025 in the multiply connected geometry case and 1% error when ε ≤ 0.1 in the simply connected geometry case. We also run simulations in the planar case, with a cut of the multiply connected geometry. This enabled us to use a much smaller space-step and thus to reach a better precision. We take a cut of the multiply connected geometry in the xy-plane and run a simulation in the interval [0, 0.25ε 2 ] with the Dirichlet boundary conditions for the velocity at the extremities of the three tubes. We use a characteristic Galerkin method, with P 3 -P 2 Taylor-Hood elements in space and BDF2 time integrator (scheme (12) in [START_REF] Bermejo | Lagrange-Galerkin methods for the incompressible Navier-Stokes equations: a review[END_REF]). It was implemented with FreeFem++ [START_REF] Hecht | New developments in FreeFem++[END_REF]. Results are reported on table 7 and match the theoretical results of [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF]. For this planar case, we also computed qε the function (not necessarily continuous) affine on each edge such (q ε -p ε )(

2 3 O ij + 1 3 O kj ) = 0, (q ε -p ε )( 1 3 O ij + 2 3
O kj ) = 0 for j ∈ {1, . . . , M }. We observed qε converges very rapidly to P on the edges where the flux is prescribed (i.e. e 1 , e 2 , e 3 ) in L 2 (0, T, H 1 dc (e 1 ∪ e 2 ∪ e 3 )) norm, with a 0.003% relative error when ε = 0.1. According to [START_REF] Panasenko | Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure[END_REF], the convergence should be exponential.

Pressure Velocity magnitude

Velocity magnitude 6: Comparison between the asymptotic model (dashed lines) and the Navier-Stokes numerical solution (blue lines) for the simply connected geometry when T = 0.875ε 2 , ε = 0.1. On the left, the pressure along tubes. On the right, the velocity magnitude across the middle of the three tubes with respect to the distance to the axis of the tube Figure 7: Comparison between the asymptotic model (dashed lines) and the Navier-Stokes numerical solution (blue lines) for the multiply connected geometry when T = 0.0875ε 2 , ε = 0.1. On the left, the pressure along tubes. On the right, the velocity magnitude across the middle of the six tubes with respect to the distance to the axis of the tube Table 5: Comparison between Navier-Stokes numerical solution and the asymptotic model in the case of the simply connected geometry. T = 0.875ε 2 . In this table, Φ ε j , Φ j is the flux accross the j-th tube according to Navier-Stokes numerical solution and the asymptotic model. The last line of the table is an estimate of the accuracy of the numerical method used to compute the Navier-Stokes solution. 7: Multiply connected two-dimensional case. Φ ε j , Φ j , Ψ ε j denotes the same fluxes as in Table 6. Notice in column 5 that qε is very accurately predicted by P on edges (j ∈ {1, 2, 3}) on which the flux is imposed.

ε P -p ε L 2 (0,T,H 1 (B))
p ε L 2 (0,T,H 1 (B))

p ε -q ε L 2 (0,T,H 1 (B))
p ε L 2 (0,T,H 1 (B))

P -q ε L 2 (0,T,H 1 (B))

q ε L 2 (0,T,H

Figure 1 :

 1 Figure 1: Example of thin structure when d = 2, N 1 = 6, N = 8, M = 7.

Figure 2 :

 2 Figure 2: ∞ -error on ∂P ∂x (e) curves obtained from the data presented in Table4(top: β = 0) and Table3(bottom: β = 1). On each graph, only one parameter varies, the two others are set by default to h = 2 -10 , k = 0.1 • 2 -14 , H = π2 -10 .

2 , 0) e 5 Figure 3 :

 253 Figure 3: The full dimension Navier-Stokes simulations were run on two geometries. On the left, the simply connected one (N 1 = 3,N = 4, M = 3), and on the right, the multiply connected one (N 1 = 3,N = 6, M = 6).

Figure 5 :

 5 Figure 5: Pressure and velocity magnitude in the transversal and longitudinal cuts of the multiply connected geometry at t = 0.0875ε 2 , ε = 0.1.

Figure

  Figure 6: Comparison between the asymptotic model (dashed lines) and the Navier-Stokes numerical solution (blue lines) for the simply connected geometry when T = 0.875ε 2 , ε = 0.1. On the left, the pressure along tubes. On the right, the velocity magnitude across the middle of the three tubes with respect to the distance to the axis of the tube

L 2

 2 (0,T,H 1 (B)) 0.0233166 0.0144485 0.69 0.0093076 0.63 0.0114437 -.30(Φ j -Φ ε j ) j L 2 ({1,...,M }×[0,T ]) (Φ ε j ) j L 2 ({1,...,M }×[0,T ])0.0054451 0.0088381 -.70 0.0166595 -.91 0.0288674 -.79

p ε L 2

 2 (0,T,H 1 (B)) 0.144626 0.103521 0.48 0.080028 0.37 0.062730 0.35p ε -q ε L 2 (0,T,H 1 (B)) p ε L 2 (0,T,H 1 (B))0.140002 0.101388 0.47 0.073540 0.46 0.056146 0.39P -q ε L 2 (0,T,H 1 (B)) q ε L 2 (0,T,H 1 (B))0.036639 0.021016 0.80 0.031650 -.59 0.028022 0.18(Φ j -Φ ε j ) j L 2 ({1,...,M }×[0,T ]) (Φ ε j ) j L 2 ({1,...,M }×[0,T ])0.035878 0.030196 0.25 0.056603 -.91 0.053661 0.08(Ψ -Ψ ε ) L 2 ({1,...,N }×[0,T ]) (Ψ ε ) L 2 ({1,...,N }×[0,T ])0.035841 0.030109 0.25 0.056649 -.91 0.053275 0.09 Table

  .32) Let us now estimate p h,k -p h,k . Using (3.19), we have:

Table 3 :

 3 ∞ -error on

				∂P ∂x (e) for the disc cross-section, β = 1, p(x, t) = e (1-t)x-1/(t+ε) .
		kernel		Numerical	Corrected	Exact
	h	k	H	∞ -error p x order	∞ -error p x order	∞ -error p x order
	2 -1			1.37116e-06	1.37126e-06	1.37063e-06
	2 -2			3.98424e-07 1.78 3.98532e-07 1.78 3.97835e-07 1.78
	2 -3			1.07928e-07 1.88 1.08041e-07 1.88 1.07306e-07 1.89
	2 -4			2.85134e-08 1.92 2.86301e-08 1.92 2.78745e-08 1.94
	2 -5 0.1 • 2 -14 2 -10 π 7.75301e-09 1.88 7.87126e-09 1.86 7.10518e-09 1.97
	2 -6			2.44720e-09 1.66 2.56627e-09 1.62 1.79485e-09 1.99
	2 -7			1.10689e-09 1.14 1.22637e-09 1.07 4.52255e-10 1.99
	2 -8			7.70502e-10 0.52 8.90184e-10 0.46 1.14715e-10 1.98
	2 -9			6.86454e-10 0.17 8.06241e-10 0.14 3.00928e-11 1.93
	2 -10			6.65557e-10 0.04 7.85396e-10 0.04 8.90762e-12 1.76
		0.1 • 2 -1		1.40603e-05	1.50528e-05	1.70035e-05
		0.1 • 2 -2		8.23443e-06 0.77 9.36902e-06 0.68 1.16079e-05 0.55
		0.1 • 2 -3		1.46155e-06 2.49 3.60481e-06 1.38 4.71718e-06 1.30
		0.1 • 2 -4		7.63423e-07 0.94 1.03297e-06 1.80 1.49344e-06 1.66
		0.1 • 2 -5		7.67132e-07 -0.01 2.91920e-07 1.82 4.20353e-07 1.83
		0.1 • 2 -6		4.14515e-07 0.89 6.10037e-08 2.26 1.11756e-07 1.91
	2 -10 0.1 • 2 -7 2 -10 π 1.82874e-07 1.18 1.23173e-08 2.31 2.88777e-08 1.95
		0.1 • 2 -8		7.31796e-08 1.32 1.70640e-09 2.85 7.35660e-09 1.97
		0.1 • 2 -9		2.75420e-08 1.41 4.54956e-10 1.91 1.86343e-09 1.98
		0.1 • 2 -10		9.73626e-09 1.50 4.85544e-10 -0.09 4.73955e-10 1.98
		0.1 • 2 -11		3.03157e-09 1.68 6.40152e-10 -0.40 1.24211e-10 1.93
		0.1 • 2 -12		5.59160e-10 2.44 7.17743e-10 -0.17 3.64187e-11 1.77
		0.1 • 2 -13		3.40868e-10 0.71 7.60524e-10 -0.08 1.44120e-11 1.34
		0.1 • 2 -14		6.65557e-10 -0.97 7.85396e-10 -0.05 8.90762e-12 0.69
			2 -1 π 3.69456e-05	3.68546e-05
			2 -2 π 3.36543e-07 6.78 3.85246e-07 6.58
			2 -3 π 2.97071e-06 -3.14 3.01664e-06 -2.97
			2 -4 π 7.43567e-07 2.00 7.86382e-07 1.94
	2 -10 0.1 • 2 -14 2 -5 π 8.39831e-08 3.15 6.60872e-08 3.57
			2 -6 π 9.45994e-08 -0.17 1.13113e-07 -0.78
			2 -7 π 3.63271e-08 1.38 4.33101e-08 1.38
			2 -8 π 1.13907e-08 1.67 1.22260e-08 1.82
			2 -9 π 3.07604e-09 1.89 3.17743e-09 1.94
			2 -10 π 6.65557e-10 2.21 7.85396e-10 2.02

Table 4 :

 4 ∞ -error on

				∂P ∂x (e) for the disc cross-section, β = 0, p(x, t) = e (1-t)x .
		kernel		Numerical		Corrected	Exact
	h	k	H	∞ -error p x order	∞ -error p x order	∞ -error p x order
	2 -1			4.65502e-02		4.65403e-02	4.64901e-02
	2 -2			1.36905e-02	1.77	1.36707e-02	1.77 1.36146e-02 1.77
	2 -3			3.77399e-03	1.86	3.75018e-03	1.87 3.69076e-03 1.88
	2 -4			2.42701e-03	0.64	1.02249e-03	1.87 9.61254e-04 1.94
	2 -5 0.1 • 2 -14 2 -10 π 3.19538e-03	-0.40	3.07484e-04	1.73 2.45304e-04 1.97
	2 -6			3.40548e-03	-0.09	1.24617e-04	1.30 6.19574e-05 1.99
	2 -7			3.46538e-03	-0.03	7.84674e-05	0.67 1.55663e-05 1.99
	2 -8			3.48384e-03	-0.01	6.69214e-05	0.23 3.89957e-06 2.00
	2 -9			3.49017e-03	-0.00	6.40570e-05	0.06 9.74829e-07 2.00
	2 -10			3.49261e-03	-0.00	6.33551e-05	0.02 2.43051e-07 2.00
		0.1 • 2 -1		3.74748e-01		2.51838e-01	1.02594e-02
		0.1 • 2 -2		2.51240e-01	0.58	1.67151e-01	0.59 3.50858e-03 1.55
		0.1 • 2 -3		1.70551e-01	0.56	5.48469e-02	1.61 1.21113e-03 1.53
		0.1 • 2 -4		1.17002e-01	0.54	1.84214e-02	1.57 4.20822e-04 1.53
		0.1 • 2 -5		8.09255e-02	0.53	4.56308e-03	2.01 1.46833e-04 1.52
		0.1 • 2 -6		5.63137e-02	0.52	1.82718e-03	1.32 5.13093e-05 1.52
	2 -10 0.1 • 2 -7 2 -10 π 3.93576e-02	0.52	6.16791e-04	1.57 1.78715e-05 1.52
		0.1 • 2 -8		2.75894e-02	0.51	2.57248e-04	1.26 6.13280e-06 1.54
		0.1 • 2 -9		1.93776e-02	0.51	8.02652e-05	1.68 2.00335e-06 1.61
		0.1 • 2 -10		1.36270e-02	0.51	1.50374e-05	2.42 5.48597e-07 1.87
		0.1 • 2 -11		9.59438e-03	0.51	3.52991e-05	-1.23 2.20105e-07 1.32
		0.1 • 2 -12		6.77321e-03	0.50	4.65737e-05	-0.40 2.35100e-07 -0.10
		0.1 • 2 -13		4.81817e-03	0.49	5.50126e-05	-0.24 2.40775e-07 -0.03
		0.1 • 2 -14		3.49261e-03	0.46	6.33551e-05	-0.20 2.43051e-07 -0.01
			2 -1 π 1.53506e+00		1.37410e+00
			2 -2 π	2.89715e-01	2.41	2.13087e-01	2.69
			2 -3 π	1.26877e-01	1.19	1.26814e-01	0.75
			2 -4 π	7.56310e-02	0.75	7.59343e-02	0.74
	2 -10 0.1 • 2 -14 2 -5 π	4.20970e-02	0.85	5.09920e-02	0.57
			2 -6 π	2.21178e-02	0.93	8.15857e-03	2.64
			2 -7 π	1.14523e-02	0.95	1.85572e-03	2.14
			2 -8 π	6.01148e-03	0.93	8.35024e-04	1.15
			2 -9 π	4.15333e-03	0.53	2.39586e-04	1.80
			2 -10 π 3.49261e-03	0.25	6.33551e-05	1.92

Table 6 :

 6 Comparison between Navier-Stokes numerical solution and the asymptotic model in the case of the multiply-connected geometry. T = 0.875ε 2 . In this table, Φ ε j , Φ j is the flux accross the j-th tube according to Navier-Stokes numerical solution and the asymptotic model, andΨ ε = O ∈ej α ,j Φ ε j .The last line of the table is an estimate of the accuracy of the numerical method used to compute the Navier-Stokes solution.

	ε	0.2	0.1	0.05	0.025
	P -p ε	L 2 (0,T,H 1 (B))

  -Φ ε j ) j L 2 ({1,...,M }×[0,T ]) (Φ ε j ) j L 2 ({1,...,M }×[0,T ]) (Ψ -Ψ ε ) L 2 ({1,...,N }×[0,T ]) (Ψ ε ) L 2 ({1,...,N }×[0,T ])

	4e-01	1.84158e-01	1.74262e-01	6.04826e-02
	3e-01	1.61135e-01 0.46	1.54638e-01 0.42	4.58469e-02 0.96
	2e-01	1.33035e-01 0.47	1.29469e-01 0.44	3.08584e-02 0.98
	1e-01	9.51649e-02 0.48	9.38936e-02 0.46	1.55718e-02 0.99
	7e-02	7.99077e-02 0.49	7.91613e-02 0.48	1.09305e-02 0.99
	5e-02	6.76983e-02 0.49	6.72472e-02 0.48	7.82038e-03 1.00
	4e-02	6.06171e-02 0.50	6.02939e-02 0.49	6.26215e-03 1.00
	3e-02	5.25627e-02 0.50	5.23528e-02 0.49	4.70006e-03 1.00
	2e-02	4.29754e-02 0.50	4.28610e-02 0.49	3.13686e-03 1.00
	7e-03	2.54107e-02 0.50	2.53871e-02 0.50	1.09576e-03 1.00
	5e-03	2.15416e-02 0.49	2.15236e-02 0.49	8.81800e-04 0.65
	ε	(Φ j P -qε qε L 2 (0,T,H 1 L 2 (0,T,H 1 dc (e1∪e2∪e3)) dc (e1∪e2∪e3))
	4e-01	1.61063e-02	1.99512e-04	1.06522e-06
	3e-01	1.04089e-02 1.52	3.58527e-05 5.97	4.56238e-09
	2e-01	6.08976e-03 1.32	3.34245e-06 5.85	8.20514e-09
	1e-01	2.71136e-03 1.17	3.20877e-06 0.06	2.40961e-08
	7e-02	1.84568e-03 1.08	4.59828e-06 -1.01	4.19010e-08
	5e-02	1.28786e-03 1.07	8.35613e-06 -1.78	1.63812e-07
	4e-02	1.02498e-03 1.02	1.21299e-05 -1.67	2.76485e-07
	3e-02	7.58649e-04 1.05	1.71698e-05 -1.21	3.57710e-07
	2e-02	5.01930e-04 1.02	2.57071e-05 -1.00	8.34215e-07
	7e-03	1.74694e-04 1.01	8.98785e-05 -1.19	1.82707e-05
	5e-03	1.37612e-04 0.71	3.15042e-04 -3.73	6.95978e-05

1

(B))

,q+1 ) are straightforward to compute. Note that the system to inverse is the same as for the block triangular approach. If we denote R the time to compute the Choleski factorization of (I -k(T (σ) H ) -1 ) and r the number of non-zero elements of the corresponding factors, then the total cost is O N S + N

+ N R + QS(N S + N 2 )r .

Acknowledgements

The authors are grateful to Konstantin Pileckas (University of Vilnius) for useful discussions. This work was partially supported by the European Social Fund according to the activity "Improvement of researchers qualification by implementing world-class R and D projects" No. 09.3.3-LMT-K-712.

Table 2: β = 0, p(x, t) = e (1-t)x . In the last three lines, the order is computed with respect to H. h k H ∞ -error p x order ∞ -error p order 0.1 • 2 -3 2.23431e-01 1.37143e-01 0.1 • 2 -4 1.42094e-01 0.65 8.75474e-02 0.65 0.1 • 2 -5 9.28145e-02 0.61 5.73905e-02 0.61 0. Since the cross-section is a disc, one could have used the radial symmetry to obtain a faster and more accurate scheme. However, the goal here is to obtain a benchmark for an arbitrary two-dimensional cross-section. Hence, we used unregular meshes for the disc.

In this section, we test the scheme (3.4.2) with three different types of choice for ( K(σ) q ): the exact (K (σ) q ), and two approximations (K (σ) num,q ) and (K (σ) cor,q ) that are defined as in Section 3.4.3. For these three approximations, we compute the ∞ -error on ∂P ∂x and the results are reported in Table 4 for the case β = 0, and Table 3 when β = 1. The corresponding error curves are given in figure 2. The values of θ(k) for the numerical and corrected approximations will be studied in detail in the forthcoming second part of this article. Let us summarize the orders with respect to h, k, H that we have estimated:

In all cases, the use of the exact kernel is the best option or close to it. When β = 1, the three schemes yields similar accuracy if the edge space step h in the cross-section space step H is small enough. When β = 0, the numerical approximation is slowly convergent and the corrected scheme offers a substantial improvement.

Comparison with Navier-Stokes equations

The aim of this section is to compare numerically the full 3D-Navier-Stokes equations (1.1) with the continuous asymptotic model on the graph (1.2).

Geometry and boundary condtions

We do it for two different configurations, one constructed from a simply connected graph and the other one from a connected graph with a hole (see Figure 3). In both cases, let Ω ε be the interior of {M ∈ R 3 |∃i, O ∈ e i , OM ⊥ e i , OM < ε}, i.e. the union of cylinders whose diameters are ε and whose axes are the e i segments. The domain Ω ε is not smooth at the junctions, but we observe the convergence to the asymptotic problem (1.2) anyway. Let us take the following boundary conditions: