Performance of the Generalized- α (G- α ) Algorithm for Discontinuous Dynamics by the Numerical Manifold Method

: Performance of the generalized- α (G- α ) time integration scheme using the numerical manifold method was examined for the continuous and discontinuous rock dynamic problems in the present paper. Inﬂuence of the generalized- α time integration scheme on the numerical stability and accuracy were studied using the dynamic equation of motion, in which various parameters such as control parameters and kinetic damping were compared by the ampliﬁcation matrix A. Furthermore, convergence productiveness of the open–close iterations for the generalized- α time integration with the case of contact was also investigated. The validations of the generalized- α algorithm were conducted by a transient response analysis of an elastic strip subjected to harmonic function loading, contact analysis of layered rock deﬂection, and mechanical behavior of dynamical tunneling of jointed rock masses, respectively. It can be concluded from the numerical results that the proposed generalized- α scheme used in the present study has superior qualities compared to the original algorithm in solving rock dynamic problems involving more nonlinear contacts using the numerical manifold method.


Introduction
For rock dynamic simulations, a step-by-step time integration scheme can be employed to obtain a numerical solution of the problems. Furthermore, to achieve expected numerical accuracy, an integration of unconditional stability and algorithmic damping is required and widely recognized [1]. Normally, the commonly used time integration schemes in numerical simulations include the family of Newmark algorithms [2], the Wilson-θ method [3], the Hilber-Hughes-Taylor method by Hilber et al. [4], the Wood-Bossak-Zienkiewicz method of Wood et al. [5], the Houbolt method [6], the ρ method [7], the θ 1 method [8], and the generalized-α (G-α) method by Chung and Hulbert [9]. These methods both have unconditional stability and possess numerical dissipation properties of high-frequency dynamics along with second-order accuracy. However, the pathological overshooting property is found in the Wilson-θ method [1], thus the Houbolt and HHT-α have been preferred and widely used in continuum-based finite element analysis. As discontinuous-based methods such as the discontinuous deformation analysis (DDA) use the Newmark method, which allows for a larger time step and have unconditional stability and are dissipative to consider the penalty formulation of the contact analogies [10,11]. Nevertheless, the calculation of complex nonlinear contacts is often very time-consuming

The Foundational Theory of the NMM
Compared with the existing numerical methods such as the finite element method and discrete element method (DEM), the most significant feature of NMM is the use of a dual coverage system, and the details of the finite coverage system and coverage segmentation technology can be found in [31,32]. In the current section, the NMM theory is briefly illustrated through the division of unified (PU) functions and the OCI criteria for structural dynamic analysis.

Dual Cover System
The NMM uses a dual cover system: the mathematical cover (MC) and physical cover (PC). Figure 1 illustrates an example of a physical domain with an internal discontinuity overlapped by two MCs denoted by the triangle M 1 and rectangle M 2 , respectively. The MCs are independently constructed and do not need to conform to the domain boundaries and the internal discontinuities (e.g., bedding, crack, joint, interlayer, etc.). The PCs are a sub-set of the MCs obtained from intersections with the physical domain. A manifold element (ME) is generated as the common area of the overlapped PCs. As plotted in Figure 1a, each patch (i.e., triangle and rectangle) is termed as a MC, denoted by M i (i = 1, 2); external boundary and internal joints or cracks may split one MC into several separate sub-patches, each one within the material domain Ω is regarded as a PC, which is denoted by P j i (I = 1, 2; j = 1, 2). Figure 1b shows that the material domain Ω is intersected by patch M 1 to generate one PC, denoted by P 1 1 . When internal discontinuities (i.e., cracks or joints) exist, each discontinuity is considered as one special physical domain to form new PCs. When a crack passes through the whole patch within the material domain, two isolated PCs are formed by the crack in M 2 , and denoted by P 1 2 and P 2 2 , respectively. When the crack cuts MC partially, only one PC is produced within the material domain such as P 1 1 . The common area of several overlapping PCs is termed as a ME. As can be seen in Figure 1c, five MEs were formed by the overlapping P 1 1 , P 1 2 , and P 2 2 , which can be represented by E 1 P 1 1 , E 2 P 1 1 ∩ P 1 2 , E 3 P 1 1 ∩ P 2 2 , E 4 P 1 2 , and E 5 P 2 2 , respectively. that algorithms such as G-α methods implemented in NMM codes can well solve the structural dynamic problem compared to known Newmark methods involving nonlinear contacts.

The Foundational Theory of the NMM
Compared with the existing numerical methods such as the finite element method and discrete element method (DEM), the most significant feature of NMM is the use of a dual coverage system, and the details of the finite coverage system and coverage segmentation technology can be found in [31,32]. In the current section, the NMM theory is briefly illustrated through the division of unified (PU) functions and the OCI criteria for structural dynamic analysis.

Dual Cover System
The NMM uses a dual cover system: the mathematical cover (MC) and physical cover (PC). Figure 1 illustrates an example of a physical domain with an internal discontinuity overlapped by two MCs denoted by the triangle 1 and rectangle 2 , respectively. The MCs are independently constructed and do not need to conform to the domain boundaries and the internal discontinuities (e.g., bedding, crack, joint, interlayer, etc.). The PCs are a sub-set of the MCs obtained from intersections with the physical domain. A manifold element (ME) is generated as the common area of the overlapped PCs. As plotted in Figure  1a, each patch (i.e., triangle and rectangle) is termed as a MC, denoted by (i = 1, 2); external boundary and internal joints or cracks may split one MC into several separate sub-patches, each one within the material domain is regarded as a PC, which is denoted by (I = 1, 2; j = 1, 2). Figure 1b shows that the material domain is intersected by patch 1 to generate one PC, denoted by 1 1 . When internal discontinuities (i.e., cracks or joints) exist, each discontinuity is considered as one special physical domain to form Since the NMM does not require MC sides to coincide with the material boundaries and discontinuities, arbitrary shape MCs can be employed. For convenience, a regularly structured mathematical mesh was employed in the NMM. As seen in Figure 2a, a regular triangular mesh was constructed, in which each MC was formed through six neighboring triangular elements sharing a common node (i.e., nodal star). Each MC had two degrees of freedom (DOFs), which was similar to the node property of the FEM. The mathematical Since the NMM does not require MC sides to coincide with the material boundaries and discontinuities, arbitrary shape MCs can be employed. For convenience, a regularly structured mathematical mesh was employed in the NMM. As seen in Figure 2a, a regular triangular mesh was constructed, in which each MC was formed through six neighboring triangular elements sharing a common node (i.e., nodal star). Each MC had two degrees of freedom (DOFs), which was similar to the node property of the FEM. The mathematical mesh covered the whole physical domain to form a PC system, where the common areas denoted by E i (i ∈ Natural number) was formed by the neighboring three hexagonal PCs. When the linear triangular element PU function is applied into the cover system, the global displacement function over a ME can be expressed as: where w i (x, y) is the PU function over the three associated PCs; FEM is a special case of the NMM when the cover function is constant and the weight function is the shape function in the FEM. In this study, we considered a 2-dimensional problem based on regular triangular grids. The finite element shape function of the triangle naturally forms a uniform division (PU) from the mathematical patch, as shown in Figure 2b. As shown in the following case, consider a ME ( ∈ ) formed by three associated PCs, three overlapping PCs of PU functions on a regular hexagonal patch, where each PC consists of six equilateral triangle elements, called flat-top PU [33]. The images and profiles of the PU function are shown in Figure 2b,c, respectively. Using the method in [34], we proved that the flat-top PU approximation built by the triangular mesh was linearly independent at the global level.

Open-Close Iteration Algorithm
In the original NMM computations, when the contact criteria of no-penetration and no-tension could not be satisfied within six iterations in one open-close iteration (OCI), the calculated time step was reduced to one-third of the initial time step, and then the OCI was restarted again. The number of OCIs also depends on the input parameters (e.g., the  In the two-dimensional NMM, a PU function w i (x, y) is defined on a PC P i such that The PU function is a partition of unity and satisfies Using the weight function w i (x, y), the global displacement function U(x, y) and V(x, y) on the whole PC system can be defined based on the following local displacement functions u i (x, y) and v i (x, y), respectively, as where n is the number of PCs equal to 3 for a 2D problem. When the triangular finite element mesh and a constant cover function are used, the global displacement on each ME can be rewritten as Thus, the FEM is a special case of the NMM when the cover function is constant and the weight function is the shape function in the FEM. In this study, we considered a 2-dimensional problem based on regular triangular grids. The finite element shape function of the triangle naturally forms a uniform division (PU) from the mathematical patch, as shown in Figure 2b. As shown in the following case, consider a ME E i (i ∈ N) formed by three associated PCs, three overlapping PCs of PU functions on a regular hexagonal patch, where each PC consists of six equilateral triangle elements, called flat-top PU [33]. The images and profiles of the PU function are shown in Figure 2b,c, respectively. Using the method in [34], we proved that the flat-top PU approximation built by the triangular mesh was linearly independent at the global level.

Open-Close Iteration Algorithm
In the original NMM computations, when the contact criteria of no-penetration and no-tension could not be satisfied within six iterations in one open-close iteration (OCI), the calculated time step was reduced to one-third of the initial time step, and then the OCI was restarted again. The number of OCIs also depends on the input parameters (e.g., the time step). Larger time steps can lead to more OCIs to achieve contact criteria [10]. Highperformance computing enables the NMM to challenge structural dynamic simulations including numerous discontinuous blocks. Nevertheless, the OCIs in NMM met the requirements of non-free and tension-free criteria. The correctness of the contact judgment criterion and the computational cost of the contact requirements at each time step control the entire computational efficiency of the NMM.
The traditional NMM uses a penalty method with parameter p and gap function of G = G(Λ), where Λ is the gap space vector for the detected contact pairs. As plotted in Figure 3, given a contact judgment distance 2ρ c to detect the potential contact boundary, where ρ c = max ∑ u l (x, y) 2 , ∀(x, y) ∈ B l , l = 1, 2, · · · , r; u l (x, y) is the global approximation to the block boundary B l ; and r is the total number of potential contact where k is the number of iterations for the current OCI and d (k) is the penetration distance at the kth iteration. K (k) = K (k) forms the global stiffness term, and F (k) = F (k) generates the global force vectors for the equilibrium equation of motion system. The next iterative displacements of D (k+1) over the iterative number k satisfy the relational expression of where . Within the number of iterations of k ≤ 6, the final condition for OCI convergence at where d (0) is the initial search penetration distance of the contact pair, and δ with the set tolerance value between 10 −4 and 10 −6 . Therefore, an appropriate contact stiffness p or time step is required to prevent poor contact penetration.
the entire computational efficiency of the NMM. The traditional NMM uses a penalty method with parameter and gap function of = ( ), where is the gap space vector for the detected contact pairs. As plotted in Figure 3, given a contact judgment distance 2 to detect the potential contact boundary, where = {√∑ ( , ) 2 }, ∀( , ) ∈ , = 1, 2, ⋯ , ; ( , ) is the global approximation to the block boundary ; and is the total number of potential contact boundaries. The contact stiffness and contact force can be determined using the penalty method as the form of where is the number of iterations for the current OCI and ( ) is the penetration distance at the th iteration. ( ) ⇢ ̿ ( ) forms the global stiffness term, and ( ) ⇢ ̿ ( ) generates the global force vectors for the equilibrium equation of motion system. The next iterative displacements of ( +1) over the iterative number satisfy the relational expression of where ̅ ( ) = ̅ ( −1) + ( ) ̿ ( ) , and ̅ ( ) = ̅ ( −1) + ( ) ̿ ( ) for the initiation of ̅ (0) = ̿ (0) and ̅ (0) = ̿ (0) . Within the number of iterations of ≤ 6, the final condition for OCI convergence at ( ) = 0 A is where (0) is the initial search penetration distance of the contact pair, and with the set tolerance value between 10 −4 and 10 −6 . Therefore, an appropriate contact stiffness or time step is required to prevent poor contact penetration.

The G-α Time Integration
The motion equation of linear dynamic problems can be expressed in the following under the consideration of the initial conditions:

The G-α Time Integration
The motion equation of linear dynamic problems can be expressed in the following under the consideration of the initial conditions: where M, C, and K are the mass, damping, and stiffness matrices, respectively; F(t) is the vector of applied load; ..       u is the increment of the acceleration. Substituting Equations (10a), (10b) and (10c) into Equation (9), the discrete form of the equation can be re-expressed as follows: (11) where F n+1 denotes the force vector at time t n+1 = (n + 1)∆t. In the traditional NMM, an implicit Newmark time integration is conducted by minimizing the potential energy with parameters β = 1/2, γ = 1. Thus, substituting the parameters into Equation (11), the equation of motion is rewritten as where ∆u = (u n+1 − u n ) is the displacement increment from the time step t n to (t n + ∆t).
The mass matrix of M and stiffness matrix of K inherits the properties of symmetry and sparsity similar to the FEM formula. Since the implicit algorithm is required to assemble the global stiffness matrix and solve the coupled system equation, the computational cost is raised, especially when the OCIs of contact treatment are applied in the NMM code. Moreover, an appropriate ∆t is determined under the requirement of numerical damping. For the best results from stability and numerical dissipation point of view, a time step ∆t is considered as [35]: where ω e max is the element maximum angle frequency of system. The G-α method in the present study is an unconditionally stable, second-order accuracy time integration, which possesses an optimal combination of high-frequency and low-frequency dissipation. Formula of the G-α method resorting to the Newmark method frame can be expressed as follows: where where n ∈ {0, 1, · · · , N − 1}; N is the total number of time step; coefficients α m and α f are the two control parameters of the algorithm; and .. u n+1 , u n and u n+1 , are the accelerations, velocities, and displacements at time t n and t n+1 , respectively. These approximations satisfy the Newmark method to be expressed as where β and γ are the two parameters determining the algorithm property and θ is the kinetic damping to discount velocity and satisfies θ ∈ [0, 1]. Accordingly, the global equilibrium equation can be obtained when θ equals to 1, and the more details of the derivation can be referred to in Appendix A. The G-α method is of second-order accuracy for which achieves optimal high frequency dissipation with a minimal low frequency effect to damp spurious oscillations when the following conditions hold: With an appropriate choice of parameters of β and γ, the G-α method simplifies to the HHT-α method [4] in the case of α m = 0 (i.e., G-α 1 hereafter), the WBZ-α method [5] in the case of α f = 0 (i.e., G-α 2 hereafter), and the Newmark family method [2] when α m = α f = 0. The original NMM code adopts the Newmark method with α m = α f = 0 and γ = 2β = 1, which can be assumed as a special simplified case of the G-α method.

Stability Analysis of the G-α Method
For a single degree of freedom (SDOF) system, the motion equation can be represented by ..
where ζ = C/ 2 √ KM is the damping ratio; ω = √ K/M is the natural frequency; and q represents the mass loading.
The G-α algorithm described in Equation (19) can be investigated using the analytical method [1,8]. For this purpose, the algorithm can be converted to another form: where A is the amplification matrix that determines algorithmic characteristics such as stability, accuracy, and numerical dissipation; L denotes the load operator; and q n and q n+1 are the loading at time t n and t n+1 , respectively. The amplification matrix A is obtained by the derivations of Equation (20), expressed as: with The eigenvalues of A are determined by the characteristic equation of where To measure the accuracy of the time integration, two solutions of the characteristic Equation (23) are considered: two complex conjugate roots (denoted as λ 1 , λ 2 ) and one real root (denoted as λ 3 ), which satisfies the condition of Specifically, a and b are the real parts and imaginary parts, i = √ −1, or the three real roots of λ 1 , λ 2 , and λ 3 .
The spectral radius ρ(A), the algorithmic damping ζ, and the relative period error E r are commonly used as criteria for the performance of an algorithm in comparisons with single step integration methods [19]. These parameters can be defined as: where Using the procedure in [8], the stability criterion |λ i | ≤ 1, i = 1, 2, 3, is satisfied if the following five conditions are constructed (i.e., space of stability in Figure 4):  (29-d)and(29-e)) all satisfy the requirement of unconditional stability of the G-α algorithm. Therefore, the G-α method is unconditionally stable when the parameters of = 0.3025 and = 0.6; since as described in [19], the NMK-α method used in the original NMM is unconditionally stable with the parameters of = 0.5 and = 1. Furthermore, if 2 < and ≥ 1 2 are both satisfied, the algorithm is conditionally stable. Under this condition, ≤ must be satisfied in the requirement of ( ) ≤ 1, where is the critical sampling frequency and can be expressed as: Suppose that θ [0, 1] and 2β ≥ γ ≥ 1 2 , the above inequalities (i.e., in Equation (29a-e)) all satisfy the requirement of unconditional stability of the G-α algorithm. Therefore, the G-α method is unconditionally stable when the parameters of β = 0.3025 and γ = 0.6; since as described in [19], the NMK-α method used in the original NMM is unconditionally stable with the parameters of β = 0.5 and γ = 1. Furthermore, if 2β < γ and γ ≥ 1 2 are both satisfied, the algorithm is conditionally stable. Under this condition, Ω ≤ Ω crit must be satisfied in the requirement of ρ(A) ≤ 1, where Ω crit is the critical sampling frequency and can be expressed as: In the present study, the G-α 0 method (e.g., see Table 1), combined with the G-α 1 and G-α 2 methods and the Newmark method of γ = 2 β = 1 (i.e., denoted as NMK-α), was employed to investigate the algorithm's performance by the NMM code; more details of the parameters used in the G-α and NMK-α methods can be referred to in Table 1. Table 1. Selected parameters for the referred methods.

The Referred Methods
As can be seen in Figure 5, when different values of kinetic damping θ are employed (i.e., θ = 0, 0.2, 0.4, 0.6, 0.8, 1), the spectral radius ρ(A) of the referred algorithms can be calculated as presented in Table 2: Table 2. Spectral radius ρ(A) calculated by the referred algorithms.

The Methods
It can be found that the G-α 0 method possessed similar numerical stability and dissipation to the G-α 1 and G-α 2 methods when the same parameters of β and γ were applied. On the other hand, when θ = 0 was used in the NMK-α method, ρ(A) → 0 could be obtained and the maximum algorithm damping ζ was considered to achieve the fasted convergence velocity, but the numerical accuracy could not be guaranteed. In terms of the period for optimization of the above-mentioned time integration schemes, the G-α method possessed more inbuilt advantage than that of the NMK-α method. Furthermore, the value of Ω bif = 4 was also computed in Figure 5d when the NMK-α method was applied to obtain the highest convergences efficiency, where Ω bif is the bifurcate sampling frequency.
To avoid the decay and saw tooth pattern oscillation, the solution of u n = ∑ 3 i=1 c i λ n i must satisfy the conditions of |λ 2 | > |λ 1 | > |λ 3 | and λ 2 < 0. Moreover, Ω < Ω bi f or θ > θ crit is hold to satisfy the stability condition, where θ crit is called critical kinetic damping. If the value of Ω is too large or θ is minor, oscillation will occur, the contact and deformation will also be distorted when the time integration of the G-α 0 method is employed into the computations. At this case, more time steps are required to achieve the solution and OCI convergence. Therefore, if Ω = Ω bi f or θ = θ crit is adopted, the convergence velocity is the fastest and computational efficiency is the highest. As plotted in Figure 6a, the critical kinetic damping θ crit of the G-α 1 method, the G-α 2 method and the G-α 0 method approach to 0.4407 when Ω → ∞ . On the other hand, the NMK-α is different since bifurcation condition has a significant effect on the solution. θ crit satisfies achieved at the case of Ω ∈ Ω bi f , ∞), where the value of Ω bi f is equal to 4. Obviously, it is more complicated to determine time step ∆t by Ω bi f than cases that considering viscous damping in the NMK-α method, ∆t = Ω bi f /ω normally is larger than the used in practical discontinuous deformation simulations. Thus, the critical kinetic damping θ crit is more preferred than Ω bi f in the actual computations. Figure 6b shows the spectral radius ρ by the θ crit when the different methods are considered. Both the G-α 1 method and the G-α 2 method possess the same spherical radius ρ → 0.54 at the cases of Ω → ∞, θ crit → 0.4407 , but the G-α 0 method is different, which satisfies the condition of Ω → ∞, ρ → 0.818 when the same θ crit is selected. For the NMK-α method, when Ω bi f = 4 is taken into account, θ crit = 1 and ρ = 0.333 can be only determined by the amplification matrix of A; Ω → ∞, ρ → 0 and Ω → 0, ρ → 1 are also obtained when different θ crit values are chosen. It is also found that the NMK-α method achieves the maximum θ crit comparing with the other three methods, and minimum spectral radius ρ to attain the fastest convergence of the solution. It can be found that the G-α 0 method possessed similar numerical stability and dissipation to the G-α 1 and G-α 2 methods when the same parameters of β and γ were applied. On the other hand, when = 0 was used in the NMK-α method, ρ(A) → 0 could be obtained and the maximum algorithm damping ζ ̅ was considered to achieve the fasted convergence velocity, but the numerical accuracy could not be guaranteed. In terms of the period for optimization of the above-mentioned time integration schemes, the G-α method possessed more inbuilt advantage than that of the NMK-α method. Furthermore, the value of bif = 4 was also computed in Figure 5d when the NMK-α method was applied to obtain the highest convergences efficiency, where bif is the bifurcate sampling frequency. → ∞, → 0 and → 0, → 1 are also obtained when different values are chosen. It is also found that the NMK-α method achieves the maximum comparing with the other three methods, and minimum spectral radius to attain the fastest convergence of the solution.

Accuracy of the G-α Method
Generally, parameters of algorithmic damping ratio ̅ and relative period error are used to measure the numerical dissipation and dispersion of the time integration scheme. In particular, the values of ̅ in the cases of = 0 and = 1 were plotted in Figure 7 to evaluate the G-α method algorithm accuracy. When = 0 was considered,

Accuracy of the G-α Method
Generally, parameters of algorithmic damping ratio ζ and relative period error E r are used to measure the numerical dissipation and dispersion of the time integration scheme. In particular, the values of ζ in the cases of θ = 0 and θ = 1 were plotted in Figure 7 to evaluate the G-α method algorithm accuracy. When θ = 0 was considered, the NMK-α method obtained ζ = 0; the G-α method possessing positive ζ implies that the computed result is convergence. ζ = 0 means that no numerical dissipation exists by the algorithm, and that spurious oscillations will occur in the computations. When θ = 1 was suggested, the NMK-α method possessed the maximum ζ so that the numerical solution approached an exact solution. The G-α 1 method and the G-α 2 method had higher ζ than that of the G-α 0 method in both cases. Thus, it can be concluded that when θ crit is employed, the optimized values of ζ can be observed to control the numerical accuracy of the G-α method. the computed result is convergence. ̅ = 0 means that no numerical dissipation exists by the algorithm, and that spurious oscillations will occur in the computations. When = 1 was suggested, the NMK-α method possessed the maximum ̅ so that the numerical solution approached an exact solution. The G-α 1 method and the G-α 2 method had higher ̅ than that of the G-α 0 method in both cases. Thus, it can be concluded that when is employed, the optimized values of ̅ can be observed to control the numerical accuracy of the G-α method.

Numerical Examples
To present the performance of the G-α method using the NMM, the aforementioned G-α 0 , G-α 1 , G-α 2 , and NMK-α methods were considered to calibrate the G-α method first, in which one transient dynamical response problem of an elastic bar was simulated. Then, a two-layer system of rock beam deflection was computed to validate the proposed G-α

Numerical Examples
To present the performance of the G-α method using the NMM, the aforementioned G-α 0 , G-α 1 , G-α 2 , and NMK-α methods were considered to calibrate the G-α method first, in which one transient dynamical response problem of an elastic bar was simulated. Then, a two-layer system of rock beam deflection was computed to validate the proposed G-α method. Finally, the mechanical behavior of dynamical tunneling was studied using the proposed G-α method to further explore the computational accuracy and efficiency of the proposed algorithm.

Transient Response of the Elastic Bar
To evaluate the accuracy of the proposed G-α method using the NMM, the transient response of an elastic rock bar was computed. As shown in Figure 8a, a 1.0 m × 2.0 m elastic bar bonded at the left end, and a harmonic function at the right end was applied, which represents the sine function form of P(t) = P 0 sin ω 0 t, as shown in Figure 8b where is the density of strip; Φ( ) is the mode shape function vector; is the bending stiffness vector; is the length of the strip; and is the Young's modulus.   The theoretical solution of deformation u(t) to the elasticity problem can be expressed as follows [36]: where η = 1 − α 2 2 + 4ζ 2 α 2 − 1 2 is the dynamic magnification factor; δ is the end deflection under unit loading and satisfies expression of δ = L 3 /3EI; and ψ and ϑ are the two angle phases, which satisfies where α = ω 0 /ω; ω = √ K * /M * is the natural frequency; ω D = ω 1 − ζ 2 is the damped natural frequency; ζ is the damping ratio; and the effective mass M * and stiffness K * are expressed as: where ρ is the density of strip; Φ(x) is the mode shape function vector; EI is the bending stiffness vector; L is the length of the strip; and E is the Young's modulus. Figure 9 shows the simulations of the elastodynamic problem by the referred time integration schemes with time steps ∆t = 0.01 s and ∆t = 0.001 s, respectively. As presented in Figure 9a, the numerical damping was found by both the G-α method and three other methods, and the NMK-α used in the NMM code possessed the strongest damping effect to the simulations. When a smaller time step (i.e., ∆t = 0.001 s) is employed, as plotted in Figure 9b, an analytical solution can be approached exactly. However, when a large time step (i.e., = 0.01 s) is adopted, stray oscillations occur, and the simulation results are dispersed into the analytical solution.
To verify the computational accuracy to the theoretical solution, the crests' displacements were calculated using the proposed methods (i.e., the G-0 , G-1 , G-2 , and NMKα methods), and the simulated results are listed in Table 3. It was found that deformation decreased from the first cases of ∆ = 0.01 s and∆ = 0.001 s. The results from the G-α methods were more accurate than that of the NMK-α method. The G-0 method achieved 18.6958 mm at the first crest compared to the NMK-α method of 17.5063 mm when ∆ = 0.01 s was used. When∆ = 0.001 s was used in the simulations, the G-0 method approached 18.9768 mm compared to 18.1625 mm of the NMK-α method, while the analytical solution of 19.3378 mm was calculated. The maximum relative errors of the NMK-α method to the analytical solution were 9.4711% and 6.0773% in the different cases.

Layered Rock Deflection Simulations
This example analyzed a two-layer system, as shown in Figure 10. The block system was considered as a composite beam with both ends clamped, consisting of a thin beam  However, when a large time step (i.e., t = 0.01 s) is adopted, stray oscillations occur, and the simulation results are dispersed into the analytical solution.
To verify the computational accuracy to the theoretical solution, the crests' displacements were calculated using the proposed methods (i.e., the G-α 0 , G-α 1 , G-α 2 , and NMK-α methods), and the simulated results are listed in Table 3. It was found that deformation decreased from the first cases of ∆t = 0.01 s and ∆t = 0.001 s. The results from the G-α methods were more accurate than that of the NMK-α method. The G-α 0 method achieved 18.6958 mm at the first crest compared to the NMK-α method of 17.5063 mm when ∆t = 0.01 s was used. When ∆t = 0.001 s was used in the simulations, the G-α 0 method approached 18.9768 mm compared to 18.1625 mm of the NMK-α method, while the analytical solution of 19.3378 mm was calculated. The maximum relative errors of the NMK-α method to the analytical solution were 9.4711% and 6.0773% in the different cases.

Layered Rock Deflection Simulations
This example analyzed a two-layer system, as shown in Figure 10. The block system was considered as a composite beam with both ends clamped, consisting of a thin beam with Young's modulus E 1 , thickness h 1 , unit weight γ 1 , and a thick beam with Young's modulus E 2 , thickness h 2 , and unit weight γ 2 . The thin beam was overlaid above the thick beam. Assuming only under gravity, if the deflection of the thin beam is larger than that of the thick beam, the thin beam will push down the bottom thick beam.  An analytical solution for the system can be referred to in [37] when only the selfweigh of the beam was considered, which is expressed as = 2 32 ℎ 2 (34) where is the maximum deflection of the beam; is the unit weight; and L and h are the length and thickness of the beam, respectively. The derived deflection of the bottom of the system is calculated by assigning it an increased unit weight , given by: In this simulation, the length of the composite beam was 36 m, and the thicknesses of the upper beam and the lower beam were ℎ 1 = 2 m and ℎ 2 = 3 m, respectively. Both beams were assumed to have the same material properties: Young's modulus 1 = 2 = 4.8 GPa, Poisson's ratio of ν 1 = ν 2 = 0.25, and unit weight γ 1 = γ 2 = 25,500 N/m 3 . Then, the time step size of 0.001 s was used in the modeling. The interface between two layers was assumed to be smooth with zero friction angle and cohesion. In the NMM model, as shown in Figure 11, the composite beam was divided into 1126 manifold elements, and 19 measured points with a spacing of 2 m were prescribed to measure the deflection of the two-layer system. According to the analytical solution in Goodman [37], it can also be seen that a good agreement was achieved. The maximum deflection of the composite beam was 0.0408 m. After checking the displacements of the measured points from the NMM model, it was found that the maximum deflections of the bottom were both close to the analytical solution with 0.0414 m by the proposed G-α 0 method and 0.0408 m using the original NMKα method. The vertical deflections of these measured points from the NMM model were also plotted in Figure 12 against those from the NMM. An analytical solution for the system can be referred to in [37] when only the self-weigh of the beam was considered, which is expressed as where D max is the maximum deflection of the beam; γ e is the unit weight; and L and h are the length and thickness of the beam, respectively. The derived deflection of the bottom of the system is calculated by assigning it an increased unit weight γ e , given by: In this simulation, the length of the composite beam was 36 m, and the thicknesses of the upper beam and the lower beam were h 1 = 2 m and h 2 = 3 m, respectively. Both beams were assumed to have the same material properties: Young's modulus E 1 = E 2 = 4.8 GPa, Poisson's ratio of ν 1 = ν 2 = 0.25, and unit weight γ 1 = γ 2 =25,500 N/m 3 . Then, the time step size of 0.001 s was used in the modeling. The interface between two layers was assumed to be smooth with zero friction angle and cohesion. In the NMM model, as shown in Figure 11, the composite beam was divided into 1126 manifold elements, and 19 measured points with a spacing of 2 m were prescribed to measure the deflection of the two-layer system.  An analytical solution for the system can be referred to in [37] when only the selfweigh of the beam was considered, which is expressed as = 2 32 ℎ 2 (34) where is the maximum deflection of the beam; is the unit weight; and L and h are the length and thickness of the beam, respectively. The derived deflection of the bottom of the system is calculated by assigning it an increased unit weight , given by: In this simulation, the length of the composite beam was 36 m, and the thicknesses of the upper beam and the lower beam were ℎ 1 = 2 m and ℎ 2 = 3 m, respectively. Both beams were assumed to have the same material properties: Young's modulus 1 = 2 = 4.8 GPa, Poisson's ratio of ν 1 = ν 2 = 0.25, and unit weight γ 1 = γ 2 = 25,500 N/m 3 . Then, the time step size of 0.001 s was used in the modeling. The interface between two layers was assumed to be smooth with zero friction angle and cohesion. In the NMM model, as shown in Figure 11, the composite beam was divided into 1126 manifold elements, and 19 measured points with a spacing of 2 m were prescribed to measure the deflection of the two-layer system. According to the analytical solution in Goodman [37], it can also be seen that a good agreement was achieved. The maximum deflection of the composite beam was 0.0408 m. After checking the displacements of the measured points from the NMM model, it was found that the maximum deflections of the bottom were both close to the analytical solution with 0.0414 m by the proposed G-α 0 method and 0.0408 m using the original NMKα method. The vertical deflections of these measured points from the NMM model were also plotted in Figure 12 against those from the NMM. According to the analytical solution in Goodman [37], it can also be seen that a good agreement was achieved. The maximum deflection of the composite beam was 0.0408 m. After checking the displacements of the measured points from the NMM model, it was found that the maximum deflections of the bottom were both close to the analytical solution with 0.0414 m by the proposed G-α 0 method and 0.0408 m using the original NMK-α method. The vertical deflections of these measured points from the NMM model were also plotted in Figure 12 against those from the NMM. It could also be seen that good agreement was achieved. The contour of the lower beam obtained by the G-α 0 and the NMK-α method are presented in Figures 13 and 14, respectively. In terms of efficiency, the G-α 0 and the NMK-α methods were employed to analyze the example on the same computer with a processor speed of 4.0 GHz and installed memory of 8 GB. The CPU time was 5.22 and 6.36 min for both methods, respectively. It can be concluded that given the same accuracy, the G-α 0 method was relatively more efficient than the original NMK-α method by the NMM.

Simulation of Behaviors of Blocks on Shaking Table
In the simulations, as plotted in Figure 15, an assembled block system consisted of three cubic concrete blocks, the size of each block was 15 cm × 15 cm × 15 cm, and the density was 2200 kg/m 3 . The numerical model is shown in Figure 16   It could also be seen that good agreement was achieved. The σ x contour of the lower beam obtained by the G-α 0 and the NMK-α method are presented in Figures 13 and 14, respectively. In terms of efficiency, the G-α 0 and the NMK-α methods were employed to analyze the example on the same computer with a processor speed of 4.0 GHz and installed memory of 8 GB. The CPU time was 5.22 and 6.36 min for both methods, respectively. It can be concluded that given the same accuracy, the G-α 0 method was relatively more efficient than the original NMK-α method by the NMM. It could also be seen that good agreement was achieved. The contour of the lower beam obtained by the G-α 0 and the NMK-α method are presented in Figures 13 and 14, respectively. In terms of efficiency, the G-α 0 and the NMK-α methods were employed to analyze the example on the same computer with a processor speed of 4.0 GHz and installed memory of 8 GB. The CPU time was 5.22 and 6.36 min for both methods, respectively. It can be concluded that given the same accuracy, the G-α 0 method was relatively more efficient than the original NMK-α method by the NMM.

Simulation of Behaviors of Blocks on Shaking Table
In the simulations, as plotted in Figure 15, an assembled block system consisted of three cubic concrete blocks, the size of each block was 15 cm × 15 cm × 15 cm, and the density was 2200 kg/m 3 . The numerical model is shown in Figure 16 and the physical properties of block and parameters used in this analysis were: Poisson's ratio of 0.25, friction angle of 36.4°, Young's modulus of 14.9 GPa, and contact stiffness of 200 GPa, respectively. Then, the time step size of 0.001 s was applied to both methods and the results of the simulation in the case of a frequency of 5.0 Hz and amplitude of 0.7138 g vibration are shown in Figure 17. It could also be seen that good agreement was achieved. The contour of the lower beam obtained by the G-α 0 and the NMK-α method are presented in Figures 13 and 14, respectively. In terms of efficiency, the G-α 0 and the NMK-α methods were employed to analyze the example on the same computer with a processor speed of 4.0 GHz and installed memory of 8 GB. The CPU time was 5.22 and 6.36 min for both methods, respectively. It can be concluded that given the same accuracy, the G-α 0 method was relatively more efficient than the original NMK-α method by the NMM.

Simulation of Behaviors of Blocks on Shaking Table
In the simulations, as plotted in Figure 15, an assembled block system consisted of three cubic concrete blocks, the size of each block was 15 cm × 15 cm × 15 cm, and the density was 2200 kg/m 3 . The numerical model is shown in Figure 16 and the physical properties of block and parameters used in this analysis were: Poisson's ratio of 0.25, friction angle of 36.4°, Young's modulus of 14.9 GPa, and contact stiffness of 200 GPa, respectively. Then, the time step size of 0.001 s was applied to both methods and the results of the simulation in the case of a frequency of 5.0 Hz and amplitude of 0.7138 g vibration are shown in Figure 17.

Simulation of Behaviors of Blocks on Shaking Table
In the simulations, as plotted in Figure 15, an assembled block system consisted of three cubic concrete blocks, the size of each block was 15 cm × 15 cm × 15 cm, and the density was 2200 kg/m 3 . The numerical model is shown in Figure 16 Figure 17.
Compared with the experimental test results [38] (i.e., Figure 18) and the proposed G-α method using the NMM (i.e., Figure 19), it was found that the G-α 0 method could be used to simulate sliding and rotation (locking) of the block seen in the test well. Figure 20 shows the results of horizontal displacement and acceleration spectrum of the top block of the numerical methods and experiments. From these data, horizontal displacement and acceleration spectrum can be well expressed by the G-α 0 method. Thus, it was proven that the proposed G-α 0 method could well simulate the vibration response of assembled blocks. On the other hand, it was also proven that the G-α 0 method could be used to simulate the vibration response of the assembled blocks accurately.   Compared with the experimental test results [38] (i.e., Figure 18) and the p G-α method using the NMM (i.e., Figure 19), it was found that the G-α0 method used to simulate sliding and rotation (locking) of the block seen in the test well. F shows the results of horizontal displacement and acceleration spectrum of the to of the numerical methods and experiments. From these data, horizontal displacem acceleration spectrum can be well expressed by the G-α0 method. Thus, it was pro the proposed G-α0 method could well simulate the vibration response of as blocks. On the other hand, it was also proven that the G-α0 method could be used ulate the vibration response of the assembled blocks accurately.   Compared with the experimental test results [38] (i.e., Figure 18) and the p G-α method using the NMM (i.e., Figure 19), it was found that the G-α0 method c used to simulate sliding and rotation (locking) of the block seen in the test well. F shows the results of horizontal displacement and acceleration spectrum of the to of the numerical methods and experiments. From these data, horizontal displacem acceleration spectrum can be well expressed by the G-α0 method. Thus, it was pro the proposed G-α0 method could well simulate the vibration response of as blocks. On the other hand, it was also proven that the G-α0 method could be used ulate the vibration response of the assembled blocks accurately.   Compared with the experimental test results [38] (i.e., Figure 18) and the proposed G-α method using the NMM (i.e., Figure 19), it was found that the G-α0 method could be used to simulate sliding and rotation (locking) of the block seen in the test well. Figure 20 shows the results of horizontal displacement and acceleration spectrum of the top block of the numerical methods and experiments. From these data, horizontal displacement and acceleration spectrum can be well expressed by the G-α0 method. Thus, it was proven that the proposed G-α0 method could well simulate the vibration response of assembled blocks. On the other hand, it was also proven that the G-α0 method could be used to simulate the vibration response of the assembled blocks accurately.    Compared with the experimental test results [38] (i.e., Figure 18) and the proposed G-α method using the NMM (i.e., Figure 19), it was found that the G-α0 method could be used to simulate sliding and rotation (locking) of the block seen in the test well. Figure 20 shows the results of horizontal displacement and acceleration spectrum of the top block of the numerical methods and experiments. From these data, horizontal displacement and acceleration spectrum can be well expressed by the G-α0 method. Thus, it was proven that the proposed G-α0 method could well simulate the vibration response of assembled blocks. On the other hand, it was also proven that the G-α0 method could be used to simulate the vibration response of the assembled blocks accurately.

Simulation of the Mechanical Behavior of Dynamical Tunneling
In this simulation, the mechanical behavior of inclined rock mass during dynamic tunnel construction was studied. Rock fractures are composed of many discontinuous, discrete blocks. The effect of the stress arch on the stress distribution and surface subsidence during dynamic tunneling was studied, and the calculation results compared with previous experiments [39]. Figure 21 presents a simplified model of jointed rock masses with a dip angle of 0 • . Both the depth of the discontinuous blocks and the width of trap door B were assumed to 200 mm. The input parameters of the modeling were: the Young's modulus of 62 MPa, Poisson's ratio of 0.15, unit density of 26.4 kN/m 3 , and the internal frictional angle of 20 • . In the simulations, the effect of the deposit blocks is denoted by 10, 11, 12, and 13 on the vertical stress around the area of rock masses bottom was studied using the G-α and the three other methods, respectively. As presented in Figure 21b, 705 NMM elements and 130 blocks were divided using the FEM meshing technique, where a 2 mm gap was prescribed between the bottom of the blocks and trap door. Furthermore, in order to achieve numerical stability and save computational time, the kinetic damping of θ crit−1 = 0.8033 was used in the G-α 0 , G-α 1 , and G-α 2 methods, and θ crit−2 = 0.8080 was applied in the NMK-α method, respectively.

Simulation of the Mechanical Behavior of Dynamical Tunneling
In this simulation, the mechanical behavior of inclined rock mass during dynamic tunnel construction was studied. Rock fractures are composed of many discontinuous, discrete blocks. The effect of the stress arch on the stress distribution and surface subsidence during dynamic tunneling was studied, and the calculation results compared with previous experiments [39]. Figure 21 presents a simplified model of jointed rock masses with a dip angle of 0°. Both the depth of the discontinuous blocks and the width of trap door B were assumed to 200 mm. The input parameters of the modeling were: the Young's modulus of 62 MPa, Poisson's ratio of 0.15, unit density of 26.4 kN/m 3 , and the internal frictional angle of 20°. In the simulations, the effect of the deposit blocks is denoted by 10, 11, 12, and 13 on the vertical stress around the area of rock masses bottom was studied using the G-α and the three other methods, respectively. As presented in Figure 21b, 705 NMM elements and 130 blocks were divided using the FEM meshing technique, where a 2 mm gap was prescribed between the bottom of the blocks and trap door. Furthermore, in order to achieve numerical stability and save computational time, the kinetic damping of θ crit−1 = 0.8033 was used in the G-α0, G-α1, and G-α2 methods, and θ crit−2 = 0.8080 was applied in the NMK-α method, respectively.  Figure 16 presents the normal stress distribution simulated by the referred methods using ∆ = 0.1 and∆ = 0.01 ms, respectively. It was found that the normal stress distributions were symmetrical, and the simulated results were in agreement with the experiments. When ∆ = 0.01 ms was used in the computations, the local yield section [40] was found by the three referred algorithms. On the other hand, when a larger time step∆ = 0.1 ms was applied to the modeling, the spurious oscillation became weaker and the simulated results were more accurate with the G-α method than that of the NMK-α method.
When a smaller step time (i.e., ∆ = 0.01 ms) was employed, as shown in Figure 22b, the simulated results of the G-α0, G-α1, and G-α2 methods were more stable than that of the NMK-α method. Figure 23 shows the simulated results of the surface subsidence by the above referred algorithms using the different time steps (i.e., ∆ = 0.1 ms and ∆ = 0.01 ms). Due to the numerical stability properties of the algorithms, the simulated results of the G-α algorithm still oscillated, but approached the experiments more than that of the NMK-α method.  Figure 16 presents the normal stress distribution simulated by the referred methods using ∆t = 0.1 and ∆t = 0.01 ms, respectively. It was found that the normal stress distributions were symmetrical, and the simulated results were in agreement with the experiments. When ∆t = 0.01 ms was used in the computations, the local yield section [40] was found by the three referred algorithms. On the other hand, when a larger time step ∆t = 0.1 ms was applied to the modeling, the spurious oscillation became weaker and the simulated results were more accurate with the G-α method than that of the NMK-α method. When a smaller step time (i.e., ∆t = 0.01 ms) was employed, as shown in Figure 22b, the simulated results of the G-α 0 , G-α 1 , and G-α 2 methods were more stable than that of the NMK-α method. Figure 23 shows the simulated results of the surface subsidence by the above referred algorithms using the different time steps (i.e., ∆t = 0.1 ms and ∆t = 0.01 ms). Due to the numerical stability properties of the algorithms, the simulated results of the G-α algorithm still oscillated, but approached the experiments more than that of the NMK-α method. The computational cost of the OCIs using the referred time integration schemes was further calculated in the present study. As listed in Table 4, the CPU cost by the G-α algorithm significantly increased from 0.51839 to 6.90864 h while the number of OCIs from 13,537 to 285,273 with the time step ∆ = 1.0 ms to ∆ = 0.01 ms. Accordingly, the NMKα algorithm had the same increase in terms of the CPU cost and the OCIs. As can be seen in Figure 24, the G-α schemes were more efficient than that of the traditional NMK-α method used in the NMM, not only in terms of the CPU time, but also the OCIs when the smaller time step (i.e., ∆ = 0.01 ms) was used. The G-α0 method took a minimum of 6.90864 h compared with 8.13725 h (i.e., the maximum relative error was 15.10%) by the NMK-α method when the time step ∆ = 0.01 ms was used. The G-α scheme was more efficient with the increase in the OCIs caused by more discrete block contacts. These results suggest that except for the original NMK-α method used in the NMM, the G-algorithm can be effectively applied to models of tunnel behavior involving more complex discontinuous rock bodies.  The computational cost of the OCIs using the referred time integration schemes was further calculated in the present study. As listed in Table 4, the CPU cost by the G-α algorithm significantly increased from 0.51839 to 6.90864 h while the number of OCIs from 13,537 to 285,273 with the time step ∆ = 1.0 ms to ∆ = 0.01 ms. Accordingly, the NMKα algorithm had the same increase in terms of the CPU cost and the OCIs. As can be seen in Figure 24, the G-α schemes were more efficient than that of the traditional NMK-α method used in the NMM, not only in terms of the CPU time, but also the OCIs when the smaller time step (i.e., ∆ = 0.01 ms) was used. The G-α0 method took a minimum of 6.90864 h compared with 8.13725 h (i.e., the maximum relative error was 15.10%) by the NMK-α method when the time step ∆ = 0.01 ms was used. The G-α scheme was more efficient with the increase in the OCIs caused by more discrete block contacts. These results suggest that except for the original NMK-α method used in the NMM, the G-algorithm can be effectively applied to models of tunnel behavior involving more complex discontinuous rock bodies.  The computational cost of the OCIs using the referred time integration schemes was further calculated in the present study. As listed in Table 4, the CPU cost by the G-α algorithm significantly increased from 0.51839 to 6.90864 h while the number of OCIs from 13,537 to 285,273 with the time step ∆t = 1.0 ms to ∆t = 0.01 ms. Accordingly, the NMK-α algorithm had the same increase in terms of the CPU cost and the OCIs. As can be seen in Figure 24, the G-α schemes were more efficient than that of the traditional NMK-α method used in the NMM, not only in terms of the CPU time, but also the OCIs when the smaller time step (i.e., ∆t = 0.01 ms) was used. The G-α 0 method took a minimum of 6.90864 h compared with 8.13725 h (i.e., the maximum relative error was 15.10%) by the NMK-α method when the time step ∆t = 0.01 ms was used. The G-α scheme was more efficient with the increase in the OCIs caused by more discrete block contacts. These results suggest that except for the original NMK-α method used in the NMM, the G-algorithm can be effectively applied to models of tunnel behavior involving more complex discontinuous rock bodies.

Conclusions
In this paper, the generalized-α (G-α) time integration scheme was studied within the NMM frame, in which a thorough investigation was carried out to reveal the effects of the G-α0, G-α1, and G-α2 methods and the Newmark method (i.e., NMK-α method) on the numerical stability and accuracy for dynamic simulations. For the first elastodynamic problem, the simulated results by the G-α scheme was in good agreement with the analytical solution. It was also found that the advantages of the G-α algorithm were apparent when the open-close iterations (OCIs) were applied for accurate contact analysis of the layered rock deflection. Comparing the simulated results of the behaviors of an assembled block system under seismic effect with the experiments, it was also clear that the G-α method could be applied to model the complicated discontinuous structural dynamic problems. Comparing the simulation results of the dynamical tunneling mechanical behavior with the experimental results, we also found that the G-α method can be used to simulate complex discontinuous structure dynamic problems. Furthermore, the convergence efficiency of the OCIs of the NMM at different time steps was investigated, with a significant decrease in the CPU time and OCIs of the NMK-α algorithm. The efficiency of the G-α algorithm implemented in the NMM code was proven, and it can be predicted that the G-α algorithm can simulate large scale discontinuous dynamic problems involving more nonlinear contacts.

Conclusions
In this paper, the generalized-α (G-α) time integration scheme was studied within the NMM frame, in which a thorough investigation was carried out to reveal the effects of the G-α 0 , G-α 1 , and G-α 2 methods and the Newmark method (i.e., NMK-α method) on the numerical stability and accuracy for dynamic simulations. For the first elastodynamic problem, the simulated results by the G-α scheme was in good agreement with the analytical solution. It was also found that the advantages of the G-α algorithm were apparent when the open-close iterations (OCIs) were applied for accurate contact analysis of the layered rock deflection. Comparing the simulated results of the behaviors of an assembled block system under seismic effect with the experiments, it was also clear that the G-α method could be applied to model the complicated discontinuous structural dynamic problems. Comparing the simulation results of the dynamical tunneling mechanical behavior with the experimental results, we also found that the G-α method can be used to simulate complex discontinuous structure dynamic problems. Furthermore, the convergence efficiency of the OCIs of the NMM at different time steps was investigated, with a significant decrease in the CPU time and OCIs of the NMK-α algorithm. The efficiency of the G-α algorithm implemented in the NMM code was proven, and it can be predicted that the G-α algorithm can simulate large scale discontinuous dynamic problems involving more nonlinear contacts.

Informed Consent Statement: Not applicable.
Data Availability Statement: The authors confirm that the data supporting the findings of this study are available within the article.