李亞男 李博文 丁潔玉,2? 潘振寬

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

引言

多體系統(tǒng)動(dòng)力學(xué)通常由微分-代數(shù)方程描述,其數(shù)值積分方法的研究是計(jì)算多體系統(tǒng)動(dòng)力學(xué)研究的重要內(nèi)容.該領(lǐng)域微分-代數(shù)方程求解的傳統(tǒng)方法由常微分方程數(shù)值求解拓展而來(lái),研究的重點(diǎn)是通過(guò)約束穩(wěn)定達(dá)到微分-代數(shù)方程求解的穩(wěn)定[1,2].

近年來(lái)逐漸發(fā)展起來(lái)的保結(jié)構(gòu)幾何數(shù)值積分方法[3]則通過(guò)保持連續(xù)系統(tǒng)盡量多的不變量以提高數(shù)值計(jì)算的穩(wěn)定性,如辛算法[4]、能量方法[5,6],及基于離散力學(xué)變分原理的能量、辛保持的變分?jǐn)?shù)值積分方法[7,8]、Lie群方法[9]等.其中,Lie群方法致力于保持物體姿態(tài)的特殊正交群特性,而其他方法主要針對(duì)物體姿態(tài)的參數(shù)化表達(dá),如用經(jīng)典的Euler角、Bryant角等,相應(yīng)的算法不能保持物體姿態(tài)固有的Lie群特性.

直接用特殊正交群表達(dá)物體姿態(tài)可有效避免參數(shù)化表達(dá)導(dǎo)致的奇異性,但由于群空間的非矢量特性,傳統(tǒng)的數(shù)值積分方法不能直接應(yīng)用.Simo等[10,11]較早提出在其矢量Lie代數(shù)上設(shè)計(jì)算法,從而保持姿態(tài)矩陣的正交性,以達(dá)到提高數(shù)值穩(wěn)定性的目的.相關(guān)基礎(chǔ)論研究是指針對(duì)矩陣微分方程的CG方法[12]和MKRK方法[13].前者基于第二類正則坐標(biāo),即離散Lie代數(shù)矢量,并已用于設(shè)計(jì)基于指標(biāo)3的微分-代數(shù)方程的廣義-α方法設(shè)計(jì)[14].后者基于第一類正則坐標(biāo),即離散計(jì)算偽轉(zhuǎn)動(dòng)角度增量,且被應(yīng)用于指標(biāo)1的微分-代數(shù)方程的Runge-Kutta方法設(shè)計(jì)[15].

本文主要針對(duì)Lie群表達(dá)的多剛體系統(tǒng)動(dòng)力學(xué)微分-代數(shù)方程指標(biāo)1、2和3形式,設(shè)計(jì)約束穩(wěn)定方法,使位移約束、速度級(jí)約束和加速度級(jí)約束同時(shí)精確保持,從而提高仿真精度.

1 Lie群表達(dá)的數(shù)學(xué)模型

多體系統(tǒng)動(dòng)力學(xué)方程可在Lie群空間表達(dá)為如下指標(biāo)3 微分-代數(shù)方程:

(1)

(2)

(3)

利用方程(2)和方程(3),Lie群表達(dá)的指標(biāo)2和指標(biāo)1微分-代數(shù)方程形式分別為:

(4)

(5)

方程(1),(4)和(5)均只包含了約束方程、速度級(jí)約束方程或者加速度級(jí)約束方程,使用不同的數(shù)值方法求解時(shí)只能精確保持一種約束方程.為了在計(jì)算過(guò)程中同時(shí)保持三種約束方程使其不發(fā)生違約,本文對(duì)上述微分-代數(shù)方程進(jìn)行改進(jìn).

2 Lie群微分-代數(shù)方程約束穩(wěn)定方法

引入Lagrange乘子參數(shù)μ,ω,構(gòu)造新的Lie群表達(dá)微分-代數(shù)方程如下:

(6)

(7)

(8)

(9)

此時(shí)方程(6)離散為:

(10)

由初始條件x1,v1,Ω1,R1使用牛頓迭代方法迭代求解可以得到xk,vk,Ωk,k=2,3,…,從而由式(9)可得Rk,k=2,3,….

3 數(shù)值算例

圖1為雙連桿空間機(jī)械臂,設(shè)桿1的質(zhì)量為m1=5kg,質(zhì)心坐標(biāo)為(x1,y1,z1),桿長(zhǎng)為l1=1m；桿2的質(zhì)量為m2=5kg,質(zhì)心坐標(biāo)為(x2,y2,z2),桿長(zhǎng)為l2=1m.

圖1 雙連桿空間機(jī)械臂Fig.1 Two-link space manipulator

Φ1(q)=x1+R1X10

Φ2(q)=x1+R1X11-x2-R2X21

(11)

圖2 兩連桿末端運(yùn)動(dòng)軌跡Fig.2 Trajectories of the ends of the two links

圖3 兩連桿質(zhì)心位置和Lagrange乘子Fig.3 Displacements of the mass center of the two links and the Lagrange multipliers

圖4 雙連桿空間機(jī)械臂能量Fig.4 Energies of the two-link space manipulator

圖5 雙連桿空間機(jī)械臂約束Fig.5 Constraints of the two-link space manipulator

圖6 約束穩(wěn)定方法得到的旋轉(zhuǎn)矩陣各分量Fig.6 Rotation matrix components obtained by using constraints stabilization method

圖2給出了仿真過(guò)程中的連桿末端軌跡,圖3為連桿質(zhì)心位置和Lagrange乘子λ的時(shí)間歷程圖.從圖中可以看到仿真過(guò)程穩(wěn)定.圖4為雙連桿空間機(jī)械臂系統(tǒng)動(dòng)能、勢(shì)能和總能量變化圖,可以看出系統(tǒng)總能量守恒.圖5為系統(tǒng)所受約束方程(11)及相應(yīng)的速度級(jí)約束和加速度級(jí)約束時(shí)間歷程圖,可以看出各級(jí)約束均精確保持,且精度較高,從而驗(yàn)證了約束穩(wěn)定方法的作用.圖6和圖7給出了約束穩(wěn)定方法得到的旋轉(zhuǎn)矩陣R的各分量變化圖及誤差eR=RTR-I,進(jìn)一步驗(yàn)證了約束穩(wěn)定方法對(duì)旋轉(zhuǎn)矩陣R的正交性的保持,即Lie群結(jié)構(gòu)的保持.

圖7 約束穩(wěn)定方法得到的旋轉(zhuǎn)矩陣正交性誤差Fig.7 Errors of the orthogonality of the rotation matrix obtained by using constraints stabilization method

表1 約束穩(wěn)定方法與指標(biāo)1、2、3DAEs結(jié)果比較(h=0.00001)Table 1 Comparison of the results obtained by constraints stabilization method and index 1,2,3 DAEs (h=0.00001)

表2 約束穩(wěn)定方法與指標(biāo)1、2、3DAEs結(jié)果比較(h=0.0001)Table 2 Comparison of the results obtained by constraints stabilization method and index 1,2,3 DAEs (h=0.0001)

從表1中可以看出,對(duì)指標(biāo)1微分-代數(shù)方程使用隱式方法求解,可以精確保持加速度級(jí)約束,但是位移約束和速度級(jí)約束誤差較大；對(duì)指標(biāo)2微分-代數(shù)方程使用隱式方法求解,可以精確保持速度級(jí)約束,而位移約束和加速度級(jí)約束產(chǎn)生違約；對(duì)指標(biāo)3微分-代數(shù)方程求解,可以精確保持位移約束,而速度級(jí)約束和加速度級(jí)約束誤差較大.但是使用約束穩(wěn)定方法求解,可以同時(shí)精確保持位移約束、速度級(jí)約束和加速度級(jí)約束.另外,從系統(tǒng)總能量相對(duì)誤差比較可以看出,約束穩(wěn)定方法的能量誤差與指標(biāo)2方法誤差接近,小于另外兩種方法.這四種方法得到的系統(tǒng)總能量相對(duì)誤差變化圖如圖8所示,可以看出約束穩(wěn)定方法在保持各級(jí)約束穩(wěn)定的同時(shí),也減小了系統(tǒng)總能量誤差. 從表2可以看出,當(dāng)步長(zhǎng)增大時(shí),上述分析仍然成立,與表1比較可以得出,步長(zhǎng)增大時(shí),各結(jié)果誤差有不同程度的增加,其中加速度約束誤差增加明顯,說(shuō)明該類方法中加速度約束對(duì)步長(zhǎng)變化較為敏感.

圖8 使用不同方法得到的系統(tǒng)總能量Fig.8 Total energies of the system obtained by using different methods

4 結(jié)論

本文針對(duì)Lie群表達(dá)的多體系統(tǒng)動(dòng)力學(xué)微分-代數(shù)方程,設(shè)計(jì)了約束穩(wěn)定的求解方法,使仿真過(guò)程中位移約束、速度級(jí)約束和加速度級(jí)約束能夠同時(shí)得到保持.數(shù)值算例表明,相比較使用同樣向后差商隱式方法求解的指標(biāo)1、指標(biāo)2和指標(biāo)3微分-代數(shù)方程,約束穩(wěn)定方法得到的結(jié)果更為精確,在保持各級(jí)約束方程的同時(shí),也保持了旋轉(zhuǎn)矩陣的正交性,是一種保Lie群結(jié)構(gòu)的數(shù)值方法.在該方法的基礎(chǔ)上可以很方便地設(shè)計(jì)高階數(shù)值方法,進(jìn)一步提高結(jié)果精度.

1Simenon B. Computational flexible multibody dynamics—A differential algebraic approach. Berlin:Springer, 2013

2Fon-Llagunes J M. Multibody dynamics-computational methods and applications. Switzerland:Springer, 2016

3Betsch P. Structure-preserving integrators in nonlinear structural dynamics and flexible multibody dynamics. Switzerland:Springer, 2016

4Jay L O. Symplectic partitioned Runge-Kutta methods for constrained Hamiltonian systems.SIAMJournalonNumericalAnalysis, 1996,33(1): 368～387

5Lens E, Cardona A. An energy preserving/decaying scheme for nonlinearly constrained multibody systems.MultibodySystemDynamics, 2007,18(3):435～470

6Betsch P, Uhlar S. Energy-momentum conserving integration of multibody dynamics.MultibodySystemDynamics, 2007,17(4):243～289

7Marsden J E, West M. Discrete mechanics and variational integrators.ActaNumerica, 2003,10(1):357～514

8Ding J Y, Pan Z K. Higher order variational integrators for multibody system dynamics with constraints.AdvancesinMechanicalEngineering, 2014,6(1):383680-1-8

9Iserles A, Munthe-Kaas H Z, N?rsett S, et al. Lie-group methods.ActaNumerica, 2016,9(2):215～365

10 Simo J C, Vu-Quoc L. On the dynamics in space of rods undergoing large motions—A geometrically exact approach.ComputerMethodsinAppliedMechanicsandEngineering, 1988,66(2):125～161

11 Lewis D, Simo J C. Conserving algorithms for the dynamics of Hamiltonian systems on Lie groups.JournalofNonlinearScience, 1994,4(1):253～299

12 Crouch P E, Grossman R. Numerical integration of ordinary differential equations on manifolds.JournalofNonlinearScience, 1993,3(1):1～33

13 Munthe-Kaas H. High order Runge-Kutta methods on manifolds.AppliedNumericalMathematics, 1999,29(1):115～127

14 Arnold M, Brüls O, Cardona A. Error analysis of generalized-α Lie group time integration methods for constrained mechanical systems.NumerischeMathematik, 2015,129(1):149～179

15 Terze Z, Müller A, Zlatar D. Lie-group integration method for constrained multibody systems in state space.MultibodySystemDynamics, 2015,34(3): 275～305