Three-Dimensional Numerical Model for Seismic Analysis of Structures

Response Spectrum Analysis (RSA) and modal combination techniques are widely used to estimate the peak response of structures subjected to earthquake vibrations. This paper presents the numerical modeling and seismic analysis of structures. The response spectra based on Eurocode 8 is implemented for a five storey moment-resisting 3D structure to estimate peak response to earthquake-induced vibration. Eigenvector analysis is carried out to determine undamped free vibration mode shapes and frequencies. The lumped mass matrix, stiffness matrix of the full 3D structure is derived and natural periods, modal shapes, modal participation factors, and effective modal mass ratios are calculated. The structural response in different modes is combined to estimate the total dynamic response. The Square Root of Sum of Squares (SRSS), Complete Quadratic Combination (CQC), and Absolute Sum (AbsSum) rules are discussed and compared. The results of the SRSS and CQC rule for the examined building are almost identical as all cross-correlation coefficients are approximately equal to zero, and the responses of the individual modes are completely uncorrelated. CQC rule is appropriate when the natural frequencies are close to each other. AbsSum method shows an upper bound estimate of the total response since it assumes that all modal peaks occur at the same time. The simulation for numerical modeling and response spectrum analysis (RSA) is implemented in MATLAB® software.


Introduction
Many structures are vulnerable to earthquake excitation.
Reliable and practical numerical simulation of structures subjected to earthquake excitation is essential. In this paper, straightforward mathematical models were used, including a model of MDOF frame systems and lumped mass at the floor level [1]. Response Spectrum Method (RSM) was implemented to approximate the peak response of a structure to earthquake excitation. RSM was established based on linear structural behavior, well-separated natural modes, and classical damping assumptions [2]. The response spectrum was generated by solving the equation of motion for a group of single degree of freedom (SDOF) systems subjected to an accelerogram or conducting a Probabilistic Seismic Hazard Analysis (PSHA) to obtain a uniform-hazard spectrum [3], a scenario spectrum [4] or a conditional mean spectrum [5]. The peak value of each mode contribution to the response or the modal static response was calculated by static analysis of the building subjected to lateral forces. Then floor displacement, story drift, base shear, and base overturning moment was calculated. The peak total response is estimated by combining its peak modal values by the SRSS, CQC, and AbsSum rule [6]. Considering that the SRSS and CQC modal combination rules are based on random vibration theory, the total response should be interpreted as the mean of the peak values of the response of many earthquake excitation. The smooth spectrum can be the mean or median of the individual response spectra, or it can be a more conservative spectrum, such as the mean-plus-one-standard-deviation spectrum.
The SRSS rule assumes that the response of the modes is statistically independent, while the coupling between modes caused by modal damping is not taken into account [7]. AbsSum method combines the modal results by the absolute sum of all values, and fully correlated modes are considered. This method is over-conservative in structural design. The CQC uses the assumption of statistical coupling among nearly spaced Modes caused by modal damping [8] and degenerates the SRSS rule when damping is zero for all Modes. Like CQC, the General Modal Combination (GMC) technique is based on the statistical coupling between nearly-spaced Modes except for the Rosenblueth correlation coefficient. Increasing the modal damping increases the coupling among nearly-spaced modes; contrarily, SRSS degenerates if modal damping is zero [9]. If the loading direction of two spectra is unknown to consider the envelope of loading at all possible angles, the CQC3 method is used, which is the extension of the SRSS method [10].
Dynamic analyses based on a particular set of load-dependent Ritz vectors estimate more precise results than when the same number of natural free-vibration mode shapes is used [11]. The reason why Ritz vectors provide more accurate results is that they are generated by the spatial distribution of the dynamic loading. In contrast, the direct use of natural mode shapes neglects this fundamental data. Besides, the Ritz-vector algorithm automatically uses the advantages of static condensation and Guyan reduction. [12].
Since the late 1960s, various formulations for estimating the peak response to earthquake excitation have been demonstrated. Most of them are different in the mathematical expressions given for the correlation coefficient, namely, research carried out by E. Rosenblueth and J. Elorduy [13] or N. M. Newmark and E. Rosenblueth [14]. Also, the work done in 1981 by Armen Der Kiureghian [15], is now widely used.
The SRSS and CQC rules are derived according to random vibration theory, and they would be more accurate for earthquake excitation with a wide range of frequencies and less accurate for short-duration impulsive ground motions [16].

Model Description
A 5-story building with plan dimensions 11.2m × 14m was modeled using MATLAB software. The average inter-story height of floor one to five was assumed 3m, and the first-floor height was 4m. All beams feature typical cross-sections 0.3m × 0.6m with the lengths represented in Figure 1. There are five types of columns with square cross-sections with decreasing order of dimensions from the first floor to the top floor (0.55m, 0.45m, 0.40m, 0.35m, and 0.30m). The elastic modulus of the columns was considered to be 16×10 9 Pascal (Pa), and the shear-type behavior of the structure was used. Rigid floor diaphragms in their planes and infinite axial deformability of columns are applied. Two translational and one rotational degree of freedom for each slab is assumed, so the 5-story building has 15 degrees of freedom. Table 1 represents the properties of the column of each floor. values represent the columns stiffness in the X and Y direction.  The dead load includes the self-weight of the slabs, beams, and columns. The contribution of gravity loads to the effective seismic weight was obtained by the combination rule established in equation 3-17 of EN 1998EN -1 (2004 -section 3.2.4 [17]. The effect of vertical loads must be combined with those of seismic actions [18]. The following combination was used for the quasi-permanent loading condition, according to Eq (1): accounting for the summary of dead loads and the combination coefficients , take into account the likelihood of the live loads ( , ). Permanent loads (G) and Live loads (Q) of floor slabs are 7.0 KN/m 2 , and 2.0 KN/m 2 , and for the roof 6.0 KN/m 2 and 1.3 KN/m 2 respectively. Permanent and Live loads coefficients are assumed 1.00 and 0.30, respectively.
Mass matrix is composed of three 5×5 diagonal matrices, , . Where = is the inertial translational masses of the ith floor, and matrix computed from the polar (rotational) moment of inertia of the slabs [19]. Since there are 3 degrees of freedom per floor, the mass matrix of the 3D structure is 15×15 diagonal. (2) Code-based elastic and design spectral accelerations were implemented for the five-story shear 3D frame to estimate peak dynamic responses subjected to earthquake excitation in the direction in which the building has uncoupled modes (ground motion in the y-direction). According to EN1998-3 Section 2.1, performance requirements refer to the state of damage in the structure defined through three limit states, namely Near Collapse (NC), Significant Damage (SD), and Damage Limitation (DL) [20]. For the SD limit state, 5.0% damping ratio, and peak ground acceleration (PGA) of 0.160g for the return period of 475 years is assumed. The soil class B of Type 1, according to EN1998-1 definition of ground types, is considered [17]. The elastic spectral accelerations are divided by the behavior factor (q) equal to 3.45 to draw the design spectral accelerations at the SD limit state.
To obtain the stiffness matrix of the 3D structure, first, we computed the frame's 2D lateral stiffness matrix in local coordinates in x and y directions. Shear-type behavior of frame is used; it means that the beams of each frame have infinite flexural stiffness, and it was assumed that all structural elements do not have axial deformability. In this case, the structure in each x and y-direction has only 5 degrees of freedom (only the horizontal displacement). Direct Stiffness Method or Unit Dummy Force Method has been used to derive the stiffness matrix stiffness. Frame stiffness matrices in x and y direction are described by Eq (3) and Eq (4): (3) Local frame's stiffness matrices ( ) were correlated with the global stiffness matrix of the 3D structure ( ) by a transformation. Transformation Matrices were adopted to relate the displacements of each floor of a jth frame to the displacements of the center of mass. The generic transformation matrix for the jth frame in x and y direction are: = � 5,5 ; 5,5 ; − 5,5 �, = � 5,5 ; 5,5 ; + 5,5 � The center of mass ( ) for each floor is assumed in the geometrical center of the slab, and this is the same for all floors. G Since there is symmetry concerning y-direction, the center of stiffness ( ) must be on the y-axis ( = 0). X and y are the K displacements between each frame and the So, the stiffness matrix of the 3D structure is equal to Eq (6): G.
After defining the and Modal analysis is performed for the 3D structure to find natural frequencies, Typically using linear superposition, one component of ground motion at a time is studied. The influence vector ( ) represents the ground motion directions in x, y, and θ. in the y-direction is = [ × × × ] , i.e., one in the excitation direction (y) and zero in rotational and translational excitation in the x-direction. When the structure is subject to a ground-motion in the y-direction, only five uncoupled translational modes in the y-direction are activated, and coupled modes that have rotational terms are not activated due to structural symmetry in the y-direction. The design equivalent static forces at the ith floor of the structure for n-modes were calculated by: Where is the mass of the ith floor, and n is the considered mode. is used to calculate the design story shear and column's bending moments for the LS limit state of each mode. For each mode the floor displacements vector (u) and Inter-storey drift ratio for the DL and LS limit state are calculated using the equation: For the LS limit state to consider the nonlinear displacements of structure, the Deflection Amplification Factor assumed equal to the Behavior Factor (q=3.45).
Since the slabs are rigid and the columns have the same displacement and stiffness in each floor, the shear diagram for each column is calculated in the y-direction. Also, it is possible to calculate the bending moment at each column, considering that it behaves as a doubly clamped column with lateral displacement.
The estimate for peak total response obtained by the CQC rule depends on the cross-correlation coefficient ( ) for the two modes. The equation for the correlation coefficient, according to Der Kiureghian is: This equation simplifies to Eq. (9) by considering each mode as the same damping ratio.

Results and Discussions
The stiffness matrix of the 3D structure calculated according to Eq. (6) is shown in the Matrix (10): We computed the Mass matrix with the assumption of lumped mass on the floors; therefore, the Mass matrix is diagonal and off-diagonal terms are equal to zero. Inertial translational mass matrices ( ) and matrix of the rotational moment of inertia ( ) are presented in matrices (11): The calculated natural frequencies ( ) and natural periods ( ) are reported in vector (12) and vector (13) After calculating the Natural frequencies and Periods, the corresponding eigenvectors or Modal shapes ( ) was calculated. The maximum amplitude was normalized to a value of one. The resulting matrix of modal shapes is presented in Eq. (14). The analyzed building is symmetric about the y-axis, therefore Modes 3, 7, 10, 12, 14 have zero terms in the x and θ directions. In other words, the examined model is 'regular' building in the y-direction in which the centers of mass and stiffness are coincident, resulting in uncoupled modes with well-separated periods when the earthquake excitation is in the y-direction. The examined building is 'irregular' with the stiffness offset from the mass center of the diaphragm resulting in coupled modes with the rotational components. So, if the building is subjected to ground motion in the y-direction, since there is symmetry concerning the mass and stiffness, modes 3, 7, 10, 12, and 14 will be activated just with translational displacements. (14) It is possible to define two separate vectors, considering the components of the ground acceleration in the x and y-direction separately. The Base Shear Effective Modal Mass Ratios for ground acceleration in x and y direction are shown by vector (15) and vector (16). The Effective Modal Mass ratio vector obtained for a ground motion in the y-direction only has five components different from zero and the first mode in this with � = 70% has the highest contribution to the total base shear, while the vector for a ground motion in x-direction has ten non zero components. The Elastic Spectral Pseudo-accelerations for DL and SD limit state and accelerations at the natural periods of the structures are presented in Figure 2. The design spectral acceleration at the SD limit state was plotted by dividing the elastic spectrum to the behavior factor q = 3.45. In DL limit state concept structure should be designed to resist the seismic force with a larger probability of occurrence, without damage experience and the associated limitations of use.
It was assumed that the structure is subject to a ground acceleration in the y-direction, and due to symmetry concerning Three-Dimensional Numerical Model for Seismic Analysis of Structures the y-axis, just five translational modes are activated without rotational components. The design equivalent static forces for each mode were presented in each column of Eq. (17). As can be seen, there are five columns that each relates to one mode. For the first Mode (first column), the seismic design force at the building roof level is 101.91 kN, and at the first-floor level is 24.37 kN. Using the equivalent lateral force of the first Mode, as described in ASCE 7-16, is a method for regular structures for which the fundamental period of vibration is dominant in mass participation, although RSA is more sophisticated and requires more knowledge of structural dynamics [21].
So far, the dynamic responses were evaluated for each mode. To capture the effect of all modes, Modal combinations, e.g., SRSS, CQC, and AbsSum, were done as the last step of the analysis.
Variation of the correlation coefficients concerning frequency ratio with a five percent damping ratio is plotted in Figure 3. Blue stars represent correlation coefficients at the natural periods of the structure, note that the two equations proposed by two different researchers provide approximately equal values. However, with an increase in the damping ratio, the difference is more noticeable. The correlation coefficient decreases or vanishes quickly as the two natural frequencies move farther apart. When the two frequencies are the same, and in a narrow range, the correlation coefficient is equal to one. In other words, Modes with closer eigenvalues may have higher correlation, and with well-separated natural frequencies, correlation coefficients vanish. Therefore, all cross-correlations in the CQC rule can be neglected, and the SRSS rule is generated.  Eurocode 8 Part 1 presents guidelines on the inter-storey drift ratio to limit the structural and non-structural damage due to frequent and rare seismic events. Also, inter-storey drifts should be limited to reduce columns P-∆ effects. Displacements produced by the seismic design action has to be calculated based on the elastic deformations of the structures multiplied by displacement behavior factor to take into account the inelastic deformations. Figure 4 shows the resulting absolute storey displacements for the first two modes at the DL and SD limit states. The inter-storey drift ratio is about less than 0.2% for the DL limit state, while the inter-storey drift ratio of the structure for the SD limit state is 0.7%. As it can be seen in Figure 5, the design equivalent static forces of different modes are combined with three different modal combination methods. At the level of floor five, there are the highest values for the static force, with about 120 kN with a decreasing trend toward the first floor. The values comparison obtained from both SRSS and CQC methods is nearly identical. This is because the cross-correlation coefficients in the CQC method are equal to zero, and the SRSS method was generated. AbsSum results with 200 kN at the top floor indicated about 70% more with respect to the other two methods in the estimation of seismic design forces. If we consider a specific performance level and spectra for the building, AbsSum overestimated the seismic forces induced by earthquake excitation. For each mode, the shear at the various floors was calculated by considering the design equivalent static forces vector related to the mode. From the roof downwards, at each floor level at the columns, the shear is the summation of the forces of all the floors above. As it is shown in Figure 6, the shear force of the floors increases from top to bottom. The storey shear force calculated by the CQC method is approximately the same as the SRSS method, while, in general, the AbsSum is higher than both methods about 10-30 percent.
Shear force calculated by the CQC and SRSS methods is approximately 352 kN and 340 kN for the first floor, respectively, while, this value in the AbsSum method is 494 kN. The storey shear of the fifth floor is about 200 kN in AbsSum method, while this value is about 121 kN for the other two methods.  Design bending moments at the LS limit state for the columns calculated with three different modal combination methods are represented in Figure 7. From the top of the structure downward, the structural members bending moment increases. In this case, in general, AbsSum rule shows 40% more compared with two other modal combination methods. The results for the two CQC and SRSS methods are very close to each other. At the top of the building on the fifth floor, the maximum bending anchor is 25.11 kN.m in the AbsSum method, while this value is 15.2 kN.m for the other two methods. The column's bending moment at the level of the first floor calculated approximately 58 kN.m in CQC and SRSS method, although this value for the AbsSum method is about 86.6 kN.m. Another observation is that CQC results may be larger or smaller than the estimate provided by the SRSS rule.

Conclusions
A crucial review and theoretical development of modal-combination methods usually used in the seismic analysis of structures have been evaluated. The main hypotheses of each technique have been discussed, and the strengths and weaknesses have been summarized. A five storey moment-resisting structure has been used to examine the relative performance of the modal combination methods. The results obtained from the SRSS and CQC rules have been found to be approximately similar to each other. The SRSS rule should be applied to structures with well-separated natural frequencies of those modes that contribute significantly to the response. The work has shown that if the periods of significant modes are different from each other by more than 10%, the SRSS and CQC combination method was equally efficient. In other words, it has been proven, the SRSS rule is generated from the CQC rule for the structures with well-separated natural frequencies. The CQC rule provides excellent response estimates for structures with closely spaced natural frequencies, such as multistory buildings with an unsymmetric plan and a broader class of structures as it overcomes the limitations of the SRSS rule. In addition, it has been shown that the assumption that all modal peaks occur at the same time, and their algebraic sign is ignored results in an upper bound estimate of the total response in AbsSum rule.