計算多體系統動力學中引入多體系統離散時間傳遞矩陣法的混合算法研究

高浩鵬，黃映云，趙建華，孫宇鵬

(海軍工程大學 船舶與動力學院，武漢 430033)

隨著計算機技術特別是計算機應用機技術的發展，人們對于復雜機械系統的多體動力學建模與求解成為可能。目前以多體系統為研究對象的方法主要有兩種:需要系統總體動力學方程的傳統多體系統動力學方法［1］和無需系統總體動力學方程的多體系統傳遞矩陣法［2］。這兩種方法都是對機械系統中的機構進行研究，雖然都從經典力學發展而來，但這在建模和求解時都具有自己的特點。

建模是多體系統動力學分析的第一步，建模的合理性與準確性直接影響計算結果，鑒于此本文擬結合計算多體系統動力學與多體系統離散時間傳遞矩陣法，形成一種混合算法。該算法主要將兩種算法進行互補，一定程度上擴展了研究對象的內容，也進一步提高了建模精度。

1 研究背景

對于復雜機械系統，若構件間的拓撲關系較為明確，在三維實體建模的基礎上，使用計算多體動力學建模分析較為方便，人為因素引起的建模差異較小，易于求解［3］。隨著計算機技術，特別是計算機應用技術的迅速發展，人們的研究對象也趨向于復雜機械系統［4］;然而復雜機械系統存在有大量非線性因素、多學科耦合，特別是存在有大量的非確定性因素，這樣使得以物理實體建模為基礎的多體動力學分析就比較困難。

柴油機結構復雜，整個工作過程涉及多部件的大空間運動，包括平動、轉動以及平動和轉動的合成運動，其適合通過計算多體動力學建模及分析;然而整個柴油機結構中存在有大量的螺栓連接結構，而大量研究結果表明螺栓連接結構復雜，影響其動力學特性因素較多，很難準確建立其模型進行特性分析［5］。柴油機結構中有很多如缸蓋螺栓連接這種受力較大的結合面，如果簡單地使用剛性連接必然會對整個柴油機特性產生影響，所以在通過多體動力學對柴油機進行建模時，螺栓連接結構就成了建模的難點。

為了解決多體動力學中螺栓結合面建模問題，須引入人為控制的數值計算方法，然而數值計算方法的引入需以結合面的特性參數為依據。對于如螺栓這種復雜結合面的動力學特性，很難通過建立其實體模型得到;為了較為準確得到其動力學特性就必須結合試驗的方法?，F代測試技術的發展，使得通過試驗得到螺栓結合面動力學特性參數成為可能。故本文的目的就是研究一種結合實體建模和數學建模的方法解決上述問題，該方法要易于通過目前應用軟件實現。本文提出計算多體系統動力學與多體系統離散時間傳遞矩陣法混合算法來解決上述問題。

2 算法研究

原則上講，計算多體動力學和多體系統傳遞矩陣法能單獨進行大多數機械系統動力學分析;然而這兩個方法的建模原理以及求解算法有一定差異，導致各自方法具有自己的特點。為了充分發揮不同方法在相關方面的優勢，便研究其混合算法，該算法要易于在目前較為流行的多體動力學分析軟件中實現。下面通過分別對這兩種算法進行論述的基礎上，引出本文所提出的混合算法。

2.1 計算多體動力學方法

計算多體動力學由經典力學發展而來，目前發展到剛柔混合多體系統動力學［1］。根據其建模方法的不同主要分為兩種方法:拉格朗日法和笛卡爾法。多體動力學笛卡爾數學建模方法起源于機械領域，由于該建模方法有利于計算機數值求解，目前在應用軟件中多采用該方法建模［3］。

下面主要對多剛體系統動力學建模及求解理論進行概述。通過應用軟件在建立多體系統模型時，每個構件具有如式(1)所列舉的15個變量和如式(2)所列舉的15個方程。

式(1)和式(2)中各參數的含義見附錄，由式(1)可知其變量為平動速度、笛卡爾位置坐標、轉動慣量、角速度及歐拉角坐標。根據牛頓第二定律、位移與速度的微分關系、帶拉格朗日乘子的拉格朗日第一類方程以及轉動慣量定律等可以得到如式(2)所示的微分－代數方程。機械系統中每個剛體構件就是通過式(1)建立動力學方程，在此基礎上通過構件間的約束力構成構件間作用關系，再建立系統的邊界條件，便構成了整個系統的動力學方程。建模與求解過程如圖1所示。該方法的優點是:基于物理實體建模，整個模型精度高，易于軟件實現可視化界面;剛柔混合多體動力學的發展對柔性體特性分析較易;整個建模過程人為控制因素少，分布式建模方便。缺點是:整個實體建模工作量大;建立的微分方程存在高階次，求解速度較慢，求解過程可能由于方程的病態而不能計算。

圖1 計算多體動力學建模與求解流程圖Fig.1 The modeling and solving flow chart of computational dynamics of multi－body systems

2.2 多體系統離散時間傳遞矩陣法

多體系統傳遞矩陣法是用傳遞矩陣研究多體系統動力學的方法，其研究對象是以各種方式相聯結的多個構件組成的多體系統，涉及兵器、航空航天、機械等多種行業。根據所研究對象的特性不同其主要有線性多體系統傳遞矩陣法和多體系統離散時間傳遞矩陣法［6］。線性多體系統傳遞矩陣法解決線性多體系統振動問題［7］。多體系統離散時間傳遞矩陣法主要解決時變、非線性、大運動、受控等一般多體系統動力學問題［8］。

雙眼矯正視力差值與遠立體視呈負相關性（r=-0.673，P＜0.001），與RDS測得的立體視銳度、交叉立體視、非交叉立體視均呈負相關性（r=-0.439、-0.466、-0.451，P＜0.001）。

下面主要論述多體系統離散時間傳遞矩陣法的求解過程。多體系統離散時間傳遞矩陣法的理論較為簡潔，即通過傳遞矩陣關聯系統輸入狀態與輸出狀態矢量，如式(3)所示;其中Z為狀態矢量，U為傳遞矩陣。傳遞矩陣由具體結構抽象得到多個力學特征元件并按規則拼裝而成，故其重點和難點主要集中在對系統傳遞矩陣的研究。狀態矢量是一個單列陣，表示研究系統中一節點的力學狀態;該矢量一般由該點的位置坐標、角坐標、內力以及內力矩組成。

本文研究重點為算法，故對傳遞矩陣的研究不作為重點，整個求解過程如圖2所示。該方法的優點是:無需建立系統的總體動力學方程，易于結合試驗數據進行數學建模;由傳遞矩陣建立的方程階次較低，計算量小，計算速度快，計算過程不易產生病態。缺點是對研究系統結構進行抽象，傳遞矩陣的建立具有一定難度，對數學知識要求高，人為因素影響較大;在動力學分析時很難對于柔性體構件的動力學特性進行分析，而且對于構件自身不同部位的特性分析較為困難。

圖2 多體系統離散時間

2.3 混合算法

多體動力學和多體系統離散時間傳遞矩陣法的研究對象都為多體系統，在對這兩種方法優缺點分析的基礎上引出混合算法，其特點主要是:在多體動力學分析軟件基礎上，引入多體系統離散時間傳遞矩陣法，實現可視化的建模，易于對系統進一步優化及特性分析，并且建模過程中人為因素較小，基于實體模型準確度高;多體系統離散時間傳遞矩陣法的引入，使得復雜子系統特別是具有不確定因素的子系統建模趨向于數學建模，再通過試驗確定傳遞矩陣的參數，其模型精度較高，并且建模和求解過程無需系統總體動力學方程。

在多體動力學軟件中引入多體系統離散時間傳遞矩陣法的建模及求解過程如圖3所示?；旌纤惴ㄊ窃诟髯宰酉到y建模的基礎上研究的，故其研究的重點和求解的可行性集中在兩個子系統建模與求解的結合處。對于多體子系統其一般將力/力矩作為邊界條件，其求解一般容易得到上文中論述的15個變量;多體系統離散時間傳遞矩陣法一般將位置和力/力矩作為狀態矢量。故混合算法是將多體子系統求解得到的位置量作為多體系統離散時間傳遞矩陣法的輸入與輸出狀態矢量的已知部分，再通過求解多體系統離散時間傳遞矩陣法所建立的代數方程得到狀態矢量中的內力和內力矩，再將內力和內力矩作為多體系統與多體系統離散時間傳遞矩陣子系統連接處的力/力矩邊界條件添加到多體系統中進行求解得到位置量，這樣循環計算便使得計算得以延續，這點便是混合算法的核心。在多體動力學中，這種混合算法可以通過用戶子程序進行實現［9］。

圖3 混合算法建模與求解流程圖Fig.3 The modeling and solving flow chart of hybrid algorithm

3 算例

螺栓連接結構大量存在于柴油機結構中，影響螺栓結合面特性的因素較多，而且很多因素具有不確定性，故本文以螺栓連接兩個剛體的結構為例來使用混合算法對其建模與求解。

3.1 建模

本文在多體動力學分析軟件中建立的混合模型如圖4所示。多體子系統1中剛體受到向上的外力，多體子系統2中剛體與大地固定，文中將螺栓結合面簡化為4個多體系統離散時間傳遞矩陣法中彈簧－阻尼器典型元件，分別位于多體子系統1的四個角上，并對稱分布。

圖4 混合模型Fig.4 The hybrid model

對每個子系統進行孤立的分析:多體子系統1的邊界條件為外力以及多體系統離散時間傳遞矩陣子系統對其作用力;多體子系統2的邊界條件為地面約束力以及多體系統離散時間子系統對其作用力;多體系統離散時間傳遞矩陣子系統主要由傳遞矩陣連接其輸入和輸出狀態矢量。由參考文獻［8］可得到單向彈簧－阻尼器的傳遞方程，拓展后得到三向彈簧－阻尼器的傳遞方程如式(4)所示，其狀態矢量為 Z=［x，y，z，qx，qy，qz，1］T，下標 OUT 表示輸出狀態矢量、IN 表示輸入狀態矢量，傳遞矩陣中各參數的含義見附錄。本文通過Fox－Euler法進行線性化，故傳遞矩陣中的參數已知或者為迭代過程中上個時刻的函數，其傳遞矩陣中相關線性化系數如式(5)所示。

根據文中提出的混合算法，由多體子系統得到多體系統離散時間傳遞矩陣子系統中狀態矢量的位置量，然后通過傳遞方程求解得到多體系統離散時間傳遞矩陣子系統中狀態矢量的內力和內力矩量，故可以將式(4)變形為式(6)，以便于線性方程組求解。所以混合算法建立模型求解的迭代算法為:求解多體子系統1和多體子系統2，分別得到式(6)中增廣矩陣系數中的參數 xIN、xOUT、yIN、、yOUT、zIN、zOUT以及線性化系數中上一時刻的位置和速度量，再加上已知的剛度和阻尼系數，即已知式(6)中所有增廣矩陣系數所涉及的參數;將這些參數代入式(6)便可以求解得到多體系統離散時間傳遞矩陣子系統中輸入和輸出的狀態矢量、和，之后將這幾個量作為多體子系統1和多體子系統2的已知邊界條件并進行求解，依次循環計算直到設定時刻。

整個建模過程中，多體子系統通過可視化實體建模實現，建立了微分－代數方程組;多體系統離散時間傳遞矩陣子系統通過用戶自定義程序實現，并嵌入到多體動力學分析軟件中，進而建立線性代數方程組?；旌纤惴ㄖ杏脩舳x子程序的編寫是難點，其主函數的格式為“VOID_FUNCTION VFOSUB(int*id，REAL*time，REAL*par，int*nPar，BOOL*dflag，BOOL*iflag，REAL*value)”，其中所涉及的主要問題有:用戶定義參數包括剛度、阻尼系數、求解步長以及子系統之間連接點的標識，這些參數通過par變量實現從多體動力學到用戶定義子程序的傳遞;由式(6)可知其增廣矩陣形式較為簡單，但為了方法的適用性以及增強混合算法的可執行性，本文選取全選主元高斯消去法求解式(6)的線性方程組;式(6)增廣矩陣系數中涉及到的上一時刻的位置和速度量通過靜態內存變量的定義來實現。

3.2 結果分析

模型建立之后，在多體動力學分析軟件中進行求解，其中求解時間1 s、求解步長0.001 s。多體子系統施加的外力如圖5所示，其幅值為1 000 N。當4個彈簧－阻尼器的剛度和阻尼系數都為500 N/mm和0.1 N·s/mm時，4個連接點處狀態矢量的力和位移值相等，分別如圖6中實線和虛線所示。由圖6可知，聯結點處力狀態矢量的幅值為261 N，大于系統施加外力幅值的四分之一，并且其有一定的振蕩衰減，這都是由于剛度和阻尼的影響。

圖5 外力邊界條件時間歷程曲線Fig.5 The boundary condition of external force

圖6 力和位移狀態矢量時間歷程曲線Fig.6 The force and displacement state vector curves

螺栓連接結合面實際結構復雜，四個螺栓的等效剛度和阻尼很難相等(如由于預緊力大小不一致，或者螺栓結構上細小區別等所導致)。文中假設其中一個剛度小于其它三個剛度20%，則該彈簧－阻尼器處的力狀態矢量如圖7中實線所示，圖中虛線為相鄰彈簧－阻尼器的力狀態矢量曲線;圖8為四個聯結點處加速度曲線。由結果分析可知，剛度較小的彈簧－阻尼器聯結點處其力狀態矢量幅值較相鄰聯結點處小，并且其大小與對角聯結點處相等，驗證了實際在擰緊螺栓時按對角順序擰緊，也驗證該方法的正確性;由圖8知，雖然只有一個彈簧－阻尼器的剛度有變化，但四個聯結點處的加速度在幅值和相位上都有一定的影響。

圖7 力狀態矢量對比曲線Fig.7 The contrast curves of force state vector

圖8 加速度對比曲線Fig.8 The contrast curves of acceleration

混合算法的研究主要是解決一些難于通過實體建模來模擬其結構的子系統，而該子系統對整個動力學特性有一定的影響，該算法的優點是繼承了計算多體動力學和多體系統離散時間傳遞矩陣法各自的優點:在可視化實體建模的基礎上，可以得到整個動力學分析模型較為精確的邊界條件;難于實體建模的子系統可以通過多體系統離散時間傳遞矩陣法實現，只需建立該子系統的傳遞方程即可。文中重點對算法進行研究，故所列舉算例較為簡單，選取了彈簧－阻尼典型元件作為混合算法研究對象;當然，在多體動力學中該典型元件也可按照彈性連接的方式添加，但其與混合算法有質的區別:多體系統離散時間傳遞矩陣法的引入是建立代數方程組實現，而彈性連接是建立微分－代數方程組實現。混合算法并不僅僅局限于將彈簧－阻尼元件引入到計算多體動力學中，其可以將所有多體系統離散時間傳遞矩陣法中典型元件以及典型元件的組合結構引入到計算多體動力學中，即混合算法的實現拓展了多體系統動力學的研究對象。

4 結論

本文以計算多體動力學和多體系統離散時間傳遞矩陣法建模及求解理論為依據，提出了在多體動力學中引入多體系統離散時間傳遞矩陣法的混合算法，并以柴油機結構中大量存在的螺栓結合面為算例通過混合算法進行建模與求解。文中算例較為簡單，目的主要是論證混合算法的可行性與正確性;在基于實際結構分析時，傳遞矩陣中的參數應基于試驗數據，以保證建模的準確性。

混合算法在多體動力學分析軟件中實現，其特點是:結合物理實體建模與數學建模，提高了模型的準確性;多體系統離散時間傳遞矩陣法的引入結合試驗可以建立一些難以通過物理實體模擬的結構，一定程度上擴展了研究對象，并且結合試驗數據提高了分析結果的可靠性;混合算法無需建立系統總體動力學方程并且便于分析傳遞矩陣中特性參數對整個系統動力學的靈敏度，便于結構優化，一定程度上減小了實體建模的數據量。

［1］洪嘉振.計算多體系統動力學［M］.北京:高等教育出版社，1999.

［2］芮筱亭，贠來峰，陸毓琪，等.多體系統傳遞矩陣法及其應用［M］.北京:科學出版社，2008.

［3］陳立平，張云清，任衛群，等.機械系統動力學分析及ADAMS應用教程［M］.北京:清華大學出版社，2006.

［4］熊曉燕.復雜機械系統動態特性分析和試驗辨識方法的研究［D］.太原:太原理工大學，2008.

［5］毛寬民，黃小磊，田紅亮，等.機床固定結合面參數識別及其擬合方法［J］.華中科技大學學報(自然科學版)，2011，39(3):18－21.MAO Kuan-min，HUANG Xiao-lei，TIAN Hong-liang，et al.Identifying dynamic parameters of fixed joints in machine tools and its fitting method［J］.Huazhong Univ.of Sci.＆Tech.(Natural Science Edition)，2011，39(3):18 －21.

［6］Rui X T，He B.Advanced in discrete time transfer matrix method formultibody system dynamics［M］.IUTAM Symposium onMultiscaleProblemsin MultibodySystem Contacts，Stuttgart，Germany，2006.

［7］芮筱亭，王國平，贠來峰，等.某遠程多管火箭振動特性研究［J］.振動與沖擊，2005，24(1):9 －12.RUI Xiao-ting，WANG Guo-ping，YUN Lai-feng，et al.Study on the vivration characteristics of LRMLRS［J］.Journal of Vibration and Shock，2005，24(1):9 －12.

［8］楊富鋒.基于多體系統傳遞矩陣法的工程系統的動力學及其應用［D］.南京:南京理工大學，2006.

［9］邢俊文，陶永忠，譯.MSC.Software.MSC.ADAMS/View 高級培訓教程［M］.北京:清華大學出版社，2004.

附 錄

符 號

V平動速度矢量 QR平動廣義力

Vx，Vy，Vz平動速度矢量的三個分量 CR平動約束反力

R位移矢量 Qr轉動廣義力

x，y，z位移矢量的三個分量 Cr轉動約束反力

Pγ廣義動量矢量 B坐標變換矩陣

Pψ，Pθ，Pφ廣義動量矢量的三個分量 z狀態矢量

ωe角速度矢量 q狀態矢量中力分量

ωψ，ωθ，ωφ角速度矢量的三個分量 K 剛度

γ角位移矢量 C阻尼

ψ，θ，φ 角位移矢量的三個分量 Cdi、DO、DI線性化系數

M質量矩陣