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.
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 aswhere is object density; is the differential of internal energy density with respect to time; 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 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:where ζ is the energy dissipation per unit volume; T is the function of empirical temperature, which is constantly positive; S is the entropy; 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: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 function T 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
Considering that the temperature function T was a constant, the differential form of (3) is
Thus, (4) and (5) could be combined:
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:
By substituting (8) into (7),could be derived.
In view of the stochasticity of the strain rate , the sufficient and necessary condition to satisfy the above inequality is
It could be derived from (10) that
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]):where is the differential of damage variable with respect to time, is the damage-evolution factor, and 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:wherewhere is the hardening parameter, “+” means tension, “－” means compression, is the initial tensile plastic free energy, and is the initial compressive plastic free energy.
The variable 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 , the first law of thermodynamics and the second law of thermodynamics could be expressed as
Since the laws of thermodynamics hold with any strain rate and the rate of temperature change, the necessary conditions for (20) and (21) are
Equation (22) could be equivalently transformed intowherewhere is the hardening function, fcs the initial yield stress of uniaxial compression, fc the uniaxial compression strength, the initial plastic accumulated strain, and and 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),
From the elastoplastic decoupling of the energy release rate,where and 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 using the following equation: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 was added when determining the plastic Helmholtz free energy, which means that ψ0p is defined as
Thus, (18) is converted into
The corrected plastic Helmholtz free energy could be obtained by combining equations (14)–(18), (19), and (29):
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
Owing to the stochasticity of ,
If equation (31) holds, the expressioncould be obtained by substituting (14) into (32).
The elastoplastic stochastic damage constitutive relation could be derived by substituting (15) into (33):wherewhere 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 iswhere is the max principal stress, is the min principal stress, is the initial elastic modulus of materials, 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:where 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; 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 bywhere is the yield function of the effective stress space and is the differential of the plastic multiplier .
Considering that the energy dissipation of the damage system is generated by two mechanisms, the damage energy dissipation and plastic dissipation , the total energy dissipation iswhere is the damage strain rate and the plastic strain rate.
The energy dissipation in the damage process is subject to the restriction of the damage criterion G (σ) = 0. According to the principle of the maximum energy-dissipation extreme, the conditions for taking the extremum should satisfywhere is the Lagrange multiplier and G is the damage criterion.
Substituting (43) into (44),where 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 σ,could be derived if D was fixed, i.e., D was independent of σ, wherewhere ß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 variable D varying with the other internal variables. The evolution rule of the damage variable isaccording 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.
The expression was obtained by dynamic expansion of the damage energy release rate (see Li et al. [7]), where is the dynamic damage energy release rate, and are dynamic effect empirical parameters, and = 2.1 × 109 N/s, = 6 × 1010 N/s, = 5.5, and = 4.5 according to Liu [12].
Equation (51) is a first-order equation regarding . The dynamic damage-evolution process under tension and compression could be solved using (52) and (53), respectively: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
When G±<0, and D± = 0; the concrete was in the damage unloading phase or neutral loading phase. When , G± = 0. was obtained based on the damage consistency conditions:where is the differential form of the damage criterion.
Considering the initial condition D0± = 0, it was obtained by substituting (55) into (49).
Therefore, the damage variable could be obtained by (56) once the expression of function G± is given: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 α iswhere 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.
The principle of this process is shown below. Letting the strain be , stress , plastic strain , and plastic hardening parameter at previous increment step tn, the predictive variable of stress , that of stress-strain , that of plastic hardening parameter , and that of damage at the next increment step tn+1 should be
Equation (60) is substituted into yield function . If the yield function satisfies (61) at the time tn+1,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 usingwhere is the increment of plastic flow factor; is the hardening-parameter modulus; and , , , and are the effective stress, damage, effective strain, and state variable of the hardening function after plastic correction, respectively.
The effective stress tensor was divided into tensile and compressive parts after plastic correction. The new damage statecould be obtained by substituting them into (51), the tensile-compressive damage energy release formula, where and 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, , , , and 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].
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:where is the joint probability density function of Z and is the partial differential of Z with respect to .
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 in the effective stress space involving the numerical simulation is taken as the linear limit strengthening of materials; = 2.04 MPa and = 27.5 MPa. With calibration, the other model parameters are = 0.0, = 1.0, = 0.20, and = 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.
As shown in Figures 3–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 Figure 6.
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 Figure 7) 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.
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 Figures 8 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.
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.
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 Figure 11, a fluid element with a length of and a height of was taken at seam . 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 at section 1-2 and at section 3-4. Fluid speed at section 1–4 was () 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 at section 1-2 and at section 3-4.
Within a microtimescale, the total mass of the net inflow control body was
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,
The speed at the upper seam surface could be obtained by taking the derivate of the displacement of the upper seam surface with respect to time:
The density of incompressible fluid was a constant. Equation (67) could be reduced by substituting (68) into it, thus obtaining
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 Figure 11, the resultant force along the direction of fluid element x is
According to Newton’s second law of motion (F = ma),along the direction of x, which could be reduced to
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 Figure 12, where is the initial crack opening, is the crack opening at 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), is the displacement of upper seam surface at after loading, is the air depth within the seam with unsaturated water at the initial moment, and is the imposed load.
Since the track slab deformed slightly within the elastic range with the effect of train load, the displacement of the upper seam surface can be denoted aswhere is the deformation coefficient associated with seam depth and the acting position of load.
It can be seen from Figure 11 that the initial volume of the air domain within the seam per unit length at the initial moment is
Owing to the balance between air and water within the seam and the absolute pressure at the outlet of the seam, the initial pressure of air and water within the seam is equal to the absolute pressure at the outlet of the seam.
The interface between air and water shifts from to at time t. According to the incompressibility of water, volumes 1 and 2 are equal; thus,
The volume variation of the air domain at this time is
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,
It could be derived from (74), (76), and (77) that the pressure in the action of train load at time t is
The pressure at the air-water interface is equal to that of the air domain at this point; thus,
At stationary point ,
Since the pressure at the air-water interface within the seam is equal to that of the air domain ,
The pressure at the seam outlet is equal to the absolute pressure of liquid, namely,
5.3.4. Analytical Expression of Pressure Distribution within Seam
The air domain pressure in area is
The pressure distribution expression of water in area 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):
Similarly, the pressure distribution expression of water in area within the seam could be obtained based on (82):
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; is the vertical displacement of point A in the numerical calculation model, and 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 Figure 10 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.
As shown in Table 5, 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 Figures 12 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 Figures 13 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 Table 5, 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.
Damage amounts Dvibration + 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].
In addition, the formula was calculated based on Palmgren–Miner cumulative fatigue-damage theory, where N denotes damage fatigue life (in years); is the critical value of cumulative concrete damage, which is set to 0.9500000 based on current structure damage diagnosis and maintenance and reinforcement; is the initial structure damage amount; and is the damage amount at point A, as shown in Figure 16.
The fatigue life (rounded to the closest integer; see Table 6) 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.
It 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.Data Availability
The data used to support the findings of this study are included within the article.
Additional Points
Highlights. (1) Dynamic response analysis of high-speed railway tunnel base with defects. (2) Fatigue-life analysis of high-speed railway tunnel base with defects. (3) Antifatigue life was 56 years in the presence of groundwater in defect, and antifatigue life was 62 years without groundwater in defect.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was financially supported by the Science and Technology Development Plan Project of China Academy of Railway Sciences Group Co., Ltd (2019YJ027 and 2020YJ170).