Instability of vibrations of an oscillator moving at high speed through a tunnel embedded in soft soil

This paper investigates the instability of vertical vibrations of an object moving uniformly through a tunnel embedded in soft soil. Using the indirect Boundary Element Method in the frequency domain, the equivalent dynamic stiffness of the tunnel-soil system at the point of contact with the moving object, modelled as a mass-spring system or as the limiting case of a single mass, is computed numerically. Using the equivalent stiffness, the original 2.5D model is reduced to an equivalent discrete model, whose parameters depend on the vibration frequency and the object’s velocity. The critical velocity beyond which the instability of the object vibration may occur is found, and it is the same for both the oscillator and the single mass. This critical velocity turns out to be much larger than the operational velocity of high-speed trains and ultra-high-speed transportation vehicles. This means that the model adopted in this paper does not predict the vibrations of Maglev and Hyperloop vehicles to become unstable. Furthermore, the critical velocity for resonance of the system is found to be slightly smaller than the velocity of Rayleigh waves, which is very similar to that for the model of a half-space with a regular track placed on top (with damping). However, for that model, the critical velocity for instability is only slightly larger than the critical velocity for resonance (of the undamped system), while for the current model the critical velocity for instability is much larger than the critical velocity for resonance due to the large stiffness of the tunnel and the radiation damping of the waves excited in the tunnel. A parametric study shows that the thickness and material damp- ing ratio of the tunnel, the stiffness of the soil and the burial depth have a stabilising effect, while the dam ping of the soil may have a slightly destabilising effect (i.e., lower critical velocity for instability). In order to investigate the instability of the moving object for velocities larger than the identiﬁed critical velocity for instability, we employ the D- decomposition method and ﬁnd instability domains in the space of system parameters. In addition, the dependency of the critical mass and stiffness on the velocity is found. We conclude that the higher the velocity, the smaller the mass of the object should be to ensure stability (single mass case); moreover, the higher the velocity, the larger the stiffness of the spring should be when a spring is added (oscillator case). Finally, in view of the stability assessment of Maglev and Hyperloop vehicles, the approach presented in this paper

Tunnel embedded in soft soil High-speed oscillator Vibration instability Indirect BEM D-decomposition method Critical velocity for instability a b s t r a c t This paper investigates the instability of vertical vibrations of an object moving uniformly through a tunnel embedded in soft soil. Using the indirect Boundary Element Method in the frequency domain, the equivalent dynamic stiffness of the tunnel-soil system at the point of contact with the moving object, modelled as a mass-spring system or as the limiting case of a single mass, is computed numerically. Using the equivalent stiffness, the original 2.5D model is reduced to an equivalent discrete model, whose parameters depend on the vibration frequency and the object's velocity. The critical velocity beyond which the instability of the object vibration may occur is found, and it is the same for both the oscillator and the single mass. This critical velocity turns out to be much larger than the operational velocity of high-speed trains and ultra-high-speed transportation vehicles. This means that the model adopted in this paper does not predict the vibrations of Maglev and Hyperloop vehicles to become unstable. Furthermore, the critical velocity for resonance of the system is found to be slightly smaller than the velocity of Rayleigh waves, which is very similar to that for the model of a half-space with a regular track placed on top (with damping). However, for that model, the critical velocity for instability is only slightly larger than the critical velocity for resonance (of the undamped system), while for the current model the critical velocity for instability is much larger than the critical velocity for resonance due to the large stiffness of the tunnel and the radiation damping of the waves excited in the tunnel. A parametric study shows that the thickness and material damping ratio of the tunnel, the stiffness of the soil and the burial depth have a stabilising effect, while the dam ping of the soil may have a slightly destabilising effect (i.e., lower critical velocity for instability). In order to investigate the instability of the moving object for velocities larger than the identified critical velocity for instability, we employ the Ddecomposition method and find instability domains in the space of system parameters. In addition, the dependency of the critical mass and stiffness on the velocity is found. We conclude that the higher the velocity, the smaller the mass of the object should be to ensure stability (single mass case); moreover, the higher the velocity, the larger the stiffness of the spring should be when a spring is added (oscillator case). Finally, in view of the stability assessment of Maglev and Hyperloop vehicles, the approach presented in this paper

Introduction
Dynamic effects are of significance for modern high-speed trains as propagating waves may be generated in the railway track and subsoil. The study of dynamic train-track-soil interactions has been of interest for researchers for decades. Popp et al. [1] gave a comprehensive review of the existing models that can be used to study the dynamic train-track-soil interaction.
In general, studies on moving trains fall into two categories. The first category is environmental vibrations induced by moving trains to assess vibration hindrance and to ensure the safety of the nearby structures. The second category is the instability of vibrations of moving trains to ensure the safety and comfort of the passengers in the trains. For the former category, the steady-state regime is assumed when investigating the dynamic amplification due to resonance [e.g., [2][3][4].
Other studies in the first category are devoted to transition radiation which occurs when a train passes an inhomogeneity [5,6] .
For studies falling into the second category, the train is usually modelled as a single-or multi-degree of freedom system [7] . When instability occurs, the free vibration (i.e., the vibration in the absence of an external force) of the train grows exponentially, resulting in an infinite displacement when time goes to infinity, which implies that a steady-state solution does not exist. This is very different from resonance, which happens when the steady-state response as induced by an external moving load is extreme (either bounded or not depending on the presence of damping). Both phenomena come with a certain critical velocity. The critical velocity for resonance is defined as the velocity at which the steady-state response induced by a moving load is extreme (i.e., resonance takes place at certain specific velocities) [8,9] , while the critical velocity for instability is defined as the velocity beyond which instability can occur (i.e., instability occurs in a range of velocities). Another crucial difference between resonance and instability is that resonance can be totally removed by increasing damping, while damping mostly shifts the instability domain, for example, to a region of larger velocities [10,11] .
The first study on instability of vibrations is that of a mass that moves uniformly along an elastically supported beam [10] . The physical explanation of instability was given by Metrikine [12] who argued that the instability is caused by the radiated anomalous Doppler waves [13] which increase the energy of the vibrating object. In addition, the physical mechanism of instability was discussed using the laws of conservation of energy and momentum [12,14] .
After the pioneering works on the instability phenomenon [10,15] , several aspects that influence the stability of an object have been discussed. For example, the effect of thermal stresses in the structure was studied considering different models of moving oscillators [11,16,17] . Other papers considered the effect of more than one contact point between the object and the structure [17] , and of contact nonlinearities [18] . Moreover, a more accurate beam model to represent the rail was considered and a comparison between the Timoshenko and Euler-Bernoulli beam models was given [19] . Four different beam and plate models were considered in [14] . Furthermore, Verichev et al. introduced other complexities in their model, i.e., a bogie model which consists of a rigid bar of finite length on two identical supports [20] . Recently, Mazilu studied the instability of a train of oscillators moving along an infinite Euler-Bernoulli beam on a viscoelastic foundation [21] . Another work focused on the stability of a moving mass in contact with a system of two parallel elastically connected beams, with one of them being axially compressed [22] . Later, the stability of vibrations of a railway vehicle moving along an infinite three-beam/foundation system has been considered, with an emphasis on the effect of the damping and stiffness of the secondary suspension of the railway vehicle [23] . For a more simple model, Dimitrovová presented a semi-analytical solution for the evolution of the beam deflection shapes and oscillator vibrations [24] . In that paper, not only the onset of instability, but also the severity is addressed.
All the above works used a one-dimensional or two-dimensional model of the railway track, which may be less accurate than three-dimensional models [25,26] to predict the instability of moving trains. However, they all convey the very important message that, in the presence of damping, the instability of moving trains may happen at speeds that exceed the critical velocity for resonance of the undamped system (which is equal to the minimum phase velocity of waves in the structure); that is, the critical velocity for instability is larger than the critical velocity for resonance. The few existing works related to instability analysis employing three-dimensional models of the railway track consider trains moving on a track founded on the ground surface [25,26] . It has been shown that the critical velocity for instability of the moving object is close to the Rayleigh wave speed in the soil. Instability of trains moving through an underground tunnel has not been analysed yet. In this paper, we therefore aim to conduct the instability analyses for an oscillator and the limiting case of a single mass moving through a tunnel embedded in soft soil. We will investigate whether the critical velocity for instability of the moving object is also close to the Rayleigh wave speed in the soil, for both a shallow and a deep tunnel. The re- sults are of practical relevance especially for contemporary high-speed railway tracks as well as upcoming ultra-high-speed transportation systems such as Maglev and Hyperloop, respectively [27][28][29] .
The paper is organised as follows. We present the model and a framework to conduct the instability analysis in Section 2 . Section 3 discusses the 2.5D Green's functions of a full-space and a half-space [30,31] , presents the Green's functions of the shell and the formulation of the indirect Boundary Element Method (BEM). Validations of the proposed indirect BEM are given in Section 4 . In Section 5 , we conduct the instability analysis of the single mass and the mass-spring oscillator. To this end, we study the equivalent dynamic stiffness to find the critical mass and stiffness, and analyse the effect of the tunnel thickness, the material damping ratios in the tunnel-soil system, the Lamé parameters of the soil and the burial depth of the tunnel on the critical velocity for instability. Moreover, the dependency of the critical mass and stiffness on the velocity is investigated. Conclusions are given in Section 6 .

Model description
In this paper, we study the vibrations of an object moving through a tunnel embedded in soft soil using a so-called 2.5D model. The soil is modelled as an elastic continuum, whereas the tunnel is modelled by the Flügge shell theory [32] . Both the soil and tunnel are assumed to be linear, elastic, homogeneous and isotropic. The soil is characterised by density ρ 1 , Poisson's ratio ν 1 , and complex Lamé parameters λ * 1 = λ 1 ( 1 + 2 i sgn (ω) ξ 1 ) and μ * 1 = μ 1 ( 1 + 2 i sgn (ω) ξ 1 ) , where i is the imaginary unit, ω is the frequency and ξ 1 the material damping ratio of the soil related to the adopted hysteretic damping model. The parameters of the tunnel are density ρ 2 , Poisson's ratio ν 2 , and complex Lamé parameters λ * 2 = λ 2 ( 1 + 2 i sgn (ω) ξ 2 ) and μ * 2 = μ 2 ( 1 + 2 i sgn (ω) ξ 2 ) , with ξ 2 being the hysteretic material damping ratio of the tunnel.
The burial depth of the tunnel is H , and its inner and outer radii are R i and R o . The object is modelled by a mass-spring oscillator (see Fig. 1 ), which is characterised by its mass M and spring stiffness K , and moves through the tunnel with a constant velocity V . Note that there is no vertical external force acting on the mass because the presence of such a force is irrelevant for the dynamic-instability analysis. Shallow and deep embedded tunnels are considered in this paper. For the shallow tunnel, the soil medium is modelled as a half-space, while for the deep tunnel, the soil is modelled as a full-space. Fig. 1 only shows the configuration of the shallow tunnel. If H → ∞ , it essentially becomes a deep tunnel.
The governing equations of the shell are presented later, in Section 3 . The current section only presents the framework to conduct instability analysis.

Method of solution
To analyse instability of vibrations of the moving object, the concept of the equivalent stiffness (also referred to as dynamic stiffness) is employed [33,34] . The procedure is illustrated in Fig. 2 and goes as follows. First, we compute the steadystate response of the system shown in Fig. 2 (a) which is subject to a uniformly moving oscillatory point load applied at the tunnel invert ( The oscillatory load has the form of P (t) = P 0 exp ( i t ) , in which P 0 is the amplitude, = 2 π f 0 is the angular frequency; f 0 is the load frequency in Hz. P ( t ) essentially represents a harmonic interaction force between the moving object and the tunnel-soil system. The steady-state radial displacement at the loading point can be expressed as U r 1 is the complex amplitude of this harmonic vibration. The indirect Boundary Element Method is employed to compute the response of the system, which is presented in detail in Section 3 . From the result, we obtain the equivalent stiffness of the tunnel-soil system at the loading point using the following relation: (1) By doing so, the original 2.5D model can be reduced to an equivalent discrete model, shown in Fig. 2 (b), consisting of a mass-spring system resting on an equivalent spring with a complex-valued stiffness K eq ( , V ), which depends on the frequency and velocity of the oscillator.
To study the instability of a moving oscillator, we apply, in accordance with previous dynamic-instability studies [11,24] , the Laplace integral transform with respect to time t ( s denotes the Laplace parameter) to the well-known governing equation of the vertical motion of the oscillator. Assuming zero initial conditions (which can be done as they do not influence the stability [11,24] ), the following characteristic equation for the free vibration of the oscillator is obtained: The roots of Eq. (3) determine the (complex) eigenfrequencies, = − i s, of the vertical motion of the oscillator as it interacts with the tunnel-soil system. If one of the roots s of the characteristic equation has a positive real part, the response will grow exponentially, which implies that the vertical vibration of the oscillator is unstable [11,19,25,26,35] . Obviously, the equivalent stiffness must be single-valued for all in order for Eq. (3) to be meaningful; see also Section 2.3 . It has been shown in [11] that the instability of a moving object may occur if and only if the imaginary part of the equivalent stiffness K eq is negative in a frequency band. In [11] , a single moving mass is considered, the motion of which is even necessarily unstable as soon as Im( K eq ) < 0 at any frequency band. The imaginary part of the equivalent stiffness can be considered to be the damping coefficient of the dashpot in the equivalent mass-spring system. A negative imaginary part of the equivalent stiffness indicates a negative damping, which makes the vibration of the moving mass unstable.
It would be very laborious to determine all the roots of the characteristic equation and check whether one of these roots has a positive real part. Alternatively, we follow a convenient method of root analysis, namely the D-decomposition method, to determine the number of 'unstable roots'. This method has been used in several papers [11,19,25,26,35] . The idea of this method is to map the imaginary axis of the complex s plane (i.e., the border between stability and instability) onto the plane of a system parameter, M or K , which is allowed to be complex. The mapped line divides the M or K plane into domains with different numbers of unstable roots. It is noted that the imaginary part of the complex system parameter has no physical meaning. Only the positive real part of the system parameter is physical, and the crucial question is whether one of the so-called instability domains overlay the positive real axis.
The procedure is as follows. Consider s = i , where serves as the parameter of the mapping, is real valued and has the meaning of frequency (same as introduced above), and has to be varied from minus to plus infinity. We discuss the following two cases in this paper. The first one is the limit case of a single mass moving through the tunnel, thus assuming K → ∞ . The characteristic equation for a single mass is reduced from Eq. (3) to Substituting s = i into Eq. (4) gives the following rule for the mapping: The second case we consider is the more general one of the moving oscillator, taking into account both the mass and the spring of the oscillator. In this case, the stiffness K will be used as the parameter for the D-decomposition assuming M to be constant because it is of practical relevance. Taking the limit case of the moving mass as the starting point, it is interesting to know what the added stiffness of the spring should be to render the oscillator vibration unstable (see also Section 5.2 ). Substituting s = i into Eq. (3) , we get the following mapping rule for the complex K plane: By replacing s by i in K eq ( s, V ), which essentially entails considering the limit case of s → i , one can use the equivalent stiffness ( Eq. (1) ) which is determined based on the steady-state response to the harmonic loading. Employing Eq. (5) or (6) as the mapping rule, one can plot the D-decomposition curve, for example, Im( M ) versus Re( M ) (as shown in Fig. 11 , for example), where is the running parameter along this curve. One side of the D-decomposition curve is shaded, and this side is related to the right-hand side of the imaginary axis in the s plane. Crossing the curve in the direction of the shading once indicates that there is an additional unstable root. Thus, one can find information on the relative number of unstable roots in domains of the complex M or K planes. The number of unstable roots in all the domains can be determined if the absolute number of those is known for any arbitrary value of the considered system parameter. By doing so, the instability domains can be found in the M or K plane, which generally allows to identify the critical velocity for instability (defined as the velocity at which an instability domain first overlays the positive real axis when increasing the velocity). The instability analysis is conducted in Section 5 .

Response to a moving oscillatory point load and derivation of equivalent stiffness
As shown in Section 2.2 , it is customary to employ the equivalent stiffness K eq to conduct the instability analysis. In the current section, we aim to derive the expression for the equivalent stiffness. To this end, we first derive the steady-state response to a moving oscillatory point load at the loading point. Additionally, the response at a fixed observation point is derived, which is needed for validating the indirect BEM in Section 4 . All the responses can be computed using the indirect BEM presented in Section 3 . Here we summarise the important steps and outcome in view of the specified aim.
The shear stresses σ r 1 θ 1 and σ r 1 x at the inner surface of the tunnel wall induced by the moving oscillatory point load (see Fig. 2 (a)) are zero. The non-zero normal stress σ r 1 r 1 (R i , θ 1 , x, t ) can be expressed as where δ(.) is the Dirac delta function.
As the considered problem is linear, we apply the Fourier Transform to derive the response of the system subject to the uniformly moving oscillatory point load in the wavenumber-frequency ( k x , ω) domain. The Fourier Transform applied with respect to time t and spatial coordinate x is defined in the following form (for an arbitrary function g ( r 1 , θ 1 , x, t )): with the inverse Fourier Transform given by The Fourier series, which is used to derive the response in the ( k x , ω) domain, of a general response quantity f ( θ 1 ) reads Expanding the term δ(θ 1 + π 2 ) in Eq. (7) into a Fourier series, the normal stress can be rewritten as Applying the Fourier Transform defined by Eq. (8) to Eq. (11) , the normal stress in the wavenumber-frequency domain is obtained as: The response induced by the auxiliary stress ˜ ˜ σ aux ( R i , θ 1 , k x , ω ) , which relates to a radial stress in the form of δ θ 1 + π 2 · δ(x ) δ(t) , can be computed using the indirect BEM ( Section 3 ) and is denoted as ˜ ˜ U 1 , aux ( r 1 , θ 1 , k x , ω ) . Thereafter, we get the expression of the actual displacement vector excited by the stress ˜ ˜ σ r 1 r 1 ( R i , θ 1 , k x , ω ) shown in Eq. (12) : We obtain the space-time domain response by applying the inverse Fourier Transform over wavenumber k x and frequency ω to Eq. (13) : In Eq. (14) , the inverse Fourier Transform over wavenumber k x has been evaluated analytically, whereas the inverse Fourier Transform over frequency ω needs to be evaluated numerically.
The radial displacement component of the steady-state response at the loading point can be obtained by substituting x = V t into Eq. (14) : In accordance with Eq. (15) , the complex amplitude of this harmonic response, which is relevant for the computation of the equivalent stiffness, is given as Using Eq. (16) , the equivalent stiffness K eq defined in Eq. (1) is obtained: We note that this result is single-valued for all as the Green's functions (see Section 3 ) used in the indirect BEM computations are uniquely defined.
We also consider the steady-state response at a fixed observation point x = 0 , which is needed for the validation of the indirect BEM. Substituting x = 0 into Eq. (14) gives the corresponding displacement vector: Eq. (18) contains the responses observed at x = 0 ; for an observation point at the tunnel invert, t < 0 indicates that x = 0 > V t, which means that the moving load has not reached the observation point yet; t = 0 indicates that the moving load is at the observation point; t > 0 indicates that x = 0 < V t, which means that the moving load has passed the observation point.
For the case of a stationary (i.e., non-moving) harmonic point load, which is also used in Section 4 for validation, an expression for the induced displacements is given in Appendix A .

Indirect boundary element method
In this paper, the indirect BEM is employed to compute the response of the tunnel-soil system in the wavenumberfrequency ( r 1 , θ 1 , k x , ω) domain. To this end, the Green's functions of the soil and tunnel are needed; note that the indirect BEM uses the Green's functions of the soil without cavity (full-space or half-space). In Sections 3.1 and 3.2 , the Green's functions of the soil and tunnel are presented. The indirect BEM is formulated in Section 3.3 .

Green's functions of the soil
The so-called two-and-a-half dimensional Green's functions of an elastodynamic full-space [30] and a half-space [31] are used in our work. The source considered in the mentioned papers is a spatially varying line load in the longitudinal direc- x indicates the direction of the load, subscript "s" indicates the coordinates of the source point, and F j is the amplitude of the source (Green's functions can be obtained by setting F j = 1 ).
The Green's functions of the half-space consist of source terms which are the same as those of the full-space, and of surface terms which are necessary to satisfy the stress-free boundary conditions at the surface of the half-space [31] . However, we found that the stress-free conditions are not satisfied using the Green's functions presented in [31] , while they are satisfied when the source terms are replaced by the ones presented in [30] that contains the full-space Green's functions. Therefore, the Green's functions of the half-space used in the current paper consist of the source terms presented in [30] and of the surface terms presented in [31] . Because the reference frames in [31] are different from that in the current paper, we have to transform the Green's functions for displacements and stresses given in [31] through the following relations ˜ ˜ where subscripts "1" and "ref" denote the responses defined in the coordinate systems of the current and reference papers, respectively. The transformation matrix reads: ˜ ˜ G u, 1 and ˜ ˜ G σ, 1 are the Green's functions for displacements and stresses of the soil without tunnel/cavity (i.e., of the fullspace or half-space), and are 3 × 3 matrices. In matrix ˜ ˜ G u, 1 , the first, second and third rows represent the displacement components ˜ ˜ u r 1 , ˜ ˜ u θ 1 and ˜ ˜ u x , while the first, second and third columns correspond to the spatially varying unit line loads acting in y, z and x directions, respectively. In matrix ˜ ˜ G σ, 1 , the first, second and third rows represent the stress components ˜ ˜ σ r 1 r 1 , ˜ ˜ σ r 1 θ 1 and ˜ ˜ σ r 1 x , while the columns also correspond to the loads in different directions.

Green's functions and responses of a cylindrical shell
The tunnel is modelled by an infinitely long cylindrical Flügge shell. The associated coordinate system is shown in Fig. 3 . The equations of motion of the shell read [32] : where ū , v and w are the mid-surface displacements in directions r 2 , θ 2 and x 2 , respectively. E 2 is the Young's modulus of the shell and h its thickness. The radii of the inner and outer surface of the shell can be expressed as R i = R − h 2 and R o = R + h 2 , respectively. q r 2 , q θ 2 and q x 2 are the net external stresses acting on the shell, namely the difference between the stresses acting at the inner and outer surfaces.
The governing equations of the shell can be rewritten into matrix form as where ū 2 = ( ū , v , w ) is the displacement vector, q 2 = q r 2 , q θ 2 , q x 2 is the net stress vector corresponding to the mid-surface of the shell, and A is an operator matrix given in Appendix B .
The stress vector q 2 is related to the stress vectors q o 2 and q i 2 corresponding to the outer and inner surfaces of the shell, respectively, through the following relations [36] : According to Love's simplification in the shell theory [37] , the longitudinal and tangential displacements vary linearly across the shell's thickness, whereas the radial displacement is independent of radial coordinate. Therefore, the mid-surface displacement vector ū 2 is related to the displacement vector u o 2 corresponding to the outer surface of the shell through the following relation After applying the Fourier Transform over time t and spatial coordinate x 2 to Eq. (23) , computing the Fourier coefficients of the circumferential harmonics (i.e., the second relation in Eq. (10) ), and considering Eqs. (24) and (25) , the governing equations of the shell can be written as which is essentially a set of algebraic equations; n denotes the number of the circumferential harmonic, and matrices ˜ ˜ The associated Green's functions of the shell can be derived by solving Eq. (26) for each of the load components, and subsequently adding the solutions for all components in the Fourier series (see Eq. (10) ): where ˜ ˜ g o n interrelates ˜ ˜ u o 2 ,n and ˜ ˜ q o 2 ,n , ˜ ˜ g i n interrelates ˜ ˜ u o 2 ,n and ˜ ˜ q i 2 ,n , and ˜ ˜ g o , ˜ ˜ g i , ˜ ˜ g o n and ˜ ˜ g i n are 3 × 3 matrices. The positive directions of the longitudinal axes in the global coordinate system (see Fig. 1 ) and the local coordinate system for the shell (see Fig. 3 ) are opposite to each other (i.e., x 2 = −x ). Therefore, the relation between the longitudinal wavenumbers k x 2 and k x (which, for the moving oscillatory point load considered in Section 2.3 , is defined as ω− V ) is as follows: k x 2 = −k x ; this relation is used below.
Using the convolution rule, the displacement vector of the shell under an arbitrary load can be obtained as The displacements in Eq. (28) are defined in the local coordinate system of the shell ( r 2 , θ 2 , x 2 , see Fig. 3 ). To satisfy the continuity of displacements and stresses at the shell-soil interface, the displacement and stress vectors ˜ ˜ u o 2 , ˜ ˜ q o 2 and ˜ ˜ q i 2 defined in the local coordinate system of the shell have to be transformed to ˜ ˜ U o 2 , ˜ ˜ Q o 2 and ˜ ˜ Q i 2 , respectively, defined in the global cylindrical coordinate system of the soil ( r 1 , θ 1 , x ), which has origin at the center of the tunnel (see Fig. 1 in which θ 1 = θ 2 and Substituting Eq. (29) into Eq. (28) , we obtain the displacements of the shell in the global cylindrical coordinate system:

Formulation of the indirect boundary element method
In this section, the formulation of the employed indirect BEM is presented.
where ˜ ˜ U 1 and ˜ ˜ U 2 are the displacement vectors of the soil and tunnel, respectively, and ˜ ˜ 1 and ˜ ˜ 2 their stress vectors.
The displacement and stress vectors of the soil and shell at the tunnel-soil interface are expressed as ˜ ˜ Note that all these displacements and stresses are defined in the global cylindrical coordinate system ( r 1 , θ 1 , x ).
According to the indirect BEM, the displacement and stress vectors in the soil are given as [36] ˜ ˜ where F ( x s ) is the yet unknown vector of the source amplitudes placed inside the fictitious cavity which is commonly used for the indirect BEM (see Fig. 4 ). Vectors x r = [ x r , y r , z r ] and x s = [ x s , y s , z s ] are coordinates of the receiver and source points, respectively. L s is the surface at which the source points are located, and the radius of the surface L s is taken as where N r denotes the number of receiver points, and N r ≥ 20 as suggested in [36] . The surface L r at which the receiver points are located lies at r 1 = R o , which is the outer surface of the actual tunnel. An expression for the displacement vector of the shell has been obtained in Eq. (31) , and can be rewritten as where ˜ ˜ G 9 Table 1 Displacement components (20 ·log 10 | U i |) at different locations ( r 1 , θ 1 , x ) for a tunnel embedded in an elastic full-space subjected to a stationary harmonic point load with excitation frequency f 0 = 10 Hz obtained using different numbers of source and receiver points ( N s , N r ). In each row, the displacements are normalised by the corresponding response obtained using (N s , N r ) = (20 , 40) .  Fig. 2 (a)).
Substituting Eqs. (34) and (36) into Eq. (32) , and considering Eqs. (35) and (37) , we derive the boundary integral equation in terms of the unknown source amplitude vector: In order to compute the source vector for every k x and ω combination, Eq. (38) is discretised (i.e., surfaces L r , L and L s ), as indicated above. Note that the size of the matrix related to ˜ ˜ K is (3 N r × 3 N s ), where N s denotes the number of source points; this implies that it is a small matrix, and the inversion of the matrix does not take much time. The most time consuming part (even in the entire procedure of instability analysis) is the evaluation of integrals over the horizontal wavenumber k y in the Green's functions of the half-space, for which we use the "quadv" routine in Matlab.

Validations
To validate the accuracy of the presented indirect BEM, we compare the results obtained by the proposed method with those calculated by Yuan et al. [38] . The first validation is performed for the case of a tunnel embedded in an elastic fullspace subject to a stationary harmonic point load at the tunnel invert. The excitation for the indirect BEM computation in this case is ˜ σ aux ( Eq. (A.3) ), with P 0 = 1 N , and the steady-state response is given in Eq. (A.5) . The elastic full-space is characterised by its longitudinal wave speed C P , 1 = 944 m / s , shear wave speed C S , 1 = 309 m / s , density ρ 1 = 20 0 0 kg / m 3 and material damping ratio ξ 1 = 0 . 03 . The elastic parameters for the tunnel are the Young's modulus E 2 = 50 GPa , Poisson's ratio ν 2 = 0 . 3 , density ρ 2 = 2500 kg / m 3 and material damping ratio ξ 2 = 0 . The inner and outer radii of the tunnel are Before showing the results, we first present a convergence test for the proposed method for the considered loading case. As shown in Eq. (A.5) , the inverse Fourier Transform over longitudinal wavenumber k x has to be evaluated to get the harmonic response in the space-time domain. The integral was computed numerically using an inverse fast Fourier Transform algorithm in Matlab. The convergence was tested regarding the discretisation of k x (i.e., k x and k max x ), the maximum number of circumferential modes of the shell N max shell in Eq. (31) and the maximum number of Fourier components N max load in Eq. (A.3) , and the number of source and receiver points ( N s , N r ). We found that it is sufficient to use k max x = 2 π , k x = 2 π 1023 , N max shell = 20 and N max load = 20 . The convergence test for the number of source and receiver points ( N s , N o ) at different locations ( r 1 , θ 1 , x ) is given in Table 1 . Responses at the tunnel invert, tunnel apex ( R o , π 2 , 0), tunnel side ( R o , π , 0) and at a point far from the load (20 m, π , 20 m) are presented; the load is characterised by P 0 = 1 N and f 0 = 2 π = 10 Hz . It is clear that converged results can indeed be obtained using ( N s , N r ) = (20, 40). Fig. 5 shows the converged vertical displacements at the tunnel invert, tunnel apex and tunnel side as a function of frequency for the first validation case. A good agreement can be observed between the results obtained by different methods, which validates the proposed method and its implementation.
The second validation case is that of a shallow tunnel embedded in an elastic half-space subject to a stationary harmonic load. The soil is characterised by its longitudinal wave speed C P , 1 = 400 m / s , shear wave speed C S , 1 = 200 m / s , density ρ 1 = 1800 kg / m 3 and material damping ratio ξ 1 = 0 . 02 . The parameters of the tunnel are the same as in the previous case, Fig. 5. Vertical displacements ( 20 · log 10 | U z 1 | ) at locations of (a) tunnel invert ( tunnel side ( r 1 = R o , θ 1 = π, x = 0 ) for the case of a tunnel embedded in an elastic full-space subject to a stationary harmonic point load.

Table 2
Velocities ( V i ( y, z, x )) and displacement (U r 1 (r 1 , θ 1 , x )) at different locations for a tunnel embedded in an elastic half-  except that ξ 2 = 0 . 015 , and the burial depth of the tunnel is H = 5 m . It is noted that in the reference paper [ 38 ], the mentioned material damping ratio should be the loss factor (there is a difference of a factor 2), which is indicated in paper [39] . This also holds for the next validation case. The displacement components (again obtained using Eq. (A.5) ) at a point on the ground surface ( y = −20 m , z = 0 , x = 20 m ) are presented in Fig. 6 , where again a good match between the results is observed. The minor differences can be attributed to the use of a continuum to model the tunnel in [38] , instead of a shell. The third validation comprises the case of a tunnel embedded in an elastic half-space subject to a uniformly moving constant point load. The excitation for the indirect BEM computation in this case is ˜ ˜ σ aux ( Eq. (12) ) and the steady-state response is given in Eq. (14) . The parameters for the soil and tunnel are as follows: μ 1 = 1 . 154 × 10 7 N / m 2 , λ 1 = 1 . 731 ×  (14) ). The convergence for the moving point load case was tested regarding the discretisation of ω (i.e., ω and ω max ), N max shell in Eq. (31) and N max load in Eq. (12) , and the number of source and receiver points ( N s , N r ). Numerical results related to two points on the ground surface and one at the tunnel invert are presented in Table 2 for the considered case of the moving load using different numbers of source and receiver points ( N s , N r ). We found that converged results can be obtained using f max = ω max 2 π = 15 Hz , f = ω 2 π = 0 . 05 Hz , N max shell = 20 , N max load = 20 and ( N s , N r ) = (20, 40). This is clear from  Table 2 which presents the responses observed at x = 0 for varying time moments: t = 0 means that the load is right below the observation point, whereas t = 1 s indicates that the load has passed that point. Fig. 7 presents the comparison between the results obtained by the proposed method and those shown in the literature. The good agreement gives confidence about the accuracy of the proposed method. The convergence requirements for the computation of the equivalent dynamic stiffness for different velocities are presented in Appendix C .

Instability of vibrations
The main framework to conduct instability analysis has been given in Section 2.2 . In the current section, we present the results for both the full-space and half-space. The base-case parameters of the tunnel-soil system are listed in Table 3 ; it is noted that the base case assumes that the burial depth H → ∞ . The parameters presented in Table 3 represent a soft soil and a concrete tunnel.
We chose the full-space case as the base case simply because the results for the half-space are pretty similar, and the computation of the two-and-a-half dimensional Green's functions of the full-space [30] is less expensive than that of the half-space Green's functions [31] . The Green's functions of the full-space are available analytically and can be evaluated very fast; however, the surface-terms part of the Green's functions of the half-space (see Section 3.1 ) are not available analytically, and integrals over the horizontal wavenumber k y need to be evaluated numerically. Table 3 Base-case parameters of the tunnel-soil system.

Critical velocity for instability of the moving object
As has been discussed in Section 2.2 , the imaginary part of the equivalent stiffness being negative indicates that the vibration of the object (i.e., moving mass or oscillator) can become unstable. Therefore, we first study the equivalent stiffness to find the critical velocity for instability of the moving object (here defined as the velocity at which Im( K eq ) < 0 first takes place). Note that the critical velocity for instability generally differs from the classical critical velocity at which the steadystate response induced by a moving load is extreme (i.e., resonance). In the general case, where the oscillator has dissipative components, the critical velocity for instability should be identified from the D-decomposition curve in the complex M or K plane (see Sections 2.2 and 5.2 ), not from the analysis of Im( K eq ) alone [19,26] . Im( K eq ) < 0 is only a necessary condition for instability. As the moving oscillator and moving mass considered in this paper do not have intrinsic dissipative components, their critical velocities for instability are the same.  Table 3 , and V inst cr = 942 m / s . In previous studies [11,19] where beam on elastic foundation models (without damping) are considered, it is shown that the critical velocity for instability V inst cr of the moving object is equal to the critical velocity for resonance (of the undamped system), which in turn is equal to the minimum phase velocity V min ph of waves in the system. For a half-space model with a regular track on top [25,26] , the critical velocity for instability in the presence of damping is slightly larger than V min ph . Additionally, V min ph , which is close to and smaller than the velocity of Rayleigh waves, is easily found from the dispersion relation of the system [12,40] . However, for the tunnel-soil system considered in this paper, it is very difficult to get the dispersion curves, as the dispersion characteristics of the system are considerably more complicated. Therefore, the minimum phase velocity cannot be easily computed. We can, however, compute the steady-state response of the tunnelsoil system subject to a uniformly moving non-oscillatory load and check the features of responses for different velocities to determine the critical velocity for resonance (for the system with damping, strictly speaking, but the influence of the damping on V res cr is small). This analysis is presented in Appendix C , and it shows that V res cr ≈ 70 m / s for the current tunnelsoil system, which is also close to and smaller than the velocity of Rayleigh waves, like for the above-mentioned half-space model.
For the system with the base-case parameters, we find that the imaginary part of the equivalent stiffness starts having a negative sign for at least a small frequency range at a velocity of V inst cr = 942 m / s . Based on this critical velocity for instability, we study the behavior of K eq ( , V ) for different velocities in the range of (0 . 9 − 1 . 2) V inst cr . The real and imaginary parts are shown in Figs. 8 and 9 , respectively. Nine different velocities were chosen to show the features of Re( K eq ) and Im( K eq ) as a function of the load frequency . Note that, if its mass is relatively small, the vibration of the object can still be stable when it moves faster than V inst cr , as will be demonstrated in the next section (see also Fig. 11 ). In Fig. 8 , we observe that the real part of the equivalent stiffness is positive, and the decaying trend of Re( K eq ) with frequency is similar for each velocity. The decaying trend may be related to the effect of inertia. The trough can probably  Table 3 , and V inst cr = 942 m / s . The read dots indicate the crossings. be interpreted as a quasi-resonance which takes place at low frequencies and is related to the wave resonance which occurs if the velocity of the moving load is the same as the group velocity of a wave excited by the load [25] .
The imaginary part of the equivalent stiffness is shown in Fig. 9 for each of the chosen velocities. For the 'sub-critical' ( V < V inst cr ) case shown in Fig. 9 (a), V = 0 . 9 V inst cr , Im( K eq ) is positive for all the frequencies , which indicates that the damping coefficient of the equivalent mass-spring system is positive, and thus, the system is always stable (see Section 2.2 ).
The frequency band considered in this study of the dynamic stiffness is limited to (0 − 40) Hz , because instability is determined by the behaviour at low frequencies [25,26] ; see also the explanation given at the end of this section. For the 'critical' and 'super-critical' ( V ≥ V inst cr ) cases, shown in Figs. 9 (b) -(i), the imaginary part of the equivalent stiffness is negative at low frequencies and becomes positive at higher frequencies. We can verify that the curves of Im( K eq ) in the high-frequency band have the same trend as that in the sub-critical case; they are not shown here since we focus on features of Im( K eq ) in the low-frequency band. There are peaks and troughs in the curves of Im( K eq ), and these are suppressed or enlarged as the velocity increases. We observe that the Im( K eq ) curve crosses the real axis 0, 4, 3, 7, 5, 3, 3, 1 and 1 times for V = (0 . 9 , 1 . 00 , 1 . 01 , 1 . 02 , 1 . 03 , 1 . 04 , 1 . 05 , 1 . 06 , 1 . 20) V inst cr , respectively. We can verify that for velocities cr , similar features of Im( K eq ) are observed (i.e., the Im( K eq ) curve crosses the real axis only once), the only difference is that the crossing occurs at higher frequency as the velocity increases (see Figs. 9 (h) and (i)). The different features of Im( K eq ) observed in the entire considered velocity range imply that in the complex M or K plane different amounts of separated domains, each having a specific number of 'unstable roots/eigenvalues', are expected for different velocities (see Section 5.2 ). At this point, it is concluded that the equivalent dynamic stiffness, especially its imaginary part, strongly depends on V .
In order to trace similarities and differences, let us now compare two different instability problems: the above mentioned model of an object moving on a track placed on the ground surface [25,26] , and the current model of an object moving through a tunnel embedded in a half-space. The critical velocity for instability for the current model with all parameters in accordance with the base case, except the burial depth H which is taken as 15 m, is found to be 891 m/s (see also Section 5.3.4 ). Clearly, V inst cr is much larger than the critical velocity for resonance ( V res cr ≈ 70 m / s ) in the model with an embedded tunnel, while V inst cr is just slightly larger than V res cr in the model with a track directly placed on the ground [25,26] ( V res cr is related to the undamped system in these studies, but the influence of the damping on V res cr is small). The difference is due to the large stiffness of the tunnel and the radiation damping/leaky character of the waves excited in the tunnel. However, there are similarities regarding ground vibrations in these two models. In the regime of V < V res cr ( ≈ V min ph ) , for both models mostly the medium in the vicinity of the load is disturbed by the eigenfield excited by the moving nonoscillating object, while in the regime of V > V res cr , both the vicinity of the moving source and the field far from the source are disturbed because waves are generated (see Appendix C ). As mentioned above, the critical velocity for resonance of the current tunnel-soil system is close to and smaller than the velocity of Rayleigh waves, like that of the other model. Therefore, from the ambient-vibration point of view, there is a clear similarity between both problems. However, instability happens only far beyond the critical velocity for resonance for the model with the tunnel, which is clearly different from the finding for the half-space with a track placed on top.
As shown in [12] , an external source has to supply a vibrating object with energy in order to maintain its uniform motion. In the case of unstable vibrations, the work done by the source is partially transferred to vibration energy of the object by the so-called anomalous Doppler waves [13] , which are waves of negative frequency. Typical dispersion curves of the current tunnel-soil system, which are similar to the ones of the beam on elastic foundation model [35] , are shown in Fig. 10 and can be used to explain why the instability only happens in the low-frequency band, as stated above. Fig. 10 also shows the so-called kinematic invariant ω = k x V + , which is essentially found in the argument of the Dirac function in the response to a moving oscillatory load (see Eq. (14) ). The kinematic invariant is a straight line indicating the relation between the load frequency , and the frequency ω and wavenumber k x of the waves that are potentially excited by the moving object; different realizations (i.e., being zero and nonzero, together with two different velocities) are shown in Fig. 10 . Intersections of the kinematic invariant with the dispersion curves represent the excited waves. We observe that intersections with negative frequency ω (i.e., anomalous Doppler waves) are only possible when the load frequency is relatively small. If the load frequency is large, the kinematic invariant will practically never intersect the dispersion curves at negative frequency, which explains why the vibration of the moving object is always stable in the high-frequency band (i.e., Im( K eq ) > 0, see Fig. 9 ); this also justifies that we restricted the analysis of K eq to the low-frequency band in Figs. 8 and 9 .

D-decomposition: complex M and K planes
In order to investigate the instability of the object vibrations for velocities larger than the corresponding critical velocity for instability identified in the previous section, we apply the D-decomposition method. We first investigate the limit case of the single mass moving through the tunnel. Considering the base case, the D-decomposition curve can be plotted in the complex M plane (i.e., Im( M ) versus Re( M )) using the mapping rule shown in Eq. (5) and is presented in Fig. 11 . For most of the considered velocities, the D-decomposition curve crosses the positive real axis, that is, one or more crossing points M * are obtained. It can be verified that the frequency at which the curve crosses the real axis corresponds to the frequency at which the imaginary part of the dynamic stiffness changes sign (see Fig. 9 ). A crossing point lying on the positive real axis can be explained by the fact that Re( K eq ) is positive when the Im( K eq ) changes its sign (see Fig. 8 ).
As it is clearly shown in Fig. 11 , the crucial difference between the D-decomposition curves in the super-critical and sub-critical cases -compare Figs. 12 (b) and (a), for example -is that there are crossing points M * 1 , M * 2 , M * 3 and M * 4 on the positive axis of Re( M ) for the super-critical case. The existence of such crossing points means that the number N of . The procedure to determine N is as follows. The relative number of unstable roots in domains of the complex M plane can be calculated by counting the number of times that one crosses the D-decomposition curve in the direction of the shading, which has been explained in Section 2.2 . To get the absolute number in all domains, the number of unstable roots for M = 0 has to be determined. M = 0 means that there is essentially no moving mass, which implies that the vibration of the mass cannot be unstable (i.e., the number of unstable roots N = 0 ). Thereafter, the absolute number of unstable roots in each domain can be determined and the result is shown in Fig. 11 . The number of unstable roots has also been validated using the Argument Principle, but this is not shown in the paper.
In Fig. 11 , we observe that the vibration of the moving mass is stable for all values of M for V = 0 . 9 V inst cr ; for V = 1 . 00 V inst cr , the vibration of the moving mass is unstable when M Note that the vibration of a mass which moves faster then the critical velocity is not necessarily unstable. For relatively small values of the mass, for example, the vibration is stable even for super-critical velocities (as illustrated in Section 5.3 ).
The question of practical relevance when studying the instability of the moving mass is whether adding flexibility (by creating a spring between the mass and the tunnel) may destabilize the system. For the mass-spring oscillator with the mass being constant, the D-decomposition curve can be plotted in the complex K plane (i.e., Im( K ) versus Re( K )) using the mapping rule shown in Eq. (6) and is presented in Fig. 12 . The mass of the moving oscillator is taken as M = 2 × 10 4 kg , which is a realistic value for a train wagon. In order to get the absolute number of unstable roots in the complex K plane with this mass, we have to connect the instability analysis of the moving oscillator to that of the single moving mass shown in Fig. 11 . From that figure, we find that N = 0 for M = 2 × 10 4 kg , which implies that the system is stable for this value of the mass for all the considered velocity cases. The single mass case corresponds to the oscillator case with K → ∞ . Therefore, knowing that the number of unstable roots at K → ∞ is zero and following the direction of the shading, the absolute number of unstable roots in domains of the complex K plane can be determined.
The following can be observed from the D-decomposition curve in the complex K plane shown in Fig. 12 . For V = 0 . 9 V inst cr , the vibration of the moving oscillator is stable for all values of the stiffness; for V = 1 . 00 V inst cr , the vibration of the moving oscillator is destabilized by the added spring when K * 1 < K < K * 2 and K * 3 < K < K * 4 ( N = 2 ); for V = 1 . 01 V inst cr , that happens when K < K * 1 and K * 2 < K < K * 3 , etc. Using these findings, it can be readily concluded that for the sub-critical case, the vibration of the oscillator is stable independently of the oscillator's stiffness, while for the super-critical cases, the stability of the oscillator depends on the stiffness of the added spring.

Parametric study
In Section 5.1 , we found the critical velocity beyond which instability of the moving object may occur. As the critical velocity for instability is the most important outcome of the instability analysis, the effects of the tunnel thickness, the material damping ratios in the tunnel-soil system, the Lamé parameters of the soil and the burial depth of the tunnel on the critical velocity are studied here. In addition, the dependency of the critical mass and stiffness (identified in Section 5.2 ) of the corresponding moving mass and moving oscillator on velocity is considered; this is only done for full-space cases because of the computational demand of the calculations for the half-space.

The effect of the thickness of the tunnel
Four different thicknesses of the tunnel are considered, and the corresponding critical velocities for instability are shown in Table 4 , where the subscript "B" (also shown in  indicates the parameters of the base case shown in Table 3 . Table 4 shows that the critical velocity for instability decreases as the tunnel thickness decreases. The reason of this reduction is the reduction of the stiffness of the tunnel. Moreover, we observe that even for the thinnest tunnel with thickness

The effect of the material damping ratios in the tunnel-soil system
We consider four different combinations of the material damping ratios of the soil and tunnel as shown in Table 5 . It demonstrates that the critical velocity for instability increases as the damping ratio of the tunnel increases and that of soil decreases. Therefore, we can conclude that the material damping of the tunnel stabilises the vibration of the moving object, while the material damping ratio of the soil may have a destabilising effect, which is similar to the finding in [26] . Table 5 Critical velocities for instability for different material damping ratios of the soil and tunnel.
Damping ratios of the soil and tunnel

The effect of the Lamé parameters of the soil
Three sets of the Lamé parameters of the soil are considered, see Table 6 . It shows that the critical velocity for instability of the moving object increases as the Lamé parameters of the soil increase. Thus, the stiffness of the soil has a stabilising effect on the vibration of the moving object, which is in line with the literature finding [26] . Table 6 Critical velocities for instability for different Lamé parameters of the soil.

The effect of the burial depth of the tunnel
It is interesting to compare the critical velocity for instability for the full-space and half-space cases. Table 7 shows that V inst cr decreases as the depth of embedded tunnel decreases. This reduction is probably and mostly because the Rayleigh wave, which is slower than the body waves in the soil, also starts to play a role, although its influence is not very large.

Dependency of the critical mass and stiffness on velocity
In this section, three cases (see Table 8 ) are considered to investigate the dependency of the critical mass and stiffness of the moving mass and moving oscillator, respectively, on the velocity in the range of V = (1 . 00 − 1 . 20) V inst cr . We chose three full-space cases to investigate the dependency relationship. As observed in Figs. 11 and 12 , there are many critical masses and stiffnesses for some velocities. For these velocities, we only consider the smallest critical mass M * 1 and the largest critical stiffness K * max because we are interested to find the regions where the vibration of the moving oscillator is stable   Table 8 , and three small intervals are considered to clearly show the trend in each interval. The region below the line relates to stable vibrations of the single mass.  Table 8 , and three small intervals are considered to clearly show the trend in each interval. The region above the line relates to stable vibrations of the oscillator. The reason is that the critical mass and stiffness can decrease and increase dramatically as the velocity increases, and by distinguishing the three sub-intervals, we clearly show the trend in each interval. In Fig. 13 , we observe that the critical mass in all three cases decreases as the velocity increases, which is in line with [11] . Fig. 13 also shows that the critical mass of the moving object (single-mass case) in case III is the largest and the one in case I is the smallest when V = V inst cr , and that the difference between the three critical masses is very large. However, the critical mass in case III becomes the smallest and the one in case I the largest when V = 1 . 20 V inst cr , and the difference between the three critical masses becomes much smaller. Fig. 14 shows that the critical stiffness of the oscillator in the three cases increases as the velocity increases, which is again in line with the literature finding in [25] . Another similarity between our result and the literature is that the instability of the moving oscillator occurs when its stiffness is in the order of 10 6 (kg/s 2 ), which is approximately the same as the value of the stiffness of springs used for conventional trains. However, the critical stiffness increases dramatically up to the order of 10 8 (kg/s 2 ) for our model while K * stays in the same magnitude as the velocity increases for the model considered in the literature [25] . Fig. 14 also shows that the critical stiffness in case III is the smallest and the one in case I is the largest when V = V inst cr , and that the difference between the three critical stiffnesses is very small. However, when V = 1 . 20 V inst cr , the critical stiffness in case III becomes the largest and the one in case I the smallest, and the difference between the three critical stiffnesses becomes much larger. Finally, Figs. 13 and 14 show that the dependency of the critical mass and stiffness on the velocity is similar in the considered velocity range for the three cases; however, as is clear from the comparison of cases II and III, the Lamé parameters of the soil have a larger effect on the curves compared to the damping ratio of the tunnel.
In Fig. 13 (14) , the region below (above) the line relates to stable vibrations of the object. In the region above (below) the line, the object vibration can be either purely unstable, or alternately either stable or unstable, which is the case when Im( K eq ) has many zero crossings (see Figs. 11 and 12 ). In case I, for velocities V = (1 . 06 − 1 . 20) V inst cr , the vibration of the moving mass in the region above the line shown in Fig. 13  in the subregions of K * 1 < K < K * 2 and K * 3 < K < K * 4 , but stable in the sub-regions K * 2 < K < K * 3 and K < K * 1 . In case II, for velocities V = 1 . 00 V inst cr and V = 1 . 02 V inst cr , the number of the critical masses and stiffnesses can be verified to be 2 and 9, respectively; for velocities V = (1 . 04 − 1 . 20) V inst cr , however, the number is 1. In case III, in the velocity range of V = (1 . 00 − 1 . 08) V inst cr , the number of the critical masses and stiffnesses varies significantly and jumps from 1 to 11 and then back to 3; for velocities cr , the number is again 1. Clearly, the precise number of sub-regions highly depends on the stiffness of the soil, and on the damping ratios of soil and shell.

Conclusions
In this paper, instability of the vibration of an object moving through a tunnel embedded in soft soil has been studied. We employed the concept of the equivalent dynamic stiffness, which reduces the original 2.5D model to an equivalent discrete model, whose parameters depend on the vibration frequency and the object's velocity. The frequency-domain indirect Boundary Element Method was used to obtain the equivalent stiffness of the tunnel-soil system at the point of contact with the moving object (i.e., the mass-spring system and the limit case of a single mass). Prior to that, the indirect BEM was validated for specific problems: the response of the system to a stationary harmonic point load and to a moving non-oscillatory load acting at the invert of a tunnel. Using the equivalent stiffness, the critical velocity beyond which the instability of the object may occur was found (it is the same for both the moving mass and the moving oscillator). The critical velocity for instability is the most important result of the instability analysis. We found that the critical velocity for instability turns out to be much larger than the operational velocity of high-speed trains and ultra-high-speed Hyperloop pods, which implies that the model adopted in this paper predicts the vibrations of these objects moving through a tunnel embedded in soft soil to be stable.
For the model of a track founded on top of the elastic half-space, considered for comparison, the critical velocity for instability in the presence of damping is just slightly larger than the critical velocity for resonance of the undamped system (which is equal to the minimum phase velocity of the system). However, for the current model, the critical velocity for instability is much larger than the critical velocity for resonance (of the damped system, strictly speaking, but the influence of the damping on the resonant velocity is small). For both models, the critical velocity for resonance is slightly smaller than the velocity of Rayleigh waves, and the fact that the critical velocity for instability is so much larger in the model with the embedded tunnel is due to the large stiffness of the tunnel and the radiation damping of the waves excited in the tunnel. Other parameters affect the instability as well. A parametric study shows that the thickness of the tunnel, the material damping ratio of the tunnel, the stiffness of the soil and the burial depth have a stabilising effect, while the damping of the soil may have a slightly destabilising effect.
In order to investigate the instability of the moving object in case the velocity exceeds the identified critical velocity for instability, we employed the D-decomposition method and found the instability domains in the space of system parameters. For a deep tunnel, the dependency of the critical mass and stiffness on the velocity was investigated. We conclude that the higher the velocity, the smaller the mass of the object (single mass case) should be to ensure its stability. Furthermore, the higher the velocity, the larger the stiffness of the spring should be when the spring is added (oscillator case). Our findings regarding the velocity dependency of the critical mass and stiffness are aligned with the conclusions obtained by Metrikine et al. [11,25] for other models.
The fact that the critical velocity for instability for the current model is much higher than the operational velocity of contemporary and future vehicles is promising for the Maglev and Hyperloop transportation systems. Furthermore, the approach presented in this paper can be applied to more advanced models with more points of contact between the moving object and the tunnel, which would resemble reality even better. Finally, as the dynamic stiffness is very important for the instability analysis for the tunnel-soil system, a refined model of the tunnel, which can potentially increase the accuracy of the response at its interior, can be considered in future work.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  In this appendix, our first aim is to find the critical velocity for resonance V res cr of the current tunnel-soil system. To this end, we analyse the response of the tunnel-soil system induced by a uniformly moving non-oscillatory ( f 0 = 2 π = 0 ) point load observed at the tunnel invert at x = 0 as derived in Section 2.3 (see Eq. (18) Table 3 . Fig. C.1 (a) shows the amplitude spectra of the radial displacements observed at the tunnel invert, at the observation point x = 0 , for one sub-Rayleigh, one super-Rayleigh and one supersonic case (the others look similarly). The Fourier transformed displacement, ˜ U r 1 (r 1 , θ 1 , x, ω) , is defined as the integrand of Eq. (18) except for the term exp ( i ωt) . It is shown that the spectra are spread around f = 0 . The time-domain responses are shown in Fig. C.1 (b). Clearly, for V = 30 m / s and V = 69 m / s , the disturbance is localized around the moving load, and there is no significant wave radiation in these    cases. For V = 70 m / s , a wave pattern emerges, which comes with significant asymmetry of the profile. For V = 75 m / s and V = 150 m / s , a more clear wave pattern can be observed; in these two cases, Rayleigh waves, and Rayleigh, shear and compressional waves are generated, respectively. Furthermore, the response is extreme for V = 70 m / s , which indicates resonance. Therefore, we conclude that for the current tunnel-soil system, V res cr ≈ 70 m / s . Clearly, a constant load moving faster than this critical speed will radiate waves.

CRediT authorship contribution statement
As explained in Subsection 2.2 , the radial displacement at the loading point, excited by a uniformly moving oscillatory load (i.e., f 0 = 0), is a key element for the instability analysis (see Eqs. (16) and (17) ). The second aim of this appendix is therefore to find the requirements that need to be met to get converged steady-state responses. Essentially, these requirements need to be defined based on U 0 ( , V ) in Eq. (16) , as that quantity is used to obtain the equivalent stiffness ( Eq. (17) ).
For illustration purposes, we consider the three cases of V = 30 m / s , V = 75 m / s and V = 150 m / s with a load frequency f 0 = 5 Hz shown in Fig. C.2 . In this figure, small frequency and time windows are shown in order to present clear features of the spectra and time-domain responses. One observes that the amplitude spectra of displacements become wider compared with the case of f 0 = 0 (compare Figs. C.1 (a) and C.2 (a)) and are spread around f = f 0 . In the time-domain responses, oscillatory patterns are observed even for the sub-critical case due to the oscillation of the load. In addition, the Doppler effect is observed in Fig. C.2 (b) [10,40] , which implies that waves are generated at frequencies different from that of the load; these frequencies are usually found from the intersections of the kinematic invariant (e.g., line 3 or 4 in Fig. 10 ) and the dispersion curves.
Based on the amplitude spectra and time-domain responses for different velocities and loading frequencies, we found the requirements (in terms of f, f max and ( N s , N r )) to obtain converged results, and they are shown in Table C.1 . Note that N shell = 20 and N load = 20 were sufficient for all the computed cases.