Development of Non-Dissipative Direct Time Integration Method for Structural Dynamics Application

Sun-Beom Kwonand Jae-Myung Lee,

Abstract: A direct time integration scheme based on Gauss-Legendre quadrature is proposed to solve problems in linear structural dynamics. The proposed method is a oneparameter non-dissipative scheme. Improved stability, accuracy, and dispersion characteristics are achieved using appropriate values of the parameter. The proposed scheme has second-order accuracy with and without physical damping. Moreover, its stability, accuracy, and dispersion are analyzed. In addition, its performance is demonstrated by the two-dimensional scalar wave problem, the single-degree-of-freedom problem, two degrees-of-freedom spring system, and beam with boundary constraints.The wave propagation problem is solved in the high frequency wave regime to demonstrate the advantage of the proposed scheme. When the proposed scheme is applied to solve the wave problem, more accurate solutions than those of other methods are obtained by using the appropriate value of the parameter. For the single-degree-offreedom system, two degrees-of-freedom system, and the time responses of beam, the proposed scheme can be used effectively owing to its high accuracy and lower computational cost.

Keywords: Structural dynamics, finite elements, direct time integration, Gauss-Legendre quadrature, non-dissipative scheme.

1 Introduction

In the finite element method, there are various sources of errors [Bathe (1996)], such as discretization [Bohinc, Brank and Ibrahimbegovics (2014); Huang and Griffiths (2015);Lee and Cangellaris (1992); Mohite and Upadgyay (2015); Payen and Bathe (2012)],numerical integration in space [Ham and Bathe (2012); Idesman, Schmidt and Foley(2011); Ja?kowiec and Pluciński (2017); Lee and Cangellaris (1992)], evaluation of constitutive relations [Lee, Kim, Park et al. (2015); Lee, Chun, Kim et al. (2009); Lee,Kim, Park et al. (2016); Park, Lee, Chun et al. (2011); Yoo, Lee, Park et al. (2011)],solution of dynamic equilibrium equations [Bathe and Baig (2005); Idesman (2011);Kwon and Lee (2017); Newmark (1959); Noh and Bathe (2013); Noh, Ham and Bathe(2013); Rougier, Munjiza and John (2004); Wen, Duan, Yan et al. (2017)], solution of finite element equations by iteration [Bogaers, Kok, Reddy et al. (2016); Golub and van Loan (1996); Varga (1962); Xu and Prozzi (2014)], and round-off [Alvarez-Aramberri,Pardo, Paszynski et al. (2012); Fried (1986); Vignes (1988)]. The finite element method can provide accurate solutions of dynamic equilibrium equations when an appropriate direct time integration method is used. Time integration schemes can be classified into two categories: explicit and implicit.

Direct time integration is primarily used for solving wave propagation and structural dynamic problems. When a wave propagation problem is solved using the time integration scheme, a number of spurious oscillations occur. Dissipative schemes,filtering, or damping schemes are used in order to reduce these oscillations. By filtering specific time and space points, errors from numerical oscillations can be reduced[Holems and Belytschko (1976); Idesman, Samajder, Aulisa et al. (2009); Idesman,Schmidt and Foley (2011)]. However, these approaches are not suitable for analyses involving global solutions [Noh, Ham and Bathe (2013)]. Damping schemes are methods of adding a viscous pressure term to the dynamic equilibrium equation [Benson (1992);Johnson and Beissel (2001)]. Their disadvantage is that the artificial viscosity depends on the problem. In the case of dissipative schemes, as there are discarded frequency modes owing to the numerical dissipation property, there are few spurious oscillations in the solutions of wave propagation problems. However, numerical errors often occur due to numerical dissipations [Bathe (1996); Noh and Bathe (2013); Noh, Ham and Bathe(2013); Payen and Bathe (2012)]. Therefore, in structural dynamics, errors can be reduced using non-dissipative schemes.

In non-dissipative schemes, the central difference method (CDM) and the trapezoidal rule(Newmark method with β =0.25 and γ =0.5) with second-order accuracy are still widely used in explicit and implicit schemes, respectively. In these schemes, dispersion errors occur in high-frequency modes due to non-dissipation. Therefore, in wave propagation problems, non-dissipative schemes result in spurious oscillations. However,in structural dynamics, CDM and trapezoidal rule are widely used owing to their computational simplicity and lower computational cost.

In this study, a new time integration scheme is proposed based on Gauss-Legendre quadrature. Gauss-Legendre quadrature is a numerical integration formula. It is defined as follows:

To use this formula, the node point xiand the corresponding weight wishould be determined. The node points and weights for several values of n have been obtained in previous research [Abramowitz and Stegun (1972); Cheney and Kincaid (2012); Stroud and Secrest (1996)]. The Gauss-Legendre quadrature formula yields reasonably accurate solutions with fewer function evaluations [Cheney and Kincaid (2012)]. The Gauss-Legendre quadrature is commonly used to calculate integrals in isoparametric finite element analysis [Mizusawa and Takami (1992); Ray, Dong and Atluri (2016)]. The Gauss-Legendre quadrature is generally more efficient, for it provides more accurate integration results for the same number of evaluations [Bathe (1996)]. However, the domain of integration is taken as [-1, +1]; thus, the time integral should be transformed.After translation, the number of node points n+1 is determined according to the required accuracy. Then, the Gauss-Legendre rule with n+1 node points is used for the time integration.

As the form of the time integration scheme depends on n, the proposed scheme is considered to be a one-parameter scheme. The parameter controls stability, accuracy, and dispersion, and its effect is analyzed in the present study. Finally, solutions are provided for the wave propagation problem and some examples in structural dynamics.

2 Time integration scheme based on Gauss-Legendre quadrature

2.1 Proposed implicit time integration scheme

The governing finite element equation for linear structures is as follows [Coppolino(2016)]:

The displacement at time t+Δt is computed as follows:

where the domain of integration is not [-1, 1] but [t, t+Δt]. Thus, in order to use Gauss-Legendre quadrature, the integral is transformed as follows:

Then, the integral in Eq. (4) is approximated by Gauss-Legendre quadrature, that is,

where p and q are parameters to be determined. Difficulty arises in the computation of u(t+ Δt) because the velocity vectors are calculated at half-step points. The best known solution to this problem is to linearly interpolate u(t) as follows [Park and Underwood(1980)]:

In order for the proposed scheme to have second-order accuracy with and without damping, the following relation holds:

To demonstrate second-order accuracy in time for the proposed scheme when Eq. (8)holds, local truncation error analyses are explained in Appendix A.

In the proposed method, the implicit and explicit schemes are considered for the displacement vector u and the velocity vector u, respectively. The proposed method for p≠1 is an implicit time integration scheme as a matrix system is solved when the proposed method is used, whereas the proposed method (p=1) is an explicit scheme as the effective stiffness matrix can be seen as the effective mass matrix owing a4=0 (see Tab. 1).

Table 1: Computation procedure of the proposed scheme

2.2 Stability analysis

Using the equations in Tab. 1, in the modal equations the recursive relationship of the proposed scheme is expressed as follows:

where A and L are the integration approximation and load operators, respectively,which are given in Appendix B. These operators depend on the direct time integration scheme. Eq. (9) can be used to provide the solution at any time t+nΔt, thereby providing stability and accuracy characteristics.

In the undamped case, the characteristic polynomial of A becomes as follows:

where

and Ω0= ω0Δt. Here, ω0is the modal natural frequency. Therefore, the eigenvalues of A are as follows:

When the spectral radius of the approximation operator is smaller than or equal to 1, the direct time integration is stable. If the stability criterion is met, Anis bounded for n→ ∞. Moreover, if ρ(A)＜ 1, Anis getting closer to0, decreasing more rapidly when ρ(A) is small [Bathe (1996)]. The spectral radius of A is defined as follows:

The bifurcation point Ωb, where the eigenvalues change from complex to real, i.e.,-A2=0, is given by following equation:

When Ω0＜ Ωb, the eigenvalues become complex and the spectral radius of A is always 1, whereas when the eigenvalues become real, the spectral radius of A is always lager than 1 (obtained with some arithmetic). Therefore, for the proposed scheme, the stability limit Ωsis the same as the bifurcation point Ωband is expressed as follows:

However, if p≤0.5, the denominator of right-hand side of Eq. (15) is zero or imaginary.When p is smaller than or equal to 0.5, the eigenvalues become complex and then the spectral radius of A is always 1.

The spectral radius at the bifurcation point is denoted by ρb. For the proposed scheme,ρbis 1 for any parameter p. This implies that the proposed scheme is non-dissipative.

As shown in Fig. 1, the parameter p plays an important role in stabilization behavior.When p≤0.5, the proposed scheme has unconditional stability, as in the case of general implicit methods. However, when p＞0.5, the parameter p affects the stability limit.Therefore, the proposed scheme is not unconditionally stable when p＞0.5.

Fig. 2 shows the spectral radius of the proposed scheme for various values of p, and other schemes. As the equations of the proposed method (p=1) and the half-step CDM are identical, the two methods yield the same results. When p≤0.5, the proposed scheme with non-dissipative properties has unconditional stability. Therefore, the proposed scheme yields the same results as the trapezoidal rule. However, the Noh-Bathe explicit method, the Bathe implicit method, and three sub-step method [Wen, Wei, Lei et al. (2017)] have numerical dissipation. Therefore, these methods are non-dissipative(resp., dissipative) for low (resp., high) frequency modes. For the Noh-Bathe method, the spectral radius of the approximation operation rapidly decreases for Ω0greater than,approximately, 1.25 (Δt/T0≈ 0.2) [Noh and Bathe (2013)], whereas for the Bathe method, ρ(A)≈ 1 for Ω0less than, approximately, 1.88 (Δt/T0≈ 0.3) [Noh, Ham and Bathe (2013)].

Figure 1: Stability limit line of the proposed scheme for various values of p

Figure 2: Spectral radius of the proposed scheme for various values of p, and other schemes considering no damping. Results for the proposed scheme with p=1 and CDM are identical, whereas results for the proposed scheme with p≤0.5, trapezoidal rule, generalized composite scheme are identical. The generalized composite scheme is described in Kim et al. [Kim and Choi (2018)] and the three sub-step scheme is described in Wen et al. [Wen, Wei, Lei et al. (2017)]

2.3 Accuracy analysis

While the stability analysis is employed to find the stability limit and discarded wave modes caused by numerical damping, the accuracy of a time integration method is not evaluated by conducting that same stability analysis. Therefore, an accuracy analysis is additionally carried out.

The accuracy of the proposed scheme is investigated by solving the following simple initial value problem [Bathe (1996)]:

The exact solution of Eq. (17) is x=cosωt. Although the damping ratio ξ should be considered, the significant characteristics of the solution can be investigated by solving the undamped problem [Bathe (1996)]. Figs. 3-4 show the period elongations (PE) and amplitude decays (AD) for the proposed and other schemes, respectively. PE and AD are defined as PE=(T -T)/T and AD=-2πξ (see Fig. 5), respectively. Here T, T ,and ξ are the exact period, numerical period, and algorithmic damping ratio,respectively. Please see Bathe [Bathe (1996)] regarding details of the accuracy analysis.In terms of PE, the proposed scheme with p=1, p=5/6, p=4/6, and p=3/6 has similar accuracy to that of CDM, the Noh-Bathe scheme, the Bathe implicit method,and the trapezoidal rule, respectively. However, in terms of AD, except for p=1, the proposed scheme has smaller solution errors as compared with those of CDM, Noh-Bathe method, and Bathe method.

Figure 3: Percentage PE of the proposed scheme for various values of p and other schemes. Results for the generalized composite scheme and three sub-step scheme are identical

For the proposed scheme, PE decreases as p increases. When p＜5/6, PE is negative,whereas it is positive when p≥5/6. AD is negative if p＜3/6, whereas it is positive elsewhere. When 3/6≤p＜4/6 and p≥5/6, AD increases with p, whereas it increases as p decreases when 4/6≤p＜5/6. Therefore, the proposed scheme(p=5/6) yields more accurate solutions. As shown in Fig. 3, the proposed scheme(p=5/6) has similar PE to that of the Noh-Bathe method. However, since the Noh-Bathe method consists of two sub-steps [Noh and Bathe (2013)], it is more timeconsuming than the proposed scheme (see Fig. 21). The proposed scheme yields highly accurate solutions and requires less time than the Noh-Bathe method.

Figure 4: Percentage AD of the proposed scheme for various values of p and other schemes. Results for the proposed scheme with p=5/6 and the generalized composite scheme are identical

Figure 5: Definition of PE and AD

2.4 Dispersion analysis

In this section, the dispersion errors are analyzed using the proposed scheme and other schemes for a two-dimensional wave propagation problem. A regular mesh of 4-node elements is used for spatial discretization.

The scalar wave governing equation is

where u is the field variable, ?2is the Laplace operator, and c0is the exact wave velocity. The corresponding finite element discretization yields

where M and K are the global mass and stiffness matrices, respectively. The matrices M and K can be obtained by summing the matrices M(k)and K(k)[Bathe (1996)],which are the local mass and stiffness matrices of the element (k), respectively, with volume V(k):

where H(k)is the displacement shape (or interpolation) matrix of the element (k).

Using the equations in Tab. 1, the linear multistep form of the proposed scheme is expressed as follows:

where γ =CFL2?2, CFL is the CFL numberand ? is the length of the side of the elements.

In the two-dimensional analysis, the general and numerical solutions of Eq. (18) areand, respectively,

where ω0, k0= ω0/c0, ω, k= ω/c, c, and θ are the frequency of the wave mode,corresponding wave number, numerical frequency, numerical wave number, numerical wave speed, and propagating angle measured from the x-axis, respectively. The numerical wave speed c and exact wave speed c0are not the same. The dispersion error is observed by their difference.

Using a regular mesh of 4-node elements with nodes equally spaced at a distance ? from each other along the x and y axes (Δx= Δy= ?) for the spatial discretization, the solution of the wave equation at time ntΔt and location nx?,ny? is obtained using the following equation [Noh and Bathe (2013); Noh, Ham and Bathe (2013)]:

For a 4-node element, the rows of Mc, Ml, and K are obtained as follows [Kwon and Lee (2017)]:

where Mcand Mlare the consistent and lumped mass matrices, respectively. It is natural to use consistent and lumped mass matrices for an implicit and explicit method,respectively. In implicit time integration, when a lumped mass matrix is used, there are no computational and storage advantages because accelerations are computed asAsis not diagonal (see Tab. 1) even if M is a lumped mass matrix, a system of equations must be solved. On the other hand, for explicit methods, a lumped mass matrix can offer significant computational advantages for calculations because effective mass matrixis used instead of effective stiffness matrixto calculate accelerations[Bathe (1996)]. However, if M is a consistent mass matrix, accelerations are obtained by solving a system of equations. Moreover, for explicit methods, the decomposition of a consistent mass matrix is generally not affordable because an implicit scheme with a much larger time step may be used at the same computational cost [Huang, Kamenski and Lang (2016)]. Therefore, for the CDM and the Noh-Bathe method, the lumped mass matrix is considered, whereas for the proposed scheme, trapezoidal rule and Bathe method, the considered mass matrix is consistent. Although the proposed scheme(0.5＜p＜1) has conditional stability, a matrix system needs to be solved. Therefore, it is natural to use the consistent mass matrix. However, as the proposed scheme (p=1) is considered to be the explicit method, it is natural to use the lumped mass matrix.

The error curves of the Noh-Bathe method [Noh and Bathe (2013)] and Bathe method[Noh, Ham and Bathe (2013)] are shown with the optimal CFL values in Figs. 6 and 7,respectively. Because the Noh-Bathe explicit method and the Bathe method are dissipative, there are discarded wave modes due to numerical damping in the numerical solutions. In the Noh-Bathe method and Bathe method, wave modes with Ω0＞ 0.6π are discarded [Noh and Bathe (2013); Noh, Ham and Bathe (2013)]. Using the definition of CFL number, Δx/(λ/2) in terms of Ω0is expressed as follows:

When we assume that the numerical and exact wave speeds are approximately the same,the discarded wave modes of the Noh-Bathe method with CFL=1.85 and Bathe method with CFL=1 are the wave modes with Δx/(λ/2)＞ 0.32 and Δx/(λ/2)＞0.6 using Eq. (28), respectively. As the high wave modes with larger dispersion errors are discarded in dissipative schemes, the dispersion errors of the Noh-Bathe method and Bathe method are smaller than those of the non-dissipative schemes (see Figs. 6 and 7).However, this property is not desirable for the high frequency wave modes as high wave modes do not participate in the total solution.

Figure 6: Relative wave speed errors of the Noh-Bathe method (CFL=1.85) for various values of the propagation angle with lumped mass matrix (dotted lines: discarded wave modes)

In the proposed method, CDM, and the trapezoidal rule, all wave modes participate in the total solution as these schemes do not have numerical damping characteristics (these schemes are non-dissipative schemes). When the high frequency wave propagation problems are solved, it is necessary to calculate all the wave modes. Therefore, nondissipative time integration methods might produce an accurate solution of high frequency wave propagation problem.

For the proposed scheme, the important point to note here is that as the linear multistep forms of the proposed scheme with p=1 and p=3/6 are the same as in CDM [Noh and Bathe (2013)] and the trapezoidal rule [Noh, Ham and Bathe (2013)], respectively,the relative wave speed errors are the same (see Figs. 8 and 9). The proposed scheme with p=1 provides no dispersion error in the one-dimensional case (θ=0). However,in actual analyses, the waves will propagate in all directions, and then the relative wave speeds are important considering the propagation angles. The maximum dispersion error of the proposed scheme with p=3/6 is ca. 17 %, whereas when p=1 it is ca. 33 %.For the proposed scheme with p=1 (resp., p=3/6), the maximum value of the dispersion error increases (resp., decreases) as θ increases. When the propagation angles are considered, the proposed scheme with p=3/6 has smaller dispersion errors than one with p=1.

Figure 7: Relative wave speed errors of the Bathe method (CFL=1) for various values of the propagation angle with consistent mass matrix (dotted lines: discarded wave modes)

Figure 8: Relative wave speed errors of the proposed scheme (p=1 and CFL=1) and CDM (CFL=1) for various values of the propagation angle with lumped mass matrix;the results for the proposed scheme (p=1) and CDM are identical

Figure 9: Relative wave speed errors of the proposed scheme (p=3/6 and CFL=0.65)and trapezoidal rule (CFL=0.65) for various values of the propagation angle with consistent mass matrix; the results for the proposed scheme (p=3/6) and trapezoidal rule are identical

Figure 10: Convergence rate of relative wave speed errors for the proposed scheme(p=5/6) with consistent mass matrix

Figure 11: Relative wave speed errors of the proposed scheme (p=4/6 and CFL=0.705) for various values of the propagation angle with consistent mass matrix

Furthermore, as the proposed scheme with p=4/6 and p=5/6 yields accurate solutions (see Section 2.3), we additionally consider two parameters (p=4/6 and p=5/6) for the dispersion analysis. For the proposed method with p=5/6, the relative wave speed errors are positive for any CFL number (see Fig. 10). The maximum value of the relative wave speed errors decreases with a decrease in the CFL number. The maximum value of the relative wave speed errors is converged after the CFL number is ca. 0.2. For the proposed method (p=5/6) with CFL=0.2, the maximum speed error is ca. 20%. However, for the proposed method with p=4/6, the maximum speed error is ca. 5% when the optimal CFL value (=0.705) is used (see Fig. 11). Although the proposed scheme with p=5/6 yields a more accurate solution than one with p=4/6,the maximum speed error of the proposed scheme (p=5/6) is approximately four times that of the proposed method (p=4/6). In addition, comparing the minimum dispersion errors in a two-dimensional analysis when the parameter p takes on the values 3/6, 4/6,5/6, and 1, we note that the proposed method with p=4/6 provides a minimum dispersion error in the two-dimensional analysis. Therefore, the parameters p=1 and p=4/6 are recommended in one- and multi-dimensional cases, respectively.

3 Numerical simulations

3.1 Wave propagation problem

To test the proposed scheme for wave propagation solutions, the pre-stressed membrane problem is considered, as can be seen in Fig. 12. We compare the solutions of the twodimensional wave propagation using the proposed scheme, CDM, the Noh-Bathe method,the trapezoidal rule, and the Bathe method. The governing wave equation of the prestressed membrane problem is expressed as follows [Noh and Bathe (2013); Noh, Ham and Bathe (2013); Wen, Duan, Yan et al. (2017)]:

where u is the transverse displacement and. In this problem, we assume that the exact wave speed c0is 1.0. Due to symmetry, the domain [0,10]×[0,10] is only considered. The load is given as follows:

The exact solution of this problem is obtained using Green’s function G(x,y,t) as follows [Wen, Duan, Yan et al. (2017)]:

Figs. 13-14 show the numerical and analytical solutions of the displacement and velocity considering a 60×60 finite element mesh and various propagating angles at time t=7.36. The numerical results for the proposed scheme with p=1, CDM, and the Noh-Bathe method are shown in Fig. 13, whereas Fig. 14 shows the wave problem solutions of the proposed scheme with p=4/6, the trapezoidal rule, and the Bathe method. As shown in Fig. 13, the solutions using the proposed scheme with p=1 and CDM are less accurate than those using the Noh-Bathe method as these schemes have high dispersion errors. However, in the Noh-Bathe method, it is difficult to describe the high frequency wave shape as there are discarded frequency modes (see Fig. 6).

Furthermore, the proposed scheme (p=4/6), the trapezoidal rule, and Bathe method yield more accurate solutions than the proposed scheme (p=1), CDM, and Noh-Bathe method. The numerical solutions of the Bathe method do not describe the high frequency wave mode as the Bathe method is a dissipative scheme, whereas the proposed scheme with p=4/6 and trapezoidal rule yield solutions including all the wave modes. In the numerical results of the non-dissipative implicit scheme, as the proposed scheme(p=4/6) has a smaller PE than the trapezoidal rule (see Section 2.3), we observe that the proposed scheme has more accuracy than the trapezoidal rule.

Figs. 15-16 show the numerical results for the displacement at t=7.36 for a 140×140 finite element mesh using the lumped and consistent mass matrices a respectively. In the proposed scheme (p=4/6), the results for θ=are more accurate than those for θ=0 in the case of the displacement and velocity. The difference can be explained by the wave speed errors. When θ=, there are few spurious oscillations as the dispersion error is almost zero, whereas the spurious oscillations occur forward of the wave front when θ=0, because the speed errors are positive (see Fig. 11). We observe that the Noh-Bathe method and Bathe method yield the solutions, including the high frequency mode, when a high spatial discretization density is used (e.g., 140×140 mesh) as the minimum wave number decrease as the element size Δx increases (see Eq. (28)).However, in the case of a low spatial discretization density (e.g., 60×60 mesh), the dissipative schemes provide less accurate solutions of the high wave mode as the minimum wave number increases.

Figure 12: Pre-stress membrane problem: initial displacement, velocity, and acceleration are 0; exact wave speed is 1.0

Figure 13: Displacement and velocity variation of the proposed scheme (p=1), CDM,and Noh-Bathe method along various propagating angles, at time t=7.36, 60×60 finite element mesh

Figure 14: Displacement and velocity variation of the proposed scheme (p=4/6),trapezoidal rule, and Bathe method along various propagating angles, at time t=7.36,60×60 finite element mesh

Figure 15: Displacement and velocity variation of the proposed scheme (p=1), CDM,and Noh-Bathe method along various propagating angles, at time t=7.36, 140×140 finite element mesh

Figure 16: Displacement and velocity variation of the proposed scheme (p=4/6),trapezoidal rule, and Bathe method along various propagating angles, at time t=7.36,140×140 finite element mesh

To evaluate the accuracy of time integration scheme, the Euclidean norm error Ej,θ(j=0,1 and θ =0,) is used. It is defined as follows:

In the implicit time integration scheme, the proposed scheme with p=4/6 provides highly accurate solutions regardless of spatial discretization density for displacement and velocity (see Fig. 18). Moreover, in the proposed scheme (p=4/6), the results for θ=are more accurate than those for θ=0 in the case of the displacement and velocity as the dispersion error for θ=is smaller than one for θ=0. In contrast, we note that for the trapezoidal rule, the numerical results for θ=0 are more accurate owing relative speed errors. Although the proposed scheme (p=4/6) and trapezoidal rule have similar average values of the error norms when the high spatial discretization density (e.g., 140×140 mesh) is used, the proposed scheme (p=4/6) yields more accurate results than other methods even when using the low spatial discretization density.Therefore, the parameter p=4/6 is recommended for wave propagation, and the proposed scheme will suffice for usual engineering use.

In Fig. 19, the computational efforts of the various methods are compared. For the explicit method, the Noh-Bathe method costs more time than the proposed method(p=1) and CDM as the Noh-Bathe method has two sub-steps and a CFL number that is 1.85 times that of the proposed scheme (p=1) and CDM. On the other hand, for the implicit method, the Bathe method with two sub-steps requires more time than the other methods. Moreover, comparing the computational efforts of the proposed scheme(p=4/6) and trapezoidal rule, the proposed scheme (p=4/6) requires less effort as it has a CFL number that is 1.08 times that of the trapezoidal rule. Therefore, the proposed methods (p=1 and p=4/6) are desirable for wave propagation problems on the explicit and implicit methods, respectively, owing to their reasonable accuracy and acceptable time cost.

Table 2: Euclidean norm errors of various schemes, at time t=7.36, 60×60 finite element mesh

Table 3: Euclidean norm errors of various schemes, at time t=7.36, 140×140 finite element mesh

Figure 17: Accuracy for various explicit schemes in wave propagation problem

Figure 18: Accuracy for various implicit schemes in wave propagation problem

Figure 19: Required computational efforts of various methods considering the wave propagation problem at time t=7.36

3.2 Standard undamped SDOF system

The free vibration of a standard undamped SDOF (single-degree-of-freedom) system is used to test the accuracy of the proposed scheme. To solve the undamped SDOF system,we consider the representative cases of the proposed scheme: p=1 (with desirable dispersion in the one-dimensional case), p=5/6 (with highest accuracy), and p=4/6(with desirable dispersion in the multi-dimensional case). The equilibrium equation of the undamped SDOF system is given by Eq. (17) [Wen, Duan, Yan et al. (2017)].

In this section, we assume that the natural frequency is ω = π (that is, the natural period is Tn=2), and the considered time duration T is 4 sec. In order to calculate the accuracy and efficiency of various schemes, the Euclidean norm error Ej(j=0,1,2) is also used as follows:

where u(j)(i.e., u,, and) andare the numerical results and the exact results at discrete time, respectively.

Figure 20: Accuracy and convergence rate for various schemes in undamped SDOF system

Fig. 20 shows the Euclidean norm error of the proposed scheme, CDM, the Noh-Bathe method, the trapezoidal rule, and the Bathe method. In general, the relation between displacement, velocity and acceleration is linear. However, for the Noh-Bathe method,the relation between the displacement and velocity, or between velocity and acceleration is not linear [Wen, Duan, Yan et al. (2017)]. Besides, for the proposed method (p=5/6),a linear relation exists between displacement and acceleration because of the applied computation procedure. For displacement, velocity, and acceleration, the proposed schemes (p=1 and p=4/6) and the Bathe method have similar error norms. The results of the proposed schemes (p=1 and p=4/6) and CDM are similar only for displacement. The proposed scheme (p=5/6) yields similar numerical results for velocity to those of the Noh-Bathe method, which provides highly accurate solutions.Moreover, the proposed scheme (p=5/6) yields more accurate results for displacement and acceleration than the other schemes. This corresponds to the results for PE (see Fig.3). As the proposed scheme has no sub-steps, the proposed scheme has a lower computational cost. Therefore, the proposed time integration method with p=5/6 is more effective for the undamped SDOF system.

Figure 21: Required computational efforts of various methods in undamped SDOF system when Δt/T=0.02 is used

3.3 Damped SDOF system

The forced vibration of a damped SDOF system is used to verify the accuracy of a new numerical scheme [Wen, Duan, Yan et al. (2017)]. In the viscously damped SDOF system, which is shown in Fig. 22, the equation of motion is given as follows [Rao(2011)]:

where x(t) and y(t) are the displacements of the mass and base from their static equilibrium positions at time t, respectively. In this study, we assume that the mass m is 10 kg/m, viscous damping coefficient c is 20 N-s/m, and spring constant k is 4000 N/m. The displacement of the base is y(t)=sin(5t) m and the initial conditions are x(0)=0.1 m and x(0)=0 m/s.

Figure 22: Damped system under the hormonic motion of the base

Figure 23: Homogenous, particular, and general solutions in the example 3.3

As seen in Fig. 23, the exact solution of this problem is x(t)=x?(t)+xp(t), where the homogenous solution x?(t) is 0.1004e-tcos(19.975t+0.0831) and the particular solution xp(t) is {53.331sin(5t)-0.0888cos(5t)}×103. Fig. 24-26 show the relative error norms, which are defined by the following equation, when the time increment Δt is 0.2Tn(the natural period

It is important to note that damped SDOF problem is used to analyze the effect of time integration schemes for harmonic and transient responses. The error curves in Figs. 24-26 start converging after the time is τ (≈ 10 s), where τ is the time when the amplitude of the transient response is c.a. 0.0015%; that is, the transient response is negligible (see Fig.23). Therefore, the physical meaning of the relative error is the accuracy of the transient responses when t≤ τ, whereas the physical meaning of the error is the accuracy of the harmonic response when t＞ τ.

As seen in Figs. 24-26, the proposed scheme with p=4/6 provides more accurate solutions for velocity than the other methods when t＞ τ. However, for the proposed scheme (p=4/6), the numerical results for velocity are less accurate than those of other methods when t≤ τ, and this scheme has low accuracy in the prediction of the displacement and acceleration. In contract, for displacement and acceleration, the proposed scheme with p=5/6 has the highest accuracy compared with the other schemes, and the numerical results for velocity are more accurate than those of other methods when t≤ τ.

Figs. 27-29 show the accuracy and convergence rate using Eq. (33) for t≤ τ, τ＜ t≤100Tn, and t≤100Tncorresponding to the transient response, harmonic response, and both transient and harmonic responses, respectively. The Noh-Bathe method has a similar accuracy as that of the proposed scheme with p=5/6 when t≤ τ, and shows a lower accuracy than the proposed scheme (p=5/6) when t＞ τ. Therefore, for the transient response, the proposed scheme (p=5/6) and the Noh-Bathe method can yield accurate solutions. For the harmonic response, the proposed scheme with p=5/6 provides more accurate solutions for displacement and acceleration than the Noh-Bathe method,whereas the p=4/6 case shows high accuracy for velocity when compared with other schemes. However, as shown in Fig. 29, for both transient and harmonic responses, the proposed scheme (p=5/6) can be used to acquire desirable solutions in engineering problems. Therefore, we can identify that the optimal parameter p is 5/6 for the damped SDOF system.

Furthermore, as the Noh-Bathe method and Bathe method have two sub-steps, the computational time is approximately twice as high as that of the proposed scheme, CDM,and trapezoidal rule. Fig. 30 shows the quantitative comparison of the computational time/cost for various schemes. As shown in Fig. 30, the proposed scheme (p=5/6) that yields small errors can be used to reduce computational efforts.

Figure 24: Accuracy of the displacement for various schemes in damped SDOF system when the time increment Δt is 0.2Tn

Figure 24: (continued)

Figure 24: (continued)

Figure 25: Accuracy of the velocity for various schemes in damped SDOF system when the time increment Δt is 0.2Tn

Figure 25: (continued)

Figure 25: (continued)

Figure 26: (continued)

Figure 26: (continued)

Figure 27: Accuracy and convergence rate of transient responses for various schemes in damped SDOF system

Figure 28: Accuracy and convergence rate of harmonic response for various schemes in damped SDOF system

Figure 29: Accuracy and convergence rate of both harmonic and transient responses for various schemes in damped SDOF system

Figure 30: Required computational cost of various methods for example 3.3

Figure 31: Two degrees-of-freedom spring system

3.4 Two degrees-of-freedom spring system

In the two degrees-of-freedom system shown in Fig. 31, we compare the solutions of the velocity and acceleration by using the proposed scheme, CDM, Noh-Bathe explicit method, trapezoidal rule, and Bathe implicit method. The governing equation of two degrees-of-freedom system is expressed by as follows:

In the present study, we assume that m1=1,m2=1,k1=10 and k2=10. Also, we consider the external force F(t) at the node 1 as follows:

with A=10 and ωp=3. In this example, the considered time increment Δt is 0.1. In this study, the reference solution and numerical solutions are investigated based on the time history of velocities and accelerations of each method at node 1 and 2. Newmark method with computational time increment of Δt=0.0001 is referred to as ‘reference solution’. As shown in Figs. 32-33, CDM, the trapezoidal rule, and Bathe method provide large period errors as time increases, while the proposed scheme and Noh-Bathe method present smaller period errors as shown in Figs. 32-33. Although the Noh-Bathe method has reasonably accurate solutions, the accuracy of the proposed method is higher than the Noh-Bathe method.

Figure 32: Numerical solutions of node 1 for two degrees-of-freedom problem

Figure 33: Numerical solutions of node 2 for two degrees-of-freedom problem

3.5 Beam with boundary springs

In this section, we consider a beam with boundary constraints (see Fig. 34), modeled with five-node element. The specific boundary conditions are as shown in the following equation:

The important point to note here is that the boundary conditions (i.e., w(0,t)=0 and w(L,t)=0) is equivalent to kw1=∞ and kθ2=∞, respectively. However, in implementing the computations, it is impractical to apply kw1=kθ2= ∞ due to floating point precision.

We consider the applied force f(L/2,t) at the center node of Fig. 34 as follows:

where the forcing frequency is set to ff=1.5(ω1/2/π). Here, ω1is the fundamental frequency of the beam vibrations.

Figs. 35-36 shows the time responses of the beam at its mid-span using the proposed scheme, trapezoidal rule, Bathe method. The step increment Δt used for the three methods is 2.28e-7sec. Newmark method with step size of Δt=2.28e-8sec is referred to as ‘reference solution’. As shown in Figs. 35-36, the solutions provided by the proposed method are more accurate than those yielded by other methods.

Figure 34: Beam with boundary springs modeled using five-node elements

Figure 35: Vertical time response at the midspan of the beam

Figure 36: Rotational time response at the midspan of the beam

4 Concluding remarks

In the present study, the direct time integration scheme based on Gauss-Legendre quadrature is proposed for structural dynamics. The proposed scheme possesses secondorder accuracy and is non-dissipative. Moreover, this scheme has one parameter that affects stability, accuracy, and dispersion. When p≠1, the proposed scheme is an implicit method, as the methods used to solve for displacement and velocity are implicit and explicit, respectively. However, when p=1, the proposed scheme is explicit method. By using the appropriate value of the parameter, the desirable stability, accuracy,and dispersion could be achieved. When p＜1, the stability limit of the proposed scheme is larger than that of CDM, whereas the proposed scheme with p≤0.5 has unconditional stability.

In order to show the uniqueness of the proposed algorithm, the solutions to the wave propagation problem are obtained using various schemes, and the numerical solutions are compared with the exact solutions. The proposed scheme (p=4/6) provides more accurate numerical results than other methods (including explicit and implicit methods)even in the case of a low spatial discretization density. This requires a large computational cost as compared to explicit methods; however, this method requires less time than other implicit methods. Furthermore, we presented the performance of the solutions of a standard undamped system and a damped system. In a standard undamped SDOF system, the proposed scheme (p=5/6) has the highest accuracy. In addition, in a damped SDOF system, two degrees-of-freedom system, and time responses of the beam,the proposed scheme (p=5/6) has higher accuracy. As there is no sub-step in the calculation of the solutions, the computational cost and time of the proposed scheme is lower. In addition, the accuracy of the proposed scheme is higher than that of the Noh-Bathe method and Bathe method consisting of two sub-steps. Therefore, the proposed scheme is desirable for linear structural dynamics owing to its accuracy and acceptable time cost.

The proposed scheme has the following technical merits and demerits in engineering. The proposed scheme has the same or a higher stability limit than CDM when p≤1; that is,the proposed scheme is good for practical use as it has good stability characteristics. As shown in Section 2.3, the method with p=5/6 has high accuracy characteristics when compared with those of the other methods, which is in good agreement with the numerical results provided in Sections 3.2-3.3. The proposed scheme is less expensive than both the Noh-Bathe method and the Bathe method, since the proposed scheme has no sub-step, making it almost 2 times cheaper than the others. On the other hand, the proposed scheme (p＞0.5) lacks unconditional stability, which has a critical time step corresponding to the stability limit. However, for the proposed scheme, the time step corresponding to the optimal CFL number is smaller than the critical time step. Therefore,in wave propagation solutions, the proposed method is stable and yields accurate solutions in engineering

For the general structure dynamics, two most important issues are the linear and nonlinear dynamic response of a finite element system. In this study, we consider only the linear equations in dynamic analysis. However, the proposed scheme discussed previously for linear problems can also be extended for nonlinear dynamic analysis.Therefore, further research on the solution of the nonlinear dynamic response by using the proposed scheme would be necessary in the future.

Acknowledgement:This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Ministry of Science and ICT (MSIT) (No.2018R1A2B6007403). This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) through GCRC-SOP (No.2011-0030013).

Appendix A. Local truncation error analyses

Local truncation error measures the order of accuracy of time integration methods at an arbitrary time step, assuming that the method was exact at the previous time step [Chung and Lee (1994); Zhang, Liu and Liu (2017); Kim and Choi (2018)]. In a single iteration of a time integration algorithm, the local truncation error is defined as the difference between the exact solution and the numerical solution. Local truncation error analyses similar to those described in Zhang et al. [Zhang, Liu and Liu (2017)] have been performed.

The equation of motion in Eq. (2) at time t and t+Δt are as follows:

Where u(t), u(t), ü(t), and F(t) are the numerical solutions of the exact displacement ut,velocity ut, acceleration üt, and external nodal force vector Ftat time t, respectively.The exact solutions satisfying Eq. (2) at time t and t+Δt yields

Substituting Eqs. (6)-(7) and half-step CDM into Eq. (41), we obtain the following relation:

As there is no previous truncation error at time t and the external nodal force vectors are known, the following relationships hold:

The acceleration local truncation error eü(t+Δt)and the corresponding weighted local truncation error eü(t+Δt)are defined as follows:

By using Eqs. (43) and (46) to express the external nodal force vector in Eq. (44) and Taylor series expansion to approximate ut+Δt,ut+Δt, and üt+Δt, substituting Eq. (44)into Eq. (48) yields

By expanding Eq. (47) using Taylor series and comparing it with Eq. (49), the acceleration local truncation error is as follows:

Similarly, for velocity and displacement, we obtain the following local truncation errors:

As shown in Eqs. (50)-(52), the local truncation errors are O(Δt3) for q=(1-p)/2.That is, the proposed scheme is second-order accurate in time for any p if q=(1-p)/2. Besides, as the damping term is included in Eqs. (50)-(52), the proposed method is always second-order accurate with and without damping.

Appendix B. Integration and load operators

The integration operator A and load operators Laand Lbare computed using the following equations:

where a=1+ ξΩ0, β =2-a, and