Research Article | Open Access
Yin Zhang, Jianguo Ding, Hui Zhuang, Yu Chang, Peng Chen, Xiangxiang Zhang, Wenhao Xie, Jin Fan, "Pounding between Adjacent Frame Structures under Earthquake Excitation Based on Transfer Matrix Method of Multibody Systems", Advances in Civil Engineering, vol. 2019, Article ID 5706015, 31 pages, 2019. https://doi.org/10.1155/2019/5706015
Pounding between Adjacent Frame Structures under Earthquake Excitation Based on Transfer Matrix Method of Multibody Systems
In this paper, the case of two adjacent frame structures is studied by establishing a mechanical model based on the transfer matrix method of multibody system (MS-TMM). The transfer matrices of the related elements and total transfer equation are deduced, combining with the Hertz-damp mode. The pounding process of two adjacent frame structures is calculated by compiling the relevant MATLAB program during severe ground motions. The results of the study indicate that the maximum error of the peak pounding forces and the peak displacements at the top of the frame structure obtained by the MS-TMM and ANSYS are 6.22% and 9.86%, respectively. Comparing the calculation time by ANSYS and MS-TMM, it shows that the computation efficiency increases obviously by using the MS-TMM. The pounding mainly occurs at the top of the short structure; meanwhile, multiple pounding at the same time may occur when the separation gap is small. The parametric investigation has led to the conclusion that the pounding force, the number of poundings, the moment of pounding, and the structural displacement are sensitive to the change of the seismic peak acceleration and the separation gap size.
Pounding will occur when the relative displacement of adjacent buildings is greater than the width of their separation gaps under the excitation of earthquake; it will directly affect the failure mode and the degree of damage of the structure. During the 1985 Mexico City earthquake, about 40% of damaged structures were subjected to pounding and about 15% of buildings collapsed due to collision . In the 1989 Loma Prima earthquake, there were about 200 impact events in San Francisco, Oakland, Santa Cruz, and Watsonville involving more than 500 buildings . In the 1995 Hanshin earthquake, the 2008 Wenchuan earthquake, and the 2010 Yushu earthquake, a considerable part of the damage caused by structural pounding was discovered [3–5]. In the past, many countries did not specify the setting of the separation gap that results in the distance between adjacent buildings being very close to or even zero in many buildings, which may lead buildings collide with each other under the excitation of earthquakes. At present, most countries have established regulations on separation gaps, but the pounding may occur in the event of a rare earthquake . Therefore, the pounding of structures during the earthquake has attracted more and more attention from earthquake-resistant workers, and they have carried out many related studies .
Nowadays, there are two main methods for structural pounding research: the classical contact method and the contact force-based method. The classical method is based on the momentum conservation of the system and the Newtonian velocity recovery coefficient . Papadrakakis et al. , Desroches and Muthukumar , and Mahmoud et al.  analyzed the structural pounding problem based on the classical method; the application of this method is greatly limited because it cannot reflect the impact factors such as impact force and impact deformation and is not easy to combine with the existing software. An analytical technique based on the contact force-based method is developed, where the contact element is activated when the structures come into contact. The contact element uses an equivalent spring element and an equivalent damper element to simulate the interaction and energy dissipation during the collision which is placed in the event of a pounding and is withdrawn when disengaged. Scholars have conducted extensive researches on the contact force-based method, such as the linear spring pounding model which uses only one spring element to simulate structural pounding . In order to reflect the nonlinear process of pounding, the Hertz model was used to simulate structural pounding responses [13, 14]. In order to further accurately simulate the structural collision response, the Kelvin model  and the Hertz-damp model [16, 17], which can consider the energy consumption of pounding, have been proposed and applied.
A number of researchers have studied the problem of pounding for adjacent structural under earthquake: Efraimiadou et al. [18, 19] performed seismic pounding response analysis of the different layer height structures. Naderpour et al.  studied the case of pounding between two adjacent buildings by the application of single degree-of-freedom structural model. Ghandil and Aldaikh  developed a series of SSSI models to accurately study the problem of SSSI-included pounding of two adjacent buildings. Karayannis and Naoum  investigated the influence of two adjacent structures with different stories, different layouts, and initial distance on torsion under earthquake. Furthermore, some more recent numerical analyses have been carried out to study the influence of different parameters in pounding of buildings [23, 24].
The multibody system is a system in which a number of rigid and flexible bodies are connected in some way. The current various multibody system dynamics methods have the following common features: it is necessary to establish the global dynamics equations of the system; the global dynamics equation of the complex system involves high-order matrices and makes the computational workload large. Rui et al. established the transfer matrix method of multibody system (MS-TMM) by combining transfer matrix method with modern calculation method in 1993 , which has the advantages as follows: without the system global dynamics equations, high programming, low order of system matrix, and high computational efficiency . The MS-TMM was mainly divided into the transfer matrix method for linear multibody systems  and discrete time transfer matrix method . The discrete time transfer matrix method is suitable for linear time-varying, nonlinear, large-motion, and general multibody systems. In the field of civil engineering, some applications have been carried out on the application of the MS-TMM. For example, Ding et al. applied this method to the vibration analysis of building structures and the dynamic response analysis of structures under earthquake action. They have studied new single-story frame bent structures , portal frame structures , frame structures , reinforced concrete shear wall structures , etc. The results show that the MS-TMM calculation results are similar to the finite element calculation results, and the calculation efficiency is significantly improved.
Many scholars have studied the pounding problem of adjacent structures under earthquakes. However, most of them use finite element software for simulation calculation, such as ANSYS, ABAQUS, LUSAS, and DRAIN-2DX, which has a large computational workload and is time-consuming. Therefore, it is an important research direction to seek a highly efficient calculation method. In this paper, the MS-TMM is introduced into the study of pounding between two adjacent frame structures under the excitation of earthquake. First, the appropriate pounding model is selected. Next, the mechanical model is established based on MS-TMM. Then, the transfer matrices of elements and the total transfer equation of the structure are derived. Finally, the corresponding MATLAB program is compiled to analyze the influence of the separation gap size and the peak seismic acceleration on pounding process of two unequal adjacent frame structures during severe ground motions. Meanwhile, the results calculated by MS-TMM and ANSYS are compared further.
2. Pounding Model Selection
Based on the contact force-based method, scholars have conducted extensive research on structural pounding problems; the structural pounding analysis model is shown in Figure 1, where and represent the masses of the two colliding bodies, and are the corresponding displacements, and are the speeds of the initial contact moments of the two colliding bodies, respectively, and are the corresponding speeds of the colliding bodies at separation moment, k represents the stiffness coefficient of the contact unit, c is the damping coefficient, and is the initial gap.
When different pounding models are used to simulate the structural pounding process, the force-deformation relationship expressions of the contact elements are as follows.(1)Linear elastic model :where is the pounding impact force and is the relative deformation of two colliding bodies during contact.(2)Hertz model [13, 14]:(3)Kelvin model :where is the relative deformation speed of the collision body during the pounding, ζ is the corresponding damping ratio, and are the derivatives of m1 and m2 for time t, and e is the Newtonian speed recovery coefficient before and after pounding.(4)Hertz-damp model [16, 17]:where is the hysteresis damping coefficient, is the contact stiffness, reinforced concrete is usually taken as , and e represents the energy recovery coefficient, which is usually taken as 0.65 in the typical concrete structure.
The linear spring model cannot simulate the energy dissipation and the changes in local stiffness with compression. The Kelvin model also cannot represent changes in the compression stiffness of the pounding, but it can represent the energy dissipation; however, when pounding is from the maximum compression to the uncompressed regression, the Kelvin model’s simulated pounding force appears as a tensile force for a period of time before the movement is about to come out of contact, which is inconsistent with the facts. The Hertz model can simulate changes in compression stiffness, but does not represent energy dissipation during pounding. The Hertz-damp model improves this shortcoming of the Hertz model; therefore, this paper selects the Hertz-damp model.
3. Mechanical Model of a Frame Structure
In this paper, a frame structure is used as a multibody system. Since the planar arrangement of the frame structure is generally regular and symmetrical, it is simplified into a planar structure for calculation and analysis. Model simplification concerns three parts: the first part is the simplification of the column, based on the discretization idea, the column is equivalent to several concentrated masses connected by vertical elastic beams; the second part is the beam-column joint, which is equivalent to a plane hinge unit; and the third part is the connection of the plate and the beam, which is simplified as a transverse elastic beam connected to the concentrated mass. The self-weight and external load of the structure are applied to the concentrated mass, ignoring the axial deformation of the elastic beam, only considering the lateral deformation. In summary, the mechanical model of the frame structure is shown in Figure 2.
4. Derivation of Element Transfer Matrix
In this paper, the discrete time transfer matrix method is used to analyze the pounding response of two structures under the excitation of earthquake. For time-varying systems, there is no linear relationship between state vectors; so, we need to introduce a linearization idea to linearize the state vector in order to transfer the force and displacement between points. When the time step is small to a certain extent, the relationship between physical parameters can be approximated, seen as linear in the physical process corresponding to each time step. In this paper, the linearization of acceleration and velocity is performed by the Newmark-β method ; the relationship between velocity, acceleration, and displacement is as follows:where
When calculating the pounding response of two frame structures under earthquake, it is necessary to consider the seismic excitation and pounding force; therefore, the concept of extending the transfer matrix is introduced here to distinguish the transfer matrix based on the linear multibody system transfer matrix method and the multibody system discrete time transfer matrix method. The state vector of elements in physical coordinates is defined aswhere x is the displacement in x direction, y is the displacement in y direction, θz is the angular displacement around z, mz is the internal moment, qx is the internal force in x direction, and qy is the internal force in y direction. The external force is , where represents the pounding force, and it is assumed in this paper that the pounding only occurs at the concentrated mass point.
4.1. Concentrated Mass at One Input End and One Output End
The concentrated mass of mass m is shown in Figure 3(a), the input end is I and the output end is O. According to its stress balance and deformation relationship, we have
4.2. Concentrated Mass at One Input End and Two Output Ends
The concentrated mass of mass m is shown in Figure 3(b), the input end is I and the output ends are O1 and O2. According to its stress balance and deformation relationship, we have
According to the transfer equation, the state vector of the output can be obtained as follows:where
4.3. Concentrated Mass at Two Input Ends and Two Output Ends
The concentrated mass of mass m is shown in Figure 3(c), the input ends are I1 and I2 and the output ends are O1 and O2. According to its stress balance and deformation relationship, we have
According to the transfer equation, the state vector of the output can be obtained as follows:
4.4. Concentrated Mass at Two Input Ends and One Output End
The concentrated mass of mass m is shown in Figure 3(d), the input ends are I1 and I2 and the output end is O. According to its stress balance and deformation relationship, we have
4.5. Massless Elastic Beam
The longitudinal elastic beam with a length of l and a flexural rigidity of k is shown in Figure 3(e), the input end is I and the output end is O. According to its stress balance and deformation relationship, we have
The transfer equation is ; combining equation (20), we can get the transfer matrix as follows:
Similarly, the transfer matrix of transverse massless beams is
4.6. Planar Elastic Hinge
A plane elastic joint with angular stiffness of , lateral spring stiffness of , and longitudinal spring stiffness of is shown in Figure 3(f), and the transfer equation iswhere
5. Derivation of the Total Transfer Equation of Frame Structure
The transmission direction in this paper is shown in Figure 4. The bottom ends of the frame column are the input ends, and is the output end. The mechanical element model number description is shown in the figure: the bottom ends of the frame columns are numbered as from left to right; the mass points from the left to the right of the first beam-column joint are numbered as ; the mass points from the left to the right of the second beam-column joint are numbered as , and so on; the mass point of the beam-column joint of the nth row of the mth layer is numbered as ; the first layer of beams is numbered as from left to right, and so on; the number of the ()th beam of the mth layer is . Take the number beam and number column as examples to illustrate the internal numbering of the beam and the inside of the column, as shown in Figure 5.
The state vector that defines the elements in physical coordinates is
Taking frame column and beam as examples, the transfer matrices are derived as follows:
For the side cross-frame column , input from the and passed up, the transfer equation can be obtained as follows:
For , the input is and the outputs are , and the transfer matrices of the concentrated mass output according to one end input are
In the same way, we continue to pass up and sort out the following equation:
Equation (29) can be organized into a matrix form:where are state vectors at the left end of beams and the transfer matrix of each beam of the first span is represented by ; according to its transfer matrix and transfer equation, the state vector at the right end of each beam is
UL1 can be defined as follows:
For the intermediate frame column , pass up from input , the transfer equation of a21 is
Element a21 is a concentrated mass at one input end and two output ends, the input ends are , the output ends are , so we have
In the same way, we continue to pass up and sort outwhich is organized into a matrix formU1A2 and U1A2 can be defined as follows:where are state vectors at the left end of the beam and the transfer matrices of the second span beam are represented by ; according to the transfer matrix and the transfer equation of a span beam, the state vector of the right end of each beam can be obtained as follows:
The transfer of the state vector of the intermediate frame is the same as that of the frame column , so it can be derived in the same way. The transfer equation of the frame column is
U1Aj and U2Aj can be defined as follows:where are state vectors at the left end of beams , and using to represent the transfer matrix of the jth span of each beam, the state vector at the right end of each beam can be obtained as follows:
For the rightmost column , are the state vectors on the right side of the span beam. Passing down from , the transfer equation is
From the transfer relationship, the following equation can be obtained in turn:
For , the input end is and the output end is ; so the transfer equation is expressed as follows:
Equation (45) can be organized into a matrix form as follows:
Substituting the above derivation formula into equation (46), the relationship between the output end and the input end is
Then, the total transfer equation of the frame structure iswhere
6. Example Analysis
Two adjacent frame structures are selected as models, one of which is a 5-layer structure (J1) and the other is a 10-layer structure (J2); the height of each layer of the two structures is 3.2 m, the concrete grade is C30, the strength grade of the longitudinal reinforcement of the column and the beam is HRB400, the strength grade of the stirrup is HPB300, and the cross-sectional dimensions of the J1 and J2 structural columns are 500 mm ∗ 500 mm and 600 mm ∗ 600 mm, respectively, and the beam cross-section dimensions are 600 mm ∗ 300 mm. The structural schematic is shown in Figure 6.
A simplified model of the structure based on the MS-TMM is shown in Figure 7. J1 simplifies each frame column into four concentrated masses and four sections of massless beams; the first span of each beam is simplified into two elastic hinges, three concentrated masses, and three sections of massless elastic beams; the second span of each beam is simplified into two planar elastic hinges, seven concentrated masses, and seven sections of massless elastic beams. J2 simplifies each frame column into four concentrated masses and four sections of massless beams; simplifies each of the first and third span of each beam into two planar elastic hinges, seven concentrated masses, and seven sections of massless elastic beams; simplifies each beam of the second span into two elastic hinges, three concentrated masses, and three sections of massless elastic beams. The specific parameters are shown in Table 1.
When the peak acceleration of seismic is large, the structure may exhibit elastoplastic deformation. Therefore, the stiffness and damping of the structural material should be considered as a function of time, that is, determines the structural restoring force model. Currently used resilience models include bilinear restoring model and trilinear restoring model. The bilinear restoring model is too simple and rough; so, this paper uses the trilinear restoring model, as shown in Figure 8.
The steps for calculating the pounding response of two structures under earthquake based on MS-TMM are shown in Figure 9. In this paper, the calculation is achieved by compiling the relevant MATLAB program.
In this paper, two natural earthquake waves (the El Centro earthquake wave and the Taft earthquake wave) and one artificial earthquake wave (the Nanjing earthquake wave) are chosen. We adjust the peak acceleration to 35 cm/s2, 70 cm/s2, and 220 cm/s2, respectively, and the pounding response of the structure under different separation gap sizes is calculated, as shown in Table 2. In order to compare and analyze the pounding response of adjacent structures under the excitation of earthquake based on MS-TMM, this paper also uses the element software ANSYS to perform modeling and calculation. In the analysis of ANSYS, because the structure is regular, BEAM161 is selected as the model for both the beam and the column elements, and SHELL63 is selected as the model for the floor elements, the contact type is defined as ASSC, as shown in Figure 10.
7. Pounding Force Analysis
The time history of the pounding force based on MS-TMM and ANSYS is shown in Figures 11 to 14; it can be seen that the results obtained by the two methods are similar. It is found that changing the peak value of seismic acceleration and the width of separation gap will have a great influence on pounding force, the number of pounding, and the moment of pounding. Under the earthquake with same peak acceleration, as the width of the separation gap decreases, the pounding force and the number of poundings increase. When the separation gap width is the same, as the peak acceleration of the earthquake increases, the pounding force and the number of poundings increase. It can be seen from Figure 11 that when the acceleration peak is 35 cm/s2 and the separation gap width is 30 mm, J1, and J2, will not pound. As shown in Figure 12, when the peak acceleration is 70 cm/s2 and the width of the gap is 60 mm, there will be no pounding between J1, and J2, . As shown in Figure 13, when the peak acceleration is 220 cm/s2 and the width of the gap is 150 mm, there will be no pounding between J1, and J2, . It can be seen from Figure 14 that when the acceleration peak is 400 cm/s2 and when the gap size is 220 mm, J1, and J2, will not pound.
According to Table 3, we can see, under the excitation of the El Centro wave with a peak acceleration of 35 cm/m2, when the widths of the separation gaps are 10 mm and 20 mm, respectively, the errors of the peak value of the pounding force calculated by the two methods are 5.59% and 3.88%, respectively. Under the excitation of the El Centro wave with a peak acceleration of 70 cm/m2, when the widths of the separation gaps are 10 mm and 40 mm, respectively, the errors of the peak value of the pounding force calculated by the two methods are 6.22% and 4.62%, respectively. Under the excitation of the El Centro wave with a peak acceleration of 220 cm/m2, when the widths of the separation gaps are 100 mm and 130 mm, respectively, the errors of the peak value of the pounding force calculated by the two methods are 5.37% and 4.39%, respectively. Under the excitation of the El Centro wave with a peak acceleration of 400 cm/m2, when the widths of the separation gap are 200 mm and 210 mm, respectively, the errors of the peak value of the pounding force calculated by the two methods are 6.25% and 5.32%, respectively. Evidently, the results of the two methods are almost the same; however, using ANSYS to calculate the pounding response of two adjacent frame structures under the excitation of earthquake takes about 8 hours, while the calculation by MS-TMM only takes about 20 minutes; its time consumption is only 1/24 of ANSYS, and the calculation speed advantage is very obvious.
8. Displacement Analysis
Select the point of J1, and the point J2, for study. The displacement time-history diagrams are shown in Figures 15 to 22. It can be seen from the comparison that the trends of the displacement time-history diagram obtained by the two methods are almost the same. As illustrated in the figure, the acceleration peak of seismic wave and the separation gap both have an effect on the structural displacement response. The main feature is that the displacement on the pounding side is limited, which makes the displacement reaction produce obvious directional difference, and the difference becomes more significant with the increase of the acceleration peak and the decrease of the separation gap. Statistical comparisons are made between the displacement peaks obtained by the two methods, as shown in Tables 4 and 5.