吉 磊， 錢林方， 陳光宋， 尹 強(qiáng)

(南京理工大學(xué) 機(jī)械工程學(xué)院，南京 210094)

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

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

本文將傳統(tǒng)的多體系統(tǒng)動(dòng)力學(xué)方程整理為顯含廣義阻尼矩陣的形式，為Bathe積分算法提供高效計(jì)算的基礎(chǔ).采用Baumgarte約束穩(wěn)定法在動(dòng)力學(xué)方程中添加了約束穩(wěn)定項(xiàng)，減小約束違約對系統(tǒng)響應(yīng)的影響.推導(dǎo)了利用Bathe積分算法求解動(dòng)力學(xué)方程的計(jì)算過程，將廣義阻尼矩陣應(yīng)用于迭代求解時(shí)的初值計(jì)算以提高效率.最后利用算例驗(yàn)證了Bathe積分算法求解多體系統(tǒng)動(dòng)力學(xué)方程的準(zhǔn)確性和使用顯含廣義阻尼矩陣形式計(jì)算的高效性.

1 多體系統(tǒng)動(dòng)力學(xué)方程一般形式

1.1 多體系統(tǒng)運(yùn)動(dòng)學(xué)方程

建立多體系統(tǒng)的運(yùn)動(dòng)學(xué)方程[24]，首先需推導(dǎo)兩個(gè)相鄰剛體的運(yùn)動(dòng)學(xué)關(guān)系.兩個(gè)相鄰剛體之間的運(yùn)動(dòng)關(guān)系如圖1所示，圖中O0x0y0z0為參考坐標(biāo)系，Oixiyizi和Ojxjyjzj為分別為Bi和Bj剛體上的連體坐標(biāo)系，piOξiOηiOζiO和pjIξjIηjIζjI為鉸Jj分別與Bi和Bj兩個(gè)剛體相連處的坐標(biāo)系.

圖1 相鄰剛體位置關(guān)系

假設(shè)體Bi和Bj的連體坐標(biāo)系Oixiyizi和Ojxjyjzj相對參考坐標(biāo)系O0x0y0z0的速度與加速度矢量分別為

(1)

根據(jù)相對坐標(biāo)法可遞推得到體坐標(biāo)系運(yùn)動(dòng)矢量與鉸坐標(biāo)系相對運(yùn)動(dòng)矢量之間的關(guān)系：

(2)

(3)

(4)

(5)

1.2 多體系統(tǒng)動(dòng)力學(xué)方程

對于任意剛體Bi，假定ci為體內(nèi)任意一點(diǎn)相對連體坐標(biāo)系Oixiyizi的矢徑，則可得該點(diǎn)的加速度為

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

(6)

式中：ai為Oi相對參考坐標(biāo)系O0x0y0z0的加速度.

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

(7)

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

(8)

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

(9)

式中：R、W、E、Qa和Qn為式(7)中各項(xiàng)在整個(gè)多體系統(tǒng)中的整合形式.

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

(10)

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

(11)

(12)

根據(jù)變分原理可得系統(tǒng)動(dòng)力學(xué)方程為

(13)

當(dāng)多體系統(tǒng)中包含完整約束時(shí)，通過引入拉格朗日乘子λ寫出帶約束的動(dòng)力學(xué)方程：

(14)

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

(15)

1.3 Baumgarte約束違約穩(wěn)定法

(16)

(17)

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

(18)

2 Bathe積分用于求解多體系統(tǒng)動(dòng)力學(xué)方程

Bathe積分的求解思路是已知t時(shí)刻的運(yùn)動(dòng)參數(shù)和計(jì)算步長Δt，待求t+Δt時(shí)刻的響應(yīng)，通過引入?yún)?shù)γ，先計(jì)算t+γΔt時(shí)刻的響應(yīng)，再利用t和t+γΔt時(shí)刻的響應(yīng)計(jì)算t+Δt時(shí)刻的響應(yīng)，從而完成一個(gè)時(shí)間步長的求解.

(19)

(20)

(21)

(22)

t+γΔt時(shí)刻的動(dòng)力學(xué)方程為

(23)

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

(24)

即為t+γΔt時(shí)刻以Pt+γΔt為自變量的動(dòng)力學(xué)方程.

獲得t+γΔt時(shí)刻的運(yùn)動(dòng)參數(shù)后即可計(jì)算t+Δt時(shí)刻的運(yùn)動(dòng)參數(shù).在t+Δt時(shí)刻的函數(shù)f的導(dǎo)數(shù)可寫為t、t+γΔt和t+Δt時(shí)刻函數(shù)的組合形式[22]：

(25)

(26)

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

c2c3Pt+γΔt+c3c3Pt+Δt

(29)

t+Δt時(shí)刻的動(dòng)力學(xué)方程為

(30)

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

(31)

即為t+Δt時(shí)刻以Pt+Δt為自變量的動(dòng)力學(xué)方程.

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

(32)

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

(33)

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

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

(36)

(37)

3 數(shù)值算例

3.1 曲柄滑塊機(jī)構(gòu)

下面以曲柄滑塊機(jī)構(gòu)為例，驗(yàn)證Bathe積分算法求解多體系統(tǒng)動(dòng)力學(xué)方程的準(zhǔn)確性和高效性.如圖2所示，圖中曲柄、連桿和接觸體均為均質(zhì)實(shí)心圓柱體，點(diǎn)O1、O2和O3分別為曲柄、連桿和滑塊的質(zhì)心.Oxy為參考坐標(biāo)系，局部坐標(biāo)系O1x1y1、O2x2y2和O3x3y3與3個(gè)構(gòu)件質(zhì)心固接，θ1和θ2為曲柄和連桿的轉(zhuǎn)動(dòng)角.滑塊末端與彈簧阻尼機(jī)構(gòu)相連，剛度為k，阻尼為c.當(dāng)θ1和θ2為0時(shí)，該彈簧處于平衡位置.圖中l(wèi)1和l2分別表示曲柄和連桿的長度，m1、m2和m3分別表示曲柄、連桿和滑塊的質(zhì)量，各參數(shù)取值如表1所示.假設(shè)系統(tǒng)只受重力作用，重力加速度g=9.806 65 m/s2.

圖2 曲柄滑塊示意圖

表1 曲柄滑塊參數(shù)列表

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

圖3 θ1隨時(shí)間變化曲線

圖4 θ2隨時(shí)間變化曲線

圖5 位移約束方程違約值

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

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

圖7 系統(tǒng)總能量隨時(shí)間變化曲線

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

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

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

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

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

圖10 計(jì)算的CPU時(shí)間隨時(shí)間步長的變化

3.2 含接觸的曲柄滑塊

以含接觸的曲柄滑塊機(jī)構(gòu)為例，驗(yàn)證Bathe積分求解算法對含接觸的多體系統(tǒng)同樣適用.如圖11所示，曲柄滑塊機(jī)構(gòu)與算例3.1中的相同，在此基礎(chǔ)上添加一接觸體，接觸體與地面固連，在系統(tǒng)運(yùn)動(dòng)過程中與連桿發(fā)生接觸，接觸體為均質(zhì)實(shí)心圓柱，點(diǎn)O4為其質(zhì)心，坐標(biāo)系O4x4y4z4固接于接觸體質(zhì)心，r4、l4和m4分別表示接觸體的半徑、長度和質(zhì)量.選擇使用純彈性接觸力模型模擬剛體之間的接觸碰撞過程.接觸力Fc的模型[27]為

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

(38)

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

表3 接觸體及接觸力參數(shù)

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

圖12 θ1隨時(shí)間變化曲線

圖13 θ2隨時(shí)間變化曲線

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

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

4 結(jié)論

本文主要研究了Bathe積分算法在多體系統(tǒng)動(dòng)力學(xué)中的應(yīng)用.將多體系統(tǒng)動(dòng)力學(xué)方程整理成顯含廣義阻尼矩陣的形式，并在動(dòng)力學(xué)方程中添加了Bamugarte違約穩(wěn)定項(xiàng)，推導(dǎo)了Bathe積分算法求解多體系統(tǒng)動(dòng)力學(xué)方程的具體流程，將廣義阻尼矩陣用于迭代計(jì)算時(shí)選取雅克比矩陣的初值.通過數(shù)值算例，得到如下結(jié)論：

(1)利用Bathe積分算法求解多體系統(tǒng)動(dòng)力學(xué)方程具有準(zhǔn)確的結(jié)果，約束違約值很小，算法穩(wěn)定性較好.Bathe積分算法求解時(shí)的數(shù)值耗散小，不同γ參數(shù)時(shí)總能量的變化規(guī)律略有區(qū)別.

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

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