陀螺系統(tǒng)時間有限元方法

隋永楓，高 強(qiáng)，鐘萬勰

(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116023;2.杭州汽輪動力集團(tuán)中央研究院 博士后科研工作站，杭州 310022;3.西安交通大學(xué) 動力工程及工程熱物理博士后流動站，西安 710049)

動力學(xué)方程的時程積分是經(jīng)常面對的課題，差分法數(shù)值積分是最常采用的方法［1］。在物理與力學(xué)中大量的保守體系可用Hamilton體系描述。馮康［2］指出，保守體系的差分格式應(yīng)當(dāng)保辛。保辛即保持保守體系結(jié)構(gòu)特性。但常用的差分格式并未考慮保辛;應(yīng)當(dāng)指出，通常的差分格式不保辛并非錯誤，而是逼近真實解不夠好。保辛的優(yōu)點就是更好的逼近。基于變分原理的有限元是自動保辛的［3］，。根據(jù)結(jié)構(gòu)力學(xué)與最優(yōu)控制的模擬理論［4］，空間坐標(biāo)有限元方法應(yīng)用廣泛［5］，它也可以用于時間坐標(biāo)。事實上，Zienkiewicz［6］已提出時間坐標(biāo)有限元法，但并未要求保辛。動力學(xué)有作用量變分原理，故其有限元可基于此變分原理。從變分原理導(dǎo)出的時間有限元矩陣也有對稱性，即也是保辛的。

陀螺系統(tǒng)的時程積分問題可用很多方法求解，如:Runge-Kutta方法、Newmark方法及Wilson方法等，這些方法中有些只有在選取特定參數(shù)下才具有保辛特點，而以上方法均未提及保辛的內(nèi)容。本文將時間有限元方法應(yīng)用于陀螺系統(tǒng)的時程積分問題，保持了陀螺系統(tǒng)的結(jié)構(gòu)特性，更好地逼近了陀螺系統(tǒng)的性態(tài)，特別對系統(tǒng)的長期時程分析具有很好的效果。

1 陀螺系統(tǒng)

如文獻(xiàn)［4］中所述，陀螺系統(tǒng)拉格朗日函數(shù)為:

其中:q為n維廣義位移向量。由此導(dǎo)出動力方程:

其對應(yīng)的變分原理見(3)式，(4)式為對應(yīng)的初值條件。

其中n×n的對稱陣M，K分別為質(zhì)量陣與剛度陣，G為反對稱陀螺陣，此為保守系統(tǒng)方程。

2 陀螺系統(tǒng)有限元列式

圖1 時間單元Fig.1 Time element

與位移有限元分析方法推導(dǎo)原理類似，在一時間區(qū)段(ta，tb)內(nèi)(圖1)的作用量函數(shù)(單元變形能)為:

其中:

式中:Kk即為時間單元剛度陣對稱陣，f1k為時間單元節(jié)點上等效外力，N為形函數(shù)矩陣。

引入動量(對偶)向量:

由最小作用量原理，聯(lián)合式(7a)、式(7b)與式(8)可導(dǎo)出:

其中:

3 穩(wěn)定性分析［7－8］

積分算法的穩(wěn)定性只依賴于逼近算子的本征值，假定λi是n維傳遞矩陣S的譜半徑，定義:

則當(dāng)且僅當(dāng)ρ(S)≤1時算法才是穩(wěn)定的。

陀螺系統(tǒng)傳遞算子矩陣S的本征問題表達(dá)為:

根據(jù)辛矩陣特點及式(12)可得:

式(13)說明1/μ仍是ST的本征值，而S與ST應(yīng)有相同的本征值，故知μ與其倒數(shù)1/μ必同時是辛矩陣S的本征值。

綜上可知，本文算法的穩(wěn)定性準(zhǔn)則為傳遞辛矩陣S的譜半徑必全為1。

4 形函數(shù)矩陣

在時間單元 (ta，tb)區(qū)段的進(jìn)行插值，每個對應(yīng)的系統(tǒng)自由度是相互獨立的，本文各自由度的插值函數(shù)相同，即:

以兩個自由度陀螺系統(tǒng)為例，討論時域有限元形函數(shù)的選取。取一次線性時域響應(yīng)模式，即一維Lagrange多項式插值，有:

其中:

則:

將(14)式代入(6)式可得時間單元的剛度陣:

其中:

5 非齊次外力項處理

同上節(jié)，以兩點線性插值形函數(shù)為例，通過(6)式計算系統(tǒng)的節(jié)點外力列式。如果外力f1=f0與時間無關(guān)，系統(tǒng)的節(jié)點外力可解析求出，對單自由度系統(tǒng)，外力式可表示為:

如果外力f1=f0(t)與時間有關(guān):① 當(dāng)外力f0(t)為可積分函數(shù)(如正弦，余弦，多項式和冪指數(shù)等形式)時，則節(jié)點外力列式可解析求解。② 當(dāng)外力f0(t)為不可積分函數(shù)時，可用數(shù)值積分方法(如高斯積分，柯茨積分等)求解［9］。

6 例題

本節(jié)給出兩自由度無阻尼自由陀螺振動系統(tǒng)，該方程能方便得到解析解，便于比較。

初值為:

本文解析解為:

Δt=0.15 s時，可得系統(tǒng)時間單元剛度陣及相應(yīng)的傳遞辛矩陣，如下。可以證明該辛矩陣本征值的模均為1，滿足本文算法收斂性條件。

圖2～圖7分別為應(yīng)用Runge-Kutta法、Newmark法及本節(jié)方法(TDT)求得3000 s內(nèi)系統(tǒng)的時間歷程曲線及對應(yīng)的絕對誤差曲線圖。其中Newmark法的參數(shù)為 α =0.25，δ=0.5，圖2、圖4 及圖6 中兩橫線分別為精確解的包絡(luò)線。從圖2和圖3可以看到，Runge-Kutta法的解在積分時間過長時超出了邊界，振幅是發(fā)散的，其絕對誤差將發(fā)散到無窮大，算法不穩(wěn)定。圖4及圖6表明Newmark法與本文方法(TDT)與精確解包絡(luò)線吻合的很好，可知這兩種方法的振幅并未發(fā)散，保證了該保守系統(tǒng)總體的能量守恒，算法上無耗散［2］。而圖5和圖7的絕對誤差曲線圖表明，Newmark法與本文方法(TDT)的絕對誤差是由于相位的誤差導(dǎo)致的，振幅不發(fā)散，絕對誤差呈周期性擴(kuò)大。在計算長時間歷程曲線時本節(jié)方法的計算精度穩(wěn)定，優(yōu)于Newmark法。

7 結(jié)論

綜上所述，由于時域有限元方法是保辛的，所以無論是長期時程分析還是短期時程分析，算法都能很好地保持原系統(tǒng)的性態(tài)。而振幅保持不變則表明了本文算法沒有產(chǎn)生能量耗散，雖然Newmark法的振幅也無明顯變化，但其絕對誤差呈相對較快速地擴(kuò)大趨勢，而本文方法的絕對誤差呈周期性變化且在相當(dāng)長的時間步長內(nèi)無明顯擴(kuò)大。通過算例證明，在同等條件下，本文方法的計算精度高于Runge-Kutta法、Newmark法，具有很好的應(yīng)用前景。

［1］Press W H，Teukolsky S A，Vetterling W T，et al.Numerical recipes in c［M］.Cambridge Univ.Press，1992.

［2］馮 康，秦孟兆.Hamilton體系的辛計算格式［M］.杭州:浙江科技出版社，2004.

［3］鐘萬勰.分析結(jié)構(gòu)力學(xué)與有限元［J］.動力學(xué)與控制學(xué)報，2004，2(4):1 －8.

［4］鐘萬勰.應(yīng)用力學(xué)對偶體系［M］.北京:科學(xué)出版社，2002.

［5］Oden J T，Belytschko T，Babuska I，et al.Research directions in computational mechanics［J］.Comput. Methods Appl.Mech.Engrg，2003，192(7 －8):913 －922.

［6］Zienkiewicz O C，Taylor R.The finite element method.5-th ed［M］.NY:McGraw-Hill，2000.

［7］Bathe K J，Wilson E L.Numerical methods in finite element analysis［M］.Prentice-Hall，Inc.，1976.

［8］Richtmyer R D，Morton K W.Difference methods for intialvalue problems(2nd Edition)［M］.New York:John Wiley and Sons ，1967.

［9］隋永楓.轉(zhuǎn)子動力學(xué)的求解辛體系及其數(shù)值計算方法［D］.大連:大連理工大學(xué)，2006.