多體系統(tǒng)動力學(xué)Lie群微分-代數(shù)方程約束穩(wěn)定方法*

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

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

引言

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

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

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

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

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

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

(1)

(2)

(3)

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

(4)

(5)

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

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

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

(6)

(7)

(8)

(9)

此時方程(6)離散為:

(10)

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

3 數(shù)值算例

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

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

Φ1(q)=x1+R1X10

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

(11)

圖2 兩連桿末端運動軌跡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 雙連桿空間機械臂能量Fig.4 Energies of the two-link space manipulator

圖5 雙連桿空間機械臂約束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給出了仿真過程中的連桿末端軌跡,圖3為連桿質(zhì)心位置和Lagrange乘子λ的時間歷程圖.從圖中可以看到仿真過程穩(wěn)定.圖4為雙連桿空間機械臂系統(tǒng)動能、勢能和總能量變化圖,可以看出系統(tǒng)總能量守恒.圖5為系統(tǒng)所受約束方程(11)及相應(yīng)的速度級約束和加速度級約束時間歷程圖,可以看出各級約束均精確保持,且精度較高,從而驗證了約束穩(wěn)定方法的作用.圖6和圖7給出了約束穩(wěn)定方法得到的旋轉(zhuǎn)矩陣R的各分量變化圖及誤差eR=RTR-I,進一步驗證了約束穩(wěn)定方法對旋轉(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)定方法與指標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)定方法與指標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中可以看出,對指標1微分-代數(shù)方程使用隱式方法求解,可以精確保持加速度級約束,但是位移約束和速度級約束誤差較大；對指標2微分-代數(shù)方程使用隱式方法求解,可以精確保持速度級約束,而位移約束和加速度級約束產(chǎn)生違約；對指標3微分-代數(shù)方程求解,可以精確保持位移約束,而速度級約束和加速度級約束誤差較大.但是使用約束穩(wěn)定方法求解,可以同時精確保持位移約束、速度級約束和加速度級約束.另外,從系統(tǒng)總能量相對誤差比較可以看出,約束穩(wěn)定方法的能量誤差與指標2方法誤差接近,小于另外兩種方法.這四種方法得到的系統(tǒng)總能量相對誤差變化圖如圖8所示,可以看出約束穩(wěn)定方法在保持各級約束穩(wěn)定的同時,也減小了系統(tǒng)總能量誤差. 從表2可以看出,當(dāng)步長增大時,上述分析仍然成立,與表1比較可以得出,步長增大時,各結(jié)果誤差有不同程度的增加,其中加速度約束誤差增加明顯,說明該類方法中加速度約束對步長變化較為敏感.

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

4 結(jié)論

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