1 Introduction

During the last several decades, a large number of researches have been conducted on the exploration of efficient time integration schemes for the linear and nonlinear analysis of structures^[1-4]. Generally, there are two major categories of time integration methods, namely, the explicit method^[5-7] and the implicit method^[8-10].

The implicit method is popular for a practical use because it can give a stable solution in the case of linear system, and the time increment can be easily selected. In the case of nonlinear analysis, it may be necessary for implicit schemes to update the stiffness matrix and conduct the iteration calculation within each time increment. Consequently, a lot of time is occupied for matrix inversion and iteration calculation. Meanwhile, a numerical divergence may be easily induced when a strong nonlinearity exists in a system^[11]. Therefore, various explicit schemes have been designed to overcome the above disadvantages of the implicit method. Compared with implicit methods, explicit methods need no effort for iteration calculation in each time step. Thus, fewer storage and computational efforts are required for explicit methods^[5]. The shortcoming of the explicit method is that almost all explicit methods are conditionally stable. Consequently, a small time increment and relatively large time steps may be required in a time history analysis, which, to a great extent, lowers the computation efficiency. However, a very small time increment can decrease the linearization errors for nonlinear systems. Actually, a small time increment can effectively capture intermediate and high frequency modes with desirable accuracy. Therefore, the explicit method is very suitable for the cases such as wave propagation and shock response problems where intermediate and high frequency modes are very significant for dynamic responses^[12-14].

Recently, Bathe and Noh^[9] developed a simple implicit time integration method. This method can remain unconditionally stable, without the use of adjustable parameters. Then, Noh and Bathe^[15] presented an explicit time integration scheme for the solution of wave propagation problems. The presented method possesses the second-order accuracy for different dynamic systems. A clear advantage of this method is that it is very effective in suppressing or filtering out spurious high-frequency modes. However, for the above two methods, each time step consists of two sub-steps, and thus requires much larger computation efforts than other conventional methods.

Recently, the B-spline method, as a powerful numerical tool, was used to develop a novel explicit time integration method by Rostami et al.^[16]. Compared with other conventional schemes, the presented B-spline method possesses higher computation efficiency. However, it is conditionally stable.

Wen et al.^[17-18] developed an explicit time integration method using septuple and quintic B-splines. With three adjustable algorithmic parameters, the presented methods can achieve both conditional and unconditional stability. Although the proposed method possesses higher computation accuracy, much more computation time is consumed, especially for the dynamic system with high degrees of freedom.

In this study, a new method is thus proposed by use of a family of uniform cubic B-spline functions. The paper is organized as follows. In Section 2, we introduce the explicit definition of B-spline and the expression of cubic B-spline interpolation. Then, the cubic B-spline interpolation is used to solve the equilibrium equation of linear single degree of freedom (SDOF) system in Section 3. Subsequently, the implementation of cubic B-spline interpolation on the multiple-degree of freedom (MDOF) system coupled with its calculation procedure is given in Section 4. The stability and accuracy analyses are conducted in Section 5. The validity of the proposed method is verified with four numerical examples in Section 6. Finally, in Section 7, we draw the conclusions of this study.

2 Overview of cubic B-splines 2.1 Explicit definition of B-spline function

There are many ways to define B-spline functions^{[17, 19]}. To simplify computer implementation, here we adopt an explicit definition to calculate B-spline functions^[19]. Thus, the ith B-spline basis function of the dth degree B_{i, d}(t), briefly denoted by B_{i, d}, can be defined by

(1)

(2)

where the knot span {t₀, t₁, …, t_i, t_i+1, …, t_d+i, t_d+i+1, …} must be a non-decreasing sequence of real numbers, that is, t_i≤ t_i+1. Because Eq. (2) can yield the quotient 0/0, we define this quotient to be zero. Without going into details, B_{i, d}(t) are linear combinations of two (d-1)-degree B-spline functions.

B_{i, d}(t) are non-negative and local supported over the interval [t_i, t_i+d+1]. As shown in Fig. 1, for cubic B-spline basis functions used in this study, the usable parameter range of B_{i, 3}(t) is [t_i, t_i+4]. Besides, we have within the usable interpolation domain [t_i, t_i+1]. Clearly, when t_i < t_i+1, B_{i, d}(t) are (d-1) times continuously differentiable.

2.2 Preliminary results of cubic B-spline interpolation

In this study, we use piecewise polynomials for the subsequent analysis. Let the time domain [a, b] be divided into n subintervals I_i[t_i, t_i+1] (i=0, 1, …, n-1) by a set of equidistant knots t_i=a+iΔt, where Δt=h= As shown in Fig. 1, to construct the cubic B-spline interpolation functions within the whole time domain [a, b], we should extend the set of knots as t_-3 < … < t_-1 < a=t₀ < t₁ < … < t_n=b < … < t_n+3.

For any knot t_j (j=-3, -2, …, t_n+3), we have t_j=t_i+(j-i)h. Thus, B_{i, 3}(t) (i=0, 1, …, n-1) can be simplified as

(3)

Equation (3) shows that B_{j, 3}(t) are the shifted instances of B_{j, 3}(t), that is, B_{j, 3}(t)=B_{i, 3}(t-(j-i)h). For any usable interpolation subinterval I_i [t_i, t_i+1] in Fig. 1, only four piecewise cubic B-spline functions are contributory to the B-spline interpolation.

In this study, the extended subintervals I'_i [t_i, t_i+φ] (1≤ φ < 2, i=0, 1, …, n-1) are used for the analysis. More specifically, the cubic B-spline interpolation function over any subinterval I'_i is used to approximate the exact solution within the subintervals I_i.

Further, let τ_i=(t-t_i)/(φΔt) (i.e., t=t_i+τ_iφΔt, 1 < φ < 2), the cubic B-spline functions in Fig. 2 and their lth derivatives with respect to the variable t can be expressed as

(4)

(5)

(6)

(7)

where l=0, 1, 2.

Thus, the cubic B-spline interpolation function within any subinterval I'_i can be defined by

(8)

where C_kⁱ are unknown real coefficients.

With Eqs. (4) -(7), Eq. (8) can be extended as

(9)

where

(10)

(11)

Clearly, for t=t_i, t=t_i+φ, and t=t_i+1, we have, respectively,

(12)

(13)

and

(14)

To solve B^(l)(0) and B^(l)(1), we tabulate B_{i, 3}^(l)(τ_i) (i=-3, -2, -1, 0, l=0, 1, 2) that evaluate at τ_i=0 (i.e., t=t_i) and τ_i=1 (i.e., t=t_i+φ) in Table 1.

3 Implementation of cubic B-spline on SDOF systems

For an SDOF system, a linear differential equation of motion can be written as

(15)

where ξ, ω, and f(t) are the damping ratio, the undamped circular natural frequency, and the modal forcing excitation of the system, respectively. u(t), u⁽¹⁾ (t), and u⁽²⁾ (t) are the displacement, velocity, and acceleration functions, respectively.

The initial value problem is to solve Eq. (15) to meet the given initial conditions of u(t₀) =a₀, u⁽¹⁾ (t₀) =a₁, and u⁽²⁾ (t₀) =a₂.

Let t₀, t₁, t₂, …, t_n-1 be n time knots over the whole time interval [0, t_f] so that t_i=iΔt, i=0, 1, 2, …, n-1, t₀=0, and t_n-1=t_f. Using Eqs. (8) and (9), the approximate solution to Eq. (15) and its lth derivatives, within any subinterval I'_i (i=0, 1, 2, …, n-1), can be expressed as

(16)

where the subscript i of u_i(t) denotes the number of time subinterval.

To simplify the subsequent deduction, we let

(17)

(18)

Substituting t=t_i (i.e., τ_i=0) and t=t_i+φ (i.e., τ_i=1) into Eq. (16) begets six equations. With Eqs. (17) and (18) and Table 1, these six equations can be expressed in the explicit forms of

(19a)

(19b)

(19c)

and

(20a)

(20b)

(20c)

Using Eqs. (19a)-(19c), C_kⁱ (k=-3, -2, -1, i=0, 1, 2, …, n-1) can be calculated by

(21a)

(21b)

(21c)

Prescribing that Eq. (15) is satisfied at t=t_i+φ (i.e., τ_i=1) renders

(22)

With Eq. (17), Eq. (22) can be written as

(23)

Substituting Eqs. (20a)-(20c) into Eq. (23) yields

(24)

where α, β, and γ are given by

(25a)

(25b)

(25c)

C_-3ⁱ, C_-2ⁱ, and C_-1ⁱ can be previously solved by Eqs. (21a)-(21c). Therefore, with Eq. (24), we can solve C₀ⁱ as

(26)

After solving C_-3ⁱ, C_-2ⁱ, C_-1ⁱ, and C₀ⁱ, we can use Eqs. (16) -(18) to determine as

(27)

After solving , we need to substitute into Eqs. (21a)-(21c) to proceed the whole recurrence calculation. Obviously, we have

After substituting into Eq. (17), we can obtain all needed numerical displacements (i.e., u_i(t_i)), velocities (i.e., u_i⁽¹⁾ (t_i)), and accelerations (i.e., u_i⁽²⁾ (t_i)).

4 Implementation of cubic B-spline on MDOF systems

The dynamic equilibrium equations governing a linear MDOF system can be expressed as

(28)

where M, D, and K are the mass, damping, and stiffness matrices, respectively. F is the vector of externally applied loads. U(t), U⁽¹⁾ (t), and U⁽²⁾ (t) are the unknown displacement, velocity, and acceleration function vectors of the discrete nodes on the structure.

The initial conditions are U(t₀) =a₀, U⁽¹⁾ (t₀) =a₁, and U⁽²⁾ (t₀) =a₂.

As previously discussed in Section 3 for the SDOF system, Eq. (16) is used to represent the approximate solution to Eq. (28). First, we let

(29)

Assume that u_{j, i}(t) are the cubic B-spline interpolation functions of the exact solution u_j(t) within any time interval I'_i (i=0, 1, 2, …, n-1). Then, by use of Eq. (16), we have

(30)

Substituting Eq. (30) into Eq. (29) gives the approximate function vectors of the exact solution U^(l)(t) within any time interval I'_i as

(31)

where C_kⁱ are the unknown coefficient vectors of size N×1. C_kⁱ can be written as

(32)

Similar to the SDOF system, here, we let

(33)

By substituting t=t_i (i.e., τ_i=0) and t=t_i+φ (i.e., τ_i=1) into Eq. (31) and using Eq. (33) and Table 1, we have

(34a)

(34b)

(34c)

and

(35a)

(35b)

(35c)

With Eqs. (34a)-(34c), C_-3ⁱC_-2ⁱ, and C_-1ⁱ can be calculated by

(36a)

(36b)

(36c)

Prescribing that Eq. (28) is satisfied at t=t_i+φ (i.e., τ_i=1) yields

(37)

Using Eq. (33), Eq. (37) can be further expressed as

(38)

Substituting Eqs. (35a)-(35c) into Eq. (38) yields

(39)

where

(40)

(41)

(42)

C_-3ⁱ, C_-2ⁱ, and C_-1ⁱ have been previously solved by Eqs. (36a)-(36c). Thus, we can solve C₀ⁱ by

(43)

Then, substituting t=t_i+1 into Eqs. (31) and (33) yields

(44)

In Eq. (44), to complete the recurrence procedure, we let . Meanwhile, after substituting into Eq. (33), we can obtain (l=0, 1, 2, i=0, 1, 2, …, n-1).

To clarify the complete algorithm of the proposed method, we give the calculation procedure of the MDOF system from the perspective of computer programming in Table 2.

Table 2 illustrates that all the C_kⁱ depend on the previous values already calculated. Therefore, the proposed algorithm is an explicit one. From Table 2, we can see that only one symmetric matrix inversion needs to be conducted within each time step. Moreover, compared with the methods of Refs. [9] and [15], all intermediate variables C_kⁱ have simpler forms and can be easily solved, which, to a great extent, can save computation time greatly. Thus, theoretically, the proposed method has the more desirable computation speed than the methods of Refs. [9] and [15].

5 Numerical stability and accuracy analyses 5.1 Stability analysis

In general, any global equation of motion can be degraded into a set of uncoupled SDOF systems by use of the modal decomposition. Meanwhile, the integration of the uncoupled equations is equivalent to the integration of the original system. Therefore, to study the stability properties of the proposed method, it suffices to only consider the SDOF system.

The proposed method transfers the state at the ith step to the (i+1) th step as elucidated in Section 3. Clearly, the recurrence formula of the proposed algorithm for any SDOF system can be intrinsically written as

(45)

where A is the amplification matrix, and are given by

(46)

(47)

In addition, we need to define as

(48)

To find the amplification matrix A, at first, we need to construct the relationships among and by use of coefficients C_kⁱ (k=-3, -2, -1, 0).

If we solve Eqs. (20a)-(20c) to reach C_-2ⁱ, C_-1ⁱ, and C₀ⁱ, we would have

(49a)

(49b)

(49c)

Similarly, with Eqs. (21a)-(21c), C_-3ⁱ, C_-2ⁱ, and C_-1ⁱ can be represented by

(50a)

(50b)

(50c)

Combining Eqs. (49a)-(49c) with Eqs. (50a)-(50c) gives

(51)

In addition, substituting Eq. (49c) and Eqs. (50b)-(50c) into Eq. (24) renders

(52)

where

(53)

Actually, Eq. (52) is the equivalent form of Eq. (22). Equations (51) and (52) can be expressed in a unified form of

(54)

where

(55)

With Eq. (27), can be written as

(56)

where

(57)

(58)

With Eqs. (49a)-(49c) and Eqs. (50a)-(50c), we have

(59)

(60)

where

(61)

Then, after substituting Eqs. (59) and (60) into Eq. (56) and using Eq. (54), we have

(62)

In Eq. (62), the amplification matrix A (i.e., ) is used to investigate stability of the proposed algorithm. Because the amplification matrix A contains too many terms, the concrete expression of A is not presented here. In fact, we can easily obtain A by virtue of some general mathematical softwares.

In this study, the stability analysis is conducted by solving the eigenvalue problem of the amplification matrix A. The eigenvalues of A are obtained by solving |A-λ I|=0, where I is the identity matrix of size 3× 3. Here, let λ_i (i=1, 2, 3) be the eigenvalues of A. To obtain a convergent algorithm, the norm of all eigenvalues should be less than unity, and thus we should have where ρ(A) is called the spectral radius.

For the proposed step-by-step method, it suffices to investigate ρ(A) only at ξ=0. At this moment, ρ(A) is a function only in terms of ω Δt and the algorithmic parameter φ. The natural frequency ω satisfies ω=2π/T, where T is the period of the considered system. To ensure that the algorithm in question is stable, the effects of the parameter φ on the spectral radius ρ(A) need to be investigated. Here, we conduct this study through numerical experimentation instead of rigorous mathematical operation.

The variation of the spectral radius as a function of Δt/T is shown in Fig. 3. For the convenience of comparison, only the spectral radius curves of some representative parameter values are shown in Fig. 3, where the spectral radius curves move downward with the increase in the parameter φ. Clearly, when φ≥ 1.37, the spectral radius achieves unconditional stability. However, increasing the value of φ would undoubtedly lower interpolation accuracy of the proposed method. Therefore, we give the suggested parameter condition as 1.0≤φ≤ 1.4. Furthermore, we provide five parameter cases of the cubic B-spline method as

Case Ⅰ φ=1.00, Δt/T≤0.550;

Case Ⅱ φ=1.20, Δt/T≤0.764;

Case Ⅲ φ=1.30, Δt/T≤1.175;

Case Ⅳ φ=1.37 (unconditionally stable);

Case Ⅴ φ=1.40 (unconditionally stable).

In Fig. 4, the radii from the cubic B-spline method are especially compared with the radii from the Noh-Bathe method^[15]. Clearly, the cubic B-spline method is more flexible in controlling the stability and stable region than the Noh-Bathe method. In Fig. 4, for the Noh-Bathe method, when the algorithmic parameter p is equal to 0.5, the spectral radius remains to be 1 for a slightly wider range than the one from any other given case. However, compared with the central difference method schemes, CaseⅠ has a much wider stable region. The proposed scheme gives very similar spectral radii as the Wilson-θ method, but possesses different calculation procedure which is formulated based on the uniform B-spline interpolation.

5.2 Accuracy analysis

In general, a standard procedure used to estimate the numerical error of a step-by-step time integration method is to quantify the difference between numerical displacements and theoretical displacements for the free vibration of an undamped SDOF system. For the free vibration problem, there are two definitions for error estimation: amplitude decay (A_D) and period elongation (P_E). A_D is sometimes reported as the "algorithmic damping ratio".

In Fig. 5, the algorithmic damping ratio and period elongation curves from the cubic B-spline method are compared with the curves from the Noh-Bathe method^[15] and the Bathe-Noh method^[9]. Figure 5(a) illustrates that the algorithmic damping ratio from the cubic B-spline method increases with the increase in the parameter φ. Figure 5(b) shows that the period elongation from the cubic B-spline method decreases with the increase in the parameter φ. Clearly, the Noh-Bathe method of Ref. [15] produces the smaller damping ratio and period elongation than the proposed method. Figure 5(b) shows that the Bathe-Noh method of Ref. [9] produces the larger period elongation than the proposed method. However, considering that all curves in Fig. 5 are close, the proposed method has the desirable algorithmic damping ratio and period elongation.

6 Numerical examples

To demonstrate validity of the proposed method for dynamic analysis, four numerical examples have been tested in this section, and five parameter cases previously given in Subse-ction 5.1 are used for computation. In addition, for comparison, the generalized-α (α_f= method^[2], the Bathe-Noh method^[9], and the Noh-Bathe method^[15] are used in this section. Here, it is important to note that, for the Noh-Bathe method, we adopt the suggested algorithmic parameter value p=0.54 for computation^[15].

6.1 Standard undamped system

To test the calculation accuracy of the proposed scheme, the free vibration of a standard undamped SDOF system is considered. This system is

(63)

where ω=π, that is, T=2.

A time duration of t=2 s is considered and calculated for accuracy and efficiency tests, and the global error norms E_l (l=0, 1, 2) are employed and defined as

(64)

in which u_i^(l) are numerical results at the time t_i, and are the corresponding exact ones.

In Fig. 6, global errors of various schemes are plotted in a log form. It can be noted that accuracy of the proposed scheme decreases with the increase in the parameter φ. CaseⅤ shows acceptable accuracy compared with other considered methods.

6.2 SDOF system

An SDOF system for numerical simulation is given by

(65)

whose exact solution is u(t)= e^-2t(cos t+2 sin t)-(8 cos (2t)-sin (2t))/65.

To investigate the quantified effects of ratio Δt/T on calculation accuracy, we select four different time increments Δt for calculation and accordingly list the global errors thereof in Table 3. Here, the global errors are defined by

(66)

where u_num⁽ⁱ⁾ and u_exact⁽ⁱ⁾ are the numerical result and the exact solution at some given time, respectively. n_g is the number of time steps. n_g is calculated by n_g=t_a/Δt, where t_a is the selected time duration. Here, we select t_a=4 s for computation.

Table 3 shows that the calculation accuracy of the proposed method decreases with the increase in the parameter φ. Clearly, CaseⅤ gives roughly the same magnitude of displacements and velocities as the generalized-α method^[2]. Nevertheless, the acceleration errors from CaseⅤ are almost one order of magnitude less than the generalized-α method. All the global errors from CaseⅠ and CaseⅡ are much smaller than those from the generalized-α method.

In Table 3, compared with the Noh-Bathe method^[15], CaseⅠ produces slightly larger displacement errors, but gives clearly smaller velocity and acceleration errors. Table 3 shows that the Bathe method^[9] possesses higher accuracy than the cubic B-spline method. However, note that the Bathe method adopts two sub-steps in each time step, and the computation time consumed by these two methods is definitely larger than other conventional methods. Thus, to conduct the computation efficiency analysis, we give the time consumption results of various methods in Table 4, where Δt=0.05 s and three different values of n_g are considered for the time consumption analysis. Table 4 shows that the proposed method spends almost the same computation time as the generalized-α method. However, for the Bathe-Noh method^[9] and Noh-Bathe method^[15], much more computation time than the proposed method is consumed.

In Table 5, to quantify the efficiency differences between the proposed method and two given novel methods, various Δt are adopted for different methods. In the table, we define the relative errors for numerical results in the given time as

(67)

Table 5 shows that the cubic B-spline method could give almost the same magnitude of accuracy when roughly the same computation time as the Noh-Bathe method^[15] is consumed. Tables 3-5 demonstrate that, compared with the Bathe-Noh method^[9] (unconditionally stable) and Noh-Bathe method^[15] (conditionally stable), the computation efficiency of the proposed method is acceptable.

6.3 Two-dimensional Howe truss under impact loads

A Howe truss under four concentrated impact loads is shown in Fig. 7. Material properties for all elements are shown in the figure. For comparison, the generalized-α method is used for computation. The least period of this system (T_min) is equal to 0.008 2 s. Thus, we select Δt=2× 10^-3 s (Δt≤0.550 T_min) for CaseⅠ, and Δt=8× 10^-3s(Δt≤1.175 T_min) for Case Ⅲ.

The time-history of vertical displacement over a time interval between 0 s and 2 s of Node 5 is illustrated in Fig. 8, and the horizontal displacement of Node 13 is plotted in Fig. 8. It is easily observed that the displacement curves from two conditional B-spline cases match very well with the curves from the generalized-α method^[2]. As for Case Ⅴ, we can select a relatively large time increment Δt for computation as shown in Figs. 8 and 9, and thus the numerical results from CaseⅠ are definitely less accurate than the other given B-spline cases. However, Figs. 8linebreak and 9 show that all cases of cubic B-spline method are valid for dynamic responses.

As a sample, the results of time consumption for a complete analysis of Subsection 6.2 are tabulated in Table 6, where all methods use Δt=4×10^-4 s for computation. For different time steps n_g, the proposed method consumes almost the same time as the generalized-α method. By contrast, the Bathe-Noh method^[9] and Noh-Bathe method^[15] consume much more time than the proposed method. Furthermore, from Tables 4 and 6, we can safely conclude that the cubic B-spline method has high computation efficiency for the SDOF and MDOF systems.

6.4 Forced vibration of simply supported uniform continuous beam

To verify validity of the proposed method for the finite element method (FEM) calculation, a simply supported uniform continuous beam under a lateral concentrated variable load at the middle point is shown in Fig. 10. The dimensions and parameters used for analysis are the overall length L=8 m, the radius of circular section R=2×10^-2 m, the cross-sectional area A=π R², the sectional inertia moment I=π R⁴/2, Young's modulus E=100 GPa, Poisson's ratio μ=0.3, the material density of the beam ρ=4×10^-4 kg/m ³, the damping ratio ξ=0, and the lateral concentrated variable load . Here, we select F₀=1 kN, and ω₀=4 rad/s. The analytical solution of the beam's bending deflection can be written as

(68)

where the natural frequency ω_r is determined by

The equilibrium equation of this problem is formulated by use of the cubic Hermite finite element with its element number N_s=8. The boundary conditions are satisfied by directly setting w=0 at x=0 and x=L. The least period of the equilibrium equation is equal to 0.005 6 s. Thus, similar to Subsection 6.2, we select Δt=3×10^-4 s for Case Ⅰ, Δt=5×10^-3 s for Case Ⅲ, and Δt=5×10^-2 s for Case Ⅳ. For simplicity, Fig. 11 just provides response values at the maximum deflection position (i.e., the middle point). In Fig. 11, although the various time increments Δt are used for different cubic B-spline cases, all curves from the proposed method are in desirable agreement with the exact curves.

7 Conclusions

In this study, a new time integration method is proposed for the structural dynamic analysis by use of the cubic B-spline interpolation. By solving the differential equations of motion, a simple algorithm is formulated to calculate the dynamic responses of SDOF and MDOF systems.

The stability analysis shows that the proposed method can achieve both unconditional and conditional stability by adjusting the algorithmic parameter. The corresponding parameter conditions for conditional stability have been given in this study. The proposed method gives the acceptable numerical damping ratio and period elongation.

The validity of the proposed method has been confirmed with four numerical simulations. Time consumption results from numerical simulations demonstrate that, compared with other well-known methods, the proposed method has an appropriate computation speed. The computation efficiency analysis shows that the proposed method still possesses desirable computation efficiency compared with other schemes. Thus, the proposed method is a good choice for the structural dynamic analysis due to its high computation efficiency.

Actually, in terms of the capability of direct integration method for dealing with the nonlinear dynamic analysis, this method has potential ability to be generalized for a nonlinear system. This can be a good topic for the further research.