**2. Model-Based Integration Algorithms**

To numerically obtain the structural responses induced by earthquakes or other dynamic actions, the following equation of motion (EOM) should be solved using a direct integration algorithm:

**.**

..

$$\mathbf{M}\mathbf{X}\_{i+1} + \mathbf{C}\mathbf{X}\_{i+1} + \mathbf{R}\_{i+1} = \mathbf{F}\_{i+1} \tag{1}$$

**.**

..

where **M** and **C** are the mass and damping matrices, respectively; **X** and **X** are the velocity and acceleration vectors, respectively; the subscripts *i* + 1 and *i* represent the next and current time steps, respectively; **R** and **F** are the restoring force and external force vectors, respectively. The restoring force **R***i*+<sup>1</sup> is generally displacement-dependent and can be degraded into **KX***i*+<sup>1</sup> if the structure is linear elastic, where **K** and **X** are the stiffness matrix and displacement vector, respectively.

The main task for an integration algorithm is to predict the velocity and displacement at the next step. For instance, the difference equations of the classical Newmark family of algorithms are expressed as:

$$\dot{\mathbf{X}}\_{i+1} = \dot{\mathbf{X}}\_{i} + \Delta t \left[ (1 - \gamma) \ddot{\mathbf{X}}\_{i} + \gamma \ddot{\mathbf{X}}\_{i+1} \right], \\ \mathbf{X}\_{i+1} = \mathbf{X}\_{i} + \Delta t \dot{\mathbf{X}}\_{i} + \Delta t^{2} \left[ (1/2 - \beta) \ddot{\mathbf{X}}\_{i} + \beta \ddot{\mathbf{X}}\_{i+1} \right] \tag{2}$$

where Δ*t* is the time step; *γ* and *β* are two integration parameters. The Newmark algorithms are explicit only if *β* = 0. Furthermore, the Newmark explicit algorithms are conditionally stable and still implicit for velocity. The Newmark algorithms are unconditionally stable when 2*β* ≥ *γ* ≥ 1/2. Therefore, it is impossible for the Newmark algorithms to achieve unconditional stability within the framework of an explicit displacement formulation. As *γ* and *β* are parameters independent of the structural model, the Newmark algorithms are called model-independent algorithms.

In contrast to model-independent algorithms, model-based integration algorithms are unconditionally stable and have an explicit displacement formulation. That is, model-based integration algorithms predict displacement based on equilibrium at the current time step and calculate the acceleration by satisfying equilibrium at the next time step. They can be used to solve general dynamic problems, i.e., the structural responses of linear and nonlinear SDOF and MDOF structures subjected to dynamic loads. Compared with implicit algorithms, explicit algorithms are computationally more efficient at solving nonlinear dynamic problems without iterations to obtain tangent stiffness. To solve the EOM of an MDOF system with a large number of DOFs, a very small time step is required for a conditionally stable algorithm to ensure stability, because the time step limit is inversely proportional to the highest natural frequency of the structural system. Therefore, the unconditionally stable and explicit model-based integration algorithms have very promising computational efficiency, particularly for the nonlinear MDOF dynamic problems.

Model-based integration algorithms are explicit for displacement, whereas they are not always explicit for velocity. Therefore, model-based integration algorithms can be classified as dual-explicit and semi-explicit according to their formulation of velocity. Specifically, dual-explicit means that the algorithm is explicit for both displacement and velocity, whereas semi-explicit indicates that the algorithm is explicit only for displacement, and implicit for velocity. For instance, the first model-based integration algorithm, which was developed by Chang [27], is semi-explicit, whereas the well-known CR algorithm [28] is dual-explicit. Dual-explicit algorithms are clearly more suitable for SSTT, particularly when the experimental substructure is velocity-dependent. In this study, a family of dualexplicit model-based integration algorithms [29], called GCR algorithms, is adopted. The formulation and characteristics of GCR algorithms are briefly summarized as follows.

Chen and Ricles [28] proposed the first dual-explicit model-based integration algorithm using the discrete control theory. The difference equations of the CR algorithm for the MDOF system are expressed as:

$$
\dot{\mathbf{X}}\_{i+1} = \dot{\mathbf{X}}\_{i} + \Delta t \mathbf{a}\_{1} \ddot{\mathbf{X}}\_{i}, \\
\mathbf{X}\_{i+1} = \mathbf{X}\_{i} + \Delta t \dot{\mathbf{X}}\_{i} + \Delta t^{2} \mathbf{a}\_{2} \ddot{\mathbf{X}}\_{i} \tag{3}
$$

where *α*1 and *α*2 are model-based integration parameter matrices, and are formulated as:

$$\mathfrak{a}\_1 = \mathfrak{a}\_2 = 4 \left( 4 \mathbf{M} + 2 \Delta t \mathbf{C} + \Delta t^2 \mathbf{K} \right)^{-1} \mathbf{M} \tag{4}$$

The CR algorithm has been demonstrated to be unconditionally stable for both linear elastic and nonlinear softening systems. In addition, the second-order accurate CR algorithm has no numerical damping. To obtain a family of model-based integration algorithms with more general and versatile numerical features, Fu et al. [29] recently developed GCR algorithms, which stands for generalized CR algorithms. The difference equations of the GCR algorithm are inherited from the CR algorithm, and the integration parameter matrices are updated by introducing two additional coefficients, *κ*1 and *κ*2, which are shown in Equation (5):

$$\mathfrak{a}\_1 = \left(\mathbf{M} + \kappa\_1 \Delta t \mathbf{C} + \kappa\_2 \Delta t^2 \mathbf{K}\right)^{-1} \mathbf{M}, \mathfrak{a}\_2 = (1/2 + \kappa\_1)\mathfrak{a}\_1 \tag{5}$$

It was found that GCR algorithms with *κ*1 and *κ*2 possess numerical properties identical to the classical Newmark algorithms with *γ* and *β*. The mapping relation is [*<sup>κ</sup>*1, *<sup>κ</sup>*2]=[*<sup>γ</sup>*, *β*]. The CR algorithm is a special type of GCR algorithm, with *κ*1 = 1/2, *κ*2 = 1/4; its counterpart includes Newmark algorithms with *γ* = 1/2, *β* = 1/4, which is the well-known constant average acceleration (CAA) algorithm. The members of the subfamily of GCR algorithms with *κ*1 = 1/2 have no numerical damping and are second-order accurate. GCR algorithms with 2*κ*2 ≥ *κ*1 ≥ 1/2 are unconditionally stable for linear elastic systems. GCR algorithms with *κ*1 > 1/2, *κ*2 ≥ (*<sup>κ</sup>*1 + 1/2)2/4 have numerical damping. The period elongation (PE) and equivalent damping ratios are two widely used indices of evaluating the accuracy of the integration algorithm [30]. The two accuracy indices of five GCR algorithms, i.e., [*<sup>κ</sup>*1, *<sup>κ</sup>*2]=[1/2, 1/4], [1/2, 1/2], [1/2, 1], [1, 1/2], [1, 1], are depicted in Figure 1. The abscissa in Figure 1 is Ω = *<sup>ω</sup>*Δ*t*, where *ω* is the circular frequency. Figure 1 shows that for a certain *κ*1, the PE increases with the increase in *κ*2. GCR algorithms with [*<sup>κ</sup>*1, *<sup>κ</sup>*2]=[1, 1/2] have PE close to that of [1/2, 1/4] (CR algorithm) and are minimum. Regarding the equivalent damping ratio, GCR algorithms with *κ*1 = 1/2 are zero, whereas GCR algorithms with *κ*1 = 1 are positive. Overall, the original CR algorithm has excellent accuracy and can be used as a reference for other GCR algorithms.

**Figure 1.** Numerical properties of GCR algorithms. (**a**) Period elongation; (**b**) equivalent damping ratio.

Figure 2 provides a general flowchart for solving the EOM for MDOF systems using GCR algorithms. The first step is to select appropriate integration coefficients *κ*1 and *κ*2 and the time step Δ*t*. Then, the structural model matrices **M**, **K**, and **C**, and the external force **F** are calculated. The model-based integration parameter matrices *α*1 and *α*2 are then obtained using Equation (5). Then, the initial conditions are calculated as follows:

$$
\ddot{\mathbf{X}}\_0 = (\mathbf{M})^{-1} \left( \mathbf{F}\_0 - \mathbf{C} \dot{\mathbf{X}}\_0 - \mathbf{K} \mathbf{X}\_0 \right) \tag{6}
$$

where **X**0, **X**0, and **X**0 are the initial velocity, displacement, and acceleration, respectively.

**Figure 2.** General flowchart of solving EOM for MDOF systems using GCR algorithms.

Then, the difference equations of the GCR algorithms, i.e., Equation (3), are adopted to predict the velocity **. X***i*+<sup>1</sup> and displacement **X***i*+<sup>1</sup> at the next time step. The displacement **X***i*+<sup>1</sup> can be used to determine the restoring force **R***i*+1, which is called state determination. If the structural system is linear elastic, the restoring force **R***i*+<sup>1</sup> can be directly calculated as **KX***i*+1; otherwise, it can be obtained using the FE method, e.g., the FE modeling in Section 3. The final step is to calculate the acceleration by rewriting the EOM, i.e., Equation (1):

$$
\ddot{\mathbf{X}}\_{i+1} = (\mathbf{M})^{-1} \left( \mathbf{F}\_{i+1} - \mathbf{C} \dot{\mathbf{X}}\_{i+1} - \mathbf{R}\_{i+1} \right) \tag{7}
$$

### **3. Finite Element Modeling of Frame Structure**

In this study, stiffness-based beam-column elements with fiber sections along with P-Δ effects are adopted to simulate the frame structure. The stiffness-based beam-column elements are used to acquire the first-order restoring force, and the P-Δ effects are considered to obtain the second-order restoring force. The following content summarizes the basic principles and procedures.

### *3.1. Stiffness-Based Beam-Column Element with Fiber Sections*
