An optimized three-sub-step composite time integration method with controllable numerical dissipation

Computers and Structures 231 (2020) 106210

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

An optimized three-sub-step composite time integration method with controllable numerical dissipation Yi Ji a,b, Yufeng Xing a,⇑ a b

Institute of Solid Mechanics, Beihang University, Beijing 100083, China Shen Yuan Honors College, Beihang University, Beijing 100083, China

a r t i c l e

i n f o

Article history: Received 4 June 2019 Accepted 16 January 2020

Keywords: Optimization Composite time integration method Stability and accuracy Dissipation and dispersion

a b s t r a c t This paper proposes an optimized three-sub-step composite time integration methods with controllable numerical dissipation, called the q1-Optimal-Trapezoidal-Trapezoidal-Backward-Interpolation-Formula (q1-OTTBIF) method. In this method, a novel Newmark-like method or four-point backward interpolation formula is employed in the third sub-step, instead of the four-point Euler backward difference method as in the Optimal-Trapezoidal-Trapezoidal-Backward-Difference-Formula (OTTBDF) method which was proposed by the present authors and co-workers. The proposed method has second-order accuracy, unconditional stability and controllable numerical dissipation, and the spectral radius q1 serves as a parameter controlling the degree of numerical dissipation. In addition, the properties of the q1-OTTBIF method can reduce to those of the OTTBDF method when q1=0. Since the low-frequency accuracy is maximized in construction, the proposed method has higher low-frequency accuracy than other methods with controllable numerical dissipation. Linear and nonlinear numerical simulations are conducted to check the advantages of the proposed method over other similar time integration methods. Ó 2020 Elsevier Ltd. All rights reserved.

1. Introduction Many advanced time integration methods have been developed since 1960s, and most of them were motivated by the pursuit of high accuracy, efficiency, stability and controllable numerical dissipation. Motivating factors mainly include following facts: (1) The dissipative Newmark algorithms are only first-order accurate; (2) Algorithms based on series expansion method are conditionally stable or unstable when algorithmic orders are more than two [1]; (3) Unconditionally stable second-order Newmark algorithms, such as the Trapezoidal Rule (TR), have no numerical dissipation and may be failure for simple nonlinear problems [2]; (4) The accuracy, efficiency, stability and numerical dissipation should be improved concurrently. To solve the first problem, several classical and second-order amethods, such as the Generalized-a method, the HHT-a method and the WBZ-a method, were developed, and they all have control-

⇑ Corresponding author. E-mail address: [email protected] (Y. Xing). https://doi.org/10.1016/j.compstruc.2020.106210 0045-7949/Ó 2020 Elsevier Ltd. All rights reserved.

lable numerical dissipation. Besides, the Three-Parameter SingleStep Method (TPSM) [3] was constructed recently by the present authors and co-workers, in which the equations of motion are exactly satisfied at time points, different from the classical amethods. Because of the fact 2, developing a high-order and unconditionally stable time integration method is a challenging work. To reach such a goal, instead of using series expansion approach, Fung developed a series of unconditionally stable higher-order methods based on the differential quadrature method [4] and the weighted residual method [5], and Kim et al. [6] constructed a new highorder method which can eliminate two additional procedures required in the high-order Fung’s methods [4]. Besides, a few conditionally stable high-order methods with large stability limits were also proposed in recent decade. Rezaiee-Pajand et al. proposed an Implicit Higher-Order Accuracy (IHOA) method [7] which stability limit was improved [8]. Xing et al. [9] proposed a time finite element method based on the differential quadrature rule and the Hamilton’s variational principle (DQTFEM). Afterwards, the improved DQTEM [10] was presented where both spatial and temporal domains were discretized via the differential quadrature rule, reducing the dimension of coefficients matrix without loss of accuracy. The fact 3 says that the trapezoidal rule, which is an unconditionally stable second-order method, may be failure for simple

2

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

nonlinear problems [2], so for the sake of conserving energy of nonlinear systems, many conserving energy methods, such as the Constraint Energy Method (CEM) [11], the Constraint Energy Momentum Algorithm (CEMA) [12], the Energy-Momentum Method (EMM) [13–14] and the Gonzalez method [15], have been proposed. In recent years, based on the EMM method, a weak form quadrature element formulation [16] of geometrically exact beam model was proposed for dynamic analysis of spatial beams undergoing large displacement and rotation, and a general momentumbased time integration scheme [17] with energy conservation were developed. For improving algorithmic performances comprehensively, as the fact 4 says, combining different type of methods in a scheme is a powerful way, which takes accuracy, efficiency and numerical dissipation together into account. In 2005, Bathe et al. [18] employed a composite method [19] to solve second-order linear equations. The Bathe method is a two-sub-step method where the TR and the Backward-Difference-Formula (BDF) are used in the first and second sub-steps respectively, having second-order accuracy, unconditional stability and strong numerical dissipation. After that, Bathe et al. investigated the conserving properties of the method for nonlinear dynamics in 2007 [20], gave insight into this method for linear dynamics in 2012 [21] and examined its performances in wave propagation problems in 2013 [22]. Based on the works of Bathe and his co-workers, a few novel composite time integration methods have been proposed in past decade, such as the three-sub-step Trapezoidal-Trapezoidal-Backward-DifferenceFormula (TTBDF) method [23], the three-sub-step Trapezoidal-Ba ckward-Difference-Formula-Houbolt (TBDFH) method [24], the optimized three-sub-step OTTBDF method [25] and the optimized four-sub-step Optimal-Trapezoidal-Trapezoidal-Trapezodial-Back ward-Difference-Formula (OTTTBDF) method [25] etc. Most of these composite methods combine TR and BDF. The present authors [26] have discussed the effects of designable parameters, including the number of sub-steps, the number of time difference points used in BDF and the methods used in sub-steps, on the properties of composite methods, and provided a few rules for the construction of this type of composite methods. In the composite methods mentioned above, the spectral radius q1 tends to zero in high-frequency range, in another word, these methods are L-stable or asymptotic annihilating. To improve numerical dissipation properties, some composite methods with controllable numerical dissipation were proposed. In 2019, Malakiyeh et al. [27] proposed the two-sub-step b1/b2-Bathe method, where q1 was smoothly controlled by changing parameters b1/b2, but there is no explicit relationship among them. Note that in the b1/b2-Bathe method, a so-called three-point trapezoidal rule with the parameters b1 and b2 was used in the second substep, instead of the three-point BDF. Afterwards, based on the b1/b2-Bathe method, Noh et al. [28] presented the two-sub-step q1-Bathe method, where the only one algorithmic parameter is explicitly expressed in terms of q1. In addition, another type of composite methods based on collocation method have also been proposed in recent years. Huang et al. [29] proposed a composite scheme [30] where the collocation method with third-order accuracy (CM3) and the collocation method with fourth-order accuracy (CM4) were used alternately. Kim et al. [31] proposed a composite scheme based on the collocation time finite element approach, and it was improved by using the time weighted residual finite element method [32]. Besides, a combination of TR and the higher-order Newton backward extrapolation functions generated a multi-sub-step higherorder implicit composite time integration method [33], but the method of more than four sub-steps is unstable and the method with more than two-steps is conditionally stable. For clarity, Table 1 lists five types of time integration methods reviewed

above, and representative methods for each type as well as their properties. As it can be seen from above review that the existing composite methods with controllable dissipation are not developed based on optimization methods, which means composite methods with controllable dissipation and better numerical performances are achieveable. In this context, this work proposes an optimized three-substep composite time integration methods with controllable numerical dissipation, call the q1-Optimal-Trapezoidal-Trapezoidal-Back ward-Interpolation-Formula (q1-OTTBIF) method, which is motivated by the two-sub-step q1-Bathe method. In the q1-OTTBIF method, a novel Newmark-like method is employed in the last sub-step, and second-order accuracy, unconditional stability and controllable numerical dissipation are reached through the optimization of algorithmic parameters. This paper is organized as follows. In Section 2, the governing formulations of the q1-OTTBIF method is presented. Then the properties of the q1-OTTBIF method, such as efficiency, memory, stability, accuracy and dissipation, are discussed in Section 3, and they are compared with those of the TPSMi (‘‘i” represents implicit) method [3], the q1-Bathe method [28] and the collocation time integration algorithm [32], which is referred to as the Kim method in this paper in Section 4. Numerical experiments are conducted in Section 5. And the conclusions are drawn in Section 6. 2. Formulations In this section, the formulations of the q1-OTTBIF method are provided, including fundamental equations, and relationships among algorithmic parameters. Assuming that the time step size is h and all variables up to time t are known, the q1-OTTBIF method is used to calculate the solutions at time t + h. The time step [t, t + h] is divided into three sub-steps as [t, t + c1h], [t + c1h, t + c2h] and [t + c2h, t + h], where 0 < c1 < c2 < 1. In the first two sub-steps, the trapezoidal rule is used as in the OTTBDF method [25] as follows


 xtþc1 h ¼ xt þ c12h x_ t þ x_ tþc1 h
 €tþc1 h €t þ x x_ tþc1 h ¼ x_ t þ c12h x

ð1Þ


 xtþc2 h ¼ xtþc1 h þ ðc2 2c1 Þh x_ tþc1 h þ x_ tþc2 h
 €tþc2 h €tþc1 h þ x x_ tþc2 h ¼ x_ tþc1 h þ ðc2 2c1 Þh x

ð2Þ

Instead of using the four-point Euler backward difference formula as in the OTTBDF method, the following Newmark-like method that can be seen as an extension of the last sub-step of the q1-Bathe method is employed in the third sub-step,


 xtþh ¼ xt þ h h3 x_ tþh þ h2 x_ tþc2 h þ h1 x_ tþc1 h þ h0 x_ t
 €tþh þ h2 x €tþc2 h þ h1 x €tþc1 h þ h0 x €t x_ tþh ¼ x_ t þ h h3 x

ð3Þ

which in effect are the backward interpolation formula (BIF), that is why the present method is called the q1-OTTBIF method, rather than the q1-OTTBDF method. Considering a linear dynamic system, the equilibrium equations at time t + c1h, t + c2h and t + h are

€tþc1 h þ C x_ tþc1 h þ Kxtþc1 h ¼ Rtþc1 h Mx

ð4Þ

€tþc2 h þ C x_ tþc2 h þ Kxtþc2 h ¼ Rtþc2 h Mx

ð5Þ

€tþh þ C x_ tþh þ Kxtþh ¼ Rtþh Mx

ð6Þ

where M, C and K are the mass, damping and stiffness matrices respectively; R is the external load vector; x is the displacement € are respectively the first and second derivatives vector; x_ and x of x with respect to time t. Since the TR is second-order accurate, the method to be used in the third sub-step is also expected to

3

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210 Table 1 The properties of some representative time integration methods. Types

Method

Accuracy

Stability

Dissipation

Self-starting

Newmark family

TR CDM IDM (d = 11/20, a = 3/10)

Second-order Second-order First-order

Unconditional Conditional Unconditional

Non-dissipation Non-dissipation Uncontrollable

Yes No Yes

a-Methods

Generalized-a (General form) HHT-a (General form) WBZ-a (General form) TPSMi [3]

Second-order Second-order Second-order Second-order

Unconditional Unconditional Unconditional Unconditional

Controllable (exactly) Controllable Controllable Controllable (exactly)

Yes Yes Yes Yes

Higher-order methods

Fung [4] (n = 1,2,. . .) Kim & Reddy [6] (n = 1,2,. . .) HIOA [7] DQTFEM[9]

(2n-1)th-order (2n-1)th-order  First-order  Third-order

Unconditional Unconditional Conditional Conditional

Controllable (exactly) Controllable (exactly) Non-dissipation Non-dissipation

Yes Yes Yes Yes

Conserving energy methods

CEM [11] CEMA [12] EMM [13–14]

Second-order Second-order Second-order

Unconditional Unconditional Unconditional

Non-dissipation Uncontrollable Uncontrollable

Yes Yes Yes

Composite methods

Bathe [18] OTTBDF [25] q1-Bathe [28] Kim [32]

Second-order Second-order Second-order Second-order

Unconditional Unconditional Unconditional Unconditional

Uncontrollable Uncontrollable Controllable (exactly) Controllable (exactly)

Yes Yes Yes Yes

be second-order accurate to make the present q1-OTTBIF method second-order accurate. To reach this goal, the first relation in Eq. (3) is rewritten via the Taylor series expansion as
 xtþh ¼ xt þ h h3 x_ tþh þ h2 x_ tþc2 h þ h1 x_ tþc1 h þ h0 x_ t

_ R 3 ¼ Rtþh þ M

h3 þ h2 þ h1 þ h0 ¼ 1

ð8Þ

h3 þ c2 h2 þ c1 h1 ¼ 1=2

Using Eqs. (1)–(8), the time–stepping equations for the present method can be formulated as

_ _ K 1 xtþc1 h ¼ R 1

ð9Þ

_ _ K 2 xtþc2 h ¼ R 2

ð10Þ

_ _ K 3 xtþh ¼ R 3

ð11Þ

where all effective stiffness matrices and load vectors are given as follows

_ 4 2 K1 ¼ Mþ CþK c1 h c21 h2 _ R 1 ¼ Rtþc1 h þ M

4

4 €t xt þ x_ þ x c1 h t c2 h2

ð12Þ !

1

_ K2 ¼

4 2 2

ðc2  c1 Þ h

_ R 2 ¼ Rtþc2 h þ M

Mþ

2 CþK ðc2  c1 Þh 4



2 þC x þ x_ t c1 h t

 ð13Þ

ð14Þ

4 € x þ þx x_ 2 2 tþc1 h ð c  c1 Þh tþc1 h tþc1 h ðc2  c1 Þ h 2   2 þC xtþc1 h þ x_ tþc1 h ðc2  c1 Þh

!

ð15Þ _ 1 1 K3 ¼ 2 2M þ CþK h3 h h3 h

ð16Þ


 h2 x_ tþc2 h þ h1 x_ tþc1 h þ h0 þ h3 x_ t !

h23 h

€tþc2 h þ h1 x €tþc1 h þ h0 x €t h2 x h3   h2 x_ tþc2 h þ h1 x_ tþc1 h þ h0 x_ t xt þC þ h3 h h3

2 €t ¼ xt þ hðh3 þ h2 þ h1 þ h0 Þx_ t þ h ðh3 þ c2 h2 þ c1 h1 Þx

Therefore, the algorithmic parameters in third sub-step should satisfy

h23 h

þ 2

þ

€t Þ þ hh2 ðx_ t þ c2 hx €t Þ þ hh1 ðx_ t þ c1 hx €t Þ þ hh0 x_ t ¼ xt þ hh3 ðx_ t þ hx

ð7Þ

xt

ð17Þ

In addition, the relationship c2 / c1 = 2 is used in the q1-OTTBIF method to make the first two sub-steps share the same effective stiffness matrix for linear systems, and reduces free algorithmic parameters by one. Then together with the two relationships in Eq. (8), the number of free parameters changes from six to three. Assuming c1, h0 and h3 are free parameters, other parameters can be expressed as

4c1 ð1  h0  h3 Þ  1 þ 2h3 2c1 2c1 ðh0 þ h3  1Þ þ 1  2h3 h2 ¼ 2c1 h1 ¼

ð18Þ

For determining the free parameters c1, h0 and h3, an undamped single degree-of-freedom (SDOF) system is used below to further investigate the properties of the q1-OTTBIF method, and the equation of the SDOM system is

€ þ x2 x ¼ 0 x

ð19Þ

for which the compact recursive formulation has the form as



xkþ1 hx_ kþ1



  xk ¼A hx_ k

ð20Þ

where the expression of the amplification matrix A can be derived from Eqs. (1)–(3) and Eq. (19), as



A¼

A11

A12

s2 A12

A11

with

0



ð21Þ




1

s6 c31 h3 1 þ 3c1 þ 2h3  3c1 h3  4c1 h0 þ ! C B C B 4 c31 þ 16c21 h3  16c21 þ 16c21 h0 þ 16c1 h0 h3 B s c1 þC 2 2 C B 36 c h þ 24 c h þ 6 c  24h þ 12h A @ 3 1 3 1 3 1 3
2  2 2 8s c1 þ 2h3  1 þ 16 A11 ¼
2 
 c1 s2 þ 4 2 h23 s2 þ 1 ð22Þ

4

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

0

A12

4c21 h3  3c21 þ 4c21 h0  18c1 h3 þ

!

1

C B þ sc C B c1 þ 16c1 h23 þ 16c1 h0 h3 þ 6h3  12h23 C B ! C B 2 2 2 C B 6 c  4 c h  4 c h  3 c þ 6 c h  2h þ 1 1 3 3 1 1 3 1 0 A @ 4s2 þ 16 2 4h3 ¼
2 
 c1 s2 þ 4 2 h23 s2 þ 1 4 2 1

ð23Þ To exactly describe the amount of numerical dissipation in high frequency range, the relationship between q1 and three free parameters is established with the spectral radius analysis of A, which is

h3 ¼

4c1 h0  3c1 þ 1 q1 c1  3c1 þ 2

ð24Þ

where

q1 ¼ slim qðAÞ; q1 2 ½0; 1 !1

ð25Þ

As q1 is a parameter which needs to be given by users, Eq. (24) in fact is a relation of three parameters c1, h0 and h3. Then one needs to find another two relations to uniquely determine these three parameters, and these two relations can be found through examining the relationships among the percentage amplitude decay, period elongation and free parameters. If q(A)  1 is satisfied, the q1-OTTBIF method is unconditionally stable. The q(A) has the form as

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s6 f 1 þ s4 f 2 þ s2 f 3 þ f 4 qðAÞ ¼ s6 g 1 þ s4 g 2 þ s2 g 3 þ g 4

Fig. 1. Percentage amplitude decay versus h0 for the q1-OTTBIF method (s = 1, q1 = 0).

ð26Þ

where fi (q1, c1, h0) and gi (q1, c1, h0) (i = 1, 2, 3, 4) are the known functions of parameters q1, c1 and h0. Letting gi  0 and using q (A)  1 generate

s6 h1 þ s4 h2 þ s2 h3 þ h4 6 0

ð27Þ

where hi = fi – gi. Therefore letting hi  0, then q(A)  1 is satisfied for all s. According to the conditions gi  0 and hi  0, the range of h0 for c12(0,1/2) can be obtained as



h0 2

c1 4c3

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2 c1 c2 2ðq1 þ 1Þc3 þ ;  2ðq1 þ 1Þc3 þ c3 4c3 c3

ð28Þ Fig. 2. Percentage period elongation versus h0 for the q1-OTTBIF method (s = 1, q1 = 0).

where c1 ¼  2 þ 5c1  3c21  q1 c1 þ q1 c21


 c2 ¼ 2 þ 2c1  11c21 þ 3c31 þ 2q1 1  3c1 þ 3c21 þ c31 þ c21 q21 1  c1
 c3 ¼ 8 2  4c1 þ c21 þ q1 c21

ð29Þ where c1 is a free parameters. In order to ensure c3 to be positive, in mathematical sense, c1 should be in the range as,

c1 <

2

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ð1  q1 Þ 2 þ 2ð1  q1 Þ and c1 –0 or c1 > 1 þ q1 1 þ q1

ð30Þ

It can be easily verified that c12(0,1/2) falls in the left range, and adjusting c1 can greatly affect the property of numerical dissipation, refer to Section 3.3. Figs. 1–4 illustrate the percentage amplitude decay and period elongation versus h0 for several c1 for s = 1, and the results for s = 2 are shown in Figs. 5–8. It follows that the amplitude decay and the period elongation both take the minimum value when h0 is equal to the left endpoint value of Eq. (28). Then for maximizing low-frequency accuracy, the left endpoint value of h0 is used here, which is

h0 ¼

c1 4c3

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c2 2ðq1 þ 1Þc3 þ c3

ð31Þ

Fig. 3. Percentage amplitude decay versus h0 for the q1-OTTBIF method (s = 1, q1 = 0.5).

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

Fig. 4. Percentage period elongation versus h0 for the q1-OTTBIF method (s = 1, q1 = 0.5).

Fig. 5. Percentage amplitude decay versus h0 for the q1-OTTBIF method (s = 2, q1 = 0).

Fig. 6. Percentage period elongation versus h0 for the q1-OTTBIF method (s = 2, q1 = 0).

5

Fig. 7. Percentage amplitude decay versus h0 for the q1-OTTBIF method (s = 2, q1 = 0.5).

Fig. 8. Percentage period elongation versus h0 for the q1-OTTBIF method (s = 2, q1 = 0.5).

Fig. 9. Percentage amplitude decay versus c1 for the q1-OTTBIF method (q1 = 0).

6

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210 Table 3 Step-by-step solution procedure of the q1-OTTBIF method for linear systems. A. Initial calculations: 1. Form mass matrix M, stiffness matrix K and damping matrix C. €0. 2. Initialize x0, x_ 0 and x 3. Select time step h, the spectral radius q1 and c1, and calculate constants: c1 ¼  2 þ 5c1  3c21  q1 c1 þ q1 c21 ;


 c2 ¼ 2 þ 2c1  11c21 þ 3c31 þ 2q1 1  3c1 þ 3c21 þ c31 þ c21 q21 1  c1 ;
 c3 ¼ 8 2  4c1 þ c21 þ q1 c21 ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
 ffi 4c h 3c þ1 c1 2 q1 þ 1 c3 þ cc23 ; h3 ¼ q 1c 0 3c 1 þ2 ; h0 ¼ 4c 3 1 1

h1 ¼

ð4c1 1h0 h3 Þ1þ2h3 ; 2c1

h2 ¼

1

2c1 ðh0 þh3 1Þþ12h3 ; 2c1

b 3: b 1 and K 4. Form effective stiffness matrices K _ _ 4 2 1 1 K 1 ¼ 2 2 M þ c h C þ K; K 3 ¼ 2 2 M þ h h C þ K c1 h

h3 h

1

3

b 1 and K b 3: 5. Triangularize K b 1 ¼ L1 D1 LT ; K b 3 ¼ L3 D3 LT K 1 3 B. For each time step: Fig. 10. Percentage period elongation versus c1 for the q1-OTTBIF method (q1 = 0).

1. Calculate effective loads at time t + c1h:  

_ €t þ C c2h xt þ x_ t R 1 ¼ Rtþc1 h þ M 24 2 xt þ c4h x_ t þ x c1 h

Eq. (31) is a relation between c1 and h0, together with Eq. (24), there are two relations for three parameters. In order to uniquely determine them, one needs to determine any one of them first. Here, as an example, c1 is determined through the optimization of the percentage amplitude decay and period elongation. Figs. 9 and 10 illustrate the percentage amplitude decay and period elongation versus c1 for several s when q1 = 0. One can find the maximum point of the amplitude decay curves and the minimum point of the percentage period elongation curves are the same, as shown in Figs. 9 and 10, and they are

ð32Þ

The values of the obtained optimal c1 corresponding to some q1 are listed in Table 2. Since it is difficult to achieve the explicit relationship between c1 and q1 even if this relationship is definite and exact, therefore, Table 2 provides some optimal c1 for corresponding q1, and the optimal c1 for any other q1 can also be obtained through the optimization of the percentage amplitude decay and period elongation. As it can been seen from Table 2 that all c1 do not differ too much, so for other q1, the corresponding optimal c1 can be accurately achieved by linear interpolations. Therefore, it follows that the q1-OTTBIF method is an optimized composite time integration method since the low-frequency accuracy and the high-frequency dissipation are optimized in the determination of the three parameters. The step-by-step procedure of the q1-OTTBIF method is presented for linear systems in Table 3. In addition, it should be noted that the optimal c1 does not satisfy 4h23 = c21 which means that the third sub-step cannot has the same effective stiffness matrix as that of the first two sub-steps, resulting in extra computational cost for linear systems. Considering that the effective stiffness matrices are calculated only once in the entire solution process for linear systems, so the extra computational cost caused by not using an identical stiffness matrix for all sub steps is generally negligible, especially for long-term simulations.

1

1

1

3. Calculate velocities and accelerations at time t + c1h: 2ðxtþc1 h xt Þ 2 x_tþ h  €t þc1 h ¼ ð c1h x_t Þ  x €t  x_ t ; x x_ tþc1 h ¼ c h c1 1

1. Calculate effective loads at time t + c2h:  _ R 2 ¼ Rtþc2 h þ M

4 2 ðc2 c1 Þ2 h

€tþc h xtþc1 h þ ðc 4c Þh x_ tþc1 h þ x 1 2 1





þ C ðc 2c Þh xtþc1 h þ x_ tþc1 h 2 1

2. Solve for displacement at time t + c2h: b2 L1 D1 LT xtþc h ¼ R 1

2

3 .Calculate velocities and accelerations at time t + c2h: _ tþc h Þ 2ðxtþc2 h xtþc1 h Þ 2 x_ h x 1 €tþc h ¼ ð tþðcc2 €tþc h  x_ tþc1 h ; x x x_ tþc2 h ¼ ðc c Þh c Þh 2 1 2

c1 ¼ 0:3608506129

1

2. Solve for displacement at time t + c1h: b1 L1 D1 LT xtþc h ¼ R

1

2

1

1. Calculate effective loads at time t + h:   _ €tþc h þh1 x €tþc h þh0 x €t h2 x_ tþc2 h þh1 x_ tþc1 h þðh0 þh3 Þx_ t h2 x 2 1 R 3 ¼ Rtþh þ M 2xt 2 þ  h3 h23 h h3 h

h2 x_ tþc2 h þh1 x_ tþc1 h þh0 x_ t þC hxth þ h 3

3

2. Solve for displacement at time t + h: b3 L3 D3 LT xtþh ¼ R 3

3. Calculate velocities and accelerations at time t + h: € € € ðxtþh xt Þhðh2 x_ tþc2 h þh1 x_ tþc1 h þh0 x_ t Þ ðx_ x_ Þh h x þh x þh x €tþh ¼ tþh t ð 2 tþhhc2 h 1 tþc1 h 0 t Þ ;x x_ tþh ¼ hh3 3

3. Properties This section aims to examine the properties of the proposed method, including efficiency, memory, stability, accuracy and dissipation. 3.1. Efficiency and memory The performances of different methods should be compared in such a way: to compare the accuracy under the same computational cost, or to compare the cost under the same accuracy. Here the former is used. With the same step size h, the computational cost of m-sub-step method is roughly m/n times that of

Table 2 The relationships between the optimal c1 and q1 for the q1-OTTBIF method.

q1

0

0.1

0.2

0.3

0.4

0.5

c1 q1

0.3608506129 0.6

0.3572389164 0.7

0.3538916132 0.8

0.3507710317 0.9

0.3478472158 0.95

0.3450959228 1

c1

0.3424972371

0.3400345835

0.3376940094

0.3354636515

0.3343865666

0.3300100000

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

Fig. 11. Spectral radius versus s for the q1-OTTBIF method etc.

Fig. 12. Percentage amplitude decay versus s for the q1-OTTBIF method etc.

Fig. 13. Percentage period elongation versus s for the q1-OTTBIF method etc.

n-sub-step method. In other words, if the time step size of m-sub-step method is h and the time step size of n-sub-step method is nh/m, they should have the same computational cost, which means that the methods with more sub-steps require fewer

7

Fig. 14. Spectral radius versus s for the q1-(O)TTBIF method for several c1 (q1 = 0).

Fig. 15. Spectral radius versus s for the q1-(O)TTBIF method for several c1 (q1 = 0.5).

Fig. 16. Spectral radius versus s for different methods.

time steps of updating and less memory to store step state vari€tþh . The sub-step state variables within a step ables xtþh , x_ tþh and x €tþci h (i = 1, 2, 3, . . .), whose storages can be include xtþci h , x_ tþci h and x reused in next step if they are not necessarily stored.

8

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

Fig. 20. Harmonic motion of the support.

Fig. 17. Spectral radius versus ds for different methods.

Fig. 21. Displacement of node 1 for various methods.

Fig. 18. Percentage amplitude decay versus ds for these methods.

Fig. 22. Velocity of node 1 for various methods.

Fig. 19. Percentage period elongation versus ds for different methods.

3.2. Stability and accuracy Considering the SDOF system (19), Figs. 11–13 show the spectral radii, percentage amplitude decay and percentage period elon-

gation versus s respectively. Fig. 11 illustrates that the q1-OTTBIF method is unconditionally stable, and the numerical dissipation can be controlled exactly and smoothly by q1. Moreover, Figs. 12 and 13 show that the amplitude and phase properties of the q1-OTTBIF method with q1 = 0 can reduce to those of the OTTBDF method, and percentage amplitude decay and percentage period elongation decrease with the increase in q1.

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

9

different c1, Figs. 14 and 15 show the relationship between the spectral radius q and s, it can be seen that there are great decreases in the low-frequency accuracy when the optimal values of c1 are not taken. Moreover, when c1 is a positive number, the amount of dissipation tends to be constant with the increase in c1, and the amount of dissipation increases with the decrease in c1 when

Fig. 25. Displacement of node 1 in the z direction for various methods (Case 1).

Fig. 23. Acceleration of node 1 for various methods: (a) the interval [0, 20]; (b) the interval [8,18].

3.3. Dissipation As mentioned in Section 2, c1 can take any values of satisfying Eq. (30) in mathematical sense, not only within (0, 1/2). But if c1 in the q1-OTTBIF method is not optimal, as shown in Table 2, the present method is no longer an optimized method. For

Fig. 26. Velocity of node 1 in the z direction for various methods (Case 1).

Fig. 24. 3D pin-jointed truss dome.

10

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

can achieve better amplitude accuracy, but it is often accompanied with greater phase errors [28]. 5. Numerical examples Here linear and nonlinear examples are investigated to demonstrate the properties of the proposed method. The linear

Fig. 28. A clamped-free bar excited by end load.

Fig. 27. Acceleration of node 1 in the z direction for various methods (Case 1).

c1 is a negative number. So it follows that the q1-TTBIF method can provide appropriate numerical dissipation by using an appropriate c1, which is useful in the solution of some special problems, such as the wave propagation problems etc. Note that, for all c1 satisfying Eq. (30), the present method is always second-order accurate and unconditionally stable; for ensuring the low-frequency accuracy, c1 should be positive.

4. Comparisons In this section, the properties of the proposed method are compared with several existing time integration methods using parameter q1 to adjust the degree of numerical dissipation, including the TPSMi method [3], the q1-Bathe method [28] and the Kim method [32]. Fig. 16 shows the curves of the spectral radii versus s for different time integration methods. It follows that for all methods, the numerical dissipation in high-frequency range can be accurately controlled by q1. Fig. 17 shows the variations of the spectral radii with ds, where ds = s/n is used to ensure that all methods have the same computational cost (‘‘n” is the number of sub-steps). It can be clearly observed that the q1-OTTBIF method has more accurate low-frequency amplitude than other methods. To further assess accuracy, Figs. 18 and 19 show the percentage amplitude decay and the percentage period elongation against ds for different methods. It can be seen that under the same computational cost, the accuracy ranks: the q1-OTTBIF method > the q1-Bathe method  the Kim method > the TPSMi method, which means the proposed method has higher amplitude accuracy and phase accuracy than other methods, especially for a small q1. In addition, the case of q1 = 0.99 is compared for overall understanding the numerical property of all methods, and q1 = 0.99 is used here is because q1 cannot be equal to 1 in the Kim method. In these comparisons, c = 0.5 is used for the q1-Bathe method. If other values of c, such as 0.01, are used, the q1-Bathe method

Fig. 29. Displacement of the middle point.

Fig. 30. Velocity of the middle point.

Table 4 The average absolute errors of different methods in example 2. TPSMi

q1-Bathe

Kim

q1-OTTBIF

Case 1

Displacement Velocity Acceleration

2.4229e-3 4.2156e + 0 8.0006e + 3

2.4382e-3 4.2209e + 0 8.0480e + 3

2.4383e-3 4.2209e + 0 8.0480e + 3

2.4283e-3 4.2186e + 0 8.0117e + 3

Case 2

Displacement Velocity Acceleration

1.8568e-2 6.0394e + 0 1.1448e + 4

1.8691e-2 7.1443e + 0 1.1487e + 4

1.8691e-2 7.1443e + 0 1.1487e + 4

1.8689e-2 6.0890e + 0 1.1471e + 4

Error

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

Fig. 31. Section A-A of Fig. 30.

Fig. 32. Section B-B of Fig. 30.

11

12

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

Fig. 36. Mass-spring system.

Fig. 33. Displacement of the softening spring problem (0  t  2 T).

Fig. 37. Displacement of m1 for different methods (Case 1).

Fig. 34. Displacement of the softening spring problem (60 T  t  63 T).

Fig. 38. Velocity of m1 for different methods (Case 1).

Fig. 35. Displacement of the softening spring problem (102 T  t  105 T).

mass-spring system [21,25,28,29,32,34,36] and the threedimensional truss structure [7,25] are used to examine the dissipation and accuracy of the proposed method. The one-dimensional wave propagation problem [27,28,35,37] is used to check the performance in eliminating oscillations. The softening spring problem [31,38,39] and the nonlinear mass-spring system [3] are employed to examine the performance in solving nonlinear problems. It should be noted that the relationships of time step sizes of different methods satisfy h(TPSMi) = h(q1-Bathe)/2 = h(Kim)/2 = h(q1-OTTBIF)/3 to ensure that the computational costs of these methods are equal roughly.

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

Example 1: Harmonic motion of the support Considering a system as shown in Fig. 20, where the support undergoes a harmonic motion. The example represents a general structure with stiff and flexible parts due to the huge difference between k1 and k2. The high stiffness values in practice mainly provide constraints. Therefore, the example is used to examine the dissipation property of the proposed method or whether or not it can effectively filter out the contribution of high-frequency oscillations. Assuming k1 = 107, k2 = 1, m1 = m2 = 1 and the initial conditions are zero. The equation of motion of the system is



m1

0

0

m2



  €1 u k1 þ k2 þ € u2 k2

k2 k2



u1 u2



 ¼

k1 y



0

ð33Þ

where

E = 2.1  1011N/m2, q = 7860 kg/m3, A1-6 = 2300 mm2, A7-12 = 1250 mm2, A13-30 = 2125 mm2, respectively. The system is motivated by the initial vertical velocity of 50 m/s at three nodes 1, 3 and 5. For this problem, numerical dissipation is meaningless, so q1 = 0.99 is used in all methods. Two time steps are considered, which are h (TPSMi) = 0.0001s  Tmin/32 (denoted by Case 1) and h (TPSMi) = 0.0002s  Tmin/16 (denoted by Case 2), where Tmin = 2pxmax = 3.284810-3s. The reference solutions are achieved with the mode superposition method, and Figs. 25–27 show the numerical results of the tip node 1 in the z direction for Case 1. Note that in Figs. 25–27, there are no observable differences among these methods, so in the following, the average absolute ‘Error’ is used to distinguish their differences, and its definition is N P

y ¼ sin1:2t

ð34Þ

The dynamic responses of this system include the transient and steady state responses. The transient responses including lowfrequency component and high-frequency component can be damped out by numerical dissipation. The time step size of the q1-Bathe method is chosen to be T / 20  0.2618 where T = 2p / 1.2 is the period of steady response. And q1 = 0 is used for all methods. The reference solutions include the low-frequency and high-frequency steady state responses and are obtained by the mode superposition method. Figs. 21–23 compare the displacements, velocities and accelerations respectively. It follows that all methods can predict accurate displacements and velocities. And it can also be seen that all methods have acceleration ‘overshoot’; the TPSMi method also has velocity ‘overshoots’; the accelerations of the Kim method are less accurate than those of other methods. Fig. 23 provides the results of the interval [8,18] to observe the differences between these methods, where one can find that the q1-Bathe method and the Kim method show observable phase lags, especially the Kim method. As well known, large amount of dissipation is beneficial for relieving ‘overshoot’, so the acceleration results obtained by the q1-TTBIF method with c1 = 5 are presented in Fig. 23, showing that the q1-TTBIF method (c1 = 5) can effectively reduce the acceleration overshoot and provide satisfactory results for the rest of time steps. Example 2: The three-dimensional truss A three-dimensional pin-jointed dome, as shown in Fig. 24, is investigated here to show the accuracy of different methods. The structure has 19 nodes, 30 elements and 21 degree-of-freedoms. The modulus of elasticity, density and sectional areas are

jxNumerical  xReference j

Error ¼ k¼1

N

ð35Þ

where N is the number of time steps with [0, 0.3s]. It can be seen from Table 4 that all methods almost have the same accuracy for two cases, and the errors of the q1-Bathe method and the Kim method are slightly higher than those of the other two methods. This conclusion is consistent with the conclusion from the spectral analysis, refer to Fig. 19. Example 3: One-dimensional wave propagation in a clamped-free bar In this example, a clamped bar is excited by a constant external load F(t) = 104 at its end, as shown in Fig. 28, and the modulus of elasticity, density, sectional area and length are E = 3  107, q = 0.00073, A = 1 and L = 200 respectively. The bar is modelled with 1000 uniform linear rod elements. The time step size pffiffiffiffiffiffiffiffiffi h = CFL  l  q=E, where CFL = 1. To achieve high dissipation, the parameter q1 = 0 is used in the TPSMi method, the Kim method and the q1-(O)TTBIF method. For the q1-Bathe method, the parameter q1 = 0.7321 and c = 1.5774 are used [37]. In the q1-(O)TTBIF method, the parameter c1 = 1000 is used to arrive at appropriate dissipation, refer to Fig. 14. The displacement and velocity of the middle point (x = 100) are shown in Figs. 29–32. The results show that the displacement can be accurately predicted using all methods, but the velocity predictions contain some oscillations. It can be seen that the proposed method can completely eliminate the oscillations, and the q1-Bathe method can rapidly eliminate the oscillations. In addition, for the q1-TTBIF method and the q1-Bathe method, c1 (or c) can serve as a parameter to adjust the amount of numerical dissipation, beneficial to filtering out oscillations Example 4: Softening spring problem Here a strong nonlinear problem is considered, and its equation of motion is

€ þ stanhðuÞ ¼ 0 u

ð36Þ

where s = 100. The initial displacement and velocity are u0 = 4 and The period T of this system is equal to 1.1419 [32,38,39]. Using a small time step size h(q1-Bathe) = T/40  0.0285, the displacements for q1 = 0 are presented in Figs. 33–34. Fig. 33 shows that all methods can accurately predict displacements at the beginning. But after some time, the errors as shown in Fig. 34, especially the phase errors of the TPSMi method, the q1-Bathe method and the Kim method are much worse than that of the present method. Then using a larger step size h(q1-Bathe) = T/10  0.1140, the displacements for q1 = 0.99 are compared in Fig. 35. It follows that the q1-OTTBIF method still has highly accurate amplitude and phase even after more than 100 periods, and the q1-Bathe method and the Kim method show noticeable phase leads, but the q1-Bathe method is better than the Kim method. Example 5: Mass-spring system

v0 = 0 respectively.

Fig. 39. Acceleration of m1 for different methods (Case 1).

13

14

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210

Table 5 The errors of different methods in example 5. Error

TPSMi

q1-Bathe

Kim

q1-OTTBIF

Case 1

Displacement Velocity Acceleration

1.6479e-3 1.1561e-2 1.0120e-1

1.6827e-3 1.1906e-2 1.0338e-1

1.6829e-3 1.1908e-2 1.0340e-1

1.6197e-3 1.1484e-2 9.9637e-2

Case 2

Displacement Velocity Acceleration

4.3203e-2 2.9179e-1 2.6703e + 0

4.3015e-2 2.9997e-1 2.6506e + 0

4.3019e-2 3.0000e-1 2.6511e + 0

3.9174e-2 2.8930e-1 2.3956e + 0

Case 3

Displacement Velocity Acceleration

1.6717e-1 1.1849e + 0 1.0366e + 1

1.7259e-1 1.2004e + 0 1.1105e + 1

1.7263e-1 1.2006e + 0 1.1108e + 1

1.6116e-1 9.6542e-1 9.1320e + 0

In this example, as shown in Fig. 36, a nonlinear two DOF massspring system subjected to a harmonic external excitation is considered. The mass of particles and stiffness of springs are m1 = m2 = 100 kg, k1 = 103  (1 + 9u21 )N/m, k2 = 103N/m, and the external load is F = 1000sin(5 t)N. Initial conditions are assumed to be zero. The equations of motion of the system have the forms as

€ 1 þ 2000u1  1000u2 þ 9000u31 ¼ 1000sinð5tÞ 100u € 2 þ 1000u2  1000u1 ¼ 0 100u

ð37Þ

The parameter q1=0.99 is used for all the methods in this example, and three different step sizes are respectively considered, including h (TPSMi) = 0.005s (Case 1), h (TPSMi) = 0.05s (Case 2) and h(TPSMi) = 0.1s (Case 3). Figs. 37–39 show the results of m1 for Case 1. It follows that there is no observable differences between these methods. To clearly see the differences, the absolute errors for the interval [0, 5s], refer to Example 2, are presented in Table 5, showing that the q1-OTTBIF method is the most accurate for three cases. As well known, in the recursive calculations of dynamic responses, errors are accumulated step by step, indicating that small accuracy advantage of a time integration method is also meaningful. 6. Conclusions The paper proposed an optimized three-sub-step composite time integration methods with controllable numerical dissipation, called the q1-OTTBIF method. The present method is motivated by two-sub-step q1-Bathe method, but the amplitude decay and period elongation were improved through using additional substep and optimizing algorithmic parameters. Compared with the OTTBDF method which is an optimized three-sub-step composite method, the q1-OTTBIF method uses the four-point backward interpolation or the Newmark-like formulation in the last sub-step rather than the four-point backward difference formulae. The q1-OTTBIF method possesses second-order accuracy, unconditional stability and controllable numerical dissipation, and q1 is used to exactly control the degree of numerical dissipation. When q1 = 0, the properties of the q1-OTTBIF method reduce to those of the OTTBDF method. Spectral analysis and several numerical experiments have validated the advantages of the q1-OTTBIF method in the amplitude decay and period elongation over other similar methods, such as the q1- Bathe method and the Kim method. However, compared with the two-sub-step methods, the present three-sub-step method is more complicated due to the additional sub-step, but only moderately since the TR is used in the first two steps. Acknowledgments The support of the National Natural Science Foundation of China (11872090, 11672019) is gratefully acknowledged.

References [1] Dahlquist G. On accuracy and unconditionally stability of linear multistep methods for second order differential equations. BIT Numer Math 1978;18 (2):133–6. [2] Xie YM, Steven GP. Instability, chaos, and growth and decay of energy of timeintegration schemes for non-linear dynamic equations. Commun Numer Meth En 1994;10:393–401. [3] Zhang HM, Xing YF. A three-parameter single-step time integration method for structural dynamics analysis. Acta Mech Sin 2019;35(1):112–28. [4] Fung TC. Solving initial value problems by differential quadrature method–Part 2: second- and higher-order equations. Int J Numer Meth Engng 2001;50:1429–54. [5] Fung TC. Weighting parameters for unconditionally stable higher-order accurate time step integration algorithms. Part 2: second-order equations. Int J Numer Meth Engng 1999;45:971–1006. [6] Kim W, Reddy JN. A new family of higher-order time integration algorithms for the analysis of structural dynamics. Int J Appl Mech 2017;84:071008. [7] Rezaiee-Pajand M, Alamatian J. Implicit higher-order accuracy method for numerical integration in dynamic analysis. J Struct Eng 2008;134(6):973–85. [8] Rezaiee-Pajand M, Sarafrazi SR, Hashemian M. Improving stability domains of the implicit higher order accuracy method. Int J Numer Meth Engng 2011;88:880–96. [9] Xing YF, Qin MB, Guo J. A time finite element method based on the differential quadrature rule and Hamilton’s variational principle. Appl Sci 2017;7(138): app7020138. [10] Qin MB, Xing YF, Guo J. An improved differential quadrature time element method. Appl Sci 2017;7(471):app7050471. [11] Hughes TJR, Caughey TK, Liw WK. Finite-element methods for nonlinear elastodynamics which conserve energy. Int J Appl Mech 1978;45:366–70. [12] Kuhl D, Ramme E. Constraint energy momentum algorithm and its application to nonlinear dynamics of shells. Comput Methods Appl Mech Engrg 1996;136:293–315. [13] Simo JC, Wong KK. Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum. Int J Numer Meth Eng 1991;31:19–52. [14] Simo J, Tarnow N. The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics. J Appl Math 1992;43:747–92. [15] Gonzalez O. Time integration and discrete Hamiltonian systems. J Nonlinear Sci 1996;6:449–67. [16] Zhang R, Zhong HZ. A quadrature element formulation of an energymomentum conserving algorithm for dynamic analysis of geometrically exact beams. Comput Struct 2016;165:96–106. [17] Krenk S. Global format for energy-momentum based time integration in nonlinear dynamics. Int J Numer Meth Engng 2014;00:1–21. [18] Bathe KJ, Baig MMI. On a composite implicit time integration procedure for nonlinear dynamics. Comput Struct 2005;83:2513–24. [19] Bank RE, Coughran WM, Fichtner W, Grosse EH, Rose DJ, Smith RK. Transient simulation of silicon devices and circuits. IEEE Trans CAD 1985;4:436–51. [20] Bathe KJ. Conserving energy and momentum in nonlinear dynamics: a simple implicit time integration scheme. Comput Struct 2007;85:437–45. [21] Bathe KJ, Noh G. Insight into an implicit time integration scheme for structural dynamics. Comput Struct 2012;98–99:1–6. [22] Noh G, Ham S, Bathe KJ. Performance of an implicit time integration scheme in the analysis of wave propagations. Comput Struct 2013;123(9):93–105. [23] Chandra Y, Zhou Y, Stanciulescu I, Eason T, Spottswood S. A robust composite time integration scheme for snap-through problems. Comput Mech 2015;55:1041–56. [24] Wen WB, Wei K, Lei HS, Duan SY, Fang DN. A novel sub-step composite implicit time integration scheme for structural dynamics. Comput Struct 2017;182:176–86. [25] Zhang HM, Xing YF. Optimization of a class of composite method for structural dynamics. Comput Struct 2018;202:60–73. [26] Xing YF, Ji Y, Zhang HM. On the construction of a type of composite time integration methods. Comput Struct 2019;221:157–78. [27] Malakiyeh MM, Shojaee S, Bathe KJ. The Bathe time integration method revisited for prescribing desired numerical dissipation. Comput Struct 2019;212:289–98.

Y. Ji, Y. Xing / Computers and Structures 231 (2020) 106210 [28] Noh G, Bathe KJ. The Bathe time integration method with controllable spectral radius: The q1-Bathe method. Comput Struct 2019;212(9):299–310. [29] Huang C, Fu HM. A composite collocation method with low-period elongation for structural dynamics problems. Comput Struct 2018;195:74–84. [30] Fung TC. Unconditionally stable higher-order accurate collocation time-step integration algorithms for first-order equations. Comput Methods Appl Mech Engrg 2000;190:1651–62. [31] Kim W, Reddy JN. An improved time integration algorithm: a collocation time finite element approach. Int J Struct Stab Dyn 2016;16(10):1–38. [32] Kim W, Choi SY. An improved implicit time integration algorithm: The generalized composite time integration algorithm. Comput Strcut 2018;196:341–54. [33] Rezaiee-Pajand M, Sarafrazi SR. A mixed and multi-step higher-order implicit time integration family. ARCHIVE Proc Inst Mech Eng Part C J Mech Eng Sci 2010;1(1):1–12.

15

[34] Bathe KJ. Finite Element Procedures. Prentice Hall; 1996. [35] Geradin M, Rixen DJ. Mechanical vibrations: theory and application to structural dynamics. John Wiley & Sons; 2015. [36] Noh G, Bathe KJ. Further insights into an implicit time integration scheme for structural dynamics. Comput Struct 2018;202:15–24. [37] Noh G, Bathe KJ. For direct time integrations: A comparisons of the Newmark and -Bathe schemes. Comput Struct 2019;225:106079. [38] Xie YM. An assessment of time integration schemes for nonlinear dynamic equations. J Sound Vib 1996;192(1):321–31. [39] Wood WL, Oduor ME. Stability properties of some algorithms for the solution of nonlinear dynamic vibration equations. Int J Numer Meth Biomed Eng 1988;4(2):205–12.