Formulation
This page summarizes the mathematical and numerical methods used by STEM to simulate railway-induced ground vibrations. STEM performs the numerical analysis of the train-track-soil system using the finite element method (FEM) and is powered by Kratos Multiphysics.
Governing equation
STEM solves the dynamic equilibrium equation, following a Total Lagrangian formulation with small strains. The governing finite element equation to be solved is:
\[\mathbf{M}\mathbf{a} + \mathbf{C}\mathbf{v} + \mathbf{K}\mathbf{u} = \mathbf{F_{ext}}\left( t \right)\]
where \(\mathbf{M}\) is the mass matrix, \(\mathbf{C}\) is the damping matrix, and \(\mathbf{K}\) is the stiffness matrix of the entire system, \(\mathbf{F_{ext}}\) denotes the vector of the external forces and \(\mathbf{a}\), \(\mathbf{v}\), \(\mathbf{u}\) are, respectively, the acceleration, the velocity and the displacement in the nodes.
In STEM it is also possible to perform quasi-static analyses, in which inertial and damping effects are neglected. In this case, the governing equation is reduced to:
\[\mathbf{K}\mathbf{u} = \mathbf{F_{ext}}\left( t \right).\]
Spatial discretisation
The spatial discretisation of the problem is performed using the finite element method (FEM). The geometry of the problem is built using gmsh and the mesh is generated and exported to Kratos mdpa format.
STEM supports several types of elements:
Beam elements (for 2D and 3D problems):
Euler-Bernoulli beam elements (2 nodes).
2D plane strain elements:
First order triangle elements (3 nodes per element)
Second order triangle elements (6 nodes per element)
3D solid elements:
First order tetrahedral elements (4 nodes per element)
Second order tetrahedral elements (10 nodes per element)
Time integration and solution strategies
The system of equations is solved in the time domain. A range of solution strategies are supported, including:
Implicit and explicit time integration schemes.
Nonlinear solution strategies (Newton-Raphson, line search, arc length).
For linear systems, it is recommended to use a custom explicit integration scheme based on the Newmark method, LinearNewtonRaphsonStrategy. This integration scheme takes advantage of the linear nature of the problem, and it is formulated explicitly. Following the Newmark method, the acceleration and velocity at time step \(n+1\) are calculated as:
where \(\Delta t\) is the time step for the numerical integration, the subscript \(\sqcup_{t + \Delta t}\) denotes kinematic variables at time \(t + \Delta t\), and \(\gamma\) and \(\beta\) are the parameters that control the integration accuracy and stability (the default values are \(\gamma = 0.5\) and \(\beta=0.25\)).
Substituting the Newmark equations into the governing equation, and taking advantage of the linearity of the problem, the system of equations yields:
where \(\mathbf{\widetilde{K}}\), \(\mathbf{\widetilde{a}}\) and \(\mathbf{\widetilde{v}}\) are the effective stiffness, acceleration and velocity, respectively, defined as:
The main advantage of using the custom Newmark explicit integration scheme is that it only requires once the assembly and factorisation of the entire system of matrices, making it computationally efficient for linear problems.
In STEM it is also possible to perform quasi-static analyses. In the quasi-static case, the governing equation is reduced to:
\[\mathbf{K}\mathbf{u} = \mathbf{F_{ext}}\left( t \right)\]
In this case, the system of equations is solved at each time step.
To solve the linear system of equations of type:
\[\mathbf{A}\mathbf{x} = \mathbf{b},\]
STEM has several linear solvers available, including:
LU decomposition
Conjugate Gradient
Algebraic Multigrid Iterative
The LU decomposition is a direct solver, which is suitable for small problems. The Conjugate Gradient and Algebraic Multigrid Iterative solvers are iterative solvers, which are suitable for large problems. The iterative solvers can use a preconditioner to improve the convergence rate, and they are recommended for large problems. In STEM the default preconditioner is the Jacobi preconditioner, which is a simple diagonal preconditioner that can be used for symmetric positive definite systems.
Boundary conditions
STEM support different types of boundary conditions, including:
Dirichlet boundary conditions (prescribed displacement or rotation)
Neumann boundary conditions (force or traction boundary conditions)
Absorbing boundaries [3]
Loads
STEM supports different types of loads, including:
Point, line, surface and moving loads
All the loads can be defined as time-dependent, and they can be applied to structural elements or soil.
For STEM it is important to be able to model the train-track interaction, which is the main source of the dynamic loading in the problem of railway-induced vibrations. To model the train-track interaction force STEM uses point moving loads, in combination with the UVEC (see Train-track model).
Constitutive models
STEM supports linear elastic material models. Non-linear material models can be implemented by the user and integrated into STEM through the UMAT interface (see Interface definitions).
Linear elastic material
STEM supports linear elastic material models. In STEM a material is defined by specifying the material formulation and the linear elastic material parameters.
Currently STEM supports the OnePhaseSoil material formulation. This formulation is used to model the soil as a single phase material, and it is defined by the following parameters:
Density of the solid phase \(\rho_s\)
Porosity \(n\)
Is_drained (boolean flag to indicate if the soil is drained or not)
The linear elastic material requires the following parameters for the OnePhaseSoil formulation:
Young’s modulus \(E\)
Poisson’s ratio \(\nu\)
In STEM the soil should be defined as OnePhaseSoil and the solution should be ‘drained’. This means that the soil is modelled as a single phase material and that the pore water pressure is not considered in the analysis. This is an assumption that is commonly used in the analysis of dynamic problems, such as railway-induced vibrations, and it is justified by the fact that the time scale of the problem is much shorter than the time scale of the pore water pressure dissipation.
Typically, the density of the bulk material is available from site investigations. The bulk density consists of the density of the solid phase and the density of the fluid phase, and it can be calculated as:
where \(\rho_f\) is the density of the fluid phase (water). Because in STEM the soil is modelled as a drained single phase material this equation simplifies to:
If the user has the bulk density and wants to use it as input, the porosity can be set to zero and the density of the solid phase can be set as the bulk density.
To model soil layers below ground water table, since the soil is modelled as a drained single phase material, the Poisson ratio of the soil layer should be set to 0.495. This means that the material is nearly incompressible, i.e. cannot experience volumetric deformation and this accurately mimics the saturated behaviour of soil, subjected to dynamic loading, in the short term (before pore water pressure dissipation occurs).
Train-track model
The train and the train-track interaction model can be modelled by a user-defined vehicle model (UVEC) as described in Interface definitions. STEM provides a default train and train-track interaction model.
Vehicle model
The train is modelled as a 2D mass-spring-damper system, with 10 degrees of freedom, that consists of a car body, bogies, suspension systems and wheels [5]. This model takes advantage of the system’s symmetry by simulating only half of the train, with symmetry boundary conditions applied to the central plane. The figure shows the schematic of the train model.
STEM train model
The railway track consists of rail, railpad and sleeper. The rail is modelled as an Euler-Bernoulli beam element. The railpad is modelled as a spring-damper element, which provides the stiffness and damping of the railpad. The sleeper can either be modelled as a concentrated mass, or a volume element, depending on the level of detail required by the user.
Interaction model
The vehicle and track systems are coupled through contact forces that occur at the interface between the wheels and the rails. These forces are determined using the non-linear Hertzian contact theory for metals [4]:
where \(k_\text{c}\) is the Hertzian stiffness coefficient and \(\delta_j\) is the indentation for wheel \(j\). The indentation \(\delta_j\) is calculated as:
where \(u_{\text{w},j}\) is the displacement of wheel \(j\), \(u_{\text{r},j}\) is the displacement of the rail at the position of wheel \(j\) and \(z_j\) is the irregularity of the rail at the position of wheel \(j\).
Irregularities model
The vertical irregularities presented in the rail are expressed following the methodology presented in [6]. The power spectral density (PSD) of the rail irregularities \(S(\Omega)\) is used to produce the samples of the irregularities.
where \(\Omega\) the wavenumber, \(A_v\) the rail irregularity parameter and \(\Omega_c\) is the critical wavenumber. These two parameters define the quality of the railway track [2].
The sample of rail irregularities can be produced by inverse Fourier transform shown as follows:
where \(\omega_n\), is a circular frequency within the interval in which the PSD function is defined, \(\theta_n\), a random phase angle uniformly distributed from 0 to \(2\pi\), and \(\Delta \omega\) is defined as the total number of frequency increments \(N\) in the range of the circular frequency.
The parameter \(A_v\) can be estimated based on the track quality [2].
Line grade |
Av Value (m² rad / m) |
|---|---|
1 (very poor) |
1.2107e-4 |
2 |
1.0181e-4 |
3 |
0.6816e-4 |
4 |
0.5376e-4 |
5 |
0.2095e-4 |
6 (very good) |
0.0339e-4 |
Rail joint model
The rail joint is modelled as a hinge along a beam element, with a rotational stiffness defined based on the properties of the rail and the joint. The rail joint might have a local geometric irregularity, commonly referred to as a dipped joint, following the approach proposed by [1]. The dipped joint represents the permanent plastic deformation that develops at rail ends due to repeated wheel-rail impact loads at insulated joints.
The dip is introduced as a prescribed vertical rail surface irregularity with a smooth cosine shape, defined by its maximum depth \(d\) and total length \(l\). The vertical irregularity profile \(r(x)\) along the rail longitudinal coordinate \(x\) is given by:
where \(x = 0\) corresponds to the start of the dipped region. This formulation ensures a continuous and smooth rail profile with zero slope at the beginning and end of the joint, avoiding artificial high-frequency excitation.
A schematic of the dipped joint is shown in the figure below.
Schematic of the dipped joint model