張冰冰 王剛 丁潔玉

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

引言

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

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

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

1 離散變分?jǐn)?shù)學(xué)模型

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

(1)

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

(2)

其中λ∈Rm為L(zhǎng)agrange乘子,Φ(q,t)∈Rm為完整約束函數(shù).

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

λiTΦ(qi,ti))

(3)

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

(4)

2 廣義坐標(biāo)和速度的Hermite插值

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

(5)

(6)

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

(7)

(8)

當(dāng)s=3時(shí),三點(diǎn)三次Hermite插值函數(shù)中,

(9)

其插值余項(xiàng)為:

(10)

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

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

(11)

(12)

其中,

k=1,…,s

(13)

(14)

(15)

當(dāng)s=2時(shí),兩點(diǎn)三次Hermite插值函數(shù)中,

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

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

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

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

(16)

其插值余項(xiàng)為:

(17)

3 數(shù)值算例

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

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

(18)

約束方程為:

(19)

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

圖2 兩連桿末端運(yùn)動(dòng)軌跡,t=20sFig.2 Trajectories of the ends of two links when t=20s

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

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

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

圖6 Runge-Kutta法兩連桿末端運(yùn)動(dòng)軌跡,t=20sFig.6 Trajectories of the ends of two links using Runge-Kutta method when t=20s

圖7 Runge-Kutta法兩連桿末端運(yùn)動(dòng)軌跡,t=100sFig.7 Trajectories of the ends of the two links using Runge-Kutta method when t=100s

圖8 三點(diǎn)三次Hermite插值兩連桿末端運(yùn)動(dòng)軌跡,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中可以看出,基于兩點(diǎn)三次和三點(diǎn)三次Hermite插值的仿真結(jié)果相差不大,三點(diǎn)三次Hermite插值的結(jié)果精度稍高一些,在加速度級(jí)約束誤差中表現(xiàn)比較明顯.在相同步長(zhǎng)情況下,Runge-Kutta方法系統(tǒng)總能量誤差較大,在較大步長(zhǎng)情況下尤為明顯,這是因?yàn)镽unge-Kutta方法對(duì)步長(zhǎng)有限制要求,適用于較小步長(zhǎng)情況.對(duì)不同步長(zhǎng)情況結(jié)果進(jìn)行比較,h=0.005時(shí),Runge-Kutta方法的能量誤差和運(yùn)行時(shí)間與h=0.01時(shí)基于Hermite插值的離散變分方法結(jié)果相當(dāng).但是,從時(shí)間歷程圖上可以進(jìn)一步看出,Runge-Kutta方法的能量誤差隨時(shí)間增加逐漸增大,而離散變分方法能量誤差能夠保持在固定范圍之內(nèi),不隨時(shí)間增加而增大,如圖9和圖10所示.

4 結(jié)論

本文基于Hermite插值,建立了含已知導(dǎo)數(shù)信息和含未知導(dǎo)數(shù)信息的多體系統(tǒng)動(dòng)力學(xué)離散變分方程,求解進(jìn)行多體系統(tǒng)動(dòng)力學(xué)仿真.數(shù)值算例表明,利用導(dǎo)數(shù)信息可以得到精確度較高的仿真結(jié)果.在步長(zhǎng)較大情況下,其能量誤差和約束誤差均優(yōu)于同步長(zhǎng)的Runge-Kutta方法.在長(zhǎng)時(shí)間仿真情況下,基于Hermite插值的離散變分方法可以保持系統(tǒng)總能量在一定范圍內(nèi)有界變化,不會(huì)產(chǎn)生誤差累積.由于該方法模型中只含有約束方程,因此結(jié)果對(duì)約束方程精確保持,在速度級(jí)約束和加速度級(jí)約束中產(chǎn)生違約誤差,可以通過(guò)約束投影方法進(jì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)時(shí)間離散變分積分法. 動(dòng)力學(xué)與控制學(xué)報(bào), 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高強(qiáng),鐘萬(wàn)勰. 非完整約束動(dòng)力系統(tǒng)的離散積分方法. 動(dòng)力學(xué)與控制學(xué)報(bào), 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白龍,戈新生. 基于球擺模型的離散變分積分子算法研究. 動(dòng)力學(xué)與控制學(xué)報(bào), 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 夏麗莉,國(guó)忠金,張偉. 基于離散Legenda變換的Hamilton系統(tǒng)的變分算法和辛結(jié)構(gòu). 華中師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 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)動(dòng)力學(xué)離散變分方法. 青島大學(xué)學(xué)報(bào)(自然科學(xué)版), 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))