Jérôme Weiss
email: jerome.weiss@ujf-grenoble.fr

Lucas Girard

Florent Gimbert

David Amitrano

Damien Vandembroucq

Finite) statistical size effects on compressive strength

Keywords: Mechanical strength, size-effect, compressive failure, critical transition, depinning

come

Significance statement

Large structures generally fail under stresses significantly lower than small ones. This is the size effect on strength, one of the oldest problems of engineering, already discussed by Leonardo da Vinci and Edmé Mariotte centuries ago. One classical explanation is the weakest link hypothesis: the largest the "chain", the larger the probability to find a weak link whose breaking will set the failure of the whole chain. We show, however, that it is irrelevant in case of compressive loading, a situation particularly crucial for e.g. geotechnical problems.

Interpreting compressive failure as a critical transition between an "intact" and a "failed" state, we quantitatively explain the size effects on compressive strength of materials such as concrete, rocks, coal, ice, or granular materials. \body

Introduction

Owing to its importance for structural design [START_REF] Bazant | Fracture and size effect in concrete and other quasibrittle materials[END_REF], the elaboration of safety regulations [START_REF]Eurocode 2: Design of Concrete Structures -Part 1-1: General rules and rules for buildings[END_REF], or the extrapolation of laboratory results to geophysical field scales [START_REF] Heuze | Scale effects in the determination of rock mass strength and deformability[END_REF], the size effects on strength of materials is one of the oldest problems in engineering, already discussed by Leonardo da Vinci and Edmé Mariotte(4) several centuries ago, but still an active field of research [START_REF] Bazant | Scaling theory for quasibrittle structural failure[END_REF][START_REF] Zapperi | Current challenges for statistical physics in fracture and plasticity[END_REF]. As early as 1686, Mariotte(4) qualitatively introduced the weakest-link concept to account for size effects on mechanical strength, a phenomenon evidenced by Leonardo da Vinci almost two centuries earlier. This idea, which states that the larger the system considered, the larger the probability to find a particularly faulty place that will be at the origin of global failure, was formalized much later by Weibull [START_REF] Weibull | A statistical theory of the strength of materials[END_REF]. Considering a chain of elementary independent links, the failure of the chain is obtained as soon as one link happens to break. By virtue of the independence between the potential breaking events, the survival probability of a chain of N links is obtained by the simple multiplication of the N elementary probabilities. Depending on the properties of the latter, the global survival probability converges toward one of the three limit distributions identified by Weibull, Gumbel [START_REF] Gumbel | Statistics of extremes[END_REF] and Fréchet respectively. Together with Fisher and Tippett [START_REF] Fisher | Limiting forms of the frequency distribution of the larges or smallest members of a sample[END_REF], these authors pioneered the field of extreme value statistics. This purely statistical argument, undoubtedly valid in 1D, was extended by Weibull [START_REF] Weibull | A statistical theory of the strength of materials[END_REF][START_REF] Weibull | A statistical distribution function of wide applicability[END_REF] to account for the risk of failure of 3D samples or structures. Beside the hypothesis of independence, it thus requires an additional hypothesis of brittleness: the nucleation of any elementary crack at the microscopic scale from a pre-existing flaw is assumed to immediately induce the failure at the macro-scale. More specifically, following linear elastic fracture mechanics (LEFM) stating that crack initiation from a flaw of size s occurs at a stress ~ / , one gets a probability of failure of a system of size L under an applied stress σ, ) , ( L P F  , that depends on the distribution of pre-existing defect sizes. Assuming a power law tail for this distribution, Weibull statistics are expected [START_REF] Weibull | A statistical theory of the strength of materials[END_REF],

                           m u d F L L L P    0 exp 1 ) , (
, whereas Gumbel statistics are expected for any distribution of defect sizes whose the tail falls faster than a power law [START_REF] Gumbel | Statistics of extremes[END_REF][START_REF] Alava | Size effects in statistical fracture[END_REF][START_REF] Sornette | Critical phenomena in natural sciences[END_REF],

                           u d F L L L P    exp exp 1 )
, ( 0, where m is the so-called Weibull's modulus, d is the topological dimension and L 0 and  u are normalizing constants. For Weibull statistics, the mean strength and the associated standard deviation (σ f ) then scale with sample size L as ~ ~ / . This approach has been successfully applied to the statistics of brittle failure strength under tension [START_REF] Weibull | A statistical theory of the strength of materials[END_REF][START_REF] Beremin | A local criterion for cleavage fracture of a nuclear pressure vessel steel[END_REF], with m in the range 6 to 30 [START_REF] Miannay | Fracture Mechanics[END_REF]. It implies a vanishing strength for L→+∞, although this decrease can be rather shallow, owing to the large values of m often reported.

Although relying on strong hypotheses, this weakest-link statistical approach was almost systematically invoked until the 70's to account for size effects on strength whatever the material and/or the loading conditions. However, as shown by Bažant [START_REF] Bazant | Fracture and size effect in concrete and other quasibrittle materials[END_REF][START_REF] Bazant | Scaling theory for quasibrittle structural failure[END_REF], in many situations the hypothesis of brittleness is not obeyed. This is in particular the case when the size of the fracture process zone (FPZ) becomes non-negligible with respect to the system size. In this so-called quasi-brittle case, an energetic, non-statistical size effect applies [START_REF] Bazant | SIZE EFFECT IN BLUNT FRACTURE -CONCRETE, ROCK, METAL[END_REF], which has been shown to account for a large variety of situations [START_REF] Bazant | Scaling theory for quasibrittle structural failure[END_REF]. In what follows, we do not consider (deterministic) energetic size effects and explore a situation, compressive failure, where both the hypotheses of brittleness (in the sense given above) and independence are not fulfilled, up to very large scales. Relaxing these initial hypotheses of the weakest-link theory, our statistical physics approach remains statistical by nature and relies on the interplay between internal disorder and stress redistributions. It is based on a mapping of brittle compressive failure onto the critical depinning transition of an elastic manifold, a class of models widely used in non-equilibrium statistical physics characterized by a dynamic phase transition [START_REF] Fisher | Collective transport in random media: from superconductors to earthquakes[END_REF]. This approach does not consider sample's shape effects [START_REF] Hustrulid | A review of coal pillar strength formulas[END_REF], only statistical size effects. The critical scaling laws associated to this phase transition naturally predict a saturation of the compressive strength at large scale and are in remarkable agreement with measurements reported for various materials such as rocks, ice, coal, or concrete.

Compressive failure cannot be captured by a weakest-link approach

Compressive loadings are particularly relevant in rocks mechanics and geophysical situations [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF] as the result e.g. of lithostatic pressure, and consequently for geotechnical problems (e.g. ( 18)). Brittle compressive failure is a complex process, as the local tensile stresses at crack tips are counteracted by the far-field compressive stresses. Consequently, Griffith-like energy balance arguments, or related tools such as fracture toughness, cannot be developed to describe the instability leading to terminal failure, thus making the weakest-link approach inoperative. Instead, brittle compressive failure involves an initiation phase, elastic interactions and stress redistributions, as well as frictional sliding along rough surfaces. In what follows, we mean by brittle failure a situation where microscopic ductile deformation processes such as creep or dislocation motion play a negligible role [START_REF] Schulson | Brittle failure of ice[END_REF]. During the initiation phase, secondary cracks nucleate from the local tensile stresses generated by the frictional sliding along pre-existing defects such as grain boundaries, small joints, or microcracks [START_REF] Nemat-Nasser | Compression-induced nonplanar crack extension with application to splitting, exfoliation, and rockburst[END_REF][START_REF] Schulson | On the initiation of shear faults during brittle compressive failure: a new mechanism[END_REF]. The propagation of these mode I secondary cracks is however rapidly stopped by the farfield compression. Instead, such nucleation events locally soften the material [START_REF] Kachanov | Elastic solids with many cracks and related problems[END_REF][START_REF] Cox | MICROCRACK FORMATION AND MATERIAL SOFTENING IN ROCK MEASURED BY MONITORING ACOUSTIC EMISSIONS[END_REF] and thus cause a redistribution of elastic stresses, which in turn can trigger other microcrackings. Then, in the course towards failure, the interaction and linking of secondary cracks is considered to be at the onset of shear fault formation, from which the macroscopic instability is thought to result [START_REF] Schulson | On the initiation of shear faults during brittle compressive failure: a new mechanism[END_REF][START_REF] Katz | Microfracturing, damage, and failure of brittle granites[END_REF]. This process is characterized by a progressive localization of damage and deformation along a fault [START_REF] Lockner | Quasi-static fault growth and shear fracture energy in granite[END_REF].

The above description shows that all the assumptions of the weakest-link theory are inappropriate in case of compressive failure. Summarizing experimental field and laboratory data obtained for 50 years, it is thus not surprising that their weakest-link predictions are poorly obeyed. When the compressive strength of brittle materials has been measured from laboratory tests over a limited scale range (generally between ~ 10 -2 m and ~ 10 -1 m), either non-significant [START_REF] Kuehn | The compressive strength of ice cubes of different sizes[END_REF][START_REF] Van Mier | Multiaxial strain-softening of concrete, Part I: Fracture[END_REF] or limited [START_REF] Mogi | The influence of the dimensions of specimens on the fracture strength of rocks[END_REF][START_REF] Del Viso | Shape and size effects on the compressive strength of high-strength concrete[END_REF] size effects on were observed, whereas, when reported, the associated variance seemed to increase towards small scales [START_REF] Kuehn | The compressive strength of ice cubes of different sizes[END_REF]. Consequently, these results do not fully constrain empirical or theoretical size effect formulations. Some studies were performed instead several decades ago over a much larger scale range (~ 10 -2 m to few m), combining laboratory and in-situ tests [START_REF] Hustrulid | A review of coal pillar strength formulas[END_REF][START_REF] Pratt | The effect of specimen size on the mechanical properties of unjointed diorite[END_REF][START_REF] Bieniawski | The effect of specimen size on compressive strength of coal[END_REF][START_REF] Brace | The effect of size on mechanical properties of rock[END_REF]. All of them reported a significant scale-dependence of at small scales, tentatively and empirically fitted as a power law decrease(18) ~ , but also a non-zero asymptotic strength at large (> 1m) scales, not explained by the weakest link approach. So far, there is no clear explanation for this non-vanishing compressive strength. Instead, empirical formulations of size effects on compressive strength of brittle materials [START_REF] Heuze | Scale effects in the determination of rock mass strength and deformability[END_REF][START_REF] Hustrulid | A review of coal pillar strength formulas[END_REF][START_REF] Hoek | Practical estimates of rock mass strength[END_REF] generally ignore such asymptotic behaviour. Following observations at small scales, they all share a common power law scaling ~ , with β varying from very small values(29) (i.e. almost no size effect), to the LEFM scaling β =1/2. The weakest-link concept has been sometimes put forth to explain this scaling for small β values( 29), although it is clear from above that this approach is irrelevant in case of compressive failure. On the other hand, a (deterministic) energy analysis of compression failure based on physical (micromechanical) considerations has been proposed [START_REF] Bazant | Fracture and size effect in concrete and other quasibrittle materials[END_REF][START_REF] Bazant | Compression failure of quasibrittle materials and size effect[END_REF]. In agreement with the scenario described above, it considers that the nucleation of microcracks roughly parallel to the principal compression axis form a band whose mechanical instability, triggered by the buckling of the microslabs separating the microcracks, leads to failure. However, the microcracks, and therefore the associated band, are assumed to nucleate suddenly, just preceding macro-failure. I.e., this approach does not consider the progressive route towards the failure, characterized by elastic interactions between cracks and progressive damage localization. In other words, the transition to failure is considered as a "first-order" transition. This, in addition to an assumed constant scaling between the band length and the size of the system, gives a vanishing strength towards large scales with an asymptotic scaling ~ / , i.e. slightly shallower than the LEFM scaling. Consequently, the observed nonvanishing strength is not explained. In addition, this deterministic approach cannot, by nature, account for a size dependence of the variability of strength. We propose instead to consider compressive failure as a critical transition and develop a mapping onto the depinning transition that allows accounting for the interplay between local disorder and long-ranged elastic interactions leading to a statistical finite size effect.

Compressive failure as a critical depinning transition

The modeling of the mechanical behavior of heterogeneous materials induced in recent years an intense research activity. From the early nineties emerged the idea that non-linear processes such as fracture, plasticity and damage could be discussed as critical phenomena [START_REF] Herrmann | Statistical models for the fracture of disordered media[END_REF]. In the context of damage, a paradigmatic example of this approach is given by the fiber bundle model [START_REF] Pradhan | Failure processes in elastic fiber bundles[END_REF]. However, the scope of this model as well as its variants [START_REF] Delaplace | Damage cascade' in a softening interface[END_REF] are restricted to the catastrophic failure occurring under tensile conditions, i.e. the transition from an initiation stage to a propagating stage triggered by the development of a critical nucleus. In contrast our interest here is the study of progressive damage under compressive conditions.

To our knowledge the first attempt of a description of compressive damage as a critical phenomenon is due to Toussaint and Pride [START_REF] Toussaint | Fracture of disordered solids in compression as a critical phenomenon. I. Statistical mechanics formalism[END_REF]. They developed a statistical mechanics formalism based on ensemble averages obtained over the rock seen as collection of disordered mesovolumes. A specific Hamiltonian [START_REF] Toussaint | Fracture of disordered solids in compression as a critical phenomenon. II. Model Hamiltonian for a population of interacting cracks[END_REF] accounted for the interaction between cracks and the traditional tools of equilibrium statistical mechanics (partition function, maximum of entropy) were used to characterize the localization transition associated with the failure of the material.

We here follow a different route. We proposed recently a numerical progressive damage model whose results are consistent with an interpretation of brittle compressive failure as critical phase transition [START_REF] Girard | Fracture as a critical phenomenon in a progressive damage model[END_REF][START_REF] Girard | Damage-cluster distributions and size effect on strength in compressive failure[END_REF]. This finite-element model( 41) considered a continuous elastic material with progressive local damage: the elastic modulus of an element decreases each time the stress state on that element exceeds a given threshold defined by a Coulomb criterion. This elastic softening simulates an increase in microcrack density at the element scale [START_REF] Kachanov | Elastic solids with many cracks and related problems[END_REF][START_REF] Cox | MICROCRACK FORMATION AND MATERIAL SOFTENING IN ROCK MEASURED BY MONITORING ACOUSTIC EMISSIONS[END_REF]. Disorder was introduced on the local stress threshold. As the result of elastic interactions, the stress redistribution following a damage event can trigger an avalanche of damage. We showed(41, 42) (i) that the size of the largest damage cluster as well as of the largest damage avalanche diverge at peak load, which just precedes failure, and (ii) the divergence of a correlation length  at failure, ~Δ / , where Δ (respectively Δ

) is the control parameter for strain-(respectively stress-) driven simulations, the applied macroscopic strain, the corresponding value at peak stress  f (failure), and =1.00.1 the correlation length exponent.

Here, in the spirit of a recent model of amorphous plasticity [START_REF] Talamali | Strain localization and anisotropic correlations in a mesoscopic model of amorphous plasticity[END_REF], we formalize this interpretation of compressive failure as a critical transition through a mapping onto a depinning model, a class of models that exhibit out-of-equilibrium phase transition. The damaging process is described as the motion of a d-dimensional elastic manifold with longrange interactions through a random field of obstacles within a space of dimension d+1. In our case, the macroscopic stress plays the role of the driving force and a local damage event corresponds to the depinning from an obstacle, with elastic redistributions in both cases.

Damage is represented by a scalar field at a mesoscopic scale λ, that typically corresponds to the correlation length of the structural disorder of the material, and it occurs whenever the stress state σ reaches the boundary of an elastic domain defined by the Coulomb criterion, | | . This criterion is of wide applicability for brittle materials under compressive stress states to define the onset of damage [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF][START_REF] Weiss | Coulombic faulting from the grain scale to the geophysical scale: Lessons from ice[END_REF]. and are respectively the shear and normal stress (sign convention positive in tension) over a plane maximizing the Coulomb stress | | , and μ tan is an internal friction coefficient ( is the socalled angle of internal friction). The heterogeneous nature of the material, i.e. the disorder, is accounted for by a statistical variability of the cohesion , .

A crucial feature is the non-independence of the local damaging events occurring in the material. Any local event is characterized by a local decrease of the elastic modulus that occurs in a small region surrounded by the remainder of the material. The latter can be seen as an elastic matrix and its reaction induces an elastic stress field all over the material. We started from the classical inhomogeneity problem of Eshelby [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | Elastic field of an ellipsoidal inclusion[END_REF] to calculate the stress field induced by a damaged inclusion. Since the damaged material is elastically disordered, the internal stress cannot be obtained by a simple superposition of the contributions of isolated inclusions. In order to partly account for interactions between inclusions, we develop a two-step strategy to compute the internal stress (see SI). The damage field is first used to obtain a self-consistent estimate of the average, macroscopic elastic behavior, . This effective value, which partly accounts for interactions between inclusions, is then used to obtain a fluctuating damage field, defined from the contrast between the actual elastic moduli at micro-scale and . The internal stress is finally obtained from the sum of the elastic contributions of the associated effective inclusions. The interplay between local disorder and elasticity is the basis for depinning models, which have proved successful in recent years to describe physical and mechanical phenomena as various as the advance or receding of triple contact line on a disordered substrate [START_REF] Joanny | A MODEL FOR CONTACT-ANGLE HYSTERESIS[END_REF][START_REF] Moulinet | Width distribution of contact lines on a disordered substrate[END_REF][START_REF] Bonn | Wetting and spreading[END_REF], the motion of a magnetic wall in a thin film [START_REF] Lemerle | Domain wall creep in an Ising ultrathin magnetic film[END_REF], the propagation of a crack front in a heterogeneous material [START_REF] Gao | A First-Order Perturbation Analysis of Crack Trapping by Arrays of Obstacles[END_REF][START_REF] Santucci | Statistics of fracture surfaces[END_REF][START_REF] Dalmas | Crack Propagation through Phase-Separated Glasses: Effect of the Characteristic Size of Disorder[END_REF][START_REF] Bonamy | Failure of heterogeneous materials: A dynamic phase transition?[END_REF], etc. The full derivation of our problem is given in the SI. It allows a complete mapping onto a depinning model, with the following equation of evolution of the damage field:

, , / ,

where denotes the positive part, µ is a mobility coefficient and λ is the characteristic length scale of the microstructural disorder. In the language of depinning models, we identify as the external forcing term, as the elastic contribution induced by the damage field D, via the knowledge of the damage dependent effective modulus , and the disorder.

This formal mapping enables us to apply to progressive damage the various analytical and numerical results obtained in the framework of the depinning transition. In the "thermodynamic" limit, i.e. for a system of infinite size, a well-defined critical threshold  th separates the static phase (limited damage) from the dynamic one (failure). For our problem, this is the expression of a non-vanishing asymptotic strength,  ∞ >0. In case of samples of finite size, fluctuations of the measured threshold  f are expected in the vicinity of  th . We thus expect, as for the threshold force of the depinning transition [START_REF] Zapperi | Current challenges for statistical physics in fracture and plasticity[END_REF][START_REF] Narayan | Threshold critical dynamics of driven interfaces in random media[END_REF], a finite-size scaling for the compressive strength of the form: ). This implies that in case of weak disorder and will be of the order of (e.g. grain size, aggregate size,..), but might be significantly larger in case of strong disorder when e.g. cracks or joints widely distributed in size are initially present in the material. The classical assumption( 55) is = FS , while the mean-field prediction( 56) is =1. Equation ( 2)

expresses the variability on strength intrinsically related to the failure process, on which experimental sources of variability should be added. Towards very small scales, ≪ , the proposed scaling (equations ( 2) and ( 3)) necessarily breaks down when  f approaches the material strength limit [START_REF] Bazant | Fracture and size effect in concrete and other quasibrittle materials[END_REF].

Application to experimental data in cohesive materials

In full qualitative agreement with experimental data (see above), this finite-size scaling implies an apparent power law decay of the mean strength at small sizes, a non-vanishing strength for L→+∞, and an increasing variability towards small sizes. Relation (2) is hardly testable from experimental data, as (σ f ) values, when reported, are based on a limited number of independent tests and include experimental-related scatter. For studies including field tests at the meter scale [START_REF] Hustrulid | A review of coal pillar strength formulas[END_REF][START_REF] Pratt | The effect of specimen size on the mechanical properties of unjointed diorite[END_REF][START_REF] Bieniawski | The effect of specimen size on compressive strength of coal[END_REF][START_REF] Brace | The effect of size on mechanical properties of rock[END_REF] and assuming that the asymptotic strength  was reached at the largest scale, we fitted the data with relation (3). The agreement is remarkable, with the bestfit  FS value ranging from 0.8 to 1.05 (Fig. 1), i.e. close to the mean-field prediction,  =1.

The corresponding length scales L B range from ~ 20 to 40 cm, a possible sign of relatively strong disorder (joints, microcracks) in these natural rock samples. For studies based only on lab tests, either:

(i) no significant size effect on strength is reported, as for fresh-water granular ice [START_REF] Kuehn | The compressive strength of ice cubes of different sizes[END_REF], limestone [START_REF] Thuro | Scale effects in rock strength properties. Part I: unconfined compressive test and Brazilian test[END_REF], granite [START_REF] Thuro | Scale effects in rock strength properties. Part I: unconfined compressive test and Brazilian test[END_REF], or concrete [START_REF] Van Mier | Multiaxial strain-softening of concrete, Part I: Fracture[END_REF]. This can be explained by a small L B in equation ( 3) and/or an insufficient dataset to properly sample size effects;

(ii) or the data can be well fitted by (3) assuming  FS =1, as shown on figure S1 for highperformance (HP) concrete [START_REF] Del Viso | Shape and size effects on the compressive strength of high-strength concrete[END_REF] and marble [START_REF] Mogi | The influence of the dimensions of specimens on the fracture strength of rocks[END_REF]. In case of HP concrete, the scale L B is close to the maximum size of the andesite aggregate (12 mm) [START_REF] Del Viso | Shape and size effects on the compressive strength of high-strength concrete[END_REF]. In agreement with our former expectation, in such initially unfractured materials, the microstructural scale (aggregate size, grain size, ..) likely sets this L B scale.

The confining pressure  3 increases the axial compressive strength  1f of rocks, ice, coal or concrete [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF][START_REF] Van Mier | Multiaxial strain-softening of concrete, Part I: Fracture[END_REF][START_REF] Hoek | Practical estimates of rock mass strength[END_REF][START_REF] Renshaw | Universal behavior in compressive failure of brittle materials[END_REF]. Up to a confining ratio  3 / 1f of about 30%, failure is brittle and occurs through microcrack initiation and interactions, followed by shear fault formation at the onset of macroscopic instability, as described above [START_REF] Renshaw | Universal behavior in compressive failure of brittle materials[END_REF]. This failure mode is sometimes called Coulombic faulting, reminiscent of the importance of solid friction in this case [START_REF] Weiss | Coulombic faulting from the grain scale to the geophysical scale: Lessons from ice[END_REF].

Consequently, one expects our mapping to the depinning transition to hold in this case. The combination of the effects of size and of confining pressure on strength has been rarely

studied, but the available data on coal [START_REF] Hoek | Practical estimates of rock mass strength[END_REF] are well explained by equation ( 2) with  FS =1 and an increasing asymptotic strength  ∞ with increasing confinement, as expected (Figure 2). For these natural samples, the scale L B is once again relatively large (several cm). It slightly decreases with increasing confinement, suggesting a secondary effect of confinement on the sensibility of L B to the variability of the local strength, . For such multiaxial compression tests, the deviatoric stress appears as the most relevant variable. Thus, the strength has been defined as . This choice for ,instead of the axial strength σ 1f , obviously doesn't change the value of the exponent  FS or of the scale L B . For confinements larger than ~ 30%, compressive failure is no more brittle, and another failure mode occurs, as mode I secondary crack nucleation is inhibited. This failure mode, called plastic faulting, involves thermal softening and an adiabatic shear instability [START_REF] Renshaw | Plastic faulting: Brittle-like failure under high confinement[END_REF]. In this case, we no longer expect elastic interactions between microcracks to occur, i.e. our size effect formalism to hold. Indeed, it has been found that for large confining pressure, size effects on compressive strength disappear [START_REF] Habib | Sur la disparition de l'effet d'échelle aux hautes pressions[END_REF]. This sets the range of applicability of our formalism.

Application to granular media

This mapping onto the depinning problem is likely not restricted to brittle cohesive materials.

As described in [START_REF] Talamali | Strain localization and anisotropic correlations in a mesoscopic model of amorphous plasticity[END_REF] and recalled in the SI, it can be extended to the macroscopic plastic instability in amorphous media. The case of a cohesionless frictional granular medium compressed under confinement can be interpreted as an intermediate case between amorphous plasticity and compressive damage. Indeed, shear-induced local rearrangements of the granular structure lead to irreversible local strains but not to a systematic degradation of local stiffness. Compared to amorphous plasticity, other complications are present, however, such as dilatancy. When compressed under confinement, these media exhibit a macroscopic flowing instability associated to strain localization [START_REF] Desrues | Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry[END_REF], which sets the yield stress, i.e. the "strength". This instability can also be considered as a critical transition [START_REF] Gimbert | Crossover from quasi-static to dense flow regime in compressed frictional granular media[END_REF]. In this case, the disorder is topological, coming from the arrangement of particles.

From this analogy, we expect finite-size scaling (relations (2) and ( 3)) to ensue. However, to our knowledge, there is so far no experimental data over a significant range of scales to check this anticipation. We therefore simulated the mechanical behaviour of frictional granular materials using the Molecular Dynamics discrete element method [START_REF] Radjai | Discrete-element modelling of granular materials[END_REF]. Two-dimensional granular assemblies made of a set of frictional circular grains were considered. The dynamic equations were solved for each grain, which interact via linear elastic laws and Coulomb friction when they are in contact [START_REF] Agnolin | Internal states of model isotropic granular packings. I. Assembling process, geometry, and contact networks[END_REF]. Neither cohesion between grains, nor rolling resistance were considered. In order to build granular assemblies with strongly different initial (before loading) characteristics, in terms of coordination number and/or packing density, specific sample preparation procedures were used. Details on the discrete element model as well as on these procedures are given in the SI. These granular assemblies were loaded under a multi-axial configuration, with the external axial stress 1 prescribed in order to impose a constant axial strain-rate, whereas the radial stress 3 , i.e. the confining pressure, was kept constant.

The 2D sample sizes varied from 100 to ~ 45000 grains.

Whatever the initial characteristics of the assemblies, finite size scaling of compressive strength was observed, in full agreement with equations ( 2) and (3) (Fig. 3),

showing the generic nature of the concept proposed here. In agreement with our expectation, the scales L A and L B were slightly larger than the average particle size, and increased for less dense, less coordinated samples.

Failure strength statistics

As noted in the introduction, the weakest-link hypothesis leads to extreme value statistics for the probability of failure under an applied stress σ. As the weakest-link theory appears irrelevant for compressive failure, we do not expect such extreme statistics for the distribution of strength in this case. Published experimental data with a sufficient number of failure tests to analyze strength distributions are rare. Results obtained on ice indeed exclude extreme statistics, either Weibull or Gumbel, and argue instead for Gaussian statistics (Fig. 4). The same is true for the discrete-element modeling of frictional granular media (Fig. S2). We

Concluding comments

This new, statistical physics interpretation of compressive failure of continuous and granular media has important practical consequences. First, when lab-scales (cm to dm) studies show no significant size effect, one expects that lab strength values give a good estimate of the asymptotic (field) strength. Extrapolation of lab-scales data to scales smaller than L A or L B will be more difficult, owing to the intrinsic variability at such scales. However, the meanfield estimate of the finite-size exponent,  FS =1, obtained from theoretical considerations, well describes the fluctuations and the associated finite-size corrections, whereas for initially unfractured materials, L A and L B are related to the characteristic microstructural scale (grain size, aggregate size, ..). Therefore, owing to its predictive potential, we believe that the proposed scaling is a useful, simple to use guidance for future structural design rules or regulations (e.g. ( 2)).

Materials and Methods

The characteristics and the simulation settings of the discrete-element model of frictional granular media are given in the Supporting Information, along with the formal derivation of the mapping of brittle compressive failure onto the depinning transition of an elastic manifold.

  FS is the finite-size exponent and  ∞ a non-vanishing asymptotic value of the strength for L→+∞. A, B (in Pa.m / ), L A = and L B (in m) are constants. These length scales define the scales below which respectively the fluctuations and the finite-size corrections become important compared to the asymptotic strength  ∞ . We expect these to scale as , ~ , where represents the associated variability on the local cohesive strength (see S.I.

  anticipate, from the criticality of the transition, the scaling form of the distribution , of the fluctuations for a system of size L as , Ψ . Such a scaling form naturally leads to the scaling relations for the mean value (relation (3)) and the standard deviation (relation (2)) of the compressive strength discussed above. However, the precise form of the statistical distribution Ψ is not prescribed by this simple scaling analysis. In particular, Ψ is not expected to obey the predictions of extreme value statistics whose hypotheses (absence of interactions) are not satisfied in the present problem. In recent results obtained in a similar framework (depinning model of amorphous plasticity (65)), Gaussian-like distributions were observed as well.Combining Gaussian statistics with equations (2) and (3) leads to the following expression for the probability of failure at scale L under a stress :

Figures

  Figures

Acknowledgments

All numerical simulations were performed at SCCI-CIMENT Grenoble. J.W. and D.V. acknowledge the hospitality of the Aspen Center for Physics, which is supported by the National Science Foundation Grant No. PHY-1066293, as the seminal ideas of this work came up during their stay at the Center. S. Zapperi and D. Bonamy are acknowledged for interesting discussions and suggestions.

for the mean strength and equation ( 2) for the standard deviation, with  FS =1.07, L A =1.68, L B =4.21 and

The best-fit exponent  FS and scale L A were obtained from the standard deviation scaling (bottom), the asymptotic strength  ∞ and scale L B were then obtained from the scaling of (inset of (top)). Inset of (top): Same data and fits, in a vs The collapse onto a single straight line, which corresponds to equation (4), argues for Gaussian statistics.


N Lestoille

C Soize

G Perrin

C Funfschilling

Long time evolution of train dynamics with respect to track geometry

Keywords: high speed train dynamics, stochastic modeling, track irregularities, stochastic dynamics, long time evolution, vehicle-track coupling

This paper aims at characterizing the long time evolution of the vehicletrack system. The knowledge of the evolution of such a system is of great concern for the railway industry, in order to maintain a high level of safety and comfort in the high speed trains. We propose a computational stochastic approach to predict the long time evolution of a given track portion. The approach is based on an adaptation of the global stochastic model of track irregularities previously identified with a large experimental data basis. The nonlinear stochastic dynamics of the train excited by track irregularities are carried out using a computational multibody dynamics model. Some indicators concerning the dynamic responses of the train are introduced in order to start off the maintenance or not of the given track portion.

Introduction

The tracks for high speed trains are submitted to more and more solicitations, because of the increase of the train traffic, the load and the speed of the trains. These solicitations induce degradations of the track geometry, making evolve track irregularities. Such degradations impact the train dynamic response in return. To guarantee a good level of safety and comfort of the train, maintenance operations of the track have to be regularly undertaken. These maintenance operations are heavy and costly. They would gain being started off by indicators on the train safety and comfort, and no more on track-irregularities measures. However, the vehicle-track system is highly nonlinear and there is no obvious correlation between track irregularities and the vehicle-track forces (see [START_REF] Li | Assessing track geometry quality based on wavelength spectra and track-vehicle dynamic interaction[END_REF][START_REF] Lonnbark | Characterization of track irregularities with respect to vehicle response[END_REF]). The knowledge of the effects of the track geometry on the train dynamic response needs to be increased in order to define robust indicators for the simulated train dynamic response and to characterize their evolution.

The goal of this study is to characterize the influence of the long time evolution of the track irregularities on the evolution of the train dynamic response. The distinction has to be done between the long time evolution, which will be denoted by τ , and the time of the train dynamics (denoted by t). Actually, the vehicle-track system is a complex system, with high nonlinearities and coupling between inputs (track geometry, track stiffness, track mass) and outputs (train responses). The track geometry is the main source of excitation for the train. The track geometry is known thanks to track measurements which are made very precisely and frequently. A global stochastic model of the track geometry has been built by Perrin et al. in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] using a very large experimental data basis concerning the French railway network for high-speed trains. The stochastic modeling is very useful to carry out nonlinear stochastic dynamic analysis of the train excited by the random track geometry. Nevertheless, we need now to analyze the long time evolution of a given track portion, in order to construct criteria which allow maintenance operations to be started off. The most relevant information we have about the geometry of the given portion is only the projection of this portion on the global stochastic model. Since this projection requires a large number of coefficients (several thousands), the use of the long time evolution of all these coefficients does not allow to study the long time evolution of this portion, because the dimension (the number of coefficients) is too high. That is why we chose to analyze the long time evolution of the train dynamics on the given portion. For the long time evolution analysis of a given track portion, for which measurements are periodically carried out, the global stochastic model is adapted to this portion in introducing a noise which allows measurements errors and variability to be taken into account. Using this adapted stochastic model of the track geometry for such a given portion, the train dynamic response is numerically simulated using a commercial software (Vampire). To assess the long time evolution of the train dynamics, we have to define indicators on the train dynamic response and to observe their evolution over a long time.

Section 2 deals with the proposed approach. Section 3 will focus on the adapted stochastic modeling of the track irregularities, taking into account the measurements for the given track portion analyzed. In Section 4, the train dynamic response will be analyzed, setting up indicators on the train dynamics to quantify the long time evolution of the train stability and comfort.

Proposed approach

Objectives

The track irregularities change during long time τ , inducing modifications of the train dynamic response. Criteria of safety and comfort are needed to assess the long time evolution of the train dynamic response. To guarantee a good level of safety and comfort, maintenance operations on the track have to be regularly undertaken. The final objective of this work is to be able to decide to start off maintenance operations of a given track portion as a function of the simulated train dynamic response, and even to be able to anticipate maintenance operations. To achieve this goal, the long time evolution of the vehicle-track system has to be characterized using measurements of the track, which are periodically performed for the track portion under consideration.

The first main stage consists in studying the sensitivity of the train stochastic dynamics to the long time evolution of track irregularities.

Numerical simulation of the train dynamic response

The inputs of the simulation are the track design, the track irregularities modeled with the adapted stochastic model evoked in Section 1 and detailed in Section 3, and a model of the train. The model of the train used for the simulation is the one presented in [START_REF] Kraft | Parameter identification for a TGV model[END_REF]. It is a multibody dynamical system whose dynamic responses are computed using a commercial software (Vampire). It is assumed that the speed of the train is constant and fixed. The simulation outputs are accelerations in the train and contact forces between the wheels and the rails. The numerical simulation of the train response is performed on the given track portion for which the length, S, has to be defined.

The decision to start off the maintenance of the given portion is taken on the basis of indicators (called dynamic indicators) which are defined as functions of the train dynamic responses (outputs). These indicators must be sensitive to the long time evolution of the track irregularities. The dynamic indicators are based on statistical quantities related to the simulation outputs. Especially, random maxima of the stochastic outputs are of great concern, since we want to analyze safety and comfort levels for the given track portion. The obtained results could be used later to define thresholds on the dynamic indicators, in order to be able to start off maintenance operations from these thresholds.

It should be noted that, since the train is assumed to have a constant speed, it is equivalent to analyze the dynamic response as a function of the time t or as a function of the spatial location of the train defined by the curvilinear abscissa s. We also use the dual variable of s which is the wave number, k, and the associated wavelength λ = 2π/k. In such a case, we can use the terminology: "wavelength in the spectrum of the train dynamic response".

Choice of the length of the studied given track portion

The length S of the given track portion has to be a compromise between the need of local information to start off the maintenance and the respect of characteristic lengths related to the train dynamics. Two characteristic lengths are taken into account to define S: the longest wavelength in the spectrum of the dynamic response and the length of the track transition section which is located between the given portion and the previous section used in the numerical simulation.

Longest wavelength in the spectrum of the train dynamic response.. The train dynamic response is simulated on a section (including the given portion) which is much longer than the train length, to be sure to observe all the wavelengths contained in the spectrum of the train response. The power spectral density functions (spectrum) of the train responses are estimated (periodogram method with a Hamming window). Figure 1 is then defined such that:

S 1 2π k min . ( 1
)
Length of the track transition section.. The length of the track transition section, S trans is defined using 25 simulations of the train dynamics for 25 tracks made up of two successive portions including the same given portion.

The length S trans is such that the 25 train dynamic responses are the same on the portion S \ S trans . Length S 2 is chosen such that:

S 2 S trans . ( 2
)
Length S of the studied given track portion.. Length S is defined as the maximum of S 1 and S 2 :

S = max{S 1 , S 2 } . ( 3
)

Stochastic modeling of track irregularities

Track measurements

The track geometry is measured by a measuring train very precisely and very frequently. The measuring train is equipped with laser cameras which measure the distances between the rails. The track is described by two data sets:

• the initial track design, which corresponds to the theoretical track (as it was planned before the construction) and which is made of perfect lines and curves.

• the irregularities of the track which appear during the track life cycle and which have to be added to the track design.

The track irregularities are modeled by a vector-valued random field X denoted by

s → X(s; τ ) = (X 1 (s; τ ), X 2 (s; τ ), X 3 (s; τ ), X 4 (s; τ )) , (4)
which is defined on S and which depends on long time parameter τ . Its components are the lateral and vertical alignment irregularities X 1 and X 2 , the cant deficiencies X 3 and the gauge irregularities X 4 . The given track portion is defined for s in Ω = [0, S], in which S the length chosen in Section 2.3. Long time τ is a discrete parameter that rises between successive measurements of the given track portion,

τ 0 < τ 1 < τ 2 < ... < τ ντ , ( 5
)
in which τ 0 is the time of the first measurement just after a maintenance operation, τ 1 , τ 2 , ..., τ ντ corresponds to the successive long times for which there are measurements of the track geometry, and τ ντ is the time of the last measurement before the next maintenance operation.

Adaptation of the global stochastic model of the track geometry to the given track portion

The global stochastic model proposed in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] is very robust with respect to measurements errors. This model has the capability to generate a track irregularity for any given portion belonging to the French railway network. Nevertheless, as explained in Sections 1 and 2, we are interested in constructing a stochastic model adapted to the given track portion which has the capability to predict the long time evolution of this portion using the local measurements x meas τ 0 , x meas τ 1 , x meas τ 2 , . . . of this portion at long times τ 0 , τ 1 , τ 2 , . . .. The objective of this adapted stochastic model which has to be constructed is (i) to take into account measurements noise associated with local measurements x meas τ 0 , x meas τ 1 , x meas τ 2 , . . . and (ii) to take into account the local variability of the given track portion in order to decrease the "statistical distance" between the global stochastic model and the local measurements.

The method proposed to construct the adapted stochastic model of the track irregularities consists in introducing a random field noise for which the spatial properties are driven by the global stochastic model and whose intensity of its statistical fluctuations is identified at long time τ 0 using measurement x meas τ 0 . This construction is detailed in Section 3.2.2. In order to simplify the reading, we give a short summary of the global stochastic modeling presented in [START_REF] Perrin | Track irregularities stochastic modeling[END_REF] whose details can be found in [START_REF] Perrin | Random fields and associated statistical inverse problems for uncertainty quantifications -Applications to railway track geometries for high-speed trains dynamical responses and risk assessment[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension form a set of realizations[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF][START_REF] Perrin | Karhunen-Love expansion revisited for vector-valued random fields: scaling, errors and optimal basis[END_REF][START_REF] Perrin | A posteriori error and optimal reduced basis for stochastic processes defined by a finite set of realizations[END_REF].

Short review on global stochastic modeling for track irregularities

The track irregularities vector X = (X 1 , X 2 , X 3 , X 4 ) is modeled by a vector-valued random field, defined on a probability space (Θ, T , P), indexed by Ω = [0, S], with values in R 4 . It has been proven that random field X is neither Gaussian nor stationnary (not homogeneous).

The global stochastic model of the track irregularities has been built and identified solving an inverse statistical problem using a very large experimental data basis related to the French railway network. Random field X is centered,

E{X(s)} = 0 , ∀s ∈ [0, S] , (6)
where E{.} is the mathematical expectation. The continuous vector-valued random field X(s), s ∈ Ω, is replaced by its spatial discretization at curvilinear abscissa s n = n h with h the spatial step and n = 0, . . . , N s , where S = N s h. Keeping the same notation for the continuous random field and its spatial discretization, the following random vector is introduced, X = (X(0), X(h), X(2h), . . . , X(N s h)) [START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data[END_REF] with values in R 4(Ns+1) .

The first step of the construction consists in using a truncated principal component representation of random vector X which is written as

X Nη k=1 λ k u k η k , ( 8
)
in which N η 4(N s + 1), where {λ k , 1 ≤ k ≤ N η } are the N η first largest eigenvalues, and where {u k , 1 ≤ k ≤ N η } are the associated eigenfunctions of the covariance matrix [C XX ] of X. The random coefficients {η k , 1 ≤ k ≤ N η } are non-correlated centered and normalized random variables such that:

η k = 1 √ λ k X T u k , 1 ≤ k ≤ N η . ( 9
)
Introducing the (4(N s + 1) × N η ) real matrix [U] defined by

[U] = u 1 . . . u Nη , [U] T [U] = [I] , (10)
in which [I] is the unity matrix, the diagonal matrix [λ] whose diagonal entries are λ 1 , . . . , λ Nη , and the random vector η such that, η = (η 1 , . . . , η Nη ) , [START_REF] Le Maitre | Spectral Methods for Uncertainty Quantification[END_REF] then random vector X can be rewritten as

X [Q] η , ( 12
)
in which the matrix [Q] is written as

[Q] = [λ] 1/2 [U] . ( 13
)
The second step consists in constructing a polynomial chaos expansion [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF][START_REF] Le Maitre | Spectral Methods for Uncertainty Quantification[END_REF][START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF] of the non-Gaussian second-order random vector η whose non-correlated components are statistically dependent. The truncated polynomial chaos representation is then written as

η η chaos (N) = N j=1 y (j) ψ j (ξ) , ( 14
)
in which ξ = (ξ 1 , ..., ξ Ng ) is a uniform second-order centered random vector whose covariance matrix is the unity matrix (consequently, the components are statistically independent and have a unit variance). The non-Gaussian dependent random variables ψ 1 (ξ), . . . , ψ N (ξ) are the multivariate normalized Legendre polynomials of random vector ξ. The parameters N g , N and N η are identified in order to minimize the errors of approximation. Using Equations ( 12) and ( 14), random vector X can be represented by,

X [Q] [y] Ψ(ξ) , ( 15
)
in which the (N η × N) real matrix [y] is defined by [y] = [y (1) . . . y (N ) ] and where the R N -valued random variable Ψ(ξ) is defined by Ψ(ξ) = (ψ 1 (ξ), ..., ψ N (ξ)). Any realization X(θ), θ ∈ Θ, of random field X, generated with the above stochastic model, is representative of the measurements x meas τ .

Local stochastic modeling

The local stochastic modeling consists in constructing an adapted stochastic model of the track irregularities related to a given track portion. Since the track irregularity is composed of four different and dependant types of irregularities, this model is adapted to each component X 1 , X 2 , X 3 and X 4 of the random field {X(s) , s ∈ Ω}. The corresponding spatial discretization yields the following decomposition of the random vector X = (X 1 , X 2 , X 3 , X 4 ). For k = 1, . . . , 4, X k is a random vector of dimension N s + 1 which is rewritten (using Equation ( 12)) as

X k [Q k ] η , ( 16
)
in which the ((N s + 1)

× N η ) real matrix [Q k ] is extracted from matrix [Q].
The proposed adapted stochastic model is written as

X k (δ k ) = [Q k ] η + δ k G k , k = 1, 2, 3, 4 , ( 17
)
in which δ = (δ 1 , δ 2 , δ 3 , δ 4 ) is the vector-valued parameter allowing the noise level to be controlled and which has to be identified using experimental data (see Section 3.2.3). For fixed k, the R Nη -valued random variable G k is a spatial discretization of a random field noise. In the model proposed,

G = (G 1 , G 2 , G 3 , G 4
) is chosen as a R 4Nη -valued Gaussian second-order centered random variable whose covariance matrix is the unity matrix and which is defined on the probability space (Θ , T , P ). From Equation ( 17), the adapted stochastic model can be rewritten as

X k (δ k ) = X k + B k (δ k ) , ( 18
)
in which the random vector B k (δ k ) which depends on δ k is such that

B k (δ k ) = δ k [Q k ] G k . ( 19
)

Identification of the parameter δ opt

The optimal parameter δ opt is estimated by using the maximum loglikelihood method with experimental data. For fixed k, the following random observation is introduced,

W k (δ k ) = X k (δ k ) E{ X k } . ( 20
)
in which X k = X k (0) is the Euclidean norm of the global stochastic model X k . We then introduce the random observation vector defined by

W(δ) = (W 1 (δ 1 ), W 2 (δ 2 ), W 3 (δ 3 ), W 4 (δ 4 )) . ( 21
)
For fixed k, the corresponding experimental measurement w meas k is defined by

w meas k = x k,meas E{ X k } , 1 ≤ k ≤ 4 . ( 22
)
Finally, the experimental observation vector w meas is introduced such that

w meas (δ) = (w meas 1 (δ 1 ), w meas 2 (δ 2 ), w meas 3 (δ 3 ), w meas 4 (δ 4 )) . ( 23
)
Let L W (w meas ; δ) = log p W (w meas ; δ) be the log-likelihood in which p W (w meas ; δ) is the value of the probability density function w → p W (w; δ) of random vector W for w = w meas . The optimal value δ opt is then identified solving the following optimization problem,

δ opt = arg max δ {L W (w meas ; δ)} . ( 24
)
The quantity p W (w meas ; δ) is computed using independent realizations of W generated with the adapted stochastic model and fitted by using the multivariate Gaussian kernel method (see for instance [START_REF] Terrell | Variable Kernel Density Estimation[END_REF][START_REF] Bowman | Applied Smoothing Techniques for Data Analysis[END_REF]).

Adapted stochastic modeling for the long time evolution of a given track portion

For a given track portion, it is assumed that ν τ +1 measurements x meas τ 0 , . . . , x meas τν τ corresponding to long times τ ∈ {τ 1 , . . . , τ ντ } are available. The optimal value δ opt of δ is identified using only time τ 0 (the first time).

δ opt = arg max δ {L W (w meas τ 0 , δ)} . ( 25
)
It is then assumed that this optimal value is representative to construct the adapted stochastic modeling allowing the long time evolution of this given track portion to be estimated.

The admissible set for δ is defined as the domain [0, 2] 4 of R 4 . For estimating the optimal value, the following deterministic search algorithm is used. The identification of δ opt is made in two steps. The admissible set is partitioned in 20 4 meshes corresponding to a constant step of 0.1 for each coordinate. For each value of δ corresponding to a node of the mesh, the loglikelihood is computed with 10, 000 realizations of W(δ). This stage allows the node of the mesh corresponding to the maximum of the log-likelihood to be identified and this node is then denoted as

δ 1 = (δ 1 1 , δ 1 2 , δ 1 3 , δ 1 4
). Since the dimension of the admissible set of parameter δ is small, and since each component δ k belongs to a given finite interval, a deterministic algorithm can be used for solving the optimization problem defined by Equation (25). Then, the subdomain Π 4 k=1 [δ 1 k -0.25, δ 1 k + 0.25] is explored for δ with a precision of 0.05 for each coordinate (10 4 meshes for the subdomain). The log-likelihood is computed with 100, 000 realizations for W(δ). The values obtained for each coordinate are coherent with Figure 2. The variation of the marginal probability density function (PDF) w 1 → p W 1 (δ 1 ) (w 1 ; δ 1 ) of random variable W 1 (δ 1 ) is plotted in Figure 3 The adapted stochastic modeling of the long time evolution for the given track portion is constructed as follows. At long time τ , the measurement x meas τ is projected on the global stochastic model defined by Equation [START_REF] Soize | Stochastic Models of Uncertainties in Computational Mechanics[END_REF]. Taking into account Equations ( 10) and ( 13) yields the projection η meas τ which is written as

η meas τ = [Q] T [λ] -1 (x meas τ ) . ( 27
)
Then, this projection is used in the adapted stochastic model defined by Equation ( 17) and gives

X k τ (δ opt k ) = [Q k ] η meas τ + δ opt k G k , k = 1, 2, 3, 4 . ( 28
)
As an illustration, x 1,meas τ 0 and the confidence region at 95% of X 1 τ 0 (δ opt 1 ) are compared in Figure 4. It can be noticed that the geometrical and physical properties of the irregularities are preserved with the identified adapted stochastic modeling.

Dynamic response of the train

Stochastic dynamic response of the train

For fixed long time τ , the stochastic model of the track irregularities is given by X k τ (δ opt k ), k = 1, 2, 3, 4, defined by Equation (28). The stochastic response of the train is then computed using the Monte-Carlo method. For each realization X k τ (δ opt k ; θ ) for θ ∈ Θ , the deterministic realization of the train response is computed. The set of these realizations allows statistical estimators to be constructed for analyzing the stochastic responses through dynamic indicators. The calculations are performed with 100 realizations and the responses are expressed in terms of the curvilinear abscissa s = v t, in which v is the constant speed of the train and t is the time of the train dynamics. The train is a TGV Duplex for which the number of wheelsets is ν w = 26, and wheelsets are numbered j = 1, . . . , ν w .

curvilinear abscissa s (m) cant irregularity (m) 5 % of X 1 τ0 95 % of X 1 τ0 x meas 1,τ0

Description of the dynamic indicators

The chosen dynamic indicators are based on criteria described in norm UIC 518 [START_REF] Uic | Testing and approval of railway vehicles from the point of view of their dynamic behaviour -Safety -Track fatigue -Running behaviour[END_REF] for the homologation of railway vehicles. These indicators C k j are related to each wheelset, j, and are functions of the dynamic outputs A k j (s) which are themselves functions of index R which stands for "right wheel" and of index L which stands for "left wheel", and of curvilinear abscissa s. The dynamic outputs are:

• the sum of the lateral forces on each wheelset, A 1 j = (Y j,R + Y j,L )(s);

• the maximum between right and left wheels of the quotient of transverse forces Y on vertical load Q at wheel-rail contact, which is defined by the quantity

A 2 j = max{(Y j,R /Q j,R )(s), (Y j,L /Q j,L )(s)};
• the lateral accelerations, A 3 j = Ÿj (s);

• the vertical accelerations, A 4 j = Zj (s); • the sum of the right and the left vertical forces on each wheelset, A 5 j = (Q j,R + Q j,L )(s);

• the difference between the right and the left vertical forces on each wheelset, A 6 j = (Q j,R -Q j,L )(s). According to [START_REF] Uic | Testing and approval of railway vehicles from the point of view of their dynamic behaviour -Safety -Track fatigue -Running behaviour[END_REF], these indicators are smoothed and filtered. The indicator C k j is defined by

C k j = max s∈[0,S] |A k j (s)| , 1 ≤ j ≤ ν w , 1 ≤ k ≤ 6 . ( 29
)
The number of indicators is reduced with respect to the number of wheelsets in analyzing the correlation between them. The indicators which are strongly correlated are removed. For instance, the correlation matrix for the sum of lateral forces restricted to 15 wheelsets among the 26 is plotted in Figure 5. Finally, 7 indicators on 6 wheelsets are kept, which yields 42 indicators to assess the train dynamic response.

Long time evolution of the dynamic indicators

We are now able to assess the train dynamic response on a given track portion and to observe its long time evolution as a function of τ using the dynamic indicators C k i (τ ), 1 ≤ i ≤ 6, 1 ≤ k ≤ 7. For each time τ , the probability density function (PDF) of C k i (τ ) is plotted in Figure 6. Its long time evolution shows an increase and a dispersion of the safety and comfort indicators of the train dynamics.

Conclusion and perspectives

This work corresponds to a first step of the development of a stochastic modeling allowing the long time evolution of a given track portion to be analyzed. The stochastic model which is proposed has been constructed using a very large experimental data basis related to all the french railway network for the high speed trains, and consequently, can be considered as very robust. The next steps which are in progress consist in introducing uncertainties in the dynamical model of the train and in developing a stochastic model for the long time evolution.

2 /s 3 ) wave number k (m - 1 )Figure 1 :

 2311 Figure 1: Power spectral density functions of two outputs of the simulated train dynamic response at a given location in the train.

4 Figure 2 :

 42 Figure 2: Sections of the log-likelihood L W (w meas τ0 , δ) for each component δ k of δ

as a function of δ 1 . 1 variation of δ 1 w meas 1 Figure 3 :

 1113 Figure 3: Graphs the PDF, w 1 → p W1(δ1) (w 1 ; δ 1 ), in function of δ 1 (the bold line is obtained for δ 1 = δ 1,opt ).

Figure 4 :

 4 Figure 4: Irregularity x 1,meas τ0and confidence region at 95% of X1 τ0 (δ opt 1 ) for the given track portion.

Figure 5 :

 5 Figure 5: Correlation between the indicators on the sum of lateral forces for 15 wheelsets.

Figure 6 :

 6 Figure 6: Evolution of the PDF of the dynamic indicators C 1 4 and C 3 4 as a function of time τ .

Acknowledgement

This work was supported by SNCF.