Dynamic Response Simulation of Lining Structure for Tunnel Portal Section under Seismic Load

Portal section is the weak link of seismic fortification for tunnel structure. Assuming that seismic wave is the vertical incident elastic plane wave, the plane wave input method for the portal section was discussed in this paper; that is, the wave input problem can be converted to the problem of calculating equivalent nodal force at artificial boundaries. Based on different damage evolution processes of concrete under tension and compression conditions, the tension and compression damage variables were defined and solved, respectively. And then a simple elastic dynamic damaged constitutive model for concrete lining was built. According to the characteristics of dynamic interaction between the lining and rock, and based on the dynamic contact force algorithm, an analytical model for joint loading between the lining and rock was built. This model can simulate lining features such as bond, separation, and slip under seismic load.The dynamic response characteristics of lining structure for the portal section under seismic load were analyzed by taking example for an exit section of Dianzhong diversion project in strong earthquake area. The results show that the relative displacement magnitudes of the lining parts are related to the vibration direction of the seismic wave, and the peak displacements decrease gradually to the fixed values from the portal to the interior. The damage coefficients of the lining parts accumulate gradually over time, and the farther the lining is away from the portal, the less serious the seismic damage is. The separation and slip zone distributions of the lining are basically consistent with its severe seismic damage area, which are mainly at haunch, spandrel, and arch foot within a certain range of distance from the portal. The seismic fortified length and key fortified parts of tunnel structure for the portal section can be determined according to the calculation results.


Introduction
It was found that the portal section is destroyed easily when an earthquake occurs according to the seismic damage investigation of highway tunnels in epicentral area after Wenchuan earthquake.The portal section is the seismic weak section of tunnel structure, second only to the fault fracture zone [1][2][3][4][5].The buried depth of the tunnel portal section is usually shallow, and the geological conditions are usually poor.The rock formations are mostly seriously weathered deposits, which can easily cause slope instability.When a strong earthquake happens in the portal section, ground surface crack, rockfall, and slope collapse occur frequently, which can easily cause lining crack, dislocation, and collapse, thus affecting the normal operation of tunnel.Therefore, study of seismic stability of the portal section is one of the urgent problems that need to be solved in tunnel.
There have been some research results for the seismic response analysis of the portal section.The main research methods are model test and numerical simulation.In the model test, Shen et al. [6] carried out a shaking table model test of dynamic response of the portal section and analyzed the failure modes of the lining by taking Ya'an-Lugu expressway mountain tunnel as the support project.Jiang et al. [7] studied the acceleration response of the soil and lining and the strain response of the lining in the portal section of Galongla tunnel through a large-scale shaking table model test.Li et al. [8] studied the acceleration response of the rock and the internal force distribution characteristics of the lining in the portal section of shallow buried biased and unbiased tunnel through a shaking table model test.Nevertheless, the application and development of the model test are subject to many restrictions, considering that it is difficult to simulate the real engineering environment and artificial boundary conditions and that it spends a lot of manpower, materials, and financial resources.In contrast, the numerical simulation is favored by researchers for its economy and convenience.In this aspect, Jiang et al. [9] analyzed the acceleration, displacement, and stress response characteristics of the lining in the portal section of Galongla tunnel using the dynamic finite element.Y. S. Li and T. B. Li [10] evaluated the seismic safety of the portal section from the portal slope, portal building, and lining structure based on FLAC3D by taking as an example a tunnel of national highway 318.Bai et al. [11] adopted the nonlinear dynamic analysis method to simulate the seismic response law of the lining and slope in the portal section by combining with a nuclear power intake tunnel project.
The key to seismic stability analysis of the portal section is simulating the dynamic response of the concrete lining, so establishing a reasonable dynamic constitutive model for the concrete and a joint loading model between the lining and rock is particularly important.However, most of the present studies assume that the lining material is linear elastic and that the lining and rock share nodes without considering the dynamic contact behavior between the two.These assumptions are obviously inconsistent with the actual situation.
This paper assumed that the seismic wave is the elastic plane wave of vertical incidence.The input methods of plane P wave and S wave for the portal section were discussed.A simple dynamic damage model for the concrete lining was built according to different damage evolution processes of the concrete under tension and compression conditions.An analytical model for joint loading between the lining and rock was built based on the dynamic contact force algorithm.The seismic response characteristics of the lining structure were simulated combining with a case study for an exit section of Dianzhong diversion project, which is expected to provide a useful reference for seismic design of the tunnel portal section.

Input Method of Plane Seismic Wave for Portal Section
2.1.Basic Ideas.The total wave field can be decomposed into internal and external wave field based on the wave field decomposition principle.There is no need to deal with the external wave field specifically, because it can be transmitted at artificial boundaries and has no influence on the finite element calculation.Therefore, the key to realizing the accurate input of seismic load is to obtain the internal wave field at artificial boundaries.The internal wave field at the bottom boundary of the model is equal to the incident wave field considering the vertical incidence of the seismic wave.Because the incident field of the lateral boundary and the wave field reflected by ground surface are both parallel to the boundary, whereas the local artificial boundary cannot simulate the wave field parallel to the lateral boundary, the internal field of the lateral boundary can be calculated as the free field at this point.When the incident wave field is known, the free field of the lateral boundary can be obtained by superimposing the incident wave field and reflected wave field.The wave input problem can be transformed into the problem of calculating equivalent nodal force acting on artificial boundaries based on the viscoelastic artificial boundary [12], so as to realize the input of the seismic wave.
Assuming that the seismic wave is the elastic plane wave of vertical incidence, the wave front is parallel to the bottom of the tunnel model, as shown in Figure 1.The following provisions are given in this paper considering different vibration directions of the seismic wave: P wave vibrates along the vertical direction ( direction), SH wave vibrates along the transverse direction ( direction), and SV wave vibrates along the tunnel axis direction ( direction).The calculation methods of seismic load for plane P wave and S wave under vertical incidence condition are stated as follows, respectively.

Plane P Wave Incidence.
Assuming that  P () is the displacement time-history of incident P wave, Δ 1 , Δ 2 , and Δ 3 are the delay time of a node  relative to zero time when P wave propagates to the left, right, and front boundaries, and Δ 4 , Δ 5 , and Δ 6 are the corresponding delay time of reflected P wave by ground surface, respectively, the displacement field of boundary nodes of the model is where where  P is the wave velocity of P wave. is the vertical coordinate of the node .
Given the displacement field, the corresponding velocity field can be solved by derivation or difference.The corresponding stress field can be calculated by the internal wave field at artificial boundaries according to the generalized Hook's law, and then the equivalent nodal force at artificial boundaries can be obtained [13,14].
For the bottom boundary: For the left boundary: For the right boundary: For the front and back boundaries: where  and  are the density and lame constant of the medium, respectively.  is the control area of the node .
The superscript of the equivalent nodal force denotes the artificial boundary surface where the node  locates, and the subscript denotes the component direction of the nodal force. T and  T are the tangential parameters of physical components at artificial boundaries, and  N and  N are the normal parameters, which are calculated by where  is the shear modulus of the medium. S is the wave velocity of S wave. is the distance between the geometric center of the tunnel model and current artificial boundary surface.

Plane S Wave
where The equivalent nodal force at artificial boundaries is as follows.
For the bottom boundary: For the left boundary: For the right boundary: For the front and back boundaries:

Dynamic Damage Model for Concrete Lining
Concrete is a kind of heterogeneous man-made stone.Its material composition determines natural micro cracks existing in the concrete.It is easy to produce damage characteristics under seismic load.Therefore, establishing a reasonable dynamic damage model for the concrete is the key to simulating the seismic response of the lining structure.

Dynamic Damaged Constitutive Model for Concrete.
Assuming that the concrete is an isotropic medium, and based on the basic theory of continuum damage mechanics, the stress-strain relationship of the concrete can be expressed as where  is the stress,  is the strain, and  is the damage variable. 0 is the initial dynamic elastic modulus of the concrete.
The concrete shows different strength and stiffness characteristics under tension and compression conditions, which is called unilateral effect [15].To well describe the influence of unilateral effect, the damage variables  t and  c are introduced to reflect the weakness of the concrete material caused by tension and compression damage, so the stressstrain relationships of the concrete under uniaxial tension and compression conditions are obtained:

Dynamic Damage Evolution Model for Concrete.
It is necessary to define and solve the damage variables  t and  c in order to perfect the dynamic damaged constitutive model for the concrete.The dynamic damage model for the concrete in the literature [16] is introduced, as shown in Figure 2.This model assumes that the concrete is in linear elastic state in the initial stage, and the stress and strain vary in proportion.
When the strain exceeds the critical strain, the concrete shows damage.The greater the strain is, the more serious the damage is, until the concrete is completely destroyed.The evolution processes of tension and compression damage of the concrete show some differences in this model, which are explained below.The tensile stress (strain) is set as positive and the compressive stress (strain) is set as negative in this paper.
The damage evolution process of the concrete can be divided into three stages under uniaxial tension condition: elastic stage, damage softening stage, and complete damage stage.Correspondingly, the damage variable  t is expressed as where  t0 is the elastic ultimate tensile strain corresponding to the tensile strength  t0 . tu is the ultimate tensile strain, and  tu =  t0 . is the ultimate tensile strain coefficient and is generally 2-5.

Compression Tension
The damage evolution process of the concrete can be also divided into three stages under uniaxial compression condition: elastic stage, damage softening stage, and residual strength stage.Correspondingly, the damage variable  c is expressed as where  c0 is the elastic ultimate compressive strain corresponding to the compressive strength  c0 . cr is the compressive strain corresponding to the residual strength  cr , and  cr =  c0 ,  cr =  c0 . is the residual compressive strain coefficient, and  is the residual strength coefficient.

Three-Dimensional Dynamic Damage Model for Concrete.
Assuming that the concrete damage is also isotropic in threedimensional stress state, when an element meets the maximum tensile strain criterion and produces tension damage, the one-dimensional damage model can be extended to threedimensional state.At this point,  t in ( 16) can be replaced by the equivalent strain . is defined as where  1 ,  2 , and  3 are the principal strains of the element in three directions, arranged in order of magnitudes.⟨ ⟩ is a functional expression and is defined as Similarly, when an element meets the Mohr-Coulomb yield criterion and produces compression damage, the threedimensional dynamic compressive damage model for the concrete can be obtained by replacing  c in (17) with the maximum principal compressive strain  3 . 3 is calculated by the following equation [17]: where  is the internal friction angle of the concrete. is the Poisson ratio. 1 ,  2 , and  3 are the principal stresses of the element in three directions.

Dynamic Strength Parameters Values of Concrete.
The concrete material generally shows obvious rate correlation property under dynamic load, and the mechanical parameters increase with the increase of the strain rate.To describe this growth property, DIF is adopted to express the dynamic increase factor of the concrete parameters.Research shows that the elastic modulus changing law of the concrete is relatively discrete as the strain rate.Therefore, the dynamic elastic modulus can be replaced by the static elastic modulus from the bias safety perspective of engineering project.However, the concrete strength shows some rules as the strain rate changes.The European International Concrete Committee gave the calculation formulas of the dynamic tensile strength increase factor DIFt and compressive strength increase factor DIFc of the concrete under different strain rates based on a large number of experimental data [18][19][20]: where  t0 and  c0 are the dynamic uniaxial tensile and compressive strengths, respectively. t and  c are the static uniaxial tensile and compressive strengths, respectively.ε d is the dynamic strain rate and is generally 3 × 10 −6 s −1 -300 s −1 .ε ts and ε cs are the static strain rates under tension and compression conditions, whose magnitudes are 3 × 10 −6 s −1 and 30 × 10 −6 s −1 , respectively. = (10 + 3  c /5) −1 ,   c is the static compressive strength of a cylinder, and log  = 7.11 − 2.33. = (5 + 3 cu /4) −1 ,  cu is the static compressive strength of a cube, and log  = 6.156 − 0.49.

Analytical Model for Joint Loading between Lining and Rock
There are complex dynamic contact behaviors in the joint seismic response process between the lining and rock.The earthquake may lead to lining damage, crack, separation, or slip.Engineering practice shows that the separation or slip failure of the lining structure is mainly partial in the tunnel without the fault fracture zone, and the large slip failure is rare.Therefore, it is assumed that the contact interface nodes  between the lining and rock are always or approximately in point to point contact state during the earthquake and that the contact node pairs are in bond contact state before the earthquake [21].
For a contact node pair   and , as shown in Figure 3, if the displacements of this node pair are u    and u   at  time, and the dynamic contact forces are R    and R   (R    + R   = 0), the displacements u +Δ   and u +Δ  of this node pair at  + Δ time can be obtained by solving the following equation using the displayed center difference method: where M, C, and K are the mass, damping, and stiffness matrices of the contact nodes, respectively.ü , u , and u are the acceleration, velocity, and displacement vectors of the contact nodes, respectively.F is the external load vector.R is the dynamic contact force vector.Thus, the relative normal displacement Δ 1 and tangential displacement Δ 2 of the contact node pair at  + Δ time can be expressed as where n  is the unit normal vector of the contact node pair, which points to the node  from the node   .Assuming that the contact node pair is still in bond contact state at  + Δ time, the dynamic contact force R +Δ  of the node  at this time is where N +Δ  and T +Δ  are the normal and tangential components of R +Δ  , respectively.  is the control area of the contact node pair. n and  t are the normal and tangential stiffness of the contact interface, respectively, which are calculated by the following equation learning from FLAC3D [22]: where  and  are the bulk modulus and shear modulus of the medium with the maximum stiffness on both sides of the contact interface, respectively.Δ min is the minimum size of the connecting area at normal direction on both sides of the contact interface.
After each time step is calculated, the contact state of the contact interface needs to be checked, and the dynamic contact force needs to be corrected.The detailed method is as follows.
(1) If ‖Δ 1 ‖ < 0, it shows that the contact node pair has a tendency to be separated from each other.And if ‖N +Δ  ‖ >   , the contact node pair is in separation state; here, where  s and  d are the static and dynamic friction coefficients of the contact interface, respectively. is the cohesion force of the contact interface.If the contact node pair is in separation or slip state at  time, the cohesion force will not be considered at  + Δ time and the later time; that is,  = 0.

Calculation Model and Parameters. Dianzhong diversion
project is mainly composed of water source project, water transmission project, auxiliary project, and so on.It is a water conservancy project with comprehensive utilization of water resources.The total length of the water transmission trunk line is 848.18 km.The project scale is huge, and the geological conditions along the route are complex.The seismic fortified intensity is high, and the tunnel stability problem is salient.An exit section along the route is selected to analyze its seismic response in this paper.The maximum buried depth of the section is 35 m, and the slope angle is about 50 ∘ .The design depth of water is 7.51 m, and the surrounding rock type is mainly IV.The excavation section is in horseshoe shape with a maximum excavation size of 8.3 m × 9.87 m.C25 reinforced concrete lining structure is adopted, and the lining thickness is 0.5 m.
A three-dimensional finite element model for the exit section is built.This model is discretized totally by 8-node   hexahedron elements, and 110,532 elements and 117,280 nodes are divided, 5,760 of which are lining elements.-axis of the model is perpendicular to the tunnel axis along horizontal direction.-axis coincides with the tunnel axis, with positive direction going along the water flow.-axis coincides with the geodetic coordinate.The detailed sizes of the finite element model are shown in Figure 4.
According to the geostress test results, the lateral pressure coefficients are valued as   = 1.0,   = 0.7, and   = 1.0.The mechanical parameters of the lining and rock are provided in Table 1.

Calculation Conditions.
The calculation program adopts the three-dimensional elastoplastic dynamic damaged displayed finite element platform [23,24], in which the seismic wave input method, dynamic damage model for the concrete, and joint loading model between the lining and rock in this paper are embedded.
The viscous-elastic artificial boundary is applied at the bottom of the model, and the free field artificial boundary is applied at the four sides.American strong earthquake record, El-Centro wave, is selected as the seismic wave, whose peak acceleration is 3.417 m/s 2 .A section of 20 s is intercepted as the incident wave, as shown in Figure 5, which is processed by filtering and baseline correction.The seismic wave is imported vertically from the bottom of the model.The actions of P wave and SH wave on the tunnel are taken into account.SH wave takes the incident wave as shown in Figure 5, and P wave is 2/3 that of SH wave.
A monitoring section is arranged every 6 m from the portal along the tunnel axis, and a total of 12 monitoring sections are arranged.Five monitoring points are arranged on the lining structure in each monitoring section, as shown in Figure 6, which are used to monitor the characteristics of stress, displacement, and so on of the lining during the earthquake.

Displacement Analysis of Lining.
Different parts of the lining are in synchronous vibration state during the earthquake.The maximum positive displacement is defined as the peak displacement; then the variation laws of and peak displacements along the tunnel axis of the 5 monitoring points in the 12 monitoring sections can be obtained, as plotted in Figures 7 and 8  As seen in Figure 7, -peak displacements of the lining parts decrease gradually along the tunnel axis.The peak displacements decrease greatly within 0 m-48 m distance from the portal.When the distance is larger than 48 m, the peak displacements vary a little.This is mainly because the portal slope has reflection action on the seismic wave.Meanwhile, the reflected wave produces waveform conversion phenomenon, which leads to the complex wave field in the 0 m-48 m section; thus the seismic response of the lining is affected [25].As seen in Figure 8, -peak displacements variation of the lining parts show a similar law as that of  directional.
The peak displacements of the spandrel and haunch of the lining are larger than those of the other parts, and the peak displacement of the bottom arch is the smallest under  SH wave action.The maximum relative displacement reaches 0.98 cm.The peak displacement of the bottom arch of the lining is larger than that of the other parts, and the peak displacement of the top arch is the smallest under P wave action.The maximum relative displacement is 0.55 cm.These indicate that the relative displacement magnitudes of the lining parts are related to the vibration direction of the seismic wave.

Stress Analysis of Lining.
The seismic damage failure of the lining is mainly caused by tensile damage, because the compressive strength of the concrete is much greater than its tensile strength.For that reason, this paper mainly analyzes the maximum principal stress variation law of the lining.Taking the sixth monitoring section as an example, the maximum principal stress time-history curves of monitoring points are plotted in Figure 9.
As seen in Figure 9, the maximum principal stresses of the lining parts are not large at 0 s-2 s.They increase rapidly to the tensile strength of the concrete after 2 s, which indicates that the lining shows tension damage.The maximum principal stresses fluctuate violently at 2 s-7 s, which may aggravate the lining damage and cause the lining fatigue failure.The maximum principal stresses after the earthquake are not much different from that at the initial time.

Damage Analysis of Lining.
Based on the dynamic damage model for the concrete in this paper, the time-history curves of the damage coefficients of monitoring points are obtained, as plotted in Figure 10 (taking the sixth monitoring section as an example).
As seen in Figure 10, there is no damage appearing on the lining parts at 0 s-2 s.The damage coefficients increase gradually over time after 2 s, and they increase obviously at 2 s-7 s when the seismic wave vibrates violently.The damage coefficients increase a little after 7 s.It can be seen that the damage evolution process of the lining is closely related to the changing law of the maximum principal stresses.Overall, the damage coefficient of the haunch is the largest, the next are the spandrel and arch foot, and the smallest are the top arch and bottom arch.
The damage coefficients' distribution of the lining structure (taking 60 m range of distance from the portal) after the earthquake is shown in Figure 11.
As seen in Figure 11, the damage degradation degrees of different parts of the lining are different.Looking horizontally, the haunch, spandrel, and arch foot are damaged seriously, while the top arch and bottom arch are damaged lightly.Looking vertically, the portal is damaged most seriously, and the farther the lining is away from the portal, the smaller the damage is.If the area where the damage coefficients are greater than 0.2 is defined as the seismic damage area, the seismic damage area is mainly distributed within 48 m range of distance from the portal.The distribution range of the severe seismic damage area ( > 0.5) is about 21 m. the dynamic contact force of the contact nodes between the lining and rock breaks through the cohesive force during the earthquake, the contact interface will separate or slip easily.And the process is irreversible.The separation and slip zone distributions of the lining structure (taking 60 m range of distance from the portal) after the earthquake are shown in Figures 12 and 13.

Separation and Slip Analysis of
As seen in Figure 12, the separation zone of the lining is mainly distributed at the haunch, top arch, and bottom arch.The separation zone of the haunch extends long from the portal to the interior of the tunnel, reaching 21 m.The separation zone distribution ranges of the top arch and bottom arch are smaller, and the extension length is limited to 9 m.As seen in Figure 13, the slip zone of the lining is mainly distributed at the spandrel and arch foot, on which the extension length of the slip zone reaches 21 m.It can be seen that the separation and slip zones of the lining are mainly distributed at the haunch, spandrel, and arch foot within a certain range of distance from the portal, which is basically consistent with the distribution law of the severe seismic damage area of the lining.

Determination of Seismic Fortified Length for Tunnel.
The displacement response of the lining is large within 48 m range of distance from the portal, where the lining damage is relatively serious according to its displacement and damage analysis.Therefore, the seismic fortified length of tunnel structure for the portal section can be set to 48 m.In addition, the separation and slip zones of the lining are mainly distributed within 21 m range of distance from the portal, which can be served as the key fortified length.

Conclusions
Based on the dynamic algorithm of the concrete lining, the dynamic response characteristics of the lining structure under seismic load were analyzed by taking as an example an exit section of Dianzhong diversion project.The following conclusions are drawn.
(1) The relative displacement magnitudes of the lining parts are related to the vibration direction of the seismic wave under seismic load.The farther the lining is away from the portal, the smaller the peak displacements are, and they tend gradually to the fixed values.
(2) The maximum principal stresses of the lining parts can easily reach the tensile strength.The violent fluctuation of the maximum principal stresses may aggravate the lining damage and cause the lining fatigue failure.
(3) The damage coefficients of the lining parts accumulate gradually over time.The seismic damage area of the lining is mainly distributed at the haunch, spandrel, and arch foot within a certain range of distance from the portal.The farther the lining is away from the portal, the less serious the seismic damage is.
(4) The separation and slip zones of the lining are mainly distributed at the haunch, spandrel, and arch foot within a certain range of distance from the portal, which are basically consistent with the distribution law of the severe seismic damage area of the lining.
(5) The seismic fortified length and key fortified parts of tunnel structure for the portal section can be determined according to the displacement, damage, separation, and slip analysis of the lining.

Figure 1 :
Figure 1: Input model for seismic wave.

Figure 2 :
Figure 2: Dynamic damage model for concrete.

Figure 3 :
Figure 3: Point to point contact.

Figure 4 :
Figure 4: Calculation model for portal section.

Figure 9 :
Figure 9: Maximum principal stress time-history of monitoring points.

Figure 10 :Figure 11 :
Figure 10: Damage coefficients variation time-history of monitoring points.

Figure 12 :Figure 13 :
Figure 12: Separation zone distribution of lining after the earthquake.
Incidence.S wave is divided into SH wave and SV wave.The calculation methods of seismic load of the two under vertical incidence condition are similar.In this paper, SH wave is taken as an example for illustration.Assuming that  S () is the displacement time-history of incident SH wave, Δ 7 , Δ 8 , and Δ 9 are the delay time of the node  relative to zero time when P wave propagates to the left, right, and front boundaries, and Δ 10 , Δ 11 , and Δ 12 are the corresponding delay time of reflected SH wave by ground surface, respectively, the displacement field of boundary nodes of the model is shows that the contact node pair has a tendency to be embedded in each other.And if ‖T +Δ ,

Table 1 :
Mechanical parameters of materials.