On the selection of the most plausible non-linear axial stress–strain model for railway ballast under different impulse magnitudes

Abstract

The effective maintenance and health monitoring of ballasted railway tracks, which involves the determination of differential settlement, track support stress and stiffness, and the strain-hardening property of ballast, is essential. The vertical stress–strain behavior of the ballast layer is primarily responsible for the irrecoverable strains and settlements in tracks, leading to further track degradation. This article reports the development of a series of applicable yet simple uniaxial models and the selection of the most plausible one for capturing the behavior of vertical stresses and strains in ballasts utilizing a set of measured vibration data of the rail–sleeper–ballast system from a Bayesian perspective. From the literature, the dynamic behavior of ballast can be divided into linear and non-linear regions. Under small amplitude vibration, the stress–strain property is linear and elastic, while the behavior becomes non-linear and inelastic once the elasticity limit is exceeded. By integrating the linear phase to some well-known non-linear engineering material laws, a list of new ballast stress–strain model classes was developed. An enhanced Markov chain Monte Carlo–based Bayesian scheme was utilized to explicitly handle the uncertainties in the model updating process, while the Bayesian model class selection method was employed to select the most plausible ballast stress–strain model class under the prevailing system conditions. The proposed methodology was verified using three sets of measured acceleration data from impact hammer tests on an in situ sleeper with simulated ballast damage. The obtained results suggest that the linear-elastic model is sufficient for small amplitude vibrations, while the modified Voce model is the most plausible amongst the investigated model classes for high impact load. The results also demonstrate the importance of the non-linear ballast model in ballast damage identification and the potential applicability of the selected ballast model in field track monitoring.

Keywords

Markov chain Monte Carlo Bayesian modeling non-linear behavior uniaxial stress–strain models ballast damage detection

Introduction

The primary concerns in the railway industry are the safety and comfort of passengers, reductions in the journey time, and consequently increased rail trips. The increase in maximum axle load and train speed results in enormous quasi-static and dynamic loads, which causes an increased rate of ballast abrasion, grain fracturing, and ballast contamination.¹ Under traffic loading, granular geomaterials possess non-linear and time-dependent elastoplastic responses such that their mechanical behavior depends mainly on the loading stress path.^2,3

A prerequisite for the favorable implementation of the mechanics-based analysis (e.g., damage detection and structural health monitoring) is that the behavior of the constituent materials is appropriately understood.^2,4 However, the true nature of the deformation mechanism of particles in granular layers is very complex and not yet fully understood. The vagueness is evidenced in the often-conflicting conclusions about the effects of some parameters on the behavior of granular materials under repeated loading, resulting in a vast number of models. Therefore, there arises a need to develop general models with robust prediction capabilities.

Although the accumulation of plastic deformation under the effect of loading is evident, the majority of permanent way engineers and researchers consider the ballast as a linear-elastic medium in rail track models.⁵ While the linear-elasticity assumption may be feasible for deep layers of the track (i.e., subgrade), the high level of loading fluctuation and stress changes in the ballast layer, which receives load directly from the sleeper, subject it to variations in strain and stiffness level. Hence, for capturing the general track mechanics, non-linear evaluation of the ballast material become almost inevitable. In the studies carried out by Paixão et al.³ and Varandas et al.⁶ using a 3D track model validated with experimental data, it was observed that the linear-elastic model was only able to provide accurate estimates of the overall track response. Meanwhile, the stress estimations (peak stress and stress evolutions) inside the ballast layer were very poor. It is thus unambiguous that generalizing the linear-elastic theory approximation undermines the analysis of ballast resiliency and deformation. Misrepresenting or neglecting the non-linearity phenomenon of the ballast layer will lead to unreliable estimation of the stiffness and stress of the ballast under a sleeper, resulting in an inefficient track health evaluation.

The most common approach in the literature is to construct models of ballast behavior based on laboratory triaxial testing using constant or varying confining pressure. Among the earliest models, Monismith et al.⁷ proposed a power law model of the resilient modulus based solely on the effect of the confining stress for cases where the deviator stress does not cause excessive plastic deformation. A more acceptable model for the stress dependence of material stiffness, known as the K-θ model, was proposed by Brown and Pell.⁸ The simple hyperbolic relationship assumes a constant Poisson’s ratio to calculate radial strain. The K-θ model was extended by other researchers like Uzan⁹ to include the deviator stress and material density effect. Other methods considering the bulk and shear moduli,¹⁰ Hertz power law,¹¹ and the hardening soil model¹² have also been proposed in the literature. The limitations attributed to models formulated by curve-fitting of large-scale experimental tests are undesirable. Thus, the numerical method of formulating ballast stress–strain is pursued further in this article.

For studying the structural behavior, optimizing the designs, and monitoring the condition of structures (railway tracks inclusive), numerical models are often adequate tools used for research and in industries. However, for such empirical models based on theories and general behaviors to provide valuable and reliable information, proper calibration and validation are necessary.³ Damage detection using vibration-based nondestructive testing (NDT) by model updating of constructed finite element (FE) models is also being increasingly studied utilizing the Bayesian theory.¹³ Although investigations in deterministic and linear-elastic model updating are mature, studies on probabilistic model updating and structural identification in non-linear models (e.g., reference 14) lag far behind. In general, the few studies involving model updating of non-linear models in the literature are usually either overly simplified or performed in the frequency domain, majorly due to the computation cost. In this article, model updating is carried out in the time domain using the track system acceleration measurements, owing to the nature of material nonlinearity considered at the local level.

The problem of ballast deformation and degradation is more of a geotechnical problem than a structural one. Hence the concentration of studies on the use of geomechanics method^12,15–18 and only a few studies from the structural viewpoint.^19–23 Although less prevalent than in the field of structural mechanics, Bayesian statistics–based methodologies have found usage in geomechanical problems related to optimal modeling and parameter identification (e.g., references 24,25). Few studies (e.g., references 26,27) have also implemented Bayesian modeling in assessing rail track geometry degradation. The uncertainty information of the damage identification results obtained by Bayesian analysis is particularly important for engineers and infrastructure managers who are trying to make judgments about remedial work, especially in situations where the damage is in inaccessible locations.²⁸ In this article, the time domain–enhanced Markov chain Monte Carlo (MCMC) simulation method²⁹ initially developed in reference 30 is employed. By drawing samples from the kernel density estimates of the posterior probability density function (PDF), the scheme can approximate the marginal PDFs of uncertain model parameters, as well as quantify model classes’ uncertainties, by which the most suitable model can be selected objectively. With the uncertainties of model parameters and models quantified, the resulting ballast model can then be extrapolated with a certain level of confidence for conditions different from those investigated herein.

The main contributions of this article are the development of non-linear models for representing the mechanics of the ballast layer based on ballasted track using in situ sleeper vibration data, as well as presenting an objective Bayesian theory–based methodology for pinpointing the most suitable non-linear model from the constructed models under varying load amplitudes, for effective ballast damage. In the Non-Linear FE Model of Ballasted Track section, the theoretical background, track model, and non-linear ballast models used in this study are presented. The MCMC-based Bayesian statistical framework for parameter identification and model class selection is introduced next in the Bayesian Statistical Framework section. Finally, in the Verification of Proposed Methodology Using In Situ Sleeper Measurements section, the proposed methodology is demonstrated using experimental data from three impact cases and five track model classes.

Non-linear FE model of ballasted track

Ballasted track FE model

In this study, a section of ballasted track across a single sleeper is modeled based on Winkler’s assumption of a beam on a single layer of elastic foundation,^31,32 considering the track substructure (ballast and subgrade) as the foundation. Although several more complicated models, which consider the foundation as multiple layers of elastic materials as well as complex continuous isotropic elastic body, are available, Winkler’s theory has proven adequate in calculating stresses and deflections in railroad tracks.³³ In this study, Winkler’s theory is used to model the ballast layer despite its simplicity. The track system consists of the rails, rail pads and fasteners, pre-stressed concrete sleepers, the ballast layer, and the subgrade, as in Figure 1. The framework of the 2-dimensional (2D) track model was developed in MATLAB^®, and the full details can be found in the authors’ previous work.³⁴ Overall, the 2D track FE model consists of 99 elements (48 Timoshenko beam elements, 49 spring elements, and 2 mass–spring elements), 100 nodes (51 moveable and 49 fixed), and 100 DOFs (51 translational and 49 rotational). The nominal values of the track model parameters are provided in Table 1. The behavior of the ballast springs in the FE model depends on the adopted engineering material law. To improve the modeling accuracy, the ballast springs are considered as possessing nonlinearity (in line with the actual behavior of geomaterials).

Table 1:

Track model parameters.

Track parameters	Nominal values
Modulus of elasticity of the sleeper	38.45 kN/m²
Density of concrete sleeper	2750 kg/m³
Total mass of the sleeper	305 kg
Poisson’s ratio of concrete	0.25
UIC60 rail mass	60 kg/m
Rail-fastening stiffness	50 MN/m
Track system damping ratio	18%

Modeling the uniaxial nonlinearity of ballast using elastoplasticity theory

The enormous number of models available for modeling the non-linear resiliency of ballast in the literature are basically from two sources: empirical models based on measured data from triaxial tests (e.g., reference 12), and theoretical models, which are first analytically derived and then validated (or verified) using data from laboratory tests (e.g., references 5,35). Dawson et al.³⁶ attribute the variations in the models available in the literature (from both sources) to the quality and constraints of the apparatus, as well as the testing procedures. Other limitations of some triaxial test-based models include the difficulty of replicating test conditions in the field, the need for high sophistication to understand granular behavior, and the large number of model parameters required in most cases. For example, the elastoplastic model developed in Sun et al.⁵ requires 15 parameters. On the other hand, empirical models are entirely based on principles that may be different from actual structural behaviors.

Among the defining factors that influence the mechanical behavior of ballast, Wang and Markine³⁷ reported that the maximal stress is the most commanding as it controls the permanent settlement of the ballast layer. Vertical stress is the most prominent form of stress from the effect of rolling wheel loads on unbounded aggregates.² Therefore, only the vertical ballast stress and strain are considered in this study. In the literature, the ballast has been demonstrated to possess work-hardening behavior.^5,12,26 The most common methods used in empirical models are the hyperbolic relationship^18,38 and the exponential/power law.^11,12,22 In conformity with previous works,³⁹ ballast models based on power law, exponential, hyperbolic, and linear functions were developed in this study.

One major development in this article is the formulation of elastoplastic models (involving both linear-elastic and plastic behaviors) for the ballast layer of the railway track. Essentially, the elastoplastic models are obtained by integrating an initial linear-elastic regime to some well-known non-linear uniaxial plasticity models found in the literature. Two categories of non-linear models were considered in the current study, viz:

Uniaxial plastic models which involve total strain: For this category of models, the stress–strain curve starts at the origin $(ε = 0, σ = 0)$ such that the total strain is used in the original non-linear model. Hence, to obtain the elastoplastic model, the part of the original model curve between the origin and the yield point $(ε = ε_{y}, σ = σ_{y})$ is cut and replaced by a chord. This approximation introduces two different stress–strain moduli along the full path of the elastoplastic model—the linear-elastic modulus for the linear part which is equivalent to the secant modulus of the original non-linear curve calculated at the yield point and the non-linear modulus which is obtained as tangent modulus in the plastic regime.

Uniaxial plastic models which involve only the plastic strain: For this category of models, the original non-linear model begins at the yield point $(ε = ε_{y}, σ = σ_{y})$ , instead of the origin $(ε = 0, σ = 0)$ as in the first category. Hence, to obtain the elastoplastic model, a straight line is inserted from the origin to the yield point to represent the linear-elastic portion of the model.

In both categories, it is demonstrated that the point which ensures the fusion of linear-elasticity and plasticity is the yield point. It is thus necessary to enforce the equality of strain and stress values at this point, even though the corresponding slope may be discontinuous. In this article, four different uniaxial elastoplastic material models with strain-hardening properties were investigated to model the general uniaxial stress–strain behavior of the ballast layer. The considered models have their origins from the Ludwik model,⁴⁰ Swift model,⁴¹ Prager model,⁴² and Voce model.⁴³ Since all of these models have been modified by introducing the linear-elastic path, they are named as the modified Ludwik, modified Swift, modified Prager, and modified Voce models in this study. For comparison purposes, the pure linear-elastic model (σ = Eε) was also included amongst the investigated ballast models, bringing the total number of models to five. The schematic representations and basic relationships for all considered models are presented in Figure 2. As shown in Figure 2, the modified Ludwik and modified Prager models are of the first category, while the modified Swift and modified Voce models are of the second category.

Modified Ludwik model

The stress–strain curve in the original Ludwik model is a function of the total strain ε and starts from the origin $(ε = 0, σ = 0)$ . In the modified Ludwik model, the initial portion of the curve up to the yield point $(ε = ε_{y}, σ = σ_{y})$ is replaced by a straight line,^29,44 as shown in Figure 2(a). Note that the slope changes discontinuously from $E$ to $n E$ at the yield point. The stress–strain equation can be expressed as

σ = {\begin{array}{l} E ε & for ε \leq ε_{y} \\ σ_{y} {(\frac{ε}{ε_{y}})}^{n} & for ε \geq ε_{y} \end{array}

(1)

where

n

is the degree of nonlinearity, and other parameters were previously defined. The fusion of the linear-elastic and the plastic portions of the curve reduces the number of parameters involved in model updating as the hardening strength

K

in Figure 2(a) has been replaced in Equation (1) by

σ_{y} / ε_{y}^{n}

. Unlike the original model, the modified Ludwik model does not suffer from possible misrepresentation in the stress–strain curve^29,44 since the modification ensures that the stress–strain curve passes through the yield point in a manner that the tangent modulus

n (σ / ε)

never surpasses the modulus of elasticity

E

Modified Swift model

The Swift model is also a power law–based stress–strain model. However, incorporating a parameter $- ε_{o}$ , which is the intercept of the curve on the $ε$ -axis (as shown in Figure 2(b)), differentiates the Swift model from the Ludwik model. $ε_{o}$ could be considered as a measure of plastic pre-strain of the material. The stress–strain equation is given by

σ = {\begin{array}{l} E ε & for ε \leq ε_{y} \\ σ_{y} {(\frac{ε_{0} + ε}{ε_{0} + ε_{y}})}^{n} & for ε \geq ε_{y} \end{array}

(2)

Indeed, Equation (2) is similar to Equation (1), and Equation (2) can be obtained by shifting the $σ$ -axis of Equation (1) along the positive $ε$ direction by $ε_{o}$ . By imposing the yield stress at the end of the linear-elastic phase to the beginning of the plastic phase (as shown in Figure 2(b)), the hardening strength $K$ can be expressed (through Equation (2)) as $σ_{y} / {(ε_{0} + ε_{y})}^{n}$ . The tangent modulus for the non-linear part of the modified Swift model is $n σ / (ε + ε_{0})$ .

Modified Prager model

The Prager model is the only hyperbolic function–based model investigated in this article. Its stress–strain curve starts at the origin with the initial slope of $σ_{s} / ε_{y} = E_{s}$ as shown in Figure 2(c), where $σ_{s}$ is the asymptotic yield stress and $ε_{y}$ is the elastic yield strain. The curve then gradually bends over and approaches the asymptotic yield stress. At the elastic yield strain, the stress is equal to $\tan h (1) σ_{s} \approx 0.76 σ_{s}$ , which is considered as the elastic yield stress $σ_{y}$ in the modified Prager model. Although the original Prager model is defined along the stress path by a single equation, a modified variant involving the linear-elastic representation of the strain up to the yield strain $ε_{y}$ is developed, in line with other considered models in Figure 2. The total stress–strain equation becomes

σ = {\begin{array}{l} E_{s} ε \tan h (1) & for ε \leq ε_{y} \\ σ_{s} \tan h (\frac{ε}{ε_{y}}) & for ε \geq ε_{y} \end{array}

(3)

In the modified Prager model, the hyperbolic curve from the origin to the yield stress of the original Prager model is replaced by the linear-elastic chord with a slope of $E = σ_{y} / ε_{y} = E_{s} \tan h (1)$ . Beyond the yield strain, the tangent modulus is $E_{s} (1 - {(σ / σ_{s})}^{2})$ .

Modified Voce model

As for the Voce model, an exponential function of the total strain and the strain-hardening parameter $N$ is involved. A definite yield point is assumed in the modified Voce model, and its formulation can be written as

σ = {\begin{array}{l} E ε & for ε \leq ε_{y} \\ \frac{σ_{y}}{1 - m} (1 - m e^{N (ε_{y} - ε)}) & for ε \geq ε_{y} \end{array}

(4)

where

e

is the exponential constant,

m

is the initial state of hardening of the material, while

N \in ℝ^{+}

is the rapidity at which the asymptotic stress

S

is approached (see Figure 2(d)). At the yield point, the slope of the curve changes from a constant value

E

to a linearly varying pattern defined as

N (S - σ)

. For reduction in model parameters, the asymptotic stress

S

(see Figure 2(d)) has been expressed as

σ_{y} / (1 - m)

in Equation (6) by incorporating the yield point. The Voce model is similar to the Prager model in a sense that both non-linear curves converge to a given stress (i.e.,

S

for the Voce model and

σ_{s}

for the Prager model). From another viewpoint, the Voce model is similar to the Swift model as the non-linear curve of both models extends to the negative strain region when the stress is very low or even zero.

Following the derivations in Adeagbo et al.,²⁹ the stiffness of the ballast

k_{B}

in both the linear-elastic and plastic vibration states of the elastoplastic models can be expressed as

k_{B} = {\begin{array}{l} \frac{b}{h} E & for ε < ε_{y} \\ \frac{b}{h} (\frac{Δ σ}{Δ ε}) & for ε \geq ε_{y} \end{array}

(5)

where

b

is the width of the sleeper in contact with the ballast and

h

is the original ballast thickness (see Figure 1). The expression

Δ σ / Δ ε

in Equation (5) is the tangent modulus which was obtained in sections 2.2.1–2.2.4 for all considered models. The nominal values of all ballast model parameters are provided in Table 2. Note that the ballast stiffness will change over time due to the accumulated traffic load and ballast degradation.⁴⁵ For unsupported or poorly supported sleepers, the reduction in ballast stiffness can be drastic.^46–49 The ballast damage detection methodology presented in this article provides a possible way to estimate the reduction in ballast stiffness based on the measured vibration of the sleeper.

Table 2.

Ballast model parameters.

Notations	Calibrating factors	Parameters	Nominal values
$ε_{y}$	$θ_{ε_{y}}$	Yield strain	1 × 10⁻⁴
$ε_{o}$	$θ_{ε_{o}}$	Measure of plastic pre-strain	1 × 10⁻⁴
$n$	$θ_{n}$	Degree of hardening	0.5
$m$	$θ_{m}$	Initial state of hardening	0.5
$N$	$θ_{N}$	Rapidity of approaching the asymptote	100
$k_{B}$	$θ_{k_{B}}$	Ballast stiffness	150 MN/m²

In all the considered models, Equations (1)–(4) and their representations in Figure 2 are only for the loading path but not the unloading one. At the removal of the load, for a ballast spring stressed into the plastic range, the unloading path is assumed to follow a linear-elastic path with modulus $E$ , thus resulting in an irrecoverable/permanent plastic strain. When the impact force is large, the vertical vibration amplitude at certain locations along the sleeper may be too large, resulting in a separation between the bottom of the sleeper and the top of the ballast. This phenomenon is implemented in the proposed model using tensionless spring to prevent negative stress, especially during the unloading process.

Bayesian statistical framework

Parameterization and description of model classes

In predicting the responses of the track system $\tilde{\ddot{z}} (t_{n})$ , resulting from input $F (t_{n})$ at a given time $t_{n}$ at all observed DOFs, it is necessary to consider the uncertainties arising from errors in the model class, parameter values, measurement errors, etc. In this study, the Bayesian method of stochastic modeling^50–52 is followed such that the system behavior can be described by a model class $M_{j} : M_{j} \in M$ and $M$ contains $N_{m}$ classes of models. To achieve this, the $j$ ^th model class $M_{j}$ is parameterized with a vector $θ_{j}$ such that the concept of uncertainty can be handled from the fundamental level of model parameters. The model class $M_{j}$ is then used to define a probability model $p (\tilde{\ddot{z}} | F, θ_{j}, M_{j})$ and prior probability density function (PDF) $p (θ_{j} | M_{j})$ , which describes the plausibility of the probability model.

In this work, all parameters which are likely to change by local conditions, load pattern, or damage are considered uncertain and must be updated, based on the measured data. For each uncertain model parameter, a dimensionless calibrating factor is assigned such that the product of the calibrating factor and the nominal value of the model parameter results in the actual value. Thus, these calibrating factors are the elements in $θ_{j}$ for the model class $M_{j}$ . The left and right rail masses in Figure 1 are calibrated with $θ_{M_{r}}$ and $θ_{M_{l}}$ , respectively. The rail pad stiffnesses at both rail seats are considered to be equal and parameterized by $θ_{k_{p}}$ . For the ballast regions created by grouping the ballast segments, a calibrating factor $θ_{k_{B i}}$ for $i$ = 1, …, $N_{R}$ calibrates the stiffness of the $i$ th discrete region in the linear-elastic regime, where $N_{R}$ is the number of ballast regions. The damping ratio for the entire track system is parameterized by $θ_{c}$ , to consider the classical system damping matrix.

The uncertain calibrating factors for the model classes are defined by the parameters of each ballast model involved, as in Table 3 and discussed in section 2.2. As a result of the varying complexities of the non-linear ballast models, these calibrating factors are all considered to be elements of a row vector

θ_{M_{j}} \in ℝ^{1 \times N_{p, j}}

for

j

= 1 to

N_{m}

, where

N_{m}

is the number of ballast models considered. Hence, the model classes are explicitly dependent on the uncertain model parameters considered for each ballast model. The total number of uncertain parameters considered for the

j

th model class is then

N_{j} = N_{p, j} + N_{R} + 4

, and the vector of calibrating factors is expressed generally as

θ_{j} = {θ_{M_{j}}, θ_{c}, θ_{M_{r}}, θ_{M_{l}}, θ_{k_{p}}, θ_{k_{B i}}}, for i = 1, \dots, N_{R} and j = 1, \dots, N_{m}

(6)

Table 3.

Model classes and the parameters of ballast models considered in this study.

$M_{j}$	Models	Parameters	$θ_{M_{j}}$
$M_{1}$	Ludwik	$ε_{y}, n$	${θ_{ε_{y}}, θ_{n}}$
$M_{2}$	Swift	$ε_{y}, n, ε_{o}$	${θ_{ε_{y}}, θ_{n}, θ_{ε_{o}}}$
$M_{3}$	Prager	$ε_{y}$	${θ_{ε_{y}}}$
$M_{4}$	Voce	$ε_{y}, N, m$	${θ_{ε_{y}}, θ_{N}, θ_{m}}$
$M_{5}$	Linear-elastic	—	${}$

In Equation (6), $θ_{j} \in ℝ^{1 \times (N_{p, j} + N_{R} + 4)}$ is a row vector. For all five model classes considered in this study, the corresponding ballast model parameters and vectors $θ_{M_{j}}$ are provided in Table 3. The calibrating factors will be referred directly to as uncertain parameters in subsequent sections for easy discussion. For simplicity, the joint prior PDF is taken as the product of the individual uniform prior PDFs for each parameter in $θ_{j}$ , assuming independency. Considering the nominal values in Table 2, $θ_{m}$ and $θ_{n}$ are bounded by [0.0 2.0], $θ_{c}$ is bounded by [1.0 1.5], $θ_{M_{r}}$ and $θ_{M_{l}}$ are bounded by [0.00 2.50], $θ_{k_{p}}$ is bounded by [0.4 2.0] according to recommendations, and $θ_{k_{B i}}$ is bounded by [0.0 1.0]. For the other model parameters, which are very susceptible to changes with impact magnitude, the bounds are obtained in each impact case by rounding up (down) the maximum (minimum) values at the final sampling level for each model class to the nearest 0.5.

Bayesian identification of model parameters in the time domain

Using a measured dataset of the system responses and input over a specified time range defined as $D = {\tilde{\ddot{z}}, F}$ , and assuming stationary covariance and independence of the responses at $N_{o}$ observed degrees of freedom (DOFs), the PDF of the dataset conditional on a model class $M_{j}$ can be expressed using a multivariate Gaussian PDF based on stochastic embedding^50–52

p (D | θ_{j}, M_{j}) = {(2 π κ^{2})}^{- N_{o} N_{t} / 2} exp [- \frac{1}{2 κ^{2}} \sum_{t = 1}^{N_{t}} ‖ \tilde{\ddot{z}} (t) - \ddot{z} {(t, θ_{j}) ‖}^{2}]

(7)

where

κ^{2}

is the prediction error variance and

\ddot{z} (t, θ_{j})

denotes the model-predicted responses at the

t

^th time step. The dependency of the responses on the model parameters,

θ_{j}

, and so as the model class,

M_{j}

, is emphasized. The joint probability model in equation (7) is referred to as the likelihood. From equation (7), the measure-of-fit function can be defined as

J (θ_{j} | D, M_{j}) = \frac{1}{N_{o} N_{t}} \sum_{t = 1}^{N_{t}} ‖ \tilde{\ddot{z}} (t) - \ddot{z} {(t, θ_{j}) ‖}^{2}

(8)

An updated PDF for the uncertain parameters $θ_{j}$ conditional on the model class, and system data can be constructed based on Bayes’ theory^51,52 using the likelihood and the prior PDF

p (θ_{j} | D, M_{j}) = c p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j})

(9)

where

p (θ_{j} | D, M_{j})

is the posterior PDF, and c ensures that the posterior PDF integrates to unity over the parameter space. By substituting equations (7) and (8) into Equation (9), the updated joint PDF of the uncertain calibrating factors can be expressed as

p (θ_{j} | D, M_{j}) = c_{0} exp [- \frac{N_{o} N_{t}}{2 κ^{2}} J (θ_{j} | D, M_{j})]

(10)

where

c_{0}

combines all the constants and prior PDF. From the updated PDF, the most probable values (MPVs) of each uncertain parameter in

θ_{j}

which maximize the posterior PDF as well as their respective uncertainties can be determined.

Model-class assessments in the time domain

Selecting the “best” ballast non-linear stress–strain model (defining the model class) from amongst competing models can be impartially achieved by utilizing the uncertainties as well as the plausibility of each model class so obtained. In general, the objective becomes that of selecting a model class $M_{j}$ from a set of candidate model classes $M$ .^53,54 By Bayes’ rule, the posterior probability of a model class $M_{j}$ conditional on the measured system data and the set of model classes is given as

P (M_{j} | D, M) = c_{1} p (D | M_{j}) P (M_{j} | M)

(11)

where

c_{1}

is a normalizing constant,

p (D | M_{j})

is the evidence of the model class

M_{j}

, and

P (M_{j} | M)

is the prior probability of the model class, which can be taken as

1 / N_{m}

if all model classes are equally probable. Thus, Equation (11) reduces to

P (M_{j} | D, M) \propto p (D | M_{j})

(12)

Hence, maximizing the evidence of the model class $M_{j}$ based on system data $D$ is, therefore, in this case, the same as maximizing the posterior probability of the model class, to obtain the most optimal model class.⁵⁵ By the total probability theorem, the evidence of the model class, $M_{j}$ , is expressed as

p (D | M_{j}) = \int^{} p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j}) d θ_{j}

(13)

The evaluation of equation (13), just like Equation (10), is cumbersome. An approach to simplify these computations is presented subsequently in the following section.

Enhanced Markov chain Monte Carlo–based sampling scheme

Toward evaluating the multi-dimensional integral of the posterior PDF in equation (10) and the evidence of a model class in equation (13), an MCMC-based sampling scheme developed in Lam et al.’s work³⁰ following the idea in reference 56 and modified in Lam et al.’s work³⁴ is employed. Readers are directed to references 30 and 34 for the detailed formulation, and only the most important equations together with a brief description are given below for the completeness of this article.

The general procedures involved, as well as the convergence criteria, are presented in Figure 3. For smooth convergence to the desired joint posterior PDF, the sampling process is divided into several levels interconnected by “bridge PDFs,” which take the form of equation (10). The popular Metropolis–Hastings algorithm^57,58 is used to draw $N_{s s}$ samples in each sampling level, with the proposal PDF being a uniform PDF (at the first level) or a kernel density estimate of the posterior PDF (in subsequent levels). The kernel density estimate at a sampling level $r$ is essentially a weighted multivariate Gaussian PDF using samples from the preceding $(r - 1)$ th level and is given by

k_{r} (θ_{j}) = \frac{1}{N_{s}} \sum_{g = 1}^{N_{s s}} {(w λ_{r - 1}^{(g)})}^{- N_{j}} N (S_{r - 1}^{(g)}, w λ_{r - 1}^{(g)} Σ_{r - 1})

(14)

where

w

is the window bandwidth,

λ_{r - 1}^{(g)}

is the

g

th local bandwidth at the

(r - 1)

th level,

N_{j}

is the total number of calibrating factors in

θ_{j}

S_{r - 1}^{(g)}

is the

g

^th vector of the MCMC samples of

θ_{j}

at the

(r - 1)

th level, and

Σ

is the covariance matrix. The prediction error variance

κ^{2}

follows a geometric sequence, with a reduction in values by a factor

α

at each level. By following the series and evaluating the estimated prediction error variance

{\hat{κ}}^{2} \approx J ({\hat{θ}}_{j} | D, M_{j})

as in Equation (8) at each sampling level, the theoretical “required” number of sampling levels,

r_{r e q}

, can thus be determined as the first stopping criterion for the MCMC²⁹

r_{req} = int (1 + \frac{\log_{10} (κ_{1}^{2} / {\hat{κ}}_{r}^{2})}{\log_{10} (α)})

(15)

where

κ_{1}^{2}

is the initial prediction error variance,

{\hat{κ}}_{r}^{2}

is the estimated prediction error variance at the

r

th sampling level, and the

int (\cdot)

function returns the integer part of a real number. The convergence of the MCMC is further diagnosed using the procedures on the left-hand side of Figure 3, by observing the relative change in the prediction error variance over a number of sampling levels

N_{c} = int ({(\log_{10} α)}^{- 1})

. In Figure 3,

ξ

is the threshold value (herein 1%), and

ψ

is the threshold reduction factor. The optimal MCMC sampling level

r_{opt}

is not taken as the final sampling level

r_{e n d}

but selected within the convergence region

[r_{end} - N_{c} + 1, \dots, r_{end}]

as the level closest to the computed

r_{req}

. With the samples generated at the optimal level

r_{opt}

, the corresponding kernel density estimate approximates the posterior PDF. For the implementation of the MCMC scheme and convergence diagnostics, readers are directed to Adeagbo et al.²⁹

In order to evaluate the evidence in Equation (13), the samples of the posterior PDF generated from the model updating process are utilized since samples from the prior PDF could be concentrated in unimportant regions. Equation (13) is then re-expressed in terms of the posterior PDF as

p (D | M_{j}) = \int \frac{p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j})}{p (θ_{j} | D, M_{j})} p (θ_{j} | D, M_{j}) d θ_{j}

(16)

Applying MCMC approximation using samples of the posterior PDF and taking the natural logarithm to prevent numerical overflow problems, the evidence becomes⁵⁹

\ln (p (D | M_{j})) \approx \frac{1}{N_{s s}} (\sum_{i = 1}^{N_{s s}} \ln (p (D | θ_{j}, M_{j})) - \sum_{i = 1}^{N_{s s}} \ln (\frac{p (θ_{j} | D, M_{j})}{p (θ_{j} | M_{j})}))

(17)

Hence, the plausibility of each model class is dependent on how well the model can replicate measured responses (i.e., the likelihood as the first term on the right-hand side) and on a penalty function. The penalty function is the second term on the right-hand side, also known as the Kullback–Leibler information (K–L information), increases exponentially with the number of model parameters (i.e., calibrating factors) and with the amount of information extracted.^50,60,61 With the penalty function, simpler models that are reasonably consistent with data should be preferred over more complex models that lead to only slightly better agreement with the data.

Verification of proposed methodology using in situ sleeper measurements

The developed classes of uniaxial models for the ballast springs’ behavior are adjudged based on real track vibration data following the Bayesian model class selection method. Their applicability to damage detection in ballasted tracks is also investigated.

Experimental set-up and procedures

The impact hammer test was carried out on a full-scale indoor test track section under laboratory conditions, as shown in Figure 4(a). The nominal values of the track section parameters are provided in section 2. Ballast damage was artificially simulated under the right-end and middle of the sleeper (two-third of the sleeper length) by replacing the normal-sized ballast (∼60 mm) with small-sized ballast (∼15 mm), as shown in Figure 4(b) and (c); thus, three ballast regions are considered in the track model. The ballast in the right one-third of the sleeper was less compacted compared to other regions. Hence, in equation (6), the ballast stiffness calibrating factors $θ_{k_{B 1}}$ , $θ_{k_{B 2}}$ , and $θ_{k_{B 3}}$ are utilized for the right, middle, and left ballast regions, respectively. The selected sleeper was within the 8-sleepers panel (see Figure 4(a)) and not at the extremes, to prevent boundary effects. Also, as the monitoring and failure probability estimation of railway infrastructures can be adversely affected by environmental factors, the temperature and relative humidity in the laboratory was kept at 26°C and 70%, respectively, throughout the whole test.

The sleeper is excited near the right end (at observation point 2 (OP2) in Figure 4(c)) using a Dytran 5805A instrumented hammer (which has a load cell fixed to the hammer head) with varying impact forces. The vibration responses of the sleeper were measured using 13 uniaxial accelerometers assembled on the top surface of the sleeper along the centerline (on square-shaped steel sheets), as in Figure 4(c). The piezoelectric accelerometers are of Type 8776A50M3 with a sensitivity of 100 ± 5 mV/g and measurement range of ±50 g. The sleeper vibration responses were measured at a sampling frequency of 6.4 kHz, for 20s with several hammer impulses. Since time-domain analysis is carried out, the measured responses from a single hammer hit are enough for model updating and model class selection. Vibration responses for ∼0.025s were utilized in the model updating since the sleeper vibration dampens out already within this period, owing to the high system damping properties. However, to ensure that the utilized data are representative of the system behavior and for the purpose of studying the effect of different levels of impulse magnitudes, several hammer hits with different impact forces were conducted. The measured data are transferred to the computer using the National Instruments analog/digital modules and chassis. Some of the experimental equipment are shown in Figure 4(a), and readers are directed to references^34,62 for further details.

Impact case 1: Hammer impulse of ∼1.40 kN

The proposed methodology was implemented for model updating and model class selection using the five considered ballast stress–strain models utilizing the in situ sleeper’s response data under a relatively low hammer impact of ∼1.40 kN, as shown in Figure 5(a).

With the 5000 MCMC samples generated at the optimal level, the Bayesian model class selection was carried out on all five model classes (using equation (17)), and the results are presented in Table 4. The best choice among model classes will be those with the higher evidence values, leading to higher posterior probabilities. As the model class with the highest log-evidence value of 2078.0561, the linear-elastic model class

(M_{5})

is selected as the most plausible model class to capture the mechanics of the track. This selection confirms that the behavior of the ballast can be approximated as linear-elastic under small amplitude vibration.³⁴ It is clear from the table that all considered model classes have very similar prediction accuracies as evidenced by the measure-of-fit values

J ({\hat{θ}}_{j} | D, M_{j})

. The K–L information is observed to increase with the complexity of the model classes. However, even though the modified Swift and modified Voce model classes are both with 10 parameters, the modified Voce model class

(M_{4})

is the most penalized. Thus, it is expected that the modified Voce model class

(M_{4})

has a higher uncertainty level compared to the modified Swift model class

(M_{2})

Table 4.

Bayesian model class selection results for all model classes in Impact case 1.

Model classes	Models	$J ({\hat{θ}}_{j} \| D, M_{j})$	Log-likelihood	K–L information	Log-evidence
$M_{1}$	Ludwik	0.007764	2112.2665	37.9612	2074.3053
$M_{2}$	Swift	0.007738	2116.6443	39.6972	2076.9470
$M_{3}$	Prager	0.007777	2111.6316	33.6116	2078.0200
$M_{4}$	Voce	0.007757	2114.7114	42.8334	2071.8779
$M_{5}$	Linear-elastic	0.007827	2104.8539	26.7978	2078.0561

The bolded values are the best values under each headings... that is, the least J(theta), the highest Log-likelihood, the least K-L information, the highest log-evidence.

The marginal posterior PDFs for all uncertain parameters are constructed by using MCMC approximation on the samples at the optimal level for all model classes. The marginal PDFs for calibrating factors, which are common to all model classes, are compared in Figure 6. Generally, the common calibrating factors from different model classes are closely distributed in the same region. The shapes of the marginal PDFs are mostly non-Gaussian (some even with multiple peaks), which showcases the benefit of the enhanced MCMC-based Bayesian technique employed. It is also observed that the parameters for the modified Ludwik model $(M_{1})$ have marginal PDFs with relatively sharp prominent peaks compared to those of other model classes. Hence, it is very likely that this model class has the lowest uncertainties for most of its model parameters.

The range of sample distribution for $θ_{ε_{y}}$ in the modified Prager model class is different from that in other model classes (except the linear-elastic model class for which this parameter does not exist). The discrepancy in the range and mean of samples is due to the hyperbolic tangent function (in the Prager model), which ensures that the stress–strain curves rapidly turn into the non-linear region toward the asymptotic yield stress $σ_{s} \approx 1.3130 σ_{y}$ .

For discussion purposes, Gaussian PDFs were fitted to the marginal distributions to estimate their MPVs and associated uncertainties in terms of COV%. The results of model updating for all uncertain parameters are presented in Table 5, where the COV% values are included in parenthesis. From Table 5, it is observed that ballast model calibrating factors

θ_{n}

and

θ_{N}

possess the highest COV in the considered model classes (50.5751% and 39.8491%, respectively). In general, the calibrating factors for the ballast parameters (i.e.,

θ_{ε_{y}}, θ_{ε_{0}}, θ_{m}, θ_{n}, and θ_{N}

) have higher uncertainties compared to those of the system parameters. For all the model classes, the updated rail pad stiffness is within the 20–100 MN/m range stipulated for the dynamic stiffness of rubber rail pads in reference 63. The updated modal damping ratio (∼22%) is comparable with the result for the first modal damping ratio of in-field condition (30.41%) reported in Kaewunruen and Remennikov⁶⁴ for AUSTRAK narrow-gauge sleeper. The difference in these values may be due to the sleeper’s properties and underlying ballast conditions. The MPVs for the ballast stiffness calibrating factors in all three regions correctly portray the simulated damage scenario with all the model classes utilized. The stiffness of the middle and right ballast regions

(θ_{k_{B 1}} and θ_{k_{B 2}})

always have lower values compared to that of the left region

(θ_{k_{B 3}})

in all model classes considered, confirming the damage location and severity prediction ability.

Table 5:

The identified MPV and COV% of the MCMC-based Bayesian model updating in Impact case 1.

	Calibrating factors	$M_{1}$	$M_{2}$	$M_{3}$	$M_{4}$	$M_{5}$
Ballast parameters	$θ_{ε_{y}}$	2.4223 (0.8342%)	2.3782 (0.3329%)	2.3153 (0.9166%)	2.4096 (0.6432%)	—
	$θ_{ε_{o}}$	—	8.8335 (12.3424%)	—	—	—
	$θ_{m}$	—	—	—	0.3895 (17.0767%)	—
	$θ_{n}$	0.4896 (11.8660%)	0.1816 (50.5751%)	—	—	—
	$θ_{N}$	—	—	—	2.2459 (39.8491%)	—
System parameters	$θ_{c}$	1.1874 (0.1043%)	1.1739 (0.3109%)	1.1852 (0.2735%)	1.1910 (0.2500%)	1.2042 (0.4734%)
	$θ_{M_{r}}$	1.9053 (0.2534%)	1.8930 (0.5111%)	1.8811 (0.6512%)	1.8842 (0.7856%)	1.8965 (0.8845%)
	$θ_{M_{l}}$	0.4434 (0.6908%)	0.4127 (1.6017%)	0.4339 (1.3192%)	0.4315 (0.7251%)	0.4263 (3.1022%)
	$θ_{k_{p}}$	0.9002 (0.1730%)	0.8850 (0.3205%)	0.8841 (0.5067%)	0.8767 (0.3887%)	0.8892 (0.7690%)
	$θ_{k_{B 1}}$	0.2700 (0.5399%)	0.2754 (0.2904%)	0.2744 (0.6243%)	0.2717 (0.2172%)	0.2701 (0.8020%)
	$θ_{k_{B 2}}$	0.4327 (0.3745%)	0.4375 (0.8817%)	0.4299 (2.2069%)	0.4165 (0.9544%)	0.4328 (1.8044%)
	$θ_{k_{B 3}}$	0.5276 (0.7946%)	0.5172 (1.1605%)	0.5229 (1.2807%)	0.5221 (1.2305%)	0.5128 (2.0572%)

After model updating, the system responses (i.e., accelerations of the sleeper) were calculated for all model classes using the updated model parameters. In Figure 7, the time-domain responses at six selected representative OPs are compared for all considered model classes as well as the measured data. The observed time histories for all models are almost overlying each other, and they are all very close to the measured data at all observation points. It can be concluded that the system responses can be reasonably predicted by model updating of any of the considered model classes. Furthermore, some deviations are observed between the model-predicted responses and the measured data at observation points located far away from the impact location (i.e., OPs 10–13). In this segment of the sleeper, higher frequencies are dominant (owing to the ballast stiffness), and therefore, the single modal damping ratio assumed for all modes may be unsuitable.

For observing the non-linear stress–strain behavior of the ballast, the ballast springs’ stress and strain values for the updated system parameters are evaluated using various stress–strain models. A comparison of the evaluated stress–strain behaviors for the ballast segment 1 (beneath sleeper element 1) is presented in Figure 8. Interestingly, the peak stresses and irrecoverable plastic strains (for the non-linear models) vary from one ballast model to the other. The non-linear models (i.e., modified Ludwik, modified Swift, modified Prager, and modified Voce) demonstrate some non-linear vibrations even with the small impact force applied. However, the value of the plastic strains (∼1.3 × 10⁻⁵ maximum) is small, relative to the total strains. Thus, the linear-elastic model is deemed sufficient to model the ballast behavior. In Figure 8, the modified Prager model showcased the highest strain-hardening property, the second-highest peak stress, and the least maximum strain (∼2.56 × 10⁻⁴), while the linear-elastic model results has the highest stress of ∼1.32 × 10⁴ N/m².

Impact case 2: Validation with hammer impulse of ∼2.50 kN

For validation purposes, another set of vibration data was collected from the in situ sleeper under a higher hammer impulse of ∼2.50 kN (shown in Figure 5(b)). The use of a higher impact force (when compared to case 1) is intended to increase the vibration amplitude and so as to clearly demonstrate the non-linear behavior of ballast.

With the 5000 MCMC samples generated at the optimal level, the Bayesian model class selection was carried out on all model classes, and the results are presented in Table 6. According to the log-evidence values, the model classes can be grouped into two: the modified Prager

(M_{3})

, modified Voce

(M_{4})

, and linear-elastic

(M_{5})

model classes with high log-evidence values near 255.0; and the other model classes (modified Ludwik and modified Swift) with lesser log-evidence values. Generally, the modified Ludwik model class

(M_{1})

has the highest log-evidence value of 255.5732 and is thus selected as the most probable model class in this impact case. This selection confirms that at this level of vibration, the behavior of the ballast can no longer be approximated as linear-elastic. Again, all non-linear model classes have better measure-of-fit values,

J ({\hat{θ}}_{j} | D, M_{j})

, and higher log-likelihood values compared to the linear-elastic model. The trend of the log-likelihood values, however, does not fit perfectly with the trend of the measure-of-fit values. This is because the measure-of-fit function is calculated from a single model (i.e., the MPVs,

{\hat{θ}}_{j}

), while the likelihood is calculated from all MCMC samples.

Table 6.

Bayesian model class selection results for the considered model classes in Impact case 2.

Model classes	Models	$J ({\hat{θ}}_{j} \| D, M_{j})$	Log-likelihood	K–L information	Log-evidence
$M_{1}$	Ludwik	0.04177	286.1618	31.4539	254.7079
$M_{2}$	Swift	0.04177	286.3822	35.6943	250.6878
$M_{3}$	Prager	0.04195	282.9948	27.8364	255.1584
$M_{4}$	Voce	0.04175	286.2158	30.6426	255.5732
$M_{5}$	Linear-elastic	0.04219	279.0196	23.7396	255.2801

The bolded values are the best values under each headings... that is, the least J(theta), the highest Log-likelihood, the least K-L information, the highest log-evidence.

The marginal posterior PDFs for all uncertain parameters are constructed using MCMC samples at the optimal level for all model classes. The marginal PDFs for calibrating factors, which are common to all model classes, are compared in Figure 9. Generally, the common calibrating factors from different model classes are closely distributed in the same region. The shapes of the marginal PDFs are in general non-Gaussian, and the modified Ludwik model class $(M_{1})$ has multiple peaks in most model parameters.

For discussion purposes, Gaussian PDFs were fitted to the marginal distributions to estimate their MPVs and associated uncertainties in terms of COV%. The results of model updating for all uncertain parameter factors are presented in Table 7, from which it is observed that the ballast model calibrating factors

θ_{m}

and

θ_{n}

possess the highest COV in the considered model classes (105.7393% and 292.1037%, respectively). In general, the calibrating factors for the ballast model parameters (i.e.,

θ_{ε_{y}}, θ_{ε_{0}}, θ_{m}, θ_{n}, and θ_{N}

), again have higher uncertainties compared to those of the system structural parameters. The MPVs for the ballast stiffness calibrating factors in all three regions are similar to those obtained in the previous Impact case 1, except for

θ_{k_{B 3}}

which increased fairly from an average of 0.52 (in Impact case 1) to 0.61 in this case, owing to the larger impact force. Nonetheless, the stiffness MPVs correctly portray the simulated damage scenario with all the model classes utilized. The stiffness of the middle and right ballast regions

(θ_{k_{B 1}} and θ_{k_{B 2}})

always has lower factors compared to the left region with normal ballast

(θ_{k_{B 3}})

in all model classes considered.

Table 7.

The identified MPV and COV% of the MCMC-based Bayesian model updating in Impact case 2.

	Calibrating factors	$M_{1}$	$M_{2}$	$M_{3}$	$M_{4}$	$M_{5}$
Ballast parameters	$θ_{ε_{y}}$	2.9044 (1.2779%)	2.9468 (0.9144%)	2.8049 (2.0107%)	2.9234 (0.9869%)	—
	$θ_{ε_{o}}$	—	0.8174 (76.3904%)	—	—	—
	$θ_{m}$	—	—	—	0.2566 (105.7393%)	—
	$θ_{n}$	0.0493 (138.6467%)	0.0561 (292.1037%)	—	—	—
	$θ_{N}$	—	—	—	6.7223 (21.6553%)	—
System parameters	$θ_{c}$	1.0650 (1.1279%)	1.0719 (0.3094%)	1.0845 (0.4421%)	1.0741 (0.5306%)	1.0977 (0.7402%)
	$θ_{M_{r}}$	1.9156 (0.7172%)	1.8901 (0.6262%)	1.9040 (0.8833%)	1.9124 (1.1787%)	1.9223 (1.4345%)
	$θ_{M_{l}}$	0.4270 (2.7311%)	0.4178 (2.2023%)	0.4048 (3.3526%)	0.4265 (3.8443%)	0.4290 (4.3301%)
	$θ_{k_{p}}$	0.8636 (0.6645%)	0.8575 (0.5184%)	0.8627 (0.6863%)	0.8645 (1.1686%)	0.8725 (1.2568%)
	$θ_{k_{B 1}}$	0.2706 (1.3007%)	0.2624 (0.5874%)	0.2665 (0.9714%)	0.2659 (1.2171%)	0.2578 (1.4311%)
	$θ_{k_{B 2}}$	0.4357 (2.1468%)	0.4329 (1.1115%)	0.4322 (1.7895%)	0.4329 (2.6900%)	0.4298 (2.8363%)
	$θ_{k_{B 3}}$	0.6116 (3.4351%)	0.6197 (1.1506%)	0.6018 (1.8535%)	0.6108 (2.0494%)	0.6154 (2.6370%)

After model updating, the system responses (i.e., accelerations of the sleeper) were calculated for all model classes using the updated model parameters. In Figure 10, the time-domain responses at six selected representative OPs are compared for all considered model classes as well as the measured data. The observed time histories for all models infer that the system responses can be reasonably predicted by the model updating of any of the considered model classes. Furthermore, some deviations are observed between the model-predicted responses and the measured data at OPs located far away from the impact location (i.e., OPs 10–13) due to higher frequencies of vibration.

For observing the non-linear stress–strain behavior of the ballast, the ballast springs’ stress and strain values for the updated system parameters are evaluated by various stress–strain models. A comparison of the evaluated stress–strain behaviors for the ballast segment 1 (beneath sleeper element 1) is presented in Figure 11 for all models. Just like in Impact case 1, the peak stresses and irrecoverable plastic strains (for the non-linear models) vary from one ballast model to the other. The non-linear models (modified Ludwik, modified Swift, modified Prager, and modified Voce) demonstrate that higher plastic strain is involved with higher impact intensity. The highest value of the plastic strain is ∼4.0 × 10⁻⁵. In this impact case, the linear-elastic model again has the highest stress (∼1.61 × 10⁴ N/m²) and total strain (∼3.35 × 10⁴) level, while the modified Prager model has the highest nonlinearity but lowest plastic strain.

Impact case 3: Revalidation with hammer impulse of ∼2.84 kN

For revalidation purpose and to investigate the effect of impact loads on the model parameters of the investigated ballast models, model updating and model class selection were carried out based on the sleeper’s response data from a hammer impact of about 2.84 kN (the highest in Figure 5(c)).

The Bayesian–MCMC model class selection results for all five model classes are presented in Table 8. Herein, the modified Voce model class

(M_{4})

, with the highest log-evidence value of 41.5750, is again selected as the most probable model class. It is observed that the trend in model class rankings, in this case, differs from that in Impact cases 1 and 2. This discrepancy is consequent of the non-uniform response of the ballast to varying loading conditions, a general characteristic of granular soils as evidenced in Akeju et al.²⁴ The modified Swift model class

(M_{2})

is the most penalized, thus relegating it to the model class with the least log-evidence value, even below the linear-elastic model

(M_{5})

, which is believed to have the lowest fitting ability. The modified Prager model is almost as good as the selected modified Voce model.

Table 8.

Bayesian model class selection results for the considered model classes in Impact case 3.

Model classes	Models	$J ({\hat{θ}}_{j} \| D, M_{j})$	Log-likelihood	K–L information	Log-evidence
$M_{1}$	Ludwik	0.05409	70.5638	31.9327	38.6311
$M_{2}$	Swift	0.05411	71.0526	35.9362	35.1164
$M_{3}$	Prager	0.05439	66.9066	27.6978	39.2088
$M_{4}$	Voce	0.05408	71.0676	31.5232	39.5444
$M_{5}$	Linear-elastic	0.05471	60.9302	23.7147	37.2155

The bolded values are the best values under each headings... that is, the least J(theta), the highest Log-likelihood, the least K-L information, the highest log-evidence.

For model updating purposes, the marginal posterior PDFs of calibrating factors were obtained from the MCMC samples at the optimal level, and Gaussian PDFs were fitted to evaluate the corresponding MPVs as well as the COVs. The marginal posterior PDFs showing the spread of values (uncertainty) for the calibrating factors, which are common to all model classes, are presented in Figure 12. It is conspicuous that the marginal PDFs of calibrating factors for all model classes are far from being Gaussian, with some calibrating factors having multimodal peaks. The calibrating factor $θ_{ε_{y}}$ for the modified Prager model class has a flat multimodal shape, depicting its higher uncertainties, compared with other model classes. The ballast stiffness marginal PDFs for the Prager model class again shift away from those of other model classes to accommodate the hyperbolic function.

From amongst the updated calibrating factors, the MPVs of those factors related to the ballast hardening model

(θ_{ε_{y}}, θ_{ε_{0}}, θ_{m}, θ_{n} and θ_{N})

and the ballast stiffness (especially

θ_{k_{B 3}}

) are the most sensitive to the changes in the impact amplitude. This behavior further reinforces the deductions in Lekarp et al.² that the behavior of ballast is affected by local factors like stress path and amplitude. Consequent to the marginal PDF plots in Figure 12, Table 9 shows the MPVs and estimated COV for the calibrating factors of all model classes. The results of model updating show that the MPVs for the stiffness of the ballast in the three regions accurately portray the simulated damage scenario regardless of the strain-hardening model utilized. Slight increments are observed in the average values of the ballast stiffness calibrating factors, in this case, due to the higher load on the ballast. The stiffness calibrating factors for the middle and right ballast regions are very similar (with damaged ballasts) and are lower compared to the left region stiffness factor (with normal ballasts) in all considered models, thus demonstrating the effectiveness of the proposed methodology in ballast damage detection.

Table 9.

The identified MPV and COV% of the MCMC-based Bayesian model updating in Impact case 3.

	Calibrating factors	$M_{1}$	$M_{2}$	$M_{3}$	$M_{4}$	$M_{5}$
Ballast parameters	$θ_{ε_{y}}$	2.8771 (0.7082%)	2.8770 (0.5004%)	2.8240 (3.2575%)	2.8881 (0.9504%)	—
	$θ_{ε_{o}}$	—	4.7309 (9.1060%)	—	—	—
	$θ_{m}$	—	—	—	0.6953 (52.9357%)	—
	$θ_{n}$	0.0284 (123.2817%)	0.3325 (31.0706%)	—	—	—
	$θ_{N}$	—	—	—	1.3893 (55.7812%)	—
System parameters	$θ_{c}$	1.0884 (0.4168%)	1.0919 (0.4275%)	1.0974 (0.6251%)	1.0869 (0.4671%)	1.1144 (0.6159%)
	$θ_{M_{r}}$	1.8785 (0.6293%)	1.8906 (0.6597%)	1.9043 (1.0238%)	1.8925 (0.8088%)	1.9167 (1.4061%)
	$θ_{M_{l}}$	0.4707 (3.9269%)	0.4516 (3.4186%)	0.4593 (4.0151%)	0.4583 (3.0353%)	0.4618 (4.1533%)
	$θ_{k_{p}}$	0.8652 (0.6611%)	0.8638 (0.7178%)	0.8699 (0.8464%)	0.8659 (0.7997%)	0.8781 (1.1592%)
	$θ_{k_{B 1}}$	0.2711 (1.1610%)	0.2702 (0.5793%)	0.2665 (1.4952%)	0.2693 (1.5039%)	0.2587 (1.2881%)
	$θ_{k_{B 2}}$	0.4305 (2.1323%)	0.4527 (1.2625%)	0.4454 (2.3325%)	0.4440 (1.8277%)	0.4453 (2.8208%)
	$θ_{k_{B 3}}$	0.6859 (2.5237%)	0.6522 (1.4423%)	0.6696 (0.4275%)	0.6706 (2.2648%)	0.6672 (2.7525%)

After updating the selected modified Voce model class, the model-predicted vibration responses at the measured DOFs were calculated and compared with the measured data. In Figure 13, the responses from both the calibrated models and experimental data are plotted at six representative sensor locations (OPs 1, 3, 6, 8, 11, and 13), and only small discrepancies were observed.

With the MPVs of the model classes, the ballast springs’ stress and strain values for the updated system parameters are evaluated, and a comparison of the ballast spring behavior for the ballast segment 1 is presented in Figure 14. Amongst the compared ballast models, it was observed that the predicted maximum stresses for the non-linear models are close to one another. Thus, the variation in the stress–strain curves follows the trend observed in the measure-of-fit values in Table 8. The selected modified Voce model class $(M_{4})$ has the highest plastic strain value (∼4.25 × 10⁻⁴), but the lowest maximal stress (∼1.46 × 10⁴ N/m²) among the model class with nonlinearity. Meanwhile, the modified Prager model class $(M_{4})$ maintains the position as the elastoplastic model with the highest degree of nonlinearity and the highest maximal stress, while the stress–strain curves of the modified Ludwik $(M_{1})$ , modified Swift $(M_{2})$ , and modified Voce $(M_{4})$ model classes are close to each other. The linear-elastic model turns out to have peak values and behaviors very different from other models, hence demonstrating the need for modelers to incorporate ballast nonlinearity in track studies. The significance of the model so selected is also demonstrated in the stress–strain plots of the monitored ballast segments 1, 24, and 48, which are representatives of the ballast segments within the right, middle, and left ballast regions, as presented in Figure 15. It can be observed that the curve for the ballast segment 48 in the normal ballast region has the highest stress–strain slope. The ballast segment 24 is in the damaged region and so has a slope that is lower than that of ballast segment 48, but higher than that of ballast segment 1 (in the damaged ballast region with lower identified stiffness). Furthermore, the stress–strain plots indicate that only ballast segment 1, which is very close to the excitation point, experiences nonlinearity. In this impact case, the maximum ballast stress and strain obtained are 14.6334 kN/m² and 3.2651 × 10⁻⁴, respectively. Since the peak stress and strain evolution plays an important role in permanent ballast settlement, necessary cautions must be taken in selecting a suitable empirical stress–strain model for reliable long-term track monitoring studies, and the proposed methodology provides a feasible way to implement this in a practical and theoretical rigorous manner.

Discussion and conclusion

In this study, a methodology is proposed for the modeling of the uniaxial strain-hardening property of railway ballasts for ballast damage detection purposes following the Bayesian framework. For generality and better representation of the track section as a beam on elastic foundation, axial strain-hardening of the ballast layer in ballasted tracks was investigated using four newly developed empirical stress–strain models with elastoplasticity as well as the linear-elastic properties. For handling the uncertainties in the model parameters, an MCMC-based Bayesian scheme was employed for model updating utilizing a set of measured vibration data.

The Bayesian model class selection method was used to select the most plausible ballast model class, conditional on the set of measured system acceleration data from impact hammer tests. By employing the proposed methodology on three different impact hammer test datasets (with low, medium, and high impact force) from the same in situ sleeper, it was observed that the general system calibrating factors generally remain the same with an increment of about 1.10 kN (from Impact cases 1 and 2) and 1.44 kN (from Impact cases 1–3) in the impact force. However, the parameters, which are direct properties of the ballast layer or function of the ballast condition, are susceptible to changes in impact level, mainly due to the variation in the rate of ballast grain contact. In general, all investigated models using the three impact cases were able to infer the health condition of ballast underneath the measured sleeper, in terms of location and severity. The high level of uncertainties involved in some calibrating factors, as well as the non-Gaussian nature of the marginal distributions, makes a strong case for employing the enhanced MCMC-based model updating technique instead of conventional methods and the assumption of unimodal Gaussianity, as common in the literature. Even though the predicted time-domain responses at all observed DOFs are almost equal using various model classes, the effects in terms of stresses and strains are quite different.

In all impact cases considered, the behavior of the ballast layer ensured that different ballast model classes were selected. For a relatively small impact (∼1.40 kN), the linear-elastic model is selected as being sufficient to model the ballast springs, even though some plastic strains exist. However, with higher impact loads (∼2.50 kN and ∼2.84 kN) and higher nonlinearity demands, the modified Voce model class is the most plausible model class, based on the measured in situ sleeper vibration. The modified Prager model class performs very well in both the low impact and high impact cases. However, the modified Swift model is the worst non-linear model in all considered impact cases, owing to the large uncertainty introduced by the model parameters. Since the peak stress and stress evolution plays an essential role in permanent ballast settlement and sleeper design, necessary cautions must be taken in selecting a suitable empirical stress–strain model for reliable long-term track monitoring and reliability studies. The main contributions of this article are

to demonstrate the importance of considering the non-linear behavior of ballast in reducing modeling errors and ensuring accurate response predictions of the rail–sleeper–ballast system under large amplitude vibrations,

to develop a series of non-linear models for representing the mechanical behavior of the ballast utilizing the in situ sleeper vibration data,

to demonstrate the use of Bayesian approach in selecting the most plausible non-linear model class from the list of constructed model classes using measured vibration data under different levels of impact loads,

to identify the modified Voce law as the “best” elastoplastic model for capturing the stress–strain behavior of ballast for impact load larger than 1.5 kN, and

to present a non-destructive vibration-based methodology for accurate prediction of the location and severity of ballast damage underneath the target sleepers.

In the future research, investigations on the applicability of Bayesian model updating for ballasted tracks using the modified Voce model will be carried out, using measurement data from ballasted railway tracks under operational conditions. In the field tests, train-induced track vibration data shall be utilized instead of impact hammer test. This phase of research will help in the development of a real-time in-service long-term monitoring system for ballasted railway tracks.

Footnotes

Acknowledgments

The help from Dr Jia-Hua Yang, Dr Jun Hu, Mr Yi-Ming Liu, and Dr Man-Tat Wong in the experiment is highly appreciated.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The work described in this article was fully supported by two grants from the Research Grants Council of the Hong Kong Special Administrative Region, China [Project No. CityU 11242716 (GRF 9042336) and R5020-18 (RIF 8799008)] and the National Natural Science Foundation of China (Project No. 51708242).

ORCID iDs

Mujib Olamide Adeagbo

Qin Hu

References

Moormann

Lehn

Aschrafi

, et al. Numerical analysis of a track-substructure system considering different modelling approaches. Proc Third Int Conf Railw Technol Res Dev Maint. 2016; 143: 1–10. DOI: 10.4203/ccp.110.38.

Lekarp

Isacsson

Dawson

. State of the art. i: resilient response of unbound aggregates. J Transportation Eng 2000; 126: 66–75.

Paixão

Varandas

Fortunato

, et al. Non-linear behaviour of geomaterials in railway tracks under different loading conditions. Proced Eng 2016; 143: 1128–1135. DOI: 10.1016/j.proeng.2016.06.147.

Wei Fan

Pizhong Qiao

. Vibration-based damage identification methods: a review and comparative study. Struct Health Monit 2011; 10: 83–111. DOI: 10.1177/1475921710365419.

Sun

Indraratna

Nimbalkar

. An elasto-plastic method for analysing the deformation of the railway ballast. Proced Eng 2016; 143: 954–960. DOI: 10.1016/j.proeng.2016.06.082.

Varandas

Paixão

Fortunato

, et al. A numerical study on the stress changes in the ballast due to train passages. Proced Eng 2016; 143: 1169–1176. DOI: 10.1016/j.proeng.2016.06.127.

Monismith

Seed

Mitry

, et al. Predictions of pavement deflections from laboratory tests. In: Proc. Second Int. Conf. Struct. Des, Ann Arbor, Michigan, 7-11 August, Asph. Pavements, 1967, 53–88.

Brown

Pell

. An experimental investigation of the stresses, strains and deflections in a layered pavement structure subjected to dynamic loads. Proc 2nd Int Conf Struct Des Asph Pavements 1967, pp. 487–504.

Uzan

. Characterization of granular material. Transp Res Rec 1985; 1022: 52–59.

10.

Dawson

Thom

Paute

. Mechanical characteristics of unbound granular materials as a function of condition. In: Correia

(ed). Flexible Pavements, Proceedings Euroflex Symposium. Rotterdam, The Netherlands: Balkema; 1993, pp. 35–44.

11.

Thompson

. The effects of local preload on the foundation stiffness and vertical vibration of railway track. J Sound Vibration 1999; 219: 881–904. DOI: 10.1006/jsvi.1998.1939.

12.

Indraratna

Salim

Rujikiatkamjorn

. Advanced Rail Geotechnology - Ballasted Track. Florida, USA: CRC Press/Balkema, 2011.

13.

Jensen

Esse

Araya

, et al. Implementation of an adaptive meta-model for Bayesian finite element model updating in time domain. Reliability Eng Syst Saf 2017; 160: 174–190. DOI: 10.1016/j.ress.2016.12.005.

14.

Marseguerra

Zio

. Monte Carlo simulation for model-based fault diagnosis in dynamic systems. Reliability Eng Syst Saf 2009; 94: 180–186. DOI: 10.1016/j.ress.2008.02.013.

15.

Lackenby

Indraratna

McDowell

, et al. Effect of confining pressure on ballast degradation and deformation under cyclic triaxial loading. Géotechnique 2007; 57: 527–536. DOI: 10.1680/geot.2007.57.6.527.

16.

Indraratna

Lackenby

Christie

. Effect of confining pressure on the degradation of ballast under cyclic loading. Géotechnique 2005; 55: 325–328. DOI: 10.1680/geot.55.4.325.65490.

17.

Sun

Indraratna

Nimbalkar

. Effect of cyclic loading frequency on the permanent deformation and degradation of railway ballast. Géotechnique 2014; 64: 746–751. DOI: 10.1680/geot.14.t.015.

18.

Sun

Indraratna

Nimbalkar

. Deformation and degradation mechanisms of railway ballast under high frequency cyclic loading. J Geotechnical Geoenvironmental Eng 2016; 142: 04015056. DOI: 10.1061/(ASCE)GT.1943-5606.0001375.

19.

Kaewunruen

Remennikov

. Effect of improper ballast packing/tamping on dynamic behaviors of on-track railway concrete sleeper. Int J Struct Stab Dyn 2007; 07: 167–177.

20.

Ngamkhanong

Kaewunruen

Baniotopoulos

. A review on modelling and monitoring of railway ballast. Struct Monit Maint 2017; 4: 195–220. DOI: 10.12989/smm.2017.4.3.195.

21.

Lam

Wong

. The Bayesian methodology for the detection of railway ballast damage under a concrete sleeper. Eng Structures 2014; 81: 289–301. DOI: 10.1016/j.engstruct.2014.08.035.

22.

Adeagbo

Lam

. A feasibility study on damage detection of nonlinear railway ballast by measured vibration of in-situ sleeper. In: Wang

JCM

Kitipornchai

(eds). ACMSM25 Lecture notes civil engineering. Singapore: Springer; 2020, Vol. 37, 351–360.

23.

Lam

Yang

, et al. Railway ballast damage detection by Markov chain Monte Carlo-based Bayesian method. Struct Health Monit 2018; 17: 706–724. DOI: 10.1177/1475921717717106.

24.

Akeju

Senetakis

Wang

. Bayesian parameter identification and model selection for normalized modulus reduction curves of Soils. J Earthquake Eng 2019; 23: 305–333. DOI: 10.1080/13632469.2017.1323051.

25.

Wang

Aladejare

. Selection of site-specific regression model for characterization of uniaxial compressive strength of rock. Int J Rock Mech Mining Sci 2015; 75: 73–81. DOI: 10.1016/j.ijrmms.2015.01.008.

26.

Chiachío

Prescott

, et al. A knowledge-based prognostics framework for railway track geometry degradation. Reliability Eng Syst Saf 2019; 181: 127–141. DOI: 10.1016/j.ress.2018.07.004.

27.

Andrade

Teixeira

. Statistical modelling of railway track geometry degradation using Hierarchical Bayesian models. Reliability Eng Syst Saf 2015; 142: 169–183. DOI: 10.1016/j.ress.2015.05.009.

28.

C-T

. Bayesian model updating approach for experimental identification of damage in beams using guided waves. Struct Health Monit 2014; 13: 359–373. DOI: 10.1177/1475921714532990.

29.

Adeagbo

Lam

H-F

. A Bayesian methodology for detection of railway ballast damage using the modified Ludwik nonlinear model. Eng Struct 2021; 236: 112047.

30.

Lam

H-F

Yang

S-K

. Bayesian model updating of a coupled-slab system using field test data utilizing an enhanced Markov chain Monte Carlo simulation algorithm. Eng Structures 2015; 102: 144–155. DOI: 10.1016/j.engstruct.2015.08.005.

31.

Winkler

. The Theory of Elasticity and Strength. Prague: Dominicus, 1867.

32.

Grassie

Cox

. The dynamic response of railway track with unsupported sleepers. Proc Inst Mech Eng D: Transport Eng 1985; 199: 123–136. DOI: 10.1243/pime_proc_1985_199_149_01.

33.

Hetenyi

. Beams on Elastic Foundation. Ann Arbor: The University of Michigan Press, 1979.

34.

Lam

H-F

Adeagbo

Yang

Y-B

. Time-domain markov chain monte carlo-based bayesian damage detection of ballasted tracks using nonlinear ballast stiffness model. Struct Heal Monit 2020. DOI: 10.1177/1475921720966950.

35.

Salim

Indraratna

. A new elastoplastic constitutive model for coarse granular aggregates incorporating particle breakage. Can Geotechnical J 2004; 41: 657–671. DOI: 10.1139/T04-025.

36.

Dawson

Correia

Jaure

, et al. Modelling resilient and permanent deflections in granular and soil pavement layers. Proc 4th Int Conf Bear Capacit Roads Airfields 1994; 2: 847–861.

37.

Wang

Markine

. Modelling of the long-term behaviour of transition zones: prediction of track settlement. Eng Structures 2018; 156: 294–304. DOI: 10.1016/j.engstruct.2017.11.038.

38.

Duncan

Chang

C-Y

. Nonlinear analysis of stress and strain in soils. J Soil Mech Foundations Division 1970; 96: 1629–1653.

39.

Adeagbo

Lam

. Bayesian modeling of axial strain-hardening behaviour in railway ballast. In: Lecture Notes in Civil Engineering. Singapore: Springer, 2021, pp. 1089–1098. DOI: 10.1007/978-981-15-8079-6_102

40.

Ludwik

. Elemente der Technologischen Mechanik. Berlin: Julius Springer, 1909.

41.

Swift

. Plastic instability under plane stress. J Mech Phys Sol 1952; 1: 1–18.

42.

Prager

. A new mathematical theory of plasticity. Rev Fac Sci Univ Istanbul 1941; 5: 215–226.

43.

Voce

. The relationship between stress and strain for homogeneous deformations. J Inst Met 1948; 74: 537–562.

44.

Chakrabarty

. Theory of Plasticity. 3rd ed. Oxford, UK: Butterworth-Heinemann, 2006.

45.

Augustin

Gudehus

Huber

, et al. Numerical model and laboratory tests on settlement of ballast track. In: System Dynamics and Long-Term Behaviour of Railway Vehicles, Track and Subgrade. Berlin, Germany: Springer Berlin Heidelberg, 2003, 317–336. DOI: 10.1007/978-3-540-45476-2_19.

46.

Lazarević

Vuĉković

Popović

. Assessment of sleeper support conditions using micro-tremor analysis. Proc Inst Mech Eng F J Rail Rapid Transit 2016; 230: 1828–1841. DOI: 10.1177/0954409715615629.

47.

Zhang

Xiao

Wen

, et al. Effect of unsupported sleepers on wheel/rail normal load. Soil Dyn Earthquake Eng 2008; 28: 662–673. DOI: 10.1016/j.soildyn.2007.08.006.

48.

Zhu

Ahmed

AKW

Rakheja

, et al. Development of a vehicle-track model assembly and numerical method for simulation of wheel-rail dynamic interaction due to unsupported sleepers. Vehicle Syst Dyn 2010; 48: 1535–1552. DOI: 10.1080/00423110903540751.

49.

Lazarević

Vučković

. Assessment of sleeper stability in ballast bed using micro-tremor sampling method. In: Adv. Intell. Syst. Comput. Berlin, Germany: Springer-Verlag, 2020, 290–299. DOI: 10.1007/978-3-030-19756-8_26.

50.

Cheung

Beck

. Calculation of posterior probabilities for Bayesian model class assessment and averaging from posterior samples based on dynamic system data. Computer-Aided Civil Infrastructure Eng 2010; 25: 304–321. DOI: 10.1111/j.1467-8667.2009.00642.x.

51.

Katafygiotis

Beck

. Updating models and their uncertainties. II: model identifiability. J Eng Mech 1998; 124: 463–467. DOI: 10.1061/(ASCE)0733-9399.

52.

Beck

Katafygiotis

. Updating models and their uncertainties. I: Bayesian Statistical Framework. J Eng Mech 1998; 124: 455–461. DOI: 10.1061/(asce)0733-9399.

53.

Yuen

K-V

. Recent developments of Bayesian model class selection and applications in civil engineering. Struct Saf 2010; 32: 338–346. DOI: 10.1016/j.strusafe.2010.03.011.

54.

Asaadi

Heyns

. A computational framework for Bayesian inference in plasticity models characterisation. Comp Methods Appl Mech Eng 2017; 321: 455–481. DOI: 10.1016/j.cma.2017.04.017.

55.

Fatahi

Moradi

. Multiple crack identification in frame structures using a hybrid Bayesian model class selection and swarm-based optimization methods. Struct Health Monit 2018; 17: 39–58. DOI: 10.1177/1475921716683360.

56.

Beck

. A new adaptive importance sampling scheme for reliability calculations. Struct Saf 1999; 21: 135–158.

57.

Hastings

. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 1970; 57: 97–109. Available at: http://www.jstor.org/stable/2334940.

58.

Metropolis

Rosenbluth

, et al. Equation of state calculations by fast computing machines. J Chem Phys 1953; 21: 1087–1092. DOI: 10.1063/1.1699114.

59.

Muto

Beck

. Bayesian updating and model class selection for hysteretic structural models using stochastic simulation. J Vib Control 2008; 14: 7–34. DOI: 10.1177/1077546307079400.

60.

Kullback

Leibler

. On information and sufficiency. Ann Math Stat 2007; 22: 79–86. DOI: 10.1214/aoms/1177729694.

61.

Cover

Thomas

. Elements of Information Theory. 2nd Editio. New York: John Wiley & Sons. 2005. DOI: 10.1002/047174882X.

62.

Lam

H-F

Alabi

Yang

J-H

. Identification of rail-sleeper-ballast system through time-domain Markov chain Monte Carlo-based Bayesian approach. Eng Structures 2017; 140: 421–436. DOI: 10.1016/j.engstruct.2017.03.001.

63.

Kaewunruen

Remennikov

. Sensitivity analysis of free vibration characteristics of an in situ railway concrete sleeper to variations of rail pad parameters. J Sound Vibration 2006; 298: 453–461. DOI: 10.1016/j.jsv.2006.05.034.

64.

Kaewunruen

Remennikov

. Field trials for dynamic characteristics of railway track and its components using impact excitation technique. NDT E Int 2007; 40: 510–519. DOI: 10.1016/j.ndteint.2007.03.004.