Study on the Seismic Response of a Portal Frame Structure Based on the Transfer Matrix Method of Multibody System

Abstract

Portal frame structures are widely used in industrial building design but unfortunately are often damaged during an earthquake. As a result, a study on the seismic response of this type of structure is important to both human safety and future building designs. Traditionally, finite element methods such as the ANSYS and MIDAS have been used as the primary methods of computing the response of such a structure during an earthquake; however, these methods yield low calculation efficiencies. In this paper, the mechanical model of a single-story portal frame structure with two spans is constructed based on the transfer matrix method of multibody system (MS-TMM); both the transfer matrix of the components in the model and the total transfer matrix equation of the structure are derived, and the corresponding MATLAB program is compiled to determine the natural period and seismic response of the structure. The results show that the results based on the MS-TMM are similar to those obtained by ANSYS, but the calculation time of the MS-TMM method is only 1/20 of that of the ANSYS method. Additionally, it is shown that the MS-TMM method greatly increases the calculation efficiency while maintaining accuracy.

1. Introduction

Because portal frame structures have many advantages, which include good aesthetics, low cost, and rapid construction [1], this type of structure is widely used in small and medium-sized industrial buildings, warehouses, and large supermarkets [2]. Additionally, China is a country in which earthquakes frequently occur, according to readily available earthquake disaster data [3]. Although steel structures have some advantages, they are still regularly damaged during an earthquake, requiring special attention and study, particularly with regard to portal frame structures [4].

With the increasing demand for a calculation method for the seismic response of a building, low calculation efficiency has become a serious hindrance to many methods, preventing them from being widely applied in the field of civil engineering [5]. Therefore, searching for a simple, fast, and accurate earthquake response analysis method is an important area of research that engineers and scholars have been working towards for many years [6]. Currently, the dynamic response with a time history for a portal frame structure under the action of an earthquake is generally calculated by a finite element software [7, 8], such as ANSYS [9] or MIDAS; unfortunately, for the complex civil structure with a large number of degrees of freedom, the computational work load is heavy, which cause long calculation times and low calculation efficiencies.

The transfer matrix method of multibody system (MS-TMM) [10] was developed by Xiaoting Rui in 1993 and was mainly divided into the transfer matrix method of linear multibody system [11] and discrete time transfer matrix method [12]. The former method mainly determined the vibration characteristics of a linear multibody system [13], while the latter mainly solved the time-varying, nonlinear, large motion, controlled, or another general multibody system dynamic problem [14]. Compared to the classical transfer matrix method, MS-TMM has the following advantages: (1) no global dynamics equations of the system are required; (2) only high-level programming is required; (3) a low order of system matrix usually results; (4) high computational efficiency is existed [15]; and (5) the “morbidity” problems inherent to eigenvalue calculations are avoided [16]. At present, the transfer matrix method of multibody system has already been applied in the field of civil engineering [17]. Additionally, the transfer matrix of some typical elements of civil engineering structures has been established.

In this paper, the MS-TMM method is used to calculate and analyze the vibration characteristics and seismic response of a portal frame structure, which are then compared to the results obtained by the traditional ANSYS method. The calculation efficiencies of these two methods are also compared.

2. Mechanical Model of Portal Frame Structure

The general idea of the transfer matrix method of a multibody system is first to divide the complex multibody system into several components by using discrete thought. The components are considered to contain two types: a “body” and a “hinge.” The “body” mainly includes the rigid body, the flexible body and a concentrated mass. The “hinge” is the connection between one “body” and the next “body”; it is massless because the mass of hinge is included in the adjacent “body.” Next, a matrix is used to express the mechanical properties of each element. Finally, “assembly” of the transfer matrix of each element is combined to obtain the total transfer equation and the total transfer matrix of the system.

In this paper, a portal frame structure is considered as a multibody system. According to the features of the multibody system for a portal frame structure, beams and columns are divided into particles (i.e., the “bodies”) and massless beams (i.e., the “hinges”); massless beams are used to link particles together; the axial deformations of the beams are ignored. Because each portal frame in the longitudinal direction is nearly symmetrical, the original spatial system can be simplified into a planar system. Now, a representative planar portal frame model of the structure is selected to study its response during an earthquake; this mechanical model is shown in Figure 1. The elements of the columns and the inclined beams include the massless beams and particles. The massless beams of the left, middle, and right columns are labeled as 1, 3, …, j–1; 1′, 3′,…, j′ – 1; and 1″, 3″,…, j″ – 1, respectively, and the particles are labeled as 2, 4,…, j–2,j; 2′, 4′,…, j′-2,j′; and 2″, 4″,…, j″-2,j″, respectively. The massless beams of the left and right inclined beams are labeled as l1,l3,…, li and r1,r3,…, ri, and the particles are labeled as l2,l4,…, li–1 and r2,r4,…, ri–1.

Figure 1:

Mechanical model of the multibody system for a portal frame structure.

3. Calculation Principles of the MS-TMM

3.1. Transfer Path and Method of the Multibody System for a Portal Frame Structure Based on Its Characteristics

Because the mechanical model of the multibody system for a portal frame structure is a bifurcated system, and the number of elements is large, the transfer equation will thus be very complex when using the traditional transfer matrix method of multibody system. Therefore, in this paper, the physical links between the beams and columns are assumed to be disconnected, and the corresponding relations at the nodes including j, j′, and j″ are replaced by state vectors based on the deformation compatibility relation of the inclined beams. As shown in Figure 2, the main research objects are the three steel columns, each of which is a chain system, and the state vector of each point is transferred in the direction shown by the arrow.

Figure 2:

Calculation principles of the transfer matrix of the multibody system for a portal frame structure.

The left column has one input (Z_0,1) and one output (Z_j,l1), the middle column has two inputs (Z_{0′, 1′}, Z_j′li) and one output (Z_{j′, r1}), and the right column has one input (Z_{j″, ri}) and one output (Z_{0″, 1″}). Among them, the column caps’ state vectors Z_j,l1 and Z_{j′, li} can be associated with the deformation compatibility relation of the left inclined beam; Z_{j′, r1} and Z_{j″, ri} can be associated with that of the right inclined beam. Finally, the relation between the three boundary state vectors Z_{0, 1}, Z_{0′, 1′}, and Z_{0″, 1″} must be established; this will allow the transfer relationship between the state vectors to yield each point state vector in this system.

The mechanical model of the portal frame structure involves only two elements: a massless beam and particle. The transfer matrix method of a linear multibody system is used to solve the natural period of the structure, and the discrete time transfer matrix method is used to analyze the seismic response of the structure to both common and rare earthquakes.

3.2. Analysis of the Natural Period of the Structure Based on the Transfer Matrix Method of Linear Multibody System

The transfer matrix method of linear multibody system is used to determine the natural period of the structure; the mechanical model of the portal frame structure used is shown in Figure 1. The state vector under physical coordinates of the element in plane motion is defined as $z={\left[x,y,{\theta }_z,m_z,q_x,q_y\right]}^T$ , where x,y,θ _z ,m _z ,q _x ,q _y denote the displacements in the x and y directions, the angular displacement around z, the internal moment, the internal force in the x direction, and the internal force in the y direction, respectively. The corresponding state vector under the modal coordinates is $\mathbf{Z}={\left[X,Y,{\Theta }_z,M_z,Q_x,Q_y\right]}^T$ . The transfer process is derived as follows. For the left column:

$\begin{array}{c}\hspace{0.5em}{\mathbf{Z}}_{\mathrm{2,1}}={\mathbf{U}}_{\mathrm{1}}{\mathbf{Z}}_{\mathrm{0,1}}\\ {\mathbf{Z}}_{\mathrm{2,3}}={\mathbf{U}}_{\mathrm{2}}{\mathbf{Z}}_{\mathrm{2,1}}\\ ⋮\\ {\mathbf{Z}}_{j,j-\mathrm{1}}={\mathbf{U}}_{j-\mathrm{1}}{\mathbf{Z}}_{j-\mathrm{2},j-\mathrm{1}}\\ {\mathbf{Z}}_{j,l\mathrm{1}}=\mathbf{H}{\mathbf{U}}_j{\mathbf{Z}}_{j,j-\mathrm{1}};\end{array}$

(1)

for the middle column:

$\begin{array}{c}\hspace{0.5em}\hspace{0.17em}{\mathbf{Z}}_{{\mathrm{2}}^{\prime },{\mathrm{1}}^{\prime }}={\mathbf{U}}_{{\mathrm{1}}^{\prime }}{\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}\\ \hspace{0.17em}{\mathbf{Z}}_{{\mathrm{2}}^{\prime },{\mathrm{3}}^{\prime }}={\mathbf{U}}_{{\mathrm{2}}^{\prime }}{\mathbf{Z}}_{{\mathrm{2}}^{\prime },{\mathrm{1}}^{\prime }}\\ ⋮\\ {\mathbf{Z}}_{j^{\prime },j^{\prime }-\mathrm{1}}={\mathbf{U}}_{j^{\prime }-\mathrm{1}}{\mathbf{Z}}_{j^{\prime }-\mathrm{2},j^{\prime }-\mathrm{1}};\end{array}$

(2)

for the right column:

$\begin{array}{c}\hspace{1em}{\mathbf{Z}}_{j^{\mathrm{\prime \prime }},j^{\mathrm{\prime \prime }}-\mathrm{1}}=\mathbf{H}{\mathbf{U}}_{j^{\prime }}{\mathbf{Z}}_{j^{\prime },ri}\\ {\mathbf{Z}}_{j^{\mathrm{\prime \prime }}-\mathrm{2},j^{\mathrm{\prime \prime }}-\mathrm{1}}={\mathbf{U}}_{j^{\mathrm{\prime \prime }}-\mathrm{1}}{\mathbf{Z}}_{j^{\mathrm{\prime \prime }},j^{\mathrm{\prime \prime }}-\mathrm{1}}\\ ⋮\\ {\mathbf{Z}}_{{\mathrm{2}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}={\mathbf{U}}_{{\mathrm{2}}^{\mathrm{\prime \prime }}}{\mathbf{Z}}_{{\mathrm{2}}^{\mathrm{\prime \prime }},{\mathrm{3}}^{\mathrm{\prime \prime }}}\\ {\mathbf{Z}}_{{\mathrm{0}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}={\mathbf{U}}_{{\mathrm{1}}^{\mathrm{\prime \prime }}}{\mathbf{Z}}_{{\mathrm{2}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}.\end{array}$

(3)

To unify the coordinate systems, H is defined as a coordinate transformation matrix:

$\begin{array}{c}\mathbf{H}=\left[\begin{array}{cc}{\mathbf{H}}_{\mathrm{1,1}}& {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}\\ {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}& {\mathbf{H}}_{\mathrm{2,2}}\end{array}\right],\\ {\mathbf{H}}_{\mathrm{1,1}}=\left[\begin{array}{ccc}\cos \theta & \sin \theta & \mathrm{0}\\ -\sin \theta & \cos \theta & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],\\ {\mathbf{H}}_{\mathrm{2,2}}=\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \cos \theta & \sin \theta \\ \mathrm{0}& -\sin \theta & \cos \theta \end{array}\right].\end{array}$

(4)

The three steel columns of the portal frame structure can be viewed as three chain systems; the resulting transfer equations are as follows:

$\begin{array}{c}{\mathbf{Z}}_{j,l\mathrm{1}}={\mathbf{U}}_{{\rm I}}{\mathbf{Z}}_{\mathrm{0,1}}\\ {\mathbf{Z}}_{j^{\prime },j^{\prime }-\mathrm{1}}={\mathbf{U}}_{{\rm I}{\rm I}}{\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}\\ {\mathbf{Z}}_{{\mathrm{0}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}={\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{Z}}_{j^{\mathrm{\prime \prime }},ri},\end{array}$

(5)

where

$\begin{array}{c}{\mathbf{U}}_{{\rm I}}={\mathbf{U}}_{\mathrm{1}}{\mathbf{U}}_{\mathrm{2}}\cdots {\mathbf{U}}_{j-\mathrm{1}}\mathbf{H}{\mathbf{U}}_j\\ {\mathbf{U}}_{{\rm I}{\rm I}}={\mathbf{U}}_{{\mathrm{1}}^{\prime }}{\mathbf{U}}_{{\mathrm{2}}^{\prime }}\cdots {\mathbf{U}}_{j^{\prime }-\mathrm{1}}\\ {\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}={\mathbf{U}}_{{\mathrm{1}}^{\mathrm{\prime \prime }}}{\mathbf{U}}_{{\mathrm{2}}^{\mathrm{\prime \prime }}}\cdots {\mathbf{U}}_{j^{\mathrm{\prime \prime }}-\mathrm{1}}\mathbf{H}{\mathbf{U}}_{j^{\mathrm{\prime \prime }}}.\end{array}$

(6)

The element j′ has two inputs (Z_{j′, li}, Z_{j′, j′ – 1}) and one output (Z_{j′, r1}). The input state vector can be defined as follows:

$\begin{array}{c}{\mathbf{Z}}_{I,j^{\prime }}=\left[X_{j^{\prime },li},Y_{j^{\prime },li},{\Theta }_{j^{\prime },li},M_{z{j}^{\prime },li},Q_{x{j}^{\prime },li},Q_{y{j}^{\prime },li},\right){\left(M_{z{j}^{\prime },j^{\prime }-\mathrm{1}},Q_{x{j}^{\prime },j^{\prime }-\mathrm{1}},Q_{y{j}^{\prime },j^{\prime }-\mathrm{1}}\right]}^T\\ \hspace{1em}={\mathbf{E}}_{\mathrm{1}}{\mathbf{Z}}_{j^{\prime },li}+{\mathbf{E}}_{\mathrm{2}}{\mathbf{Z}}_{j^{\prime },j^{\prime }-\mathrm{1}},\end{array}$

(7)

where

${\mathbf{E}}_{\mathrm{1}}=\left[\begin{array}{c}{\mathbf{I}}_{\mathrm{6}}\\ {\mathbf{O}}_{\mathrm{3}\times \mathrm{6}}\end{array}\right],\hspace{1em}\hspace{1em}{\mathbf{E}}_{\mathrm{2}}=\left[\begin{array}{c}{\mathbf{O}}_{\mathrm{6}\times \mathrm{3}}\\ {\mathbf{O}}_{\mathrm{3}\times \mathrm{3}}\end{array}\begin{array}{c}{\mathbf{O}}_{\mathrm{6}\times \mathrm{3}}\\ {\mathbf{I}}_{\mathrm{3}}\end{array}\right].$

(8)

The output state vector is

${\mathbf{Z}}_{j^{\prime },r\mathrm{1}}={\mathbf{U}}_{j^{\prime }}{\mathbf{Z}}_{I,j^{\prime }}.$

(9)

According to the deformation compatibility relation of the left and right inclined beams, respectively, the following is derived:

${\mathbf{Z}}_{j^{\prime },li}={\mathbf{U}}_L{\mathbf{Z}}_{j,l1},\hspace{1em}\hspace{1em}{\mathbf{Z}}_{j,ri}={\mathbf{U}}_R{\mathbf{Z}}_{j^{\prime },r1}.$

(10)

For this model, UL, UR are the products of the transfer matrices of the elements in left and right inclined beam, respectively, which are

$\begin{array}{c}{\mathbf{U}}_L={\mathbf{U}}_{l\mathrm{1}}{\mathbf{U}}_{l\mathrm{2}}\cdots {\mathbf{U}}_{li-\mathrm{1}}{\mathbf{U}}_{li}\\ {\mathbf{U}}_R={\mathbf{U}}_{r\mathrm{1}}{\mathbf{U}}_{r\mathrm{2}}\cdots {\mathbf{U}}_{ri-\mathrm{1}}{\mathbf{U}}_{ri}\end{array}$

(11)

Therefore,

$\begin{array}{c}{\mathbf{Z}}_{j^{\prime },r\mathrm{1}}={\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{1}}{\mathbf{Z}}_{j^{\prime },li}+{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{2}}{\mathbf{Z}}_{j^{\prime },j^{\prime }-\mathrm{1}}\\ \hspace{1em}={\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{1}}{\mathbf{U}}_L{\mathbf{U}}_I{\mathbf{Z}}_{\mathrm{0,1}}+{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{2}}{\mathbf{U}}_{\text{II}}{\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}\end{array}$

(12)

Because the element j′ is a particle, all displacements are equivalent and the following is true:

$\begin{array}{c}{\mathbf{Z}}_{{\mathrm{0}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}={\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{Z}}_{j^{\mathrm{\prime \prime }},ri}={\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{U}}_R{\mathbf{Z}}_{j^{\prime },r\mathrm{1}}\\ \hspace{1em}={\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{U}}_R{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{1}}{\mathbf{U}}_L{\mathbf{U}}_{{\rm I}}{\mathbf{Z}}_{\mathrm{0,1}}+{\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{U}}_R{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{2}}{\mathbf{U}}_{{\rm I}{\rm I}}{\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}.\end{array}$

(13)

where

${\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{U}}_R{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{1}}{\mathbf{U}}_L{\mathbf{U}}_{{\rm I}}{\mathbf{Z}}_{\mathrm{0,1}}+{\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{U}}_R{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{2}}{\mathbf{U}}_{{\rm I}{\rm I}}{\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}-{\mathbf{Z}}_{{\mathrm{0}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}=\mathrm{0}.$

(14)

This yields

${\mathbf{E}}_{\mathrm{3}}{\mathbf{Z}}_{j^{\prime },li}={\mathbf{E}}_{\mathrm{3}}{\mathbf{Z}}_{j^{\prime },j^{\prime }-\mathrm{1}}={\mathbf{E}}_{\mathrm{3}}{\mathbf{Z}}_{j^{\prime },r\mathrm{1}},$

(15)

Combining equations (14) and (17) yields

${\mathbf{E}}_{\mathrm{3}}=\left[\begin{array}{cc}{\mathbf{I}}_{\mathrm{3}}& {\mathbf{O}}_{\mathrm{3}\mathrm{\hspace{0.17em}\hspace{0.17em}}\times \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{3}}\end{array}\right].$

(16)

where

${\mathbf{E}}_{\mathrm{3}}{\mathbf{U}}_L{\mathbf{U}}_I{\mathbf{Z}}_{\mathrm{0,1}}-{\mathbf{E}}_{\mathrm{3}}{\mathbf{U}}_{\text{II}}{\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}=\mathrm{0}.$

(17)

The relationship between the state vectors under the physical and modal coordinates is z = Ze ^iωt . Thus, the transfer matrix of the massless beam can be derived as follows:

${\mathbf{Z}}_{\text{all}}{\mathbf{U}}_{\text{all}}=\mathbf{0}\mathbf{,}$

(18)

The transfer matrix of a particle is expressed as follows:

$\begin{array}{c}{\mathbf{U}}_{\text{all}}=\left[\begin{array}{ccc}{\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{U}}_R{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{1}}{\mathbf{U}}_L{\mathbf{U}}_{{\rm I}}& {\mathbf{U}}_{{\rm I}{\rm I}{\rm I}}{\mathbf{U}}_R{\mathbf{U}}_{j^{\prime }}{\mathbf{E}}_{\mathrm{2}}{\mathbf{U}}_{{\rm I}{\rm I}}& -{\mathbf{I}}_{\mathrm{6}}\\ {\mathbf{E}}_{\mathrm{3}}{\mathbf{U}}_L{\mathbf{U}}_{{\rm I}}& -{\mathbf{E}}_{\mathrm{3}}{\mathbf{U}}_{{\rm I}{\rm I}}& {\mathbf{O}}_{\mathrm{3}\mathrm{\hspace{0.17em}\hspace{0.17em}}\times \mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{6}}\end{array}\right],\\ {\mathbf{Z}}_{\text{all}}=\left[\begin{array}{c}{\mathbf{Z}}_{\mathrm{0,1}}\\ {\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}\\ {\mathbf{Z}}_{{\mathrm{0}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}\end{array}\right].\end{array}$

(19)

with the following boundary conditions:

$\mathbf{U}=\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& l& \frac{l^{\mathrm{2}}}{\mathrm{2}EI}& \mathrm{0}& \frac{l^{\mathrm{3}}}{\mathrm{6}EI}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \frac{l}{EI}& \mathrm{0}& \frac{l^{\mathrm{2}}}{\mathrm{2}EI}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& l\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right].$

(20)

Zero elements can be removed from Z_all, and the elements in the corresponding column also can be removed. This yields the equation ${\overline{\mathbf{Z}}}_{\text{all}}{\overline{\mathbf{U}}}_{\text{all}}=\mathbf{0}$ . Based on the characteristic equation $\Delta =\left|{\overline{\mathbf{U}}}_{\text{all}}\right|=\det ({\overline{\mathbf{U}}}_{\text{all}})=\mathrm{0}$ , the natural frequency ω can be obtained using the MATLAB programming language.

3.3. Analysis of the Seismic Response of the Structure Based on the Discrete Time Transfer Matrix Method

The discrete time transfer matrix method of a multibody system is used to determine seismic response of the structure during common and rare earthquakes; an elastic-plastic restoring force model is also established to analyze the response to a rare earthquake. Numerical solutions are typically used because analytic solutions are usually not available when the change in external load is complex. Generally, during the physical process at each time step, the relationship between the physical variables can be approximately regarded as linear as long as the time step is small. In this paper, the equations are linearized on the basis of the Newmark-β method, which provides an efficient solution.

3.3.1. Extended Transfer Matrix of a Massless Beam and a Particle in a Portal Frame Structure

An external force f should be considered in the transfer matrices of elements. The state vector of the element in plane motion is thus defined as $z={\left[x,y,,{\theta }_z,m_z,q_x,q_y,\mathrm{1}\right]}^T$ , where x,y,θ _z ,m _z ,q _x ,q _y denote the displacements in the x and y directions, the angular displacement around z, the internal moment, the internal force in the x direction, and the internal force in the y direction, respectively.

The extended transfer matrix of the massless beam is as follows:

$\mathbf{U}=\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& m{\omega }^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]$

(21)

The extended transfer matrix of the particle is as follows:

$\begin{array}{c}{\mathbf{Z}}_{\mathrm{0,1}}={\left[\mathrm{0,0},\mathrm{0},M_z,Q_x,Q_y\right]}^T,\\ {\mathbf{Z}}_{{\mathrm{0}}^{\prime },{\mathrm{1}}^{\prime }}={\left[\mathrm{0,0},\mathrm{0},M_z,Q_x,Q_y\right]}^T,\\ {\mathbf{Z}}_{{\mathrm{0}}^{\mathrm{\prime \prime }},{\mathrm{1}}^{\mathrm{\prime \prime }}}={\left[\mathrm{0,0},\mathrm{0},M_z,Q_x,Q_y\right]}^T.\end{array}$

(22)

where m is the mass of the particle; c_ζ is the damping of the structure; $f_{y,c}=-m\stackrel{••}{x_g}(t_{i-\mathrm{1}})$ when the earthquake acceleration is $\stackrel{••}{x_g}(t_{i-\mathrm{1}})$ at t = t_i–1. Moreover,

$\mathbf{U}=\left[\begin{array}{ccccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& l& \frac{l^{\mathrm{2}}}{\mathrm{2}EI}& \mathrm{0}& \frac{l^{\mathrm{3}}}{\mathrm{6}EI}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \frac{l}{EI}& \mathrm{0}& \frac{l^{\mathrm{2}}}{\mathrm{2}EI}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& l& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right].$

(23)

3.3.2. Calculation Process of Solving the Response of the Structure

The calculation process is shown as Figure 3; MATLAB programming software was used to calculate the solutions shown.

Figure 3:

The process of solving the seismic response of a structure using MS-TMM.

4. Engineering Example

In this section, a single-story double-span portal frame structure is used as an example. The span is 18 m × 2; the wall and roof panel are both made of colored profiled sheets whose theoretical areal masses are 6.5 kg/m². The grade of the loads is I, and the snow load is 0.35 kg/m². For a steel structure, E = 2.06 × 10¹¹ N/m².

The seismic fortification intensity is set to 7, the site type is II, and the designed earthquake is designated as a 1st classification. According to the Code for Seismic Design of Buildings of China, when the response of the structure is studied under common and rare earthquakes, the peak value of the acceleration versus time data is adjusted to 35 cm/s² and 220 cm/s², respectively. In this study, a natural earthquake wave (i.e., the El Centro earthquake wave) is chosen for the analysis of the seismic response of the structure. Due to the earthquake duration time, which is required to be from 5 to 10 times the structure's natural period, according to the Code for Seismic Design of Buildings of China, the earthquake duration time is set to 10's in this study.

The calculation model of the portal frame structure is shown as Figure 4, and the calculation parameters are shown as follows:

$\mathbf{U}=\left[\begin{array}{ccccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -(mA+c_{\zeta }C)& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -mB+c_{\zeta }D+f_{y,c}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],$

(24)

$\begin{array}{c}A=\frac{\mathrm{1}}{\alpha \Delta t^{\mathrm{2}}},\hspace{1em}\hspace{1em}B=A\left[-y_{t_{i-1}}-{\stackrel{•}{y}}_{t_{i-\mathrm{1}}}\Delta t-\left(\frac{\mathrm{1}}{\mathrm{2}}-\alpha \right){\stackrel{••}{y}}_{t_{i-\mathrm{1}}}\Delta t^{\mathrm{2}}\right],\\ C=\beta \Delta tA,\hspace{1em}\hspace{1em}D={\stackrel{•}{y}}_{t_{i-\mathrm{1}}}+\left(\mathrm{1}-\beta \right)\Delta t{\stackrel{••}{y}}_{t_{i-\mathrm{1}}}+\beta \Delta tB\left(\alpha =\frac{\mathrm{1}}{\mathrm{2}},\beta =\frac{\mathrm{1}}{\mathrm{4}}\right).\end{array}$

(25)

$\begin{array}{c}{m}_{\mathrm{2}}=m_{{\mathrm{2}}^{\mathrm{\prime \prime }}}=\mathrm{150.741}\hspace{0.17em}\text{kg},\hspace{1em}\hspace{1em}{m}_{\mathrm{4}}=m_{{\mathrm{4}}^{\mathrm{\prime \prime }}}=\mathrm{153.387}\hspace{0.17em}\text{kg},\\ m_{\mathrm{6}}=m_{{\mathrm{6}}^{\mathrm{\prime \prime }}}=\mathrm{156.393}\hspace{0.17em}\text{kg},\\ m_{{\mathrm{2}}^{\prime }}=m_{{\mathrm{4}}^{\prime }}=m_{{\mathrm{6}}^{\prime }}=\mathrm{182.512}\hspace{0.17em}\text{kg},\\ m_{l\mathrm{2}}=m_{l\mathrm{4}}=m_{r\mathrm{2}}=m_{r\mathrm{4}}=\mathrm{3500}\hspace{0.17em}\text{kg},\\ I_{\mathrm{1}}=I_{{\mathrm{1}}^{\mathrm{\prime \prime }}}=\mathrm{1.81}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}{m}^{\mathrm{4}},\hspace{1em}\hspace{1em}{I}_{\mathrm{3}}=I_{{\mathrm{3}}^{\mathrm{\prime \prime }}}=\mathrm{2.13}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}{m}^{\mathrm{4}},\\ I_{\mathrm{5}}=I_{{\mathrm{5}}^{\mathrm{\prime \prime }}}=\mathrm{2.49}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}{m}^{\mathrm{4}},\\ I_{{\mathrm{1}}^{\prime }}=I_{{\mathrm{3}}^{\prime }}=I_{{\mathrm{5}}^{\prime }}=\mathrm{1.51}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}{m}^{\mathrm{4}},\\ I_{l\mathrm{1}}=I_{r\mathrm{5}}=\mathrm{1.89}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}{m}^{\mathrm{4}},\hspace{1em}\hspace{1em}{I}_{l\mathrm{3}}=I_{r\mathrm{3}}=\mathrm{1.65}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}{m}^{\mathrm{4}}\\ I_{l\mathrm{5}}=I_{r\mathrm{1}}=\mathrm{2.52}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}{m}^{\mathrm{4}}.\end{array}$

(26)

Figure 4:

The calculation model of the portal frame structure.

To analyze the seismic response of the portal frame structure obtained by the transfer matrix method of a multibody system, the finite element software ANSYS was used to model the system to produce results for comparison with the new proposed method. When modeling, beams and columns were established by Beam188 elements, the I-type and circular tube sections were defined, the left and right columns and the inclined beams were both set as tapered members, and the sections were set as I-type sections. The beam-column joint sections were enlarged, the middle column section was set as a circular tube, and the column foot was fixed.

4.1. Vibration Characteristics of the Portal Frame Structure

The first five orders of the natural period of the portal frame structure obtained by MS-TMM and ANSYS are listed in Table 1, and the errors between these two methods are also listed.

Table 1:

Natural period T (s) of the structure obtained by the two methods.

	1st	2nd	3rd	4th	5th
MS-TMM	0.5414	0.2550	0.1375	0.1012	0.0845
ANSYS	0.5384	0.2496	0.1305	0.0925	0.0735
Error	0.58%	2.16%	5.36%	9.41%	14.87%

As listed in Table 1, the natural periods obtained by these two methods were similar, yielding a maximum error of 14.87%; this is within the limits of engineering precision. The calculation time of the MS-TMM method was approximately 1/7 that of the ANSYS method; therefore, MS-TMM greatly improved the calculation efficiency of this problem.

4.2. Analysis of the Seismic Response of the Portal Frame Structure during a Common Earthquake

4.2.1. Analysis of the Displacement for the Top of the Column of the Portal Frame Structure

The displacements for the tops of columns A and B during a common earthquake based on MS-TMM and ANSYS are shown in Figures 5 and 6, respectively.

Figure 5:

The displacement time history curves for the top of column A during a common earthquake.

Figure 6:

The displacement time history curves for the top of column B during a common earthquake.

The maximum displacements and elastic story drift angles for the tops of columns A and B during a common earthquake obtained by MS-TMM and ANSYS are listed in Tables 2 and 3, respectively, and the error between these two methods is also listed.

Table 2:

Maximum displacements for the top of column A during a common earthquake.

Methods	Maximum displacement (mm)	Time (s)	Story drift angle	Error
MS-TMM	7.3480	5.80	1/735	1.67%
ANSYS	7.2275	5.88	1/747

Table 3:

Maximum displacements for the top of column B during a common earthquake.

Methods	Maximum displacements (mm)	Time (s)	Story drift angle	Error
MS-TMM	6.9761	5.78	1/946	2.31%
ANSYS	6.8188	5.84	1/966

As listed in Tables 2 and 3, the trends of the displacement time history curves obtained by MS-TMM and ANSYS are similar. The maximum elastic story drift angle also did not exceed 1/250; this indicates that the portal frame structure satisfies the requirement of the Code for Seismic Design of Buildings of China. The maximum displacements for the top of columns obtained by these two methods are similar and show a maximum error of 2.31%, which is in agreement with typical engineering precision. The calculation time based on MS-TMM was approximately 1/15 that of the ANSYS method, reflecting the high calculation efficiency of MS-TMM.

In addition, the shapes of the displacement history curves for the tops of columns A and B are similar and almost reach their respective maximums at the same time; the latter value is slightly smaller than the former, and their ratio is approximately 0.95, which indicates that the deformation of beam is primarily an integral translation caused by the horizontal seismic force.

4.2.2. Analysis of the Shear Force on the Bottom of the Column of the Portal Frame Structure

The shear forces on the bottom of columns A and B during a common earthquake based on MS-TMM and ANSYS are shown in Figures 7 and 8, respectively.

Figure 7:

The shear force time history curves for the bottom of column A during a common earthquake.

Figure 8:

The shear force time history curves for the bottom of column B during a common earthquake.

The maximum shear forces on the bottom of column A and B during a common earthquake obtained by MS-TMM and ANSYS are listed in Tables 4 and 5, respectively, and the error between these two methods is also listed.

Table 4:

Maximum shear forces on the bottom of column A during a common earthquake.

Methods	Maximum shear forces (kN)	Time (s)	Error
MS-TMM	6.3140	5.78	2.87%
ANSYS	6.5004	4.50

Table 5:

Maximum shear forces on the bottom of column B during a common earthquake.

Methods	Maximum shear forces (kN)	Time (s)	Error
MS-TMM	6.3695	2.76	2.17%
ANSYS	6.5105	2.72

As listed in Tables 4 and 5, the trends of the shear force time history curves obtained by MS-TMM and ANSYS are similar, their maximum values are similar, and the maximum error is a low 2.87%. The calculation time of the MS-TMM method was approximately 1/15 that of the ANSYS method, showing a significant improvement in calculation efficiency. However, the shear force on the bottom of column B is shown to be slightly larger than that of column A.

4.3. Analysis of the Seismic Response of the Portal Frame Structure during a Rare Earthquake

4.3.1. Elastic-Plastic Restoring Force Model of the Structure

Because the structure exhibits elastic-plasticity during rare earthquakes, the stiffness varies with time; thus, it is important to establish the relations between the restoring force and the displacement. According to the commonly used restoring force models, the model of the bilinear restoring force is selected for this paper.

As shown in Figure 9, two broken lines are used for this bilinear restoring force instead of positive and reverse skeleton curves; the fold point is used as the yielding point, and the stiffness is not degraded but rather kept at the initial stiffness k1 during the unloading process.

Figure 9:

Bilinear restoring force model.

The bilinear restoring force model has the following features.

1. Two broken lines are used to represent the positive and reverse skeleton curves;

2. the unloading stiffness is equal to the elastic stiffness k1 and does not degrade;

3. the positive and reverse skeleton broken lines are antisymmetric.

The bilinear model curve can be determined according to the following parameters:

1. the elastic stiffness k1,

2. the yielding load F _y ,

3. the stiffness k2 after yielding.

For the steel structure in this paper, k1 = 200 GPa, F _y = 235 MPa, and k₂ = 6100 MPa.

4.3.2. Analysis of the Displacement of the Top of the Column of the Portal Frame Structure

The displacements for the tops of columns A and B during a rare earthquake based on MS-TMM and ANSYS are shown in Figures 10 and 11, respectively.

Figure 10:

Displacement time history curves for the top of column A during a rare earthquake.

Figure 11:

Displacement time history curves for the top of column B during a rare earthquake.

The maximum displacements and elastic-plastic story drift angles for the tops of columns A and B during a rare earthquake obtained by MS-TMM and ANSYS are listed in Tables 6 and 7, respectively, and the errors between these two methods are also listed.

Table 6:

Maximum displacements for the top of column A during a rare earthquake.

Methods	Maximum displacements (mm)	Time (s)	Story drift angle	Error
MS-TMM	48.77	5.58	1/111	9.69%
ANSYS	54.00	6.72	1/100

Table 7:

Maximum displacements for the top of column B during a rare earthquake.

Methods	Maximum displacements (mm)	Time (s)	Story drift angle	Error
MS-TMM	47.00	5.58	1/140	4.08%
ANSYS	49.00	6.72	1/135

As listed in Tables 6 and 7, the trends of the displacement time history curves obtained by MS-TMM and ANSYS are similar, and the maximum elastic-plastic story drift angles are approximately 1/100 in both methods, which meet the requirements of seismic resistance. The maximum displacements for the tops of the columns obtained by these two methods are similar with a maximum error of 9.69%, which is in accordance with the engineering precision. The calculation time of the MS-TMM method was approximately 1/30 that of the ANSYS method, which reflects the high calculation efficiency of the MS-TMM method.

Additionally, the shapes of the displacement history curves for the tops of columns A and B are similar and nearly reach their maximums at the same time: the latter value is shown to be slightly smaller than the former.

4.3.3. Analysis of the Shear Force on the Bottom of the Column of the Portal Frame Structure

The shear forces on the bottom of columns A and B during a rare earthquake based on MS-TMM and ANSYS are shown in Figures 12 and 13, respectively.

Figure 12:

Shear force time history curves for the bottom of column A during a rare earthquake.

Figure 13:

Shear force time history curves for the bottom of column B during a rare earthquake.

The maximum shear forces on the bottoms of columns A and B during a rare earthquake obtained by MS-TMM and ANSYS are listed in Tables 8 and 9, respectively, and the errors between these two methods are also listed.

Table 8:

Maximum shear forces on the bottom of column A during a rare earthquake.

Methods	Maximum shear forces (kN)	Time (s)	Error
MS-TMM	112.5	5.80	8.17%
ANSYS	104	5.32

Table 9:

Maximum shear forces on the bottom of column B during a rare earthquake.

Methods	Maximum shear forces (kN)	Time (s)	Error
MS-TMM	120	5.80	7.14%
ANSYS	112	5.32

As listed in Tables 8 and 9, the trends of the shear force time history curves obtained by MS-TMM and ANSYS are similar, the maximum values are similar, and the maximum error is 8.17%. The calculation time of the MS-TMM method was approximately 1/30 that of the ANSYS method, showing a greatly improved calculation speed.

5. Conclusion

In this paper, the transfer matrix method of a multibody system was established as a mechanical model of a portal frame structure, where its characteristics were used to derive the transfer matrices of the elements in the model and obtain the total transfer equation. As a practical example, the structure's natural period and seismic response during common and rare earthquakes were analyzed. The results from the analysis based on the MS-TMM method were then compared to those based on the traditional ANSYS method, and the following conclusions were drawn.

1. The first five order natural periods of the portal frame structure obtained by the two methods are similar, showing a maximum error of 14.87%. The calculation time of the MS-TMM method was shown to be approximately 1/7 that of the ANSYS method.

2. During a common earthquake, the maximum elastic story drift angles of the top of a column obtained by the two methods both meet the requirements of the code for seismic design. The maximum error of the displacements for the top of column was 2.31%, the maximum error of the shear force for the bottom of column was 2.63%, and the calculation time based on the MS-TMM method was approximately 1/15 that based on the ANSYS method.

3. During a rare earthquake, the maximum elastic-plastic drift angles of the top of a column obtained by the two methods both meet the requirements of the code for seismic design. The maximum error of the displacements for the top of column was 9.69%, the maximum error of the shear force for the bottom of the column was 8.17%, and the calculation time based on the MS-TMM method was approximately 1/30 that based on the ANSYS method.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.