# Dynamic Response and Fatigue-Life Analysis of High-Speed Railway Tunnel Base with Defects

Authors: Zhiwei Wang; Weibin Ma
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Based on the repeated “tensile-compressive-tensile-compressive” stress characteristics of tunnel-base concrete structure, double scalar damage variables were defined using the concrete stochastic damage mechanics. The elastoplastic Helmholtz free energy was corrected by introducing the hardening parameter. A stochastic dynamic damage constitutive model with the nonlinear and strain-rate effect behavior of concrete was derived via dynamic expansion of damage energy based on the principle of energy dissipation. An “elastic prediction-plastic correction” numerical analysis algorithm was developed based on the solution of the probability density evolution equation. Secondary development of the algorithm was achieved using the Universal Distinct Element Code (UDEC) for numerical calculation of the proposed constitutive model. A comparison was made between the calculation result and the result of laboratory rapid uniaxial compression test to verify the model. Dynamic response and damage features of the tunnel-base structure in three working conditions, i.e., filled layer without seam, filled layer with seam and ground water, and filled layer with seam but without water, were analyzed based on engineering practice. According to the analysis, structural vibration response was intensified in the presence of a seam. With and without groundwater, the vertical dynamic stress attenuation at both sides of the seam was 41.07% and 47.13%, respectively; vertical vibration acceleration was attenuated by 91.17% and 91.73%, respectively; and the acceleration amplitude at the upper structure of the seam increased by 724.67% and 765.02%, respectively. Groundwater in the seam would aggravate damage accumulation. It could be seen from the analysis that the current design parameters satisfied the antifatigue requirements within the design reference period at a train speed of 300 km/h when there was no seam in the tunnel-base concrete structure with IV-class surrounding rocks. When there was a seam in the tunnel-base concrete structure, however, antifatigue life was 56 years in the presence of groundwater and 62 years without groundwater, which suggested that current design parameters failed to satisfy the antifatigue requirements within the design reference period.


---

## Body

## 1. Introduction

Currently, high-speed railways have been constructed at a furious pace in China, which makes trips more convenient and has even fundamentally changed the conventional lifestyle. Yet, some challenging problems, especially those related to safe operation and disaster prevention and mitigation, have also arisen in high-speed railway construction and operation. In most cases, the high-speed railway base is constructed using a ballast-less track slab with a large rigidity. The concrete structure under it bears strong vibration, while the high ride comfort during high-speed train-operation relies on high dynamic stability of the concrete structure.As a multiphase composite, concrete naturally features microscopic damage defects, such as stochastically distributed microvoids and microfractures. In addition, some initial macroscopic defects, such as a seam in the filled layer [1], lining crack, and sedimentation deformation, can be observed at some tunnel-base structures during initial operation due to different construction techniques, such as tunnel-base concrete structure stratification and casting in sections, and various factors, such as surrounding rocks and soil and groundwater environment. With the coupling dynamics of long-term periodic dynamic load of high-speed trains and groundwater, the aforementioned microscopic damage and initial macroscopic defects will inevitably result in a defect-accumulation effect, which leads to local resistance fading in the structure. Macroscopic concrete failure can occur as the defects develop to a certain level over time, thus affecting the train-operation safety.Current studies on the bearing behavior and disaster-causing mechanism of high-speed railway tunnel-base structure mostly focus on the mechanical performance nonlinearity and load response mechanism of concrete (Huang [2]; Ding [3]; Cai [4]; Li [5]; and Wang [6]), in which the constitutive models are targeted at the nonlinearity and strength-reduction characteristics during concrete failure. In terms of microscopic damage defect distribution, investigation suggests that the destruction of base concrete structure shows stochastic statistical characteristics along the direction of the tunnel axis. The reason is that there are stochastically distributed microfissures and microvoids within the concrete materials that conceivably impact the evolution process of concrete damage and destruction. Macroscopically, this is manifested by the significant stochasticity of concrete’s mechanical property. Li et al. [7] also pointed out that nonnegligible stochastic distribution characteristics could be observed in the initial damage distribution and subsequent evolution process of concrete due to the stochastic distribution of concrete material components.In addition, most research on the mechanical behavior of high-speed railway tunnel-base has been focused on the dynamic response of lining structure in an ideal condition (Huang et al. [8]). However, they have failed to probe the damage dynamics and fatigue issues of the concrete structure with initial macroscopic defects in the coupling effects of dynamic load of high-speed trains and groundwater while studying the microscopic damage defects. Considering the characteristics of high-speed railway lining structure, such as high maintenance difficulty, high cost, and irreversibility, it is urgent that relevant studies be carried out.Concrete damage is an irreversible energy-dissipation process, the damage mechanics of which involves energy and irreversible thermodynamics. In the thermodynamic framework, damage affects not only elastic deformation but also plastic deformation. Normally, for the sake of numerical calculation, elastoplastic and damage deformation is decoupled, which means that the plasticity has a direct influence on damage in a current load step, while the influence of damage on plasticity will be reflected in the next load step. However, four parts, including the plastic-hardening parameter, can be derived from the elastoplastic decoupling of the Helmholtz free-energy equation, in which plastic Helmholtz free energy is also associated with the hardening parameter. To fully reflect the thermodynamic compatibility conditions, the influence of the hardening parameter should not be ignored. In this study, the hardening parameter was also taken into consideration during derivation, from which the corrected elastoplastic Helmholtz free energy was derived.In this study, a stochastic dynamic damage constitutive model for high-speed railway tunnel-base concrete was derived from the improved elastoplastic Helmholtz free energy and energy-dissipation principle using the stochastic damage theory according to the energy variable theory and the laws of thermodynamics. The energy equivalent strain concept was introduced. Damage-evolution rules and damage criteria were built in the effective stress space, which were entered into the calculation software by means of secondary programming. The validity of the proposed constitutive model was verified by comparing the calculation result with the laboratory test result based on which the dynamic response and damage characteristics of the defective high-speed tunnel bases in different working conditions under load were analyzed by the way of engineering practice. The fatigue life of tunnel-base structure was predicted based on the calculation result.


## 2. Modified Elastoplastic Helmholtz Free-Energy-Based Stochastic Dynamic Damage Constitutive Model

### 2.1. Stress State of Tunnel-Base Concrete Structure

The concrete-failure pattern under the effect of complex stress can be generally into three types: tensile failure, compression failure, and crushing failure under high hydrostatic pressure. A strain-stress test of the typical tunnel base by Yin et al. [9] suggested that the stress of concrete at the invert of the tunnel bottom was subject to repeated changes of “tension-compression-tension-compression.” Thus, concrete damage and failure primarily originated from the tensile damage and compressive damage without considering the strain hardening caused by high hydrostatic pressure. Therefore, tensile damage and compressive damage should be distinguished in the selection of variables appropriate for tunnel-base concrete structure damage.


### 2.2. Thermodynamic Fundamentals and Basic Equation of Damage Constitutive Model

Considering that concrete microcracks cannot be healed and are followed by inevitable energy dissipation, the constitutive model could be built based on the irreversibility of damage-incurred energy dissipation.According to the first law of thermodynamics, the rate of change of kinetic energy and internal energy within a system with time is equal to the sum of the power of work done by external force and the rate of change of thermal energy transferred to the system from outside; its differential form can be expressed as(1)ρe˙=ρq˙+σ:ε˙−divh,where ρ is object density; e˙ is the differential of internal energy density with respect to time; q˙ is the heat generated by the internal heat source per unit mass in unit time; σ is the effective stress tensor; ε˙ is the differential of the effective strain tensor ε; the operator “:” is the second-order contraction product of the tensor, i.e., the power σ:ε˙=σijε˙ij of work by external force; div() is the divergence operator; and h is the heat outflow across the unit area in unit time.There are single-valued functions of the thermodynamic temperature and entropy in the irreversible thermodynamics system (the second thermodynamics system). The relationship between energy dissipation and temperature gradient is usually reflected by the Clausius–Duhem dissipated-energy inequality:(2)ζ=ρTS˙−e˙+σ:ε˙−1Th∇T≥0,where ζ is the energy dissipation per unit volume; T is the function of empirical temperature, which is constantly positive; S is the entropy; S˙ is the change in entropy; ▽T is the gradient, ▽T=(∂T/∂xi)ei; xi is the space coordinate; and ei is the corresponding unit vector.Entropy changes of the system and environment should be considered when the entropy function serves as the basis of free-energy judgement. However, considering that the thermodynamic reaction usually takes place in the isothermal and isobaric condition or in the isothermal and isovolumetric condition, it would be much more convenient to judge the direction and limit of thermodynamic reaction using the change values of the system state function in this condition. Helmholtz (Von Helmholtz, H.L.P) defined a new state function:(3)ψ≡U−TS=ρe−TS,where ψ is the Helmholtz free energy, and U is the internal energy of system, which, in the physical sense, means the max work that the system does equal to the decrement of the Helmholtz free energy of the system in the isothermal irreversible process.In an isothermal pure mechanical process, temperature functionT is a constant regardless of the influence of the internal heat source q and heat flow h. (1) of the first law of thermodynamics and that, (2), of the second law of thermodynamics were reduced to(4)ρe˙=σ:ε˙,(5)ζ=ρTS˙−e˙+σ:ε˙≥0.Considering that the temperature functionT was a constant, the differential form of (3) is(6)ψ˙=ρe˙−TS˙.Thus, (4) and (5) could be combined:(7)ζ=σ:ε˙−ψ˙≥0.From (7), some of the work by external force was converted into Helmholtz free energy, and some was used for energy dissipation during deformation, which was a nonnegative value.Definitions:D is the damage variable that signifies the degree of the loss of Helmholtz free energy; D+ the degree of the loss of Helmholtz free energy when the concrete was under tension, which is called the tensile damage variable; and D－the degree of the loss of Helmholtz free energy when the concrete was under compression, which is called the compressive damage variable. Since concrete damage is an irreversible thermodynamic process and there was strain with energy dissipation during concrete failure process, the Helmholtz free energy was the state function of effective strain tensor ε and damage variable:(8)ψ=ψε,D.By substituting (8) into (7),(9)ζ=σ−∂ψ∂ε:ε˙−∂ψε,D∂DD˙≥0,could be derived.In view of the stochasticity of the strain rateε˙, the sufficient and necessary condition to satisfy the above inequality is(10)σ−∂ψε,D∂ε=0,ζ=−∂ψε,D∂DD˙≥0.It could be derived from (10) that(11)σ=∂ψε,D∂ε.Equation (11) is the fundamental equation for the damage mechanical constitutive model since it contains the damage variable D and the stress-strain relationship of the effective strain ε. The damage-evolution equation could be inferred from the max-dissipation-energy principle (Li et al. [7]):(12)D˙=λ˙D∂G∂β,where D˙ is the differential of damage variable with respect to time, λ˙D is the damage-evolution factor, and G is the damage trial function or dissipation function that could be obtained by experience.


### 2.3. Corrected Elastoplastic Helmholtz Free Energy

When the effective strain tensorε was decomposed into the elastic strain tensor εe and plastic strain tensor εp, ε  = εe + εp. Supposing that there was no coupling between the elastic Helmholtz free energy ψe and plastic Helmholtz free energy ψp of the material in the isothermal condition, the total elastoplastic Helmholtz free energy ψ is their sum:(13)ψεe,κ,D+,D−=ψeεe,D+,D−+ψpκ,D+,D−,where(14)ψeεe,D+,D−=ψe+εe,D++ψe−εe,D−,(15)ψe+εe,D+=1−D+ψe+ε=121−D+σ+:εe,(16)ψe−εe,D−=1−D−ψe−ε=121−D−σ−:εe,(17)ψpκ,D+,D−=ψp+κ,D++ψp−κ,D−,(18)ψp+κ,D+=1−D+ψ0p+κ=1−D+∫0εpσ+:dεp,(19)ψp−κ,D−=1−D−ψ0p−κ=1−D−∫0εpσ−:dεp,where κ is the hardening parameter, “+” means tension, “－” means compression, ψ0p+ is the initial tensile plastic free energy, and ψ0p− is the initial compressive plastic free energy.The variablevk was defined to characterize the change of the internal state of thermodynamic system and fk the corresponding generalized force variable (adjoint variable). If ε and T signify the state variables, the corresponding adjoint variables are force σ and entropy S. If elastic strain εe, plastic strain εp, accumulated plastic strain q, and damage variable D are taken as the internal state variables, the corresponding conjugate variable, radius of the yield surface, and damage energy release rate are σ, R, and β, respectively. According to Li et al. [7], when ψ is the function of ε, D, and vk, the first law of thermodynamics and the second law of thermodynamics could be expressed as(20)σ−ρ∂ψ∂ε:ε˙−S+∂ψ∂TT˙−ρTS˙−ρ∂ψ∂vk:vk−ρ∂ψ∂DD˙+ρr−divq=0,(21)σ−ρ∂ψ∂ε:ε˙−S+∂ψ∂TT˙−ρTS˙−ρ∂ψ∂vk:vk−ρ∂ψ∂DD˙−qTgradT≥0.Since the laws of thermodynamics hold with any strain rate and the rate of temperature change, the necessary conditions for (20) and (21) are(22)σ=ρ∂ψ∂ε,s=−∂ψ∂T,fk=ρ∂ψ∂vk,β=−ρ∂ψ∂D.Equation (22) could be equivalently transformed into(23)σ±=∂ψe±∂εe,β±=−∂ψ±∂D±,lq=−∂ψp∂q,where(24)lq=fc1+α1+2αqq0−αq2q02,q≤q0,fc,q>q0.q0T=q10+q10−q20sin232T,α=fcfcs,q10=q0π3,q20=q00,where lq is the hardening function, fcs the initial yield stress of uniaxial compression, fc the uniaxial compression strength, q0 the initial plastic accumulated strain, and q10 and q20 the equivalent plastic strain values corresponding to the peak stress under uniaxial compression and biaxial equal-strength compression, respectively.By derivation of (13) and substituting the three equalities of (23) into (13),(25)σ+−ρ∂ψe+∂εe:εe+σ−−∂ψe−∂εe:εe+σ:ε˙p−∂ψ+∂D+D˙+−∂ψ−∂D−D˙−−∂ψp∂qq˙≥0.From the elastoplastic decoupling of the energy release rate,(26)β±=βe±+βp±=ψe0±+ψp0±,where ψe0± and ψp0± are the initial elastic Helmholtz free energy and initial plastic Helmholtz free energy that are numerically equal to the elastic part and plastic part of the energy release rate, respectively.Faria et al. [10] proposed calculation of ψe0± using the following equation:(27)ψe0±=12σ±:εe=12σ±:C0:σ,where C0 is the initial flexibility tensor.Wu et al. [11] argued that ψ0p+ = 0 in (15) and (16) reduces the computational complexity since the evident quasibrittleness of concrete could be observed under tension, and plastic deformation was minimal after unloading. However, the third term of (23) suggests that plastic Helmholtz free energy was associated with the hardening parameter. To fully reflect the thermodynamic compatibility conditions, the influence of hardening parameter should not be ignored. The term −∫ldq was added when determining the plastic Helmholtz free energy, which means that ψ0p is defined as(28)ψ0pκ=ψ0p+κ+ψ0p−κ=−∫ldq+ψ0p−κ.Thus, (18) is converted into(29)ψp+κ,D+=1−D+ψ0p+κ−∫ldq=−∫ldq.The corrected plastic Helmholtz free energy could be obtained by combining equations (14)–(18), (19), and (29):(30)ψεe,κ,D+,D−=121−D+σ+:εe+121−D−σ−:εe+1−D−∫0εpσ−:dεp−∫ldq.


### 2.4. Stochastic Dynamic Damage Constitutive Model

The material damage constitutive relationship could be derived from the uniform (7) of thermodynamics laws. The differential of (30) is taken with respect to time and substituted into (7), thus obtaining(31)σ−∂ψe∂εe:ε˙e+−∂ψ∂D+D˙++−∂ψ∂D−D˙−+σ:ε˙p−−∂ψp∂κ¯κ¯˙≥0.Owing to the stochasticity ofε˙,(32)σ=∂ψeεe,D+,D−∂εe.If equation (31) holds, the expression(33)σ=1−D+∂ψ0e+∂εe+1−D−∂ψ0e−∂εe,could be obtained by substituting (14) into (32).The elastoplastic stochastic damage constitutive relation could be derived by substituting (15) into (33):(34)σ=1−D+σ++1−D−σ−=1−D:σ=1−DE0:ε−εP,where(35)D=D+P++D−P−,(36)P+=∑iHσ^iPi,(37)P−=I,P+,where D is the fourth-order tensor, and P+ and P− are the positive and negative fourth-order projection tensors, respectively, of stress tensor σ.Damage evolution results in the anisotropy of concrete under two-dimensional plane stress. Therefore, it is better to symmetrize the corresponding secant stiffness matrix. According to (34), the stochastic damage stress-strain relation in the principal stress direction under biaxial tension and compression is(38)σ1=E01−v021−D+ε1e+1−D+1−D−v0ε3e,(39)σ3=E01−v021−D+1−D−v0ε31e+1−D−ε3e,where σ1 is the max principal stress, σ3 is the min principal stress, E0 is the initial elastic modulus of materials, v0 is the initial Poisson’s ratio of materials, ε1e is the elastic strain in the direction of max principal stress, and ε3e is the elastic strain in the direction of min principal stress.The microscopic failure of tensile damage and compressive damage of concrete at the tunnel bottom was stochastic and uncontrolled. Mathematically, the damage variables were stochastic variables with probability description, such as the probability density distribution and homogeneity and the standard deviation. Since the damage evolution was complex, the relationship between stochastic damage and deterministic damage is determined by taking the mean value of the damage. Li et al. [7] developed the uniform expression of the homogeneity and standard deviation of concrete tensile and compressive damage:(40)μDεej=FΔθiεej,(41)VD2εej=2∫011−ηFΔθiεej,θiεej;ηdη−FΔ2θiεej,where μD is the mean function of damage D; subscript j (j = 1, 2) of εe is the elastic strain in the direction of max or min principal stress; F is the one-dimensional function of stochastic field (x); θ is the equivalent coefficient of energy transformation, and the subscript i (i = 1, 2) of θ is the equivalent transformation coefficient of tensile and compressive energy, respectively; VD2 is the variance of damage D; and η is the absolute value of the relative distance between independent variables of the probability density function in the one-dimensional stochastic field.The damage criterion is the basis for determining whether there is damage, and the process of damage failure is a process in which balance and stability are broken and reestablished at any time. Therefore, the principle of minimum energy dissipation also applies to the establishment of the damage criterion. When the energy dissipation rate of the energy dissipation system and the required constraint conditions for energy dissipation are known, the conditions for the stationary value of the system energy dissipation can be identified, thus establishing the uniform damage criterion. The flow rule of plastic theory can be expressed by(42)ε˙p=λ˙p∂f∂σ,where f is the yield function of the effective stress space and λ˙p is the differential of the plastic multiplier λp.Considering that the energy dissipation of the damage system is generated by two mechanisms, the damage energy dissipationΘD and plastic dissipation Θp, the total energy dissipation Θ is(43)Θ=ΘD+Θp=σ:ε˙D+σ:ε˙p,where ε˙D is the damage strain rate and ε˙p the plastic strain rate.The energy dissipation in the damage process is subject to the restriction of the damage criterionG (σ) = 0. According to the principle of the maximum energy-dissipation extreme, the conditions for taking the extremum should satisfy(44)∂Θ+ϑG∂σ=0,where ϑ is the Lagrange multiplier and G is the damage criterion.Substituting (43) into (44),(45)∂G∂σ=−2D˙ϑ1−D2σ:C0−λ˙pϑ∂f∂σ−λ˙pϑσ:∂2f∂σ⊗∂σ,where C0 is mathematically an integration constant when compared with variable σ and its physical meaning is the influence of damage D on G. It is a parameter related to the damage threshold.When taking the differential of (45) with respect to σ,(46)Gσ=β−β0,could be derived if D was fixed, i.e., D was independent of σ, where(47)β=σ:C0:σ2+f+σ:∂2f∂σ∂σ:σλ˙p1−D22D˙,(48)β0=ϑ1−D2C02D˙,where ß0 is the threshold of the damage factor. Physically, this means controlling the evolution of the damage surface.The key to damage theory is to determine the damage-evolution rule, i.e., the rule of damage variableD varying with the other internal variables. The evolution rule of the damage variable is(49)D˙±=λ˙pD±∂G±∂β±,(50)β˙0±=λ˙pD±,according to the rule of orthogonal flow based on the damage criterion shown in (42), where the superscripts “+” and “–” are the states of tension and compression, respectively.(51)β˙n±=μ±β0±βv±βn±−1n±.The expression was obtained by dynamic expansion of the damage energy release rateβ (see Li et al. [7]), where βv± is the dynamic damage energy release rate, μ± and n± are dynamic effect empirical parameters, and μ+  = 2.1 × 109 N/s, μ−  = 6 × 1010 N/s, n+  = 5.5, and n−  = 4.5 according to Liu [12].Equation (51) is a first-order equation regarding βn±. The dynamic damage-evolution process under tension and compression could be solved using (52) and (53), respectively:(52)β˙n+=μ+β0+βv+βn+−1n+,D+=1−β0+βn+A+1−βn+β0+,(53)β˙n−=μ−β0−βv−βn−−1n−,D−=1−β0+βn+1−A−−A−B+1−βn−β0−,where A+, A−, B+, and B− are model parameters. Their physical meanings are the softening effects of concrete before and after uniaxial peak stress under tension and compression. They could be calibrated by the stress-strain curves of uniaxial tension and compression obtained in the test [12].The damage loading and unloading condition, i.e., the Kuhn–Tucker relationship, is(54)λ˙pD±≥0,G±≤0,λ˙pD±G±=0.WhenG±<0, λ˙pD±=0 and D± = 0; the concrete was in the damage unloading phase or neutral loading phase. When λ˙pD±>0, G± = 0. λ˙pD± was obtained based on the damage consistency conditions:(55)G±β±,β0±=0⟹β0±=β±,λ˙D±=β˙0±=β±≥0,where G˙ is the differential form of the damage criterion.Considering the initial conditionD0± = 0, it was obtained by substituting (55) into (49).(56)D˙±=β∂G±β±∂β±=G±β±=G±β0±⟺D±=G±β0±.Therefore, the damage variable could be obtained by (56) once the expression of function G± is given:(57)D+=G+β0±=1−β0initial+β0±1−A+expB+1−β0±β0initial++A+,β0±≥β0initial+,(58)D−=G−β0−=1−β0initial−β0−1−A−−A−expB+1−β0+β0initial+,β0−≥β0initial−,where ß0initial± is the threshold of the initial energy damage rate that physically signifies that there will be damage when the damage energy release rate is higher than this threshold. ß0initial+ is the linear limit strength of concrete under tension, which is the uniaxial tensile strength; ß0initial-is the linear limit strength of concrete under compression, ß0initial-= (1 − α) f0, f0 = (0.3 – 0.5) fc, where fc is the uniaxial compressive strength; and α is the function of the ratio Δ between biaxial compressive yield strength fby- and uniaxial compressive yield strength fy-. The calculation formula for α is(59)α=Δ−12Δ−1,where the range of Δ is 1.1–1.2.


## 2.1. Stress State of Tunnel-Base Concrete Structure

The concrete-failure pattern under the effect of complex stress can be generally into three types: tensile failure, compression failure, and crushing failure under high hydrostatic pressure. A strain-stress test of the typical tunnel base by Yin et al. [9] suggested that the stress of concrete at the invert of the tunnel bottom was subject to repeated changes of “tension-compression-tension-compression.” Thus, concrete damage and failure primarily originated from the tensile damage and compressive damage without considering the strain hardening caused by high hydrostatic pressure. Therefore, tensile damage and compressive damage should be distinguished in the selection of variables appropriate for tunnel-base concrete structure damage.


## 2.2. Thermodynamic Fundamentals and Basic Equation of Damage Constitutive Model

Considering that concrete microcracks cannot be healed and are followed by inevitable energy dissipation, the constitutive model could be built based on the irreversibility of damage-incurred energy dissipation.According to the first law of thermodynamics, the rate of change of kinetic energy and internal energy within a system with time is equal to the sum of the power of work done by external force and the rate of change of thermal energy transferred to the system from outside; its differential form can be expressed as(1)ρe˙=ρq˙+σ:ε˙−divh,where ρ is object density; e˙ is the differential of internal energy density with respect to time; q˙ is the heat generated by the internal heat source per unit mass in unit time; σ is the effective stress tensor; ε˙ is the differential of the effective strain tensor ε; the operator “:” is the second-order contraction product of the tensor, i.e., the power σ:ε˙=σijε˙ij of work by external force; div() is the divergence operator; and h is the heat outflow across the unit area in unit time.There are single-valued functions of the thermodynamic temperature and entropy in the irreversible thermodynamics system (the second thermodynamics system). The relationship between energy dissipation and temperature gradient is usually reflected by the Clausius–Duhem dissipated-energy inequality:(2)ζ=ρTS˙−e˙+σ:ε˙−1Th∇T≥0,where ζ is the energy dissipation per unit volume; T is the function of empirical temperature, which is constantly positive; S is the entropy; S˙ is the change in entropy; ▽T is the gradient, ▽T=(∂T/∂xi)ei; xi is the space coordinate; and ei is the corresponding unit vector.Entropy changes of the system and environment should be considered when the entropy function serves as the basis of free-energy judgement. However, considering that the thermodynamic reaction usually takes place in the isothermal and isobaric condition or in the isothermal and isovolumetric condition, it would be much more convenient to judge the direction and limit of thermodynamic reaction using the change values of the system state function in this condition. Helmholtz (Von Helmholtz, H.L.P) defined a new state function:(3)ψ≡U−TS=ρe−TS,where ψ is the Helmholtz free energy, and U is the internal energy of system, which, in the physical sense, means the max work that the system does equal to the decrement of the Helmholtz free energy of the system in the isothermal irreversible process.In an isothermal pure mechanical process, temperature functionT is a constant regardless of the influence of the internal heat source q and heat flow h. (1) of the first law of thermodynamics and that, (2), of the second law of thermodynamics were reduced to(4)ρe˙=σ:ε˙,(5)ζ=ρTS˙−e˙+σ:ε˙≥0.Considering that the temperature functionT was a constant, the differential form of (3) is(6)ψ˙=ρe˙−TS˙.Thus, (4) and (5) could be combined:(7)ζ=σ:ε˙−ψ˙≥0.From (7), some of the work by external force was converted into Helmholtz free energy, and some was used for energy dissipation during deformation, which was a nonnegative value.Definitions:D is the damage variable that signifies the degree of the loss of Helmholtz free energy; D+ the degree of the loss of Helmholtz free energy when the concrete was under tension, which is called the tensile damage variable; and D－the degree of the loss of Helmholtz free energy when the concrete was under compression, which is called the compressive damage variable. Since concrete damage is an irreversible thermodynamic process and there was strain with energy dissipation during concrete failure process, the Helmholtz free energy was the state function of effective strain tensor ε and damage variable:(8)ψ=ψε,D.By substituting (8) into (7),(9)ζ=σ−∂ψ∂ε:ε˙−∂ψε,D∂DD˙≥0,could be derived.In view of the stochasticity of the strain rateε˙, the sufficient and necessary condition to satisfy the above inequality is(10)σ−∂ψε,D∂ε=0,ζ=−∂ψε,D∂DD˙≥0.It could be derived from (10) that(11)σ=∂ψε,D∂ε.Equation (11) is the fundamental equation for the damage mechanical constitutive model since it contains the damage variable D and the stress-strain relationship of the effective strain ε. The damage-evolution equation could be inferred from the max-dissipation-energy principle (Li et al. [7]):(12)D˙=λ˙D∂G∂β,where D˙ is the differential of damage variable with respect to time, λ˙D is the damage-evolution factor, and G is the damage trial function or dissipation function that could be obtained by experience.


## 2.3. Corrected Elastoplastic Helmholtz Free Energy

When the effective strain tensorε was decomposed into the elastic strain tensor εe and plastic strain tensor εp, ε  = εe + εp. Supposing that there was no coupling between the elastic Helmholtz free energy ψe and plastic Helmholtz free energy ψp of the material in the isothermal condition, the total elastoplastic Helmholtz free energy ψ is their sum:(13)ψεe,κ,D+,D−=ψeεe,D+,D−+ψpκ,D+,D−,where(14)ψeεe,D+,D−=ψe+εe,D++ψe−εe,D−,(15)ψe+εe,D+=1−D+ψe+ε=121−D+σ+:εe,(16)ψe−εe,D−=1−D−ψe−ε=121−D−σ−:εe,(17)ψpκ,D+,D−=ψp+κ,D++ψp−κ,D−,(18)ψp+κ,D+=1−D+ψ0p+κ=1−D+∫0εpσ+:dεp,(19)ψp−κ,D−=1−D−ψ0p−κ=1−D−∫0εpσ−:dεp,where κ is the hardening parameter, “+” means tension, “－” means compression, ψ0p+ is the initial tensile plastic free energy, and ψ0p− is the initial compressive plastic free energy.The variablevk was defined to characterize the change of the internal state of thermodynamic system and fk the corresponding generalized force variable (adjoint variable). If ε and T signify the state variables, the corresponding adjoint variables are force σ and entropy S. If elastic strain εe, plastic strain εp, accumulated plastic strain q, and damage variable D are taken as the internal state variables, the corresponding conjugate variable, radius of the yield surface, and damage energy release rate are σ, R, and β, respectively. According to Li et al. [7], when ψ is the function of ε, D, and vk, the first law of thermodynamics and the second law of thermodynamics could be expressed as(20)σ−ρ∂ψ∂ε:ε˙−S+∂ψ∂TT˙−ρTS˙−ρ∂ψ∂vk:vk−ρ∂ψ∂DD˙+ρr−divq=0,(21)σ−ρ∂ψ∂ε:ε˙−S+∂ψ∂TT˙−ρTS˙−ρ∂ψ∂vk:vk−ρ∂ψ∂DD˙−qTgradT≥0.Since the laws of thermodynamics hold with any strain rate and the rate of temperature change, the necessary conditions for (20) and (21) are(22)σ=ρ∂ψ∂ε,s=−∂ψ∂T,fk=ρ∂ψ∂vk,β=−ρ∂ψ∂D.Equation (22) could be equivalently transformed into(23)σ±=∂ψe±∂εe,β±=−∂ψ±∂D±,lq=−∂ψp∂q,where(24)lq=fc1+α1+2αqq0−αq2q02,q≤q0,fc,q>q0.q0T=q10+q10−q20sin232T,α=fcfcs,q10=q0π3,q20=q00,where lq is the hardening function, fcs the initial yield stress of uniaxial compression, fc the uniaxial compression strength, q0 the initial plastic accumulated strain, and q10 and q20 the equivalent plastic strain values corresponding to the peak stress under uniaxial compression and biaxial equal-strength compression, respectively.By derivation of (13) and substituting the three equalities of (23) into (13),(25)σ+−ρ∂ψe+∂εe:εe+σ−−∂ψe−∂εe:εe+σ:ε˙p−∂ψ+∂D+D˙+−∂ψ−∂D−D˙−−∂ψp∂qq˙≥0.From the elastoplastic decoupling of the energy release rate,(26)β±=βe±+βp±=ψe0±+ψp0±,where ψe0± and ψp0± are the initial elastic Helmholtz free energy and initial plastic Helmholtz free energy that are numerically equal to the elastic part and plastic part of the energy release rate, respectively.Faria et al. [10] proposed calculation of ψe0± using the following equation:(27)ψe0±=12σ±:εe=12σ±:C0:σ,where C0 is the initial flexibility tensor.Wu et al. [11] argued that ψ0p+ = 0 in (15) and (16) reduces the computational complexity since the evident quasibrittleness of concrete could be observed under tension, and plastic deformation was minimal after unloading. However, the third term of (23) suggests that plastic Helmholtz free energy was associated with the hardening parameter. To fully reflect the thermodynamic compatibility conditions, the influence of hardening parameter should not be ignored. The term −∫ldq was added when determining the plastic Helmholtz free energy, which means that ψ0p is defined as(28)ψ0pκ=ψ0p+κ+ψ0p−κ=−∫ldq+ψ0p−κ.Thus, (18) is converted into(29)ψp+κ,D+=1−D+ψ0p+κ−∫ldq=−∫ldq.The corrected plastic Helmholtz free energy could be obtained by combining equations (14)–(18), (19), and (29):(30)ψεe,κ,D+,D−=121−D+σ+:εe+121−D−σ−:εe+1−D−∫0εpσ−:dεp−∫ldq.


## 2.4. Stochastic Dynamic Damage Constitutive Model

The material damage constitutive relationship could be derived from the uniform (7) of thermodynamics laws. The differential of (30) is taken with respect to time and substituted into (7), thus obtaining(31)σ−∂ψe∂εe:ε˙e+−∂ψ∂D+D˙++−∂ψ∂D−D˙−+σ:ε˙p−−∂ψp∂κ¯κ¯˙≥0.Owing to the stochasticity ofε˙,(32)σ=∂ψeεe,D+,D−∂εe.If equation (31) holds, the expression(33)σ=1−D+∂ψ0e+∂εe+1−D−∂ψ0e−∂εe,could be obtained by substituting (14) into (32).The elastoplastic stochastic damage constitutive relation could be derived by substituting (15) into (33):(34)σ=1−D+σ++1−D−σ−=1−D:σ=1−DE0:ε−εP,where(35)D=D+P++D−P−,(36)P+=∑iHσ^iPi,(37)P−=I,P+,where D is the fourth-order tensor, and P+ and P− are the positive and negative fourth-order projection tensors, respectively, of stress tensor σ.Damage evolution results in the anisotropy of concrete under two-dimensional plane stress. Therefore, it is better to symmetrize the corresponding secant stiffness matrix. According to (34), the stochastic damage stress-strain relation in the principal stress direction under biaxial tension and compression is(38)σ1=E01−v021−D+ε1e+1−D+1−D−v0ε3e,(39)σ3=E01−v021−D+1−D−v0ε31e+1−D−ε3e,where σ1 is the max principal stress, σ3 is the min principal stress, E0 is the initial elastic modulus of materials, v0 is the initial Poisson’s ratio of materials, ε1e is the elastic strain in the direction of max principal stress, and ε3e is the elastic strain in the direction of min principal stress.The microscopic failure of tensile damage and compressive damage of concrete at the tunnel bottom was stochastic and uncontrolled. Mathematically, the damage variables were stochastic variables with probability description, such as the probability density distribution and homogeneity and the standard deviation. Since the damage evolution was complex, the relationship between stochastic damage and deterministic damage is determined by taking the mean value of the damage. Li et al. [7] developed the uniform expression of the homogeneity and standard deviation of concrete tensile and compressive damage:(40)μDεej=FΔθiεej,(41)VD2εej=2∫011−ηFΔθiεej,θiεej;ηdη−FΔ2θiεej,where μD is the mean function of damage D; subscript j (j = 1, 2) of εe is the elastic strain in the direction of max or min principal stress; F is the one-dimensional function of stochastic field (x); θ is the equivalent coefficient of energy transformation, and the subscript i (i = 1, 2) of θ is the equivalent transformation coefficient of tensile and compressive energy, respectively; VD2 is the variance of damage D; and η is the absolute value of the relative distance between independent variables of the probability density function in the one-dimensional stochastic field.The damage criterion is the basis for determining whether there is damage, and the process of damage failure is a process in which balance and stability are broken and reestablished at any time. Therefore, the principle of minimum energy dissipation also applies to the establishment of the damage criterion. When the energy dissipation rate of the energy dissipation system and the required constraint conditions for energy dissipation are known, the conditions for the stationary value of the system energy dissipation can be identified, thus establishing the uniform damage criterion. The flow rule of plastic theory can be expressed by(42)ε˙p=λ˙p∂f∂σ,where f is the yield function of the effective stress space and λ˙p is the differential of the plastic multiplier λp.Considering that the energy dissipation of the damage system is generated by two mechanisms, the damage energy dissipationΘD and plastic dissipation Θp, the total energy dissipation Θ is(43)Θ=ΘD+Θp=σ:ε˙D+σ:ε˙p,where ε˙D is the damage strain rate and ε˙p the plastic strain rate.The energy dissipation in the damage process is subject to the restriction of the damage criterionG (σ) = 0. According to the principle of the maximum energy-dissipation extreme, the conditions for taking the extremum should satisfy(44)∂Θ+ϑG∂σ=0,where ϑ is the Lagrange multiplier and G is the damage criterion.Substituting (43) into (44),(45)∂G∂σ=−2D˙ϑ1−D2σ:C0−λ˙pϑ∂f∂σ−λ˙pϑσ:∂2f∂σ⊗∂σ,where C0 is mathematically an integration constant when compared with variable σ and its physical meaning is the influence of damage D on G. It is a parameter related to the damage threshold.When taking the differential of (45) with respect to σ,(46)Gσ=β−β0,could be derived if D was fixed, i.e., D was independent of σ, where(47)β=σ:C0:σ2+f+σ:∂2f∂σ∂σ:σλ˙p1−D22D˙,(48)β0=ϑ1−D2C02D˙,where ß0 is the threshold of the damage factor. Physically, this means controlling the evolution of the damage surface.The key to damage theory is to determine the damage-evolution rule, i.e., the rule of damage variableD varying with the other internal variables. The evolution rule of the damage variable is(49)D˙±=λ˙pD±∂G±∂β±,(50)β˙0±=λ˙pD±,according to the rule of orthogonal flow based on the damage criterion shown in (42), where the superscripts “+” and “–” are the states of tension and compression, respectively.(51)β˙n±=μ±β0±βv±βn±−1n±.The expression was obtained by dynamic expansion of the damage energy release rateβ (see Li et al. [7]), where βv± is the dynamic damage energy release rate, μ± and n± are dynamic effect empirical parameters, and μ+  = 2.1 × 109 N/s, μ−  = 6 × 1010 N/s, n+  = 5.5, and n−  = 4.5 according to Liu [12].Equation (51) is a first-order equation regarding βn±. The dynamic damage-evolution process under tension and compression could be solved using (52) and (53), respectively:(52)β˙n+=μ+β0+βv+βn+−1n+,D+=1−β0+βn+A+1−βn+β0+,(53)β˙n−=μ−β0−βv−βn−−1n−,D−=1−β0+βn+1−A−−A−B+1−βn−β0−,where A+, A−, B+, and B− are model parameters. Their physical meanings are the softening effects of concrete before and after uniaxial peak stress under tension and compression. They could be calibrated by the stress-strain curves of uniaxial tension and compression obtained in the test [12].The damage loading and unloading condition, i.e., the Kuhn–Tucker relationship, is(54)λ˙pD±≥0,G±≤0,λ˙pD±G±=0.WhenG±<0, λ˙pD±=0 and D± = 0; the concrete was in the damage unloading phase or neutral loading phase. When λ˙pD±>0, G± = 0. λ˙pD± was obtained based on the damage consistency conditions:(55)G±β±,β0±=0⟹β0±=β±,λ˙D±=β˙0±=β±≥0,where G˙ is the differential form of the damage criterion.Considering the initial conditionD0± = 0, it was obtained by substituting (55) into (49).(56)D˙±=β∂G±β±∂β±=G±β±=G±β0±⟺D±=G±β0±.Therefore, the damage variable could be obtained by (56) once the expression of function G± is given:(57)D+=G+β0±=1−β0initial+β0±1−A+expB+1−β0±β0initial++A+,β0±≥β0initial+,(58)D−=G−β0−=1−β0initial−β0−1−A−−A−expB+1−β0+β0initial+,β0−≥β0initial−,where ß0initial± is the threshold of the initial energy damage rate that physically signifies that there will be damage when the damage energy release rate is higher than this threshold. ß0initial+ is the linear limit strength of concrete under tension, which is the uniaxial tensile strength; ß0initial-is the linear limit strength of concrete under compression, ß0initial-= (1 − α) f0, f0 = (0.3 – 0.5) fc, where fc is the uniaxial compressive strength; and α is the function of the ratio Δ between biaxial compressive yield strength fby- and uniaxial compressive yield strength fy-. The calculation formula for α is(59)α=Δ−12Δ−1,where the range of Δ is 1.1–1.2.


## 3. Numerical Implementation of Damage Constitutive Model

As a two-dimensional universal distinct element code, UDEC is a tool providing geotechnical engineers with accurate and effective analysis using explicit problem-solving ability. The explicit problem-solving scheme has been widely used since it provides unstable physical processes with stable solutions and simulates the failure process of the object [13].UDEC could be deemed as a rigid body that satisfies the Mohr–Coulomb criterion in modeling. Model deformation or failure is primarily controlled by the weak planes between blocks. Since the rigidity and strength of the weak plane were lower than those of the block, stress concentration could occur near the weak plane under disturbance, and the model could macroscopically collapse, slide, or crack due to expansion, slippage, fracture, and so on. The weak plane satisfies a certain constitutive relationship and yield function. Concrete cracks and erosive wear initially occur at the junctions between aggregates [7], which means that concrete failure or damage is associated only with the junctions. Considering that the junction strength and rigidity are lower than those of the aggregate, the aggregates in the concrete could be regarded as the blocks in UDEC and the junctions as the weak planes. Thus, the aggregate satisfies the Mohr–Coulomb plastic model and the junction satisfies the stochastic dynamic damage model developed in this study. The two key functions, Initialize() and Run(), could be reloaded using C++ based on the “user-defined joint constitutive models” interface of UDEC. Since yield criterion and tensile-compressive damage combination coefficient are expressed as the principal stress or principal stress invariant, the element principal stress at each circulation step should be obtained, which is performed by the two important functions of UDEC—Resoltopris() and Resoltoglob(). The stochastic dynamic damage constitutive model is embedded into UDEC by means of elastic prediction, plastic correction, and damage correction, thereby establishing the joint constitutive model code “User stochastic,” as shown in Figure 1. The calculation of “elastic prediction-plastic correction-damage correction” is shown in Figure 2.Figure 1
Flow chart of damage constitutive program.Figure 2
Numerical calculation flow.The principle of this process is shown below. Letting the strain beεn, stress σn, plastic strain εnp, and plastic hardening parameter hn at previous increment step tn, the predictive variable of stress σn+1predicted, that of stress-strain εn+1p,predicted, that of plastic hardening parameter hn+1predicted, and that of damage Dn+1±,predicted at the next increment step tn+1 should be(60)σn+1predicted=E0εn+1−εn+1p,predictedεn+1p,predicted=εnphn+1predicted=hnDn+1±,predicted=Dn±.Equation (60) is substituted into yield function f. If the yield function satisfies (61) at the time tn+1,(61)fn+1predicted=fσn+1predicted,hn+1predicted≤0,there would be plastic flow within increment step [tn,tn+1], which suggests that the yield criterion is not satisfied and that plastic correction is required. Since the plastic state of the effective stress space is unrelated to damage, it could be implemented using the return mapping algorithm [13]. Plastic modification is conducted using(62)εn+1p=εn+Δλ∂fσ¯n+1,hn+1∂σ,hn+1=hn+ΔλEhσn+1,hn+1,where Δλ is the increment of plastic flow factor; Eh is the hardening-parameter modulus; and σn+1, Dn+1±, εn+1p, and hn+1 are the effective stress, damage, effective strain, and state variable of the hardening function after plastic correction, respectively.The effective stress tensorσn+1 was divided into tensile and compressive parts σn+1+,σn+1− after plastic correction. The new damage state(63)σn+1=1−Dn+1+σn+1++1−Dn+1−σn+1−,could be obtained by substituting them into (51), the tensile-compressive damage energy release formula, where Dn+1+ and Dn+1− are the tensile and compressive damage variables after plastic correction, respectively.The influence of the stochasticity of concrete initial damage on the nonlinearity of its failure involves the propagation of stochasticity in the physical system. Li et al. [14] investigated the principle of probability conservation during a stochastic event state transition. They built a generalized probability density evolution equation for general stochastic systems and a general theory of stochastic nonlinear analysis, which is also followed in this study. Let Z, G, ϒ, Γ, and x characterize the probability system, probability function, stochastic damage variable reflecting structure material property, structure loading parameter, and spatial location variable determining the location of structure response, respectively. In general loading conditions, the nonlinear structure response could be described by [7].(64)Z=Gϒ,Γ,x.There is a stochastic damage constitutive relation for (64), as shown in (65), which could be used to calculate the probability density evolution process of stress with respect to strain:(65)∂px,ϒ,Γ∂Γ+∑x=1nz˙x∂px,ϒ,Γ∂x=0,where px,ϒ,Γ is the joint probability density function of Z and z˙x is the partial differential of Z with respect to x.When the generalized probability density evolution equation is an explicit expression, an analytical solution exists in the given initial condition. However, it is difficult to set the given explicit expression for the elastoplastic damage constitutive relationship. Therefore, the numerical solution is used in normal cases. A numerical solution could be found using the finite-difference and finite-volume methods. The finite-difference method enjoys distinct advantages since finding the solution of the differential equation could be converted into a recursive algebraic calculation at the increment node, thus significantly reducing the problem complexity. The solving process includes the probability subdivision of the stochastic parameter vector, structure response obtained using the nonlinear analysis, and numerical computation based on the differential. Li and Chen [15] introduced the progress of research in the generalized probability density evolution equation. Gopalaratna and Shah [16] proposed the steps to solve for the nonlinear stochastic response.


## 4. Numerical Verification of Stochastic Dynamic Damage Constitutive Model

Dube et al. [17] modified the classical Lemaitre model using the overstress theory in viscoplasticity mechanics and developed a power-function expression that took the influence of loading rate into consideration. In this section, the loading test proposed by previous researchers is simulated in the same conditions using the proposed stochastic dynamic damage constitutive model. The yield function f¯ in the effective stress space involving the numerical simulation is taken as the linear limit strengthening of materials; f¯+  = 2.04 MPa and f¯−  = 27.5 MPa. With calibration, the other model parameters are A+  = 0.0, A−  = 1.0, B+  = 0.20, and B−  = 0.213. Because the damage variable is a random variable, its mean value and variance should be given to describe the probability. The experimental results of the uniaxial compression fast loading, Dube model results, and numerical analysis results are shown in Figures 3–5.Figure 3
Comparison of calculated results with strain rate of 2 × 10−4/s.Figure 4
Comparison of calculated results with strain rate of 2 × 10−4/s.Figure 5
Comparison of calculated results with strain rate of 12 × 10−2/s.As shown in Figures3–5, the calculation result of the proposed damage constitutive model is consistent with the test result; the concrete strain-softening process under compression could be favorably simulated by the model. However, the test point falls within the range of mean stress plus and minus 1 or 2 times the mean-square error. Thus, the proposed constitutive model not only reflected the stress-strain relationship of concrete under compression but it was also able to simulate its discrete range and fully reflect the nonlinearity, stochasticity, and rate-dependent effect of the mechanical behavior of concrete.


## 5. Dynamic Response and Damage Characteristics of Defective High-Speed Tunnel-Base Structure

### 5.1. Engineering Overview

A double-track high-speed railway tunnel located in Western China was taken as the research object. The formation in the tunnel-site area was composed of IV-class surrounding rocks with joint fissures. Average joint spacing and max depth of the tunnel were 0.5 and 270 m, respectively. The tunnel section is shown in Figure6.Figure 6
Cross section of tunnel lining (unit: mm).According to the periodic reinspection by the railway engineering department using a track-inspection car during tunnel operation, a 1.268-km-long track bed slab bulge was observed along k1679 + 139–+562 of the upline and k1680 + 044–+889 of the downline. Compared with the design values, the max difference of elevation reached 12 cm. An on-site drill-hole survey suggested seam defects (see Figure7) between the base filled layer and the leveling layer, with a max seam gap of 1.2–2.6 cm. There was also high-pressure groundwater, for which groundwater upwelled quickly to the track bed after drilling.Figure 7
Crack in the fill layer.


### 5.2. UDEC Numerical Calculation Model and Material Parameters

To analyze the law of dynamic response and damage characteristics of a base filled layer with seam defects, numerical models of the filled layer with and without a seam were built based on UDEC. The seam was 2 cm wide and 20 cm from the top surface of the filled layer. The space at the left, right, top, and bottom of the model was 4 times the tunnel span. In this paper, the viscous boundary is adopted in the calculation model. The deficiency of model depth was compensated by load. The local enlarged views of the models are shown in Figures8 and 9. Train load depends on axle weight, suspended mass, operating speed, track structure, track regularity, and so on. In this study, the loading method developed by Liang et al. [18] was used. Train load was applied along the left line during simulation. As shown in Figure 10, five monitoring points, A, B, C, D, and E, were arranged at depths of 20, 22, 52, 112.2, and 173.3 cm, respectively, below the top surface of the filled layer along the center line of the left rail line, points A and B of which were located on the upper interface of the seam; D was on the interface between the invert and filled layer; and E was on the interface between the primary support of the invert and surrounding rock.Figure 8
Local mesh division of the numerical model for the filled layer with no crack.Figure 9
Local mesh division of the numerical model for the filled layer with a crack.Figure 10
Schematic diagram of monitoring points (Unit: mm).The yield of blocks in the formation satisfied the Mohr–Coulomb criterion, while the joint yield conformed to the Coulomb sliding criterion of surface contact (see Wang [19]). According to the data specified in the Handbook for Railway Engineering Geology and Code for Design of Railway Tunnel, the physical parameters of formation blocks are shown in Table 1. Supposing that the concrete structure was comprised of coarse aggregates (mean particle size of 2 cm) that combined each other by “corner-corner” contact, “corner-side” contact, and “side-side” contact, the yield of coarse aggregates satisfied the Drucker–Prager criterion. The coarse aggregate interface was simulated using the proposed stochastic dynamic damage constitutive model. Two joints are available in UDEC: a statistical joint generator that generates the joints using the parameters defined by conventional rock mechanics and a Voronoi tessellation generator that creates randomly sized polygonal blocks. The Voronoi tessellation generator was used to build the stochastic aggregate model in this study. To limit block size, the max average edge length of polygons was limited to 2 cm. Similar coarse aggregates were used at different positions of the lining structure; its physical parameters are shown in Table 2. Aggregate junctions at different positions were made up of cement, coal ash, sand, and water, with their physical parameters measured by the laboratory test (Table 3). The parameters for calculating the damage model were calibrated by uniaxial tensile and uniaxial compressive tests [20], as shown in Table 4.Table 1
Physical and mechanical parameters of rock blocks and joints.
TypeElastic modulus (GPa)Volume modulus (GPa)Shear modulus (GPa)Poisson’s ratioCohesion (MPa)Internal friction angle (°)Unit weight (kg/m3)Block20.6015.109.520.245.0045.0023.00JointN/AN/AN/AN/A0.505.00N/ATable 2
Physical and mechanical parameters of coarse aggregate.
MaterialDensity (kN/m3)Elastic modulus (GPa)Poisson’s ratioCohesion (MPa)Internal friction angle (°)Coarse aggregate23.0032.000.252.3055.00Table 3
Physical and mechanical parameters of joint of coarse aggregate in different positions.
Position of coarse aggregate junctionDensity (kN/m3)Elastic modulus (MPa)Poisson’s ratioCohesion (MPa)Internal friction angle (°)Primary lining19.7015.800.181.0422.00Secondary lining20.4016.100.181.3223.00Filled layer16.9014.300.181.0420.90Concrete foundation21.8016.600.181.2023.40Track slab23.8017.200.181.5224.50Table 4
Calculation parameters of damage model of coarse aggregate joint with different linings.
Position of coarse aggregate junctionf¯+ (MPa)f¯− (MPa)A+A−B+B−Secondary lining1.41315.0300.0000.5850.9900.153Filled layer1.14310.7100.0000.4500.7200.126Concrete foundation1.28712.8700.0000.5400.9000.135Track slab1.70120.7900.0000.9001.2600.162


### 5.3. Dynamic Water-Pressure Calculation within Seam

UDEC can be used to simulate both confined groundwater flow and free water flow. UDEC is provided with three fluid-flow-analysis calculation algorithms, basic, steady-state flow, and experimental algorithms, of which the steady-state flow algorithm enhances calculation speed since it ignores the influence of domain volume and its change on fluid pressure and, thus, does not consider the influence of fluid rigidity on time step. Therefore, the steady-state flow algorithm was used in this study. Before calculation, determination of the dynamic water pressure in the seam was required.According to an on-site survey, periodic interface clapping was observed in the water-bearing seam in the filled layer when the train approached and departed. Solid-liquid-gas three-phase coupled vibration was formed by the seam, groundwater, and air in the seam in the loading process of train load since the seam was connected with the side drain. To simplify the derivation, the following hypotheses were made: (1) water was an ideal incompressible fluid without viscosity; (2) the air in the seam was ideal gas; (3) the upper and lower bounds of the seam were impermeable; and (4) the water within the seam migrated along the direction of the tunnel cross section.


#### 5.3.1. Speed Distribution Equation of Water within Seam

As shown in Figure11, a fluid element with a length of dx and a height of ωx was taken at seam x. Since the boundary area of the fluid element was small, the speed along each boundary could be deemed to be uniformly distributed. The speed was υx at section 1-2 and υx+∂υx/∂xdx at section 3-4. Fluid speed at section 1–4 was υy (ωx) and consistent with the speed at the upper seam surface; section 2-3 was the solid-fluid boundary where fluid speed was 0. There was no shear stress but there was pressure for an ideal fluid. The pressure was px at section 1-2 and px+∂px/∂xdx at section 3-4.Figure 11
Fluid element in the crack.Within a microtimescale, the total massm of the net inflow control body was(66)m=ρυy−∂ρυx∂xωxdxdt.It could be inferred from the law of conservation of mass that the total mass of the net inflow control body was equal to the increment of fluid mass in the control body, namely,(67)ρυy−∂ρυx∂xωxdxdt=h∂ρ∂tdxdt.The speedυy at the upper seam surface could be obtained by taking the derivate of the displacement of the upper seam surface with respect to time:(68)υy=ux′=βxdFdt.The densityρ of incompressible fluid was a constant. Equation (67) could be reduced by substituting (68) into it, thus obtaining(69)∂υx∂x=βxωxdFdt.Equation (69) is the differential equation of water-speed distribution in the seam with the action of train load according to the law of conservation of mass.


#### 5.3.2. Pressure Distribution Differential Equation of Water within Seam

As shown in Figure11, the resultant force along the direction of fluid element x is(70)∑Fx=−∂px∂xωxdx.According to Newton’s second law of motion (F = ma),(71)−∂px∂xωxdx=ρωxdxdυxdt,along the direction of x, which could be reduced to(72)∂px∂x=−ρdυxdt.Equation (72) is the differential equation of water-pressure distribution in the seam with the action of train load based on Newton’s second law of motion.


#### 5.3.3. Boundary Conditions

The enlarged crack-deformation diagram is shown in Figure12, where h is the initial crack opening, ωx is the crack opening at x after loading (for simplicity, the subscript t was omitted since it meant the fluid state at any time in this study; the same as below), ux is the displacement of upper seam surface at x after loading, Lai is the air depth within the seam with unsaturated water at the initial moment, and Fx is the imposed load.Figure 12
Schematic diagram of crack deformation under train load.Since the track slab deformed slightly within the elastic range with the effect of train load, the displacementux of the upper seam surface can be denoted as(73)ux=βxF,where βx is the deformation coefficient associated with seam depth and the acting position of load.It can be seen from Figure11 that the initial volume Vai of the air domain within the seam per unit length at the initial moment is(74)Vai=Laih.Owing to the balance between air and water within the seam and the absolute pressurepai at the outlet of the seam, the initial pressure of air and water within the seam is equal to the absolute pressure pai at the outlet of the seam.The interface between air and water shifts fromLai to Lap at time t. According to the incompressibility of water, volumes 1 and 2 are equal; thus,(75)∫LapLaiωxdx=∫LaiLstauxdx.The volume variationΔV of the air domain at this time is(76)ΔV=∫0Lstauxdx=F∫0Lstaβxdx.It can be seen from the ideal-gas state equation that the product of gas pressure P and its volume V is a constant when the amount and temperature of gas remain unchanged, namely,(77)pV=C.It could be derived from (74), (76), and (77) that the pressure pair in the action of train load at time t is(78)pair=paiVaiVai+ΔV=paiLaihLaih+F∫0Lstaβxdx.The pressurepap at the air-water interface is equal to that of the air domain at this point; thus,(79)pap=pair.At stationary pointx=Lsta,(80)υxLsta=0.Since the pressurepap at the air-water interface within the seam is equal to that of the air domain pair,(81)pxLap=pap=pair.The pressure at the seam outlet is equal to the absolute pressure of liquid, namely,(82)pxL=pai.


#### 5.3.4. Analytical Expression of Pressure Distribution within Seam

The air domain pressure in areax=0∼Lap is(83)pair=paiLaihLaih+F∫0Lstaβxdx.The pressure distribution expression of water in areax=Lap∼Lsta of the seam could be obtained based on the differential equation of water current speed within the seam (69), the differential equation of water-pressure distribution (72), and the boundary-condition expressions (80) and (81):(84)px=ρF∫Lapxdx∫xLstaβxωxdx+pair.Similarly, the pressure distribution expression of water in areax=Lsta∼L within the seam could be obtained based on (82):(85)px=−ρF∫xLdx∫Lstaxβxωxdx+pai.When the initial volume of air in the seam was determined, the water-pressure distribution in different areas of the seam could be obtained based on equations (83)–(85). According to calculation, the proportion of the initial air volume in the total seam volume is 5%; Lsta = 1.5 m and L = 4.7 m; βx is the vertical displacement of point A in the numerical calculation model, and pai is the standard atmospheric pressure, 101.325 kPa.


### 5.4. Calculation Conditions and Steps

The train speed was set to 300 km/h, and the calculation conditions included the following: (1) there was no seam in the filled layer; (2) there was a seam with groundwater in the filled layer; and (3) there was a seam without groundwater in the filled layer. Dynamic water pressure was not considered in conditions 1 and 2.Calculation steps were as follows: modeling⟶meshing⟶static boundary definition (horizontal displacement at both sides and vertical displacement at the bottom were fixed)⟶excavation (by upper and lower benches)⟶support ⟶dynamic water-pressure addition (if any)⟶dynamic boundary⟶dynamic calculation.


### 5.5. Dynamic Response of Tunnel-Base Concrete Structure

The amplitudes of vertical dynamic stress and vibration accelerations of the five monitoring points in Figure10 in the three working conditions, as well as the comparison between them, are shown in Figures 13 and 14. The trends of vertical dynamic stress and vibration accelerations varying with depth below the top surface of the filled layer in the three working conditions are shown in Figures 15 and 16. With the values of point A (at a depth of 0.2 m) shown in Figure 10 taken as the reference values in the three working conditions, the amplitude attenuation of vertical dynamic stress and vibration accelerations is given in Table 5. It can be seen that the vertical dynamic stress and vertical vibration acceleration decreased with depth in all three working conditions, but the presence of a seam in the filled layer in working conditions 2 and 3 led to a sharp decrease of vertical dynamic stress and vertical vibration acceleration. Compared with the upper seam surface (at a depth of 0.2 m), the attenuation values of vertical dynamic stress and vertical vibration acceleration at the monitoring point on the lower seam surface (at a depth of 0.22 m) were (41.07%; 47.13%) and (91.17%; 91.73%), respectively. Such attenuation values were far larger than the vertical dynamic stress attenuation (1.55%) and vertical vibration acceleration attenuation (0.00%) in working condition 1.Figure 13
Comparison of vertical dynamic stress amplitude at different depth monitoring points under three conditions.Figure 14
Comparison of vertical vibration acceleration amplitude at different depth monitoring points under three conditions.Figure 15
Variation of dynamic stress amplitude at different depths in the three conditions.Figure 16
Variation trend of vertical vibration acceleration at different depths in the three conditions.Table 5
Comparison of vertical dynamic stress and acceleration attenuation at different depth-monitoring points under three conditions.
Monitoring pointDepth (m)Attenuation of vertical dynamic stress amplitude (%)Attenuation of vertical acceleration amplitude (%)Working condition 1Working condition 2Working condition 3Working condition 1Working condition 2Working condition 3B0.221.5541.0747.130.0091.1791.73C0.5220.1555.7763.9724.5092.7693.94D1.12238.8373.4181.4337.7994.2795.34E1.73349.5683.8690.9249.5095.2496.30As shown in Table5, the dynamic stress at point E on the invert bottom in working conditions 2 and 3 was attenuated by 90.92% and 83.86%, respectively. Thus, the upward pressure applied by surrounding rocks could not be effectively reduced, which went against the overall stress of the tunnel lining structure.According to the comparison between Figures12 and 14, the seam had little influence on vertical dynamic stress transmission from the top surface of the filled layer to point A. Compared with the case in working condition 1, vertical dynamic stress at point A increased by 1.15% in working condition 2 and by 1.24% in working condition 3 because the concrete on both sides of the seam was not completely separated. The propagation rule also remained the same when contact area and load were invariant.The comparison between Figures13 and 15 indicates that the seam has significant influence on vertical vibration acceleration from the top surface of the filled layer to point A. Compared with the case in working condition 1, vertical vibration acceleration at point A increased by 724.67% in working condition 2 and by 765.02% in working condition 3 because the mass of the concrete structure involved in forced vibration decreased in the presence of a seam. Thus, acceleration was bound to increase when load remained the same.


### 5.6. Damage Fatigue Properties of Concrete within Filled Layer of Tunnel Base

As shown in Table5, the vertical dynamic stress attenuation at point B on the lower seam surface compared with point A was lower in working condition 2 than in working condition 3. It followed that the groundwater in the seam was conducive to train load transmission in the tunnel bottom, but groundwater migration within the seam with the action of train load would inevitably cause abrasion damage to concrete aggregate junctions on the upper and lower seam surfaces. Thus, the damage amounts at points A and B were the heaviest. Therefore, the damage at both points should be calculated; the damage fatigue life of the point with greater damage should be considered.Figure 17
Damage points of A and B caused by train load at the 1 time of the three conditions.Damage amountsDvibration + at points A and B with the action of train load for once in three working conditions are shown in Figure 17. It can be seen that the damage amount at point A was greater in the three working conditions; it turned out to be the greatest and reached 0.4436992 × 10−6 when there was a seam with groundwater in the filled layer. It was far smaller than 1, suggesting that the damage amount resulting from the single pass of the train was too tiny to noticeably influence the concrete structure safety.However, the total traffic flow can reach 365 × 104 times if there were 100 16-wagon high-speed trains passing through the tunnel every day during its 100-year service period, which is a high-cycle frequency. Thus, damage-fatigue-life analysis is required. To facilitate comparison with similar studies, the simulation was carried out at an initial damage amount of 0.0501300 [2].(86)N=Dtotal−DinitialDvibration+×1365×100.In addition, the formula was calculated based on Palmgren–Miner cumulative fatigue-damage theory, whereN denotes damage fatigue life (in years); Dtotal is the critical value of cumulative concrete damage, which is set to 0.9500000 based on current structure damage diagnosis and maintenance and reinforcement; D初 is the initial structure damage amount; and D振+ is the damage amount at point A, as shown in Figure 16.The fatigue life (rounded to the closest integer; see Table6) at different train speeds could be obtained using (63). As shown in Table 6, the relative difference between the result calculated using the proposed constitutive model in working condition 1 and that obtained by Du et al. [1] is 2.08%. This not only demonstrates the applicability of the proposed constitutive model but also indicates that current high-speed railway tunnel parameters satisfied the antifatigue requirements within the design reference period when there was no seam in the tunnel base with IV-class surrounding rocks, and train speed was less than or equal to 300 km/h. In working conditions 2 and 3, however, damage fatigue life did not satisfy the antifatigue requirements within the design reference period; the minimum service life was only 56 years; and the presence of groundwater in the seam would aggravate structure damage.Table 6
Comparison of fatigue life under different working conditions.
Working conditionDamage fatigue life (years)Relative difference (%)Result obtained using proposed constitutive modelResult calculated by Huang [2]Working condition 1There is no seam in the filled layer2402493.61Working condition 2There is a seam with groundwater in the filled layer56N/AN/AWorking condition 3There is a seam without groundwater in the filled layer62N/AN/AIt is worth noting that structural failure was not only caused by damage but was also highly related to the characteristics of surrounding rocks. In the engineering case investigated by Du et al. [1], the bottom structure would be damaged soon once a seam occurred in the filled layer since the surrounding rocks were swelling rocks, which resulted in rapid bottom deformation and affected traffic safety. Therefore, countermeasures must be taken immediately. Considering that structure failure and service life are subject to the combined influence of various factors, the result obtained in this study is presented for reference only.


## 5.1. Engineering Overview

A double-track high-speed railway tunnel located in Western China was taken as the research object. The formation in the tunnel-site area was composed of IV-class surrounding rocks with joint fissures. Average joint spacing and max depth of the tunnel were 0.5 and 270 m, respectively. The tunnel section is shown in Figure6.Figure 6
Cross section of tunnel lining (unit: mm).According to the periodic reinspection by the railway engineering department using a track-inspection car during tunnel operation, a 1.268-km-long track bed slab bulge was observed along k1679 + 139–+562 of the upline and k1680 + 044–+889 of the downline. Compared with the design values, the max difference of elevation reached 12 cm. An on-site drill-hole survey suggested seam defects (see Figure7) between the base filled layer and the leveling layer, with a max seam gap of 1.2–2.6 cm. There was also high-pressure groundwater, for which groundwater upwelled quickly to the track bed after drilling.Figure 7
Crack in the fill layer.


## 5.2. UDEC Numerical Calculation Model and Material Parameters

To analyze the law of dynamic response and damage characteristics of a base filled layer with seam defects, numerical models of the filled layer with and without a seam were built based on UDEC. The seam was 2 cm wide and 20 cm from the top surface of the filled layer. The space at the left, right, top, and bottom of the model was 4 times the tunnel span. In this paper, the viscous boundary is adopted in the calculation model. The deficiency of model depth was compensated by load. The local enlarged views of the models are shown in Figures8 and 9. Train load depends on axle weight, suspended mass, operating speed, track structure, track regularity, and so on. In this study, the loading method developed by Liang et al. [18] was used. Train load was applied along the left line during simulation. As shown in Figure 10, five monitoring points, A, B, C, D, and E, were arranged at depths of 20, 22, 52, 112.2, and 173.3 cm, respectively, below the top surface of the filled layer along the center line of the left rail line, points A and B of which were located on the upper interface of the seam; D was on the interface between the invert and filled layer; and E was on the interface between the primary support of the invert and surrounding rock.Figure 8
Local mesh division of the numerical model for the filled layer with no crack.Figure 9
Local mesh division of the numerical model for the filled layer with a crack.Figure 10
Schematic diagram of monitoring points (Unit: mm).The yield of blocks in the formation satisfied the Mohr–Coulomb criterion, while the joint yield conformed to the Coulomb sliding criterion of surface contact (see Wang [19]). According to the data specified in the Handbook for Railway Engineering Geology and Code for Design of Railway Tunnel, the physical parameters of formation blocks are shown in Table 1. Supposing that the concrete structure was comprised of coarse aggregates (mean particle size of 2 cm) that combined each other by “corner-corner” contact, “corner-side” contact, and “side-side” contact, the yield of coarse aggregates satisfied the Drucker–Prager criterion. The coarse aggregate interface was simulated using the proposed stochastic dynamic damage constitutive model. Two joints are available in UDEC: a statistical joint generator that generates the joints using the parameters defined by conventional rock mechanics and a Voronoi tessellation generator that creates randomly sized polygonal blocks. The Voronoi tessellation generator was used to build the stochastic aggregate model in this study. To limit block size, the max average edge length of polygons was limited to 2 cm. Similar coarse aggregates were used at different positions of the lining structure; its physical parameters are shown in Table 2. Aggregate junctions at different positions were made up of cement, coal ash, sand, and water, with their physical parameters measured by the laboratory test (Table 3). The parameters for calculating the damage model were calibrated by uniaxial tensile and uniaxial compressive tests [20], as shown in Table 4.Table 1
Physical and mechanical parameters of rock blocks and joints.
TypeElastic modulus (GPa)Volume modulus (GPa)Shear modulus (GPa)Poisson’s ratioCohesion (MPa)Internal friction angle (°)Unit weight (kg/m3)Block20.6015.109.520.245.0045.0023.00JointN/AN/AN/AN/A0.505.00N/ATable 2
Physical and mechanical parameters of coarse aggregate.
MaterialDensity (kN/m3)Elastic modulus (GPa)Poisson’s ratioCohesion (MPa)Internal friction angle (°)Coarse aggregate23.0032.000.252.3055.00Table 3
Physical and mechanical parameters of joint of coarse aggregate in different positions.
Position of coarse aggregate junctionDensity (kN/m3)Elastic modulus (MPa)Poisson’s ratioCohesion (MPa)Internal friction angle (°)Primary lining19.7015.800.181.0422.00Secondary lining20.4016.100.181.3223.00Filled layer16.9014.300.181.0420.90Concrete foundation21.8016.600.181.2023.40Track slab23.8017.200.181.5224.50Table 4
Calculation parameters of damage model of coarse aggregate joint with different linings.
Position of coarse aggregate junctionf¯+ (MPa)f¯− (MPa)A+A−B+B−Secondary lining1.41315.0300.0000.5850.9900.153Filled layer1.14310.7100.0000.4500.7200.126Concrete foundation1.28712.8700.0000.5400.9000.135Track slab1.70120.7900.0000.9001.2600.162


## 5.3. Dynamic Water-Pressure Calculation within Seam

UDEC can be used to simulate both confined groundwater flow and free water flow. UDEC is provided with three fluid-flow-analysis calculation algorithms, basic, steady-state flow, and experimental algorithms, of which the steady-state flow algorithm enhances calculation speed since it ignores the influence of domain volume and its change on fluid pressure and, thus, does not consider the influence of fluid rigidity on time step. Therefore, the steady-state flow algorithm was used in this study. Before calculation, determination of the dynamic water pressure in the seam was required.According to an on-site survey, periodic interface clapping was observed in the water-bearing seam in the filled layer when the train approached and departed. Solid-liquid-gas three-phase coupled vibration was formed by the seam, groundwater, and air in the seam in the loading process of train load since the seam was connected with the side drain. To simplify the derivation, the following hypotheses were made: (1) water was an ideal incompressible fluid without viscosity; (2) the air in the seam was ideal gas; (3) the upper and lower bounds of the seam were impermeable; and (4) the water within the seam migrated along the direction of the tunnel cross section.


### 5.3.1. Speed Distribution Equation of Water within Seam

As shown in Figure11, a fluid element with a length of dx and a height of ωx was taken at seam x. Since the boundary area of the fluid element was small, the speed along each boundary could be deemed to be uniformly distributed. The speed was υx at section 1-2 and υx+∂υx/∂xdx at section 3-4. Fluid speed at section 1–4 was υy (ωx) and consistent with the speed at the upper seam surface; section 2-3 was the solid-fluid boundary where fluid speed was 0. There was no shear stress but there was pressure for an ideal fluid. The pressure was px at section 1-2 and px+∂px/∂xdx at section 3-4.Figure 11
Fluid element in the crack.Within a microtimescale, the total massm of the net inflow control body was(66)m=ρυy−∂ρυx∂xωxdxdt.It could be inferred from the law of conservation of mass that the total mass of the net inflow control body was equal to the increment of fluid mass in the control body, namely,(67)ρυy−∂ρυx∂xωxdxdt=h∂ρ∂tdxdt.The speedυy at the upper seam surface could be obtained by taking the derivate of the displacement of the upper seam surface with respect to time:(68)υy=ux′=βxdFdt.The densityρ of incompressible fluid was a constant. Equation (67) could be reduced by substituting (68) into it, thus obtaining(69)∂υx∂x=βxωxdFdt.Equation (69) is the differential equation of water-speed distribution in the seam with the action of train load according to the law of conservation of mass.


### 5.3.2. Pressure Distribution Differential Equation of Water within Seam

As shown in Figure11, the resultant force along the direction of fluid element x is(70)∑Fx=−∂px∂xωxdx.According to Newton’s second law of motion (F = ma),(71)−∂px∂xωxdx=ρωxdxdυxdt,along the direction of x, which could be reduced to(72)∂px∂x=−ρdυxdt.Equation (72) is the differential equation of water-pressure distribution in the seam with the action of train load based on Newton’s second law of motion.


### 5.3.3. Boundary Conditions

The enlarged crack-deformation diagram is shown in Figure12, where h is the initial crack opening, ωx is the crack opening at x after loading (for simplicity, the subscript t was omitted since it meant the fluid state at any time in this study; the same as below), ux is the displacement of upper seam surface at x after loading, Lai is the air depth within the seam with unsaturated water at the initial moment, and Fx is the imposed load.Figure 12
Schematic diagram of crack deformation under train load.Since the track slab deformed slightly within the elastic range with the effect of train load, the displacementux of the upper seam surface can be denoted as(73)ux=βxF,where βx is the deformation coefficient associated with seam depth and the acting position of load.It can be seen from Figure11 that the initial volume Vai of the air domain within the seam per unit length at the initial moment is(74)Vai=Laih.Owing to the balance between air and water within the seam and the absolute pressurepai at the outlet of the seam, the initial pressure of air and water within the seam is equal to the absolute pressure pai at the outlet of the seam.The interface between air and water shifts fromLai to Lap at time t. According to the incompressibility of water, volumes 1 and 2 are equal; thus,(75)∫LapLaiωxdx=∫LaiLstauxdx.The volume variationΔV of the air domain at this time is(76)ΔV=∫0Lstauxdx=F∫0Lstaβxdx.It can be seen from the ideal-gas state equation that the product of gas pressure P and its volume V is a constant when the amount and temperature of gas remain unchanged, namely,(77)pV=C.It could be derived from (74), (76), and (77) that the pressure pair in the action of train load at time t is(78)pair=paiVaiVai+ΔV=paiLaihLaih+F∫0Lstaβxdx.The pressurepap at the air-water interface is equal to that of the air domain at this point; thus,(79)pap=pair.At stationary pointx=Lsta,(80)υxLsta=0.Since the pressurepap at the air-water interface within the seam is equal to that of the air domain pair,(81)pxLap=pap=pair.The pressure at the seam outlet is equal to the absolute pressure of liquid, namely,(82)pxL=pai.


### 5.3.4. Analytical Expression of Pressure Distribution within Seam

The air domain pressure in areax=0∼Lap is(83)pair=paiLaihLaih+F∫0Lstaβxdx.The pressure distribution expression of water in areax=Lap∼Lsta of the seam could be obtained based on the differential equation of water current speed within the seam (69), the differential equation of water-pressure distribution (72), and the boundary-condition expressions (80) and (81):(84)px=ρF∫Lapxdx∫xLstaβxωxdx+pair.Similarly, the pressure distribution expression of water in areax=Lsta∼L within the seam could be obtained based on (82):(85)px=−ρF∫xLdx∫Lstaxβxωxdx+pai.When the initial volume of air in the seam was determined, the water-pressure distribution in different areas of the seam could be obtained based on equations (83)–(85). According to calculation, the proportion of the initial air volume in the total seam volume is 5%; Lsta = 1.5 m and L = 4.7 m; βx is the vertical displacement of point A in the numerical calculation model, and pai is the standard atmospheric pressure, 101.325 kPa.


## 5.3.1. Speed Distribution Equation of Water within Seam

As shown in Figure11, a fluid element with a length of dx and a height of ωx was taken at seam x. Since the boundary area of the fluid element was small, the speed along each boundary could be deemed to be uniformly distributed. The speed was υx at section 1-2 and υx+∂υx/∂xdx at section 3-4. Fluid speed at section 1–4 was υy (ωx) and consistent with the speed at the upper seam surface; section 2-3 was the solid-fluid boundary where fluid speed was 0. There was no shear stress but there was pressure for an ideal fluid. The pressure was px at section 1-2 and px+∂px/∂xdx at section 3-4.Figure 11
Fluid element in the crack.Within a microtimescale, the total massm of the net inflow control body was(66)m=ρυy−∂ρυx∂xωxdxdt.It could be inferred from the law of conservation of mass that the total mass of the net inflow control body was equal to the increment of fluid mass in the control body, namely,(67)ρυy−∂ρυx∂xωxdxdt=h∂ρ∂tdxdt.The speedυy at the upper seam surface could be obtained by taking the derivate of the displacement of the upper seam surface with respect to time:(68)υy=ux′=βxdFdt.The densityρ of incompressible fluid was a constant. Equation (67) could be reduced by substituting (68) into it, thus obtaining(69)∂υx∂x=βxωxdFdt.Equation (69) is the differential equation of water-speed distribution in the seam with the action of train load according to the law of conservation of mass.


## 5.3.2. Pressure Distribution Differential Equation of Water within Seam

As shown in Figure11, the resultant force along the direction of fluid element x is(70)∑Fx=−∂px∂xωxdx.According to Newton’s second law of motion (F = ma),(71)−∂px∂xωxdx=ρωxdxdυxdt,along the direction of x, which could be reduced to(72)∂px∂x=−ρdυxdt.Equation (72) is the differential equation of water-pressure distribution in the seam with the action of train load based on Newton’s second law of motion.


## 5.3.3. Boundary Conditions

The enlarged crack-deformation diagram is shown in Figure12, where h is the initial crack opening, ωx is the crack opening at x after loading (for simplicity, the subscript t was omitted since it meant the fluid state at any time in this study; the same as below), ux is the displacement of upper seam surface at x after loading, Lai is the air depth within the seam with unsaturated water at the initial moment, and Fx is the imposed load.Figure 12
Schematic diagram of crack deformation under train load.Since the track slab deformed slightly within the elastic range with the effect of train load, the displacementux of the upper seam surface can be denoted as(73)ux=βxF,where βx is the deformation coefficient associated with seam depth and the acting position of load.It can be seen from Figure11 that the initial volume Vai of the air domain within the seam per unit length at the initial moment is(74)Vai=Laih.Owing to the balance between air and water within the seam and the absolute pressurepai at the outlet of the seam, the initial pressure of air and water within the seam is equal to the absolute pressure pai at the outlet of the seam.The interface between air and water shifts fromLai to Lap at time t. According to the incompressibility of water, volumes 1 and 2 are equal; thus,(75)∫LapLaiωxdx=∫LaiLstauxdx.The volume variationΔV of the air domain at this time is(76)ΔV=∫0Lstauxdx=F∫0Lstaβxdx.It can be seen from the ideal-gas state equation that the product of gas pressure P and its volume V is a constant when the amount and temperature of gas remain unchanged, namely,(77)pV=C.It could be derived from (74), (76), and (77) that the pressure pair in the action of train load at time t is(78)pair=paiVaiVai+ΔV=paiLaihLaih+F∫0Lstaβxdx.The pressurepap at the air-water interface is equal to that of the air domain at this point; thus,(79)pap=pair.At stationary pointx=Lsta,(80)υxLsta=0.Since the pressurepap at the air-water interface within the seam is equal to that of the air domain pair,(81)pxLap=pap=pair.The pressure at the seam outlet is equal to the absolute pressure of liquid, namely,(82)pxL=pai.


## 5.3.4. Analytical Expression of Pressure Distribution within Seam

The air domain pressure in areax=0∼Lap is(83)pair=paiLaihLaih+F∫0Lstaβxdx.The pressure distribution expression of water in areax=Lap∼Lsta of the seam could be obtained based on the differential equation of water current speed within the seam (69), the differential equation of water-pressure distribution (72), and the boundary-condition expressions (80) and (81):(84)px=ρF∫Lapxdx∫xLstaβxωxdx+pair.Similarly, the pressure distribution expression of water in areax=Lsta∼L within the seam could be obtained based on (82):(85)px=−ρF∫xLdx∫Lstaxβxωxdx+pai.When the initial volume of air in the seam was determined, the water-pressure distribution in different areas of the seam could be obtained based on equations (83)–(85). According to calculation, the proportion of the initial air volume in the total seam volume is 5%; Lsta = 1.5 m and L = 4.7 m; βx is the vertical displacement of point A in the numerical calculation model, and pai is the standard atmospheric pressure, 101.325 kPa.


## 5.4. Calculation Conditions and Steps

The train speed was set to 300 km/h, and the calculation conditions included the following: (1) there was no seam in the filled layer; (2) there was a seam with groundwater in the filled layer; and (3) there was a seam without groundwater in the filled layer. Dynamic water pressure was not considered in conditions 1 and 2.Calculation steps were as follows: modeling⟶meshing⟶static boundary definition (horizontal displacement at both sides and vertical displacement at the bottom were fixed)⟶excavation (by upper and lower benches)⟶support ⟶dynamic water-pressure addition (if any)⟶dynamic boundary⟶dynamic calculation.


## 5.5. Dynamic Response of Tunnel-Base Concrete Structure

The amplitudes of vertical dynamic stress and vibration accelerations of the five monitoring points in Figure10 in the three working conditions, as well as the comparison between them, are shown in Figures 13 and 14. The trends of vertical dynamic stress and vibration accelerations varying with depth below the top surface of the filled layer in the three working conditions are shown in Figures 15 and 16. With the values of point A (at a depth of 0.2 m) shown in Figure 10 taken as the reference values in the three working conditions, the amplitude attenuation of vertical dynamic stress and vibration accelerations is given in Table 5. It can be seen that the vertical dynamic stress and vertical vibration acceleration decreased with depth in all three working conditions, but the presence of a seam in the filled layer in working conditions 2 and 3 led to a sharp decrease of vertical dynamic stress and vertical vibration acceleration. Compared with the upper seam surface (at a depth of 0.2 m), the attenuation values of vertical dynamic stress and vertical vibration acceleration at the monitoring point on the lower seam surface (at a depth of 0.22 m) were (41.07%; 47.13%) and (91.17%; 91.73%), respectively. Such attenuation values were far larger than the vertical dynamic stress attenuation (1.55%) and vertical vibration acceleration attenuation (0.00%) in working condition 1.Figure 13
Comparison of vertical dynamic stress amplitude at different depth monitoring points under three conditions.Figure 14
Comparison of vertical vibration acceleration amplitude at different depth monitoring points under three conditions.Figure 15
Variation of dynamic stress amplitude at different depths in the three conditions.Figure 16
Variation trend of vertical vibration acceleration at different depths in the three conditions.Table 5
Comparison of vertical dynamic stress and acceleration attenuation at different depth-monitoring points under three conditions.
Monitoring pointDepth (m)Attenuation of vertical dynamic stress amplitude (%)Attenuation of vertical acceleration amplitude (%)Working condition 1Working condition 2Working condition 3Working condition 1Working condition 2Working condition 3B0.221.5541.0747.130.0091.1791.73C0.5220.1555.7763.9724.5092.7693.94D1.12238.8373.4181.4337.7994.2795.34E1.73349.5683.8690.9249.5095.2496.30As shown in Table5, the dynamic stress at point E on the invert bottom in working conditions 2 and 3 was attenuated by 90.92% and 83.86%, respectively. Thus, the upward pressure applied by surrounding rocks could not be effectively reduced, which went against the overall stress of the tunnel lining structure.According to the comparison between Figures12 and 14, the seam had little influence on vertical dynamic stress transmission from the top surface of the filled layer to point A. Compared with the case in working condition 1, vertical dynamic stress at point A increased by 1.15% in working condition 2 and by 1.24% in working condition 3 because the concrete on both sides of the seam was not completely separated. The propagation rule also remained the same when contact area and load were invariant.The comparison between Figures13 and 15 indicates that the seam has significant influence on vertical vibration acceleration from the top surface of the filled layer to point A. Compared with the case in working condition 1, vertical vibration acceleration at point A increased by 724.67% in working condition 2 and by 765.02% in working condition 3 because the mass of the concrete structure involved in forced vibration decreased in the presence of a seam. Thus, acceleration was bound to increase when load remained the same.


## 5.6. Damage Fatigue Properties of Concrete within Filled Layer of Tunnel Base

As shown in Table5, the vertical dynamic stress attenuation at point B on the lower seam surface compared with point A was lower in working condition 2 than in working condition 3. It followed that the groundwater in the seam was conducive to train load transmission in the tunnel bottom, but groundwater migration within the seam with the action of train load would inevitably cause abrasion damage to concrete aggregate junctions on the upper and lower seam surfaces. Thus, the damage amounts at points A and B were the heaviest. Therefore, the damage at both points should be calculated; the damage fatigue life of the point with greater damage should be considered.Figure 17
Damage points of A and B caused by train load at the 1 time of the three conditions.Damage amountsDvibration + at points A and B with the action of train load for once in three working conditions are shown in Figure 17. It can be seen that the damage amount at point A was greater in the three working conditions; it turned out to be the greatest and reached 0.4436992 × 10−6 when there was a seam with groundwater in the filled layer. It was far smaller than 1, suggesting that the damage amount resulting from the single pass of the train was too tiny to noticeably influence the concrete structure safety.However, the total traffic flow can reach 365 × 104 times if there were 100 16-wagon high-speed trains passing through the tunnel every day during its 100-year service period, which is a high-cycle frequency. Thus, damage-fatigue-life analysis is required. To facilitate comparison with similar studies, the simulation was carried out at an initial damage amount of 0.0501300 [2].(86)N=Dtotal−DinitialDvibration+×1365×100.In addition, the formula was calculated based on Palmgren–Miner cumulative fatigue-damage theory, whereN denotes damage fatigue life (in years); Dtotal is the critical value of cumulative concrete damage, which is set to 0.9500000 based on current structure damage diagnosis and maintenance and reinforcement; D初 is the initial structure damage amount; and D振+ is the damage amount at point A, as shown in Figure 16.The fatigue life (rounded to the closest integer; see Table6) at different train speeds could be obtained using (63). As shown in Table 6, the relative difference between the result calculated using the proposed constitutive model in working condition 1 and that obtained by Du et al. [1] is 2.08%. This not only demonstrates the applicability of the proposed constitutive model but also indicates that current high-speed railway tunnel parameters satisfied the antifatigue requirements within the design reference period when there was no seam in the tunnel base with IV-class surrounding rocks, and train speed was less than or equal to 300 km/h. In working conditions 2 and 3, however, damage fatigue life did not satisfy the antifatigue requirements within the design reference period; the minimum service life was only 56 years; and the presence of groundwater in the seam would aggravate structure damage.Table 6
Comparison of fatigue life under different working conditions.
Working conditionDamage fatigue life (years)Relative difference (%)Result obtained using proposed constitutive modelResult calculated by Huang [2]Working condition 1There is no seam in the filled layer2402493.61Working condition 2There is a seam with groundwater in the filled layer56N/AN/AWorking condition 3There is a seam without groundwater in the filled layer62N/AN/AIt is worth noting that structural failure was not only caused by damage but was also highly related to the characteristics of surrounding rocks. In the engineering case investigated by Du et al. [1], the bottom structure would be damaged soon once a seam occurred in the filled layer since the surrounding rocks were swelling rocks, which resulted in rapid bottom deformation and affected traffic safety. Therefore, countermeasures must be taken immediately. Considering that structure failure and service life are subject to the combined influence of various factors, the result obtained in this study is presented for reference only.


## 6. Conclusions

(1)
A damage-evolution equation is introduced based on the law of thermodynamics and energy-dissipation principle. The elastoplastic Helmholtz free energy is corrected by introducing the hardening parameter. A corrected elastoplastic Helmholtz free-energy-based stochastic dynamic damage constitutive model is derived by means of elastoplastic decoupling and dynamic expansion of the damage energy-release rate.(2)
The elastic prediction-plastic modification-damage correction numerical analysis and solution framework for the proposed model is established. Programmed calculation of the proposed constitutive model is implemented based on UDEC. The validity of the proposed model is verified by comparisons with the laboratory uniaxial rapid loading test, which also suggests that the proposed constitutive model is able to accurately simulate the strain softening, stochasticity, nonlinearity, and strain-rate effect of the mechanical behavior of concrete.(3)
Tunnel-base dynamic response in three working conditions—namely, there is no seam in the filled layer, there is a seam with water in the filled layer, and there is a seam without water in the filled layer—is analyzed based on engineering application, which indicates that dynamic stress and vertical vibration acceleration decrease with depth.(4)
In the working conditions with a seam, vertical dynamic stress and vertical vibration acceleration at both sides of the seam decline sharply by (41.07%; 47.13%) and (91.17%; 91.73%) with and without groundwater, respectively. Both values are far higher than the vertical dynamic stress attenuation (1.55%) and vertical vibration acceleration attenuation (0.00%) in the working condition without a seam. Vertical dynamic stress is less affected by the seam within the 20 cm range from the top surface of the filled layer to the upper seam surface, whereas the seam has a more significant influence on vertical vibration acceleration within the same range. In the two working conditions with and without groundwater, vertical dynamic stress increases by 1.15% and 1.24%, while vertical vibration acceleration increases by 724.67% and 765.02%, respectively.(5)
Groundwater in the seam is beneficial for train-load transmission in the tunnel base, but it is inevitable that the groundwater migration in the seam causes abrasion damage to the concrete aggregate junctions on the upper and lower seam surfaces, which aggravates accumulative damage. The damage fatigue life of the lining structure of the high-speed railway tunnel is preliminarily analyzed using the proposed model without considering crustal stress and special surrounding rocks. According to the results of the analysis, the current high-speed railway tunnel parameters can satisfy the antifatigue requirements within the design reference period when there is no seam in the tunnel base with IV-class surrounding rocks, and train speed is 300 km/h. When there is a seam in the tunnel base, however, antifatigue life is 56 years in the presence of groundwater and 62 years without groundwater, which suggests that current design parameters fail to satisfy the antifatigue requirements within the design reference period. Therefore, this type of defect should be given adequate consideration in order to identify deficiencies and repair them in time.



---

*Source: 6622998-2022-01-20.xml*