*4.1. Fundamental Theory of the Dynamic Contact Force Method*

According to the basic theory of kinematics, the dynamic balance equation for the contact nodes can be obtained after the tunnel system is discretized using finite elements:

$$\mathbf{M}\ddot{\mathbf{U}} + \mathbf{F}\_{\text{damp}} + \mathbf{F}\_{\text{int}} = \mathbf{F} + \mathbf{R} \tag{13}$$

where *M* is the mass matrix, *U* is the displacement vector, *F*damp is the damping force vector, *F*int is the internal force vector, which is obtained by the stress integration of the elements, *F* is the load vector of the seismic wave, and *R* is the contact force vector, *R* = *N* + *T*.

When the movement states of the contact nodes are known at time *t*, the integration schemes of the movement variables at time *t* + Δ*t* can be obtained, using the explicit central difference method to solve Equation (13):

$$\mathcal{U}^{t+\Delta t} = \mathcal{U}^{t+\Delta t} + \Delta \mathcal{U}^{t+\Delta t} \tag{14}$$

$$
\dot{\boldsymbol{\mathcal{U}}}^{t+\Delta t} = \frac{2(\boldsymbol{\mathcal{U}}^{t+\Delta t} - \boldsymbol{\mathcal{U}}^t)}{\Delta t} - \dot{\boldsymbol{\mathcal{U}}}^t \tag{15}
$$

$$\overset{\text{L}}{\mathcal{U}}^{t+\Delta t} = \mathcal{U}^t + \Delta t \dot{\mathcal{U}}^t + \frac{\Delta t^2}{2} \mathcal{M}^{-1} (\mathcal{F}^t - \mathcal{F}\_{\text{int}}^t - \mathcal{F}\_{\text{damp}}^t) \tag{16}$$

$$
\Delta \mathbf{U}^{t+\Delta t} = \frac{\Delta t^2}{2} \mathbf{M}^{-1} \mathbf{R}^t \tag{17}
$$

*t*+Δ*t*

where ¯ *U* is the predicted displacement when ignoring the contact force *R<sup>t</sup>* and Δ*Ut*+Δ*<sup>t</sup>* is the modified displacement caused by *R<sup>t</sup>* .

From the above integration schemes, we can see that the key to determining the movement states of the contact nodes at time *t* + Δ*t* is solving for their contact force *R<sup>t</sup>* . *R<sup>t</sup>* should be determined based on the contact states and the contact conditions between the rock mass and the fault within time *t*–*t* + Δ*t*.

#### *4.2. Solving for the Contact Force and Judging of the Contact State*

The point-to-point contact is the most widely used contact type [14] and is very effective in simulating the small sliding problem, but it is difficult to simulate the large sliding problem well using this method. Based on the point-to-point contact type, the calculation formula for the contact force was derived, and methods of judging the contact state were designed. In order to simulate the large sliding problem between the rock mass and the fault, we needed to extend the point-to-point contact method to include the point-to-surface contact.

#### (1) Point-to-Point Contact Type

When the contact point pair has no large relative sliding, it can be considered that the contact points are of the point-to-point type (as shown in Figure 2). In order to solve the contact force *R<sup>t</sup>* for this contact type, it can be assumed that the contact points are in the bond contact state. Then, they should satisfy the contact displacement conditions (Equations (5) and (6)). The normal and tangential contact gaps of the contact point pair are denoted Δ<sup>n</sup> and Δt, and hence

$$
\Delta\_{\mathbf{n}} = \mathbf{n}\_1^T (\mathbf{U}\_{l'} \quad -\mathbf{U}\_l \quad \text{)}, \\
\Delta\_{\mathbf{t}} = \mathbf{t}\_1^T [(\mathbf{U}\_{l'} \quad -\mathbf{U}\_l \quad \mathbf{t}) - (\mathbf{U}\_{l'}^t \quad \mathbf{t}) - (\mathbf{U}\_{l'}^t - \mathbf{U}\_l^t)] \tag{18}
$$

**Figure 2.** Point-to-point contact.

Combining Equations (5), (6), (14), (17), and (18), and according to the contact force conditions (Equation (7)), the expression for the contact force *R<sup>t</sup> <sup>l</sup>* can be derived in the form of a normal component *N<sup>t</sup> <sup>l</sup>* and a tangential component *<sup>T</sup><sup>t</sup> l* :

$$\Delta \mathbf{N}\_l^t = \frac{2M\_l M\_{l'}}{(M\_l + M\_{l'})\Delta t^2} \Delta\_{\mathbf{n}} \mathbf{u}\_{1\prime} \ T\_l^t = \frac{2M\_l M\_{l'}}{(M\_l + M\_{l'})\Delta t^2} \Delta\_{\mathbf{t}} \mathbf{t}\_1 \tag{19}$$

where *Ml* and *Ml* are the lumped masses of *l* and *l*, respectively.

If Δ<sup>n</sup> > 0, this indicates that the contact points are in a tension state. If *N<sup>t</sup> <sup>l</sup>* satisfies Equation (8), this assumption is correct; otherwise, the contact points enter a separation state, and then the contact force is modified according to Equation (10). If Δ<sup>n</sup> ≤ 0, this indicates that the contact points are in a compression state. If *T<sup>t</sup> <sup>l</sup>* satisfies Equation (9), this assumption is correct; otherwise, the contact points enter a sliding contact state, and then the contact force is modified according to Equation (11).

#### (2) Point-to-Surface Contact Type

When the contact point pair has large relative sliding, it can be considered that the contact points are of the point-to-surface type (as shown in Figure 3), and the cohesive force is no longer considered in this case. In order to solve for the contact force *R<sup>t</sup>* , it can be assumed that the contact points are in a static contact state, so they should satisfy the contact displacement conditions (Equations (5) and (6)). Assuming that *l* is the projection point of *l* on the contact interface, the displacement of *l* is given by

$$\mathbf{U}\_{l'}^{t} = \sum\_{p} \Phi\_{p} \mathbf{U}\_{p}^{t} \tag{20}$$

where *p* is the node number of the element to which the contact interface belongs, Φ*<sup>p</sup>* is the shape-function value of *p* and *Ml* is the equivalent lumped mass of *l* , as shown below [24].

$$M\_{l'} = \sum\_{p} m\_{p'} \ m\_p = \frac{M\_p \Phi\_p}{\sum\_{q} \Phi\_q^2} \tag{21}$$

where *mp* is the mass contribution of *p*, *Mp* is the lumped mass of *p*, and *q* has the same meaning as *p*.

**Figure 3.** Point-to-surface contact.

According to the contact conditions, the expression for the contact force *R<sup>t</sup> <sup>l</sup>* can also be derived in the form of Equation (19). The contact force of *p* is calculated via

$$\mathbf{N}\_p^t = \frac{m\_p}{M\_{l'}} \mathbf{N}\_{l'}^t \; \; T\_p^t = \frac{m\_p}{M\_{l'}} T\_{l'}^t \tag{22}$$

*t*+Δ*t*

If Δ<sup>n</sup> > 0, this indicates that the contact points enter a separation state, and then the contact force is modified according to Equation (10). If Δ<sup>n</sup> ≤ 0, this indicates that the contact points are in a compression state. If *T<sup>t</sup> <sup>l</sup>* satisfies Equation (12), this assumption is correct; otherwise, the contact points enter a sliding contact state, and then the contact force is modified according to Equation (11).

#### *4.3. Calculation Flow of the Improved Dynamic Contact Force Method*

Before the seismic action is applied, the discontinuity of the finite element model can be realized by adding common nodes on the contact interface between the rock mass and the fault. The main calculation procedures for the improved dynamic contact force

method are: (1) calculate the predicted displacement ¯ *U* , ignoring the contact force *R<sup>t</sup>* ; (2) determine the contact types of all the contact point pairs; (3) use the improved method to calculate and modify *R<sup>t</sup>* ; (4) calculate the modified displacement Δ*Ut*+Δ*<sup>t</sup>* caused by *R<sup>t</sup>* ; (5) calculate the vibration deterioration coefficient *Dt*+Δ*<sup>t</sup>* . The dynamic finite element calculation flow for the tunnel contact system at time *t* + Δ*t* is shown in Figure 4.

**Figure 4.** Calculation flow for the dynamic contact system of tunnel at time *t* + Δ*t*.
