闞子云 彭海軍 陳飆松

(大連理工大學(xué)工程力學(xué)系工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧大連116024)

考慮彈簧阻尼作動(dòng)器解析雅可比矩陣的多剛體動(dòng)力學(xué)分析1)

闞子云 彭海軍2)陳飆松

(大連理工大學(xué)工程力學(xué)系工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧大連116024)

彈簧--阻尼--作動(dòng)器(spring-damper-actuator，SDA)是多體系統(tǒng)中常見(jiàn)的力元，在工程領(lǐng)域中有著廣泛的應(yīng)用.采用絕對(duì)坐標(biāo)方法建立的多體系統(tǒng)動(dòng)力學(xué)控制方程通常是復(fù)雜的非線性微分--代數(shù)方程組.為了保證數(shù)值解的精度和穩(wěn)定性，通常需要采用隱式算法求解動(dòng)力學(xué)方程，而雅可比矩陣的計(jì)算在隱式數(shù)值求解過(guò)程中至關(guān)重要.對(duì)于含有SDA的多體系統(tǒng)，SDA造成的附加雅可比矩陣是與廣義坐標(biāo)和廣義速度相關(guān)的高度非線性函數(shù).目前的很多研究工作專(zhuān)注于廣義力向量的計(jì)算，然而對(duì)附加雅克比矩陣的計(jì)算則少有關(guān)注.針對(duì)含SDA的多剛體系統(tǒng)進(jìn)行動(dòng)力學(xué)分析，首先基于Newmark算法研究其在動(dòng)力學(xué)方程求解中的雅可比矩陣的構(gòu)成形式;然后推導(dǎo)SDA的廣義力向量對(duì)應(yīng)的附加雅可比矩陣，其中包括廣義力向量對(duì)廣義坐標(biāo)和對(duì)廣義速度的偏導(dǎo)數(shù)矩陣.最后通過(guò)兩個(gè)數(shù)值算例研究附加雅可比矩陣對(duì)動(dòng)力學(xué)分析收斂性的影響;數(shù)值分析表明：當(dāng)SDA的剛度、阻尼和作動(dòng)力數(shù)值較大時(shí)，SDA導(dǎo)致的附加雅可比矩陣對(duì)數(shù)值解的收斂性有重要影響；當(dāng)考慮SDA對(duì)應(yīng)的附加雅可比矩陣時(shí)，動(dòng)力學(xué)分析可以以較少的迭代步實(shí)現(xiàn)收斂，從而減少分析時(shí)間.

彈簧阻尼作動(dòng)器，雅可比矩陣，多體系統(tǒng)，隱式算法，張拉整體結(jié)構(gòu)

引言

對(duì)復(fù)雜機(jī)械系統(tǒng)進(jìn)行動(dòng)力學(xué)分析前需要建立其多體系統(tǒng)動(dòng)力學(xué)模型，建模的實(shí)質(zhì)是對(duì)系統(tǒng)中4個(gè)要素進(jìn)行定義[1]：物體、鉸、主動(dòng)力、力元.其中力元反映了物體與物體間的相互作用.彈簧--阻尼--作動(dòng)器(spring-damper-actuator，SDA)系統(tǒng)是多體系統(tǒng)中的常見(jiàn)力元,可以有效地模擬工程中機(jī)構(gòu)或結(jié)構(gòu)的某些元器件的靜力和動(dòng)力特性,在機(jī)械、汽車(chē)、土木等領(lǐng)域得到廣泛的應(yīng)用[26].

采用絕對(duì)坐標(biāo)方法對(duì)多體系統(tǒng)進(jìn)行動(dòng)力學(xué)分析得到的控制方程通常是復(fù)雜的非線性微分--代數(shù)方程組.在數(shù)學(xué)上，微分方程的求解通?？煞譃轱@式算法和隱式算法兩類(lèi).顯式算法不需要迭代求解，算法穩(wěn)定域較小，因此需要較小的時(shí)間步長(zhǎng)才能滿足長(zhǎng)時(shí)間仿真的精度要求，通常不用于非線性程度較高且呈現(xiàn)剛性特點(diǎn)的多體系統(tǒng)[7]；而隱式算法采用迭代求解策略，在每一增量步內(nèi)都需要對(duì)離散非線性平衡方程進(jìn)行迭代求解，可以采用較大的時(shí)間步長(zhǎng)同時(shí)滿足解的精度和穩(wěn)定性要求.另一方面，多體動(dòng)力學(xué)方程中包含代數(shù)方程，其求解不同于一般的常微分方程組，根據(jù)約束方程的不同處理方式，又可以將微分--代數(shù)方程組的求解歸納為增廣法和縮并法兩大類(lèi).直接時(shí)間積分法是增廣法中的一類(lèi)重要求解算法，經(jīng)典的有Newmark法[8]、HHT法[910]、廣義α法[11]以及近年來(lái)出現(xiàn)的Bathe積分方法[1213]、祖沖之類(lèi)保辛算法[1415]等.這些算法在數(shù)學(xué)上均屬于隱式算法，在求解過(guò)程中同時(shí)對(duì)連續(xù)時(shí)間域的動(dòng)力學(xué)控制方程(包含約束方程)進(jìn)行離散，得到離散時(shí)間點(diǎn)處的非線性平衡方程，通過(guò)滾動(dòng)求解非線性平衡方程來(lái)完成整個(gè)時(shí)間域的分析.在每一時(shí)間步的非線性平衡方程的迭代求解中，需要計(jì)算系統(tǒng)的雅可比矩陣，類(lèi)似于非線性有限元領(lǐng)域的切線剛度陣或結(jié)構(gòu)動(dòng)力學(xué)中的等效剛度陣.SDA作為多體系統(tǒng)中的力元可以理解成系統(tǒng)中的內(nèi)力源之一，影響著多體動(dòng)力學(xué)方程中的廣義力向量的集成，并對(duì)系統(tǒng)的總體雅可比矩陣產(chǎn)生貢獻(xiàn).對(duì)應(yīng)SDA的雅可比矩陣則是與廣義坐標(biāo)以及廣義速度相關(guān)的復(fù)雜非線性函數(shù).而計(jì)算多體系統(tǒng)動(dòng)力學(xué)領(lǐng)域中的多部經(jīng)典專(zhuān)著[1,1619]在進(jìn)行相關(guān)介紹時(shí)，一般專(zhuān)注于廣義力向量的計(jì)算，而對(duì)SDA導(dǎo)致的附加雅可比矩陣的計(jì)算尚未見(jiàn)有完整系統(tǒng)的報(bào)道.雅可比矩陣的計(jì)算在多體動(dòng)力學(xué)分析中至關(guān)重要，當(dāng)不易解析地計(jì)算非線性方程的雅克比矩陣時(shí)，一般可以采用兩種處理方法：一是采用數(shù)值差分[20]，計(jì)算出離散非線性平衡方程的近似雅可比矩陣，這種方式簡(jiǎn)單易行但對(duì)于大規(guī)模問(wèn)題，效率十分低下；另一種折衷的辦法是選取易于解析計(jì)算的部分項(xiàng)參與總體雅可比矩陣的組裝.當(dāng)系統(tǒng)中含有SDA時(shí)，則忽略SDA產(chǎn)生的附加雅克比矩陣.然而在多體動(dòng)力學(xué)求解中不計(jì)入SDA造成的雅可比矩陣貢獻(xiàn)時(shí)，會(huì)使系統(tǒng)的總體雅可比矩陣不精確，最終可能會(huì)導(dǎo)致非線性求解過(guò)程中收斂速度變慢，甚至完全不能收斂.

鑒于此，本文首先以多體動(dòng)力學(xué)方程求解的Newmark算法為例，說(shuō)明其雅可比矩陣的具體構(gòu)成形式，以及 SDA對(duì)應(yīng)的廣義力向量.進(jìn)一步推導(dǎo)了SDA在動(dòng)力學(xué)分析中附加的雅可比矩陣的精確形式，包括廣義力向量對(duì)廣義坐標(biāo)以及廣義速度的偏導(dǎo)數(shù)矩陣.并在附錄中以姿態(tài)角的四元數(shù)描述為例，給出了雅可比矩體組成關(guān)鍵項(xiàng)的具體表達(dá)式.當(dāng)姿態(tài)角采用歐拉角或卡爾丹角等其他類(lèi)型描述時(shí)，只需計(jì)算出其對(duì)應(yīng)的關(guān)鍵項(xiàng)的表達(dá)式，進(jìn)行整體代入即可.最后通過(guò)數(shù)值算例研究了SDA產(chǎn)生的附加雅可比矩陣對(duì)多體系統(tǒng)動(dòng)力學(xué)分析收斂性的影響.

1 多體動(dòng)力學(xué)方程及其數(shù)值求解

典型的多體系統(tǒng)動(dòng)力學(xué)方程是如下形式的指標(biāo)-3的微分代數(shù)方程組

其中M為系統(tǒng)的廣義質(zhì)量矩陣，Φ和Φq分別為系統(tǒng)的約束方程及約束雅可比矩陣，Qi,Qa和Qe分別為系統(tǒng)的耦合慣性力、外加力、內(nèi)力.對(duì)于含有SDA的多體系統(tǒng)，其對(duì)系統(tǒng)總體廣義力向量的影響體現(xiàn)在內(nèi)力Qe中.此處以及后文中黑斜體表示矩陣或向量，斜體均表示標(biāo)量.

以經(jīng)典的Newmark算法為例，考察其雅可比矩陣組成.在時(shí)間區(qū)間[tn,tn+1]內(nèi)，Newmark算法采用如下假設(shè)公式

式中?t為時(shí)間步長(zhǎng)，α和δ為Newmark參數(shù)，帶頂標(biāo)“～”的量為當(dāng)前時(shí)間步的待求未知量.為了書(shū)寫(xiě)方便，將式(1)中第1項(xiàng)的等號(hào)右端項(xiàng)統(tǒng)稱(chēng)為總體廣義力向量Qt.動(dòng)力學(xué)方程離散后得到的系統(tǒng)的非線性代數(shù)方程組為

其中

當(dāng)多體系統(tǒng)中含SDA時(shí)，其將對(duì)總體廣義力向量Qt產(chǎn)生貢獻(xiàn)，具體通過(guò)式(4)中 P3項(xiàng)具體影響雅可比矩陣J的構(gòu)成.上述雖僅針對(duì)Newmark算法說(shuō)明其雅可比矩陣的構(gòu)成，對(duì)其他時(shí)間積分算法，雅可比矩陣的具體形式雖然略有不同，但是都不可避免地需要對(duì)矩陣?Qe/?q和進(jìn)行計(jì)算.此外，在多體系統(tǒng)的靜平衡分析中，當(dāng)忽略速度和加速度相關(guān)項(xiàng)后，獲取精確的雅可比矩陣仍需要計(jì)算?Qe/?q.

2 彈簧--阻尼--作動(dòng)器系統(tǒng)的動(dòng)力學(xué)模型

如圖1所示,物體I和物體J通過(guò)SDA連接，連接點(diǎn)分別為Pi和Pj.當(dāng)前時(shí)刻兩物體相對(duì)于全局坐標(biāo)系的位形坐標(biāo)(廣義坐標(biāo))為

廣義速度

圖1 SDA示意圖Fig.1 Schematic diagram of SDA

式 (8)和式 (9)中 Ri∈ R3×1和 Rj∈ R3×1分別表示固結(jié)在物體 I和 J上的局部坐標(biāo)系的原點(diǎn)相對(duì)于全局坐標(biāo)系的矢量在全局坐標(biāo)系的投影.Θi和Θj分別表示各自物體的局部坐標(biāo)系相對(duì)于全局坐標(biāo)系的姿態(tài)角.此種建模方式即為多體動(dòng)力學(xué)中最為經(jīng)典和工程中普遍使用的笛卡爾絕對(duì)建模方式.當(dāng)采用歐拉角或卡爾丹角描述姿態(tài)角Θi∈R3×1,Θj∈R3×1，一個(gè)物體共有6個(gè)廣義坐標(biāo).當(dāng)采用四元數(shù)描述姿態(tài)角時(shí)，Θi∈ R4×1，Θj∈ R4×1，一個(gè)物體共有7個(gè)廣義坐標(biāo).值得說(shuō)明的是：由于四元數(shù)之間并非相互獨(dú)立，此種情況下動(dòng)力學(xué)方程中需要添加相應(yīng)的約束方程.為統(tǒng)一起見(jiàn)，以下將單個(gè)物體對(duì)應(yīng)的廣義坐標(biāo)個(gè)數(shù)記為N.

記Pi和Pj點(diǎn)在全局坐標(biāo)系的坐標(biāo)分別為RPi∈R3×1和 RPj∈ R3×1，則有

將式(14)對(duì)時(shí)間求導(dǎo)，根據(jù)數(shù)學(xué)關(guān)系知

令

其中，Bi∈ R3×N,Bj∈ R3×N，根據(jù)剛體運(yùn)動(dòng)學(xué)的相關(guān)知識(shí),Bi和Bj的表達(dá)式亦可以簡(jiǎn)潔寫(xiě)成如下形式

彈簧的當(dāng)前長(zhǎng)度為

定義沿矢量h方向的單位矢量為?h∈R3×1，其相對(duì)與全局坐標(biāo)系中的投影為

彈簧當(dāng)前長(zhǎng)度的變化率為

設(shè)彈簧的剛度系數(shù)為k，阻尼器的阻尼系數(shù)為c，作動(dòng)器的主動(dòng)作用力為 fa.SDA的作用力始終沿Pi和Pj的連線，其大小為

物體I和J上所受到的作用力在全局坐標(biāo)系下的列陣為

由式(16)和式(17)知

另一方面，SDA的廣義力向量Qei∈RN×1和Qej∈RN×1與虛功有如下關(guān)系

將式(25)，式(26)，式(28)，式(29)代入式(27)并與式(30)比較可知，作用在物體I和物體J上的廣義力向量為

和有限元分析中的節(jié)點(diǎn)力向量的組裝類(lèi)似，對(duì)于含有SDA的多體系統(tǒng)，在靜平衡或動(dòng)力學(xué)分析中，按上述公式計(jì)算其對(duì)應(yīng)的廣義外力陣，累加到全局廣義坐標(biāo)陣對(duì)應(yīng)的位置.

3 彈簧--阻尼--作動(dòng)器系統(tǒng)的雅可比矩陣

3.1 雅可比矩陣之一：?Q/?q

第2節(jié)給出了動(dòng)力學(xué)方程中廣義外力陣和廣義坐標(biāo)之間的關(guān)系，本節(jié)推導(dǎo)廣義力向量對(duì)廣義坐標(biāo)的偏導(dǎo)數(shù)項(xiàng).方便起見(jiàn)，根據(jù)廣義坐標(biāo)所屬物體的不同，將此矩陣寫(xiě)成分塊形式，即

在總體雅可比矩陣中，?Qei/?qi和 ?Qej/?qj分布在主對(duì)角線附近，而 ?Qei/?qj和 ?Qej/?qi為耦合項(xiàng)，其分布在遠(yuǎn)離主對(duì)角線位置.為推導(dǎo)方便，首先進(jìn)行以下中間量的計(jì)算

綜合式(24)、式(33)～式(38)得

同理可得

結(jié)合式(31)和式(32)，對(duì)于主對(duì)角元項(xiàng)有

前述已經(jīng)提及Bi的表達(dá)式和qj無(wú)關(guān)聯(lián)，Bj和qi無(wú)關(guān)聯(lián).對(duì)于任意的V向量，V∈R3×1，有

故對(duì)于耦合項(xiàng)，有

綜合式(39)～式(46)，即得SDA對(duì)應(yīng)的廣義力向量對(duì)廣義坐標(biāo)的偏導(dǎo)數(shù)矩陣.注意到在式(39)～式(42)中涉及到Bi和Bj(或其轉(zhuǎn)置矩陣)與某一向量乘積后對(duì)廣義坐標(biāo)求導(dǎo)的矩陣，即

理論上這些項(xiàng)可以由矩陣Bi和Bj(每一個(gè)元素)對(duì)廣義坐標(biāo)qi，qj的偏導(dǎo)數(shù)構(gòu)成的三維矩陣和待乘向量復(fù)合而成，或?qū)懗蓮埩啃问剑@樣做并不利于程序的設(shè)計(jì).附錄中直接以四元數(shù)為例給出了其具體的表達(dá)式.對(duì)于其他形式的姿態(tài)描述方式，可直接根據(jù)Bi(Bj)矩陣的每一項(xiàng)經(jīng)過(guò)簡(jiǎn)單求導(dǎo)得出.

3.2 雅可比矩陣之二：?Q/?˙q

當(dāng)SDA中的阻尼器的阻尼不為0時(shí)，廣義力向量和廣義速度相關(guān).本節(jié)給出廣義力向量對(duì)廣義速度的偏導(dǎo)數(shù).在動(dòng)力分析中，此項(xiàng)對(duì)于雅可比矩陣的精確計(jì)算，亦必不可少.注意到

故有

根據(jù)式(24)，有

綜合式(49)～式(54)有

同理

可以看出廣義力向量對(duì)廣義速度的偏導(dǎo)數(shù)相比于廣義力向量對(duì)廣義坐標(biāo)的偏導(dǎo)數(shù)公式更為簡(jiǎn)潔，且當(dāng)寫(xiě)成整體形式時(shí)，此矩陣具有對(duì)稱(chēng)性.

3.3 考慮雅可比矩陣的計(jì)算流程

3.1和3 .2節(jié)中給出的SDA對(duì)應(yīng)附加的雅克比矩陣是完全精確的，未引入任何的近似.從相關(guān)公式可以看出其涉及到的計(jì)算量較大，且表達(dá)式復(fù)雜，呈現(xiàn)出高度的非線性.為了便于有序、高效地進(jìn)行數(shù)值計(jì)算，下面給出了SDA對(duì)應(yīng)的廣義力向量和雅可比矩陣的計(jì)算流程：

(4)根據(jù)式(24)計(jì)算內(nèi)力Fs,根據(jù)式(31),式(32)計(jì)算廣義外力陣Qei,Qej.

(8)根據(jù)式(41),式(42),式(45),式(46)計(jì)算

(9)根據(jù)式(55)～式(58),計(jì)算

4 數(shù)值算例

為了說(shuō)明上述雅可比矩陣的有效性，本節(jié)針對(duì)兩個(gè)算例，進(jìn)行多體動(dòng)力學(xué)計(jì)算，并著重從雅可比矩陣的構(gòu)成對(duì)多體動(dòng)力學(xué)分析的收斂性以及計(jì)算效率的角度進(jìn)行分析.多體系統(tǒng)中物體的姿態(tài)均采用四元數(shù)進(jìn)行描述，以避免歐拉角或卡爾丹角建?？赡茉斐傻钠娈悊?wèn)題[1].程序的運(yùn)行環(huán)境為MATLAB R2016a，處理器：Intel Core i5-4300U@1.90GHz、內(nèi)存4GB.

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

如圖2所示的曲柄滑塊機(jī)構(gòu)，1號(hào)連桿與地面，2號(hào)連桿與1號(hào)連桿，2號(hào)連桿與滑塊3之間分別通過(guò)旋轉(zhuǎn)鉸連接，滑塊3可沿地面水平滑動(dòng)，其上作用大小為1kN，水平向右的力.連桿1和2的之間連接有SDA，連接點(diǎn)在分別在兩連桿的質(zhì)心處，SDA的剛度為 k(單位:N/m)，阻尼為 c(單位:N·s/m)，作動(dòng)力為 fa(單位:N).連桿1和連桿2桿長(zhǎng)均為2m，質(zhì)量均為50kg，垂直于桿長(zhǎng)方向的轉(zhuǎn)動(dòng)慣量均為15 kg·m2，沿桿長(zhǎng)方向的轉(zhuǎn)動(dòng)慣量均為 0.02kg·m2.滑塊的質(zhì)量為1kg，繞慣性主軸方向的轉(zhuǎn)動(dòng)慣量均為1kg·m2，不計(jì)重力，無(wú)初始速度.

圖2 曲柄滑塊機(jī)構(gòu)Fig.2 Slider-crank mechanism

采用Newmark算法對(duì)上述算例進(jìn)行動(dòng)力學(xué)分析.時(shí)間步長(zhǎng)取?t=0.01s，仿真步數(shù)：100.Newmark參數(shù)取α=0.25，δ=0.5，在此種參數(shù)選擇下Newmark算法為二階精度且無(wú)能量耗散.每一時(shí)間步的非線性代數(shù)方程組采用牛頓--拉夫遜迭代算法計(jì)算，初始迭代值取為上一時(shí)間步末端的廣義加速度和拉氏乘子向量，收斂誤差按殘差的2范數(shù)定義，取ε=1×10?7.

圖3～圖5分別給出了在參數(shù)c=0，fa=0下變化彈簧剛度k、不計(jì)入和計(jì)入SDA的雅可比矩陣對(duì)應(yīng)的收斂曲線、以及兩種情況下的分析時(shí)間對(duì)比.可以看出，當(dāng)不計(jì)入SDA造成的附加雅可比矩陣項(xiàng)時(shí)，每一時(shí)間步需要的收斂步數(shù)隨著剛度k的增加而增加.相比之下，計(jì)入SDA造成的附加雅可比矩陣可以保證每一時(shí)間步更快的收斂，在當(dāng)前情況，各時(shí)間步的非線性方程組求解均只需要3步即可收斂到指定精度.

圖3 不計(jì)SDA的雅可比矩陣，每一時(shí)間步的迭代步數(shù)Fig.3 Iteration number of each time step excluding the Jacobian matrix of SDA

圖4 計(jì)入SDA的雅可比矩陣，每一時(shí)間步的迭代步數(shù)Fig.4 Iteration number of each time step including the Jacobian matrix of SDA

圖5 分析時(shí)間對(duì)比Fig.5 Comparison of analysis time

圖6 不計(jì)SDA的雅可比矩陣，每一時(shí)間步的迭代步數(shù)Fig.6 Iteration number of each time step excluding the Jacobian matrix of SDA

上述結(jié)果只考慮了SDA的參數(shù)c=0的情況.根據(jù)式(55)～式(58)，此種情況下SDA的廣義力向量對(duì)廣義速度的偏導(dǎo)數(shù)為零矩陣，相當(dāng)于只考察了SDA的廣義力向量對(duì)廣義坐標(biāo)的偏導(dǎo)數(shù)對(duì)收斂性的影響.為了考察廣義力向量對(duì)廣義速度的偏導(dǎo)數(shù)對(duì)收斂性的影響，選取固定的k=2×104，fa=0，c分別取 1×101，1×102，1×103，1×104等4種情況進(jìn)行動(dòng)力學(xué)計(jì)算.分別考慮在總體雅克比矩陣中不計(jì)和計(jì)入SDA的廣義力向量對(duì)廣義速度的偏導(dǎo)數(shù)對(duì)雅可比矩陣的附加項(xiàng)的情況.每一時(shí)間步的迭代步數(shù)和分析時(shí)間如圖6～圖8所示.從圖6可以發(fā)現(xiàn)當(dāng)不計(jì)入此附加修正項(xiàng)時(shí)，每一時(shí)間步需要的收斂步數(shù)隨著c的增加而增加.由圖7和圖8可知計(jì)入本文推導(dǎo)的精確雅克比修正項(xiàng)后，所需迭代步更少，節(jié)省計(jì)算時(shí)間，特別是對(duì)于SDA的阻尼較高的情況.采用控制變量法對(duì)作動(dòng)器的作動(dòng)力fa進(jìn)行分析，可以得到類(lèi)似的結(jié)論.另外當(dāng)嘗試進(jìn)行靜平衡分析時(shí)，不計(jì)入總體雅可比矩陣中SDA對(duì)應(yīng)的部分，則完全不能收斂.

圖7 計(jì)入SDA的雅可比矩陣，每一時(shí)間步的迭代步數(shù)Fig.7 Iteration number of each time step including the Jacobian matrix of SDA

圖8 分析時(shí)間對(duì)比Fig.8 Comparison of analysis time

4.2 Tensegrity的找形分析

Tensegrity是由多個(gè)不連續(xù)的壓桿和連續(xù)的拉索構(gòu)成的空間機(jī)構(gòu)，最早來(lái)自于藝術(shù)領(lǐng)域，后被應(yīng)用到航空航天、建筑結(jié)構(gòu)等工程領(lǐng)域[2324].找形分析[2530]是Tensegrity設(shè)計(jì)過(guò)程中的重要環(huán)節(jié).本算例來(lái)自文獻(xiàn)[31]，考察初始位形如圖9所示的由6根壓桿，21根拉索的兩層Tensegrity結(jié)構(gòu)，粗線表示壓桿，細(xì)線表示拉索.采用多體動(dòng)力學(xué)方法，將壓桿考慮為物體，拉索由SDA模擬.借鑒動(dòng)力松弛法[3233]的思想，在初始的不平衡力的作用下構(gòu)件發(fā)生運(yùn)動(dòng)，SDA的阻尼器耗散能量最終使系統(tǒng)趨于平衡狀態(tài)，從而完成找形分析.壓桿的質(zhì)量均為0.3kg，垂直于桿長(zhǎng)方向的轉(zhuǎn)動(dòng)慣量均為0.1kg·m2，沿桿長(zhǎng)方向的轉(zhuǎn)動(dòng)慣量均為0.01kg·m2.拉索的原始長(zhǎng)度和連接信息詳見(jiàn)文獻(xiàn)[31].

圖9 系統(tǒng)的初始位形Fig.9 Initial position of the system

通過(guò)此算例考察，當(dāng)系統(tǒng)中含有多個(gè)SDA時(shí)，雅可比矩陣對(duì)收斂的影響，不考慮作動(dòng)器的作用，即fa=0.所有彈簧的剛度k(單位:N/m)，以及阻尼c(單位:N·s/m)保持一致，無(wú)初始速度.不對(duì)系統(tǒng)做外在約束，以體現(xiàn)Tensegrity結(jié)構(gòu)“自平衡”的特點(diǎn).Newmark參數(shù)取α=0.3，δ=0.6，步長(zhǎng)取?t=0.01s，收斂誤差取ε=1×10?7，每一時(shí)間步的最大迭代步數(shù)設(shè)為1000.將SDA對(duì)應(yīng)的?Qe/?q和項(xiàng)分別簡(jiǎn)記為M2和M3，精確雅可比矩陣的其他部分記為M1.

表1給出了在若干SDA參數(shù)下，選取不同的部分參與雅可比矩陣組裝，仿真前100步的平均收斂步數(shù)和和計(jì)算時(shí)間.從中可以看出動(dòng)力學(xué)分析的收斂情況與參數(shù)的選取有著緊密的關(guān)系.對(duì)于此算例，采用完全精確的雅可比矩陣相比于只考慮部分項(xiàng)的雅可比矩陣，每一步所需要的迭代步更少.對(duì)于彈簧剛度k較小，阻尼c較大的情況，M3項(xiàng)對(duì)收斂速度的影響較大.當(dāng)彈簧的剛度k較大時(shí)，各種情況的收斂性均變差.對(duì)于第8種參數(shù)，甚至完全精確的雅可比矩陣也不收斂，這是因?yàn)?，此種情況下彈簧的剛度大且無(wú)阻尼，在較短的時(shí)間內(nèi)，系統(tǒng)的位形變化大，當(dāng)前的時(shí)程分析步長(zhǎng)?t=0.01s取得過(guò)大.對(duì)于第5種參數(shù)可以觀察到似乎反常的情況，只計(jì)入M1項(xiàng)和計(jì)入M1+M2或M1+M3項(xiàng)相比，平均迭代步變多而總體分析時(shí)間卻變少，這是因?yàn)橛?jì)算SDA的雅可比矩陣需要耗費(fèi)一定的時(shí)間，迭代步減少所節(jié)省的時(shí)間并不一定能夠以抵消矩陣計(jì)算所帶來(lái)的額外開(kāi)銷(xiāo)，但是對(duì)完全精確的雅可比矩陣(M1+M2+M3)相比于只計(jì)入M1項(xiàng)在收斂速度和時(shí)間上仍有明顯的優(yōu)勢(shì).對(duì)于含有較多的SDA的多體系統(tǒng)，精確計(jì)算其造成的雅可比矩陣是必要的.

表1 不同參數(shù)下的平均迭代步數(shù)與計(jì)算時(shí)間對(duì)比Table 1 The average iteration number and computing time by di ff erent parameters

圖10給出了彈簧剛度k=200N/m，阻尼系數(shù)c=10N·s/m時(shí)，系統(tǒng)在0～0.3s區(qū)間的能量變化曲線，其中殘余能量定義為：物體的總動(dòng)能+SDA的總彈性勢(shì)能，SDA耗散的能量即為阻尼器做的功，算法耗散的能量定義為：系統(tǒng)的初始能量–殘余能量–SDA耗散的能量.首先可以發(fā)現(xiàn)系統(tǒng)的殘余能量迅速被耗散到較低的水平并緩慢衰減,這說(shuō)明系統(tǒng)在較短的時(shí)間內(nèi)即運(yùn)動(dòng)到平衡狀態(tài)附近.此外可以觀察到Newmark算法自身耗散能量遠(yuǎn)小于SDA阻尼器做的功.隨著仿真時(shí)間的增加，各能量將趨于穩(wěn)定.圖11給出了在上述參數(shù)選取下，仿真至20s后的結(jié)果，以及文獻(xiàn)[31]的結(jié)果，可以看出，二者在形態(tài)上比較一致，從而側(cè)面說(shuō)明了此算例分析的正確性.進(jìn)一步的分析表明，在平衡構(gòu)型下，桿件全部處于受壓狀態(tài)，繩索均處于受拉狀態(tài)，其中上層的3根側(cè)邊繩索所受的拉力最大，約為4.92N，說(shuō)明此算例構(gòu)造的合理性.

圖10 系統(tǒng)的能量變化曲線Fig.10 System energy variation curve

圖11 找形的結(jié)果對(duì)比Fig.11 Comparison of the form- fi nding result

5 結(jié)論

采用隱式算法對(duì)多體系統(tǒng)動(dòng)力學(xué)微分--代數(shù)方程求解，需要計(jì)算系統(tǒng)的雅可比矩陣.對(duì)于含有SDA的多體系統(tǒng)，其造成的雅可比矩陣是廣義坐標(biāo)的復(fù)雜非線性函數(shù).本文以多體動(dòng)力學(xué)微分-代數(shù)方程的Newmark求解算法為例，說(shuō)明其雅可比矩陣的構(gòu)成形式，以及SDA對(duì)雅可比矩陣的具體的影響項(xiàng).詳細(xì)推導(dǎo)了多體系統(tǒng)動(dòng)力學(xué)中SDA的精確的雅可比矩陣，包括廣義力向量對(duì)廣義坐標(biāo)的偏導(dǎo)數(shù)矩陣和廣義力向量對(duì)廣義速度的偏導(dǎo)數(shù)矩陣.以姿態(tài)角的四元數(shù)描述為例，給出了其組成關(guān)鍵項(xiàng)的具體表達(dá)式.當(dāng)姿態(tài)角采用歐拉角或卡爾丹角等其他類(lèi)型描述時(shí)，只需計(jì)算出其對(duì)應(yīng)的關(guān)鍵項(xiàng)的表達(dá)式，進(jìn)行整體代入即可.在此基礎(chǔ)上，進(jìn)一步研究了SDA引起的雅可比矩陣對(duì)動(dòng)力學(xué)分析收斂性的影響，結(jié)果表明隨著SDA的彈簧剛度，阻尼器的阻尼以及作動(dòng)器的作用力的增大，SDA造成的附加雅可比矩陣對(duì)多體動(dòng)力學(xué)分析收斂性影響明顯.采用不考慮SDA造成的附加雅可比矩陣會(huì)動(dòng)力學(xué)分析導(dǎo)致收斂速度變慢，甚至不能收斂.利用本文推導(dǎo)的公式，計(jì)入SDA影響的精確的雅可比矩陣無(wú)論是從迭代步數(shù)還是從最終的仿真時(shí)間上來(lái)看，都能夠達(dá)到很好的分析效果.本文中考慮的SDA為線性彈簧和線性阻尼器，且SDA的連接物體均為剛體的情況，對(duì)于非線性彈簧或非線性阻尼器的SDA或SDA關(guān)聯(lián)物體為柔性體的情況可按照類(lèi)似的思路導(dǎo)出.

1洪嘉振.計(jì)算多體系統(tǒng)動(dòng)力學(xué).北京:高等教育出版社,1999(Hong Jiazhen,Computational Dynamics of Multibody Systems.Beijing:Higher Education Press,1999(in Chinese))

2石奇端,馬履中.六自由度并聯(lián)機(jī)構(gòu)組合彈簧阻尼減振裝置.農(nóng)業(yè)機(jī)械學(xué)報(bào),2007,38(8):128-131(Shi Qiduan,Ma Luzhong.Research on the device of 6-dof parallel mechanisms combined with elastic damping system.Transactions of the Chinese Society for Agricultural Machinery,2007,38(8):128-131(in Chinese))

3王其東,陳無(wú)畏,張炳力.基于多體模型的汽車(chē)半主動(dòng)懸架控制方法研究.機(jī)械工程學(xué)報(bào),2004,40(1):104-108(Wang Qidong,Chen Wuwei,Zhang Bingli.Study on the control method of automobile semi-active suspension based on the multibody model.Journal of Mechanical Engineering,2004,40(1):104-108(in Chinese))

4劉穎,王瑞.基于計(jì)算多體系統(tǒng)動(dòng)力學(xué)的機(jī)織卡爾丹角多體動(dòng)力學(xué)模型.機(jī)械工程學(xué)報(bào),2014,50(5):102-107(Liu Ying,Wang Rui.A cardan angular model of rigid bodies for woven fabrics based on the dynamics of multibody systems.Journal of Mechanical Engineering,2014,50(5):102-107(in Chinese))

5譚蔚，賈占斌，聶清德.彈簧阻尼器在塔器防振中的應(yīng)用.化學(xué)工程,2013,41(12):16-19,34(Tan Wei,Jia Zhanbin,Nie Qingde.Application of spring damper in antivibration of columns.Chemical Engineering,2013,41(12):16-19,34(in Chinese))

6王庚祥，劉宏昭.多體系統(tǒng)動(dòng)力學(xué)中關(guān)節(jié)效應(yīng)模型的研究進(jìn)展.力學(xué)學(xué)報(bào)，2015,47(1):31-50(Wang Gengxiang,Liu Hongzhao.Research progress of joint e ff ects model in multibody system dynamics.Chinese Journal of Theoretical and Applied Mechanics,2015,47(1):31-50(in Chinese))

7王國(guó)平.多體系統(tǒng)動(dòng)力學(xué)數(shù)值解法.計(jì)算機(jī)仿真,2006,23(12):86-89(Wang Guoping.Numerical algorithms of multibody system dynamics.Computer Simulation,2006,23(12):86-89(in Chinese))

8 Newmark N.A method of computation for structural dynamics.ASCE Journal of the Engineering Mechanics Division,1959,85:67-94

9 Negrut D,Rampalli R,Ottarsson G,et al.On an implementation of the Hilber-Hughes-Taylor method in the context of index 3 di ff erential-algebraic equations of multibody dynamics.Journal of Computational and Nonlinear Dynamics,2007,2(1):73-85

10馬秀騰，翟彥博，謝守勇.多體系統(tǒng)指標(biāo)2運(yùn)動(dòng)方程HHT方法違約校正.力學(xué)學(xué)報(bào),2017,49(1):175-181(Ma Xiuteng,Zhai Yanbo,Xie Shouyong.HHT method with constraints violation correction in theindex2equationsofmotionformultibodysystem.ChineseJournal of Theoretical and Applied Mechanics,2017,49(1):175-181(in Chinese))

11 Chun J,Hulbert G.A time integration algorithm for structural dynamics with improved numerical dissipation:the generalized-α method.ASME Journal of Applied Mechanics,1993,60(2):371-375

12 Bathe KJ,Baig MMI.On a composite implicit time integration procedure for nonlinear dynamics.Computers and Structures,2005,83(31-32):2513-2524

13 Bathe KJ.Conserving energy and momentum in nonlinear dynamics:A simple implicit time integration scheme.Computers and Structures,2007,85(7-8):437-445

14鐘萬(wàn)勰,高強(qiáng).約束動(dòng)力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分．動(dòng)力學(xué)與控制,2006,4(3):193-200(Zhong Wanxie,Gao Qiang.Integration of constrained dynamical system via analytical structural mechanics.JournalofDynamicsandControl,2006,4(3):193-200(inChinese))

15吳鋒,高強(qiáng),鐘萬(wàn)勰.基于祖沖之類(lèi)方法的多體動(dòng)力學(xué)方程保能量保約束積分.計(jì)算機(jī)輔助工程,2014,23(1):64-68(Wu Feng,Gao Qiang,Zhong Wanxie.Energy and constraint preservation integration for multibody equations based on Zu Chongzhi method.Computer Aided Engineering,2014,23(1):64-68(in Chinese))

16 Shabana AA.Dynamics of Multibody Systems(3rd edition).New York:Cambridge University Press,2005

17 Bayo E,Garc′?a de Jal′o J.Kinematics and Dynamics Simulation of Multibody System,the Real-time Challenge.New York:Springer,1994

18張雄,王天舒.計(jì)算動(dòng)力學(xué).北京:清華大學(xué)出版社,2007(Zhang Xiong,Wang Tianshu.Computational Dynamics.Beijing:Tsinghua University Press,2007(in Chinese))

19劉延柱,潘振寬,戈新生.多體系統(tǒng)動(dòng)力學(xué)(第2版).北京:高等教育出版社,2014(LiuYanzhu,PanZhenkuan,GeXinsheng.Dynamics of Multibody Systems(2nd edition).Beijing:Higher Education Press,2014(in Chinese))

20 Hussein BA,Shabana AA.Sparse matrix implicit numerical integration of the sti ffdi ff erential/algebraic equations:implementation.Nonlinear Dynamics,2011,65:369-382

21田強(qiáng).基于絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的柔性多體系統(tǒng)動(dòng)力學(xué)研究與應(yīng)用.[博士論文].武漢:華中科技大學(xué)，2009(Tian Qiang.Research and application of fl exible multibody system dynamics based on the absolute nodal coordinate method.[PhD Thesis].Wuhan:HuazhongUniversityofScienceandTechnology,2009(inChinese))

22 Bottasso CL,Dopico D,Trainelli L.On the optimal scaling of index three DAEs in multibody dynamics.Multibody System Dynamics,2007,19:3-20

23 Motro R.Tensegrity systems:the state of the art.International Journal of Space Structures,1992,7:75-82

24 Skelton R,de Oliveka M.Tensegrity Systems. New York：Springer，2009

25 Pagitz M,Mirats Tur JM.Finite element based form- fi nding algorithm for tensegrity structures.International Journal of Solids and Structures,2009,46:3235-3240

26 Gasparini D,Klinka KK,Arcaro VF.A fi nite element for formfi nding and static analysis of tensegrity structures.Journal of Mechanics of Materials&Structures,2011,6(9-10):1239-1254

27 Koohestani K.Form- fi nding of tensegrity structures via genetic algorithm.International Journal of Solids and Structures,2012,49:739-747

28 Tran HC,Lee J.Form- fi nding of tensegrity structures using double singular value decomposition.Engineering with Computers,2013,29:71-86

29 Koohestani K,Guest S.A new approach to the analytical and numerical form fi nding of tensegrity structures.International Journal of Solids and Structures,2013,50(19):2995-3007

30 Zhang LY,Li Y,Cao YP,et al.Sti ff ness matrix based form- fi nding method of tensegrity structures.Engineering Structures,2014,58:36-48

31 Massimo C,Josep M,Mirats T.A comprehensive dynamic model for class-1 tensegrity systems based on quaternions.International Journal of Solids and Structures,2011,48:785-802

32 Barnes MR.Form fi nding and analysis of tension structures by dynamic relaxation.International Journal of Space Structures,1999,14(2):89-104

33 Ali NBH,Rhode-Barbarigos L,Smith IFC.Analysis of clustered tensegrity structures using a modi fi ed dynamic relaxation algorithm.InternationalJournalofSolidsandStructures,2011,48(5):637-647

RIGID BODY SYSTEM DYNAMIC WITH THE ACCURATE JACOBIAN MATRIX OF SPRING-DAMPER-ACTUATOR1)

Kan Ziyun Peng Haijun2)Chen Biaoshong
(State Key Laboratory of Structural Analysis for Industrial Equipment，Department of Engineering Mechanics，Dalian University of Technology，Dalian 116024，China)

The spring-damper-actuator(SDA)is a common force element in multibody system and widely used in the fi eld of engineering.The governing equations of multibody dynamic system established by absolute coordinate methods are di ff erential-algebraic equations which are usually nonlinear and complex.To ensure the stability and accuracy of the numerical solutions,the implicit algorithms are commonly used to solve the dynamic equations.While the calculations of Jacobian matrices are the crucial process in implicit algorithms.For a multibody system containing the SDA,the additional Jacobian matrices induced by the SDA are highly nonlinear functions of the generalized coordinates and generalized velocities.A lot of current research works focus on the calculation of generalized force vector,however the calculations of additional Jacobian matrices are less concerned.This paper focuses on dynamic analysis of multi-rigid-body systems containing the SDA.Firstly,the construction of the accurate Jacobian matrices in solving the dynamic equations is investigated based on the Newmark algorithm.Then,the additional Jacobian matrices relating to the generalized force vector of the SDA are analytically derived.These matrices consist of the partial derivative of generalized force vector with respect to the generalized coordinates and the generalized velocities.Finally,the in fl uence of additional Jacobian matrices on the convergence of dynamic analysis is investigated via two numerical examples.The numerical results indicate that when the values of sti ff ness,damping and active force are large,the additional Jacobian matrices induced by the SDA have a signi fi cant in fl uence on the convergence of dynamic analysis.When the additional Jacobian matrices induced by the SDA are taken into account,the dynamic analysis can achieve convergence with less iteration steps and the computational time thus can be reduced.

spring-damper-actuator,Jacobian matrix,multibody systems,implicit algorithm,tensegrity

O313.7

A

10.6052/0459-1879-17-030

2017–01–23收稿，2017–05–23 錄用,2017–05–24 網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金(11472069,11772074,91648204)和國(guó)家重點(diǎn)研發(fā)計(jì)劃(2016YFB0200702)資助項(xiàng)目.

2)彭海軍，副教授，主要研究方向:動(dòng)力學(xué)與控制.E-mail:hjpeng@dlut.edu.cn

闞子云，彭海軍，陳飆松.考慮彈簧阻尼作動(dòng)器解析雅可比矩陣的多剛體動(dòng)力學(xué)分析.力學(xué)學(xué)報(bào)，2017,49(5):1103-1114

Kan Ziyun,Peng Haijun,Chen Biaoshong.Rigid body system dynamic with the accurate Jacobian matrix of spring-damper-actuator.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1103-1114

附錄A 本文中約定的矩陣的求導(dǎo)運(yùn)算法則

附錄B 關(guān)鍵矩陣的表達(dá)式

設(shè)物體局部坐標(biāo)系原點(diǎn)在全局坐標(biāo)系中的坐標(biāo)為 R=?x y z?T，四元數(shù)描述下的姿態(tài)坐標(biāo)為Θ=[l0l1l2l3]T.該.SDA與該物體的連接點(diǎn)在該物體局部坐標(biāo)系的矢量為：[v1v2···vn]T為任意n維列向量.下面給出SDA對(duì)應(yīng)雅可比矩陣的關(guān)鍵項(xiàng).(未給出指標(biāo)索引的項(xiàng)均為0)