A second-order accurate three sub-step composite algorithm for structural dynamics

Applied Mathematical Modelling 77 (2020) 1391–1412

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

A second-order accurate three sub-step composite algorithm for structural dynamics Jinze Li, Kaiping Yu∗, Haonan He Department of Astronautic Science and Mechanics, Harbin Institute of Technology, No.92 West Dazhi Street, Harbin 150001, China

a r t i c l e

i n f o

Article history: Received 7 February 2019 Revised 1 July 2019 Accepted 19 August 2019 Available online 22 August 2019 MSC: 00-01 99-00 Keywords: Three sub-step algorithm Bathe algorithm Structural dynamics Implicit integration algorithms

a b s t r a c t In this paper, a novel three sub-step composite algorithm with desired numerical properties is developed. The proposed method is a self-starting, unconditionally stable and second-order accurate implicit algorithm without overshoot. Particularly, the second-order accuracy in time is achieved in its final form, but it is not required in each sub-step. Its unique algorithmic parameter is analyzed to achieve the unconditional stability and it shares the identical effective stiffness matrix inside three sub-steps to save the computational cost in linear analyses. The same as the Bathe algorithm, the proposed algorithm is always L-stable, meaning that the spurious high-frequency modes can be effectively eliminated. Three numerical examples are simulated to illustrate the superiority of the proposed algorithm over some existing implicit algorithms. The first numerical simulation, solving a linear single-degree-of-freedom system, shows less period elongation errors and the second-order accuracy of the present scheme. The second one, a clamped-free bar excited by the end load, shows the ability of effectively damping out the unexpected highfrequency modes. The last example solves the nonlinear mass-spring system with variable degree-of-freedoms and illustrates that the composite sub-step algorithm can save more computational cost than the traditional implicit algorithm when the integration step size is selected appropriately. © 2019 Elsevier Inc. All rights reserved.

1. Introduction 1.1. Background In the linear transient analysis, the following second-order equations of motion [1,2] can be obtained by the finite element discretization.

MU¨ (t ) + CU˙ (t ) + KU (t ) = F (t )

(1)

where the M, C and K represent the mass, damping and stiffness matrix, respectively; F(t) denotes the load vector, as a known function of time t; U (t ), U˙ (t ) and U¨ (t ) are the displacement, velocity and acceleration vector, respectively. One of the procedures for solving the equation of motion (1) is the direct time integration algorithm. When solving the linear structural dynamical problems, an alternative procedure, the modal superposition technique [2], can be employed to obtain ∗

Corresponding author. E-mail addresses: [email protected] (J. Li), [email protected] (K. Yu), [email protected] (H. He).

https://doi.org/10.1016/j.apm.2019.08.022 0307-904X/© 2019 Elsevier Inc. All rights reserved.

1392

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

sufficiently accurate numerical solutions. Hence, the direct time integration algorithm is generally used to solve the nonlinear structural problems, which the modal superposition technique cannot effectively deal with. The common direct time integration algorithms include the Newmark-β family of algorithms [3], Wilson-θ algorithm [4], Houbolt method [5], the optimal collation scheme [6], WBZ-α algorithm [7], HHT-α algorithm [8], CH-α algorithm [9], and the GSSSS algorithm [10]. The detailed analysis about these classical integration algorithms can be found in the literature [11]. The above mentioned integration algorithms are all implicit, meaning that a system of linear equations needs to be solved in each time step, and unconditionally stable, meaning that their integration steps are not limited by the stability condition. For the points from [6], there are some shortcomings for these implicit algorithms. For example, these traditional dissipative algorithms, including Wilson-θ algorithm [4], WBZ-α algorithm [7], HHT-α algorithm [8], and CH-α algorithm [9], exhibit the unexpected overshoot behavior because of using the weighted combination of the equations of motion or solving the inaccurate initial accelerations [12]; The GSSSS algorithm [10] without overshoot cannot achieve the maximum of numerical dissipations; And the Newmark-β algorithm [3] achieves numerical dissipations at the expense of accuracy. Recently, some novel direct integration algorithms, such as the composite sub-step algorithm [13,14], the (semi-)explicit structure-dependent algorithm [15,16], the time finite element method [17], the finite difference scheme [18] and other higher order algorithms [19,20], are proposed to improve numerical accuracy and/or stability properties. The composite sub-step algorithm divides an integration step into several sub-steps, and different numerical integration schemes could be used inside each step. As an example, the first composite sub-step algorithm proposed by Bathe et al. [13], named as the (standard) Bathe algorithm in this study, is a single step but two sub-step composite method, where the trapezoidal rule [1] and the three points backward differential scheme [21] are employed in the first and second sub-step, respectively. By generalizing the above procedure, they also proposed a three sub-step algorithm, where the trapezoidal rule is used twice in the first two sub-steps and the Houbolt method [5] is employed in the last sub-step. In addition, various composite sub-step algorithms [22–26] using two or three sub-step strategy have been developed. A simple summary about these existing composite sub-step algorithms can be found in literature [23]. Particularly, Bathe and co-workers have recently developed two new composite sub-step algorithms: the ρ ∞ -Bathe algorithm [25] and β 1 /β 2 -Bathe algorithm [26]. In the case of ρ∞ = 0, the ρ ∞ -Bathe algorithm [25] reduces to the standard Bathe algorithm [13,14], while in the case of ρ∞ = 1 it corresponds to the composite two sub-step trapezoidal rule. Although the β 1 /β 2 -Bathe algorithm [26] presents two algorithmic parameters (β 1 and β 2 ) intuitively, it is essentially a single-parameter algorithm since these two parameters should be appropriately selected to achieve the accuracy condition. In the case of β2 = 2β1 , it only achieves the first-order accuracy and L-stability, and hence solves the wave propagation problems successfully. On the other hand, in the case of β2 = 1 − β1 it possesses the second-order accuracy and controllable numerical dissipations, but this controllably dissipative and secondorder accurate β 1 /β 2 -Bathe algorithm is essentially the ρ ∞ -Bathe algorithm with γ = 0.5. Therefore, among three Bathetype algorithms, only the standard Bathe algorithm [13,14] is L-stable and second-order accurate; the β 1 /β 2 -Bathe algorithm [26] with β2 = 2β1 is L-stable and first-order accurate. Moreover, they need to select the appropriate algorithmic parameters to achieve the identical effective stiffness matrices inside each sub-step. In addition, the composite sub-step explicit algorithm is also developed to significantly reduce the computational cost while preserving the advantages of the composite sub-step strategy. Generally, the composite sub-step explicit method only achieves the conditional stability. Considering the fact that the central difference method [1] is the most popular explicit method, but its integration step size h is controlled by h ≤ 2/ωmax and it is zero-dissipative. Hence, some composite sub-step explicit methods [27–30] are proposed to achieve the larger integration step size and/or numerical dissipations. These explicit methods succeed solving wave propagation and some nonlinear problems. In terms of the computational cost of the implicit scheme, time that the composite n sub-step algorithm elapses is generally n times of that of the traditional implicit single step and single sub-step algorithm, such as the CH-α algorithm [9], when solving the linear elastic problems using the same integration step. But generally, the implicit composite sub-step algorithm can achieve the unconditional stability and hence can employ larger time step to significantly save the computational cost [31]. It is necessary to point out that for the traditional implicit algorithms, using the larger integration step, can also reduce the computational cost, but contribute to two additional side effects. One is introducing more unexpected numerical dissipations and dispersions in the important low-frequency region, resulting in the undesired numerical errors. The other is requiring more nonlinear iterations to converge to the satisfactory solutions within each step, leading to more computational time. The desired numerical property of the composite sub-step algorithm is that the identical effective stiffness matrices can be achieved inside each sub-step to further save the computational cost. This ensures that in the case of solving linear problems, only one matrix decomposition of effective stiffness matrices is required during the whole simulation.√For example, the Bathe algorithm [13,14] achieves this numerical property only when its parameter γ is equal to be 2 ± 2 while other composite sub-step algorithms [22,24,32] fail to possess it. 1.2. Objective and scope In this paper, particular attention is paid to the composite sub-step implicit algorithm used to effectively solve structural dynamical problems (1). The proposed composite sub-step algorithm is a single step but three sub-step implicit algorithms, where the whole integration interval [t, t + h] is split into three sub-step: [t, t + γ1 h], [t + γ1 h, t + γ2 h] and [t + γ2 h, t + h]. As pointed out in [33], the algorithmic parameters γ 1,2 are allowed to be greater than one. Further, the proposed algorithm is designed to achieve the following desired numerical characteristics:

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412 •

•

•

• •

•

1393

The second-order accuracy in the displacement, velocity and acceleration is obtained in its final form, but it is not required in each sub-step. The unconditional stability is achieved by analyzing its algorithmic parameters. It ensures that the proposed algorithm can use a larger integration step to solve structural dynamical problems to reduce the computational cost. The same as the Bathe algorithm [13], the proposed algorithm is also L-stable, which indicates that the spurious highfrequency modes can be damped out effectively and quickly. There is no overshoot for the proposed algorithm when any nonzero initial displacements and/or velocities are used. Independent of the algorithmic parameter, the identical effective stiffness matrices inside three sub-step are always ensured to save the computational cost in the linear analysis. The proposed algorithm is self-starting. This numerical property is not embedded into the composite sub-step algorithm presented in [22].

The remaining text is organized as follows. The newly proposed three sub-step algorithm is presented in Section 2, where numerical schemes used in each step are given in detail. The algorithmic parameter γ 1 is analyzed to achieve the unconditional stability in Section 3. The numerical accuracy including the amplitude decay and the period elongation is presented in Section 4. The analysis of the overshoot behavior is carried out in Section 5 by solving a linear single-degreeof-freedom system. In Section 6, three numerical examples are simulated to show the advantage of the proposed algorithm over other implicit algorithms. Results and discussions are presented in Section 7. Conclusions are drawn in Section 8. 2. Composite three sub-step algorithm In this section, a new composite three sub-step algorithm is derived in detail, where the whole integration step h is divided into three sub-step: [t, t + γ1 h], [t + γ1 h, t + γ2 h] and [t + γ2 h, t + h]. In the first sub-step t ∈ [t, t + γ1 h], the backward Euler scheme [34] is employed

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

(2)

Ut+γ1 h = Ut + γ1 hU˙ t+γ1 h

(3)

which yields

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

(4)

The equation of motion in linear systems at time t + γ1 h is used to obtain the acceleration solutions, that is

MU¨ t+γ1 h + CU˙ t+γ1 h + KUt+γ1 h = Ft+γ1 h

(5)

The above Eqs. (2)–(5) can be rewritten as follows:

K 1U¨ t+γ1 h = F 1

(6)

K 1 = (γ1 h )2 K + (γ1 h )C + M

(7)

F 1 = Ft+γ1 h − CU˙ t − K (Ut + γ1 hU˙ t )

(8)

where

The acceleration solutions at time t + γ1 h are obtained by solving Eq. (6). Then, the velocity and displacement responses are given by substituting U¨ t+γ1 h into Eqs. (2) and (3), respectively. In the second sub-step t ∈ [t + γ1 h, t + γ2 h], the following backward difference scheme is employed, which is

hU˙ t+γ2 h = c1Ut + c2Ut+γ1 h + c3Ut+γ2 h

(9)

hU¨ t+γ2 h = c1U˙ t + c2U˙ t+γ1 h + c3U˙ t+γ2 h

(10)

with the equation of motion at time t + γ2 h

MU¨ t+γ2 h + CU˙ t+γ2 h + KUt+γ2 h = Ft+γ2 h

(11)

where the coefficients ci , i = 1, 2, 3 are determined by

c1 =

γ2 − 2γ1 γ1 − γ2 c2 = γ12 γ12

c3 =

1

γ1

(12)

Similar to the first sub-step, the following form can be obtained in the second sub-step.

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

(13)

1394

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

where

 2

h K2 = c3

 

h C+M c3

K+

F 2 =Ft+γ2 h − C



−K

c

1

c3

U˙ t +

c2 U˙ c3 t+γ1 h

(14)



1 h (c1Ut + c2Ut+γ1 h ) + 2 (c1U˙ t + c2U˙ t+γ1 h c3 c3

 (15)

In the last sub-step t ∈ [t + γ2 h, t + h], the similar numerical scheme is used to update the numerical responses.

hU˙ t+h = d1Ut + d2Ut+γ1 h + d3Ut+γ2 h + d4Ut+h

(16)

hU¨ t+h = d1U˙ t + d2U˙ t+γ1 h + d3U˙ t+γ2 h + d4U˙ t+h

(17)

with the equation of motion at time t + h

MU¨ t+h + CU˙ t+h + KUt+h = Ft+h

(18)

Using Eqs. (16) and (18), the following equation can be obtained:

K 3U¨ t+h = F 3 where

 K3 =

h d4

(19)

2

 K+

 F 3 =Ft+h − C +



h C+M d4

d1 d2 d3 U˙ t + U˙ t+γ1 + U˙ t+γ2 d4 d4 d4

(20)



 −K

d1 d1 h Ut + 2 U˙ t d4 d4

d2 d2 h d3 d3 h U + 2 U˙ t+γ1 h + Ut+γ2 h + 2 U˙ t+γ2 h d4 t+γ1 h d4 d4 d4



(21)

In order to save the computational cost in linear analysis, three effective stiffness matrices (i.e. K i , i = 1, 2, 3) in the new composite three sub-step algorithm should be designed as the identical one [35]. Therefore, the d4 is considered to be

d4 =

1

(22)

γ1

The remaining coefficients di , i = 1, 2, 3 in the third sub-step are determined by

d1 =

−6γ12 + 6γ1 − 1 2γ13

(23)

d2 =

γ13 − (2γ2 + 1 )γ12 + (3γ1 − 1/2 )γ2 γ13 (γ1 − γ2 )

(24)

d3 =

γ12 − 2γ1 + 1/2 γ12 (γ1 − γ2 )

(25)

to achieve the second-order accuracy in time. The numerical test about its second-order accuracy is given in Section Numerical experiments. The parameter γ 2 is calculated by

γ2 = 2γ1

(26)

to reduce the number of unknown parameters. Obviously, there is a parameter γ 1 for the new three sub-step algorithm. In next section, the only parameter γ 1 will be analyzed to achieve the unconditional stability of the composite algorithm. In the end of this section, the pseudo code of the new three sub-step algorithm for solving nonlinear problems (27) in the finite element model is summarized in Algorithm 1.

MU¨ (t ) + R(U˙ (t ), U (t )) = F (t )

(27)

where R(U˙ (t ), U (t )) is the vector of nodal forces determined by the displacement and velocity. If the modified Newton–Raphson iteration scheme is employed, some mild modifications are required in Algorithm 1, and (1 ) the value of Ut+ predicted in step 31 is obtained by using some other numerical schemes, such as the Hermite interpolation h [34].

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

1395

Algorithm 1 Three sub-steps composite algorithm for nonlinear problems Set: parameter γ1 and calculate the corresponding γ2 using (26). Compute: the initial acceleration U¨ 0 from U˙ 0 and U0 . 3: for t = t0 to t = tend do //———————–In the first sub-step——————————– 4: ( γ h )2 5: Predict: U (1 ) = Ut + γ1 hU˙ t + 1 U¨ t . 1:

2:

t+γ1 h

6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34:

2

for i = 1 to max − iter do (i ) (i ) Compute: U˙ t+ γ h = (Ut+γ

− Ut )/(γ1 h ).

1h

1

(i ) ˙ (i ) ˙ Compute: U¨t+ γ1 h = (Ut+γ1 h − Ut )/ (γ1 h ). (i ) (i )  = 1 M + 1 C (i ) + K (i ) . ˙ Compute: R(Ut+γ h , Ut+γ h ) and K γ1 h t+γ1 h t+γ1 h γ12 h2 1 1 (i ) (i ) ( i ) ( i )  ¨ (i ) . U Solve: U using K U = Ft+γ1 h − R(U˙ t+ , U ) − M γ1 h t+γ1 h t+γ1 h

if convergence then (i ) Compute: Ut+γ1 h = Ut+ γ

1h

Break. end if (i+1 ) (i ) Update: Ut+ γ h ← Ut+γ

1h

1

and U˙ t+γ1 h , U¨t+γ1 h .

+ U (i ) .

end for //———————–In the second sub-step——————————– ( γ h )2 Predict: U (1 ) = Ut+γ h + γ1 hU˙ t+γ h + 1 U¨ t+γ h . t+γ2 h

1

2

1

1

for i = 1 to max − iter do (i ) (i ) Compute: U˙ t+ γ h = (c1Ut + c2Ut+γ1 h + c3Ut+γ h )/h. 2

2

(i ) ˙ ˙ ˙ (i ) Compute: U¨t+ γ2 h = (c1Ut + c2Ut+γ1 h + c3Ut+γ2 h )/h. 2  = c3 M + c3 C ( i ) Compute: R(U˙ (i ) , U (i ) ) and K t+γ2 h

Solve: U (i ) using

t+γ2 h U (i ) K

if convergence then (i ) Compute: Ut+γ2 h = Ut+ γ Break. end if (i+1 ) (i ) Update: Ut+ γ h ← Ut+γ 2

2h

(i ) + Kt+ h t+γ2 h γ2 h . ( i ) ( i ) (i ) = Ft+γ2 h − R(U˙ t+γ h , Ut+γ h ) − MU¨ t+ γ2 h . 2 2 h2

2h

and U˙ t+γ2 h , U¨ t+γ2 h .

+ U (i ) .

end for //———————–In the last sub-step——————————– (1 ) Predict: Ut+ = Ut+γ2 h + (1 − γ2 )hU˙ t+γ2 h + [(1 − γ2 )h]2U¨ t+γ2 h /2. h for i = 1 to max − iter do (i ) (i ) Compute: U˙ t+ = (d1Ut + d2Ut+γ1 h + d3Ut+γ2 h + d4Ut+ )/h. h h (i ) (i ) Compute: U¨ t+ = (d1U˙ t + d2U˙ t+γ1 h + d3U˙ t+γ2 h + d4U˙ t+ )/h. h h 2

35: 36: 37: 38: 39: 40: 41: 42: 43: 44:

d (i ) (i ) (i ) (i ) Compute: R(U˙ t+ , Ut+ ) and K = h42 M + dh4 Ct+ + Kt+ . h h h h (i ) (i ) (i ) ( i ) ( i )  ˙ Solve: U using K U = Ft+h − R(Ut+h , Ut+h ) − MU¨ t+ . h if convergence then (i ) Compute: Ut+h = Ut+ and U˙ t+h , U¨ t+h . h Break. end if (i ) (i ) Update: Ut+ ← Ut+ + U (i ) . h h end for //———————Next integration step——————————— end for

3. The analysis of γ 1 The stability study of direct integration algorithms designed for structural dynamics is inevitably linked to the following single-degree-of-freedom (SDOF) system

x¨ (t ) + 2ξ ωx˙ (t ) + ω2 x(t ) = 0

(28)

1396

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

Fig. 1. Stability region in the (I1 , I2 ) plane [2], where ρ = max(|λi | ), i = 1, 2, 3 represents the spectral radius of the integration algorithm.

where the ξ and ω are the viscous damping ratio and natural frequency, respectively. Applying the new three sub-step algorithm to this SDOF system yields



xt+h hx˙ t+h h2 x¨t+h





xt = A hx˙ t h2 x¨t



(29)

where A is the amplification matrix of the three sub-step algorithm, whose explicit derivation is given in Appendix A. The characteristic polynomial of matrix A is given by calculating [1]

det(λI − A ) = 0

(30)

which gives

λ3 − 2I1 λ2 + I2 λ = 0

(31)

In the case of ξ = 0, the explicit expressions of I1 and I2 are

I1 = −

I2 =

(−12γ14 + 16γ13 − 3γ12 )4 + (1 − 6γ12 )2 − 2 2(1 + γ12 2 )3

(9γ14 − 18γ13 + 12γ12 − 3γ1 + 1/4 )4 + 3γ12 2 + 1 (1 + γ12 2 )3

(32)

(33)

In order to achieve the unconditional stability, the I1 and I2 should satisfy the following conditions [1]



−(I2 + 1 ) ≤ 2I1 ≤ I2 + 1, −1 < I1 < 1,

−1 ≤ I2 < 1 I2 = 1

(34)

The above stability conditions describe in (I1 , I2 ) space three straight lines forming a triangular region, as shown in Fig. 1. Substituting Eqs. (32) and (33) into conditions (34) and after some algebraic manipulations yields

f := −6γ14 + 18γ13 − 12γ12 + 3γ1 − 1/4 ≥ 0

(35)

The condition (35) requires that the value of γ 1 should satisfy

γ1 ∈ [0.181, 2.186]

(36)

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

1397

Fig. 2. Spectral radii (ρ ) of various direct integration algorithms.

d AD

Exact Numerical PE

0

T0/2

T0

T

t

Fig. 3. Definition of the period elongation error (PE) and amplitude decay (AD).

to achieve the unconditional stability. Fig. 2a presents spectral radii of the proposed three sub-step algorithm with different γ 1 ∈ [0.181, 2.186]. Fig. 2a indicates that the proposed composite algorithm achieves the unconditional stability when the parameter γ 1 satisfies the condition (36). Further, it is shown that the proposed composite algorithm is asymptotically annihilating (L-stable) since its spectral radii tend to zero as the  increases into infinity. It means that the new composite algorithm can effectively eliminate the spurious high-frequency modes. Fig. 2b gives√ the comparison of spectral radii for different L-stable integration algorithms, including the Bathe algorithm with γ = 2 − 2, CH-α algorithm with ρ∞ = 0.0, and the Park algorithm [36]. Fig. 2 reveals that the numerical accuracy in the low-frequency region can be adjusted by changing γ 1 to some extent. Although this desirable numerical property can be embedded into the Bathe algorithm by controlling its parameter γ [33], the effective stiffness matrices inside two sub-step are not the same when the γ varies. 4. Numerical accuracy In this section, the amplitude decay [1,26,33] and period elongation error [26], are used to show the numerical characteristics of the proposed three sub-step algorithm. The mathematical expressions of the amplitude decay (AD) and period elongation error (PE), which are illustrated in Fig. 3, can be summarized as follows [33]

PE =

T − T0  = −1 T0 d

AD = 1 − exp(−2π ξ

 ) d

(37)

(38)

1398

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

a

b

Fig. 4. Numerical accuracy of the proposed algorithm with different γ 1 : (a) Amplitude decay and (b) Period elongation.

a

b

Fig. 5. Comparison of numerical accuracy for various algorithms: (a) Amplitude decay and (b) Period elongation.

in which the numerical damping ratio ξ and the d are determined by

ξ =−

1



ln(ρ ),

d = tan−1 ( I2 − I12 /I1 )

(39)

In the above equation, the I1 and I2 are from Eqs. (32) and (33), respectively. Fig. 4 plots amplitude decays and period elongation errors of the proposed three sub-step algorithm as the γ 1 varies. From these two figures, for the given smaller γ 1 , the desired numerical accuracy in the low-frequency region can be obtained. For example, in the case of γ1 = 0.181 or 0.2, there are the small amplitude decays and period elongation errors compared with the other cases. In fact, when the γ 1 is less than 0.5, which makes γ 2 less than one, the desired numerical accuracy can be obtained. On the other hand, the significant amplitude decays and period elongation errors can be found in the low-frequency region for the case of larger γ 1 , such as γ1 = 2.0. The comparative study about amplitude decays and period elongations is presented in Fig. 5, where the Bathe algorithm √ [13,14] with γ = 2 − 2, CH-α algorithm [9] with ρ∞ = 0.0, the classical Newmark algorithm [1,3] with 2β = γ > 0.5, and the β 1 /β 2 -Bathe algorithm [26] with β2 = 2β1 are included for comparison. Some basic observations from Fig. 5 can be summarized as follows, •

The proposed algorithm with smaller γ 1 provides the most accurate results in terms of the amplitude decay and period elongation in the low-frequency region.

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

a

1399

b

Fig. 6. Displacement and velocity errors in the first step for different integration algorithms when initial conditions d0 = 1, v0 = 0 is used.

•

•

•

For larger γ 1 , the proposed algorithm could give some acceptable numerical accuracy compared with the traditional implicit algorithm CH-α algorithm with ρ∞ = 0.0. As pointed out in literature [1,25], the β 1 /β 2 -Bathe algorithm with β2 = 2β1 and the Newmark algorithm with 2β = γ > 0.5 only achieve the first order accuracy. Hence, there are very significant amplitude decays in the low-frequency region for these two methods, as shown by the dotted lines in Fig. 5a. In addition, the fact [37] has been confirmed in Fig. 5a that the amplitude decay curves of other second-order accurate methods are tangent to the -axis at  = 0.0. As shown in Fig. 4b, the curves of period elongations for the proposed algorithm with γ 1 > 0.5 firstly decrease into the negative minimum and then increases continuously in the positive direction.

In a word, the proposed algorithm can provide the desired numerical accuracy in the case of γ 1 ∈ [0.181, 0.5]. For the case of γ 1 ∈ (0.5, 2.186], the proposed algorithm could be suitable for solving some structural dynamical problems, see Section 6.2.

5. Overshoot Most traditional implicit dissipative algorithms, such as the CH-α [9], HHT-α [8], WBZ-α [7] and Wilson-θ [4] algorithm, suffer from overshoot that is typically observed when solving problems with large eigenfrequencies and simulations with relatively large integration steps. Overshoot describes the tendency of an integration algorithm to strongly overestimate the numerical response of the structural problem in the first time steps. The first ones to systematically analyze overshoot behavior are Hilber and Hughes [6] who showed that the large elements in powers of the amplification matrix could be a factor leading to this phenomenon. The overshoot behavior of the above implicit dissipative algorithms has been eliminated by using some techniques, such as introducing more algorithmic parameters [38] and using the structure-dependent parameter [39]. In this section, a linear SDOF system [38] is solved to numerically analyze the overshoot behavior of the proposed algorithms. The natural frequency ω and physical damping ratio ξ are considered to be 2π rad/s and 0.5, respectively. First, considering the first case of initial conditions d0 = 1, v0 = 0, Fig. √ 6 shows displacement and velocity errors in the first step, where numerical results from the Bathe algorithm with γ = 2 − 2 and CH-α algorithm with ρ∞ = 0.0 are given for comparison. Fig. 6 indicates that there exists no overshoot tendency for the proposed composite algorithm with any γ 1 and the Bathe algorithm in the first case of initial conditions, while the very significant errors can be found for CH-α algorithm as the time step h increases, revealing its overshoot behaviors. The same conclusion can be found for these algorithms in the second case of initial conditions d0 = 0, v0 = 1, as shown in Fig. 7. The reason why CH-α algorithm with ρ∞ = 0.0 exists very significant overshoot tendency can be explained as follows. Considering  → ∞, the displacement and velocity of CH-α algorithm at the first step can be written as [6]

where

d1 ≈ c11 d0 + c12 v0 h

(40)

v1 ≈ c21 ωd0 + c22 v0

(41)

1400

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

a

b

Fig. 7. Displacement and velocity errors in the first step for different integration algorithms when initial conditions d0 = 0, v0 = 1 is used.

c11 = c21

( ρ∞ − 1 ) 2 ξ

,

2 ( ρ∞ − 1 ) 2 = , 4

c12 = (ρ∞ − 1 )2 ξ 2 c22 = −

( ρ∞ − 1 ) 2 ξ 2



(42)

Hence, in the case of ξ = 0, CH-α algorithm exhibits a tendency to overshoot linearly in  in the velocity because of the nonzero initial displacement and/or velocity. Moreover, the d1 of CH-α algorithm also overshoots linearly in  due to initial displacement and linearly in h due to initial velocity. Note that for the case of ρ∞ = 1, the coefficients (42) are not suitable and CH-α algorithm, essentially reducing to the trapezoidal rule in this case, exists no overshoot tendency. 6. Numerical experiments To check the numerical performance of the proposed three sub-step algorithm, some linear and nonlinear examples are solved in this section. 6.1. A linear SDOF system In this section, a simple linear SDOF system [26] is considered to highlight some observations mentioned and test the convergence ratio of the proposed algorithm.

u¨ (t ) + 4u(t ) = 0,

u(0 ) = 1, u˙ (0 ) = 0

(43)

whose exact displacement solution √ is u(t ) = cos(2t ). Herein, two new three sub-step algorithm (γ1 = 0.181 and 2.186), two Bathe algorithms with γ = 2 ± 2, trapezoidal rule and CH-α algorithm with ρ∞ = 0.0 are employed to solve this SDOF system. The integration step h is different for different algorithms to ensure the identical computational cost for the same integration interval. Figs. 8 and 9 show displacement responses of this linear SDOF system, where the√proposed algorithm with γ1 = 0.181 can provide more accurate numerical solutions than the Bathe algorithm with γ = 2 − 2 and CH-α algorithm. The unexpected numerical dissipations in the low-frequency region can be confirmed for the proposed algorithm with γ1 = 2.186 and the √ Bathe algorithm with γ = 2 + 2. This numerical test also indicates that the trapezoidal rule [1] is a good choice to solve such simple problems without spurious high-frequency modes. Fig. 10 plots the results about convergence tests for different algorithms, where the global error norms [23] are employed to calculate the errors in the displacement, velocity and acceleration, which are



Edis =

N i=0

 Evel =

N i=0

 Eacc =

[ui − u(ti )]2

[vi − u˙ (ti )]2

N i=0

N i=0

N

[ai − u¨ (ti )]2

i=0

[u(ti )]2

(44)

[u˙ (ti )]2

(45)

[u¨ (ti )]2

(46)

N i=0

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

Fig. 8. Displacement responses of the linear SDOF system for t ∈ [0, 8].

Fig. 9. Displacement responses of the linear SDOF system for t ∈ [45, 50].

Fig. 10. Convergence rates of different algorithms for solving the linear SDOF system: (a) Displacement, (b) Velocity and (c) Acceleration.

1401

1402

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

Fig. 11. Convergence rates of the β 1 /β 2 -Bathe algorithm [26] for solving the linear SDOF system: (a) Displacement, (b) Velocity and (c) Acceleration.

Fig. 12. A clamped-free bar excited by the end load.

where ui , vi and ai are the numerical displacement, velocity and acceleration responses of the linear SDOF system. The second order accuracy of the proposed algorithm with γ 1 ∈ [0.181, 2.186] can be confirmed in Fig. 10. Moreover, the first order accuracy in acceleration of CH-α algorithm is easily found in Fig. 10c. Fig. 11 presents the convergence ratio of the β 1 /β 2 -Bathe algorithm [26] with β2 = 2β1 , where the first order accuracy in the displacement, velocity and acceleration can be observed. 6.2. A clamped-free bar excited by end load In this section, a clamped-free bar [2,26] subjected to an stepwise end load F(t), as shown in Fig. 12, is considered. The numerical model is made up of N = 100 elements of equal length l = L/N, and the lumped mass matrix is used in this model. The lumped mass and stiffness matrices are

⎡

2 ⎢−1

⎢ ⎢ EA ⎢ K= ⎢ l ⎢ ⎢ ⎣

⎤

−1 2

−1

−1

0 ..

2 .. .

.

..

. −1

0

−1 2 −1

⎡

2

⎥ ⎥ ⎥ ⎥ ⎥ ⎦

0 2

ml ⎢ M= ⎢ 2 ⎢

⎣

1

(47)

N×N

⎤

2

⎢ ⎢

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ −1⎦

.. 0

. 2 1

(48)

N×N

where the physical parameters E = 3 × 107 , A = 1, L = 200 and m = 7.3 × 10−4 are used in this model. The applied end load is considered as F (t ) = 1 × 104 . The integration step size h used by all algorithms shown herein satisfies

h = (9.88 × 10−7 ) × CFL

(49)

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

1403

a

b

Fig. 13. Numerical displacements and velocities at node 50 using CH-α algorithm with ρ∞ = 0.0.

Fig. 14. Section A-A of Fig. 13.

Firstly, the CH-α algorithm [9] with ρ∞ = 0.0 is employed to solve the above model. The displacement and velocity responses at node 50 are plotted in Fig. 13, where the exact displacement solution of this bar is [2]

u(x, t ) =



∞

8F (t )L F (t ) x− 2 EA π EA s=1



(−1 )s−1 ( 2s − 1 )π x (2s − 1 )π ct sin cos 2L 2L ( 2s − 1 )2



(50)

in which c = EA/m denotes the wave propagation velocity in the bar. Figs. 13 and 14 indicate that for the CH-α algorithm, there are very significant numerical oscillations in the displacement and velocity compared with the exact solutions. Next, the standard Bathe algorithm [13,14] with γ = 1.9 and the β 1 /β 2 -Bathe algorithm [26] with β1 = 0.75 and β2 = 2β1 are used to solve this model. As shown in Figs. 15 and 16, the β 1 /β 2 -Bathe algorithm can provide more accurate numerical responses than the standard Bathe algorithm with γ = 1.9. Particularly, very significant numerical oscillations in velocity can be found for the Bathe algorithm with γ = 1.9 since it introduces less numerical dissipations in the low-

1404

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

a

b

Fig. 15. Numerical displacement and velocities at node 50 using the Bathe algorithm with γ = 2 −

a

b

Fig. 16. Numerical displacement and velocities at node 50 using the β 1 /β 2 -Bathe algorithm.

√ 2.

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

1405

a

b

Fig. 17. Numerical displacements and velocities at node 50 using the first-order accurate Newmark algorithm.

frequency region than the β 1 /β 2 -Bathe algorithm. It is necessary to point out that although the Bathe-β 1 /β 2 algorithm with β2 = 2β1 predicts better numerical solutions, there are two undesired shortcomings: (i) it only achieves the first-order accuracy in the displacement, velocity and acceleration, as shown in Fig. 11; (ii) it does not share the identical effective stiffness matrices inside two sub-step. Note that it is the first-order accuracy that results in its larger amplitude decays in the low-frequency region ( ≤ 0.1), as shown in Fig. 5a. More amplitude decays in the low-frequency region makes it be capability of accurately solving this kind of wave propagation problems. It is well-known that the classical Newmark algorithm [1,3] only achieves the first-order accuracy when numerical dissipations (namely, amplitude decays) are required. In this case, two algorithmic parameters of the Newmark algorithm should fulfill the requirement: 2β ≥ γ > 0.5. Herein, the first-order accurate Newmark algorithm with γ = 1.3 and β = γ /2 is used to solve this model. Numerical results are presented in Fig. 17. Obviously, Figs. 16 and 17 illustrate that two first-order accurate algorithms are more suitable for solving this model due to their more amplitude decays in the low-frequency region, as shown in Fig. 5a. It necessary to point out that in the case of using the same integration step, the Newmark algorithm needs less computational cost than the β 1 /β 2 -Bathe algorithm due to less calculations in each step. Last, the numerical model can be resolved by using the proposed composite algorithm with γ1 = 0.75. Fig. 18 shows numerical displacements and velocities at node 50. Obviously, in the case of CFL=1, there are very significant numerical oscillations in velocity for the proposed algorithm with γ1 = 0.75. More accurate numerical solutions can be obtained as the CFL increases or the value of γ 1 is appropriately chosen. In fact, numerical solutions also exist very significant oscillations in velocity if the proposed algorithm with γ1 = 0.181 is used (although its result is not presented herein) because, similar to CH-α algorithm with ρ∞ = 0.0 shown in Fig. 13b, there is less numerical dissipations in the low-frequency region. In order to select an appropriate parameter γ 1 , users need to carry out more than one numerical experimentations. The parameter γ1 = 0.181 is firstly suggested simulating structural dynamical problems by default. Then, according to numerical results obtained, users could appropriately increase the value of γ 1 to yield the satisfactory numerical solutions. 6.3. An n-degree-of-freedom nonlinear mass-spring system An n-degree-of-freedom nonlinear stiffness system [40] shown in Fig. 19 is considered in this subsection to compare the computational cost between the composite sub-step algorithm and the traditional implicit algorithm. All physical properties of this model are considered as

mi = 1

ki =

k

i = 1, . . . , n

(51) i=1

k[1 + α (ui − ui−1 )2 ] 2 ≤ i ≤ n

where the parameter α = −2 is used in this test.

(52)

1406

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

Fig. 18. Numerical displacements and velocities at node 50 using the proposed composite algorithm.

Fig. 19. Model of the n-degree-of-freedom nonlinear mass-spring system.

Firstly, the case of n = 50 is considered to test the numerical accuracy of various integration algorithms for solving nonlinear systems, where time step h of the three sub-step algorithm, two sub-step algorithm, and traditional implicit algorithm is taken to be 0.03, 0.02 and 0.01, respectively. Fig. 20 presents numerical displacement and velocity responses obtained from the proposed algorithm with γ1 = 0.181 and 0.4, the Bathe algorithm, the β 1 /β 2 -Bathe algorithm algorithm, the trapezoidal rule, and CH-α algorithm with ρ∞ = 0.0. Note that the reference solutions are calculated by using higher order integration scheme [20] for comparison. Figs. 20 and 21 reveal that the proposed algorithm can provide more accurate numerical solutions by choosing the appropriate γ 1 . Particularly, there are very significant amplitude errors for the β 1 /β 2 -Bathe algorithm [26] and CH-α algorithm since they introduce more numerical dissipations in the low-frequency region. Hence, the integration method with less dissipations in the low-frequency region can provide more accurate numerical results for this model, such as the proposed algorithm with γ1 = 0.181 and the Bathe algorithm. Next, three values of n (50 0, 10 0 0 and 1500) in this model are considered to simulate three different nonlinear stiffness mass-spring systems in order to test the computational cost. The time steps used by various algorithms are the same as the above case. Fig. 22 presents displacement responses at the largest degree-of-freedom of three nonlinear stiffness systems, which demonstrates that although two composite √ sub-step algorithms (i.e. the proposed three sub-step algorithm with γ1 = 0.181 and the Bathe algorithm with γ = 2 − 2) use larger integration step, the almost same accurate numerical solutions as two traditional implicit algorithms can be provided. Besides, the CPU time elapsed by four algorithms is recorded and presented in Fig. 23. Two primary facts can be found in Fig. 23. One is that the CPU time used by the proposed algorithm is almost identical to that of the Bathe algorithm since their step sizes are proportional to the number of sub-step. Another is that as the number of degree-of-freedoms increases, the CPU time elapsed by the traditional implicit algorithms is significantly more than that of the composite sub-step methods. When solving nonlinear systems, each

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

1407

a

b

Fig. 20. Displacement and velocity responses of the nonlinear mass-spring system using various algorithms in the case of n = 50.

Fig. 21. Section A-A of Fig. 20a.

(sub-)step requires the nonlinear iterative schemes to converge to the expected solutions. The number of iterations in each (sub-)step is generally unknown for different integration algorithms, which is the main reason why the computational time elapsed by the traditional implicit algorithms and composite sub-step schemes is significantly different. Combined with Fig. 22, Fig. 23 indicates that the composite sub-step algorithms can employ larger time steps to reduce the computational cost without sacrificing accuracy.

1408

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

Fig. 22. Displacement responses of three nonlinear mass-spring systems: (a) n = 500, (b) n = 1000 and (c) n = 1500.

Fig. 23. The CPU time elapsed by four algorithms when solving three nonlinear mass-spring systems.

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

1409

Last, it is not difficult to imagine that if the same integration step is used for these four algorithms, the computational time elapsed by composite sub-step algorithms is definitely more than that of traditional implicit algorithms since the former requires more computational operations in the whole integration interval. As a result, more accurate numerical results can be provided by the composite algorithms. 7. Results and discussions 7.1. Result statements In this paper, the present three sub-step algorithm includes the unique algorithmic parameter γ 1 , which should satisfy the stability condition γ 1 ∈ [0.181, 2.186] to achieve the unconditional stability (L-stability). The spectral radii, amplitude decays and period elongation errors are calculated for the present scheme with different γ 1 , and their comparisons with other existing algorithms are given. The spectral analysis shows that the present scheme indeed achieves the L-stability. When the γ 1 decreases from 2.186 into 0.181, the numerical accuracy in the low-frequency region increases. In other word, in the case of smaller γ 1 , smaller amplitude decays and period elongation errors can be found in the important low-frequency region and larger γ 1 indicates more amplitude decays and the unexpected period elongation errors. Hence, the algorithmic parameter γ1 = 0.181 is first recommended for inexperienced users to use. Compared with two L-stable Bathe algorithms √ (the standard Bathe algorithm [13,14] with γ = 2 − 2 and the β 1 /β 2 -Bathe algorithm [26] with β2 = 2β1 ), the present scheme can provide less amplitude decays and period elongation errors in the low-frequency region. The complete overshoot analysis shows that the present scheme exists no overshoot in both displacement and velocity when any nonzero initial conditions are used, while the traditional implicit dissipative CH-α algorithm [9] presents very significant overshoot results. Three numerical examples are simulated to confirm the above theoretical results. In detail, (i) the first numerical example, a standard SDOF system, shows the second-order accuracy and the smaller period error of the proposed algorithm and confirms that the β 1 /β 2 -Bathe algorithm [26] with β2 = 2β1 is only the first-order accurate in the displacement, velocity and acceleration. In this example, the parameter γ1 = 0.181 can be chosen in default to solve the structural dynamical problems because of its smaller amplitude decays and less period elongation errors in the low-frequency region. This numerical example also indicates that the trapezoidal rule [1] is also a good choice to solve this kind of simple systems without spurious high-frequency modes. (ii) The numerical model of a clamped-free bar excited by the end load shows the ability of the present scheme in terms of effectively eliminating the spurious numerical oscillations. In this test, the first-order accurate β 1 /β 2 -Bathe algorithm [26] with β2 = 2β1 and Newmark algorithm [3] with 2β = γ > 0.5 provide more accurate numerical solutions than other second-order algorithms due to their larger amplitude decays in the low-frequency region. This test also illustrates that users maybe need to solve more than one numerical experimentations in order to select an appropriate parameter (such as, the γ 1 of the proposed algorithm, the β 1,2 of the β 1 /β 2 -Bathe algorithm [26] and the γ of the standard Bathe algorithm [13,14]). (iii) An n-degree-of-freedom nonlinear stiffness mass-spring system indicates that in case of the identical integration step, the composite sub-step algorithm can significantly save more computational cost than the traditional implicit algorithm when the almost same solutions are provided. 7.2. Discussions The present three sub-step algorithm fully implements the six desired numerical characteristics mentioned in the Section 1.2 and its unique algorithmic parameter γ 1 should be chosen in the range of γ 1 ∈ [0.181, 2.186], where the γ1 = 0.181 is recommended by default. As for these six numerical properties, the implicit integration method should generally achieve two most common properties: the second-order accuracy in time and unconditional stability, because the implicit algorithms without these two properties are very limited in use. For example, the first-order β 1 /β 2 -Bathe algorithm with β2 = 2β1 may only be considered for use in solving wave propagation problems, but in this case the inexpensive first-order Newmark algorithm could be considered as an alternative, as shown in the second numerical example. Fortunately, the almost all implicit composite sub-step algorithms easily possess these two properties except for the one developed in [41]. The numerical dissipation is also regarded as the desired properties of the integration method since the dissipative algorithm can effectively eliminate the spurious high-frequency modes and stabilize the nonlinear iterative process. The traditional dissipative algorithms exist significant overshoot behavior due to using the weighted equilibrium equations. The present scheme introduces the maximum of numerical dissipations but exists no overshoot. In fact, all implicit composite sub-step algorithms do not have the overshoot tendency. Therefore, the superiority of the composite sub-step technique in terms of developing dissipative algorithms without overshoot is very obvious. Herein, it is necessary to point out that for structural dynamical problems without spurious high-frequency modes, the non-dissipative trapezoidal rule is still the best choice. Besides, a well-designed integration algorithm should ensure its self-starting, which is not embedded into the algorithm developed in [22]. In contrary to three Bathe-type algorithms and other algorithms [22,24], the present scheme shares the identical effective stiffness matrices within each sub-step, independent of its algorithmic parameter γ 1 . Due to the use of sub-step strategy, the implicit composite sub-step algorithm is undoubtedly more expensive than the traditional implicit method in the case of using the same integration step. As a result, in this case, the implicit composite

1410

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

sub-step algorithm can provide more accurate numerical solutions. In the last numerical simulation, the total computational cost of the composite sub-step algorithms can be significantly reduced by selecting the larger time step. Although the traditional implicit algorithm can use larger integration step to reduce the computational cost, larger time step is not conducive to nonlinear iterative convergence, which in turn requires more iterations to achieve satisfactory errors. As a result, the total computational time has not decreased or even increased. 8. Conclusion In this paper, a second-order accurate and L-stable direct integration algorithm is proposed for structural dynamical problems. The scheme has been formulated by employing the strategy of three sub-steps within a integration step to achieve some desired numerical properties, such as the second-order accuracy in time and the unconditional stability. In contrast to common composite sub-step algorithms, each sub-step in the proposed algorithm does not achieve the second-order accuracy. As an example, the first-order accurate backward Euler scheme is used in the first sub-step. But, the second-order accuracy in its final form can be obtained by designing the reasonable integration schemes in the second and third sub-step. The unique algorithmic parameter γ 1 is analyzed to achieve the unconditional stability, requiring that the γ 1 should satisfy γ 1 ∈ [0.181, 2.186]. The analysis of the amplitude decay and the period elongation error indicates that, compared with the Bathe algorithm, smaller period errors and less numerical dissipations in the low-frequency region can be provided by the proposed algorithm when the value of γ 1 lies in the small neighborhood of 0.181. On the other hand, when the γ 1 tends to 2.186, the proposed algorithm introduces more numerical dissipations in the low-frequency region, as well as the period elongation error. Besides, there is no overshoot for the proposed algorithm for any value of γ 1 . In the end, three numerical examples are solved to verify the theoretical analysis mentioned in this paper. Acknowledgments 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 and Li Xiangyang have led to the improvements of this paper; The authors gratefully acknowledge this assistance. Appendix A. The amplification matrix of the proposed algorithm The amplification matrix A of the proposed three sub-step algorithm is explicitly given as follows: First, applying the backward Euler scheme in the first sub-step to the SDOF system yields



xt+γ1 h hx˙ t+γ1 h h2 x¨t+γ1 h





= A10

xt hx˙ t h2 x¨t



(A1)

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



A10 =

1

β1

1 + 2γ1 ξ  −γ1 2 −2



γ1

0 0 0

1 −2ξ  − γ1 2

(A2)

where β1 = 1 + 2γ1 ξ  + γ12 2 . ω = /h and ξ are the undamped natural frequency and viscous damping ratio of the SDOF system, respectively. Second, the similar calculation can be done in the second sub-step. The matrix-vector form can be given as



xt+γ2 h hx˙ t+γ2 h h2 x¨t+γ2 h





xt+γ1 h hx˙ t+γ1 h h2 x¨t+γ1 h

= c2 A20





xt hx˙ t h2 x¨t

+ c1 A20

where the matrix A20 is



A20 =

1

β2

−2ξ  − c3



2

c 3 2

−1 −c3 2 c 3 ξ  + 2

0 0 0



(A3)

 (A4)

β2 = 2 + 2c3 ξ  + c32

(A5)

The coefficients ci , i = 1, 2, 3 are determined by Eq. (12). Then, the calculation in the third sub-step yields



xt+h hx˙ t+h h2 x¨t+h





= d3 A30

xt+γ2 h hx˙ t+γ2 h h2 x¨t+γ2 h





+ d2 A30

xt+γ1 h hx˙ t+γ1 h h2 x¨t+γ1 h



 + d1 A30

xt hx˙ t h2 x¨t

 (A6)

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

where the matrix A30 can be expressed as



A30 =

1

β3

−2ξ  − d4



2

d 4 2

−1 −d4 2 d 4 ξ  + 2

1411



0 0 0

(A7)

where β3 = 2 + 2d4 ξ  + d42 . The coefficients d1 , d2 , d3 and d4 are determined by Eqs. (22)–(25). Finally, substituting Eqs. (A1) and (A3) into Eq. (A6) yields



xt+h hx˙ t+h h2 x¨t+h





xt = A hx˙ t h2 x¨t



(A8)

where the amplification matrix A of the new sub-step algorithm is given as follows

A = d3 A30 (c2 A20 A10 + c1 A20 ) + d2 A30 A10 + d1 A30

(A9)

Hence, substituting Eqs. (A2), (A4), and (A7) into Eq. (A9) yields the explicit expression of the amplification matrix A. References [1] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Civil and Mechanical Engineering, Dover Publications, 20 0 0. [2] M. Géradin, D.J. Rixen, Mechanical Vibrations: Theory and Application to Structural Dynamics, 3, Wiley, 2015. [3] N.M. Newmark, A method of computation for structural dynamics, J. Eng. Mech. Div. 85 (3) (1959) 67–94. [4] 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 (1968). [5] J.C. Houbolt, A recurrence matrix solution for the dynamic response of elastic aircraft, J. Aeronaut. Sci. 17 (9) (1950) 540–550. [6] H.M. Hilber, T.J.R. Hughes, Collocation, dissipation and ‘overshoot’ for time integration schemes in structural dynamics, Earthq. Eng. Struct. Dyn. 6 (1) (1978) 99–117. [7] W.L. Wood, M. Bossak, O.C. Zienkiewicz, An alpha modification of Newmark’s method, Int. J. Numer. Methods Eng. 15 (10) (1980) 1562–1566. [8] H.M. Hilber, T.J.R. Hughes, R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthq. Eng. Struct. Dyn. 5 (3) (1977) 283–292. [9] J. Chung, G.M. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method, J. Appl. Mech. 60 (2) (1993) 371–375. [10] X. Zhou, K.K. Tamma, Design, analysis, and synthesis of generalized single step single solve and optimal algorithms for structural dynamics, Int. J. Numer. Methods Eng. 59 (5) (2004) 597–668. [11] J. Har, K.K. Tamma, Advances in Computational Dynamics of Particles, Materials and Structures, Wiley, 2012. [12] J.M. Benítez, F.J. Montáns, The value of numerical amplification matrices in time integration methods, Comput. Struct. 128 (5) (2013) 243–250. [13] K.J. Bathe, M.M.I. Baig, On a composite implicit time integration procedure for nonlinear dynamics, Comput. Struct. 83 (31–32) (2005) 2513–2524. [14] K.J. Bathe, Conserving energy and momentum in nonlinear dynamics: a simple implicit time integration scheme, Comput. Struct. 85 (78) (2007) 437–445. [15] S.Y. Chang, Explicit pseudodynamic algorithm with unconditional stability, J. Eng. Mech. 128 (9) (2002) 935–947. [16] J. Li, K. Yu, Noniterative integration algorithms with controllable numerical dissipations for structural dynamics, Int. J. Comput. Methods 15 (3) (2018) 52. [17] L. Wang, H. Zhong, A time finite element method for structural dynamics, Appl. Math. Model. 41 (2017) 445–461. [18] K. Hejranfar, K. Parseh, Numerical simulation of structural dynamics using a high-order compact finite-difference scheme, Appl. Math. Model. 40 (3) (2016) 2431–2453. [19] W. Kim, J.N. Reddy, Effective higher-order time integration algorithms for the analysis of linear structural dynamics, J. Appl. Mech. 84 (7) (2017) 071009. [20] W. Kim, J.N. Reddy, A new family of higher-order time integration algorithms for the analysis of structural dynamics, J. Appl. Mech.-ASME 84 (7) (2017) 071008–071017. [21] C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, 1971. [22] S. Dong, BDF-like methods for nonlinear dynamic analysis, J. Comput. Phys. 229 (8) (2010) 3019–3045. [23] J. Li, K. Yu, An alternative to the Bathe algorithm, Appl. Math. Model. 69 (2019) 255–272. [24] W.B. Wen, K. Wei, H.S. Lei, S.Y. Duan, D.N. Fang, A novel sub-step composite implicit time integration scheme for structural dynamics, Comput. Struct. 182 (2017) 176–186. [25] G. Noh, K.-J. Bathe, The Bathe time integration method with controllable spectral radius: The ρ ∞ -Bathe method, Comput. Struct. 212 (2019) 299–310. [26] M.M. Malakiyeh, S. Shojaee, K.-J. Bathe, The Bathe time integration method revisited for prescribing desired numerical dissipation, Comput. Struct. 212 (2019) 289–298. [27] G. Noh, K.-J. Bathe, An explicit time integration scheme for the analysis of wave propagations, Comput. Struct. 129 (2013) 178–193. [28] W. Kim, J.H. Lee, An improved explicit time integration method for linear and nonlinear structural dynamics, Comput. Struct. 206 (2018) 42–53. [29] W. Kim, A simple explicit single step time integration algorithm for structural dynamics, Int. J. Numer. Methods Eng. 119 (5) (2019) 383–403. [30] W. Kim, A new family of two-stage explicit time integration methods with dissipation control capability for structural dynamics, Eng. Struct. 195 (2019) 358–372. [31] J. Li, K. Yu, X. Li, A novel family of controllably dissipative composite integration algorithms for structural dynamic analysis, Nonlinear Dyn. 96 (4) (2019) 2475–2507. [32] W. Kim, J.N. Reddy, An improved time integration algorithm: A collocation time finite element approach, Int. J. Struct. Stab. Dyn. 17 (2) (2016) 1750024. [33] G. Noh, K.-J. Bathe, Further insights into an implicit time integration scheme for structural dynamics, Comput. Struct. 202 (2018) 15–24. [34] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, 3rd edition, Wiley, 2016. [35] H.M. Zhang, Y.F. Xing, Optimization of a class of composite method for structural dynamics, Comput. Struct. 202 (2018) 60–73. [36] K.C. Park, An improved stiffly stable method for direct integration of nonlinear structural dynamic equations, J. Appl. Mech. 42 (2) (1975) 464. [37] C. Kolay, J.M. Ricles, Development of a family of unconditionally stable explicit direct integration algorithms with controllable numerical energy dissipation, Earthq. Eng. Struct. Dyn. 43 (9) (2014) 1361–1380. [38] K. Yu, A new family of generalized-α time integration algorithms without overshoot for structural dynamics, Earthq. Eng. Struct. Dyn. 37 (12) (2008) 1389–1409.

1412

J. Li, K. Yu and H. He / Applied Mathematical Modelling 77 (2020) 1391–1412

[39] S. Mohammadzadeh, M. Ghassemieh, Y. Park, Structure-dependent improved Wilson-θ method with higher order of accuracy and controllable amplitude decay, Appl. Math. Model. 52 (2017) 417–436. [40] J. Li, K. Yu, X. Li, A generalized structure-dependent semi-explicit method for structural dynamics, J. Comput. Nonlinear Dyn. 13 (11) (2018) 111008–111020. [41] S.-B. Kwon, J.-M. Lee, A non-oscillatory time integration method for numerical simulation of stress wave propagations, Comput. Struct. 192 (Supplement C) (2017) 248–268.