一種基于Hamel 形式的無條件穩定動力學積分算法1)

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

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

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

** (北京應用物理與計算數學研究所,北京 100088)

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

*** (北京理工大學數學與統計學院,北京 100081)

??? (復雜信息數學表征分析與應用北京市重點實驗室,北京 100081)

引言

時間積分算法被廣泛應用于動力學系統的求解中.動力學方程的時間積分經常會出現數值不穩定現象,有限元空間離散也通常會造成偽高頻振蕩,因此,針對上述問題的數值積分算法的研究一直備受關注.本文將基于活動標架法和離散力學,發展對動力學系統無條件穩定的數值積分格式.

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

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

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

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

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

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

1 Hamel 形式與Hamel 變分積分子

1.1 Hamel 場形式

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

給定兩個向量 ξ,η ∈W,定義一個反對稱括號算子[·,·]q:W×W→W

其中 [·,·] 是流形上兩個向量場的Jacobi-Lie 括號.對于任意 ξ,η ∈W和 α ∈W*,括號算子 [·,·]q的對偶算子通過下式給出

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

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

令力學系統位形空間為Q,TQ是切叢,T*Q為余切叢.令q∈Q,(q,v)∈TQ且 (q,p)∈T*Q.令L:TQ→R為Lagrange 函數,l是關于活動標架下局部坐標(q,ξ)的Lagrange 函數,則下列陳述等價.

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

關于獨立變分 δq,δv,δp的臨界點,即

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

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

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

關于獨立變分 δq=Ψqη ,δ ξ 和 δμ 的臨界點

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

對于任意點q∈Q,曲線 (q,ξ,μ) 滿足隱式Hamel 方程組的強形式

1.2 Hamel 變分積分子

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

式中參數 α ∈[0,1],Δt是時間步長,令

為關于活動標架 Ψqk+α的速度分量.

通過連續Lagrange 函數,定義其離散形式

相應離散作用和為

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

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

其中Dil為Lagrange 函數對第i個參數的偏導數.

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

本節考慮向量空間中力學系統,基于活動標架法,構造力學系統的Hamel 廣義 α 方法.

2.1 Hamel 方程

考慮一個受控系統,其位形空間為向量空間Q,狀態空間為TQ,T*Q表示相空間.對 ?q∈Q取標架算子 Ψq為恒算子,定義系統的Lagrange 函數L:TQ→R

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

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

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

關于獨立變分 δq,δv,δp的臨界點,即

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

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

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

2.2 Hamel-廣義 α 方法

給定時間步長 Δt,區間[ 0,T] 等分為N個子區間,為簡單起見,令區間 [tn,tn+1] 為 [0,Δt],為構建數值積分算法,定義離散形式

并且滿足Mv0=p0,最后根據加速度函數定義相應的離散加速度

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

因此,可知力f所對應的離散形式為

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

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

式中參數分別為

其次,tn+1時刻離散加速度函數為

為了得到更穩定的算法,令

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

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

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

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

可得到更高階數值積分算法(33),其中參數μ0,μ1,γ0,γ1變為

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

3 算法性質

3.1 算法分析

考慮一個力學系統

相應的Euler-Lagrange 方程組為

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

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

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

其中 λi為矩陣的第i個特征值.若 ρ ≤1,則對于線性問題算法是無條件穩定的.

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

其中a和b是實數,只要

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

滿足這個條件,Hamel 廣義 α 方法便可實現高頻耗散最大化.此外,當A=0.5 時,可得A4=0,無法確定A3的值,故而設定A3=-0.5.

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

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

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

3.2 算法穩定性

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

圖1 展示了Hamel 廣義 α 方法的譜半徑,發現是無條件穩定的,并且在高頻區域內可通過 ρ∞精確控制數值耗散程度.

圖1 譜半徑Fig.1 Spectral radius

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

圖2 數值阻尼比Fig.2 Numerical damping ratio

圖3 周期延長率Fig.3 Period elongation

3.3 算法收斂性

本節通過數值算例測試算法的收斂速度.收斂速度是特定時間內隨著時間步長遞減,絕對誤差的下降率.令系統方程(38)式中 ω=π.初始條件為

該系統的精確解為

通過Hamel-廣義 α 方法t=1 時刻的絕對誤差,如圖4可以發現本文方法關于位移是二階精度,而且隨著ρ∞增大,算法準確度也在提升.

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

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

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

4 數值算例

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

考慮一個兩自由度系統,其Lagrange 函數為

系數矩陣為

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

相應Euler-Lagrange 方程組為

采用如下初值

圖6 和圖7 分別給出了系統1 和系統2 三種方法的運行結果,通過位移、速度和加速度曲線,可以看到Hamel 廣義 α 方法和廣義 α 方法在穩定性方面要比Newmark 方法更具優勢.圖8 和圖9 分別給出了系統1 和系統2 Hamel 廣義 α 方法和廣義 α 方法兩種方法的數值誤差,通過與經典算法對比發現,本文提出的方法具有良好的穩定性和精度性質.

圖6 系統1 各種方法數值結果對比圖Fig.6 Comparison of numerical results of various methods for system 1

圖7 系統2 各種方法數值結果對比圖Fig.7 Comparison of numerical results of various methods for system 2

圖8 系統1 Hamel 廣義 α 方法與廣義 α 方法數值結果誤差Fig.8 Errors of numerical results between Hamel generalized-α method and generalized-α method for system 1

圖9 系統2 Hamel 廣義 α 方法與廣義 α 方法數值結果誤差Fig.9 Errors of numerical results between Hamel generalized-α method and generalized-α method for system 2

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

本節考慮李群空間中力學系統的Hamel 形式,并提出相應的運動方程—–Hamel 方程,在此基礎上,構造李群空間上力學系統的時間積分算法.

5.1 Hamel 方程

考慮李群空間上的力學系統,其位形流形為李群G,狀態空間為TG,T*G表示相空間,定義系統的Lagrange 函數L:TG→R.對 ?g∈G令標架算子Ψg為左不變群TLg,因此,定義活動標架下的Lagrange函數l:G×g→R

給定非保守力

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

λ為Lagrange 乘子,Φ(g,ξ)=0 為系統的約束方程.由前文定理1,可得出下面的表述是等價的.

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

關于獨立變分 η ,δ ξ 和 δμ 的臨界點,即

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

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

證明參見文獻[32].

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

設計控制力fτ可獲得穩定的時間積分數值算法.

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

給定時間步長 Δt,區間[ 0,T] 等分為N個子區間,為簡單起見,令區間 [tn,tn+1] 為 [0,Δt],為構建李群上數值積分算法,定義G×g*上局部坐標的離散形式

其中變量 Δgn表示時間區間 [tn,tn+1] 的平均速度場,同時令

并且滿足 μ0=Mξ0,最后根據加速度函數給出相應的離散加速度

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

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

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

結合方程(55)可得出

其次,對于無約束系統tn+1時刻離散加速度函數為

若為帶約束系統,則

為了得到更穩定的算法,令

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

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

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

5.3 數值算例

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

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

圖10 柔性雙擺系統的初始構型Fig.10 Initial configuration of a flexible double pendulum

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

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

6 結論

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

目前,上述方法和算法的研究仍需進一步完善,進而提出具無條件穩定性,高階精度和可以控制數值耗散以及可變時間步長的數值積分方法,同時發展出相應的李群形式.

致謝

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