Efficient numerical approach to unbounded systems subjected to a moving load

The present paper solves numerically the problem of vibrations of infinite structures under a moving load. A velocity formulation of the space–time finite element method was applied. In the case of simplex shaped space–time finite elements, the ‘steady state’ dynamic behaviour of the system was obtained. A properly performed discretization allowed of propagating information in a given direction at a limited velocity. The solutions were obtained under the assumption that the deformation is quasi-stationary, i.e., stationary in the coordinate system that moves with the load. The unbounded Timoshenko beam subjected to a distributed moving load was used as a test example. The dynamical system is placed on an elastic foundation. The matrices describing an infinite dynamical system subjected to a moving load are derived and the stability of the numerical scheme is analysed. The numerical results are compared with the analytical solutions in the literature and the classical numerical method.


Introduction
Simulations of a moving load travelling along a periodic structures are extensively performed in modern transportation, robotics, industry, etc. The velocities of moving objects, such as a train, can grow rapidly. As a result, we can observe a new phenomenon that has not been properly investigated. Particularly important are wave problems when the train B. Dyniewicz (B) Institute of Fundamental Technological Research, Polish Academy of Sciences, Pawińskiego 5b, 02-106 Warsaw, Poland e-mail: bdynie@ippt.gov.pl velocity gets close to the characteristic wave velocities generated in a railroad track. This results in an increasing stress in the rail, which can cause a buckling of the track [13,17,18]. Waves phenomena are also important in catenarypantograph dynamic interaction problems [2,11,15]. Broad overview of the numerical methods is presented in the paper [25]. It discuss in detail the ways of modeling of pantographs and overhead contact lines. The results of numerical simulations and experimental studies are also presented.
In literature we can find different models of generation and propagation of vibration in an unbounded structure. In the case of simple structures for the track or catenary and the moving loads, analytical solutions of the problem are known. The solution can be found in a moving coordinate system, and then the complex Fourier transform must be applied. Finally, we obtain the structural response in the form of standing waves [10,21] or travelling waves [1,8,9]. Broad review of analytical solutions of the load moving along an infinite beam on an elastic foundation was presented in [16]. In more complex models numerical methods can be used. Numerical description of the loads travelling along a string, beams or plates have been widely presented in the literature [3,14,27]. This concerns both the classical FEM with spatial discretization as well as space-time discretization. Recent paper of numeric calculations concern mainly inertial loads. In [12] the multiplex shaped space-time finite element approach applied to general description of a moving inertial load.
In numerical simulations an infinite structure is replaced by a finite one. In order to minimize the influence of wave reflections from the ends of the structure, the numerical methods require the simulation of large sections of a track or a catenary. This approach results in a task with a large number of degrees of freedom and this significantly increases the computational time. Different methods of eliminating adverse phenomena have been used [19,20,28]. Artificial boundary Fig. 1 Information flow in methods with full inertia matrix conditions applied to the edges of the structure have been designed to suppress the reflected waves. These methods are acceptable in problems where we can predict the shape of the reflected wave. Presented in mainstream literature numerical simulations of dynamic systems are based on the explicit and implicit methods of the integration at time. Due to the form of the inertia matrix various methods of calculation have important feature in wave problems. It is the way of the information flow in space-time mesh during the step-by-step calculations. Propagation of information is a kind of wave unrelated to physical system but associated with the numerical system. According to [22] the schematic information flow in the classical schemes of the numerical integration is presented in Fig. 1.When the matrix of inertia is consistent or band, the information is transferred from node to node in the successive time step with an infinitely high speed. This mechanism of information flow is presented on Fig. 1. A disturbance F at one point of the structure is influences to every node of the system. An impulse applied to a node of the mesh results in non-zero values in the entire space-time mesh in successive time steps. This is not a physical feature in the case of wave problems. The infinite information speed in space-time mesh occurs both in explicit and implicit methods. This feature is characteristic for multiplex shaped space-time finite element approaches. If the inertia of the structure can be provided as a diagonal matrix of inertia, the calculation scheme is simplified. We obtain a system of the separated algebraic equations. According to Fig. 2 in each successive time step an initial single impulse of the external force F propagates to the neighboring nodes in the spatial mesh at a rate of one node in one step. Therefore, the information propagates at a finite speed in both directions. Thus the type of the spacetime mesh and the time integration scheme influences results. The limited velocity of the information flow is characteristic for hyperbolic problems while the infinite speed is typical for parabolic problems. This parabolic-like behaviour of the wave solutions is considered as a serious disadvantage, since wave reflections from boundaries of the mesh influence sig- In order to avoid these disadvantages, we will use the simplex shaped space-time finite element method. This special type of space-time finite elements has an interesting property. It allows the flow of numerical information in one direction at a limited speed. It results in triangular matrices of resulting systems of algebraic equations. Examples of onedimensional simplex elements are depicted in Fig. 3. In the case (a) the information flows from left to right while in the case of (b)-oppositely. In the further sections of the paper we consider the case (a). The property of one-way information flow can be successfully used in moving load problems. If we use triangular elements with its slope sides directed to the moving load, the specific selection of the space-time element size (b/ h ratio) allows us to coincide the characteristic line of the speed of numerical information wave with the speed of external moving load in the structure. Therefore a properly performed discretization allows the propagation of information in a given direction at a reduced velocity (Fig. 4). Nodes of the mesh before the front of the moving load are not disturbed by non-physical informations in the successive time steps, as in conventional numerical methods. Thus information, and consequently the reflected waves following the subjected nodes, never reach the moving load. The area of disturbances caused by the load allows us to limit the observations to a region located in a closed neighbourhood of the load position (Fig. 5). Isolation of the subsystem of finite elements allows us to bind the reference system and a moving load. Therefore, if in the physical system a stationarity of a phenomenon occurs, we can simply assume it in our considerations. The steadystate problem can not be solve by conventional numerical methods without artefacts.
In this paper, we focus on numerical approach with the propagation of information in a given direction at a reduced velocity. Considered solution of a simplified model of some systems that can be directly related to engineering problems. An unbounded Timoshenko beam resting on elastic supports was taken as a reference problem and solved numerically. The stationary solutions are considered. Elementary matrices describing the beam, the boundary conditions, and the distributed moving load are derived. Then a comparison with the literature examples and the discussion of numerical and analytical solutions are given. The next section includes notes on the stability of the numerical approach. Finally, some concluding remarks will summarize all these considerations.

Formulation of the problem and analytical solution
Let us consider an unbounded Timoshenko beam resting on an elastic foundation (Fig. 6). The equations of motion are written in following form where dot and prime denote the differentiation with respect to t and with respect to x. The coordinate x is fixed to the beam, t is time variable, and c is the elastic coefficient of the foundation. A uniformly distributed load travelling along the beam with a constant speed of v 0 is described by the Heaviside step function H . The remaining parameters of a Timoshenko beam with a constant cross-sectional area A, moment of inertia I , mass density ρ, Young's modulus E, shear modulus G, and shear coefficient k are known. The motion of the structure is completely described by the set of Eq. (1). This problem was solved analytically in [8] and examples given there were used to verify our numerical solutions. The analytical solution describes the stationary case. According to moving coordinate system fixed with the front of the load the equations of motion (1) can be written in the ordinary differential form with respect to spatial variable. Let us introduce the following notation: .
v * 1 and v * 2 are the bending wave speed and shear wave speed in the Timoshenko beam, respectively. The set of Eq. (1) is reduced to a fourth order differential equation in terms of displacements or rotations. In the case of rotation we obtain in which Contrary to the displacement case, the form of (3) enables the Fourier transformation and simple return to the spatial domain. According to (1) (3) can assume hyperbolic or parabolic type of the differential equation. So, these characteristic dimensionless velocities V 1 and V 2 are important for qualitative properties of the solutions. Two cases: They result in particular forms of the solution of (1). Distinction between these two types of the solutions can be important for calculation of engineering structures such as railway tracks. When the simplified model is taken into account, the partial weight of the ballast is associated with the mass of the rail. Properly designed classic railway track fulfills the hyperbolic type of the differential equation. The increase of the ballast weight, that can occur in the case of hydration, the transverse wave speed increases. The type of the problem starts to be parabolic. Both cases of analytical solution exhibit different responses of the system to the same velocities of the moving load. In the space-time finite element solution the split of the differential equation of Timoshenko beam to hyperbolic or parabolic type has no direct relevance. However, both cases were considered since numerical results are compared with analytical solutions. More detailed discussion of the presented analytical solutions can be found in [8].

Simplex shaped space-time approach
The conception of symplectic elements was proposed for the first time in [23]. In [4], a displacement variant of the triangular space-time finite element method was successfully applied to beams. A velocity variant was presented in [3]. Simplex shaped elements have been used in contact problems [5] and tested in the context of advection-diffusion Eq. [7].
This continuous Galerkin method discretizes the spatial variables and time variable simultaneously. Therefore, we can postulate a balance of energy over the interval of time, not only at some moments. In the formulation of the method, we integrate the physical quantities analytically in the time interval rather than numerically, as in the classical finite element method. The numerical method presented will be used to solve the evolutionary processes of the vibrations. In this paper, the velocity formulation of the space-time finite element method is used. The equations of motion are discretized both in space and time. This means that in the present numerical scheme the velocity is distributed in a finite space-time element according to interpolation functions and nodal velocities in two adjacent time layers used as parameters. The analytic form of the velocity function in space and time allows the integration and differentiation with respect to these variables. Functions of displacement and acceleration determined in this way, however, still depend on the nodal velocities. Thus, as the result we have both inertia and stiffness matrices multiplied by velocities. The integration of the velocity function with respect to time results in the displacement function and contains the term with initial displacements in time layer. Initial displacements in each time step of time stepping scheme in the final equation of force equilibrium express nodal forces at the beginning of time interval. According to the equations of motion (1) calculated acceleration and displacement functions allow the analytical determination of the energy of the system. The energy of the external forces is derived from the right hand side of the Eq. (1). Classical energy minimization and assembly of the global system leads to the following matrix solution scheme Here, i and i + 1 denote the known and calculated state, respectively. The problem is reduced to the numerical solution of the system of algebraic Eq. (5). The vector v contains the velocity of the nodal displacements and angles of rotation, and the vector w contains the nodal displacements and the angles of rotation. M g , K g , and E g are the global matrices of inertia, stiffness, and nodal forces. F g is the global vector of external forces. The global matrices are assembled from the local matrices M, K, and E, which will be derived further. These elemental matrices are merged in appropriate locations of the global matrices, based on the topology of the mesh. Similarly, the global load vector is assembled from elemental force vectors F. This vector can assume zero values or values describing the distributed external load in finite elements. According to the current position of the moving load, the vector F g varies in each time step. The vector F g has zero components at the beginning and in each successive time step gains new non-zero values contributed by vectors F. The current vector of displacements and rotation angles can be computed using The partition of the space-time area into elements of simplex shape allows us to obtain the stationary solution numerically. In order to demonstrate the properties of the method, let us consider a single simplex finite element as depicted in Fig. 3a.
In this case, two space-time subdomains We assume a linear distribution of the velocity displacements v =ẇ and the velocity rotation angles ψ =θ inside the triangular element: v =ẇ(x, t) = a 1 x + a 2 t + a 3 , ψ =θ(x, t) = a 1 x + a 2 t + a 3 .
In the finite element shown in Fig. 3a, the interpolation of nodal velocities can be written The same shape functions were assumed for the velocities of the nodal rotations. Below, calculations for part A of the space-time finite element are presented. We must emphasize here that in the case of part B, the procedure will be totally analogous. Nodal displacements and rotation angles can be written as the integrals of the velocities: In order to determine the virtual energy of the problem, we multiply the equations of motion (1), by the virtual functions v * A and ψ * A and integrate the resulting power over the spacetime domain A . We consider a distributed load and assume that the load is moved in each time step from one finite element to another, without intermediate steps. The Heaviside step function describing the external load is replaced by a constant function. Thus, virtual energy of the part A of the space-time is given by the following form The proper choice of virtual functions is a fundamental question of the space-time approach. Various functions in time result in solution schemes of different accuracy and stability. Virtual function review in the case of multiplex shaped space-time finite element approach is presented in [6]. In this case the virtual linear shape function is assumed. After integration by parts, we have The classical energy minimization of (11) leads us to the matrices M A , K A , and E A for part A of the simplex element. We do the same for part B. After aggregating the matrices of parts A and B, we get the local matrices of inertia, stiffness, and nodal forces: The stiffness matrix K (13) and nodal forces matrix E (14) can be split into three parts: a part related to bending, a part related to foundation, and a part related to shear. Assuming the linear shape function (7) in virtual energy (11), we indicate that the bending strain component is constant and the shear strain component varies linearly. In this case the exact integration of the shear strain components in virtual energy (11) results in the element which is too stiff. This over-stiffness of the element is known and is called locking. The shear components in relation with the bending component are high. A reduced integration of the shear part is a classical literature remedy. Bending terms are integrated exactly while shear terms are integrated with only one point of the Gauss quadrature. More details can be found in [24,26], e. g. As mentioned previously, the matrices (12), (13), and (14) were assembled into the global matrices M g , K g , and E g . As a result of the space-time integration of the right hand side of (11), the local vector of external force takes the following form Finally, according to the current position of the travelling distributed load, we fill up the global external load vector F g with the local vectors (15), based on the subjected nodes. These nodes alter from step to step.

Numerical results
As mentioned in Sect. 2 the numerical results was compared with analytical solution from [8]. We present two examples that depend on the dimensionless shear wave speed V 1 and dimensionless bending wave speed V 2 . Each example consideres three characteristic ranges of the dimensionless velocity V of the moving load with qualitatively different solutions. The resulting dimensionless displacements W are presented in a dimensionless coordinate system, X, fixed to the beam. According to the assumptions the stationary case is considered. The shape of the Timoshenko beam can be observed in the moving frames of reference, i. e. the moving coordinate system. The number of steps of calculation relates to the velocity of the moving load. Properly selected number of time steps allows us to plot accurately the curve of displacement wave of the deformed Timoshenko beam. This number of steps is required for the formation of displacement wave in a region located in the neighbourhood of the front of the load. In our case for V 2 = 4 we used 40 steps while for V 2 = 55 we used 16 steps.

The case of
The characteristic form of the displacement wave at dimensionless velocity V 2 = 4 was captured and depicted in Fig. 7. This velocity of the moving load is in the range described by the condition V 2 < V 2 1 . The solution tends monotonically to the asymptotes W = 0 and W = 1. In the range V 2 1 < V 2 < V 2 2 , the dimensionless speed V 2 = 8 was used. The solution is illustrated by the curve in Fig. 8. The solution vanishes monotonically before the load front. If the load is moving faster than the bending wave in the beam, i.e., V 2 > V 2 2 , the characteristic form of the displacement wave at dimensionless velocity V 2 = 55 was presented in Fig. 9. The displacements are equal to zero before the load front.

The case of
Within the range V 2 < V 2 1 , the solution at dimensionless speed V 2 = 4.5 is presented in Fig. 10. Some differences between the analytical and numerical solutions are can be seen. In this range the solution substantially changes its quantitative feature. In the analytical case, the solution consists then of two noticeable waves. One wave with a low amplitude and a short wavelength is placed in front of the load. The second wave is characteristic of higher amplitude and longer wavelength. It appears behind the load front. Prepared numerical solution revealed deficiencies of plots of analytical solution in the literature. The wave with a constant amplitude behind and ahead the load front occurs but with a significantly lower amplitude. In fact these waves are not visible in practice. Perhaps the author of the figure presenting the analytical solution in [8] was going to highlight qualitative features of his solution. The solutions for the range V 2 1 < V 2 < V 2 2 and V 2 > V 2 2 exhibit similar features to the solutions in subsection 4.1. The characteristic form of the displacement wave at dimensionless speeds V 2 = 6 and V 2 = 15 are depicted in Figs. 11 and 12, respectively.
The resulting numerical solutions were compared with the classical implicit Newmark method. Characteristic forms of the displacement wave are depicted in Figs. 13 and 14. Figure  13 corresponds to the hyperbolic case 4.1, whereas Fig. 14 corresponds to the parabolic case 4.2. A comparison between the simplex shaped space-time finite element approach and the classical Newmark approach revealed non-physical response of the classically solved problem. We can observe waves behind and ahead the load front which are results of reflections from the end of the beam. The properties of the presented numerical method allow the simulation of the infinite structure. The results coincides with the analytical solutions.

Notes on stability
This stability analysis is based on the assumption that the numerical method can not transfer the error with increasing amplitude from step to step. In order to investigate the stability of the method, the characteristic set of equations binding the nodal displacements with the velocities was built. Using (5) and (6) Fig. 15 Influence of the space-time finite element size on the spectral radius where T is the transition matrix. The Neumann necessary condition is used as a stability criterion: The eigenvalues of the transition matrix T must be within the unit circle. Their moduli satisfy the inequality ρ ≤ 1.
The analysis of a single space-time element of the Timoshenko beam (Fig. 3a) was carried out. According to (16), the element with four degrees of freedom leads to a matrix T of size 8 x 8. Thus, a characteristic polynomial of high order was obtained, and its roots can be calculated numerically. Generally, they are complex numbers. For the case V 2 2 > V 2 1 (V 2 1 + 1) discussed in subsection 4.1, the analysis of the influence of the ratio size of the space-time element on the spectral radius ρ was computed. The results of the calculations for the nodal displacements are shown in Fig. 15. The horizontal axis shows the dimensionless length of the finite element b which is related to the characteristic length l * . The vertical axis describes the dimensionless time step h, related to the characteristic length l * and the velocity v * . Both parameters are given using the notations in Sect. 2. According to Fig. 15, the numerical method exhibits damping at higher time step h. However, the stability region of the solution can provide propagation of information at a limited velocity with a wide range.

Conclusions
The space-time simplex shaped finite element approach simulates an unbounded system under a moving load. We considered the feature of the limited speed of information propagation in the finite element mesh. A special arrangement of the mesh allows us to avoid wave reflections from boundaries. This novelty can be used in practical engineering problems, especially in transportation. The time anisotropic property of the method allows us to reduce a large scale problem with a multi-node mesh to a small system that in a few steps settles down to a quasi-stationary solution. The reduced infinite systems, tested without the numerical tricks that others apply to suppress the reflected waves from the boundaries of the structure, exhibit very good coincidence of the numerical results with the analytical ones. Moreover, in the numerical analysis, we can observe properties of the solutions identical with those previously published that were obtained analytically. Thus, the space-time anisotropy affects on the correctness of the results of the numerical wave problems. The phenomenon of a limited speed of information propagation in the space-time mesh can be used in multidimensional problems of structures subjected to moving loads.