一種基于Hamel 形式的無條件穩(wěn)定動(dòng)力學(xué)積分算法1)

顧崴 劉 鋮?, 安志朋 史東華,???

* (北京理工大學(xué)信息與電子學(xué)院,北京 100081)

? (中國(guó)工程物理研究院高性能數(shù)值模擬軟件中心,北京 100088)

** (北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100088)

?? (北京理工大學(xué)宇航學(xué)院,北京 100081)

*** (北京理工大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,北京 100081)

??? (復(fù)雜信息數(shù)學(xué)表征分析與應(yīng)用北京市重點(diǎn)實(shí)驗(yàn)室,北京 100081)

引言

時(shí)間積分算法被廣泛應(yīng)用于動(dòng)力學(xué)系統(tǒng)的求解中.動(dòng)力學(xué)方程的時(shí)間積分經(jīng)常會(huì)出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象,有限元空間離散也通常會(huì)造成偽高頻振蕩,因此,針對(duì)上述問題的數(shù)值積分算法的研究一直備受關(guān)注.本文將基于活動(dòng)標(biāo)架法和離散力學(xué),發(fā)展對(duì)動(dòng)力學(xué)系統(tǒng)無條件穩(wěn)定的數(shù)值積分格式.

Newmark 方法[1]是廣泛使用的時(shí)間積分方法,它包括許多格式,比如中心差分法(CDM)和梯形規(guī)則(TR)等,其中隱式TR 具有二階精度和無條件穩(wěn)定性.此外,對(duì)于線性保守系統(tǒng),TR 可以長(zhǎng)時(shí)間保持能量.但對(duì)于波傳播和剛性系統(tǒng),TR 不能過濾掉不需要的高頻.為解決此問題,研究人員開發(fā)了更理想的方法,如廣義 α 方法[2]、WBZ-α 方法[3]和HHT-α方法[4].

廣義 α 方法[2]是求解結(jié)構(gòu)動(dòng)力學(xué)問題的一種隱式方法,可以通過單個(gè)參數(shù)來控制數(shù)值耗散,使得在耗散系統(tǒng)高頻響應(yīng)的同時(shí)較好地保持系統(tǒng)的低頻響應(yīng).廣義 α 方法在時(shí)間上具有二階精度且無條件穩(wěn)定,已成功用于廣泛的工程應(yīng)用[5-8].為應(yīng)用于多體動(dòng)力學(xué),Yen 等[9]將其擴(kuò)展到約束系統(tǒng).文獻(xiàn)[10]分析了廣義 α 格式和Lagrange 乘子的二階收斂性.張雄和王天書[11]和田強(qiáng)[12]比較了上述算法的精度和耗散特性,研究了力學(xué)系統(tǒng)守恒量的變化特性.此外,文獻(xiàn)[13]將二階廣義 α 方法推廣到三階,并且可以通過單個(gè)參數(shù)控制數(shù)值耗散,同時(shí)也證明了該方法是無條件穩(wěn)定的.季奕等[14-15]提出了非線性動(dòng)力問題的單步以及兩步可控耗散時(shí)間積分法,這些方法是二階精度和無條件穩(wěn)定的,并且不需要初始加速度矢量.

為了研究具有李群結(jié)構(gòu)的柔性多體系統(tǒng),Brüls等[16]提出了李群廣義 α 時(shí)間積分算法,并證明了該方法具有全局二階精度.Arnold 等[17]詳細(xì)分析了李群廣義 α 方法的收斂性.文獻(xiàn)[18]給出了SE(3)群的廣義 α 方法.文獻(xiàn)[19] 基于BDF(backward differentiation formulae)方法提出了更高精度的李群積分BLieDF 方法.

在幾何力學(xué)領(lǐng)域,變分積分子[20-21]也是一種數(shù)值積分算法.它源于離散力學(xué)中離散Hamilton 原理,能如實(shí)反映動(dòng)力學(xué)系統(tǒng)性質(zhì),目前已被廣泛應(yīng)用于Lagrange 系統(tǒng)[22]和約束系統(tǒng)[23]等各種問題,并取得良好的效果.變分積分子的核心思想是先對(duì)系統(tǒng)相應(yīng)的變分原理進(jìn)行離散,然后通過離散變分運(yùn)算,得到連續(xù)Euler-Lagrange 方程的離散形式.Lew 等[21]指出,在有限維情形下,由于直接從離散的變分原理出發(fā),變分積分子能夠如實(shí)反映連續(xù)系統(tǒng)的動(dòng)力學(xué)性質(zhì),如保守系統(tǒng)的保動(dòng)量和保辛結(jié)構(gòu),保守和非保守系統(tǒng)的離散能量守恒,因而避免了傳統(tǒng)數(shù)值方法存在的數(shù)值耗散問題,在數(shù)值上體現(xiàn)出更好的性質(zhì).隨著離散力學(xué)的發(fā)展,研究人員提出了辛能量–動(dòng)量變分積分子[24]和李群變分積分子[25]等.文獻(xiàn)[26]結(jié)合變分積分子,提出了二階數(shù)值時(shí)間積分方法RATTLie,并將它們應(yīng)用于構(gòu)型空間為李群的約束力學(xué)系統(tǒng).同時(shí),研究人員又發(fā)展了高階變分積分子,如Galerkin 變分積分子[27]和譜變分積分子[28]等,這些高階變分積分子既能保辛和保動(dòng)量,又能實(shí)現(xiàn)任意階精度或具有指數(shù)收斂性.

Bloch 等[29]將活動(dòng)標(biāo)架法融入變分原理,推導(dǎo)了有限維力學(xué)系統(tǒng)的Hamel 動(dòng)力學(xué)方程,并提出了Hamel 變分積分子[30],其框架可統(tǒng)一描述李群和李代數(shù)變分積分子.Hamel 變分積分子產(chǎn)生了一種更簡(jiǎn)單的算法,目前已應(yīng)用于機(jī)械臂最優(yōu)控制問題[31],展現(xiàn)了良好性能.Shi 等[32]提出了無窮維力學(xué)系統(tǒng)的Hamel 形式,進(jìn)而研究了經(jīng)典場(chǎng)論的Hamel 形式[33-34],并在此基礎(chǔ)上構(gòu)建了Hamel場(chǎng)變分積分子.王亮等[35]利用Hamel 場(chǎng)變分積分子,給出了幾何精確梁的動(dòng)力學(xué)建模和數(shù)值仿真,該算法能長(zhǎng)時(shí)間保持系統(tǒng)的能量、動(dòng)量和幾何結(jié)構(gòu),高山等[36]對(duì)其保持動(dòng)量給出了理論證明.因此,Hamel場(chǎng)變分積分子可以對(duì)一類無窮維力學(xué)系統(tǒng)進(jìn)行準(zhǔn)確刻畫,是力學(xué)系統(tǒng)長(zhǎng)時(shí)間數(shù)值模擬的理想選擇.

Newmark 方法和廣義 α 方法等一系列時(shí)間積分算法最初均是通過經(jīng)驗(yàn)推導(dǎo)或總結(jié)出來的,后來賦予了各種理論基礎(chǔ),如泰勒級(jí)數(shù)和Galerkin 弱形式等,最后利用數(shù)值算法的分析手段確定其收斂性和穩(wěn)定性.這些方法都不是為了保持能量和動(dòng)量而設(shè)計(jì)的,在一些特殊參數(shù)情況下,它們可以保持動(dòng)量[37].文獻(xiàn)[38]從Hamilton 方法出發(fā),指出Newmark 族以及相關(guān)的積分算法,在離散力學(xué)的Veselov 形式下才是變分的[38].迄今為止,在變分積分子的框架下一直未出現(xiàn)一種能耗散系統(tǒng)高頻響應(yīng)同時(shí)又能較好保持系統(tǒng)低頻響應(yīng)的數(shù)值積分算法.因此,在該領(lǐng)域中,構(gòu)造一種無條件穩(wěn)定的時(shí)間積分算法是必要的,這樣既可以拓展變分積分子的理論體系,同時(shí)又可以豐富時(shí)間積分算法的理論基礎(chǔ).在這種背景下,本文基于離散力學(xué)框架和活動(dòng)標(biāo)架法,采用特殊的變分形式[39-40],提出了一種Hamel 廣義 α 方法,并推廣到了一般李群空間,理論分析和數(shù)值測(cè)試均顯示出了本文方法在精度、耗散和穩(wěn)定性方面的優(yōu)勢(shì).從而為進(jìn)一步研究無條件穩(wěn)定的高階數(shù)值積分算法提供參考依據(jù).

1 Hamel 形式與Hamel 變分積分子

1.1 Hamel 場(chǎng)形式

令Q為向量空間W上的無窮維光滑流形.Lagrange函數(shù)為L(zhǎng):TQ→R,其中TQ是位形空間Q的切叢.系統(tǒng)的動(dòng)力學(xué)由Euler-Lagrange 方程組決定.若選擇合適的標(biāo)架,可獲得系統(tǒng)動(dòng)力學(xué)的更簡(jiǎn)潔形式——Hamel 方程[34].令U是q∈Q的鄰域,定義一族可逆的有界線性算子 Ψq:W→TqQ為活動(dòng)標(biāo)架算子,?q∈U,同時(shí)也是可逆的有界線性算子.因而,對(duì)于 ?ξ ∈W定義向量場(chǎng)

給定兩個(gè)向量 ξ,η ∈W,定義一個(gè)反對(duì)稱括號(hào)算子[·,·]q:W×W→W

其中 [·,·] 是流形上兩個(gè)向量場(chǎng)的Jacobi-Lie 括號(hào).對(duì)于任意 ξ,η ∈W和 α ∈W*,括號(hào)算子 [·,·]q的對(duì)偶算子通過下式給出

給定非保守力F∈Γ(T*Q),其中Γ (T*Q) 表示T*Q的所有截面,則虛功記為

定理1 Hamilton-Pontryagin 原理的Hamel 形式[34]

令力學(xué)系統(tǒng)位形空間為Q,TQ是切叢,T*Q為余切叢.令q∈Q,(q,v)∈TQ且 (q,p)∈T*Q.令L:TQ→R為L(zhǎng)agrange 函數(shù),l是關(guān)于活動(dòng)標(biāo)架下局部坐標(biāo)(q,ξ)的Lagrange 函數(shù),則下列陳述等價(jià).

(1) 令 (q,v,p) 為Pontryagin 叢TQ⊕T*Q的局部坐標(biāo),則曲線 (q(t),v(t),p(t)),0 ≤t≤T為其相應(yīng)的作用積分

關(guān)于獨(dú)立變分 δq,δv,δp的臨界點(diǎn),即

且δq(0)=δq(T)=0.

(2) 曲線 (q(t),v(t),p(t)) 滿足隱式Euler-Lagrange 方程組

(3) 曲線 (q,ξ,μ)∈U×(W⊕W*) 是作用泛函

關(guān)于獨(dú)立變分 δq=Ψqη ,δ ξ 和 δμ 的臨界點(diǎn)

(4) 曲線 (q,ξ,μ) 滿足隱式Hamel 方程組的弱形式

對(duì)于任意點(diǎn)q∈Q,曲線 (q,ξ,μ) 滿足隱式Hamel 方程組的強(qiáng)形式

1.2 Hamel 變分積分子

假設(shè)Q是向量空間,連續(xù)曲線q:[0,T]→Q用離散序列近似表示.定義

式中參數(shù) α ∈[0,1],Δt是時(shí)間步長(zhǎng),令

為關(guān)于活動(dòng)標(biāo)架 Ψqk+α的速度分量.

通過連續(xù)Lagrange 函數(shù),定義其離散形式

相應(yīng)離散作用和為

同理,變量 η 的離散形式 ηk+α可以線性展開

通過離散變分原理,可得隱式Hamel 方程組強(qiáng)形式(12)的離散形式—–Hamel 變分積分子

其中Dil為L(zhǎng)agrange 函數(shù)對(duì)第i個(gè)參數(shù)的偏導(dǎo)數(shù).

2 向量空間上的Hamel 廣義 α 方法

本節(jié)考慮向量空間中力學(xué)系統(tǒng),基于活動(dòng)標(biāo)架法,構(gòu)造力學(xué)系統(tǒng)的Hamel 廣義 α 方法.

2.1 Hamel 方程

考慮一個(gè)受控系統(tǒng),其位形空間為向量空間Q,狀態(tài)空間為TQ,T*Q表示相空間.對(duì) ?q∈Q取標(biāo)架算子 Ψq為恒算子,定義系統(tǒng)的Lagrange 函數(shù)L:TQ→R

其中fτ,fext分別為反饋控制力和外力.

由前文定理1,可知下面的表述是等價(jià)的.

(1) 定義在TQ⊕T*Q上的曲線(q(t),v(t),p(t)) 是作用泛函

關(guān)于獨(dú)立變分 δq,δv,δp的臨界點(diǎn),即

(2) 曲線 (q(t),v(t),p(t)) 滿足隱式Hamel 方程組的弱形式

曲線 (q(t),v(t),p(t)) 滿足隱式Hamel 方程組的強(qiáng)形式

此時(shí)Hamel 方程組為隱式Euler-Lagrange 方程組.所以由式(24)可定義加速度函數(shù)

2.2 Hamel-廣義 α 方法

給定時(shí)間步長(zhǎng) Δt,區(qū)間[ 0,T] 等分為N個(gè)子區(qū)間,為簡(jiǎn)單起見,令區(qū)間 [tn,tn+1] 為 [0,Δt],為構(gòu)建數(shù)值積分算法,定義離散形式

并且滿足Mv0=p0,最后根據(jù)加速度函數(shù)定義相應(yīng)的離散加速度

此外,令外力和反饋控制律離散形式為

因此,可知力f所對(duì)應(yīng)的離散形式為

通過控制律fτ可獲得無條件穩(wěn)定的數(shù)值積分算法.定義輔助未知量 δq,δp和 δv的離散形式

其中,A1,A2,A3,A4為四個(gè)參數(shù).將離散形式(26)～式(32)帶入隱式Hamel 方程組的弱形式(23),求解關(guān)于的駐點(diǎn)方程,可得到比廣義 α 方法更一般的數(shù)值積分算法

式中參數(shù)分別為

其次,tn+1時(shí)刻離散加速度函數(shù)為

為了得到更穩(wěn)定的算法,令

綜上,方程(33)～(35)稱之為Hamel 廣義 α 方法.

注1.Hamel 廣義 α 方法可以退化為一些經(jīng)典算法,如Newmark、HHT-α 和廣義 α 方法等.如表1 所示,令A(yù)3=2(3β-1),A4=-3(4β-1),表格中參數(shù)α,αm,αf的含義參見文獻(xiàn)[1-2,4].

表1 Hamel 廣義 α 方法與經(jīng)典算法比較Table 1 Comparison between Hamel generalized-α method and classical algorithm

注2.對(duì)離散形式(26)進(jìn)行改進(jìn)如下

可得到更高階數(shù)值積分算法(33),其中參數(shù)μ0,μ1,γ0,γ1變?yōu)?/p>

其余參數(shù)不變,其中I為單位陣,K=?2V(?qi?qj) 為彈性矩陣.選取特定的參數(shù),該數(shù)值格式會(huì)得到更高的收斂速率,具體數(shù)值算例及其收斂速率見3.3 節(jié).

3 算法性質(zhì)

3.1 算法分析

考慮一個(gè)力學(xué)系統(tǒng)

相應(yīng)的Euler-Lagrange 方程組為

Hamel 廣義 α 方法應(yīng)用于上述方程,令f=0,則離散方程可寫成矩陣形式

故本文所提方法為二階精度算法.

算法的穩(wěn)定性和數(shù)值耗散取決于其輔助矩陣的特征值.譜半徑是數(shù)值耗散的一種度量,譜半徑越小,數(shù)值耗散越大.算法的譜半徑由矩陣的特征值決定,定義為

其中 λi為矩陣的第i個(gè)特征值.若 ρ ≤1,則對(duì)于線性問題算法是無條件穩(wěn)定的.

理想算法應(yīng)該是嚴(yán)格控制高頻,但是低頻區(qū)域不引起過度耗散,在低頻域的譜半徑值接近于1.隨著 Ω 的增加,譜半徑值平滑地減小.若輔助矩陣的主根形式為

其中a和b是實(shí)數(shù),只要

成立,則高頻耗散最大化.由(41)可知

滿足這個(gè)條件,Hamel 廣義 α 方法便可實(shí)現(xiàn)高頻耗散最大化.此外,當(dāng)A=0.5 時(shí),可得A4=0,無法確定A3的值,故而設(shè)定A3=-0.5.

令ρ∞表示高頻極限處的譜半徑,變量可表示為

為方便描述方程(44),借助 ρ∞,定義B1和B2

綜上所述,同時(shí)滿足方程式(40)、式(42)和式(45),則Hamel 廣義 α 方法是無條件穩(wěn)定、二階精度的算法,并且可實(shí)現(xiàn)高頻耗散最大,低頻耗散最小.

3.2 算法穩(wěn)定性

譜分析是分析時(shí)間積分算法穩(wěn)定性與精度的常用方法.下文通過力學(xué)系統(tǒng)(37)來探討算法性質(zhì).令式(38)中 ω=2π.若 0 ≤ρ ≤1,則Hamel 廣義 α 方法無條件穩(wěn)定;在此基礎(chǔ)上,若 ρ ＜1,該方法是數(shù)值耗散的,當(dāng)ρ=1 時(shí),該方法是非耗散的.定義數(shù)值阻尼比和周期延長(zhǎng)率為

圖1 展示了Hamel 廣義 α 方法的譜半徑,發(fā)現(xiàn)是無條件穩(wěn)定的,并且在高頻區(qū)域內(nèi)可通過 ρ∞精確控制數(shù)值耗散程度.

圖1 譜半徑Fig.1 Spectral radius

在0 ≤ρ∞≤1 情況下,當(dāng) Ω ≤0.1 時(shí),譜半徑接近于1,保持了基本的低頻模式,隨后快速降至 ρ∞,過濾掉不需要的高頻信息.圖2 和圖3 分別給出了該方法的數(shù)值阻尼比和周期延長(zhǎng)率,可以發(fā)現(xiàn)振幅和相位精度隨著 ρ∞的增加而提高.

圖2 數(shù)值阻尼比Fig.2 Numerical damping ratio

圖3 周期延長(zhǎng)率Fig.3 Period elongation

3.3 算法收斂性

本節(jié)通過數(shù)值算例測(cè)試算法的收斂速度.收斂速度是特定時(shí)間內(nèi)隨著時(shí)間步長(zhǎng)遞減,絕對(duì)誤差的下降率.令系統(tǒng)方程(38)式中 ω=π.初始條件為

該系統(tǒng)的精確解為

通過Hamel-廣義 α 方法t=1 時(shí)刻的絕對(duì)誤差,如圖4可以發(fā)現(xiàn)本文方法關(guān)于位移是二階精度,而且隨著ρ∞增大,算法準(zhǔn)確度也在提升.

圖4 Hamel 廣義 α 方法t=1 時(shí)刻的誤差Fig.4 Errors of Hamel generalized-α method at time t=1

若A=0.5,A3=-1,A4=1,B1=B2=0,通過Hamel廣義 α 方法t=1 時(shí)刻絕對(duì)誤差,如圖5,發(fā)現(xiàn)注2 的數(shù)值積分格式關(guān)于位移接近于四階精度.

圖5 Hamel 廣義 α 方法t=1 時(shí)刻的誤差Fig.5 Errors of Hamel generalized-α method at time t=1

4 數(shù)值算例

本節(jié)通過數(shù)值算例,在 ρ∞=0.8 的情況下將Hamel 廣義 α 方法與廣義α 方法,Newmark 方法(β=0.25,γ=0.5)的結(jié)果做對(duì)比,來驗(yàn)證其性質(zhì).

考慮一個(gè)兩自由度系統(tǒng),其Lagrange 函數(shù)為

系數(shù)矩陣為

其中m1=m2=1 和k1=104,k2=1.

相應(yīng)Euler-Lagrange 方程組為

采用如下初值

圖6 和圖7 分別給出了系統(tǒng)1 和系統(tǒng)2 三種方法的運(yùn)行結(jié)果,通過位移、速度和加速度曲線,可以看到Hamel 廣義 α 方法和廣義 α 方法在穩(wěn)定性方面要比Newmark 方法更具優(yōu)勢(shì).圖8 和圖9 分別給出了系統(tǒng)1 和系統(tǒng)2 Hamel 廣義 α 方法和廣義 α 方法兩種方法的數(shù)值誤差,通過與經(jīng)典算法對(duì)比發(fā)現(xiàn),本文提出的方法具有良好的穩(wěn)定性和精度性質(zhì).

圖6 系統(tǒng)1 各種方法數(shù)值結(jié)果對(duì)比圖Fig.6 Comparison of numerical results of various methods for system 1

圖7 系統(tǒng)2 各種方法數(shù)值結(jié)果對(duì)比圖Fig.7 Comparison of numerical results of various methods for system 2

圖8 系統(tǒng)1 Hamel 廣義 α 方法與廣義 α 方法數(shù)值結(jié)果誤差Fig.8 Errors of numerical results between Hamel generalized-α method and generalized-α method for system 1

圖9 系統(tǒng)2 Hamel 廣義 α 方法與廣義 α 方法數(shù)值結(jié)果誤差Fig.9 Errors of numerical results between Hamel generalized-α method and generalized-α method for system 2

5 李群空間上的Hamel-廣義 α 方法

本節(jié)考慮李群空間中力學(xué)系統(tǒng)的Hamel 形式,并提出相應(yīng)的運(yùn)動(dòng)方程—–Hamel 方程,在此基礎(chǔ)上,構(gòu)造李群空間上力學(xué)系統(tǒng)的時(shí)間積分算法.

5.1 Hamel 方程

考慮李群空間上的力學(xué)系統(tǒng),其位形流形為李群G,狀態(tài)空間為TG,T*G表示相空間,定義系統(tǒng)的Lagrange 函數(shù)L:TG→R.對(duì) ?g∈G令標(biāo)架算子Ψg為左不變?nèi)篢Lg,因此,定義活動(dòng)標(biāo)架下的Lagrange函數(shù)l:G×g→R

給定非保守力

其中fτ,fext,fcon分別為反饋控制力,外力和約束力

λ為L(zhǎng)agrange 乘子,Φ(g,ξ)=0 為系統(tǒng)的約束方程.由前文定理1,可得出下面的表述是等價(jià)的.

(1) 定義在G×(g⊕g*) 上的曲線 (g(t),ξ(t),μ(t)) 是作用泛函

關(guān)于獨(dú)立變分 η ,δ ξ 和 δμ 的臨界點(diǎn),即

(2) 曲線 (g(t),ξ(t),μ(t)) 滿足隱式Hamel 方程組的弱形式

曲線 (g,ξ,μ) 滿足隱式Hamel 方程組的強(qiáng)形式

證明參見文獻(xiàn)[32].

由式(53)可定義加速度函數(shù)

設(shè)計(jì)控制力fτ可獲得穩(wěn)定的時(shí)間積分?jǐn)?shù)值算法.

5.2 李群空間上的Hamel 廣義 α 方法

給定時(shí)間步長(zhǎng) Δt,區(qū)間[ 0,T] 等分為N個(gè)子區(qū)間,為簡(jiǎn)單起見,令區(qū)間 [tn,tn+1] 為 [0,Δt],為構(gòu)建李群上數(shù)值積分算法,定義G×g*上局部坐標(biāo)的離散形式

其中變量 Δgn表示時(shí)間區(qū)間 [tn,tn+1] 的平均速度場(chǎng),同時(shí)令

并且滿足 μ0=Mξ0,最后根據(jù)加速度函數(shù)給出相應(yīng)的離散加速度

式中外力,約束力和反饋控制律離散形式為式(31).

定義輔助未知量 ηs(t) 和 μs(t) 的離散形式

將離散形式(55)～(59)代入方程組(51),求解關(guān)于的駐點(diǎn)方程,其中關(guān)于的駐點(diǎn)方程可推導(dǎo)出

結(jié)合方程(55)可得出

其次,對(duì)于無約束系統(tǒng)tn+1時(shí)刻離散加速度函數(shù)為

若為帶約束系統(tǒng),則

為了得到更穩(wěn)定的算法,令

綜上所述,方程(60)～(65)為一般李群上的Hamel廣義 α 方法.

注: Hamel 廣義 α 方法可以退化為李群廣義 α 方法等經(jīng)典算法(見表2).令A(yù)3=2(3β-1),A4=-3(4β-1),表格中參數(shù)αm,αf參見文獻(xiàn)[16].

表2 Hamel 廣義 α 方法與經(jīng)典算法比較Table 2 Comparison between Hamel generalized-α method and classical algorithm

5.3 數(shù)值算例

本節(jié)通過空間雙擺模型[18,40]的數(shù)值仿真來說明Hamel 廣義 α 方法的性能,模型采用的是一階SE(3)幾何精確梁?jiǎn)卧?最終使用Matlab2016 b 實(shí)現(xiàn)了Hamel-廣義 α 方法.

空間雙擺模型初始放置于XOY水平面上,如圖10 所示,僅在重力作用下運(yùn)動(dòng).兩個(gè)擺臂的長(zhǎng)度分別為0.3 m 和0.5 m,截面均為寬與高0.5 mm 的矩形,物理參數(shù)為:ρ=2700 kg/m3,E=2.1×1010Pa,v=0.3,重力加速度為g=9.81 m/s2.我們采用廣義α方法作為參考,步長(zhǎng)均為5×10-4s,仿真時(shí)長(zhǎng)1 s,算法參數(shù)均為 ρ∞=0.8,并且對(duì)兩種方法的數(shù)值結(jié)果做了對(duì)比,發(fā)現(xiàn)與傳統(tǒng)廣義 α 方法結(jié)果保持一致,如圖11 和圖12,說明本方法適用于李群空間上的系統(tǒng).

圖10 柔性雙擺系統(tǒng)的初始構(gòu)型Fig.10 Initial configuration of a flexible double pendulum

圖11 柔性雙擺系統(tǒng)兩桿連接處軌跡Fig.11 Trajectory at the joint of a flexible double pendulum

圖12 柔性雙擺系統(tǒng)第二桿末端軌跡Fig.12 Trajectory at the end of a flexible double pendulum

6 結(jié)論

本文基于活動(dòng)標(biāo)架法和Hamel 場(chǎng)變分積分子,討論如何發(fā)展新的數(shù)值積分算法的構(gòu)造方法,并提出了一種無條件穩(wěn)定的Hamel 廣義 α 方法.該方法具有如下特點(diǎn): 無條件穩(wěn)定,二階精度和可以控制數(shù)值耗散;選擇特殊算法參數(shù),可以退化成傳統(tǒng)的廣義α方法、Newmark 等;由于活動(dòng)標(biāo)架法的特性,該方法既可以用于處理一般線性空間上的問題,也能夠用于處理李群空間上的問題;通過構(gòu)造特殊的變分形式,該方法可以演化出更高精度的數(shù)值格式.本文通過理論分析,給出了算法滿足精度,數(shù)值耗散等條件,指出該算法繼承了傳統(tǒng)數(shù)值積分算法的特點(diǎn),可以通過單個(gè)參數(shù)來控制數(shù)值耗散,可實(shí)現(xiàn)高頻損耗并同時(shí)最小化低頻損耗的目的;甚至可以提出更高階精度的數(shù)值格式.算法也可以從線性空間拓展到李群空間.通過和Newmark 方法、廣義 α 方法比較的結(jié)果表明,本文方法具有理想的精度和穩(wěn)定性,并且可以實(shí)現(xiàn)更高的精度.與李群廣義 α 方法相比,本文方法可以實(shí)現(xiàn)同等效果,這些數(shù)值算例驗(yàn)證了該方法的上述特點(diǎn)和可行性.

目前,上述方法和算法的研究仍需進(jìn)一步完善,進(jìn)而提出具無條件穩(wěn)定性,高階精度和可以控制數(shù)值耗散以及可變時(shí)間步長(zhǎng)的數(shù)值積分方法,同時(shí)發(fā)展出相應(yīng)的李群形式.

致謝

感謝北京理工大學(xué)宇航學(xué)院季奕博士在廣義α方法方面提供的幫助.