浮放物體平面多剛體動力學建模與算法研究1)

張潤森 王 琪

(北京航空航天大學航空科學與工程學院,北京100191)

動力學與控制

浮放物體平面多剛體動力學建模與算法研究1)

張潤森 王 琪2)

(北京航空航天大學航空科學與工程學院,北京100191)

采用非光滑多體系統(tǒng)動力學的方法研究浮放物體與基礎平臺組成的多體系統(tǒng)，建立其非光滑接觸的動力學方程與數(shù)值算法.浮放物體由主體部分和支撐腿組成，其間通過含黏彈性阻力偶的轉(zhuǎn)動鉸連接.支撐腿與基礎平臺間的接觸力簡化為接觸點的法向接觸力和摩擦力，采用擴展的赫茲接觸力模型描述接觸點的法向接觸力，采用庫倫干摩擦模型描述其摩擦力.采用笛卡爾坐標系下的位形坐標作為系統(tǒng)的廣義坐標.首先，將基礎平臺運動看作非定常約束，用第一類拉格朗日方程建立系統(tǒng)的動力學方程，并采用鮑姆加藤約束穩(wěn)定化的方法解決違約問題.隨后給出基于事件驅(qū)動法和線性互補方法的數(shù)值算法.當相對切向速度為零時，構(gòu)造靜滑動摩擦力的正負余量和正、負向加速度的互補關系，從而將接觸點黏滯--滑移切換的判斷以及靜滑動摩擦力的計算轉(zhuǎn)化為線性互補問題進行求解，并采用Lemke算法求解線性互補問題.最后，通過數(shù)值仿真選擇合適的步長；通過仿真結(jié)果說明浮放物體運動中存在的黏滯--滑移切換現(xiàn)象以及基礎平臺運動、質(zhì)心位置對浮放物體運動的影響.

浮放物體，非光滑，庫倫干摩擦，線性互補，接觸力

引言

海上回收的火箭(如圖1所示)、文物、家具、儀器設備等以浮放方式放置于基礎平臺上的物體統(tǒng)稱為浮放物體，其受到基礎平臺運動的影響而產(chǎn)生滑移、搖晃、傾覆、甚至墜落等現(xiàn)象而受到損壞.為此，國內(nèi)外學者對浮放物體的動力學特性展開了大量研究.

圖1 獵鷹9號火箭Fig.1 Rocket Falcon IX

Housner[1]將浮放物體簡化為為對稱均質(zhì)的塊狀剛體，假設浮放物體與基礎平臺之間為點接觸且無相對滑動，采用經(jīng)典碰撞理論(沖量與動量方法)建立了浮放物體的動力學方程，并分析了浮放物體傾覆時的最小脈沖加速度.

學者們在Housner模型的基礎上對浮放物體的搖晃運動進行了大量理論和實驗的研究，并取得了階段性的成果[2-4].Zhang等[5]對基礎平臺不同運動規(guī)律下浮放物體的搖晃運動進行了研究，得到了浮放物體不同的傾覆方式；Shi等[6]對非對稱浮放物體的搖晃運動進行了數(shù)值模擬，研究了其動力學行為.Shenton等[7-8]指出基礎平臺激勵下浮放物體存在著搖晃、滑移、搖晃--滑移、相對基礎平臺靜止等運動形式，并推導了浮放物體搖晃、滑移、搖晃--滑移3種運動形式下的動力學方程，指出了考慮浮放物體與基礎平臺之間相對滑移的重要性.周乾[7]，Roussis等[8]，Taniguchi等[9]，Contento 等[10]，Pompei等[11]在Shenton研究的基礎上展開了進一步的研究.研究了考慮浮放物體與基礎平臺間的相對運動時，基礎平臺不同運動形式下浮放物體的動力學行為.這些研究基于Housner關于碰撞的假設，碰撞在瞬間完成，且為完全塑性碰撞.

不同于 Housner模型，一些文獻考慮了浮放物體和基礎平臺的變形，采用離散單元方法 (distinct/discrete element method)研究了浮放物體的動力學行為[12-13].當浮放物體與基礎平臺發(fā)生接觸時，該方法將法向接觸力和切向接觸力表示為接觸點相對位移和相對速度的函數(shù)，然后將法向接觸力和切向接觸力代入動力學方程進行求解.

在某些情況下，浮放物體和基礎平臺可以視為一個多體系統(tǒng)[14].然而與行進的被動行走器[15]、考慮間隙的非理想鉸[16-17]、行駛的車輛[18]、著陸的飛機[19]等非光滑多體系統(tǒng)相同，浮放物體和基礎平臺之間存在的含摩擦接觸問題給動力學方程的求解帶來了困難.Zhang等[20]和Yilmaz等[21]分別將非光滑多體系統(tǒng)動力學的方法應用于浮放物體的研究中，研究了浮放物體搖晃運動時的動力學行為，但均未考慮浮放物體相對于基礎平臺的滑移運動.Zhuang等[17]，Pfei ff er等[22]將相對切向速度為零時接觸點黏滯--滑移的判斷以及靜摩擦力的求解等非光滑問題轉(zhuǎn)化為線性互補問題進行求解，與試算法相比減少了對非光滑事件的判斷，更易于編程計算，為浮放物體和基礎平臺組成的多體系統(tǒng)求解提供了思路.

本文采用非光滑多體系統(tǒng)動力學理論，研究了含支撐腿浮放物體與基礎平臺組成的平面多剛體系統(tǒng)，建立了其動力學方程與數(shù)值算法.通過數(shù)值仿真，揭示了接觸點摩擦狀態(tài)黏滯--滑移的切換現(xiàn)象，說明了基礎平臺運動規(guī)律、質(zhì)心的偏移等對浮放物體運動的影響.不同于Housner采用的完全塑性碰撞假設或離散單元法，本文將考慮浮放物體和基礎平臺間接觸部位的局部變形，用擴展的赫茲接觸力模型描述接觸點的法向接觸力[23]，用庫倫干摩擦模型描述摩擦力.采用第一類拉格朗日方程建立系統(tǒng)的動力學方程，并用鮑姆加藤違約修正方法改善數(shù)值仿真過程中的違約問題[24].利用事件驅(qū)動算法和文獻[16-17]的思路，將接觸點黏滯--滑移狀態(tài)的判斷轉(zhuǎn)化為線性互補問題進行求解.但不同于文獻[16-17]中的方法，本文方法不需要將法向接觸力的計算轉(zhuǎn)化為線性互補問題，從而降低了線性互補問題的互補維數(shù)，提高了線性互補問題的計算效率.最后選擇合適的步長進行數(shù)值仿真，通過算例反映了浮放物體運動過程中的黏滯--滑移現(xiàn)象以及基礎平臺運動規(guī)律、質(zhì)心的偏移對浮放物體運動的影響，定性說明了方法的正確性.

1 系統(tǒng)的描述及接觸力模型

1.1 動力學模型

將含支撐腿的浮放物體簡化為由主體部分、支撐腿構(gòu)成，并與基礎平臺組成平面多體系統(tǒng)，其模型如圖2所示.將支撐腿與主體的連接簡化為含有黏彈性阻力偶的轉(zhuǎn)動鉸A和B，其等效扭轉(zhuǎn)剛度系數(shù)為k，等效扭轉(zhuǎn)阻尼系數(shù)為c.考慮基礎平臺的平移運動，其轉(zhuǎn)動也可以按照本文給出的方法進行研究.浮放物體主體部分為非均質(zhì)剛體，質(zhì)心位置如圖2中C1所示.將兩條支撐腿為均質(zhì)細桿，C2和C3分別表示兩條支撐腿的質(zhì)心.

圖2 浮放物的多體模型Fig.2 Model of multi-body system of free-standing body

為敘述方便，用剛體0，1，2，3依次代表基礎平臺，浮放物體主體部分以及左右兩條支撐腿.剛體1質(zhì)心C1到鉸鏈A和B的距離分別記作r1和r2.剛體2和剛體3具有相同的長度2l0，則質(zhì)心C2和C3到鉸鏈的距離為l0.令 2b0=r1sinα1+r2sinα2，h0=r2cosα2=r1cosα1，e=r1sinα1?b0，則e為剛體1質(zhì)心相對于幾何對稱軸的偏移量.用mi(i=1,2,3)表示剛體i的質(zhì)量，Ji(i=1,2,3)表示剛體i相對質(zhì)心Ci的轉(zhuǎn)動慣量.g為重力加速度.

建立如圖2中所示的直角坐標系，其中x軸指向水平向右，y軸指向鉛錘向上.用(xi,yi)(i=0,1,2,3)表示剛體i的質(zhì)心坐標.用θi(i=1，2,3)表示剛體i的轉(zhuǎn)角，以逆時針方向為正.則該系統(tǒng)的廣義坐標可表示為

1.2 法向接觸力模型

不同于Hounser所采用的經(jīng)典碰撞理論，擴展的赫茲接觸力模型可以更好地反映碰撞過程中法向接觸力的大小及其變化規(guī)律[23].當物體發(fā)生接觸時，將碰撞力表示為嵌入量和嵌入速度的非線性函數(shù)[25-27]，本文采用文獻 [26]提出的連續(xù)碰撞力模型，其表達式為

其中,FNi表示作用于接觸點i的法向接觸力.δi和i分別表示接觸點i對應的法向嵌入量和嵌入速度.指數(shù)n和廣義剛度系數(shù)K1可以通過理論分析得到也可以直接通過實驗擬合獲得，而廣義阻尼系數(shù)c2與法向嵌入速度無關，往往通過實驗擬合得到[26].法向嵌入量δi可表示為系統(tǒng)廣義坐標的函數(shù)

將上式對時間求導可得法向嵌入速度

1.3 切向接觸力模型

學者們根據(jù)不同的摩擦現(xiàn)象提出了多種摩擦模型，如庫倫干摩擦模型，庫倫 +黏性摩擦模型，Stribeck摩擦模型.庫倫干摩擦模型簡潔地反映了摩擦的庫倫現(xiàn)象以及黏滯--滑移切換現(xiàn)象，因而本文采用了庫倫干摩擦模型[28].庫倫干摩擦模型是關于相對切向速度的多值函數(shù)，會給數(shù)值計算帶來一定的困難，因而學者們提出了關于相對切向速度連續(xù)的修正的庫倫摩擦模型，但是其無法反映相對切向速度為零時的摩擦特性.為了反映摩擦黏滯--滑移切換現(xiàn)象以及相對切向速度為零時的摩擦特性，本文將采用庫倫干摩擦模型描述切向接觸力.

根據(jù)庫倫干摩擦模型，剛體i(i=2,3)受到的摩擦力與法向接觸力的關系可以表示為

其中，F(xiàn)fi和FNi(i=2,3)分別表示剛體i與基礎平臺接觸時受到的摩擦力和法向接觸力.和分別表示接觸點i的相對切向速度和相對切向加速度.μi和μ0i(i=2,3)分別表示接觸點i與基礎之間的動、靜滑動摩擦系數(shù).sgn()是符號函數(shù)，其函數(shù)表達式如下

Sgn()是集值函數(shù)，其定義如下

2 系統(tǒng)的動力學方程及其算法

2.1 動力學方程

在圖2所示系統(tǒng)中，設fi(q)=0(i=1,2,3,4)分別為鉸鏈A和B的約束方程，其表達式如下

將基礎平臺的水平方向和鉛錘方向的平移運動看作非定常約束，其約束方程可以表示為

將約束方程(8)和(9)用矩陣形式表示為

由第一類拉格朗日方程可得系統(tǒng)的動力學方程

其中，M為系統(tǒng)的廣義質(zhì)量矩陣；為約束矩陣Φ關于q的雅可比矩陣的轉(zhuǎn)置；λ為對應于約束方程(10)的拉格朗日乘子列向量；QN為法向接觸力的廣義力列向量；Qf為摩擦力的廣義力列向量；Q為重力以及鉸鏈A和B處黏彈性力偶矩對應的廣義力列向量，表達式為

將法向接觸力和摩擦力的廣義力分別表示為

其中，WN和WT是與廣義坐標有關的系數(shù)矩陣.FN和Ff分別由基礎平臺作用于剛體i(i=2,3)的法向接觸力FNi和摩擦力Ffi組成，可表示為

將式(13)和式(14)代入動力學方程中，可將式(11)中的第1式改寫為

2.2 摩擦力的互補關系

庫倫干摩擦模型在零點的不連續(xù)性給動力學方程(17)的求解帶來了困難，其難點在于接觸點摩擦狀態(tài)黏滯--滑移切換(非光滑事件)的判斷.本文借鑒前人的工作，將非光滑事件的判斷轉(zhuǎn)化為線性互補問題進行求解.

接觸點相對切向加速度的正向加速度和負向加速度分別定義為[28]

從而上述物理量的互補關系可以表述為

2.3 摩擦力的廣義力化簡

為便于動力學方程的求解,在此將摩擦力的廣義力進行化簡.由式(14)可得

將式(22)代入式(21)可得

對于矩陣S，存在矩陣B，滿足

當矩陣S存在非零元素時，矩陣B由矩陣S中的非零列向量組成.當矩陣S不存在非零元素時，即為零矩陣時，取B=S.

將式(24)代入式(23)可得

則將式(25)代入式(17)中可得

2.4 動力學方程的算法

本文采用鮑姆加藤約束穩(wěn)定化的方法解決數(shù)值求解動力學方程中的違約問題.將式(11)中的約束方程改寫為[24]

其中，α和β是大于零的常數(shù)，其大小與計算步長有關，可以根據(jù)具體問題通過試算得到合理的取值[24].令

則式(27)可表示為

由式(26)可得

將式(29)代入式(28)可得

將式(30)代入式(29)可得

其中

接觸點的相對切向加速度列向量可表示為

其中，列向量wT是廣義坐標和廣義速度的函數(shù).

將式(33)代入式(32)可得

由式(19)可得

將式(35)和式(31)代入式(34)可得

其中

由式(18)可得

結(jié)合式(36)和式(37)可得

因為在式(38)中，滿足線性互補條件

由此將接觸點相對切向速度為0時，接觸點摩擦狀態(tài)的判斷和靜滑動摩擦力大小的求解轉(zhuǎn)化為線性互補問題，可以采用萊姆克(Lemke)算法，轉(zhuǎn)軸算法等最優(yōu)化方法進行求解.給出第n步到第n+1步的計算流程圖，如圖3所示.

圖3 算法流程圖Fig.3 Simulation fl ow chat

3 數(shù)值仿真結(jié)果

3.1 步長的選取

Wrinkle[13]指出當使用顯式算法(如本文采用的四階龍格--庫塔方法)，需采用合適的時間步長以保證數(shù)值穩(wěn)定性.

如圖4所示，等效扭轉(zhuǎn)剛度系數(shù)為k=1.0×105N，等效扭轉(zhuǎn)阻尼系數(shù)為c=0N·s，各剛體的質(zhì)量分別為：m1=1.2kg，m2=m3=0.1kg.動/靜滑動摩擦系數(shù)分別為：μ02=μ03=0.6，μ2=μ3=0.48.法向接觸力的參數(shù)分別為K1=2.0×107N/m1.5，c2=0N·m/s2.圖5為一對稱浮放物體從高度0.4m處自由下落時，剛體1質(zhì)心縱坐標y1的時間歷程圖，與Wrinkle的計算結(jié)果保持一致.在系統(tǒng)運動過程中，摩擦力不做功，系統(tǒng)機械能守恒，因而每次彈起高度保持一致.隨計算步長h減小，計算結(jié)果趨于收斂，步長h=1.0×10?6s與h=1.0×10?7s時計算結(jié)果幾乎重合，故選計算步長h=1.0×10?6s.

圖4 數(shù)值穩(wěn)定性測試算例Fig.4 The case of numerical stability test

圖5 y1時間歷程圖Fig.5 Time history of y1

3.2 數(shù)值仿真結(jié)果

表1為系統(tǒng)數(shù)值仿真過程中各參數(shù)取值.其中θ20和θ30分別表示鉸鏈處力偶矩為零時θ2和θ3的大小.令xxd=x1?x0?x1(0)，則xxd表示剛體1相對于基礎平臺的水平坐標.動、靜摩擦系數(shù)滿足μi=0.8μ0i.

情形1

靜滑動摩擦系數(shù)μ02=μ03=0.6，基礎平臺做水平運動，X0=0.8m，ωx= πrad/s，αx=0.初始時，浮放物體放置于基礎平臺上，相對基礎平臺靜止.圖6分別給出了剛體1質(zhì)心偏移量e=?0.915m,e=0 m時xxd的時間歷程圖.當e=?0.915m時，由于浮放物體非對稱，其左右搖晃幅度不同，將會沿基礎平臺產(chǎn)生單方向運動；而對于浮放物體對稱的情況，浮放物體在基礎平臺水平激勵下做周期運動.圖7和圖8分別給出質(zhì)心偏移量e=?0.915m時接觸點相對切向速度τi、法向接觸力FNi和切向接觸力Ffi的時間歷程圖.從圖7可以看到接觸點的摩擦狀態(tài)發(fā)生黏滯--滑移切換現(xiàn)象，圖8所展示的法向接觸力和切向接觸力變化規(guī)律與摩擦點狀態(tài)的切換保持同步.

表1 系統(tǒng)參數(shù)取值Table 1 Parameters’value of the system

圖6 xxd時間歷程圖Fig.6 Time history of xxd

圖7 τi時間歷程圖Fig.7 Time history of τi

圖8 FNi和Ffi時間歷程圖Fig.8 Time history of FNiand Ffi

情形2

靜滑動摩擦系數(shù)μ02=μ03=0.2，基礎平臺做鉛錘運動，Y0=0.99mm，ωy=30.0πrad/s，αy=0.h0=3.536m，b0=7.208m，其余系統(tǒng)參數(shù)如表1所示.初始時浮放物體放置于基礎平臺，且有鉛錘方向的相對速度，i=0?1.0(i=1，2，3).圖9分別給出了剛體1質(zhì)心偏移量e=?3.673m和e=0m時xxd的時間歷程圖.當e=?3.673m時，由于浮放物體非對稱，其左右搖晃幅度不同，在基礎平臺鉛錘激勵下有水平方向運動；當e=0m時，對稱浮放物體始終未產(chǎn)生水平方向運動.圖10給出了e=?3.673m時接觸點法向接觸力FNi的時間歷程圖.從圖10可以看出浮放物體與基礎平臺發(fā)生碰撞后，法向接觸力呈現(xiàn)周期性的變化規(guī)律.

圖9 xxd時間歷程圖Fig.9 Time history of xxd

圖10 FNi時間歷程圖Fig.10 Time history of FNi

情形3

靜滑動摩擦系數(shù)μ02=μ03=0.8，基礎平臺既做水平運動，又做鉛錘運動，X0=1.0m，ωx=ωy=πrad/s，αx=αy=0，b0=2.5m，e=0m，其余系統(tǒng)參數(shù)如表1所示.初始時浮放物體放置于基礎平臺，且相對基礎平臺靜止.圖11分別給出了基礎平臺鉛錘激勵Y0=0.9m和Y0=0m時xxd的時間歷程圖.當Y0=0.9m時，由于基礎平臺鉛錘運動和水平運動的耦合，浮放物體朝單方向運動；當Y0=0m時，基礎平臺做水平運動，xxd呈現(xiàn)周期的變化規(guī)律.

圖11 xxd時間歷程圖Fig.11 Time history of xxd

4 結(jié)論

本文基于非光滑多體系統(tǒng)動力學方法，建立了含支撐腿浮放物體與基礎平臺所組成的平面多剛體系統(tǒng)的動力學方程與數(shù)值算法.用擴展的赫茲接觸力模型描述支撐腿與基礎平臺之間的法向接觸力，用庫倫干摩擦模型描述摩擦力.將第一類拉格朗日方程與包姆加藤約束穩(wěn)定化方法結(jié)合建立系統(tǒng)的動力學方程.由于摩擦的存在使動力學方程不連續(xù)，為此建立正負摩擦余量與正、負向加速度，將靜滑動摩擦力的求解和黏滯--滑移的判斷轉(zhuǎn)化為線性互補問題的求解.

最后在數(shù)值仿真時，采用Wrinkle的方法選擇合適的步長以保證數(shù)值仿真的穩(wěn)定性.通過仿真算例分析了浮放物體運動中存在的黏滯--滑移問題，定性地說明了方法的正確性與特定系統(tǒng)參數(shù)下基礎平臺運動規(guī)律和質(zhì)心的偏移對浮放物體運動的影響.當基礎平臺做水平運動時，對稱浮放物體水平方向上做周期運動，而存在質(zhì)心偏移量時，其將沿基礎平臺單方向運動；當基礎平臺做鉛錘運動時，對稱浮放物體無水平方向上相對運動，而存在質(zhì)心偏移量時，其將沿基礎平臺單方向運動；若基礎平臺做水平方向和鉛錘方向的耦合運動，對稱浮放物體水平方向上將沿基礎平臺單方向運動.

1 Housner GW.The behavior of inverted pendulum structures during earthquakes.Bulletin of the Seismological Society of America,1963,53(2):403-417

2 Caliò I,Marletta M.Passive control of the seismic rocking response of art objects.Engineering Structures,2003,25(8):1009-1018

3 Zhang J,Makris N.Rocking response of free-standing blocks under cycloidal pulses.Journal of Engineering Mechanics,2001,127(5):473-483

4 Shi B,Anooshehpoor A,Zeng Y,et al.Rocking and overturning of precariously balanced rocks by earthquakes.Bulletin of the Seismological Society of America,1996,86(5):1364-1371

5 Shenton HW,Jones NP.Base excitation of rigid bodies.I:Formulation.Journal of Engineering Mechanics,1991,117(10):2286-2306

6 Shenton HW,Jones NP.Base excitation of rigid bodies.II:Periodic slide-rock response.Journal of Engineering Mechanics,1991,117(10):2307-2328

7 周乾,閆維明,紀金豹.地震作用下浮放文物滑移及搖晃響應仿真.華北地震科學,2016,34(1):13-20(Zhou Qian,Yan Weiming,Ji Jinbao.Simulation of oscillation and sliding responses of freestanding cultural relics under earthquakes.North China Earthquake Sciences,2016,34(1):13-20(in Chinese))

8 Roussis P,Odysseos S.Multi-mode response of base-isolated rigid blockstogroundexcitation//Proceedingofthe10thNationalConference in Earthquake Engineering,Earthquake Engineering Reasearch Institude,Anchorage,AK,2014

9 Taniguchi T.Non-linear response analyses of rectangular rigid bodies subjected to horizontal and vertical ground motion.Earthquake Engineering&Structural Dynamics,2010,31(8):1481-1500

10 Contento A,Egidio AD.Investigations into the bene fi ts of base isolation for non-symmetric rigid blocks.Earthquake Engineering&Structural Dynamics,2010,38(7):849-866

11 Pompei A,Scalia A,Sumbatyan MA.Dynamics of rigid block due to horizontal ground motion.Journal of Engineering Mechanics,1998,124(7):713-717

12 Psycharis IN,Lemos JV,Papastamatiou DY,et al.Numerical study of the seismic behaviour of a part of the parthenon pronaos.Earthquake Engineering&Structural Dynamics,2003,32(13):2063-2084

13 Winkler T,Meguro K,Yamazaki F.Response of rigid body assemblies to dynamic excitation.Earthquake Engineering&Structural Dynamics,1995,24(10):1389-1408

14 Zhang H,Brogliato B,Liu C.Dynamics of planar rocking-blocks with Coulomb friction and unilateral constraints:Comparisons between experimental and numerical data.Multibody System Dynamics,2014,32(1):1-25

15 段文杰,王琪,王天舒.圓弧足被動行走器非光滑動力學仿真研究.力學學報,2011,43(4):765-774(Duan Wenjie,Wang Qi,Wang Tianshu.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))

16 Zhuang F,Wang Q.Modeling and analysis of rigid multibody systems with driving constraints and frictional translation joints.Acta Mechanica Sinica,2014,30(3):437-446

17 Zhuang F,Wang Q.Modeling and simulation of the nonsmooth planar rigid multibody systems with frictional translational joints.Multibody System Dynamics,2013,29(4):403-423

18 范新秀,王琪.車輛縱向非光滑多體動力學建模與數(shù)值算法研究.力學學報，2015,47(2):301-309(Fan Xinxiu,Wang Qi.Research on modeling and simulation of longitudinal vehicle dynamics based on non-smooth dynamics of multibody systems.Chinese Journal of Theoretical and Applied Mechanics,2015,47(2):301-309(in Chinese))

19 徐梓堯,王琪.含單邊非完整約束飛機滑跑的建模與仿真方法.北京航空航天大學學報,2015,41(5):835-840(Xu Ziyao,Wang Qi.Method for modeling and simulation of aircraft taxiing with unilateral and non-holonomic constraints.Journal of Beijing University of Aeronautics and Astronautics,2015,41(5):835-840(in Chinese)).

20 Zhang H,Brogliato B.The planar rocking-block:analysis of kinematic restitution laws,and a new rigid-body impact model with friction.Industrial&Commercial Training.2011,43(7):451-459

21 Yilmaz C,Gharib M,Hurmuzlu Y.Solving frictionless rocking block problem with multiple impacts.Proceedings of the Royal Society A Mathematical Physical&Engineering Sciences,2009,465(2111):3323-3339

22 Pfei ff er F,Glocker C.Multibody Dynamics with Unilateral Contacts.New York:John Wiley&Sons.Inc,1996

23 董富祥,洪嘉振.多體系統(tǒng)動力學碰撞問題研究綜述.力學進展,2009,39(3):352-359(Dong Fuxiang,Hong Jiazhen.Review of impact problem for dynamics of multibody system.Advances in Mechanics,2009,39(3):352-359(in Chinese))

24 Flores P,Machado M,Seabra E,et al.A parametric study on the baumgartestabilization methodfor forwarddynamics of constrained multibody systems.Journal of Computational&Nonlinear Dynamics,2011,6(1):11019

25 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(2):165-190

26 Yigit AS,Ulsoy AG,Scott RA.Spring-dashpot models for the dynamics of a radially rotating beam with impact.Journal of Sound&Vibration,1990,142(3):515-525

27 Flores P,Ambr′osio J,Claro JP.Dynamic analysis for planar multibody mechanical systems with lubricated joints.Multibody System Dynamics,2004,12(1):47-74

28 劉麗蘭,劉宏昭,吳子英等.機械系統(tǒng)中摩擦模型的研究進展.力學進展,2008,38(2):201-213(Liu Lilan,Liu Hongzhao,Wu Ziying,et al.An overview of friction models in mechanical systems.Advances in Mechanics,2008,38(2):201-213(in Chinese))

RESEARCH ON MODELING AND NUMERICAL METHOD OF FREE STANDING BODY ON PLANAR RIGID MULTIBODY DYNAMICS1)

Zhang Runsen Wang Qi2)
(School of Aeronautic Science and Engineering，Beihang University，Beijing100191，China)

A multibody system composed of free standing body and the basic platform is investigated by non-smooth dynamics of multibody system.Dynamic equations and numerical method of the system with non-smooth contacts are proposed.The free standing body consists of main body and supporting legs,which are connected by revolute joints with viscoelastic moments.The contact forces between free standing body and the basic platform are simpli fi ed as normal forces and frictional forces of contact points.Moreover,the modi fi ed Hertz contact model and Coulomb’s law for dry friction are employed to describe normal forces and frictional forces,respectively.And the con fi guration coordinates of Cartesian coordinate system are used as the generalized coordinates.Firstly,the system’s dynamic equations are established by Lagrange’s equations of the fi rst kind and the motion of basic platform is regarded as a rheonomous constraint.The problem of constraints violations is solved by Baumgarte stabilization method.Secondly,the numerical method of the multibody system are proposed,which is based on the event-driven schemes and linear complementarity formulations.The complementary formulations of friction saturations and the relative accelerations in the tangential are given,while the relative velocities in the tangential of the contact points are equal to zero.Therefore the judgements of stick-slip transitions for contact points and the solutions of frictional forces in stick situation could be solved as a linear complementarity problem.And the linear complementarity problem is solved by Lemke’s algorithm.Finally,an appropriate step is chosen by the simulation.Then the numerical simulations denote the stick-slip phenomenon and the in fl uence of basic platform as well as mass centre’s position.

free standing body,non-smooth,Coulomb dry friction,linear complementarity problem,contact force

O313.3

A doi：10.6052/0459-1879-17-235

2017–06–25 收稿，2017–09–29 錄用,2017–09–29 網(wǎng)絡版發(fā)表.

1)國家自然科學基金資助項目(11372018，11772021).

2)王琪,教授,主要研究方向:多體系統(tǒng)動力學與控制.E-mail:bhwangqi@sina.com

張潤森,王琪.浮放物體平面多剛體動力學建模與算法研究.力學學報,2017,49(6):1370-1379

Zhang Runsen,Wang Qi.Research on modeling and numerical method of free standing body on planar rigid multibody dynamics.Chinese Journal of Theoretical and Applied Mechanics,2017,49(6):1370-1379