A review of automatic time-stepping strategies on numerical time integration for structural dynamics analysis

Engineering Structures 80 (2014) 118–136

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Review article

A review of automatic time-stepping strategies on numerical time integration for structural dynamics analysis Diogo Folador Rossi a, Walnório Graça Ferreira a,⇑, Webe João Mansur b, Adenilcia Fernanda Grobério Calenzani a a b

Department of Civil Engineering, Universidade Federal do Espírito Santo, Vitória, ES, Brazil Program of Civil Engineering, Instituto Alberto Luiz Coimbra de Pós Graduação e Pesquisa de Engenharia, Universidade Federal do Rio de Janeiro, Rio de Janeiro, RJ, Brazil

a r t i c l e

i n f o

Article history: Received 10 November 2013 Revised 7 August 2014 Accepted 8 August 2014

Keywords: Dynamic structural analysis Automatic time stepping Direct time integration methods

a b s t r a c t This paper discusses some of the algorithms available for the automatic adaptive selection of time step size, applied to the step-by-step direct time integration methods of structural dynamics problems. Three adaptive strategies based on different concepts are explored and compared: the algorithm of Bergan and Mollestad (1985), which is based on the ‘current characteristic frequency’; the strategy of Hulbert and Jang (1995), which uses a ‘local error estimator’; and the method of Lages et al. (2013), which is based on the ‘geometric indicator of displacements history curvature’. The reviewed strategies are applied to the Newmark integration scheme to solve various numerical examples of linear dynamic systems, which are presented to compare the performance between the three algorithms that are tested. To conclude, a brief analysis about the considerations of the computational cost is made. Ó 2014 Elsevier Ltd. All rights reserved.

Contents 1. 2.

3.

4.

5. 6.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Basic structural dynamics concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Dynamic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Newmark formulae . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Automatic time-stepping strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Strategy of Bergan and Mollestad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Strategy of Hulbert and Jang . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1. Normalised error and error tolerance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2. Algorithm for automatic time step size control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3. Application to problems with quiescent start conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. Strategy of Lages et al. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1. Indicator of curvature and time step strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2. Curvature regularisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical examples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Single-degree-of-freedom systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. Shear building with three floors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. Simply supported beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. A three-bar plane frame structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Computational cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

119 119 120 120 120 121 121 121 122 122 122 123 123 124 124 126 129 132 134 135

⇑ Corresponding author. Address: Av. Fernando Ferrari, 514, Goiabeiras, 29075-910 Vitória, ES, Brazil. Tel.: +55 27 40092666, cell: +55 27 99699645. E-mail addresses: [email protected] (D.F. Rossi), [email protected] (W.G. Ferreira), [email protected] (W.J. Mansur), [email protected] (A.F.G. Calenzani). http://dx.doi.org/10.1016/j.engstruct.2014.08.016 0141-0296/Ó 2014 Elsevier Ltd. All rights reserved.

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

119

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

1. Introduction It is well known that there are two general approaches that can be used to analyse the dynamic responses of structural systems. The first approach employs mode-superposition methods in which the total response is expressed as the sum of individual responses in the various normal modes of vibration. This approach could be useful for finding analytical solutions for the dynamic problem, but it cannot be applied in non-linear systems because of the superposition scheme involved, which typically requires that the system remains linear during the process. The second general approach that is applicable in the analysis of an arbitrary set of nonlinear dynamic equations, which is also effective for addressing coupled linear modal equations, is the numerical step-by-step direct time integration. These schemes are often used in practice because they can be easily applied to solve the cases of both nonlinear and linear systems numerically and because the analytical solutions for most of the real dynamic problems do not exist or have an excessive computational cost. To compute the numerical solution at specific time ti (where i is an integer that denotes the number of discrete time steps taken to reach ti), most direct time integration methods require the solution to be specified in the previous instant, ti1. Also required is the specified time step size, Dti = ti  ti1, among other parameter specifications for each specific method. It is also well known that it is difficult to choose the appropriate time step size for the direct time integration process. Selecting a small time step increases the number of evaluations in the dynamic analysis and provides a more detailed description of the response; however, selecting a small time step increases the computational cost and the time spent on the integration process. The correct choice of the time step size must consider the conflict between the demands of solution accuracy and computational cost (the stability is also important in some cases). The ideal time step size maximises the accuracy while minimising the computational effort. However, the optimal time step may change during the computation due to changes in the forcing function and/or system nonlinearities. Thus, although a fixed time step size is typically used for time domain analyses and its choice is frequently based on intuition and/or experience, the best approach to address this parameter is employing an algorithm for automatic time step selection. The algorithm should typically seek the largest possible time step while maintaining a prescribed accuracy. It should also be efficient and economical and satisfy some specific criteria. Several automatic time-stepping algorithms for dynamic problems have been studied in past research and have been described in the literature. Some of the early methods were proposed by Hibbitt and Karlsson [1], who developed a strategy based on a measure of the residual force vector computed at the midpoint of the time interval, and by Oughourlian and Powell [2], who proposed a cheaper method to obtain the mid-step force residual through the use of the product between the stiffness matrix and the incremental velocity vector. However, both methods appeared to have too significant additional computational cost. Felippa and Park [3], Park and Underwood [4], and Underwood and Park [5] proposed an algorithm for step size control that was applied to the central difference method based on the ‘apparent highest frequency’, which is computed from the highest ratio between incremental accelerations and displacements for all degrees of freedom. Following the concept of apparent highest

frequency, Bergan and Mollestad [6] proposed an automatic time-stepping algorithm based on a ‘current characteristic frequency’ that is estimated using an expression similar to the Rayleigh quotient. Their scheme adjusted the time step according to a period measure that is tied to the current stiffness and response. To prevent the step size from changing too frequently, they also introduced a ‘tuning function’ that determines the sensitivity of the time-stepping algorithm. Zienkiewicz et al. [7] proposed a different approach to time step length selection based on a simple expression as an indicator of local error, which was obtained at little computational cost. Observing that this indicator does not provide an accurate measure of the real local error, Zienkiewicz and Xie [8] proposed another method that has a more accurate error estimate. Zeng et al. [9] presented an identical error indicator, which has a slightly simple and intuitive derivation scheme. Li et al. [10] developed a new error estimator that was derived from the difference in the solutions between an ordinary integration method and another higher-order method, which assumes that the derivatives of the accelerations vary linearly within each time step. The last three methods cited here were applied on the ordinary Newmark [11] time integration process using a posteriori error estimators and strategies that control the time step size by decreasing, maintaining, or increasing it in such a way that the local error fell inside or outside the lower and upper limits of a specified tolerance interval. Based on the same approach of local error estimates, Hulbert and Jang [12] proposed an automatic time-stepping algorithm that was applied to the generalised-a method integration process, developed by Chung and Hulbert [13], which includes the Newmark [11], HHT-a [14], and WBZ-a [15] integration methods as special cases. Later, Chung et al. [16] developed an a priori error estimator for the generalised-a method, which is computed using only information in the previous and current time steps (without a feedback process, which is required in most conventional a posteriori error estimators). Recently, Lages et al. [17] introduced a new automatic time step strategy that employs an estimator based on the geometric indicator of displacement history curvature. The developed estimator is obtained at little computational effort and is easily applicable to the various existing direct time integration methods. In this paper, three specific automatic time-stepping strategies are explored in more detail and compared with each other. Specifically, the algorithms of Bergan and Mollestad [6], Hulbert and Jang [12], and Lages et al. [17] are discussed when applied to the time integration process of the Newmark scheme. These strategies were chosen because they are based on three different concepts: ‘current characteristic frequency’, ‘local error estimator’, and ‘indicator of curvature’, respectively. The reviewed strategies are implemented with some necessary modifications, and the compared results are presented in numerical examples of linear structures, including single-degree-of-freedom and multi-degree-of-freedom systems.

2. Basic structural dynamics concepts This section describes the basic concepts of structural dynamics adapted for step-by-step procedures and the equations of the Newmark integration method.

120

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

Tj Dt 6 Dtcr ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p 1  4b

2.1. Dynamic equations After spatial discretisation, a structural dynamic system can be represented by the following matrix system of second-order linear ordinary differential equations:

€ þ CðtÞDðtÞ _ þ KðtÞDðtÞ ¼ PðtÞ MDðtÞ

ð1Þ

where M is the mass matrix, C(t) and K(t) are the viscous damping and stiffness matrices at time t (which might be constants in a lin_ € and DðtÞ are the displacement, velocity, and ear system), DðtÞ, DðtÞ, acceleration vectors at time t, respectively (superposed dots indicate differentiation with respect to time), and P(t) is the load vector at time t. The initial value problem consists of finding a function that satisfies (1) for all t 2 [0, tN], tN > 0, and the known initial _ _ 0. conditions: D(0) = D0, Dð0Þ ¼D The purpose of all numerical direct time integration algorithms _ i , and D € i , which are given approximations of D(ti), is to find Di, D € i Þ, respectively, in discrete instants ti. The most gen_ i Þ, and Dðt Dðt eral procedure to reach this aim is to consider the following incremental equation in a generic step i on the step-by-step time integration [18,19]:

€ i þ Ci DD_ i þ Ki DDi ¼ DPi M DD

ð2Þ

_ i , and DD € i are the displacement, velocity, and accelwhere DDi, DD eration vector increments, respectively, that occur from the beginning to the end of the time step, DPi = P(ti)  P(ti1), and Ci and Ki are the current incremental (initial tangent) damping and stiffness matrices, respectively, which are assumed to be constants for the short increment of time or deformation within the time step for nonlinear systems analysis. Linear systems can also be represented by this equation, which becomes simplified when the matrices Ci and Ki remain constant over the entire analysis. After the incremental quantities have been determined with some direct integration procedure, the displacement and velocity vectors at the end of the time increment can be evaluated by using the following:

Di ¼ Di1 þ DDi

_ i ¼ D_ i1 þ DD_ i D

ð3Þ

The acceleration at the end of the time step could also be found € i , but to avoid an accumulation of errors in using the increment DD a nonlinear analysis, the initial acceleration vector is calculated directly from the equilibrium condition at the beginning of the following step:

€ i1 ¼ M1 ðPi1  Ci D _ i1  Ki Di1 Þ D

ð4Þ

where the inverse of the mass matrix, M1, is calculated once because it is typically constant. 2.2. Newmark formulae The integration method of Newmark [11] is one of the most popular and widely used time integration algorithms in dynamic analysis. This method uses the following implicit expressions as assumptions over the displacement and velocity increments:

€ i1 þ cDti D €i DD_ i ¼ ð1  cÞDti D

ð5Þ

  1 € i1 þ bDt2 D €i DDi ¼ Dti D_ i1 þ  b Dt2i D i 2

ð6Þ

where c and b are parameters chosen by the analyst to control the stability and accuracy. The method is unstable for c < 1/2 [20]. When c = 1/2, the method is unconditionally stable for b P 1/4 [21]. For b < 1/4 and c = 1/2, it is only conditionally stable, and the condition for stability is given by the critical time step size [22]:

ð7Þ

where Tj is the period of natural vibration for the jth normal mode. This stability condition must be satisfied for each mode in the system; thus, the minimum natural period, Tmin, is the critical value. When the incremental dynamic equilibrium Eq. (2) is considered with Eqs. (5) and (6), the following system of simultaneous linear equations (which has the form of an incremental pseudostatic equilibrium equation) can be obtained:

KDDi ¼ DP

ð8Þ

where K and DP are the effective stiffness matrix and incremental effective load vector, respectively, given by

K¼

1 c Mþ Ci þ Ki bDt i bDt2i

  1 _ 1 € Di1 þ Di1 DP ¼ DP þ M bDt i 2b    c € i1 þ c D _ i1  1 Dt i D þ Ci 2b b

ð9Þ

ð10Þ

The solution of system (8) provides the displacement increment, and with this result, the velocity increment may be found using the following equation (which does not require the € i ), which is obtained by combining (5) and (6): value of D

DD_ i ¼

  c € i1  c D_ i1 DDi þ 1  Dt i D bDt i 2b b

c

ð11Þ

3. Automatic time-stepping strategies Three adaptive time-stepping procedures will be discussed in the following sections, which were chosen because they are based on three different concepts. First, it is important to define some characteristics that make an automatic time-stepping algorithm efficient, economic, and applicable in real structural dynamic analyses. Bergan and Mollestad [6] suggested a set of criteria to evaluate step size control strategies. Briefly, these criteria are: 1. The optimised time step size should not be influenced by the initial step size estimate. 2. The time step should remain constant during a linear stationary response. 3. The time step should not be influenced by the selection of physical units or by the number of degrees of freedom in the dynamic equations. 4. All of the input parameters should be easy to describe. 5. The additional step size selection computational cost should be as small as possible. 6. The step size should react immediately to sudden changes in a dynamic response. 7. The time step should not change without necessity. The last point above has a special importance in the integration process. Indeed, a change in the time step of a linear system analysis means that a new effective stiffness matrix K (given by Eq. (9)) must be formed and factorised to solve Eq. (8) and to find DDi, which is an expensive operation. Thus, the time step should be changed only when necessary, and only major changes in the response should lead to adjustments in the time step size for linear problems. Too frequent changes in the time step might not introduce additional computational cost only in nonlinear analysis, in which the tangent damping and stiffness matrix (Ci and Ki) are recomputed for almost every step because of nonlinear effects,

121

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

thus changing the matrix K. In this paper, only linear cases are compared in numerical examples. The three strategies described in the following sections have procedures to avoid frequent changes in the time step size when it is not necessary. 3.1. Strategy of Bergan and Mollestad The strategy of Bergan and Mollestad [6] proposes the use of an expression that is similar to the Rayleigh quotient, which is defined as a ‘current characteristic angular frequency’, given by the first of the following equations for step i:

 2i ¼ x

DDTi Ki DDi DDTi MDDi

2p T i ¼  1=2 x  2

ð12Þ

 2i ¼ x  2i1 when kDDi k < ekDDi1 k x

i

Based on this frequency, a corresponding measure of the ‘current characteristic period’ is calculated by the second of the previ 2i was taken because it might ous equations. The absolute value of x  2i is not be negative and T i must be real. The defined frequency x any specific natural eigenfrequency, but it reflects the incremental response for all eigenmodes, being an estimate for the dominant frequency of the response. The period T i expresses a dominant increment response period for the total system, and it is tied to the current stiffness and response, but it is not generally a measure of the time elapsed between zero amplitudes. Using this characteristic period, a preliminary estimate of the time increment for the next time step is calculated by

Dt iþ1 ¼ kT i

ð13Þ

where the asterisk indicates that this equation is only an estimate and not the actual new time step size. Here, k is an input constant that is denoted as the ‘step-length parameter’, which is selected according to the desired accuracy and could lie between 0.001 (high accuracy) and 0.1 (low accuracy). A typical value of k would be 0.05, which corresponds to 20 time steps for one characteristic period. This number is the fixed value used in the computing software ANSYS 11.0 [23], in which this method is applied. Using the estimate Dtiþ1 , the ‘current time-step ratio’,  ni , is defined by the first of the following Eq. (14), and the actual new time increment of the next time step is calculated by the second of these equations. 

ni ¼ Dt iþ1 Dt i

 Dt iþ1 ¼ f ni Dti

step length does not change unless a major change in the response has taken place. Outside this plateau, the function f  ni is simply  ni , which makes the setting of the time step in accordance with Eq. (13) except that a maximum increase has been set at  nm . Bergan and Mollestad reported that  np should be specified as an input parameter (as well as the step-length parameter, k, and the first time step size, Dt1) and that typical values of  np would be 1.3 or 1.4. The maximum increase in the time step,  nm , can be kept fixed (not an input parameter), and its value would be 2.0. Finally, it is worth noting that, in some steps, the values of DDi could become very small or vanish, such as when passing through  2i maximum amplitudes. In these situations, it was proposed that x could be assigned to its value in the previous step:

ð14Þ

 ni is a ‘tuning function’ of  ni that determines the sensitivwhere f  ity of the time-stepping algorithm. Bergan and Mollestad discussed some examples of tuning functions in their work and concluded that the best type of such functions has the form of the function shown in Fig. 1. ni ¼ 1, which This type of tuning function has a plateau around  is defined by a value  np and serves as a protection against making unnecessary adjustments of the time step, what means that the

ð15Þ

where e is a fixed small number, which is specified to be 0.1. Other comments were made in the original work about various aspects of the algorithm, but the description given here is necessary for implementation. 3.2. Strategy of Hulbert and Jang The adaptive strategy of Hulbert and Jang [12] proposes to use an algorithm based on the local error estimate approach. The local error ei is a measure of the error in the numerical solution for each time step. This error is obtained by the difference: ei ¼ Di  Dðti Þ, where DðtÞ is the solution of the local problem, i.e., it is an analytical function that satisfies the following equations:

€ þ Ci DðtÞ _ þ Ki DðtÞ ¼ PðtÞ MDðtÞ Dðti1 Þ ¼ Di1

_ i1 Þ ¼ D_ i1 Dðt

ð16Þ € i1 Þ ¼ D € i1 Dðt

ð17Þ

Using Taylor series expansions, Hulbert and Jang developed a simple expression to define the local error, which is applied to the generalised-a method integration process. This time integration algorithm uses the same Eqs. (5) and (6) of the Newmark method combined with a modified balance equation, in which two additional parameters are introduced: am and af. When am = af = 0, the algorithm reduces to the original Newmark method. Therefore, by making these two parameters equal to zero, the local error expression developed by Hulbert and Jang, applied to the Newmark integration process, is given by

  1 €i ei ¼ Dt 2i b  DD 6

ð18Þ

which is identical to the expression presented by Zienkiewicz and Xie [8] and Zeng et al. [9]. This equation is invalid only if b = 1/6 (ei is equal to zero). However, for this case, an error estimate can be found by using high-order terms in the Taylor series expansion [10]. In this paper, only Eq. (18) will be used to define the local error of this strategy. 3.2.1. Normalised error and error tolerance To specify an appropriate tolerance on the local error for choosing the time step sizes, Hulbert and Jang proposed to relate it to an input quantity, (Dt/T)target, for which the relationship with the error tolerance is obtained by considering an unforced, undamped single-degree-of-freedom model, where T denotes the period of vibration. To obtain such a relation, the norm of the local error was scaled using the following equation:

RLi ¼

Fig. 1. Tuning function for the time-stepping algorithm of Bergan and Mollestad.

kei k sclfaci

  9 sclfaci ¼ max kDDi k; sclfaci1 10

ð19Þ

where sclfac0 = 0. The resulting discussion of this strategy applied to the Newmark scheme demonstrates that given the user specified

122

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

parameter (Dt/T)target, the corresponding tolerance for this normalised error can be computed by

tol ¼ C d

 2 Dt T

   1 C d ¼ ð2pÞ2 b   6

ð20Þ

Common practice dictates that there should be a minimum of 10 time steps per period of the maximum frequency of interest in the response [12], which suggests that typical values for (Dt/T)target would be 0.1. The original work of Hulbert and Jang also suggests another scheme for the normalised error RLi. However, the expression of the tolerance in this case was derived numerically by curve fitting, and the resulting equation is a function of the ‘high-frequency dissipation parameter’ of the generalised-a method, q1, which has no relationship to the Newmark method. For this reason, only the normalisation scheme of Eq. (19) will be used in this paper. 3.2.2. Algorithm for automatic time step size control After calculating the error tolerance by Eq. (20) and the normalised error RLi for the current interval, the first procedure of the time step control algorithm is to request, at each step, that

lb  tol 6 RLi 6 tol

ð21Þ

where lb is a multiplication factor for the lower-bound tolerance, for which a typical fixed value would be 0.75. Three situations can occur from this point: (1) if the condition (21) is satisfied, then the solution is accepted and the integration process proceeds to the next time increment without changing the step size; (2) if RLi < lb  tol, then the time step size can increase because the error is considered to be too small. However, this increase in step size shall occur only for the next time increment, and the solution for the current step is taken to be acceptable. Due to the periodic nature of the local error, the time step length would also be periodic, and because it is not desirable to change the time step too frequently, a counter ‘cont’ was introduced to register the number of steps in which the situation ‘RLi < lb  tol’ occurs consecutively. If this counter exceeds a specified value, ‘lcount’, only then the step size is increased for the next time increment, in accordance with the following equations:

Dtiþ1 ¼ f inc  Dt i

f inc ¼



tol RLmax i

1=pinc ð22Þ

where pinc is an input parameter that is related to the rate of convergence of the local error and RLmax is the maximum value of RL that i has occurred since the counter was reset to zero. The counter is reset every time that a step size change occurs (an increase or decrease) or when the condition (21) is satisfied. The remaining situation that could occur is: (3) RLi > tol. When this situation occurs, where the error is computed using a current step size, Dt old i , then the time increment can be reduced: the current solution is discarded, and a new solution is computed using a smaller step size, Dtnew . Two cases can occur when this step size reduci tion is required. In the first case, if Dt was increased in the previous step, then the new trial Dt is taken as its value before the step size increase has occurred. In the other case, i.e., if the step size did not increase in the previous step (e.g., if it did not change or decreased), then the new step size is computed using the following equations:

Dtnew ¼ f dec  Dtold i i

 f dec ¼

tol RLi

Finally, it is important to expose that the user-specified value (Dt/T)target can be used to determine the number lcount. Hulbert and Jang applied a criterion where this number is defined as the reciprocal of (Dt/T)target, which can be stated as follows:

1 lcount ¼ int ðDt=TÞtarget

!

3.2.3. Application to problems with quiescent start conditions The time-step control strategy described previously cannot be used to compute the step size for the first time increment when the initial conditions and load are zero (quiescent conditions). The normalised error for the first time step would be independent of Dt for a quiescent start. It can be shown that under such conditions, the normalised error is given by

RL1 ¼

b  16 b

ð25Þ

Since b is an algorithm parameter, the time-step control strategy may not be utilised. Thus, the first step increment, Dt1, is not changed in a quiescent start problem, and this time step size should be given as an input parameter. For the second time step, Dt1 should be used as initial trial step size. Then, if RL2 > tol, the step size is decreased. However, due to quiescent start, the value ||DD2||, obtained with the decreased time step, may not be large enough to be used as a normalisation factor. For this reason, the normalisation factor is kept constant in the second time step with the value kD02  D1 k, where D02 is the displacement computed for time t2 using the initial step size estimate Dt1. 3.3. Strategy of Lages et al. The automatic adaptive strategy proposed by Lages et al. [17] is based on a geometric indicator defined as the curvature of the displacement history. Using the definitions of elementary vector calculus, the parametric representation of a curve in space, described by a given parameter t, can be expressed as a vector function rðtÞ (Fig. 2). The vector v ðt a Þ ¼ dr=dtjt¼ta – 0 is the ‘geometric velocity vector’ of the curve r at time ta, and it is tangent to the curve at the point r(ta). The corresponding ‘unit tangent vector’ is defined by TðtÞ ¼ v ðtÞ=kv ðtÞk, and the ‘arc length function’ s(t) is the integral of the norm ||v(t)|| taken from a specified instant point t0 to a generic instant t. Based on these concepts, the ‘curvature of a curve’ at time t, which is denoted by j(t), is defined by



dT

ðtÞ

ds

jðtÞ ¼

ð26Þ

Using definition (26), it is possible to show that the curvature can be expressed by the following equation:

jðtÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðv  v Þða  aÞ  ðv  aÞ2 ðv  v Þ3=2

1=pdec ð23Þ

where pdec has the same definition as pinc and the fixed value used for these parameters was pinc = pdec = 2. As cited before, the solution at time t i ¼ t i1 þ Dt old is discarded, and a new solution is obtained i by using Dt new , thereby defining a new value of RLi and restarting i the adaptive strategy for that step.

ð24Þ

Fig. 2. Parametric representation of a curve in space.

ð27Þ

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

where a(t) = dv/dt = d2r/dt2 is the ‘geometric acceleration vector’ and the symbol ‘‘ ’’ indicates the dot product of two vectors. The benefit of using Eq. (27) is related to the computational economy, which is the reason why it is recommended. The function jðtÞ gives a positive measure about the geometry of rðtÞ that indicates how much it is curved. The curvature will be close to zero at points where rðtÞ is similar to a straight line. At the other points, jðtÞ will increase or decrease in accordance with the deviation of the function rðtÞ from a straight line. This information is simple, but it can be used in an efficient time adaptive strategy.

3.3.1. Indicator of curvature and time step strategy Let r(t) describe the displacement history of a structure by the first of the following equations. Thereby, r(t) will be a vector of n + 1 elements, where n is the number of degrees of freedom in the structure. Thus, v(t) and a(t) will be given by the corresponding temporal derivatives:

rðtÞ ¼



t



vðtÞ ¼

DðtÞ



1 _ DðtÞ



aðtÞ ¼



0 € DðtÞ

 ð28Þ

The indicator of ‘displacement history curvature’ should then be obtained by an expression that is similar to Eq. (27). Using the vector functions (28), the curvature should be given by (for a discrete time ti)

ji ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _iD _ i ÞðD €i  D € i Þ  ðD _iD € i Þ2 ð1 þ D 3=2 ð1 þ D_ i  D_ i Þ

ð29Þ

This curvature could be used to control the time step size. As cited before, ji will be small when the graphs of D are similar to a straight line and larger when the graphs are more curved. It is well known that plotting a response history in locations that are similar to a straight line requires few points to describe the graph, whereas more points are required in parts that have greater curvature. This fact suggests that a smaller step size should be used in parts of the response history that have a greater curvature, where a more detailed solution is needed to guarantee the accuracy. Similarly, a larger step size can be used in parts that have a smaller curvature. Because the curvature is a positive measure, some relationships between it and the time step size can be established. Lages et al. [17] suggested a simple exponential expression to find the time step for the next increment, which is given by

Dt iþ1 ¼ Dtmax expðcp ji Þ

123

Fig. 3. Intervals of curvature regularisation.

3.3.2. Curvature regularisation The value of the time increment obtained by Eq. (30) or Eq. (31) should not be used in the next step solutions without some necessary adjustments. In many cases, the values of ji could present a highly oscillatory behaviour, which causes many changes in the time step size, and as discussed earlier, this scenario is not typically desirable. To prevent this situation from occurring, Lages et al. introduced a method of ‘curvature regularisation’, which makes the curvature values change in regular plateaus during the step-by-step analysis. Some techniques of curvature regularisation were suggested, but the most appropriate technique was the strategy based on the ‘maximum interval value’: during numerical integration, for each time instant, the curvature value is investigated and compared with the maximum value that occurs in the interval that corresponds to that instant; the highest value among them is chosen to be used in time step computation with Eq. (31). It was suggested that the regularisation interval (Dtreg) be a function of the critical time step size (Dtcr), which is a constant quantity. The next figures present the scope of the regularisation technique. Fig. 3 illustrates the step iterations, and Fig. 4 illustrates the programming process used in this paper.

ð30Þ

where cp is a positive constant and Dtmax is an upper value to the time step size, which can be set equal to the critical value in conditionally stable integration methods. This paper proposes another equation to relate the step size and curvature, which has the form

Dt iþ1

Dt max ¼ 1 þ cp ji

ð31Þ

In Eqs. (30) and (31), the time step varies from an upper specified value (Dtmax) to zero in the limit case when the curvature is infinity. Eq. (31) appeared to be more appropriate than Eq. (30) because the exponential function varies more than this ratio form for the same variation in curvature. Moreover, the exponential form yields a zero time step for a curvature value that is not overly high, i.e., the number exp(ji) is close to infinity in numerical calculations when ji is close to 700, which causes a zero time step size in Eq. (30). For this reason, Eq. (31) was then chosen for use in the numerical examples of this paper.

Fig. 4. Algorithm for regularisation in a maximum interval value, where: Dt reg ) Regularisation interval, ct ) Positive constant, max jðt 0j ; ti Þ ) Maximum value of ji between times t0j and ti , a ) Weighting constant.

Fig. 5. Single-degree-of-freedom system.

124

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

displacements; a simply supported beam with a harmonic load; and a plane frame structure composed of three bars with different applied loads.

Table 1 Set of parameters. (1)

(2)

m k c D0 D_ 0

1.0 1.0 0.6 0 0

1.0 p2/4 0 1.0 0

p

1.0

0

4.1. Single-degree-of-freedom systems

The presented algorithm refers to the maximum curvature value that has occurred in the previous time interval (jreg i1 on the eighth line of Fig. 4). This process is used to minimise time step changes in the first instants of the regularisation interval, Dtreg, where the maximum value between the times t0j and ti, max jðt0j ; t i Þ, does not accurately represent the behaviour of function j(t) when ti is close to t 0j , which could cause undesired oscillating changes in the time step size. This technique also guarantees that the regularised curvature values are always greater than the non-regularised values, which is convenient. The values of the parameters a, ct, and cp will be chosen for each case separately because the values suggested in the original work that proposes this strategy did not yield reasonable results when applied in the dynamic problems of frame structures, which is treated in this paper. The first time step for problems with quiescent start conditions, Dt1, could also not be changed because the initial curvature, j0, is zero in such cases, which would make the first time step be set equal to the maximum value, Dtmax. Normally, for transient problems, the most important dynamic events occur at the very beginning; as a result, it is not convenient to choose an excessively long first time step. Therefore, this first step size should be given as an input parameter in quiescent start problems in the same way as in the strategy of Hulbert and Jang. 4. Numerical examples This section presents a comparison of the results that are given by the three previous related strategies when applied to simple numerical analyses of dynamic structural systems. The same integration process (Newmark) was used with the three techniques. The first example studies the model of a single-degree-of-freedom system with two sets of parameters. Then, three different multi-degree-of-freedom problems are presented: a shear building structure that assumes a lumped mass system with only lateral

The spring-mass system shown in Fig. 5 is analysed with two different sets of parameters, which are given in Table 1. The first set of parameters is a forced vibration system, whereas the second is an undamped free vibration system. The analytical solutions for both cases are well known from the elementary theory of mechanical vibration. For the Newmark integration process, c = 1/2 was taken in both cases. The value of b was chosen as 1/8 and 1/5 for the first and second sets of parameters, respectively. For the damped system, the first time step was chosen as Dt1 = 0.5 s, the value of (Dt/T)target in the strategy of Hulbert and Jang was taken as 1/10, and the parameters cp, ct, a for the method of Lages et al. were set equal to 40, 1.0, and 0.2, respectively. For the undamped system, the first time step was chosen as Dt1 = 0.01 s, the value of (Dt/T)target was chosen as 1/100, and the parameters cp, ct, and a were set to 40, 1.0, and 1.0, respectively. The algorithm of Bergan and Mollestad is not compared here because it does not apply to single-degreeof-freedom systems. The value Dtmax used in the method of Lages et al. was taken to be equalptoffiffiffi the critical time step pffiffiffi (Dtcr) in both cases: this increment is 2 2 ¼ 2:8284 s and 4 5=p ¼ 2:8470 s for the first and second cases, respectively. To compare the accuracy of the numerical solutions, the ‘absolute error’ (er) was calculated from the difference between the numerical response ui and the exact solution uex(ti) found by analytical equations. The absolute error history is shown in Fig. 6(a), and the time step size histories for the different adaptive strategies are shown in Fig. 6(b), using the first case for the parameters. The numerical error is greater when a constant time step is used compared to when the adaptive methods are applied on the integration process. The highest values of this error occurred in the first stages of the analysis, which are exactly when the most important dynamic events occur, before the solution is dissipated by damping. Both strategies of Hulbert and Jang and Lages et al. start by adjusting the first time increment, Dt1, in accordance with the local error tolerance and initial curvature, respectively. The strategy of Lages et al. yields the smallest error at the beginning, which grows as the solution advances, becoming greater than the other errors in the final instants. The strategy of Hulbert and Jang presents a greater error at the beginning. This behaviour is 2.5

0.012 Constant Δt Hulbert & Jang Lages et al.

2

Time step size history (Δt)

0.010

0.008

Error (er)

Constant Δt Hulbert & Jang Lages et al.

0.006

0.004

0.002

0.000

1.5

1

0.5

0 0

4

8

12

16

20

0

4

8

12

Time (t)

Time (t)

(a)

(b)

Fig. 6. Absolute error and time step size history for the first set of parameters.

16

20

125

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136 1.000

0.018

RL

0.001

0.015

0.0008

0.012

0.0006

0.009

0.0004

0.006

0.0002

0.003

0.900

Non-Regularized Curvature Regularized Curvature

0.800 0.700

Curvature (κ)

|| e ||

Normalized local error history (RL)

Norm of the local error history (|| e ||)

0.0012

0.600 0.500 0.400 0.300 0.200 0.100

0

0.000

0

0

4

8

12

16

0

20

4

8

12

16

20

Time (t)

Time (t)

(a)

(b)

Fig. 7. Normalised local error history and curvature history for the first set of parameters.

3.0

0.06

0.05

Regularized Curvature

2.5

0.04

Curvature (κ)

Time step size history ( Δt)

Non-Regularized Curvature

Constant Δt Hulbert & Jang Lages et al.

0.03

2.0

1.5

0.02

1.0

0.01

0.5

0.0

0 0

5

10

15

20

0

5

10

15

Time (t)

Time (t)

(a)

(b)

20

Fig. 8. Time step size history and curvature history for the second set of parameters.

1.8x10 RL

-4

1x10

-5

8x10

-6

1.2x10

6x10

-6

9.0x10

4x10

-6

1.5x10

-4

-5

-5

6.0x10

-5

-6

0x10

3.0x10

0

0

0.0x10 0

5

10

15

20

Constant Δt Hulbert & Jang Lages et al.

0.006 0.005

Error (er)

|| e ||

2x10

0.007

-4

-5

Normalized local error history (RL)

Norm of the local error history (|| e ||)

1x10

0.004 0.003 0.002 0.001 0

0

5

10

Time (t)

Time (t)

(a)

(b)

15

20

Fig. 9. Normalised local error and absolute error for the second set of parameters.

in accordance with the variation in the time steps. As shown in Fig. 6(b), for the strategy of Hulbert and Jang, the time steps start at close to half of the initial value (Dt1) and tend to increase as the analysis advances, which is in accordance with the behaviour of the normalised local error. For a short interval of time, the value of Dt has been set to be equal to the initial estimate (Dt = 0.5 s)

before it is increased again, searching for stabilisation at a target value given by the input (Dt/T)target, which is approximately Dt = 0.628 s. In many steps, the normalised error RL lies below the lower tolerance limit lb  tol = 0.0123 in condition (21), as shown in Fig. 7(a). This observation explains why the time step increases with time. Fig. 7(a) also illustrates that the error RL never

126

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

Table 2 Summary of results of single-degree-of-freedom systems. Set

T

Dtcrit

(1) (2)

6.283 4.000

2.828 2.847

Hulbert and Jang

Table 3 Set of parameters.

Lages et al.

Dt min

Dtmax

NTOT

Dtmin

Dtmax

NTOT

0.2782 0.0100

0.8062 0.0400

47 576

0.0690 0.0286

2.2898 0.0286

85 701

exceeds the tolerance tol = 0.0164 that is predicted by the algorithm. This graph still shows the history of the local error norm, e, which is used in the strategy of Hulbert and Jang (the left vertical axis shows the values of e and the right axis shows the values of RL). Still observing Fig. 6(b), the time steps of the strategy of Lages et al. are initially considerably lower than the others, becoming greater when the solution advances over time. This pattern occurs because the curvature has a high initial value, which reduces substantially as the response is dissipated by the damping. Fig. 7(b) presents the curvature history, with and without regularisation, applied in the strategy of Lages et al. The regularised curvature history in a constant plateau is very sensitive to changes in the parameter a, which was taken as 0.2 in this case. The analysis of the second case also illustrates some disadvantages of the strategy of Lages et al. in this type of system. Fig. 8(a) presents the time step size history for the second set of parameters. The time steps in the strategy of Lages et al. remained constant and were equal to the step obtained with the maximum curvature, i.e., Dt = Dtcr/(1 + cpjmax). This pattern occurred because of the choice a = 1.0 in the curvature regularisation, which makes the regularised curvature assume a constant value equal to the maximum obtained from the displacement history, which is reached at the beginning, as shown in Fig. 8(b). This choice was made because the actual curvature assumes a permanent oscillatory behaviour, which results from the system being in an undamped free vibration (c = 0). Unless the regularisation interval was taken to be equal to half of the natural period, no value of a between zero and one would yield a reasonable behaviour of the regularised curvature that would transfer this permanent oscillation to the time step sizes. In the strategy of Hulbert and Jang, the time step was adjusted rapidly to reach the constant value given by the target ratio (Dt/T)target. On the first steps, the increment is relatively small, and the computed error RL is less than the lower bound of the tolerance; as a result, the step size remains small until the counter cont reaches the value lcount = 1/(Dt/T)target = 100, when the time step increases to the permanent value Dt = 0.04 s. This behaviour can also be seen in Fig. 9(a), which provides the error history of the algorithm. After reaching the permanent time step size, the norm of the local error (||e||) becomes oscillatory, and the normalised error (RL) is controlled within the tolerance interval, assuming small values periodically for a few steps and becoming equal to the upper bound value tol = 1.3  104 in almost every step.

m1 m2 m3 k1 k2 k3

(2)

150 t 150 t 150 t 100,000 kN/m 200,000 kN/m 300,000 kN/m

180 t 270 t 360 t 105,000 kN/m 210,000 kN/m 315,000 kN/m

Fig. 11. Applied loads on the shear building systems.

The absolute errors for this set of parameters are shown in Fig. 9(b). The strategy of Hulbert and Jang provided the largest error, whereas the method of Lages et al. assumed smaller error values. These results are in accordance with the time step size behaviour, where the strategy that used larger values of Dt reached larger errors. Table 2 summarizes the results of analyzes of the two singledegree-of-freedom systems with respect to maximum and minimum time increments found with adaptive strategies, and the total number of points (NTOT) plotted with the different algorithms. The higher the number NTOT, the more expensive is the strategy and the more detailed is the dynamic response. It is concluded that the strategy of Lages et al. provided more points in both analyzes. 4.2. Shear building with three floors To illustrate the solution of a dynamic analysis in a simple multi-degree-of-freedom system, the structure will be considered, in this example, to be a ‘shear building’ with three floors. In this type of structure, the supporting columns provide only horizontal displacements in each floor. Thus, the shear building has only three degrees of freedom, as shown in Fig. 10. The matrix M of this structure is obtained by the lumped mass model, where the masses of the floors are disposed in the first diagonal of the matrix. Two cases of shear buildings were tested in this example using the parameters provided in Table 3. In the first case, there is a transient horizontal load that acts on the lower floor with a constant value for 0.38 s, becoming zero after that instant, as shown in Fig. 11(a). In the second case, there is a horizontal load acting on the upper floor, with a triangular impulsive load history that has a duration of 0.05 s, which is illustrated in Fig. 11(b). With this information, the load vectors for the two cases can be given, respectively, by

pðtÞ ¼ ½ 0 0 p3 ðtÞ T

Fig. 10. Shear building model with three floors.

(1)

and pðtÞ ¼ ½ p1 ðtÞ 0 0 T

ð32Þ

When solving the linear analysis of these systems, we considered Rayleigh’s proportional damping with damping ratios of 15% in the 1st and 3rd modes of vibration for the first case and of 5% for the second case, also in the 1st and 3rd modes. For the first case, the mass matrix, stiffness matrix, and natural eigenfrequencies are formed and found, giving the following results (in tons, kN/m and rad/s, respectively):

127

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136 0.015

0.025

First Floor

Bergan & Mollestad

Time step size history - Δt (s)

Second Floor Third Floor

Displacement (m)

0.01

0.005

0

Hulbert & Jang Lages et al.

0.02

0.015

0.01

0.005

-0.005

0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

Time (s)

Time (s)

(a)

(b)

0.7

0.8

0.9

1.0

Fig. 12. Dynamic response history and variation of Dt for the first case of the shear building.

Table 4 Parameters of the adaptive strategies used in the shear building examples. Case

Bergan and Mollestad

(1) (2)

2

Hulbert and Jang

Lages et al.

k

 np

(Dt/T)target

cp

ct

a

0.05 0.05

1.4 1.3

1/10 1/20

2 25

1.0 2.0

0.6 0.5

150

6 m¼4 0

0 150

0

0

0

3

2

1

1

0

3

2

16:65

3

7 6 7 6 7 0 5 k ¼ 1  105 4 1 3 2 5 x ¼ 4 39:11 5 150 0 2 5 64:76 ð33Þ

The damping matrix c can then be formed with Rayleigh’s damping using the damping ratios, eigenfrequencies, and matrices m and k. The displacement responses for the three floors, obtained from the dynamic analysis using a constant time step Dt = 0.001 s, are shown in Fig. 12(a). The parameters of the numerical integration procedure were set equal to c = 1/2 and b = 1/8, and the parameters for the different adaptive strategies are provided in Table 4. The value Dtmax of the method of Lages et al. was set equal to the critical time step, Dtcr, in both of the p cases ffiffiffi on the shear building. In the first case, this value is Dtcr ¼ T J 2=p. Fig. 12(b) presents the variation in the time step sizes with all of the adaptive

procedures discussed earlier, for which the first time step was chosen as Dt1 = 0.01 s in the three strategies. All of the algorithms predicted an increase in the time step size value during the first 0.35 s, sharply reducing it near time t = 0.4 s and increasing it again in the subsequent moments. Because the value of RL1 is greater than the bound tol, the strategy of Hulbert and Jang initiates the integration process by adjusting the initial time step, and the same is done in the strategy of Lages et al. The algorithm of Bergan and Mollestad uses the given value Dt1 in the first step integration and reduces it immediately afterward. We observed a slight tendency of the time step to reach the value Dt = 0.015 s at the end of all of the analyses. The strategy of Lages et al. provided the lowest values of Dt, whereas the other strategies provided values that were closer to each other. Fig. 13(a) presents the curvature sequences used in the strategy of Lages et al. Fig. 13(b) presents the variations of the characteristic period T calculated with the strategy of Bergan and Mollestad, as well as the ratio Dt=k, which has plateaus of constant values. The three horizontal traced lines represent the values of the natural modal periods of vibration of the structure, which are equal to T1 = 0.377 s, T2 = 0.161 s, and T3 = 0.097 s. The values of T vary between the values of the fundamental period T1 and the mode period of the highest frequency, T3, passing through the value T2. To confirm that the algorithms are independent of the initial step size Dt1, the same analysis is conducted with a different value 0.4

15.0 Non-Regularized Curvature

T

(s)

Regularized Curvature

Current characteristic period -

Curvature -

κ (1/m)

12.0

9.0

6.0

3.0

0.3

0.2

0.1 T

Δt /λ 0

0.0 0.0

0.1

0.2

0.3

0.4

0.5

Time (s)

0.6

0.7

0.8

0.9

1.0

0.0

0.1

0.2

0.3

0.4

0.5

Time (s)

(a) Fig. 13. Curvatures and characteristic periods for the first case of the shear building.

(b)

0.6

0.7

0.8

0.9

1.0

128

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

Bergan & Mollestad

Hulbert & Jang

0.025

0.025

Δt1 = 0.01 Δt1 = 0.001

Δt1 = 0.001

Time step size history - Δt (s)

Time step size history - Δt (s)

Δt1 = 0.01

0.02

0.015

0.01

0.005

0.02

0.015

0.01

0.005

0

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.0

1.0

0.1

0.2

0.3

Time (s)

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Time (s)

(a)

(b) Fig. 14. Variation in the time step size using different values of Dt1.

0.035 Bergan & Mollestad

Norm of the local error history (|| e ||)

Time step size history - Δt (s)

8x10

Hulbert & Jang Lages et al.

0.03 0.025 0.02 0.015 0.01

0.005 0

|| e || RL

6x10

-7

3.00x10-3

4x10

-7

2.00x10

2x10

-7

10

0x10 0.0

0.5

1.0

1.5

4.00x10-3

-7

-3

0.00x100

0

0.0

2.

-3

0.5

1.0

1.5

Time (s)

Time (s)

(a)

(b)

2.0

Fig. 15. Time step size history and normalised error for the second case of the shear building.

1.8

0.5

Non-Regularized Curvature 1.6

Current characteristic period - T (s)

Regularized Curvature

Curvature - κ (1/m)

1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0

0.4

0.3

0.2

0.1 T

Δt /λ 0

0.0

0.5

1.0

1.5

2.0

Time (s)

0.0

0.5

1.0

1.5

2.0

Time (s)

(b)

(a) Fig. 16. Curvatures and characteristic periods for the second case of the shear building.

Dt1 = 0.001 s using the strategies of Bergan and Mollestad and Hulbert and Jang. The compared results are shown in Fig. 14. The behaviour of the values of Dt are generally insensitive toward selection of the initial time step in both methods. The strategy of

Lages et al. was not compared because it does not employ the value Dt1 unless the initial start conditions are quiescent. Fig. 14 indicates that the initial optimised time step would be between 0.005 s and 0.01 s in both strategies. This optimal time

129

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

oscillatory until approximately the same value was reached. This number is 1/20 of the average value T = 0.3 s, and the characteristic period of Bergan and Mollestad oscillates about this value, as shown in Fig. 16(b) (this result is compatible with the choice k = (Dt/T)target = 1/20). Fig. 16(b) also indicates that the values of T vary in the region defined by the fundamental period T1 = 0.439 s, and the period of the second mode T2 = 0.205 s, oscillating around the average value, which is close to 0.3 s. In the method of Lages et al., the time step starts smaller and increases in most parts of the analysis, following the inverse behaviour of the regularised curvature shown in Fig. 16(a). Fig. 17. Model of a simply supported beam under sinusoidal load.

4.3. Simply supported beam

increment is reached in the method of Bergan and Mollestad immediately after the first step with the two choices of Dt1. In the strategy of Hulbert and Jang, the initial value Dt1 = 0.01 s is adjusted before the analysis starts, as was noted before, but the step Dt1 = 0.001 s remains with its small value in the first steps until the value of RL stays below the limit lb  tol during lcount steps. Despite this different initial behaviour, the adjusted time step sizes stayed closer in the algorithm of Hulbert and Jang than in the strategy of Bergan and Mollestad. The procedure for defining the structural matrices and vectors are identical when solving the second case of the Shear Building, in which are applied the load p1(t) on the upper floor. The parameters for the integration method were chosen as c = 1/2 and b = 1/5. The results are shown in Figs. 15 and 16. The parameters for the different adaptive strategies are provided in Table 4. This case is a problem that has quiescent start conditions, and for this reason, all of the algorithms employ the initial increment Dt1 = 0.001 s in the integration of the first step. The error RL1 in the strategy of Hulbert and Jang is not used on the adaptive procedure and has a value that exceeds the limits of Fig. 15(b), which are obtained directly from Eq. (25). The initial curvature used in the method of Lages et al. is zero, and for this reason, the maximum time increment (Dtmax) is avoided in the first step by using the user-specified value Dt1. However, this procedure did not avoid the small value of the curvature on the second step, which made the increment Dt2 assume a relatively high value. In the strategy of Bergan and Mollestad, the time step size oscillated close to the value Dt = 0.015 s. In the algorithm of Hulbert and Jang, the increment increased in stages that were not very

A simply supported beam is excited in transverse motion with a harmonic point load acting at the midpoint of its length. Fig. 17 provides the problem parameters and geometry. The beam is modelled by 12 equal beam elements over its length. A Rayleigh proportional damping is introduced with damping ratios in the first and third modes of vibration set equal to 10% of the critical. Two cases were studied using different forcing frequencies and amplitudes. In the first case, the forcing frequency was equal to the fundamental frequency x1 = 2p rad/s and the loading amplitude was taken as p0 = 450 N; in the second case, the forcing frequency was equal to the third mode frequency x3 = 18p rad/s and the amplitude was chosen as p0 = 4,500 N. The beam was analysed with the Newmark procedure, choosing c = 1/2 and b = 1/4 and the initial time step Dt1 equal to 0.001 s. The parameters k and  np were set equal to 0.05 and 1.3, respectively, for the strategy of Bergan and Mollestad. The value (Dt/T)target was taken as 0.05 for the algorithm of Hulbert and Jang. The parameters cp, ct, and a for the method of Lages et al. were set equal to 1.0, 1.0, and 0.9, respectively, and the maximum value Dtmax = 0.1 s was employed, which corresponds to 1/10 of the fundamental period T1 = 1 s. The regularisation interval was calculated as Dt reg ¼ ct Dtmax because Dtcr = + 1 in this case. Along with the above-mentioned adaptive strategies, the same transient analysis was executed with the computer program ANSYS 11.0 [23], applying the automatic time-stepping procedure of the program. Fig. 18 presents the obtained results for the first mode-forcing frequency case. The responses are almost identical for all of the analysis. In the strategy of Bergan and Mollestad, the parameter k corresponds to an ideal time step of 5% of the fundamental period. Fig. 18(b)

Displacement history of the central node

Time step size history

0.08

0.10 Constant Δt

0.09

Time step size history - Δt (s)

Displacements (m)

0.06 0.04 0.02 0.00 -0.02 -0.04

Constant Δt

-0.06

Ansys Bergan & Mollestad Hulbert & Jang Lages et al.

Ansys Bergan & Mollestad Hulbert & Jang Lages et al.

0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01

-0.08

0.00 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0

0.5

1

1.5

Time (s)

2

2.5

Time (s)

(a)

(b) Fig. 18. Results of the beam analysis for the first case of loading.

3

3.5

4

4.5

5

130

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

Curvatures for

7.00

Curvatures for

7.00

Non-Regularized Curvature

Regularized Curvature

6.00

5.00

Curvature - κ (1/m)

Curvature - κ (1/m)

Non-Regularized Curvature

Regularized Curvature

6.00

4.00

3.00

5.00

4.00

3.00

2.00

2.00

1.00

1.00

0.00

0.00 0

0.5

1

1.5

2.5

2

3

3.5

4

4.5

5

0

0.5

1

1.5

2

Time (s)

2.5

3

3.5

4

4.5

5

Time (s)

(a)

(b) Fig. 19. Variation in the curvature for different values of a.

Bergan and Mollestad

Lages et al. 0.10

0.10 α = 0.9

Time step size history - Δt (s)

Time step size history - Δt (s)

Δt1 = 0.001

0.09

α = 1.0

0.08

0.06

0.04

0.02

Δt1 = 0.1

0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01

0.00

0.00 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0

0.5

1

1.5

2

Time (s)

2.5

3

3.5

4

4.5

5

Time (s)

(a)

(b) Fig. 20. Time step size history for different values of the input data.

Hulbert and Jang

Ansys 11.0 0.14

0.05

Δt1 = 0.001

0.13

Δt1 = 0.05

0.12

Time step size history - Δt (s)

Time step size history - Δt (s)

0.06

0.04

0.03

0.02

0.01

Δt1 = 0.001 Δt1 = 0.01 Δt1 = 0.1

0.11 0.10 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01

0.00

0.00

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0

0.5

1

1.5

Time (s)

(a)

2

2.5

Time (s)

(b) Fig. 21. Time step size history for different values of the input data.

3

3.5

4

4.5

5

131

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136 0.049

0.035

Displacements (m)

0.028

Constant Δt

0.040

Ansys Bergan & Mollestad Hulbert & Jang Lages et al.

0.035

Constant Δt

Time step size history - Δt (s)

0.042

0.021 0.014 0.007 0.000 -0.007 -0.014

Ansys Bergan & Mollestad Hulbert & Jang Lages et al.

0.030 0.025 0.020 0.015 0.010

-0.021 0.005

-0.028 0.000

-0.035 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0

5

0.5

1

1.5

2

2.5

Time (s)

Time (s)

(a)

(b)

3

3.5

4

4.5

5

Fig. 22. Results of the beam analysis for the second loading case.

illustrates that for this strategy, the time step increases rapidly from the initial low value of 0.001 s up to nearly 0.05 s and has a brief fall that is close to the time t = 0.5 s before a new increase to an unchanged value close to the ideal value. In the strategy of Hulbert and Jang, the target value given by (Dt/T)target is the same as the method of Bergan and Mollestad, which was almost reached after some intermediary plateaus and remained stable. The time step history obtained with ANSYS 11.0 is close the history given by the method of Bergan and Mollestad. The same behaviour shown in the previous example with the strategy of Lages et al. in quiescent start conditions can be observed in this case: the size of the second time step (Dt2) is close to the maximum value Dtmax = 0.1 s due to an initial curvature that is close to zero; however, the final sizes of the time increment oscillated in values that were much smaller than those obtained with the other strategies. Fig. 19(a) illustrates the variation in the curvature for this analysis, with and without regularisation, illustrating that the result of jreg is very oscillatory between local peak values. This behaviour leads to an undesirable oscillation in the time step sizes, which is eliminated by taking a = 1.0 in the adaptive procedure, leading to the curvature history of Fig. 19(b). Fig. 20(a) presents a comparison between the time step size histories obtained with this adaptive strategy while choosing the different values a = 0.9 and a = 1.0. In the second case, the time

step is less oscillatory and decreases in constant plateaus. In general, the value a = 1.0 should be chosen in problems where the peak values of the curvature grow over time, such as in the last example. An alternative choice to this procedure would be to take a value of ct that makes the regularisation interval Dtreg close to the oscillation period of the curvature, which is difficult to predict. The same analysis of this loading case was conducted with the strategies of Bergan and Mollestad and Hulbert and Jang, taking different values for the initial value Dt1 equal to 0.1 s and 0.05 s, respectively. The comparison results for variations in the time step are shown in Figs. 20(b) and 21(a). In the strategy of Bergan and Mollestad, the time step decreases rapidly from the high initial value to the same value as in Fig. 18(b), assuming a similar behaviour. In the strategy of Hulbert and Jang, the initial step size decreased after the first step due to a quiescent start, and the time steps reached the target value after some intermediary plateaus, remaining stable. This result proves that both strategies are not dependent on the choice of the initial time step. The variation in the time step was also analysed in the program ANSYS 11.0, with different initial values Dt1 = 0.01 s and Dt1 = 0.1 s, for which the results are shown in Fig. 21(b). However, in this case, the time step behaviour changed with the choice of the initial step size, and the software made no changes in the value Dt1 = 0.1 s.

Curvatures for

90.00

Curvatures for

90.00

Non-Regularized Curvature

Non-Regularized Curvature

80.00

80.00

Regularized Curvature

70.00

Curvature - κ (1/m)

70.00

Curvature - κ (1/m)

Regularized Curvature

60.00 50.00 40.00 30.00

60.00 50.00 40.00 30.00

20.00

20.00

10.00

10.00

0.00

0.00 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0

0.5

1

1.5

2

2.5

Time (s)

Time (s)

(a)

(b) Fig. 23. Variations in the curvature for different values of a.

3

3.5

4

4.5

5

132

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

Fig. 24. Plane frame steel structure.

Fig. 25. Loading and models of the plane frame structure.

The results for the third mode forcing frequency case are shown in Fig. 22. Fig. 22(a) shows the response history at the midpoint of the beam’s length. The transient part of the vibration is gradually damped out and the particular part of the solution dominates the response. The characteristic period calculated by Bergan and Mollestad would become nearly equal to the period of the loading function, which is 1/9 = 0.111. . .s. In both strategies of Bergan and Mollestad and Hulbert and Jang the ideal time step size would assume nearly 5% of this period, what is successfully obtained in the second strategy, as shown in Fig. 22(b). The result using algorithm of Bergan and Mollestad is very oscillatory, but it remained close to the ideal plateau in most part of the analysis. The time steps obtained with ANSYS 11.0 stayed nearly close to 0.01 s, but assumed so high values in some moments, that may have caused a loss of precision in the response seen in Fig. 22(a), compared with the others methods. In the strategy of Lages et al., the curvature increases and oscillates considerably, as shown in Fig. 23(a). The regularised curvature with a = 0.9 yields an oscillating time step that is close to a

nearly constant value of 0.0015 s. This behaviour can be fixed by taking a = 1.0, which leads to the variation shown in Fig. 23(b). In the same manner, the last example illustrates that the value a = 1.0 should be chosen when the average behaviour of the peak curvatures are increasing. As noted before, this procedure is useful in linear analyses, where an oscillating time step can increase the computational cost. 4.4. A three-bar plane frame structure The following example illustrates the application of the reviewed time adaptive techniques to a problem of a plane frame steel structure, with the geometric and physical properties shown in Fig. 24. The objective is to find the lateral displacement history Dx(t) of the top-right beam–column joint node, as illustrated in the last figure. Some cases of loading were considered. In the first case, a lateral transient load F1(t) is applied at the same target node, with the loading function illustrated in Fig. 25(a). To compare the results

Non-subdivided model

Subdivided model 0.022

0.011

Bergan & Mollestad

0.021

Hulbert & Jang Lages et al. Ansys

0.010

Time step size history - Δt (s)

Time step size history - Δt (s)

Bergan & Mollestad

0.004

0.003

0.002

0.001

0.000

Hulbert & Jang Lages et al. Ansys

0.020 0.004

0.003

0.002

0.001

0.000 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Time (s)

Time (s)

(a)

(b) Fig. 26. Time step size histories.

0.4

0.45

0.5

0.55

0.6

133

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

Non-subdivided model

4000

Subdivided model

400000

Non-Regularized Curvature

Non-Regularized Curvature

Regularized Curvature

Regularized Curvature

Curvature - κ (1/m)

Curvature - κ (1/m)

3000

2000

1000

300000

200000

100000

0

0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

0

0.05

0.1

0.15

0.2

Time (s)

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

Time (s)

(a)

(b) Fig. 27. Curvature history for both models of the frame structure.

3000 Non subdivided Model

Curvature - κ (1/m)

Subdivided Model

2000

Fig. 29. Loading and model of the plane frame structure. 1000

0 0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

Time (s) Fig. 28. Curvatures for the frame structures from time t = 0.05 s. Fig. 30. First multiple-loading case.

obtained with different finite-element beam meshes, the analysis is conducted on two different models, which are illustrated in Fig. 25(b) and (c). In the first model, there is no internal element division over the bars, which are described only by the initial and final nodes, yielding only four nodes and three bars in the entire structure. In the second model, the finite element mesh was refined in 24 equal beam elements of length equal to 0.5 m, yielding a total of 25 nodes. Because the support nodes of the structure are completely restrained, the model of Fig. 25(b) has six degrees of freedom and the model of Fig. 25(c) has 69 degrees of freedom. The Rayleigh proportional damping is introduced with damping ratios that are set to 5% of the critical amount in the first and third natural modes of vibration in both models. The first initial step in the dynamic analysis was chosen as Dt1 = 0.001 s. np and k were set equal to 1.4 and 0.05, respectively. (Dt/T)target was set equal to 0.1. The parameters ct and a were set equal to 1.0 and 0.7, respectively, and the value Dtmax in the algorithm of Lages et al. was chosen as 0.03 s, which is close to the fundamental period T1  0.0346 s. The parameter cp was set equal to 0.1 and 0.001 in the first and second models, respectively. The Newmark integration process was applied with c = 1/2 and b = 1/4, and in the method of Lages et al., the expression Dtreg = ctDtmax was used again because there is no Dtcr in this case. The step size history results obtained for the divided and nondivided models are shown in Fig. 26. The responses for the lateral

Fig. 31. Second multiple-loading case.

displacement Dx appear to be identical in all of the analyses and are not shown. The dominant period T calculated by Bergan and Mollestad stayed close to the fundamental period T1, yielding an ideal time step size for this strategy that is equal to 1/20 of T1, which is close to 0.00175 s. The target time step of Hulbert and Jang was close to twice that value because (Dt/T)target is twice k. ANSYS 11.0 yielded time step size values that are between the two plateaus. When applying the strategy of Lages et al., a characteristic situation of multi-degree-of-freedom systems has occurred, as cited in Section 3.3.1: the curvature takes such high values that the time step would be very small if the parameter cp did not reduce the value of j. The curvature value at the initial time step increases with an increasing number of degrees of freedom, which

134

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

justifies the choice of cp = 0.1 in the non-divided bars case and cp = 0.001 in the divided case because the magnitude of the initial curvature in the second model is approximately 100 times greater than in the first model, as shown in Fig. 27. However, this fact occurs only in the initial time steps because in the next instants, the curvatures in both models are approximately the same magnitude, as shown in the interval plotted in Fig. 28, which could explain why the time step in the method of Lages et al. starts with close to the same value in both models and reaches much higher values at the end of the second model, which results from the reduced cp. However, in the first model, the time step size ends with values that are much greater than in the other strategies (Fig. 26(a)) because the initial curvature already has a value that is much greater than in the subsequent instants, which leads to a small initial assumed value for cp in the final step calculations. This phenomenon occurs because when the initial velocity vec_ 0 ) is equal to 0, the curvature at time t = 0 is given by the tor (D € 0 k), which will have a greater norm of the acceleration vector (kD value when there are more degrees of freedom in the structure. Even if this behaviour does not occur in the subsequent moments, the curvature regularisation algorithm continues to use a high percentage of the values that were originally calculated. This phenomenon does not occur in quiescent start condition problems because the initial acceleration is zero in this case. To explore the results of the automatic time-stepping procedures for multiple transient loadings, the same plane frame steel structure was analysed with two acting loads, Fa(t) and Fb(t), which are applied on the nodes shown in Fig. 29. In this case, only the subdivided model was considered. The proportional damping is introduced with a damping ratio of 10% for the first and third modes. Two groups of loading were

considered, which are given by F2 and F3, as illustrated in Figs. 30 and 31, respectively. The variations in the time steps for both loading combinations are shown in Fig. 32. The same parameter values applied in the previous example were used, changing only the values of a and Dtmax from the strategy of Lages et al., to 0.8 and 0.05 s, respectively. The previous values of cp, 0.1 and 0.001 were used on the cases of F2 and F3, respectively. In the same way as in the last example, the time step sizes calculated by Bergan and Mollestad remained close to 1/20 of the fundamental period in both cases of loading. ANSYS yielded results with similar behaviour, but with values that were slightly greater, staying closer in the second case. From using the strategy of Hulbert and Jang, the time step was rather oscillatory in the final instants in both cases because the normalised error RL stays lower than and greater than the tolerance limits, sequentially. From using the algorithm of Lages et al., the time step was increased due to the continuous curvature reduction, staying in the same range of values that was obtained with the other strategies. The results were almost the same in both loading cases, even with very different values of cp, because the first case has quiescent start conditions, which lead to zero initial curvature, and the second case is similar to the previous example, with a very large initial curvature.

5. Computational cost Similarly to Table 2, the following Table 5 summarizes the results of analyzes of the previous multi-degree-of-freedom systems, in an attempt to compare the computational cost of the algorithms. The maximum and minimum time increments found with

First loading case

Second loading case

0.0035

0.0030

Bergan & Mollestad

Time step size history - Δt (s)

Bergan & Mollestad

Time step size history - Δt (s)

Hulbert & Jang Lages et al. Ansys

0.0030 0.0025 0.0020 0.0015 0.0010

Hulbert & Jang Lages et al. Ansys

0.0025

0.0020

0.0015

0.0010

0.0005

0.0005

0.0000

0.0000 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

0

0.05

0.1

0.15

0.2

0.25

Time (s)

0.3

0.35

0.4

0.45

0.5

0.55

0.6

Time (s)

(a)

(b)

Fig. 32. Time step size history results for the two cases of applied loads in the frame structure.

Table 5 Summary of results of multi-degree-of-freedom systems. Example

Shear build. (1) Shear build. (2) Beam (1) Beam (2) Frame, load F1 Frame, load F2 Frame, load F3

Bergan and Mollestad

Hulbert and Jang

Lages et al.

Dt multip. factors

Dtmin

Dt max

NTOT

Dtmin

Dt max

NTOT

Dtmin

Dtmax

NTOT

5.4994 1.0000 1.0000 1.0000 1.9931 2.8048 0.3907

18.3063 21.8786 48.5011 19.8014 16.2555 17.2158 17.2371

103 131 120 801 376 403 451

4.4521 0.3445 0.3546 0.3549 0.0017 1.8334 0.0009

18.4874 14.0845 44.8030 5.5313 33.3489 30.7454 13.3440

93 258 227 1035 284 948 1712

1.5787 1.0000 1.0000 1.0000 0.9327 2.9520 3.0195

18.2792 33.3247 98.0206 35.5056 211.381 26.6504 19.2714

247 306 249 3234 1093 812 831

103 seg.

104 seg.

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

adaptive strategies are presented, as well as the total number of points (NTOT) plotted with the different algorithms. The higher the number NTOT, the more expensive is the strategy and the more detailed is the dynamic response. For the simply supported beam systems, only results of the analysis with the same data input for all strategies are shown, and for the three-bar plane frame systems, the results obtained only with the subdivided models are presented. The last column of Table 5 brings the multiplication factors of the values of Dt presented in the table. Observing the summary results, it can be noted that the strategy of Bergan and Mollestad lead to the smallest values of NTOT in almost all analyses (with higher values of Dt), except in two cases, where the strategy of Hulbert and Jang reached smaller values. The strategy of Lages et al., on the other hand, gave bigger values of NTOT in many of the numerical examples, or a value near to the higher one, except for the last case, where the strategy of Hulbert and Jang gave a value of NTOT greater than twice the value obtained with Lages et al. By only this comparison, one would conclude that the strategy of Bergan and Mollestad leads to a less expensive analysis, and that an analysis using the strategy of Lages et al. has a higher computational cost. However, the cost of operation in each time step size selection should be also considered (which should be small, as it was noted in the statement 5 in the beginning of Section 3), as well as the accuracy obtained with the numerical analyses. The computational cost of the strategies in the time step size selection can be compared by counting the number of numeric operations performed by each algorithm. A common technique for estimating the efficiency of an algorithm is to analyse it according to the size of the problem, given by the number of processed elements ‘‘n’’. By calculating the number of operations performed on the n elements, a complexity function f(n) can be found for the algorithm, and the asymptotic behaviour of the efficiency can be estimate by taking big values of n. For a structural analysis problem, n would be the number of degrees of freedom of the system. To compare the three automatic time-stepping strategies with this approach of ‘‘complexity function’’, it should be considered only the routine of time step size selection in each problem, because the numerical integration process is the same for all algorithms. The most expensive operation executed in the routine of Bergan and Mollestad is a matrix–vector multiplication that appears on the terms of Eq. (12), which has complexity of n2. The second most expensive operation of this strategy is the norm-vector calculation of Eq. (15) (which has the cost of n) followed by minor fixed simple arithmetic operations in each step. Thus, it is found that the strategy of Bergan and Mollestad has complexity of n2 , because its complexity function has the form f(n) = An2 + Bn + C. Both strategies of Hulbert and Jang and Lages et al., in the other hand, has complexity of n, since their complexity function has the form f(n) = Bn + C, because the most expensive operations in the routines of these algorithms are the norm-vector calculations of Eqs. (19) and (29), followed by minor fixed operations. It is concluded that the strategy of Bergan and Mollestad is the most expensive when only the additional step size selection computational cost is considered. However, it should be noted that the expensive operation of factorising and solving the matrix Eq. (8) is always performed in the numerical analysis, because it is part of the numerical integration process. This operation may have a computational cost of complexity greater than n2 in some cases, which leads the whole algorithm analysis to this order of asymptotic complexity anyway. Moreover, the gain in not make unnecessary changes the value of Dt could be better than decrease the cost of the individual time step size selection, because it avoid a new factorization of the effective stiffness matrix in Eq. (9) to solve the system.

135

The global analysis that would be made to compare the computational cost of the automatic time-stepping strategies, and to give a final conclusion about the less expensive technique, must consider these three criteria together: the total number of steps taken (NTOT, or the number of plotted points), the cost of individual time step size selection routine (which could be given not only by the asymptotic behaviour of efficiency, but by the explicit complexity function f(n), to be more specific) and the frequency of chances in the value of Dt (which leads to changes in the matrix K). The global comparison depends on these aspects, which are affected by various parameters and become very complex. The accuracy obtained with the numerical analyses by different automatic time-stepping strategies should also be considered for completeness, what could be made only by extensive simulation of various systems and comparison with analytical or more accurate responses, in the same way it was done in some numerical examples in the previous section. The strategy of Bergan and Mollestad is very similar to that used by the software ANSYS, which has been used and which reached accurate results in practical analyses. The strategy of Hulbert and Jang has itself an error estimate that intends to control the accuracy results. The strategy of Lages et al. should be more applied in future simulations to evaluate accuracy, because it was developed recently. 6. Conclusions Observing the formulation and numerical example results, one can draw some conclusions about the automatic adaptive timestepping procedures. The strategy of Bergan and Mollestad [6], which is based on the current characteristic frequency, can easily be implemented in any numerical integration procedure, facilitating programming. This approach is the oldest method among the three studied approaches and has been successfully applied in commercial finite element software. However, it was observed that the oscillations of the time increments were too large in some examples when using this strategy (Fig. 22b) despite remaining close to a constant plateau  is close to zero value. This pattern occurs because the frequency x when the response passes through maximum amplitudes, as the norm ||DD|| becomes relatively small in these situations (a fact that has been observed in other research). Even if verifying the  i , the technique condition (15) for discarding some frequency x was not very effective, even when using values for e that were different from 0.1. In some cases, this fact leads to a very frequent oscillation in the time step sizes. This method also yields no additional gain in the analysis of single-degree-of-freedom systems, which have only one characteristic period. The suggested values of the parameters k,  np ;  nm , and e were used, but how to define these values remains unclear. The strategy of Hulbert and Jang [12] is very intuitive when using integration errors, which is a concept that has been suggested in various other studies (e.g., [7,10]) and appears to be a more suitable technique for future research, perhaps using more recent error estimators [16]. The algorithm was very effective in the solution of single-degree-of-freedom systems and also led to time step sizes that were more regular when solving multidegree-of-freedom systems. However, the last example of a frame structure in Section 4.4 displayed values of Dt that were very oscillatory and that were calculated by this strategy. Moreover, the formulation of the local error presented by Hulbert and Jang is attached to the integration algorithm of the generalised-a method and must be adapted to any other numerical procedure, which is not necessary in the other strategies. Even with the intuitive use of Eq. (21) for tolerance intervals, the time step changing Eqs. (22) and (23) are not very intuitive and should be studied in more depth. Again, the suggested values of the parameters

136

D.F. Rossi et al. / Engineering Structures 80 (2014) 118–136

(Dt/T)target, lb, pinc, and pdec were used, but how they are defined remains rather unclear. The strategy of Lages et al. [17] is the latest strategy that was developed, and it has an interesting and intuitive derivation, which uses the geometric indicator of the displacement history curvature. The approach is not attached to the specific choice for the numerical integration procedure and can be applied to any of them. However, the implementation had some difficulties in the relation between the curvature indicator and time step size. Despite the regularisation procedure used in processing the curvature, the curvature is not stabilised to any constant value, which causes a continuous change in the time increment even when the system has a stationary response, which violates the second criterion shown in Section 3. Thus, the strategy needs a better development on an application to structural vibration, particularly to solve the problems observed when analysing the quiescent start condition systems (which could have initial curvatures that are close to zero, which leads to large time step sizes) and systems with many degrees of freedom (which could have greater initial curvatures, leading to unworkable time step sizes that are close to zero). The computational cost analysis of the three strategies should consider various aspects when choosing the one with better performance. For accuracy optimisation only, for example, it could be adopted the procedure of calculate the time step size with all strategies and take the smaller value of Dt to be used in the step solution. Acknowledgments The authors would like to thank the institutions CNPq, CAPES, and FAPES, which provided support for this work. References [1] Hibbitt HD, Karlsson BI. Analysis of Pipe Whip. Electric Power Research Institute – rept. no. EPRI NP-1208, Palo Alto; November 1979. [2] Oughourlian CV, Powell GH. ANSR-III: general purpose computer program for nonlinear structural analysis. Earthquake Engineering Research Center – report no. UCB/EERC-82/21, Berkeley; November 1982.

[3] Felippa CA, Park KC. Direct time integration methods in nonlinear structural dynamics. Comput Methods Appl Mech Eng 1979;17–18(2):277–313. [4] Park KC, Underwood PG. A variable-step central difference method for structural dynamics analysis – part 1. Theoretical aspects. Comput Methods Appl Mech Eng 1980;22(2):241–58. [5] Underwood PG, Park KC. A variable-step central difference method for structural dynamics analysis – part 2. Implementation and performance evaluation. Comput Methods Appl Mech Eng 1980;23(3):259–79. [6] Bergan PG, Mollestad E. An automatic time-stepping algorithm for dynamic problems. Comput Methods Appl Mech Eng 1985;49(3):299–318. [7] Zienkiewicz OC et al. A unified set of single step algorithms. Part 1: general formulation and applications. Int J Numer Methods Eng 1984;20(8):1529–52. [8] Zienkiewicz OC, Xie YM. A simple error estimator and adaptive time stepping procedure for dynamic analysis. Earthq Eng Struct Dynam 1991;20(9):871–87. [9] Zeng LF et al. A posteriori local error estimation and adaptive time-stepping for Newmark integration in dynamic analysis. Earthq Eng Struct Dynam 1992;21(7):555–71. [10] Li XD, Zeng LF, Wiberg N-E. A simple local error estimator and an adaptive time-stepping procedure for direct integration method in dynamic analysis. Commun Numer Methods Eng 1993;9(4):273–92. [11] Newmark NM. A method of computation for structural dynamics. ASCE J Eng Mech Div 1959;85(3):67–94. [12] Hulbert GM, Jang I. Automatic time step control algorithms for structural dynamics. Comput Methods Appl Mech Eng 1995;126(1–2):155–78. [13] Chung J, Hulbert GM. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-a method. J Appl Mech 1993;60(2):371–5. [14] Hilber HM, Hughes TJR, Taylor RL. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthq Eng Struct Dynam 1977;5(3):283–92. [15] Wood WL, Bossak M, Zienkiewicz OC. An alpha modification of Newmark’s method. Int J Numer Methods Eng 1980;15(10):1562–6. [16] Chung J, Cho E-H, Choi KA. A priori error estimator of the generalized-a method for structural dynamics. Int J Numer Methods Eng 2003;57(4):537–54. [17] Lages EN et al. An adaptive time integration strategy based on displacement history curvature. Int J Numer Methods Eng 2013;93(12):1235–54. [18] Clough RW, Penzien J. Dynamics of structures. 3rd ed. Berkeley: Computers & Structures Inc.; 1995. [19] Chopra AK. Dynamics of structures: theory and applications to earthquake engineering. Englewood Cliffs: Prentice Hall; 1995. [20] Cook RD, Malkus DS, Plesha ME. Concepts and applications of finite element analysis. 3rd ed. Madison: John Wiley & Sons; 1989. [21] Hughes TJR. A note on the stability of Newmark’s algorithm in nonlinear structural dynamics. Int J Numer Methods Eng 1977;11(2):383–6. [22] Hughes TJR. The finite element method: linear static and dynamic finite element analysis. Englewood Cliffs: Prentice Hall; 1987. [23] ANSYS, INC. Theory reference for ANSYS and ANSYS workbench; 2007.