domain: math spi

text:
Jean Auriol

Federico Bribiesca-Argomedo

Silviu-Iulian Niculescu

Jeanne Redaud

Stabilization of a hyperbolic PDEs-ODE network using a recursive dynamics interconnection framework

In this paper, we design a state-feedback law that exponentially stabilizes an underactuated network of scalar hyperbolic systems coupled at the end of the chain with a finite-dimensional system modeled by an ordinary differential equation. Our approach uses a recursive dynamics interconnection framework. More precisely, for each subsystem, we solve a stabilization problem, a tracking problem and design a statepredictor. Then, we are able to recursively design a control law that stabilizes the whole chain. Some illustrative examples complete the presentation.

I. INTRODUCTION

Interconnected systems represented by Partial Differential Equations (PDEs) naturally appear when modeling industrial processes. Some well-known examples include traffic network systems with different types of vehicles [START_REF] Burkhardt | Stop-and-go suppression in twoclass congested traffic[END_REF], [START_REF] Yu | Simultaneous stabilization of traffic flow on two connected roads[END_REF], ventilation in buildings [START_REF] Rodriguez-Diaz | Control strategies for ventilation networks in smallscale mines using an experimental benchmark[END_REF], density-flow systems [START_REF] Bastin | Stability of linear densityflow hyperbolic systems under PI boundary control[END_REF], [START_REF] Hayat | Exponential stability of density-velocity systems with boundary conditions and source term for the h 2 norm[END_REF], open channels [START_REF] Gu | Exact boundary controllability of nodal profile for quasilinear hyperbolic systems in a tree-like network[END_REF], [START_REF] Gu | Exact boundary controllability of nodal profile for unsteady flows on a tree-like network of open canals[END_REF], [START_REF] Hayat | PI controller for the general saint-venant equations[END_REF], communication networks [START_REF] Espitia | Fluid-flow modeling and stability analysis of communication networks[END_REF]. Additional couplings with Ordinary Differential Equations (ODEs) arise when there is a lumped element coupled to the distributed dynamics such as in heavy chain systems [START_REF] Petit | Flatness of heavy chain systems[END_REF], Rijke tubes [START_REF] De Andrade | Backstepping stabilization of a linearized ODE-PDE rijke tube model[END_REF] (where the ODE is sandwiched between two PDEs) or mechanical vibrations in drilling applications [START_REF] Saldivar | A control oriented guided tour in oilwell drilling vibration modeling[END_REF]. In the last case case, the hyperbolic PDE system models axial and torsional stress propagation (with potential discontinuities due to the junctions between different types of pipes), while the ODE represents the dynamics of the Bottom hole Assembly (BHA). Several approaches have been proposed in the literature to stabilize such interconnected systems. Among them, we can cite PI boundary controllers in the case of fully actuated networks [START_REF] Bastin | Stability of linear densityflow hyperbolic systems under PI boundary control[END_REF], flatness-based designs [START_REF] Schmuck | Feed-forward control of an HVDC power transmission network[END_REF], or the method of characteristics for quasilinear hyperbolic systems in a treelike networks [START_REF] Gu | Exact boundary controllability of nodal profile for unsteady flows on a tree-like network of open canals[END_REF], [START_REF] Wang | Exact boundary controllability of nodal profile for quasilinear wave equations in a planar tree-like network of strings[END_REF]. Recently, this last approach has been used in [START_REF] Strecker | Output feedback boundary control of series interconnections of 2× 2 semilinear hyperbolic systems[END_REF] to design output-feedback laws for networks of interconnected scalar semilinear PDEs. However, the proposed state-feedback control law requires solving a set of PDEs online which is computationally expensive. Regarding networks of ODEs and PDEs, the backstepping approach [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF] has been shown to be a powerful tool for the design of explicit stabilizing control laws. It was first applied to hyperbolic equations in [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF] to express the delays acting on an ODE as PDE states of transport equations. This led to a new interpretation of the delay-compensating Finite Spectrum Assignment controller. Based on extensions of the backstepping technique, more complex stabilization problems have been progressively solved: non-linear ODEs with delays [START_REF] Bekiaris-Liberis | Nonlinear control under nonconstant delays[END_REF], ODEs coupled with a beam [START_REF] Wu | Static output feedback control via pde boundary and ode measurements in linear cascaded ode-beam systems[END_REF], or velocity ing author: jean.auriol@centralesupelec.fr 2 Federico Bribiesca-Argomedo is with Univ Lyon, INSA Lyon, CNRS, Ampère, F-69621, Villeurbanne, France. recirculation in a wave equation [START_REF] Su | Boundary stabilization of wave equation with velocity recirculation[END_REF]. The backstepping approach has been successfully used in [START_REF] Auriol | Delay-robust stabilization of an n + m hyperbolic pde-ode system[END_REF], to stabilize a system of hyperbolic PDEs coupled with an ODE in a PDE-ODE structure. In the case of a scalar PDE system, [START_REF] Irscheid | Trajectory tracking control for a class of 2× 2 hyperbolic pde-ode systems[END_REF] solved a similar problem (with a potentially non-linear ODE) by solving a tracking problem. The ODE-PDE-ODE structure has been considered in [START_REF] Wang | Control of a 2× 2 coupled linear hyperbolic system sandwiched between 2 odes[END_REF], [START_REF] Deutscher | Output feedback control of general linear heterodirectional hyperbolic ode-pde-ode systems[END_REF], [START_REF] Gehring | A systematic design of backstepping-based state feedback controllers for ode-pde-ode systems[END_REF] and recently in [START_REF] Bou Saba | Strictly proper control design for the stabilization of 2 × 2 linear hyperbolic ODE-PDE-ODE systems[END_REF] with the design of a strictly proper control design. In recent contributions [START_REF] Aarsnes | Delay robust control design of under-actuated PDE-ODE-PDE systems[END_REF], the authors have stabilized systems with a PDE-ODE-PDE structure. Finally, the general problem of the stabilization of an arbitrary number of interconnected scalar PDEs has been solved in [START_REF] Auriol | Output feedback stabilization of an underactuated cascade network of interconnected linear pde systems using a backstepping approach[END_REF], introducing a complex backstepping transformation.

So far, one major limitation of the backstepping approach for the stabilization of general networks of interconnected ODE-PDE systems is its lack of genericity, i.e., the involved backstepping transformations are specific to the system under consideration. More precisely, adding a new subsystem into the network implies designing a completely new backstepping transformation to encompass this new subsystem. In other research fields, such as electronics [START_REF] Araujo | A modular approach for modeling and simulation of semiconductor power devices[END_REF] or bioengineering [START_REF] Simon | Peptoids: A modular approach to drug discovery[END_REF], [START_REF] Sunder | Hyperbranched polyether polyols: A modular approach to complex polymer architectures[END_REF], as well as in the process industry [START_REF] Ladiges | Integration of modular process units into process control systems[END_REF], modular approaches have been developed in order to efficiently model complex interconnected systems.

In this paper, we develop a new recursive interconnected dynamics framework to design an explicit state-feedback control law that stabilizes a network of interconnected scalar linear hyperbolic systems coupled with an ODE at the end of the chain. Roughly speaking, the proposed control law is recursively obtained by considering stabilizing virtual inputs for each subsystems and by ensuring the outputs of each subsystem converge to these desired virtual inputs. The control design becomes simpler and is based on simple assumptions that can be independently verified for each subsystem. This new framework allows for a "plug-and-play"-like approach to control design since additional subsystems, satisfying similar conditions, can be added to the network using the same procedure.

The approach followed in this paper consists on solving three independent problems for each subsystem: (i) a stabilization problem, (ii) a delayed-trajectory tracking problem and (iii) building a predictor for the system. If we consider the specific case of the last subsystem (the one furthest in the chain from the control action), it means in particular that we know an input that stabilizes it. Then we are able, using the second-to last subsystem, to track this desired signal (with inevitable delay due to finite speeds of propagation, which is compensated using the predictor of the last subsystem). Iterating this procedure, at each step we know a control input that solves the stabilization problem for all downstream subsystems and we use the current subsystem to track this trajectory, compensating the inevitable delay using the predictor for the downstream subsystem, which is constructed using a composition of the subsystem predictors.

This paper is structured as follows: in Section II, the control problem is presented, as well as the general ideas behind our new recursive dynamics interconnection framework. In Section III, we present necessary conditions to apply this framework and show that for the class of systems under consideration, they are always satisfied. We design the stabilizing control law in Section IV and present some simulation results in Section V.

II. PROBLEM DESCRIPTION A. System under consideration

We consider in this paper a system composed of a chain of n > 1 subsystems interconnected through their boundaries. Each subsystem consists of two coupled linear hyperbolic PDEs, except for the last one that can be an ODE. This class of system may appear when considering oil production systems made of networks of pipes (whose principal line is known as the manifold) [START_REF] Jadid | Performance evaluation of virtual flow metering models and its application to metering backup and production allocation[END_REF]. More precisely, the lower part of the drill-string is usually made up of drill collars that can have a great impact on the global dynamics due to their inertia [START_REF] Aarsnes | Linear stability analysis of selfexcited vibrations in drilling using an infinite dimensional model[END_REF]. In particular, these pipes may have different lengths, density, inertia or Young's modulus. These spatial variations in the characteristic line impedance may cause reflections to appear in the junctions. Each subsystem (except for the last one) is defined by the following set of

PDEs (i ∈ {1, • • • , n -1}) ∂ t u i (t, x) + λ i ∂ x u i (t, x) = σ + i (x)v i (t, x), (1)
∂ t v i (t, x) -µ i ∂ x v i (t, x) = σ - i (x)u i (t, x), (2)
evolving in {(t, x) s.t. t > 0, x ∈ [0, 1]}, where λ i > 0 and µ i > 0 are the constant transport velocities respectively associated to equations ( 1) and (2), and where the in-domain coupling terms σ + i and σ - i are continuous functions. The associated initial conditions are denoted u 0 i (•) = u i (0, •) and v 0 i (•) = v i (0, •) and are defined in L 2 ([0, 1], R). These subsystems are connected through their boundary conditions, which satisfy

u i (t, 0) = q i,i v i (t, 0) + q i-1,i u i-1 (t, 1) + δ i 1 V (t), (3)
and

v i (t, 1) = ρ i,i u i (t, 1) + ρ i,i+1 v i+1 (t, 0) + δ i n-1 CX(t) (4
) where the different coupling and ρ ij are constant and where V is a real-valued control input. By convention q 0,1 = 0 and ρ n-1,n = 0. The notation δ i j stands for the Kronecker symbol (i.e. δ i j = 1 if i = j, and δ i j = 0 if i = j). The last subsystem is an ODE defined by

Ẋ(t) = AX(t) + Bu n-1 (t, 1), (5)
where the matrices

A ∈ R p×p , B ∈ R p×1 , C ∈ R 1×p are constant. The initial condition X 0 = X(0) belongs to R p .
The well-posedness of system (1)-( 5) is guaranteed by [9, Theorem A.6, page 254]. The system (1)-( 5) is schematically pictured in Figure 1. The advantage of such a representation is that it highlights that the interaction between the different subsystems only occurs at the boundaries (black arrows). For instance, the output of the second subsystem entering the first subsystem could be seen as a kind of disturbance acting on the first subsystem (even if such a disturbance signal is indirectly modified by the corresponding output of the first subsystem). For a subsystem i, we will call the subsystem i + 1 the downstream subsystem and i -1 the upstream subsystem.

V (t) System 1 System 2 System n

Fig. 1. Schematic representation of the system under consideration (1)-( 5)

B. Problem formulation

The objective of this paper is to exponentially stabilize system (1)-( 5) using the available control V (t). The stabilization is done in the sense of the Ξ-norm defined by

||(u, v, X)|| Ξ = ( n-1 i=1 ||u i || 2 L 2 + ||v i || 2 L 2 + ||X|| 2 R p ) 1 2
. This problem has already been solved in [START_REF] Auriol | Output feedback stabilization of an underactuated cascade network of interconnected linear pde systems using a backstepping approach[END_REF] for the case of a chain of PDEs (no ODE at the end of the chain). The method presented in [START_REF] Auriol | Output feedback stabilization of an underactuated cascade network of interconnected linear pde systems using a backstepping approach[END_REF] uses a triangular integral transformation to move most of the in-domain couplings at the actuated boundary. We believe that such an approach could be adjusted to the presence of the ODE. However, this would require deriving and solving a new set of complex kernel equations. Generally speaking, the approach proposed in [START_REF] Auriol | Output feedback stabilization of an underactuated cascade network of interconnected linear pde systems using a backstepping approach[END_REF] cannot be straightforwardly extended when adding new blocks inside the chain that have a different structure (i.e. adding ODEs in the chain or PDEs that have a different number of states). This explains why here, we develop a new recursive dynamics interconnection framework. Based on simple conditions placed on each independent subsystem, we derive a recursive control law that stabilizes the interconnected system (1)-( 5). The main advantage of the proposed methodology is its potential extension to arbitrary chains of ODEs and PDEs. Moreover, the conditions derived in this paper are less restrictive than those used in [START_REF] Auriol | Output feedback stabilization of an underactuated cascade network of interconnected linear pde systems using a backstepping approach[END_REF].

C. Description of the method

The recursive dynamics interconnection framework we develop in this paper is based on a recursive approach. Considering the last subsystem as an independent system, we design a theoretical input that stabilizes it. This input is a state-feedback law that only depends on the last subsystem's state. Then, we consider the penultimate subsystem as an independent system only disturbed by the downstream subsystem's action. This penultimate subsystem acts on the last subsystem through its right output. For this system, we can design a (virtual) control input that stabilizes the system in the absence of the disturbance and, more importantly, that guarantees the convergence of its right output to the theoretical control law that stabilizes the last subsystem. Such a control law guarantees the stabilization of the last subsystem and finally of the penultimate subsystem. Due to the hyperbolic systems' transport structure, such a tracking control law may require a prediction of the downstream subsystem's state. Iterating such a procedure, it is possible to design a stabilizing control law V (t) for the whole system.

One of the main advantages of such an approach is that, provided we can solve each subsystem's tracking problem, and correctly predict each subsystem's state (which is possible for ODEs and hyperbolic PDEs), it becomes straightforward to add new subsystems. This results in a simpler and more generic control design procedure than the one developed in [START_REF] Auriol | Output feedback stabilization of an underactuated cascade network of interconnected linear pde systems using a backstepping approach[END_REF].

III. CONDITIONS TO DEVELOP A RECURSIVE DYNAMICS INTERCONNECTION FRAMEWORK

In this section, we introduce several assumptions that are necessary to apply our new control framework. They correspond to a set of simple and natural conditions defined for each subsystem of (1)- [START_REF] Auriol | Delay-robust stabilization of an n + m hyperbolic pde-ode system[END_REF]. We show that these conditions are always satisfied for the class of systems we consider. In the next section, we will show that under these assumptions it becomes straightforward to stabilize the system (1)- [START_REF] Auriol | Delay-robust stabilization of an n + m hyperbolic pde-ode system[END_REF]. We start giving the following definition to qualify the interactions between the different subsystems.

Definition 1: For a subsystem i < n, we denote Vi (t) the action of the (i -1) th subsystem on this subsystem. This function will be called virtual input acting on the i th subsystem. We also denote χ i (t) the action of the (i + 1) th subsystem on the i th subsystem. This function will be called virtual disturbance acting on the i th subsystem. This terminology is justified by the recursive dynamics interconnection framework we develop since each subsystem is stabilized using its leftward neighbor. With these notations, we have for instance

V2 (t) = q 2,1 u 1 (t, 1), χ 2 (t) = ρ 2,3 v 3 (t, 0), χ n-1 (t) = CX(t).
For the last subsystem (ODE) we do not have any virtual output or virtual disturbance and we denote Vn (t) = u n-1 (t, 1). These different notations are summarized in Figure 2. We are now ready to give the different assumptions for each subsystem and show that they are always satisfied.

System i

Vi (t)

χ i-1 (t) Vi+1 (t) χ i (t)

A. Delay-robustness assumption

The first assumption we make aims to guarantee the existence of robustness margins for the closed-loop system.

Assumption 1: The open loop system (1)-(5) (i.e. V (t) ≡ 0) in the absence of in-domain coupling terms (i.e. σ • i ≡ 0) and of the ODE (i.e X ≡ 0) is exponentially stable (in the sense of the L 2 norm). It has been shown in [START_REF] Logemann | Conditions for robustness and nonrobustness of the stability of feedback systems with respect to small delays in the feedback loop[END_REF] that a necessary condition to guarantee the existence of robustness margins for an arbitrary closed-loop system is that the open-loop transfer function must have a finite number of poles on the closed right halfplane. For a deeper discussion of the sensitivity wrt the delay parameter, the reader is referred to [START_REF] Michiels | Stability and stabilization of timedelay systems: an eigenvalue-based approach[END_REF]. It is worth mentioning that for the system (1)-( 5), [START_REF] Auriol | Delay-robust stabilization of an n + m hyperbolic pde-ode system[END_REF], [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF] proved that this implies that the open-loop system without in-domain couplings and without the ODE must be exponentially stable. Some explicit conditions to verify such an assumption can be found in [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF] and [START_REF] Auriol | Output feedback stabilization of an underactuated cascade network of interconnected linear pde systems using a backstepping approach[END_REF]. Consequently, Assumption 1 is not restrictive as it is necessary for the existence of robustness margins for the closed-loop system.

B. Stabilizability of each subsystem

The second assumption we make guarantees the possibility to independently stabilize each subsystem.

Assumption 2: In the absence of the virtual disturbance χ i (i.e. χ i (t) ≡ 0), the i th subsystem subject to the virtual actuation Vi must be stabilizable by a state-feedback law. More precisely, if i < n (resp. i = n), there exist a linear operator K i : (L 2 ([0, 1])) 2 → R (resp. K n : R p → R) such that Vi = K i (u i , v i ) (resp. Vn = K n (X)) stabilizes the system in the absence of the virtual disturbance. This assumption is natural, since to stabilize our interconnected system for any initial condition, we need to have each subsystem independently stabilizable. For the ODE to be stabilizable, we must have the existence of a matrix K such that A + BK is Hurwitz. (i.e. has all its eigenvalues in the complex left half plane). Indeed, if Vn (t) = KX(t), then the ODE system is stabilized. Let us now consider the i thsubsystem (i ∈ {1, 2, • • • , n -1}) with the virtual actuation Vi (t). Assuming that χ i (t) ≡ 0, this PDE system rewrites

∂ t u i (t, x) + λ i ∂ x u i (t, x) = σ + i (x)v i (t, x) (6) ∂ t v i (t, x) -µ i ∂ x v i (t, x) = σ - i (x)u i (t, x), (7)
with the boundary conditions

u i (t, 0) = q i,i v i (t, 0) + Vi (t), (8)
v i (t, 1) = ρ i,i u i (t, 1). ( 9
)
Let us consider the following backstepping transformation adjusted from [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF], [START_REF] Coron | Local exponential H 2 stabilization of a 2× 2 quasilinear hyperbolic system using backstepping[END_REF], [START_REF] Vazquez | Local exponential H 2 stabilization of a 2× 2 quasilinear hyperbolic system using backstepping[END_REF] and defined as

α i (t, x) β i (t, x) = u i (t, x) v i (t, x) - 1 x K uu i (x, y) K uv i (x, y) K vu i (x, y) K vv i (x, y) u i (t, y) v i (t, y) dy, (10)
where the kernels K •• i are L ∞ functions defined on the upper triangular part of the unit square T = {(x, y) ∈ [0, 1] 2 , x ≤ y}. They satisfy the following set of equations

λ i ∂ x K uu i (x, y) + λ i ∂ y K uu i (x, y) = -σ - i (y)K uv i (x, y), (11)
λ i ∂ x K uv i (x, y) -µ i ∂ y K uv i (x, y) = -σ + i (y)K uu i (x, y), (12)
µ i ∂ x K vu i (x, y) -λ i ∂ y K vu (x, y) = σ - i (y)K uv i (x, y), ( 13
) µ i ∂ x K vv i (x, y) + µ i ∂ y K vv (x, y) = σ + i (y)K uv i (x, y), (14)
K uv i (x, x) = - σ + i λ i + µ i , K uv i (x, x) = σ - i λ i + µ i , (15)
K uu i (x, 1) = ρ i,i µ i λ i K uv i (x, 1), K vv (x, 1) = 0. ( 16
)
It has been proved in [START_REF] Coron | Local exponential H 2 stabilization of a 2× 2 quasilinear hyperbolic system using backstepping[END_REF], [START_REF] Vazquez | Local exponential H 2 stabilization of a 2× 2 quasilinear hyperbolic system using backstepping[END_REF] that this set of equations admits a unique solution. Defining •)) (and therefore defining the operators K i in Assumption 2) as

V BS i . = K i (u i (t, •), v i (t,
V BS i (t) = -q i,i v i (t, 0) + 1 0 K uu i (0, y)u i (t, y)dy + 1 0 K uv i (0, y)v i (t, y)dy, (17)
the transformation (10) maps the system (6)-( 9) to the target system [START_REF] Espitia | Fluid-flow modeling and stability analysis of communication networks[END_REF] with the boundary conditions

∂ t α i (t, x) + λ i ∂ x α i (t, x) = 0, ( 18
) ∂ t β i (t, x) -µ i ∂ x β i (t, x) = λ i K vu (x, 1)ρ i α i (t, 1),
α i (t, 0) = 0, β i (t, 1) = ρ i,i α i (t, 1). ( 20
)
This target system is exponentially stable. This in turns implies the exponential stability of the original system due to the invertibility of the transformation [START_REF] Bastin | Stability of linear densityflow hyperbolic systems under PI boundary control[END_REF]. However, it is worth mentioning that the control law V BS i (t) cancels the boundary reflection terms. As shown in [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF], this may lead to vanishing delay margins. To avoid this problem and make the control law strictly proper, we can combine it with a well-tuned low pass filter as done in [START_REF] Auriol | Stabilization of an underactuated 1+ 2 linear hyperbolic system with a proper control[END_REF].

C. Output tracking in presence of a disturbance signal

The proposed approach requires each subsystem to be able to track any arbitrary function. More precisely, we have the following assumption Assumption 3: Consider a subsystem i < n and define Φ i ∈ L 2 loc (R + ) an arbitrary known function. Let us assume that the virtual disturbance χ i acting on this subsystem is known. Then, there exists a control law Vi that exponentially tracks the function Vi+1 (t) to the desired function Φ i (t). Moreover, if χ i (t) ≡ Φ i (t) ≡ 0, then, such a control law stabilizes the i th subsystem. More precisely, there exist linear causal operators

L i , such that Vi (t) = K i (u i (t, •), v i (t, •)) + L i (Φ i (t + 1 λi ), χ i ([t, t + 1 λi ]
)), with the operators K i being defined in Assumption 2 and ( 17), and the operators L i defined later in this section.

The function Vi+1 (t) (virtual input of the (i+1) th subsystem) corresponds to the right output of the subsystem i. The fact that the operators L i require future values of the Φ i and χ i is related to the transport delay to go from the left boundary of the PDE (where is the actuation) to the right boundary (where is defined the output we want to track). This assumption does not have to be fulfilled for the last subsystem (i.e. for the ODE), as explained in Section II-C. Since the function Vi+1 is defined by Vi+1 (t) = q i,i+1 u i (t, 1) (with the convention q n-1,n = 1), we are led to the following condition: for all 1 < i ≤ n, we have q i-1,i = 0. Indeed, if it is not the case, we cannot track any arbitrary function since the output will always be zero. Regarding our class of PDEs subsystems we can adjust the approach presented in [START_REF] Hu | Control of homodirectional and general heterodirectional linear coupled hyperbolic PDEs[END_REF] to solve this tracking problem. More precisely, in the presence of a disturbance χ i (t) acting on the i th system, the transformation (10) maps the system ( 6)-( 9) to the target system

∂ t α i + λ i ∂ x α i = -µ i K uv i (x, 1)χ i (t) (21) ∂ t β i -µ i ∂ x β i = λ i K vu i (x, 1)α i (t, 1), (22)
with the boundary conditions

α i (t, 0) = V tr i (t), (23)
β i (t, 1) = ρ i,i α i (t, 1) + χ i (t). ( 24
)
where we have defined

V tr i (t) = Vi (t) -V BS i (t)
. The backstepping transformation [START_REF] Bastin | Stability of linear densityflow hyperbolic systems under PI boundary control[END_REF] does not modify the function u i (t, 1). Indeed, we have Vi+1 (t) = q i+1,i u i (t, 1) = q i+1,i α i (t, 1) (with the convention q n,n-1 = 1). Applying the method of characteristics on equation [START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF], we have

α i (t, 1) = V tr i (t - 1 λ i ) -µ i 1 λ i 0 K uv i (x -λ i s, 1)χ i (t -s)ds. (25)
Thus, we can simply choose

V tr i . = L i (Φ i (t + 1 λi ), χ i ([t, t + 1 λi ]
)) (fixing the operators L i in Assumption 3) as

V tr i (t) = 1 q i+1,i Φ i (t + 1 λ i ) + µ i 1 λ i 0 K uv i (1 -λ i s, 1)χ i (t + 1 λ i -s)ds, (26)
to guarantee the convergence of y i (t) = u i (t, 1) = α i (t, 1) to Φ i (t). If the functions Φ i and χ i converge to zero, so does the state α i (t, 1) and so does the whole i th -subsystem due to its transport structure.

D. Prediction of the state

The tracking operators L i require future values of the functions Φ i and χ i . Since these functions depend on the states of the downstream subsystems, we must be able to predict future values of these states. This leads us to the following assumption.

Assumption 4: For all x ∈ [0, 1], it is possible to obtain a i-1 j=1

1 λj + x
λi units of time ahead prediction of the PDE states u i (t, x) and v i (t, x), and a n-1 j=1

1 λj + x
λi units of time ahead prediction of the ODE X. More precisely there exist predictor functions P vi , P ui and P X such that for all t > 0, P vi (t, x) = v i (t + i-1 j=1

1 λj + x λi , x), P ui (t, x) = u i (t + i-1 j=1 1 λj + x λi , x), P X (t) = X(t + n-1 j=1 1 λj
). Let us explicitly design these predictors. In what follows, to avoid useless case disjunctions, we use the convention ρ n-1,n = C and v n (t, 0) = X. Applying the transformation [START_REF] Bastin | Stability of linear densityflow hyperbolic systems under PI boundary control[END_REF] to each subsystem ( 6)-( 9), we obtain the target system

∂ t α i (t, x) + λ i ∂ x α i (t, x) =f + i (x)v i+1 (t, 0) (27) ∂ t β i (t, x) -µ i ∂ x β i (t, x) =f - i (x)u i (t, 1), (28)
with the boundary conditions

α i (t, 0) = q i,i v i (t, 0) + q i,i-1 u i-1 (t, 1) + δ i 1 V (t) + 1 0 L αα i α i (t, ξ) + L αβ i β i (t, ξ), (29)
β i (t, 1) = ρ i,i α i (t, 1) + ρ i,i+1 v i+1 (t, 0) + δ i n-1 CX. ( 30
)
where we recall that the kernels L •• i are the inverse kernels of the kernels K •• i (see [START_REF] Coron | Local exponential H 2 stabilization of a 2× 2 quasilinear hyperbolic system using backstepping[END_REF], [START_REF] Vazquez | Local exponential H 2 stabilization of a 2× 2 quasilinear hyperbolic system using backstepping[END_REF] for details), and where we have defined

f - i (x) = λ i K vu i (x, 1) and f + i (x) = -µ i ρ i,i+1 K uv i (x, 1
). With these definitions, equations ( 27)-( 28) correspond to the tracking system ( 21)- [START_REF] Gu | Exact boundary controllability of nodal profile for quasilinear hyperbolic systems in a tree-like network[END_REF]. To these PDEs, we add the ODE [START_REF] Auriol | Delay-robust stabilization of an n + m hyperbolic pde-ode system[END_REF]. Note that we decided to preserve the terms u i (t, 1) and v i (t, 0) inside the target systems to avoid adding integral terms in α and β. We will first design a predictor for the states u i (t, 1), v i (t, 0), α i (t, 0) and X(t) and show that these predictors give us the wanted predictor for the whole original states. Using the method of characteristics, we have for i < n + 1

u i (t, 1) =α i (t - 1 λ i , 0) + 1 λ i 0 f + i (1 -νλ i )v i+1 (t -ν, 0)dν, (31)
(for i = n -1, the term v i+1 (t, 0) has to be replaced by X(t)). Following the approach given in [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF] (we choose to skip the complete proof here, due to space restrictions), we obtain

v i (t, 0) = ρ i,i u i (t - 1 µ i , 1) + ρ i,i+1 v i+1 (t - 1 µ i , 0) + τi 0 g 1 i (ν) α i (t -ν, 0) + g 2 i (ν)u i (t -ν, 1) + g 3 i (ν)v i+1 (t -ν, 0)dν. (32)
where τ i = 1 λi + 1 µi and where the functions g 1 i , g 2 i and g 3 i are defined by

g 1 i (ν) = -1 [0, 1 λ i ] (ν)λ i L βα i (0, λ i ν), (33)
g 2 i (ν) = 1 [0, 1 µ i ] (ν)(f - i (µ i ν) -µ i ρ i,i L ββ i (0, 1 -µ i ν) - 1-µiν 0 L ββ i (0, ξ)f - i (ξ + µ i ν)dξ) (34)
g 3 i (ν) = -1 [0, 1 µ i ] (ν)(µ i ρ i,i+1 L ββ i (0, 1 -µ i ν)) -1 [0, 1 λ i ] (ν)( 1 λiν L βα i (0, ξ)f + i (ξ -λ i ν)dξ), (35)
where 1 Ω denotes the characteristic function of the set Ω. Similarly, we obtain

α i (t, 0) = q i,i v i (t, 0) + q i,i-1 u i-1 (t, 1) + δ 1 i V (t) + τi 0 k 1 i (ν)α i (t -ν, 0) + k 2 i (ν)u i (t -ν, 1)dν + τi 0 k 3 i (ν)v i+1 (t -ν, 0)dν. (36)
where the functions k 1 i , k 2 i and k 3 i are defined by

k 1 i (ν) = 1 [0, 1 λ i ] (ν)λ i L αα i (0, λ i ν), (37)
k 2 i (ν) = 1 [0, 1 µ i ] (ν)(µ i ρ i,i L αβ i (0, 1 -µ i ν) + 1-µiν 0 L αβ i (0, ξ)f - i (ξ + µ i ν)dξ) (38)
k 3 i (ν) = 1 [0, 1 µ i ] (ν)(µ i ρ i,i+1 L αβ i (0, 1 -µ i ν)) + 1 [0, 1 λ i ] (ν)( 1 λiν L αα i (0, ξ)f + i (ξ -λ i ν)dξ). (39)
We can now define the state predictors P αi (t, s), P vi (t, s), P vi (t, s) and P X (t, s) (see [START_REF] Bekiaris-Liberis | Simultaneous compensation of input and state delays for nonlinear systems[END_REF], [START_REF] Bresch-Pietri | Prediction-based control of linear input-delay system subject to state-dependent state delay-application to suppression of mechanical vibrations in drilling[END_REF]). Let us define τ = i τ i as the total transport delay. Inspired by [START_REF] Auriol | Stabilization of an underactuated 1+ 2 linear hyperbolic system with a proper control[END_REF], let us define for t ≥ 0 and s ∈ [t -τ -i-1 j=1 1 λj , t],P αi (t, s) (resp. P vi (t, s)) as the state prediction of α i (t, 0), (resp. v i (t, 0)), ahead a time 40)- [START_REF] Wang | Control of a 2× 2 coupled linear hyperbolic system sandwiched between 2 odes[END_REF]. From these definitions, we immediately have

i-1 j=1 1 λj , for s ∈ [t-τ - i j=1
P αi (t, s) = α i (s + i-1 j=1 1 λ j , 0), s ∈ [t -τ - i-1 j=1 1 λ j , t], (44)
P vi (t, s) = v i (s + i-1 j=1 1 λ j , 0), s ∈ [t -τ - i-1 j=1 1 λ j , t], (45)
P ui (t, s) = u i (s + i j=1 1 λ j , 1), s ∈ [t -τ - i j=1 1 λ j , t], P X (t, s) = X(s + n-1 j=1 1 λ j ), s ∈ [t -τ - n-1 j=1 1 λ j , t]. (46)
From these predictions, it is straightforward to obtain the predictions of the states α i (t, x) and β i (t, x), using the fact that α i (t, x) = α i (t -x λi , 0) +

x λ i 0 f + i (x -λ i ν)v i+1 (t - ν, 0)dν, and
β i (t, x) = ρ i,i u i (t - 1 -x µ i , 1) + ρ i,i+1 v i+1 (t - 1 -x µ i , 0) + 1-x µ i 0 f - i (x + µ i ν)u i (t -ν, 1)dν.
Finally, using the invertibility of the backstepping transformation [START_REF] Bastin | Stability of linear densityflow hyperbolic systems under PI boundary control[END_REF], we obtain the desired predictions for the original PDE and ODE states.

IV. STATE FEEDBACK CONTROL LAW We now have all the tools to apply our recursive dynamics interconnection framework.

Theorem 1: For i ∈ [1, n -1], let us define the following sequences

Vn (t) = K n (P X (t - n-1 j=1 1 λ j )), (47)
χ n-1 (t) = CP X (t - n-1 j=1 1 λ j ), (48)
Vi (t) = K i (P ui (t - i-1 j=1 1 λ j - x λ i , x), P vi (t - i-1 j=1 1 λ j - x λ i , x)) + L i ( Vi+1 (t + 1 λ i ), χ i (t + 1 λ i )), (49)
χ i-1 (t) = ρ i-1,i P vi (t - i-1 j=1 1 λ j , 0), i > 1, (50)
where the operators K i are defined in Assumption 2, the operators L i are defined in Assumption 3 and the predictors P ui , P vi and P X are defined in Assumption 4. Then, the control law V (t) = V1 (t) exponentially stabilizes the system (1)-( 5) in the sense of the Ξ-norm.

Proof: The control law V (t) is well-defined and causal, due to the definition of the different predictors (Assumption 4.) Then, applying Assumption 3 on each subsystem, we obtain that u n-1 (t, 1) exponentially converges to K n (X(t)). Consequently, X(t) exponentially converges to zero. Using Assumption 3, we can recursively show that each subsystem exponentially converges to zero starting from i = (n -1). This concludes the proof. One major advantage of the proposed framework and of the recursive design proposed in Theorem 1 is that it can easily be extended to different classes of subsystems (ODEs, non-scalar PDEs), as long as Assumptions 1, 2, 3, 4 can be verified. One must be aware that Theorem 1 completely neglects the robustness aspects of the system. More precisely, the operators K we have designed in Section 2 cancel all the boundary reflection terms and may not be strictly proper

Pα i (t, s) =                    αi(s + i-1 j=1 1 λj , 0) if s ∈ [t -τ - i-1 j=1 1 λj , t - i-1 j=1 1 λj ] qi,iPv i (t, s) + qi,i-1Pu i-1 (t, s) + δ 1 i V (s) + τ i 0 k 1 i (ν)Pα i (t, s -ν) + k 2 i (ν)Pu i (t, s -ν - 1 λi ) + dν + τ i 0 k 3 i (ν)Pv i+1 (t, s -ν - 1 λi )dν otherwise, (40)
Pv i (t, s) =                    vi(s + i-1 j=1 1 λj , 0) if s ∈ [t -τ - i-1 j=1 1 λj , t - i-1 j=1 1 λj ] ρi,i+1Pv i+1 (t, s - 1 λi - 1 µi ) + ρi,iPu i (t, s - 1 λi - 1 µi ) + τ i 0 g 1 i (ν)Pα i (t, s -ν) + g 2 i (ν)Pu i (t, s -ν - 1 λi )dν + τ i 0 g 3 i (ν)Pv i+1 (t, s -ν - 1 λi )dν otherwise, (41)
Pu i (t, s) =            ui(s + i j=1 1 λj , 0) if s ∈ [t -τ - i j=1 1 λj , t - i j=1 1 λj ] Pα i (t, s) + 1 λ i 0 f + i (1 -νλi)Pv i+1 (t, s -ν)dν otherwise, (42)
PX (t, s) =              X(s + n-1 j=1 1 λj ) if s ∈ [t -τ - n-1 j=1 1 λj , t - n-1 j=1 1 λj ] e A n-1 j=1 1 λ j (X(s) + s+ n-1 j=1 1 λ j s e A(s-ν) BPu n-1 (t, ν - n-1 j=1 1 λj )dν) otherwise. (43)
(that may lead to zero robustness margins as shown in [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF]). However, due to Assumption 1, it is possible to low-pass filter these operators and make them strictly proper, as done in [START_REF] Auriol | Stabilization of an underactuated 1+ 2 linear hyperbolic system with a proper control[END_REF].

V. APPLICATIONS A. Simulation results

The proposed control law has been tested in simulations using Matlab. The PDE systems are simulated using a classical finite volume method based on a Godunov scheme [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. We used 61 spatial discretization points (and a CFL number of 1). The predictors are implemented using a backward Euler approximation of the integral terms. The numerical values used are: n = 3, λ 1 = 1, λ 2 = 2, µ 1 = 1.3, µ 2 = 1.8, σ + 1 = 1, σ - 1 = 0.4, σ + 2 = -0.3, σ - 2 = 0.7, q 11 = 0.5, q 12 = 0.3, q 22 = 0.2, ρ 11 = 0.275, ρ 12 = 0.4, ρ 22 = 0.4, A = 0.1, B = 0.1, C = 0.1 These coefficients are chosen such that the whole system is unstable in open-loop. We have pictured in Figure 3 the evolution of the Ξ-norm of the system in open-loop and in closed-loop. As expected, the proposed control law guarantees the exponential stability. Note that the control law has been filtered by a low-pass filter to guarantee the existence of robustness margins. The important amplitude of the norm of the closed-loop system during the transient is due to the important initial value of the ODE (compared to the values of the PDEs states). Compensating this state implies propagating it through the different PDEs, which may create an artificial high norm during the transient.

B. Illustrative example

Since the control law resulting from our recursive dynamics framework requires the knowledge of all the state, it is yet unusable for practical applications. However, combined with 5) and of the closed-loop system using the control law defined in Theorem 1 a state observer, it could be applied to industrial problems. More precisely, using the mathematical model proposed by [START_REF] Aarsnes | Linear stability analysis of selfexcited vibrations in drilling using an infinite dimensional model[END_REF], [START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF], our framework could help controlling drilling devices.

VI. CONCLUDING REMARKS

In this paper, we have designed a full-state feedback control law that stabilizes a network of interconnected scalar linear hyperbolic systems coupled with an ODE at the end of the chain. We propose a recursive dynamics interconnection framework that consists of solving simpler tracking and stateprediction problems for each subsystem. Then, it becomes possible to design a stabilizing control law for the whole system recursively. The proposed control law can be lowpass filtered in order to guarantee the existence of robustness margins. We believe that this new approach can be easily adapted to stabilize more general classes of interconnected systems due to its relative simplicity (as opposed to starting a new design with an ad hoc transformation). In future works, we will consider the design of a state-observer. It is a crucial step to obtain an output-feedback law. Such an observer should encompass the prediction of the state. We will also consider a more general class of networks, for which we may have non-scalar PDE subsystems or for which we have ODE subsystems between two PDE subsystems.

Fig. 2 .

 2 Fig. 2. Schematic representation of the interactions between the different subsystems

Fig. 3 .

 3 Fig. 3. Evolution of the Ξ-norm of the open-loop system (1)-(5) and of the closed-loop system using the control law defined in Theorem 1

Jean Auriol, Silviu-Iulian Niculescu and Jeanne Redaud are with Université Paris-Saclay, CNRS, Inria, CentraleSupélec, Laboratoire des Signaux et Systèmes, 91190, Gif-sur-Yvette, France. Correspond-


J W Gao

T Miteva

Y Wu

J G Wang

A Dubois

N Sisourat
email: nicolas.sisourat@sorbonne-universite.fr

Single and double ionization processes using Gaussian type orbitals: benchmark on antiproton-helium collisions in the keV energy range

We present a simple and efficient method for computing single and double ionization cross sections in ion-atom and ion-molecule collisions using L 2 Gaussian basis sets. Gaussian functions are widely employed to compute bound states of ions, atoms and molecules. However, the description of continuum states, and therefore ionization phenomena, remains a theoretical challenge. Our approach is tested on the benchmark system antiproton-helium collisions in the so-called intermediate energy range. A good agreement with numerically exact calculations is observed. The proposed method is general and can thus be employed in any collisional systems in the challenging non-perturbative regime. Our work opens the way to investigate multiple ionization processes by ion impact in multi-center poly-electronic systems.

I. INTRODUCTION

Since the pioneer work of Boys and coworkers [START_REF] Boys | Proc. R. Soc. Lond. A[END_REF] Gaussian functions are widely used in electronic structure calculations. Although they cannot reproduce the cusp at the nucleus of the true eigenstates [START_REF] Klopper | [END_REF], they provide a fairly good description of bound states and allow the efficient computation of integrals [3]. Gaussian functions, also called Gaussian type orbitals (GTOs), are therefore used in most quantum chemistry packages (see e.g. the special issue in [4]).

Because GTOs are square-integrable L 2 functions they cannot describe the exact boundary conditions of continuum states. Although there have been large efforts to find optimal GTOs to approximate these states [5][6][7][8][START_REF] Woźniak | A systematic construction of gaussian basis sets for the description of laser field ionization and high-harmonic generation[END_REF] or to provide corrections to the GTOs finite boundary conditions (see e.g. [START_REF] Riss | [END_REF][11][12][13][14]), it is still challenging to describe molecular states above the ionization thresholds.

An accurate description of continuum states is of particular relevance to describe ionization processes. In this work, we are interested in single and double ionization phenomena in ion-atom and ion-molecule collisions in the so-called intermediate energy range. In the context of these collisions, GTOs have been widely employed to describe bound-bound transitions (i.e. charge transfer and excitation) (see for example [15][16][17][18][19][20][21] and references therein).

Ionization processes have also been considered but with less success. Using GTOs, the continuum states are approximated by the pseudo continuum states with energy above the ionization thresholds obtained by the diagonalization of the electronic Hamiltonian within the basis set. Single ionization cross sections have thus been rather accurately computed for small systems (mostly one-active electron atomic systems). Computing these cross sections for larger systems and considering multiple ionization processes remain, however, challenging and relevant tasks. In this work, we report on a simple and efficient approach for computing single and double ionization cross sections in ion-atom and ion-molecule collisions in the non-perturbative regime using GTOs. Using unsupervised machine learning classification techniques combined with pertinent properties of the pseudo continuum states, we sort these states into single and double electron continuum states. While we test this approach with the benchmark antiproton-helium collisional system, our proposed method can straightforwardly be applied to multicenter and polyelectronic systems.

The present paper is organized as follows. In the next section we outline the methods and computational details employed to compute the total and state-resolved single ionization and total double ionization cross sections. The results of the clustering analysis used to classify the pseudo continuum states into single and double electron continuum states are then reported (section III A). In section III B and III C, we present and discuss the single and double ionization cross sections in antiproton-helium collisions, respectively. They are followed by the conclusions in section IV. Atomic units are used throughout, unless explicitly indicated otherwise.

II. METHODS AND COMPUTATIONAL DETAILS

We use the well established straight-line impact parameter method (see [START_REF] Bates | Proc. R. Soc. Lond. A[END_REF][START_REF] Bransden | Charge exchange and the theory of ion-atom collisions[END_REF] and references therein) in which the nuclei follow classical trajectories while a quantum treatment describes the electronic dynamics. The position vector of the projectile relative to the target is R(t) = b + vt where b and v are the impact parameter and projectile velocity, respectively.

For a two-electron system, the time-dependent electronic Schrödinger equation reads

H e -i ∂ ∂t Ψ( r 1 , r 2 , R(t)) = 0, (1)
where H e is the electronic Hamiltonian,

H e = i=1,2 - 1 2 2 i + V T ( r i ) + V P ( r p i ) + 1 | r 1 -r 2 | , (2)
and r i , r i p = r i -R(t) are the position vectors of the electrons with respect to the target and the projectile respectively.

In our work, the electronic wavefunction is written as

Ψ( r 1 , r 2 , t) = N k=1 a k (t)χ k ( r 1 , r 2 )e -iE k t (3)
where χ k ( r 1 , r 2 ) are singlet helium pseudo states and E k the corresponding energies. These states are obtained by diagonalization of the electronic Hamiltonian matrix which is built using Slater determinants of Hartree-Fock (HF) orbitals. All configurations (i.e. HF ground state, single and double excitations) are included. These full Configuration Interaction (CI) calculations are performed using the CIPPRES [START_REF] Sisourat | CIPPRES[END_REF] plugin of the QP2 quantum chemistry package [START_REF] Garniron | [END_REF]26].Note that in Eq. 3, a single generic index k is used. However, it may refer to bound states as well as single and double electron pseudo continuum states (see below).

The insertion of Eq. 3 into Eq. 1 results in a system of first-order coupled differential equations, written in matrix form as,

i d dt a(t) = M( b, v, t)a(t) (4)
In the above equation, a(t) is the column vector of the time-dependent expansion coefficients and M is the coupling matrix. These equations are solved for a set of initial conditions (initial state i, b and v) using the Adams-Bashforth-Moulton predictor-corrector integrator. The probabilities of the transition i → f for a given impact parameter and velocity read

P i→f (b) = lim t→∞ |a f (t)| 2 . ( 5
)
The cross section for a given electronic process is then calculated as

σ f = 2π +∞ 0 bP i→f (b)db. ( 6
)
The cross sections for single and double ionization processes are obtained by plugging in the above equations (Eqs. 5 and 6) the expansion coefficients a f (t) (Eq. 3) corresponding to single and double electron pseudo continuum states. The procedure to separate these two kinds of states is discussed below.

In our computations, we used even-tempered Gaussian basis sets [27]: the Gaussian exponents α p are expressed as a geometric sequence according to α p = αβ 1-p (p = 1, 2, ...).

Three different basis sets (namely 12s9p6d3f1g, 13s10p7d4f and 14s11p7d4f) have been used

in order to evaluate the convergence of the cross sections. We report below the mean and standard deviation values of the set of cross sections obtained with these basis sets (the cross sections obtained with each individual basis set are given in the supplementary materials).

The notation 12s9p6d for example indicates that 12 Gaussian functions s ( =0), 9 functions p ( =1) and 3 functions d ( =2) were included in the basis set. In our calculations, the lowest and largest values of α were fixed, respectively, at 0.001 and 400 for s functions, and 0.001 and 40 for higher angular momentum functions. The parameter β was then chosen to reproduce the geometric sequence for a given number of functions. Depending on the basis set, between 9800 and 12 000 configurations are thus used to describe the pseudo states.

After full diagonalization of the Hamiltonian matrix, pseudo states with energy up to 5 a.u.

are kept in the close-coupling expansion (i.e. between 3000 and 5000 states).

The pseudo states of helium can be grouped into three classes: bound states which describe elastic and excitation processes, single and double electron pseudo continuum states used to model single and double ionization phenomena, respectively. The first class of states have energy below the lowest ionization potential of the system (i.e. -2.0 a.u. for helium) and can thus be clearly identified. States with higher energies are pseudo continuum states because Gaussian functions do not have the correct outgoing-wave boundary conditions.

Pseudo states with energy between the lowest ionization threshold and the double ionization one can be considered as single electron continuum states. Note that for these pseudo states, the cationic states of the remaining bound electron cannot be clearly defined because it should involve the true continuum of the cationic channels. Furthermore, above the double ionization threshold a distinction between single and double electron continuum states cannot be rigorously made. We propose below a ad hoc procedure to i) obtain the states of the remaining cation after single ionization, allowing thus to compute state-resolved excitation-ionization cross sections and ii) to separate single and double electron continuum states.

To sort the pseudo continuum states into single and double electron continuum states, one needs to find some properties, so-called features below, that are different for the two kind of states. The choice of suitable features is a difficult task of any machine learning algorithm and has thus fed the literature for the last two decades [28]. In this work, we use two features for each pseudo continuum state: the sum of the norm of the Dyson orbitals and the electron-electron repulsion energy. Hereafter, we use the notation SND for the former and EER for the latter quantity. The SND for the state k is defined as

SND k = j d r |φ kj ( r)| 2 . (7)
In this equation, φ kj ( r) are the so-called Dyson orbitals,

φ kj ( r) = d r 2 χ He + j ( r 2 )χ He k ( r 1 , r 2 ) ( 8
)
and the sum j is performed over all bound states of He + . The EER for the state k is computed as

EER k = d r 1 d r 2 χ He k ( r 1 , r 2 ) 1 | r 1 -r 2 | χ He k ( r 1 , r 2 ) ( 9
)
While there is no rigourous justification for the use of these features, our results suggest that they are relevant to sort the pseudo continuum states. Furthermore, their relevance in sorting the states may be rationalized as follows: the SND of a true single electron continuum state should be one since one electron remains bound after the ionization process while for a true double electron continuum state the SND should go to zero. The EER of both kinds of true continuum states should be close to zero and this feature should not be relevant.

However, as shown and discussed below, in a finite basis set the EER of the double electron pseudo continuum states are generally larger that those of single electron pseudo continuum states. The use of the EER is therefore pertinent.

Using these features, we have performed a clustering analysis (see III A) to group the pseudo continuum states into two "clusters" corresponding to single and double electron continuum states. K-means [START_REF] Arthur | Proceedings of the Twenty-Second Annual Symposium on Computational Geometry, SCG '06[END_REF] and spectral clustering [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] approaches, which are two wellestablished methods for unsupervised machine learning classification, have been employed.

Same clusters are obtained with both methods and only results from the K-means clustering approach are reported below. These analysis have been done with the Scikit-learn package [START_REF] Pedregosa | [END_REF].

The partial single ionization cross sections σj (i.e. after the collision He + is in the ground or excited state j) are obtained as:

σj = k P kj j P kj σ k ( 10
)
where σ k is the cross section of the pseudo continuum state k (classified as a single continuum state in the clustering analysis) and P kj is the norm of the Dyson orbital φ kj . The total double ionization cross sections are the sum of the cross sections of the pseudo continuum states belonging to the double electron continuum cluster.

III. RESULTS AND DISCUSSIONS

A. Single versus double electron continuum states

Before reporting the single and double ionization cross sections in antiproton helium collisions, we present the results of our approach to separate single and double electron continuum states. The energy of the helium pseudo states obtained with the 12s9p6d3f1g basis set are shown in Fig. 1. There are 43 bound states, 1351 states whose energies are between -2.0 and 0.0 a.u. (i.e. assigned to single electron continuum states) and 989 states with energies above the double ionization threshold. In order to illustrate the relevance of the SND and the EER as features for the clustering analysis, we plot them for every state in Fig. 2: the bound states exhibit a SND close to 1 since the helium cation states (within the basis set) represent a complete basis for the helium states. Except for the ground state of helium, the EER of the bound states are small. The SND for the higher states with energy below the double ionization threshold are also close to one. However, as the energy of the states rises the SND decreases owing to the finite basis set. The states with energy above 0 a.u. are the most interesting ones for our analysis. In order to separate these states into two groups (i.e. single and double electron continuum states) we applied the K-means clustering analysis as implemented in the Scikit-learn package. The results are shown in The antiproton-helium collisional system is a prototype and has thus been widely investigated [32][33][34][35][36][37][38][39][40] (see also [41,[START_REF] Schulz | Ion-Atom Collisions: The Few-Body Problem in Dynamic Systems[END_REF] for recent reviews on antiproton collisions and [37] for a detailed discussion on electron-electron correlation effects in the antiproton-helium system).

For this two-electron system, a full numerical solution can be performed as in [38]. In the following, we use the results of [38] as references to compare with the total cross sections computed within our L 2 Gaussian basis set approach. Total single and double ionization cross sections have also been measured, and computed with different approaches. The results of [38] agree well with the available measurements and a detailed comparison can be found in [38]. To our knowledge, state-resolved single ionization cross sections have only been reported in [START_REF] Nagy | [END_REF] and references therein. Note that these cross sections were computed with perturbative approaches. We use the results of [START_REF] Nagy | [END_REF]44] below for comparison.

-

B. Single ionization cross sections

The mean total and state-resolved single ionization cross sections computed with our three basis sets are shown in Fig. 4. The total cross sections thus obtained agree well with the reference data of [38]. The corresponding standard deviations are fairly small, showing a good convergence of the total cross sections with respect to the basis set. Our state-resolved single ionization cross sections are in general larger than those reported in [START_REF] Nagy | [END_REF], especially for the cross sections of ionization-excitation into He + (2p) (i.e. about a factor of 3 at the lowest overlapping energy and up to an order of magnitude at higher energies). We mention that in [START_REF] Nagy | [END_REF], the authors have employed the same perturbative method to compute the state-resolved single ionization cross sections for proton and electron impact too. The cross sections they obtained are overall smaller than the experimental ones (see [START_REF] Nagy | [END_REF] for a detailed discussion on the possible reasons for such disagreement). The perturbative approach may therefore underestimate the cross sections for antiproton collisions in the energy domain

C. Double ionization cross sections

The total double ionization cross sections computed with our approach are shown with the reference data of [38] in Fig. 5. For comparison, we report the cross sections obtained after the clustering analysis (labeled as "From analysis" in the figure) and the sum of the cross sections over all pseudo continuum states of positive energies ("Upper limit" in the figure). The latter represent an upper limit of the double ionization cross sections. The difference between the "Upper limit" and "From analysis" cross sections thus represent the single ionization cross sections of states having positive energies. Without the clustering analysis, the double ionization cross sections are far too large because part of the single tot Tot [38] 2s [44] 2p [START_REF] Nagy | [END_REF] FIG. 4. Total and state-resolved single ionization cross sections. The reference calculations of [38] and of [START_REF] Nagy | [END_REF]44] are shown in full lines. ionization cross sections are included. Furthermore, these upper limits depend strongly on the basis sets as shown by the large standard deviations. The cross sections obtained after the clustering analysis agree well with the reference results and depend less than the upper limits on the basis sets. It is therefore clear that the clustering analysis is crucial to obtain accurate double ionization cross sections.

IV. CONCLUSIONS

In this work, we have theoretically investigated single and double ionization processes in antiproton-helium collisions in the non-perturbative regime. The impact parameter method and a close-coupling approach using Gaussian-type orbitals are employed. Gaussian functions are widely employed to compute bound states of ions, atoms and molecules. However,

From analysis

Upper limit [38] FIG. 5. Double ionization cross sections. The cross sections obtained after the clustering analysis are labeled "From analysis" while the sum of the cross sections over all pseudo continuum states of positive energies are marked "Upper limit". The reference calculations of [38] are shown in full black line.

the description of continuum states with such L 2 functions, and therefore ionization phenomena, remain a theoretical challenge. A simple and efficient method for computing single and double ionization cross sections in ion-atom and ion-molecule collisions using these basis sets is reported and tested. Good agreement with available reference data is observed.

Our proposed method has been tested on collisions with antiproton but it can straightforwardly be applied to other collision systems. The clustering analysis method reported here is used to sort the pseudo continuum states into single and double electron continuum states and can be employed to any systems. In many relevant ion-atom and ion-molecule collisions, electron transfer can take place, which is not the case with the antiproton. In our close-coupling expansion, the electron transfer channels can be accounted for thanks to the pseudo continuum states. As in any truncated expansion, convergence of the cross sections with respect to the set of states must however be checked. Providing that convergence is reached the overall approach presented here is general. Our results thus open the way to investigate multiple ionization processes by ion impact in multi-center poly-electronic systems.

Fig. 3 :

 3 Fig.3: single electron continuum states present in general a SND close to 0.5 and small EER, while the SND for the double electron continuum states are near zero and the EER are larger. The latter observation might be counter-intuitive since two electrons in a real continuum should have nearly zero electron-electron repulsion. However, owing to the finite character of our basis set, the two electrons in the pseudo continuum are actually spatially closer than when one electron is bound (i.e. close to the nuclei) and one is "free". Using the results of this classification and the ad hoc procedure presented above, we compute the state-resolved single ionization and the double ionization cross sections. Results are shown and discussed hereafter.

FIG. 1 .

 1 FIG.1. The Helium pseudo energies obtained with the 12s9p6d3f1g basis set.

FIG. 2 .

 2 FIG.2. Sum of the norm of the Dyson orbitals (SND) plotted as a function of the electron-electron repulsion energy (EER) for every Helium states: red triangles are used for states with energy below the lowest single ionization potential, blue diamonds represent the states with energy between -2.0 and 0.0 a.u. and green dots show the features for states with positive energy. These quantities are reported in atomic units.

FIG. 3 .

 3 FIG.3. Sum of the norm of the Dyson orbitals (SND) plotted as a function of the electron-electron repulsion energy (EER) for Helium states with positive energy: those identified as single and double electron continuum states using the K-means clustering analysis are shown in red diamonds and black dots, respectively. These quantities are reported in atomic units.


Claire Silvani
email: claire.silvani@insa-lyon.fr

Julien Réthoré

| Sylvie Nicaise

Measuring coarse grain deformation by digital image correlation

Keywords: coarse granular materials, DIC, oedometric compression, mesoscopic scale, strain

This work presents results from oedometric compression of coarse granular material. Coarse granular media exhibit significant deformations making it complicated to predict the settlement of structures. In this paper, a measurement technique was developed for the analysis of two-dimensional (2D) images of a deforming coarse granular medium to investigate its deformation. This was achieved by realizing grainbased image correlation to measure the grain transformation in gravel with the use of a Digital Image Correlation technique. The 2D displacement fields enable us to explore the behavior of granular media at different scales: microscopic, mesoscopic, and macroscopic scales. The mesoscopic scale is defined from branches that connect the centers of three neighboring grains, using a Delaunay triangulation to account for an equivalent continuum media. While the consistency of the macroscopic strain and the average mesoscopic strain is assessed, it is shown that a deviation from the normalized microscopic vertical displacement is an indicator of the heterogeneity of the mesoscopic strain field. The proposed mesoscopic analysis allows us to investigate these heterogeneities. Another important result is that even if the amplitude of the microscopic strain is small (approximately 100 times smaller) compared with the other strain measures, it confirms that the grains are not rigid and that their ultimate strain can be estimated using the proposed approach.

INTRODUCTION

This study focuses on the behavior of coarse granular media under compression. Rockfill media are coarse media composed of quarry rock debris or crushed rock fragments with a characteristic size of up to 1 m. This kind of medium is used in various civil engineering works, such as roads, retaining walls, railway engineering works, and in a large number of embankment dams for the coarsest material.

The mechanical behavior of granular materials is complex essentially due to their discrete nature. To improve our knowledge of granular media it is necessary to understand the evolution of the internal state of this heterogeneous material. Here, experimental work was conducted to offer insights into the strain behavior of granular media. Many authors have explored the strain behavior of granular material via theoretical, numerical, and experimental approaches (Satake, 1992 [28], Kruyt and Rothenburg, 1996 [12], 2016 [START_REF] Kruyt | A micromechanical study of dilatancy of granular materials[END_REF], 2018 [START_REF] Kruyt | A strain-displacement-fabric relationship for granular materials[END_REF], Kruyt 2003 [14], Bagi, 1996 [4], Kuhn, 1999 [16], Cambou et al., 2000 [6], Kuhn and Daouadji 1999 [START_REF] Kuhn | Multi-directional behavior of granular materials and its relation to incremental elastoplasticity[END_REF] for example for theoretical and numerical approaches). for example, for theoretical and numerical approaches). Furthermore, extensive experimental work has been performed to study the kinematics mostly of shear bands in granular materials (such as sands) in triaxial loading paths with X-ray micro-tomography, and with experiments carried out in situ [START_REF] Desrues | Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography[END_REF] [START_REF] Desrues | Void ratio evolution inside shear bands in triaxial sand specimens studied by computed tomography[END_REF], Viggiani et al. 2004, [31], [START_REF] Hall | Strain localisation in sand under triaxial loading: characterisation by x-ray micro tomography and 3d digital image correlation[END_REF] [START_REF] Hall | Strain localisation in sand under triaxial loading: characterisation by x-ray micro tomography and 3d digital image correlation[END_REF]) or during plane strain compression (Rechenmacher 2006 [24] and Rechenmacher et al. 2011 [START_REF] Rechenmacher | Characterization of mesoscale instabilities in localized granular shear using digital image correlation[END_REF]) with the use of Digital Image Correlation (DIC) to obtain the grain kinematics. These investigations provided important data concerning shear band features such as orientation and width, patterns of shear band formation, and their evolution (Viggiani et al., 2001 [32]). Similarly, [START_REF] Sibille | A numerical photogrammetry technique for measuring microscale kinematics and fabric in schneebeli materials[END_REF] [START_REF] Sibille | A numerical photogrammetry technique for measuring microscale kinematics and fabric in schneebeli materials[END_REF] ] measured the microscale kinematics and fabric properties of cylindrical wooden rods such as two-dimensional (2D) analog granular materials. In the past few years (Andò et al. 2012 [START_REF] Andò | Experimental micromechanics: grain-scale observation of sand deformation[END_REF] and [START_REF] Andò | Grain-scale experimental investigation of localised deformation in sand: a discrete particle tracking approach[END_REF]), the mechanisms involved in strain localization (in triaxial tests) were experimentally investigated at the microscale and the quantification of individual grain kinematics (displacements and rotations) of large quantities of sand grains was realized (in a test sample undergoing loading). However, there is still a challenge to extract quantitative descriptors of grain-scale processes, e.g., kinematics and spatial interactions (Vlahinic et al. 2014 [START_REF] Vlahinić | Towards a more accurate characterization of granular media: extracting quantitative descriptors from tomographic images[END_REF]) to explain macroscopic strain. Concerning the theoretical works, the micro-mechanical formulation of the stress tensor has been well established for a long time, whereas the corresponding formulation for the strain tensor has proved to be much more evasive and is the subject of much discussion (Duràn et al., 2010 [START_REF] Durán | Analysis of three-dimensional micro-mechanical strain formulations for granular materials: evaluation of accuracy[END_REF]). To improve the models, some authors have introduced an intermediate scale, called "the mesoscale," at which strain can be defined using the tools of continuum mechanics (Bagi, 2006 [5]). Establishing a link between particle-level displacements and macro-level deformations is important for researchers who intend to develop micromechanically based constitutive theories (Nicot and Darve 2011 [START_REF] Nicot | The h-microdirectional model: accounting for a mesoscopic scale[END_REF], Chaze and Cambou 2014 [START_REF] Chaze | Analysis of internal state and strains in granular material at meso-scale: influence of particle shape[END_REF], Xiong et al. 2017 [START_REF] Xiong | A three-dimensional micromechanically based model[END_REF]). Such an intermediate local scale is expected to allow the local structure of granular materials to be more precisely described. Hence, micro-mechanical models embedding both microscale and mesoscale can usually describe a variety of constitutive features in a natural way.

The aim of the present study is to contribute to the theoretical and numerical works on the introduction of an intermediate scale for representing macro strain. To do so, an experimental study is carried out, based on the development of a DIC tool to extract the evolution of key properties directly from the grain-scale mechanics of granular media under oedometric compression. The objective is to measure the displacement of each individual grain and then to compute mesoscopic strain from grain displacements. As theoretical and numerical approaches, we introduce this intermediate scale to link grain motion to the macroscopic deformation of the media. To our knowledge, there are few experimental data avalaible on this topic (Houda et al., 2016 [11]), ]), with most papers focusing mainly on discrete element simulations (Nguyen et al., 2009 [19], 2012 [START_REF] Nguyen | Local stress analysis in granular materials at a mesoscale[END_REF] and 2016 [START_REF] Nguyen | Towards a new approach for modeling the behavior of granular materials: A mesoscopic-macroscopic change of scale[END_REF], Zhu et al., 2016 [38] and 2017 [START_REF] Zhu | On the mechanics of meso-scale structures in two-dimensional granular materials[END_REF]), since they represent solutions that are more affordable to implement.

A series of experiments on gravel are conducted in this work. The first sections of this paper describe the experimental set-up and materials, the DIC technique, and the developments. The following sections examine the results of the experiments through different scales (macro-, micro-, and mesoscales) and conclude with a comparison between these different scales and the possible relationships between them.

EQUIPMENT, MATERIALS, AND METHODS

Experimental set-up and materials

The study of coarse granular materials requires suitable laboratory facilities whose dimension must be in proportion to the maximum size of the grains. To this end, oedometric tests are performed. The oedometric cell used here was developed by the IRSTEA laboratory in Aix-en-Provence (France) and has a square surface of 30 cm x30 cm and is 40 cm high. Some of the cell walls are in Plexiglas (Figure 1), to guarantee the optic accessibility of surfaces. Only the upper part of one face of the box is transparent; the bottom is opaque. DIC is then used to measure the grain displacements in 2D on the projected surface of the grains. Stress is applied via a loading plate attached to a frame that is connected to a hydraulic cylinder to impose the load. The frame used to impose the load does not allow for the use of a single camera to follow the whole upper part (see Figure 2). Therefore, two cameras are used to follow the whole upper part, and each camera points toward half of the upper part (Figure 2). Lighting of the surfaces of interest is made with LED or with more diffuse ligthing. Several tests on limestone were carried out: The gravel is from the Ain region near Lyon (France) and is composed of grains whose diameters range between 6 mm and 12 mm. Grains are dropped off randomly into the cell. Granulometry was performed before and after the tests. The principal characteristics of the material are summed up in Table 1. Several of the tests performed have the same initial granulometry charateristics (Figure 4). Cameras with ultra-high resolution (29 megapixels) are used to analyze precisely the grain field. Pictures were taken every 60 s. Tests with imposing load velocity ̇ are performed (2 kN/min or 26 kPa/min). The displacement of the loading plate is monitored via a linear variable differential transformer (LVDT) sensor. Parts of the Plexiglas walls were sometimes damaged by the material, and scratches appeared that made DIC difficult, if not impossible. Thin sheets of removable Plexiglas were then placed behind the thicker Plexiglas cell walls so that several experiments could be performed.

Digital image correlation

In the past decade, image correlation techniques have gained attention in experimental geomechanics. Today, the use of this technique, which was introduced approximately 30 years ago by Peters and Sutton (Peters et al., 1982 [23], Sutton et al., 1983 [30]), is widespread. Image-based analysis methods are widely used in particulate media in several laboratory investigations in the field of geomechanics: cone penetration testing (e.g, Paniagua et al. 2013 [START_REF] Paniagua | Soil deformation around a penetrating cone in silt[END_REF], Arshad et al. 2014 [START_REF] Arshad | Experimental study of cone penetration in silica sand using digital image correlation[END_REF]), load transfer in granular platforms supported by piles in soft soil (Houda et al. 2016 [11]), and fluctuations in shear tests (Richefeu et al., 2012 [26]). DIC was chosen to follow the surface displacement (vertical and horizontal) of each grain individually. This technique assumes that the surface texture (gray-level distribution) does not change during the test. Let us call the reference image and the deformed image, the displacement field, the pixel coordinates in the reference image frame, and * = + ( ) the advected coordinates of these pixels in the deformed image. It is possible to write (Equation ( 1)):

( ) = ( + ( )) = ( * ) (1)
Usually, it is necessary to apply a speckle on the object studied to make the machting between and easier. In our case, the color gradation of the granular media enables us to work without applying a speckle pattern. We used the software UFreckles, described in Réthoré 2018 [START_REF] Réthoré | Ufreckles[END_REF]. All the DIC parameters are summarized in Appendix B. We consider only 2D displacements (vertical and horizontal). We disregarded the out-of-plane component of the displacement, since grains were placed behind the window, and the media is very constrained in such an oedometric compression test. To solve the ill-posed problem (Equation 1), the reference image is cut into small zones Ω, and attached to each individual grain. Since hundreds of grains are considered, it is not possible to define these zones with a manual grain contour capture. However, the center of each grain is picked up manually (see Figure 5). From this set of points (assumed to be the center of the particles), a Delaunay triangulation is derived (see Figure 6). Voronoi cells are then computed from this triangulation and each analyzed grain surface Ω is defined as a scaling of its Voronoi cell (Figure 6). The displacement of each individual grain is then searched for as an affine transformation (constant strain and rigid body motion), meaning that for each grain a full micro strain tensor (later referred to as ) is obtained as well as a rigid translation vector ( , ) and an in-plane rotation .

FIGURE 5

Methodology for the image processing: identification of the grain center.

From this assumtion, the discretized displacement field over the grain surface Ω reads as (2):

( ) = ∑ ∈N ( ) (2)
where N is the set of degrees of freedom, and the basis functions of the considered displacement description. The basis functions can be any basis of a linear function of the pixel coordinates. The corresponding degrees of freedom are related to , ( , ) and which are thus obtained independently for each grain. To find the transformation of each zone Ω, the goal is to obtain the best approximation of the field displacement . This consists in the resolution of an inverse problem. The best solution is the one that minimizes the function as defined in Equation ( 4).

( ) = ( 2 ) (3) 2 = ∫ ∫ Ω ( ( ) -( * )) 2 (4)

MULTI-SCALE ANALYSIS

The purpose of this section is to explain how each grain of the granular media behaves individually and how it interacts with its neighbors. More than 250 grains are followed behind the transparent glass. The ratio of cell size/grain size is equal to 26 to limit the influence of the lateral boundary conditions. No breakage appears for the range of macroscopic loads considered here. Some grains near the edges are not taken into account to avoid convergence problems. Different scales and parameters have been studied herein to compare their order of magnitude:

• is the strain derived from the loading plate displacement and the height of the box,

• will be computed from the displacements of the grain center ( , ) based on the Delaunay triangulation,

• is the microscopic strain of the grains obtained directly from the DIC analysis at the grain scale.

An other quantity ∕ 0 is also computed for each grain, with 0 the initial position of the grain center. This normalized microscopic displacement is intrisically related to the grain motion and it will thus be classified as a microscopic parameter. This parameter has been defined because if the set of grains is deformed as a homogeneous continuum submitted to a constant strain then its value should also be constant and its averaged value should match the one of the macroscopic strain (see Appendix A). Thus, by computing the averaged normalized displacement ⟨ ∕ 0 ⟩ , any deviation of the grain set kinematics compared with a homogeneous continuum can be detected.

The different quantities enumerated will be explored and compared thereafter and are presented in the following scheme (Figure 7). Note that compression will be considered positive and expansion will be considered negative.

Macroscopic scale

Figure 8 shows some of the experiments conducted. Tests were performed several times to check the macroscopic reproducibility. The curves (Figure 8) are obtained by the LVDT sensor. The two curves are quite close, which confirms the reproducibility of the test In the rest of the paper, the macroscopic strain is evaluated from a DIC measurement of the vertical displacement of the loading plate. This allows us to have exactly the same reference state (defined by the reference image for DIC) for the different strain measures presented.

Microscopic scale

The horizontal and vertical ⟨ ⟩ are compared (Figure 9). The horizontal mean micro-strain is low compared with the vertical one, which increases more with the loading. Note that the grains are compressed in the two directions and their volume change is mainly governed by the sum of these two strains. A maximum relative volume change of approximately 0.0125% is obtained.

Normalized vertical displacements ∕ 0 are shown in Figure 10 as a function of the macroscopic vertical strain. Grains move downward in a continuous way (the vertical displacement is nearly linear with respect to

). No remarkable rearrangement is visible. Note that the slope of the plots in Figure 10 should be the same for all the grains and should be 1 if the set of grains behaves as a continuum submitted to a homogeneous strain transformation. The fact that this slope is not the same for all the grains is proof that the deformation of the set of grains differs from that of a continuum submitted to the same loading.

The mean ⟨ ∕ 0 ⟩ is now compared with the mean micro-vertical strain ⟨ ⟩ . ⟨ ⟩ is very low compared with the mean normalized vertical displacement of the whole grains ⟨ ∕ 0 ⟩ (Figure 11), as expected. Note that, again, this figure is an illustration of the non-homogeneous deformation of the granular medium we observed. Indeed, in the case of a continuum, the average ⟨ ∕ 0 ⟩ should match the macoscopic strain and the slopes in Figure 11 should be 1 (see Appendix A). ⟨ ∕ 0 ⟩ will be compared with mesoscopic strains hereafter.

Mesoscopic scale

Recent work (e.g., Bagi, 2006 [5], Chaze and Cambou, 2014 [START_REF] Chaze | Analysis of internal state and strains in granular material at meso-scale: influence of particle shape[END_REF], Zhu et al., 2016 [38], Nguyen et al., 2016 [20]) showed that an intermediate scale is necessary to link the grain displacements to macroscopic strain, i.e., the meso-scale, which is defined hereafter.

To compute strain, the assembly is replaced by a continuous domain, and displacement of the characteristic points of the equivalent continuum media is equal to the displacement of the particles centers. To define local strain, a meso-scale is defined, whose branches connect the centers of three particles in contact, using Delaunay triangulation.

From the Delaunay triangulation built using the grain centers as nodes, the gradient of the displacement field over each triangle is determined on the basis of particle displacements and geometrical characteristics, as in the finite element method with shape function on each node. The displacement field is linear within a triangle, and thus a constant gradient can be assigned separately to each cell, as done by Bagi, 1996 [4], Cambou et al., 2000 [6] (for example). Plane strain is assumed. The meso-strain for each cell c is then computed (Equation ( 5)):

= ̄ { } (5)
with the strain vector for a cell c, ̄ the shape functions tensor derivative for each triangular element, and { } the nodal displacement vector for a cell c. A total 445 Delaunay cells are defined in the proposed example. Note that the Delaunay triangulation is defined on the initial configuration and is unique during the whole test. Then, an average strain is computed as follows (Equation 6) and defines the mean meso-strain of the sample at each step.

= 1 ∑ =1 ∑ =1 (6)
Figure 12 shows that we have a random mesh with the Delaunay triangulation. A criterion is given to eliminate the too elongated cells (in blank in the figure) that are mainly located on the boundary of the sample. The strain level of the mesoscopic cells seems to be randomly distributed and some strain localization appears for a few cells. The horizontal and vertical mean mesoscopic strains are presented in Figure 13. The mean vertical mesoscopic strain evolution is a regular function of the macroscopic strain (Figure 13). The average horizontal strain is not exactly zero. Two reasons can explain this result. First, all the grains in the experiment have not been taken into account. As mentionned earlier, some grains near the interfaces are not taken into account. Secondly, at the beginning of the test the grains do not fill the entire available space in the test box which allows for an expansion in the horizontal direction while the vertical compression is applied.

Comparison between scales

Figure 14 compares the four vertical strains already defined. The mean microscopic strain is very low compared with the macro-strain. However, with the proposed strategy it is far above the measurement uncertainty. It can thus be concluded that the micro-strain is indeed negligible compared to the macro-strain but that the grains are not rigid. The level of strain measured inside the grain is indeed related to stress which may lead to grain failure for loads of higher amplitude. is very close to , which is consistent with its definition. This implies that the major events are captured by the DIC measurements behind the window and that is a good indicator of even if the rear of the sample is not explored. One may imagine that some major rearrangements that occur (inside the loading cell or close to the observation window) have an influence limited to their FIGURE [START_REF] Kruyt | Micromechanical definition of the strain tensor for granular materials[END_REF] Map of mesoscopic strain (magnitude) at the end of the test.

FIGURE 13

Mean meso vertical strain as a function of vertical macro-strain.

local vicinity only and have no major impact on the displacement measurement of the plate and thus on . However, the average mesoscopic strain over the planes in the vicinity of the rearrangements may deviate significantly from . For example, if the rearrangements induce a rotation of the loading plate from its initial orientation, the average of the mesoscopic strain can deviate from the macroscopic strain. The mean normalized vertical displacement of the grains ⟨ ∕ 0 ⟩ also has the same order of magnitude as but it deviates more. The slope of the evolution of ⟨ ∕ 0 ⟩ as a function of is observed to be lower than 1. This indicates the presence of deviations of the local mesoscopic strain (confirmed in Figure 12) from the macroscopic strain, which can be seen as the signature of the discrete nature of the media. Plotting the evolution of ⟨ ∕ 0 ⟩ for each grain allows us then to appreciate the heterogeneity of the meso-strain field experienced by the granular material. It is also concluded from Figure 14 that this deviation from the deformation of a continuum increases continuously as the prescribed deformation increases. This confirms that no remarkable rearrangement is visible but that the relative motion between the grains and their own deformation (micro-strain) producing heterogeneities in the mesoscopic strain field increases as well.

CONCLUSION

In this paper, we present an experimental exploratory work to understand the mechanisms of coarse granular media (gravel) deformation under oedometric compression. A DIC technique in 2D is used to explore the granular media. The challenge was then to extract quantitative information from 2D images to understand the mechanics at play, and then to elucidate and validate the pertinent scale to accomodate for the macroscopic strain. A method is developed to process the digital images: An area is defined from the center of the gravel particles proportional to the Voronoi cell deduced from the set of grain centers. The tool developed enables us to measure displacement fields in an unfriendly environment: a cell with several Plexiglas windows with dust and scratches.

The application of this tool reveals different scales of interest. Several scales are investigated (microscopic, macroscopic, mesoscopic) to understand better the granular media deformation. The DIC technique offers the possibility of monitoring the phenomena at a microscopic scale, allowing to extract information of the evolution of key properties directly from grain-scale kinematics, and thus enabling a better comprehension of the behavior of the material. This microscopic strain represents only about one thousandth of the macro-strain but it confirms that the grains are not perfectly rigid and that the stress levels corresponding to the measured micro-strain can be significant. Analyzing the micro-strain could be of interest for estimating the load failure of grains in-situ (not in an experiment involving unrealistic loadings compared with what grains are submitted to in-situ). A mesoscopic scale is introduced based on a Delaunay triangulation of the set of grain centers. The gradient of the displacement field of each triangular cell is determined and defines the mesoscopic strain field. It appears that the average meso-strain equals the macroscopic strain by construction but does not allow to appreciate the heterogeneity of the strain field experienced by the analyzed material. This can be appreciated by analyzing the mesoscopic strain field. Another microscopic parameter is computed for each grain (the normalized microscopic displacement ∕ 0 ). It has been shown that the mean normalized vertical displacement deviates from the macroscopic strain, with this deviation being a consequence of the heterogeneity of the strain distribution at the mesoscopic scale. The experimental quantification of the scatter between this strain measure and the macro-strain allows to also appreciate the heterogeneous nature of the transformation of the granular material. The work presented here is a first step toward validating theoretical works and discrete element simulations concerning this issue. The key contribution of the paper is the development of a grain-kinematics measurement technique for the determination of links between deformation at the continuum macro-scale and at the grain scale, focusing on the specific parameters to be calculated so as to account for the heterogeneous nature of the transformation of the granular material. Further, with this method, any deviation from the idealized case of a uniformly applied strain can a priori be detected from a deviation between the macroscopic strain and the average mesoscopic strain.

Like many other works, this one has certain limitations. This exploratory study is carried out in 2D whereas the media is essentially 3D, which introduces artifacts disturbing the field measurements with out-of-plane motions. However, the out-ofplane motions are limited since the Plexiglass wall prevents motion and thus particles are in an intermediate case (between 2D and 3D). Performing tomography should help to overcome this limitation. Moreover, due to the nearly "1D compression" global deformation state, the rotation of the grains was not analyzed. To better understand granular media behavior, with more motions and rearrangements, triaxial tests whould also be more adequate; therefore triaxial tests under tomography may be the best and most elaborate way to proceed and should represent a major step forward in this issue. Advanced works on triaxial tests using X-ray tomography have been done in recent years by several researchers. For example, Vlahinić et al. 2017 [START_REF] Vlahinić | From computed tomography to mechanics of granular materials via level set bridge[END_REF], extracted particle morphology from 3D X-ray tomography images, but also reported with high accuracy the spatial interactions between particles (contact statistics) as well as grain-level kinematics (particle rotations and translations) as a function of load state. Wiebicke and coworkers (Wiebicke et al. 2017 [35] and [START_REF] Wiebicke | Measuring the evolution of contact fabric in shear bands with x-ray tomography[END_REF] 2020) linked the metrology of inter-particle contacts to the investigation of the contact fabric evolution. However, to date, this cutting-edge work does not include multi-scale analysis to study the deformation of granular media.

Finally, this displacement field characterization can be used to provide a concrete, experimental basis for discrete element simulations in order to obtain force distribution in the media and thereby understand force transmission and breakage criteria. This work is in progress.

where is the final vertical displacement of the plate and the initial vertical position of the plate (height of the sample) with the following definition for and (Equation A2):

= 1 ∫ 0 (A2)
The objective of this part is to compute final ⟨ ∕ ⟩ for the whole 1D system as follows (Equation A3):

⟨ ⟩ = 1 ∫ 0 (A3)
and compare it with the final . Three functions (one linear and two non-linear functions that lead to = ∕ ) are choosen for the study. Table A1 summarizes the results. When the function is quadratic, the quantity ⟨ ⟩ deviates from the macro-strain. From this simple analysis, it appears that the average meso-strain equals the macroscopic strain by construction. Thus the

Quantity Function

Function Function average meso-strain does not allow to appreciate the heterogeneity of the strain field experienced by the analyzed material. The heterogeneity of the deformation experienced by the material can only be observed from the meso-strain maps, and there is no global index derived from usual strain estimates quantifying the heterogeneity of the deformation. Alternatively, for any case when the mesoscopic strain is not homogeneous, the strain measure ⟨ ⟩ differs from the macroscopic strain. Hence, this is a way to quantify the heterogeneity of the mesoscopic strain field from a global point of view. Once one detects that ⟨ ⟩ differs from the macroscopic strain , a deeper inspection of the mesoscopic strain field is required.

B DIC HARDWARE AND ANALYSIS PARAMETERS

All the DIC parameters are summurized in Tables B2 andB3.

FIGURE 1

 1 FIGURE 1 Experimental devices: two cameras follow a quarter of the sample each (only the upper part of the cell is transparent).

FIGURE 2

 2 FIGURE 2 Principle of field measures with the oedometric cell using DIC (top view of the experiment). Each camera follows a quarter of the sample, with only the upper part of the cell being transparent.

FIGURE 3 TABLE 1

 31 FIGURE 3 Gravel with 6 mm < D < 12 mm.

FIGURE 4

 4 FIGURE 4 Typical initial granulometry for the gravel tests.

FIGURE 6

 6 FIGURE 6 Methodology for image processing: definition of grain area (dark gray) as a scaling of the Voronoi cells (black dashed lines) and Delaunay triangulation (red) computed from the grain center.

FIGURE 7

 7 FIGURE 7 Schematic view of the different quantities computed.

FIGURE 8

 8 FIGURE 8Macroscopic behavior of the gravel (LVDT sensor measures for strain).

FIGURE 9

 9 FIGURE 9Mean micro-strain as a function of vertical macro-strain.

FIGURE 10

 10 FIGURE 10 Normalized grain vertical displacement as a function of vertical macro-strain. A different color is assigned to each grain. The dashed line indicates a slope equal to 1.

FIGURE 11

 11 FIGURE 11Mean micro-strain ( and ∕ 0 ) as a function of vertical macro-strain.

FIGURE 14

 14 FIGURE 14Comparison of various strains as a function of vertical macro-strain.

TABLE A1

 A1

	Correspondences between	⟨ ⟩ ,	and	.

ACKNOWLEDGMENTS

The authors thank IRSTEA Aix en Provence for the technical assistance for performing oedometric tests.

Conflict of interest

The authors declare no potential conflict of interests.

APPENDIX A ESTIMATION OF THE NORMALIZED VERTICAL DISPLACEMENT IN 1D COMPRESSION

To have an order of magnitude of the normalized vertical grain displacement, an assumption is made to deal with a 1D compression test in the vertical y direction (low horizontal displacements). Let us call U this vertical displacement. Several grain displacement functions can be assumed. Let us assume different functions for grain displacements (linear and non linear since is not constant for each element). The final of the sample can be written with the following notations (Equation A1): Small strain Post-Filtering of Strains NA Displacement Noise-Floor* 0.05 px Strain Noise-Floor* Emicro 7.5e-6, Emeso 0.003 (*) Before the test started, four images were acquired. The displacement and strain computed from the analysis of these images are used to estimate the noise-floor of the proposed methodology. The overall noise-floor for the displacement, Emicro and Emeso are obtained by computing the standard deviation of the sought field including all its components at all points and for all four images.

TABLE B3 DIC Analysis Parameters