Numerical error assessment in nonlinear dynamic analysis of structures

ABSTRACT The complex nature of exciting forces in earthquake engineering applications mandates employing numerical methods to obtain structural response. Several numerical methods are used in earthquake engineering to solve the ordinary differential equation of motion. Previous studies assumed the appropriate time step as a constant fraction of natural period of the system disregarding period range. This study is concerned with assessing, by numerical experiment; the numerical error of commonly used schemes in nonlinear dynamic analyses and assessing their appropriateness for different excitations with different structural periods. The current investigation involved testing Newmark, HHT and HHT1 methods. The three methods were tested for nonlinear single degree of freedom system representing a typical concrete bridge structure. Three different exciting forces were used to test the schemes; half-sine pulse, harmonic force, and actual ground motion record. Three natural periods were used to conduct the experiment, representing short, medium and long period systems including the resonant condition. Three structural damping ratios, representing lightly, moderately and heavily damped systems, were used to assess the damping effect on the accuracy of the schemes. The results of this investigation indicated that commonly used assumption of time step as a constant fraction of natural period of the system disregarding the period range could result in significant numerical errors. The study also showed the significant effect of damping ratio of the system on the accuracy of the numerical schemes. The study presents a recommendation matrix of the most appropriate time step for each scheme for different applications and structural conditions.


Introduction
The complex nature of dynamic exciting forces in earthquake engineering applications mandates the use of numerical methods to obtain the response of structures to ground motions. Many numerical methods are in use in structural dynamic analysis to solve the ordinary differential equation of motion. This study is concerned with assessing the effect of time step selection on the numerical error of some commonly used numerical schemes in nonlinear CONTACT Wael M. Hassan wmhassan@alaska.edu University of Alaska, Anchorage, USA dynamic analysis of structures and evaluating the suitability of each method to different excitations and structural characteristics. Several numerical methods for solving ordinary differential equations are already in use by structural analysis software to tackle the dynamics of structures associated with earthquake excitations. Some of these methods are well-researched for accuracy and stability for linear systems while the accuracy of some other methods still needs more investigation. The research on assessment of time step selection in these methods in analyzing nonlinear structures is limited.

Scope and objectives
For both linear and nonlinear systems, the available research studies suggest selecting the analysis time step as a constant fraction of natural period of the system regardless of the period range. The validity of this assumption for all period ranges and for nonlinear systems is questionable. This study is concerned with assessing the validity of that assumption for nonlinear systems. The current study is aiming to evaluate the numerical error and stability of several numerical methods when applied to solve the ordinary differential equation of motion for nonlinear SDOF systems. The following objectives are of particular interest: exploring the available literature on the numerical solutions of the ordinary differential equation of motion of SDOF systems with emphasis on numerical methods used for earthquake engineering applications; conducting an error analysis and stability study on select numerical methods to assess the effect of mathematical and structural parameters of the accuracy of each numerical method; identifying and recommending the most accurate numerical method for each earthquake engineering application and finally recommending the most accurate and practical time step for each application and structural condition. The current study will explore the available literature body on the numerical methods commonly used in earthquake engineering applications with an emphasis on those methods implemented in commercial and research software packages. The proposed analysis matrix and analytical procedure to assess the accuracy of different numerical schemes will be presented along with methods of error calculation. The results of the error analysis performed as well as the discussion of analysis results will be presented as well as suggestions for time step and numerical method usage for each application.

Numerical integration methods for earthquake engineering applications
The numerical schemes used in solving the dynamic equilibrium equation of motion, or any other ODE, can be classified into explicit or implicit methods.
In the explicit methods, the differential equation at time t is used to predict the solution at time t + h where h is the time step. The explicit methods are conditionally stable with respect to the size of time step. On the contrary, implicit methods satisfy the differential equation at time t after the solution at time t-h is obtained. Implicit methods can be conditionally or unconditionally stable. A significant number of accurate, higher-order, multi-step methods are available for numerical solution of ordinary differential equations for several applications. The solution in these multi-step methods is assumed to be a smooth function with continuous higher derivatives. This assumption could work well with linear structures. However, the nonlinear structures usually produce a non-smooth discontinuous higher derivative, namely the acceleration of the second derivative of displacement. The underlying reasons for that discontinuity are the nonlinear nature of the hysteresis response and the buckling and contact of structural elements, [1]. This makes the multi-step methods not very good candidates for nonlinear analysis. Wilson [1] states, after significant experience in this field, that only one-step; implicit, unconditionally stable methods can be used for numerical integration purposes in seismic analysis of practical structures. Derived from some of the methods mentioned in the previous section, earthquake engineers commonly use particular time-stepping methods such as Central Difference Method (CDM), HHT-α method; including HHT and HHT-1, and Newmark-β method; including linear acceleration method (LAM) and average acceleration method (AAM). Wilson-θ method is another less popular method used in structural dynamics applications. Forty years after the development of Wilson-θ method, Wilson [1] stated that he can no longer recommend his method for the direct integration of the dynamic equilibrium equations. The reason is that a peculiar property of this method was discovered that is a significant overshoot of the exact solution in the first few time steps (Hughes [2]). A good candidate numerical method for structural dynamics is a controversial topic. However, it is generally agreed that the method should possess the following properties to be competitive, [2]: be unconditionally stable when applied to linear systems, no more than one set of implicit equations to be solved at each time step, possess secondorder accuracy and controllable algorithmic dissipation (numerical damping) in the higher modes and is self-starting. The commonly used CDM, Newmark-β LAM and AAMs schemes are presented elsewhere [3].
The HHT-α method was presented for the first time by the finite element scientists Hilber, Hughes and Taylor in 1977, [4]. The HHT-α method is a generalized form of the Newmark-β method. This method reduces to Newmark-β method for α = 0. The HHT method is a dissipative one-step implicit method that allows for numerical energy dissipation for the higher modes maintaining second order accuracy. Depending on the choice of input parameters, the method can be unconditionally stable. The HHT-α method can be second-order accurate and unconditionally stable, for linear systems, if the following method parameters are satisfied: However, it is possible to use α higher than 1/3 but no more than 1/2. HHT method is a numerical scheme that attempts to increase the numerical damping without degrading the order of the accuracy of the method. The smaller the α, the higher the numerical damping. The HHT-α method is particularly useful in earthquake engineering simulations incorporating many degrees of freedom MDOF. It also helps when it is desirable to numerically attenuate or dampen the response at high frequencies. If β and γ are selected according to Equation 3 and Equation 4, increasing α decreases (dampens) the displacement response at frequencies above 1/ (2h). It is worth mentioning that structural damping terms in HHT-α equation accommodate non-classical damping case.

Existing recommendations for selection of time step h
Newmark originally suggested a constant time step from h = T/6 to h = T/5 to achieve a minimum error in the numerical solution; where T is the fundamental natural period. Other researchers, such as Humar [5] and Chopra [3], suggested h = T/10 for acceptable numerical solution. However, the latter value seems very un-conservative and hence will be investigated. Bathe and Wilson [6] suggested that a satisfactory numerical result can be only obtained using the very small time step of h = T/100, which is unaffordable and computationally demanding for most applications. Because of this controversy, a medium value of h/T = 0.05 is typically used as a default value in finite element software packages, [7]. Pavic et al. [7] indicated that the use of the standard h/T = 0.05 time step to integrate numerically a damped SDOF system under a resonant harmonic excitation may lead to a considerable underestimation of the calculated steady-state response. This underestimation is more significant if the damping ratio is 2% or less. They recommended that the numerical error could be tolerated up to damping ratio of 1.3%. A time step refinement is necessary for damping less than 1.3%. They selected time step refinement of T/30 to T/50 based on damping ratio. The numerical error increased with decreasing damping ratio for that particular case of resonant condition. The underlying reason for this numerical error was the period elongation. Systems with damping higher than 5% are less prone to such great error. The main sources of numerical errors are the period distortion (elongation or shortening) and the numerical damping. Many researchers investigated thoroughly these effects for linear systems. Chang [8] indicated that the period distortion significantly differs between different numerical schemes. He also confirmed that the period distortion is reduced with decreasing time step. Consequently, the error in amplitudes of sine and cosine terms increases with increasing time step. He tested the numerical error for linear system with harmonic excitation. Chang concluded that the steady-state response for linear systems can be exactly obtained, with a very small time step, for any member of Newmark method if a linearly varying load is assumed within each time step, such as an earthquake ground acceleration record. However, that does not apply to transient response. He also found that period and amplitude distortions are cumulative for the Newmark method except for AAM. He also concluded that numerical damping is effective only for the transient response and not for steadystate response of harmonic excitation.

Computational efficiency of numerical schemes
Besides defining the appropriate time step size for accurate results, it is important also not to achieve this accuracy on the expense of the efficiency of the method. In other words, the time step could be the same for equal accuracy for two algorithms but the computing efficiency reflected by computing time is smaller in one algorithm than the other. In this case, the latter algorithm will be the best candidate because of it is accuracy and efficiency. In an attempt to address this topic, Negrut et al. [9] tested several low order numerical integration formulas for rigid and flexible multi-body dynamics as an application for mechanical engineering. The authors identified the importance of low order methods since the majority of large real-life models contain discontinuities, friction, and contacts that make low-order integration schemes the 'only viable robust candidates capable of handling these classes of problems'. They studied the efficiency of the numerical scheme represented by order of convergence and the time of processing needed by the computer central processing unit CPU for different methods. The results of the study were in favor of Newmark family methods. The HHTI3 method also proved to be very efficient.

Model structure
The current investigation is concerned with estimating numerical error for numerical schemes commonly used to solve the equation of motion of SDOF oscillator. The SDOF oscillator is the primary system that represents dynamic response of structures, even in MDOF, since the response spectrum analysis is based on equivalent SDOF system. The model structure selected in this study is a typical California Department of Transportation highway overpass concrete bridge [10]. The bridge is five span overpasses with 150 ft. typical interior span. The bridge pier height is varied in this study to vary the structural period as an influence parameter on the accuracy of numerical scheme. The pier cross section is 60 × 60 inches with 2% Gr. 60 steel. Seismic weight is 2000 kips.

Analysis matrix
The proposed analysis program comprises performing nonlinear numerical analysis simulation of the SDOF oscillator described in Section 4.1 using the criteria depicted in Figure 1. The analysis is performed using OpenSees [11,12] and Matlab software platforms. The analysis criteria comprise testing numerical error and stability under the effect of different loading parameters, structural parameters, and numerical analysis parameters. The loading parameters include testing the effect of different loading excitations, namely half-sine pulse excitation, harmonic excitation and ground motion earthquake excitation. The structural parameters include changing bridge period to represent short, medium and long period ranges for different structural damping ratios. Details of the parameters are presented next.

Experiment parameters
Loading function Three different loading functions are used in the current study. The first is a half-sine pulse wave; the second is the harmonic excitation while the third is an actual ground motion acceleration record. The pulse excitation is defined as a sudden applied force with short excitation duration relative to structural period. When the excitation duration is very short relative to the structural period, the force can be considered impulsive. Pulse and impulse forces have numerous applications in the field of modeling applied natural and man-made forces on structures. Near-fault earthquakes, external explosion shock wave, tsunami and flood shock wave, lateral shock of vehicles and trains, and machine impact forces are examples of such forces. The applied pulse excitation in the current investigation is a half-sine pulse that can be expressed as: where p o is the amplitude of the sine pulse and 2π is the forcing frequency ω. The pulse duration t d is 0.50 second. The half-sine pulse excitation is particularly chosen among other shapes of pulse functions, such as rectangular and triangular pulses, since it best matches the shape of the most energetic component of a near-fault ground motion.

Harmonic excitation
The harmonic excitation has many engineering applications. This includes, for example, forced vibration tests to determine modal parameters for structures, the effect of machine vibrations in mechanical design, cyclic loading tests and fatigue related cases of loading. This excitation has a period of 1 second and can be expressed as: The total duration of the harmonic excitation is 6 seconds. After that, the equation of motion's right hand side will be zero implying a free vibration response of the structure with decaying amplitudes of motion until reaching zero displacement.

Ground motion
The ground motion adopted in this study is a classical ground motion record that is used extensively in research studies (El Centro 1940 record) obtained during the Imperial Valley, California earthquake of 18 May 1940. This is one of the earliest full records obtained from a significant earthquake in the U.S.
(El Centro 1940 north-south component record of Imperial Valley Irrigation District) [13]. The peak ground acceleration of this record is 0.319 g. In structural engineering design and assessment of buildings and bridges, usually many ground motion records are used to represent the expected variation of earthquakes at a certain site. These records are statistically averaged for the purpose of design. Seven to twenty records are normally used. However, in the current study, the purpose is to estimate the comparative numerical error rather than delivering a sound structural design. Thus, only one record is used to obtain the relative error between the different numerical schemes.

Structural parameters
Natural period of the structure The dynamic response of structures to applied excitation varies according to the natural period of the structure. For example, short period structures are more affected by a typical ground motion for a stiff soil condition than longer period structures. Another example is the response of the structure to harmonic motion. The period of that motion is 1 second. If the period of the structure is close to the period of excitation, the amplified displacement resonance response takes place. In this case, the effect of damping is significant. If the two periods are far apart, the deformation of the structure is much less intense and the effect of damping is less significant. It is intended in the current study to investigate the numerical error associated with resonant condition as well as other non-resonant responses of the structure. Accordingly, three period values are selected as following; T ¼ 0:50 sec to represent short period structures, with an equivalent pier height of 242.5 in.,T ¼ 1 sec to represent medium period structures, with an equivalent pier height of 384 in. and T ¼ 2 sec to represent long period structures, equivalent pier height of 611.4 in. The selection of the 0.50 sec. period represents a common period of short stiff structures and bridges, buildings up to 5 stories, and buildings designed to respond elastically to earthquakes. The 1 second period is common in flexible bridges and higher than 5 story structures. It represents buildings that can dissipate the earthquake energy through inelastic deformation. The long period range of 2 sec. represents tall structures, base isolated structures, and very slender bridge piers. It also implies structures founded on soft soils.

Structural damping
The structural damping is the phenomenon of the decay of displacement response amplitudes with time. Damping represents internal energy dissipation of the structure. Damping results from many sources such as internal heat generation, friction with air, soil radiation damping, and friction at crack and joint interfaces of the structural elements. Most structures possess a damping ratio ζ less than 20%. The most common damping ratio for concrete structures within normal operating conditions is 2%-5%. This ratio becomes about 7%-10% if the concrete structure is close to or at the nonlinear range [15]. In an 'uncracked' elastic concrete structure, the damping ratio is about 2%. The maximum practical range of damping is 15%-20% which can exist in a wood structure close to its inelastic limit. Damping is beneficial during earthquakes to quickly dissipate the earthquake energy and dampen out the motion. Therefore, damping is artificially enhanced in the structural design process by adding supplemental damping devices. In the case of using supplemental damping devices like viscous dampers, the damping ratio could easily top 10-15%. Based on the above-mentioned facts, it is decided to select three damping ratios in the current analysis to reflect the commonly used values in practice as following; ζ ¼ 2% to represent elastic stiffer concrete structures, ζ ¼ 5% to represent more flexible general case of concrete structures and ζ ¼ 20% to represent upper bound of practical damping range or structures with supplemental damping devices.

Structural analysis type
The structural analysis performed in the current investigation is nonlinear analysis. A strength-deformation relationship has to be defined in the inelastic structure for the purpose of nonlinear analysis. The concrete compressive strength f c = 4 ksi and modulus of elasticity E c equals 57000 ffiffiffi ffi f c p . The yield strength of steel reinforcement is f y = 60 ksi. An elastic-strain hardening concentrated plasticity element is used for nonlinear modeling. The post yield tangent is 0.01 the elastic slope. The moment curvature analysis was performed to obtain moment-curvature diagram that can be related to force-displacement diagram through column height functions. The yield strength is 130,000 kip.in and the yield curvature is 0.65 × 10 −4 1/in. The OpenSees library Concrete02 and Steel01 were used.

Numerical analysis parameters
As mentioned earlier, two of the most commonly used numerical methods in earthquake engineering applications are to be evaluated. These methods are Newmark-β and HHT-α method.

Numerical scheme parameters
In Newmark-β method and HHT-α method, the choice of method parameters γ, β, and α may have an effect on the numerical scheme accuracy and stability, especially for nonlinear systems. Therefore, a preliminary study was performed varying these parameters confirmed that the most common values used in the current dynamic analyses practice are the best choice. This selection implies the investigation of the following schemes: Constant average Acceleration Method, AAM; where γ = 1/2, β = 1/4, and α = 1 HHT Method; where γ = 1/2, β = 1/4, and α = 1/2 HHT1 Method [16]; where γ = 1/2, β = 1/4, and α = 1/2. The main difference between HHT and HHT1 is in the trial step. HHT only updates the velocity and acceleration to be the predicted ones at t+ αh, however, HHT1 in addition to updating the derivatives updates the displacement too.

Time step
The effect of time step selection on the accuracy of numerical methods for nonlinear dynamic analysis is the main deliverable of the current study. This effect is directly informative to both engineers and mathematicians to suggest suitable time step relative to structure's period. The selection of time step is usually tied in the literature as a constant fraction of the natural period of the structure. For example, it is recommended to use h = 0.10T for reasonable accuracy for Newmark-β and central difference methods for linear analysis. However, this value is restricted to linear system and may not hold accurate in nonlinear analysis with different types of excitations, different damping ratios and different period ranges of the structure. Accordingly, the analysis is performed in this study using time step of 0.005, 0.01, 0.02, 0.025, 0.04, 0.05, 0.10 and 0.20 seconds. These values can be related to the structural periods used in Section 4.2.1 to give h/T choice recommendations. Time steps up to 0.025 second are considered, in the earthquake engineering community, appropriate for ground motion record analyses. Larger time steps, if accurate enough could be used for other types of excitations.

Exact solution
The exact solution for displacement response to a half-sine pulse or harmonic excitation for a linear system can be easily derived from first principles of dynamics using classical or Duhamel's integral methods, which can be found in Clough and Penzein [15]. However, the 'exact' solution for nonlinear system is simulated by solving the problem numerically with a very small time step (0.001 sec.) compared to the practical time-0step range for routine analyses. The exact solution of the equation of motion for a nonlinear SDOF subjected to ground acceleration record is practically impossible due to the jagged random nature of the ground acceleration along with material nonlinearity of the system. Therefore, the 'exact' solution adopted is the current study is the numerical solution with an extremely small time step that can be considered impractical for routine analyses. The time step used for 'exact' analysis is 0.001 second, which is ten times smaller than the least time step used in professional earthquake engineering practice for ground motion analysis.

Numerical error calculation
The numerical error is estimated as the absolute value of the difference between peak positive (and negative) displacement response obtained using exact solution and that obtained using numerical analysis. Because of its practical importance, the peak displacement error is particularity chosen since it is the main input used for design purposes. Details of error computation is found in a companion study [17]. The peak displacement positive error E +ve and negative error E −ve can be expressed as: where u ext. and u num. are the positive or negative peak displacement responses of the exact and numerical solutions, respectively. The error at time instances other than peak displacement is graphically presented if it has significant implication such as describing numerical damping or period elongation errors of numerical methods. From a structural viewpoint, the error at time instances other than peak displacement could be important for design purposes for calculation of post-peak response rotations to determine inelastic rotation demand and decide the rotation capacity-demand relationship for some structural elements. The error over the full displacement response history is generally described in a subsequent section. Other applications that require knowing the error over the full history is estimating vibration characteristics to check structures serviceability.

Newmark-β method (AAM)
Half-sine pulse excitation Figure 2 presents displacement response history of nonlinear SDOF system under half-sine pulse excitation obtained using AAM for 2% and 20% damping ratios. The analysis is performed using OpenSees simulation using concentrated plasticity model. The numerical solution tends in general to underestimate the peak inelastic displacement and the free vibration response after the end of pulse forcing phase. The error in peak inelastic displacement and free vibration response are similar. The error increases with increasing natural period of vibration. Period elongation of the numerical scheme is evident. This period elongation decreases with increasing natural period of the structure. This may suggest constant systematic numerical scheme error that becomes more apparent as the natural period gets shorter. It is evident too that the period elongation is proportional to the time step selected. Numerical damping is not evident in the current numerical scheme. It is note mentioning that the damped period of free vibration does not coincide with the natural period of the system. Instead, it vibrates at much longer period. The numerical solution has nothing to do with this period elongation. This elongation is structural and is introduced because of the yielding of the system. Yielding   structures inherently have additional damping of the response due to yielding, which is reflected by additional period elongation. This structural period elongation is more significant in longer period systems. It is also noteworthy that the system vibrates in its free vibration phase about a shifted position from the original equilibrium position due to the permanent inelastic displacement of a yielding system. Figure 3 shows the effect of time step on the accuracy of AAM for nonlinear SDOF system under half-sine pulse excitation. It can be clearly observed that the error is proportional to the natural period of the system for the same h/T. This confirms that is it is inappropriate to assume the time step as a constant fraction of natural period disregarding period range. As opposed to the inverse proportionality generally observed in linear systems [17], [18], the numerical error is proportional to damping ratio in nonlinear systems. For heavy damping (20%), the time step should be limited to 0.01T for short period systems and 0.005T for medium and long period systems for acceptable numerical error (<1%). In medium damping range, the time step can be relaxed to 0.01T for all period values. For lightly damped systems, the time step is limited to 0.04T, 0.025T, 0.02T for acceptable  numerical solution for short, medium, and long period inelastic systems, respectively. Figure 4 displays the full displacement response history of nonlinear SDOF system subjected to harmonic excitation for 2% and 20% damping. The numerical solution for harmonic force tends in general to underestimate both the peak inelastic displacement and the free vibration response after the end of forcing phase for systems with light damping. The amount of error in both phases is not close like the case of pulse excitation. In the free vibration phase, the error is proportional to the natural period of the system for the same h/T  value. However, the relationship is not as clear for the peak inelastic displacement in the forced vibration phase. For heavily damped systems (20%), the numerical solution underestimates the peak inelastic response in forced vibration phase; however, the error in the free vibration response is not as systematic. Slight period elongation of the numerical scheme is evident in the free vibration phase. Period elongation decreases with increasing natural period of the structure and damping ratio. It is also evident that period elongation error decreases with decreasing time step. Figure 12 shows the error trends in estimating peak displacement response of nonlinear SDOF system subjected to harmonic excitation using Newmark-AAM. As opposed to the inverse proportionality generally observed in linear systems [17], [18], the numerical error generally increases with damping ratio in nonlinear systems ( Figure 5). It is apparent that the originally suggested h/T = 0.2 by Newmark gravely errs in estimating both the peak response and residual displacement moderately damped systems. Moreover, the commonly used h/T = 0.10 develops considerable error in estimating the response of short and long period systems and is generally not accepted for nonlinear systems. This emphasizes the inappropriateness of using a universal approach of time-step selection as a constant fraction of the natural period of the structure. For heavy damping (20%), the time step should be limited to 0.01T for short period systems and 0.005T for medium and long period systems for acceptable numerical error (<1%). In medium damping range, the time-step should be limited to 0.005T, 0.01T and 0.02T for long, medium, and short natural period, respectively. For lightly damped systems, the time step is limited to 0.01T for acceptable numerical solution for all period ranges.

Earthquake excitation
In performing response history analysis for systems subjected to earthquake excitation, a very fine time step is typically used. The time step of the input acceleration record is very small (usually 0.01 or 0.02 sec.) to enable describing the highly jagged and irregular ground motion. This value normally restricts the time step used in performing the numerical analysis of the ODE of motion.
In the input ground motion record used in this manuscript, namely El Centro 1940 record, the time step of the input record is 0.01 sec. This value corresponds to 0.02, 0.01, and 0.005 times the natural period of short, medium, and long period, respectively. However, it was decided to investigate a wider range of time steps (0.005 to 0.20 sec.) to assess the associated numerical error and explore the possibility of using different time steps. Figure 6 displays the displacement response history of nonlinear SDOF system subjected to El Centro ground motion for 2% damping, obtained using Newmark-AAM. Period elongation is not as evident in the case of ground motion in nonlinear systems, especially in forced vibration phase. For lightly damped systems, the numerical error increases with increasing the fundamental period of the structure for the same time step-period ratio. However, the relationship is not as clear for damping ratios other than 2%. In general, the numerical error tends to increase with increasing damping ratio for all period ranges.
The typical time step used in earthquake simulation, which is equal to or smaller than input record time step, yields very satisfactory results as can be clearly observed in Figure 6. Using a larger time step than the input ground motion one seems logically unacceptable due to the chance of missing a peak acceleration value, especially for nonlinear systems. However, slightly relaxing the time step in linear systems was possible [18]. For inelastic systems, a more stringent time step assumption seems to be necessary, especially with the increasing error due to high damping. The h/T corresponding to input ground motion record is 0.02, 0.01, and 0.005 for short, medium and long natural period, respectively. Accordingly, it is recommended, based on Figures 13, 14, and 15, to limit the time step to 0.005T for 20% damping for all period ranges, to 0.01T for 5% damping for all period ranges, and relax it to 0.02T for medium and short period ranges for lightly damped systems.

HHT method
Half-sine pulse excitation Figure 7 presents displacement response history of nonlinear SDOF system under half-sine pulse excitation obtained using HHT method for 2% damping ratio. Similar to AAM, the numerical solution tends in general to underestimate the response. The error in peak inelastic displacement and free vibration response are close in value. Consistent with HHT method results for linear system [17], [18]      relationship of the numerical error and the fundamental natural period of the systems is not clear as can be observed from Figure 7. Period elongation of the numerical scheme is evident. This period elongation decreases with increasing natural period of the structure. Period elongation is proportional to the time step selected. Figure 12 shows the effect of time step on the accuracy of HHT method for nonlinear SDOF system under half-sine pulse excitation. It can be observed that the numerical error increases with increasing damping ratio in nonlinear systems. Time step should be limited to 0.05T, 0.02T, and 0.01T for light, medium, and heavy damping nonlinear systems, respectively. Figure 8 displays the displacement response history of nonlinear SDOF system subjected to harmonic excitation for 2% damping and 20% damping using HHT method. The numerical error increases with increasing natural period in the free vibration response. However, the relationship is not very clear in the forced vibration phase, although it is noticed that the short period system has the least error in the heavy damping case. The numerical error in both forced and free vibration phases increase with increasing damping ratio. Period elongation error is more pronounced than that in Newmark-AAM for nonlinear systems. The period elongation is clearer in the free vibration phase. Figures 9, 13, 14 and 15 show the numerical error in estimating peak displacement response of nonlinear SDOF system subjected to harmonic excitation using HHT method. It is clear that commonly used time step of 0.1T in linear systems is inappropriate for nonlinear systems. The figure suggests using a time step of 0.005T for heavily damped systems and 0.01T for moderately and lightly damped systems for all period ranges. Figure 10 displays the displacement response history of nonlinear SDOF system subjected to El Centro ground motion for 2% damping, obtained using HHT method. Consistent with nonlinear system harmonic and pulse excitation, period elongation is not as evident in the case of ground motion in nonlinear systems, especially in forced vibration phase. However, period elongation is more pronounced than that of Newmark-AAM for nonlinear systems. The numerical error increases with increasing the fundamental period of the structure for the same time step-period ratio. In general, the numerical error tends to increase with increasing damping ratio of the system for medium and long period ranges. Figure 16, 17 and 18 depict the errors developed using HHT method to solve the equation of motion for inelastic systems subjected to earthquake excitation. The time step of the input record is 0.01 sec, which corresponds to 0.02, 0.01, and 0.005 times the natural period of short, medium, and long period, respectively. For short period range, the h/T value of 0.01 seems appropriate for light and medium damping and the value of 0.005 is reasonable for the heavily damped systems. For the medium period range, h/T should be limited to 0.01 for heavily damped systems and 0.005 for light and medium damping. For long period range, h/T ratio of 0.01, which is twice that of input ground motion, is appropriate for all damping ratios for accurate numerical results.   Figure 11 shows the displacement response of  and 2% damping earthquake excitation for long period range. Given these two exceptional cases and the very fine time step used, it is recommended that HHT1 method not used for inelastic systems and should be limited to linear systems.

Analysis of numerical schemes errors and usage recommendations
The purpose of this section is to suggest, based on observations on numerical error trends obtained throughout the numerical experiment performed in this study, the most suitable numerical scheme(s) for each loading condition, damping ratio and natural period for nonlinear SDOF systems. Figure 12 through 18 exhibit comparisons between numerical errors developed using the Newmark AAM and HHT methods. Based on observing error trends, Table 1 was derived to suggest the recommended maximum time step for an acceptable numerical error less than 1% in earthquake engineering applications.  Figure 10. Displacement response history for HHT method of nonlinear SDOF: El Centro earthquake.

Conclusions
The numerical experiment performed in the current study confirmed some previously known aspects of the accuracy of commonly used numerical schemes in the field of nonlinear structural dynamic analysis in earthquake engineering applications. It also unveiled some new findings regarding the effect of natural period, damping ratio, and nature of exciting force on the numerical accuracy of the investigated numerical schemes. Accordingly, the most appropriate numerical scheme and time step choice for reasonable accuracy were recommended for each application. Within the boundaries of the investigated system and parameters,   Figure 11. Displacement response history for HHT1 method of nonlinear SDOF.
forces, materials and methods, the following concluding remarks can be drawn: (1) The results of the current investigation clearly show the inappropriateness of the universal approach of selecting time step as a function of natural period independent of the value of that period. The time step of 0.1T recommended in many references is suitable only for some exceptional cases and cannot be generalized due to the significant associated error. The numerical error was proven to be sensitive to natural period range and structural damping as well as the exciting force and the numerical scheme itself.  Figure 13. Assessment of numerical method accuracy for NL SDOF, 2% damping: harmonic excit.
complete out of phase response. This elongation is more pronounced with short period and lightly damped systems. The period elongation of HHT method is more pronounced than that of Newmark AAM method for forced harmonic excitation. (5) Period elongation of HHT method is more pronounced than that of Newmark AAM in nonlinear systems.  Figure 15. Assessment of numerical method accuracy for NL SDOF, 20% damping: harmonic excit.
cases with extremely small time step. Accordingly, the use of HHT1 in nonlinear system is not recommended. (8) In AAM and HHT methods for nonlinear system, peak displacement error increases with increasing natural period of the system even for the same h/T ratio. However, the increase is less significant than that of linear systems. (9) For earthquake excitation in nonlinear systems, the choice of time step not larger than that of the input record is strictly required to obtain accurate results. Few exceptions exist where twice the time step of the input record yielded accurate results.  Figure 17. Assessment of numerical method accuracy for NL 5% damping SDOF: El Centro GM.