New Approach to Nonlinear Dynamic Analysis of Reinforced Concrete 3D Frames; An Accurate and Computational Efﬁcient Mathematical Model

to this work. Abstract: The structural engineering community often deals with the issue of inelastic incursions of the structural response. Although buildings situated in seismic regions are usually designed using elastic analysis, most encounter signiﬁcant inelastic deformations when major events occur. In general, material nonlinearity is the most important source of nonlinearity considered in the dynamic analysis commonly performed in structural design. In this paper, the recent concept of the force analogy method for 3D structures is developed and integrated into an accurate computational efﬁcient algorithmic routine for nonlinear dynamic analysis of reinforced concrete frames. Moreover, a unique straightforward mathematical model for the numerical implementation of degrading cyclic behavior of structural elements is proposed and further used to simulate the response of a 10 story reinforced concrete frame structure. A set of nonlinear dynamic analyses are performed using the proposed algorithm in order to assess the structural damage in case of different peak ground accelerations seismic recordings. The seismic structural damage is evaluated using both structural response parameters expressed in terms of displacement and energy concepts.


Introduction
Urban development through accelerated expansion of large cities, correlated with the appearance of rather daring architectural concepts, requires the use in structural design, of increasingly complex calculations, and modeling procedures.
The majority of practitioners performing in the civil engineering domain are familiar with the concept of the response spectrum method (RSM) to evaluate the structural response of the commonly designed structures. For nearly 60 years, the RSM was used to produce a maximum peak value for one force or one displacement at an unspecified point in time [1]. Consequently, in order to take into account the actual behavior, the structural response is reduced by specific factors that are correlated with well-defined geometrical and ductility requirements imposed by the design code prescriptions. The technology evolution and its current availability oblige the structural engineers to move forward and to consider more advanced analysis methods, which take into account the full nonlinear dynamic behavior of the structural elements as well as more accurate models in terms of seismic excitation. In this respect, a simple, accurate, and fast procedure to obtain the nonlinear dynamic structural response is one of the challenges and is still being a hot topic in this domain. Among the methods available to determine the nonlinear dynamic structural response [2], the force analogy method (FAM) has proved to be both accurate and numerical efficient [3].
The key idea in FAM is to retain the initial stiffness matrix and to modify the displacement vector such that the same level of restoring forces are produced. Consequently, the displacement vector is decomposed in its elastic and inelastic components. The inelastic component of the displacement is moved to the right-hand side of the equation of motion and regarded as some perturbation forces. Coupling FAM with an explicit solver to integrate the equation of motion has proved to be both accurate and numerically efficient not only for regular problems of obtaining the nonlinear structural response of isolated structures [4] but also in the more complex case of structural pounding due to earthquake action [5]. In recent papers, Safaei S. et al. [6] studied the accuracy of the FAM for nonlinear static analysis. They concluded that the results are in good agreement with those obtained by other analysis techniques commonly used for seismic performance assessment of structures. J.-T. Qu et al. [7] concluded that FAM is effective and accurate to estimate the plastic energy dissipation of multi-story and high-rise structures. They conducted several numerical simulations on different types of structures in order to compare the plastic energy dissipation of MDOF systems and their equivalent SDOF systems. Using the force analogy method, Toloue I. et al. [8] investigated the seismic performances of two plane portal frames, one modeled with Timoshenko and one using Bernoulli beam finite element. They showed that the resulting differences in structural response are not significant for engineering purposes. Moreover, Ningthoukhongjam et al. [9] used the FAM to investigate the effects of location and intensity of single-floor mass irregularity, on seismic response of mid-rise moment-resisting steel frames.
A major breakdown of the FAM has been achieved when the concept of FAM 3D was introduced [10]. The 3D case of FAM was addressed in [11] for the analysis of structural models subjected to long-duration time-histories of wind loads and for the assessment of the structural damage state. The concept of axial force-bending moment interaction surface, coupled with an efficient return mapping algorithm was used [12]. Although the main features of the FAM 3D modeling technique were introduced in [10], there are several topics that were not fully addressed. An obvious issue in [10] is that the usage of material bilinear model is arguable for modeling cyclic behavior of reinforced concrete (RC) elements. It must to be noted that no previous paper tackles the implementation of hysteresis models for 3D structures using FAM. In other papers, Li. G et al. [13] presented a general framework for seismic damage analysis of RC plane frame structures using FAM. The inelastic behavior was modeled using two types of plastic hinges in order to account for damage behavior caused by concrete cracking and material nonlinearity. Moreover, a damage index is computed and discussed. Chao et al. [14] used a modified force analogy method to simulate the nonlinear inelastic response of reinforced concrete plane structures. Beam-column elements with three different plastic mechanisms were utilized to simulate inelastic responses caused by moment and shear force. On the other hand, Gang Li et al. [15] implemented an existing shear wall macroscopic model for use in the force analogy method. The shear wall model is composed of two vertical sliding hinges and one horizontal sliding hinge that are assigned to capture the relationship of lateral deflection or rotation versus the RC shear wall force. Moreover, shear hysteresis functions are implemented and discussed.
Having said that, the main scope of this paper is to introduce a facile computational efficient algorithm for nonlinear dynamic analysis of reinforced concrete frames. The entire model is built in a non-commercial code developed in Matlab/Simulink. In order to accomplish the defined scope of this paper, the parameter inter-dependence when building such a code is also illustrated. The numerical efficiency and nonlinear insights are further evaluated through the concept of energy balance. Moreover, a damage assessment analysis is conducted using a model based on a well known damage index. At the same time, another addressed topic is the modeling of slab contribution. Consequently, a condensation technique is used to account for rigid diaphragm behavior.
The authors' main contribution is the development of a FAM 3D-based algorithm that can be used to conduct, in a simple manner, nonlinear dynamic analysis for RC momentresisting 3D frame structures. A related important contribution is that the paper is the first to present a detailed methodology for the implementation of polygonal hysteresis models in FAM 3D, which is validated against the results obtained from commercial software. Moreover, the entire algorithm is implemented in a finite element-based numerical routine developed in Matlab/Simulink.

Force Analogy Method Theoretical Background
The force analogy method distinguishes from classical nonlinear analysis methods implemented in commercial software through the fact that throughout the entire dynamic analysis, the initial stiffness matrix of the structure K 0 is used. The basic idea of the FAM is to use a change in displacement vector in order to obtain the same level of force one will obtain using a variable stiffness model. Therefore, the equation of motion considered by the FAM is: where M, C, and K 0 are the mass, damping, and stiffness matrices, respectively, u(t) is the total displacement vector, h is a influence vector, and a g (t) represents the horizontal components of the ground acceleration. The total displacement vector u(t) is split into the elastic displacementũ(t) and the plastic displacement u (t) In Equation (3) the physical nonlinearity is treated as a perturbation of the excitation forces. The last term on the right-hand side of Equation (3) K P · u (t) has the physical meaning of "the forces producing inelastic displacement in all structural coordinates".
In the civil engineering common approach, the plasticity is concentrated in zerolength hinges with moment-rotation model parameters. Although many sophisticated models that use distributed plasticity formulations have been developed, simplified models approaches proved more practical due to their relatively condensed and efficient numerical formulations. In addition, the simplified models are able to capture more effectively the relevant features with same or insignificantly lower level of accuracy [16]. Therefore considering that the beam end rotation is the only source of nonlinearity, the term K P · u (t) in Equation (3) can be replaced by unbalanced moments K P · θ (t) generated by the inelastic beam end rotation θ (t), where K P represents the member force recovery matrix. For the sake of brevity, in the FAM, the classical nonlinear differential equation is changed with a linear equation (Equation (4)) and an associated inequality (Equation (5)) which mimic the original nonlinear differential equation.
Analyzing Equation (4) it is obvious that a second equation must be added to determine the two unknown vectors u(t) and θ (t). This equation is obtained by limiting the bending moment values using the cross-section capacity functions. Moreover, the bending moment m(t) is calculated using Equation (6).
where K R represents the member restoring force matrix. The member force recovery matrix K P and member restoring force matrix K R are assembled from beam element stiffness matrix in local axis coordinates. The detailed procedure is illustrated in the authors' previous work [10].
Moreover, the second-order differential equation (Equation (3)) can be reduced to a first-order differential equation using the state-space formulation. Furthermore, in order to obtain an accurate and efficient solution when solving Equation (7), a variable time-step integrator can be used, for example, ode23 available in Matlab/Simulink. whereż

Development of Hysteresis Models
The mathematical laws commonly called hysteresis functions govern the nonlinear behavior of the structural elements. As the inelastic deformations increase, the structural response becomes more and more sensitive to the model used to account for cyclic behavior. The hysteresis relationships are highly dependent on the loading history and properties of the member and are usually determined from laboratory experimental tests.
Although many complex models have been proposed, the bilinear hysteretic function has been largely used for both reinforced concrete and steel structures. Its simplicity resides in the fact that the behavior is described by only three rules: elastic, yielding, unloading, and load reversal slopes. The main characteristic of this model is that for high amplitude deformations a large amount of energy is dissipated, while for low amplitude cycles, zero hysteresis energy is consumed. This leads to the fact that early dissipation of the seismic energy does not allow an accurate build-up of the structural deformations.
In term of reinforced concrete members, the usage of the bilinear model has the drawback of overestimating the energy dissipation capacity. Moreover, experimental tests conducted for beam-to-column connections showed that the behavior of reinforced concrete members under cyclic loading is complex and cannot be accurately described by simple bilinear hysteresis functions [17].
Although the numerical implementation of the dynamic nonlinear mathematical is usually a complex task, one of the main features that make the force analogy method attractive is the ability to account for the nonlinear behavior in a simple and direct way. This is also the case for the implementation of degrading hysteresis models that better suit the cyclic behavior of the RC structural elements. The FAM's main advantage is that it facilitates a simple path towards implementation of any polygonal hysteresis model.
The procedure for the implementation of cyclic models in FAM 3D consists of two main steps. Firstly, a mathematical hysteresis model must be defined by establishing, for every possible case, the rules that govern the behavior. Following, the classical FAM equations are modified to account for complex cyclic of the structural elements and a new matrix K d , that contains on the main diagonal the history variables defining for every plastic hinge the slopes of the branches characterizing the hysteresis model, is introduced. This matrix is utilized to determine the capacity function variation in time. Consequently, the key feature for the implementation of complex hysteresis models in the FAM is the usage of the capacity functions. They actually represent the branches of the cyclic model and are established using certain rules that define the hysteresis model. The general form of a capacity function is presented in Equation (9).
where m inter is the vector containing the conveniently chosen intercepts of the hysteresis model. Moreover, the capacity function can be used to define constrains between consecutive bending moment values (Equation (11)).
Coupling the Equations (5) and (6) leads to a system of equations in which m k+1 and θ are unknowns, u k+1 being obtained from an explicit integration scheme using previous step information u k+1 = f (u k , a g k , u k ). Moreover, the values in K d depend on the prior history and are determined using the hysteresis model (m − θ ) rules Replacing m k+1 in the first equations yields: where K d is a K R size matrix that contains on the main diagonal the history variables k br i .
k br i are the slopes of the hysteresis model branches and nph is the number of plastic hinges. A representative feature of FAM is the usage of trial bending moment values above the elastic limit which are reduced to capacity. Thus, the first step in simulating the structural nonlinear behavior is to define the trial values based on the previous step's information.
In the original approach of the FAM, the trial bending moments are calculated in each time step assuming that no plastic rotation increment ∆θ trial occurs. For the implementation of hysteresis functions a non-zero plastic rotation trial increment ∆θ trial must be introduced. Taking into account that trial bending moments and trial plastic rotation increment are determined using the last and current step data, trial K d values must also be assumed.
It is important to mention that as long as the structural model exhibits elastic behavior, every k br i value in K d is set to zero. Consequently, the ∆θ trial is a null vector and the bending moment trial values m trial k+1 , are considered to be the final values. On the other hand, after the bending capacity is exceeded, plastic deformations occur and the hysteresis rules are engaged for the determination of final m k+1 and ∆θ vectors. For every plastic hinge, the rules that define the cyclic model are implemented using the values of the K d matrix which are updated accordingly in each time step.
As stated before, the procedure for the determination of nonlinear structural response at any time step starts with an initial K d trial matrix whose values k br i describe the unloading slopes for all the plastic hinges. Next, trial plastic rotation increment and trial bending moment are determined. Furthermore, in order to update the K d matrix, for every plastic hinge, the values in m trial k+1 are compared with both the previous corresponding values in m k and the capacity function defined for that case. If the bending decreased, the unloading assumption made for k br i is correct and the trial bending moment value is set back to the previous step value. Else, if bending increases, the trial value is compared to the capacity function Equation (16). If the capacity function is exceeded, the k br i is updated to the slope of the capacity function and the trial bending value is capped to capacity (Figure 1). On the other hand, if the trial bending values are within capacity function limits, the reloading on the unloading path case is considered and therefore k br i remains unmodified and the trial bending moment value is set back to the previous step value.
Further, using the updated values of m trial k+1 and K d , the actual change in plastic rotation is calculated.
Moreover Equation (10) is used to determine the actual m k+1 . It is important to point out that in each time step, certain elements encounter plastic behavior while the rest exhibit elastic behavior. Therefore, it must be taken into account that, when determining m k+1 , it is possible, due to the bending moment redistribution, certain elements exhibit the transition from elastic to inelastic behavior. In this case, the K d corresponding values are updated while bending moments are capped to capacity. Consequently, the actual plastic rotation ∆θ and bending moments are recalculated m k+1 . This simple routine provides the framework for the nonlinear dynamic analysis of the RC frame structures.
In the next sections some well-known key concepts are brought to the reader's attention. The energy approach to structural response and cumulative damage indices are further used to demonstrate the capabilities of the new FAM 3D algorithm endowed with a hysteresis function that models the degrading cyclic behavior of the RC elements in a spatial frame structure.

Energy Approach to Structural Response
To analyze the insights of the inelastic behavior, energy formulation is used in this paper. Although the energy balance has little practical significance to most of the practitioners in the field of civil engineering, energy is a fundamental parameter for all hydrocodes used in the analysis of special problems like impact analysis, fluid-structure interaction, blast analysis, and in general in continuum mechanics. Thus, the energy balance will give an insight into the energy dissipation capacity of a civil engineering structure. The energy balance is obtained by a simple integration of the equation of motion (Equation (1)), over the path of the total displacement vector u(t). The terms given in Equation (18) are rearranged to illustrate the kinetic energy KE k , the energy dissipated by viscous damping DE k , and the input energy IE k . The third term in the left hand side of the equation of motion can be separated in accordance with the total displacement decomposition given in Equation (2) for du = dũ + du , to obtain Equation (23). The significance of the terms given in Equation (23) is the strain energy SE k and hysteresis energy, as they are detailed in Equations (24) and (25).

Damage Quantification
Worldwide there is an increasing concern regarding the resilience of the existing buildings in order to avoid major damage [18] or even collapse [19] during next significant seismic events. Therefore, damage quantification is a central subject in current structural design. In recent papers, ref. [20] investigated the cumulative damage of a 20-story highrise steel building equipped with buckling-restrained braces subjected to earthquake and wind events using Miner's damage rule, ref. [21] introduced a multiple-input convolutions neural network model in order to estimate the seismic damage buildings. Furthermore, ref. [22] studied the post seismic behavior of a reinforced concrete building which presented structural damage due to flexural and axial actions in one of its walls on the first floor. Several solutions to improve the seismic behavior were discussed. At the same time, ref. [23] proposed a novel definition of damage states for reinforced concrete members to be used for structural vulnerability assessment. Moreover, ref. [24] utilized the seismic fragility method based on incremental dynamic analysis in order to evaluate the damage levels of high-rise structures under long-period ground motions.
Furthermore, some theoretical concepts used in this paper for the implementation of structural damage assessment are presented. In order to quantify the damage induced by earthquake excitation and sustained by the structural elements through plastic deformations, a damage index approach can be used. A well-documented and highly used damaged index was given in [25] and later in [26]. The general relation for the damage index of a plastic hinge model was given by Equation (26).
where θ m is the maximum rotation obtained during loading history, θ r is the recoverable rotation at unloading, θ u is the ultimate rotation capacity of the cross-section, β is a strength degrading parameter, M y is the yield moment of cross-section, and Eh is the dissipated hysteresis energy. Having been introduced in the above section, it is important to say that the hysteresis model further utilized in this paper is based on the dependence between plastic rotation and bending moment. Therefore, the damage index evolution for a plastic hinge is given by Equation (27) where Eh PH i (t) is the cumulative dissipated hysteresis energy up to time step t for the ith plastic hinge. Furthermore, in order to quantify the damage record at each of the stories the following expression can be used.
with n being the number of plastic hinges corresponding to jth story while p represents the number of stories. Moreover λ j i (t) is the normalized dissipated hysteresis energy of plastic hinge ith corresponding to jth story with respect to the dissipated hysteresis energy of the jth story (Equation (29)).
where Eh j,storey PH i (t) is the cumulative dissipated energy up to time step t for the ith plastic hinge corresponding to jth story, and ∑ n i=1 Eh j,storey PH i (t end,curr ) is the maximum (or current) value attained during loading history as a sum of all the cumulative dissipated hysteresis energy by the plastic hinges corresponding to jth story.
In a similar manner, the damage index for the whole building can be defined by Equation (30).
where λ j (t) is the normalized cumulative dissipated energy by the jth story with respect to the cumulative dissipated energy of the building overall with its formula given by Equation (31) λ j (t) = Eh j,storey (t)

Case Study
In this section, the implementation of a peak-oriented hysteresis function is presented and further used for modeling the nonlinear behavior of beams in a 3D reinforced concrete moment resisting frame structure. A dynamic nonlinear simulation is conducted, the structural response being obtained using numerical analysis algorithm developed in Matlab/Simulink which is validated in terms of displacements and bending moments using a commercial software.
Although the plastic behavior of three-dimensional structures is very complex, the common seismic design procedures consider that beam end rotation is the only desired source of nonlinearity while shear and axial deformations are considered to be elastic. These assumptions are in close connection with the ability of structural members to dissipate high quantities of input seismic energy through sustained cyclic rotation.
In this paper, in order to account for the nonlinear behavior of structural elements, the concentrated plasticity approach is used [27]. Considering that most of the time, for civil engineering structures, beams are incorporated in slabs, it can be accepted that axial force is negligible and plastic rotations occur only around the horizontal axis of the crosssection. Therefore, bending moment plastic hinges can be adopted for beams. In the case of columns, the plastic hinge mechanism is far more complex. The bending moment values are significant on both principal the direction of the cross-sections and depend on the axial force variation. Consequently, axial force-bending moments (P-m-m) interaction plastic hinges can be utilized. Moreover, the novel algorithm proposed by [10] is used in this paper to implement this type of plastic hinges in FAM dynamic nonlinear procedure.
Moreover, laboratory tests conducted on the response of structural members [28] indicate that strength deteriorates as the number and amplitude of cycles increase, even if the monotonic displacement associated with the maximum strength has not been reached. Consequently, the cyclic model implemented herein, similar to the peak-oriented model in [29] that had been developed firstly by Clough and Johnston and modified by Mahin and Bertero, takes into account both strain hardening and softening behavior. The general bending moment-rotation hysteresis function are presented in Figures 2 and 3. Furthermore, as stated earlier, plastic rotation is considered to be the only source of nonlinearity, and therefore, a bending moment-plastic rotation hysteresis model is a better fit for the implementation in FAM equations. It is important to note that the cyclic model proposed here is only implemented for the beams of the model. Figure 4 together with Table 1 present in detail the five general rules upon which the cyclic behavior model is built. The rules describe elastic, plastic, and softening behavior together with loading, unloading, and reloading principles. The elastic and plastic behavior including loading and unloading are defined only by their corresponding fixed slope values. Softening behavior is also described by a fixed slope, but it occurs only after θ s is reached. Furthermore being a peak-oriented cycle model, the reloading branch is then computed using variable slopes aiming at the overall maximum (m − θ ) point in the I and I I I quadrant in each corresponding case. Furthermore, it must be mentioned that when designing the numerical algorithm the most changeling task in the implementation hysteresis models is the determination, for every plastic hinge, in each time step, of the proper k br i . Being a history variable, it is dependent on the hysteresis path and therefore several key parameters must be memorized throughout the entire analysis in order for the appropriate contextual value to be assigned. Firstly, taking into account that the cyclic model is a peak-oriented type model, the overall maximum and minimum bending moments together with the corresponding plastic rotation values must be stored. At the same time, two tracking parameters must be used in order to identify the current and previous quadrant position of (m i , θ i ) pair. This is necessary because different rules apply for the same position depending on the last reloading point. For example, a point in the first quadrant can be on the reloading branch originating in the first or second quadrant. Moreover, taking into account that the reloading branch has variable slopes, the last plastic rotation associated zero bending moment must also be stored. Table 2 presents the memorized variables and their role in the algorithm.    In order to conduct the nonlinear dynamic analysis, a spatial moment-resisting frame structural model was chosen as illustrated by Figure 5: a 10-story building with 2 spans (L = 6 m), 3 bays (B = 5 m), and the floor height is H = 4 m. The cross-section of the columns is 1.1 m × 1.1 m, and for the beams 0.6 m × 0.3 m, while the Young modulus E = 30,000,000 kN/m 2 and density ρ = 2.5035 tons/m 3 . At the same time, a live load of 11 kN/m is considered acting on all beams. Furthermore, assuming a proportional Rayleigh type damping model, the coefficients α and β are 0.2339 and 0.0107, respectively, and correspond to a 5% damping ratio on the first two eigenmodes. The structural elements are modeled using beam finite elements as presented in [10]. Table 2. Peak-oriented algorithm memorized variables.

Memorized Variable
Role in Algorithm θ -the plastic rotation Used to determine whether the structural behavior is elastic or plastic.
The peak positive pair (θ , m(t)) in Used to define the reloading slopes originated in the the first quadrant first and second quadrant The peak negative pair (θ , m(t)) in Used to define the reloading slopes originated in the the third quadrant third and fourth quadrant The current quadrant location of (θ , m(t)) pair Used to distinguish between The previous quadrant location of (θ , m(t)) pair different reloading situations The last plastic rotation associated zero Used to define the reloading slopes bending moment The structural model is constructed by considering 6 degrees of freedom (DOFs) per joint, 2 plastic hinges for each of the 170 beam elements, and 4 plastic hinges-interacting 2 by 2-for each of the 120 column elements. The total number of plastic hinges gives the number of plastic rotation θ to be computed in each time step, namely 820 elements. Additionally, in each time step the state vector z given in Equation (7) is evaluated. If the full model is considered, there are 720 displacement elements and additional 720 velocity elements included in the state vector. However, by imposing a diaphragm constraint at each of the 10 stories of the structure, the number of displacement components reduces to 390, corresponding to the translations TX, TY, and rotation RZ computed in the mass centers of the slabs, and the remaining rotations RX, RY and translation TZ computed at each joint. The use of a diaphragm constraint is considered to be a good balance between numeric efficiency and accuracy. Additional coordinates condensation can be performed as it is suggested in [3], but the reader is suggested to carefully consider such techniques since the member restoring force matrix K R may become singular once a static condensation is performed. It must be mentioned that the diaphragm constraint implemented here will not produce singular member restoring force matrices K R . As a conclusion, the numerical model used in this case study is described by 1600 unknowns (390 displacements, 390 velocities, and 820 plastic rotations) that require evaluation in each time step. Regarding the structural elements, beams bending capacity was considered 520 kN and the column yielding surface f (x, y, z) ( Figure 6) was defined as previously proposed by G. Powell et al. [30]: where m yp = m zp = 2000 kNm define the maximum bending moments, P p = −11,500 kN is the vertical semi-axis, P b = −8500 kN is the balance point axial force, and (m y , m z , P) are the member forces varying along y, z, and x axis respectively. Furthermore, P b + P p is the cross-section compression capacity and P b − P p is the cross-section tension capacity. Simulation of the dynamic nonlinear structural model was carried out using the North-South and East-West Vrancea 1977 seismic recordings. It must be noted that the defining characteristic of these motions is that they exhibit pulse-like behavior, most of the seismic input energy being concentrated in a short amount of time [31].

Validation of the Dynamic Nonlinear Model
The numerical algorithm elaborated in Matlab/Simulink is presented in (Figure 7), where (a) represents the seismic and scale factor blocks, (b) is the integrator (ode23ttrapezoidal rule and backward differentiation formula with variable-step), (c) multiplies the state variable with the state matrix and feeds back the result, (d) multiplies the plastic displacement vector and feeds the result into (b), while (e) and (f) are the main routines that simulate the nonlinear behavior of the structure. This diagram describes the exact data flow of the numerical routine. In the first instance using blocks (a), (b), (c), and (d), information from previous step is fed into the integrator Equation (7) in order to determine the state vector. Further, the displacements are used in block ( f ) to determine the trial plastic rotation increments (∆θ trial ) and consequently the trial bending moments (m trial k+1 ), as described by Equations (14) and (15). At the same time, the trial K d matrix is populated with slope values describing unloading Equation (13). Additionally, the axial force is established in order to determine columns bending moment capacity as described by [10]. Next the data is fed into the main nonlinear routine embedded in block (e). In this block, using a while loop, the K d matrix and (∆θ trial ) vector are updated to the actual values ( Figure 1, Equations (16) and (17)) in order to determine the real bending moments (m k+1 ) Equations (10) and (11). Furthermore the final plastic rotation (theta_p) and consequently the plastic displacement (x_p) are determined. The interested reader will find in [10] more insights regarding the FAM 3D nonlinear routine in block (e), especially for the return mapping algorithm used. Furthermore, before putting the algorithm to work, the performance of the routine was evaluated using a comparison of the dynamic nonlinear structural response obtained, with the one conducted on an equivalent model developed in a commercial software SAP2000. It must be noted that SAP2000 capabilities of modeling hysteresis behavior is limited to certain default options. Therefore, a similar model to Takeda was used for the FAM 3D validations. Beam (#125 and #651) and column (#1 and #2) bending moments together with the top story displacements time histories response are compared using as seismic input the Vrancea 1977 recordings scaled by a scale factor of 1.4 (Figure 8). The results shown in

Results
In this section the nonlinear dynamic structural model is simulated using the Vrancea 1977 seismic input. Energy response and structural damage quantification are presented and discussed. The energy balance is given in Figure 13, which illustrates a significant percentage of input energy being dissipated through inelastic incursions. These inelastic incursions have the significance of damage, distributed throughout the entire height of the structure. Another observation that ought to be made is related to the energy dissipation mode. The pulse-like nature of the seismic inputs leads to early dissipation of the induced energy through high amplitude deformations as illustrated in Figure 12. At the same time, it can be seen that significant values are also reached by the energy dissipated through viscous damping while low values are encountered by the kinetic and strain energy time variations. The second picture of Figure 13 presents the energy balance error (EBE) time history. Taking into account that each energy type is computed based on the displacement, velocity and acceleration vectors obtained from numerical integration in each time step, the EBE can be seen as a criterion for the evaluation of the global accuracy of the time integration algorithm. The EBE can be calculated using the following relation.
As stated by [32], a tolerance situated between 2% and 5% is reasonable for practical applications.
Furthermore, the concept of damage index is used in this paper to give a numerical insight on the element story and overall damage of the structure.
The Matlab/Simulink model introduced in the above sections has been used to make an incremental dynamic analysis (IDA). A series of scale factors were applied to the time history loading, spanning from 0.6 to 1.4, and a 0.1 incremental factor. Therefore, in Figures 14 and 15 is presented the development of damage indices for the plastic hinges highlighted in Figure 5, across the loading history, with respect to the scale factors used in IDA.  The results shown in these figures highlight that damage index values recorded for columns are higher than those for beams. This fact is also highlighted by the percentage of energy dissipated by the base columns, with respect to dissipated energy for the whole building, in contrast with the rest of the plastic hinges. In Figure 16 the hysteresis energy sum of all base columns plastic hinges together with the rest of the plastic hinges hysteresis energy sum is presented as a percentage from the total dissipated energy of the building (Eh b ). In all the analyses, a minimum of 30% from the total dissipated energy is due to the base columns' contribution.
An obvious second observation is that, as the demand approaches the capacity of the columns, more and more beams are engaged in the energy dissipating process, implying that even though the plastic hinges of base columns are approaching peak-capacity, there are still some beams yet to be involved in this dissipating process. The damage index is shown for the base columns, even though higher than the conventionally understood value for collapse state of 1, should not be perceived as the limit state of the building.
Damage occurring in base columns dominates also the global damage index value of the building, as displayed in Figure 17 and in correlation with Figure 18 where the evolution of 1st story maximum drift at each time step and the story damage index (DI s ) is represented.  The next Figures 19-21, display the maximum drift and damage indices for 3rd, 6th and 9th stories. It can be seen that there exists a nonlinear correlation between maximum drift and damage indices for some of the stories. This fact indicates that for stories with such correlation a higher damage degree may be encountered. Moreover, it can be seen that drift development is mainly linear during IDA but the inherited nonlinearity shown in Figure 22 is due to maximum damage index evolution which is indeed nonlinear. At the same time, in Figure 22 it can be seen that stories with a high root mean squared error value are correlated with higher values of damage.

Conclusions
The main scope of this paper is to introduce a straightforward numerical routine developed in Matlab/Simulink for nonlinear dynamic analysis of reinforced concrete (RC) frames using the insights of the force analogy method coupled with an explicit integration scheme.
The accuracy of the proposed numerical algorithm is further validated through a comparison of the structural response with the results obtained by commercial software commonly used by structural engineers in structural design.
Furthermore, an energetic evaluation of the structural response was conducted. The energy balance time history was calculated using a seismic input scaled with a scale factor of 1.4. In addition, the energy balance error was computed in each time step in order to evaluate the performance of the time integration scheme.
Moreover, in order to showcase the capabilities of the algorithm, an incremental dynamic analysis was conducted on a structural model using Vrancea 1977 seismic inputs scaled with scale factors ranging from 0.6 to 1.4. The structural response data was processed and a damage evaluation procedure, using damage indices, was conducted. The results reveal damage distribution among different plastic hinges and different stories of the model, together with the time evolution of story drifts. The overall building hysteresis dissipated energy and damage index are also displayed.
Compared to the well known classical approaches, the proposed routine has the following crucial advantages: Throughout the entire dynamic analysis, the initial stiffness matrix of the structure is used, and therefore no other numerical algorithm is necessary for the determination of stiffness matrix in each time step; the classical equation of motion described by a second-order differential system is reduced to a first-order differential system which is further transformed into an explicit time integration scheme that uses adaptive time steps; the technique used for implementation of hysteresis functions that describe the degrading behavior of concrete elements during cyclic loading is rather uncomplicated.
To point out some of the important results, the nonlinear dynamic procedure proposed in this paper is the first to present a detailed methodology for the implementation of polygonal hysteresis models in the FAM 3D. Moreover, validation of the results using commercial software shows very good compliance of the structural response as illustrated in Figures 9-12. Furthermore, the performance of the time integration procedure is also evaluated using the energy balance error (Figure 13). The step error reaches a maximum value of 0.7%, which is situated under the limit accepted by civil engineering literature (2-5%). Another interesting result aroused form damage analysis is presented in Figure 16, where the percentage of dissipated hysteresis energy by plastic hinges of base columns and by the rest of the plastic hinges with respect to the total dissipated energy is illustrated. The total dissipated energy by the base columns varies linearly over IDA from 50% to 30%. Therefore, it can be concluded that the overall amount of energy dissipated through plastic deformation is mostly influenced by the cyclic behavior of the base columns. Consequently, as the seismic scale factor is raised the plastic energy amount dissipated by beams increases.
To conclude, the authors' main contribution is the development of a force analogybased algorithm that can be used to conduct, in a simple manner, nonlinear dynamic analysis for RC moment-resisting 3D frame structures.