Stochastic heterogeneous material modeling for wave propagation in a ballast layer

This paper discusses the dynamical behavior of a randomly-fluctuating heterogeneous continuum model of the ballast. The Young’s modulus is modeled as a random field parameterized by its average, its variance and a correlation model. A numerical model of the ballast and the surrounding soil is then constructed based on an explicit spectral element solver. This model allows to describe numerically the wave field generated, as well as to construct dispersion curves for the ballast-soil model. The influence of heterogeneity is discussed by comparison with a similar model where the ballast is assumed homogeneous.


Introduction
The competition with other means of transportation has increased the demand for performance in the railway industry.The dynamical loads caused by the passage of high-speed trains accelerate track deterioration and damage neighbor buildings [1].In order to study this problem two classes of numerical models are used to try to predict this behavior: (1) discrete approaches, in which each grain of the ballast is represented by a rigid body and interacts with its neighbors through nonlinear contact forces (using molecular dynamics or non-smooth contact dynamics); and (2) continuum approaches, in which the ballast is replaced by a homogenized continuum and the classical Finite Element Method (FEM, or similar) can be used to solve the problem.The discrete approaches are today capable of solving only a few meters-length of ballast.The absence of an efficient scheme for parallelization over large clusters of computers [2] and the coupling between the grains and the soil [3] remains open problems.On the other hand, the methods based on continuum mechanics (like FEM) can solve very large models thanks to efficient parallelization schemes.However, a homogeneous continuum cannot represent the complex fields of stress and strain present in the ballast layer.
The dynamical behavior of ballast can be completely described by its dispersion curve.The numerical simulation of dispersion curves in railway tracks is usually performed on homogeneous media [4][5][6].The most ad-vanced models consider 3D finite elements and model reduction [7].Unfortunately, these homogeneous models cannot reproduce accurately the experimental observations [8].The discrete nature of the material and the complex interaction between the grains produce a "non trivial" mechanical behavior.It is necessary to consider the granular character of the material.The simulation of wave propagation in granular media for simple geometries has been treated using discrete models, and has succeeded in reproducing a wide range of experiments [9,10].The influence of pre-stress on the wave velocity in particular can be simulated [11] as well as the influence of the pressure in the grain network [12].The simple addition of random "heterogeneity" in a mass-spring network allows to reproduce the complex frequency response seen in measurements of sound propagation in a granular system [13].Note that [9] observed the coexistence of an attenuated coherent signal and the "coda", which is a classical feature of seismic signals [14].
The objective of this work is to use a statistical approach to model the mechanical properties of the ballast layer.This statistical approach aims to reproduce in a continuum model the heterogeneity of the stress field in a discrete granular model.The mechanical properties (in practice the Young's modulus) are modeled as stochastic random field.This approach was introduced in [15] and we apply it here to the construction of dispersion curves of ballasted railway tracks.

Numerical Framework
The key point of this model is the description of the mechanical properties.To build this stochastic continuum heterogeneous medium we need to define two things: (1)first-order marginal law; and (2) a correlation model.The first-order marginal law for the Young Modulus was identified using a optimization process when one discrete model was compared with one stochastic continuous one.The correlation model chosen here is based in the literature.
The second point is the wave propagation in this stochastic heterogeneous model.The construction of the dispersion diagrams was made using 3 different models, they were described in the Sec.2.2.The methodology adopted was based on a 2D-FFT.We applied the procedure explained by Alleyne [16] to recover the dispersion curves with time-displacement signals.

Discrete Model
The model consists of 29 cylinders with a radius of 35 cm and a height of 39 cm.Samples were taken from this cylinders as 48X48X32 cm 3 orthogonal parallelepiped.The total number of particles in each cylinder is around 2700.The particles are convex polyhedra with diameters between 2.5 cm and 5 cm, whose shapes have been obtained by digitalization from real ballast, one example was presented in the Fig. 1.An isotropic confinement pressure of 60 kPa is applied on all samples, a vertical load on the upper face of 63 kN, and the gravity is considered within all samples.The model was calculated using non-smooth contact dynamics code [17].

Continuum Model
The continuum model has one cubic sample of 48X48X32 cm 3 , taken from cylinders with radius of 35 cm and 39 cm high.A radial pressure of 60 kPa is applied, a pressure of 163 kPa was applied on the upper face.Gravity is considered, the density was assumed constant, 1500 kg/m 3 .In terms of first-order marginal law, we will model Young's modulus field with a log-normal distribution.This model requires the definition of a mean value Ē, a variance σ 2 , and a correlation model.We took Ē = 80 MPa given by the INNOTRACK Project report [18], variance (estimate below) and correlation length.A theoretical correlation model was used here to estimate the correlation function, based in the Percus-Yevick model [19] , the characteristic diameter was took as 3.7 cm.The value of the Poisson ratio was taken from the report [18], as ν = 0.23.

Identification Process
This task requires the comparisons between the discrete and the continuum media.We choose the stress (σ zz ) as comparison variable.This choice try to represent the same statistical behavior between the equivalent stress in discrete medium and the continuum one.
It is interesting to define a generalization of the notion of stress tensor for these granular media.Several proposals have been made by [20][21][22] for the static regime.Given a network of contact forces f c α (at contact points c and with coordinates α) in a granular medium, the equivalent stress field were used by [23,24], and can be defined as: where c is the vector linking the centers of the two particles in contact at c, the sum is on the N c contact points in the averaging volume V, see [25] for more details.Obviously, this quantity depends on the size of the averaging volume (with more variability over the whole sample for smaller volumes).One realization of the random Young Modulus was inserted in a classical FEM model.The stress distribution found in this model is compared using the L2 norm with the discrete one.At this moment the optimization process take place.We used a Nelder-Mead simplex algorithm [26].This process stopped when the tolerance was reached.The process found the minimum distance is at σ 2 E / Ē2 ≈ 15.6.For this variance a comparison of the stress distribution in log-log scale was presented in the Fig. 2, where a good approximation for the values around the mean was noted, the exponential behavior was well fitted, however the match in the tails was not perfect.
A small remark must be made here: This approach can be applied for others discrete models and contact models.Models based on molecular dynamics [27], for example, can add the influence of the grain stiffness.Other contact mechanics, like the inclusion of the rolling contact [28,29] adding rolling resistance, they will change the equivalent stress distribution.This modifications would be captures by the statistical properties of the continuum model.

Geometry and boundary conditions
Three models were built, one of those models can be seen in the Fig. 3, it was based in a 250000 elements, with 375 degree of freedoms in each element, given ≈ 94x10 6 degrees of freedom.This model is 32 meters long.One prescribed point force was impose as excitation.A rickers  pulse with amplitude 1 and central frequency 210 Hz, for the other models the central frequency were modified.The Fig. 3 show also the region where vertical load is applied in the top of the ballast layer.The other two models had 320000 and 380000 elements with 120x10 6 and 142x10 6 degrees of freedom.They only difference is the length of the railway, one is 40 and the other 48 meters long.The excitation for those models were 100 and 30Hz.Due the large number of degree of freedoms present in the model we choose a parallelizable and scalable code.This code, called SEM3D, it solve the dynamic equation in a time domain, using a explicit time integration [30,31].
The Fig. 3 show these position of the captures used to record the displacement field in the model.Three rows of captures were put inside of the ballast co-linear to the longitudinal axis, placed in the top, half, and the bottom of the ballast layer.
The mechanical properties used in the soil are the same of the ballast the follow: V s = 150m/s, V p = 380m/s, and ρ = 1500kg/m 3 .Two kind of simulations were performed in each model, homogeneous and heterogeneous for the ballast material.The statistical mean of the mechanical properties was keep the same in this two models.

Influence of the heterogeneity
We start with the homogeneous model of the ballast.The dispersion curves are plotted in Fig. 4, here we present just the result of the top line, the other had similar results.As expected, we retrieve approximately the dispersion relation for a Rayleigh wave in a homogeneous medium, with a velocity close to that of a shear wave (indicated by the black dotted line on Fig. 4).Although most of the displacement field corresponds to shear along the line that we are measuring, a small amount of P wave is also encountered at extremely low frequency, probably due a coupling between P and S waves in the low frequency.
The dispersion relations corresponding to the heterogeneous ballast are plotted in Fig. ??.Two important effects can be observed: (i) the propagating behavior is dispersive in the low frequency range, with a very strong slowing of the wave velocity with increasing frequency, similar to that observed in the experimental work of [32] on unconsolidated granular packings, and (ii) the dispersion relation vanishes at higher frequencies.Note that it does not mean that there is no energy in the system, but only that the energy does not propagate.We also observed some energy along the P-wave dispersion curve, which is probably due to a strong coupling of P and S wave on the heterogeneities of the ballast.

Conclusions
In this paper, we constructed large scale randomlyfluctuating continuum models of ballasted railway tracks.These models allowed to construct dispersion curves for the ballasted layers.The importance of the heterogeneity was stressed very strongly.In particular, it was shown to create relatively slow waves at medium frequencies and a filter effect in the higher frequencies.
Because of the geometry of the ballast, and its heterogeneity, it is absolutely impossible to solve the dynamical equations analytically.Computer simulations are therefore an irreplaceable and unavoidable tool to understand the dynamical behavior of the soil-ballast system.

Figure 1 .
Figure 1.Example of one of the 29 samples of granular medium simulated with the NSCD.

Figure 2 .
Figure 2. Relative distance between the probability distribution function σ zz with V = 8 3 cm 3 .Comparison of the probability distribution function with normalized variance σ 2 E / Ē2 ≈ 15.6.The black line with circles is the equivalent stress distribution in the granular media and the blue line with cross is the stress distribution in the FEM model.

Figure 3 .
Figure 3. Cut view of the geometry and position of the captures lines.The red lines corresponds to the position of the captures lines.Three different sets where plotted at the top of the ballast, in the half of the layer, and on the interface between soil and ballast.The force is represented by the black solid arrow.

Figure 4 .
Figure 4. Dispersion curves for the homogeneous ballast resting on soil, upper figure, and dispersion curves for the heterogeneous ballast, lower figure.Dispersion relations for the shear wave (black dotted line) and pressure wave (black solid line) in the ballast are also represented.