Megdouda Benamara

Marjorie Grzeskowiak

Antoine Diet

Gaëlle Lissorgues

Yann Le Bihan

Equivalent Input and Output Impedances in HF RFID System Including Resonator

Keywords: antenna, impedance matrix, RFID, magnetic coupling

High Frequency Radio Frequency Identification (HF RFID) system based on Magnetically Coupled Reader Resonator Coils (MCRRC) is reported. The proposed system consists of reader antenna including small resonant coil operating by magnetic coupling with the tag coil. In the proposed system, the reader and tag impedances are modified. The equivalent electrical model is used to express the equivalent impedance matrix and used to express the equivalent input and output impedances of the system. The formulas are confirmed by comparison between High Frequency Structure Simulator (HFSS) results and measures.

INTRODUCTION

RFID (Radio Frequency Identification) HF (High Frequency, 13.56 MHz) is based on inductive coupling effect between the reader and the tag coils. According to the distance and the positioning (planar and angular) between the two coils, the impedance of the tag chip is modified, the shift of the impedance is interpreted by the reader antenna, this is the principle of the load modulation [START_REF] Finkenzeller | RFID Handbook[END_REF].

In several applications like as tracking, the tag is misaligned and has small size coil compared to the reader coil. In this case, the reader antenna has to be changed in multi-loop reader antenna to ensure detection [START_REF] Grzeskowiak | Pebbles Tracking Thanks to RFID LF Multiloops Inductively Coupled Reader[END_REF][3] [START_REF] Diet | Improvement of RFID HF tags detection with a distributed diameter coil, Antennas and Wireless Propagation Letters[END_REF]. In this study the proposed concept is to changed reader antenna on magnetically coupled reader resonator coils (MCRRC), the link between the reader including resonator and the tag depends on the equivalent mutual impedance [START_REF] Benamara | Calculation of the Equivalent Mutual Impedance in Complex HF RFID Systems[END_REF][6] [START_REF] Waffenschmidt | Homogeneous Magnetic Coupling for free positioning in inductive wireless power system[END_REF] [START_REF] Zhang | Quantitative Analysis of Mutual Inductance for Optimal Wireless Power Transfer via Magnetic Resonant Coupling[END_REF]. Preliminary results of the proposed design confirm its interest in improving both surface (in the case of lateral misalignment of the tag above the reader antenna surface) and distance of detection [START_REF] Diet | Improvement of RFID HF tags detection with a distributed diameter coil, Antennas and Wireless Propagation Letters[END_REF] [START_REF] Benamara | Calculation of the Equivalent Mutual Impedance in Complex HF RFID Systems[END_REF]. The difficulty is the need of adjusting the impedance of the matching circuit of the reader antenna according to the electrical properties of the resonator in order to maintain the resonance frequency at 13.56MHz [START_REF] Benamara | Calculation of the Equivalent Mutual Impedance in Complex HF RFID Systems[END_REF].

The proposed structure seen in Fig. 1 represents the interaction between the two reader coils (coil 1 and resonator) and the coil of the tag [START_REF] Grzeskowiak | Pebbles Tracking Thanks to RFID LF Multiloops Inductively Coupled Reader[END_REF]. Resonator coil has a radius Rres of 15 mm. Coils 1 and 2 present respectively a radius of 50 mm and 15 mm. In [START_REF] Benamara | Calculation of the Equivalent Mutual Impedance in Complex HF RFID Systems[END_REF] the calculation of the equivalent impedance parameters is given, the work concerns the improvement of the equivalent mutual impedance (as link between the reader multi-coil and the tag coil) at 13.56MHz.

In this paper, all the equivalent impedance parameters of the system are studied according to the lateral misalignment of the tag above the reader antenna in the case of conventional and MCRRC RFID system with MCRRC. The calculated equivalent parameters are confirmed by comparison between HFSS simulations, calculations and measurements in section 2. In section 3, the influence of the added resonator on the equivalent input and output impedances of the system is investigated. New formulas defining the equivalent input and output impedances are given and confirmed by simulation.

A. Equivalent impedance matrix calculation

The voltages V 1 , V res and V 3 in the circuit (b) satisfy the following relations (1), with M ij =M ji being the mutual inductances between coils 1, 2, 3, (1<i,j<3)) The equivalent impedance matrix (Z IO matrix) of the MCRRC system (b) is given by (the complete calculation of the equivalent matrix is developed in [START_REF] Benamara | Calculation of the Equivalent Mutual Impedance in Complex HF RFID Systems[END_REF]):

                    ω + ω Μ 21 ω ω ω + ω Μ 12 ω ωΜ ω + =                     =           I I res I L j r M res j j M res j L res j r res M res j j j L j r I I res I Z Z res Z Z res Z resres Z res Z Z res Z V V res V 2 1 2 2 2
Z C L j M res r Z oo M res M res M j Z oi M res M res M j Z io Z C L j M res r Z ii 2 2 2 2 2 2 1 2 2 21 2 1 2 12 1 1 1 2 2 1 + + - = - = - = + + - = ω α ω α ω ω α ω ω ω α ω (2)

With

( )

r res L res j C res j + - = ω ω α
The parameters of the Z IO matrix can be expressed as:

( ) ( ) r M res M res C L j Z M res M M res M res M j eq Z M r M r M res M res M j eq Z r M res M res C L j Z eq res eq 2 2 1 ) ( 2 2 2 ) ( 2 2 2 2 1 2 22 1 2 1 2 21 2 1 2 2 1 12 . ) ( 1 2 2 ) ( 1 2 2 1 1 1 11 + +         + - = + + = + + = + +         + - = δ ω γ ω ω ω δ ω γ ω ω δ ω γ ω ω δ ω γ ω ω ω (3) With ( ) ( ) r res L res C res r res r res C res L res L res C res + ω - ω = + ω - ω ω - ω = 1 2 , 1 2
1 δ γ
In MCRRC system, the resonator influence is depicted in both reader and tag impedances by the M 1res and M 2res parameters but also by impedance coupling (in this case, the mutual impedance has a real and an imaginary part).

B. Validation of the calculated equivalent impedance matrix

The equivalent impedance matrix parameters are analyzed in this part according to the tag misalignment in both cases with and without resonator. These parameters are studied in this part by calculation, HFSS simulations (using the analytical model in equation) and confirmed by measurements. The distance between the reader and the tag coils is chosen to 10mm: at this distance, the impact of the resonator is seen clearly on the mutual impedance in [START_REF] Benamara | Calculation of the Equivalent Mutual Impedance in Complex HF RFID Systems[END_REF].

The simulated results concern structure where the coils are designed with perfect conductor to reduce the computation time under HFSS, while in measure, the results are presented for realized coils with cooper. The different materials used respectively in simulation (without substrate) and measure (with FR4 substrate) generate different electrical properties (consequently different resonance frequencies). The aim of this part is to confirm the developed electrical model by comparison of the evolution of the equivalent impedance parameters around the resonance frequency.

B.1.

Impedance of the reader (Z11 and Z 11eq )

For three positions of the misaligned tag: y=-50 mm, y=-30 mm and y=0mm in Fig. 1 corresponding respectively to the positions where the tag is above the edges of the reader antenna, the center of the resonator and the center of the reader antenna, we reported in fig. 3 In the case without resonator (Fig. 3 a,b), the resonance frequency is shifted according to the position of the tag, but the frequency shift, in this case, can be negligible, the resonance frequency is around 13.55MHz. The added resonator shifted the frequency; the resonance frequency of the Z 11eq impedance is around 13.27 MHz for the different lateral positions of the tag. The amplitude of the real part is modified according to the tag misalignment. Its maximum in the case without resonator is around 0.29 Ω. This maximum is improved with the added resonator: at y=0 mm the variation of the real part is negligible compared to the case of y=-50 mm and y=-30mm where the real part of the Z 11eq impedance is strongly influenced by the added resonator. Also, the results are confirmed by measurements. In measured results the frequency shift according to the tag misalignment for the structure with resonator is confirmed.

B.2.

Impedance of the reader (Z 22 and Z 22eq ) In the same way, fig. 4 reports, the variation of the real and the imaginary part of the self-impedance of the tag in conventional (Z 22 ) and MCRR system (Z 22eq ). In conventional system (without resonator), the resonance frequency of the tag is seen around 13.6MHz, according to the tag positioning, both resonance frequency (null of the imaginary part) and the amplitude of the real part are changed. This variation, is negligible compared to the variation in the case of MCRR system, the maximum variation of the real and the imaginary parts of the tag impedance is seen for y=-30 mm (when the tag is above the resonator). The results of the calculated equations are conformed by measurements, where the resonance frequency of the measured tag is around 15MHz.

B.3.

Impedance of the reader (Z12 and Z 12eq ) In this part, we are also interesting on the mutual impedance parameter. In conventional system, the imaginary part is used to calculate the mutual inductance between the reader and the tag coils, while the real part is very low and close to zero. A few variations of these parameters are seen according to the tag misalignments. While in MRCC system, both real and imaginary parts of the Z 12 impedance are changed, in simulation, the maximum of this variation is seen around 13 MHz for y=-30 mm. in measurements, the variation has maximum values at 25 MHz. These correspond to the resonance frequencies of the reader and the tag coils in MCRR system. In this part, we confirm that: -The added resonator modifies all the parameters of the impedance matrix.

-The developed calculation in equation ( 3) is validated by comparison between simulations and measures. -The interest of the added resonator is depicted on the improvement of the mutual impedance of the system.

-Both the equivalent self-impedances of the reader and the tag coils are modified; consequently a frequency shift is generated -According to the lateral positioning of the tag, the variation of the impedance parameters is more important when the tag is above the resonator (y=-30 mm).

III. EFFECT ON INPUT AND OUTPUT IMPEDANCES

The MCRRC system is benefic for improving generated magnetic field by the reader coil, but the main constraint of such system is the influence on the impedance parameters. In this part, the input and output impedances are studied by calculation and simulation. The input and the output impedances are defined using the equivalent impedance matrix parameters in equation [START_REF] Zierhofer | Geometric approach for coupling enhancement of magnetically coupled coils[END_REF].

A. Calculation of equivalent input and output impedances

In conventional structure (a), the input and output impedances of the system, respectively the impedance of the reader and of the tag, are depending on the self-impedances of the reader and the tag and the mutual inductance between them:

r L j r L j r Z Z Z Z Z out Z L j r L j r Z Z Z Z Z in L L chip chip M M + + + + = + - = + + + + = + - =
( ) ( )         + + + + +         + + + + - - = 2 2 (5) With: ( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( )( )
Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z re Z Z Z ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )         + + + + - +         + + + + - = 2 2
Z Z Z Z Z Z Z Z Z Z Z Z Z Z re Z Z Z
Also, the equivalent output impedance is expressed by

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )         + + + + +         + + + + - - ) 22
Im( Im L eq L eq L eq eq eq L eq L eq L eq eq eq eq eq out

Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z re Z Z Z ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( )  
      + + + + - +         + + + + - = 2 2 2 2 2
Z Z Z Z Z Z Z Z Z Z Z Z Z Z re Z Z Z
From the calculation, both input and output impedances are modified by the added resonator.

B. Simulation results

The equivalent input and output impedances are studied in this part by HFSS simulations versus the lateral misalignment at resonance frequency of 13.56MHz (using the analytical model in equations 4, 5 and 6) (Fig. 6). In this part, the distance between the reader and the tag coils is 10 mm. The study is done for lateral misalignment of the tag coil above the reader coil (-50 mm<Y<50 mm). The magnitudes of the input and output impedances are changed with the added resonators. This is due to the transformed impedances in equation (5 and 6). In both reader and tag coils, the variation of the impedances due to the mutual inductance is negligible compared to the variation due to the added resonator: the improvement of the reader and the tag impedances corresponds respectively to 5Ω and 3Ω.

In conventional reader antenna, the phase shift is depending on the tag positioning: the mutual inductance between the reader and the tag coils is the only cause of the phase shift. While, in MCRR structure, the evolution of the phase is depending on the mutual coupling between the reader and the tag coil (M in equation 4) but also on the electrical properties of the resonator (γ) and the added equivalent mutual inductances (M1res and M 2res ). The phase shifts in this case are higher than conventional structure. This means, a frequency shift due to the transformed impedances.

IV. CONCLUSION

In this paper, the improvement of mutual impedance, in HF RFID system, is achieved thanks to the addition of a coplanar resonator, made of coil with a size in the range of those of the tag coil (with geometrical similarity).

Analytical formulas for equivalent impedance parameters of the system are developed using electrical model of the structure, constituted of dual-coils reader and tag. The calculation is confirmed by comparison between HFSS simulations and VNA measurements. These parameters are used to express the equivalent input and output impedances of the system.

The added resonator is benefic to improving mutual coupling between the reader and the tag coils, but it creates a phase shift on both equivalent input and output impedance. The use of such system demands the modification of the matching circuit of the reader and the tag according to the electrical properties of the resonator and its positioning inside the reader antenna surface.

The work in progress relates to automatic matching circuit based on equations (5 and 6). The shift frequency in complex HF RFID will be studied and rectified using automatic switched capacitances.

Fig. 1 .

 1 Fig.1. Fabricated tag coil 2 (left) and MCRRC (right) II. EQUIVALENT CIRCUIT FOR RFID SYSTEM Fig.2 shows the equivalent circuit of conventional and complex (MCRRC) RFID system where: L 1 = 0.3 uH, L 2 = 0.058 uH and L res = 0.058 uH are respectively the selfinductances of coils 1, 2 and the resonator. r 1 = 0.2Ω, r 2 = 0.1 Ω and r res = 0.1 Ω are their internal resistances. C 1 =0.45 nF, C 2 =2.32 nF and C res =2.32 nF are the series tuning capacitors at 13.56 MHz.

Fig. 2 .

 2 Fig.2. Electrical model for conventional (a) and MCCRC (b) RFID system

3 .

 3 the evolution of the imaginary and the real parts self-impedance of the reader in the two cases without (Z 11 ) and with resonator (Z 11eq ). Imaginary and real part of simulated Z 11 (a, b), simulated Z 11eq (c, d) and measured Z 11eq (e, f) parameters

Fig. 4 .

 4 Fig.4. Imaginary and real part of simulated Z 22 (a, b), simulated Z 22eq (c, d) and measured Z 22eq (e, f) parameters

5 .

 5 imaginary and real part of simulated Z 12 (a, b), simulated Z 12eq (c, d) and measured Z 11eq (e, f) parameters

  the equation (4), we define the equivalent input and output impedances by replacing the impedance parameters by the equivalent impedance parameters. The equivalent input impedance can be express by:

Fig. 6 .

 6 Fig.6. modules of input and equivalent input impedances (a), modules of input and equivalent input impedances (a), phases of input and equivalent input impedances (b), modules of output and equivalent output impedances (b), phases of output and equivalent output impedances (b)


Yibiao Wang

Azita Ahmadi

Didier Lasseux

On the Inertial Single Phase Flow in 2D Model Porous Media: Role of Microscopic Structural Disorder

Keywords: Inertial one-phase flow, Non-Darcy flow, Flow regimes, Darcy-Forchheimer

In this work, single-phase incompressible laminar flow in 2D model porous media is studied and the influence of microscopic structural disorder on the flow is thoroughly investigated. Emphasis is laid upon the onset of the deviation from Darcy's law and the identification of different inertia regimes observed before the flow becomes unsteady. For this purpose, six globally disordered pore structures were generated and the values of the critical Reynolds number at which the flow becomes unsteady corresponding to the first Hopf bifurcation were determined. Numerical simulations of steady laminar single-phase flow were then carried out to investigate the effects of the microstructures on the inertial correction to Darcy's law. Different flow regimes, namely weak inertia, strong inertia and the regime beyond strong inertia, are identified. Comparisons are made with results presented in the literature which were restricted to ordered and locally disordered structures. The critical Reynolds number decreases and inertia intensity increases as more disorder is introduced into the pore structure. Results on flow inertia widely extend some previous studies on the subject and show that it is mainly influenced by the shape of the obstacles (either circular or square), slightly affected by the inclination of the square cylinders and hardly disturbed by the size distribution of the obstacles.

Introduction

Inertial flow in porous media corresponds to flow regimes where the Stokes or creeping model is not valid at the microscopic scale so that Darcy's law is not applicable at the macroscopic scale. Such flows are encountered in many practical situations as in aquifers [START_REF] Liu | A generalized forchheimer radial flow model for constant-rate tests[END_REF], in filters [START_REF] Payri | Experimental-theoretical methodology for determination of inertial pressure drop distribution and pore structure properties in wall-flow diesel particulate filters (dpfs)[END_REF], near wellbores in hydrocarbon recovery or in column reactors [START_REF] Enayati | Numerical simulations of porous medium with different permeabilities and positions in a laterally-heated cylindrical enclosure for crystal growth[END_REF]. A macroscopic empirical law classically employed to take into account inertial effects is the popular Darcy-Forchheimer momentum equation [START_REF] Forchheimer | Wasserbewgung durch Boden[END_REF] given by B Didier Lasseux didier.lasseux@u-bordeaux.fr

∂ p β β ∂ x = α 1 v β + ρ β α 2 v β v β (1)
Note that in this expression, the pressure and velocity are understood in the sense of their intrinsic and superficial average, respectively, the definitions of which will be made clear in Sect. 2.1. Some fair agreement with the quadratic nature of the correction to the linear Darcy's law was found with numerical results obtained in 2D structures made of parallel and regularly arranged circular cylinders for flow orthogonal to the cylinders axes [START_REF] Lee | Modeling of Darcy-Forchheimer drag for fluid flow across a bank of circular cylinders[END_REF]. Reasonable agreement was also found in the turbulent regime for some special 2D and 3D geometries using a general RANS model (Soulaine and Quintard 2014). In addition, the form of Eq. ( 1) was assumed to be a valid macroscopic steady model which was used as the basis for a second upscaling in order to obtain a model operating in heterogeneous porous media [START_REF] Auriault | Upscaling Forchheimer law[END_REF]). However, validity of this steady model, in particular regarding the quadratic correction term, has been the object of an active research over the past 30 years (see for instance discussions in [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], [START_REF] Lasseux | On the developments of Darcy's law to include inertial and slip effects. A century of fluid mechanics: 1870-1970[END_REF] and references therein) yielding to the identification of several regimes. Nevertheless, the macroscale description of the physics which should formally result from an upscaling procedure keeps reminiscence of the microscale flow process through the inherent dependence of the macroscopic coefficients [α 1 and α 2 in Eq. ( 1)] upon the microstructure of the porous material [see an example in [START_REF] Yazdchi | On the transition from creeping to inertial flow in arrays of cylinders[END_REF]]. The analysis of the effect of the microstructure on the inertial correction, in particular its ordered character or not, was investigated in [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF] and the present work is a thorough continuation on the same lines. In addition, if a time-dependent flow is under concern, the question on the dependence of the transition to unsteadiness upon the microsctructure is even more acute. Most works reported in the literature on this subject have adopted simple ordered or slightly disordered geometries, like, in 2D, a periodic pattern of parallel cylinders of circular or square cross section to determine the critical Reynolds number at which the flow becomes unsteady [START_REF] Kawaguti | The critical Reynolds number for the flow past a sphere[END_REF][START_REF] Jackson | A finite-element study of the onset of vortex shedding in flow past variously shaped bodies[END_REF][START_REF] Kumar | Prediction of the critical Reynolds number for flow past a circular cylinder[END_REF][START_REF] Lankadasu | Onset of vortex shedding in planar shear flow past a square cylinder[END_REF], to investigate the wake behavior [START_REF] Gera | CFD analysis of 2D unsteady flow around a square cylinder[END_REF][START_REF] Kelkar | Numerical prediction of vortex shedding behind a square cylinder[END_REF], to study the interactions of wakes and evaluate the critical Reynolds number [START_REF] Mizushima | Vortex shedding from a row of square bars[END_REF] or to study the effects of cylinder spacing on flow [START_REF] Kumar | Simulation of flow around a row of square cylinders[END_REF]). There are also many studies using more complex pore-like geometries to investigate the flow properties. In the work by Amaral Souto and Moyne (1997), a regular square tessellation of the 2D unit cell used for the computational domain and a random distribution of the phase indicator (solid or fluid) on the resulting grid was performed to compute the inertial flow, yielding the confirmation of the cubic correction to Darcy's law in a weak inertia regime identified earlier [START_REF] Barrère | Modélisation des écoulements de Stokes et Navier-Stokes en milieu poreux[END_REF][START_REF] Mei | The effect of weak inertia on flow through a porous medium[END_REF]. Similarly, a disordered structure was considered by [START_REF] Andrade | Inertial effects on fluid flow through disordered porous media[END_REF] that was obtained by removing square obstacles from a 64 × 64 square lattice until a pore space with a prescribed porosity was generated. It was used to statistically characterize the transition from linear to nonlinear flow in terms of the distribution of kinetic energy. More recently, [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF] subdivided the periodic 2D domain under consideration into equal-sized unit cells with one solid inclusion of eventually random size placed randomly inside each unit cell. In this way, they studied the effects of locally disordered (LD) structures on the inertial correction to Darcy's law. Similar structures were adopted by [START_REF] Agnaou | Une étude numérique des écoulements mono et diphasique inertiels en milieux poreux[END_REF], [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF] to study the first Hopf bifurcation for one-phase flow or to investigate the origin of inertial effects through a detailed investigation of the flow structure.

Although many studies have been dedicated to the analysis of flow in different types of pore structures, a more systematic study on the effects of the microstructure on the onset and on the role of inertia, as well as on the occurrence of unsteadiness, is still lacking.

Based on the above considerations, the present study is a continuation of previous investigations [START_REF] References Abbasian Arani | Sur quelques aspects des écoulements inertiels mono-et diphasique en milieu poreux[END_REF][START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | Une étude numérique des écoulements mono et diphasique inertiels en milieux poreux[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF] taking into account six different types of 2D globally disordered (GD) structures with solid inclusions of circular or square cross sections. The focus is laid upon 2D configurations, keeping in mind that 3D cases are also of major interest, in particular with the objective of comparing the two situations and analyzing whether the fact that the inertial correction departs from the quadratic (Forchheimer) form can be reasonably attributed to 3D effects as suggested in [START_REF] Fourar | On the non-linear behavior of a laminar single-phase flow through two and three-dimensional porous media[END_REF]. Compared to ordered and locally disordered (LD) structures studied in the above-cited references, the accent is laid upon much stronger disorder, considering random size distributions and orientations of the solid obstacles. Comparisons are made between GD structures and simple geometries mentioned above and a thorough analysis is carried out in order to investigate the effects of the microscopic pore structures on the flow regimes considering inertia.

The incipient departure from Darcy's law due to the contribution of inertia has been observed by many researchers [START_REF] Edwards | The influence of Reynolds number upon the apparent permeability of spatially periodic arrays of cylinders[END_REF][START_REF] Koch | Moderate Reynolds number flows through periodic and random arrays of aligned cylinders[END_REF][START_REF] Andrade | Inertial effects on fluid flow through disordered porous media[END_REF][START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF] in the laminar flow regime prior to the appearance of turbulence. The flow becomes unsteady for a Reynolds number larger than a critical value in the laminar regime. However, in some works on inertial effects [START_REF] Andrade | Inertial effects on fluid flow through disordered porous media[END_REF][START_REF] Chai | Non-Darcy flow in disordered porous media: a lattice-Boltzmann study[END_REF][START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], the flow was just assumed to be stationary in a limited Reynolds number range without a pre-determination of the critical Reynolds number. Therefore, in order to ensure a more rigorous and complete analysis in terms of the domain of validity of steady flow, the critical Reynolds number was first determined for the structures considered in this paper. The component of the average velocity orthogonal to the mean flow was monitored to capture flow unsteadiness. Then, steady numerical simulations were carried out to identify the different inertia regimes corresponding to weak, strong and beyond strong inertia [START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF] as well as to further assess the effects of pore structure disorder on the inertial contributions to the deviations from Darcy's law.

The paper is organized as follows. Section 2 is dedicated to a presentation of the microand macroscopic models for incompressible flow and the different pore structures used to carry out the simulations together with the numerical procedure are detailed. Results are reported in Sect. 3 with, successively, the critical Reynolds number (Sect. 3.1), the influence of the microstructure on the permeability (Sect. 3.2.1) and on the flow regimes and inertia (Sect. 3.2.2). Conclusions are drawn in Sect. 4.

Physical Models: Methodology and Tools

In this section, the governing equations for stationary isothermal flow of an incompressible Newtonian fluid (β-phase) in a porous medium at the microscale and at the macroscale are first recalled, together with the closure problem relating the two, and the inertial correction vector on which discussions will be based is introduced. The six types of pore structures used to analyze the flow are also listed, with a specific nomenclature in order to clearly distinguish them from those considered in previous studies. In addition, numerical algorithms adopted for the two softwares employed in this work are briefly introduced.

Governing Equations: Microscopic and Macroscopic Models

Model porous media under consideration are periodic 2D structures so that the analysis can be restricted to a periodic unit cell, V , representative of the process. For convenient comparison with existing reported results, the governing mass (continuity) and momentum (Navier-Stokes) conservation equations, together with the periodic boundary conditions, are written in a dimensionless form, yielding the following initial boundary value problem in this periodic domain [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF])

Re * ∂v * β ∂t * + v * β • ∇v * β = -∇ p * β + ∇ 2 v * β -∇ p * β β in V β (2a) ∇ • v * β = 0 in V β (2b) v * β = 0 at βσ (2c) v * β = 0 at t * = 0 in V β (2d) p * β β = 0 t * > 0 in V β (2e) v * β (r * + l * i ) = v * β (r * ) p * β (r * + l * i ) = p * β (r * ) i = 1, 2, 3 ( 2 f )
where, for the sake of simplicity and consistent comparison with existing results, no body force was taken into account. Here,V β and βσ represent the volume of the β-phase contained in V and the fluid-solid (β-σ ) interface, respectively; v * β , p * β are the dimensionless β-phase velocity and pressure. Time and position variables are denoted t * and r * and l * i (i = 1, 2, 3) represents the dimensionless periodic lattice vectors of the unit cell. The pressure p * β in Eq. (2a) has been replaced by its expression given by Gray's decomposition [START_REF] Gray | A derivation of the equations for multi-phase transport[END_REF]) in terms of the intrinsic average pressure, p * β β , and its deviation, p * β , namely

p * β = p * β β + p * β (3)
where the intrinsic average of any quantity ψ in V β is defined as

ψ β = 1 V β V β ψ dV = -1 ψ (4)
denoting the porosity of the structure and ψ the superficial average of ψ. The decomposition in Eq. (3) allows to impose periodicity for the pressure deviation, while considering ∇ p * β β as a driving force. Dimensionless pressure and velocity are defined from their dimensional counterparts using the reference values

v ref = l 2 μ β |∇ p β β | and p ref = l|∇ p β β |, yielding ∇ p * β β
= 1 so that any flow can be obtained from an appropriate independent choice of the Reynolds number and orientation of ∇ p * β β [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF]. The Reynolds number Re * in Eq. (2a) is defined by:

Re * = ρ β v ref l μ β = ρ β l 3 μ 2 β |∇ p β β | (5)
In the above equations, l represents the size of the periodic geometrical unit cell in the case of ordered structures (see Fig. 1a andb) or, in the case of locally and globally disordered structures, the size of the periodic unit cell divided by the square root of n, n being the number of obstacles in this unit cell (see Figs. 1c andd, 2). In addition, ρ β and μ β are the density and dynamic viscosity of the β-phase.

Under the assumptions of quasisteady flow, a scale hierarchy (l L or √ nl L, L being the size of the macroscopic medium considered as a periodic structure in its homogeneous part) and

ρ β | v β β |d p μ β L/d p
where d p is the characteristic pore size [START_REF] Whitaker | Advances in theory of fluid motion in porous media[END_REF][START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF][START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], the pore-scale boundary value problem can be upscaled to yield the macroscopic model using the volume averaging method [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF][START_REF] Whitaker | The Method of Volume Averaging: Theory and Applications of Transport in Porous Media[END_REF]. In the dimensionless form, it is given by [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF])

∇ • v * β = 0 in V (6a) v * β = -K * • ∇ p * β β -F • v * β in V (6b)
valid far enough from the macroscopic boundaries nearby which the scale hierarchy fails.

The dimensionless intrinsic permeability, K * = K/l 2 , and the Darcy's law correction tensor, F, can be determined by solving the associated closure problem for Re * = 0 and for any prescribed value of Re * respectively. In its dimensionless form, this closure problem is given by

Re * v * β • ∇M * = -∇m * + ∇ 2 M * + I in V β (7a) ∇ • M * = 0 in V β (7b) M * = 0 at βσ (7c) m * β = 0 (7d) m * (r * + l * i ) = m * (r * ) M * (r * + l * i ) = M * (r * ) i = 1, 2, 3 ( 7 e )
where m * and M * are closure variables that map v * β onto p * and ṽ * respectively [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF]. It must be noticed that an alternative form of the inertial term in Eq. ( 7a) can be used with the advantage of being fully closed in contrast with the expression Re * v * β • ∇M * which seems to indicate that v * β is required to solve this closure problem. In fact, as shown in Valdés-Parada et al. ( 2016), the left-hand side of Eq. (7a) can be written as Re * ∇ p β β • M T • ∇M * , with the drawback, however, of introducing a strong nonlinearity which makes the closure problem more complex to solve. From this closure, it follows that

M * = H * (8)
with M * = K * when Re * = 0, and

F = K * • H * -1 -I (9)

Determination of the Inertial Correction

The inertial correction to Darcy's law is given by -F • v * β as can be seen in Eq. (6b). For convenience, a renormalization is applied leading to express the correction as a dimensionless vector, f c , as [see [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF] for the details and the physical interpretation]

f c = -F • v * β | v * β | (10)
An alternative expression can be employed, avoiding the computation of the tensor F, and this is obtained by making use of Eq. ( 6b), which yields Under this form, the determination of f c only requires the solution of the flow within the unit cell resulting from ∇ p * β β and a given Reynolds number, together with the permeability tensor. Rather than Re * , a Reynolds number, depending only on macroscopic quantities is used in the present work that is given by:

f c = v * β + K * • ∇ p * β β | v * β | (11)
Re k = ρ β | v β | √ k μ β = | v * β | √ k * Re * (12)
where k * = k/l 2 is the dimensionless intrinsic permeability considering an isotropic structure. It should be noted that a Reynolds number based on the grain size of the porous medium and the average velocity is also widely employed in the literature (see a discussion on the different definitions of the Reynolds number in [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], section II.B), namely

Re d = ρ β | v β |d μ β = | v * β |d * Re * (13)
where d * = d/l is the dimensionless grain size of the porous structure (Fig. 1). For the random structures used in this work, d was estimated as the quadratic mean of the obstacle sizes, i.e., d = i=n i=1 d 2 i /n. From the above expression of f c , it appears that when the correction to Darcy's law,

-F • v * β , scales as | v * β | 3 (respectively | v * β | 2 )
as expected in the weak (respectively strong) inertia regime [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], then f c depends linearly on Re 2 k (respectively Re k ). Therefore, in this paper, the interval of validity of each regime will be identified from the relationship between f c and Re k .

Pore Structures

The nomenclature of the disordered structures considered in this work shall be first briefly introduced in order to clearly distinguish them. The random porous media studied in previous works [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF] are referred here to as locally weakly disordered (LWD, Fig. 1c) and locally strongly disordered (LSD, Fig. 1d) porous media.

Four other globally disordered (GD) 2D porous media, generated with square solid inclusions, are considered in this study, namely globally weakly disordered structures (GWD, Fig. 2a), globally weakly disordered structures with flipped squares (GWDf, Fig. 2b), globally strongly disordered structures (GSD, Fig. 2c), globally strongly disordered structures with flipped squares (GSDf, Fig. 2d). Two additional ones with circular cylinders are also considered, namely globally weakly and strongly disordered structures with circular cylinders (respectively, GWDC, Fig. 2e and GSDC, Fig. 2f). Weakly disordered structures refer to cases for which inclusions are of uniform size but are randomly placed within the unit cell while, for strong disorder, inclusion size is also randomly chosen. In contrast with the local disorder envisaged so far for which random placement was spatially confined within a region of size l, global disorder allows a random distribution of the inclusions within the unit cell of size √ nl. All the GD structures have a dimensionless size of 30×30, which means n = 900.

Indeed, in accordance with a previous study [START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF], this size (and number of inclusions) is sufficiently large to reduce artifacts related to periodic boundary conditions [START_REF] Hill | Moderate-Reynolds-number flow in a wall-bounded porous medium[END_REF]. For comparison purposes, the value of the porosity ( = 0.75) is kept the same as in previous works [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF]. In order to avoid overlapping of the cylinders and extremely narrow channels that would locally require an exceedingly large number of grid blocks, the following constraint was imposed during generation: inclusions were randomly distributed imposing a minimum pore throat size of 0.1 times the average edge size (for square inclusions) or diameter (for circular inclusions).

Numerical Tools and Procedure

In order to identify the critical Reynolds number corresponding to the first Hopf bifurcation, transient flow was first simulated in all the geometries generated in this work. This was performed by solving Eq. ( 2) using the transient solver icoFoam, which is based on the Pressure Implicit Split Operator (PISO) algorithm in the OpenFOAM package [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]). The linear system solver for the velocity is based on the Preconditioned Bi-Conjugate Gradient (PBiCG) with a tolerance of 10 -10 , while, for the pressure, the Preconditioned Conjugate Gradient (PCG) solver was chosen with a tolerance of 10 -7 . These tolerances were selected to accurately determine the critical Reynolds number intervals and are in agreement with those employed in a previous work [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF]). The time step was calculated according to the mesh size and the velocity magnitude in order to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition. The Courant number was always smaller than 0.5 in this work.

Once the critical Reynolds number delimiting the steady laminar flow regime was identified, numerical simulations of the steady incompressible flow were carried out in all the generated pore geometries, allowing a detailed analysis of the inertial correction to Darcy's law. To do so, steady-state computations were performed within the Reynolds number range that ensures steady laminar flow for which the average model in Eq. ( 6) is applicable. This was achieved by making use of the partial differential equation (PDE) module in the COM-SOL multiphysics commercial software package (COMSOL Multiphysics 2012). The system of nonlinear equations was discretized and linearized around an initial estimate using the COMSOL stationary solver. The system of linearized equations was solved with the direct solver MUMPS (Multifrontal Massively Parallel Direct Sparse Solver) [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF], which makes use of a LU factorization. The computational domain was discretized employing unstructured triangular elements with local mesh refinement in constricted zones as illustrated in Fig. 3. The average pressure gradient was such that ∇ p * β β = e x and the mean flow is expected to be aligned with the driving force due to the isotropic character of the random structures.

Results and Analysis

The hydrodynamics within the spatially periodic porous media under consideration is analyzed in order to investigate the effect of the microscopic structure on the inertial contribution to the deviation from Darcy's law. In order to ensure that results are representative of the physics at play in the corresponding structure of infinite extent and circumvent the fact that the size of the random unit cell may not be a representative domain of the medium, an ensemble average of results obtained for ten realizations of each structure was considered, except for Before proceeding any further, it is of importance to study mesh convergence for the numerical solutions [START_REF] Chai | Non-Darcy flow in disordered porous media: a lattice-Boltzmann study[END_REF]). This analysis is reported in the "Appendix."

Determination of the Critical Reynolds Number

Transition from steady to time-periodic flow within the laminar regime is marked by a Hopf bifurcation [START_REF] Marsden | The Hopf Bifurcation and Its Applications[END_REF] characterized by a critical value, Re * c (or Re dc or Re kc ), of the Reynolds number. For flow past an isolated circular cylinder, the wake behind the cylinder becomes periodically oscillating at the bifurcation which occurs at Re dc ≈ 50, leading to the formation of a Karman vortex street [START_REF] Jackson | A finite-element study of the onset of vortex shedding in flow past variously shaped bodies[END_REF]. In many engineering fields, the occurrence of the periodic vortex shedding is believed to result in the failure of structures exposed to this type of flow [START_REF] Franke | Numerical calculation of laminar vortex-shedding flow past cylinders[END_REF]). Here, the determination of this critical value is aimed at restricting our investigations of inertial effects to the steady regime.

Capturing the occurrence of unsteadiness in the random structures under consideration is not an easy task a priori. The best indicator was found to be the time evolution of the y-component of the average dimensionless velocity, v * β y , and this was further checked to be in agreement with the occurrence of unsteadiness of the pointwise y-velocity component at particular locations in the pore space. Simulations of unsteady flow were carried out over a sufficiently long time for successive increasing Reynolds number values until oscillatory flow was observed. It should be emphasized that, unlike other works in the literature where a disturbance was introduced into the simulations [START_REF] Kelkar | Numerical prediction of vortex shedding behind a square cylinder[END_REF][START_REF] Chen | Bifurcation for flow past a cylinder between parallel planes[END_REF], no perturbation was added in this study and the periodic oscillation is purely due to the physical instability of the problem, depending on the pore geometries and Reynolds number. To keep the simulation time reasonable, the critical Reynolds number values were estimated on only two distinct realizations of each structure.

As an illustration, evolutions of v * β y obtained on a GWDf realization for Re * = 12,000 (Re k = 5.48, Re d = 26.74) and Re * = 14,000 (Re k = 6, Re d = 29.31) are reported in Fig. 4a andc respectively. For Re * = 12,000, v * β y remains constant over time after an initial transient period, showing that the flow stabilizes in a steady regime. Conversely, for Re * = 14,000, v * β y exhibits a permanent oscillatory behavior, typical of a post-bifurcation regime. This is further confirmed by the Fourier spectra represented versus the Strouhal number, St, defined as

St = f l v ref = μ β f l|∇ p β β |
, f being the frequency of the oscillation, in Fig. 4b andd, respectively. These two figures show that v * β y (t) has a zero frequency at Re * = 12,000 and is characterized by a unique frequency for Re * = 14,000. As Re * increases, additional nonzero frequencies are observed in the spectra and the amplitudes of the corresponding signals also increase in accordance with previously reported results [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF]). As a consequence, the critical Reynolds number for the particular case of Fig. 4 can be identified to be Re * c = 13,000 ± 1000 (Re kc = 5.75 ± 0.26, Re dc = 28.05 ± 1.28). Critical Reynolds numbers for the different structures under consideration in this work are listed in Table 1. Values of Re dc and Re kc are computed from the values of Re * c after the simulation is completed using Eqs. ( 12) and ( 13). Several conclusions on the effect of the geometry on Re dc can be drawn. Firstly, as has been pointed out by [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF], the critical Reynolds number, Re dc , decreases drastically as more disorder is introduced into the pore structure. This can be observed from the values of Re dc obtained for all globally disordered structures (GWD, GWDf, GSD and GWDC) compared to those reported for locally disordered (LWD an LWDC) and ordered (OS and OC) media. Secondly, the effect of the size distribution of solid squares or circular obstacles is quite limited, as can be seen from the values of Re dc obtained for GWD and GSD structures. For this reason, critical values of the Reynolds number are not listed for GSDf and GSDC structures as they are expected to be similar to those obtained for GWDf and GWDC structures. Thirdly, unlike the study of [START_REF] Sohankar | On the use of a Darcy-Forchheimer like model for a macro-scale description of turbulence in porous media and its application to structured packings[END_REF], in which the angle of incidence of the flow with respect to the square obstacle has a great influence on the critical Reynolds number, there is little effect of the square inclinations as can be seen from Re dc values for GWD and GWDf structures. This contrast may be attributed to the fact that, in the present study, randomness of the angle of incidence from 0 • to 45 • applied on the squares may lead to strong interactions of the wakes behind neighboring obstacles which may weaken the effects of the local flow incidence. Fourthly, Re dc values for structures composed of circular obstacles are slightly larger than those for structures with solid squares, as can be observed from a comparison of Re dc results obtained for GWD, LWD and OS structures on the one hand and for GWDC, LWDC and OC structures on the other hand. Finally, although the observation is restricted to only two realizations, a significant dispersion of the critical Reynolds number over the different realizations of the same structure seems to occur. A more detailed analysis, however beyond the scope of this work, would be necessary to quantify this effect. It should be noted that the above conclusions drawn for Re dc cannot be straightforwardly extended to Re kc or Re * c , since for instance, Re d and Re k may have opposite variations, as can be seen from the results between GSD(2) and GWDC(2) or between OS and OC. Indeed, different definitions of the Reynolds number may lead to opposite conclusions as pointed out by [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF]. In the following sections, analysis of the inertial correction is restricted to Re k up to the lower bound value obtained for each structure for which steady flow is ensured.

Influence of the Microstructure

The microstructure has an impact not only on the permeability of the medium but also on the inertial regimes and the ranges of Re k on which they can be identified.

On the Permeability

The permeability tensor, which represents an important physical parameter for the process under study as it is further used to compute f c , is analyzed in this section. The attention is focussed on the mean diagonal term k * = locally disordered structures [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF]. Two major remarks must be made on these results. First, the extra diagonal term remains approximately 2 to 3 orders of magnitude smaller than the diagonal term, confirming that the structures are isotropic and justifying the use of k * to characterize their permeability. Secondly, as mentioned in the appendix, the standard deviation, σ , on the diagonal term is roughly 100 times smaller than k * . This confirms that the ensemble average provides an accurate representation of the infinite medium for each structure. In addition, results in Table 2 also show that the two GD structures containing circular cylinders (GWDC and GSDC) are more permeable than those containing square cylinders, confirming the same tendency observed on ordered structures [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF]. This result can be explained by the fact that local variations in the direction of flow streamlines are more abrupt in the presence of square obstacles than for circular inclusions, yielding a stronger resistance to flow. Moreover, k * for GD structures made up of square cylinders is slightly smaller than that obtained on corresponding less disordered geometries, namely OS, LWD and LSD. This results from the fact that the randomly distributed obstructions reduce the possibility of the formation of 'preferential channels,' confirming an observation reported elsewhere [START_REF] Clavier | Etude expérimentale et modélisation des pertes de pression lors du renoyage d'un lit de débris[END_REF]. As a final remark, one may observe that structures made up of square obstacles have a permeability which, in the increasing order, can be classified as follows: GWDf, GWD, GSDf, GSD, OS, LWD and LSD.

On the Flow Regimes and Inertia

In Fig. 5, the two components, f cx and f cy , of the correction vector to Darcy's law (see Eq. 11) obtained with ∇ p * β β = e x are plotted as functions of Re k for the GWD and GWDf structures. A similar behavior was obtained for the other structures and is not represented here for the sake of brevity. It must be noted that f cy ≈ 0 and is hence completely negligible compared to f cx , confirming that the structures remain isotropic in the inertial regime. As a consequence the analysis of the different flow regimes can be carried out on f cx only, in agreement with results for simpler geometries reported in the work of [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF]. The variations of f cx versus Re k are represented in Fig. 6 for the different GD structures under consideration along with results on ordered and LD structures reported earlier [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF]. Data points are average values over all the realizations for each structure.

The dependence of f cx on Re k allows to distinguish the so-called weak inertia and strong inertia regimes [START_REF] Skjetne | New insights on steady, non-linear flow in porous media[END_REF][START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF] for which linear correlations versus Re 2 k and Re k are applicable, respectively, corresponding to a cubic and quadratic dependence of the correction to the Darcy velocity in a dimensional form. This is highlighted in Fig. 7 representing these correlations on Reynolds number intervals where the Pearson coefficient, R 2 , remains ≥ 0.9999. Data points used for the fits and represented in this figure were obtained from an average of all the numerical results obtained on the ten realizations for each structure. Reynolds values intervals corresponding to the two regimes are indicated in Table 3 together with those obtained for LWD and ordered structures reported in [START_REF] Agnaou | Une étude numérique des écoulements mono et diphasique inertiels en milieux poreux[END_REF]. Obviously, the interval of the weak inertia regime shrinks as more disorder is introduced, confirming an observation already made [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF]. This may consistently explain why this regime is not easily identified during experimental observations on natural Following a nomenclature used in a previous work [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], namely

f cx = ARe 2 k + B in the weak inertia regime (14a) f cx = C Re k + D in the strong inertia regime (14b)
the above-mentioned linear fits yield the coefficients A, B, C and D in the two regimes as reported in Table 4. From these results, it must be noted that the values of B and D are extremely small and have a negligible contribution. More importantly, it is striking to observe that A and C are significantly larger for GD than for LD structures, equivalent to an increase of the inertia intensity, extending a conclusion reached earlier between LD and ordered structures [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF].

The Reynolds number interval corresponding to the transitional flow between weak and strong inertia can be characterized by a single Reynolds number, Re kt , which is the one minimizing the error between the weak and strong inertia correlations in Eqs. (14a) and (14b) hence given by The values of Re kt for the different structures are listed in Table 3. As can be seen from this table, for both media made up of circular and square inclusions, Re kt decreases as disorder is introduced from ordered to LD and, further, to GD structures, without any impact of weak or strong disorder in the latter case. Values of Re kt for GSDf and GWDf structures are intermediate between those for GWD (which is the same as for GSD) and LWD media. Moreover, Re kt is always larger for the structure made up of circular inclusions compared to the corresponding one including square obstacles.

Re kt = C 2 A (15)
Beyond the strong inertia regime, the dependency of f cx with respect to Re k can be fitted according to the following polynomial correlation [START_REF] Agnaou | Une étude numérique des écoulements mono et diphasique inertiels en milieux poreux[END_REF]; [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF]):

f cx = E Re 2 k + F Re k + G (16)
and an illustration of the corresponding coefficients is given in Table 5 for GWDf, GWDC and GSDC together with some recalls for ordered and LD structures reported in [START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF]. As can be seen from Table 3, the range of Reynolds numbers where this regime is observed on GD structures is significantly reduced compared to ordered and LD structures. In particular, this regime was not identified for Re k ≤ Re kc for GWD and GSD structures, which means the bifurcation just follows strong inertia in these two pore geometries. Moreover, it must be noted that the values of E remain quite small compared to F and G. Because the maximum value of Re k where the correlation in Eq. ( 16) applies is about 8.5, the contribution of the quadratic term in this correlation remains quite small with respect to the linear term, suggesting that a linear correlation, close to that in the strong inertia regime, is robust in the Reynolds number range beyond the strong inertia regime. Indeed, with a less restrictive constraint on the Pearson coefficient to fit the data, regime beyond strong inertia identified above would be merged into the strong inertia regime. This remark does not apply to the LWD or ordered structures for which the maximum Reynolds number in the regime beyond strong inertia is about 20 or more. These observations may be better highlighted while comparing the Re k dependence of ∂ f cx ∂ Re k , normalized by the corresponding maximum value, over the whole steady regime for the different structures under consideration. This is represented in Fig. 8, except for the GSD structures (see the justification below). For all structures, at low Reynolds numbers, a linear relationship is obtained in this graph, as expected in the weak inertia regime. The quadratic correction to Darcy's law in the strong inertia regime at larger values of Re k , corresponding to a classical Forchheimer correction widely used in the literature, is much more robust for GD (in particular for GWD and GWDC) (Fig. 8a) than for ordered or LD structures (see Fig 8b). As mentioned above, for the GWDf structure, no significant change with respect to the strong inertia regime is observed in the regime identified beyond strong inertia, confirming that a quadratic correction remains quite relevant in this Reynolds number range (3.98 ≤ Re k ≤ 5.59).

As final important remarks, it should be noted that:

-For ordered structures, inertia intensity is larger for media with circular inclusions (OC) than for those with square obstacles (OS) -Inertia is more intense for GD (either weakly or strongly disordered) structures made up of square obstacles than for those with circular inclusions. This might be physically explained by the fact that, since the wake downstream of a square cylinder (see illustrations in [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF]) is much wider than behind a circular cylinder [START_REF] Yuce | A numerical analysis of fluid flow around circular and square cylinders[END_REF][START_REF] Zaki | Numerical and experimental investigation of flow past a freely rotatable square cylinder[END_REF], interaction of the wakes in structures with squares become stronger than that with circular cylinders, thus increasing inertia intensity. This effect is not present in ordered structures as wakes are well separated in the direction orthogonal to the flow. -As can be seen from Table 4 and Figs. 6 and 7, inertia is not significantly altered by polydispersity for globally disordered media, i.e., GS or GW structures feature almost the same behavior. This contrasts with LD structures as previously observed in [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF]. This might be attributed to the fact that global disorder prevails over size distribution of the inclusions. From that, it may be concluded that one-phase flow analysis can be restricted to monodisperse GD structures. This justifies why attention was only focused on these types of media in Fig. 8. -Random inclinations of square obstacles leads to slightly larger inertial effects in GD structures. A similar phenomenon was observed for ordered structures (OS and OC) in the work of [START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF], in which an average pressure gradient oriented at an angle θ = 45 • with respect to the lattice vectors of the geometrical unit cell yields more intense inertia compared to the case where θ = 0 • . This phenomenon needs to be further investigated in random structures as flow topologies differ a lot between ordered and GD structures.

Conclusions

Effect of the microscopic structure on single-phase incompressible flow in 2D model porous media has been thoroughly investigated in this paper, widely extending previous analyses on this subject. Firstly, six types of spatially periodic porous structures were generated to introduce more structural disorder, compared to the ordered and locally disordered media studied so far in the literature [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF]). The critical Reynolds number at which laminar flow becomes unsteady was determined. On the basis of this critical value, numerical simulations of steady laminar flow were performed in order to thoroughly investigate the inertial correction to Darcy's law. It was found that inertial flow is significantly affected by the microstructure. The critical Reynolds number based on the obstacle size becomes much smaller as more disorder is introduced into the pore structures. In contrast with locally disordered structures considered in previous works, the critical Reynolds number for globally disordered structures significantly varies from one realization to another, suggesting that, for the size of the unit cell employed in this work, an ensemble average is required to capture a value representative of the same structure of infinite extent. The critical Reynolds number is significantly affected by the cylinder shape, but hardly depends on the size distribution and inclination of the cylinders.

The Reynolds number interval corresponding to the weak inertia regime shrinks when pore structures become more disordered. Compared to the weak inertia regime, the strong inertia regime is much more robust for GD structures than for less disordered structures studied in previous works. The intensity of inertia increases with structural disorder and is mainly influenced by the shape of the obstacles but is slightly affected by the inclination of the square cylinders and is barely disturbed by the size distribution of the obstacles.

While the present analysis has been dedicated to the homogeneous part of a porous medium, a lot remains to be done to investigate inertial effects close to the macroscopic boundaries of the system. A derivation of a macroscopic model valid everywhere, including in the neighborhood of the macroscopic boundaries, was proposed in [START_REF] Paéz-García | Macroscopic momentum and mechanical energy equations for incompressible single-phase flow in porous media[END_REF]. However, specific numerical simulations must be performed to determine the modifications of the inertial correction expected close to these boundaries.

1

  Arts et Métiers, CNRS, I2M, Esplanade des Arts et Métiers, 33405 Talence Cedex, France 2 CNRS, I2M, UMR 5295, Esplanade des Arts et Métiers, 33405 Talence Cedex, France

-

Fig. 1

 1 Fig. 1 Ordered (a, b) and locally disordered (c, d) pore structures[START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF]. OS, OC, LWD and LSD, respectively, refer to ordered squares, ordered circles, locally weakly disordered and locally strongly disordered structures.

Fig. 2

 2 Fig. 2 Examples of generated periodic unit cells of globally disordered pore structures including 900 solid obstacles

Fig. 3

 3 Fig. 3 Example of locally refined mesh

Fig. 4

 4 Fig. 4 Time evolution of the y-component of the average velocity (a, c) and corresponding Fourier spectra (b, d) for a GWDf realization. Re * = 12,000 (Re k = 5.48, Re d = 26.74) (a, b) and Re * = 14,000 (Re k = 6, Re d = 29.31) (c, d)

Fig. 5

 5 Fig. 5 Variations of the x-and y-components of f c with Re k for the ten realizations and their averages for the GWD and GWDf structures

Fig. 7

 7 Fig. 7 Variation of f cx versus Re 2 k (resp. Re k ) for the weak (resp. strong) inertia regime. Data points results from an average over the ten realizations for each structure

Fig. 8

 8 Fig. 8 Variation of ∂ f cx ∂ Re k , normalized by the corresponding maximum value, with Re k in the steady flow regime for GD structures with cylinders of uniform size (a) and comparison with ordered and locally disordered structures (b)

Table 1

 1 Critical Reynolds numbers (Re dc , Re kc , Re * c ) for GD structures and comparisons with results reported in previous works for ordered and LWD structures[START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF]

	Structure	Re dc	Re kc	Re * c
	GWD(1)	36.12 ± 1.26	7.67 ± 0.27	17,000 ± 1000
	GWD(2)	31.96 ± 1.37	6.79 ± 0.29	14,000 ± 1000
	GWDf(1)	35.7 ± 1.08	7.39 ± 0.22	19,000 ± 1000
	GWDf(2)	28.05 ± 1.28	5.75 ± 0.26	13,000 ± 1000
	GSD(1)	27.31 ± 1.95	6.28 ± 0.45	9000 ± 1000
	GSD(2)	37.11 ± 1.18	8.24 ± 0.26	18,000 ± 1000
	GWDC(1)	41.16 ± 2.28	8.57 ± 0.47	11,000 ± 1000
	GWDC(2)	38.21 ± 2.34	7.73 ± 0.47	11,000 ± 1000
	LWD (Agnaou et al. 2016)	8 3 .51 ± 2.45	19.15 ± 0.56	30,000
	LWDC (Agnaou et al. 2017)	93.46	19.87	30,000
	OS (Agnaou et al. 2016)	100.82 ± 2.5	2 3 .02 ± 0.57	17,250 ± 500
	OC (Agnaou et al. 2017)	101.39	21.09	17500

  * yx , the overbar indicating an average over the ten realizations. Values of these two quantities are reported in Table2together with some previously reported results obtained on ordered and

	k * xx +k *

yy 2 and on the extra diagonal term k * xy = k

Table 2

 2 Average values of the components of K * obtained over ten realizations of GD structures and comparisons with previous works[START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Agnaou | Origin of the inertial deviation from Darcy's law: an investigation from a microscopic flow analysis on two-dimensional model structures[END_REF]

	Structure	k * (σ )	k * xy (σ )
	GWD	0.01123 (3.7380e-4)	-1.50624e-5 (1.82847e-4)
	GWDf	0.01039 (2.5727e-4)	2.90673e-5 (2.05526e-4)
	GSD	0.01285 (4.71718e-4)	3.08100e-4 (2.30627e-4)
	GSDf	0.01231 (3.47639e-4)	8.133e-5 (2.28609e-4)
	GWDC	0.01367 (5.20535e-4)	3.85605e-4 (2.33563e-4)
	GSDC	0.01547 (4.4480e-4)	3.47906e-4 (2.334833e-4)
	LWD (Lasseux et al. 2011) a	0.01308 (3.82e-4)	-1.01e-4 (4.75e-4)
	LWD (Agnaou et al. 2017) b	0.01321
	LSD (Lasseux et al. 2011) a	0.01350 (7.72e-4)	-3.43e-4 (6.69e-4)
	LSD (Agnaou et al. 2017) b	0.01440
	OS (Lasseux et al. 2011)	0.01304
	OS (Agnaou et al. 2017)	0.01303
	OC (Agnaou et al. 2017)	0.01377

a Results were obtained with n = 100 and 10 realizations b Results were obtained with n = 900 and one realization

Table 3

 3 Re k intervals for the different flow regimes and Re kt values

	Structure	Weak (Re k )	Re kt	Strong (Re k )	Beyond (Re k )
	GWD 0.35 porous materials. The weak and strong inertia regimes are separated by a transitional flow [1.19e-5, 0.1]
	characterized again by a Reynolds number interval.

Table 4

 4 Values and standard errors (σ ) of the coefficients in correlations (14a) and (14b). Note that for LWD[START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF] and LSD[START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], coefficients are the average values obtained on the fits performed on 10 different realizations and σ is the standard deviation of all the values

	Structures	A(σ )	B(σ )	C(σ )	D(σ )
	GWD	0.291 (5.3e-4) -8.40e-5 (1.6e-6) 0.203 (4.5e-4) -0.06 (1.3e-3)
	GWDf	0.294 (6.9e-4) 4.71e-6 (2.5e-6)	0.248 (9.2e-4) -0.12 (3e-3)
	GSD	0.299 (7.3e-4) 6.57e-5 (3.3e-6)	0.210 (4.1e-4) -0.05 (1.2e-3)
	GSDf	0.309 (6.2e-4) 4.57e-6 (2.5e-6)	0.253 (8.8e-4) -0.12 (3.2e-3)
	GWDC	0.193 (4.1e-4) 4.30e-6 (2.3e-6)	0.160 (4.6e-4) -0.05 (1.2e-3)
	GSDC	0.198 (5.3e-4) 5.30e-6 (2.9e-6)	0.163 (5e-4)	-0.05 (1.3e-3)
	LWD (Lasseux et al. 2011) 0.070 (6.4e-3) 7.10e-6 (2.1e-6)	0.10 (8e-3)	-4.8e -2 (5.9e-3)
	LWD (Agnaou et al. 2017)	0.0683		0.0989
	LSD (Lasseux et al. 2011)	0.20 (4.9e-2) 6.8e -6 (1.6e-6)	0.15 (1.2e-2) -4.1e -2 (6.6e-3)
	LWDC (Agnaou et al. 2017) 0.0631		0.1042
	OS (Lasseux et al. 2011)	0.0103	1.3e-6	0.0153	-6.7e-3
	OS (Agnaou et al. 2017)	0.0102		0.0153
	OC (Agnaou et al. 2017)	0.0212		0.0411

Table 5

 5 Coefficients of the polynomial fit for the regime beyond strong inertia

	Structure	E(σ )	F(σ )	G(σ )

Acknowledgements

The financial support from the program of China Scholarship Council (No. 2015068300-49) and computational resources provided by the computing facility MCIA (Mésocentre de Calcul Intensif Aquitain) of the University of Bordeaux are gratefully acknowledged.

Appendix

In this appendix, the mesh convergence analysis is summarized. It was performed by evaluating the impact of the mesh size on k * = k * xx +k * yy 2 , k * x x and k * yy being the dimensionless diagonal permeability components averaged over the ten realizations, the off-diagonal components, k * xy , remaining negligible (see Sect. 3.2.1). The variation of k * versus the number of grid blocks is reported in Fig. 9. It must be noted that the standard deviation on diagonal term of the permeability is of the order of 10 -4 , a value which is roughly two orders of magnitude smaller than the typical value of k * , showing that each of the ten realizations provides a quite accurate value of the permeability (see Sect. 3.2.1 for additional comments).

Convergence was assessed by considering a relative variation of k * smaller than 0.5%, taking the reference value as the one obtained for the finest mesh size (i.e., ∼ 3.2 × 10 6 grid blocks). From this criterion, a mesh made up of around 2.2 × 10 6 grid blocks was found adequate for all the structures. To make sure this convergence remains for larger Reynolds numbers, the relative error on the x-component of f c , f cx , for these two grid block numbers was computed on a realization of a GWD structure for the maximum value of the Reynolds number considered for this structure (see Sect. 3.1), namely Re * = 18,000. The relative error was found to be around 0.22%, confirming the choice for the mesh size.


Six months after the launch of a national consultation, the Libault Report (issued in March 2019) finally paves the way for the muchneeded reform of the French long-term care system. The Report reflects the diverse opinions of the stakeholders concerned and contains 175 measures. It presents "loss of autonomy" as a social risk to be integrated into social security funding legislation and will be used as a background document for the legislative reform of long-term care policy scheduled for the end of 2019.

LEGAL NOTICE

This document has been prepared for the European Commission. However, it reflects the views only of the authors, and the Commission cannot be held responsible for any use which may be made of the information contained therein. The Report starts with a precise overview of policy targeted at elderly people. It highlights the high level of public funding -€23.7 billion, i.e. 1.2% of national wealth (Gross Domestic Product) -and identifies "an unclear system despite increasing public efforts and numerous innovative initiatives" (Press Pack, March 2019, p. 20). Reform is thus necessary to handle the predicted increase in public expenditure in favour of elderly people, evaluated at 1.6% of GDP in 2030. While maintaining older people in their homes and increasing freedom of choice in the organisation of care are put forward as priorities, the 175 measures presented in the Report concern both the development of quality support in the home and care in institutions. Concretely, three main strands can be identified. or identifying difficulties. In the absence of sufficient funding to implement the recommendations put forward in the Report, the main danger is that the country will slide into a public intervention policy that primarily focuses on support to family carers.

Future trends in

Further reading

  2019, six months after the launch of a major national consultation process to tackle the challenges of ageing, Dominique Libault, Councillor of State and President of the High Council of Financing for Social Protection, handed over his Report on autonomy and old age to Agnès Buzyn, Minister of Solidarity and Health. The Report outlines the foundations of "a strong, new policy on old age in France". It builds on a major national consultation carried out between September 2018 and February 2019, and involving all actors concerned with old age along with experts in the field.

  Libault, D., Grand âge, le temps d'agir, Concertation Grand âge et Autonomie, Ministère des solidarités et de la santé, March 2019. Grand âge, le temps d'agir, Concertation Grand âge et Autonomie, Dossier de presse, ministère des solidarités et de la santé, 28 March 2019. IGAS, 2017, Evaluation de la mise en oeuvre de la loi d'adaptation de la société au vieillissement pour le volet domicile, Tome 1, September 2017. Le Bihan, B., France anticipates ageing society through new piece of legislation, ESPN Flash Report 2016/18, April 2016.

  The first of these consists of a series of measures aimed at reorganising the various types of existing financial support, which include:

		allowances received have been
		deducted.	This	estimate	is
		questionable since it only takes
		into account services defined in the
		care plan (drawn up as part of the
		Personal		autonomy	allowance
		delivery procedure) and does not
		include varying service prices,
		additional services, or the cost of
		informal care.
		In sum, while the Libault Report is
		comprehensive,	arguments	in
		favour of high public investment
		should	not	be	limited	to
	BLANCHE LE BIHAN -EUROPEAN SOCIAL POLICY NETWORK	highlighting existing funding and MAY 2019 measures,
		 creating a new home-based cash
		benefit (to replace the current
		"Personal	autonomy	allowance"
		[Allocation			personnalisée
		d'autonomie])	with	three
		components: personal assistance,
		technical support and respite (for
		carers);
		 merging healthcare and social care
		expenditure in residential homes to
		reduce the remaining amount which
		is payable by residents, evaluated at
		€1,850.
		The second strand concerns the
		development of services in the home and
		in residential homes through:
		 renovating residential homes and
		making them more open to the
		outside world: a renovation plan
		worth	€3	billion	has	been
		announced;
		 increasing the supervision rate in
		residential homes by 25%, by
		recruiting 80,000 employees at an
		estimated cost of €1.2 billion;
		 creating 60,000 places in residential
		homes;
		 developing alternatives to residential
		homes or home-based care, i.e.
		temporary	accommodation	and
		collective housing;
		 creating indicators to measure the
		quality of services available in
		residential homes.

Author

Blanche Le Bihan (blanche.lebihan@ehesp.fr) and Alis Sopadzhiyan, Université de Rennes, Ecole des Hautes études en Santé Publique, Arènes (UMR 6051).

Finally, the third strand aims to increase support to the 3.9 million informal carers who provide care to older relatives. This will be done by:  simplifying procedures and access to information;  providing financial aid to support informal carers;

 facilitating informal carers' work/life balance.

The Report argues in favour of financing long-term care (LTC) policy through national solidarity, by recognising "loss of autonomy" as a genuine social protection risk, and including this in social security funding legislation. Recourse to private funding is presented as additional to public funding. One possibility could be to take account of a share of property assets when calculating the level of the benefit received, in order to support funding of home-based and residential care.

Outlook and commentary

In the light of demographic trends, the estimated number of elderly people eligible for the "Personal autonomy allowance" is likely to increase considerably in the coming decades: from 1,265,000 in 2015 to 1,582,000 in 2030 and 2,235,000 in 2050. The French LTC policy reform is therefore urgently needed. Announced in the early 2000s, it was launched with the 2015 Act on Adapting Society to an Ageing Population, but with a budget considered to be highly insufficient (Le Bihan 2016, IGAS 2017). The Libault Report continues in the vein of this legislation by aiming to change the negative representation of population ageing by developing a more inclusive society that fully integrates older people. Commissioned by the Minister of Solidarity and Health herself, it puts a strong emphasis on the involvement of public authoritiesat both central and local levels -in the loss of autonomy sector. The Report is a response to the current crisis in residential homes and to requests for support from families; it is thus generally the object of consensus.

The Report is conceived as a background document in preparation for legislation to reform LTC policy, scheduled for the end of 2019. This reform will have to address three key challenges.

First, the set of measures announced in the Report will require massive public funding; although difficult to evaluate, the estimated additional amount is €9.2 billion by 2030. The favoured scenario is that an existing pay deduction (the Contribution to Reducing the Social Debt [CRDS]), which will have been fully paid by 2024, will be converted into funding for loss of autonomy. Nevertheless, political and social developments due to the "yellow jacket" movement and multiple costly reform projects resulting from the "Great Debate" organised in response to this movement could call into question the availability of this funding. Alternative scenarios, such as creating new mandatory pay deductions, or extending working time by cancelling a national holiday, have so far been ruled out.

Secondly, the governance of LTC policy needs to be clarified, as it involves two main institutional actors -regional health agencies and local councils -relations between which are currently very tense.

Thirdly, the declared priority of maintaining elderly people in their home is based on the assumption that the average amount to be paid by the care user (and/or his/her family) to stay in the home is €60 per month once possible