L De Abreu Corrêa

R Cottereau
email: cottereau@lma.cnrs-mrs.fr

B Faure

Dispersion analysis in ballasted railway tracks and Anderson localization in granular media

Keywords: Ballasted railway tracks, Granular media, Heterogeneous media, Random media, Anderson localization 00-01, 99-00

This paper analyzes the dynamic behavior of a randomly-heterogeneous continuum model of a ballasted railway track through the numerical construction of dispersion curves along the axis of the track. At intermediate frequencies (50-100 Hz), the dispersion curves display a strong decrease of the group velocity. At higher frequencies (above 100 Hz), the wave are strongly localized and cannot propagate along the track. This behavior is coherent with experimental observations reported in the literature on ballasted railway tracks and other granular media. Additional simulations are performed to prove that the localization is not a consequence of the geometry of the track or the impedance mismatch between the ballast and the underlying soil, but is indeed Anderson localization due to the heterogeneity of the ballast.

Introduction

The ballast is the uppermost layer of the railway track superstructure, on which the sleepers stand. It is made of coarse crushed stones, whose size distribution lies in the range 25 -50 mm. This ballast layer plays an important role in the transmission and repartition of static and dynamic train loads, in the absorption of mechanical and acoustical vibrations, and in the drainage of rain water [START_REF] Indraratna | Advanced rail geotechnology. Ballasted track[END_REF]. The noise of the passage of a train and its vibration impact on the environment can be important for all types of trains: heavy freight coaches, whose low frequency impact in the soil carries far away; highspeed trains, that generate high amplitude excitations; and even tramways (although more rarely on ballasted tracks), which stand very close to surrounding buildings. This vibration impact may sometimes be strong enough to force the train companies to locally reduce the velocity of the passing trains or force the companies change the track.

The modeling of the ballast aims at improving the understanding of the mechanical behavior of the track system and providing clues to mitigate this issue [START_REF] Heckl | Structure-borne sound and vibration from rail traffic[END_REF][START_REF] Connolly | Benchmarking railway vibrations -track, vehicle, ground and building effects[END_REF]. It represents a challenging task because granular materials exhibit various distinct behaviors (such as transition between fluid-like and solid-like behavior) depending on the applied stresses and strains.

Classical numerical models of granular media include discrete models and continuum homogeneous models. Unfortunately, they are both inappropriate for dynamic analyses of ballasted railway tracks. Indeed, the latter produce vibration levels away from the track that are much larger than experimentally measured, while the former are too numerically expensive for such applications. Recently, a randomly-fluctuating heterogeneous continuum model of the ballast was proposed [START_REF] De Abreu Corrêa | Randomly-fluctuating heterogeneous continuum model of a ballasted railway track[END_REF]. This model is based on continuum mechanics, so that it is parameterized with Lamé parameters and density, which are constant within the domain in classical models. In the proposed model, however, they fluctuate at the scale of the grains, which allows to recover the strong heterogeneity of the stress field, characteristic of granular media and very influential on its dynamic behavior [START_REF] Leibig | Model for the propagation of sound in granular materials[END_REF][START_REF] De Abreu Corrêa | Randomly-fluctuating heterogeneous continuum model of a ballasted railway track[END_REF]. Finally, the fluctuations of the mechanical parameters are modeled as random fields and parameterized with statistical parameters (average, variance, correlation length). More details on that model will be recalled in Section 3.1.

The model was shown to be coherent with previously unexplained experimental observations, and to reproduce stress distributions obtained by discrete simulation of assemblies of grains. Being continuum-mechanics based, it is amenable to very efficient numerical approximation of large scale problems involving ballasted railway tracks.

The main objective of this paper is to analyze precisely the dynamical behavior of this randomly fluctuating model and confirm its validity for ballasted railway tracks. This analysis will be performed through the numerical construction of dispersion curves for that model of the ballast.

Dispersion curves are among the most classical tools used in the analysis of wave propagation. These curves link the values of pairs of wave numbers and frequencies (k, ω) for which a plane wave can actually propagate in the medium. These simple 2D diagrams allow for a quick understanding of the properties of the wave propagation medium, independently of the source that may excite it. For a given excitation frequency, the curves chart which waves will propagate in the medium considered, and with what wave number. They also allow to understand the way a given wave will evolve with frequency, the so-called dispersion characteristics of that wave. In some rare cases, the dispersion curves can be constructed analytically. For instance, in linear elastic isotropic infinite media, the dispersion curves are two straight lines, for the pressure and shear waves, respectively. For half-spaces and thin media, respectively, the dispersion curves for Rayleigh waves [START_REF] Rayleigh | On waves propagated along the plane surface of an elastic solid[END_REF] and Lamb waves [START_REF] Lamb | On waves in an elastic plate[END_REF] can likewise be constructed analytically. In three dimensions, the wave number is a vector, with each coordinate describing propagation along a given direction of space, and the dispersion curve is actually a volume in a 4D space. To simplify the presentation, this paper concentrates on the construction of the dispersion relation for wave numbers along the axis of the track, hence exposing classical plot, with frequency on one axis, and a scalar wave number along the other axis.

In granular media, among which the ballast, wave propagation phenomena are more difficult to analyze [START_REF] Langer | Lectures in the theory of pattern formation[END_REF][START_REF] Duran | Sands, powders, and grains: an introduction to the physics of granular materials, Partially Ordered Systems[END_REF]. For instance, dispersion curves can be constructed experimentally for waves at the surface of unconsolidated granular media in dry condition [START_REF] Jacob | Acoustic probing of the jamming transition in an unconsolidated granular medium[END_REF][START_REF] Bodet | Elasticity profile of an unconsolidated granular medium inferred from guided waves: Toward acoustic monitoring of analogue models[END_REF] and present a very complex behavior. Waves in granular media are carried by complex force chains, which create a heterogeneous and non-isotropic behavior, as observed with photo-elastic imaging [START_REF] Shukla | Dynamic photoelastic studies of wave propagation in granular media[END_REF][START_REF] Geng | Footprints in sand: the response of a granular material to local perturbations[END_REF]. Also, the contact chains are easily broken, deformed and rearranged [START_REF] Liu | Sound and vibration in granular materials[END_REF], creating a nonlinear behavior of the medium.

Even in the purely linear case, granular media are also dispersive, which is difficult to capture through classical continuum models [START_REF] Buckingham | Theory of acoustic attenuation, dispersion, and pulse propagation in unconsolidated granular materials including marine sediments[END_REF]. Through the construction of numerical dispersion curves for the randomly-fluctuating continuum model of the ballast, this paper will reproduce some of the features of these experimental dispersion curves.

After the presentation of a general methodology for the numerical construction of dispersion curves in Section 2, this construction will be the focus of Section 3. There, it will be shown that high-frequency waves (above 100 Hz) cannot propagate in this ballasted railway track. This will be analyzed as Anderson localization [START_REF] Anderson | Absence of Diffusion in Certain Random Lattices[END_REF][START_REF] Anderson | The question of classical localization: a theory of white paint ?[END_REF], which is a very specific physical phenomenon that appears in disordered media, such as granular packings. When localized, the waves cannot propagate over long distances because of the heterogeneity of the parameters and remain in the vicinity of the source. Although more classical in quantum systems, experimental evidence of Anderson localization exist for classical waves [START_REF] Weaver | Anderson localization of ultrasound[END_REF][START_REF] Liu | Sound in sand[END_REF][START_REF] Lobkis | Anderson localization of ultrasound in plates: Further experimental results[END_REF], and in particular in granular media [START_REF] Hu | Localization of ultrasound in a three-dimensional elastic network[END_REF]. To make sure that the localization observed in Section 3 is not a simple wave guide effect, or a consequence of the impedance mismatch between the ballasted railway track and the underlying soil, additional simulations will be performed in Section 4.

Construction of the dispersion curve using numerical solutions in time-space

In this section, we present a methodology to obtain dispersion curves using numerical solutions of the wave equation in space-time.

Definition of dispersion curves

Let us consider the equilibrium equation, neglecting the influence of the gravity forces, states:

∇ • σ + f = ρ ü (1)
where u(x,t) is the displacement field, σ represents the stress tensor, ε = (∇u + ∇u T )/2 is the strain, f represents external volume forces, and ∇ is the differential operator. We assume that the material is heterogeneous, linear, and isotropic, so that σ = λ (x)TrεI + 2µ(x)ε with λ (x) and µ(x) the Lamé coefficients, and I the identity second-order tensor. The density is denoted ρ(x). The parameters are related to the wave velocities through V p = (λ + 2µ)/ρ and V s = µ/ρ. This problem is closed with boundary and initial conditions. For instance, in an infinite medium, the Sommerfeld radiation conditions state:

lim r→∞ r ∂ u ∂ r -iku = 0, (2)
where r = x .

Assuming that the displacement field in Eq. ( 1) is expressed as a monochromatic plane wave:

u = ũ exp i(k • x -ωt),
where ω is the angular frequency, k the wave number and ũ the amplitude vector, we obtain in the case with no volume force:

(Γ(k) -ρω 2 I) ũ = 0 (3)
where

Γ(k) = (λ + µ)k ⊗ k -i(∇λ ⊗ k + k ⊗ ∇µ) + (µ|k| 2 -ik • ∇µ) is the Christoffel tensor.
The dispersion equations are defined as det(Γ(k)ρω 2 I) = 0 and indicate the (k, ω) pairs that allow for the existence of non-vanishing plane wave solutions of the wave equation. The values of k can be real, complex or imaginary. With the convention used here (as in [START_REF] Auld | Acoustic fields and waves in solids[END_REF]), the real values correspond to propagative modes, pure imaginary values refer to evanescent modes, and complex values correspond to propagative modes with attenuation.

For an elastic, isotropic and 3-D homogeneous medium, the propagative solutions are such that ω = V p |k| or ω = V s |k|. In the latter case, two different polarizations (direction of vector ũ) exist. For most other media, dispersion relations have to be computed numerically [START_REF] Lagasse | Higher order finite element analysis of topographic guides supporting elastic surface waves[END_REF][START_REF] Knopoff | A matrix method for elastic wave problems[END_REF].

General methodology to obtain dispersion curves using time-space data

A first approach consists in diagonalizing directly the Christoffel tensor, when available. In this paper, because the problem is too large, we rather use a methodology based on a 2D Fourier transform of numerical approximations of Eq. ( 1) to obtain an approximation of the pairs (k, ω) that satisfy the dispersion relation [START_REF] Alleyne | A two-dimensional fourier transform method for the measurement of propagating multimode signals[END_REF][START_REF] Jacob | Acoustic probing of the jamming transition in an unconsolidated granular medium[END_REF][START_REF] O'donovan | Micromechanics of wave propagation through granular material[END_REF]. The general methodology to construct the dispersion relation consists in computing a displacement field, u(x,t) and performing Fourier transforms in both the time and space domains to obtain a 2D-variate function

H k(k, ω) in a specific direction k = k/k with k = k : H k(k, ω) = +∞ -∞ +∞ -∞ u(x,t)e -i(k•x+ωt) dx k dt. ( 4
)
Note that x k = x • k represents the coordinate of the space position along the selected direction. The locations in (k, ω) where the amplitude H k(k, ω) is large match the dispersion equation [START_REF] Alleyne | A two-dimensional fourier transform method for the measurement of propagating multimode signals[END_REF]. Note that advantage is taken of the property of the method to resolve multimodal signals (more than one wave mode in the same time signal).

In order to compute an approximation of Eq. ( 4), we use u(x,t) estimated in positions {x m } 1≤m≤M and aligned along direction k, and record the quantity of interest in those sensors at different time steps {t n } 1≤n≤N . Assuming that the sensors are uniformly distributed along a segment of length L, and that the time instants are uniformly distributed in an interval of length T , we obtain the following map:

H pq = M ∑ m=1 N ∑ n=1 w(x m • k,t n )u(x m ,t n )e -2πi( pm M + qn N ) = DFT 2D [w(x m • k,t n )u(x m ,t n )] (5)
where p = kL is the normalized wave-number, q = f T is the normalized frequency and DFT 2D represents the 2D Discrete Fourier Transform. To avoid the presence of side lobes and reduce leakage, a bi-dimensional Tukey (tapered cosine) window function w(x k ,t) was applied [START_REF] Harris | On the use of windows for harmonic analysis with the discrete fourier transform[END_REF][START_REF] Nuttall | Some windows with very good sidelobe behavior[END_REF].

In many examples described in this paper, the maps cited above are computed as averages of such functions using P lines of sensors, and N MC realizations of the random medium (to be discussed later). In that case, the formulas above are computed for each line and each realization independently, and the average is computed dividing by P × N MC . Also note that the source used for the construction of the maps has a very flat spectrum so the solution u(x,t) can be used directly. However, all results will be de-convolved (actually normalized by the source in the Fourier space of frequencies) before computing the maps.

Spectral Element solver

In this paper, Eq. ( 1) is approximated with the Spectral Element Method (SEM).

The SEM is a high-order Finite Element Method that uses Gauss-Lobatto-Legendre (GLL) quadrature rule and Lagrange polynomials based on the nodes of that same quadrature. This ensures that the mass matrix is diagonal and allows to use explicit time integration schemes and very efficient parallelization (for a complete description of the method, see for instance [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF], as well as the original papers [START_REF] Maday | Spectral element methods for the incompressible Navier-Stokes equations[END_REF][START_REF] Bernardi | Approximations spectrales de problèmes aux limites elliptiques[END_REF]). The space discretization of the variational form of the wave propagation Eq. ( 1) gives:

M V = F ext -F int (U, V) (6)
where U and V are vectors containing the components of the displacement and velocity at the nodes, M is the (diagonal) mass matrix, and the vectors F ext and F int are vectors of external loads and internal forces. Using an explicit second-order finite-difference scheme in time leads to:

1 ∆t M [V n+1 -V n ] = F ext n+1/2 -F int (U n+1/2 , V n+1/2 ), (7)
U n+1 = U n + 1 2∆t [V n + V n+1 ] , (8)
A n+1 = 1 ∆t [V n+1 -V n ] , (9)
where

U n+1/2 = (U n+1 + U n )/2 and F ext n+1/2 = (F ext n+1 + F ext n )/2.
Even though the stability condition requires to use very small time steps [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF][START_REF] Sevilla | Influence of periodically fluctuating material parameters on the stability of explicit high-order spectral element methods[END_REF][START_REF] Cottereau | Stability of an explicit high-order spectral element method for acoustics in heterogeneous media based on local element stability criteria[END_REF], the construction of the solution at each time step is very cheap because the inversion of the mass matrix is instantaneous. Our implementation uses hexahedral non-structured meshes and Perfectly Matched Layers (PML) to account for infinite domains. PMLs are fictitious layers of materials added at the outer edges of the computational domain to absorb outgoing waves with as little reflection as possible at the interface [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Collino | Optimizing the perfectly matched layer[END_REF][START_REF] Festa | The Newmark scheme as velocity-stress time-staggering: An efficient PML implementation for spectral element simulations of elastodynamics[END_REF][START_REF] Modave | Optimizing perfectly matched layers in discrete contexts[END_REF]. All simulations presented in this paper ran at Moulon Mesocentre facility in France, jointly managed by CentraleSupélec and ENS Paris-Saclay. The cluster uses Intel Xeon processors E5-2670 v3 @ 2.30 GHz (Haswell).

Verification of the methodology in the infinite isotropic case

This last part of the section is devoted to a verification of the methodology introduced in Section 2.2 in the isotropic case, for which the dispersion curve is known. Indeed, as indicated in Section 2.1, the dispersion curve in infinite isotropic media is composed of two lines, with equations ω = V p |k| and ω = V s |k|. This last section given by an Ormsby function:

A(t) = (π f 4 sinc(π f 4 t)) 2 π( f 4 -f 3 ) - (π f 3 sinc(π f 3 t)) 2 π( f 4 -f 3 ) - (π f 2 sinc(π f 2 t)) 2 π( f 2 -f 1 ) - (π f 1 sinc(π f 1 t)) 2 π( f 2 -f 1 ) , (10)
with The mesh is composed of hexahedral elements of edge λ min /2 = 0.3125 m. Each 135 element uses 4th-order Lagrange polynomials in each direction, for a total of 375 de- grees of freedom per element and 4'096'000 elements in the mesh. On all sides, 5 layers of PMLs are added, with a total thickness 0.255λ max . The damping evolution in the PMLs is parameterized by a polynomial degree p = 3 and amplitude α = 9.87 [START_REF] Drozdz | Efficient finite element modelling of ultrasound waves in elastic media[END_REF].

f 1 = 20
The model was run on 144 processors, and required around 140 h of CPU time.

The methodology proposed in Sec. 2.2 was then applied. In order to evaluate Eq. ( 5), sensors are placed along different lines and record the displacement field.

Three groups of 5 lines of sensors were placed aligned with the cube axes. Each line has a total length of 3.5λ max , with equally spaced sensor by λ min /10. Finally, the distance between two parallel lines is λ min /2. The general setting of sensors is represented on Fig. 2. Fig. 3 presents the time history in the numerical domain. As the PML domain is fictitious, its local solution is not relevant to the mechanical problem and is not displayed in the images. There is no apparent reflection coming out of the PML, which indicates that it is accurately playing its role of absorbing outgoing waves. We can clearly see the P and S waves with different velocities. The P wave, faster then the S wave, presents large amplitudes along the x-axis and vanishing amplitude along the y and z axes. This is due to the polarization of the source along the x-axis. As the P-wave is a pressure wave, it is indeed strongly excited by the source. On the other hand, such a source does not generate pressure perpendicular to its polarization. For the S-wave, snapshots of the solution in time and analyze them. Finally, the construction of the dispersion curve itself is presented and discussed.

Randomly heterogeneous continuum model of a ballasted railway track

The ballasted railway track is modeled over a length of 84.375 m as an elastic halfspace topped by a wedge-shaped structure, whose geometry is invariant along one direction (see Fig. 5). The structure represents the ballast while the elastic half-space represents the soil that supports it. The soil is modeled as a bounded box with dimensions Note that this model was already used in [START_REF] Innotrack D | First phase on the modelling of poor quality sites[END_REF], although with a different mesh.

A series of point forces polarized in the [1; 1; 1] direction, are positioned at the top of the ballast layer, 9 m away from the border of the model in a non-symmetric arrange- 10)). The displacement is recorded in sensors located along six parallel lines at the top of the ballast layer and aligned with the axis of the track (see Fig. 5). The lines are 72 meters-long and the sensors are separated by 0.0732 m ≈ λ min /10. This sampling is sufficient to be able to compute the Fourier transform in space accurately.

The ballast is modeled as a continuum-based randomly-fluctuating elastic medium, as introduced in [START_REF] De Abreu Corrêa | Randomly-fluctuating heterogeneous continuum model of a ballasted railway track[END_REF]. Following the polar decomposition of the strain-stress tensor [START_REF] Ta | Modeling of random anisotropic elastic media and impact on wave propagation[END_REF]), the density, ρ, and the shear, µ, and bulk moduli, κ, are taken as fully correlated stationary random fields indexed in space. The first-order marginal density for each of these random fields is taken as gamma law with averages ρ = 1570 kg/m 3 , κ = 179.60 MN/m 2 and µ = 35.32 MN/m 2 , and standard deviations σ ρ = 3925 kg/m 3 , σ κ = 449.02 MN/m 2 and σ µ = 88.31 MN/m 2 . This corresponds to a ballast layer with (arithmetic) average velocities V s = 150 m/s and V p = 380 m/s [START_REF] Innotrack D | First phase on the modelling of poor quality sites[END_REF]. In such heterogeneous setting, the apparent velocities of the homogenized waves are expected to be close to the harmonic averages, which are here approximately V p = 280 m/s and V s = 110 m/s for the pressure and shear waves, respectively. Note that the standard deviations are equal to 2.5 times the averages, so that the fluctuations are very large. The correlation model [START_REF] Torquato | Microstructure of two-phase random media. III. The n-point matrix probability functions for fully penetrable spheres[END_REF][START_REF] Quintanilla | Microstucture functions for random media with impenetrable particles[END_REF][START_REF] Torquato | Random heterogeneous materials. Microstructure and macroscopic properties[END_REF] corresponds to that of a dense packing of impenetrable spheres of uniform diameter d = 3.9 cm in void, with volume ratio φ = 0.583. More details, including identification and comparison with simulations with the discrete element method, can be found in [START_REF] De Abreu Corrêa | Randomly-fluctuating heterogeneous continuum model of a ballasted railway track[END_REF]. Realizations of the random fields can be obtained using a spectral representation method, specifically tailored for large scale parallel applications [START_REF] De Carvalho Paludo | Scalable parallel scheme for sampling of Gaussian random fields over large domains[END_REF], and made available as an open-source software [START_REF]Random field library[END_REF]. One sample of the random field is provided in Fig. 6.

Analysis of the wave patterns

Fig. 7 displays snapshots at different times of the normalized displacement field in the ballasted railway track, with a normalization similar to that of Fig. 3, that is to say a normalization by the maximum displacement at time t = 0.25 s. The first observation is that a large part of the energy seems to remain localized within the vicinity of the source. This is not what is expected from a simulation in a homogeneous domain [START_REF] De Abreu Corrêa | Impact of the heterogeneity of the ballast on the dynamical behavior of the ballast-soil system[END_REF][START_REF] De Abreu Corrêa | Randomly-fluctuating heterogeneous continuum model of a ballasted railway track[END_REF], where waves should propagate non-dispersively with the velocity of a Rayleigh wave. In a setting including a homogeneous ballast layer, the situation would be more complex, and probably involve dispersive waves in the waveguide [START_REF] De Abreu Corrêa | Impact of the heterogeneity of the ballast on the dynamical behavior of the ballast-soil system[END_REF]. How-ever, the fact that part of the energy remains localized close to the source is a specific feature of the heterogeneous model of the ballast. On the snapshots at longer times (at t = 0.585 s for instance), there seems to be a separation between waves with longer wavelengths that actually propagate away from the source, and shorter wavelengths that remain localized around the source. This effect therefore seems to be dependent on frequency. A possible explanation for this phenomenon, that seems to be stationary in time, and strongly selective in frequency, is Anderson localization [START_REF] Anderson | Absence of Diffusion in Certain Random Lattices[END_REF][START_REF] Leibig | Model for the propagation of sound in granular materials[END_REF][START_REF] Lawney | Frequency filtering in disordered granular chains[END_REF]. First discovered in 1958 by P. W. Anderson [START_REF] Anderson | Absence of Diffusion in Certain Random Lattices[END_REF] for the Schrödinger equation, it explains why pure metals are conductors while strongly heterogeneous metals -that is to say metals with high density of defects in the cristalline structure -are insulators. Later on, it was acknowledged as ubiquitous to all classical wave equations (acoustic, electromagnetic, elastic) [START_REF] Anderson | The question of classical localization: a theory of white paint ?[END_REF], and explains that energy cannot propagate away from the source in strongly heterogeneous non-periodic media. Anderson received the Nobel Prize in physics in 1977 for that discovery. The absence of propagative features can be observed in the dispersion relation through the absence of (ω, k) pairs in the considered frequency range, that is to say a disappearance of the dispersion curve. Such a disappearance will be discussed in the next section.

Construction of the dispersion curves

We now turn to the construction of the dispersion curves for our model, using the methodology presented in Sec. 2.2 and based on displacement field acquired along the sensors lines indicated in Fig. 5. The maps H pq are obtained for each of the 6 lines independently and for 5 different random samples of the ballast layer. The final map is obtained as an average over these 30 individual maps and is plotted in Fig. 8.

In this figure, three main frequency regimes are identified. Above about 100 Hz, the figure appears blank, which means that no energy propagates. Of course, this does not mean that no energy exists in that frequency range, since the source contains frequencies up to 400 Hz. It just does not propagate along the axis of the track (the direction that is being monitored here), which is consistent with the observations at short wavelengths made on Fig. 7 in the previous section. Energy in that high-frequency range actually remains mostly localized in the vicinity of the source, as predicted by Ander- . son localization. Below about 50 Hz, there seems to be mainly an S-wave (or more probably a Rayleigh wave, whose velocity is slightly smaller) propagating in the soil.

At such long wavelengths (compared to the thickness of the ballast), the wave does not really interact with the ballast layer and propagates as if in a homogeneous soil without ballast. Between about 50 Hz and 100 Hz, the soil S/Rayleigh wave starts interacting with the ballast, and its apparent velocity is strongly reduced, approximating that of the ballast layer. A P wave also seems to appear in the ballast (faster than in the soil), bending rapidly with increasing frequency. This general behavior is consistent with experimental observations in unconsolidated granular packings [START_REF] Jacob | Acoustic probing of the jamming transition in an unconsolidated granular medium[END_REF].

To get an idea of the fluctuations among the different independent maps (before averaging), the maps corresponding to two different lines for the first realization of the random ballast are plotted in Fig. 9, along with the maps corresponding to the average over the 6 lines for two different random realizations of the ballast layer. The few maps presented here are representative of the entire set. Although the average on the lines and the realizations reduces noise, the main features discussed above are present in each and every individual curve. These features are therefore not a chance effect obtained from averaging but really a characteristic behavior of ballasted railway tracks.

Figs. 8 and9 clearly state that no wave can propagate along the direction of the track with frequency above 100 Hz (wave number above 1 m-1), and that waves with frequency between 50 Hz and 100 Hz propagate more slowly than wave below 50 Hz.

This of course reads more precisely in the dispersion curves but was already qualitatively seen at later times on Fig. 7. What cannot be read on the dispersion curves is what happens to the energy of these higher frequency waves since there is no damping in this model. Observation of Fig. 7 indicates that at least part of this energy remains in the vicinity of the source. Although it is widely known that high frequency waves get dissipated within the ballast, the most novel finding here is a proposal for the mechanism that leads to the trapping of high frequency waves in that ballast: Anderson localization. Note that Anderson localization does not involve any non-linear mechanism, so that it is advocated that energy becomes trapped in the ballast layer by heterogeneity, before it is eventually dissipated (in particular by friction between grains, which is not modeled here).

Strong localization analysis

In the previous section, appearance of Anderson localization in ballasted railway tracks was hinted at, with its classical stationary character and a clear frequency dependence. However, there is a possibility that this localization might be some type of guided wave effect, due to the particular geometry of the ballast. This section aims at removing the ambiguity by showing that the same feature would be present in an infinite domain, composed only of ballast material.

Description of the numerical model

To remove the effect of the ballast geometry and its interface with the underlying soil, we consider an infinite domain, as in Section 2.4. This domain is modeled as a cube of edge 10.98 m, with PMLs all around. The mesh is discretized with a uniform grid of ≈ 7.

Analysis of the wave patterns

For each of the two cases considered, Figs. 11 and 12 display snapshots of the displacement magnitude at different time instants. In the homogeneous case, the expected symmetrically spherical waves are observed. As in Section 2.4, P and S waves are observed with different velocities, and with amplitudes depending on the orientation with respect to the polarization of the source. On the other hand, in the heterogeneous cases, there is no clear coherent pulse propagating away from the source and the energy remains localized in the vicinity of the source.

Dispersion curve for an infinite medium of ballast material

Eq. ( 5) was evaluated as an average over the 91 lines monitored for both the homogeneous and heterogeneous cases and Fig. 13 presents the obtained maps H pq . The reference value for the decibel scale was picked as the maximum value from the homogeneous case for both maps. As expected, a good agreement is obtained in the homogeneous case with both the analytical shear wave and pressure wave modes. The amplitude of the dispersion curve for the P wave appears smaller because of the average among 91 lines who are mostly oriented to record S-polarized displacements. Indeed, only 1 line is perfectly oriented to record the P wave (along the axis of polarization of the source), while 10 lines are perfectly oriented to record S waves (in the plane going through the source and normal to its polarization).

Turning to the heterogeneous model, the filtering effect at 100 Hz is clearly observed as in the ballasted railway track of the previous section. This shows that it is indeed an intrinsic property of the heterogeneous material, and not a wave guide effect due to the geometry of the ballasted track. Below 100 Hz, however, the material behaves as a purely isotropic material, which shows that the velocity reduction observed in the ballasted railway track was indeed the consequence of the combination of ballast and soil at different wavelengths. Note that there is no apparent P-wave in this heterogeneous model for the same reason as discussed in the homogeneous case, but amplified by the heterogeneity and the weaker averaging for lines recording P waves.

Characterization of Anderson localization

In this last section, we provide additional evidence that the filtering effect above 100 Hz is indeed Anderson localization. Anderson localization is characterized by its strong frequency dependence, which is already apparent in the dispersion curves constructed previously. It is also characterized by the stationarity in time of the localized energy, which means that the energy remains in the vicinity of the source rather than slowly diffusing away. To illustrate this, the displacement magnitude is therefore averaged over all recording stations at the same distance from the source and plotted in Figs. 14 and 15 as a function of time and distance from the source, filtered in three frequency bands (using a second order bandpass Butterworth filter). The normalization factor for the decibel scale is the maximum displacement magnitude for complete spectra of frequencies. In the homogeneous case, the energy clearly propagates along a straight line, whose slope is controlled by the velocity of the medium. As expected in a non-dispersive medium, this effect is independent of frequency. On the other hand, in the heterogeneous medium, the behavior is strongly dispersive. Below 50 Hz, an isotropic homogenized behavior is retrieved. Above 100 Hz, the energy seems to propagate along vertical lines (zero velocity), which clearly indicates that it remains localized in time at the same distance from the source.

Conclusion

In this paper, we constructed and analyzed the dispersion curve of a ballasted railway track set on a homogeneous soil. A continuum-based randomly-fluctuating heterogeneous model was used to model the ballast. The dispersion curve of the ballasted railway track shows a tremendous influence of the heterogeneity. The behavior is dispersive with higher frequency waves propagating slower than lower frequency waves.

Below 50 Hz, the waves do not interact with the ballast layer and we retrieve the behavior of the homogeneous soil. In the frequency range between 50 Hz and 100 Hz, waves seem to slow down, but still propagate at that slower pace. Above 100 Hz, Anderson localization sets in, where the higher frequency waves remain trapped in the vicinity of the source. Although railway engineering applications would really consider only the [0, 150] Hz range, this paper considered a larger range to draw definite conclusions about the lack of propagation of the higher frequency waves.

Although very widely studied in the physics literature, Anderson localization has been studied in the acoustical literature mostly from the experimental angle, trying to discriminate between dissipation and localization in measurements [START_REF] Hu | Localization of ultrasound in a three-dimensional elastic network[END_REF]. The conclusions of this paper seem to indicate that the higher part of the energy generated by the passage of trains on a ballasted railway track remains trapped in the ballast rather than being dissipated in the surrounding soil. This is opposite to the common knowledge in the field of numerical modeling of railway-generated noise in the environment [START_REF] Paolucci | Numerical prediction of low-frequency ground vibrations induced by high-speed trains at[END_REF],

although railway engineers also experience strong ballast degradation which indicates large dissipation of energy in that layer. This paper therefore proposes a physical explanation for this experimental observation, which may be a first building block towards improving ballasted railway tracks.

Figure 1 :

 1 Figure 1: Time history of the amplitude of the point source (a) and frequency spectrum (b) for the verification example in Section 2.4.

  Hz, f 2 = 25 Hz, f 3 = 160 Hz, and f 4 = 165 Hz. This function displays a quasi-plane spectrum between frequencies f 2 and f 3 , with a low-cut frequency f 1 and a high-cut frequency f 4 . Fig. 1 displays the time history and frequency spectrum of the loading.

Figure 2 :

 2 Figure 2: General setting of the verification example in Sec. 2.4: Monitors positions (5 blue lines along each direction of space), source position (red dot). The exterior cube represents the numerical domain without the PMLs.

Figure 3 : 3 .Figure 4 :

 334 Figure 3: Normalized displacement fields in an infinite elastic isotropic homogeneous model at times: (a) t = 0.000 s; (b) t = 0.141 s; (c) t = 0.282 s; (d) t = 0.423 s; (e) t = 0.564 s; (f) t = 0.705 s.

3 × 8 ×

 8 84.375 m3 and with the following properties: V s = 180 m/s, V p = 1100 m/s, and ρ = 1900 kg/m 3 (see further down for the description of the ballast). The box of soil is surrounded on five sides by a layer of PMLs, that allow the wave to exit the box with little reflection, mimicking an infinite medium. The total thickness of the PML layer is 2.45 m, and the damping evolution in the PMLs is parameterized by a polynomial degree p = 3 and amplitude α = 9.87. The mesh, shown as insert in Fig.5, is discretized by 2.4 millions of non-structured elements, with 375 degrees of freedoms (DOFs) in each element, which amounts to a total of ≈ 160 millions DOFs in the entire model.

Figure 5 :

 5 Figure 5: Geometry, loading (black arrows), and sensors (lines of red dots) for the model. The insert displays the mesh. The ballast appears in yellow, the soil in blue, and the PML region in green.

Figure 6 :

 6 Figure 6: One realization of the density field in a ballasted railway track (PML not represented).

Figure 7 :

 7 Figure 7: Normalized displacement fields in the ballasted railway track at times: (a) t = 0.135 s; (b) t = 0.225 s; (c) t = 0.315 s; (d) t = 0.405 s; (e) t = 0.495 s; (f) t = 0.585 s.

Figure 8 :Figure 9 :

 89 Figure 8: Dispersion curve for a heterogeneous ballast resting on a homogeneous soil. Dispersion relations for the homogenized ballast shear wave (black solid line) and soil shear wave (black dashed line) are also indicated. The color scale is the amplitude of the average H pq in Eq. (5), averaged over 6 sensor lines and 5 random realizations of the ballast layer. The color scale is limited at -40dB.

  2 millions cube elements with edge 0.06 m. Polynomials of 4th order are used, inducing close to 1.4 billion DOFs in the entire model. The PML uses 3 layers of elements with power degree 2, amplitude 3.71 and total thickness of the PML equal to 1.69 m. A point force is positioned at [0.66, 0.66, 0.66] m, and polarized in the x direction. The time history of that source is that of Eq. (10) with a plane spectrum between f 2 = 40 Hz and f 3 = 140 Hz, and cutoff frequencies at f 1 = 35 Hz and f 4 = 145 Hz. Two different media will be considered in this section, and the corresponding two dispersion curves will be constructed. The first medium is homogeneous, with mechanical parameters corresponding to the (arithmetic) average of the ballast in the previous section, with κ = 179.60 MN/m 2 , µ = 35.32 MN/m 2 , and ρ = 1570 kg/m 3 . The second medium is the same heterogeneous ballast as in the previous section and in Section 3.1, with average values of κ = 179.60 MN/m 2 , µ = 35.32 MN/m 2 , and ρ = 1570 kg/m 3 , and standard deviations σ ρ = 3925 kg/m 3 , σ κ = 449.02 MN/m 2 and σ µ = 88.31 MN/m 2 . The correlation model is identical to the one used above to describe the ballast layer and the fields correspond to average velocities of V p = 380 m/s and V s = 150 m/s, and apparent velocities close to V p = 280 m/s and V s = 110 m/s. Since the geometry is this time isotropic, the sensor lines for the construction of the dispersion curves can be set in all directions before averaging. Fig. 10 shows the position of the sensors, along 91 lines separated by angles of 10 o in both the horizontal and vertical directions, and with sensors along the lines separated by distances of 0.066 m, from 0 m to 10.98 m from the source. The monitors record the amplitude of the displacement.

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: Position of the monitors (blue line) and the source (red point).

Figure 13 :

 13 Figure 13: Dispersion curves for infinite homogeneous (figure (a)) and heterogeneous media (figure (b)). The color scale is the amplitude of H pq evaluated in Eq. (5), limited to values above -40 dB. The dotted black line is the analytical solution for the homogenized pressure waves dispersion relation. The dashed black line is the analytical solution for the homogenized shear wave.

Figure 14 :Figure 15 :

 1415 Figure 14: Displacement magnitude in an infinite homogeneous medium over different frequency ranges: (a) 0 -50 Hz; (b) 50 -100 Hz; (c) 100 -150 Hz. The color scale is the amplitude of the displacement field, limited to values above -80 dB.

Acknowledgements

The spectral element software used for the simulations in this paper is developed jointly by CentraleSupélec, CEA Commissariat à l' Énergie Atomique and Institut de Physique du Globe de Paris. Within the SINAPS@ project, this development benefited from French state funding managed by the National Research Agency under program RNSR Future Investments bearing reference No. ANR-11-RSNR-0022-04. This work was performed using HPC resources from the "Mésocentre" computing center of Cen-traleSupélec and École Normale Supérieure Paris-Saclay, supported by CNRS and the French Ministère de l'Enseignement Supérieur, de la Recherche, et de l'Innovation (http://mesocentre.centralesupelec.fr/).


H N Vo

M F Pucci

S Corn

N Le Moigne

W Garat

S Drapier

P J Liotier
email: liotier@emse.fr

Capillary wicking in bio-based reinforcements undergoing swelling -Dual scale consideration of porous medium

Keywords: A, Biocomposite Natural fibre E, Liquid composite moulding Capillary wicking

Capillary wicking of liquids in natural fibrous reinforcements is a significant phenomenon in fibre-reinforced composites manufacturing through Liquid Composite Moulding (LCM) processes. Such phenomenon is however difficult to analyse due to heterogeneous, multiscale and variable fibrous medium during flow. Taking into account the fibre swelling, a modified Washburn's model is proposed to predict the capillary rise of liquids in flax fibres. Particularly, swelling and wicking are studied at two different scales, i.e. elementary fibres and individual yarns. Swelling effects have been considered in the model and a modified Washburn's theory has been developed. The proposed model satisfactory fits experimental results from capillary wicking tests, and hold for a 60% fibre volume fraction (V f ) of fabric on a relatively short duration while it fits well on a longer duration for lower V f .

Introduction

The use of natural fibres as reinforcements in composites has grown in the last decades. Indeed, natural fibres have interesting mechanical and damping properties, low density, and better sustainability [START_REF] Baley | Analysis of the flax fibres tensile behaviour and analysis of the tensile stiffness increase[END_REF][START_REF] Duigou | Environmental impact analysis of the production of flax fibres to be used as composite material reinforcement[END_REF][START_REF] Tran | Wetting analysis and surface characterisation of coir fibres used as reinforcement for composites[END_REF][START_REF] Fuentes | Wetting behaviour and surface properties of technical bamboo fibres[END_REF][START_REF] Placet | Diameter dependence of the apparent tensile modulus of hemp fibres: A morphological, structural or ultrastructural effect? Compos Part A[END_REF]. In some cases, e.g. flax fibres, the specific modulus is comparable or even higher than for glass fibres [START_REF] Wambua | Natural fibres: can they replace glass in fibre reinforced plastics?[END_REF]. The potential of flax fibres as reinforcements in structural composites can only be realized by taking advantage of twistless long fibre yarns. Besides, structural composite parts with high mechanical performances can only be obtained if high fibre volume fraction is reached. Therefore, Liquid Composite Moulding (LCM) processes are considered as promising and effective to manufacture structural composite parts reinforced with natural fibres.

In LCM processes, such as resin transfer moulding and resin infusion, the impregnation of fibrous reinforcements by resin is a crucial step contributing to the quality and performances of composite parts [START_REF] Park | Modeling and simulation of voids and saturation in liquid composite molding processes[END_REF][START_REF] Salvatori | 3D spacers enhance flow kinetics in resin transfer molding with woven fabrics[END_REF][START_REF] Caglar | Assessment of capillary phenomena in liquid composite molding[END_REF]. Those processes allow the manufacturing of large parts suitable for transportation industry applications [START_REF] Umer | The effect of yarn length and diameter on permeability and compaction response of flax fibre mats[END_REF][START_REF] Francucci | Processing of plant fiber composites by liquid molding techniques: An overview[END_REF]. The interest of composites reinforced by natural fibres in the transportation industry is due to their relatively high specific properties [START_REF] Charlet | Characteristics of Hermès flax fibres as a function of their location in the stem and properties of the derived unidirectional composites[END_REF]. Fibrous reinforcements are commonly considered as porous media, and they can be ideally approximated as a capillary tube network [START_REF] Schuchardtl | Liquid transport in composite cellulose -superabsorbent fiber networks[END_REF][START_REF] Masoodi | Darcy's law-based model for wicking in paper-like swelling porous media[END_REF]. Hence, wicking in those media can be described by the Washburn's theory [START_REF] Washburn | Note on a method of determining the distribution of pore sizes in a porous material[END_REF], where the medium porosity is considered as constant, e.g. for synthetic fabrics. In the case of natural fibres fabrics, fibres can be subjected to swelling during wicking, especially with polar liquids, leading to changes in the porous medium configuration [START_REF] Nguyen | Permeability of natural fiber reinforcement for liquid composite molding processes[END_REF][START_REF] Nguyen | Modeling of resin flow in natural fiber reinforcement for liquid composite molding processes[END_REF]. Natural fibres also have some limitations compared to synthetic fibres such as intrinsic hydrophilic character, low thermal stability, large variability in terms of morphology, chemical composition, and surface properties [START_REF] Charlet | Scattering of morphological and mechanical properties of flax fibres[END_REF]. The absorption of liquids during LCM processes results in an increase of fibre diameter [START_REF] Nguyen | Permeability of natural fiber reinforcement for liquid composite molding processes[END_REF]. Subsequently, the pore size in between fibres decreases which changes the flow characteristics of resin during the impregnation process. Therefore, the swelling of natural fibres should take into account in the Washburn's equation to better describe the wicking phenomenon in these specific porous media. The behaviour of fibres in terms of swelling will be different in resin then in water. Even if dynamic wetting of cellulosic substrates by liquid polymer has been studied recently [START_REF] Pucci | Temperature effect on dynamic wetting of cellulosic substrates by molten polymers for composite processing[END_REF], it has been chosen to use water as test liquid because it represents the worst possible case for swelling of natural fibres.

Swelling of natural fibres has been studied extensively using different measurement techniques such as optical microscopy, automated laser scanning, or X-ray micro-computed tomography [START_REF] Nguyen | Permeability of natural fiber reinforcement for liquid composite molding processes[END_REF][START_REF] Masoodi | A study on moisture absorption and swelling in bio-based jute-epoxy composites[END_REF][START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF][START_REF] Pucci | Wetting and swelling property modifications of elementary flax fibres and their effects on the Liquid Composite Molding process[END_REF][START_REF] Joffre | A Method to measure moisture induced swelling properties of a single wood cell[END_REF][START_REF] Garat | Swelling of natural fibre bundles under hygro-and hydrothermal conditions: Determination of hydric expansion coefficients by automated laser scanning[END_REF]. The effect of swelling on permeability of the fibrous medium has also been studied [START_REF] Nguyen | Permeability of natural fiber reinforcement for liquid composite molding processes[END_REF][START_REF] Francucci | Study of saturated and unsaturated permeability in natural fiber fabrics[END_REF]. Several natural fibres, used as reinforcements in composites, have been studied such as flax, hemp, sisal, jute, and wood. Swelling measurements were performed on elementary fibres, fibres bundles or yarns (case of flax) immersed in water (in optical

= ⎡ ⎣ ⎤ ⎦ h t cr γ cosθ η t ( ) ( ¯) 2 L a 2 (1)
where c is a parameter taking into account tortuosity, r ¯is a mean porous radius, θ a is an apparent advancing contact angle, γ L and η are the liquid surface tension and viscosity, respectively.

Eq. ( 1) can be rewritten as a function of the square mass gain over time:

= ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ m t cr ε πR ρ γ cosθ η t ( ) ( ¯) ( ) 2 L a 2 2 2 2 2 (2)
where ε, the porosity, is 1 -V f (V f : the fibre volume fraction), R is the radius of a cylindrical sample holder, and ρ is the liquid density. The squared-brackets term is related to both morphology of the porous medium and geometry of sample holder. Eq. ( 2) fits very well to the experimental data for synthetic fabrics because the medium porosity (related to geometric product cr ¯) remains constant over time [START_REF] Pucci | Capillary wicking in a fibrous reinforcement -Orthotropic issues to determine the capillary pressure components[END_REF].

Model development

For natural fibres, the porous medium changes over time during wicking of water due to fibre swelling. In a first approach, assuming that tortuosity is not affected by liquid and that the volume of lumen is neglected, hence the capillary wicking only takes place in the spaces between individual yarns and between elementary fibres. The dual scale capillary radius is then proposed as follows:

= -∅ +∅ cr c R r ¯((1 ) ¯¯')
yarn yarn [START_REF] Tran | Wetting analysis and surface characterisation of coir fibres used as reinforcement for composites[END_REF] where R ¯and r ¯' are respectively the mean capillary radii between individual yarns, and between elementary fibres (Fig. 1), and ∅ yarn is the surface ratio of yarns for a cross-section of the composite.

During wicking in natural porous media, the mean capillary radius was assumed to decrease linearly over time t due to swelling [START_REF] Schuchardtl | Liquid transport in composite cellulose -superabsorbent fiber networks[END_REF][START_REF] Masoodi | Darcy's law-based model for wicking in paper-like swelling porous media[END_REF][START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF][START_REF] Stuart | Modelling of wicking and moisture interactions of flax and viscose fibres[END_REF]. In line with this approach, the swelling effects on the evolution of the mean capillary radii R ¯and r ¯' should also be described. Eq. ( 3) can therefore be expressed as follow:

= -∅ - +∅ - cr t c R bt r at ( ¯)( ) ((1 )( ¯) ( ¯) ) yarn yarn 0 0 ' (4)
where R ¯0 and r ¯0' are the mean capillary radii before swelling between the individual yarns, and between the elementary fibres, respectively. Constants b and a express the swelling rates of a yarn and of an elementary fibre, respectively. Inserting Eq. (3) (before swelling: cr ( ¯)ini corresponding to R ¯0 and r ¯0' ) into Eq. (4), the following equation is obtained:

= - -∅ + ∅ cr t cr ct b a ( ¯)( ) ( ¯) ( (1 )
)

ini yarn yarn [START_REF] Placet | Diameter dependence of the apparent tensile modulus of hemp fibres: A morphological, structural or ultrastructural effect? Compos Part A[END_REF] which expresses for the evolution of the geometric product cr ¯during capillary wicking in natural fibres taking into account the fibre swelling at a double scale. In this equation, parameters such as ∅ yarn , b, a, R ¯0, and r ¯0' can be determined experimentally.

It must be noted that modifications of porosity ∊ t ( ) related to the variation of the geometric product cr t ( ¯)( ) are taken into account in the microscopy, automated laser scanning device) or placed in high relative humidity environment (in tomography machine). Optical microscopy and automated laser scanning directly give apparent diameter measurements while X-ray micro-computed tomography require image post analysis, but allow 3D reconstruction of fibre morphology. Particularly, with optical microscopy, fibre diameter was measured manually by picking two points at the edge of a single elementary fibre, and monitored over time during swelling. As a result, the fibre diameter as a function of time was obtained and the swelling rate could be estimated [START_REF] Nguyen | Permeability of natural fiber reinforcement for liquid composite molding processes[END_REF][START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF]. Medium porosity is defined as spaces between fibres, that is to say all spaces that are not occupied by fibres. It can be calculated as the complement of a ratio of the volume of fibres on the overall volume (ε, the porosity, is 1 -V f ; V f being the fibre volume fraction). Depending on scales, such spaces can be considered as in between yarns, fibre bundles, and elementary fibres [START_REF] Park | Modeling and simulation of voids and saturation in liquid composite molding processes[END_REF][START_REF] Nguyen | Permeability of natural fiber reinforcement for liquid composite molding processes[END_REF][START_REF] Nguyen | Modeling of resin flow in natural fiber reinforcement for liquid composite molding processes[END_REF][START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF][START_REF] Lee | Analysis and minimization of void formation during resin transfer molding process[END_REF]. As pointed out above, the conventional Washburn's equation could not describe wicking for natural fibres absorbing water [START_REF] Schuchardtl | Liquid transport in composite cellulose -superabsorbent fiber networks[END_REF][START_REF] Masoodi | Darcy's law-based model for wicking in paper-like swelling porous media[END_REF][START_REF] Stuart | Modelling of wicking and moisture interactions of flax and viscose fibres[END_REF], because the porous medium changes over time due to fibre swelling. More precisely, it was shown that the geometric product (cr ¯) in Washburn's equation related to tor- tuosity c and mean porous radius r ¯ changed during wicking time. Therefore, the characteristic product (cr ¯) of the porous medium should take into account swelling effects and the Washburn's equation is modified accordingly.

In a previous work [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF], a model accounting for the swelling of elementary fibres (related to the geometrical product) was proposed. However, both swelling of yarns and elementary fibres, and then spaces in between elementary fibres and individual yarns, should be considered to describe the porous structure of fibre network and its evolution during wicking [START_REF] Nguyen | Modeling of resin flow in natural fiber reinforcement for liquid composite molding processes[END_REF]. Therefore, in the present work, a modified Washburn's model is proposed, taking into account the porous media and the swelling effects. Particularly, the swelling effects on wicking are evaluated at two scales, i.e. the scale of elementary fibres and the scale of individual yarns, based on optical microscopy and X-ray tomography measurements. Furthermore, the evolution of the spaces in between elementary fibres and individual yarns are also accounted in the model. Experimental wicking data were obtained with a tensiometer in the warp direction of quasi-unidirectional carbon and flax fabrics with 60% V f . It should be noticed that V f of 60% was chosen in this study because structural composites with high mechanical performance are commonly processed at high V f . In the experiment of wicking, such V f corresponds to the ratio of the volume of fibres to the internal volume of the cylindrical sample holder. Then, the swelling behaviour of flax fibres in water was characterized at both scales, i.e. elementary fibres and yarns, as input experimental data to establish a modified model of wicking adapted to natural fabrics and evaluate its validity at high fibre volume fraction.

Theory

Conventional Washburn's equation for fabrics

A liquid rise in porous media is usually described by Washburn's equation [START_REF] Washburn | Note on a method of determining the distribution of pore sizes in a porous material[END_REF]. This equation models the complex porous structure of fibre network as bundles of capillary tubes in a cylindrical holder [START_REF] Schuchardtl | Liquid transport in composite cellulose -superabsorbent fiber networks[END_REF]. The model predicts the flow front position h of liquid as a function of time:

> if cr t cr ( ¯)( ) ( ¯)fin = - -∅ + ∅ cr t cr ct b a ( ¯)( ) ( ¯) ( (1 )
)

ini yarn yarn [START_REF] Wambua | Natural fibres: can they replace glass in fibre reinforced plastics?[END_REF] else

= cr t cr ( ¯)( ) ( ¯)fin (7)

Washburn's relation with swelling

Feeding the proposed modification of porosity, and of pore radius diminution (Section 2.2 and Section 2.3), together with using initial wicking parameter ( cr ( ¯)ini is determined experimentally) into the Washburn's equation (Eq. ( 2)), it is noted that the apparent advancing contact angle θ a and the parameter accounting for tortuosity c in Eq. ( 2) still need to be identified. These two parameters are obtained when fitting the model to the experimental data of wicking (presented in Section 4.4).

The proposed model implicitly considers an important assumption. Namely, the entire volume of the fibre sample begins to swell when the wicking starts. In other words, no differential or incremental swelling is specified in each position along the elementary fibres and yarns when liquid starts to be in contact with the sample. This assumption is supposed to be valid due to the similar time scales of the fibre swelling and the wicking in the present work. For more complex bio-based preforms, deeper analysis requiring a numerical model will be needed. It will be the main focus of further studies.

Materials and methods

Experiments were designed to investigate capillary effects coupled with swelling during wicking of water in flax fabrics. Wicking test was first performed at the macro-scale of synthetic and natural fibre fabrics, i.e. carbon fabrics and flax fabrics in the warp direction. Carbon fabrics were used in this test to validate the method at high V f . Swelling tests were then carried out at the micro-scale on yarns and elementary fibres for natural fibres only.

Materials

The reinforcements are quasi-unidirectional fabrics of carbon and flax. Carbon fabrics were provided by Hexcel with commercial designation 48580®. The fabrics contain 3% of glass yarns in volume, and have areal weight of 540 g/m 2 and density of 1.8 g/cm 3 . Flax fabrics were supplied by Libeco with material code FLAXDRY UD 180®. Their areal weight and density are 180 g/m 2 and 1.45 g/cm 3 , respectively.

The cylindrical sample holder in which fibres were placed is made of poly(methyl methacrylate) polymer, it has an inner radius of 6 mm and a height of 20 mm. The tube was cleaned with sulfochromic acid (50% of a saturated solution of potassium dichromate and 50% of concentrated sulfuric acid), rinsed with water and dried at ambient condition before testing.

n-Hexane and distilled water were selected as reference liquids for this study. n-Hexane is a totally wetting liquid, meaning that the advancing contact angle θ a is supposed to be zero. This chemical is sui- table for determination of initial non-swollen geometric parameter of the fibrous medium cr ( ¯)ini . Distilled water was used to determine the apparent advancing contact angle θ a from wicking tests, and was also used for swelling tests on elementary flax fibres and yarns. Some characteristics of n-Hexane and distilled water such as surface tension, density, and viscosity are presented in Table 1.

Methods

Wicking tests

The details of procedure were well described in previous works [START_REF] Pucci | Capillary wicking in a fibrous reinforcement -Orthotropic issues to determine the capillary pressure components[END_REF][START_REF] Pucci | Capillary effects on flax fibers -Modification and characterization of the wetting dynamics[END_REF]. In this study, the amount of carbon and flax fabrics were determined so as to get a global fibre volume fraction V f of 60% which is higher than in previous studies [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF][START_REF] Pucci | Capillary wicking in a fibrous reinforcement -Orthotropic issues to determine the capillary pressure components[END_REF]. These fabrics were then rolled along the transverse direction, and inserted into the cylindrical tube. Capillary wicking was performed with a tensiometer (DCAT11, Data-Physics Instrument GmbH, Filderstadt, Germany). The tube was clamped to the tensiometer and the vessel with test liquid moved with a speed of 0.5 mm/s up to the detection of the holder. The holder is linked to an electronic microbalance with a resolution of 10 µg. The data were then recorded with the square of mass gain as a function of time. After measurement with n-Hexane, each sample was dried at 60 °C for at least 2 h. Their mass loss related to evaporation of n-Hexane was checked with microbalance until the weight of sample did not change before measuring with water. All the wicking tests were performed at 20 °C. For each fabric, at least nine samples were tested.

To determine the initial geometric product cr ( ¯)ini , the procedure is as follows. In constant porous medium, the evolution of square mass gain as a function of time was observed linearly with both n-Hexane and water as described in conventional Washburn's equation. An example of linear fits for carbon fabrics is presented in Fig. 2. The geometric product cr ¯is determined in Eq. ( 2) using n-Hexane, for which the apparent advancing contact angle θ a is supposed to be 0°. Because the porous medium remains constants during wicking, the cr ¯obtained is considered as cr ( ¯)ini . From linear fits of wicking with water and knowing cr ( ¯)ini , the apparent advancing contact angle with water is then de- termined.

Mean capillary radius between yarns by X-ray tomography

In order to evaluate the fabric internal structure in 3D within the entire tube area, nanofocus X-ray computed tomography (nano-CT) was used. Measurements were performed with a Phoenix NanoTom S (GE Measurements and Control solutions, Wunstorf, Gernany). The device was equipped with a tungsten target and the applied voltage and current were 100 kV and 100 µA, respectively. After scanning, 3D images were reconstructed from the radiographs using "Phoenix datos | x 2.5" reconstruction software.

The images obtained by tomography were then analysed by the image analysis software ImageJ to determine the surface ratio of yarns ∅ yarn as well as the mean capillary radius between the yarns before With ρ, η, γ L , γ L p , and γ L d representing respectively the density, viscosity, surface tension and its polar and dispersive components, of liquid which are found in literature [START_REF] Garat | Swelling of natural fibre bundles under hygro-and hydrothermal conditions: Determination of hydric expansion coefficients by automated laser scanning[END_REF].

proposed model the same way as in the previous model [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF].

Pore radius diminution

It should be noted that fibres and yarns swell until their maximal dimensions in a definite time. Assuming tortuosity as a constant, the geometric product cr ¯ thus varies over time until its final lower value (cr ¯)fin is reached [21]. Therefore, (cr ¯)fin can be defined from minimal mean capillary radius. In the previous model [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF], this parameter was limited to a decrease of 30% because fibres swelled of around this value and also the model could not hold further diminution of mean capillary radius. The present model allows diminution down to 50% of initial mean capillary radius value, which is the chosen minimal value to extend the time range of validity of the model.

Eq. ( 7) gives the boundary condition for cr ¯, taking into account the maximal swelling of fibres: swelling R ¯0. This software was used to determine the pore volume percentage. Flax fabrics were studied in dry condition and saturated with water (Fig. 3).

Swelling rate and mean capillary radius between elementary fibres by optical microscopy

In order to determine the swelling rate of flax elementary fibres and individual yarns, an optical microscope (Laborlux 11 POLS, Leitz, Germany) equipped with a "mon-CDD Leica DFC420" digital camera (resolution 1728 × 1296 pixels) and an image acquisition software Archimed® (Microvision Instruments, France) has been used. Elementary fibres (extracted from yarns in dry condition) and yarns were positioned in a container made of two glass plates separated with double-sided tape. Distilled water was introduced by capillary forces between the two plates until the water filled the container and fibre swelling phenomenon occurred. The space between the two plates was sufficient for the fibres to freely move in water. No tension was applied to the elementary fibres and yarns. Four to six images per second were recorded from the time of water introduction, until fibres reached a stabilized swollen state. The ImageJ software was used to measure fibre diameter changing over time (assuming a circular section). Because diameter of elementary fibres and yarns vary along their length and greatly change from one to another, diameters were measured at least at five different positions (Figs. 4 and5). It should be pointed out that in the case of flax yarns, diameter measurements were made more difficult due to fibrillation as seen in Fig. 5. At least five samples of elementary fibres and yarns were analyzed for the swelling test.

Optical microscope was also chosen to determine the mean capillary radius between elementary fibres before swelling r ¯0' . It would be very complicated to measure the r ¯0' with tomography since the minimal voxel size (around 800 μm 3 ) obtained in a single yarn is much higher than the diameter of elementary fibres and spaces in between them. Therefore, images taken from an epoxy based composite reinforced by flax quasi-UD fabrics (Fig. 1) were analyzed with ImageJ in order to estimate r ¯0' . Those samples were manufactured in dry conditions under vacuum, assuming that the resin had no significant swelling effect on fibres. It is also assumed that the r ¯0' does not change much with the V f .

Results and discussion

Validation from wicking in carbon fabrics

To validate that this method can be applied for reinforcements with V f of 60%, the wicking tests were carried out with n-Hexane on carbon fabrics, following the procedure detailed in Section 3.2.1. The evolution of square of mass gain as a function of time was monitored, as presented in Fig. 6a. The linearity confirms the validity of Washburn's equation with the totally wetting liquid for carbon fabrics. As expected, the equilibrium mass gain in this study at 60% V f is lower than the one in a previous study at 40% V f [START_REF] Pucci | Capillary wicking in a fibrous reinforcement -Orthotropic issues to determine the capillary pressure components[END_REF]. This is logical since the mass of liquid inside the sample holder decreases with the void space available. Table 2 presents the results of linear fitting with the correlation coefficient (R 2 ), and subsequently the determined initial geometric product cr ( ¯)ini for each test. The mean result of geometric product cr ( ¯)ini was then calculated as 2.34 ± 0.35 µm. This mean value is lower than the one calculated at 40% V f (12.1 ± 1.5 µm) [START_REF] Pucci | Capillary wicking in a fibrous reinforcement -Orthotropic issues to determine the capillary pressure components[END_REF]. This is also expected since the r ¯decreases as V f increases [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF].

Similar tests were also performed with distilled water to determine an apparent advancing contact angle θ a . The results of capillary wicking with water are presented in Fig. 6b. The linear slope is also observed, and hence the Washburn's equation can be applied. The advancing contact angle is calculated based on both the slope of the linear fit and the geometric product determined with n-Hexane. The result of each test is presented in Table 2, and the mean calculated value of the advancing contact angle was 55.4°± 3.8°. The linearity of wicking and the resulting parameters validate that this method can be used for V f up to 60%.

Wicking in flax fabrics

Using the same method as for carbon fabrics, the initial geometric product cr ( ¯)ini was determined with capillary wicking of flax fibres with n-Hexane. Fig. 7a shows the results of wicking curves for five different samples. The linearity was also observed, supporting that Washburn's equation can be used. The results of wicking tests for flax fabrics are presented in Table 3 and a mean cr ( ¯)ini of 3.77 ± 0.37 µm was ob- tained.

Regarding capillary wicking with water, non-linear trends were observed with a large variation between the samples (Fig. 7b). Consequently, Washburn's equation cannot be applied. As pointed out in the introduction, the chemical composition and surface properties of natural fibres and the variations of diameter along their length lead to heterogeneous wetting as well as complex tortuosity of pore between the fibres (in both individual yarns and elementary fibres), influencing the capillary wicking behaviour [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF][START_REF] Pucci | Capillary wicking in a fibrous reinforcement -Orthotropic issues to determine the capillary pressure components[END_REF][START_REF] Pejic | The effects of hemicelluloses and lignin removal on water uptake behavior of hemp fibers[END_REF]. Furthermore, the swelling of flax fibres during wicking causes an increase of fibre diameter within the cylindrical tube, and hence a decrease of porosity, and of the mean capillary radius r ¯. As a result, capillary rise speed evolves during the wicking test (i.e. it is decreased), and the equilibrium weight is more slowly reached. It should be pointed out that the equilibrium weight is also determined by water sorption within the flax fibre cell walls which occurs concomitantly with fibre swelling.

Measurement of capillary radii, surface ratio of yarn, and swelling

Before swelling, R ¯0 and r ¯0' are considered as the largest mean ca- pillary radii between individual yarns and between elementary fibres, respectively. ∅ yarn is defined as the surface ratio of yarns at the scale of composites, taking into account the largest mean radius between individual yarns R ¯0 and the largest mean capillary radius between elementary fibres r ¯0' . Their values (presented in Table 4) were determined by images taken by X-ray tomography (for R ¯0 and ∅ yarn ) and optical microscopy (for r ¯0' ) and then image analysed with ImageJ as described in Section 3.2.2 and Section 3.2.3.

As mentioned earlier, the diameter of an elementary flax fibre as well as the one of a yarn are not uniform. Table 5 shows results of diameter distribution before wetting. This is in line with several studies which show that flax fibres exhibit a large scatter of diameter along their length [START_REF] Nguyen | Permeability of natural fiber reinforcement for liquid composite molding processes[END_REF][START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF][START_REF] Baley | Average tensile properties of French elementary flax fibers[END_REF][START_REF] Garat | Analysis of the morphometric variations in natural fibres by automated laser scanning: Towards an efficient and reliable assessment of the cross-sectional area[END_REF]. After wetting in water, elementary flax fibres and yarns have significantly swelled as clearly seen in Figs. 4 and5. Fig. 8 presents the results of swelling tests with water for elementary fibres and yarns. It is noticed that faster swelling rates were related to larger initial diameters for both elementary fibres and yarns. The mean values of their swelling rate are given in Table 4, respectively coefficients b and a in Eq. (4).

Comparison between experiment and model

Capillary wicking is well described by conventional Washburn's equation for invariant porous structure morphology. However, wicking curves presented in Fig. 7b evidence that changes in the porous medium during wicking in case of natural fibres require modification of the Washburn's equation. Therefore, the modified Washburn's model proposed in Section 2.2 was applied to fit experimental data.

As mentioned in Section 2.4, the apparent advancing contact angle θ a and the parameter accounting for tortuosity c in Eq. ( 2) require to be identified to fit the experimental data. As observed in Fig. 9a, based on the modified Washburn's equation, the values of θ a and c could be , b and a representing respectively the mean capillary radius before swelling between yarns, the mean capillary radius before swelling between elementary fibres, the surface ratio of yarns at the scale of composites, the swelling rate of a yarn, and the swelling rate of an elementary fibre.

* Only one sample was scanned by nano-CT due to complexity of the set-up, hence no standard deviation is given. It must be noted that ∅ yarn (as defined in Section 2.2) only considers a surface ratio of yarns on an overall section of sample, and it thus does not consider spaces in between elementary fibres. Therefore, the value of ∅ yarn obtained does not correspond to 60% V f and is thus necessary higher.

Table 5

Diameter distribution of an elementary flax fibre (Fig. 4a) and of a flax yarn (Fig. 5a) before wetting. determined, and the model satisfactory fits the experimental data with values of θ a and c of 67.5°and 0.166, respectively. Another originality of the proposed model compared to the previous one is that it allows to determine the value of the factor c related to tortuosity that becomes an independent and adjustable parameter from the geometric product cr ¯. This is the first time that c and r ¯are considered separately.

Moreover, the modified Washburn's model holds well at the first three seconds of wicking. The wicking of water in natural fibre fabrics at 60% V f is more complex than at lower V f . This is in agreement with the previous work [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF] which showed that the effect of swelling is more significant at higher fibre volume fraction. This is also why the previous model [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF] cannot describe phenomenon occurring at V f higher than 40%. In addition, at 40%, the previous model could represent experimental data for four seconds whereas the proposed model holds for eight seconds as seen in Fig. 9b (twice longer). The present model can also satisfactory predict experimental data at V f of 30% and 35% [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF], which highlights its versatility. Even if the previous model was predicting the wicking at low V f (30-35%), it was not able to cover the full V f range (30-60%). The present one has shown its reliability from 30% to 60% of V f . Therefore, the proposed modified Washburn's model that considers a double scale of swelling (elementary fibres and yarns) and setting measured swelling rates, expands the range of validity of the previous model [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF] in terms of V f , and also improves it for lower V f .

Conclusions

The understanding and prediction of capillary wicking in natural fibre fabrics is determinant for improving manufacturing with LCM processes and the quality of composite parts. An experimental methodology was developed to evaluate capillary rise, taking into account the multi-scale swelling effects in natural fibres. In the present work, swelling was studied at two scales: elementary fibres and individual yarns. Wicking tests were performed on quasi-unidirectional carbon and flax fabrics at 60% V f . Carbon fabrics showed a linear trend following the Washburn's equation, while a non-linear trend was observed with flax fabrics for wicking tests in water. This highlights the need to take into account swelling effects in natural fibres to better describe capillary wicking. Based on a modified Washburn's equation that considers a double scale of swelling and experimental swelling rates, the capillary parameters have been determined experimentally for both reinforcements. The proposed model has shown to fit very accurately experimental data during the first three seconds of wicking, at V f up to 60%. At lower V f , i.e. 40%, the proposed model can fit experimental data for eight seconds which shows the improvement to twice longer duration compared to the standard model [START_REF] Pucci | Capillary wicking in flax fabrics -Effects of swelling in water[END_REF]. A numerical implementation of this model would allow taking into account the differential swelling of fibres depending on the duration that they were immersed in liquids. Particularly, elementary fibres and yarns at the bottom of the sample start swelling before the ones at the middle. This effect, which would be addressed in further works, is not included in the proposed model and would extend drastically its domain of validity.

Fig. 1 .

 1 Fig. 1. Illustration of the two-scale pore radii in the cross-section of an epoxy based composite reinforced by quasi-unidirectional flax fabrics. Broken lines are added to visualize individual yarns.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Linear fit of a wicking test for carbon fabrics at 60% V f with (a) n-Hexane and (b) distilled water. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. Optical measurements of diameters for an elementary flax fibre before (a) and after (b) swelling in distilled water. The white lines marked on the fibres represent the positions where fibre diameters were measured. Caution was paid to measure the diameters at the same position along the fibre axis before and after swelling.

Fig. 5 .

 5 Fig. 5. Optical measurements of diameters for a flax yarn before (a) and after (b) swelling in distilled water. The white lines marked on the yarns represent the positions where yarns diameters were measured. Caution was paid to measure the diameters at the same position along the yarn axis before and after swelling.

Fig. 6 .

 6 Fig. 6. Wicking curves in warp direction of carbon fabrics with (a) n-Hexane and (b) distilled water at 60% V f . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .

 7 Fig. 7. Wicking curves in warp direction of flax fabrics with (a) n-Hexane and (b) distilled water at 60% V f . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .

 8 Fig. 8. Swelling rate of (a) elementary fibre and (b) yarn versus initial diameter. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

  CRediT authorship contribution statement H.N. Vo: Writing -original draft, Methodology, Investigation, Formal analysis. M.F. Pucci: Conceptualization, Validation, Investigation, Writing -review & editing. S. Corn: Conceptualization, Methodology. N. Le Moigne: Conceptualization, Methodology. W. Garat: Methodology. S. Drapier: Writing -review & editing. P.J. Liotier: Conceptualization, Project administration, Supervision, Formal analysis, Writing -review & editing.

Fig. 9 .

 9 Fig. 9. Fitting of experimental with the modified Washburn's equation for flax fabrics with water at (a) 60% V f and (b) 30, 35, 40% V f [17].

Table 1

 1 Characteristics of test liquids at 20 °C.

		ρ (g/cm 3 )	η (mPa.s)	γ L (mN/m)	γ L p (mN/m)	γ L d (mN/m)
	n-Hexane	0.659	0.32	18.4	0.0	18.4
	Water	0.998	1.00	72.8	51.0	21.8

Table 2

 2 Results of wicking tests for carbon fabrics with n-Hexane and water at 60% V f .

	Carbon	n-Hexane			Water
		Linear slope	R 2	cr ( ¯)ini (µm)	Linear slope	R 2	θ a (°)
	Sample 1	0.057	0.9973	2.24	0.097	0.9996	54.2
	Sample 2	0.052	0.9974	2.05	0.093	0.9968	52.0
	Sample 3	0.075	0.9988	2.95	0.103	0.9997	61.8
	Sample 4	0.059	0.9985	2.30	0.097	0.9984	55.2
	Sample 5	0.055	0.9977	2.15	0.095	0.9987	53.5

Table 3

 3 Results of wicking tests for flax fabrics with n-Hexane at 60% V f .

	Flax	Linear slope	R 2	cr ( ¯)ini (µm)
	Sample 1	0.106	0.9963	4.14
	Sample 2	0.093	0.9953	3.64
	Sample 3	0.104	0.9957	4.09
	Sample 4	0.095	0.9953	3.73
	Sample 5	0.082	0.9972	3.23

Table 4

 4 Results of measurement related to swelling parameters of flax fibres.

	R ¯0 (µm)	r ¯0' (µm)	∅ yarn	b (µm/s)	a (µm/s)
	108.84 ± 27.31	6.63 ± 2.55	0.66*	4.76 ± 2.66	2.11 ± 0.66
	With	∅ ¯, ¯, yarn R r 0 0 '

Acknowledgments

The authors would like to thank O. Valfort for performing the computed tomography machine of reinforcement materials.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.