王曉軍 王琪 莊方方

(1.常州工學(xué)院機(jī)電學(xué)院，常州 213002)(2.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191)

引言

非光滑多體系統(tǒng)動力學(xué)是在光滑多體系統(tǒng)動力學(xué)研究的基礎(chǔ)上逐步發(fā)展起來的新的研究領(lǐng)域，主要研究非光滑因素(如摩擦與碰撞)對多體系統(tǒng)動力學(xué)行為的影響.如步行機(jī)器人在地面上行走、航天器的空中對接、機(jī)械手抓取工件和具有非理想約束鉸鏈(考慮間隙與摩擦的滑移鉸和轉(zhuǎn)動鉸)的機(jī)械系統(tǒng)等都存在物體間的接觸與分離、滑移與粘滯等現(xiàn)象(稱為非光滑事件).由于這些非光滑事件的存在，導(dǎo)致系統(tǒng)的動力學(xué)方程不連續(xù)或分段連續(xù)，給非光滑事件的判斷和動力學(xué)方程的求解帶來了新的困難.

上個世紀(jì)末，Pfeiffer研究了具有單邊約束非光滑多體系統(tǒng)動力學(xué)［1］，通過引入互補概念，有效地解決了非光滑事件的判斷;隨著研究的不斷深入，逐步形成了事件驅(qū)動法和時間步進(jìn)法［2］.段文杰應(yīng)用非光滑動力學(xué)的時間步進(jìn)法研究了被動行走器足地間的庫倫摩擦系數(shù)和碰撞恢復(fù)系數(shù)對被動行走器動力學(xué)行為的影響［3］.Flores研究了具有非理想鉸鏈多體系統(tǒng)的運動學(xué)與動力學(xué)問題［4］，研究了鉸鏈的間隙、摩擦和潤滑劑等因素對多體系統(tǒng)動力學(xué)行為的影響，數(shù)值結(jié)果表明，隨著間隙的減小，由于接觸碰撞引起的法向約束力的突變也隨之減小，然而當(dāng)采用時間步進(jìn)法研究含非光滑鉸鏈多體系統(tǒng)動力學(xué)時，隨著間隙的減小，違約問題將逐步凸顯［5］.文獻(xiàn)［6］用非光滑多體系統(tǒng)動力學(xué)方法研究了滑移鉸含摩擦多體系統(tǒng)的建模與數(shù)值計算方法，當(dāng)滑移鉸的間隙充分小時，將機(jī)械系統(tǒng)中滑移鉸的滑塊視為質(zhì)點，滑塊與滑道間的幾何約束視為定常的雙邊約束，建立了雙邊約束法向力的互補關(guān)系，應(yīng)用水平線性互補［7］和事件驅(qū)動法給出了非光滑事件判斷的計算方法，應(yīng)用Baumgarte約束穩(wěn)定化方法［8］在一定程度上解決了約束的違約問題.若能將文獻(xiàn)［6］的方法推廣到非定常約束(如驅(qū)動約束)的多體系統(tǒng)，且滑移鉸的滑道不是固定的而是運動的(如曲柄搖桿機(jī)構(gòu)中的搖桿)，則可使非光滑多體系統(tǒng)動力學(xué)方法的應(yīng)用領(lǐng)域更加廣泛.

本文將研究具有非定常約束(驅(qū)動約束)及滑移鉸含摩擦的非光滑多體系統(tǒng)動力學(xué)的建模方法與數(shù)值計算方法，應(yīng)用庫侖干摩擦模型作為滑移鉸間的摩擦模型，應(yīng)用線性互補方法建立滑移鉸法向約束力的互補關(guān)系，將約束分為幾何的定常約束(滑移鉸約束)和非定常約束(驅(qū)動約束)，應(yīng)用具有約束穩(wěn)定化的增廣法和事件驅(qū)動法建立該系統(tǒng)的動力學(xué)方程，最后通過算例說明本文給出方法的有效性.

1 非光滑滑移鉸力學(xué)模型

1.1 滑移鉸幾何約束模型

滑移鉸是機(jī)械系統(tǒng)中常見的運動副.在有些機(jī)構(gòu)中，當(dāng)間隙充分小時，滑移鉸被視為雙邊約束且其中的滑塊被視為質(zhì)點［6］，此時滑移鉸的力學(xué)模型如圖1所示，其中是滑道兩側(cè)分別作用在滑塊上的法向約束力.設(shè)系統(tǒng)中的第i個滑移鉸的約束方程為

式中，q=［q1，…，qk］為系統(tǒng)的廣義坐標(biāo)，n*為滑移鉸的個數(shù).若用遞推法［9］或用距離函數(shù)列寫滑移鉸的約束方程，則約束方程(1)對應(yīng)的Lagrange乘子為滑道作用在滑塊i上的法向約束力.

圖1 滑移鉸模型Fig.1 The model of translational joints

當(dāng)λNi=0時，滑道與滑塊無接觸，如圖1(a)所示;當(dāng)λNi＞0時，滑道的一側(cè)與滑塊接觸，如圖1(b)所示;當(dāng)λNi＜0時，滑道的另一側(cè)與滑塊接觸，如圖1(c)所示.滑道作用于滑塊上的兩個法向約束力，與Lagrange乘子λNi具有下列關(guān)系［6］

且滿足下列互補條件

設(shè)

則由式(2)可得

1.2 滑移鉸摩擦模型

機(jī)械系統(tǒng)中常用的摩擦模型有多種［10］，其中庫侖摩擦模型又可分為庫侖干摩擦模型和修正的庫侖摩擦模型［4］，前者是相對速度的非連續(xù)函數(shù)(給數(shù)值計算帶來一定的困難)，后者是相對速度的連續(xù)函數(shù)(不易反映庫侖摩擦的靜動態(tài)特性)［6］.本文將采用庫侖干摩擦模型作為滑移鉸的摩擦模型.

應(yīng)用非光滑動力學(xué)方法，庫侖干摩擦模型可表示為

式中，F(xiàn)fi為滑道作用在滑塊上的摩擦力在切向上的投影;μi，μ0i分別為滑道與滑塊間的動、靜摩擦因數(shù);|λNi|為作用于滑塊上法向約束力的大小;Sgn()為符號函數(shù);vri，˙vri分別為滑塊相對滑道的相對速度和相對切向加速度;Sgn()為集值函數(shù)［11，12］，可表示為

由式(4)和式(5)，可以看出，當(dāng)滑塊的相對速度和相對切向加速度均為零時，摩擦力的取值是一個范圍，即當(dāng)滑道內(nèi)的滑塊處于粘滯狀態(tài)時，摩擦力的取值在該范圍內(nèi).

2 非光滑動力學(xué)方程

2.1 多體系統(tǒng)動力學(xué)方程

第一類Lagrange方程是建立多體系統(tǒng)動力學(xué)方程的有效方法之一.設(shè)滑移鉸的約束方程由方程(1)表示，將其用向量形式表示為

式中，Φ=［φ1，…，φn*］T.設(shè)系統(tǒng)的驅(qū)動約束方程為

則由第一類Lagrange方程可得到系統(tǒng)的動力學(xué)方程為

式中，P中各元素為系統(tǒng)廣義坐標(biāo)及其對時間一階和二階導(dǎo)數(shù)的函數(shù)=［|λN1|，…，|λNn*|］T.

將式(9)代入式(8)，并應(yīng)用Baumgarte約束穩(wěn)定化方法，方程(8)可表示為

式中，

方程(10)～(12)為具有約束穩(wěn)定化的動力學(xué)方程.當(dāng)系統(tǒng)是光滑時，P=0，方程組(10)～(12)是關(guān)于，λN的線性代數(shù)方程組，可用相關(guān)的數(shù)值計算方法求解，但是對于非光滑系統(tǒng)，P≠0中含有|λNi|(i=1，…，n*)，則該方程組不是關(guān)于，λN的線性代數(shù)方程組，不能用線性代數(shù)方程組的數(shù)值計算方法求解.

2.2 線性互補算法

利用式(3)，可將方程組(10)～(12)表示成線性互補方程.為便于推導(dǎo)，將方程(10)～(12)表示成

式中，

將(13)式代入(15)式，可求得

式中，

將式(16)代入(13)得

將上式代入式(14)得

再將上式表示為

式中，

將式(3)代入式(17)可得

將上式表示成

且

式中，

式(18)和式(19)為標(biāo)準(zhǔn)的水平線性互補問題(HLCP)［6］，應(yīng)用線性互補的計算方法可求出，;然后將其代入式(3)可得到λN，;再將其代入式(16)，則可求出驅(qū)動約束力(或力偶)~λ.將相關(guān)的約束力代入方程(10)，應(yīng)用常微分方程的數(shù)值計算方法即可求出該系統(tǒng)的運動特性.

3 算例

3.1 力學(xué)模型

本文以曲柄搖桿機(jī)構(gòu)為例，如圖2所示，其中搖桿OA和曲柄O1B均視為剛體，滑塊B視為質(zhì)點.

該系統(tǒng)參數(shù)分別為:搖桿對O軸的轉(zhuǎn)動慣量為J1=16/3kg·m2，質(zhì)量為m1=4.0kg，其質(zhì)心到O軸的距離為L1=1.0m;曲柄對O1軸的轉(zhuǎn)動慣量為J2=1/6kg·m2，質(zhì)量為m2=2.0kg，其質(zhì)心到O1軸的距離為L2=0.25m;滑塊的質(zhì)量為m3=0.5kg，到O1軸的距離為L3=0.5m;O軸到O1軸的距離為L4=1.0m;滑塊與滑道間的動滑動摩擦因數(shù)為μ=0.2;曲柄的角速度ω=3.0rad/s.

圖2 曲柄-搖桿機(jī)構(gòu)Fig.2 Crank-rocker mechanism

系統(tǒng)的廣義坐標(biāo)為q=［θ1，θ2］T，分別為搖桿和曲柄與鉛垂軸y的夾角;γ=θ2-θ1為搖桿與曲柄間的夾角.設(shè)滑移鉸B的約束方程和曲柄的驅(qū)動約束方程分別為

由于該系統(tǒng)在運動過程中，滑塊相對搖桿無相對靜止?fàn)顟B(tài)(當(dāng)滑塊的相對速度為零時，其相對加速度不為零)，因此滑塊相對滑道無粘滯狀態(tài)，則摩擦力的廣義力可表示為

式中，vr為滑塊相對搖桿的相對速度.

設(shè)搖桿在運動過程中受到線性阻力矩的作用(Mf=-c)，其中為阻尼系數(shù).當(dāng)約束方程約束的是線位移時，對應(yīng)的Lagrange乘子是約束力，當(dāng)約束方程約束的是角位移時，對應(yīng)的Lagrange乘子為約束力偶矩.

3.2 數(shù)值分析

根據(jù)約束方程(20)和驅(qū)動約束方程(21)可知，曲柄勻角速度轉(zhuǎn)動，搖桿往復(fù)擺動.應(yīng)用本文給出的計算方法對該系統(tǒng)進(jìn)行數(shù)值仿真，圖3給出了搖桿的角速度(實線)和曲柄的角速度(虛線)的時間歷程.圖4給出了搖桿的角加速度與曲柄轉(zhuǎn)角θ2的對應(yīng)關(guān)系，從圖中可以看出，當(dāng)θ2=0rad或θ2=πrad時，搖桿的角加速度為零(這與定性分析的結(jié)果相吻合).

圖3 和的時間歷程圖Fig.3 The time history of and

圖4 搖桿OA的角加速度圖Fig.4 Angular acceleration of rocker OA

圖5 λN和Ff圖(無阻尼)Fig.5 λN and Ff with c=0

圖6給出了當(dāng)搖桿的阻尼系數(shù)不為零時(即:c=8.0N·m·s/rad)，作用于滑塊上的摩擦力Ff與廣義坐標(biāo)θ2的對應(yīng)關(guān)系.當(dāng)θ2=0rad或θ2=πrad時，有θ1=0°=0.0rad/s2，λN≠0且滑塊相對速度的方向發(fā)生改變，導(dǎo)致Ff的值發(fā)生突變.

圖6 Ff～θ2(c=8.0N·m·s/rad)Fig.6 Ff～θ2(c=8.0N·m·s/rad)

上述仿真結(jié)果與用牛頓-歐拉方法建立的動力學(xué)方程求得的數(shù)值結(jié)果完全吻合.

在本算例中，若α＞1.0，β＞1.0，在數(shù)值計算時，設(shè)計算步長為h，當(dāng)計算步長滿足0＜h≤0.001時，約束方程滿足若α=0，β=0，h=0.001，約束方程不能被滿足，其計算結(jié)果發(fā)散.

4 結(jié)論

本文研究了非光滑滑移鉸平面多剛體系統(tǒng)驅(qū)動力的數(shù)值計算方法.應(yīng)用第一類Lagrange方程建立了該系統(tǒng)的動力學(xué)方程，將滑移鉸視為雙邊約束，分別建立驅(qū)動約束方程和滑移鉸的運動約束方程，與驅(qū)動約束方程對應(yīng)的Lagrange乘子為驅(qū)動力或驅(qū)動力偶矩，與滑移鉸的運動約束方程對應(yīng)的Lagrange乘子為作用在滑移鉸的法向約束力;滑移鉸的摩擦模型采用庫侖摩擦模型，為便于計算摩擦力的廣義力，建立了滑移鉸法向約束力的互補關(guān)系，將其法向約束力的計算轉(zhuǎn)化為線性互補方程的求解;結(jié)合Baumgarte約束穩(wěn)定化方法和常微分?jǐn)?shù)值計算方法，給出了該系統(tǒng)驅(qū)動力及系統(tǒng)動力學(xué)響應(yīng)的數(shù)值計算方法.最后以曲柄-搖桿機(jī)構(gòu)為例，通過數(shù)值仿真說明了本文給出方法的有效性.

1 Pfeiffer F，Glocker C.Multibody dynamics with unilateral contacts.New York:John Wiley＆Sons.Inc.，1996

2 Pfeiffer F，F(xiàn)oerg M，Ulbrich H.Numerical aspects of nonsmooth multibody dynamics.Computer Methods in Applied Mechanics and Engineering，2006，195:6891～6908

3 段文杰，王琪，王天舒.圓弧足被動行走器非光滑動力學(xué)仿真研究.力學(xué)學(xué)報，2011，43(4):765～774(Duan W J，Wang Q，Wang T S.Simulation research of a passive dynamic walker with round feet based on non-smooth method.Chinese Journal of Theoretical and Applied Mechanics，2011，43(4):765～774(in Chinese))

4 Flores P，Ambrosio J，Pimenta Claro J C，et al.Kinematics and dynamics of multibody systems with imperfect joints.Berlin:Springer-Verlag，2008

5 Flores P，Leine R，Glocker C.Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach.Multibody System Dynamics，2010，23:165～190

6 Wang Q，Peng H L，Zhuang F F.A constraint-stabilized method for multibody dynamics with friction-affected translational joints based on HLCP.Discrete and Continuous Dynamical Systems Series B，2011，16(2):589～605

7 韓繼業(yè)，修乃華，戚厚鐸.非線性互補理論與算法.上海:上?？茖W(xué)出版社，2006(Han JY，Xiu N H，Qi H D.Theory and algorithm of nolinear complementarity.Shanghai:Shanghai Scientific＆Technical Publishers，2006(in Chinese))

8 Flores P，Machado M，Seabra E，Silva M T.A parametric study on the Baumgarte stabilization method for forward dynamics of constrained multibody systems.ASME Journal of Computational and Nonlinear Dynamics，2011，6:011019

9 洪嘉振.計算多體系統(tǒng)動力學(xué).北京:高等教育出版社，2001(Hong J Z.Calculation of dynamics of multibody system.Beijing:Higher Education Press，2006(in Chinese))

10 劉麗蘭，劉宏昭，吳子英等.機(jī)械系統(tǒng)中摩擦模型的研究進(jìn)展.力學(xué)進(jìn)展，2008，38(2):200～213(Liu L L，Liu H Z，Wu Z Y，et al.An overview of friction models in mechanical systems.Advances in Mechanics，2008，38(2):200～213(in Chinese))

11 Zhuang F F，Wang Q.Modeling and simulation of the nonsmooth planar rigid multibody systems with frictional translational joints.Multibody System Dynamics，2013，29:403～423

12 Acary Vincent，Brogliato Bernard.Numerical methods for nonsmooth dynamical systems.Berlin:Springer-Verlag，2008