Enhanced studies on the composite sub-step algorithm for structural dynamics: The Bathe-like algorithm

Enhanced studies on the composite sub-step algorithm for structural dynamics: The Bathe-like algorithm

Journal Pre-proof

Enhanced studies on the composite sub-step algorithm for structural dynamics: The Bathe-like algorithm Jinze Li, Xiangyang Li, Kaiping Yu PII: DOI: Reference:

S0307-904X(19)30715-2 https://doi.org/10.1016/j.apm.2019.11.033 APM 13163

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

12 April 2019 6 October 2019 19 November 2019

Please cite this article as: Jinze Li, Xiangyang Li, Kaiping Yu, Enhanced studies on the composite sub-step algorithm for structural dynamics: The Bathe-like algorithm, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.11.033

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights • Further studies on the Bathe algorithm is given and hence a novel class of the Bathe-like algorithm is presented. • It provides a new family of algorithms with the optimal numerical characteristics. • It gives a novel family of algorithms, which can be regarded as the alternative to the original Bathe algorithm. • This study shows that the Bathe algorithm can reduce to the trapezoidal rule and the backward Euler formula.

1

Enhanced studies on the composite sub-step algorithm for structural dynamics: The Bathe-like algorithm Jinze Lia , Xiangyang Lia , Kaiping Yua,∗ a Department

of Astronautic Science and Mechanics, Harbin Institute of Technology, No.92 West Dazhi Street, Harbin 150001

Abstract The Bathe algorithm is superior to the trapezoidal rule in solving nonlinear problems involving large deformations and long-time durations. Generally, the √ parameter γ = 2 − 2 is highly recommended due to its optimal numerical properties. This paper further studies this implicit composite sub-step algorithm and thus presents a class of the Bathe-like algorithm. It not only gives a novel family of composite algorithms whose numerical properties are the ex√ actly same as the original Bathe algorithm with γ = 2 − 2, but also provides the generalized alternative to the original Bathe algorithm with any γ. In this study, it has been shown that the Bathe-like algorithm, including the original Bathe algorithm, can reduce to two common single-step algorithms: the trapezoidal rule and the backward Euler formula. Besides, a new parameter called the algorithmic mode truncation factor is firstly defined to describe the numerical property of the Bathe-like algorithm and it can estimate which modes to be damped out. Finally, numerical experiments are provided to show the superiority of the Bathe-like algorithm over some existing methods. For example, the novel Bathe-like algorithms are superior to the original Bathe algorithm when solving the highly nonlinear pendulum. Keywords: Sub-step composite algorithm, Bathe algorithm, Structural dynamics, Implicit time integration, Algorithmic mode truncation factor ∗ Corresponding

author Email addresses: [email protected] (Jinze Li), [email protected] (Xiangyang Li), [email protected] (Kaiping Yu)

Preprint submitted to Applied Mathematical Modelling

November 21, 2019

2010 MSC: 00-01, 99-00

1. Introduction Transient analysis in structural dynamics requires to obtain the accurate and stable numerical solutions of the following equation ¨ ˙ MU(t) + CU(t) + KU(t) = F(t)

(1)

where M, C and K are mass, damping and stiffness matrices, respectively. F(t) is the external load vector, known as a function of time t; U(t) is the displacement vector and superposed dots imply differential with respect to time t. The Houbolt method [1] is one of the earliest algorithms to solve equation (1), and the Newmark-β algorithm [2] is the widely used direct integration algorithm since almost all open source finite element programs, such as the FEAPpv[3] and Calculix [4], achieve it well. The Newmark algorithm includes the Newmark explicit method and constant average acceleration method as the special cases. In other literatures [5, 6], these two algorithms are also called as the central difference method and trapezoidal rule, respectively. So far, a number of researchers [6–23] have developed various direct integration methods based on different design ideas. In general, these numerical integration schemes can be split into two groups: explicit and implicit algorithms. The central difference method [5] is the oldest and widely used explicit algorithm. The main advantage of the explicit method is significantly saving the computational cost because of requiring no solution of linear equations. Further, it is more significant than the implicit algorithm for nonlinear problems due to no nonlinear iterations within each step. However, the traditional explicit method is generally conditionally stable. Its integration step must be limited by its stability condition. As an example, the time step h of the central difference method should satisfy the condition ω · h ≤ 2 where the ω is the largest natural frequency of the system analyzed. Hence, for the explicit methods, their integration steps are generally chosen to be very small to obtain stable

3

numerical responses. Recently, the so-called explicit/semi-explicit algorithms [6, 16, 20, 24–26] have been proposed to remedy the conditional stability of the traditional explicit method. These algorithms employ the structure-dependent algorithmic parameters to achieve the unconditional stability for the linear elastic and nonlinear stiffness softening systems. But compared with the standard explicit method, structure-dependent algorithms still need the factorization of an effective stiffness matrix [27]. The trapezoidal rule [5] and CH-α algorithm [10] are the widely used nondissipative and dissipative implicit integration methods, respectively. In contrast to the explicit method, the implicit algorithm can generally achieve unconditional stability. As a result, they can use a relatively larger time step to save the computational cost. However, there are two shortcomings when using a larger integration step. One is that a larger time step for the common implicit algorithms could present more undesirable numerical dissipations and dispersions in the important low-frequency region, leading to lower solution accuracy; the other one is that more nonlinear iterations are required to converge to the satisfactory solution in each step. Recently, the novel implicit composite sub-step algorithms [12–14, 21–23, 28–32] are proposed to overcome the above two problems well. The original intention of the composite sub-step algorithm [12, 13] is to effectively solve some nonlinear problems which the non-dissipative trapezoidal rule cannot solve well. The Bathe algorithm [12, 13] is the most representative of these composite sub-step algorithms. The Bathe algorithm is a single step but two sub-step composite algorithm where the trapezoidal rule and three-points backward differential formula [33] are employed in the first and second sub-step, respectively. The Bathe algorithm achieves the maximum of numerical dissipations (i.e. Lstability)and another two Bathe-type algorithms [21, 22] with controllable numerical dissipations have been developed. Note that the second-order accurate β1 /β2 -Bathe algorithm [21] is essentially the ρ∞ -Bathe algorithm [22] with γ = 0.5 and the original Bathe algorithm [12, 13] corresponds to the ρ∞ -Bathe algorithm [22] with ρ∞ = 0.0. In other words, three Bathe-type algorithms 4

[13, 21, 22] (the original Bathe algorithm, the β1 /β2 -Bathe and ρ∞ -Bathe algorithm) will reduce to the identical integration scheme (i.e. the original Bathe algorithm) when only considering the maximum of numerical dissipations and second-order accuracy. Hence, the L-stable Bathe algorithm [12, 13] is the focus of attention in this study. The property of controllable numerical dissipations [34] is always considered as the desirable numerical characteristic of the direct integration algorithm. But this paper shows that the algorithm with controllable numerical dissipations only controls how fast the spurious modes are eliminated by adjusting its algorithmic parameter (generally ρ∞ ). For the given time step, the controllably dissipative algorithm cannot depend on which modes to be damped out by adjusting ρ∞ . The only approach to control which modes to be damped out is to select the appropriate integration step, which is obviously impractical in engineering applications. In this study, further studies on the original Bathe algorithm are presented and hence a novel composite sub-step algorithm, called as the Bathe-like algorithm, is given by introducing a new parameter θ. The present scheme uses the generalized trapezoidal rule [5] in the first sub-step and the three-points backward differential formula with different coefficients is employed in the second sub-step. By the spectral analysis, the present scheme can provide an optimal family of composite sub-step algorithms, whose numerical characteristics are √ the exactly same as the Bathe algorithm with γ = 2 − 2. In other words, this optimal family of algorithms can be considered as the alternative to the Bathe √ algorithm with γ = 2 − 2. Moreover, for given θ 6= 0, the present scheme can achieve the same numerical properties as the original Bathe algorithm with any γ, which is the reason why the present scheme is named as the Bathe-like algorithm. Particularly, in the case of θ = 1/2, the present scheme reduces to the original Bathe algorithm and in the case of θ = 1, it can be regarded as the generalized version of the algorithm proposed in [23]. Importantly, it has been shown in this paper that the Bathe-like algorithm, including the original Bathe algorithm, can reduce to the trapezoidal rule [5] and the backward Euler 5

formula [33]. In order to describe numerical properties of the Bathe-like algorithm well, a new parameter, algorithmic mode truncation factor, is also defined in this paper. The remaining part of this paper is organized as follows. The Bathe-like algorithm and its accuracy analysis are presented in Subsection 2.1 and 2.2, respectively. A new parameter is defined as the algorithmic mode truncation factor and its discussion is presented in Section 3, where the algorithmic mode truncation factors of traditional implicit dissipative algorithms are also calculated. In Section 4, four different cases of the algorithmic parameter are analyzed in terms of the spectral radius, numerical dissipation and dispersion. The special cases and members from the Bathe-like algorithm are presented in Section 5, where the Bathe-like algorithm can reduce to the trapezoidal rule, backward Euler formula, and original Bathe algorithm. Section 6 gives the implementation of the Bathe-like algorithm for solving nonlinear problems. Four numerical examples are carried out in Section 7 to demonstrate the superiority of the Bathe-like algorithm over other existing methods. Conclusions are drawn in Section 8.

2. Composite sub-step algorithm and accuracy analysis 2.1. The Bathe-like algorithm In the Bathe-like algorithm, the entire time step h is still divided into two sub-steps as [t, t + γh] and [t + γh, t + h]. In the first sub-step [t, t + γh], the generalized trapezoidal rule [5] is used to obtain numerical responses at time t + γh, which is ˙ t + θU ˙ t+γh ] Ut+γh = Ut + γh[(1 − θ)U ˙ t+γh = U ˙ t + γh[(1 − θ)U ¨ t + θU ¨ t+γh ] U

(2)

and in the second sub-step [t + γh, t + h], the three-points backward differential formula with different coefficients is employed, ˙ t+h = c2 Ut+h + c1 Ut+γh + c0 Ut hU ¨ t+h = c2 U ˙ t+h + c1 U ˙ t+γh + c0 U ˙t hU 6

(3)

where coefficients ci , i = 0, 1, 2 are determined by c0 = −

2γ 2 θ − 2γ + 1 , 2γ 2 θ − γ

c1 =

1 , 2γ 2 θ − γ

c2 =

2γθ − 2 2γθ − 1

(4)

to ensure second-order accuracy in its final form. Note that when the parameter θ = 1/2 is chosen, the Bathe-like algorithm obviously reduces to the original Bathe algorithm [12–14]. In the case of θ = 1, the first-order accurate backward Euler formula [33] is used in the first sub-step but coefficients (4) still ensure the second-order accuracy of numerical solutions at time t + h. In this case, the present scheme can be regarded as the generalized version of the algorithm proposed in [23]. Of course, the motion of equation (1) at time t + γh and t + h is also needed to obtain the acceleration responses. ¨ t+γh + CU ˙ t+γh + KUt+γh = Ft+γh MU ¨ t+h + CU ˙ t+h + KUt+h = Ft+h MU

(5) (6)

With Eq.s (2)-(6), the Bathe-like algorithm can be rewritten as ˜ 1U ¨ t+γh = F ˜1 K

(7)

˜ 2U ¨ t+h = F ˜2 K

(8)

˜ 1 = (γθh)2 K + γθhC + M K  2 h ˜2 = h K K+ C+M c2 c2

(9)

where

(10)

and ˜ 1 =Ft+γh − C[U ˙ t + (1 − θ)γhU ¨ t ] − K[Ut + γhU ˙ t + θ(1 − θ)(γh)2 U ¨ t ] (11) F   1 1 1 ˙ ˙ ˜ ˙ ˙ F2 =Ft+h + C(c0 Ut + c1 Ut+γh ) + K (c0 Ut + c1 Ut+γh ) + 2 h(c0 Ut + c1 Ut+γ ) c2 c2 c2 (12) ˜ 1 and K ˜ 2) It is necessary to point out that two effective stiffness matrices (K become identical when γθc2 = 1 is satisfied. In this paper, Further analyses of parameters γ and θ will be presented in the subsequent sections. 7

2.2. Consistency and accuracy For simplicity, the free vibration of an undamped single-degree-of-freedom (SDOF) system is considered, which is u ¨(t) + ω 2 u(t) = 0

(13)

where ω is the natural frequency of the SDOF system. Numerical responses at time t + h can be written as the following matrixvector form when the Bathe-like algorithm is applied to the SDOF system (13), 

ut+h





ut



         hu˙ t+h  = A  hu˙ t      h2 u ¨t+h h2 u ¨t

(14)

where matrix A is the amplification matrix of the Bathe-like algorithm and its explicit expression is given in appendix A. The characteristic polynomial of matrix A is obtained by calculating det(λI − A) = 0 as λ3 − 2A1 λ2 + A2 λ = 0

(15)

where λ is the eigenvalue of matrix A and coefficients A1,2 read 4Ω2 (γθ)4 − 8Ω2 (γθ)3 + (6Ω2 + 4)(γθ)2 − 8γθ − Ω2 + 4 (Ω2 (γθ)2 + 1)[(4Ω2 + 4)(γθ)2 − (4Ω2 + 8)γθ + Ω2 + 4] 4Ω2 (γθ)4 − 8Ω2 (γθ)3 + (8Ω2 + 4)(γθ)2 − (4Ω2 + 8)γθ + Ω2 + 4 A2 = (Ω2 (γθ)2 + 1)[(4Ω2 + 4)(γθ)2 − (4Ω2 + 8)γθ + Ω2 + 4]

A1 =

(16)

where Ω = ω · h. It is well-known that numerical characteristics of direct integration algorithms, including accuracy and stability conditions, depend only on the eigenvalues of the amplification matrix A [5]. Hence, equations (15)-(16) indicate that numerical properties of the Bathe-like algorithm are determined only by the product of γ and θ. Hence, another parameter µ is introduced to represent the product of γ and θ, which is µ=γ·θ 8

(17)

In Eq.(15), the velocity and acceleration can be eliminated to obtain a difference equation with respect to displacements [34], ut+h − 2A1 ut + A2 ut−h = 0

(18)

where coefficients A1,2 are given in Eq.(16). The local truncate error (l.t.e.) of displacement ut+h can be defined as l.t.e. =

1 [u(t + h) − 2A1 u(t) + A2 u(t − h)] h2

(19)

where u(t + jh), j = −1, 0, 1 are the exact displacements at time t + jh, j = −1, 0, 1, respectively. Substituting the Taylor series expansion of displacements u(t + h) and u(t − h) at time t and coefficients A1,2 into the above yields l.t.e. =

1 {4(µ − 1)2 [¨ u(t) + ω 2 u(t)] + O(h2 )} (Ω2 µ2 + 1)[(2µ − 1)2 Ω2 + 4(µ − 1)2 ] (20)

Using the equation of motion at time t, equation (20) shows that the Bathe-like algorithm achieves the second-order accuracy in time for any µ. In fact, in the case of µ = 1, the Bathe-like algorithm will reduce to the first-order accurate backward Euler formula, which is analyzed in Section 5. Considering the following linear SDOF system is to examine the numerical accuracy of the Bathe-like algorithm. u ¨(t) + 4u(t) ˙ + 5u(t) = sin(2t),

u(0) =

57 2 , u(0) ˙ = 65 65

(21)

whose exact displacement solution is u(t) = exp(−2t)(cos(t) + 2 sin(t)) −

1 (8 cos(2t) − sin(2t)) 65

(22)

Time duration of t = 2T ≈ 5.62s is considered and calculated for the convergence test. The global error norm [28] is used to compute errors in the displacement, which is Edis

v uN N uX X = t [ui − u(ti )]2 / [u(ti )]2 i=0

(23)

i=0

where ui and u(ti ) are the numerical and exact solution at time ti , respectively, and N is the total number of time steps. Similar expressions of Evel and Eacc can be easily obtained for the velocity and acceleration. 9

E rro r in d is p la c e m e n t

µ= 0 . 2 9 & µ= 0 . 5 5 & 1 0

-2

1 0

-3

1 0

-4

1 0

-5

1 0

-6

1 0

-7

θ= 0 . 5 θ= 0 . 5

µ= 0 . 2 9 & µ= 0 . 5 5 &

θ= 1 . 0 θ= 1 . 0

µ= 0 . 2 9 & µ= 0 . 5 5 &

θ= 1 . 5 θ= 1 . 5

µ= - 0 . 1 5 & θ= 0 . 5 B a t h e γ= 2 - √2

µ= - 0 . 1 5 & θ= 1 . 0 T ra p e z o id a l ru le

µ= - 0 . 1 5 & θ= 1 . 5 C H - α. ρ∞ = 0 . 0

2

1 0

-3

1 0

-2

1 0

-1

h

Figure 1: Accuracy and convergence rates of various algorithms in displacement.

E rro r in v e lo c ity

µ= 0 . 2 9 & µ= 0 . 5 5 &

1 0

-2

1 0

-3

1 0

-4

1 0

-5

1 0

-6

θ= 0 . 5 θ= 0 . 5

µ= 0 . 2 9 & µ= 0 . 5 5 &

θ= 1 . 0 θ= 1 . 0

µ= 0 . 2 9 & µ= 0 . 5 5 &

θ= 1 . 5 θ= 1 . 5

µ= - 0 . 1 5 & θ= 0 . 5 B a t h e γ= 2 - √2

µ= - 0 . 1 5 & θ= 1 . 0 T ra p e z o id a l ru le

µ= - 0 . 1 5 & θ= 1 . 5 C H - α. ρ∞ = 0 . 0

2

1 0

-3

1 0

-2

1 0

-1

h

Figure 2: Accuracy and convergence rates of various algorithms in velocity.

Three error norms for different algorithms versus the step size h on the log-log scale are presented in Fig.s 1-3. It is obvious to find that the Bathe-like algorithm achieves the second-order accuracy and it produces the smallest errors √ in displacement, velocity and acceleration in the case of µ = 0.29 ≈ 1 − 2/2. Furthermore, for given µ, the Bathe-like algorithm with larger θ can give the mildly smaller global errors in displacement, velocity and acceleration than the present scheme with θ = 1/2 (i.e. the original Bathe algorithm), although these differences are vary small and ignorable. It means that the integration scheme in the first sub-step is not necessary to limit into the second-order accurate trapezoidal rule and can be replaced by other numerical schemes, such as the

10

E rro r in a c c e le ra tio n

µ= 0 . 2 9 & µ= 0 . 5 5 &

1 0

-1

1 0

-2

1 0

-3

1 0

-4

1 0

-5

1 0

-6

θ= 0 . 5 θ= 0 . 5

µ= 0 . 2 9 & µ= 0 . 5 5 &

θ= 1 . 0 θ= 1 . 0

µ= 0 . 2 9 & µ= 0 . 5 5 &

θ= 1 . 5 θ= 1 . 5

µ= - 0 . 1 5 & θ= 0 . 5 B a t h e γ= 2 - √2

µ= - 0 . 1 5 & θ= 1 . 0 T ra p e z o id a l ru le

µ= - 0 . 1 5 & θ= 1 . 5 C H - α. ρ∞ = 0 . 0

1

2

1 0

-3

1 0

-2

1 0

-1

h

Figure 3: Accuracy and convergence rates of various algorithms in acceleration.

backward Euler scheme [23].

3. Algorithmic mode truncation factor Herein, a new parameter is firstly defined as the algorithmic mode truncation factor Ωtr to describe the numerical dissipative property of the Bathe-like algorithm. Definition 1 (Algorithmic Mode Truncation Factor). The algorithmic mode truncation factor of an unconditionally stable direct integration algorithm is defined as the dimensionless frequency Ωtr = ωtr · h, where the second derivative of the spectral radius is equal to zero, that is d2 ρ(Ω) =0 dΩ2 Ω=Ωtr

(24)

where ρ(Ω) is the spectral radius, as a function of the dimensionless frequency Ω, of the integration algorithm in the absence of physical damping ratio ξ. Some additional informations indicated by Definition 1 are listed as follows: (i) Equation (24) indicates that the algorithmic mode truncation factor Ωtr is essentially the inflection point of its spectral curve ρ(Ω). (ii) In general, there exists the algorithmic mode truncation factor Ωtr only for the unconditionally stable integration algorithm since the dimensionless 11

frequency Ω of the conditionally stable method is limited by its stability condition. For example, Ωmax = ωmax h ≤ 2 for the central difference method [5]. (iii) In Definition 1, it is required that the integration algorithm should possess the smooth spectral radius curve to fulfill the mathematical requirement of differential twice of ρ(Ω). Hence, there is no algorithmic mode truncation factor Ωtr for the Wilson-θ algorithm [7] with θ = 1.4 since there exists a point in its spectral curve where its first derivative is not continuous (see Fig.9.3.1 in [5]). (iv) The algorithmic mode truncation factor Ωtr is defined as the infinity for the zero-dissipative integration algorithms, such as the trapezoidal rule, since their spectral radii are always equal to one for any value of Ω. Note that the algorithmic mode truncation factor Ωtr is introduced herein to determine which modes to be damped out, while the classical controllably dissipative parameter ρ∞ only controls how fast the modes to be eliminated. As an example, the CH-α [10] and G-α algorithm [11] are such algorithms that they possess the controllable numerical dissipations in the high-frequency range but the constant algorithmic mode truncation factor Ωtr . On the contrary, the Bathe-like algorithm achieves the maximum numerical dissipation in the highfrequency region and the controllable algorithmic mode truncation factor Ωtr . Hence, when solving multi-degree-of-freedom systems, the unique way to depend on which modes to be damped out for the CH-α and G-α algorithm is to adjust their time step h, which is not actually practical in engineering applications. But, the Bathe-like algorithm can achieve this goal by simply adjusting µ. 3.1. The algorithmic mode truncation factor of the Bathe-like algorithm Herein, the functional relationship between the parameter µ and Ωtr is analyzed in order to compute the value of µ for given Ωtr . With Eq.(15), the spectral radius ρ(Ω) of the Bathe-like algorithm can be given by ρ(Ω) = max(|λi |, i = 1, 2, 3) = 12

p

A2

(25)

where A2 is from the second line of Eq.(16). Squaring on both sides of Eq.(25) and differentiating it with respect to Ω twice yield ρ(Ω)2 = A2 dA2 dρ(Ω) = dΩ dΩ 2  d2 ρ(Ω) d2 A2 dρ(Ω) + 2ρ(Ω) = 2 dΩ dΩ2 dΩ2 2ρ(Ω)

(26) (27) (28)

Reorganizing these equations yields dρ(Ω) 1 dA2 = √ dΩ 2 A2 dΩ "  2 # 2 d ρ(Ω) 1 d2 A2 1 dρ(Ω) = √ − dΩ2 2A2 dΩ 2 A2 dΩ2

(29) (30)

Hence, the algorithmic mode truncation factor Ωtr can be given by solving " 2 #  d2 A2 1 dρ(Ω) − =0 (31) dΩ2 2A2 dΩ Ω=Ωtr

Substituting the second line of (16) into the above gives

a8 Ω8tr + a6 Ω6tr + a4 Ω4tr + a2 Ω2tr + a0 = 0

(32)

where coefficients aj are determined by 1 2 µ (2µ − 1)2 (2µ2 − 2µ + 1)4 56 1 a6 = − (4µ4 − 20µ3 + 26µ2 − 10µ + 1)(4µ4 + 128 a8 =

(33a)

4µ3 − 10µ2 + 2µ + 1)(2µ2 − 2µ + 1)2

(33b)

− 44µ3 + 26µ2 − 8µ + 1)(µ − 1)2

(33c)

1 a4 = − (16µ8 − 64µ7 + 104µ6 − 88µ5 + 58µ4 8 9 a2 = − (2µ2 − 2µ + 1)2 (µ − 1)4 8

(33d)

a0 = − 3(µ − 1)6

(33e)

Obviously, it is almost impossible to obtain the analytical expression of µ on Ω. But, numerical solutions of Eq.(32) can be easily calculated for given Ωtr , see Table 1. 13

10

3

10

2

10

1

10

0

10

-1

10

-2

10

-3

10

-4

10

-5

10

-6

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

tr

Figure 4: The relationship between parameter µ and Ωtr in the case of µ > 0.

-10 -6 -10 -5 -10 -4 -10 -3

µ

-10 -2 -10 -1 -10 0 -10 1 -10 2 -10 3 10-2

10-1

100

101

102

103

104

105

106

Ωtr

Figure 5: The relationship between parameter µ and Ωtr in the case of µ < 0.

14

Figures 4 and 5 show numerical solutions of Eq.(32) when the Ωtr ranges from 10

−2

to 106 . These figures show that the Bathe-like algorithm can control the

algorithmic mode truncation factor Ωtr from 0.74 (approximately) to infinity by changing µ. Particularly, as the parameter µ tends to 0 or 0.5, the larger value of Ωtr can be obtained. But, when the absolute value of µ continuously increases to infinity, there is no significant change for Ωtr near 0.74 (approximately). Therefore, the interval of µ is limited into µ ∈ [−10, 1] in the following analysis. The spectral radii and their first derivative curves of the Bathe-like algorithm with various µ are plotted in Fig.6. It is observed in Fig.6 that the first derivative of ρ(Ω) achieves the minimum at Ωtr . It is interesting to find that in Fig.6 the spectral radii ρ(Ω) give the almost same value (about 0.8) at Ωtr , independent of µ. Note that there are still some numerical dissipations in the small region Ω ≤ Ωtr , since spectral radii in this small region are less than unity. 1.0

1.0

(a)

1.0

= -10

(b)

= -0.001

(c)

= 0.0001

0.8 0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.6

0.4

0.2

0.0

-0.2

0.0

0.0

-0.4 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

1.0

1.0

(d)

= 1-

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

1.0

2/2

(e)

(f)

= 0.45

= 1.00

0.8 0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.6

0.4

0.2

0.0

-0.2 0.0

0.0

-0.4 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

Figure 6: Spectral radii and their first derivative curves of the Bathe-like algorithm with various µ. The solid line indicates the spectral radii and the dotted line presents its first derivative.

15

10

5

10

6

Table 1: Numerical solutions of equation (32).

µ ∈ (0, 1 −

√

µ ∈ [−10, 0)

µ ∈ (0.5, 1]

Ωtr

0.74

-8.14941

0.94535

4.13

0.22878

0.35168

1

-1.86140

0.82526

5

0.15980

0.40490

2

-0.50258

0.66724

8

0.09008

0.45050

4

-0.19178

0.58046

16

0.04424

0.47685

8

-0.08926

0.54097

25

0.02829

0.48544

16

-0.04423

0.52118

50

0.01414

0.49283

50

-0.01414

0.50697

100

0.00707

0.49644

100

-0.00707

0.50351

200

0.00354

0.49823

300

-0.00236

0.50118

300

0.00236

0.49882

2/2]

µ ∈ (1 −

√

Ωtr

2/2, 0.5)

Table 1 presents numerical results of µ for some given values of Ωtr , where the interval of µ is divided into four subintervals (Note that the µ = 0 and µ = 0.5 are not included in four subintervals due to Eq.(4). In fact, these two cases correspond to the non-dissipative trapezoidal rule). In each subinterval of µ, the controllable algorithmic mode truncation factor Ωtr can be achieved. Minimums of Ωtr are approximately 0.74 and 3.92 in the cases of µ ∈ [−10, 0) ∪ (0.5, 1] and µ ∈ (0, 0.5), respectively. These results can be observed in Fig.s 4 and 5. In Section 4, numerical characteristics of the Bathe-like algorithm in four subintervals of µ will be analyzed and compared in detail. 3.2. The algorithmic mode truncation factor of traditional implicit dissipative algorithms In general, the amplification matrix of a self-starting integration algorithm can be expressed as the 3-dimensional square matrix. Note that if the direct integration algorithm, such as the Newmark-β algorithm [2], uses equilibrium at step ends, then its 3-dimensional amplification matrix can be reduced to the 2-dimensional square matrix [35]. In this case, the above procedure used in the Bathe-like algorithm to calculate Ωtr still works. However, if the integration algorithm, such as the G-α algorithm [11], employs the weighted combinations 16

of equilibrium, then a 3 × 3 amplification matrix is adequate such that the Eq.(31) is not suitable. The work presented in this subsection focuses on this case. Consider the fact that the G-α algorithm [11] not only includes three wellknown traditional implicit dissipative algorithms (i.e. the CH-α [10], HHT-α [8], and WBZ-α algorithm [9]) as the special cases, but also develops their corresponding no-overshoot versions (i.e. the NOCH-α, NOHHT-α, and NOWBZ-α algorithm). Hence, particular attention is paid to the G-α algorithm to show that their algorithmic mode truncation factors are not adjusted by users. There are three nonzero eigenvalues for the amplification matrix of the G-α algorithm [11] so that it is very difficult and complex to explicitly derivate the similar expression (31). Hence, some numerical methods are used to calculate its algorithmic mode truncation factor, as shown in Fig. 7. Obviously, when the unique parameter ρ∞ varies, their algorithmic mode truncation factors of these six controllably dissipative algorithms almost don’t move in the Ω-direction. It means that for these six algorithms, their algorithmic mode truncation factors are not adjusted by choosing the ρ∞ . 1.00

1.00 = 1.00

= 1.00

= 0.75

= 0.75

1.00 = 1.00

(a) NOCH/CH0.75

0.75

= 0.75

(b) NOHHT/HHT-

0.75 (c) NOWBZ/WBZ= 0.50

= 0.50

0.50

= 0.50

0.50

= 0.00

0.50

= 0.00

= 0.25

= 0.25 = 0.25

= 0.50

0.25

= 0.75 = 1.00

= 1.00

= 0.00

0.00

= 0.00

0.00

-0.25

0.00

-0.25 -2

10

-1

10

0

10

1

10

2

10

3

10

4

10

= 0.25

0.25

= 0.75

= 1.00

10

= 0.50

= 0.50

0.25 = 0.75

-0.25 -2

10

-1

10

0

10

1

10

2

10

3

10

4

10

-2

10

-1

10

0

10

1

10

Figure 7: Spectral radii and their first derivative curves of the G-α algorithm. The solid line indicates the spectral radii and the dotted line presents its first derivative. Note that NOHHT/HHT-α algorithm are unconditionally stable only for ρ∞ ≥ 0.5.

17

2

10

3

10

4

4. Numerical characteristics Actually, the reason why the interval of µ can be divided into four subintervals is further given in Fig.8a, where the Bathe-like algorithm with any µ achieves the unconditional stability since spectral radii are always less than one. Similar results can be found for damped systems, see Fig.8b. Therefore, the in√ √ terval of µ is divided into four subintervals: [−10, 0), (0, 1− 2/2], (1− 2/2, 0.5) and (0.5, 1] in order to better explore numerical characteristics of the Bathe-like algorithm. The cases of µ = 0 and µ = 1/2 will be analyzed in Section 5. 1.0

1.0

0.4

0.2

0.0 -2.0

0.8

-1.5

-1.0

)

h/T=0.05 h/T=0.10 h/T=0.20 h/T=0.30 h/T=0.50 h/T=0.75 h/T=1.00 h/T=2.00 h/T=10.00 h/T=100.0

0.6

(

0.6

(

)

0.8

0.4

0.2

0.0 -0.5

0.0

0.5

1.0

1.5

2.0

(a)

-2.0

h/T=0.05 h/T=0.10 h/T=0.20 h/T=0.30 h/T=0.50 h/T=0.75 h/T=1.00 h/T=2.00 h/T=10.00 h/T=100.0 -1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

(b)

Figure 8: Spectral radii of the Bathe-like algorithm with different values of µ: (a) ξ = 0 and (b) ξ = 0.1.

4.1. Case one When considering µ ∈ [−10, 0), spectral radii of the Bathe-like algorithm with various Ωtr are plotted in Fig.9. The spectral radii of CH-α algorithm [10] √ with ρ∞ = 0.0 and the original Bathe algorithm [12, 14] with γ = 2 − 2 are also presented for comparison. Figure 9 indicates that the algorithmic mode truncation factor Ωtr does control the horizontal movement of spectral radii curves in the Ω direction. That is, the newly defined parameter Ωtr can control effectively the accuracy of low-frequency region and hence adjust numerical dissipations in the low-frequency range. As a result, the parameter Ωtr can be used by users to estimate which modes to be damped out. It is necessary to point out that the larger Ωtr > 300 is not given in this paper since these 18

modes are considered as the spurious high-frequency responses and hence must be eliminated. It is interesting to plot spectral radii of the Bathe-like algorithm when the physical damping ratio ξ is assumed to be 0.1, see Fig. 10. Obviously, there are some numerical dissipations in the middle-frequency range due to the effect of ξ. 1.0

tr

0.8

tr

tr

)

tr

0.6

(

tr

tr

tr

0.4 tr

tr

=0.75 =1.00 =2.00 =4.00 =8.00 =16.0 =50.0 =100 =300

CH-

=0.0

0.2 Bathe. =2-

0.0 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

Figure 9: Spectral radii of the Bathe-like algorithm with µ ∈ [−10, 0) when ξ = 0.0.

Figure 11 shows that amplitude decays and period elongations of the Bathelike algorithm with µ ∈ [−10, 0), and numerical results from CH-α algorithm √ with ρ∞ = 0.0 and the Bathe algorithm with γ = 2 − 2 are also included for comparison. Figure 11 illustrates that the smaller values of Ωtr introduce more amplitude decays in the low-frequency region to eliminate the corresponding modes (which is also indicated in Fig.9). Compared with CH-α algorithm, most values of Ωtr can provide acceptable numerical accuracy in terms of amplitude decays and period errors. It may be observed that period errors of the Bathelike algorithm decrease into the specific curve (it belongs to the trapezoidal rule) with the increase of Ωtr . In this case µ ∈ [−10, 0), the Bathe algorithm with 19

1.0

tr

0.8

tr

tr

)

tr

0.6

(

tr

tr

tr

0.4 tr

tr

=0.75 =1.00 =2.00 =4.00 =8.00 =16.0 =50.0 =100 =300

CH-

=0.0

0.2 Bathe. =2-

0.0 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

Figure 10: Spectral radii of the Bathe-like algorithm with µ ∈ [−10, 0) when ξ = 0.1.

γ =2−

√

2 gives the smallest period error among all algorithms shown herein. √ Later, the same period errors as the Bathe algorithm with γ = 2 − 2 can be found for the present integration scheme. 4.2. Case two √ When considering µ ∈ (0, 1− 2/2], spectral radii of the Bathe-like algorithm with different Ωtr are plotted in Fig. 12. Similarly, the parameter Ωtr can adjust the accuracy of low-frequency region from 4.13 to 300. Figure 13 gives amplitude decays and period elongations of the Bathe-like √ algorithm with µ ∈ (0, 1 − 2/2]. Obviously, all amplitude decay curves of the Bathe-like algorithm are below that of the Bathe algorithm, and all period elongation curves are above that of the Bathe algorithm. As the Ωtr increases from 4.13 to 300, period errors also increase into the constant curve, which essentially belongs to the trapezoidal rule. The CH-α algorithm with ρ∞ = 0.0 achieves the largest amplitude decay and period elongations among all algorithms provided herein. 20

tr

=0.75

tr

=1.00

tr

=2.00

tr

=4.00

tr

=8.00

tr

=16.0

0.8

tr

=50.0

tr

=100

tr

=300

CH-

=0.0

Bathe. =2-

0.6

0.7 0.5

Period elongation

Amplitude decay

0.6

0.5

0.4

0.3

0.4

0.3

0.2

0.2 0.1 0.1

0.0

0.0

0.00

0.25

0.50

0.75

1.00

1.25

1.50

0.00

0.25

0.50

0.75

1.00

1.25

Figure 11: Amplitude decays and period elongations of the Bathe-like algorithm with µ ∈ [−10, 0).

1.0

tr

0.8

tr

tr

)

tr

0.6

(

tr

tr

tr

0.4 tr

tr

=4.13 =5.00 =8.00 =16.0 =25.0 =50.0 =100 =200 =300

CH-

=0.0

0.2 Bathe. =2-

0.0 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

Figure 12: Spectral radii of the Bathe-like algorithm with µ ∈ (0, 1 −

21

10

√

6

2/2].

1.50

tr

=4.13

tr

=5.00

tr

=8.00

tr

=16.0

tr

=16.0

tr

=25.0

tr

=100

tr

=200

tr

=300

CH-

=0.0

Bathe. =2-

0.06 0.14

0.05

Period elongation

Amplitude decay

0.12

0.04

0.03

0.02

0.10

0.08

0.06

0.04 0.01

0.02 0.00 0.00 0.00

0.25

0.50

0.75

1.00

1.25

1.50

0.00

0.25

0.50

0.75

1.00

1.25

Figure 13: Amplitude decays and period elongations of the Bathe-like algorithm with µ ∈ √ 2/2].

(0, 1 −

Figure 14 presents comparisons of spectral radii between some integration algorithms. For the same Ωtr , the Bathe-like algorithm in two different cases provides the identical spectral radii, but gives different amplitude decays and period errors, as shown in Fig.15. Generally, in the case of same Ωtr , the √ Bathe-like algorithm with µ ∈ (0, 1 − 2/2] can produce fewer period errors and amplitude decays than with µ ∈ [−10, 0), such as the case of Ωtr = 8. √ Hence, the Bathe-like algorithm with µ ∈ (0, 1− 2/2] is a better choice to solve most structural dynamical problems compared with µ ∈ [−10, 0). Figure 15 also reveals that period errors and amplitude decays of the Bathe-like algorithm tend to those of the trapezoidal rule [5] as the Ωtr increases. Importantly, Figures 14 and 15 show that in the case of µ = 1 −

√

2/2

(i.e. Ωtr = 3.92), the Bathe-like algorithm possesses the same spectral radius, √ amplitude decay and period error as the Bathe algorithm with γ = 2 − 2. In this case, both the Bathe-like and original Bathe algorithm achieve the identical effective stiffness matrices inside two sub-steps, which saves the computational cost for the analysis of linear systems. Recall that when θ = 1/2, the Bathe-like algorithm reduces to the original Bathe algorithm [12, 13], but when θ is not equal to be 1/2, these two composite algorithms are not the identical, where

22

1.50

1.0

0.6

(

)

0.8

(

tr

(

0.4

tr

(

tr

(

0.2

(

tr

tr

=300)

=300) =8)

=8)

=3.92)

Trapezoidal rule Bathe. =2-

0.0 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

Figure 14: Comparisons of spectral radii from some integration algorithms.

0.05

0.150

(

0.04

Amplitude decay

( ( (

0.03

tr

tr

tr

(

=300)

tr

tr

=300)

(

0.125

=8)

(

=8)

Period elongation

(

=3.92)

Trapezoidal rule Bathe. =2-

0.02

(

0.100 (

tr

tr

tr

tr

tr

=300)

=300) =8)

=8)

=3.92)

Trapezoidal rule

0.075

Bathe. =2-

0.050

0.01 0.025

0.00 0.00

0.000 0.25

0.50

0.75

1.00

1.25

1.50

0.00

0.25

0.50

0.75

1.00

1.25

Figure 15: Comparisons of some integration algorithms with respect to amplitude decays and period elongations

23

1.50

the Bathe-like algorithm can provide more possibility of choosing θ(or γ), but still achieves the same numerical characteristics as the original Bathe algorithm √ with γ = 2 − 2. The Bathe-like algorithm with the larger θ could produce more reasonable numerical solutions when solving some nonlinear systems, as shown in Section 7.3. Further analysis about the θ will be given in Section 5. 4.3. Case three √ When considering the case of µ ∈ (1− 2/2, 0.5), spectral radii of the Bathelike algorithm with different Ωtr are presented in Fig.16. The amplitude decays and period elongations are plotted in Fig.17. From these figures, there are the √ same conclusions as the case of µ ∈ (0, 1 − 2/2]. Hence, a more detailed analysis is not given herein for brevity. 1.0

tr

0.8

tr

tr

)

tr

0.6

(

tr

tr

tr

0.4 tr

tr

=4.13 =5.00 =8.00 =16.0 =25.0 =50.0 =100 =200 =300

CH-

=0.0

0.2 Bathe. =2-

0.0 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

Figure 16: Spectral radii of the Bathe-like algorithm with µ ∈ (1 −

5

√

10

6

2/2, 0.5).

4.4. Case four When considering the last case of µ ∈ (0.5, 1], spectral radii of the Bathelike algorithm with different Ωtr are presented in Fig.18. Figure 19 gives their 24

tr

=4.13

tr

=5.00

tr

=8.00

tr

=16.0

tr

=16.0

tr

=25.0

tr

=100

tr

=200

tr

=300

CH-

=0.0

Bathe. =2-

0.06 0.14

0.05

Period elongation

Amplitude decay

0.12

0.04

0.03

0.02

0.10

0.08

0.06

0.04 0.01

0.02 0.00 0.00 0.00

0.25

0.50

0.75

1.00

1.25

1.50

0.00

0.25

0.50

0.75

1.00

Figure 17: Amplitude decays and period elongations of the Bathe-like algorithm with µ ∈ √ 2/2, 0.5).

(1 −

amplitude decays and period elongations. It is easy to find that there are the same results as the case of µ ∈ [−10, 0). Hence, further analysis of this case µ ∈ (0.5, 1] is also not given herein. Instead, the simple comparison of these four cases is given in the next subsection. 4.5. Comparisons of four cases Figure 20 presents comparisons of spectral radii for the Bathe-like algorithm when three values of Ωtr = 8, 16 and 50 are considered in four subintervals of µ. Obviously, for each of given Ωtr , the Bathe-like algorithm presents the identical spectral radii in four different subintervals of µ. However, the identical spectral radii don’t mean the same amplitude decays and period errors, as shown in Fig. 21. Figure 21 indicates that for any value of Ωtr , the µ ∈ [−10, 0) provides the same amplitude decays and period elongations as the µ ∈ (0.5, 1]. Similar results √ √ can be found for the µ ∈ (0, 1 − 2/2] and µ ∈ (1 − 2/2, 0.5). Importantly, √ √ for the same Ωtr , µ ∈ (0, 1 − 2/2] ∪ (1 − 2/2, 0.5) can provide less amplitude decays and period errors than µ ∈ [−10, 0) ∪ (0.5, 1]. Hence, it is generally recommended by authors to choose µ ∈ (0, 0.5) for given Ωtr . Based on the analysis of four cases. Some conclusions can be summarized

25

1.25

1.50

1.0

tr

0.8

tr

tr

)

tr

0.6

(

tr

tr

tr

0.4 tr

tr

=0.75 =1.00 =2.00 =4.00 =8.00 =16.0 =50.0 =100 =300

CH-

=0.0

0.2 Bathe. =2-

0.0 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

Figure 18: Spectral radii of the Bathe-like algorithm with µ ∈ (0.5, 1].

tr

=0.75

tr

=1.00

tr

=2.00

tr

=4.00

tr

=8.00

tr

=16.0

0.8

tr

=50.0

tr

=100

tr

=300

CH-

=0.0

Bathe. =2-

0.6

0.7 0.5

Period elongation

Amplitude decay

0.6

0.5

0.4

0.3

0.4

0.3

0.2

0.2 0.1 0.1

0.0

0.0

0.00

0.25

0.50

0.75

1.00

1.25

1.50

0.00

0.25

0.50

0.75

1.00

1.25

Figure 19: Amplitude decays and period elongations of the Bathe-like algorithm with µ ∈ (0.5, 1].

26

1.50

1.0 =-0.08926( = 0.09008( = 0.45050(

0.8

= 0.54097(

= 0.04424(

0.6

= 0.47685(

(

)

=-0.04423(

= 0.52118( =-0.01414(

0.4

= 0.01414( = 0.49283(

0.2

= 0.50697(

tr

tr

tr

tr

=8) =8)

tr

tr

tr

tr

tr

tr

=16)

=16) =16) =16)

tr

tr

=8)

=8)

=50)

=50) =50) =50)

0.0 10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

10

5

10

6

Figure 20: Comparisons of spectral radii for the Bathe-like algorithm with three different values of Ωtr in four cases.

0.016

0.20

= 0.09008( = 0.45050(

Amplitude decay

0.012

= 0.54097( =-0.04423(

0.010

0.008

= 0.04424( = 0.47685( = 0.52118(

0.006

=-0.01414( = 0.01414(

0.004

= 0.49283( = 0.50697(

tr

tr

tr

tr

tr

tr

=8)

tr

tr

0.16

= 0.09008( = 0.45050( = 0.54097(

=8) =16)

=16) =16) =16)

tr

tr

0.18

=8)

tr

tr

=-0.08926(

=8)

Period elongation

=-0.08926(

0.014

=50)

=50)

0.14

0.12

= 0.04424( = 0.47685(

0.10

= 0.52118( =-0.01414(

0.08 = 0.01414(

0.06

=50)

=-0.04423(

= 0.49283( = 0.50697(

=50)

tr

tr

tr

tr

=8) =8)

tr

tr

tr

tr

tr

tr

=16)

=16) =16) =16)

tr

tr

=8)

=8)

=50)

=50) =50) =50)

0.04

0.002

0.02 0.000 0.00 0.00

0.25

0.50

0.75

1.00

1.25

1.50

0.00

0.25

0.50

0.75

1.00

Figure 21: Comparisons of amplitude decays and period elongations for the Bathe-like algorithm with three different values of Ωtr in four cases.

27

1.25

1.50

as follows: (i) The Bathe-like algorithm always achieves the second-order accuracy and unconditional stability (the L-stability). (ii) The interval of parameter µ can be divided into four subintervals accompanied by four subfamilies of algorithms, each of which achieves the controllable algorithmic mode truncation factor Ωtr . √ (iii) In the case of µ = 1 − 2/2, numerical properties of the Bathe-like algorithm, including the spectral radius, amplitude decay and period elonga√ tion, are exactly identical to those of Bathe algorithm with γ = 2 − 2. In this case, the condition γθc2 = 1 is automatically satisfied to achieve the identical effective stiffness matrices inside two sub-steps. Hence, the √ present scheme with µ = 1 − 2/2 is considered as the optimal family of algorithms in this study. (iv) The larger values of µ ∈ [−10, 0) ∪ (0.5, 1] give very significant amplitude decays and period elongations which could be useful for some specific problems. Its common values within µ ∈ (0, 0.5) give the analogous results with √ the Bathe algorithm with γ = 2 − 2, which indicates that the Bathe-like algorithm with µ ∈ (0, 0.5) is suitable for solving most structural dynamic problems. (v) Among all algorithms shown herein, the Bathe-like algorithm with µ = √ √ 1 − 2/2 and the Bathe algorithm with γ = 2 − 2 achieve the smallest relative period errors. (vi) As for the choices of three parameters µ, θ and γ, the parameter µ can be determined by selecting the algorithmic mode truncation factor Ωtr based on Table 1 or Eq.(32). Then, the θ is considered to be 0.5, 1 or 1.5 (which means that the trapezoidal rule, backward Euler formula or extrapolation scheme is employed in the first sub-step, respectively). Lastly, the value of γ can be given by calculating γ = µ/θ. The simple estimate of the Ωtr will be given in Section 6.

28

5. Some algorithms from the Bathe-like algorithm In this section, the different cases of θ are considered to develop some novel time integration algorithms. Expect for the trapezoidal rule (i.e. µ = 0 and 1/2) and the backward Euler formula (i.e. µ = 1), the Bathe-like algorithms with any θ achieve the aforementioned numerical properties because numerical characteristics of the present scheme depend only on the algorithmic parameter µ instead of either θ or γ. An equivalent form in the second sub-step for the Bathe-like algorithm is given before continuing to discuss the algorithmic parameters. Substituting the generalized trapezoidal rule (2) in the first sub-step into (3) and rewriting the ˙ t+h as Ut+h and U ˙ t+h + d1 U ˙ t+γh + d0 U ˙ t) Ut+h = Ut + h(d2 U ˙ t+h = U ˙ t + h(d2 U ¨ t+h + d1 U ¨ t+γh + d0 U ¨ t) U

(34)

where d0 =

θ−1 , 2γθ − 2

d1 =

θ , 2 − 2γθ

d2 =

2γθ − 1 2γθ − 2

(35)

Obviously, when using (34) in the second sub-step, the algorithmic parameter γ and/or θ can be set to zero and the µ = γθ can be set to 1/2, not 1. On the other hand, when using (3) in the second sub-step, the γ cannot be set to zero and the µ = γθ can be taken to 1, not 1/2. Hence, there are not limitations introduced for the algorithmic parameters when the appropriate numerical scheme in the second sub-step are considered. Note that numerical schemes (34) have the same form as those of second sub-step of the ρ∞ -Bathe algorithm, and both the Bathe-like algorithm with θ = 1/2 and the ρ∞ -Bathe algorithm with ρ∞ = 0.0 correspond to the original Bathe algorithm. Hence, a new equivalent representation of the Bathe algorithm can be obtained, which is  γh  ˙ ˙ t+γh Ut + U 2   γh ˙t+ ¨t +U ¨ t+γh =U U 2

Ut+γh = Ut + ˙ t+γh U

29

(36)

and ˙ t+h + c1 U ˙ t+γh + c0 U ˙ t) Ut+h = Ut + h(c2 U ˙ t+h = U ˙ t + h(c2 U ¨ t+h + c1 U ¨ t+γh + c0 U ¨ t) U

(37)

where three coefficients c0 , c1 and c2 are determined by c0 = c1 =

1 , 4 − 2γ

c2 =

γ−1 γ−2

(38)

Note that equation (38) requires γ 6= 2 while the original study [12, 13, 15] can use γ = 2; similarly, the original works [12, 13, 15] cannot employ γ = 0 and 1 while they can be used herein in the equivalent form. Hence, it is said that the parameter γ of the Bathe algorithm is completely free when using the appropriate numerical scheme. 5.1. The trapezoidal rule: µ = 0 or 1/2 The µ = γθ = 0 indicates γ = 0 and/or θ = 0. Herein, the case of θ = 0 is firstly considered. In this case, the first sub-step uses the forward Euler formula, which is ˙t Ut+γh = Ut + γhU

(39a)

˙ t+γh = U ˙ t + γhU ¨t U

(39b)

and the numerical scheme in the second sub-step can be expressed as   ˙ t+h = 2Ut+h − 1 Ut+γh + 1 − 2 Ut hU γ γ   1 ¨ t+h = 2U ˙ t+h − U ˙ t+γh + 1 − 2 U ˙t hU γ γ

(40a) (40b)

Substituting (39) into (40) and reorganizing the results yield  h˙ ˙ t+h Ut + U 2   h ˙t+ ¨t +U ¨ t+h =U U 2

Ut+h = Ut +

(41a)

˙ t+h U

(41b)

Numerical schemes (41) can also be directly obtained by using (34). Obviously, the numerical scheme (41) is the common trapezoidal rule [5] in the interval [t, t+h]. This fact also explains why, in the first two subintervals of µ, numerical 30

characteristics of the Bathe-like algorithm tend to the trapezoidal rule as the µ changes to zero. In the case of γ = 0. the composite two sub-step algorithm actually reduces to the single-step algorithm. In the first sub-step, numerical solutions at time t + γh are directly equal to those at time t, i.e. ˙ t+γh = U ˙ t, U

Ut+γh = Ut ,

¨ t+γh = U ¨t U

(42)

In the second sub-step with (34), coefficients d0 , d1 and d2 read d0 =

1−θ , 2

d1 =

θ , 2

d2 =

1 2

(43)

Using Eq.(42) and (43), numerical scheme (34) corresponds to the trapezoidal rule (41). When considering the case of µ = γθ = 1/2. the numerical scheme (34) in the second sub-step should be used, which is simplified as   ˙ t+γh + (1 − θ)U ˙t Ut+h = Ut + h θU   ˙ t+h = U ˙ t + h θU ¨ t+γh + (1 − θ)U ¨t U

Combined with the numerical scheme in the first sub-step,   1 ˙ 1 ˙ Ut+γh = Ut + h (γ − )Ut + Ut+γh 2 2   1 ¨ 1¨ ˙ ˙ Ut+γh = Ut + h (γ − )Ut + Ut+γh 2 2

(44)

(45)

numerical scheme (44) can be rewritten as Ut+γh − Ut h γh ¨ ¨ ˙ t + Ut+γh − Ut h =U γh

Ut+h = Ut + ˙ t+h U

(46)

Note that the numerical scheme (44) in the second sub-step is explicit while (45) is implicit in the first sub-step. Hence, its computational cost is the same as the common traditional implicit algorithm, such as CH-α algorithm [10] and the trapezoidal rule. Consider the fact that the amplification matrix given in

31

Appendix A is not suitable for the case of µ = 1/2, the amplification matrix of numerical schemes (44)-(45) is presented in Appendix B. It is interesting to find that although numerical schemes (45) and (46) differ from the common trapezoidal rule in form, it achieves the same spectral properties and the same computational cost as the common trapezoidal rule, independent of the choice of γ. Hence, the present scheme with µ = 1/2 is also referred to the trapezoidal rule for brevity. 5.2. The Backward Euler formula: µ = 1 Considering the case of µ = 1, the numerical scheme in the second sub-step should use (3), not (34). The numerical schemes in the first and second sub-step are

and

  1 ˙ 1 ˙ Ut+γh = Ut + γh (1 − )U t + Ut+γh γ γ   1¨ 1 ¨ ˙ ˙ Ut+γh = Ut + γh (1 − )Ut + Ut+γh γ γ

(47)

˙ t+h = 1 (Ut+γh − Ut ) hU γ ¨ t+h = 1 (U ˙ t+γh − U ˙ t) hU γ

(48)

Similar to the Bathe-like algorithm with µ = 1/2, the computational cost of numerical schemes (47) and (48) is still the same as the single-step implicit algorithm. Particularly, for the case of γ = 1, the above scheme reduces to the first-order accurate backward Euler formula [33]. For general γ, although numerical schemes (47) and (48) are different from the backward Euler formula in form, the same spectral properties can be observed, independent of γ. Remark (a) The Bathe-like algorithm with three special values of µ (0,1/2 and 1) is presented in Fig.22, where some algorithms are given when the appropriate parameters γ and θ are given. In Fig.22, the optimal algorithm can be √ obtained by solving γθc2 = 1 (i.e. µ = 1 ± 2/2). In the case of µ = √ 1 − 2/2, the smallest period error and identical effective stiffness matrices can be ensured. 32

(b) For each given θ 6= 0, the Bathe-like algorithm can provide the same numerical characteristics as the original Bathe algorithm [12, 13]. As an example, not only the algorithm developed in [23] but also the present scheme with √ θ = 3/2 and γ = (2 − 2)/3 can be regarded as the alternative to the Bathe √ algorithm with γ = 2 − 2. In terms of the original Bathe with general γ := γb , the present scheme can provide the same numerical properties only if its γ is set to be γ = γb /(2θ). This is the reason why the present scheme is named as the Bathe-like algorithm. (c) For µ = 0 or 1/2, the Bathe-like algorithm corresponds to the trapezoidal rule. Figure 23 also illustrates this conclusion numerically. Besides, some first-order accurate integration algorithms, including the first-order accurate β1 /β2 -Bathe algorithm [21] and Newmark algorithm [2] with γ > 1/2, are compared with respects to the amplitude decay and period elongation. Among the algorithms presented herein, the backward Euler formula provides the largest numerical dissipations. It is necessary to point out that due to the use of the sub-step strategy, the first-order accurate β1 /β2 -Bathe algorithm presents the smallest period elongation errors.

5.3. The original Bathe algorithm: θ = 1/2 When the algorithmic parameter θ = 1/2 is considered, the present integration scheme reduces obviously to the original Bathe algorithm [12, 13]. Note that the original Bathe algorithm also achieves the controllable algorithmic mode truncation factor by µ = γθ = γ/2. Similar conclusions have been analyzed by directly considering γ in [15], but the cases of γ = 1, 2 and γ ≤ 0 have not been considered there. As analyzed in the above subsection, when γ = 0 or 1 is used, numerical properties of the Bathe algorithm tend to those of the trapezoidal rule while in the case of γ = 2, the Bathe algorithm tends to the backward Euler formula. These results also further explain the numerical phenomenon presented in [15] when the γ is close to 0, 1 or 2.

33

3.0

2.5 Backward Euler Formula (

2.0

c

=1)

=1

2

Bathe-like algorithm with

=3/2

1.5 Bathe-like algorithm with

1.0

=1

0.5

Bathe-like algorithm with

0.0

Trapezoidal rule

0.0

0.5

=1/2

=1/2

(the original Bathe algorithm)

1.0

1.5

2.0

2.5

Figure 22: The members of the Bathe-like algorithm when different algorithmic parameters are given.

5.4. The newly developed algorithm: θ = 1 If the case of θ = 1 is chosen, the backward Euler formula [33] is used in the first sub-step. The new composite sub-step algorithm, whose numerical dissipations are introduced in the first sub-step, is developed. The similar design has been used in [30] to give the generalized composite sub-step algorithm with controllable numerical dissipations. For the case of θ = 1, the coefficients in the second sub-step can be rewritten as c0 = −

2γ 2 − 2γ + 1 , 2γ 2 − γ

c1 =

1 , 2γ 2 − γ

c2 =

2γ − 2 2γ − 1

(49)

To ensure the identical effective stiffness matrices insider two sub-steps, the requirement γθc2 = 1 yields γ =1±

√

2 2

(50)

With Eq.(50), the newly developed algorithm reduces to the composite sub-step algorithm proposed in [23]. Note that the composite sub-step scheme developed

34

/

1

/

1

-Bathe alg.

2

=1.00 &

=2

2

1

= /2

-Bathe alg.

=0.39 &

2

Newmark alg.

1

=0.65 &

=2

2

1

= /2

/

1

-Bathe alg.

2

Newmark alg.

=0.55 &

1

=0.85 &

=2

2

1

= /2

Trapezoidal rule

1.0

0.5

0.8

0.4

Period elongation

Amplitude decay

Newmark alg.

=0.85 &

1

0.6

0.4

0.3

0.2

0.2 0.1

0.0 0.0 0.0

0.3

0.6

0.9

1.2

1.5

0.0

0.3

0.6

0.9

1.2

Figure 23: Amplitude decays and period elongations for the trapezoidal rule and some firstorder accurate algorithms.

in [23] achieves the same numerical properties as the Bathe algorithm with √ γ = 2 − 2. Hence, the newly developed algorithm (i.e. the Bathe-like algorithm with θ = 1) can be considered as the generalized version of the algorithm developed in [23]. Furthermore, it is interesting to find that if the γ = 1 is used, leading to µ = 1, the present scheme with θ = 1 will reduce to the first-order accurate backward Euler formula [33], as shown in Fig.22. 5.5. The newly developed algorithm: θ = 3/2 When the θ = 3/2 is considered, another newly developed algorithm is given. In this case, none of existing composite sub-step algorithms included as the special cases has been proposed. Particularly, if the requirement γθc2 = 1 is achieved, the algorithmic parameter γ will be given as √ 2± 2 γ= 3

(51)

√ to ensure the identical effective stiffness matrices. With γ = (2 − 2)/3, the √ corresponding value of µ is 1 − 2/2 so that the same desired numerical prop√ erties as the Bathe algorithm with γ = 2 − 2 can be obtained due to the same µ. 35

1.5

In the case of γ = 1/θ = 2/3, the proposed algorithm also corresponds to the first-order accurate backward Euler formula, as indicated in Fig.22. Remark (a) In the Bathe-like algorithm, when the θ is given, the corresponding new composite sub-step algorithm can be developed. Four values of θ = 0, 1/2, 1 and 3/2 have been considered in detail. For each of given θ 6= 0, the desired numerical properties, such as the smallest period error, can be obtained √ when considering µ = γθ = 1 − 2/2, in which case the requirement γθc2 = 1 is also satisfied to ensure the identical effective stiffness matrices inside two sub-steps. (b) For three values of θ = 1, 1/2 and 3/2, three algorithms with the corre√ √ √ sponding γ = 1 + 2/2, 2 + 2 and (2 + 2)/3 also achieves the same numerical properties because of the same µ. (c) When the θ 6= 0 is given, the present scheme can provide a composite sub-step algorithm with the same numerical characteristics as the original Bathe algorithm [12, 13]. Hence, the Bathe-like algorithm with given θ can be considered as the alternative to the original Bathe algorithm with any γ.

6. Implementation for nonlinear problems After the numerical analysis of the Bathe-like algorithm based on linear systems, it is important to give the procedure of solving nonlinear problems. In fact, equations (7)-(12) present the procedure of solving linear elastic problems. The pseudo code for the application of the Bathe-like algorithm to nonlinear problems (52) is given in Algorithm 1. ¨ ˙ MU(t) + R(U(t), U(t)) = F(t)

(52)

˙ where R(U(t), U(t)) is the vector of nodal forces depending on velocity and displacement vectors. Algorithm 1 Bathe-like algorithm for nonlinear problems 1:

Set: parameter θ, µ and calculate the corresponding γ using (17). 36

2:

Set: the constant coefficients ci , i = 0, 1, 2 using (4).

3:

¨ 0 from U ˙ 0 and U0 . Compute: the initial acceleration U

4:

for t = t0 to t = tf do

5:

//———————–In the first sub-step——————————–

6:

˙t+ Predict: Ut+γh = Ut + γhU

7:

for i = 1 to max-iter do

8: 9: 10: 11: 12:

(1)

(γh)2 ¨ 2 Ut .

(i) 1 1−θ ˙ µh (Ut+γh − Ut ) − θ Ut . ¨ (i) = 1 (U ˙ (i) − U ˙ t ) − 1−θ U ¨ t. Compute: U t+γh t+γh µh θ (i) (i) (i) (i) 1 1 ˙ e Compute: R(U t+γh , Ut+γh ) and K = µ2 h2 M + µh Ct+γh + Kt+γh . (i) e ˙ (i) , U(i) ) − MU ¨ (i) . Solve: ∆U(i) using K∆U = Ft+γh − R(U t+γh t+γh t+γh

˙ (i) = Compute: U t+γh

if convergence then

13:

(i) ˙ t+γh , U ¨ t+γh using Steps 8 and 9. Compute: Ut+γh = Ut+γh , U

14:

Exit iteration loop.

15:

end if

16:

Update: Ut+γh ← Ut+γh + ∆U(i) .

(i+1)

(i)

17:

end for

18:

//———————–In the second sub-step——————————– h i (1) ˙ t−U ˙ t+γh ] . Predict: Ut+h = γ13 (γ 3 − 3γ + 2)Ut + (3γ − 2)Ut+γh + (γ − 1)γh[(γ − 1)U

19: 20:

for i = 1 to max-iter do

21:

˙ (i) = 1 (c2 U(i) + c1 Ut+γh + c0 Ut ). Compute: U t+h t+h h

22:

˙ (i) + c1 U ˙ t+γh + c0 U ˙ t ). ¨ (i) = 1 (c2 U Compute: U t+h t+h h

23: 24: 25:

˙ (i) , U(i) ) and K e = Compute: R(U t+h t+h

c22 h2 M

(i) e Solve: ∆U(i) using K∆U = Ft+h −

(i) c2 (i) h Ct+h + Kt+h . ˙ (i) , U(i) ) − MU ¨ (i) . R(U t+h t+h t+h

+

if convergence then

26:

(i) ˙ t+h , U ¨ t+h using Steps 21 and 22. Compute: Ut+h = Ut+h and U

27:

Exit iteration loop.

28:

end if

29:

Update: Ut+h ← Ut+h + ∆U(i) .

(i+1)

(i)

30:

end for

31:

// Next integration step.

32:

end for 37

(1)

Note that the initial predicted value of Ut+h (i.e. Ut+h in Step 19) is determined by constructing the cubic Hermite interpolation based on displacements and velocities at time t and t + γh. Of course, other interpolation techniques, such as the Taylor series truncation [36], can be used to obtain the similar ini(1)

(i)

(i)

tial guess Ut+h . In Step 10, matrices Ct+γh and Kt+γh are determined by, respectively (i) Ct+γh

(i)

Kt+γh

˙ ∂R(U(t), U(t)) = ˙ (i) ∂ U(t) ˙ ˙ (i) U U(t)=Ut+γh ,U(t)= t+γh ˙ ∂R(U(t), U(t)) = (i) ∂U(t) ˙ ˙ (i) U(t)=U ,U(t)= U t+γh

(i)

(53a) (53b)

t+γh

(i)

Similar computations for matrices Ct+h and Kt+h can be obtained in Step 23. In Step 1, it is essential to appropriately select two algorithmic parameters θ and µ. Since the parameter θ has no influence on numerical properties of linear systems, any of three choices θ = 0.5, 1.0 and 1.5 analyzed in this study can be employed. For each given θ, the Bathe-like algorithm corresponds to a singleparameter (µ) integration scheme. The µ is associated with the algorithmic mode truncation factor Ωtr . Assume that (i) there are n > 1 modes in the system being analyzed, which are sorted in ascending order as ω1 < ω2 < · · · < ωj < · · · < ωn−1 < ωn and (ii) the modes more than and equal to ωj are unexpected and should be damped out effectively. The initial estimate of Ωtr is easily given as the product of ωj and time step h. If its value is less than 0.74, the finial Ωtr is recommended as 0.74 or the slightly larger value. When the initial estimate exceeds 0.74, the value slightly larger than this initial estimate is suggested in order to decrease the influence of dissipations on the modes close to ωj . Once the Ωtr is given, the corresponding µ can be obtained based on Table 1 or Eq.(32). On the other hand, if users want to determine which modes to be damped out by adjusting time step h or users feel confused about √ selecting appropriate µ, the default value of µ = 1− 2/2 is recommended, since the Bathe-like algorithm achieves optimal numerical properties in this case. It should be pointed out that if users replace the smallest unexpected mode ωj by 38

a larger mode ωk , k > j, the final Ωtr can be directly given by its initial estimate. In this case, there are still some dissipations to eliminate modes ω ∈ [ωj , ωk ]. 7. Numerical examples In this section, some numerical examples are presented to show the advantages of the Bathe-like algorithm. If not explicitly stated, the algorithmic parameter θ is set to one by default. Herein, numerical responses from the CH-α √ algorithm with ρ∞ = 0 and the original Bathe algorithm with γ = 2 − 2 are presented for comparison. 7.1. The free-vibration of the three-DOF mass-spring system In this numerical example, the three-DOF system presented in Fig.24 is considered to show the superiority of the Bathe-like algorithm over the controllably dissipative algorithm with respect to the controllable algorithmic mode truncation factor. The governing equation can be expressed as

m1

m2

m3

Figure 24: The three-DOF mass spring model.

 m  1   0  0

0 m2 0

   u ¨ k + k2   1  1    0  u ¨  +  −k2   2  m3 u ¨3 0 0

−k2 k2 + k3 −k3

    u 0   1       −k3  u2  = 0     k3 u3 0 0

(54)

where all masses are assumed to be unity and three stiffness coefficients k1 , k2 and k3 are set to be 105 , 102 and 1, respectively. Three natural frequency of this system are ω1 = 0.9950, ω2 = 10.0454 and ω3 = 316.3860. The initial velocities are assumed to be zero while the initial displacements read h iT = φ1 + 0.005φ2 + 0.005φ3 u1 u2 u3 t=0

39

(55)

where the φi , i = 1, 2, 3 are the mode shape corresponding to the natural frequency ωi , i = 1, 2, 3. The time integration algorithm with controllable numerical dissipations is firstly used to solve this model in order to show the fact that the parameter ρ∞ only controls how fast the spurious modes to be damped out and cannot determine which modes to be eliminated. Herein, the ρ∞ -Bathe algorithm [22] with √ optimal γ = (2 − 2 + 2ρ∞ )/(1 − ρ∞ ) is selected as a representative of the controllably dissipative algorithms. Figure 25 plots numerical displacements of this -5

10 2

= 0.9

= 0.6

= 0.3

= 0.0

u1

1

0

-1

-2 0.02

= 0.9

= 0.6

= 0.3

= 0.0

= 0.9

= 0.6

= 0.3

= 0.0

u2

0.01

0.00

-0.01

-0.02

1.0

u3

0.5

0.0

-0.5

-1.0 0

5

10

15

20

25

time

Figure 25: Numerical displacements of three-DOF system from the ρ∞ -Bathe algorithm using h = 0.02s.

system using the ρ∞ -Bathe algorithm with h = 0.02s. Note that the numerical displacement u3 only presents the first mode while the u2 reflects the super40

30

position of the first two modes. The third mode ω3 in the u1 is successfully eliminated by introducing numerical dissipations (i.e. decreasing ρ∞ ). Obviously, for any values of ρ∞ 6= 1, the ρ∞ -Bathe algorithm only eliminates the third mode and hence retains the first two modes. The smaller ρ∞ is, the faster the spurious third mode is filtered out. If users require to filter out the last two modes (ω2 and ω3 ) in this system, the ρ∞ -Bathe algorithm with h = 0.02s cannot work well in this case. The controllably dissipative algorithm has to change its integration step size to determine the unexpected modes to be damped out. Herein, the ρ∞ -Bathe algorithm with larger step h = 0.2s is used to filter out the last two modes. Figure 26 presents numerical displacements when using the 10

-5

2 = 0.9

= 0.6

= 0.3

= 0.0

= 0.9

= 0.6

= 0.3

= 0.0

= 0.9

= 0.6

= 0.3

= 0.0

u1

1

0

-1

-2 0.02

u2

0.01

0.00

-0.01

-0.02

1.0

u3

0.5

0.0

-0.5

-1.0 0

5

10

15

20

25

30

time

Figure 26: Numerical displacements of three-DOF system from the ρ∞ -Bathe algorithm using h = 0.2s.

larger time step size h = 0.2s. As the ρ∞ decreases from 0.9 into 0.0, numeri41

cal displacements u1 and u2 only reflect the first mode and hence the last two modes have been successfully eliminated. It is necessary to point out that the ρ∞ -Bathe algorithm with other γ would provide similar numerical results. Next, the Bathe-like algorithm with µ ∈ (0.5, 1] is employed to solve this model in order to illustrate the advantage of the controllable algorithmic mode truncation factor. First, the only third mode ω3 is considered as the spurious mode and hence should be effectively damped out. According to the Table 1, the value of µ should belong to appropriately µ ∈ [0.50118, 0.54] to filter out the ω3 . Herein, the µ = 0.53 is selected. Numerical displacements obtained from the 10

-5

2 = 0.53

u1

1

0

-1

-2 0.02 = 0.53

u2

0.01

0.00

-0.01

-0.02 1.0 = 0.53

u3

0.5

0.0

-0.5

-1.0 0

5

10

15

20

25

30

time

Figure 27: Numerical displacements of three-DOF system from the Bathe-like algorithm with µ = 0.53 using h = 0.02s.

Bathe-like algorithm with µ = 0.53 are plotted in Fig.27, where the unexpected third mode has been filtered out. In fact, the unexpected third mode can be 42

eliminated faster by selecting larger µ, such as µ = 0.56. Assume that the second mode ω2 is also unexpected and should be filtered out. Unlike the controllably dissipative algorithm, the Bathe-like algorithm can only adjust its algorithmic parameter to achieve the elimination of the last two modes. Theoretically, the initial Ωtr is estimated as ω2 · h ≈ 0.2, which is larger than 0.74. Hence, the finial Ωtr can be given as 0.74 or the sightly larger value. The Ωtr ≈ 0.93 is chosen in this example and the corresponding µ is given to be 0.85 based on Eq.(32). Figure 28 presents numerical displacements in this case. 10

-5

2 = 0.85

u1

1

0

-1

-2 0.02 = 0.85

u2

0.01

0.00

-0.01

-0.02 1.0 = 0.85

u3

0.5

0.0

-0.5

-1.0 0

5

10

15

20

25

30

time

Figure 28: Numerical displacements of three-DOF system from the Bathe-like algorithm with µ = 0.85 using h = 0.02s.

The last two modes are successfully damped out by the Bathe-like algorithm by adjusting µ from 0.53 to 0.85. Note that both the ρ∞ -Bathe algorithm with ρ∞ = 0.0 and the Bathe-like algorithm with µ = 0.85 can filter out the last two 43

modes, but the former presents more numerical errors in u3 due to using larger integration step, as shown in Fig.29. In fact, the larger µ can faster eliminate the second mode, but it could, in turn, lead to the lower solution accuracy in the long-time simulation. Hence, the choice of µ is so important that more than one numerical experiments could be required to obtain the appropriate µ. 0.05

Absolute errors of

u

3

-Bathe 0.04

Bathe-like

=0.0 & =0.85 &

h=0.2s h=0.02s

0.03

0.02

0.01

0.00 0

5

10

15

20

25

30

time

Figure 29: Absolute errors of u3 of the Bathe-like algorithm and the ρ∞ -Bathe algorithm.

7.2. A clamped-free bar excited by end load A clamped-free bar [21, 22] subjected to the stepwise end load F (t), as shown in Fig.30, is considered. The physical model is discretized into N = 100 elements of equal length l = L/N , and the lumped mass matrix is employed in this numerical model. The lumped mass and stiffness matrices are x, u F(x) L

Figure 30: A clamped-free bar excited by the end load: E = 3 × 107 , A = 1, L = 200 and m = 7.3 × 10−4 .

44



2

 −1    EA   K= l      

−1 2

−1

−1

2

..

.

..

..

.

0

.

0

−1

−1 2 −1



           −1  1

N ×N

 2      ml   M= 2      

2



0 2 ..

0

. 2 1 (56)

The applied end load is assumed to be F (t) = 1 × 104 . The time step size h used by all algorithms provided herein satisfies h = (9.88 × 10−7 ) × CFL

(57)

It has been shown in [21] that the first-order accurate dissipative algorithms are more advantageous than the second-order methods when solving this kind of wave propagation problems. The β1 /β2 -Bathe algorithm [21] with β2 = 2β1 (the first-order accuracy) is firstly used to solve this model. Figure 31 plots numerical displacements and velocities at node 50 using the β1 /β2 -Bathe algorithm with β2 = 2β1 and β1 = 0.85. Note that the β1 = 0.85 is used herein while the β1 = 0.39 is employed in [21] since the total number of linear finite elements used in this physical model is different. Figure 31 indicates that the firstorder β1 /β2 -Bathe algorithm can accurately predict displacement and velocity responses. It is well-known that the classical Newmark algorithm only achieves the firstorder accuracy and hence introduces numerical dissipations when the parameter γ > 1/2. Herein, the first-order accurate Newmark algorithm is also employed to solve this problem. The corresponding numerical results are presented in Fig. 32. Similar to the β1 /β2 -Bathe algorithm, the first-order Newmark algorithm also provides pretty accurate numerical solutions. As pointed out previously, the Bathe-like algorithm will reduce to the firstorder backward Euler formula when the µ is set to be one. Numerical results from the Bathe-like algorithm with µ = 1 are plotted in Fig.33, where the Bathelike algorithm with µ = 1 also provides very accurate numerical responses. 45

            

N ×N

0.07

Displacement

0.06

0.05

0.04

Exact =0.85 CFL=1

1

0.03

=0.85 CFL=4

1

0.02

0.01

0.00 80 60

Velocity

40 20 0 -20 Exact -40 =0.85 CFL=1

1

-60

=0.85 CFL=4

1

-80 0.000

0.004

0.008

0.012

Time

Figure 31: Numerical displacements and velocities at node 50 using the β1 /β2 -Bathe algorithm with β2 = 2β1 .

46

0.016

0.07

Displacement

0.06

0.05

0.04

Exact =1.3 &

= /2 CFL=1

=1.3 &

= /2 CFL=4

0.03

0.02

0.01

0.00 80 60

Velocity

40 20 0 -20 -40 -60

Exact =1.3 &

= /2 CFL=1

=1.3 &

= /2 CFL=4

-80 0.000

0.004

0.008

0.012

Time

Figure 32: Numerical displacements and velocities at node 50 using the first-order Newmark algorithm.

47

0.016

Compared with numerical results from three first-order accurate algorithms, although they all provide the acceptable numerical solutions, the computational cost of the β1 /β2 -Bathe algorithm is about two times that of the Newamrk and Bathe-like algorithm with µ = 1. Hence, when the first-order accuracy is required in some solutions, the first-order Newmark algorithm and the Bathe-like algorithm with µ = 1 could be a better choice due to involving less computational cost. 0.07

Displacement

0.06

0.05

0.04

Exact =1 CFL=1

0.03

=1 CFL=4

0.02

0.01

0.00 80 60

Velocity

40 20 0 -20 -40

Exact =1 CFL=1

-60 =1 CFL=4 -80 0.000

0.004

0.008

0.012

Time

Figure 33: Numerical displacements and velocities at node 50 using the Bathe-like algorithm with µ = 1.

Finally, three second-order Bathe-like algorithm (i.e. µ = −5, 0.6 and 0.95) are used to resolve this model, see Fig.34. Since the Bathe-like algorithm with µ = 0.5 is spectrally equivalent to the trapezoidal rule, very significant numerical

48

0.016

oscillations can be found for the Bathe-like algorithm with µ = 0.6 (near to 0.5). On the other hand, although the µ = 0.95 is close to one (in which case the Backward Euler formula is obtained), spurious numerical oscillations are also found easily. It is necessary to point out that although numerical responses can be further improved by adjusting µ, very accurate numerical solutions are still from the case of µ = 1 (i.e. the first-order accurate Bathe-like algorithm). 0.070

90

Exact = -5

60 0.068

CFL=1

= 0.6 CFL=1

Velocity

Displacement

= 0.95CFL=1 30 0.066

0

0.064

-30 Exact 0.062

= -5

CFL=1

-60

= 0.6 CFL=1 = 0.95CFL=1 0.060

-90 0.0132

0.0136

0.0140

0.0144

0.013

Time

0.014

0.015

Time

Figure 34: Numerical displacements and velocities at node 50 using the Bathe-like algorithm with µ = −5, 0.6 and 0.95.

7.3. A simple nonlinear pendulum The simple nonlinear pendulum problem, shown in Fig.35, has been used

Figure 35: A simple nonlinear pendulum.

49

to test numerical performances of various integration methods in the previous studies [23, 30, 37]. Its governing differential equation can be given as ϕ(t) ¨ +

g sin(ϕ(t)) = 0, L

ϕ(0) = ϕ0 ,

ϕ(0) ˙ = ϕ˙ 0

(58)

where the g and L represent the acceleration of gravity and length of pendulum, respectively. The exact solution to the nonlinear pendulum can be obtained by using the elliptic integral of the first kind [30] and can be naturally used as the good reference. In this paper, the g = 1 and L = 1 are employed and the initial conditions θ0 = 0 and θ˙0 = 1.999999238456499 are chosen to describe the highly nonlinear pendulum. With the given data, this pendulum will oscillate between two peak points (θ = ±π) instead of making a complete rotation since the initial total energy is mildly less than the minimum total energy required to make complete turns. Figure 36 presents numerical displacements of the pendulum in ϕ direction by using three Bathe-like algorithms when the same integration step h = 0.002s is used. Firstly, the original Bathe algorithm with different γ, which is actually the Bathe-like algorithm with θ = 0.5, is employed to solve this model. Figure 36a reveals that there are very significant period errors for the original Bathe algorithm with different γ. Particularly, the Bathe algorithm with γ = 1.2(µ = 0.6) cannot predict the reasonable numerical solutions since the complete rotation is found. Another two Bathe-like algorithms with θ = 1.0 and 1.5 can provide reasonable numerical results although there are still noticeable period errors. In contrast to the original Bathe algorithm (i.e. the Bathe-like algorithm with θ = 0.5), two Bathe-like algorithms with θ = 1.0 and 1.5 can provide more accurate solutions by adjusting µ. Figure 37 plots numerical displacements of this pendulum by using three Bathe-like algorithms in two cases of µ = 0.3 and 0.35, where the Bathe-like algorithms with θ = 1.0 and 1.5 can accurately predicts numerical solutions. Recall that for given µ, the Bathe-like algorithms with different θ all achieve the same spectral properties. Hence, it is revealed in Fig.37 that the newly proposed algorithms (i.e. the Bathe-like algorithm with θ = 1.0 and 1.5) are superior to the original Bathe algorithm when 50

direction Displacement in direction

9.42

a

( ) The original Bathe algorithm (The Bathe-like algorithm with

Reference =0.5)

6.28

=0.1(

=0.05)

=0.3(

=0.15)

=0.5(

=0.25)

=0.8(

=0.40)

=1.2(

=0.60)

3.14

0.00

-3.14

b) The Bathe-like algorithm with

(

=1.0

3.14

direction

Displacement in

Reference 0.00

-3.14

=0.05 &

=1

=0.15 &

=1

=0.25 &

=1

=0.40 &

=1

=0.60 &

=1

c

( ) The Bathe-like algorithm with

=1.5

3.14

Displacement in

Reference 0.00

-3.14

0.00

=0.05 &

=1.5

=0.15 &

=1.5

=0.25 &

=1.5

=0.40 &

=1.5

=0.60 &

=1.5

16.86

33.72

50.58

67.44

Time

Figure 36: Numerical displacements of the nonlinear pendulum in ϕ direction by using three Bathe-like algorithms with different µ.

51

solving this nonlinear pendulum. This numerical example also indicates that it is not always a good choice to use the trapezoidal rule (i.e. θ = 0.5) within the first sub-step. The Bathe-like algorithm proposed in this paper includes all possible numerical schemes used in the first sub-step.

Displacement in

direction

Reference =0.3 & =0.35 &

=0.5

=0.3 &

=0.5

=1.0

=0.35 &

=0.3 &

=1.0

=0.35 &

=1.5 =1.5

3.14

0.00

-3.14 0.00

16.86

33.72

50.58

67.44

Time

Figure 37: Numerical displacements of the nonlinear pendulum in ϕ direction by using three Bathe-like algorithms when u = 0.3 and 0.35 are used.

7.4. The cantilever beam A cantilever beam [13, 38] in the case of plane strain shown in Fig. 38 is considered to shown the ability of reliable stability. This cantilever beam 500 f N /m2

0.001 m 0.4 m f E=70∗109 N /m2 ρ =2700 kg/m3 μ =0.33 plane strain

1

time [s] 0

0.02

0.04

Figure 38: Cantilever beam modeled using 8-node two-dimensional elements.

52

is modeled with 400 eight-node two-dimensional elements, and subjected to pressure loading. The beam is supported to undergo large displacements [13, 38]. √ Herein, two primary facts are recalled again that (i) in the case of µ = 1 − 2/2, the Bathe-like algorithm with θ = 1/2 is essentially the Bathe algorithm with √ γ = 2− 2 [12, 13]; and (ii) the trapezoidal rule fails to provide stable numerical

d is p la c e m e n t [m ]

responses of this cantilever beam according to the results provided in [13, 38].

-

-

-

-

0 .2 0 .1 0 .1 0 .0 0 .0 0 .0 0 .1 0 .1 0 .2

0

θ= 0 .5 θ= 1 .0 θ= 1 .5 5 0 5 0 5 0 5 0

8

θ= 0 .5 θ= 1 .0 θ= 1 .5

v e lo c ity [m /s ]

6 4 2 0 -2 -4 -6 -8

a c c e le ra tio n [m /s 2]

3 0 0

θ= 0 .5 θ= 1 .0 θ= 1 .5

2 0 0 1 0 0 0 -1 0 0 -2 0 0 -3 0 0

0 .0

0 .5

1 .0

1 .5

2 .0

tim e [s ]

Figure 39: Numerical responses of the cantilever beam from the Bathe-like algorithm with √ µ = 1 − 2/2 using h = 0.002s.

Figure 39 presents numerical vertical displacements, velocities and accelera√ tions at the tip of this beam from the Bathe-like algorithm with µ = 1 − 2/2 using h = 0.002s. The results show that three curves of numerical responses using θ = 0.5, 1.0 and 1.5 are the exactly same due to the same spectral properties. Hence, the Bathe-like algorithm with different θ can still predict the same √ accurate numerical solutions as the Bathe algorithm with γ = 2− 2. The same conclusions can be obtained for the Bathe-like algorithm with other µ, as shown

53

in Fig.s 40 and 41. They indicate that the Bathe-like algorithm with different 0 .2 0

θ= 0 .5 θ= 1 .0 θ= 1 .5

d is p la c e m e n t [m ]

0 .1 5 0 .1 0 0 .0 5 0 .0 0 - 0 .0 5 - 0 .1 0 - 0 .1 5 - 0 .2 0 8

θ= 0 .5 θ= 1 .0 θ= 1 .5

v e lo c ity [m /s ]

6 4 2 0 -2 -4 -6 -8

θ= 0 .5 θ= 1 .0 θ= 1 .5

a c c e le ra tio n [m /s 2]

3 0 0 2 0 0 1 0 0 0 -1 0 0 -2 0 0 -3 0 0 0 .0

0 .5

1 .0

1 .5

2 .0

tim e [s ]

Figure 40: Numerical responses of the cantilever beam from the Bathe-like algorithm with µ = −0.001 using h = 0.002s.

values of µ can also provide stable and accurate numerical results. These figures illustrate that numerical characteristics of the Bathe-like algorithm depend only on the value of µ, not either θ or γ. Hence, for any given θ 6= 0, the Bathe-like algorithm can be considered as an alternative to the original Bathe algorithm. When increasing integration step size h into 0.01s, the Bathe-like algorithm with different values of µ and θ = 1.5 also gives the reliable numerical solutions √ in displacement and velocity, as presented in Fig.42. the µ = 1 − 2/2 obviously eliminates the high-frequency modes in acceleration, since it possesses the smaller algorithmic mode truncation factor than other two cases. In fact, Figure 42 also reveals that an alternative way to damp out the high-frequency modes is to adjust the parameter µ, not to increase time step h.

54

d is p la c e m e n t [m ] -

0 .2 0 .1 0 .1 0 .0 0 .0 0 .0 0 .1 0 .1 0 .2

0

θ= 0 .5 θ= 1 .0 θ= 1 .5 5 0 5 0 5 0 5 0

8

θ= 0 .5 θ= 1 .0 θ= 1 .5 6

v e lo c ity [m /s ]

4 2 0 -2 -4 -6 -8

a c c e le ra tio n [m /s 2]

3 0 0

θ= 0 .5 θ= 1 .0 θ= 1 .5

2 0 0 1 0 0 0 -1 0 0 -2 0 0 -3 0 0 0 .0

0 .5

1 .0

1 .5

2 .0

tim e [s ]

Figure 41: Numerical responses of the cantilever beam from the Bathe-like algorithm with

d is p la c e m e n t [m ]

µ = 0.51 using h = 0.002s.

-

-

-

-

0 .2 0 .2 0 .1 0 .1 0 .0 0 .0 0 .0 0 .1 0 .1 0 .2

5 0

µ = 1 - √2 / 2

µ= 0 .5 1

µ= -0 .0 0 1

µ = 1 - √2 / 2

µ= 0 .5 1

µ= -0 .0 0 1

µ = 1 - √2 / 2

µ= 0 .5 1

0 5 0 5 0 5 0 8 6

v e lo c ity [m /s ]

µ= -0 .0 0 1 5

4 2 0 -2 -4 -6 -8 4 0 0

a c c e le ra tio n [m /s 2]

3 0 0 2 0 0 1 0 0 0 -1 0 0 -2 0 0 -3 0 0 0 .0

0 .5

1 .0

1 .5

2 .0

tim e [s ]

Figure 42: Numerical responses of the cantilever beam from the Bathe-like algorithm with different µ and θ = 1.5.

55

8. Conclusion In this paper, a class of Bathe-like algorithm is proposed and analyzed comprehensively by reconstructing the Bathe algorithm. The Bathe-like algorithm uses the generalized trapezoidal rule and three-points backward differential formula in the first and second sub-step, respectively. The second-order accuracy and unconditional stability are still obtained in its final form. A new parameter, named as the algorithmic mode truncation factor Ωtr , is firstly introduced to better describe the numerical property of the Bathe-like algorithm. The Ωtr can characterize the capability of eliminating spurious modes. Unlike the traditional dissipative parameter ρ∞ , the Ωtr can give an estimate on which modes to be damped out while the ρ∞ only can control how fast the spurious modes to be filtered out. The present Bathe-like algorithm achieves the controllable algorithmic mode truncation factor Ωtr by simply adjusting µ based on Eq.(32), which can help users estimate which modes to be eliminated when solving large multi-degree-of-freedom systems. The spectral analysis of the Bathe-like algorithm shows its controllable ability with respect to the algorithmic mode truncation factor, and reveals that (a) when µ = 1 −

√

2/2, the Bathe-like algorithm can achieve (i) the identical

effective stiffness matrices inside two sub-steps for any values of γ and θ; (ii) the exactly same numerical characteristics as the Bathe algorithm with √ γ = 2 − 2, but θ = 1/2 is not required in the Bathe-like algorithm. Hence, in this case, the Bathe-like algorithm with given θ 6= 0 can be regarded as √ an alternative to the Bathe algorithm with γ = 2− 2; and (iii) the smallest period elongation errors among all possible values of µ. (b) when µ = 0 or 1/2, the Bathe-like algorithm reduces to the non-dissipative trapezoidal rule. In terms of the original Bathe algorithm, this case also holds when its γ is set to be zero or unity. (c) when µ = 1, the Bathe-like algorithm corresponds to the first-order accurate backward Euler formula. For the original Bathe algorithm, the Bathe algorithm also reduces to the backward Euler formula in the case of γ = 2. 56

(d) the Bathe-like algorithm with any given non-zero θ can be considered as the alternative to the original Bathe algorithm since the exactly same numerical properties can be achieved. Importantly, these novel algorithms could provide more accurate numerical solutions than the original Bathe algorithm when some nonlinear problems are solved. In the end, four numerical experiments show the superiority of the Bathelike algorithm over other algorithms, such as the CH-α, trapezoidal rule and three Bathe-type algorithms. Appendix A: the amplification matrix of the Bathe-like algorithm Note that in the case of µ = 1/2, the corresponding amplification matrix is presented in Appendix B. Generally, it is very difficult for composite substep algorithms to directly derivate the explicit expression of the amplification matrix A, which can be divided into the following steps: First, applying the generalized trapezoidal rule in the first sub-step on the SDOF system yields



ut+γh





ut



         hu˙ t+γh  = A11  hu˙ t      h2 u ¨t+γh h2 u ¨t

(A1)

where the amplification matrix A11 in the first sub-step can be written explicitly as A11

 2γθξΩ + 1 1   =  −γθΩ2 β1  −Ω2

γ + 2θγχξΩ

θγχ

1 − θγχΩ2

χ

−γΩ2 − 2ξΩ

−χ(γθΩ2 + 2ξΩ)

    

(A2)

where β1 = γ 2 θ2 Ω2 + 2γθξΩ + 1 and χ = (1 − θ)γ. ω = Ω/h and ξ are the natural frequency and viscous damping ratio of the SDOF system, respectively. Then, the similar calculation can be done in the second sub-step, which yields



ut+h





ut+γh





ut



             hu˙ t+h  = A22  hu˙ t+γh  + A21  hu˙ t        h2 u ¨t+h h2 u ¨t+γh h2 u ¨t 57

(A3)

where two iteration matrices A22 and A21 can be expressed as   −c2 − 2ξΩ −1 0  c1    A22 =  Ω2 −c2 0 β2   c2 Ω2 Ω2 + 2c2 ξΩ 0 A21 =

c0 A22 c1

(A4)

(A5)

where β2 = Ω2 + 2c2 ξΩ + c22 . c0 , c1 and c2 are determined by Eq.s (4). In the end, Substituting Eq.(A1) into (A3) yields     ut ut+h          hu˙ t+h  = A  hu˙ t      h2 u ¨t h2 u ¨t+h

(A6)

where the amplification matrix A is presented as follows A = A22 · A11 + A21

(A7)

Hence, substituting Eq.s (A2), (A4) and (A5) into Eq.(A7) gives the amplification matrix A.

Appendix B: the amplification matrix of the Bathe-like algorithm with µ = 1/2 In the case of µ = 1/2, the corresponding amplification matrix can be directly given without any difficulties, which is   a11 (1 − 2θγ)Ω2 + 4(1 − θ)ξΩ + 4 2(2γ − 1)θ  1    A= −4θΩ2 (1 − 4γθ)Ω2 − 4(2θ − 1)ξΩ + 4  a23 β2   a31 a32 a33

58

(B1)

where a11 = (1 − 2θ)Ω2 + 4ξΩ + 4 a23 = (1 − 2γθ)Ω2 − 4(2θ − 1)ξΩ + 4 − 4θ a31 = (2θ − 1)Ω4 + 4(2θ − 1)ξΩ3 − 4Ω2 a32 = (2γθ − 1)Ω4 + ((8γ + 4)θ − 6)ξΩ3 + 2(−2 + 4(2θ − 1)ξ 2 )Ω2 − 8ξΩ a33 = (4γθ − 2)ξΩ3 + 8(2θγ − 1)ξ 2 Ω2 + 2θ(1 − 2γ)Ω2 + 8(θ − 1)ξΩ

(B2)

and β2 = Ω2 + 4ξΩ + 4

(B3)

Acknowledgment This work is supported by the National Natural Science Foundation of China (Grant No. 11372084). This support is gratefully acknowledged. The helpful and constructive comments by the referees have led to the improvements of this paper; the authors gratefully acknowledge this assistance.

References [1] J. C. Houbolt, A Recurrence Matrix Solution for the Dynamic Response of Elastic Aircraft, Journal of the Aeronautical Sciences 17 (9) (1950) 540– 550. [2] N. M. Newmark, A method of computation for structural dynamics, Journal of Engineering Mechanic Division 85 (3) (1959) 67–94. [3] Feappv,

a

finite

element

analysis

program:

Personal

version,

http://projects.ce.berkeley.edu/feap/feappv/. [4] Calculix:

A three-dimensional structural finite elemente program,

http://www.calculix.de/.

59

[5] T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Civil and Mechanical Engineering, Dover Publications, 2000. [6] J. Li, K. Yu, Noniterative Integration Algorithms with Controllable Numerical Dissipations for Structural Dynamics, International Journal of Computational Methods 15 (3) (2018) 1850111. [7] E. L. Wilson, A computer program for the dynamic stress analysis of underground structure, sesm report no.68-1, division of structural engineering and structural mechanics. University of California, Berkeley. [8] H. M. Hilber, T. J. R. Hughes, R. L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Engineering & Structural Dynamics 5 (3) (1977) 283–292. [9] W. Wood, M. Bossak, O. Zienkiewicz, An alpha modification of Newmark’s method, International Journal for Numerical Methods in Engineering 15 (10) (1980) 1562–1566. [10] J. Chung, G. M. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method, Journal of Applied Mechanics 60 (1993) 371–375. [11] K. Yu, A new family of generalized-α time integration algorithms without overshoot for structural dynamics, Earthquake Engineering & Structural Dynamics 37 (12) (2008) 1389–1409. [12] K. J. Bathe, Conserving energy and momentum in nonlinear dynamics: A simple implicit time integration scheme, Computers & Structures 85 (7 – 8) (2007) 437–445. [13] K. J. Bathe, M. M. I. Baig, On a composite implicit time integration procedure for nonlinear dynamics, Computers & Structures 83 (31 – 32) (2005) 2513–2524.

60

[14] K.-J. Bathe, G. Noh, Insight into an implicit time integration scheme for structural dynamics, Computers & Structures 98 – 99 (2012) 1 – 6. [15] G. Noh, K.-J. Bathe, Further insights into an implicit time integration scheme for structural dynamics, Computers & Structures 202 (2018) 15– 24. [16] C. Kolay, J. M. Ricles, Assessment of explicit and semi-explicit classes of model-based algorithms for direct integration in structural dynamics, International Journal for Numerical Methods in Engineering 107 (1) (2016) 49–73. [17] M. Rezaiee-Pajand, M. Karimi-Rad, More accurate and stable time integration scheme, Engineering with Computers 31 (4) (2015) 791–812. [18] C. Hoff, P. J. Pahl, Development of an implicit method with numerical dissipation for time integration algorithms in structural dynamics, Computer Methods in Applied Mechanics & Engineering 67 (3) (1988) 367–385. [19] C. Hoff, R. L. Taylor, Higher derivative explicit one step methods for non-linear dynamic problems. Part II: Practical calculations and comparisons with other higher order methods, International Journal for Numerical Methods in Engineering 29 (2) (1990) 291–301. [20] J. Li, K. Yu, X. Li, A generalized structure-dependent semi-explicit method for structural dynamics, Journal of Computational and Nonlinear Dynamics 13 (11) (2018) 111008. [21] M. M. Malakiyeh, S. Shojaee, K.-J. Bathe, The Bathe time integration method revisited for prescribing desired numerical dissipation, Computers & Structures 212 (2019) 289–298. [22] G. Noh, K.-J. Bathe, The Bathe time integration method with controllable spectral radius: The p∞-Bathe method, Computers & Structures 212 (2019) 299–310.

61

[23] J. Li, K. Yu, An alternative to the Bathe algorithm, Applied Mathematical Modelling 69 (2019) 255–272. [24] S. Y. Chang, Explicit Pseudodynamic Algorithm with Unconditional Stability, Journal of Engineering Mechanics 128 (9) (2002) 935–947. [25] S.-Y. Chang, A new family of explicit methods for linear structural dynamics, Computers & Structures 88 (1112) (2010) 755 – 772. [26] C. Kolay, J. M. Ricles, Development of a family of unconditionally stable explicit direct integration algorithms with controllable numerical energy dissipation, Earthquake Engineering & Structural Dynamics 43 (9) (2014) 1361–1380. [27] D. J. Maxam, K. K. Tamma, Discussion of “development of a family of unconditionally stable explicit direct integration algorithms with controllable numerical energy dissipation” by Chinmoy Kolay and James M. Ricles, Earthquake Engineering & Structural Dynamics (2018) 6. [28] W. B. Wen, K. Wei, H. S. Lei, S. Y. Duan, D. N. Fang, A novel substep composite implicit time integration scheme for structural dynamics, Computers & Structures 182 (2017) 176–186. [29] S. Dong, BDF-like methods for nonlinear dynamic analysis, Journal of Computational Physics 229 (8) (2010) 3019–3045. [30] W. Kim, S. Y. Choi, An improved implicit time integration algorithm: The generalized composite time integration algorithm, Computers & Structures 196 (2018) 341–354. [31] J. Li, K. Yu, X. Li, A novel family of controllably dissipative composite integration algorithms for structural dynamic analysis, Nonlinear Dynamics 96 (4) (2019) 2475–2507. [32] J. Li, K. Yu, H. He, A second-order accurate three sub-step composite algorithm for structural dynamics, Applied Mathematical Modelling (2019) 42doi:10.1016/j.apm.2019.08.022. 62

[33] C. W. Gear, Numerical initial value problems in ordinary differential equations, Prentice Hall, 1971. [34] H. M. Hilber, T. J. R. Hughes, Collocation, dissipation and ’overshoot’ for time integration schemes in structural dynamics, Earthquake Engineering & Structural Dynamics 6 (1) (1978) 99–117. [35] J. M. Bentez, F. J. Montns, The value of numerical amplification matrices in time integration methods, Computers & Structures 128 (5) (2013) 243– 250. [36] W. T. M. Silva, L. M. Bezerra, Performance of Composite Implicit Time Integration Scheme for Nonlinear Dynamic Analysis, Mathematical Problems in Engineering 2008 (4) (2008) 267–290. [37] W. Kim, H. L. Jin, An improved explicit time integration method for linear and nonlinear structural dynamics, Computers & Structures 206 (2018) 42– 53. [38] T. Liu, C. Zhao, Q. Li, L. Zhang, An efficient backward Euler timeintegration method for nonlinear dynamic analysis of structures, Computers & Structures 106-107 (5) (2012) 20–28.

63