Abstract

An energy conservation algorithm for numerically solving nonlinear multidegree-of-freedom (MDOF) dynamic equations is proposed. Firstly, by Taylor expansion and Duhamel integration, an integral iteration formula for numerically solving the nonlinear problems can be achieved. However, this formula still includes a parameter that is to be determined. Secondly, through some mathematical manipulations, the original dynamical equation can be further converted into an energy conservation equation which can then be used to determine the unknown parameter. Finally, an accurate numerical result for the nonlinear problem is achieved by substituting this parameter into the integral iteration formula. Several examples are used to compare the current method with the well-known Runge-Kutta method. They all show that the energy conservation algorithm introduced in this study can eliminate algorithm damping inherent in the Runge-Kutta algorithm and also has better stability for large integral steps.

1. Introduction

Numerical stability and algorithmic damping have been long recognized as two important aspects that need to be carefully handled in time integration algorithms for solving dynamic problems. Indeed, many works have been done in this area. For instance, to improve the stability, the Generalized-𝛼method of Chung and Hulbert [1], the HHT-𝛼method of Hilber et al. [2], and the WBZ-𝛼 method of Wood et al. [3] all demonstrate very good dissipation property either at low-frequency or high-frequency ranges. In references [4–8], Fung presents a series of time-step algorithms that are based on different mathematical and mechanical principles and can be used to deal with linear dynamical problems. Recent works on numerical methods have been focusing on modeling long-term qualitative properties and stabilities in the numerical solutions of nonlinear dynamic problems. In the past, the feature of the energy conservation of a system has been widely used in various numerical integration methods to achieve satisfactory results. Since the property of exact energy conservation enables the numerical scheme to be stable without resorting to high-frequency numerical dissipation [9], the feature of the energy conservation appears particularly attractive in numerically solving nonlinear dynamic systems. Consequently, much effort has gone into the development of energy conserving time-stepping schemes. By adding an additional constraint through the energy conservation equation between adjacent time steps, Bui proposed a modified Newmark family for nonlinear dynamic analysis [10]. Similar works are also reported by LaBudde and Greenspan [11], Hughes et al. [12], Greespan [13], Simo and Tarnow [14], Simo et al. [15], Greenspan [16], and Fung and Chow [17]. Interestingly, all works cited above derived their algorithms starting with the finite difference method, and most of them were based on Hamilton’s canonical equations of motion.

Though energy conservative methods have showed some advantages, they also might suffer from some drawbacks in practical applications. In an attempt to obtain a stable large-step integration, Simo and Gonzalez [9] used the energy-momentum algorithm, which was obtained from the modification of midpoint scheme. By doing this, they risked to wrongly approximate slowly varying solution quantities in highly oscillatory systems, which was especially significant when fast and low modes are tightly coupled [18]. Moreover, with the addition of the constraint on the energy conservation, one has an overdetermined system that the exact solution has to satisfy. However, once the system is discretized and approximated in order to solve it numerically, the overdetermined system may not have a solution anymore. In other words, difficulties may occur in the multidimensional root finding of the corresponding nonlinear system of equations [19].

In this paper, the authors proposed a time integration formula and scheme which can maintain the system energy conservation constraint automatically. However, this method is also different from the energy conservative methods developed previously in the following four aspects. Firstly, the solution of the nonlinear dynamic equation was presented analytically by Duhamel integral in the current method. Secondly, the current method used the Taylor expansion to approximate the exact solution of the nonlinear equation. During this process, an undetermined parameter was introduced. Thirdly, substituting the approximate solution into the analytical solution, an iterative formula with the undetermined parameter was derived. Finally, the energy conservation equation was established and the undetermined parameter was obtained. In practice, how to introduce and obtain the undetermined parameter should be carefully considered case by case, especially when the right hand term of the nonlinear equation includes functions of velocity and time. Another merit of the current algorithm was that it behaved very stable under large time steps. Comparisons with the Rounge-Kutta method showed that the proposed method had much better stabilities in solving different types of nonlinear equations.

The remainder of this paper is arranged as follows. In Section 2, a detailed process of deriving the integral iteration formula is presented. In the iteration formula, we will show that how an undetermined parameter is introduced into the formula. In Section 3, the energy conservation equation is established for obtaining the algebraic expression of the undetermined parameter. Section 4 focuses on the calculation of the multinomial interpolation used in Section 2. Section 5 shows some representative numerical examples which compare the current method with the popular Runge-Kutta method in terms of algorithmic damping and stability. Finally, some discussions and conclusions are given in Section 6.

2. Derivation of the Integral Iteration Formula

Step-by-step time-integration algorithms are commonly used to solve dynamic equations which mostly come from actual engineering problems. By spatial discretization using the finite element method, a nonlinear system may be represented by a second order nonlinear ordinary differential equation asπŒΜˆΜ‡π±+𝐊𝐱=𝐟(𝐱,𝐱,𝑑),(2.1) where 𝐌,𝐊are 𝑛-by-𝑛 constant mass and stiffness matrices, respectively. Μ‡Μˆπ±,𝐱,𝐱,𝐟 are vectors with rank 𝑛 representing the displacement, velocity, and acceleration, respectively. 𝑑 is time. The right hand term in (2.1), Μ‡πŸ(𝐱,𝐱,𝐭), is the force vector that includes all external forces such as the damping forces and the nonlinear forces. Using the matrix decomposition, the mass matrix 𝐌 can be expressed as 𝐌=π‹β‹…π‹π‘‡πŒβˆ’1=π‹βˆ’π‘‡β‹…π‹βˆ’1.(2.2)

Substituting (2.2) into (2.1) and multiplying by π‹βˆ’1, we can obtainπ‹π‘‡Μˆπ±+π‹βˆ’1πŠπ‹βˆ’π‘‡β‹…π‹π‘‡π±=π‹βˆ’1𝐟.(2.3) Knowing that 𝐋,π‹βˆ’1,𝐋𝑇,π‹βˆ’π‘‡ are all constant matrices, a variable substitution can be executed as𝐲=𝐋𝑇̇𝐱,𝐲=π‹π‘‡Μ‡Μˆπ±,𝐲=π‹π‘‡Μˆπ±.(2.4) Substituting (2.4) into (2.3), a new dynamic system which is equivalent to the original system can be obtained as̈𝐲+𝐊𝐲=𝐅(𝑑),𝐅(𝑑)=π‹βˆ’1𝐟(𝑑),𝐊=π‹βˆ’1πŠπ‹βˆ’π‘‡.(2.5) Let𝐊=𝐊0+𝐃,(2.6) where 𝐃 is a diagonal matrix and its diagonal elements are the diagonal elements of matrix 𝐊.𝐊0 is a matrix whose diagonal elements are zero, and other elements are equal to those in 𝐊. Note that 𝐷𝑖𝑖β‰₯0 for 𝐌 which is a positive definite matrix and 𝐊 is a semipositive definite matrix. Let𝐷𝑖𝑖=𝑑2𝑖.(2.7) Then (2.5) can be rewritten as̈𝐲+𝐃𝐲=𝐅(𝑑)βˆ’πŠ0𝐲.(2.8) The separate form of the above matrix equation can be expressed as follows:Μˆπ‘¦π‘–+𝑑2𝑖𝑦𝑖=𝐹𝑖(𝑑)βˆ’π€π‘–0𝐲,𝑖=1,2,…,𝑛,(2.9) where 𝐀𝑖0 is the 𝑖th row of the matrix 𝐊0. From (2.9), it is easy to see that the analytical solution of the displacement and the velocity can be obtained by the Duhamel integral as𝑦𝑖(𝑑)=π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έcosπ‘‘π‘–ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ+Μ‡π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έsinπ‘‘π‘–ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έπ‘‘π‘–+1π‘‘π‘–ξ€œπ‘‘π‘‘π‘˜πΉπ‘–(𝜍)sin𝑑𝑖1(π‘‘βˆ’πœ)π‘‘πœβˆ’π‘‘π‘–ξ€œπ‘‘π‘‘π‘˜π€π‘–0𝐲(𝜍)sin𝑑𝑖(π‘‘βˆ’πœ)π‘‘πœ,(2.10)̇𝑦𝑖(𝑑)=βˆ’π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έπ‘‘π‘–sinπ‘‘π‘–ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ+Μ‡π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έcosπ‘‘π‘–ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ+ξ€œπ‘‘π‘‘π‘˜πΉπ‘–(𝜍)cos𝑑𝑖(ξ€œπ‘‘βˆ’πœ)π‘‘πœβˆ’π‘‘π‘‘π‘˜π€π‘–0𝐲(𝜍)cos𝑑𝑖(π‘‘βˆ’πœ)π‘‘πœ.(2.11)

In order to derive the time integral formula, let 𝑑=π‘‘π‘˜+𝜏in (2.10) and (2.11), where 𝜏 is the integral time step, then we haveπ‘¦π‘–ξ€·π‘‘π‘˜ξ€Έ+𝜏=π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έcosπ‘‘π‘–πœ+Μ‡π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έsinπ‘‘π‘–πœπ‘‘π‘–+1π‘‘π‘–ξ€œπ‘‘π‘˜π‘‘+πœπ‘˜πΉπ‘–(𝜍)sinπ‘‘π‘–ξ€·π‘‘π‘˜ξ€Έβˆ’1+πœβˆ’πœπ‘‘πœπ‘‘π‘–ξ€œπ‘‘π‘˜π‘‘+πœπ‘˜π€π‘–0𝐲(𝜍)sinπ‘‘π‘–ξ€·π‘‘π‘˜ξ€Έ+πœβˆ’πœπ‘‘πœ,(2.12)Μ‡π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έ+𝜏=βˆ’π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έπ‘‘π‘–sinπ‘‘π‘–πœ+Μ‡π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έcosπ‘‘π‘–ξ€œπœ+π‘‘π‘˜π‘‘+πœπ‘˜πΉπ‘–(𝜍)cosπ‘‘π‘–ξ€·π‘‘π‘˜ξ€Έβˆ’ξ€œ+πœβˆ’πœπ‘‘πœπ‘‘π‘˜π‘‘+πœπ‘˜π€π‘–0𝐲(𝜍)cosπ‘‘π‘–ξ€·π‘‘π‘˜ξ€Έ+πœβˆ’πœπ‘‘πœ.(2.13)

Inspecting (2.12) and (2.13), it can be found that there are still some unknown parameters that need to be identified before the time integral formula can be carried out numerically. They are the right hand side terms consisting of the undetermined variables 𝐲(𝑑) and 𝐹𝑖(𝑑). The latter one may also be a function of 𝐲(𝑑) and ̇𝐲(𝑑). To proceed, the Taylor expansion formula is used to expand 𝐲(𝑑) on the interval (π‘‘π‘˜β‰€π‘‘β‰€π‘‘π‘˜+𝜏) as𝑑𝐲(𝑑)=π²π‘˜ξ€Έ+ξ€·π‘‘βˆ’π‘‘π‘˜ξ€ΈΜ‡π²ξ€·π‘‘π‘˜ξ€Έ+ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ22Μˆπ²ξ€·π‘‘π‘˜ξ€Έ+ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ36𝐚.(2.14)

In (2.14), 𝐲(𝑑) is expanded to the third order term following exactly the Taylor expansion process, while in calculating the fourth order term, a new variable vector 𝐚is introduced. It should be pointed out that although the Taylor expansion is an approximation to the original variable 𝐲(𝑑), (2.14) is still an exact expression of the variable 𝐲(𝑑). This is because the last term in (2.14) which includes the newly introduced vector 𝐚 can be interpreted to compensate for the difference between 𝐲(𝑑) and the summation of the first three terms in (2.14). In order to express the vector 𝐚 by the variable 𝐲(𝑑), let 𝐲1=𝐲(𝑑=π‘‘π‘˜+𝜏),𝐲0=𝐲(𝑑=π‘‘π‘˜). From (2.14), it can be obtained that6𝐚=𝜏3𝐲1βˆ’π²0Μ‡π²βˆ’πœ0βˆ’πœ22̈𝐲0ξ‚Ά,𝐲(𝑑)=𝑁1(𝑑)𝐲0+𝑁2̇𝐲(𝑑)0+𝑁3̈𝐲(𝑑)0+𝑁4(𝑑)𝐲1,(2.15) where:𝑁1ξ€·(𝑑)=1βˆ’π‘‘βˆ’π‘‘π‘˜ξ€Έ3𝜏3,𝑁2(𝑑)=π‘‘βˆ’π‘‘π‘˜βˆ’ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ3𝜏2,𝑁3(𝑑)=π‘‘βˆ’π‘‘π‘˜ξ€Έ22βˆ’ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ32𝜏,𝑁4(𝑑)=π›½π‘‘βˆ’π‘‘π‘˜ξ€Έ3ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ+ξ€·βˆ’πœπ‘‘βˆ’π‘‘π‘˜ξ€Έ3𝜏3,(2.16) where an undetermined parameter 𝛽 has been introduced in 𝑁4(𝑑) to regulate the stability of the algorithm and will be determined by the energy conservation equation in the next section. By 𝑑=π‘‘π‘˜,𝑑=π‘‘π‘˜+𝜏 in (2.14), we can obtainπ²ξ€·π‘‘π‘˜ξ€Έ=𝐲0,Μ‡π²ξ€·π‘‘π‘˜ξ€Έ=̇𝐲0,Μˆπ²ξ€·π‘‘π‘˜ξ€Έ=̈𝐲0𝑑,π²π‘˜ξ€Έ+𝜏=𝐲1.(2.17)

By means of multinomial interpolation, 𝐟(𝑑)(π‘‘π‘˜β‰€π‘‘β‰€π‘‘π‘˜+𝜏) can be written as𝐟(𝑑)=𝐫0+ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έπ«1+ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ2𝐫2+ξ€·π‘‘βˆ’π‘‘π‘˜ξ€Έ3𝐫3ξ‚€ξ€·+π‘œπ‘‘βˆ’π‘‘π‘˜ξ€Έ4.(2.18)

Here we use the third order interpolation, and generally one can choose the order of interpolation discretionarily based on solely the algorithm accuracy order that is needed. Different interpolation order will lead to different integration formulas. Now, we will derive the integration formulas first, and the discussion of the interpolation will be addressed at Section 4. Using (2.5) and (2.14), the last two terms of the right hand of (2.11) can be expressed separately as1π‘‘π‘–ξ€œπ‘‘π‘˜π‘‘+πœπ‘˜π€π‘–0𝐲(𝜍)sinπ‘‘π‘–ξ€·π‘‘π‘˜ξ€Έ+πœβˆ’πœπ‘‘πœ=𝛼𝑖0𝐀𝑖0𝐲0+𝛼𝑖1𝐀𝑖0̇𝐲0+𝛼𝑖2𝐀𝑖0̈𝐲0+𝛼𝑖3𝐀𝑖0𝐲1,1π‘‘π‘–ξ€œπ‘‘π‘˜π‘‘+πœπ‘˜π‹π‘–βˆ’1𝐟(𝜍)sinπ‘‘π‘–ξ€·π‘‘π‘˜ξ€Έ+πœβˆ’πœπ‘‘πœ=𝛾𝑖0π‹π‘–βˆ’1𝐫0+𝛾𝑖1π‹π‘–βˆ’1𝐫1+𝛾𝑖2π‹π‘–βˆ’1𝐫2+𝛾𝑖3π‹π‘–βˆ’1𝐫3,(2.19) where π›Όπ‘–π‘˜,π›Ύπ‘–π‘˜,π‘˜=0,1,2,3 are scalar and can be obtained by follow polynomials𝛼𝑖0ξ€·=βˆ’6sinπ‘‘π‘–πœ+𝑑3π‘–πœ3cosπ‘‘π‘–πœβˆ’6π‘‘π‘–πœξ€Έπ‘‘5π‘–πœ3,𝛼𝑖1=ξ€·6π‘‘π‘–πœβˆ’π‘‘2π‘–πœ2sinπ‘‘π‘–πœβˆ’6sinπ‘‘π‘–πœξ€Έπ‘‘5π‘–πœ2,𝛼𝑖2=ξ€·2π‘‘π‘–πœβˆ’3sinπ‘‘π‘–πœ+π‘‘π‘–πœcosπ‘‘π‘–πœξ€Έπ‘‘5π‘–πœ,𝛼𝑖3=ξ€·6sinπ‘‘π‘–πœ+𝑑3π‘–πœ3βˆ’6π‘‘π‘–πœξ€Έπ‘‘5π‘–πœ3βˆ’π›½ξ€·π‘‘2π‘–πœ2βˆ’4+4cosπ‘‘π‘–πœ+π‘‘π‘–πœsinπ‘‘π‘–πœξ€Έπ‘‘6π‘–πœ4,𝛾𝑖0=ξ€·1βˆ’cosπ‘‘π‘–πœξ€Έπ‘‘2𝑖,𝛾𝑖1=ξ€·π‘‘π‘–πœβˆ’sinπ‘‘π‘–πœξ€Έπ‘‘3𝑖,𝛾𝑖2=𝑑2π‘–πœ2βˆ’2+2cosπ‘‘π‘–πœξ€Έπ‘‘4𝑖,𝛾𝑖3=𝑑3π‘–πœ3βˆ’6π‘‘π‘–πœ+6sinπ‘‘π‘–πœξ€Έπ‘‘5𝑖.(2.20) Substituting (2.20) into (2.12), we haveπ‘¦π‘–ξ€·π‘‘π‘˜ξ€Έ+𝜏=π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έcosπ‘‘π‘–πœ+Μ‡π‘¦π‘–ξ€·π‘‘π‘˜ξ€Έsinπ‘‘π‘–πœπ‘‘π‘–+𝛾𝑖0π‹π’βˆ’1𝐫0+𝛾𝑖1π‹π’βˆ’1𝐫1+𝛾𝑖2π‹π’βˆ’1𝐫2+𝛾𝑖3π‹π’βˆ’1𝐫3βˆ’ξ€·π›Όπ‘–0𝐀𝑖0𝐲0+𝛼𝑖1𝐀𝑖0̇𝐲0+𝛼𝑖2𝐀𝑖0̈𝐲0+𝛼𝑖3𝐀𝑖0𝐲1ξ€Έ(𝑖=1,2,…,𝑛).(2.21) Knowing that 𝐲1=𝐲(𝑑=π‘‘π‘˜+𝜏), so we can write the above equation in a matrix form:𝐲1=𝐔0𝐲0+𝐔1̇𝐲0+𝜸0π‹βˆ’1𝐫0+πœΈπŸπ‹βˆ’1𝐫1+𝜸2π‹βˆ’1𝐫2+𝜸3π‹βˆ’1𝐫3βˆ’πœΆ0𝐊0𝐲0βˆ’πœΆ1𝐊0̇𝐲0βˆ’πœΆ2𝐊0̈𝐲0βˆ’πœΆ3𝐊0𝐲1,(2.22) where 𝐔0, 𝐔1,πœΆπ‘˜,πœΈπ‘˜(π‘˜=0,1,2,3) are diagonal matrices and their diagonal elements are cosπ‘‘π‘–πœ, sinπ‘‘π‘–πœ/𝑑𝑖, π›Όπ‘–π‘˜, π›Ύπ‘–π‘˜(𝑖=1,2,…,𝑛,π‘˜=0,1,2,3), respectively. From (2.22), the iterative solution 𝑦1 can be expressed asξ€·πˆ+𝜢3𝐊0𝐲1=𝐔0βˆ’πœΆ0𝐊0𝐲0+𝐔1βˆ’πœΆ1𝐊0̇𝐲0βˆ’πœΆ2𝐊0̈𝐲0+𝜸0π‹βˆ’1π‘Ÿ0+𝜸1π‹βˆ’1π‘Ÿ1+𝜸2π‹βˆ’1π‘Ÿ2+𝜸3π‹βˆ’1π‘Ÿ3.(2.23) Multiplying (2.23) by π‹πœΆ3βˆ’1, we haveπ‹ξ€·πœΆ3βˆ’1+𝐊0𝐲1=π‹πœΆ3βˆ’1𝐔0βˆ’πœΆ0𝐊0𝐲0+π‹πœΆ3βˆ’1𝐔1βˆ’πœΆ1𝐊0̇𝐲0βˆ’π‹πœΆ3βˆ’1𝜢2𝐊0̈𝐲0+π‹πœΆ3βˆ’1𝜸0π‹βˆ’1𝐫0+𝜸1π‹βˆ’1𝐫1+𝜸2π‹βˆ’1𝐫2+𝜸3π‹βˆ’1𝐫3.(2.24)

Note that 𝐊0 is a matrix whose diagonal elements are zero, and other elements are equal to those in 𝐊 and πœΆπ‘˜ is a diagonal matrix, so 𝐋(𝜢3βˆ’1+𝐊0) is a symmetrical matrix. The left hand side of (2.24) can be written asπ‹ξ€·πœΆ3βˆ’1+𝐊0ξ€Έπ‹π‘‡π‹βˆ’π‘‡π²1=ξ€·π‹πœΆ3βˆ’1𝐋𝑇+π‹πŠ0π‹π‘‡ξ€Έπ‹βˆ’π‘‡π²1.(2.25) From (2.5) and (2.6), we haveπ‹πŠ0𝐋𝑇=π‹β‹…ξ‚πŠβˆ’πƒβ‹…π‹π‘‡=π‹β‹…π‹βˆ’1πŠπ‹βˆ’π‘‡β‹…π‹π‘‡βˆ’π‹πƒπ‹π‘‡=πŠβˆ’π‹πƒπ‹π‘‡.(2.26) From (2.3) and (2.26), the left hand of (2.24) can be written asπ‹ξ€·πœΆ3βˆ’1+𝐊0𝐲1=ξ€Ίξ€·πœΆπŠ+𝐋3βˆ’1ξ€Έπ‹βˆ’πƒπ‘‡ξ€»β‹…π±1.(2.27) All the terms in the right hand side of (2.24) can be expressed separately asπ‹πœΆ3βˆ’1𝐔0βˆ’πœΆ0𝐊0𝐲0=ξ€Ίπ‹ξ€·πœΆ3βˆ’1𝐔0+𝜢3βˆ’1𝜢0πƒξ€Έπ‹π‘‡βˆ’π‹πœΆ3βˆ’1𝜢0π‹βˆ’1πŠξ€»β‹…π±0,π‹πœΆ3βˆ’1𝐔1βˆ’πœΆ1𝐊0̇𝐲0=ξ€Ίπ‹ξ€·πœΆ3βˆ’1𝐔1+𝜢3βˆ’1𝜢1πƒξ€Έπ‹π‘‡βˆ’π‹πœΆ3βˆ’1𝜢1π‹βˆ’1πŠξ€»β‹…Μ‡π±0,π‹πœΆ3βˆ’1𝜢2𝐊0̈𝐲0=π‹πœΆ3βˆ’1𝜢2π‹βˆ’1π‹πŠ0π‹π‘‡π‹βˆ’π‘‡Μˆπ²0=π‹πœΆ3βˆ’1𝜢2π‹βˆ’1ξ€·πŠβˆ’π‹πƒπ‹π‘‡ξ€Έβ‹…Μˆπ±0.(2.28) After substituting proper variables, the iteration formula in terms of the original variables can be expressed asξ€Ίξ€·πœΆπŠ+𝐋3βˆ’1ξ€Έπ‹βˆ’πƒπ‘‡ξ€»β‹…π±1=ξ€Ίπ‹ξ€·πœΆ3βˆ’1𝐔0+𝜢3βˆ’1𝜢0πƒξ€Έπ‹π‘‡βˆ’π‹πœΆ3βˆ’1𝜢0π‹βˆ’1πŠξ€»β‹…π±0+ξ€Ίπ‹ξ€·πœΆ3βˆ’1𝐔1+𝜢3βˆ’1𝜢1πƒξ€Έπ‹π‘‡βˆ’π‹πœΆ3βˆ’1𝜢1π‹βˆ’1πŠξ€»β‹…Μ‡π±0βˆ’π‹πœΆ3βˆ’1𝜢2π‹βˆ’1ξ€·πŠβˆ’π‹πƒπ‹π‘‡ξ€Έβ‹…Μˆπ±0+π‹πœΆ3βˆ’1𝜸0π‹βˆ’1𝐫0+π‹πœΆ3βˆ’1𝜸1π‹βˆ’1𝐫1+π‹πœΆ3βˆ’1𝜸2π‹βˆ’1𝐫2+π‹πœΆ3βˆ’1𝜸3π‹βˆ’1𝐫3.(2.29)

In (2.29), there is a term consisting of a double derivative. According to (2.3), the term with the double derivative ̈𝐱0 can be replaced by βˆ’π‹βˆ’π‘‡π‹βˆ’1𝐊𝐱0+π‹βˆ’π‘‡π‹βˆ’1𝐟0. Finally, the displacement iteration formula can be obtained asξ€Ίξ€·πœΆπŠ+𝐋3βˆ’1ξ€Έπ‹βˆ’πƒπ‘‡ξ€»β‹…π±1𝜢=𝐋⋅3βˆ’1𝐔0+𝜢3βˆ’1𝜢0πƒξ€Έπ‹π‘‡βˆ’πœΆ3βˆ’1𝜢0π‹βˆ’1𝐊+𝜢3βˆ’1𝜢2π‹βˆ’1ξ€·πŠβˆ’π‹πƒπ‹π‘‡ξ€Έπ‹βˆ’π‘‡π‹βˆ’1πŠξ€»β‹…π±0𝜢+𝐋⋅3βˆ’1𝐔1+πœΆπŸ‘βˆ’1𝜢1πƒξ€Έπ‹π‘‡βˆ’πœΆ3βˆ’1𝜢1π‹βˆ’1πŠξ€»β‹…Μ‡π±0ξ€ΊπœΆ+𝐋⋅3βˆ’1𝜸0π‹βˆ’1βˆ’πœΆ3βˆ’1𝜢2π‹βˆ’1ξ€·πŠβˆ’π‹πƒπ‹π‘‡ξ€Έβ‹…π‹βˆ’π‘‡π‹βˆ’1𝐫0+π‹πœΆ3βˆ’1𝜸1π‹βˆ’1𝐫1+π‹πœΆ3βˆ’1𝜸2π‹βˆ’1𝐫2+π‹πœΆ3βˆ’1𝜸3π‹βˆ’1𝐫3.(2.30)

Substituting (2.14) and (2.15) into (2.13) and through some mathematical manipulations, the velocity iteration formula can be obtained aṡ𝐱1=π‘πœπ±1+𝐑0𝐱0+𝐑1̇𝐱0+π‹βˆ’π‘‡ξ€ΊπœΌ0π‹βˆ’1βˆ’ξ€·πœ2π‹βˆ’1πŠβˆ’πœ2πƒπ‹π‘‡ξ€Έπ‹βˆ’π‘‡π‹βˆ’1𝐫0+π‹βˆ’π‘‡πœΌ1π‹βˆ’1𝐫1+π‹βˆ’π‘‡πœΌ2π‹βˆ’1𝐫2+π‹βˆ’π‘‡πœΌ3π‹βˆ’1𝐫3,(2.31) where𝜼0, 𝜼1, 𝜼2, 𝜼3 are diagonal matrices with diagonal elements πœ‚π‘–0, πœ‚π‘–1, πœ‚π‘–2, πœ‚π‘–3, respectively. πœπ‘š is also a diagonal matrix and its diagonal elements are π‘π‘–π‘š(π‘š=0,1,2,3). Every term in the right hand side of (2.31) is given as bellow: 𝐑0=π‹βˆ’π‘‡π•ξ€Ίξ€·0𝐋𝑇+𝐜0πƒπ‹π‘‡βˆ’πœ0π‹βˆ’1πŠξ€Έ+πœπŸπ‹βˆ’1πŠβˆ’πœπŸπƒπ‹π“ξ€Έπ‹βˆ’π“π‹βˆ’1π‘β‹…πŠ,1=π‹βˆ’π‘‡ξ€·π•1𝐋T+𝐜1πƒπ‹π‘‡βˆ’πœ1π‹βˆ’1πŠξ€Έ,π‘πœ=βˆ’π‹βˆ’π‘‡πœ3ξ€·π‹βˆ’1πŠβˆ’πƒπ‹π‘‡ξ€Έ,𝑐𝑖0=𝑑3π‘–πœ3sinπ‘‘π‘–πœβˆ’6cosπ‘‘π‘–πœβˆ’3𝑑2π‘–πœ2ξ€Έ+6𝑑4π‘–πœ3,𝑐𝑖1𝑑=βˆ’2π‘–πœ2cosπ‘‘π‘–πœ+6cosπ‘‘π‘–πœ+2𝑑2π‘–πœ2ξ€»βˆ’6𝑑4π‘–πœ2,𝑐𝑖2(ξ€·=βˆ’1/2)2π‘‘π‘–πœsinπ‘‘π‘–πœ+6cosπ‘‘π‘–πœ+𝑑2π‘–πœ2ξ€Έβˆ’6𝑑4π‘–πœ,𝑐𝑖3=ξ€·6cosπ‘‘π‘–πœ+3𝑑2π‘–πœ2ξ€Έβˆ’6𝑑4π‘–πœ3βˆ’π›½ξ€·6π‘‘π‘–πœcosπ‘‘π‘–πœβˆ’24sinπ‘‘π‘–πœβˆ’π‘‘2π‘–πœ2ξ€Έ+18𝑑5π‘–πœ4,πœ‚π‘–0=sinπ‘‘π‘–πœπ‘‘π‘–πœ‚π‘–1=1βˆ’cosπ‘‘π‘–πœπ‘‘2π‘–πœ‚π‘–2=2π‘‘π‘–πœβˆ’2sinπ‘‘π‘–πœπ‘‘3𝑖,πœ‚π‘–3=ξ€œπ‘‘π‘˜π‘‘+πœπ‘˜ξ€·πœβˆ’π‘‘π‘˜ξ€Έ3cosπ‘‘π‘–ξ€·π‘‘π‘˜ξ€Έ+πœβˆ’πœπ‘‘πœ,(2.32) where 𝐕0, 𝐕1 are diagonal matrices with diagonal elements βˆ’π‘‘π‘–sinπ‘‘π‘–πœ, cosπ‘‘π‘–πœ, respectively.

3. Energy Conservation Equation

One reason of expressing the dynamic equation in the form of (2.1) is to establish the energy conservation equation more conveniently and more directly. The following steps illustrate the construction of the energy conservation equation. Multiplying ̇𝐱𝑇 to both sides of (2.1), we haveΜ‡π±π‘‡πŒΜˆΜ‡π±π±+π‘‡Μ‡π±πŠπ±=π‘‡Μ‡πŸ(𝐱,𝐱,𝑑).(3.1)

Integrating (2.32) from π‘‘π‘˜ to π‘‘π‘˜+1, we can obtainξ€œπ‘‘π‘˜+1π‘‘π‘˜Μ‡π±π‘‡πŒΜˆξ€œπ±π‘‘π‘‘+π‘‘π‘˜+1π‘‘π‘˜Μ‡π±π‘‡ξ€œπŠπ±π‘‘π‘‘=π‘‘π‘˜+1π‘‘π‘˜Μ‡π±π‘‡πŸπ‘‘π‘‘π‘˜=0,1,2,3,…,(3.2) namely,12Μ‡π±π‘‡πŒΜ‡π±|||π‘‘π‘˜+1π‘‘π‘˜+12𝐱𝑇|||πŠπ±π‘‘π‘˜+1π‘‘π‘˜=ξ€œπ‘‘π‘˜+1π‘‘π‘˜Μ‡π±π‘‡πŸπ‘‘π‘‘π‘˜=0,1,2,3,….(3.3)

Simplifying (3.3), we can obtain an energy conservation equation between π‘‘π‘˜ and π‘‘π‘˜+1as follows:π‘‡π‘˜+1βˆ’π‘‡π‘˜+π‘‰π‘˜+1βˆ’π‘‰π‘˜=ξ€œπ‘‘π‘˜+1π‘‘π‘˜Μ‡π±π‘‡πŸπ‘‘π‘‘π‘˜=0,1,2,3,…,(3.4) whereπ‘‡π‘˜+1=12Μ‡π±π‘‡π‘˜+1πŒΜ‡π±π‘˜+1,π‘‡π‘˜=12Μ‡π±π‘‡π‘˜πŒΜ‡π±π‘˜,π‘‰π‘˜+1=12π±π‘‡π‘˜+1πŠπ±π‘˜+1,π‘‰π‘˜=12π±π‘‡π‘˜πŠπ±π‘˜.(3.5)

Substituting (2.30) and (2.31) into the left hand side of (3.4), a polynomial of the undermined parameter𝛽can be easily achieved. For the right hand side of (3.4), the integral term can be firstly decomposed into two parts as follows:𝑓1=ξ€œπ‘‘π‘˜+1π‘‘π‘˜Μ‡π±π‘‡πͺ1(𝐱)𝑑𝑑,𝑓2=ξ€œπ‘‘π‘˜+1π‘‘π‘˜Μ‡π±π‘‡πͺ2(̇𝐱,𝑑)𝑑𝑑,(3.6) where the term 𝑓1 is an integral of an autonomous system and can be integrated easily. The ̇𝐱(𝑑) in term 𝑓2 can be expressed as a polynomial of time using the relationship in (2.4) and taking the derivative of (2.15) with respect to time. Two predictive methods are recommended for determining the unknown term 𝐲1 in (2.15). One is to let 𝐲1=𝐲𝟎+Μ‡π²πŸŽπ‰ and the other is to let 𝛽=1 in (2.23).

Then through (3.5), an algebraic equation with an undetermined parameter 𝛽 can be established and 𝛽 can be numerically obtained by the Newton iteration method or other algebraic methods. Finally, substituting 𝛽 into (2.30) and (2.31), a numerical result can then be achieved.

4. Calculations of the Interpolation

Before giving some numerical examples, choosing the proper interpolation form of (2.18) must be discussed because it will affect the accuracy and stability of the proposed algorithm. In the current study, the authors use the Hermite interpolation to approximate the 𝐫1,𝐫2,𝐫3,𝐫4 in (2.18), that is,⎧βŽͺβŽͺβŽͺβŽͺ⎨βŽͺβŽͺβŽͺβŽͺβŽ©π‘Ÿπ‘–0π‘Ÿπ‘–1π‘Ÿπ‘–2π‘Ÿπ‘–3⎫βŽͺβŽͺβŽͺβŽͺ⎬βŽͺβŽͺβŽͺβŽͺ⎭=βŽ‘βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ’βŽ£βˆ’11000𝜏01𝜏0βˆ’3𝜏2βˆ’2𝜏3𝜏2βˆ’1𝜏2𝜏31𝜏2βˆ’2𝜏31𝜏2⎀βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ₯βŽ¦β‹…βŽ§βŽͺβŽͺβŽͺ⎨βŽͺβŽͺβŽͺβŽ©π‘“π‘–ξ€·0,π‘₯𝑖(0),Μ‡π‘₯𝑖𝑓(0)′𝑖0,π‘₯𝑖(0),Μ‡π‘₯𝑖𝑓(0)𝑖𝑑,π‘₯𝑖(𝜏),Μ‡π‘₯𝑖𝑓(𝜏)ξ…žπ‘–ξ€·π‘‘,π‘₯𝑖(𝜏),Μ‡π‘₯π‘–ξ€ΈβŽ«βŽͺβŽͺβŽͺ⎬βŽͺβŽͺβŽͺ⎭(𝜏).(4.1)

It should be noted that there are unknowns in the right hand term of (4.1) which are π‘₯𝑖(𝜏),Μ‡π‘₯𝑖(𝜏). The prediction of the two unknowns is shown below. For example, we can letπ‘₯𝑖(𝜏)=π‘₯𝑖(0)+Μ‡π‘₯𝑖(0)β‹…πœ,Μ‡π‘₯𝑖(𝜏)=Μ‡π‘₯𝑖(0)+̈π‘₯𝑖(0)β‹…πœ.(4.2)

Then at every iteration of (3.4), the parameter 𝜷 can be updated. Submitting 𝜷 into (2.30) and (2.31), a new prediction of the displacement and velocity can therefore be obtained.

5. Numerical Examples

In this section we give some numerical examples to verify the effectiveness of the proposed algorithm, in particular, the advantage in stability of the proposed algorithm. Since (2.14) is a fourth order Taylor expansion, the energy conservation algorithm has fourth order accuracy. So we choose the Rounge-Kutta method as a numerical comparison. The numerical results show the advantages of the proposed energy conservation algorithm in terms of its integration stability and the ability to eliminate the algorithm damping inherent in the Rounge-Kutta method.

5.1. The Oscillation of a Nonlinear Simple Pendulum

The dynamic equation of a nonlinear single pendulum without damping can be written as̈π‘₯+πœ”20sinπ‘₯=0,πœ”20=1.0,π‘₯(0)=1.57,(5.1) where π‘₯ denotes the angular displacement. The numerical solutions are shown in Figure 1. From the figure we can see that the proposed energy conservation method (ECM) can keep the numerical stability and have no computing damping under large-step comparing with the Rounge-Kutta (RK) method. The numerical result of parameter 𝛽 is shown in Figure 2. Table 1 gives the comparison of the computing efficiency. The efficiency of the proposed method is not as good as the RK method due to the iteration of parameter 𝛽 and the time needed to compute the associated matrices. Figure 3 gives the error analysis between the ECM and the RK under time step 1.0 s.

5.2. The Unforced Linear Vibration of the Cuboid Rigid Body with Two DOF

The structural diagram of the system is shown in Figure 4. The mass of the rigid body is π‘š and the length of the hemline is π‘Ž. The center of mass is collocated at the geometry center (point C). The mass moment of inertia around the center of mass is 𝐽 and the stiffness of the spring is π‘˜. The deformations of the two springs are π‘₯1,π‘₯2. The displacement in vertical direction of the center of mass is π‘₯𝑐. The angular displacement of the rigid body about the mass center is πœ™. Using the above parameters, the equation of motion of the system can be written asβŽ‘βŽ’βŽ’βŽ’βŽ£π‘š4+π½π‘Ž2π‘š4βˆ’π½π‘Ž2π‘š4βˆ’π½π‘Ž2π‘š4+π½π‘Ž2⎀βŽ₯βŽ₯βŽ₯⎦⎧βŽͺ⎨βŽͺ⎩̈π‘₯1̈π‘₯2⎫βŽͺ⎬βŽͺ⎭+⎑⎒⎒⎣⎀βŽ₯βŽ₯⎦⎧βŽͺ⎨βŽͺ⎩π‘₯π‘˜00π‘˜1π‘₯2⎫βŽͺ⎬βŽͺ⎭=⎧βŽͺ⎨βŽͺ⎩00⎫βŽͺ⎬βŽͺ⎭.(5.2) Let π‘š=8,π‘Ž=1,π‘˜=2,𝐽=1,π‘₯1=1,andπ‘₯2=βˆ’1. Figures 5 and 6 compare the displacement (π‘₯1) and velocity (π‘₯2-dot) results predicted by the proposed method and the RK method. It can be seen that even with a big time step 1.0 s, the proposed method still has an accurate numerical solution but the RK method does not. Table 2 gives the comparison of the computing efficiency. Again, the efficiency of the proposed method is lower than that of the RK method in calculating these two degrees of freedom problem. Furthermore, it is noticed that the RK method almost keeps the same efficiency in Sections 5.1 and 5.2.

Figure 7 gives the error analysis between the ECM and the RK under time step 1.0 s. As Figure 7 already shows that the accuracy of the ECM is almost same as the result of RK with a 0.001 time step, the comparison does not use the RK with a small time step.

5.3. The Unforced Nonlinear Oscillation of a Spring Pendulum with Two DOF

The dynamic equation of the spring pendulum can be written as̈π‘₯1+2𝑐1Μ‡π‘₯1+πœ”21π‘₯1βˆ’π‘1π‘₯1π‘₯2=0,̈π‘₯2+2𝑐2Μ‡π‘₯2+πœ”22π‘₯2βˆ’π‘2π‘₯21=0.(5.3)

Figures 8 and 9 show the numerical solution under different damping. Parameters and initial condition are given as follows:πœ”1=1.0,πœ”2=1.5,𝑏1=𝑏2=1.0,π‘₯1=π‘₯2=0.1,Μ‡π‘₯1=Μ‡π‘₯2=0.0.(5.4)

Figure 10 shows the comparison of the numerical results between ECM and RK methods under large time steps. It is obvious that the proposed method can eliminate algorithm damping better and provides better stability than the RK scheme. Parameters and initial conditions used in the calculation are as follows:𝑐1=𝑐2=0.0,πœ”1=1.0,πœ”2𝑏=1.5,1=𝑏2=1.0,π‘₯1=π‘₯2=0.1,Μ‡π‘₯1=Μ‡π‘₯2=0.0.(5.5)

Figure 11 shows long-time response of (5.2). Parameter and initial condition is as same as (5.4). From this figure we can see that after long-term iteration the proposed method still keeps numerical stability under a relatively large time step. But the algorithm damping makes the RK lose accuracy. Table 3 gives the comparison of the computing efficiency between RK and the ECM algorithms. It is shown that, in solving the two degrees of freedom nonlinear problem, the efficiency of the EMC is lower than the RK method. Moreover, comparing the elapsed time by the ECM in Table 2 and Table 3, it also can be seen that the computation efficiency of the ECM is worse for calculating nonlinear problems than for linear problems. Figure 12 shows the error analysis between the ECM and the RK under time step 1.0 s.

6. Conclusion

(1) The energy conservation algorithm has the advantage in stability and time step compared with some numerical means because the numerical solution has been corrected by the energy conservation equation.

(2) All examples have shown that the energy conservation method can eliminate algorithm damping. It is also an effective means for calculating the long-term characteristics of nonlinear dynamic systems.

(3) The proposed method conserves the angular momentum automatically. Although the efficiency of the energy conservation method is not as good as the RK algorithm as well as some other numerical methods discussed in the literature, the integration step is large enough to implement long-term integration with good numerical stability.

(4) The reason of the low efficiency of the proposed method is because the iterations need to calculate the parameter 𝛽 and the time consumed in matrix computing needed by the algorithm. The efficiency of the EMC is lower in dealing with nonlinear problems compared with linear problems.

Acknowledgments

This work was funded by the β€œ973” National Basic Research Project of China (no. Q10110919), Key Project of the National Natural Science Foundation of China (no. 10932003), β€œ863” Project of China (no. 2009AA04Z101), and β€œ973” National Basic Research Project of China (no. 2010CB832700). These supports are gratefully acknowledged. Many thanks are due to the reviewers for their valuable comments.