Research on Dynamic Response Characteristics for Basement Structure of Heavy Haul Railway Tunnel with Defects

Abstract

Based on the basic principle of thermodynamics, an elastoplastic damage constitutive model of concrete is constructed in this paper. The model is realized and verified in FLAC3D, which provides a solid foundation for the study of dynamic response and fatigue damage to the base structure of a heavy haul railway tunnel. The dynamic response and damage distribution of the base structure of a heavy-duty railway tunnel with defects were numerically simulated by the concrete elastic-plastic damage constitutive model. Then, by analyzing the response characteristics of the tunnel basement structure under different surrounding rock softening degrees, different foundation suspension range and different foundation structure damage degree are determined. The results show the following: (1) The elastoplastic damage constitutive model of concrete can well describe the stress–strain relationship of materials, especially with the simulation results of post peak softening being in good agreement with the test results, and the simulation effect of the unloading–reloading process of the cyclic loading and unloading test also meet the requirements. (2) The initial stress field and dynamic response of the tunnel basement structure under the action of train vibration load are very different from the ideal state of the structure design when the surrounding rock of the base is softened, the base is suspended, or the basement structure is damaged. With the surrounding rock softening, basement hanging, or basement structure damage developing to a certain extent, the basement structure will be damaged. (3) The horizontal dynamic stress amplitude increases with the increase in the softening degree of the basement surrounding rock. The horizontal dynamic stress of the measuring point increases with the increase in the width of the hanging out area when the hanging out area is located directly below the loading line. When the degree of damage to the basement structure is aggravated, the horizontal dynamic tensile stress of each measuring point gradually decreases. (4) The maximum principal stress increment increases with the increase in the fracture degree of the basement structure, while the minimum principal stress increment decreases with the increase in the fracture degree of the basement structure, but the variation range of the large and minimum principal stress increments is small. The research results have important theoretical and practical significance for further analysis of the damage mechanism and control technology of the foundation structure of a heavy haul railway tunnel with defects.

1. Introduction

Through the analysis of the field investigation results of tunnel defects at home and abroad, it was found that the basement structure of the tunnel is the most common site of tunnel structure defects. After long operation periods, many tunnel basement structures will show damage, subsidence, and other defects, and the track irregularity will be aggravated, even leading to rail fracture in serious cases, which poses a threat to the safety of train operation. At present, the academic understanding of the causes of the defects in the tunnel basement structure is as follows: under the coupling effect of train vibration and groundwater and other external environmental factors, the bedrock at the bottom of the tunnel softens and breaks into fine particles, which are carried away by groundwater erosion. The bedrock gradually forms cavities, and the contact state between the tunnel basement structure and the bedrock gradually deteriorates, which eventually leads to the deterioration of the tunnel basement structure. The defect is now present. It can be seen that the condition of the tunnel basement is one of the key factors affecting the service condition of its structure [1,2,3,4,5,6,7,8,9,10,11,12,13,14].

In the analysis of the impact of train load on dangerous roads, subgrades, and the surrounding environment, field test technology is the most direct, obvious, and reliable research method. The early vibration prediction is based on a large number of field test results, and an empirical formula is devised to analyze it. Koch et al. [15] carried out a field test and analysis of the vibration and noise caused by subway operation in view of its influence on surface buildings and obtained a mathematical expression of the variation of noise wave velocity in a tunnel with the thickness of lining shed at a driving speed of 60 km/h. The research team of Zach and Rutishauser analyzed the vibration frequency and acceleration characteristics of the tunnel structure through field tests [16]. Based on the European Union scientific research project, Degrende conducted a series of field tests on the vibration of 35 test trains during operation of the London Metro at a speed of 20–50 km/h [17]. The test objects included train, rail, tunnel side wall, team invert, surrounding soil layer, surface, and two buildings 70 m away from the tunnel, and the time history curve of speed or acceleration at each measuring point was recorded. The test results showed that the dynamic responses of the rail, tunnel side wall, and invert increased with the increase in driving speed.

The field testing of tunnel vibration under train load was started late by domestic scholars. Pan and Xie obtained the dynamic response characteristics of several points on the tunnel lining gun structure through field tests of the Beijing Metro. Combined with the track acceleration data and the vibration of the vehicle system, the mathematical expression of train load was obtained [18]. Zhang and Bai obtained the train vibration signals through a field test and analyzed them. It was found that these signals could be regarded as a kind of periodic vibration combination [19]. Furthermore, the mathematical expression of track vibration acceleration caused by a train passing was obtained through Fourier transform spectrum analysis. Li and Gao obtained the acceleration of the concrete bottom layer 0.55 m below the rail surface through a field test of train vibration in the Jinjiayan Tunnel and pointed out that the vibration caused by the train load occurred on the lining, which can cause defects in the bottom structure [20]. Wang and Yang conducted a field test on the vibration of the Zhuting Tunnel structure on the Beijing to Guangzhou Railway and analyzed the measured data by power spectral density [21]. On this basis, the wheel rail vibration analysis model was established, and the mathematical expression of train vibration load was obtained. Xue et al. obtained the vertical vibration velocity of the tunnel invert and foundation soil as well as the attenuation law of soil pressure along the depth by simulating the cyclic dynamic load test of high-speed train vibration in a loess tunnel and further analyzed the variation law of foundation soil pressure and excess hydrostatic pressure with vibration frequency [22]. Peng et al. analyzed the dynamic response of the bottom structure of the Wushiling Tunnel under train load through a field test of vertical static and dynamic stress and acceleration [23]. Shi et al. conducted a long-term test and analysis of the dangerous deformation of a bottom structure and the stress of invert concrete and reinforcement of the Xichongzi Tunnel of the Wuhan Guangzhou high-speed railway and concluded that the tunnel bottom structure experienced the repeated change process of uplift subsidence, and the stress of invert concrete and reinforcement experienced the repeated change process of tension compression. In the process of tunnel construction, the invert is in an unfavorable stress state, so it should be closed as soon as possible to ensure operation safety [24]. Fu et al., through a field test of the Sanjiacun Tunnel on the Shuozhou to Huanghua Railway, statistically analyzed the dynamic stress of the filling layer in different surrounding rock sections, converted the dynamic force into static value using the soil column method, and obtained the impact coefficient of the base structure. Based on this, the dynamic response and impact coefficient of the basement structure of the dangerous track under the action of a heavy-duty train were analyzed by establishing a comprehensive analysis model of train tunnel track structure [25]. Yu et al. studied the virtual hybrid simulation method for underground structures subjected to seismic loadings [26,27]. Based on the Least Square Method, some scholars have analyzed different cases [28,29,30,31,32,33,34,35,36,37,38,39,40,41].

Based on the above research, the theoretical analysis of the vibration response of tunnel structure under train load can only analyze some simplified problems, and there are many assumptions, which is difficult to reflect the real state of the analysis object.

The reasonable description of concrete damage evolution process under the static and dynamic coupling of surrounding rock pressure and train load is the key to effectively evaluate the base condition of heavy haul railway tunnel, and it is the basis to reasonably and accurately reflect the health state of base structure of heavy haul railway tunnel. In order to reasonably and accurately reflect the dynamic response characteristics of the base structure of the defective heavy haul railway tunnel, the following studies are carried out in this paper: Section 1 introduces the previous research results; Section 2 establishes the elastic-plastic damage model based on the previous research; Section 3 verifies the rationality, correctness, and effectiveness of the elastic-plastic damage model by numerical tests; Section 4 carries out the dynamic response simulation test of the foundation structure of a heavy-duty railway tunnel by using the elastic-plastic damage model; Section 5 analyzes the test data of dynamic response simulation test of foundation structure of a heavy haul railway tunnel; Section 6 draws the research conclusions and further observations.

2. Establishment of Elastoplastic Damage Model

The tensile damage variable and shear damage variable are used to phenomenally reflect the influence of micro damage on the degradation of macro strength and stiffness of concrete. The elastic Helmholtz free energy of material is defined in the effective stress space; the development and evolution processes of material damage and plastic deformation are determined, and the elastic-plastic damage constitutive model in accordance with the basic principle of thermodynamics is constructed based on the non-associated flow law.

2.1. Lastoplastic Damage Constitutive Relation

(1)

Stress–strain relationship

According to the plastic increment theory, the total strain tensor of materials can generally be divided into elastic strain tensor ${\epsilon }^e$ and plastic strain tensor ${\epsilon }^p$.

$$\epsilon ={\epsilon }^e+{\epsilon }^p$$

(1)

The elastic part can be obtained by Hooke’s law.

$$\sigma =C{\epsilon }^e$$

(2)

where $\sigma$ is Cauchy stress tensor;

$C$ is the elastic stiffness tensor.

According to Lemaitre’s equivalent stress assumption, the elastoplastic constitutive relationship of nondestructive materials in effective stress space can be defined.

$$\overline{\sigma }=C_0\left(\epsilon -{\epsilon }^p\right)$$

(3)

where $\overline{\sigma }$ is the effective stress tensor; $C_0$ is the initial elastic stiffness tensor.

In order to clearly distinguish the different effects of tensile and compressive stress on materials, the effective stress tensor $\overline{\sigma }$ is decomposed into positive stress ${\overline{\sigma }}^{+}$ and negative stress ${\overline{\sigma }}^{-}$.

$$\overline{\sigma }={\overline{\sigma }}^{+}+{\overline{\sigma }}^{-}$$

(4)

$${\overline{\sigma }}^{+}={\displaystyle \sum _i\langle {\overline{\sigma }}_i\rangle }\left(p_i\otimes p_i\right)$$

(5)

where $\langle x\rangle =\{\begin{array}{cc}0& x<0\\ x& x\ge 0\end{array}$ is the Macaulay Function;

$p_i$ is the Eigenvector of the $\overline{\sigma }$.

(2)

Elastic Helmholtz free energy

Faria et al. proposed the elastic Helmholtz free energy function based on the positive and negative decomposition of effective stress in the effective stress space.

$$\psi =\left(1-D^{+}\right){\psi }_0^{+}+\left(1-D^{-}\right){\psi }_0^{-}=\frac{1}{2}\left(1-D^{+}\right){\overline{\sigma }}^{+}:C_0^{-1}:{\overline{\sigma }}^{+}+\frac{1}{2}\left(1-D^{-}\right){\overline{\sigma }}^{-}:C_0^{-1}:{\overline{\sigma }}^{-}$$

(6)

Inside:

$$C_0^{-1}=\frac{1}{E}\left[\frac{1+\upsilon }{2}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)-\upsilon {\delta }_{ij}{\delta }_{kl}\right]$$

(7)

where $E$ is the elastic modulus;

$\upsilon$ is the Poisson’s ratio;

${\delta }_{ij}$ is the Kronecker symbol.

Then, it is deduced that

$$\sigma =\frac{\partial \psi }{\partial \epsilon }=\left(1-D^{+}\right){\overline{\sigma }}^{+}+\left(1-D^{-}\right){\overline{\sigma }}^{-}$$

(8)

(3)

Damage criterion

The tensile damage energy release rate $Z^{+}$ and shear damage energy release rate $Z^{-}$ corresponding to the tensile e damage variable $D^{+}$ and compressive damage variable $D^{-}$ can be expressed as

$$Z^{+}=-\frac{\partial {\psi }^{+}}{\partial D^{+}}={\psi }_0^{+}$$

(9)

$$Z^{-}=-\frac{\partial {\psi }^{-}}{\partial D^{-}}={\psi }_0^{-}$$

(10)

According to the assumption of damage energy release rate by Li et al. [42], this paper adopts a similar definition of damage energy release rate.

$${\overline{\tau }}^{+}=g\left(Z^{+}\right)=\sqrt{E_0\left({\overline{\sigma }}^{+}:C_0^{-1}:{\overline{\sigma }}^{+}\right)}$$

(11)

$${\overline{\tau }}^{-}=\alpha {\overline{I}}_1+\sqrt{{\overline{J}}_2}$$

(12)

where ${\overline{\tau }}^{+}$ and ${\overline{\tau }}^{-}$ are tensile and compressive forces, respectively;

$E_0$ is the initial Young’s modulus of the material.

Based on the above definition of damage energy release rate, the following damage criteria are established

$${\overline{\tau }}^{+}-r^{+}\le 0$$

(13)

$${\overline{\tau }}^{-}-r^{-}\le 0$$

(14)

where $r^{+}$ and $r^{-}$ are the threshold values of tensile and compressive forces, respectively, and the initial values are $r_0^{+}$ and $r_0^{-}$, respectively, controlling the range of linear elastic region, and the relationship between tensile and compressive strength is

$$r_0^{+}=f_0^{+}$$

(15)

$$r_0^{-}=\left(\frac{\sqrt{3}}{3}-\alpha \right)f_0^{-}$$

(16)

where $f_0^{+}$ and $f_0^{-}$ are the linear upper limit strength under uniaxial tension and compression, which are generally taken according to the following two formulas

$$f_0^{+}=\left(0.6~0.8\right)f_t$$

(17)

$$f_0^{-}=\left(0.3~0.5\right)f_c$$

(18)

where $f_t$ and $f_c$ are uniaxial tensile strength and uniaxial compressive strength, respectively.

2.2. Internal Variable Evolution Process

(1)

Damage variable

Based on the above damage criteria and according to the orthogonal flow law, the damage change rate can be obtained by the following formula

$${\dot{D}}^{\pm }={\dot{r}}^{\pm }\frac{\partial G^{\pm }\left(r^{\pm }\right)}{\partial r^{\pm }}$$

(19)

The damage loading and unloading can be judged by Kuhn–Tucker Relationship

$${\dot{r}}^{\pm }\ge 0G^{\pm }\le 0{\dot{r}}^{\pm }{G}^{\pm }=0$$

(20)

When $G^{\pm }\le 0$, ${\dot{r}}^{\pm }=0$ then $D^{\pm }=0$, it is in the stage of damage unloading or neutral load change, and the material damage will no longer increase. When ${\dot{r}}^{\pm }\ge 0$, $G^{\pm }=0$, ${\dot{r}}^{\pm }={\overline{\tau }}^{\pm }\ge 0$, get Formula (21).

$${\dot{D}}^{\pm }={\dot{G}}^{\pm }\left(r^{\pm }\right)\ge 0$$

(21)

According to the following constraints of tensile and shear damage variables, both tensile and shear damage variables adopt the expressions proposed by Mazars et al.

$$0\le G^{\pm }\left(r^{\pm }\right)\le 1{\dot{G}}^{\pm }\left(r^{\pm }\right)\ge 0G^{\pm }\left(r_0^{\pm }\right)=0$$

(22)

$$D^{+}=1-\frac{r_0^{+}\left(1-A^{+}\right)}{r^{+}}-A^{+}\exp \left[B^{+}\left(1-\frac{r^{+}}{r_0^{+}}\right)\right]$$

(23)

$$D^{-}=1-\frac{r_0^{-}\left(1-A^{-}\right)}{r^{-}}-A^{-}\exp \left[B^{-}\left(1-\frac{r^{-}}{r_0^{-}}\right)\right]$$

(24)

where $A^{+}$, $B^{+}$, $A^{-}$ and $B^{-}$ are constitutive model parameters, $A^{+}$ and $B^{+}$ can be calibrated by uniaxial tensile test, $A^{-}$ and $B^{-}$ can be calibrated by uniaxial compression test.

(2)

Plastic deformation

The evolution process of material plastic deformation is described by the Empirical Formula [43] proposed by Wu.

$${\dot{\epsilon }}^p={\xi }^p{E}_{0}H\left({\dot{D}}^{-}\right)\frac{\langle {\epsilon }^e:\dot{\epsilon }\rangle }{\overline{\sigma }:\overline{\sigma }}:\overline{\sigma }$$

(25)

where ${\xi }^p$ is the empirical parameter of the model, which can be roughly calibrated by uniaxial compression test. $H(•)$ is the Heaviside step function.

3. Numerical Test Verification of Elastoplastic Damage Model

The reasonable description of the damage evolution process in concrete under the static and dynamic coupling of surrounding rock pressure and train load was the core issue in the establishment of the concrete elastic–plastic damage model in this chapter. It was the key to effectively evaluating the base condition of the heavy haul railway tunnel and the basis for reasonably and accurately reflecting the health state of the base structure of the heavy haul railway tunnel.

Using the elastic-plastic damage model established in Section 2, the monotonic loading and cyclic loading and unloading tests were numerically simulated in FLAC3D. The loading mode was consistent with the prototype test. The numerical simulation results and test results were compared to verify the correctness of the constitutive model and the effectiveness of the model program. However, because the path dependence of unloading stiffness is not considered in this model, the hysteretic characteristics after unloading and reloading cannot be reproduced.

3.1. Uniaxial Tensile Test

The uniaxial tensile test simulated the Gopalarantam–Shah test [42]. The initial Young’s modulus $E_0$ measured in the literature was 3.1 × 10⁴ MPa, and the tensile strength $f_0^{+}$ was 3.48 MPa. The numerical simulation parameters are shown in

Table 1. Simulation parameters of the uniaxial tensile test.

${\mathit{E}}_0/\mathbf{MPa}$	${\mathit{f}}_0^{+}/\mathbf{MPa}$	${\mathit{A}}^{+}$	${\mathit{B}}^{+}$	${\mathit{\xi }}^{\mathit{p}}$
3.1 × 104	3.48	1.0	0.518	0.52

. The comparison between the numerical simulation results and test results is shown in Figure 1. It can be seen that the model in this paper better described the whole process curve of stress–strain in concrete materials under a uniaxial tensile stress state, including the softening section.

3.2. Uniaxial Compression Test

The uniaxial compression test simulated the uniaxial compression test of Karsan and Jirsa [43]. The initial Young’s modulus $E_0$ measured was 3.17 × 10⁴ MPa, and the compressive strength $f_0^{-}$ was 27.6 MPa. The numerical simulation parameters are shown in

Table 2. Simulation parameters of uniaxial compression test.

${\mathit{E}}_0/\mathbf{MPa}$	${\mathit{f}}_0^{-}/\mathbf{MPa}$	${\mathit{A}}^{-}$	${\mathit{B}}^{-}$	${\mathit{\xi }}^{\mathit{p}}$
3.17 × 104	10.2	1.0	0.16	0.42

. The comparison between the numerical simulation results and test results is shown in Figure 2. The model in this paper also well simulated the stress–strain curve, including the softening section under a uniaxial compression stress state, and the program operation results were satisfactory.

3.3. Biaxial Stress Test

The biaxial stress test simulated the test done by Kurfer et al. [44]. In their biaxial stress test, the measured Young’s modulus was 3.10 × 10⁴ MPa and the compressive strength was 32.0 MPa. The numerical simulation parameters are shown in

Table 3. Simulation parameters of biaxial stress test.

${\mathit{E}}_0/\mathbf{MPa}$	${\mathit{f}}_0^{-}/\mathbf{MPa}$	${\mathit{A}}^{-}$	${\mathit{B}}^{-}$	${\mathit{\xi }}^{\mathit{p}}$
3.10 × 104	15.0	1.0	0.19	0.37

. The working conditions where the stress ratios of σ₁/σ₂ comprised −1:0 and −1:−1 were simulated, respectively. The comparison between the numerical calculation results and the test results is shown in Figure 3. It can be seen that the simulation effect of the model for the material stress–strain relationship under biaxial stress also met the requirements.

3.4. Cyclic Tensile Test

The cyclic tensile test simulated Taylor’s test [45]. In Taylor’s uniaxial repeated tensile test, the measured Young’s modulus was 3.17 × 10⁴ MPa and the tensile strength was 3.47 MPa. The numerical simulation parameters are shown in

Table 4. Simulation parameters of the cyclic tensile test.

${\mathit{E}}_0/\mathbf{MPa}$	${\mathit{f}}_0^{+}/\mathbf{MPa}$	${\mathit{A}}^{+}$	${\mathit{B}}^{+}$	${\mathit{\xi }}^{\mathit{p}}$
3.17 × 104	3.47	1.0	0.518	0.52

. The comparison between the numerical calculation results and test results is shown in Figure 4. It can be seen that the model not only well described the post peak strength softening stage, but also effectively simulated the elastic unloading–reloading process.

3.5. Cyclic Compression Test

The cyclic tensile test simulated the cyclic compression test conducted by Karsan and Jirsa [43]. The measured elastic modulus was 3.0 × 10⁴ MPa, and the compressive strength was 27.6 MPa. The numerical simulation parameters are shown in

Table 5. Simulation parameters of cyclic compression test.

${\mathit{E}}_0/\mathbf{MPa}$	${\mathit{f}}_0^{-}/\mathbf{MPa}$	${\mathit{A}}^{-}$	${\mathit{B}}^{-}$	${\mathit{\xi }}^{\mathit{p}}$
3.0 × 104	13.8	1.0	0.16	0.42

. The comparison between numerical calculation results and test results is shown in Figure 5. The model had a good simulation effect on the loading process of cyclic compression test. The elastic unloading–reloading process shown by the numerical calculation was slightly different from the test results and could not describe the hysteretic phenomenon in the unloading–reloading process in the test.

4. Simulation Experiment on Dynamic Response of the Foundation Structure of a Heavy Haul Railway Tunnel

4.1. Introduction of Finite Element Calculation Model

The three-dimensional dynamic numerical analysis model was established to simulate the response characteristics of the tunnel basement structure under dynamic load and to analyze the dynamic stress and damage distribution of each part of the tunnel basement structure.

4.1.1. Numerical Model

As shown in Figure 6, the model size is 80 m × 1 m × 80 m, and the unit thickness is taken along the tunnel axis. The tunnel axis is located in the center of the model, and the length and width of the model are both 80 m. The developed elastoplastic damage model is used for the surrounding rock, invert, filling layer, primary support, secondary lining, and other tunnel structures.

4.1.2. Calculation Condition

The standard cross-section of railway double line composite lining tunnel is taken as the typical calculation cross-section. As the measured field data revealed that the running speed of a heavy train has little influence on the dynamic stress of the tunnel foundation structure, the numerical analysis mainly focused on the influence of different design parameters of foundation structure and different conditions of the tunnel foundation on the dynamic response and damage characteristics of the tunnel foundation structure. The calculation conditions under different design parameters of the tunnel basement structure are shown in

Table 6. Calculation cases with different design parameters of the tunnel base structure.

Influence Factor	Grade of Surrounding Rock	Thickness Invert	Rise Span Ratio of Invert
Specific parameters	III	30 cm	1/15
	IV	40 cm	1/12
	V	50 cm	1/10
	VI	60 cm	1/8

.

The benchmark calculation example of the above conditions was as follows: the axle load of the train is 30 t, the surrounding rock grade is grade V, the invert thickness is 50 cm, the rise span ratio is 1/12, and the train load acts on the side of the heavy train line. Other conditions only need to change the corresponding parameters.

4.1.3. Calculation Parameters

The basic physical and mechanical parameters of surrounding rock were adopted as shown in

Table 7. Basic mechanical parameters of concrete and surrounding rock.

Structure	Elasticity Modulus/GPa	Poisson’s Ratio	Weight/(kN·m−3)	Cohesion/kPa	Friction Angle/°
Ballast	32.0	0.18	25.0	/	/
Filling layer	29.0	0.18	25.0	/	/
Invert	30.0	0.18	25.0	/	/
Second lining	30.0	0.18	25.0	/	/
First branch	27.5	0.18	25.0	/	/
Surrounding rock	1.5	0.36	19.5	200	42

. The mechanical damage parameters of concrete are shown in

Table 8. Mechanical damage parameters of concrete.

Structure	${\mathit{f}}_0^{-}/\mathbf{MPa}$	${\mathit{f}}_0^{+}/\mathbf{MPa}$	${\mathit{A}}^{-}$	${\mathit{B}}^{-}$	${\mathit{A}}^{+}$	${\mathit{B}}^{+}$	${\mathit{\xi }}^{\mathit{p}}$
Ballast	8.44	1.26	1.15	0.23	1.0	0.58	0.35
Filling layer	4.76	0.89	0.98	0.16	1.0	0.51	0.41
Invert	6.68	1.10	1.05	0.18	1.0	0.54	0.39
Second lining	6.68	1.10	1.05	0.18	1.0	0.54	0.39
First branch	3.84	0.77	0.94	0.14	1.0	0.46	0.43

.

4.1.4. Simulation Method of Train Dynamic Load

Under the condition of a 30 t axle load and 80 km/h running speed, the Single Side Static Wheel Weight P₀ = 150 kn and the Unsprung Mass M₀ = 1200 kg were adopted. The time history curve of the artificial excitation force as shown in Figure 7 was obtained by Formula (26). The load acts on the rail simulated by Beam element.

$$F(t)=k_1{k}_2(P_0+P_1\sin {\omega }_{1}t+P_2\sin {\omega }_{2}t+P_3\sin {\omega }_{3}t)$$

(26)

where $P_0$ is the static load of the wheel. $P_1$, $P_2$, $P_3$ are typical vibration load values related to ride irregularity, dynamic additional load, and wave abrasion, respectively, corresponding to a typical value of control condition ①–③ in Table 4. $k_1$ is the wheel rail force superposition coefficient. $k_2$ is the rail dispersion transfer coefficient.

The wavelength management values of track geometric irregularity correspond to ①, ②, and ③ in

Table 9. Track geometry irregularity management value in the UK.

Control Conditions	Wavelength/m	Normal Vector/mm
According to ride comfort ①	50.00	16.00
	20.00	9.00
	10.00	5.00
According to the dynamic additional load acting on the line ②	5.00	2.50
	2.00	0.60
	1.00	0.30
Wave abrasion ③	0.50	0.10
	0.05	0.005

.

4.1.5. Boundary Conditions

The setting of artificial boundary conditions and vibration wave input are the key links to realize dynamic numerical analysis, which directly affect the accuracy and credibility of the calculation results. In the engineering rock mass in semi-infinite space, when the vibration wave propagates around, it will transmit and refract when it meets the structure or rock mass structural plane and interfere with the incident wave to form a complex wave field. In the near-field dynamic numerical analysis, in order to effectively simulate the stray scattered wave passing through the truncated boundary of the model and make the vibration wave energy escape to infinity, it is necessary to introduce artificial boundary conditions to eliminate the reflection and oscillation of vibration wave on the artificial cutting boundary of the calculation area. In the FLAC3D dynamic analysis of this chapter, the free field boundary is used as the boundary condition of dynamic analysis.

In the numerical simulation of the vibration response of near surface structures such as slopes and shallow buried underground structures, the free field motion without structures must be considered at the lateral boundary of the numerical model. In the numerical simulation, the free field motion is applied to the model boundary, so that the artificial boundary remains non reflective. The lateral boundary of the discrete model simulates the viscous boundary through the coupling of damper and free field, and the unbalanced force on the free field node is applied to the boundary node of the discrete model. The unbalanced force on a plane parallel to the axis in the normal direction can be expressed as

$$\{\begin{array}{l}{F}_x=-\rho C_p\left(v_x^m-v_x^{ff}\right)A+F_x^{ff}\\ F_y=-\rho C_s\left(v_y^m-v_y^{ff}\right)A+F_y^{ff}\\ F_z=-\rho C_s\left(v_z^m-v_z^{ff}\right)A+F_z^{ff}\end{array}$$

(27)

where $\rho$ is the material medium density. $C_p$ and $C_s$ represent the wave velocities of P-wave and S-wave in material medium, respectively; $A$ is the influence area of free field nodes. $v_x^m$, $v_y^m$, $v_z^m$, $v_x^{ff}$, $v_y^{ff}$, $v_z^{ff}$ represent the three components along the coordinate axis of the particle vibration velocity caused by fluctuation on the main node on the artificial boundary and the corresponding free field node respectively. $F_x^{ff}$, $F_y^{ff}$, $F_z^{ff}$ represent the contribution of normal stress and tangential stress of free field node influence element to node force, respectively.

If the calculation model is a homogeneous medium without near surface structures, the calculation area is consistent with the motion law of free field, and the lateral damper does not work. When the two move relatively, the damper absorbs the energy of the external traveling wave. In FLAC3D, the free field boundary condition requires that the bottom surface of the model be a horizontal plane and the normal direction be an axial direction. The side must be vertical and the normal direction axis or axis direction.

4.1.6. Layout of Measuring Points

The arrangement of measuring points is shown in Figure 8. Measuring points 1 and 5 are respectively located at the joint of side invert and side wall of the loaded line and the empty line, and measuring points 2, 3, and 4 are respectively located at the right rail of the loaded line, the center of the tunnel, and just below the center of the empty line.

5. Analysis of Simulation Experiment Test Data

5.1. Horizontal Dynamic Stress

5.1.1. Influence of Rock Softening on Basement

The amplitude of the horizontal dynamic stress at measuring points 2, 3, and 4 under different softening degrees of the basement surrounding rock is shown in

Table 10. Horizontal dynamic stress amplitudes of measuring points 2, 3, and 4 under different softening degrees of the basement surrounding rock.

Measuring Point	Horizontal Dynamic Stress Amplitude/kPa under Different Softening Degrees of the Basement Surrounding Rock
	30% Reduction	50% Reduction	80% Reduction	90% Reduction
2	−177.9	−186.6	−214.6	−240.8
3	−93.6	−101.7	−117.4	−139.1
4	−37.6	−38.7	−47.9	−58.8

and Figure 9.

Under different softening degrees of the basement surrounding rock, the horizontal dynamic stress of measuring points 2, 3, and 4 was tensile stress. The results showed that the horizontal dynamic stress of the measuring point under the heavy vehicle line (measuring point 2) was the largest, the measuring point under the tunnel center (measuring point 3) was the second largest, and the measuring point under the empty vehicle line (measuring point 4) was the smallest; the horizontal dynamic stress amplitude increased with the increase in the softening degree of the basement surrounding rock, and the horizontal dynamic stress amplitudes of the measuring points 2, 3, and 4 were −240.5, −138.4, and −59.3 kPa, respectively, when the mechanical parameters of the basement surrounding rock were reduced by 90%.

5.1.2. The Influence of Base Hanging

The time history curves of the horizontal dynamic stress of measuring points 2, 3, and 4 under different basement hanging ranges are shown in Figure 3, and the amplitudes of horizontal dynamic stress at measuring points 2, 3, and 4 under different basement hanging ranges are shown in

Table 11. Horizontal dynamic stress amplitudes of measuring points 2, 3, and 4 under different basement hanging ranges.

Measuring Point	Horizontal Dynamic Stress Amplitude/kPa under Different Basement Hanging Ranges
	1/4 Width Loading Line	1/2 Width Empty Line	1/2 Width Loading Line	3/4 Width Loading Line
2	−211.8	−184.9	−264.3	−332.2
3	−93.7	−123.3	−119.8	−216.5
4	−37.6	−82.5	−31.6	−55.3

and Figure 10.

The results showed that the horizontal dynamic stress of the measuring point increased with the increase in the width of the hanging out area when the hanging out area was located directly below the loading line; when the hanging out area was located at different positions on the base but with the same width, the horizontal dynamic stress of the measuring point was quite different.

5.1.3. Influence of Basement Structure Damage

The amplitudes of horizontal dynamic stress at measuring points 2, 3, and 4 under different degrees of damage to the basement structure are shown in

Table 12. Horizontal dynamic stress amplitudes of measuring points 2, 3, and 4 under different degrees of damage to the foundation structure.

Measuring Point	Horizontal Dynamic Stress Amplitude/kPa under Different Degrees of Damage to the Basement Structure
	10% Reduction	20% Reduction	50% Reduction	80% Reduction
2	−147.8	−138.3	−106.8	−51.2
3	−76.4	−71.9	−47.9	−16.7
4	−31.2	−27.5	−17.7	−6.2

and Figure 11.

The results showed that the horizontal dynamic tensile stress of measuring points 2, 3, and 4 was also the tensile stress under different degrees of damage to the basement structure. When the degree of damage to the basement structure was aggravated, the horizontal dynamic tensile stress of each measuring point gradually decreased.

5.2. Principal Stress

5.2.1. Influence of Rock Softening on Basement

The maximum amplitudes of principal stress at measuring points 1, 2, 3, and 5 under different softening degrees of the basement surrounding rock are shown in

Table 13. Maximum principal stress amplitude of measuring points 1, 2, 3, and 5 under different softening degrees of basement surrounding rock.

Measuring Point	Maximum Principal Stress Amplitude/kPa under Different Softening Degrees of Basement Surrounding Rock
	Maximum Principal Stress				Maximum Principal Stress Increment
	30%	50%	80%	90%	30%	50%	80%	90%
1	421	387	346	329	−2.7	−1.9	1.4	2.1
2	109	138	162	168	35.8	31.8	28.5	24.9
3	111	96	98	103	21.6	20.2	18.8	16.9
5	437	386	341	321	0.5	−1.2	−1.1	−0.6

, Figure 12 and Figure 13.

The minimum amplitudes of principal stress at measuring points 1, 2, 3, and 5 under different softening degrees of the basement surrounding rock are shown in

Table 14. Minimum principal stress amplitude of measuring points 1, 2, 3, and 5 under different softening degrees of basement surrounding rock.

Measuring Point	Minimum Principal Stress Amplitude/kPa under Different Softening Degrees of Basement Surrounding Rock
	Minimum Principal Stress				Minimum Principal Stress Increment
	30%	50%	80%	90%	30%	50%	80%	90%
1	1706	1615	1471	1407	−55.5	−52.3	−49.2	−53.6
2	1478	1748	2206	2342	−171.2	−187.6	−215.8	−240.1
3	1382	1576	2047	2226	−95.3	−98.4	−117.9	−136.2
5	1727	1639	1498	1425	−11.5	−10.1	−8.4	−6.5

, Figure 14 and Figure 15.

The results showed that when the mechanical parameters of bedrock were reduced by 30, 50, 80, and 90%, the maximum principal stress and minimum principal stress of each measuring point were compressive stress, that is, after the train load, the position of each measuring point was still in a compression state; the increment of principal stress, especially the increment of minimum principal stress, was smaller than the initial value of principal stress; the maximum principal stress and minimum principal stress of measuring points 2 and 3 increased with the increase in softening degree of the basement surrounding rock. The large and minimum principal stresses of measuring points 1 and 5 decreased with the increase in the softening degree of the surrounding rock. The influence of dynamic load on the stress of measuring points 2 (directly below the heavy vehicle line) and 3 (directly below the tunnel center) was greater than that of measuring points 1 and 5 at the junction of the side wall and invert.

5.2.2. Influence of Basement Hanging Ranges

The maximum amplitudes of principal stress measuring points 1, 2, 3, and 5 under different basement hanging ranges are shown in

Table 15. Maximum principal stress amplitude of measuring points 1, 2, 3, and 5 under different basement hanging ranges.

Measuring Point	Maximum Principal Stress Amplitude/kPa under Different Basement Hanging Ranges
	Maximum Principal Stress				Maximum Principal Stress Increment
	1/4 Width Loading Line	1/2 Width Empty Line	1/2 Width Loading Line	3/4 Width Loading Line	1/4 Width Loading Line	1/2 Width Empty Line	1/2 Width Loading Line	3/4 Width Loading Line
1	415	436	357	346	−2.2	−6.4	4.6	3.9
2	11	212	−9	−16	8.7	40.3	8.9	−20.5
3	246	508	286	−17	35.2	−17.6	3.7	−1.3
5	437	353	435	287	−0.8	4.5	−7.4	−31.6

, Figure 16 and Figure 17.

The minimum amplitudes of principal stress measuring points 1, 2, 3, and 5 under different basement hanging ranges are shown in

Table 16. Minimum principal stress amplitude of measuring points 1, 2, 3, and 5 under different basement hanging ranges.

Measuring Point	Minimum Principal Stress Amplitude/kPa under Different Basement Hanging Ranges
	Minimum Principal Stress				Minimum Principal Stress Increment
	1/4 Width Loading Line	1/2 Width Empty Line	1/2 Width Loading Line	3/4 Width Loading Line	1/4 Width Loading Line	1/2 Width Empty Line	1/2 Width Loading Line	3/4 Width Loading Line
1	1638	2009	1259	1279	−58.1	−57.2	14.6	19.2
2	1590	2106	1397	176	−206.3	−186.4	−265.8	−296.7
3	1547	1947	2332	492	−94.5	−135.3	−94.9	−198.4
5	1819	1235	2041	2154	−11.8	−18.2	−5.8	13.3

, Figure 18 and Figure 19.

The results showed that the main stress amplitude of each measuring point in the basement hanging empty range did not exceed the tensile strength and compressive strength of concrete material. When the hanging empty area was directly below the heavy vehicle line, the main stress amplitude of measuring point 2 changed from 12 to −11 kPa when the hanging empty area width increased from 1/4 to 1/2, and the main stress amplitude of measuring point 3 changed from 293 to −18 kPa when the hanging empty area width increased from 1/2 to 3/4. The results showed that the minimum principal stress amplitudes of measuring points 2 and 3 decreased sharply from 1408 and 2344 kPa to 182 and 500 kPa, respectively; when the hanging out area was located at different positions of the tunnel bottom and the hanging out width was the same, the stress of each measuring point was very different; when the hanging out area was located under the heavy vehicle line and the width was 1/2, the minimum principal stress amplitudes of measuring points 1 and 5 were 2020 and 1243 kPa, respectively, while the hanging out area was located under the empty vehicle line and the width was 1/2. The minimum principal stresses of measuring points 1 and 5 were 1276 and 2054 kPa, respectively.

5.2.3. Influence of Basement Structure Damage

The maximum principal stress amplitudes at measuring points 1, 2, 3, and 5 under different degrees of damage to the basement structure are shown in

Table 17. Maximum principal stress amplitude of measuring points 1, 2, 3, and 5 under different degrees of damage to the basement structure.

Measuring Point	Maximum Principal Stress Amplitude/kPa under Different Degrees of Damage to the Basement Structure
	Maximum Principal Stress				Maximum Principal Stress Increment
	10%	20%	50%	80%	10%	20%	50%	80%
1	425	428	426	426	0.7	0.7	2.2	3.5
2	96	101	105	114	38.9	39.8	45.1	56.2
3	115	116	114	118	22.4	22.6	22.5	27.9
5	427	429	424	423	−2.1	−1.9	1.0	1.7

, Figure 20 and Figure 21.

The minimum principal stress amplitudes at measuring points 1, 2, 3, and 5 under different degrees of damage to the basement structure are shown in

Table 18. Minimum principal stress amplitude of measuring points 1, 2, 3, and 5 under different degrees of damage to the basement structure.

Measuring Point	Minimum Principal Stress Amplitude/kPa under Different Degrees of Damage to the Basement Structure
	Minimum Principal Stress				Minimum Principal Stress Increment
	10%	20%	50%	80%	10%	20%	50%	80%
1	1766	1764	1765	1764	−44.5	−45.6	−47.5	−47.8
2	1317	1317	1332	1361	−151.3	−148.2	−107.1	−51.1
3	1279	1279	1298	1315	−77.6	−70.2	−51.3	−16.5
5	1782	1776	1780	1769	−14.3	−11.2	−13.6	−17.9

, Figure 22 and Figure 23.

The results showed that the maximum principal stress increment increased with the increase in the fracture degree of the basement structure, while the minimum principal stress increment decreased with the increase in the fracture degree of the basement structure, but the variation range of the large and minimum principal stress increments was small.

6. Conclusions and Further Observations

6.1. Conclusions

Based on the basic principle of thermodynamics, the elastoplastic damage constitutive model of concrete is constructed in this paper. The model is realized and verified in FLAC3D, which provides a solid foundation for the study of dynamic response and fatigue damage to the base structure of a heavy haul railway tunnel. The dynamic response and damage distribution of the base structure of a heavy-duty railway tunnel with defects were numerically simulated by the concrete elastic–plastic damage constitutive model. Then, by analyzing the response characteristics of tunnel basement structure under different surrounding rock softening degree, different foundation suspension range and different foundation structure damage degree are determined. The research results have reference value for analyzing the dynamic response characteristics of the basement structure of a heavy haul railway tunnel with different defects.

The following conclusions can be made.

(1)

The elastoplastic damage constitutive model.

The elastoplastic damage constitutive model of concrete can well describe the stress–strain relationship of materials, especially the simulation results of post peak softening, which were in good agreement with the test results. Meanwhile, the simulation effect of the unloading–reloading process of the cyclic loading and unloading test also met the requirements, so as to verify the rationality of the model and the correctness of the model program.

(2)

Health status of tunnel basement structure.

The initial stress field and dynamic response of the tunnel basement structure under the action of train vibration load are very different from the ideal state of the structure design when the surrounding rock of the base is softened, suspended, or damaged. When the surrounding rock softens, basement hanging or basement structure damage develops to a certain extent, and the basement structure suffers damage.

(3)

Horizontal dynamic stress.

The horizontal dynamic stress amplitude increases with the increase in the softening degree of the basement surrounding rock. The horizontal dynamic stress of the measuring point increases with the increase in the width of the hanging out area when the hanging out area is located directly below the loading line. When the hanging out area is located at different positions of the base but has the same width, the horizontal dynamic stress of the measuring point is quite different. With the aggravation of the degree of damage to the basement structure, the horizontal dynamic tensile stress of each measuring point gradually decreases.

(4)

Principal stress

When the mechanical parameters of bedrock are reduced by 30, 50, 80, and 90%, the maximum principal stress and minimum principal stress of each measuring point are compressive stress, that is, after the train load, the position of each measuring point is still in a compression state; the increment of principal stress, especially the increment of minimum principal stress, is smaller than the initial value of principal stress. The maximum principal stress increment increases with the increase in the fracture degree of the basement structure. The minimum principal stress increment decreases with the increase in the fracture degree of the basement structure. But the variation range of the large and minimum principal stress increments is small.

6.2. Further Observations

Through theoretical analysis and numerical simulation, this paper studied the dynamic response characteristics, vibration characteristics, and damage distribution characteristics of the basement structure of a heavy haul railway tunnel with defects. Although some conclusions were drawn, there are still some imperfections in the content of this paper due to the limitations of time and conditions, and there are many problems that need further study:

(1)

This paper only performed a numerical analysis of the damage to the tunnel base structure and did not study the damage to the base surrounding rock. Mechanical properties such as strength and stiffness of the surrounding rock are weaker than that of the base structure and are more prone to damage. The elastic–plastic damage model of surrounding rock can be introduced into a later numerical analysis to study the damage to the surrounding rock.

(2)

In this paper, the dynamic response and damage distribution of the basement structure of the heavy-duty railway tunnel were analyzed only through numerical calculation. It is hoped that in future research, the damage development and defect generation mechanism of the basement of the heavy-duty railway tunnel at the mesoscale level will be analyzed with the indoor test method.