Dynamic analysis of multilayer-reinforced concrete frame structures based on NewMark-β method


 The analysis method of the simplified structure formation model provides the basis for the analysis of the reinforced concrete (RC) structure under earthquake and dynamic load, which has important significance for seismic analysis of RC structure. In this paper, the three-layer RC frame structure is simulated and analyzed by MATLAB based on the NewMark-β method, considering the influence of time-varying simple harmonic loads and seismic waves on acceleration, displacement, and velocity of RC structure. The vibration response of the RC structure is analyzed by introducing the stiffness reduction coefficient. The results show that NewMark-β method provides a new idea for the seismic response of RC frame structures, making the seismic analysis of frame structures more practical; the variation range of its elastic modulus is obtained through the analysis of the constitutive model of RC, which provides the basis for the value of the stiffness coefficient; the application of the top load and the bottom load has different structural responses to the RC frame structure, and the impact of the load on the structure is more adverse when the load acts on the bottom; with the change of time, the binding stiffness coefficient will also change, and the stability of the structure will decrease greatly; the function relationship between the acceleration of the third floor and the reduction coefficient of rigidity is obtained by taking different values of the reduction coefficient of rigidity.


Introduction
The progress of the construction industry is particularly prominent [1]. Building forms are changing with each passing day, especially many high-rise buildings. In addition to the height increase, the plane layout, and vertical shape being more and more complex, the structural system is also increasingly diversified, which puts forward higher requirements for the seismic design of building structures [2,3]. The purpose of seismic design is to make the building have the corresponding resistance ability to the earthquake of different frequency and intensity during the service life [4]. The distribution of structural stiffness is accompanied with certain unevenness. Therefore, to a large extent, it will lead to the overall deformation of the building, and in serious cases, it can cause local cracking of the building. Therefore, the bottom frame structure is difficult to be used effectively in some areas with high fortification intensity. In order to further improve the seismic design value of building structures, the influence of structural stiffness must be considered [5].
Unlike linear analysis, incremental iterative formulas for structural nonlinear analysis have not been developed to the extent that every researcher follows exactly the same logical steps [6]. In most cases, only the Newton-Raphson method or other algorithms are used as tools for incremental iterative analysis [7]. However, there are many varieties in each stage of the calculation structure, especially in the calculation of the element node forces (referred to the correction stage), which will affect the accuracy of the calculation results and the convergence of the iterative process [8]. For the NewMark-β method, when the control parameters α = 1/2 and β = 1/4, the mean constant acceleration method has second-order accuracy (to meet the engineering requirements). NewMark-β is unconditionally stable and is widely used. By controlling the NewMark-β, the results can be more accurate and converging, and the velocity and acceleration responses can be calculated more easily than other methods. Cheynet [9] used the nonstationary Kanai-Tajimi model to record the ground acceleration in MATLAB, simulating the idealized model of the frame structure. La Spina [10] investigated the seismic input single-degree of freedom (SDOF) structure and found that the acceleration changed linearly with the step. The response spectrum is used to analyze the response of the structure linear system and seismic input, without considering the nonlinear factors. Three constitutive models for nonlinear analysis of structures were proposed, as follows: linear elasto-elasticity, complete plasticity, and Armstrong Frederick cyclic hardening plasticity. Al-Rumaithi [11] calculated the dynamic response multi-degree-of-freedom nonclassical damping linear system by MDOF based on idealized shear stiffness matrix construction and assuming that all floors have the quality. It provides the seismic acceleration time history of the shock response spectrum program, which can be used to calculate peak ground acceleration and velocity based on the known shock response spectrum of the acceleration time history of the earthquake. Tazarv [12] investigated the response of linear SDOF earthquakes to ground motion through NewMark method. The response of linear single degree of freedom to seismic ground motion provides a simplified calculation method, which is easy to be implemented by MATLAB program. The nonlinear (or linear) response of a single degree of freedom damped mass-spring system under external force was predicted. The dynamic analysis of multi-structure underground excitation was presented. However, the idealized model based on the foregoing does not consider the effect of structural stiffness reduction and it is not used by considering the RC structure. Geometrical nonlinearity of materials is an important characteristic of RC structures. The design of RC structures is generally based on elastic theory. However, the stiffness coefficient of RC structure was not constant due to the plastic development of RC. Due to the effect of seismic force, the response of structural stiffness reduction became more and more obvious. Hu et al. [13] investigated the influence of nanomaterials on reinforcement plasticity, which provided a reference for the research on structural stiffness. Meng et al. [14] investigated the influence of nano-strengthening on the properties and microstructure of recycled concrete, providing another idea for the study of concrete stiffness. The review of the existing literature showed that the scholars have provided different assumptions on the seismic modes and structural mode shapes. He conducted reasonable numerical simulations based on these assumptions and proposed different ideas for structural seismic modal analysis. But none of the above researches considered time factor caused by reinforced concrete (RC) structural stiffness reduction. Therefore, the objective of this study is to introduce the stiffness reduction according to the empirical formula, compare the displacement, velocity, and acceleration responses before and after the stiffness reduction is applied at different positions, and obtain the functional relationship between the stiffness reduction factor and the acceleration of the third floor of the RC frame structure.
Only when the constitutive models of RC structures of different materials are well-mastered can the stiffness coefficient be analyzed effectively. It is of great significance for different scholars to study the material characteristics through experiments [15,16]. According to the different theoretical basis of mechanics, the existing constitutive models can be roughly divided into the following types: linear elastic and nonlinear elastic constitutive models based on elastic theory; elastoplastic and elastoplastic hardening constitutive models based on the classical plasticity theory; and the constitutive model of concrete described by inner time theory. Later Ni et al. also established a series of constitutive models under different conditions through numerical simulation and experiments [17]. Zhan et al. [18] analyzed the stiffness and deflection of RC members after cracking and established constitutive models at different stages. Vecchio and Collins [19] conducted shear and tension-compression experiments on RC plates and obtained the stressstrain relationship of the corresponding RC plates. Wang and Zhou [20] used the hypothesis of strain coordination and strength equivalence to construct the RC structure plastic damage model. This model can be used in the ABAQUS finite element damage plasticity numerical analysis of RC wall structures and provides another way of thinking for the analysis of RC wall structure seismic damage model. Zhang et al. [21] used ABAQUS to conduct finite element simulation in mass concrete and obtained the nonlinear constitutive relation and failure criterion for concrete. Fang and Sun [22] proposed a three-dimensional constitutive model for concrete and obtained good failure characteristics of RC structures in accordance with the model through experiments. Zhan et al. [23] used the William-Wanke plastic constitutive model to describe the nonlinear behavior of concrete, taking into account the interaction between steel bar and concrete, and verified it through experiments. Desai et al. [24] and Desai [25] classified the RC structure constitutive models, applying the elastoplastic model of the flow law, and took the softening behavior into consideration and used the damage model to correct it. Zhang et al. [26] conducted numerical simulation on RC beams, taking into account the influence of stiffness degradation during loading and unloading of RC structures, and verified its constitutive model. Shang and Zhuang [27] simulated the hysteretic relationship of RC columns based on the uniaxial tension and compression of concrete provided by China Code for Design of Concrete Structures [28]. These constitutive relations can well reflect the variation trend of elastic modulus in the loading process, so that the change rate of stiffness can be analyzed [28]. The analysis of the constitutive model provides a reasonable idea for the reduction of the stiffness coefficient and a reference for the reasonable value of the stiffness coefficient. The process of flowchart is shown in graphical abstract.

Introduction to NewMark-β method
NewMark-β method is a method to unify the linear acceleration method, being widely used in finite element analysis. Acceleration is modified by NewMark-β method, as follows: If α = 1/2, only keep β, which is the NewMark-β method. If α = 1/2, β = 1/6, that is equivalent to the acceleration varying linearly over t Δ , which is the Wilson-θ method with θ = 1; if the acceleration is constant in t Δ , it is the average acceleration method [29].
It can be obtained from Equations (1) and (2): Substitute the Equations (3) and (4) into the dynamic equation at + t t Δ : Where, The accuracy of the NewMark-β calculation results is mainly determined by the time step t Δ . The determination of t Δ needs to consider the load change and the length of the natural vibration period T. In general, it is required that t Δ be less than 1/7 of the natural vibration period of the minimum structure that has an important effect on the response. When ≥ β 1 4 , NewMark-β method is unconditionally stable; when < β 1 4 , it is conditionally stable [30].
Considering the vibration differential equation at time + t t Δ , it can be obtained as: Substituting Equation (10) into the structural dynamic equation： and Select the integral step size t, parameter α and β, and calculate the integral constant: Forming effective stiffness matrix: (2) Calculating the time step: Calculating the effective load at time + t t Δ : Getting the displacement at time + t t Δ : Calculating the velocity and acceleration at time + t t Δ :

Analysis and discussion
Taking the multilayer frame structure as an example for analysis, the structural calculation model is shown in Figure 1 and the stiffness information of each floor is shown in Table 1.

Structural stiffness reduction of RC structure
Stiffness degradation occurs in plastic hinge areas of RC frame structures under repeated loads. The accumulation of structural fatigue will increase the plastic zone of the structure and affect the resistance and dynamic characteristics of the structure, resulting in local or overall performance loss [32]. The experiment of high strength steel provides an idea of reinforcing steel, and the combination effect with concrete is more obvious [33]. Scholars have proposed some models, which well-reflected the force characteristics of stiffness degrading components, and obtained some hysteretic cycle rules, such as IMK hysteretic cycle rule. The MODLMK bending moment Angle model proposed through OpenSees simulation experiment well-reflected the relationship between column top force and displacement, which validates the effectiveness and stability of the model [34]. The stiffness reduction coefficient of the scaffold column took into account an overall comprehensive stiffness reduction coefficient. In the elastic second-order analysis, the overall displacement of the structure and the horizontal displacement between layers were approximately the same as the results of the nonlinear analysis [35]. In American (ACI 318-14) code [36], the structural effect was calculated, and the reduction coefficient of beam member was 0.35 and column member was 0.7. The cracked wall was 0.35 and plate was 0.25. In New Zealand (NZS3101) code [37], the reduction coefficient of beam members was different due to different sections. The rectangular section was 0.4, and the T-type and L-type components were 0.35. Column components are valued according to different axial compression ratios. China Code for Design of Concrete Structures GB50010-2010 [38] also considered the influence of stiffness reduction, and the reduction coefficient of beam members was 0.4, column members  were 0.6, and wall members were 0.45. In this paper, a harmonic load calculation example was applied. The calculation example is considered as the stiffness coefficient is 0.75. During the nonlinear analysis, a stiffness reduction coefficient of 0.75 was introduced to reduce the elastic stiffness K 0 to K 0 ⁎ 0.75, so that under the original load level of the structure, the elastic analysis and nonlinear analysis could produce approximately the same structural response [39]. The internal forces and deformations of the structure under nonlinear conditions were simulated by a unified stiffness reduction factor, rather than giving different reduction factors for different stiffness degradation sections. The unified reduction factor was called the comprehensive stiffness reduction factor. The analysis of the constitutive model provides a reasonable way to reduce the stiffness coefficient [40]. Stiffness reduction coefficient was introduced to analyze the structural acceleration response with different reduction coefficients of stiffness α in the ranges of (0.5, 1). By analyzing the experimental results made by Zhan et al. [18], the constitutive relation of RC is derived, and the variation trend of elastic modulus in the loading process is obtained by fitting the constitutive relation of RC beam, as shown in Figures 2 and 3. It is found that the variation of stiffness coefficient keeps the descending section and has obvious polynomial characteristics with the gradient of variation at 0.05. So that under the original load level of the structure, the elastic analysis and nonlinear analysis produced approximately the same structural response [41]. Using this principle, the same nonlinear constitutive model S in the structure of different input levels corresponded to different reduction factor, but with the engineering practicability, through a universal applicable reduction factor. However, in view of engineering practicability, through a universally applicable reduction factor, it was used in the elastic second-order analysis method. The elastic second-order method and the nonlinear finite element method are used to calculate the interlayer displacement Angle and interlayer displacement. The results were equivalent, instead of just constraining the stiffness reduction coefficient α of a certain state [42]. In addition, carbon nanotube cement mortar also improves the workability of concrete and has important research value in mitigating stiffness degradation response [43].

Structural response analysis
The three-layer RC frame structure was subjected to a harmonic load at the top or bottom for 5 s, and the response time was 100 s, which was completed in 2,000 steps. The damping matrix was constructed by Rayleigh damping. By comparing the loads applied at different positions before and after the reduction of stiffness, MATLAB was used to simulate and analyze the vibration shapes under different conditions. The stiffness matrix of the structure was constructed, and the beam element was analyzed in the local coordinate system. The material parameters (including elastic modulus E = 30000 Mpa, element length L = 1 m, moment of inertia of Section I = 1.33 × 10 8 mm 4 and Section A = 40000 mm 2 ) were known, and the stiffness matrix of the beam element in the local coordinate system was obtained, which could be converted into the stiffness matrix in the global coordinate system. As shown in the Figure 1  , the harmonic load F = 50 ⁎ sin(0.8 ⁎ π ⁎ t) was applied at the top or bottom, the force action time was 5 s, the response time was 100 s, and the step size was 2,000 steps. The damping matrix was constructed by Rayleigh damping. The stiffness reduction coefficient was 0.75.
(1) A harmonic load of F = 50 ⁎ sin(0.8 ⁎ π ⁎ t) was applied to the top of the frame structure, and the structural response is shown in Figure 4. The responses of displacement, velocity, and acceleration before and after stiffness reduction were compared under the condition of harmonic load applied to the top of three-layer frame structure. The variation range of displacement before reduction with time was within the range of 0.6. The displacement of the first layer and the third layer acted in different directions, and its displacement with time was approximate to the change of harmonic function. It reached the peak value in about 6 s, and then the fluctuation of its displacement gradually leveled off as the load stopped being applied. In the case that the reduction coefficient of stiffness was 0.75, the range of displacement that changed with time expands and changes within the range of (−1, 0.8), and the maximum absolute value of displacement was larger than that before the reduction. The displacement difference between different layers also increased. Compared with the velocity images, the velocity range before reduction varies within the range of (−1.5, 1.5), and the curve features were similar to those of the displacement images. When the reduction coefficient was 0.75, the velocity changed before (−2, 2.3), and the peak value was higher than that before the reduction. Compared with the peak acceleration before and after reduction, the second-and third-layer accelerations had similar frequencies. The acceleration was greatly increased after the reduction. To sum up, after the reduction of the load applied on the top, the stability decreased and the frequency approached due to the decline of the structure, so the response of vibration to the structure was more obvious.
(2) A harmonic load of F = 50 × sin(0.8 * π * t) was applied to the bottom of the frame structure, and the structural response is shown in Figure 5.
When the load was applied to the bottom, the displacement, velocity, and acceleration images before and after the reduction were compared. The curve differed greatly from the harmonic curve. The displacement trends of the first and second layers before reduction were the same and close to each other. After the reduction, the displacement difference of each layer increased, and the absolute value of its maximum displacement had a great change compared with that before. However, its acceleration had no significant change, with a range of about 5%. After the reduction, the speed increased, and the speed showed multiple peaks, which varied dramatically in 4-8 s. To sum up, under the action of bottom load, the displacement and velocity had a great change with time due to the reduction of stiffness, while the change trend of acceleration was not obvious and almost maintains the original shape. Stiffness reduction also had an obvious effect on the stability of the structure.
(3) The random seismic wave response was generated by MATLAB, and the geodetic acceleration curve was generated [9], as shown in Figure 6. The displacement, velocity, and acceleration responses of each layer of the structure were obtained, as shown in Figure 7.

Influences of the stiffness reduction factor leading to regression analysis
At the top of the frame structure applying F = 50 ⁎ sin(0.8 * PI * t) harmonic load, the stiffness coefficient in range (0.5, 1), each group of data of stiffness reduction factor was 0.02. Through different stiffness coefficient values, the maximum acceleration response of the third layer was obtained by 26 MATLAB nonlinear analysis. Through curve fitting, the acceleration change law of the three-layer frame structure under different stiffness reduction coefficients was fitted. The value was 1 for the first time and then the value was successively decreased by 0.02. The first ten simulation data are shown in Table 2, the 26 simulation results were shown in Figure 8, and the fitted polar coordinates were shown in Figure 9.
The change law of acceleration of three-layer frame structure under different stiffness reduction coefficients was obtained by fitting the function relationship between the maximum acceleration and the number of experiments (the interval of stiffness coefficient goes to 0.02): y = 1.8201 * ln(x) + 1.6345.
This indicates that the influence of the composition of RC materials on the stiffness of the structural analysis should be considered. In addition, nanoparticles in cement have an important effect on cement-based composites, which have an impact on the properties of concrete [44,45]. Therefore, it is necessary to consider the influence of nanomaterial concrete on building structural analysis in future research.

Conclusion
(1) NewMark-β method provides a new idea for the seismic response of RC frame structures, making the seismic analysis of frame structures more practical. (2) Through the analysis of the constitutive model of RC, the variation range of its elastic modulus is obtained, which provides the basis for the value of the stiffness coefficient.
The application of the top load and the bottom load has different structural responses to the RC frame structure, and the impact of the load on the structure is more adverse when the load acts on the bottom. (4) With the change of time, the binding stiffness coefficient will also change, and the stability of the structure will decrease greatly. Through the reduction of the stiffness of different coefficients, the function relationship between the acceleration of the third layer and the reduction coefficient is y = 1.8201 * ln (x) + 1.6345.   Dynamic analysis of multilayer-reinforced concrete frame structures  575