滿淑敏， 高 強(qiáng)， 鐘萬(wàn)勰

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,工程力學(xué)系，大連 116023)

1 引 言

哈密頓系統(tǒng)的一個(gè)主要性質(zhì)是其相流保持辛結(jié)構(gòu)[1,2]，因此其近似積分的數(shù)值方法也應(yīng)該具有保辛的性質(zhì)[3]。在完整約束哈密頓動(dòng)力系統(tǒng)中，約束的存在給數(shù)值積分帶來(lái)了許多問(wèn)題。如一些傳統(tǒng)算法在處理完整約束系統(tǒng)時(shí)，約束不是直接施加，而是由等價(jià)的加速度約束所代替。在這個(gè)過(guò)程中，數(shù)值積分帶來(lái)的誤差會(huì)使得結(jié)果無(wú)法滿足系統(tǒng)原有的約束條件，導(dǎo)致硬曲面下沉這種不真實(shí)現(xiàn)象的產(chǎn)生。為解決這類(lèi)誤差，引入了復(fù)原力的概念，但這類(lèi)任意添加的參數(shù)必須要根據(jù)每一個(gè)計(jì)算過(guò)程進(jìn)行調(diào)整，并要引入人工的能量損失。對(duì)于高度受限的系統(tǒng)，這類(lèi)損失在最終所生成的動(dòng)力系統(tǒng)中將占據(jù)支配作用[4]。而保辛算法能保持系統(tǒng)的固有特性，其數(shù)值結(jié)果具有能量和動(dòng)量守恒等特點(diǎn)[5,6]。因此，研究完整約束系統(tǒng)的保辛算法具有重要意義。

變分積分法是構(gòu)造保辛算法的一類(lèi)重要手段[7,8]。在一類(lèi)變量表示的拉格朗日體系內(nèi)，利用哈密頓原理，文獻(xiàn)[7,9]通過(guò)對(duì)作用量進(jìn)行近似，給出了構(gòu)造保辛算法的通用方法。文獻(xiàn)[10]討論了積分點(diǎn)個(gè)數(shù)與廣義位移近似多項(xiàng)式階數(shù)不相同時(shí)算法的收斂階數(shù)。在完整約束系統(tǒng)中，為提高算法精度，Wenger等[11]同時(shí)使用Gauss積分和Lobatto積分，得到一種高精度的Galerkin變分積分法。

上述保辛算法都是基于拉格朗日函數(shù)構(gòu)造的，但實(shí)際有些問(wèn)題只能獲得哈密頓函數(shù)，而無(wú)法給出拉格朗日函數(shù)。此時(shí)，保辛算法的構(gòu)造基于對(duì)偶變量變分原理。在無(wú)約束哈密頓系統(tǒng)中，文獻(xiàn)[12-15]基于兩種不同的對(duì)偶變量變分原理，通過(guò)選取不同的獨(dú)立變量，分別構(gòu)造了四種不同類(lèi)型的高階保辛算法。這類(lèi)保辛算法只需哈密頓函數(shù)便可實(shí)施，拓寬了保辛算法的應(yīng)用領(lǐng)域。

本文從哈密頓函數(shù)表示的對(duì)偶變量變分原理出發(fā)，研究完整約束哈密頓動(dòng)力系統(tǒng)的保辛算法。首先，引入拉格朗日乘子，得到包含完整約束方程的對(duì)偶變量變分原理；然后，選取積分區(qū)間兩端位移為獨(dú)立變量，基于變分原理的一種離散化形式，得到能夠在位移插值點(diǎn)處高精度滿足完整約束的高階保辛算法。

2 完整約束哈密頓系統(tǒng)的基本方程和對(duì)偶變量變分原理

設(shè)哈密頓動(dòng)力系統(tǒng)位置坐標(biāo)由d維向量q= {q1,q2,…,qd}T表示，d維位置坐標(biāo)受到c

g(q) = {g1(q),g2(q),…,gc(q)}T=0

(1)

哈密頓原理表明，系統(tǒng)的運(yùn)動(dòng)方程可由變分方程(2)導(dǎo)出，

(2)

(3)

則由哈密頓函數(shù)表示的對(duì)偶變量變分原理為

(4)

完成方程(4)的變分得到

pT(t1)δq(t1)-pT(t0)δq(t0)

(5)

式中G(q) = ?g(q)/?q。根據(jù)式(5)，如果積分區(qū)間兩端位移給定，則變分原理導(dǎo)出系統(tǒng)的運(yùn)動(dòng)方程為

(6)

如果系統(tǒng)的運(yùn)動(dòng)方程提前滿足，則作用量?jī)H僅是時(shí)間區(qū)間左右兩端位置坐標(biāo)的函數(shù)，即

δS=pT(t1)δq(t1)-pT(t0)δq(t0)

(7)

方程(7)對(duì)應(yīng)了q(t0),p(t0)和q(t1),p(t1)之間的一個(gè)正則變換

p(t0) =- ?S[q(t0),q(t1)]/[?q(t0)]

p(t1) = ?S[q(t0),q(t1)]/[?q(t1)]

(8)

式(8)可用于構(gòu)造數(shù)值算法。本文將以方程(8)為出發(fā)點(diǎn)，構(gòu)造完整約束哈密頓動(dòng)力系統(tǒng)的保辛算法。

3 完整約束哈密頓系統(tǒng)的保辛算法

(9)

(k= 0,1,…,n)(10)

(k= 0,1,…,m)(11)

(k= 0,1,…,z)(12)

將式(9)代入式(4)，可得

(13)

(14)

根據(jù)式(7)，近似作用量?jī)H為積分區(qū)間兩端位移q0和qn的函數(shù)，對(duì)近似作用量變分可得

(i= 1,2,…,n-1)(15)

(i= 0,1,…,m)(16)

(i= 0,1,…,z)(17)

同時(shí)，根據(jù)正則變換(8)可得

?Sd/?q0+p0=0, ?Sd/?qn-pm=0

(18)

(19)

則式(14)的最后一項(xiàng)為

(20)

(i= 1,2,…,n-1)(21)

方程(21)為約束方程在位移插值點(diǎn)處的表達(dá)式。需要注意的是，由于初值q0已知，初始點(diǎn)處的約束方程g(q0) =0自動(dòng)滿足，因此，需要另外補(bǔ)充c個(gè)方程。本文使用積分區(qū)間右端點(diǎn)處的速度約束作為補(bǔ)充方程。將約束方程(1)對(duì)時(shí)間求導(dǎo)，并利用式(6)可得右端點(diǎn)處的速度約束為

(22)

至此，求解方程組(15,16,18,21,22)便可得到完整約束哈密頓動(dòng)力系統(tǒng)的解。本文采用牛頓迭代法求解。由于上述非線性方程組由對(duì)偶變量變分原理導(dǎo)出，因此算法為保辛算法。

4 數(shù)值算例

通過(guò)兩個(gè)算例分析保辛算法的數(shù)值性質(zhì)。數(shù)值算例在64位計(jì)算機(jī)和64位操作系統(tǒng)下完成，程序使用Fortran實(shí)現(xiàn)，采用四倍精度。

算例1考慮三維球面擺在非線性勢(shì)能場(chǎng)中的運(yùn)動(dòng)。三維球面擺擺長(zhǎng)l=1.0 m，質(zhì)量m=1.0 kg，位移與動(dòng)量分別為q= {xyz}T和p= {pxpypz}T，系統(tǒng)哈密頓函數(shù)和約束方程為

(23)

g(q) =x2+y2+z2-l2= 0

(24)

系統(tǒng)初始狀態(tài)為

(25)

大量數(shù)值結(jié)果表明，當(dāng)位移和動(dòng)量的插值階次分別為n和m時(shí)，算法所需最少Gauss積分點(diǎn)個(gè)數(shù)為g>max. (n,m)，而當(dāng)g≤max. (n,m)時(shí)，由于在牛頓迭代求解過(guò)程中雅可比矩陣奇異，因此不能求解。

(26)

(27)

由表1可知，算法收斂階數(shù)s滿足如下規(guī)律，

(28)

可見(jiàn)，位移和動(dòng)量插值階次的最佳搭配方案為n=m+1，此時(shí)算法的收斂階數(shù)和計(jì)算效率最高。表2 給出了m和g取不同值，且n=m+1時(shí)，算法的收斂階數(shù)。由表2可知，僅改變g的值不會(huì)改變算法的收斂階數(shù)。綜合表1和表2，位移的插值階次n、動(dòng)量的插值階次m以及Gauss積分點(diǎn)個(gè)數(shù)g的最佳搭配方案為m=n-1 =g-2。此時(shí)，算法的收斂階數(shù)為s= 2m+2。

表1 g= max. (n,m)+1，且n和m取不同值時(shí)保辛算法的收斂階數(shù)

表2 n= m+1，且n，m和g取不同值時(shí)保辛算法的收斂階數(shù)

算例2考慮一彈簧連桿裝置。兩個(gè)彈性系數(shù)為k，可以在光滑水平面上自由振動(dòng)的彈簧分別將一端固定，另一端連接半徑為r、質(zhì)量為M的光滑圓環(huán)形軌道，軌道上各放置一質(zhì)量為m的小球，兩球在軌道中運(yùn)動(dòng)，同時(shí)又由一長(zhǎng)為D的輕質(zhì)細(xì)桿相連。兩彈簧原長(zhǎng)均為l，振動(dòng)方向在同一直線上，固定端相距2L。建立圖3所示坐標(biāo)系，系統(tǒng)位移為q= {x1,y1,x2,y2,O1,O2}T， 動(dòng)量為p= {px1,py1,px2,py2,pO1,pO2}T。 其中(x1,y1)和(x2,y2)分別為左右兩側(cè)小球的位置坐標(biāo)，O1和O2分別為左右兩側(cè)軌道圓心的x軸坐標(biāo)。系統(tǒng)的哈密頓函數(shù)及約束方程為

圖1 三維球面擺數(shù)值積分結(jié)果

圖2 三維球面擺哈密頓函數(shù)相對(duì)誤差

k(L-l-r+O1)2/2+k(L-l-r-O2)2/2

(29)

(30)

系統(tǒng)參數(shù)取M=1.5 kg，m=1.0 kg，k=0.6 N/m，L=1.0 m，l=0.4 m，r=0.2 m，D=0.8 m，初值選取

q(0) = {-0.1,0,0.7,0,-0.3,0.5}T

p(0) = {0,0.05,0,-0.05,0,0}T

(31)

數(shù)值結(jié)果表明，算法參數(shù)的最佳搭配方案為m=n-1 =g-2，此時(shí)，算法階數(shù)為s= 2m+2。

圖3 彈簧連桿裝置

圖4 彈簧連桿裝置中兩球運(yùn)動(dòng)軌跡

圖5 彈簧連桿裝置完整約束之和

5 結(jié) 論

本文選取積分區(qū)間兩端位移為獨(dú)立變量，將位移、動(dòng)量及拉格朗日乘子分別采用拉格朗日多項(xiàng)式近似，對(duì)作用量中包含約束的積分項(xiàng)使用Lobatto積分近似，對(duì)作用量中不包含約束的積分項(xiàng)使用Gauss積分近似，從而得到近似作用量。然后，在此近似作用量的基礎(chǔ)上，基于對(duì)偶變量變分原理，構(gòu)造了求解約束哈密頓系統(tǒng)的高階保辛算法。通過(guò)數(shù)值算例得到如下結(jié)論。

(1) 如果位移、動(dòng)量和拉格朗日乘子的插值階次分別為n，m和z，Gauss積分點(diǎn)個(gè)數(shù)為g，則本文算法所需最少Gauss積分點(diǎn)個(gè)數(shù)為g>max. (n,m)。

(2) 位移插值階次、動(dòng)量插值階次以及Gauss積分點(diǎn)個(gè)數(shù)的最佳搭配方案為m=n-1 =g-2。此時(shí)，算法的收斂階數(shù)為s= 2m+2。

(3) 算法在位移插值點(diǎn)處高精度地滿足完整約束方程。