基于Bathe積分算法的機械系統多體動力學方程求解方法

吉 磊， 錢林方， 陳光宋， 尹 強

(南京理工大學 機械工程學院，南京 210094)

多體系統動力學有多種建模手段[1-3]，如絕對坐標法[4]、相對坐標法[5]、凱恩方法[6]和遞推方法[7]，很多學者[8-10]對這些建模方法進行了研究.對于一般含完整約束的多體系統，采用上述任一方法并引用拉格朗日乘子，均可建立一組微分代數方程(DAEs)，從而利用各種數值方法進行求解.然而DAEs不同的指標或形式有可能對數值方法的求解與效率產生影響，因此當采用不同的數值方法時，需要考慮DAEs的形式以適應算法，從而高效準確求解.

多體系統動力學方程的數值求解算法一直是多體問題中的關鍵內容，國內外學者研究了各種數值算法用于多體動力學方程的求解：Negrut等[11-12]提出了求解指標-3的多體系統動力學方程的Newmark方法和HHT-I3方法；Jay等[13]對廣義-α法進行了擴展研究，使其能夠對包含變質量矩陣、完整約束和非完整約束的多體系統進行求解運算，同時還提出了新的變步長公式；Lunk等[14]擴展了Newmark方法和廣義-α法類的方法用于求解含約束的機械系統，通過聯立位移和速度的穩定項從而保證約束正確，使用的α-RATTLE方法具有二階精度；Shabana等[15]提出了使用Newmark法對獨立坐標積分的雙循環隱式積分方法.馬秀騰等[16-17]利用θ1方法分別求解了指標-3和指標-2的多體動力學方程，表明該方法具有2階精度，并對HHT-SI2方法進行了改進，將其校正項通過兩種不同方式進行改進，數值算例表明改進方法都具有二階精度且數值阻尼可控；郭晛等[18]對HHT-α法在求解含接觸約束的柔性多體系統動力學方程時的數值特性進行了研究；姚廷強等[19]運動廣義-α方法計算分析了不同工況下圓柱滾子軸承的動態特性和保持架的穩定性，獲得了不同工況下軸承動態響應的變化規律；丁潔玉等[20]基于廣義-α法，結合約束投影方法構造了廣義-α-S投影法，該方法能較好地保持系統總能量和滿足約束，計算效率較高；劉穎等[21]提出了一種基于離散零空間理論的Newmark積分算法，結果表明該算法能夠實現系統降維并提升效率.上述這些常用的求解算法已被用于多體系統動力學分析，而Bathe積分策略[22-23]作為一種二階隱式積分算法，最初主要針對結構動力學常微分方程組(ODEs)的積分求解，而其亦可推廣到對DAEs的求解.Bathe積分算法的優點在于求解長時間歷程的非線性動力學問題時比較穩定，數值耗散較小，在較大積分步長的情況下也能獲得較好的精度.

本文將傳統的多體系統動力學方程整理為顯含廣義阻尼矩陣的形式，為Bathe積分算法提供高效計算的基礎.采用Baumgarte約束穩定法在動力學方程中添加了約束穩定項，減小約束違約對系統響應的影響.推導了利用Bathe積分算法求解動力學方程的計算過程，將廣義阻尼矩陣應用于迭代求解時的初值計算以提高效率.最后利用算例驗證了Bathe積分算法求解多體系統動力學方程的準確性和使用顯含廣義阻尼矩陣形式計算的高效性.

1 多體系統動力學方程一般形式

1.1 多體系統運動學方程

建立多體系統的運動學方程[24]，首先需推導兩個相鄰剛體的運動學關系.兩個相鄰剛體之間的運動關系如圖1所示，圖中O0x0y0z0為參考坐標系，Oixiyizi和Ojxjyjzj為分別為Bi和Bj剛體上的連體坐標系，piOξiOηiOζiO和pjIξjIηjIζjI為鉸Jj分別與Bi和Bj兩個剛體相連處的坐標系.

圖1 相鄰剛體位置關系

假設體Bi和Bj的連體坐標系Oixiyizi和Ojxjyjzj相對參考坐標系O0x0y0z0的速度與加速度矢量分別為

(1)

根據相對坐標法可遞推得到體坐標系運動矢量與鉸坐標系相對運動矢量之間的關系：

(2)

(3)

(4)

(5)

1.2 多體系統動力學方程

對于任意剛體Bi，假定ci為體內任意一點相對連體坐標系Oixiyizi的矢徑，則可得該點的加速度為

aci=ai+εi×ci+ωi×(ωi×ci)

(6)

式中：ai為Oi相對參考坐標系O0x0y0z0的加速度.

根據虛功率原理，可得到該體的虛功率方程為

(7)

δvci=δvi+δ(ωi×ci)

(8)

將式(1)、(2)、(6)分別代入式(7)，并將單體的虛功率方程運用至多體系統中，整理成緊湊形式得

(9)

式中：R、W、E、Qa和Qn為式(7)中各項在整個多體系統中的整合形式.

將式(4)和(5)代入式(9)可得

(10)

由于約束力和力矩不做功，所以其虛功率為0，將式(10)寫成緊湊形式得

(11)

(12)

根據變分原理可得系統動力學方程為

(13)

當多體系統中包含完整約束時，通過引入拉格朗日乘子λ寫出帶約束的動力學方程：

(14)

式中：Φ(q,t)=0為系統約束方程；Φq為Φ(q,t)對q求導.將式(14)中的約束方程Φ(q,t)=0對時間t求二階導數，可將式(14)寫為如下指標-1的微分-代數方程組：

(15)

1.3 Baumgarte約束違約穩定法

(16)

(17)

式中：λ″為λ對時間的二次積分項(λ″并未參與計算，引入λ″同樣是為了保持方程形式的一致性)，則式(16)可整理為

(18)

2 Bathe積分用于求解多體系統動力學方程

Bathe積分的求解思路是已知t時刻的運動參數和計算步長Δt，待求t+Δt時刻的響應，通過引入參數γ，先計算t+γΔt時刻的響應，再利用t和t+γΔt時刻的響應計算t+Δt時刻的響應，從而完成一個時間步長的求解.

(19)

(20)

(21)

(22)

t+γΔt時刻的動力學方程為

(23)

將式(21)和(22)代入式(23)，可得

(24)

即為t+γΔt時刻以Pt+γΔt為自變量的動力學方程.

獲得t+γΔt時刻的運動參數后即可計算t+Δt時刻的運動參數.在t+Δt時刻的函數f的導數可寫為t、t+γΔt和t+Δt時刻函數的組合形式[22]：

(25)

(26)

將式(27)代入式(28)，可得

c2c3Pt+γΔt+c3c3Pt+Δt

(29)

t+Δt時刻的動力學方程為

(30)

將式(27)和(29)代入式(30)可得

(31)

即為t+Δt時刻以Pt+Δt為自變量的動力學方程.

式(24)和(31)都可通過Newton迭代方法進行求解，為提高計算效率采用Broyden擬牛頓法，迭代過程如下：

(32)

式中：上標k和k+1分別表示第k次和第k+1次迭代；下標tCalu表示當前計算的時刻(tCalu=t+γΔt或tCalu=t+Δt)；yk和sk都為迭代的中間變量；J表示雅克比矩陣；Y的表達式如下：

(33)

根據泰勒展開原理，式(24)的迭代初值可設置為

同理，式(31)的迭代初值可設置為

(36)

(37)

3 數值算例

3.1 曲柄滑塊機構

下面以曲柄滑塊機構為例，驗證Bathe積分算法求解多體系統動力學方程的準確性和高效性.如圖2所示，圖中曲柄、連桿和接觸體均為均質實心圓柱體，點O1、O2和O3分別為曲柄、連桿和滑塊的質心.Oxy為參考坐標系，局部坐標系O1x1y1、O2x2y2和O3x3y3與3個構件質心固接，θ1和θ2為曲柄和連桿的轉動角.滑塊末端與彈簧阻尼機構相連，剛度為k，阻尼為c.當θ1和θ2為0時，該彈簧處于平衡位置.圖中l1和l2分別表示曲柄和連桿的長度，m1、m2和m3分別表示曲柄、連桿和滑塊的質量，各參數取值如表1所示.假設系統只受重力作用，重力加速度g=9.806 65 m/s2.

圖2 曲柄滑塊示意圖

表1 曲柄滑塊參數列表

由于本算例無解析解，所以本文利用Adams軟件中的WSTIFF積分器SI2算法作為求解算法，誤差設置為10-10，步長設置為10-6，以此情況下的結果作為參考解.下面給出了Bathe算法與Adams計算出的參考解對比，圖3和圖4分別為θ1、θ2隨時間變化的曲線圖.

圖3 θ1隨時間變化曲線

圖4 θ2隨時間變化曲線

圖5 位移約束方程違約值

為了研究Bathe積分算法對系統總能量(包括動能、重力勢能以及彈簧的彈性變形勢能)的影響，假設彈簧阻尼c=0，圖7顯示了不同γ取值(0.1、0.3、0.5、0.7和0.9)時系統總能量保持的情況，從圖中可看出γ=0.1與γ=0.9時較為一致，γ=0.3與γ=0.7時較為一致.不論γ取值，總能量的變化范圍都較小，Bathe積分算法較穩定且數值耗散較小，在本文中，如不特別說明，γ取值為0.5.

圖6 速度約束方程違約值曲線

圖7 系統總能量隨時間變化曲線

表2 雅克比矩陣的初值J0中是否含有廣義阻尼矩陣對Bathe算法總迭代次數的影響

使用Adams軟件以10-6為時間步長時的結果作為參考值，對比Bathe算法與龍格庫塔45階算法(RK45)對θ1和θ2的計算結果在不同時間步長Δt(0.05～10-5s)情況下的相對誤差e情況以及相應CPU時間(tCPU)，對比結果如圖8～10所示(圖中均為雙對數坐標).

圖8 θ1的相對誤差隨時間步長的變化

圖9 θ2的相對誤差隨時間步長的變化

從圖8和9中可以看出：在相同的時間步長下，Bathe算法計算結果的相對誤差較小，收斂效果更好；且隨著步長的減小，Bathe算法相對誤差下降的趨勢更加顯著.從圖10中可以看出：相同時間步長下Bathe算法的計算耗時比RK45算法稍長，這是由于Bathe算法每個時間步長都需要進行迭代計算至收斂的原因；且隨著時間步長的減小，兩種算法計算的CPU時間有相互接近的趨勢.若兼顧精確度與計算CPU時間，選擇Bathe算法進行動力學計算對計算耗時的略微犧牲能獲得更精確的計算結果.

圖10 計算的CPU時間隨時間步長的變化

3.2 含接觸的曲柄滑塊

以含接觸的曲柄滑塊機構為例，驗證Bathe積分求解算法對含接觸的多體系統同樣適用.如圖11所示，曲柄滑塊機構與算例3.1中的相同，在此基礎上添加一接觸體，接觸體與地面固連，在系統運動過程中與連桿發生接觸，接觸體為均質實心圓柱，點O4為其質心，坐標系O4x4y4z4固接于接觸體質心，r4、l4和m4分別表示接觸體的半徑、長度和質量.選擇使用純彈性接觸力模型模擬剛體之間的接觸碰撞過程.接觸力Fc的模型[27]為

圖11 含接觸的曲柄滑塊示意圖

(38)

式中：δ為穿透深度；Kc為接觸剛度；nc為接觸指數.算例中參數取值如表3所示(其余參數與表1一致)，假設系統只受重力作用.

表3 接觸體及接觸力參數

采用本文的建模方法建立動力學方程，并采用Bathe算法計算0～2 s內該系統的響應，Bathe積分策略參數λ=0.5，步長取10-5s，初始狀態時θ1和θ2為0.下面給出了Bathe算法與Adams軟件計算的結果進行對比，圖12和13分別為θ1和θ2隨時間變化的曲線圖.同算例3.1一樣，從圖12和13中可以看出Bathe積分算法求解的結果與Adams軟件求解的結果是一致的.圖14和15給出了接觸力在x方向和y方向的曲線圖，顯示兩種算法結果一致，每次接觸開始和結束的時刻以及力的大小都一致，驗證了Bathe積分策略用于求解含接觸多體動力學方程同樣是有效和準確的.

圖12 θ1隨時間變化曲線

圖13 θ2隨時間變化曲線

圖14 接觸力在x方向上隨時間變化曲線

圖15 接觸力在y方向上隨時間變化曲線

4 結論

本文主要研究了Bathe積分算法在多體系統動力學中的應用.將多體系統動力學方程整理成顯含廣義阻尼矩陣的形式，并在動力學方程中添加了Bamugarte違約穩定項，推導了Bathe積分算法求解多體系統動力學方程的具體流程，將廣義阻尼矩陣用于迭代計算時選取雅克比矩陣的初值.通過數值算例，得到如下結論：

(1)利用Bathe積分算法求解多體系統動力學方程具有準確的結果，約束違約值很小，算法穩定性較好.Bathe積分算法求解時的數值耗散小，不同γ參數時總能量的變化規律略有區別.

(2)由于將動力學方程整理為顯含廣義阻尼矩陣的形式，改進了雅克比矩陣迭代初值的選擇，提高了Bathe積分算法的計算效率，且隨著步長的增大計算效率相對越高.

(3)Bathe算法在相同步長情況下的計算結果更加準確，收斂性更好.雖然Bathe算法計算時的CPU時間會略微變長，但計算結果的準確性會大幅度提高.