王曉軍 王琪 莊方方

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

引言

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

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

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

1 非光滑滑移鉸力學模型

1.1 滑移鉸幾何約束模型

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

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

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

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

且滿足下列互補條件

設(shè)

則由式(2)可得

1.2 滑移鉸摩擦模型

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

應用非光滑動力學方法，庫侖干摩擦模型可表示為

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

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

2 非光滑動力學方程

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

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

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

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

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

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

式中，

方程(10)～(12)為具有約束穩(wěn)定化的動力學方程.當系統(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)表示成線性互補方程.為便于推導，將方程(10)～(12)表示成

式中，

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

式中，

將式(16)代入(13)得

將上式代入式(14)得

再將上式表示為

式中，

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

將上式表示成

且

式中，

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

3 算例

3.1 力學模型

本文以曲柄搖桿機構(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 曲柄-搖桿機構(gòu)Fig.2 Crank-rocker mechanism

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

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

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

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

3.2 數(shù)值分析

根據(jù)約束方程(20)和驅(qū)動約束方程(21)可知，曲柄勻角速度轉(zhuǎn)動，搖桿往復擺動.應用本文給出的計算方法對該系統(tǒng)進行數(shù)值仿真，圖3給出了搖桿的角速度(實線)和曲柄的角速度(虛線)的時間歷程.圖4給出了搖桿的角加速度與曲柄轉(zhuǎn)角θ2的對應關(guān)系，從圖中可以看出，當θ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給出了當搖桿的阻尼系數(shù)不為零時(即:c=8.0N·m·s/rad)，作用于滑塊上的摩擦力Ff與廣義坐標θ2的對應關(guān)系.當θ2=0rad或θ2=πrad時，有θ1=0°=0.0rad/s2，λN≠0且滑塊相對速度的方向發(fā)生改變，導致Ff的值發(fā)生突變.

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

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

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

4 結(jié)論

本文研究了非光滑滑移鉸平面多剛體系統(tǒng)驅(qū)動力的數(shù)值計算方法.應用第一類Lagrange方程建立了該系統(tǒng)的動力學方程，將滑移鉸視為雙邊約束，分別建立驅(qū)動約束方程和滑移鉸的運動約束方程，與驅(qū)動約束方程對應的Lagrange乘子為驅(qū)動力或驅(qū)動力偶矩，與滑移鉸的運動約束方程對應的Lagrange乘子為作用在滑移鉸的法向約束力;滑移鉸的摩擦模型采用庫侖摩擦模型，為便于計算摩擦力的廣義力，建立了滑移鉸法向約束力的互補關(guān)系，將其法向約束力的計算轉(zhuǎn)化為線性互補方程的求解;結(jié)合Baumgarte約束穩(wěn)定化方法和常微分數(shù)值計算方法，給出了該系統(tǒng)驅(qū)動力及系統(tǒng)動力學響應的數(shù)值計算方法.最后以曲柄-搖桿機構(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 段文杰，王琪，王天舒.圓弧足被動行走器非光滑動力學仿真研究.力學學報，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出版社，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)動力學.北京:高等教育出版社，2001(Hong J Z.Calculation of dynamics of multibody system.Beijing:Higher Education Press，2006(in Chinese))

10 劉麗蘭，劉宏昭，吳子英等.機械系統(tǒng)中摩擦模型的研究進展.力學進展，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