# A New Elastoplastic Time-History Analysis Method for Frame Structures

Authors: Qianqian Liang; Chen Zhao; Jun Hu
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

This study aimed to analyze the formation and application of the time-domain elastoplastic response spectrum. The elastoplastic response spectrum in the time domain was computed according to the trilinear force-restoring model. The time-domain elastoplastic response spectrum corresponded to a specific yield strength coefficient, fracture stiffness, and yield stiffness. However, the force-restoring models corresponding to different structural systems and the states of the structural systems at different moments were not the same. Therefore, the dynamic characteristics of a particular periodic point corresponding to a particular structure were meaningful for the elastoplastic response spectrum. In addition, the curve in the time-domain dimension along the periodic point truly reflected the real-time response of the structure when the structure encountered a seismic load.


---

## Body

## 1. Introduction

The elastic response spectrum provides the maximum structural response to an earthquake in the elastic stage. However, it cannot be directly used for structural seismic internal stress analysis. Therefore, understanding the time-domain elastoplastic response spectrum of the structure is necessary to understand the nonlinear performance and development history of the structure [1]. Elastoplastic response analysis of structures is generally based on displacement. At present, the commonly used elastoplastic seismic response calculation methods are the dynamic elastoplastic time-history analysis method, static elastoplastic analysis method, and simplified elastoplastic analysis method [2–6]. The abovementioned methods are widely used as traditional methods to analyze elastoplastic deformation of structures, but the simplified elastoplastic analysis method is an empirical amplification method. Because of the complex and changeable forms of structures and the various influencing factors of ground motion, it is difficult to consider the influence of different factors by only one amplification coefficient, and the accuracy of the analysis is poor. The current spectrum curve cannot well reflect the actual cyclic stress process of the structure. Therefore, this paper presents a method to solve the development process of elastoplastic seismic force of structures by means of mode decomposition method using time-domain elastoplastic response spectrum. The force-restoring model of the structure needs to be determined first to solve the time-domain elastoplastic response spectrum of the structure. Then, the elastoplastic dynamic differential equation of the structure can be derived. Finally, the structural response spectrum can be analyzed by the numerical integration method [7, 8].


## 2. Elastoplastic Force-Restoring Model

The force-restoring model is a function between internal forces and deformations obtained by mathematical induction based on numerous experiments [9, 10]. The force-restoring curve is the embodiment of the structure or component strength, stiffness, and comprehensive ductility. It also reflects the energy absorption ability of the structure or component. Thus, it is the basis for studying the elastoplastic properties of the structure [11–15]. For the reinforced concrete structure, the force-restoring performance can be understood from materials and components. The force-restoring model of the material reflects the constitutive relationship between steel and concrete. The force-restoring model of the components reflects the relationship between internal force and deformation hysteresis. Two rules are applied to select the force-restoring model of reinforced concrete structures: to ensure the accuracy of the model, which can reflect the structural stress state, and to express in a mathematical form that is able to analyze the structure using nonlinear mathematics. Recently, the most common force-restoring models for reinforced concrete structural components include the bilinear model (Figure 1) and the degenerate trilinear model (Figure 2).Figure 1
Bilinear model.Figure 2
Degenerate trilinear model.The degenerate trilinear model is determined using parameters such as fracture loadPc, yield load Py, elastic stiffness, fracture stiffness, and postyield stiffness. The stiffness of the model during unloading remains unchanged, which is always the tangential stiffness between the yield point and the displacement origin [16]. The degenerate trilinear model is more reasonable and accurate in simulating the variation in internal forces and displacements in reinforced concrete structural components in an elastoplastic state during large earthquakes.


## 3. Numerical Expression of Degenerate Trilinear Restoration Force Model

Not only the mechanical properties of structural components but also the deformation properties are included in the force-restoring curve model [17]. This study used a degenerate trilinear force-restoring model to describe the structural characteristics such as stiffness degradation and ductility. The force-restoring equation is assumed to be(1)fx=k¯x+b,where k¯is the slope of the force-restoring equation and brepresents different possible stages of the restoring force and displacement of the particle during earthquakes. Different stages represented the initial stiffness as ke, the fracture stiffness as kc, and the postyield stiffness as ky. The ratio of the fracture stiffness kc to the initial stiffness ke was denoted as α, and the ratio of the postyield stiffness ky to the fracture stiffness kc was denoted as β. The load when the structural component was fractured was denoted as Pc, and the load at the time of yield was denoted as Py.For a structural component whose initial stiffnesske, fracture stiffness kc, postyield stiffness ky, fracture load Pc, and the yielding load Py were all known, its degeneration trilinear force-restoring model was then determined. As shown in Figure 2, the strain or displacement corresponding to the fracture of the component was assumed as u1, and the strain or displacement when the component yielded was u2. The numerical expressions of different stages were given as follows:(i)
Stages 0–1:fx=kex(ii)
Stages 1–2:fx=αkex+1−α⋅Pc(iii)
Stages 2–9:fx=αβkex+1−αβPc+1−β⋅Py(iv)
Stages 0–2 and stages 0–5:fx=αkePy/Py−1−αPc⋅x(v)
Stages 3–4:fx=αkePy/Py−1−αPc⋅x+Py−αkePyu3/Py−1−αPc(vi)
Stages 7–8:fx=αβkeu3+1−αβPc+1−β⋅Py/u3−u7⋅x−u7(vii)
Stages 9–10:fx=αkePy/Py−1−αPc⋅x+Py−αkePyu9/Py−1−αPcThe reactions of the structural components in different stages can be simplified using the mathematical expressions of the aforementioned stages, and the stress state during the whole process can be expressed. In this study, the maximum seismic response of the single-degree-of-freedom system with different periods was calculated based on the aforementioned force-restoring model when calculating the elastoplastic response spectrum. In the computing procedure, the maximum elastic response of the structure was first calculated. Then, the force-restoring curve was determined using the yield strength coefficient of the structure and setting different fracture strength ratios and yield strength ratios. Thus, the corresponding reaction analysis was performed according to the recordings of different earthquake ground dynamics.


## 4. Elastoplastic Dynamic Differential Equation

The elastoplastic dynamic differential equation for a single-degree-of-freedom system during horizontal earthquakes is as follows:(2)md2xdt2+ctdxdt+fx=−m⋅d2ugtdt2,where m is the mass; c is the damping; x is the particle displacement; t is the particle displacement; ug is the ground motion;fx is the elastoplastic restoring force; dx/dtis the relative velocity of the system;d2x/dt2 is the relative acceleration of the system; and d2ugt/dt2is the ground movement acceleration.During the deformation of single-degree-of-freedom systems, the variation in restoring force is nonlinear [18]. This study used the trilinear hysteresis curve for the force-restoring model and assumed that the damping remained constant during the deformation of the system. After substituting the restoring force into the single-degree-of-freedom elastoplastic dynamic differential equation of the system, the following expression was obtained:(3)md2xdt2+cdxdt+k¯x+b=−m⋅d2ugtdt2.Dividing both sides of the equation bym yielded the following:(4)d2xdt2+cmdxdt+k¯mx+bm=−d2ugtdt2.Therefore, the dynamic differential equation of the single-degree-of-freedom system was established based on the degenerate trilinear model:(5)d2xdt2+2ξωdxdt+rω2x=−d2ugtdt2−bm,where ξ is the damping ratio; ω is the natural frequency of the single-degree-of-freedom system; and r is 1 in stages 0–1, αin stages 1–2, αβ in stages 2–9, αPy/Py−1−αPc in stages 3–4, and stages 9–10, and 1−αu1+1−βαu1+αβu3/u3−u7 in stages 7–8.


## 5. Numerical Integration Method of Elastoplastic Response Spectrum

The earthquake acceleration and structure response of the elastoplastic system are both nonlinear [19, 20], which cannot be described by an analytical solution solved directly from the dynamic differential equation, such as an elastic system. Therefore, it can only be computed by the step-by-step integration method. This study used the Newmark-β method to solve the increments based on the quasistatic equation.The Newmark-β method is used to calculate the structural dynamic response according to the assumed acceleration variation within the time increment. It is a stepwise integration method. Its advantage is that it avoids any superposition application and can be well adapted to nonlinear analysis. Owing to various ways of acceleration variation in the time increment, the basic assumptions of the Newmark-β method are as follows:(6)dui+1dt=duidt+1−βd2uidt2+βd2ui+1dt2⋅Δt,(7)ui+1=ui+duidtΔt+12−γd2uidt2+γd2ui+1dt2⋅Δt2.This study used the average acceleration method (i.e.,β is 0.5 and γ is 0.25). The acceleration was assumed as the average value during time interval ti – ti+1:(8)d2utdt2=12d2uidt2+d2ui+1dt2.Thus, the speed and displacement were as follows:(9)dui+1dt=duidt+βd2uidt2+d2ui+1dt2⋅Δt,(10)ui+1=ui+duidtΔt+γd2uidt2+d2ui+1dt2⋅Δt2.Because of a sudden change in stiffness in the force-restoring model, a cumulative error occurs during the reciprocal integral calculation if the breakpoint is not accurately processed, reducing the reliability of the computing result. The breakpoint often appears within a certain integration step. A basic assumption of the Newmark-β method is that the stiffness of the structure in each calculation section is linear. Therefore, extra efforts are required to reach the stiffness breakpoint. This study used the precise breakpoint processing method [21]. The polygonal line in the degenerate trilinear model in Figure 2 indicates two types of breakpoints in the model: the loading point (Figure 3), in which the symbol of speed does not change on both sides of the breakpoint, and the unloading point (Figure 4), in which the sign of speed changes on different sides of the breakpoint. The aforementioned characteristics were used to determine the loading point based on the increase in absolute displacement and determine the unloading point using the sign change of speed.Figure 3
Loading point (first species).Figure 4
Unloading point (second species).The elastoplastic dynamic differential equation in increments was expressed as follows:(11)m⋅d2Δudt2+c⋅dΔudt+k⋅Δu=−m⋅d2Δugdt2.Assuming that the breakpoint appeared at a certain timet + Δt0 in the (t, t + Δt) time period, the displacement at the loading breakpoint was expressed as follows:(12)ut+Δt0=ut+Δut,which was converted into the following incremental form:(13)Δut=dutdt⋅Δt0+12d2utdt2Δt02+β⋅d2Δutdt2⋅Δt02,(14)dΔutdt=d2utdt2Δt0+αd2Δutdt2Δt0.Equation (13) could be used to obtain the following expression:(15)d2Δutdt2=1β⋅ΔutΔt02−1βdut/dtΔt0+12βd2utdt2.The incremental form of the elastoplastic dynamic differential equation (11) was used to obtain the following:(16)d2Δutdt2=−d2Δugtdt2−cm⋅dΔutdt−km⋅Δut.By substituting equation (16) into equation (13), the displacement increment in t to t + △t0 period was obtained:(17)Δut=dutdt⋅Δt0+12d2utdt2⋅Δt02−β⋅d2Δugtdt2⋅Δt02−βcm⋅dΔutdt⋅Δt02−βkm⋅Δut⋅Δt02.By substituting equation (15) into equation (13), the following expression was obtained:(18)dΔutdt=d2utdt2⋅Δt0+αβ⋅ΔutΔt0−αβ⋅dutdt−α2β⋅d2utdt2⋅Δt0.Substituting equation (18) into equation (17) yielded the following expression:(19)Δut=dutdt⋅Δt0+12d2utdt2⋅Δt02−β⋅d2Δugtdt2⋅Δt02−βcm⋅d2utdt2⋅Δt0+αβ⋅ΔutΔt0−αβ⋅dutdt−α2β⋅d2utdt2⋅Δt0⋅Δt02−βkm⋅Δut⋅Δt02.Equation (19) was organized into the standard form as follows:(20)A1⋅Δt03+B1⋅Δt02+C1⋅Δt0+D1=0.The coefficients ofA1, B1, C1, and D1 in the aforementioned equation were calculated as follows:(21)A1=−d2Δugt/dt2Δt+α2β−1⋅c⋅d2ut/dt2m,B1=−12βd2utdt2+α2β⋅c⋅d2ut/dt2m−k⋅Δutm,C1=−1βut−cα⋅Δutβm,D1=−Δutβ.For the unloading breakpoint, the velocity at the breakpoint was zero:(22)dut+Δt0dt=0,which was converted into the following incremental form:(23)dutdt+dΔutdt=0.Equation (23) yielded the following:(24)dΔutdt=−dutdt,dut+Δt0/dt was expressed using the incremental Newmark-β method:(25)dut+Δt0dt=dutdt+d2utdt2⋅Δt0+α⋅d2Δutdt2⋅Δt0=0.Equation (25) yielded the following:(26)d2Δutdt2=−1α⋅Δt0dutdt+d2utdt2⋅Δt0.By substituting equations (13), (24), and (26) into equation (11), the following expression was obtained:(27)m⋅−1α⋅Δt0dutdt+d2utdt2⋅Δt0+c⋅−dutdt+k⋅dutdt⋅Δt0+12d2utdt2⋅Δt02+β⋅d2Δutdt2⋅Δt02=−m⋅d2Δugdt2.The equation (27) is organized into the standard form as follows:(28)A2⋅Δt03+B2⋅Δt02+C2⋅Δt0+D2=0.The coefficients ofA2, B2, C2, and D2 in the aforementioned equation were calculated as follows:(29)A2=km12−βα−1⋅d2utdt2,B2=d2Δugt/dt2Δt+k⋅d2ut/dt2m+kβ⋅dut/dt2m,C2=c⋅dut/dtm+d2ut/dt2α,D2=−dut/dtα.The coefficients of the loading and unloading breakpoints were related to displacement, velocity, acceleration, and structural stiffness at timet, as shown by equations (20) and (28). Therefore, four coefficient values were determined based on the state of the structure at time t and the breakpoint type. Then, the time of breakpoint was accurately solved using equations (20) and (28).For any seismic wave, the elastic absolute acceleration response spectra were first analyzed to get the maximum absolute accelerationsa of the elastic system. Then, the elastic seismic response analysis of the structure and actual reinforcement structure of each component was performed by calculating the elastic seismic response of the structure. After that, the yield strength coefficient ζ of each floor was obtained by comparing the floor shear-bearing capacity calculated using real component reinforcement and standard material strength with floor elastic seismic shear stress calculated using rare earthquake interaction (i.e., the ratio of floor elastic seismic shear stress calculated using yield shear-bearing capacity of different floors and rare earthquake seismic interaction). The corresponding yield load was obtained as follows:(30)Fy=ξ⋅m⋅sa.Using the initial stiffness of the structure and setting the ratio of the postfracture stiffness to the elastic stiffness of the structure asβ and the ratio of stiffness after yielding to fracture stiffness as p, the postfracture stiffness kc and stiffness after yielding ky were obtained as follows:(31)kc=ke⋅β,ky=kc⋅p.Using the aforementioned known structural floor yield load, elastic stiffness, and postfracture stiffness and by setting the ratio of the fracture displacement and yield displacement of the structure asα, the calculation formula of the floor displacement when the structure was fractured was obtained as follows:(32)uc=Fyke+1/α−1⋅kc.The corresponding yield displacement calculation formula is as shown in equation (33):(33)uy=ucα.The corresponding trilinear force-restoring model was determined by substituting equation (30) into equation (33). Then, the displacement, velocity, and acceleration of the system at any time under seismic wave excitation were determined by numerical integration of the elastoplastic dynamic differential equation using the Newmark-β method [22] according to the displacement u0, speed du0/dt, and acceleration d2u0/dt2 of the single-degree-of-freedom system for any seismic wave recording. The seismic wave is divided as Figure 5 to the same time interval △T, and then the full period seismic response spectrum is calculated for the time period after each segmentation. The maximum value of seismic influence coefficient of each period point in tm−1—tm time slot is am1am2 … amn. When the natural vibration period of the single-degree-of-freedom system is Tn, the maximum value of seismic influence coefficient in each time slot (0—t1, t1—t2, …, tm−1—tm) is a1na2n … amn. The maximum value of seismic influence coefficient corresponding to natural vibration period Tn in tm−1—tm time slot is amn, thus the time domain response spectrum matrix of each period point in each time period is obtained. Hence, the seismic force of the time domain and the elastoplastic displacement of the layer of weak floors of the structural system were analyzed using the aforementioned data.Figure 5
Seismic wave splitting indication.The elastoplastic acceleration response spectra were computed using the lab-made computational programs based on the earthquake recordings of the level 9.0 EW direction earthquake in East Japan AOM021 on March 11, 2011 (300-s seismic recording length; recording seismic peak acceleration amplitude was 220 cm/s2 during 7 degrees 0.1 g earthquake; damping ζ was 0.05). The calculation example used in this paper is a frame with two spans and three floors with both directions of 6 m, and floor height is 3.2 m. The frame column section is 400  mm × 400  mm, the frame beam section is 200 mm × 500 mm, the floor dead load is 2.0 kN/m2, the live load is 3.5 kN/m2, the roof dead load is 2.0 kN/m2, the live load is 0.5 kN/m2, and the line load of side beam is 6.0 kN/m). The elastic analysis in this study showed that the yield strength of the elastic factor was set as 0.8 and the ratio of fracture stiffness to initial elastic stiffness of the hysteretic curve of concrete structures was set as 0.5. The ratio of postyielding stiffness to fracture stiffness was 0.1, and the ratio of fracture displacement to yield displacement was 0.5. Figure 6 is a time-domain elastoplastic response spectrum curve generated by calculating the time-domain seismic response spectrum with a selected time interval of 10 s and a period interval of 0.02 s. The corresponding matrix comprised 600 rows and 30 columns. Figure 7 is a time-domain elastoplastic seismic impact coefficient surface with a time interval of 0.01 s and a period interval of 0.02 s. The corresponding seismic impact coefficient matrix comprised 600 rows and 30,000 columns.Figure 6
Time-domain absolute acceleration elastoplastic influence coefficient surface (time interval 10 s).Figure 7
Time-domain absolute acceleration elastoplastic influence coefficient surface (time interval 0.01 s).Figures6 and 7 show the time-domain absolute acceleration seismic influence coefficient surface generated using the same seismic acceleration record but different time intervals. A comparison of the two figures showed that the period and time of the maximum peak of seismic influence coefficient were basically the same, and the calculated peaks are all 0.38. However, more peaks appeared on the coefficient surface in the time domain due to the very short time interval selected in Figure 7. Horizontal seismic force analysis in the time domain using this response spectrum surface was accurate enough to capture the change in 0.01 s time interval, but the corresponding computation cost was quite large. For example, if the response spectrum curve generated using this seismic record was used for calculating mode decomposition calculation, the number of time-domain seismic forces under each mode reached 30,000; the total number of time-domain seismic forces reached 30,000 n (n is the mode number). Moreover, the vibration modes of 30,000 n seismic forces were recombined and calculated. Such large-scale computations were obviously quite inefficient; even using computers to perform the computation took a long time. Therefore, the time-domain response spectrum matrix generated at intervals of 1, 5, and 10 s was recommended in this study to be used to analyze the time-domain seismic force according to the stage of the project and required calculation accuracy, so as to meet the design requirements of the corresponding project stage. For example, as shown in Figure 6, the time interval is 10 s and the corresponding elastoplastic response spectrum matrix has only 30 seismic forces in each mode. Although it could not reflect the trend of seismic force and deformation in the 10 s time period, the peak value of seismic force in the 10 s time period was not lost, roughly reflecting the change in seismic force in the time domain of 300 s seismic record.The premise of studying the elastoplastic response spectrum is that the load corresponds to the yield strength of the structure system. When the load response of the particle exceeds its cracking load for the first time, the stiffness of the system will be reduced, and the natural vibration period will be correspondingly extended. Therefore, the peak acceleration corresponding to each period of the elastoplastic response spectrum curve is smaller than that of the elastic response. As the seismic vibration intensity increased, the stiffness of the system underwent a second reduction when the load response of the particle exceeded the set yield loadfy. Then, the stiffness of the system was reduced for the second time, and the corresponding natural vibration period was extended for the second time, which again reduced the acceleration response of the system subjected to a seismic dynamic load. Using EPDA (Elastoplastic Dynamic Analysis) software developed by the Chinese Academy of Architectural Sciences and using the same seismic records and structures as those of Figures 6 and 7, the peak value of elastoplastic dynamic time-history analysis is 0.39, which shows that the time-domain elastoplastic seismic response surface calculated in this paper is effective and practical. The aforementioned analysis showed that the elastoplastic response spectrum curve corresponded to a specific yield strength coefficient, fracture stiffness, and yield stiffness reduction factor. Therefore, the corresponding periodic point of the specific structure was meaningful for the elastoplastic response spectrum curve. In addition, the curve corresponding to the time-domain dimension along the periodic point truly reflected the time-dependent reaction when the specific structure encountered the seismic dynamic load. With response spectrum curves corresponding to time-domain dimensions, the development history of structural seismic force can be obtained by mode decomposition response spectrum method, and then the time-domain values of interfloor shear force and displacement can be obtained. Because the maximum elastoplastic seismic influence coefficients of interfloor shear force and displacement calculated by ordinary mode decomposition response spectrum method are adopted under each mode of vibration, while the elasticolastic seismic shear force and displacement mining in time-domain are adopted. The maximum elasticolastic seismic impact coefficient in each period is used, and the corresponding calculation accuracy is greatly improved.The maximum elastoplastic interstory displacement angle of the frame structure is 1/78. The maximum interstory displacement angle of the structure is 1/76 obtained by elastoplastic time history analysis under the same seismic wave. Therefore, the calculation results of the two methods are basically consistent. The calculation theory in this paper meets the needs of engineering practice.


## 6. Conclusions

This study analyzed the formation and application of time-domain elastoplastic seismic response spectrum in depth. First, the elastoplastic force-restoring model, the numerical expression of the degenerate trilinear force-restoring model, and the treatment of the inflection point for reinforced concrete structures were analyzed. Then, the elastoplastic response spectrum curve of any seismic wave for a certain yield strength coefficient was solved using the Newlab-β numerical integration method using Matlab according to the basic principles of the elastoplastic dynamic differential equation. The practicality and effectiveness of time-domain elastoplastic seismic response spectrum were illustrated using computation examples. The main conclusions of this study were as follows:(1)
For the time-domain elastoplastic seismic response spectrum surface generated by different time intervals using the same seismic acceleration record, the period and time of occurrence of the maximum peak of the response spectrum were basically the same. For a shorter time interval, the time-domain range displayed a more common trend, and more peaks appeared in the surface figure. However, the minimum time interval could not be less than the recording step size of the seismic acceleration.(2)
The elastoplastic response spectrum curve corresponded to the specific yield strength coefficient, fracture stiffness, and yield stiffness reduction factor. Therefore, the corresponding periodic point of the specific structure was meaningful for the elastoplastic response spectrum curve. The curve corresponding to the time-domain dimension along the periodic point truly reflected the time-dependent reaction when the specific structure encountered the seismic dynamic load.(3)
The time-domain values of shear force and displacement between floors were analyzed using the time-domain elastoplastic response spectrum matrix. The maximum elastoplastic seismic influence coefficient of each mode was used for calculating the interlayer shear force and displacement by the ordinary mode decomposition reaction spectrum method. Also, the maximum elastoplastic seismic influence coefficient in each period was used to calculate the time-domain elastoplastic seismic shear force and displacement. Hence, the corresponding computation accuracy was greatly improved. The accuracy of the new method meets the needs of engineering. The calculation efficiency is obviously higher than that of the elastoplastic dynamic time-history analysis method.



---

*Source: 8818187-2020-09-30.xml*