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

闞子云 彭海軍 陳飆松

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

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

闞子云 彭海軍2)陳飆松

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

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

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

引言

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

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

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

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

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

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

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

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

其中

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

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

如圖1所示,物體I和物體J通過SDA連接，連接點分別為Pi和Pj.當(dāng)前時刻兩物體相對于全局坐標(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)系的原點相對于全局坐標(biāo)系的矢量在全局坐標(biāo)系的投影.Θi和Θj分別表示各自物體的局部坐標(biāo)系相對于全局坐標(biāo)系的姿態(tài)角.此種建模方式即為多體動力學(xué)中最為經(jīng)典和工程中普遍使用的笛卡爾絕對建模方式.當(dāng)采用歐拉角或卡爾丹角描述姿態(tài)角Θi∈R3×1,Θj∈R3×1，一個物體共有6個廣義坐標(biāo).當(dāng)采用四元數(shù)描述姿態(tài)角時，Θi∈ R4×1，Θj∈ R4×1，一個物體共有7個廣義坐標(biāo).值得說明的是：由于四元數(shù)之間并非相互獨立，此種情況下動力學(xué)方程中需要添加相應(yīng)的約束方程.為統(tǒng)一起見，以下將單個物體對應(yīng)的廣義坐標(biāo)個數(shù)記為N.

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

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

令

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

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

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

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

設(shè)彈簧的剛度系數(shù)為k，阻尼器的阻尼系數(shù)為c，作動器的主動作用力為 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é)點力向量的組裝類似，對于含有SDA的多體系統(tǒng)，在靜平衡或動力學(xué)分析中，按上述公式計算其對應(yīng)的廣義外力陣，累加到全局廣義坐標(biāo)陣對應(yīng)的位置.

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

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

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

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

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

同理可得

結(jié)合式(31)和式(32)，對于主對角元項有

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

故對于耦合項，有

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

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

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

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

故有

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

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

同理

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

3.3 考慮雅可比矩陣的計算流程

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

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

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

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

4 數(shù)值算例

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

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

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

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

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

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

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

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

圖5 分析時間對比Fig.5 Comparison of analysis time

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

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

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

圖8 分析時間對比Fig.8 Comparison of analysis time

4.2 Tensegrity的找形分析

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

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

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

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

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

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

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

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

5 結(jié)論

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

1洪嘉振.計算多體系統(tǒ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é)報,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王其東,陳無畏,張炳力.基于多體模型的汽車半主動懸架控制方法研究.機(jī)械工程學(xué)報,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劉穎,王瑞.基于計算多體系統(tǒng)動力學(xué)的機(jī)織卡爾丹角多體動力學(xué)模型.機(jī)械工程學(xué)報,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)動力學(xué)中關(guān)節(jié)效應(yīng)模型的研究進(jìn)展.力學(xué)學(xué)報，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王國平.多體系統(tǒng)動力學(xué)數(shù)值解法.計算機(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運動方程HHT方法違約校正.力學(xué)學(xué)報,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鐘萬勰,高強(qiáng).約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分．動力學(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),鐘萬勰.基于祖沖之類方法的多體動力學(xué)方程保能量保約束積分.計算機(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張雄,王天舒.計算動力學(xué).北京:清華大學(xué)出版社,2007(Zhang Xiong,Wang Tianshu.Computational Dynamics.Beijing:Tsinghua University Press,2007(in Chinese))

19劉延柱,潘振寬,戈新生.多體系統(tǒ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).基于絕對節(jié)點坐標(biāo)方法的柔性多體系統(tǒ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)國家自然科學(xué)基金(11472069,11772074,91648204)和國家重點研發(fā)計劃(2016YFB0200702)資助項目.

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

闞子云，彭海軍，陳飆松.考慮彈簧阻尼作動器解析雅可比矩陣的多剛體動力學(xué)分析.力學(xué)學(xué)報，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)運算法則

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

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