闞子云 彭海軍 陳飆松

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

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

闞子云 彭海軍2)陳飆松

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

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

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

引言

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

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

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

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

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

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

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

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

其中

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

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

如圖1所示,物體I和物體J通過SDA連接，連接點分別為Pi和Pj.當前時刻兩物體相對于全局坐標系的位形坐標(廣義坐標)為

廣義速度

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

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

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

將式(14)對時間求導，根據數(shù)學關系知

令

其中，Bi∈ R3×N,Bj∈ R3×N，根據剛體運動學的相關知識,Bi和Bj的表達式亦可以簡潔寫成如下形式

彈簧的當前長度為

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

彈簧當前長度的變化率為

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

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

由式(16)和式(17)知

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

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

和有限元分析中的節(jié)點力向量的組裝類似，對于含有SDA的多體系統(tǒng)，在靜平衡或動力學分析中，按上述公式計算其對應的廣義外力陣，累加到全局廣義坐標陣對應的位置.

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

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

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

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

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

同理可得

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

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

故對于耦合項，有

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

理論上這些項可以由矩陣Bi和Bj(每一個元素)對廣義坐標qi，qj的偏導數(shù)構成的三維矩陣和待乘向量復合而成，或寫成張量形式，但這樣做并不利于程序的設計.附錄中直接以四元數(shù)為例給出了其具體的表達式.對于其他形式的姿態(tài)描述方式，可直接根據Bi(Bj)矩陣的每一項經過簡單求導得出.

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

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

故有

根據式(24)，有

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

同理

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

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

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

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

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

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

4 數(shù)值算例

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

4.1 曲柄滑塊機構

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

圖2 曲柄滑塊機構Fig.2 Slider-crank mechanism

采用Newmark算法對上述算例進行動力學分析.時間步長取?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的雅可比矩陣對應的收斂曲線、以及兩種情況下的分析時間對比.可以看出，當不計入SDA造成的附加雅可比矩陣項時，每一時間步需要的收斂步數(shù)隨著剛度k的增加而增加.相比之下，計入SDA造成的附加雅可比矩陣可以保證每一時間步更快的收斂，在當前情況，各時間步的非線性方程組求解均只需要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

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

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

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

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

表1給出了在若干SDA參數(shù)下，選取不同的部分參與雅可比矩陣組裝，仿真前100步的平均收斂步數(shù)和和計算時間.從中可以看出動力學分析的收斂情況與參數(shù)的選取有著緊密的關系.對于此算例，采用完全精確的雅可比矩陣相比于只考慮部分項的雅可比矩陣，每一步所需要的迭代步更少.對于彈簧剛度k較小，阻尼c較大的情況，M3項對收斂速度的影響較大.當彈簧的剛度k較大時，各種情況的收斂性均變差.對于第8種參數(shù)，甚至完全精確的雅可比矩陣也不收斂，這是因為，此種情況下彈簧的剛度大且無阻尼，在較短的時間內，系統(tǒ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)在較短的時間內即運動到平衡狀態(tài)附近.此外可以觀察到Newmark算法自身耗散能量遠小于SDA阻尼器做的功.隨著仿真時間的增加，各能量將趨于穩(wěn)定.圖11給出了在上述參數(shù)選取下，仿真至20s后的結果，以及文獻[31]的結果，可以看出，二者在形態(tài)上比較一致，從而側面說明了此算例分析的正確性.進一步的分析表明，在平衡構型下，桿件全部處于受壓狀態(tài)，繩索均處于受拉狀態(tài)，其中上層的3根側邊繩索所受的拉力最大，約為4.92N，說明此算例構造的合理性.

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

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

5 結論

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

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

2石奇端,馬履中.六自由度并聯(lián)機構組合彈簧阻尼減振裝置.農業(yè)機械學報,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王其東,陳無畏,張炳力.基于多體模型的汽車半主動懸架控制方法研究.機械工程學報,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)動力學的機織卡爾丹角多體動力學模型.機械工程學報,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譚蔚，賈占斌，聶清德.彈簧阻尼器在塔器防振中的應用.化學工程,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)動力學中關節(jié)效應模型的研究進展.力學學報，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)動力學數(shù)值解法.計算機仿真,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)指標2運動方程HHT方法違約校正.力學學報,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鐘萬勰,高強.約束動力系統(tǒng)的分析結構力學積分．動力學與控制,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吳鋒,高強,鐘萬勰.基于祖沖之類方法的多體動力學方程保能量保約束積分.計算機輔助工程,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張雄,王天舒.計算動力學.北京:清華大學出版社,2007(Zhang Xiong,Wang Tianshu.Computational Dynamics.Beijing:Tsinghua University Press,2007(in Chinese))

19劉延柱,潘振寬,戈新生.多體系統(tǒng)動力學(第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田強.基于絕對節(jié)點坐標方法的柔性多體系統(tǒng)動力學研究與應用.[博士論文].武漢:華中科技大學，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 網絡版發(fā)表.

1)國家自然科學基金(11472069,11772074,91648204)和國家重點研發(fā)計劃(2016YFB0200702)資助項目.

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

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

附錄B 關鍵矩陣的表達式

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