A Time Integration Method Based on Galerkin Weak Form for Nonlinear Structural Dynamics

This paper presents a step-by-step time integration method for transient solutions of nonlinear structural dynamic problems. Taking the second-order nonlinear dynamic equations as the model problem, this self-starting one-step algorithm is constructed using the Galerkin finite element method (FEM) and Newton–Raphson iteration, in which it is recommended to adopt time elements of degree m = 1,2,3. Based on the mathematical and numerical analysis, it is found that the method can gain a convergence order of 2m for both displacement and velocity results when an ordinary Gauss integral is implemented. Meanwhile, with reduced Gauss integration, the method achieves unconditional stability. Furthermore, a feasible integration scheme with controllable numerical damping has been established by modifying the test function and introducing a special integral rule. Representative numerical examples show that the proposed method performs well in stability with controllable numerical dissipation, and its computational efficiency is superior as well.


Introduction
The step-by-step time integration method is the most commonly used numerical method for transient analysis of nonlinear dynamics, and can in general be classified into an explicit and implicit scheme. For the explicit integration algorithm, it is not required to solve coupled equations in each time step, and this method is more suitable for determining solutions of wave propagation problems. While for the implicit algorithm, although the requirement for solutions of coupled equations in each time step makes it expensive, it still has been widely used in solid mechanics and structural engineering because of the superior numerical stability. In fact, it may often occur that an implicit algorithm of unconditional stability in linear problems fails to give a stable response for some nonlinear problems. For example, significant instability may appear in the trapezoidal rule [1] when it is used for the evaluation of large deformation problems and long time range responses [2][3][4][5].
For linear dynamics, the stability of an integration algorithm can be estimated by spectral analysis, whereas in nonlinear dynamics spectral analysis is only one of the requirements to remain stable, and the conservation or decay of total energy within a time step interval should be the very sufficient condition for effective solutions [3]. From this point of view, the time integration method can be classified into three categories: One of numerical dissipation, one of enforced conservation of energy and one of algorithmic conservation of energy. Firstly, the algorithm belonging to the numerical dissipation category uses numerical dissipation to achieve energy conservation or decay. For example, numerical damping of high-order-frequency can be introduced into the Newmark method with γ > 1/2, however, in this case, the method only possesses an accuracy of the first order. For another example, in the HHT-α [6] method and generalized-α [7] method, the numerical accuracy has been improved to the second order and controllable numerical damping has been introduced as well. For such methods, the suitable selection of parameters turns out to be the key to achieving effective energy decay. Secondly, For nonlinear dynamic problems in solid mechanics and structural engineering, discretization in the space dimension generally results in the following governing equations in the time domain, which are named dynamic equations hereinafter:  (1) is a typical IVP of nonlinear ordinary differential equations (ODEs). A common strategy for solving such ODEs is to transform them into linear equations and then perform iteration. Thus the Newton-Raphson method is first applied to Equation (1) which results in the following linearized ODEs: (3) u (k) ) + (f s (k) Here u (k) , , tangential stiffness matrix -K (k) and internal force vector f (k) int are functions of time. For convenience, we omit the superscript 'k' in Equation (2) temporarily and denote f = f ext − f int . Then, Equation (2) can be written into the following simplified form, which is the main solution target of the next subsection:

Galerkin Finite Element Method
To solve the linear IVP of ODEs in Equation (5), a time-domain FEM based on the Galerkin weak form [25,34,42] was used in the present paper. Let v ∈ H 1 v and u h ∈ H 1 E denote the test function and the trial function, respectively. The weighted residual weak form of Equation (5) might be established in time interval [t n , t n+1 ] as follows: where H 1 E and H 1 v represent the trial function space and the test function space, respectively. Taking [t n , t n+1 ] as a typical time element of the FEM with ∆t = t n+1 − t n , as shown in Figure 1, the m-degree polynomials are adopted as the trial function and test function defined on this typical element, i.e.,: where the shape functions N i (t) (i = 1, 2, . . . , m + 1) are Lagrange polynomials, the shape function matrix N = [N 1 I, N 2 I, · · · , N m+1 I], and I is the unit matrix, T is the vector of nodal test displacements.  Figure 1, the m-degree polynomials are adopted as the trial function and test function defined on this typical element, i.e.,: where the shape functions () is the vector of nodal test displacements. Substituting Equation (7b) into Equation (6), because of the arbitrariness of test displacement vector v d , the following equation can be derived: After implementing integration by parts for the first term of Equation (8), which is: where +1 n u and n u represent the velocity at +1 n t and n t , respectively. Substituting Equation (7a) into Equation (8), and then taking a time coordinate mapping from 1 [ , ] nn t t t +  to [ 1,1]  − , the following algebraic equations are yielded: Equation (10) is equivalent to the equation: Substituting Equation (7b) into Equation (6), because of the arbitrariness of test displacement vector d v , the following equation can be derived: After implementing integration by parts for the first term of Equation (8), which is: where . u n+1 and . u n represent the velocity at t n+1 and t n , respectively. Substituting Equation (7a) into Equation (8), and then taking a time coordinate mapping from t ∈ [t n , t n+1 ] to ζ ∈ [−1, 1], the following algebraic equations are yielded: with Equation (10) is equivalent to the equation: which could further be written into the following tensor notation form: with δ i m+1 and δ i1 being the Kronecker delta function and: As a result, Equation (15) is the step-by-step solution formula derived from the Galerkin FEM to solve the time-dependent problem in Equation (5).

Numerical Integration
In Equation (14), there exists a number of integral terms. Since the integral accuracy of these terms might affect the numerical stability of the whole algorithm, the scheme of Gauss numerical integration is briefly described and discussed in this subsection.
With ω l and ζ l denoting the Gaussian integral coefficients and locations of n l Gaussian points in [−1, 1], Equations (16) and (17) might be expressed as: where .

Iterative Algorithm
Equation (15) is a general step-by-step solution formula for the linear IVP of Equation (5) with time finite elements of degree m, where m can be any positive integer, in theory. However, in practice, examples have shown that the formulae with m = 1,2,3, i.e., linear, quadratic, and cubic time elements, are more pragmatic and sufficiently effective, with no need that m to be much higher. Therefore, the detailed computation schemes of Equation (15) for m = 1,2,3 are given in this subsection, and then the entire iterative algorithm of the Galerkin time FEM is established for the solution of the nonlinear IVP in Equation (1).
For the linear Lagrange time element, there are only two nodes in the single element in Figure 1, i.e., t 1 = t n and t 2 = t n+1 . In this case, Equation (14), the step-by-step solution formula of the linear IVP, can be expressed as: where u n and . u n are the initial displacement and velocity at time t n , respectively. Furthermore, the final solution formula with the omitted iteration symbol "k" being added can be expressed as: Similarly, for the quadratic and cubic Lagrange time elements, the solution formula from Equation (14) can be respectively transformed into: and It is obvious that the coupled equations above need to be solved in each time step, i.e., the proposed algorithm is an implicit step-by-step time integration method. The iterative algorithm of this time integration method can be summarized into the following steps, as shown in Algorithm 1. Solve the dynamic equation at each time point t n + ∆t.

4.1
For the k-th iterative step, input the initial data: d Form the component matrix: Form the component vector: Solve the following equation: with Equations (20) and (21), let k = k + 1, and then return to step 4.2; otherwise go to the next step.
If t n + ∆t < t f , let t n = t n + ∆t and return to step 4; otherwise go to the next step.

Analysis of Numerical Stability and Accuracy
Numerical stability and accuracy are the two essential indicators to evaluate the properties of a numerical integration method. In this section, these two features of the present algorithm are analyzed and discussed, especially focusing on the schemes with a reduced Gaussian integral.

Stability
To analyze the stability of the proposed algorithm, the following single-degree-of-freedom system is considered: where ω is the circular frequency of free vibration. When the proposed algorithm is applied to solve this linear dynamic equation, Equations (16) and (17) can be expressed as: Taking the linear element as an example, all the integral terms in Equation (27) could be integrated exactly by more than two Gaussian points. If we only use one Gauss integral point, i.e., the reduced Gauss integration, the recursive difference equation at any two successive time points would be derived as follows: In this case, the present algorithm automatically degrades into the so-called 'GW (Galerkin Weak form) method' which has been already proposed [34] for linear elastodynamic problems, and turns out to be an unconditionally stable time integration method with the spectral radii ρ(A) = 1, similar to quadratic and cubic elements. Borri [25] discovered that, when reduced Gauss integration was used in calculation of the stiffness matrix, the GW method for linear elastodynamic equations is unconditionally stable, otherwise, it is conditionally stable. This conclusion might be understood as the loss of accuracy earns the improvement of stability.

Accuracy
Also taking the single-degree-of-freedom system in Equation (26) in consideration, the corresponding initial conditions of the problem are set to be: The exact solution of IVP Equations (26) and (29) is u = cos ωt. Comparing the numerical solution obtained from the proposed method with this exact solution, the errors can be measured in terms of period elongation (PE) and amplitude decay (AD). The percentage period elongation with respect to the ratio ∆t/T is plotted in Figure 2a for the proposed method, the trapezoidal rule, and Bathe method, where T denotes the period of free vibration in Equation (28). To more clearly show the results with quadratic and cubic elements, the two curves in Figure 2a are redrawn in Figure 2b,c, correspondingly. It is worth mentioning that since the proposed algorithm has no numerical dissipation (ρ(A) = 1), there is no amplitude decay when ∆t/T varies. Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 23

A Feasible Scheme with Numerical Dissipation
To be sufficiently effective in solutions of structural dynamics, the time integration scheme should possess controllable numerical dissipation, such that the false vibration components of higher orders could be removed and the accuracy of results could be guaranteed. However, as discussed in Section 3, the integral scheme of Equation (15) does not have such numerical dissipation. In order to get over this difficulty, a feasible time integration scheme with numerical dissipation is proposed using the Galerkin weak form method. The detailed construction procedure is as follows.
Considering the Galerkin weak form shown in Equation (6) with linear time elements, let the trial function remain as the common interpolation of linear polynomials while the test function is as follows: where  i N (i = 1,2) is the following constructed linear polynomial with a parameter '  ' being introduced, i.e.,:

A Feasible Scheme with Numerical Dissipation
To be sufficiently effective in solutions of structural dynamics, the time integration scheme should possess controllable numerical dissipation, such that the false vibration components of higher orders could be removed and the accuracy of results could be guaranteed. However, as discussed in Section 3, the integral scheme of Equation (15) does not have such numerical dissipation. In order to get over this difficulty, a feasible time integration scheme with numerical dissipation is proposed using the Galerkin weak form method. The detailed construction procedure is as follows.
Considering the Galerkin weak form shown in Equation (6) with linear time elements, let the trial function remain as the common interpolation of linear polynomials while the test function is as follows: where N αi (i = 1,2) is the following constructed linear polynomial with a parameter 'α' being introduced, i.e.,: That is to say, the test function in Equation (32) doesn't satisfy the C 0 continuity between elements any more. Substituting such a trial function u h and test functionṽ h into Equation (6), in consideration of the arbitrariness of d v , the following solution formula for Equation (5) may be obtained: In Equation (32), there are many integration terms in both K α and P α which are also with respect to the parameter α. Again, for the target to introduce numerical damping, the following integral formulae are used instead of the common Gauss integral in calculations of all the integral terms, i.e.,: Finally the following Newton-Raphson iterative scheme is yielded: This is the proposed time integration scheme with numerical dissipation for linear time elements based on the Galerkin weak form.
Here, we give some analysis and verification of the introduced numerical dissipation. For the integration scheme shown in Equations (34)- (38), the transfer matrix A is: When α = 0, Equation (39) is the same as Equation (28), thus the scheme with a linear element in Equation (25) is a particular case of this improved scheme. The curves of the spectral radius obtained from A are shown in Figure 3 when α is set to be different values. From Figure 3, it can be seen that the spectral radius of the algorithm is unconditionally stable, and the numerical damping is easily controllable, since following an increase of α the numerical damping increases correspondingly. The curves of the spectral radius of the trapezoidal rule, Wilson θ-method (θ = 1.4) [41], and Bathe method (γ= 1/2) [5] are shown in Figure 3 as well. The percentage of period elongation and amplitude decay are shown in Figure 4.

Numerical Examples
In this section, three representative numerical examples are given to show the efficiency of the proposed method. For each example, the results of the present method are compared with those of the following four algorithms: Tshe trapezoidal rule, Newmark [41], and Bathe method (  = 1/2 ) [5], which are commonly-used time integration algorithms in dynamics. All the results are obtained from the programming codes of MATLAB 2017a under the same computational circumstance.

The Simple Pendulum
The governing equation [43] for the nonlinear pendulum shown in Figure 5 is,

Numerical Examples
In this section, three representative numerical examples are given to show the efficiency of the proposed method. For each example, the results of the present method are compared with those of the following four algorithms: Tshe trapezoidal rule, Newmark method (  = 1/2, = 1/6 ) [1], Wilson  -method (  = 1.4 ) [41], and Bathe method (  = 1/2 ) [5], which are commonly-used time integration algorithms in dynamics. All the results are obtained from the programming codes of MATLAB 2017a under the same computational circumstance.

The Simple Pendulum
The governing equation [43] for the nonlinear pendulum shown in Figure 5

Numerical Examples
In this section, three representative numerical examples are given to show the efficiency of the proposed method. For each example, the results of the present method are compared with those of the following four algorithms: Tshe trapezoidal rule, Newmark method (γ = 1/2, β = 1/6) [1], Wilson θ-method (θ = 1.4) [41], and Bathe method (γ = 1/2) [5], which are commonly-used time integration algorithms in dynamics. All the results are obtained from the programming codes of MATLAB 2017a under the same computational circumstance.

The Simple Pendulum
The governing equation [43] for the nonlinear pendulum shown in Figure 5 is, .. u + g/L sin(u) = 0 with initial conditions u(0) = π/2 rad and . u(0) = 0 rad/s, where g is the gravitational acceleration, L is the length of a massless suspension, and u is the angle between the cycloid and vertical plane. The kinetic and potential energy of the system are T E = 0.5L 2 . u 2 and V E = g(L − L cos u), respectively.
Moreover, the total energy is the sum of them. The solutions obtained from the trapezoidal rule with ∆t = 0.0001 s are approximately regarded as the exact solutions to evaluate the accuracy and efficiency of the present method here. This example aims to test the stability and accuracy of the proposed algorithm and the scheme without using numerical damping.   Figure 8, which are all stable. It can be seen that, for Newmark method, the error of response gets bigger as the time goes on. For the Wilson method and Bathe method, noteworthy amplitude decay takes place in displacement responses since there is ineradicable numerical damping in these two methods, especially in the Wilson method. For this pendulum example, using methods without numerical damping may obtain more suitable responses, but there is no way for the Bathe method and Wilson method to control or remove the numerical damping, which will significantly affect the accuracy of the solution. The comparisons of energy responses shown in Figure 9 can further demonstrate the numerical dissipation property of the Wilson and Bathe methods. Both the displacement and velocity results of the present method with m = 2,3 are shown in Figure 10, which are stable and of high accuracy. Let g/L = 4π 2 s −2 and take the time step to be ∆t = 0.15 s. The pendulum motion is calculated with the trapezoidal rule, Newmark method, Bathe method, Wilson θ-method, and the proposed method. The displacement, velocity, and total energy histories from the trapezoidal rule are presented in Figures 6 and 7, from which it can be seen that the results become considerably large when the time approaches 30 s, i.e., an unstable response is given from the trapezoidal rule for this example. The displacement responses of the exact solution and the other four methods are shown in Figure 8, which are all stable. It can be seen that, for Newmark method, the error of response gets bigger as the time goes on. For the Wilson method and Bathe method, noteworthy amplitude decay takes place in displacement responses since there is ineradicable numerical damping in these two methods, especially in the Wilson method. For this pendulum example, using methods without numerical damping may obtain more suitable responses, but there is no way for the Bathe method and Wilson method to control or remove the numerical damping, which will significantly affect the accuracy of the solution. The comparisons of energy responses shown in Figure 9 can further demonstrate the numerical dissipation property of the Wilson and Bathe methods. Both the displacement and velocity results of the present method with m = 2,3 are shown in Figure 10, which are stable and of high accuracy.  Figures 6 and 7, from which it can be seen that the results become considerably large when the time approaches 30 s, i.e., an unstable response is given from the trapezoidal rule for this example. The displacement responses of the exact solution and the other four methods are shown in Figure 8, which are all stable. It can be seen that, for Newmark method, the error of response gets bigger as the time goes on. For the Wilson method and Bathe method, noteworthy amplitude decay takes place in displacement responses since there is ineradicable numerical damping in these two methods, especially in the Wilson method. For this pendulum example, using methods without numerical damping may obtain more suitable responses, but there is no way for the Bathe method and Wilson method to control or remove the numerical damping, which will significantly affect the accuracy of the solution. The comparisons of energy responses shown in Figure 9 can further demonstrate the numerical dissipation property of the Wilson and Bathe methods. Both the displacement and velocity results of the present method with m = 2,3 are shown in Figure 10, which are stable and of high accuracy.           Tables 1 and 2, and the convergence curves are shown in Figure 11 correspondingly. The present algorithm with the linear element, quadratic element, and cubic elemen,t respectively, possess convergence orders of 2, 4 and 6 (i.e., 2 m) for both displacement and velocity, while all the other four methods are of a convergence order of 2. In this case, the accuracy of our method with m = 2,3 is obviously higher than that of others. For m = 1, the accuracy of the present algorithm is similar to that of the Newmark and Bathe methods, and is better than that of the Wilson  -method.
where u h i and u i are the approximate displacement and the exact displacement at step i and N is the total number of time steps. Error ε and convergence order ρ from these seven schemes are given in Tables 1  and 2, and the convergence curves are shown in Figure 11 correspondingly. The present algorithm with the linear element, quadratic element, and cubic element respectively, possess convergence orders of 2, 4 and 6 (i.e., 2 m) for both displacement and velocity, while all the other four methods are of a convergence order of 2. In this case, the accuracy of our method with m = 2,3 is obviously higher than that of others. For m = 1, the accuracy of the present algorithm is similar to that of the Newmark and Bathe methods, and is better than that of the Wilson θ-method.

Long Span Planar Skeletal Structure
A large span truss model [44] with a length of 25 m and the height of 1m is shown in Figure 12. The elastic modulus is E = 2.0 × 10 5 MPa, the density is 7800 kg/m 3 , and the cross-section area of each bar is 100 mm 2 . In FEM, the large displacement Saint Venant Kirchhoff material model and complete Lagrange description [45] are adopted for the truss elements. There are 101 member bars in total, and each member bar is a single element with hinged constraints at both ends. As shown in Figure 9, Points A and B are located near the middle of the span and are set to gain an initial velocity v = 100 m/s. This example is used to test both the numerical stability and numerical dissipation of the algorithm, and the integration scheme in Equations (34)- (38) with   0 . Firstly, setting the time step = 0.01s t , this problem was solved with trapezoidal rule and the proposed method with  = 0 . The vertical responses at Point A are shown in Figure 13. It can be obviously seen that after about 2 s, the numerical errors of the acceleration start to accumulate significantly for the trapezoidal rule, while the proposed method with  = 0 still performs well. Enlarging the time step size to be = 0.02s t and = 0.04s t , the responses resulted from the proposed method with  = 0 remain stable, as shown in Figure 14, i.e., the proposed method performs better in numerical stability than the trapezoidal rule for this example.

Long Span Planar Skeletal Structure
A large span truss model [44] with a length of 25 m and the height of 1m is shown in Figure 12. The elastic modulus is E = 2.0 × 10 5 MPa, the density is 7800 kg/m 3 , and the cross-section area of each bar is 100 mm 2 . In FEM, the large displacement Saint Venant Kirchhoff material model and complete Lagrange description [45] are adopted for the truss elements. There are 101 member bars in total, and each member bar is a single element with hinged constraints at both ends. As shown in Figure 9, Points A and B are located near the middle of the span and are set to gain an initial velocity v = 100 m/s. This example is used to test both the numerical stability and numerical dissipation of the algorithm, and the integration scheme in Equations (34)

Long Span Planar Skeletal Structure
A large span truss model [44] with a length of 25 m and the height of 1m is shown in Figure 12. The elastic modulus is E = 2.0 × 10 5 MPa, the density is 7800 kg/m 3 , and the cross-section area of each bar is 100 mm 2 . In FEM, the large displacement Saint Venant Kirchhoff material model and complete Lagrange description [45] are adopted for the truss elements. There are 101 member bars in total, and each member bar is a single element with hinged constraints at both ends. As shown in Figure 9, Points A and B are located near the middle of the span and are set to gain an initial velocity v = 100 m/s. This example is used to test both the numerical stability and numerical dissipation of the algorithm, and the integration scheme in Equations (34)- (38) with   0 . Firstly, setting the time step = 0.01s t , this problem was solved with trapezoidal rule and the proposed method with  = 0 . The vertical responses at Point A are shown in Figure 13. It can be obviously seen that after about 2 s, the numerical errors of the acceleration start to accumulate significantly for the trapezoidal rule, while the proposed method with  = 0 still performs well. Enlarging the time step size to be = 0.02s t and = 0.04s t , the responses resulted from the proposed method with  = 0 remain stable, as shown in Figure 14, i.e., the proposed method performs better in numerical stability than the trapezoidal rule for this example. Firstly, setting the time step ∆t = 0.01 s, this problem was solved with trapezoidal rule and the proposed method with α = 0. The vertical responses at Point A are shown in Figure 13. It can be obviously seen that after about 2 s, the numerical errors of the acceleration start to accumulate significantly for the trapezoidal rule, while the proposed method with α = 0 still performs well. Enlarging the time step size to be ∆t = 0.02 s and ∆t = 0.04 s, the responses resulted from the proposed method with α = 0 remain stable, as shown in Figure 14, i.e., the proposed method performs better in numerical stability than the trapezoidal rule for this example. Secondly, in order to verify the efficiency in numerical dissipation of the proposed algorithm, we set the time step interval = 0.001s t and solve the problem with  = 0, 0.01, 0.05, 0.1, 0.2 respectively, the resulting responses are shown in Figures 15 and 16. It can be seen that, when  = 0 , the high-frequency oscillation of the system is very intense. With an increase of  , i.e., numerical damping is introduced, the high-order-frequency response of the system obviously decays. The Secondly, in order to verify the efficiency in numerical dissipation of the proposed algorithm, we set the time step interval = 0.001s t and solve the problem with  = 0, 0.01, 0.05, 0.1, 0.2 respectively, the resulting responses are shown in Figures 15 and 16. It can be seen that, when  = 0 , the high-frequency oscillation of the system is very intense. With an increase of  , i.e., numerical damping is introduced, the high-order-frequency response of the system obviously decays. The Secondly, in order to verify the efficiency in numerical dissipation of the proposed algorithm, we set the time step interval ∆t = 0.001 s and solve the problem with α = 0, 0.01, 0.05, 0.1, 0.2 respectively, the resulting responses are shown in Figures 15 and 16. It can be seen that, when α = 0, the high-frequency oscillation of the system is very intense. With an increase of α, i.e., numerical damping is introduced, the high-order-frequency response of the system obviously decays. The higher α is, the faster the high frequency response decays. In addition, the low frequency response remains intact while no obvious amplitude decay and no period elongation appear. Therefore, the proposed method can effectively filter out the high frequency response of the system while retaining the low frequency response well through the controllable numerical damping.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 16 of 23 higher  is, the faster the high frequency response decays. In addition, the low frequency response remains intact while no obvious amplitude decay and no period elongation appear. Therefore, the proposed method can effectively filter out the high frequency response of the system while retaining the low frequency response well through the controllable numerical damping.
(a) (b) Thirdly, the vertical displacement responses at Point A calculated with the Wilson  -method, Bathe method, and the proposed method are shown and compared in Figures 17 and 18, where 0.01s t = and the computational times are correspondingly listed in Table 3. It is found that, for the Wilson method the computational time is short, but the numerical damping is a little large, which results in excessive dissipation of low frequency modes and significant amplitude decay. For Bathe method, the low frequency modes are retained well and the high frequency modes of the response are quickly filtered out, but the computational time is longer. Nevertheless, with the setting of different α, the present method performs better than the Wilson method in numerical dissipation and better than Bathe method in computational efficiency. higher  is, the faster the high frequency response decays. In addition, the low frequency response remains intact while no obvious amplitude decay and no period elongation appear. Therefore, the proposed method can effectively filter out the high frequency response of the system while retaining the low frequency response well through the controllable numerical damping.
(a) (b) Thirdly, the vertical displacement responses at Point A calculated with the Wilson  -method, Bathe method, and the proposed method are shown and compared in Figures 17 and 18, where 0.01s t = and the computational times are correspondingly listed in Table 3. It is found that, for the Wilson method the computational time is short, but the numerical damping is a little large, which results in excessive dissipation of low frequency modes and significant amplitude decay. For Bathe method, the low frequency modes are retained well and the high frequency modes of the response are quickly filtered out, but the computational time is longer. Nevertheless, with the setting of different α, the present method performs better than the Wilson method in numerical dissipation and better than Bathe method in computational efficiency. Thirdly, the vertical displacement responses at Point A calculated with the Wilson θ-method, Bathe method, and the proposed method are shown and compared in Figures 17 and 18, where ∆t = 0.01 s and the computational times are correspondingly listed in Table 3. It is found that, for the Wilson method the computational time is short, but the numerical damping is a little large, which results in excessive dissipation of low frequency modes and significant amplitude decay. For Bathe method, the low frequency modes are retained well and the high frequency modes of the response are quickly filtered out, but the computational time is longer. Nevertheless, with the setting of different α, the present method performs better than the Wilson method in numerical dissipation and better than Bathe method in computational efficiency. Appl. Sci. 2019, 9,

Cantilever Beam
Consider the cantilever Euler-Bernoulli beam under a uniformly distributed pulse load as the example [5], which is shown in Figure 19. A linear elastic constitutive relation and degenerate Green-Lagrange strain tensor [45] are adopted to describe the finite elements, and 400 elements are used in spatial discretization. Linear Lagrange element interpolation is used for axial displacement, and cubic Hermite element interpolation is used for lateral and rotational displacements.

Cantilever Beam
Consider the cantilever Euler-Bernoulli beam under a uniformly distributed pulse load as the example [5], which is shown in Figure 19. A linear elastic constitutive relation and degenerate Green-Lagrange strain tensor [45] are adopted to describe the finite elements, and 400 elements are used in spatial discretization. Linear Lagrange element interpolation is used for axial displacement, and cubic Hermite element interpolation is used for lateral and rotational displacements.

Cantilever Beam
Consider the cantilever Euler-Bernoulli beam under a uniformly distributed pulse load as the example [5], which is shown in Figure 19. A linear elastic constitutive relation and degenerate Green-Lagrange strain tensor [45] are adopted to describe the finite elements, and 400 elements are used in spatial discretization. Linear Lagrange element interpolation is used for axial displacement, and cubic Hermite element interpolation is used for lateral and rotational displacements.   Solving the beam with different numbers of spatial elements, the corresponding computational time taken in each method is listed in Table 4, and then the variation of computational time versus the number of degrees of freedom (DOFs) is plotted in Figure 23. The computational speed of the proposed method is as good as the Wilson method and is obviously better than the Bathe method.  Solving the beam with different numbers of spatial elements, the corresponding computational time taken in each method is listed in Table 4, and then the variation of computational time versus the number of degrees of freedom (DOFs) is plotted in Figure 23. The computational speed of the proposed method is as good as the Wilson method and is obviously better than the Bathe method.

Conclusions
Taking the nonlinear dynamic Equation (1) as the model problem, a step-by-step time integration method is proposed based on the Galerkin finite element approximation and Newton-Raphson iteration. In summary, the following four conclusions of this method should be noted: 1. For the proposed method based on the Galerkin FEM, usually time elements of degree m = 1,2,3 are used in practical computation. Numerical examples show that such elements can gain a convergence order of 2m for both displacement and velocity results. 2. When a reduced Gauss integral method is applied in computation, the present method is unconditional stability. 3. For linear elements, a feasible time integration scheme with controllable numerical damping has been constructed by modifying the test function and introducing a special integral rule. 4. This is a self-starting and one-step method, which can effectively solve both linear and nonlinear structural dynamic problems. Representative numerical examples have shown its reliability and efficiency, and according to comparisons of results with the trapezoidal rule, Wilson  -method, and Bathe method, the present method performs well in stability with controllable numerical dissipation, and its computational efficiency turns out to be superior as well.
The proposed time integration method is of comprehensive capabilities in nonlinear dynamics. Furthermore, since it originates from FEM, self-adaptive FE techniques can be directly introduced into it, and then an adaptive integration algorithm can be established with the size of time step being automatically produced, leading to more efficient solutions of nonlinear dynamic problems.