Comparison between recent implicit time integration methods with frequency dissipation for nonlinear structural applications

The present paper aims to test recent (Truly self-starting two sub-step method and three-parameter single-step implicit method) and classical (Generalized-α, HHT - α, and WBZ - α methods) time integration methods using the geometrically nonlinear Positional Finite Element Method (PFEM). The numerical formulation is based on the total Lagrangian approach and uses the Hessian matrix to obtain the response. The mixed hardening inelastic model applied to PFEM is also presented. Two examples validate the time integration algorithms and the inelastic model. In the first example, the mixed hardening inelastic model is compared with the the bilinear stress-strain model and the elastic-perfectly plastic hinge model, and aspects such as amplitude decay and period elongation are discussed. In the second example, the implemented algorithms are verified in a severe geometrically nonlinear example, considering the influence of numerical dissipation, time interval, and the number of elements in the response. Results show the relevance of numerical damping for numerical stabilization and the good performance of the Generalized-α algorithm.


INTRODUCTION
The response of the equation of motion in structural systems under arbitrary forces can be obtained by modal superposition techniques and methods of direct time integration.The superposition techniques use linear modal analyses and express the response in a series expansion of the vibration modes and modal participation contribution.Its effectiveness is remarkable since the fundamental modes are predominant; however, they are limited to linear cases (Géradin and Rixen, 2015).
In this context, direct integration methods are an alternative, as they allow the consideration of the contribution of high frequencies in the response directly.In a complementary way, they can be applied in linear and nonlinear analyses (Géradin and Rixen, 2015).According to Clough (2003), direct integration methods provide the only complete general approach for the analysis of the nonlinear response and are equally valid in the linear cases, regardless of the structural behavior.These methods are so effective and convenient that time-domain analyses almost always are done by some form of step-by-step analysis, regardless of whether or not the response behavior is linear.
The term "direct" means that no transformation of the equations of motion is performed before the process of numerical integration.According to Bathe (2014), direct numerical integration is based on two ideas.In the first one, the equation of motion need not be satisfied at any time t but at discrete time intervals.The equilibrium equation (including inertial forces and damping) must be satisfied at discrete points in the solution interval.In the second idea, a variation must be assumed for the displacements (or positions), velocities, and accelerations within each time interval.This assumption is determinant in the precision, stability, and computational resource used in the solution procedure.
In the last two decades, researchers proposed improvements in time integration methods or even presented new ones.Chung and Hulbert (1993) presented the family of algorithms called Generalized-α with the numerical dissipation being controlled by the user.The parameter allows the algorithm to dissipate the high-frequencies and minimize the low-frequency dissipation.Besides, the proposed algorithm with the appropriate choice of parameters results in use of the following methods: Newmark, Hilber-Hughes-Taylor-α (HHT-α), and Wood-Bossak-Zienkiewicz-α (WBZ-α).Kontoe et al. (2004) compared the Generalized-α, Newmark, HHT-α, and WBZ-α time integration schemes for a boundary value problem with a deep foundation subjected to earthquakes.The comparison considered precision, CPU time, and the numerical dissipation control of high frequencies.Their preliminary results demonstrate the superiority of the Generalized-α algorithm in earthquake transient analyses.A simple implicit scheme of time integration for systems under large deformations and long duration of time was proposed by Bathe (2007) for cases in which the Newmark method does not conserve the system energy and is unstable.The method is about twice as expensive as the trapezoidal rule per time step, but less steps can frequently be used, and the method remains stable when the trapezoidal rule fails.Kuo et al. (2012) developed a direct integration algorithm for a nonlinear dynamic analysis of structures whose accuracy is of the fourth-order time and self-starting; this is so that errors caused by the initial estimate of acceleration are eliminated.Still, according to the authors, it is possible to use large time steps.Rostami et al. (2013) presented an algorithm where the cubic B-spline method is developed for systems with multiple degrees of freedom.The algorithm uses periodic cubic base functions over time and is conditionally stable.These authors performed analyses of stability and precision, and concluded that the results of the proposed method are coincident with those from the linear acceleration method; however, they are obtained in a faster and more efficient form since the displacement, velocity, and acceleration vectors are calculated independently.
A new family of self-initializing algorithms for dynamic analysis was proposed by Soares Jr. ( 2015) with the implicit formulation being unconditionally stable and the explicit formulation having conditioned stability.Moreover, the method has a second-order precision and a parameter for numerical dissipation control.Noh and Bathe (2019a) sought to improve the ρ ∞ -Bathe method based on the rate of division of the time step and the spectral radius over long time intervals using the least possible number of parameters.The scheme is effective to prescribe in a smooth manner from no amplitude decay to very large amplitude decays, with a correspondingly small period elongation to very large period elongations while maintaining second-order accuracy.Malakiyeh et al. (2019) considered the Bathe method in the wave propagation problem.They also pointed out that this method is unconditionally stable and often used without parameter adjustment.However, one can make use of these parameters although detailed numerical experiments are necessary.For Zhang and Xing (2019), time integration methods with three parameters and with a solution step (like the Generalized-α) lose precision when interpolating the external force vectors; this is despite improving numerical dissipation when modifying the equilibrium equation at discrete time points.To solve this problem and still allow the application of these methods in the solution of secondorder nonlinear differential equations, the authors propose a new method (also with three parameters and a solution step) called the Three-Parameter Single-Step Implicit Method with an additional variable inserted in the update of equations.structural applications William Luiz Fernandes et al.Solids and Structures, 2022, 19(3), e441 3/22

Latin American Journal of
Zhang (2020) tried to carry out a comprehensive study about the Optimal A-Stable Linear Two-Step (OALTS) time integration method.This scheme is implicit with second-order precision for displacement, velocity, and acceleration simultaneously with controllable numerical dissipation.Li and Yu (2020) developed an implicit, self-initializing, unconditionally stable, and second-order precision algorithm for solving nonlinear dynamic problems.The response is obtained through two substeps and uses two parameters to control the numerical dissipation at high and low frequencies.Borrvall and Lawson (2020) pointed out that the Newmark method is prone to numerical instability when the models are submitted to rotational movements despite its unconditional stability and energy conservation.They used algorithms implemented in LS-DYNA®, such as Bathe, HHT-α, and Finite Rotational Dynamics (FRD), to avoid these instabilities and indicate the most suitable algorithms.Behnoudfar et al. (2020) extended the Generalized-α method whose precision is in the range of second-order to third-order but controlled by a single numerical dissipation parameter.The authors state that the method is unconditionally stable and provide a direction for generalizing the method to higher-order schemes.Kim (2020) reviewed and critically analyzed explicit and implicit numerical integration algorithms, including an explicit algorithm implemented by the author himself.Computational aspects, similarities, and orientations on their applications in the transient analysis were detailed from examples.
Researches have been developed using mixed (or combined) hardening models to consider the behavior of materials based on the linear or nonlinear back-stress approach.Seisenbacher et al. (2018) performed a study of parameters for the combined model of isothermal components mechanically loaded (based on Chaboche's model), which are validated by low cycle fatigue (LCF) tests.The authors state that the model allows simulations independent of the load cycles number and the chosen strain amplitude.In the study by Geng et al. (2018), the mixed hardening model is adopted to analyze the residual stresses of butted weld joint of a huge cylinder with an ultra-thick wall via Lemaitre-Chaboche's formulation.Jiang et al. (2021) carried out a similar investigation with Chaboche's formulation but to determine the fatigue life prediction of the 316L stainless steel weld joint.Zhao et al. (2019) investigated the low cycle fatigue (LCF) and LCF-creep behaviors of P92 steel with strain-control mode and fully inverse loading waveforms at a specific temperature using the mixed hardening model.Shi et al. (2020) adopted three hardening models to investigate the size effect on metal sheets on the subsequent yield.They identified that the mixed hardening model presented better results compared to the pure isotropic and pure kinematic hardening.The mixed hardening model was also used by Liu et al. (2021) to verify the influence of lubrication conditions in the metal forming process of tubes by rotary swaging.Yang et al. (2020) introduced a methodology to calibrate Chaboche's mixed hardening formulation parameters.This model combines isotropic hardening with kinematic hardening, which is widely used to describe the change in the yield surface.Overall, mixed models are often considered for describing both the expansion of the yield surface and its translation, which, for example, allows the capturing of the Bauschinger effect.Thus, a more accurate characterization of the simulations is guaranteed.
The present work proposes a comparison of numerical time integration algorithms implemented using Positional Finite Element Method (PFEM) formulation, which is based on the total Lagrangian approach and uses the Hessian matrix (instead of stiffness matrix) to obtain the positional response.The formulation's main advantage in this case is that the chain rule is not explicitly applied to establish the strain gradient (Coda and Paccola, 2008;Coda, 2009;Coda and Paccola, 2010) occurring as a simple numerical matrix inversion.As a result, PFEM can be used to solve any geometrically nonlinear problem (Coda and Paccola, 2010;Coda and Paccola, 2011;Nogueira et al., 2013;Silva and Coda, 2012;Pascon and Coda, 2013).The mixed hardening model is presented in detail and adapted for the proposed formulation to consider the physical nonlinearity.Two highly nonlinear transient examples are used to demonstrate that the PFEM formulation, with the time integration algorithms and the mixed hardening model, can be a significant and viable alternative to geometrical and physical nonlinear analysis of structural plane frames.

THE POSITIONAL FINITE ELEMENT METHOD FORMULATION
The PFEM formulation applied to dynamic problems have been developed by Greco (2004), Coda and Greco (2004), and Greco and Coda (2006).Figure 1 presents the generic parameterized curve representing the element, and its geometry is described as a function of the dimensionless variable ξ (which varies from 0 to 1) as well as the global coordinate system.
The conservation of energy in a mechanical system occurs if the input and output of energy are at balance.If there is some kind of dissipation, the total energy of the system changes along time.The total potential energy of a system can be described as follows: Latin American Journal of Solids and Structures, 2022, 19(3), e441 4/22 where Q(X,t) is the quantity of energy withdrawn from the simple conservative idealized energy Π 0 , Π is the remaining (actual) mechanical energy of the system, X is the coordinate parameter, and t is the time.For a structural problem associated with a fixed reference system (Figure 1), the ideal potential energy function can be written as the composition of the strain energy U t , the potential energy of applied forces P, the kinetic energy K, and dissipation Q as follows: Fig. 1 Generic parameterized curve representing the plane frame finite element.
The strain energy function of the frame U t is considered stored in the initial volume of the body V 0 and is written as an integral of a specific strain energy value u e : where c 1 , c 2 and c 3 are functions of θ 1 , θ 2 , l x (=X 2 -X 1 ) and l y (=Y 2 -Y 1 ); 1/r is the curvature, E is the Young modulus, r is the curvature radius, z is the height of equivalent cross-section, ε med is the engineering strain, ε is the strain measure for geometric nonlinearity, ε p is the plastic strain, and ξ is a dimensionless variable (from 0 to 1).
In the initial position (non-deformed), the strain energy is assumed to be zero.The work of the applied conservative concentrated forces P is written as follows: where F i represents forces (or moments) applied in i direction and X i is the ith coordinate parameter of the point where the load is applied.The kinetic energy is given by the following: (5) structural applications William Luiz Fernandes et al.
Latin American Journal of Solids and Structures, 2022, 19(3), e441  5/22 where i X  is the ith velocity and ρ 0 is the specific mass.The dissipative term is written in its differential form as follows: where q(X,t) is the specific dissipative functional, λ m is a proportionality constant, and X i is the position of any specific point (for PFEM, it is a nodal position).The potential energy of the applied forces may not be zero in the reference configuration (Coda and Greco, 2004).Eq. ( 2) can be rewritten as follows: The energy function can be evaluated by the following approximation: The principle of minimum potential energy consists of the real values that satisfy the equilibrium equations and make the potential energy stationary for all admissible displacements in a certain system.This principle is used in Π0 by differentiating Eq. ( 8) regarding a generic nodal position X s , which results in the following: In Eq. ( 9), s X  is a generic nodal velocity.In a vector form, Eq. ( 9) can be rewritten as follows: where the internal forces vector ∂U t /∂X, the inertia forces F inert. , and the damping forces F damp. are given, respectively, by the following: In Eq. ( 12)-( 13), M, C, X  , and X  are mass matrix, damping matrix, velocity vector, and acceleration vector, respectively.Then, Eq. ( 10) becomes: Latin American Journal of Solids and Structures, 2022, 19(3), e441 6/22 In Eq. ( 14), time t and position vector X are the variables, and F is the external loads vector.It is necessary to discretize the system to perform the numerical integration in the time domain.Rewriting for time t+Δt: which represents the condition of geometric nonlinear dynamic equilibrium equation of motion, whose solution can be obtained using an iterative process.

NONLINEAR POSITIONAL FINITE ELEMENT FORMULATION FOR SOME TIME INTEGRATION METHODS
The following algorithm schemes use the spectral radius ρ ∞ [0; 1] as a numerical dissipation control parameter for the contribution of high frequencies to the response.In the lower limit, there is a complete dissipation of the highfrequency contribution response.On the other hand, there is no dissipation in the upper limit.

The Generalized-α, HHT-α, and WBZ-α methods
These methods have the following characteristics: implicit, second-order time precision, unconditional stability for the linear case, and numerical dissipation control of high frequencies (Chung and Hulbert, 1993).Eq. ( 15) can be rewritten as follows: The position, velocity, acceleration, and external force vectors are given, respectively, by the following: with acceleration and velocity at t+Δt, respectively, equal to: The equation of motion can be described as follows: Latin American Journal of Solids and Structures, 2022, 19(3), e441 7/22 where D M and D C are the vectors related to the mass and damping matrices, respectively, with the known variables at time t.The Hessian matrix is given by the following: with the following parameters: Eq. ( 27) presents the algorithm parameters: α m corresponds to the mass terms, and α f is related to the other terms.In addition, β and γ are complementary parameters of the method.The HHT-α and WBZ-α methods, which can be found in the research of Hilber et al. (1977) and Wood et al. (1980), are obtained considering α m =0 and α f =0, respectively.

Truly Self-Starting Two Sub-Step method
The method proposed by Li and Yu (2020) has significant features.Acceleration is not necessary to initialize the algorithm (truly self-start), but it can be achieved at each time step.Besides, the displacement (or position), velocity, and acceleration are obtained with second-order precision.The method has unconditional stability for linear elastic problems.

First sub-step
The first variation of the total potential energy, Eq. ( 14), can be approximated by the following: and position, velocity and external force vectors can be written, respectively, as follows: Thus, the equation of motion is obtained by the following: ( ) ( ) and the second variation of the total potential energy leads to the Hessian Matrix: ( ) Similarly, the Eqs.( 28)-( 35) for the second sub-step are given by the following: ( ) ( ) ( )

Updating the variables
After the sub-steps, the updated vectors for position and velocity are given, respectively, by the following: and the acceleration equation assumes the form: ( ) The parameters of the method are presented in Eq. ( 47): Latin American Journal of Solids and Structures, 2022, 19(3), e441 9/22 ( ) ( ) where γ and ρ ∞ are the coefficients to control numerical dissipation at low and high frequencies, respectively.The dependent parameters β, χ, α and γ 1 are used to ensure accuracy requirements of the present scheme.

Three-Parameter Single-Step Implicit Method
The method proposed by Zhang and Xing (2019), called the Three-Parameter Single-Step method, with an implicit form strictly satisfies the equilibrium equation in each time step and has a desirable performance according to the authors.The position, velocity, acceleration, and the arbitrary vector θ t+Δt are defined as follows: ( ) Δt 1 The arbitrary vector can be defined as θ t =0 for t=0 and must be updated at each step.The equation of motion is defined to be: The Hessian Matrix is such that: and the parameters are obtained by: The spectral radius ρ ∞ is the only independent parameter, with δ, α, and γ as a dependent parameters to optimize the implicit method.

Research Methodology Flowchart
The following flowchart algorithm (Figure 2) shows the research methodology of the present work.The standard Newton-Raphson is used as an iterative process and the Euclidean norm as a convergence criterion.

MIXED HARDENING MODEL APPLIED TO PFEM
A mixed hardening model combines isotropic hardening (expansion of the elastic surface) and kinematic hardening (translation of the elastic domain center).It is simple and has proximity to the experimental occurrence, particularly for cases of plasticity in metals.Greco (2004) presented the implementation of this inelastic model in the PFEM formulation.Figure 3 shows a graphical representation for a one-dimensional element.
The yield function f and yield criterion used are presented, respectively, by the following: ( ( ( ) In Eq. ( 57), K is the plastic module.The measures α and q (back-stress) are internal variables that measure the expansion and change of the center of the elastic domain, respectively.The yield stress of the material is σ Y .A state is said to be admissible when Eq. ( 58) is verified.
By assumption, considering a strain decomposition, the stress is given by the following: The evolution of the plastic strain (or plastic flow), which is described by Eq. ( 60), and of the internal hardening variable, Eq. ( 61), are given by the following: where γ is called plastic multiplier.The evolution of the internal variable q is given by Ziegler's Law: and H is the kinematic hardening module (constant).Two conditions govern the admissibility of stress states and the evolution of plastic deformation.The first one, called the Kuhn-Tucker condition, is given by the relationships in Eq. ( 63): which guarantees an admissible state.In this case, a plastic deformation occurs only when ( ) f ,q, 0 σ α = , i.e., 0 γ > .The second one, called consistency condition, is described as follows: and requires for p ε that 0 γ > .It must persist on the yield surface so that ( ) Based on the two conditions, Eq. ( 63) and ( 64), and the yield criterion function, Eq. ( 57), it is possible to write: When K=H=0, the model becomes perfectly plastic.For K>0 and H=0, pure isotropic hardening is obtained.

Solution based on incremental terms of the variables
The problem solution involves the consistency condition and the Kuhn-Tucker condition.The insertion of a trial stress state (supposedly elastic) allows writing the solution in incremental terms: ( ( with the initial conditions p n ε , n α , and n q to solve the differential equation problem.The following relations are based on an incremental trial state and define whether plastic correction is necessary as is shown below in Eq. ( 72): trial n 1 0, 0 then elastic step f 0, 0 then plastic step Comparison between recent implicit time integration methods with frequency dissipation for nonlinear structural applications William Luiz Fernandes et al.Solids and Structures, 2022, 19(3), e441 12/22

Latin American Journal of
Chart 1, which is based on Eqs. ( 66)-( 72), shows the step-by-step solution for the incremental plastic correction.For details regarding the algorithm, see Greco (2004) and Simo and Hughes (2006).This incremental plastic correction is updated at each time step.

Chart 1
Step-by-step algorithm for mixed hardening model.

Given the initial variables
2. Given the elastic trial stress (Eq.58) in a point, calculate by PFEM and Gauss Quadrature:

EXAMPLES
The purpose of the two numerical examples is to show the limitations of recent time integration schemes to solve dynamic nonlinear problems.Thus, the first numerical example is essentially physically nonlinear, and the second one is geometrically nonlinear only.In both examples, the spectral radii were considered as ρ ∞ =0 (complete numerical dissipation) and ρ ∞ =0.9 (to avoid convergence problems for most algorithms).Specifically, γ=0.1 and γ=0.4 were applied to the Li and Yu (2020) method since, for γ=0.0, the D M and D C vectors would be undefined, because γ appears in the divisor of these expressions, Eqs (32)-( 35), in the first sub-step.Similarly, for γ=0.5, α assumes the value 0.0, and the vectors D M and D C become undefined again, because α is in the divisor of these expressions, Eqs (40)-( 42), in the second sub-step.Both situations are independent of the ρ ∞ , even for cases without physical damping (C=0).In all cases, the γ and β parameters for the Newmark method are, respectively, 0.25 and 0.5.

Toridis-Khozeimeh plane frame
The first example consists of a portal steel plane frame with fixed ends (Figure 4) submitted to a lateral step load (P 0 =444.82kN) at the top left corner.Toridis and Khozeimeh (1971), Marur and Kant (1994), Toi and Isobe (1996), Chan and Chui (2000), and Silva et al. (2011), among others, studied this case.Beam and columns are W16x36 and W12x120 steel I profiles, respectively, discretized in three finite elements per bar.Young modulus and yield stress are equal to 207GPa and 470MPa, respectively.The mass of steel profile members is lumped at element nodes and are multiplied by 625.In this work, tolerance for the iterative process and time step is considered 10 -5 and Δt=0.001s, respectively.(Toi and Isobe, 1996).
Results presented in Figure 5, position variation X(t) on applied force degree of freedom, indicate the need of use of at least three finite elements per bar.Amplitude range of position variation indicates inelastic effects.Both classical and recent direct time integration schemes presented practically same results.No significant effect of numerical dissipation (based on parameter ρ ∞ ) was observed in these responses.The amplitude of vibration with the bilinear stressstrain model with an inelastic branch slope equal to 1/10 of the elastic branch (Toi and Isobe, 1996, 4 elem.) is less than the mixed hardening model (present work, 3 elem.)combined with other time integration methods (Figure 5).It should be emphasized that Toi and Isobe (1996) employed the Adaptively Shifted Integration (ASI) technique to the cubic element based on the Bernoulli-Euler hypothesis.Thus, there is the possibility of the formation of a plastic hinge at Gauss point or even at the ends of the element, which occurred in this example.Therefore, the formulation used by Toi and Isobe (1996) allows more energy dissipation than the one used in the present work.This fact justifies the difference in the results for the stationary vibration period.The other responses coincide with Newmark (Figures 5a, c, and d).For ρ ∞ =0.9, the Li and Yu (2020) scheme has a period elongation about Newmark (Figure 5b).The WBZ-α and Zhang and Xing (2019) presented convergence problems at t~0.8s (Figure 5a) and t~0.5s (Figure 5b), respectively, for the same ρ ∞ .
Similar behavior is observed in Figure 6.Results show a close agreement with Chan and Chui (2000) and Silva et al. (2011); see Figures 6a, c, and d.However, the contribution of the high modes in response (ρ ∞ =0.9) increased the period for the Li and Yu (2020) algorithm (Figure 6b).It is worth mentioning that inelastic effects predominate over geometrical nonlinearity in this example.Therefore, introducing the contribution of high frequencies to the response (ρ ∞ =0.9) did not significantly change the amplitude of positions or even the number of iterations per time step, compared to ρ ∞ =0.0.perfectly plastic hinge model (Chan and Chui, 2000;Silva et al., 2011).
Two aspects must be highlighted in Figures 7 and 8.The amplitude decay of the Generalized-α and the HHT-α (Figures 7e and f) is caused by the contribution of high frequencies in the response.In addition, there is a period elongation for Li and Yu (2020) response, with a half wavelength lag in t~17s (Figure 8a).The other time integration methods do not present these behaviors.Figures 8b, c, and d indicate an oscillatory behavior of the response obtained from the Li and Yu (2020) and Newmark algorithms, apparently indicating some instability for the physical nonlinearity.However, in the Generalized-α and HHT-α schemes, the amplitude decay of the response over time results in an elliptical spiral (Figures 8c and d) similar to those found in systems with physical damping.

Three-spring plane frame
The structural system consists of a steel plane frame with four steel bars (Figure 9).The bars have a circular hollow section with 100mm external diameter and 5mm thickness.The system is pinned at the ends and has three rotational springs.Young modulus, specific mass, and spring stiffness are E=200GPa, ρ=7850kg/m 3 , and k=16881.15177N/m,respectively.One finite element per bar was used, and the magnitude of the periodic triangular load is P 0 =30kN.The convergence tolerance for the iterative process is 10 -7 .It is intended to evaluate the influence of Δt in the time evolution of the algorithms without viscous damping.Results presented in Figure 10 are related to the vertical position variation on applied force degree of freedom.Initially, only one finite element per bar was used in this numerical example.Amplitude range of position variation indicates strong nonlinear geometrical effects.Figure 10b indicates the need of use of numerical dissipation (based on parameter ρ ∞ ).It is worth to remember that the smaller is the value of parameter ρ ∞ , the higher is the numerical dissipation.It can be seen in Figure 10a that the Generalized-α, Li and Yu (2020), and WBZ-α methods provided a complete time-history response for ρ ∞ =0, considering Δt=5x10 -3 s.The amplitude of vibration for WBZ-α is the smallest in the reloading interval.For the Li and Yu (2020) method, different range values provided different amplitudes in the reloading stage.For ρ ∞ =0.9, no method provided a complete time-history (Figure 10b).On the other hand, the Generalized-α, Li and Yu (2020), γ=0.1, and WBZ-α methods exhibited time-history up to 10s (ρ ∞ =0.0; Δt=1x10 -3 s).The response curves for Generalized-α and WBZ-α are superimposed over almost the entire analysis time (Figure 11a).For ρ ∞ =0.9,only Generalized-α and HHT-α provided a complete time-history response (Figure 11b).The numerical damping in the reloading stage concerning the previous cases is remarkable along with that in the free vibration interval (9s<t<10s).Figure 12a shows that only Generalized-α and WBZ-α provided the complete time-history response and are superimposed before the snap-through (t<2.2s) and at the start of reloading (t>5s).However, for ρ ∞ =0.9,only Zhang and Xing (2019) and WBZ-α did not show complete time-history (Figure 12b).It is noteworthy that in free vibration (9s<t<10s), the Generalized-α does not present large amplitudes like the other methods (Figure 12b).For the shortest time interval of this example (Δt=1x10 -4 s), two algorithms provided the complete time history for ρ ∞ =0: Generalized-α and WBZ-α (Figure 13a).For ρ ∞ =0.9, the Generalized-α and HHT-α responses showed a lag between them in reloading, despite the amplitude of the position in both responses being equivalent (Figure 13b).Comparing Figures 10-13, higher modes contribute significantly to the response by decreasing Δt (ρ ∞ =0.0).However, the amplitude increases, diverging from the response for the refined mesh model (5<t<9s).Now, considering a mesh refinement (5 and 10 elements per bar) for the Newmark algorithm, the time-history responses are presented in Figure 14.The spectral radius reduces the response amplitude (5<t<9s) by bringing it closer to the refined mesh model despite the need for a greater number of iterations per time step.Table 1 presents the time march progress for the implemented algorithms.Newmark, widely used in practice, did not present a complete time history in any of the evaluated Δt (considering one finite element for each bar), reaching a maximum time march progress of 57.06%.The least efficient method in this case was that of Zhang and Xing (2019),with no complete time history in any of the proposed Δt and with a maximum time march progress of 5.42%.The most efficient method was the Generalized-α, presenting a complete time history in all cases evaluated, except for Δt equal to 5x10 -3 s, with the contribution of high frequencies in the response.As the system is geometrically nonlinear, subjected to forced vibration, and numerically damped, the phase space presents characteristics of systems with chaotic behavior (Figure 15).For ρ ∞ =0.0, the range of velocities and positions (concerning the undeformed configuration) indicate accumulation of kinetic and strain energy for the Generalized-α and HHT-α algorithms, without amplitude decay effect.The same occurs with Newmark.There is a significant reduction in velocities and positions, indicating energy loss caused by numerical damping (ρ ∞ =0.9).The phase space diagrams demonstrate that the responses are equal for Generalized-α and WBZ-α (ρ ∞ =0.0) and for Generalized-α and HHT-α (ρ ∞ =0.9).Therefore, the Generalized-α algorithm seems to incorporate the best numerical dissipation characteristics of the HHT-α and WBZ-α methods for the evaluated conditions.

Fig. 3
Fig. 3 Evolution of the yield surface and yield stress.

Fig. 6
Fig. 6 Position versus time (top left corner) for steel plane frame with the mixed hardening model (present work) and elastic-perfectly plastic hinge model(Chan and Chui, 2000;Silva et al., 2011).

Fig. 7 22 Fig. 8
Fig. 7 Amplitude decay (top left corner) for steel plane frame with the mixed hardening model.
Comparison between recent implicit time integration methods with frequency dissipation for nonlinear structural applications William Luiz Fernandes et al.

Table 1
Time march progress for time integration algorithms (t f =10s) -three spring plane frame.