3.1.1. Initialization

**.**

..

The model-based integration parameter matrices, i.e., *α*1 and *α*2, should be predetermined before using model-based integration algorithms. As shown in Equation (5), model-based integration parameter matrices are functions of model matrices of structure,

i.e., **M**, **C**, and **K**, and two controlling coefficients *κ*1 and *κ*2. Therefore, the first step is to find an appropriate combination of *κ*1 and *κ*2 based on the numerical properties of the GCR algorithms. Then, the model matrices of the structures should be determined. Among the three model matrices of structure, the Rayleigh damping assumption is widely adopted for establishing the damping matrix **C**:

$$\mathbf{C} = a\_{\rm m}\mathbf{M} + a\_{\rm k}\mathbf{K} = \frac{2\xi\omega\_{\rm i}\omega\_{\rm j}}{\omega\_{\rm i} + \omega\_{\rm j}}\mathbf{M} + \frac{2\xi}{\omega\_{\rm i} + \omega\_{\rm j}}\mathbf{K} \tag{8}$$

where *am* and *ak* are the combination coefficients of **M** and **K**, respectively; *ξ* is the damping ratio of structure; *ωi* and *<sup>ω</sup>j* are the *i*th and *j*th circular frequencies, respectively, which can be obtained by conducting modal analysis of the structure. Therefore, determination of the three model matrices is reduced to two matrices: **M** and **K**.

The formulation of the mass matrix can be classified into two types: lumped mass and consistent mass. The detailed construction process can be found in [31].

As a stiffness-based beam-column element with fiber sections is used, the initial stiffness matrix **K** of the structure should be sequentially established from different levels. There are four levels from bottom to top: fiber, section, element, and structure. Regarding the section with fibers, the stiffness matrix *ks*(*x*) is:

$$\mathbf{k}\_s(\mathbf{x}) = \begin{bmatrix} \begin{bmatrix} \sum\_{j=1}^{N\_f} E\_j A\_j & -\sum\_{j=1}^{N\_f} E\_j A\_j y\_j\\ -\sum\_{j=1}^{N\_f} E\_j A\_j y\_j & \sum\_{j=1}^{N\_f} E\_j A\_j y\_j^2 \end{bmatrix} \end{bmatrix} \tag{9}$$

where *Ej*, *Aj*, and *yj* are the elastic modulus, area, and centroid coordinate of the *j*th fiber, as shown in Figure 3; *Nf* represents the total number of fibers for a section.

**Figure 3.** Fiber section.

> Then, the element stiffness matrix in the local coordinate system can be calculated as:

$$\overline{\mathbf{k}}\_{\mathfrak{k}} = \int\_{0}^{\mathbf{L}\_{\mathfrak{k}}} B\_{d}(\mathbf{x})^{T} \mathbf{k}\_{\mathfrak{k}}(\mathbf{x}) B\_{d}(\mathbf{x}) d\mathbf{x} \tag{10}$$

where *Le* is the element length; *Bd*(*x*) is the differential of the displacement interpolation function with the expression of:

$$B\_d(\mathbf{x}) = \begin{bmatrix} -\frac{1}{L\_r} & 0 & 0 & \frac{1}{L\_r} & 0 & 0\\ 0 & \frac{12\mathbf{x}}{L\_r^3} - \frac{6}{L\_r^2} & \frac{6\mathbf{x}}{L\_r^2} - \frac{4}{L\_r} & 0 & \frac{6}{L\_r^2} - \frac{12\mathbf{x}}{L\_r^3} & \frac{6\mathbf{x}}{L\_r^2} - \frac{2}{L\_r} \end{bmatrix} \tag{11}$$

It is impractical to integrate Equation (10) directly, so the Gauss–Legendre integration method is typically used to indirectly solve Equation (10). In this study, the Gauss–Legendre

integration method with five integration points is used; then, Equation (10) can be rewritten as:

$$\overline{\mathcal{K}}\_{\mathfrak{c}} = \sum\_{k=1}^{5} B\_{d}(\mathbf{x}\_{k})^{T} k\_{\mathfrak{s}}(\mathbf{x}\_{k}) B\_{d}(\mathbf{x}\_{k}) \omega\_{k} \mathbf{L}\_{\mathfrak{c}} \tag{12}$$

where *xk* and *ωk* are the coordinate and weight at the *k*th integration point. The coordinates and weights of five integration points are listed in Table 1.

**Table 1.** Coordinates and weights of five integration points of the Gauss–Legendre integration method.


As the initial stiffness is elastic, direct use of the element stiffness matrix of the elastic beam-column element without fiber sections is also applicable and simpler (Equation (13)). The fiber section is adopted to formulate the element stiffness in order to make it consistent with the state determination in Section 3.1.2. It should be noted that the element stiffness matrix shown in Equation (13) is a symmetric matrix:

$$
\overline{\mathbf{k}}\_{\varepsilon} = \begin{bmatrix}
EA/L\_{\varepsilon} & 0 & 0 & -EA/L\_{\varepsilon} & 0 & 0 \\
0 & 12EI/L\_{\varepsilon}^{3} & 6EI/L\_{\varepsilon}^{2} & 0 & -12EI/L\_{\varepsilon}^{3} & 6EI/L\_{\varepsilon}^{2} \\
0 & 6EI/L\_{\varepsilon}^{2} & 4EI/L\_{\varepsilon} & 0 & -6EI/L\_{\varepsilon}^{2} & 2EI/L\_{\varepsilon} \\
0 & -12EI/L\_{\varepsilon}^{3} & -6EI/L\_{\varepsilon}^{2} & 0 & 12EI/L\_{\varepsilon}^{3} & -6EI/L\_{\varepsilon}^{2} \\
0 & 6EI/L\_{\varepsilon}^{2} & 2EI/L\_{\varepsilon} & 0 & -6EI/L\_{\varepsilon}^{2} & 4EI/L\_{\varepsilon}
\end{bmatrix}
\tag{13}
$$

where *E* is the elastic modulus of the material; *A* and *I* are the area and inertia moment of the cross-section, respectively. It should be noted that only flexural deformation (without shear deformation) is considered in Equation (13).

Then, the element stiffness matrix *ke* in the global coordinate system can be transformed from the element stiffness matrix *ke* in the local coordinate system:

$$\mathbf{k}\_{\mathbf{c}} = T\_{\mathbf{c}}^{\mathrm{T}} \overline{\mathbf{k}}\_{\mathbf{c}} T\_{\mathbf{c}} \tag{14}$$

where *Te* is the coordinate transformation matrix. After calculating *ke* for all elements, the initial stiffness **K** of the structure can finally be assembled by mapping the degrees of freedom (DOFs).
