基于球擺模型的離散變分積分子算法研究*

白龍 戈新生

(北京信息科技大學(xué)機電工程學(xué)院，北京 1 00192)

基于球擺模型的離散變分積分子算法研究*

白龍?戈新生

(北京信息科技大學(xué)機電工程學(xué)院，北京 1 00192)

在動力學(xué)系統(tǒng)長時間的仿真計算中，力學(xué)系統(tǒng)固有的結(jié)構(gòu)將影響到計算精度及穩(wěn)定性．離散變分積分子能夠保持力學(xué)系統(tǒng)的能量，動量及辛結(jié)構(gòu)的守恒．結(jié)合離散變分原理，通過對系統(tǒng)的拉格朗日函數(shù)進行離散化以及求變分和積分的過程，可以得到力學(xué)系統(tǒng)的離散變分積分子算法．該算法是一種遞歸算法，給定初始條件便可得到系統(tǒng)的動力學(xué)參數(shù)的時間歷程．使用該原理可以構(gòu)造具有完整約束的拉格朗日系統(tǒng)的辛－動量積分子方法．與連續(xù)算法相比，離散變分積分子算法能夠直接在離散拉格朗日函數(shù)的基礎(chǔ)上得到姿態(tài)與角速度的遞推公式，而不需要復(fù)雜的迭代計算．本文研究是基于第一類拉格朗日函數(shù)的離散變分積分子算法．球擺模型是一個具有完整約束的拉格朗日系統(tǒng)．仿真結(jié)果表明，系統(tǒng)的能量值在長時間的仿真中得到保持，且計算的精度與步長的數(shù)量級呈現(xiàn)二次方的關(guān)系，系統(tǒng)角速度和姿態(tài)的仿真結(jié)果都符合球擺的運動規(guī)律．

離散變分原理， 力學(xué)積分子， 拉格朗日函數(shù)， 能量守恒

引言

在動力學(xué)系統(tǒng)的實際計算中，能量守恒與動量守恒與算法能否保持系統(tǒng)的結(jié)構(gòu)具有很大關(guān)系．通過拉格朗日函數(shù)得到系統(tǒng)的運動微分方程的解法在長時間的仿真中不具有保持系統(tǒng)結(jié)構(gòu)的性質(zhì)，而力學(xué)積分子能夠保持力學(xué)系統(tǒng)結(jié)構(gòu)．因此從理論上來講，它能夠保持系統(tǒng)的能量，動量及辛結(jié)構(gòu)的守恒．通過離散變分原理，得到離散拉格朗日方程，進而得到力學(xué)系統(tǒng)的離散變分積分子．與傳統(tǒng)算法比較，該方法既能保證計算的效率，又具有較好的長時間仿真的性質(zhì)．

離散變分積分算法經(jīng)過近半個世紀的發(fā)展，已經(jīng)成為一個較完整的算法體系．其中Verlet算法［1］可以用來處理無約束的力學(xué)系統(tǒng)，Shake算法［2］是Verlet算法的延伸，用來處理帶有完整約束的力學(xué)系統(tǒng)．為處理完整約束系統(tǒng)因離散而帶來的速度約束，Rattle算法［3］作為Shake算法的延伸較好的處理了該問題．

近年來，基于李群的離散變分積分子算法得到廣泛研究［4］，并且在3D剛體擺的動力學(xué)計算中取得了良好的效果．該理論擺脫了傳統(tǒng)離散變分積分子算法的約束理論，利用李群理論保證了系統(tǒng)的結(jié)構(gòu)，使系統(tǒng)的能量和動量在長時間的仿真中都保持守恒．

本文研究的是帶有約束的拉格朗日系統(tǒng)的保能量離散變分積分子方法．該方法通過離散變分積分子方法，得到一個離散的拉格朗日函數(shù)，使得約束在每個時間步長都得到滿足．在長時間的仿真中，離散拉格朗日方程保證了辛結(jié)構(gòu)和動量的守恒，進而保證了能量的守恒．球擺是3D擺的一種簡化模型，在離散變分積分子算法下，其角速度、姿態(tài)是否符合其運動規(guī)律將為離散變分積分子算法應(yīng)用到3D擺的研究中提供理論保證．仿真結(jié)果表明，在離散變分積分子算法下，系統(tǒng)的能量維持恒定，并隨著時間步長的減小，計算的精度也逐步提高．

1 離散變分原理

給定一個結(jié)構(gòu)空間Q，離散拉格朗日函數(shù)是一個映射L:Q×→．

對于固定的正整數(shù)N，對L(qk+1，qk)求和運算是映射S:QN+1→ ，定義為

其中，qk∈Q，k∈ 是離散時間．對 L(qk+1，qk)求和是對L(q，q)積分的一個離散近似．通過對S在q1，…，qN－1求極限，得到離散形式的歐拉—拉格朗日方程

其中k∈{1，…，N－1}，Φ:Q×Q→Q×Q定義為 Φ(qk，qk－1)=(qk+1，qk)．如果D2L 是可逆的，那么方程(3)定義了離散映射Φ，使得系統(tǒng)在離散時間下向前演化．定義映射

通過使用T*Q上的規(guī)范二階形式，定義Q×Q上的二階量ω

選擇Q上的坐標(biāo)qi和T*Q上的正則坐標(biāo)(qi，pi)，該坐標(biāo)滿足 ΩCAN=dqi∧dpi以及 ΘCAN=pidqi，那么離散歐拉－拉格朗日方程是:

根據(jù)映射 Φ(qk+1，qk)=(qk+1，qk)得到

對方程(5)計算得到:

由于∧=0，則

由式(7)，對ω求Φ運算

由式(9)可知:Φ－1ω=ω，即 Φ 能夠保持系統(tǒng)的辛流形．

令離散拉格朗日函數(shù)在Q上的李群G的對角線作用下恒定，令ξ∈g，其中g(shù)是G的李代數(shù)．則

對式(11)求微分，令s=0得到:

其中，ξQ是無窮小生成元．考慮對L(qk+1，qk)的求和運算，(0 ＜k＜N)，通過qk+1(s)=exp(sξ)qk+1將qk+1變s∈ 的形式．根據(jù)qk+1(0)對S求極限得到:

由方程(7)，(12)還可以寫為以下形式:

從方程(14)減去方程(12)得到:

如果動量映射J:Q×Q→g*定義為:

則方程(16)說明動量映射由Φ:Q×Q→Q×Q保持．

2 離散拉格朗日方程

球擺模型如圖1所示，質(zhì)心在慣性坐標(biāo)系中的投影可以表示為

圖1 球擺模型Fig．1 The model of spherical pendulum

球擺的動能為

勢能可以寫為

由拉格朗日函數(shù)L=T－V可得

約束條件為

定義表達式:

將q和q的離散表達式分別代入連續(xù)拉格朗日函數(shù)中，得到離散拉格朗日函數(shù):

則連續(xù)拉格朗日函數(shù)(21)的離散形式為:

由離散勒讓德變換可知:

分別在qk，qk+1對約束條件g(q)=0進行離散化可以得到

使用拉格朗日乘子，將約束條件添加到離散拉格朗日函數(shù)后，求變分可以得到

在考慮因離散化而對拉格朗日系統(tǒng)產(chǎn)生不顯含速度約束的情況時，離散系統(tǒng)在不顯含速度的作用下，會使系統(tǒng)產(chǎn)生耗散現(xiàn)象，導(dǎo)致仿真結(jié)果不具有一般性，因此，要解決此問題，就要在離散動量遞推公式后再添加一個約束項，對系統(tǒng)的速度進行約束．

由 Rattle算法［3］的定義

將式(31)代入式(30)得到qk+1的遞推公式．再將得到的遞推公式代入約束條件g(qk+1=0)可求出式(31)中常量γk．

由式(26)可以得到系統(tǒng)的動量的遞推公式表示為

對約束條件g(qk+1)=0求導(dǎo)得到包含速度的約束條件

將式(30)、(32)代入式(33)便可求出未知量λk．

3 仿真實驗

球擺模型是3D擺的一種簡化模型，設(shè)擺的重量為 1 kg，g=9．8N/kg，擺長 l =1m，選取初始位移坐標(biāo)q=()以及初始速度q=(0，0，0)，時間步長為 0．001s．仿真結(jié)果如圖所示．

圖2和圖3表明，球擺系統(tǒng)隨時間呈現(xiàn)周期性運動，由于q和q表征的是系統(tǒng)在慣性坐標(biāo)系下的運動，根據(jù)式(17)和(18)得出θ和φ以及θ·和 φ隨時間的變化圖，如圖4、圖5所示．

圖2 球擺坐標(biāo)的變化規(guī)律Fig．2 The coordinate variation of the spherical pendulum

圖3 球擺速度的變化規(guī)律Fig．3 The velocity variation of the spherical pendulum

圖4 球擺角度的變化規(guī)律Fig．4 The angle variation of the spherical pendulum

圖5 球擺角速度的變化規(guī)律Fig．5 The angle velocity variation of the spherical pendulum

圖4表明，球擺的θ角度呈現(xiàn)周期性變化;φ角保持恒定，符合球擺的運動規(guī)律．

圖5表明，球擺的角速度也呈現(xiàn)周期性變化，且 φ在較高的精度下保持恒定．

圖6 球擺動能勢能及總能量圖Fig．6 The kinetic potential and total energy of spherical pendulum

由圖6可知，球擺系統(tǒng)的能量值呈現(xiàn)周期性變化，符合球擺運動時能量的變化規(guī)律，而且系統(tǒng)的總能量在宏觀上維持恒定，未出現(xiàn)隨時間增加而逐漸偏離理論值的現(xiàn)象．

圖7 球擺的總能量比較Fig．7 the comparison of the total energy of the spherical pendulum

由圖7可知，與四階龍格庫塔算法相比，離散變分積分子算法能夠保持系統(tǒng)的總能量維持在一恒定的范圍波動，不會出現(xiàn)隨仿真時間的增加而逐漸發(fā)散的現(xiàn)象．

表1為在不同數(shù)量級的步長下，系統(tǒng)的能量的波動范圍及誤差隨時間步長的變化．通過對比時間步長與誤差可以得出，誤差的數(shù)量級與時間步長的數(shù)量級呈現(xiàn)二次方的關(guān)系．即誤差隨步長的減小而逐漸減小．

表1 系統(tǒng)總能量變化Table 1 The fluctuate of the energy

4 總結(jié)與展望

離散變分積分子算法作為一種保結(jié)構(gòu)的算法，在力學(xué)系統(tǒng)長時間的仿真中有著明顯的優(yōu)勢．本文表明，在簡單非線性系統(tǒng)的計算中，傳統(tǒng)的離散變分積分子算法取得了長時間保能量的良好的效果，如何將其應(yīng)用到更加復(fù)雜的諸如3D剛體擺等非線性系統(tǒng)的計算中將是下一步研究的重點．

傳統(tǒng)離散變分積分子算法是基于類似于帶有約束條件的第一類拉格朗日方程的形式來達到保能量的要求，3D剛體擺是具有三個轉(zhuǎn)動自由度的完整系統(tǒng)，歐拉四元數(shù)的建模方法為將離散變分積分子算法應(yīng)用到3D剛體擺中提供了可能，該方法通過四個變量實現(xiàn)對3D剛體擺的描述，而四個變量滿足的約束條件可以充當(dāng)3D剛體擺在應(yīng)用離散變分積分子算法時所需的約束條件．

除傳統(tǒng)離散變分積分子算法以外，基于李群的離散變分積分子算法通過系統(tǒng)運動學(xué)方程中的反對稱結(jié)構(gòu)來保持系統(tǒng)的辛流形，從而達到保能量的要求．該算法因其不需要單獨尋找約束條件而更加簡潔有效．國內(nèi)外相關(guān)研究也證明，基于李群的離散變分積分子算法也可應(yīng)用到3D剛體擺的動力學(xué)計算中．

1 Anderson H C．Rattle:A velocity version of the shake algorithm for molecular dynamics calculations．Journal of Computational Physics，1983，52:24 ～34

2 Wendlandt J M，Marsden J E．Mechanical integrators derived from a discrete variational principle．Physica D:Nonlinear Phenomena，1997:223～246

3 Bacconi，B．Spacecraft attitude dynamics and control［PhD Thesis］．Michigan:Department of Aerospace Engineering University of Michigan．2006:25 ～33

4 Wendlandt J M．Control and simulation of multibody systems［PhD Thesis］．Berkeley:University of California．1997:18～40

5 Shen J L．Dynamics and control of a 3 D pendulum．Proceedings of the IEEE Conference on Decision and Control，2004 43rd IEEE Conference on Decision and Control(CDC)，2004:323 ～328

6 Chaturvedi N A．Nonlinear dynamics of the 3D pendulum．Journal of Nonlinear Science，2011:3 ～32

7 丁潔玉，潘振寬．非完整約束多體系統(tǒng)時間離散變分積分法．動力學(xué)與控制學(xué)報，2011，9(4):289～292(Ding J Y，Pan Z K．Time discrete variational integrator for multibody dynamic systems with nonholonomic constraints．Journal of Dynamics and Control，2011，9(4):289 ～ 2 92)(in Chinese)

8 Marsden J E，West M．Discrete mechanics and variational integrators．Acta Numerica，2001:357 ～514

9 Lee，Taeyoung．A lie group variational integrator for the attitude dynamics of a rigid body with applications to the 3D pendulum．Proceeding of the IEEE International Conference on Control Applications，2005:962 ～ 9 67

10 Lee T Y．Lie group variational integrators for the full body problems in orbital mechanics．Celestial Mech Dyn Astr，2007，98:121 ～144

11 丁希侖，劉穎．用李群李代數(shù)分析具有空間柔性變形桿件的機器人動力學(xué)．機械工程學(xué)報，2007，43(12):184～189(Ding X L，Liu Y．Dynamic analysis of robot with spatial compliant links using lie group and lie algebra．Chinese Journal of Mechanical Engineering，2007，43(12):184～189(in Chinese))

12 高強，鐘萬勰．非完整約束動力系統(tǒng)的離散積分方法．動力學(xué)與控制學(xué)報，2012，10(3):193～198(Gao Q，Zhong W X．Numerical algorithms for dynamic system with non －h(huán)olonomic constraints．Journal of Dynamics and Control，2012，10(3):193 ～198(in Chinese))

13 Mclachlan R，Perlmutter M．Integrators for nonholonomic mechanical systems．Journal of Nonlinear Science，2006，16:283～332

14 鐘萬勰，高強．約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分．動力學(xué)與控制學(xué)報，2006，4(3):193～200(Zhong W X，Gao Q．integration of constrained dynamical system via analytical structual mechanicals．Journal of Dynamics and Control，2006，4(3):193 ～200(in Chinese))

*The project supported by the National Natural Science Foundation of China(11072038)the B category of State key development program of the Natural Science Foundation of Beijing(KZ201110772039)

? Corresponding author E－mail:bailong0316jn@126．com

THE DISCRETE VARIATIONAL INTEGRATORS METHOD OF THE SPHERICAL PENDULUM*

Bai Long?Ge Xinsheng
Beijing Information Science＆Technology University，Beijing100192，China

The mechanical system＇s intrinsic structure may influence the long time computation＇s accurate and stability．The discrete variation integrators can conserve the energy momentum and the symplectic structure of the system．Combined with the discrete variation principle，the discrete variation integrator method can be obtained through the process of discretization varaiton and integration．This is a recursive algorithm that the time history of the parameters only need the initial condition．According this theory，a sympletic－momentum integrator can be formulated for the holonomic constraint Lagrange system．This method can get the recursive formula of the attitude and the angle velocity direct form the discrete Lagrange function and don＇t need complicated iterative computation．The discrete variation integrator method explored in this paper is based on the first Lagrange function．The spherical pendulum is a Lagrange system with holonomic constraints．The simulation result states that the energy be conserved in a long time simulation，and the accuracy of the computation presents a quadratic relation with the time step．The angle velocity and the attitude also present different character under two different algorithm．

discrete variation principle， mechanical integrator， Lagrange function， energy conservation

21 August 2012，

29 May 2013．

10．6052/1672－6553－2013－077

2012－08－21 收到第 1 稿，2013－05－29 收到修改稿．

*國家自然科學(xué)基金資助項目(11072038)、北京市自然科學(xué)基金重點項目B類(KZ201110772039)

E－mail:bailong0316jn@126．com