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

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

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

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

1 研究背景

對于復(fù)雜機(jī)械系統(tǒng)，若構(gòu)件間的拓?fù)潢P(guān)系較為明確，在三維實體建模的基礎(chǔ)上，使用計算多體動力學(xué)建模分析較為方便，人為因素引起的建模差異較小，易于求解［3］。隨著計算機(jī)技術(shù)，特別是計算機(jī)應(yīng)用技術(shù)的迅速發(fā)展，人們的研究對象也趨向于復(fù)雜機(jī)械系統(tǒng)［4］;然而復(fù)雜機(jī)械系統(tǒng)存在有大量非線性因素、多學(xué)科耦合，特別是存在有大量的非確定性因素，這樣使得以物理實體建模為基礎(chǔ)的多體動力學(xué)分析就比較困難。

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

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

2 算法研究

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

2.1 計算多體動力學(xué)方法

計算多體動力學(xué)由經(jīng)典力學(xué)發(fā)展而來，目前發(fā)展到剛?cè)峄旌隙囿w系統(tǒng)動力學(xué)［1］。根據(jù)其建模方法的不同主要分為兩種方法:拉格朗日法和笛卡爾法。多體動力學(xué)笛卡爾數(shù)學(xué)建模方法起源于機(jī)械領(lǐng)域，由于該建模方法有利于計算機(jī)數(shù)值求解，目前在應(yīng)用軟件中多采用該方法建模［3］。

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

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

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

2.2 多體系統(tǒng)離散時間傳遞矩陣法

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

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

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

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

圖2 多體系統(tǒng)離散時間

2.3 混合算法

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

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

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

3 算例

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

3.1 建模

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

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

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

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

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

3.2 結(jié)果分析

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

［3］陳立平，張云清，任衛(wèi)群，等.機(jī)械系統(tǒng)動力學(xué)分析及ADAMS應(yīng)用教程［M］.北京:清華大學(xué)出版社，2006.

［4］熊曉燕.復(fù)雜機(jī)械系統(tǒng)動態(tài)特性分析和試驗辨識方法的研究［D］.太原:太原理工大學(xué)，2008.

［5］毛寬民，黃小磊，田紅亮，等.機(jī)床固定結(jié)合面參數(shù)識別及其擬合方法［J］.華中科技大學(xué)學(xué)報(自然科學(xué)版)，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］芮筱亭，王國平，贠來峰，等.某遠(yuǎn)程多管火箭振動特性研究［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］楊富鋒.基于多體系統(tǒng)傳遞矩陣法的工程系統(tǒng)的動力學(xué)及其應(yīng)用［D］.南京:南京理工大學(xué)，2006.

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

附 錄

符 號

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

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

R位移矢量 Qr轉(zhuǎn)動廣義力

x，y，z位移矢量的三個分量 Cr轉(zhuǎn)動約束反力

Pγ廣義動量矢量 B坐標(biāo)變換矩陣

Pψ，Pθ，Pφ廣義動量矢量的三個分量 z狀態(tài)矢量

ωe角速度矢量 q狀態(tài)矢量中力分量

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

γ角位移矢量 C阻尼

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

M質(zhì)量矩陣