Matrix equations models for nonlinear dynamic analysis of two-dimensional and three-dimensional RC structures with lateral load resisting cantilever elements

In control engineering and structural dynamics, mathematical models such as the state-space representation, equation of motion, and the phase plane are matrix equations describing the system equilibrium. This paper develops novel matrix equations models for linear/nonlinear dynamic analysis of reinforced concrete (RC) buildings with cantilever elements lateral load resisting systems (e.g., RC shear wall, RC core). The models offer a new approach for introducing two-dimensional and three-dimensional cantilever structures to control the theory’s state-space representation and structural dynamics’ equation of motion. The development primarily addresses the stiffness and mass matrices. The proposed displacement-related stiffness matrix of cantilever elements satisfies the necessary conditions of symmetricity and elemental boundary conditions. The nonlinear matrix structural analysis employs a smooth hysteretic model for deteriorating inelastic structures, referring to the relation between the bending moment and the bending curvature through the bending stiffness. The parameters controlling the cyclic behavior regard a composite RC cross section subject to gravitational load and bending simultaneously. The paper includes four examples that exemplify the practical utilization of the matrix equations models in analyzing two-dimensional and three-dimensional structures of linearly elastic and inelastic properties. The four examples demonstrated the idealized applicability of the matrix equations models for modal analysis, pushover analysis, and inelastic earthquake response analysis.

Abstract In control engineering and structural dynamics, mathematical models such as the statespace representation, equation of motion, and the phase plane are matrix equations describing the system equilibrium. This paper develops novel matrix equations models for linear/nonlinear dynamic analysis of reinforced concrete (RC) buildings with cantilever elements lateral load resisting systems (e.g., RC shear wall, RC core). The models offer a new approach for introducing two-dimensional and three-dimensional cantilever structures to control the theory's state-space representation and structural dynamics' equation of motion. The development primarily addresses the stiffness and mass matrices. The proposed displacement-related stiffness matrix of cantilever elements satisfies the necessary conditions of symmetricity and elemental boundary conditions. The nonlinear matrix structural analysis employs a smooth hysteretic model for deteriorating inelastic structures, referring to the relation between the bending moment and the bending curvature through the bending stiffness. The parameters controlling the cyclic behavior regard a composite RC cross section subject to gravitational load and bending simultaneously. The paper includes four examples that exemplify the practical utilization of the matrix equations models in analyzing two-dimensional and three-dimensional structures of linearly elastic and inelastic properties. The four examples demonstrated the idealized applicability of the matrix equations models for modal analysis, pushover analysis, and inelastic earthquake response analysis. Distance between the center of the rth ceiling center-of-mass and the DOFs origin in the x-direction x dis l;n Distance between the center of the 'th bending wall and the DOFs origin in the x-direction y CM r Distance between the center of the rth ceiling center-of-mass and the DOFs origin in the y-direction y dis l;n Distance between the ceiling centerof-mass and the DOFs origin in the Pinching length factor R j z; M R ; u À Á Stiffness degradation factor ROI r Rth ceiling radius-of-inertia about the DOFs origin ROI Dz;l N ? 1 9 N ? 1 Diagonal matrix whose components are the 'th element's radius-of-inertia at the mesh lines ROI l;n N = z/Dz ? 1 Component of ROI Dz;l T f !M Integration transformation matrix from lateral force into bending moment   Various nonlinear dynamic analysis methods have been developed to evaluate the response of buildings under cyclic loading. Wong and Wang [1] modified the classical force analogy method for momentresisting frame models by introducing static condensation, which reduced the degrees of freedom (DOFs) and resolved in a veritable mass matrix-essential to the force analogy method. Li and Li [2] utilized the force analogy method to analyze the dynamic response of moment-resisting frame structures with energy dissipation devices. Sun et al. [3] proposed their numerical substructure method for seismic analysis and determining the locally yielded lateral load resisting elements. Chang et al. [4] presented the simplified generalized conforming mixed quadrilateral for membrane elements. Pozo et al. [5] suggested using the Bouc-Wen model to determine the position and acceleration of dynamic systems under cyclic loads. Liu and Yu [6] addressed the previously studied quasi-zero stiffness vibration isolator and enhanced the dynamic model by referring to the damping effects. This paper aims to develop matrix equations models to analyze linear/nonlinear RC cantilever structures. The models are ideal for control theory's state-space representation, and the structural matrices apply to numerical evaluation methods that regard the equation of motion.
Introducing the matrix models to equation-ofmotion numerical evaluation methods yields an alternative to finite element simulations that may require extensive computation time and effort. This paper adds to existing alternatives, such as the knowledgeenhanced deep learning algorithm of Wang and Wu [7], which refers to physical equations and semiempirical formulas. Nevertheless, while the knowledgeenhanced deep learning algorithm is efficient in computation, it requires knowledge of neural network theory. Another difference is that the nonlinear response estimation in [7] employs the long-term memory approach, while the matrix equations models utilize cyclic mechanical behavior models.
The mechanical understanding of structural elements under cyclic loading has significantly improved in recent decades. The paper of Fragiadakis and Papadrakakis [8] alone overviews dozens of finite element and computation methods for addressing the inelastic response of solely beam-column members. The review paper by Younesian et al. [9] described methods developed in recent decades for analyzing nonlinear dynamic systems of various model types (e.g., discrete, multilayered, discontinuous). Mata et al. [10] showcase dozens of numerical approaches developed during a 20 years gap for evaluating the response of RC structures with energy dissipation devices. Very recently, Castaldo et al. [11] and Gino et al. [12] addressed the uncertainties of epistemic nature in cyclically loaded reinforced concrete systems and proposed using a partial safety factor corresponding to the resistance model uncertainties in the use of nonlinear finite element analyses. Ismail et al. [13] presented a broad survey on the development and utilization of the Bouc-Wen-based models for the nonlinear analysis of structural elements, mechanical systems, soil behavior, and more. The paper exemplifies the extent to which the Bouc-Wenbased models are popular and suitable for addressing structural elements' nonlinearities and quantifying their yielding energy under cyclic loading.
Smooth hysteretic models, that stem from the Bouc-Wen equation, have been developed to provide a continuous change of stiffness or functions that address the bending moment-curvature relationships (e.g., Sivaselvan and Reinhorn [14]; Wang et al. [15]; Charalampakis and Koumousis [16]; Charalampakis [17]). The smooth hysteretic models are utilized to analyze the cyclic response of a wide variety of mechanical systems such as steel framed-tubes [18], shape-memory-alloy bars [19], wood joints with yielding plates/nails/bolts [20], bolted joint's dynamic frictional contact [21], connection fracture in low-rise steel buildings [22], base-isolation system [23], dampers that interconnect two adjacent structures [24], and more. Smooth hysteretic models have also been incorporated in approximating seismic demands for buildings using cyclic pushover analysis [25] and evaluating their collapse fragility curve [26],27). This paper follows suit and employs the smooth hysteretic model of Sivaselvan and Reinhorn [14] for modeling the cyclic mechanical behavior of deteriorating inelastic RC cantilever elements, which function as the lateral load resisting system of the building.

Two-dimensional element
This section develops the matrix equations model for the two-dimensional cantilever element. The partial differential equation (PDE) governing the cantilever element's dynamic equilibrium (with boundary conditions) is described. Then, the displacement-related stiffness matrix is defined, and the model's numerical accuracy is examined.

Two-dimensional model and assumptions
A general elevation scheme of a slender cantilever element is depicted in Fig. 1. The element is characterized by distributed mass per unit length m z ð Þ, distributed elasticity EI z ð Þ, and is subjected to a distributed lateral load function pðz; tÞ. It is assumed the element bends about the strong/weak axes and, accordingly, this paper's models refer to either of the cross sections shown in Fig. 1 or similar types where the strong/weak axes merge with the primary axes.
The partial differential equation (PDE) governing the lateral motion of the undamped system is: where f R z; t ð Þ is the resisting force at time t at coordinate z, u w z; t ð Þis the lateral displacement at time t at coordinate z, and its second time derivative € u w z; t ð Þ is the horizontal acceleration at time t at coordinate z, the term M R z; t ð Þ identifies the bending moment at time t at coordinate z, and V R z; t ð Þ is the shear force at time t at coordinate z.
The Euler-Bernoulli assumption is adopted regarding the lateral movement of the structure. That means the cross section movement is constrained to its horizontal plane and its dimensions, B z ð Þ and L z ð Þ, remain stationary. Accordingly, the angular accelerations and shear deformations are assumed to be neglected, which leads to the following relationships: and the bending relationship: where u z; t ð Þ is the bending curvature. The term of f R z; t ð Þ in Eq. (1) is determined by substituting Eqs. (4) into (3), which yields: Using the separated variables technique for solving Eq. (1), we have u w z; t ð Þ ¼ / z ð Þq t ð Þ. The function / z ð Þ caters to the homogeneous boundary conditions (BC), and q v t ð Þ governs the initial-time conditions and dynamic response. The solution procedure starts with addressing / z ð Þ, which leads to modal superposition due to multiple solutions (see chapter 16 in [28]). That is: Fig. 1 Cantilever element of distributed mass and elasticity Substituting Eqs. (5) and (6) into Eq. (1) gives: where f R/ v z ð Þ is defined herein as the vth modal resisting force spatial function.
The analytical solution of a cantilever element typically involves defining the shape function that , determining the shape frequencies b v¼1;2;...;1 that satisfy the boundary conditions, calculating the shape function's coefficients, deriving the modal angular frequencies x v¼1;2;...;1 , and evaluating q v t ð Þ. This paper, however, takes an alternative approach and looks to define the mass and stiffness matrices of the structure. These matrices (i.e., matrix equations model) can be employed in numerical evaluation methods that regard the equation of motion.

Elemental stiffness matrix
The mass and stiffness matrices are the fundamental structural matrices since they dictate the structure's modal properties and define the inherent damping matrix (e.g., Rayleigh classical damping, Caughey classical damping). Since the PDE's inertia force relates to horizontal accelerations (i.e., the second derivative of the lateral displacements), the mass matrix is composed of the mass per unit length. The resisting force, however, relates to the bending curvature and, thus, the stiffness matrix has to apply transformation between curvature coordinates and displacement coordinates while catering to the boundary conditions. The cantilever element's stiffness matrix is equivalent to the stiffness distribution operator k z ð Þ, which provides the relationship between the vth mode shape function / v z ð Þ and the vth modal resisting force f R/ v z ð Þ: That implies the stiffness distribution operator is: Note that the modal terms in Eq. (8) are used to express the linear relation between lateral displacements and resisting force through k z ð Þ. The developed matrix equations model does not employ modal superposition analysis.
Equation (9) suggests that the stiffness matrix applies four numerical differential transformations. In addition, the stiffness matrix must be symmetric so that the modal frequencies are real numbers (to suit undamped systems) and must cater to the elemental boundaries specified in Eq. (7). Developing the stiffness matrix becomes more applicable by implementing numerical integration rather than numerical differential. As shown herein, numerical integration can address all the necessary matrix properties rather than numerical differentiation. Define the flexibility distribution operator sF z ð Þ-the inverse of the stiffness distribution operator. That is: Accordingly, the operator F z ð Þ consists of four integral operators in z such that: The vertical coordinate variable z is bounded by the structure's dimensions z 2 0 H ½ -setting the integration limits. The start point of the integration (i.e., 0 or H) plays a crucial role in the stiffness matrix development by obtaining a symmetric matrix form. Let us observe the first integration operator applied to f R/ v z ð Þand examine the two options for integration limits: option one : Option two is preferred because EI H ð Þ/ 3 ð Þ v H ð Þ ¼ 0 and, thus, it provides the distributed modal shear force V R v z ð Þ without a tail. That concept is repeated in the other three integration operators as follows: Accordingly, the flexibility distribution operator F z ð Þ is: The operator F z ð Þ is now addressed in matrix form. Figure 2 depicts a discretized vertical cantilever element in z. The lateral DOFs are uniformly distributed with mesh size Dz ¼ H=N and indexed as n ¼ 0; 1; . . .; N, where n ¼ 0 refers to the ground level and n ¼ N refers to the top level.
Referring to the mesh depicted in Fig. 2, Eq. (14) is represented in matrix terms as: where F Dz is the symmetric flexibility transformation matrix, / Dz;v is the vth mode shape vector and is the approximation of / v z ð Þ at the mesh lines z ¼ 0; Dz; . . .; N=Dz: and f R/ Dz;v is the vth modal resisting force vector and the approximation of f R/ v z ð Þ at the mesh lines: It is noted that the matrix/vector terms notated with lower case Dz refer to the discretized cantilever element's mesh lines throughout this paper. Define F Dz by the following matrix transformation: is an integration transformation matrix from ¼ 0 to ¼ z, and EI Dz is the N ? 1 9 N ? 1 diagonal elasticity matrix defined as: The left Riemann sums (LRS) numerical integration method is employed for defining T The idealized scheme of the LRS method is employed because it produces a symmetric F Dz due to its idealized form. The LRS integration transformation matrix is given by: Substituting Eqs. (20a) and (20b) into Eq. (18) yields F Dz of the general form: where F Dz is a symmetric matrix referred to herein as the N-1 9 N-1 flexibility sub-matrix, and 0 is a matrix of zeros.
The inverse form of F Dz gives the N ? 1 9 N ? 1stiffness matrix k Dz , but F Dz ð Þ À1 is singular due to the zeros on its diagonal. Therefore, the following equivalent form is proposed: Here, k Dz contains the inverse of the invertible part of the flexibility matrix F Dz ð Þ À1 , and K is a 2 9 2 diagonal matrix whose components K 1 and K 2 are large-valued numbers. The matrix K resembles inverse flexibility matrix singularity. The components K 1 and K 2 correspond to two irrelevant modal systems of mode shapes Due to the large-valued numbers, the associated modal frequencies are unreasonably large-preventing these modes from practically participating in the analysis (similar to very stiff modal systems of zero displacements).
Referring to Eqs. (15)- (22), the cantilever element's dynamic equilibrium of Eq. (7) is represented by the following matrix equations model: where m Dz is the N ? 1 9 N ? 1 mass matrix of the cantilever element and p Dz t ð Þ is the N ? 1-dimensional lateral load vector. Both refer to the mass and loads applied at the mesh lines as follows: In addition, respectively. That is: Lastly, q Dz t ð Þ is an N ? 1-dimensional vector containing the v = 1,…, N ? 1 time governing functions: The modal properties of the cantilever element are evaluated using m Dz and k Dz by the eigenvalues and eigenvectors of m Dz ð Þ À1 k Dz . The matrix equations model numerical accuracy is thus examined in terms of modal properties while referring to the exact mathematical solution.

Mesh-size efficiency
The following subsection exemplifies the relation between the mesh size and the numerical accuracy of the matrix equations model. The examination refers to a cantilever element of uniform mass per unit length mðzÞ ¼ m and constant elasticity function EIðzÞ ¼ EI. The PDE governing undamped free vibrations is then: Using the analytical solution provided in chapter 16 in [28], we have the following formulas for x v and / v ðzÞ: where C v is an arbitrary coefficient, and the quantity of e v H stems from: The coefficient C v is taken so that the Euclidean norm of k/ v 0 ð Þ; / v Dz ð Þ; . . .; / v ðHÞk equals 1.0. The modal frequencies and mode shapes of the matrix equations model stem from: where I is the N ? 1 9 N ? 1 identity matrix. For a certain Dz, the mass and stiffness matrices are: are defined in Eq. (20). The relative numerical error between the analytical solution's modal properties and the matrix equations model's modal properties is quantified as:  Figure 3 shows that for Dz=H\0:01 the relative error decreases linearly with Dz=H while Fig. 4 shows the same for Dz=H\0:002. It can also be  The relative error of the modal roof displacement approximation decided that for errðÞ\0:01 a mesh size of Dz=H\0:01 is recommended. That concludes the comparing the matrix equations model and the analytical solution. It is important to note that the comparison also reflects the numerical accuracy of three-dimensional structures and inelastic elements. Practically speaking, in the case of threedimensional structures, each cantilever element is solved individually and integrated into the global structure. Hence, the global numerical error magnitude is of an individual cantilever multiplied by the number of elements. In the case of inelastic nonlinear analysis, the elasticity function EIðzÞ shifts along the mesh lines z ¼ 0; Dz; . . .; N=Dz between elastic-loading/plastic-loading/unloading states during simulation. For an individual time step, the numerical error is as exemplified above for a two-dimensional structure. Regarding the stability of the numerical simulation, an unconditionally stable method should be employed-as shown in Example 4.

Three-dimensional structure
The three-dimensional structure's matrix equations model employs the cantilever element's mass and stiffness matrices from the previous section. The three-dimensional model refers to the Euler-Bernoulli assumption so that the ceilings are considered rigid diaphragms of zero mutuality between the horizontal x and y directions, and the elements have no movement in the z-direction. The cantilever element is now considered a part of the lateral load resisting system and is notated henceforth as the 'th element (within the three-dimensional structure). Figure 5 shows the cantilever element scheme with the new parametric notations, and Fig. 6 illustrates the three-dimensional structure. In Fig. 6, r is the ceilings' index so that r ¼ 1; . . .; R. Also, H r denotes the height of the rth ceiling above ground level so that H R is the roof level. For mesh size Dz ¼ H R =N, the DOFs, index is n ¼ 0; 1; . . .; N. A cantilever element can have less than N-DOFs when its height is H l \H R and, in that case, that element is of N l ¼ H l =Dz DOFs as depicted in Fig. 7.
The lth cantilever element's matrix equations model from the previous section is now given the following notations: Each of the l ¼ 1; 2; . . .; l cantilever elements within the three-dimensional structure has its matrix equations model, and all are added to the global dynamic equilibrium equation.
The equation of damped motion governs the dynamic equilibrium. In the case of a three-dimensional structure subjected to lateral and torsional loading, the equation-of-damped-motion formulation is notated as: The DOFs e u x;n t ð Þ and e u y;n t ð Þ denote lateral deformations at the z ¼ nDz level in the x-direction and ydirection, respectively, and e u h;n t ð Þ is the rotation about the DOF origin. The equation-of-damped-motion formulation consents with the global DOF origin to be located anyplace within the x-y plane and not necessarily at the floor plan's center-of-mass (CM). The 3N ? 3-dimensional load vector e p Dz t ð Þ is composed of three N ? 1-dimensional sub-vectors e p Dz;x t ð Þ, e p Dz;y t ð Þ, and e p Dz;h t ð Þ applied in the directions of e u Dz;x t ð Þ, e u Dz;y t ð Þ, and e u Dz;h t ð Þ, respectively: The global mass matrix of the three-dimensional structure is of the diagonal form: where m Ã Dz is the N ? 1 9 N ? 1 diagonal mass matrix of the structure, and I mÃ Dz is the N ? 1 9 N ? 1 diagonal moment-of-inertia matrix. The nth diagonal component of m Ã Dz and I mÃ Dz refers to the z ¼ nDz level. The global mass matrix e m Dz is formed using transformation matrices a m1;l for the respected 'th element and a m2;r for the rth ceiling, so that: In Eq. (42), the zeros matrices (0 N l ÂNÀN l , 0 NÀN l ÂN l , and 0 NÀN l ÂNÀN l ) treat the case of H l \H R , m r is the rth ceiling mass, a m1;l is a 3N ? 3 9 3N ? 3 transformation matrix defined according to Eq. (43), and a m2;l is a 3N ? 3 9 3 transformation matrix defined according to Eq. (44). Furthermore, in Eq. (43), 0 is a N ? 1 9 N ? 1 matrix of zeros, ROI Dz;l is a diagonal matrix whose components are the 'th element's radius-of-inertia at the mesh lines: ROI Dz;l ¼ diag ROI l;0 ; . . .; ROI l;n ; . . .; ROI l;N È É so that: The parameters x dis l;n and y dis l;n are the horizontal distances between the center of the 'th element and the DOFs origin-as depicted in Fig. 8. In Eq. (44), e 0 is a N ? 1-dimensional vector of zeros, e r denotes a vector whose elements are all zero except for the r ¼ H r =Dz þ 1 component whose value is ''1.0,'' and ROI r is the rth ceiling radius-of-inertia about the DOFs origin, which is given by: where ROI CM r is the rth ceiling radius-of-inertia about its CM, and the parameters x CM r and y CM r are the horizontal distances between the ceiling's CM and the DOFs origin-as depicted in Fig. 9.  The global stiffness matrix e k Dz is the combination of the element stiffness matrices k Dz;l¼1;2;... transformed into global displacement coordinates using the direct stiffness method (DSM): In Eq. (48), the zeros matrices (0 N l ÂNÀN l , 0 NÀN l ÂN l , and 0 NÀN l ÂNÀN l ) treat the case of H l \H R . The transformation matrix a k;l of Eq. (49a) transforms an element whose strong\primary moment of inertia is in the x-z plane, and q xyl;n is the stiffness mutuality from the x-direction to the y-direction at the nth DOF level. The transformation matrix a k;l of Eq. (49b) transforms an element whose strong moment of inertia is in the y-z plane, and q yxl;n is the stiffness mutuality from the y-direction to the xdirection at the nth DOF level. If an element has significant moment-of-inertia quantities in its primary perpendicular directions, it is transformed twice for each direction.
The 3N ? 3 9 3N ? 3 global inherent damping matrix e c Dz is determined according to e m Dz and e k Dz using Caughey classical damping: Equations (38)-(50) conclude the three-dimensional structure's matrix equations model. The matrix equations models of the two-dimensional cantilever element and three-dimensional structure are now enhanced with the material's inelastic mechanical properties.

Smooth hysteretic model
The smooth hysteretic model for deteriorating inelastic structures in Sivaselvan and Reinhorn [14] is employed to implement the mechanical behavior of a cross section subjected to gravitational load and timevarying bending moment. The smooth hysteretic model defines the relation between the first-order time derivatives of the bending moment and the bending curvature through the bending stiffness: where j z; u; _ u ð Þis the inelastic bending stiffness. The bending stiffness considers nonsymmetrical yielding between the positive and negative yield bending moments, stiffness degradation, strength degradation, and the possibility of a pinching effect.
In the case where the pinching effect is not regarded, j z; u; _ u ð Þis a combination of the permanent portion of the initial bending stiffness and the hysteretic bending-stiffness portion: where j 0 z ð Þ denote the initial bending stiffness (without deterioration), and since shear deformations are neglected, due to the introduction of the Euler-Bernoulli assumption, it is equal to the elasticity distribution: Also, in Eq. (54), j Ã z; M RÃ ; _ u À Á is the hysteretic bending-stiffness portion, and a is the ratio of postyielding to initial stiffness. The function M RÃ ðz; tÞ is the hysteretic portion of m R (z,t) and is given by In the case where the pinching effect is regarded, the total bending stiffness is calculated as follows: where j P M RÃ ; £ À Á is the pinching stiffness and is connected in ''series'' to j Ã z; M RÃ ; _ u À Á . The equation of j Ã z; M RÃ ; _ u À Á is formulated so that the hysteretic portion is related to the ratio between the bending moment and the yielding moment: In Eq. (58), m is the parameter controlling the smoothness of the transition from elastic to plastic region, g 1 and g 2 ¼ 1 À g 1 are parameters controlling the shape of the unloading curve, R j z; M R ; u À Á is the stiffness degradation factor in z whose behavior is adjusted by the degradation parameter a as follows: Also, in Eq. (58), M yldÃ z ð Þ is the hysteretic portion of the yielding moment M yld z ð Þ and is given by: In Eqs. (59) and (60), in a case where M yld z ð Þ is nonsymmetric, it is defined as: in which M yld;þ z ð Þ and M yld;À z ð Þ denote the positive/negative yielding bending moments. Equation (61) is added with a correction suggested by Wang et al. [15]. The yielding bending-moments' degraded quantity (i.e., degraded strength) is calculated using the following continuous terms: where M y0;þ ðzÞ/M y0;À ðzÞ are the initial positive/ negative yield bending moments, u þ ðz;tÞ/u À ðz;tÞ are the positive/negative bending curvatures, u ult;þ ðzÞ/ u ult;À ðzÞ are the ultimate positive/negative bending curvatures, b 1 is the ductility-based strength degradation parameter, E yld ðz;tÞ is the energy dissipated by yielding, b 2 is the energy-based strength degradation parameter, and E ult ðzÞ is the ultimate dissipated energy by yielding, which is calculated by loading the cross section monotonically to reach the ultimate curvature without any degradation: The pinching stiffness j P M RÃ ;u À Á is usually a result of cracks closure and is added in series to the hysteretic bending stiffness to model this effect. The pinching stiffness is modeled as the inverse of the probability density function of the normal distribution with an expected value M RÃ z ð Þ¼kM yldÃ z ð Þ (i.e., the mean moment level about which pinching occurs) and with variance M rÃ z ð Þ¼rM yldÃ z ð Þ (i.e., the measure of the moment range over which slip occurs): Here, s is the pinching length and is defined according to the pinching length factor R s as follows: That sums up the smooth hysteretic model. Equations (53)-(66) define the inelastic bending-stiffness behavior and are now implemented in the twodimensional cantilever element's matrix equations modeling.

Inelastic matrix equations model
The inelastic analysis follows the smooth hysteretic model and addresses the cantilever element in bending-curvature coordinates. The bending curvature relates to the bending moment. Accordingly, integrating both sides of Eq. (1) by R z H R z H dd yields the following bending-moment equilibrium: where M I z; t ð Þ is the bending-moment function applied by the lateral inertia force m z ð Þ € u w z; t ð Þ, and M p z; t ð Þ is the bending-moment function applied by the lateral load p z; t ð Þ. In the boundary conditions, V I z; t ð Þ is the distributed shear force by the lateral inertia force and V p z; t ð Þ is the distributed shear force by the lateral load.
Equation (68) represents Eq. (67) in matrix equations form by referring to the bending-moments' equilibrium at the n ¼ 0; 1; . . .; N lateral DOF levels. BC: The vector M p Dz t ð Þ is determined using the LRS integration transformation matrix: and the bending-stiffness matrix is of the diagonal form: The smooth hysteretic model for deteriorating inelastic structures presented in 4.1 determines the components of the j Dz u Dz ; Þ_ u Dz t ð Þ, the bending-moment vector M R Dz t ð Þ stems from an implicit term. Consequently, the numeric simulation should utilize an iterative method for implicit numercal schemes (see chapter 5 in [28]).
The displacement vector u w Dz t ð Þ is determined using the LRS integration transformation matrix. For determining the modal properties of the structure, one requires the N ? 1 9 N ? 1 bending-stiffness matrix j Dz u Dz ; _ u Dz ð Þand the N ? 1 9 N ? 1 bending-curvature-related-mass matrix, denoted as Dz , which is given by: That ends the matrix equations model formulation of the inelastic cantilever element. Next, the inelastic three-dimensional structure model is developed.

Inelastic three-dimensional structure
The equation of damped motion of the inelastic threedimensional structure, subjected to lateral and torsional loadings, is written as: where e m Dz , e c Dz , e p Dz t ð Þ, _ e u Dz t ð Þ, and € e u Dz t ð Þ are given by the matrix equations model in 3 and e f H Dz is the global hysteretic resisting force vector. Since this paper implements the smooth hysteretic model in 4.1, the three-dimensional structure is analyzed in bending-coordinates. Accordingly, Eq. (73) is modified using the transformation matrix from lateral force into bending moment T f!M , based on the the LRS integration transformation matrix, which is defined as: wheree Dz is the global bending-curvature-related-mass matrix and e u Dz t ð Þ is the 3N ? 3-dimensional global bending-curvature vector composed of the sub-vectors       the global bending-curvature-related-inherent-damping matrix is defined using Caughey damping: In Eq. (79), the global mode shape functions and modal frequencies e / Dz;v and e x Dz;v stem from: Performing matrix structural analysis using the above equations assesses the global bending curvature and its time derivatives, e u Dz t ð Þ, _ e u Dz t ð Þ, and € e u Dz t ð Þ. The global displacement vector, and its time derivatives are calculated using the transformation matrix T f!M : Also, the elemental bending-curvature vector, lateral deformation vector, and their time derivatives are determined using the transformations: That concludes the matrix equations model of the inelastic three-dimensional structure and the models' development. The first example employs the matrix equations model developed in 2.2 for analyzing a two-dimensional linearly elastic model and determining its modal properties. In example 2, the cantilever wall is placed in a three-dimensional structure and studied using the equations model in 3. Examples 3 and 4 define and introduce the two-dimensional and three-dimensional structures' inelastic mechanical properties. In Example 3, pushover analysis is applied to the cantilever element, and in Example 4, earthquake response analysis is performed for the three-dimensional structure.
6.1 Example 1: two-dimensional model The first example calculates the cantilever wall's mass and stiffness matrices and illustrates the first eight mode shapes. The number of lateral DOF levels is set as N ¼ H l =Dz ¼ 2000 and, thus, the size of all matrices provided in the current example is 2001 9 2001. A stage-by-stage description for determining the mass and stiffness matrix is provided herein. Stage-b. Calculate the mass matrix m Dz;l -using Eq. (24) Given that m l z ð Þ ¼ 3:375 ton=m, the mass matrix is:
Stage-d. Construct the stiffness matrix k Dz;l -using Eqs. (22) The matrices F Dz;l À1 and K are used to construct k Dz;l . The components of K are set as: and the stiffness matrix is: That concludes the mass and stiffness matrices calculation.
The modal frequencies and mode shapes are approximated (see Eqs. 33 and 34) and compared to those of the analytical solution (see Eqs. 30-32). Table 1 shows the relative error between the modal frequencies, and Figs. 11 and 12 illustrate the absolute error between the first eight mode shapes in z 2 0 H l ½ . It is shown that the more significant errors are located where the mode shape displacement is nearing zero (i.e., / v z ð Þ % 0). Nevertheless, the absolute error does not exceed the value of 4 9 10 -3 .

Example 2: three-dimensional model
The RC wall in Example 1 is used to carry the gravitational and lateral loads of a 20-story building. Figures 13 and 14 depict its floor plan and elevation scheme, respectively. The RC walls are located asymmetrically about the global DOFs origin, chosen to be positioned at the center of the circumscribed square. The ceiling thickness is 20 cm, and the story heights are 3.0 m (i.e., H r¼1 ¼ 3:0m, H r¼2 ¼ 6:0m,…,H R¼20 ¼ 60:0m). Considering the structure's geometry and the cantilever wall's stiffness and mass matrices from the previous example, the current example first calculates the mass and stiffness matrices of the three-dimensional structure according to the stage-by-stage description below.
Stage-a. Determine the transformation matrices a m1;l for the ' = 1,2,…,8 elements-using Eq. (43) The matrix a m1;l transforms the 'th cantilever wall mass matrix into the global coordinates and adds its contribution to the global mass matrix e m Dz . The elements' distance coordinates about the DOF origin are:   Using the distance coordinates in Eqs. (45) and (46) gives the cantilever walls' radius-of-inertia matrices: and, closing Stage-a, we determine a m1;l of 6003 9 6003 size for all the RC walls: The matrix a m2;r is used to transform the rth ceiling mass matrix into the global coordinates and add it to the global mass matrix. The distances between the ceilings' CM and the DOFs origin are: . . .; 20 The ceilings' radius-of-inertia about the DOFs origin is: where e r is a vector whose components are all zero except for the H r =0:03 þ 1 component whose value is ''1.0'' (H r is the ceiling height above ground).
Stage-c. Determine the DSM transformation matrices a k;l for ' = 1,…,8 elements-using Eq. (49) The matrix a k;l transforms the 'th cantilever wall stiffness matrix into the global coordinates and adds its contribution to the global stiffness matrix e k Dz . The employed RC wall is slender (i.e., L l z ð Þ ) B l z ð Þ). Therefore, the stiffness mutuality between the lateral directions can be neglected and q xy;n and q yx;n are set as: elements in x À z plane : q xy;n ¼ 0 elements in y À z plane : q yx;n ¼ 0 The DSM transformation matrices are: Calculate e m Dz using the transformation matrices a m1;l and a m2;r -using Eq. (42) The ceiling mass is m r ¼ 196:91ton when considering its self-load and the assumed amount of additional load 0:5ton=m 2 . Given a m1;l , a m2;r , m Dz;l , and m r , the global mass matrix is calculated using the term: Stage-e. Calculate e k Dz using the transformation matrices a k;l -using Eq. (48) The global stiffness matrix is calculated using the DSM transformation: That concludes example 2's stage-by-stage description. The three-dimensional structure's modal frequencies and mode shapes are determined using Eqs. (51) and (52). The 16 highest modal periods are specified in Table 2, and their respected mode shapes are depicted in Fig. 15 (modes 1-8) and Fig. 16 (modes 9-16). The  18 Axial force bending-moment interaction diagrams-strong direction Fig. 19 Axial force bending-moment interaction diagrams-weak direction modal systems' seismic contributions are also quantified and specified in Table 2. The seismic contribution refers to the participation in base shear using modal static response subjected to earthquake load distribution with the additional eccentricity of 5% of the floor-dimension perpendicular to the direction of the seismic action (i.e., AE0:75m in x and y directions). The participation factors in Table 2 indicate that the modes that bend in the x and y directions are of neglectable contribution, while the torsional modes are more effective. That indication feat the building's asymmetric plan, which is of minor lateral resistance in the diagonal direction. The stress-strain relationship between the concrete and the reinforcing steel is taken according to EC 2 [29]. The three types of RC cross section designs, depicted in Fig. 17, are employed. Each design is provided with a different amount of steel reinforcement distribution to address the stories' varying gravitational load levels. The elastic/plastic axial force bending-moment interaction diagrams for the cross section designs are depicted in Figs. 18 and 19. Figure 18 refers to bending in the longitudinal (strong) direction, and Fig. 19 refers to bending in the transverse (weak) direction.
The smooth hysteretic model's parameters stem from the axial force bending-moment interaction diagrams. The gravitational load (including the ceiling's load and wall's self-weight) is applied to all stories' walls as specified in Table 3. For each gravitational load quantity, the initial positive/negative yield bending moments (M y0;þ ðzÞ and M y0;À ðzÞ) and the ultimate positive/negative bending moments (M ult;þ ðzÞ and M ult;À ðzÞ) are determined, and the consequent positive/negative ductility demand (u ult;þ and u ult;À ) is defined. Table 4 shows the calculated values in the strong direction, and Table 5 shows the calculated values in the weak direction. The other smooth hysteretic model parameters are defined as g 1 ¼ g 2 ¼ 0:5 (unloading stiffness equals to aEI), m ¼ 10 (fast transition from elastic to the plastic range), a ¼ 200 (stiffness deterioration is neglected), and R s ¼ 1 (no pinching effect is regarded). Also, according to Table 6.5 in ASCE/SEI 41-06 [30], the post-yielding ratio to initial stiffness is considered a = 0.5. Incremental load analysis is applied to RC wall ' = 4. The wall's mass and stiffness matrices in bending-coordinates are determined using the stageby-stage description below.
Stage-a. Define the integration matrices T Stage-b. Calculate the bending-curvature-relatedmass matrix Dz;4 -using Eq. (72) The mass matrix of the two-dimensional wall is m Dz;4 ¼ 3:375I ton m , and the bending-curvature-relatedmass matrix Dz;4 is: Stage-c. Define the bending-stiffness matrix components' term-see 4.1 Eq. (71)    The earthquake response analysis shows that RC wall ' = 8 is the closest to its ultimate curvature at z = 48.5 m (i.e., DOF n = 810), reaching about 82% of ductility demand; thus, it is addressed. Figure 22 depicts the corresponding area sections' cyclic behavior and the cyclic hysteretic portion at z = 48.5 m. Also, Fig. 23 illustrates the time-variation of the normalized bending moment and degradation yielding bending moment. The plots exemplify the cyclic smoothness throughout the analysis.
Showcasing the global structural response, Fig. 24 depicts the roof displacements in the x and y directions (also indicating their maximum) and the absorbed yielding energy at the base of RC wall ' = 8. Figure 25 refers to the maximum roof displacements and illustrates the whole structure deformation. The inelastic three-dimensional matrix equations model can easily provide these plots and more.

Conclusions
The paper develops novel matrix equations models for the linear/nonlinear dynamic analysis of two-dimensional and three-dimensional RC buildings with cantilever elements lateral load resisting systems. The matrix models are suitable for control theory's state-space representation and structural dynamics' equation of motion. The development starts with discretizing the PDE governing the dynamic equilibrium of a cantilever element subjected to lateral load while introducing the Euler-Bernoulli assumption. The LRS numerical integration scheme evaluates the two-dimensional element's displacement-related stiffness matrix while adhering to the structure's boundary conditions and yields a symmetric stiffness matrix. The mesh-size efficiency in terms of numerical accuracy shows that a mesh size of less than 1% of the structure's height results in less than 10 -2 relative errors in reference to the accurate analytical solution.
The three-dimensional structure's matrix equations employ the elements' stiffness matrices and combine them into the global stiffness matrix using the DSM. The elements' mass matrices and the ceilings' mass quantities are integrated into the global mass matrix using the proposed transformation matrices. The matrix equations model considers the element's lateral resistance in both directions and the element's lateral stiffness mutuality. The nonlinear RC mechanical behavior uses a smooth hysteretic model whose Fig. 25 Illustration of the structural deformation at the time of maximum roof displacements parameters refer to an RC wall cross section subjected to gravitational load and bending moment simultaneously. The matrix equations model refers to the bending-moment equilibrium to implement the smooth hysteretic model for deteriorating inelastic structures. Accordingly, in the nonlinear case, the LRS integration method defines the mass matrix in bending-coordinates rather than the bending-stiffness matrix. Working in bending coordinates is beneficial since it exempts matrix inversion-saving computational efforts.
The paper presents four examples for utilizing the developed matrix equations models, which provide a stage-by-stage description for defining the stiffness and mass matrices. Examples 1 and 2 address two-dimensional and three-dimensional structures of linearly elastic properties and illustrate their mode shapes. Example 3 defines the inelastic properties and performs pushover analysis. Lastly, Example 4 deals with modeling a three-dimensional inelastic structure and performs earthquake response analysis. The four examples demonstrate the practicability of the developed matrix equations models, which offer an idealized approach for nonlinear dynamics and control theory.
Data availability statement Data will be made available on reasonable request to the author.

Declarations
Conflict of interest The author declares that there is no conflict of interest regarding the publication of this paper.