基于Hermite插值的多體系統(tǒng)動力學離散變分方法*

張冰冰 王剛 丁潔玉

(1.青島大學計算機科學技術(shù)學院, 青島 266071) (2.青島大學數(shù)學與統(tǒng)計學院, 青島 266071) (3.青島大學計算力學與工程仿真研究中心, 青島 266071)

引言

多體系統(tǒng)動力學仿真目前已成為工程仿真的重要內(nèi)容之一,在機械、航天、車輛、機器人學及生物力學等領(lǐng)域均有廣泛應(yīng)用.目前多體系統(tǒng)動力學仿真算法主要分為兩個方向,一是對動力學方程進行離散求解,如Runge-Kutta方法、Newmark方法、HHT方法、廣義-α方法和辛算法等,二是對Lagrange函數(shù)進行離散,利用離散力學變分原理,得到離散方程進行求解,稱為離散變分方法.

離散變分方法基本框架最早由Marsden等[1]建立,對于受完整約束的保守系統(tǒng),該方法能使辛-能量-動量均得到保持,其中能量是在一定誤差范圍內(nèi)有界變化.這類方法提出之后很快在非光滑系統(tǒng)、非完整約束系統(tǒng)、隨機系統(tǒng)等動力學領(lǐng)域得到應(yīng)用[2-4].Ober-Bl?baum等[5]基于Galerkin方法給出了無約束動力學系統(tǒng)高階離散變分方法.潘振寬和丁潔玉等[6,7]給出了完整約束和非完整約束多體系統(tǒng)動力學高階離散變分方法.高強等[8]采用了微分-代數(shù)方程求解的雷同策略,只在離散區(qū)段內(nèi)滿足無約束的積分,而在離散節(jié)點處則用約束沖量(Lagrange乘子向量)讓軌道發(fā)生轉(zhuǎn)折,給出了非完整約束動力系統(tǒng)的離散積分格式.白龍等[9]研究了基于球擺模型的離散變分積分子算法.夏麗莉等[10]給出了離散Hamilton系統(tǒng)的差分方程和辛格式.目前該類方法的研究主要為基于李群方法以及不連續(xù)Galerkin方法的離散變分方法[11,12].

離散變分方法的精度依賴于變量離散及相關(guān)的數(shù)值求積方法,變量離散主要采用Lagrange插值多項式,數(shù)值積分主要采用Gauss求積公式.拉格朗日插值在插值節(jié)點較多時存在Runge現(xiàn)象,影響插值精度,潘坤等[13]采用重心拉格朗日插值,討論了均勻節(jié)點與非均勻節(jié)點情況.本文考慮導數(shù)信息,對廣義坐標和廣義速度進行Hermite插值離散,討論基于Hermite插值的多體系統(tǒng)動力學離散變分方法.

1 離散變分數(shù)學模型

通常情況下,多體系統(tǒng)動力學Lagrange函數(shù)可以表示為:

(1)

對于完整約束多體系統(tǒng),其Hamilton作用量可以表達為:

(2)

其中λ∈Rm為Lagrange乘子,Φ(q,t)∈Rm為完整約束函數(shù).

在多體系統(tǒng)動力學數(shù)值仿真中,離散變分方法是先對Hamilton作用量進行數(shù)值離散,然后利用離散變分原理得到離散歐拉-拉格朗日方程.

λiTΦ(qi,ti))

(3)

其中Aj,j=0,1,…,M為數(shù)值積分系數(shù),tij為區(qū)間[ti,ti+1]的求積節(jié)點,i=0,1,…,N,j=0,1,…,M.

(4)

2 廣義坐標和速度的Hermite插值

2.1 含已知導數(shù)值的Hermite插值

(5)

(6)

其中αk(t),k=1,…,s和β(t)可以由s點牛頓插值多項式加上一個s次項,由已知導數(shù)值求解出待定系數(shù)得到.

(7)

(8)

當s=3時,三點三次Hermite插值函數(shù)中,

(9)

其插值余項為:

(10)

使用s點s次Hermite插值函數(shù),求解離散歐拉-拉格朗日方程(4),可以求得qi,2,…,qi,s,如果qi,s=q(ti+1),則直接可得q(ti+1),否則由式(7)可得q(ti+1)=qd(1).

2.2 含未知導數(shù)值的Hermite插值

(11)

(12)

其中,

k=1,…,s

(13)

(14)

(15)

當s=2時,兩點三次Hermite插值函數(shù)中,

α1(τ)=(1+2τ)(τ-1)2

α2(τ)=(3-2τ)τ2

β1(τ)=hτ(τ-1)2

β2(τ)=h(τ-1)τ2

(16)

其插值余項為:

(17)

3 數(shù)值算例

圖1 雙連桿平面機械臂Fig.1 Two-link planar manipulator

選取廣義坐標為q=[x1y1θ1x2y2θ2]T,則Lagrange函數(shù)為:

(18)

約束方程為:

(19)

設(shè)仿真時間為20s,將時間區(qū)間[0,20]按步長h=0.01等分,在每個小區(qū)間上對廣義坐標q(t)進行三點三次Hermite插值,仿真結(jié)果如圖2-5所示.

圖2 兩連桿末端運動軌跡,t=20sFig.2 Trajectories of the ends of two links when t=20s

圖2和圖3給出了使用三點三次Hermite插值得到的兩連桿末端運動軌跡以及位移時間歷程,可以看到連桿運動過程穩(wěn)定,沒有出現(xiàn)偏移誤差.而如果使用4階Runge-Kutta方法對系統(tǒng)動力學方程進行數(shù)值仿真,在時間步長h=0.01時的運動軌跡如圖6所示,可以看出桿1的軌跡已經(jīng)出現(xiàn)偏移.這種現(xiàn)象在長時間仿真情況中表現(xiàn)比較明顯,如圖7所示,當仿真時間為100s時,誤差積累已經(jīng)非常嚴重,以致于運動軌跡已經(jīng)嚴重失真.而使用基于三點三次Hermite插值的離散變分方法,在仿真時間100s時運動軌跡仍保持穩(wěn)定,如圖8所示.

圖3 兩連桿末端位移Fig.3 Displacements of the ends of two links

圖4 雙連桿平面機械臂能量Fig.4 Energies of the two-link planar manipulator

圖6 Runge-Kutta法兩連桿末端運動軌跡,t=20sFig.6 Trajectories of the ends of two links using Runge-Kutta method when t=20s

圖7 Runge-Kutta法兩連桿末端運動軌跡,t=100sFig.7 Trajectories of the ends of the two links using Runge-Kutta method when t=100s

圖8 三點三次Hermite插值兩連桿末端運動軌跡,t=100sFig.8 Trajectories of the ends of the two links using Hermite interpolation method when t=100s

表1 Hermite插值離散變分方法和Runge-Kutta方法 結(jié)果比較,h=0.01Table 1 Comparison of the results from Hermite interpolation discrete variational method and Runge-Kutta method when h=0.01

表2 Hermite 插值離散變分方法和Runge-Kutta 方法 結(jié)果比較,h=0.005Table 2 Comparison of the results from Hermite interpolationdiscrete variational method and Runge-Kutta methodwhen h=0.005

從表1和表2中可以看出,基于兩點三次和三點三次Hermite插值的仿真結(jié)果相差不大,三點三次Hermite插值的結(jié)果精度稍高一些,在加速度級約束誤差中表現(xiàn)比較明顯.在相同步長情況下,Runge-Kutta方法系統(tǒng)總能量誤差較大,在較大步長情況下尤為明顯,這是因為Runge-Kutta方法對步長有限制要求,適用于較小步長情況.對不同步長情況結(jié)果進行比較,h=0.005時,Runge-Kutta方法的能量誤差和運行時間與h=0.01時基于Hermite插值的離散變分方法結(jié)果相當.但是,從時間歷程圖上可以進一步看出,Runge-Kutta方法的能量誤差隨時間增加逐漸增大,而離散變分方法能量誤差能夠保持在固定范圍之內(nèi),不隨時間增加而增大,如圖9和圖10所示.

4 結(jié)論

本文基于Hermite插值,建立了含已知導數(shù)信息和含未知導數(shù)信息的多體系統(tǒng)動力學離散變分方程,求解進行多體系統(tǒng)動力學仿真.數(shù)值算例表明,利用導數(shù)信息可以得到精確度較高的仿真結(jié)果.在步長較大情況下,其能量誤差和約束誤差均優(yōu)于同步長的Runge-Kutta方法.在長時間仿真情況下,基于Hermite插值的離散變分方法可以保持系統(tǒng)總能量在一定范圍內(nèi)有界變化,不會產(chǎn)生誤差累積.由于該方法模型中只含有約束方程,因此結(jié)果對約束方程精確保持,在速度級約束和加速度級約束中產(chǎn)生違約誤差,可以通過約束投影方法進行修正.

圖9 不同方法系統(tǒng)總能量,h=0.01Fig.9 Total energies of the system using different methods when h=0.01

圖10 不同方法系統(tǒng)總能量,h=0.005Fig.10 Total energies of the system using different methods when h=0.005

1Wendlandt J M, Marsden J E. Mechanical integrators derived from a discrete variational principle.PhysicaD, 1997,106(3-4):223～246

2McLachlan R I, Perlmutter M. Integrators for non-holonomic mechanical systems.JournalofNonlinearSciences, 2006,16(4):283～328

3Leyendecker S, Marsden J E, Ortiz M. Variational integrators for constrained dynamical systems.JournalofAppliedMathematicsandMechanics, 2008,88(9):677～708

4Bou-Rabee N, Owhadi H. Stochastic variational integrators.IMAJournalofNumericalAnalysis, 2009,29(2):421～443

5Ober-Bl?baum S, Saake N. Construction and analysis of higher order Galerkin variational integrators.AdvancesinComputationalMathematics, 2015,41(6):955～986

6丁潔玉,潘振寬. 非完整約束多體系統(tǒng)時間離散變分積分法. 動力學與控制學報, 2011,9(4):289～292 (Ding J Y, Pan Z K. Time-discrete variational integrator for multibody dynamic systems with nonholonomic constraints.JournalofDynamicsandControl, 2011,9(4):289～292 (in Chinese))

7Ding J Y, Pan Z K. Higher order variational integrators of multibody system dynamics with constraints. Advances in Mechanical Engineering, 2014:383680-1-8

8高強,鐘萬勰. 非完整約束動力系統(tǒng)的離散積分方法. 動力學與控制學報, 2012,10(3):193～198 (Gao Q, Zhong W X. Numerical algorithms of dynamic system with nonholonomic constraints.JournalofDynamicsandControl, 2012,10(3):193～198 (in Chinese))

9白龍,戈新生. 基于球擺模型的離散變分積分子算法研究. 動力學與控制學報, 2013,11(4):295～300 (Bai L, Ge X S. The Discrete variational integrators method of the spherical pendulum.JournalofDynamicsandControl, 2013,11(4):295～300 (in Chinese))

10 夏麗莉,國忠金,張偉. 基于離散Legenda變換的Hamilton系統(tǒng)的變分算法和辛結(jié)構(gòu). 華中師范大學學報(自然科學版), 2017,51(4):449～454 (Xia L L, Guo Z J, Zhang W. Variational calculation and symplectic structure of Hamiltonian systems based the discrete Legenda transformation.JournalofCentralChinaNormalUnversity(NatureScienceEdition), 2017,51(4):449～454 (in Chinese))

11 Hall J. Convergence of Galerkin variational integrators for vector spaces and Lie groups[Ph.D Thesis]. San Diego: Department of Mathematics University of California, 2013

12 Muehlebach M, Heimsch T, Glocker C. Variational integrators—A continuous time approach.InternationalJournalforNumericalMethodsinEngineering, 2017,109:1549～1581

13 潘坤,丁潔玉,董賀威等. 基于重心插值的多體系統(tǒng)動力學離散變分方法. 青島大學學報(自然科學版), 2017,30(2):77～82 (Pan K, Ding J Y, Dong H W, et al. Discrete variational method of multi-body system dynamics based on center of gravity interpolation.JournalofQingdaoUniversity(NatureScienceEdition), 2017,30(2):77～82 (in Chinese))