敬 謙，劉宏昭，王庚祥

(1.西安理工大學(xué) 機(jī)械與精密儀器工程學(xué)院，西安 710048；2.隴東學(xué)院 機(jī)械工程學(xué)院，甘肅 慶陽(yáng) 745000)

機(jī)構(gòu)系統(tǒng)關(guān)節(jié)反力是約束反力的一種，關(guān)節(jié)反力計(jì)算的準(zhǔn)確性對(duì)機(jī)構(gòu)動(dòng)力學(xué)分析、工作特性分析、結(jié)構(gòu)設(shè)計(jì)、零部件選取以及機(jī)器工作壽命預(yù)測(cè)等意義重大。近年來(lái)，機(jī)構(gòu)系統(tǒng)運(yùn)動(dòng)中如何引入拉格朗日乘子并明確其物理意義仍然是國(guó)內(nèi)外研究的熱點(diǎn)[1]。關(guān)于拉格朗日乘子與關(guān)節(jié)反力的關(guān)系問(wèn)題，丁光濤[2]從理論分析的角度討論了完整約束和非完整約束力學(xué)系統(tǒng)中引入待定乘子兩種不同的途徑，著重研究了變分原理?xiàng)l件極值中引入待定乘子修正系統(tǒng)的拉格朗日函數(shù)的方式，給出了拉格朗日乘子與理想約束反力之間的關(guān)系表達(dá)式。

關(guān)節(jié)反力常見(jiàn)的計(jì)算方法有牛頓歐拉方程、達(dá)朗貝爾原理以及拉格朗日方程。趙燕等[3]利用牛頓歐拉法列出所有構(gòu)件的力和力矩平衡方程，確定了驅(qū)動(dòng)力和平臺(tái)受到的約束反力，但同時(shí)也指出當(dāng)機(jī)構(gòu)結(jié)構(gòu)復(fù)雜或需要求解的主動(dòng)力和運(yùn)動(dòng)副反力的數(shù)目更多時(shí)，這種方法很不方便。Marek 等[4-5]討論了包含有冗余約束和奇異約束的剛體系統(tǒng)約束反力的求解方法，給出了在包含非獨(dú)立約束條件下能夠確定單一關(guān)節(jié)反力或一組關(guān)節(jié)反力的前提條件。通過(guò)整車系統(tǒng)動(dòng)力學(xué)分析,劉剛同樣選用牛頓歐拉方程對(duì)工程機(jī)構(gòu)鉸接系統(tǒng)鉸點(diǎn)動(dòng)態(tài)約束反力進(jìn)行了研究[6]。文獻(xiàn)[7]建立了多對(duì)均布傳力桿空間RCCR機(jī)構(gòu)模型，利用達(dá)朗貝爾原理計(jì)算出了在動(dòng)力條件下機(jī)構(gòu)慣性力引起約束反力的理論表達(dá)式。拉格朗日增廣法將加速度約束方程與動(dòng)力學(xué)方程耦合進(jìn)行積分，求解時(shí)必須考慮速度和位移約束的違約修正,如果修正不當(dāng),往往會(huì)使得求解發(fā)散。為了避免違約修正，原亮明等[8]將完整約束方程進(jìn)行泰勒展開(kāi)后與動(dòng)力學(xué)方程聯(lián)立求解并以一個(gè)七桿機(jī)構(gòu)為例，驗(yàn)證了算法的正確性。彭慧蓮等[9]以多體系統(tǒng)中的樹(shù)形結(jié)構(gòu)為例，詳細(xì)介紹了選用局部坐標(biāo)建立約束方程的方法和優(yōu)點(diǎn)，指出對(duì)具有光滑固定面約束和光滑柱鉸鏈約束的平面多剛體系統(tǒng)的關(guān)節(jié)約束反力與Lagrange 乘子一一對(duì)應(yīng)，并以曲柄滑塊機(jī)構(gòu)為例利用矢量方程法加以驗(yàn)證。上述文獻(xiàn)對(duì)拉格朗日增廣法中待定乘子的引入原理進(jìn)行了詳細(xì)介紹，并從理論分析的角度詳述了機(jī)構(gòu)系統(tǒng)在受完整約束或非完整約束下關(guān)節(jié)反力的數(shù)值計(jì)算方法。另外，當(dāng)前利用增廣法求解機(jī)構(gòu)系統(tǒng)動(dòng)力學(xué)的目的主要集中在對(duì)機(jī)構(gòu)運(yùn)動(dòng)特性的分析[10-12]，而對(duì)關(guān)節(jié)反力本身的物理含義及求解方法有效性問(wèn)題研究較少，這對(duì)研究桿件壽命、關(guān)節(jié)磨損特性等問(wèn)題形成一定阻礙。

針對(duì)以上問(wèn)題，本文在對(duì)拉格朗日方程中約束反力與拉格朗日乘子關(guān)系的基礎(chǔ)上，指出乘子項(xiàng)是廣義位形空間中的廣義約束力，其實(shí)質(zhì)是廣義約束力中各個(gè)約束所占的權(quán)重。利用拉格朗日增廣法計(jì)算機(jī)構(gòu)的關(guān)節(jié)反力，詳細(xì)說(shuō)明每個(gè)關(guān)節(jié)對(duì)應(yīng)關(guān)節(jié)反力的計(jì)算方法和篩選步驟。為了證明計(jì)算方法的正確性，利用牛頓歐拉方程對(duì)所選模型關(guān)節(jié)反力再次求解，并應(yīng)用ADAMS對(duì)同一模型進(jìn)行動(dòng)力學(xué)仿真，將所得的三種結(jié)果進(jìn)行了對(duì)比。本文雖然以雙搖桿機(jī)構(gòu)為例，所受約束為理想約束，但由于求解過(guò)程對(duì)約束方程與結(jié)構(gòu)特性沒(méi)有特殊要求，保證了該方法對(duì)空間任意機(jī)構(gòu)系統(tǒng)關(guān)節(jié)反力求解的適應(yīng)性和準(zhǔn)確性。

1 拉格朗日乘子的含義

第一類拉格朗日方程是應(yīng)用數(shù)學(xué)分析中的乘子法建立的動(dòng)力學(xué)方程。對(duì)于機(jī)構(gòu)系統(tǒng)中任意構(gòu)件其完整的表達(dá)式為[13]

(1)

(2)

式中：F為廣義外力；Qc為約束反力。比較式(1)和式(2)明顯可得約束力與拉格朗日乘子之間的關(guān)系

(3)

2 關(guān)節(jié)反力計(jì)算與篩選

約束方程是在建立系統(tǒng)模型時(shí)，系統(tǒng)的狀態(tài)變量必須滿足一定條件所構(gòu)成的方程。機(jī)械系統(tǒng)關(guān)節(jié)約束方程數(shù)等于該處關(guān)節(jié)反力個(gè)數(shù)且與該關(guān)節(jié)自由度數(shù)之和為3(平面)或?yàn)?(空間)[14]。設(shè)q為機(jī)構(gòu)系統(tǒng)廣義坐標(biāo)，C為機(jī)構(gòu)約束方程，對(duì)C關(guān)于廣義坐標(biāo)求一階和二階偏導(dǎo)可得機(jī)構(gòu)的速度和加速度約束方程

(4)

(5)

(6)

當(dāng)約束方程中包含驅(qū)動(dòng)約束方程時(shí)，雅克比矩陣為方陣，且通常為非奇異矩陣，即m=n，當(dāng)約束方程中不包含驅(qū)動(dòng)約束方程時(shí)，雅克比矩陣不為方陣，即m≠n。將式(6)與動(dòng)力學(xué)方程式(1)聯(lián)立可得到系統(tǒng)完整的動(dòng)力學(xué)模型

(7)

將剛體質(zhì)量矩陣M、雅克比矩陣Cq、合外力矢量Qe以及Qd代入式(7)即可求得拉格朗日乘子λm×1，其標(biāo)量形式可以寫(xiě)為

λ=[λ1，λ2,...,λm]T

(8)

假設(shè)機(jī)構(gòu)總運(yùn)動(dòng)副個(gè)數(shù)為A，構(gòu)件數(shù)為b，第a個(gè)運(yùn)動(dòng)副自由度為Fa(a的取值規(guī)則與關(guān)節(jié)反力計(jì)算結(jié)果無(wú)關(guān))，則機(jī)構(gòu)總的約束方程個(gè)數(shù)m和機(jī)構(gòu)廣義坐標(biāo)個(gè)數(shù)n分別為

(9)

圖1 機(jī)構(gòu)系統(tǒng)運(yùn)動(dòng)副Fig.1 Kinematic pairs in mechanism system

同樣，與該運(yùn)動(dòng)副相關(guān)的拉格朗日乘子(λ)k為機(jī)構(gòu)系統(tǒng)乘子項(xiàng)λ中的第l+1行至第l+3-Fk行，即

(10)

將式(9)轉(zhuǎn)置與式(10)代入式(3)得

(11)

(12)

再次取所有廣義關(guān)節(jié)反力項(xiàng)N中與該運(yùn)動(dòng)副相關(guān)的關(guān)節(jié)反力，即N中的第l+1行至第l+3-Fk行

(13)

上述計(jì)算方法和篩選步驟的整個(gè)求解過(guò)程對(duì)系統(tǒng)的結(jié)構(gòu)特性并無(wú)特殊要求，適用于空間機(jī)構(gòu)系統(tǒng)，因此可以得出選用拉格朗日增廣法求解任意機(jī)構(gòu)系統(tǒng)關(guān)節(jié)反力的一般步驟：

步驟1根據(jù)機(jī)構(gòu)運(yùn)動(dòng)學(xué)已知條件求解約束方程C，并對(duì)約束方程關(guān)于廣義坐標(biāo)求偏導(dǎo)得到雅克比矩陣Cq，提取雅克比矩陣Cq中與待求關(guān)節(jié)相關(guān)的分量并轉(zhuǎn)置；

步驟2同樣方法，提取拉格朗日乘子λ中與代求關(guān)節(jié)相關(guān)的拉格朗日乘子分量；

圖2 關(guān)節(jié)反力計(jì)算流程圖Fig.2 Calculated flow chart of reaction force

3 數(shù)值算例

以平面四桿機(jī)構(gòu)OABC為例，簡(jiǎn)圖如圖3所示。表1給出各均質(zhì)桿件結(jié)構(gòu)參數(shù)，轉(zhuǎn)動(dòng)慣量相對(duì)質(zhì)心。由于各桿長(zhǎng)關(guān)系滿足l1+l4≤l2+l3，因此該機(jī)構(gòu)為雙曲柄機(jī)構(gòu)，即驅(qū)動(dòng)桿和從動(dòng)桿均可做整周回轉(zhuǎn)。桿l2為驅(qū)動(dòng)桿，逆時(shí)針勻速轉(zhuǎn)動(dòng)，轉(zhuǎn)速50π rad/s，計(jì)算時(shí)間t從0～0.4 s即一個(gè)周期內(nèi)各構(gòu)件的位置、速度和加速度。取桿件質(zhì)心建立局部坐標(biāo)系，系統(tǒng)坐標(biāo)系O1XY，原點(diǎn)位于O，X軸與桿l1重合，方向指向C點(diǎn)。

圖3 四桿機(jī)構(gòu)簡(jiǎn)圖Fig.3 Diagram of planar four-bar mechanism

表1 四桿機(jī)構(gòu)參數(shù)Tab.1 Parameters of four bar mechanism

為了證明上述拉格朗日增廣法求解各個(gè)關(guān)節(jié)對(duì)應(yīng)關(guān)節(jié)反力選擇步驟和方法的正確性，利用牛頓歐拉方程對(duì)所選模型各關(guān)節(jié)反力重新計(jì)算。設(shè)桿件質(zhì)量為m，選用質(zhì)心坐標(biāo)為廣義坐標(biāo)q=(Rx,Ry,θ)，根據(jù)桿件受到與其相鄰構(gòu)件i和j的力和力矩平衡條件，可列出如下矢量表達(dá)式

(14)

同樣選用四桿機(jī)構(gòu)OABC為對(duì)象，利用ADAMS軟件建立三維仿真模型如圖4所示。并進(jìn)行動(dòng)力學(xué)仿真分析，求解各個(gè)關(guān)節(jié)接觸反力。

ADAMS簡(jiǎn)化模型在時(shí)間t從0～0.4 s即一個(gè)周期內(nèi)進(jìn)行仿真，得到各個(gè)關(guān)節(jié)的關(guān)節(jié)反力及x向和y向的分力，將結(jié)果與前兩節(jié)的數(shù)值結(jié)果進(jìn)行對(duì)比如圖5～圖8所示。

圖4 ADAMS仿真模型Fig.4 ADAMS simulation model

圖5 O關(guān)節(jié)反力及分量Fig.5 Reaction force and components at point O

圖7 B點(diǎn)關(guān)節(jié)反力與分量Fig.7 Reaction force and components at point B

分析圖5～圖8中數(shù)值計(jì)算與ADAMS仿真結(jié)果曲線可得以下結(jié)論：

(1)拉格朗日增廣法對(duì)所選雙搖桿機(jī)構(gòu)各關(guān)節(jié)處關(guān)節(jié)反力計(jì)算結(jié)果與其他兩種方法計(jì)算結(jié)果相比略有誤差，但不影響本次對(duì)關(guān)節(jié)反力計(jì)算方法與篩選步驟的討論。

(2)拉格朗日增廣法計(jì)算關(guān)節(jié)反力產(chǎn)生誤差的原因主要是由于初始坐標(biāo)位置精度不同、給定計(jì)算結(jié)果誤差等級(jí)不夠以及動(dòng)力學(xué)方程組求解過(guò)程中出現(xiàn)的約束違約等原因造成。

(3)三種計(jì)算方法所得各關(guān)節(jié)處關(guān)節(jié)反力變化曲線基本重合，符合四桿機(jī)構(gòu)實(shí)際的運(yùn)動(dòng)規(guī)律和動(dòng)力學(xué)特性，可以證明文中拉格朗日增廣法求解關(guān)節(jié)反力計(jì)算方法和篩選步驟的正確性。

4 結(jié) 論

(1)針對(duì)長(zhǎng)期以來(lái)對(duì)拉格朗日乘子法中涉及拉格朗日乘子是否為系統(tǒng)真正關(guān)節(jié)反力的問(wèn)題，明確了拉格朗日乘子中的乘子項(xiàng)是廣義位形空間中的廣義約束力而非真正的關(guān)節(jié)約束反力。

(2)利用拉格朗日增廣法計(jì)算機(jī)構(gòu)的關(guān)節(jié)反力，詳細(xì)說(shuō)明每個(gè)關(guān)節(jié)對(duì)應(yīng)關(guān)節(jié)反力的計(jì)算步驟和方法，并以雙搖桿機(jī)構(gòu)為例對(duì)一個(gè)周期內(nèi)的所有關(guān)節(jié)的關(guān)節(jié)反力進(jìn)行數(shù)值求解。

圖8 C點(diǎn)關(guān)節(jié)反力及分量Fig.8 Raction force and components at point C

(3)為了驗(yàn)證計(jì)算方法和選擇步驟的正確性，利用牛頓歐拉方程與ADAMS軟件對(duì)同一模型進(jìn)行動(dòng)力學(xué)數(shù)值仿真，通過(guò)使用不同方法求得的各關(guān)節(jié)反力相同，且誤差在允許范圍之內(nèi)。

(4)關(guān)節(jié)反力計(jì)算方法與篩選步驟雖然以雙搖桿機(jī)構(gòu)運(yùn)動(dòng)一周各個(gè)關(guān)節(jié)所受關(guān)節(jié)反力的計(jì)算與篩選為例，但通過(guò)第二部分論述，該方法可使用于空間任意機(jī)構(gòu)系統(tǒng)的關(guān)節(jié)反力計(jì)算。