動(dòng)力學(xué)建模方法與解法總結(jié)

1、目錄1 剛體系統(tǒng) 12 彈性系統(tǒng)動(dòng)力學(xué) 63 高速旋轉(zhuǎn)體動(dòng)力學(xué) 101 剛體系統(tǒng)一般力學(xué)研究的對(duì)象， 是由兩個(gè)或兩個(gè)以上剛體通過(guò)鉸鏈等約束聯(lián)系在一起 的力學(xué)系統(tǒng)， 為一般力學(xué)研究對(duì)象。 自行車(chē)、萬(wàn)向支架陀螺儀通?？煽闯啥鄤傮w 系統(tǒng)。人體在某種意義上也可簡(jiǎn)化為一個(gè)多剛體系統(tǒng)。現(xiàn)代航天器、機(jī)器人、人 體和仿生學(xué)中關(guān)于動(dòng)物運(yùn)動(dòng)規(guī)律的研究都提出了多剛體系統(tǒng)的一系列理論模型 作為研究對(duì)象。 多剛體系統(tǒng)按其內(nèi)部聯(lián)系的拓?fù)浣Y(jié)構(gòu)， 分為樹(shù)型和非樹(shù)型 (包含 有閉鏈)；按其同外界的聯(lián)系情況，則有有根和無(wú)根之別。利用圖論的工具可以 一般地分析多剛體系統(tǒng)的構(gòu)造， 建立系統(tǒng)的數(shù)學(xué)模型和動(dòng)力學(xué)方程組。 也可從分 析力學(xué)

2、中的高斯原理出發(fā)， 用求極值的優(yōu)化算法直接求解系統(tǒng)的運(yùn)動(dòng)和鉸鏈反力。 依照多剛體系統(tǒng)動(dòng)力學(xué)的理論和方法， 廣泛采用電子計(jì)算機(jī)對(duì)這些模型進(jìn)行研究， 對(duì)于精確地掌握這些對(duì)象的運(yùn)動(dòng)規(guī)律是很有價(jià)值的。1.1 自由物體的變分運(yùn)動(dòng)方程任意一個(gè)剛體構(gòu)件i，質(zhì)量為mi，對(duì)質(zhì)心的極轉(zhuǎn)動(dòng)慣量為Ji，設(shè)作用于剛體的所有外力向質(zhì)心簡(jiǎn)化后得到外力矢量Fi和力矩ni，若定義剛體連體坐標(biāo)系 xo y 的原點(diǎn) o 位于剛體質(zhì)心，則可根據(jù)牛頓定理導(dǎo)出該剛體帶質(zhì)心坐 標(biāo)的變分運(yùn)動(dòng)方程：riTmiri Fi iJi i ni 0(1-1)其中， ri 為固定于剛體質(zhì)心的連體坐標(biāo)系原點(diǎn) o 的代數(shù)矢量， i 為連體坐標(biāo) 系相對(duì)于全局

3、坐標(biāo)系的轉(zhuǎn)角，r與i分別為ri與i的變分。定義廣義坐標(biāo)：qi riT, iT(1-2)廣義：Qi FiT,ni T(1-3)及質(zhì)量矩陣：體坐標(biāo)系原點(diǎn)固定于剛體質(zhì)心時(shí)用廣義力表示的剛體變分運(yùn)動(dòng)方程：(1-5)qiT(M iqi Qi) 01.2 束多體系統(tǒng)的運(yùn)動(dòng)方程考慮由nb個(gè)構(gòu)件組成的機(jī)械系統(tǒng)，對(duì)每個(gè)構(gòu)件運(yùn)用式(1-5)，組合后可得到系統(tǒng)的變分運(yùn)動(dòng)方程為：nbi1qiTM iqiQi 0(1-6)若組合所有構(gòu)件的廣義坐標(biāo)矢量、質(zhì)量矩陣及廣義力矢量，構(gòu)造系統(tǒng) 的廣義坐標(biāo)矢量、質(zhì)量矩陣及廣義力矢量為：T T T Tq q1 ,q2,.,qnb(1-7)Mdiag (M 1, M 2 ,., M n

4、b )(1-8)Q Q1T,Q2T,.,QnTbT(1-9)系統(tǒng)的變分運(yùn)動(dòng)方程則可緊湊地寫(xiě)為：qTMq Q 0 (1-10)對(duì)于單個(gè)構(gòu)件，運(yùn)動(dòng)方程中的廣義力同時(shí)包含作用力和約束力，但在 一個(gè)系統(tǒng)中，若只考慮理想運(yùn)動(dòng)副約束，根據(jù)牛頓第三定律，可知作用在 系統(tǒng)所有構(gòu)件上的約束力總虛功為零，若將作用于系統(tǒng)的廣義外力表示為：QA Q1AT,Q2AT,.,QnAbTT(1-11)其中：QiA FiAT,nAT， i 1,2,., nb(1-12)則理想約束情況下的系統(tǒng)變分運(yùn)動(dòng)方程為：TAqTMq QA 0 (1-13)式中虛位移q 與作用在系統(tǒng)上的約束是一致的。系統(tǒng)運(yùn)動(dòng)學(xué)約束和驅(qū)動(dòng)約束的組合如式 (1-

5、10)，為：(q,t) 0(1-14)對(duì)其微分得到其變分形式為：q q 0(1-15)3式(1-13)和(1-15) 組成受約束的機(jī)械系統(tǒng)的變分運(yùn)動(dòng)方程。 為導(dǎo)出約束機(jī)械系統(tǒng)變分運(yùn)動(dòng)方程易于應(yīng)用的形式，運(yùn)用拉格朗日乘子定 理對(duì)式 (1-13) 和(1-15) 進(jìn)行處理。拉格朗日乘子定理：設(shè)矢量 b Rn，矢量x Rn，矩陣A Rmn為常數(shù) 矩陣，如果有：bT x 0(1-16)對(duì)于所有滿(mǎn)足式(1-84 )的x條件都成立。Ax 0(1-17)則存在滿(mǎn)足式( 1-85)的拉格朗日乘子矢量Rm 。bT x TAx 0(1-18)其中x為任意的。在式(1-13)和(1-15)中，q Rn，M Rnn，

6、QA Rn， qR” n，運(yùn)用拉格朗日乘子定理于式 (1-13) 和 (1-15) ，則存在拉格朗日乘子矢量Rm ，對(duì)于任意的 q 應(yīng)滿(mǎn)足：Mq QAT q T q q MqqTQAT q 0(1-19)由此得到運(yùn)動(dòng)方程的拉格朗日乘子形式：TAMqqTQA(1-20)式(1-20) 還必須滿(mǎn)足式 (1-10) 、(1-12) 和 (1-14) 表示的位置約束方程、速度約束方程及加速度約束方程，如下：(q,t)0(1-21)(q,q,t) q(q,t)q0，t(q,t)(1-22)(q, q, q, t)q(q,t)q (q,q,t)0，( qq) qq 2 qtqtt (1-23)以上三式其維

7、數(shù)同式 (1-14)式(1-20) 、(1-21) 、(1-22) 和(1 -23) 組成約束機(jī)械系統(tǒng)的完整的運(yùn)動(dòng)方程。 將式 (1-20) 與(1-23) 聯(lián)立表示為矩陣形式：QA(1-24)式(1-24)即為多體系統(tǒng)動(dòng)力學(xué)中最重要的動(dòng)力學(xué)運(yùn)動(dòng)方程，式(1-24)還必須滿(mǎn)足式(1-22)和(1-23)。它是一個(gè)微分一一代數(shù)方程組，不同于單純的常微 分方程組問(wèn)題，其求解關(guān)鍵在于避免積分過(guò)程中的違約現(xiàn)象，此外，還要 注意DAE問(wèn)題的剛性問(wèn)題。如果系統(tǒng)質(zhì)量矩陣是正定的，并且約束獨(dú)立，那么運(yùn)動(dòng)方程就有唯一解。實(shí)際中的系統(tǒng)質(zhì)量矩陣通常是正定的，只要保證約束是獨(dú)立的，運(yùn)動(dòng) 方程就會(huì)有解。在實(shí)際數(shù)值迭代求

8、解過(guò)程中，需要給定初始條件，包括位置初始條件q(to)和速度初始條件q(to)。此時(shí)，如果要使運(yùn)動(dòng)方程有解，還需要滿(mǎn)足初 值相容條件，也就是要使位置初始條件滿(mǎn)足位置約束方程，速度初始條件 滿(mǎn)足速度約束方程。對(duì)于由式(1-24)及(1-21)、(1-22)確定的系統(tǒng)動(dòng)力學(xué)方程,初值相容條件為：(q(to),to)0(1-25)(q(to),q(to),to)q (q(to),to)q(to)(q(to),t。) 0(1-26)1.3正向動(dòng)力學(xué)分析、逆向動(dòng)力學(xué)分析與靜平衡分析對(duì)于一個(gè)確定的約束多體系統(tǒng)，其動(dòng)力學(xué)分析不同于運(yùn)動(dòng)學(xué)分析，并不需要系統(tǒng)約束方程的維數(shù)m等于系統(tǒng)廣義坐標(biāo)的維數(shù)n , m n。

9、在給定外力的作用下，從初始的位置和速度，求解滿(mǎn)足位置約束式(1-22)及速度約束式(1-23)的運(yùn)動(dòng)方程式(1-24)，就可得到系統(tǒng)的加速度和相應(yīng)的速度、位 置響應(yīng)，以及代表約束反力的拉格朗日乘子，這種已知外力求運(yùn)動(dòng)及約束 反力的動(dòng)力學(xué)分析，稱(chēng)為正向動(dòng)力學(xué)分析。如果約束多體系統(tǒng)約束方程的維數(shù)m與系統(tǒng)廣義坐標(biāo)的維數(shù)n相等，m n，也就是對(duì)系統(tǒng)施加與系統(tǒng)自由度相等的驅(qū)動(dòng)約束，那么該系統(tǒng)在運(yùn)動(dòng)學(xué)上就被完全確定，由2.2.3節(jié)的約束方程、速度方程和加速度方程可求解系統(tǒng)運(yùn)動(dòng)。在此情況下，雅可比矩陣是非奇異方陣，即：展開(kāi)式 (1-24) 的運(yùn)動(dòng)方程，為：(1-28)TA Mq q Q(1-29)由式(1-

10、29)可解得q，再由式(1-28)可求得 ，拉格朗日乘子 就唯一地確定 了作用在系統(tǒng)上的約束力和力矩(主要存在于運(yùn)動(dòng)副中) 。這種由確定的運(yùn) 動(dòng)求系統(tǒng)約束反力的動(dòng)力學(xué)分析就是逆向動(dòng)力學(xué)分析。q0如果一個(gè)系統(tǒng)在外力作用下保持靜止?fàn)顟B(tài)，也就是說(shuō)，如果：(1-30)那么，就說(shuō)該系統(tǒng)處于平衡狀態(tài)。將式(1 -30)代入運(yùn)動(dòng)方程式 (1-20)，得到5AQA由平衡方程式 (1-21) 及約束方程式(1-13) 可求出狀態(tài) q 和拉格朗日乘子平衡方程：(1-31)。這種求系統(tǒng)的平衡狀態(tài)及在平衡狀態(tài)下的約束反力的動(dòng)力學(xué)分析稱(chēng)為 靜)平衡分析。1.4 約束反力對(duì)于約束機(jī)械系統(tǒng)中的構(gòu)件 i ，設(shè)其與系統(tǒng)中某構(gòu)件

11、 j 存在運(yùn)動(dòng)學(xué)約束或驅(qū)動(dòng)約束，約束編號(hào)為 k 。除連體坐標(biāo)系 xo y 外，再在構(gòu)件 i 上以某點(diǎn) P為原點(diǎn)建立一個(gè)新的固定于構(gòu)件上的坐標(biāo)系 x Py ，稱(chēng)為運(yùn)動(dòng)副坐標(biāo)系，設(shè) 從坐標(biāo)系 x Py 到坐標(biāo)系 x oy 的變換矩陣為 Ci ，從坐標(biāo)系 xo y 到坐標(biāo)系xoy的變換矩陣為 A，則可導(dǎo)出由約束 k產(chǎn)生的反作用力和力矩分別為:TiFi(siPTT AT k T k i Ai riBiTkT rikT kT)(1-32)(1-33)以上兩式中，k為約束k對(duì)應(yīng)的拉格朗日乘子，反作用力F k和力矩Ti k均為運(yùn)動(dòng)副坐標(biāo)系 x Py 中的量。2彈性系統(tǒng)動(dòng)力學(xué)由于工業(yè)機(jī)器人、機(jī)械手、彈性聯(lián)動(dòng)裝

12、置、帶柔性附件人造衛(wèi)星、直升飛機(jī) 的旋翼等工程結(jié)構(gòu)發(fā)展的需求，使運(yùn)動(dòng)中的彈性結(jié)構(gòu)的動(dòng)力學(xué)分析得到了很大 的進(jìn)展。運(yùn)動(dòng)彈性體的動(dòng)力學(xué)分析屬于多體系統(tǒng)動(dòng)力學(xué)的范疇。而導(dǎo)出其有限元 格式的動(dòng)力學(xué)方程并研究其數(shù)值解法則是計(jì)算多體系統(tǒng)動(dòng)力學(xué)的任務(wù)。由于彈性變形與剛體運(yùn)動(dòng)的耦合導(dǎo)致了運(yùn)動(dòng)彈性體的動(dòng)力學(xué)方程為時(shí)變的或非線(xiàn)性的，因此運(yùn)動(dòng)中的彈性體會(huì)出現(xiàn)諸多非線(xiàn)性效應(yīng)。運(yùn)動(dòng)中彈性體的動(dòng)力分析問(wèn)題可分為兩類(lèi)，其一是具有給定剛體運(yùn)動(dòng)的彈性體的動(dòng)力分析，這類(lèi)問(wèn)題僅討論彈性體的剛體運(yùn)動(dòng)對(duì)其彈性變形的影響，比如機(jī)械手的彈性終端桿的振動(dòng)分析一般可歸于此類(lèi)。第二類(lèi)問(wèn)題是多體系統(tǒng)中之剛體運(yùn)動(dòng)與其中的彈性體的彈性變形的相互耦合的

13、動(dòng)力分析，在這類(lèi)問(wèn)題中，彈性體的變形會(huì)受到系統(tǒng)剛體運(yùn)動(dòng)的影響，反之彈性體的變形也會(huì)影響系統(tǒng)的剛體運(yùn)動(dòng)。下面采用運(yùn)動(dòng)參考系方法并用 Jourdain動(dòng)力學(xué)普遍方程導(dǎo)出了具有空間一 般運(yùn)動(dòng)的彈性體之通用的有限元?jiǎng)恿W(xué)方程， 其最大的優(yōu)點(diǎn)在于推導(dǎo)簡(jiǎn)單并適用 于各類(lèi)結(jié)構(gòu)及各種單元形式。對(duì)系統(tǒng)的動(dòng)力學(xué)方程的數(shù)值求解，一般可以采用直 接積分法。下面給出了對(duì)時(shí)變的運(yùn)動(dòng)彈性的動(dòng)力學(xué)方程的Neumann級(jí)數(shù)2直接積分解法，該方法可以在保證計(jì)算精度的前提下很大程度地節(jié)省機(jī)時(shí)。圖2-1所示為一運(yùn)動(dòng)的彈性體B，選用兩個(gè)坐標(biāo)系來(lái)定義彈性體B的剛體運(yùn)動(dòng)與彈性變形：靜系一oxXoX3，簡(jiǎn)記o系；原點(diǎn)在B上的01點(diǎn),固連于B

14、上的動(dòng)系一Oix0iXoix0i ,簡(jiǎn)記為01系。B的剛體移動(dòng)由01點(diǎn)對(duì)于o點(diǎn)的矢量訂，定義B的 空間轉(zhuǎn)動(dòng)則用0i系對(duì)o系的轉(zhuǎn)動(dòng)來(lái)定義，而B(niǎo)內(nèi)任意點(diǎn)P的彈性變形則用在Oi系 內(nèi)的彈性變形位移矢量u來(lái)表示。由圖可見(jiàn)B發(fā)生彈性變形后，其上任意一點(diǎn)P對(duì)o系的位置矢量可以表示ru(2-1)9(2-2)其中r是B未產(chǎn)生彈性變形時(shí)P點(diǎn)在Oi系中的位置矢量，u則表示P點(diǎn)的彈性變 形位移矢量。把（2-2）式代入（2-1）式并向o系投影，且采用矩陣形式表示為：r；r。：A"1 r。1u。1（ 2-3）其中r；和ro：分別表示r；和打向o系的投影列陣；A001表示01系向o系轉(zhuǎn)移的方向余弦矩陣。把（3-

15、3）式中u01的用有限元的格式，表達(dá)為：u01N（2-4）其中N為單元形函數(shù)矩陣，為P點(diǎn)所在單元的有限元結(jié)點(diǎn)位移列陣。把(2-4)式代入（2-3）式，并利用公式：?00 0000.A 11 A 1（2-5）其中001是:系相對(duì)于0系轉(zhuǎn)動(dòng)角速度在0系上投影的斜對(duì)稱(chēng)陣1由（2-3）式對(duì)時(shí)間分別求一次導(dǎo)數(shù)和二次導(dǎo)數(shù)可得 P點(diǎn)的速度v；和加速度a；，進(jìn)而可得到P點(diǎn)的虛速度 v； 于是P點(diǎn)鄰域之微元體的Jourdain動(dòng)力學(xué)普遍方程可以寫(xiě)作：v； 丁 A001 fmpdv aP0（2-6）其中：m；為彈性體在P點(diǎn)的質(zhì)量密度；f是作用于P點(diǎn)微元體上的全部力在。1系上的投影 對(duì)于 vpo T Aoo1 f

16、可利用常規(guī)有限元的格式將它寫(xiě)作 :?T2-7)vpo T Aoo1 f N T F K C其中：K和C分別為單元?jiǎng)偠汝嚭蛦卧枇﹃囋?P點(diǎn)的值；F為作用在P點(diǎn)微元體上的外力在Oi系的列陣，把求得的P點(diǎn)的虛速度和加速度以及(2-7)式代 ?入(2-6)式，并考慮到中諸元素之獨(dú)立性，可得P點(diǎn)微元體的動(dòng)力學(xué)方程為：mpdv Vpo apo 02-8)將(2-8) 式對(duì)單元積分便可得運(yùn)動(dòng)的彈性體的單元?jiǎng)恿W(xué)方程?MeCeKeFe2-9)式中 ：MeTmp N N dvCeNTN dv2 mp NT Aoo1oo1oo1A 1 N dvCsCdKeBTD B dvmp N TAoo1 Too1oo1 o

17、o1Aoo1N dvKsKdFeF dv?Too1 T om p N T Aoo1roo1oo1dvT oo T mp NAoooo1?oo1oo1oo1Aoo1 ro1 dvFsFd其中 Cs，KsFs 分別是常規(guī)有限元法中的單元阻力陣、剛度陣和外力向量 ,而 Cd ，Kd ，F(xiàn)d 則分別是由于剛體運(yùn)動(dòng)與彈性變形的耦合而產(chǎn)生的附加單元?jiǎng)恿ψ枘彡?、?dòng)力剛度陣和動(dòng)力力向量。 而且由于它們的表達(dá)式中含有表示彈性?體空間運(yùn)動(dòng)量 roo1 和 ， 因此 ,通常這些動(dòng)力附加項(xiàng)是時(shí)變的。當(dāng)彈性體的剛 體運(yùn)動(dòng)速度特別是轉(zhuǎn)動(dòng)速度較大時(shí) , 彈性體受到較大的慣性力作用 , 會(huì)產(chǎn)生變形 的耦合效應(yīng)。例如轉(zhuǎn)動(dòng)的梁,

18、 由于離心慣性力產(chǎn)生的軸向拉力會(huì)增大梁的抗彎剛 度, 即所謂的 “剛化效應(yīng) ”。這時(shí)在 (2-10) 式中的常規(guī)剛度陣 K s 中需計(jì)入結(jié)構(gòu) 的幾何剛度陣 , 關(guān)于各類(lèi)單元的幾何剛度陣可參閱有關(guān)非線(xiàn)性有限元的書(shū)籍。 而 結(jié)構(gòu)的幾何剛度陣往往是未知內(nèi)力的函數(shù) , 這時(shí)方程 (2-9) 式就是一個(gè)非線(xiàn)性的 動(dòng)力方程。但對(duì)于簡(jiǎn)單的彈性體 , 如梁 , 由于剛體運(yùn)動(dòng)的慣性力產(chǎn)生的軸力容易 求得 , 這時(shí)的幾何剛度陣就變?yōu)闀r(shí)變陣。 本文只討論幾何剛度陣為時(shí)變陣的情況 , 即方程 (2-9)式為時(shí)變動(dòng)力學(xué)方程時(shí)的數(shù)值解法。顯然, 若彈性體沒(méi)有剛體運(yùn)動(dòng) , 則方程 (2-9)式退化為常規(guī)的有限單元?jiǎng)恿W(xué) 方

19、程。把(2-9)式按常規(guī)有限元的組集方法進(jìn)行組集 , 便可得到對(duì)于運(yùn)動(dòng)彈性體的 具有時(shí)變特性的、通用的有限元?jiǎng)恿W(xué)方程 :? ?M C K F113高速旋轉(zhuǎn)體動(dòng)力學(xué)咼速旋轉(zhuǎn)體通常是由是由二個(gè)剛體卜環(huán)、內(nèi)環(huán)、轉(zhuǎn)子互相約束在一起而成，可使陀螺儀轉(zhuǎn)子具有空間轉(zhuǎn)動(dòng)的三個(gè)自由度。 過(guò)去曾長(zhǎng)期認(rèn)為，高速自轉(zhuǎn)的平衡 對(duì)稱(chēng)卡登陀螺儀和單剛體陀螺儀的理論模型沒(méi)有本質(zhì)區(qū)別，具有所謂定軸性。但實(shí)際上，理論研究和精密的實(shí)驗(yàn)研究都已證明這個(gè)想法是錯(cuò)誤的。平衡對(duì)稱(chēng)卡 登陀螺儀的空間定向大都具有里雅普諾夫意義下的不穩(wěn)定性（見(jiàn)運(yùn)動(dòng)穩(wěn)定性）?？ǖ峭勇輧x和單剛體陀螺儀模型有本質(zhì)區(qū)別，只有通過(guò)多剛體系統(tǒng)模型的研究才 能正確解釋卡登

20、陀螺儀的動(dòng)力學(xué)特征。圖3-1如圖3-1所示，對(duì)于外徑D與長(zhǎng)度I的比值D/I 5的轉(zhuǎn)子，如多缸內(nèi)燃機(jī)的 曲軸、機(jī)床主軸等，這些轉(zhuǎn)子的不平衡質(zhì)點(diǎn)不是集中在同一平面內(nèi)， 而是分布在 垂直于軸線(xiàn)的各個(gè)平面內(nèi)。對(duì)于這種轉(zhuǎn)子動(dòng)平衡問(wèn)題，一般都采用矢量法來(lái)求校 正質(zhì)量mb、mb的重徑積mbh和mbg。但是這種方法所帶來(lái)問(wèn)題是力多邊形不 易求解以及圖解法不夠精確。假如采用平面解法，不僅簡(jiǎn)單正確，而且對(duì)于沒(méi)有 動(dòng)平衡機(jī)的工廠(chǎng)無(wú)疑有一定的實(shí)用價(jià)值。上述轉(zhuǎn)子質(zhì)量分布簡(jiǎn)圖如圖3-2所示，不平衡質(zhì)量 葉、m2、m3分別分布在與回轉(zhuǎn)軸線(xiàn)垂直的三個(gè)平面1、2、3內(nèi)，各質(zhì)點(diǎn)距回轉(zhuǎn)軸線(xiàn)的矢徑分別為 幾、r2、r3。當(dāng)轉(zhuǎn)子以等角速度?；剞D(zhuǎn)時(shí)，各質(zhì)點(diǎn)所產(chǎn)生的離心慣性力分別為(3-1)(3-2)P2m2ri(3-3)P3m3h2圖3-2方向如圖所示。若選擇轉(zhuǎn)子左、右二端面T （過(guò)點(diǎn)A與軸線(xiàn)垂直的平面）、T （過(guò) 點(diǎn)B與軸線(xiàn)垂直的平面）作為校正平面，在T、T平面內(nèi)分別加上校正質(zhì)量mb、mb，矢徑為h、rb ，則校正質(zhì)量所產(chǎn)生的離心慣性力為Pb mbb 2和Pbmbrb2，Pi、P2、P3、Pb 和Pb組成了空間力系選取三坐標(biāo)軸x、y、z軸如圖所示，并將作用在轉(zhuǎn)子上的所有力向 YAZ平面和XAY平面投影，