張欣剛，齊朝暉，王 剛，國樹東，吳志剛

(1. 青島理工大學(xué) 理學(xué)院，青島 266525; 2. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024;3. 大連理工大學(xué) 海洋科學(xué)與技術(shù)學(xué)院，遼寧 盤錦 124221)

隨著智能機(jī)器人、航空航天等領(lǐng)域的發(fā)展，輕質(zhì)化、柔性化、復(fù)雜化、快速化的機(jī)械系統(tǒng)得到了廣泛的使用，同時也促進(jìn)了多柔體系統(tǒng)動力學(xué)學(xué)科的蓬勃發(fā)展[1]。以航天科技為例，其空間結(jié)構(gòu)往往以折疊狀態(tài)固定在艙內(nèi)，入軌后再逐步展開，如大型桁架式索網(wǎng)天線，大型太陽能電池陣,柔性抓取機(jī)械臂等[2]。為了控制航天器的重量，這些結(jié)構(gòu)多為大型輕質(zhì)柔性結(jié)構(gòu)。這些柔性構(gòu)件在展開過程中不僅呈現(xiàn)出大轉(zhuǎn)動、大變形剛?cè)狁詈系膭恿W(xué)特性，其部件也不可避免的要與周圍環(huán)境或自身發(fā)生接觸/碰撞[3]，由此產(chǎn)生的非光滑振動不僅給航天器的姿態(tài)控制帶來較大的困難，甚至導(dǎo)致航天器的失穩(wěn)[4]。因此，深入分析結(jié)構(gòu)柔性以及碰撞、摩擦對系統(tǒng)動力學(xué)特性的影響具有切實(shí)的意義。

對碰撞的描述可分為兩類：① 忽略碰撞的細(xì)節(jié)，認(rèn)為碰撞在瞬間完成[5-7]；② 用彈簧模擬物體對變形的抵抗能力，用阻尼模擬能量耗散[8-10]。第一類方法又被稱為離散分析方法，優(yōu)點(diǎn)是簡化了接觸過程且計(jì)算效率高，缺點(diǎn)在于不能計(jì)算碰撞過程中碰撞力的大小以及作用過程，其線性補(bǔ)償策略可能會違反摩擦接觸時的能量守恒定律[11]。第二類方法又被稱為連續(xù)分析方法，允許接觸體在接觸區(qū)域相互嵌入，接觸力可表示為侵入量的明確函數(shù)[12]。其中Hertz接觸理論是使用范圍最廣研究范圍最悠久的接觸力模型。文獻(xiàn)[13]系統(tǒng)總結(jié)了Hertz接觸模型的發(fā)展，并著重討論了模型中阻尼系數(shù)的發(fā)展過程。

將多剛體系統(tǒng)中應(yīng)用最為廣泛的連續(xù)接觸力模型直接引入到多柔體接觸/碰撞動力學(xué)分析中并非完全適用[14]。連續(xù)接觸力模型將接觸力描述為變形量的明確函數(shù)，相當(dāng)于引入一個本構(gòu)關(guān)系。而柔體自身有其本構(gòu)關(guān)系。引入連續(xù)接觸力模型刻畫柔體間的碰撞，相當(dāng)于給系統(tǒng)疊加了本構(gòu)關(guān)系，從而使得兩者相互干擾：一方面，模型中剛度和阻尼參數(shù)的確定含有很大的不確定性[15-17]；另一方面，由于一次碰撞歷程中相對速度可在較大的范圍內(nèi)取值，加之受到兩類本構(gòu)關(guān)系的干擾，往往會計(jì)算出負(fù)的接觸力，從而與接觸力的性質(zhì)相違背。基于線性互補(bǔ)理論的非光滑動力學(xué)模型能夠保證接觸力的非負(fù)性，且能夠高精度的求解持續(xù)接觸狀態(tài)下的接觸力，不足之處在于需要區(qū)分不同的接觸狀態(tài)，且無法同時定性和定量的分析碰撞時刻的接觸力演化過程[18]。

結(jié)合以上研究現(xiàn)狀，本文首先分析了連續(xù)接觸力模型在多柔體接觸/碰撞分析中的不足，隨后在非光滑動力學(xué)的體系下，將多柔體系統(tǒng)模型降噪[19]的思想引入柔體間動態(tài)接觸力的建模與計(jì)算當(dāng)中，建立了平均化縫隙與接觸力之間的標(biāo)準(zhǔn)線性互補(bǔ)方程。這種互補(bǔ)關(guān)系可以綜合考慮碰撞和平順接觸，在接觸狀態(tài)發(fā)生改變時不需切換模型。不僅保證了接觸力的非負(fù)性，也在很大程度上簡化了非光滑互補(bǔ)問題的分析和求解過程。

1 柔體間動態(tài)接觸力模型

1.1 連續(xù)化模型

接觸力可視為單面約束提供的約束反力，其重要力學(xué)特性是約束反力的單向性：沿縫隙函數(shù)gi增加方向的接觸力值τi必須非負(fù)，而當(dāng)縫隙函數(shù)gi大于零時τi必須為零。因此縫隙函數(shù)和法向接觸力必然滿足以下關(guān)系式

0≤gi⊥τi≥0

(1)

盡管上述關(guān)系符合客觀事實(shí)，但由于縫隙函數(shù)本身與接觸力并不存在明確的函數(shù)關(guān)系，因此我們僅能通過它定性的確定接觸力是否為零，而難以確定其幅值。

為了刻畫柔體間碰撞力的演化過程，一些學(xué)者將連續(xù)接觸力引入到柔性多體系統(tǒng)接觸/碰撞動力學(xué)中，例如應(yīng)用很廣泛的Lankarani和Nikravesh模型

(2)

式(1)的一個特例是，當(dāng)單面約束正在起作用時，由于縫隙函數(shù)及其一階變化率均為零，因此退化為平順接觸時的線性互補(bǔ)關(guān)系

(3)

由動力學(xué)普遍原理可知，縫隙函數(shù)的二階變化率與接觸力存在函數(shù)關(guān)系，因此接觸力可以被確定，大量算例也證實(shí)了這種關(guān)系的真實(shí)性。由此可見，平順接觸條件下不需人為設(shè)定接觸力模型也可以高精度的求解接觸力。事實(shí)上，連續(xù)接觸力模型將圖1(a)物體之間的瞬時碰撞抽象為圖1(b)所示的等效系統(tǒng)。由圖1可見，人為設(shè)定柔體之間的接觸力模型相當(dāng)于給系統(tǒng)疊加兩個本構(gòu)關(guān)系，從而帶來了一定程度的干擾。

(a)

(b)圖1 柔體間接觸力模型Fig.1 Contact force model between two flexible bodies

1.2 互補(bǔ)問題的光滑化方法

類比庫倫摩擦定律可見：一些物理上正確的數(shù)學(xué)描述往往無法用于數(shù)值分析。當(dāng)所研究的系統(tǒng)無stick-slip運(yùn)動狀態(tài)的切換時，人們往往采用修正的庫倫摩擦模型來代替。將間斷函數(shù)連續(xù)化也是處理非光滑問題的一類被驗(yàn)證可行的方法。

類似地，將縫隙函數(shù)在短時間[t,t+h]內(nèi)做泰勒展開，并保留到2階,則ξ時刻的縫隙函數(shù)可表示為

(4)

(5)

由此可將接觸力與縫隙函數(shù)的互補(bǔ)關(guān)系近似為

(6)

由于平均化縫隙函數(shù)表達(dá)式(5)中含有縫隙函數(shù)二階變化率項(xiàng)，因此可確定接觸力。經(jīng)此近似以后的互補(bǔ)方程可以綜合考慮接觸與碰撞，且不必將接觸過程人為的劃分為多個狀態(tài)。與式(6)等價的非線性方程形式為

(7)

2 接觸/碰撞的線性互補(bǔ)方程

建立均勻化縫隙與接觸力之間的標(biāo)準(zhǔn)線性互補(bǔ)方程，還需推導(dǎo)縫隙加速度與廣義加速度之間以及廣義加速度與接觸力之間的函數(shù)關(guān)系。建立多柔體系統(tǒng)動力學(xué)方程的物理依據(jù)是虛功率原理，可表述為彈性體外力虛功率等于彈性體慣性力虛功率與彈性體變形虛功率之和

(8)

式中：r為物體內(nèi)一點(diǎn)的矢徑；σ和ε分別為該點(diǎn)處的廣義應(yīng)力和廣義應(yīng)變；f為作用在該點(diǎn)的外力。當(dāng)給系統(tǒng)增加外力時，僅需計(jì)算該作用力的虛功率并疊加到系統(tǒng)虛功率方程中即可。

不考慮接觸力時，系統(tǒng)的虛功率方程可由廣義坐標(biāo)q表示為

(9)

式中:M為廣義質(zhì)量陣;F為廣義力陣。

接觸力的虛功率可表示為

(10)

將接觸點(diǎn)虛速度用廣義速度表示，則式(10)可改寫為

(11)

(12)

所有接觸力虛功率之和為

(13)

式中

Kf=[Tf1e1,Tf2e2,…,Tfnen]

(14)

f=[f1,f2,…,fn]T

(15)

接觸力對廣義力的貢獻(xiàn)即為Kff，如果式(9)所選的廣義坐標(biāo)獨(dú)立，則系統(tǒng)動力學(xué)方程可由此寫為

(16)

由式(16)，可將廣義加速度由接觸力線性表示

(17)

根據(jù)運(yùn)動學(xué)關(guān)系，我們也總能將縫隙函數(shù)的變化率表示為廣義坐標(biāo)變化率的函數(shù)

(18)

(19)

綜合式(7)以及式(17)～(19),并將其代入式(5),我們可以通過標(biāo)準(zhǔn)線性互補(bǔ)問題

(20)

并通過求解其等價的非線性方程組(7)來求解接觸力f，其中

(21)

(22)

求解式(20)得到接觸力后，代入式(17)即可得到廣義加速度。

用均勻化的縫隙函數(shù)與法向接觸力建立互補(bǔ)方程，實(shí)際上是將突變的接觸力進(jìn)行了光滑化處理，當(dāng)物體臨近接觸時接觸力已經(jīng)開始起作用，而當(dāng)縫隙函數(shù)及其變化率均為零時，方程又可退化為平順接觸時的接觸力方程。一方面能夠解出接觸力隨時間的演化，另一方面嚴(yán)格滿足互補(bǔ)條件，因此不會解出負(fù)接觸力。

3 數(shù)值算例

例1：圖2所示的彈性質(zhì)量系統(tǒng)，各質(zhì)量點(diǎn)均為mi=m=1 kg，彈簧剛度ki=k=104N/m。靜止?fàn)顟B(tài)下間隙為零。左側(cè)3號質(zhì)點(diǎn)在外力作用下向左移動0.03 m，達(dá)到靜平衡后釋放并與右側(cè)系統(tǒng)發(fā)生碰撞。分別采用彈簧阻尼模型以及所提線性互補(bǔ)方法進(jìn)行數(shù)值仿真。彈性阻尼模型的等效剛度取140×106N/m3/2。

圖2 彈簧質(zhì)量系統(tǒng)Fig.2 Spring-mass system

為了進(jìn)一步討論連續(xù)接觸力模型對柔體本構(gòu)關(guān)系的影響，我們以Flores模型為例，分別計(jì)算不同恢復(fù)系數(shù)條件下的系統(tǒng)響應(yīng)。圖3分別給出了不同恢復(fù)系數(shù)、相同時刻不同恢復(fù)系數(shù)兩種情況下的系統(tǒng)機(jī)械能時程曲線，圖中cr表示恢復(fù)系數(shù)。

由圖3(a)可見，恢復(fù)系數(shù)為1時為純彈性碰撞，系統(tǒng)機(jī)械能守恒。但隨著恢復(fù)系數(shù)的降低，在相同時刻，系統(tǒng)機(jī)械能并沒有隨著恢復(fù)系數(shù)的降低而降低，而是呈現(xiàn)一定的波動。圖3(b)給出了相同時刻不同恢復(fù)系數(shù)條件下系統(tǒng)機(jī)械能的變化,更為直觀了展示了這種波動。

除以上情況外，計(jì)算得到的恢復(fù)系數(shù)與預(yù)設(shè)的恢復(fù)系數(shù)也呈現(xiàn)了較大的差異，圖4給出了不同接觸力模型設(shè)定恢復(fù)系數(shù)與計(jì)算恢復(fù)系數(shù)的比較。

(a) 時程曲線

(b) 相同時刻不同恢復(fù)系數(shù)圖3 系統(tǒng)機(jī)械能Fig.3 System energy

圖4 預(yù)設(shè)與計(jì)算恢復(fù)系數(shù)關(guān)系Fig.4 Relation between the post and pre-restitution coefficients

對比圖4以及文獻(xiàn)[8]可見，F(xiàn)lores模型在計(jì)算剛體間碰撞時恢復(fù)系數(shù)與預(yù)設(shè)值吻合較好，但計(jì)算柔體間碰撞時就產(chǎn)生了較大的偏差。其余兩種經(jīng)典接觸力模型[9,10]在計(jì)算剛體間高恢復(fù)系數(shù)碰撞問題時吻合較好，但計(jì)算柔體碰撞時也產(chǎn)生了一定的偏差。

除了以上兩種情況外，將接觸力模型引入到柔體間接觸問題往往得到負(fù)的接觸力。其主要源于兩個方面：① 相對速度可在較大范圍內(nèi)取值；② 阻尼項(xiàng)中分母若含有恢復(fù)系數(shù)，則在低恢復(fù)系數(shù)時阻尼項(xiàng)絕對值被放大。如圖5所示。其中各接觸力模型及其適用范圍可參照文獻(xiàn)[20]的歸納總結(jié)。

由圖5可見，對于分母中含有恢復(fù)系數(shù)的遲滯阻尼因子，在低恢復(fù)系數(shù)下其幅值就會被放大，加之相對速度在較大范圍內(nèi)取值，在低恢復(fù)系數(shù)情況下很容易計(jì)算出小于負(fù)1的阻尼項(xiàng)，如圖6所示(Flores模型)。一次碰撞過程中接觸力與侵入深度的關(guān)系如圖7所示。當(dāng)物體間相對速度為正(脫離)時，相對速度與初始相對速度之商為負(fù)，由于模型中阻尼項(xiàng)系數(shù)被分母中的恢復(fù)系數(shù)放大，則在物體脫離對方的時刻計(jì)算出了如圖7所示的負(fù)接觸力。

圖5 阻尼因子與恢復(fù)系數(shù)關(guān)系Fig.5 Relation between damping factor and restitution coefficients

圖6 阻尼項(xiàng)時程曲線Fig.6 Time plot of damping item

圖7 接觸力/嵌入量關(guān)系Fig.7 Contact force-deformation relation

采用所提線性互補(bǔ)方法時，實(shí)際縫隙函數(shù)與均勻化縫隙函數(shù)的關(guān)系如圖8所示。由均勻化縫隙函數(shù)的表達(dá)式(5)可知，物體相互靠近時，縫隙函數(shù)大于零，其一階變化率小于零，因此均勻化縫隙比真實(shí)縫隙要小一些，則在物體實(shí)際碰撞前的一個較短時間內(nèi)接觸力已經(jīng)開始作用。并且h值越大，接觸力作用的時間越早，接觸力的峰值也隨之降低。如圖9所示。

圖8 縫隙函數(shù)與均勻化縫隙函數(shù)關(guān)系Fig.8 Relation between normal gap and time-averaged gap

圖9 不同參數(shù)h下的接觸力時程曲線Fig.9 Time lapse of contact force with different h

圖10給出了圖9相應(yīng)的接觸力與實(shí)際縫隙函數(shù)的關(guān)系。與連續(xù)接觸力模型相似，兩者均可視為局部碰撞的連續(xù)化。但由于采用均勻化縫隙與接觸力建立互補(bǔ)關(guān)系，因此物體在發(fā)生碰撞時也會有一定的嵌入量。另外，與沖量形式的摩擦碰撞線性互補(bǔ)模型[21-23]相比，優(yōu)點(diǎn)在于可直接求解接觸力的大小，且碰撞前后的狀態(tài)依賴于柔體自身的本構(gòu)關(guān)系，不必人為的設(shè)定沖量恢復(fù)系數(shù)。

圖10 接觸力與縫隙函數(shù)的關(guān)系Fig.10 Contact force-gap relation

縫隙函數(shù)均勻化區(qū)間長度h對系統(tǒng)動力學(xué)響應(yīng)以及接觸力有著很大的影響，圖11給出了不同光滑參數(shù)h下滑塊3的位移時程曲線。

結(jié)合圖11以及圖9可見，均勻化區(qū)間長度h取值越大，接觸力作用時間越長，碰撞效應(yīng)將被磨平，位移曲線較為光滑。隨著h值的減小，位移曲線逐漸吻合，更趨近于真實(shí)響應(yīng)，但過低的h值也會使得所求接觸力失真。建議h的取值應(yīng)與碰撞時間在相同數(shù)量級。

圖11 滑塊3位移時程曲線Fig.11 Time lapse of the displacement of slider 3

實(shí)際上，除了碰撞力本身外，人們更為關(guān)心的是碰撞前后系統(tǒng)的動力學(xué)行為。圖12給出了分別采用所提LCP方法以及連續(xù)接觸力模型(CCF)計(jì)算的滑塊3位移時程曲線以及接觸力時程曲線，其中連續(xù)接觸力模型中恢復(fù)系數(shù)取0.3。從中可見，除了連續(xù)接觸力模型會計(jì)算出負(fù)接觸力外，兩種模型所求接觸力略有差異，但系統(tǒng)動力學(xué)響應(yīng)基本吻合。

(a) 位移

(b) 接觸力圖12 兩種模型對比Fig.12 Comparison between CCF and LCP

例2本算例將給出所提方法在柔性多體系統(tǒng)接觸/碰撞動力學(xué)中的應(yīng)用?；谖墨I(xiàn)[19]所提的模型降噪方法以及文獻(xiàn)[24]所提的大變形空間梁的共旋坐標(biāo)法，將梁的廣義應(yīng)力修正為短時間內(nèi)的平均應(yīng)力

(23)

(24)

由于單元應(yīng)變在共旋坐標(biāo)系中為小量，因此式(8)可用單元剛度矩陣表示，并按照式(24)修改為

(25)

因此，修正后的節(jié)點(diǎn)內(nèi)力為

(26)

式中：Ke為單元剛度矩陣；ue為單元局部參數(shù)，描述截面在單元坐標(biāo)系下的位移和轉(zhuǎn)動。

對于圖13所示的柔性雙擺機(jī)構(gòu)，材料密度設(shè)為ρ=7 850 kg/m3，彈性模量E=2.1×109Pa，兩根擺長均為L=2 m，截面積均為A=3×10-4m2。機(jī)構(gòu)在重力加速度g=9.81 m/s2作用下，從水平靜止位置落下并于剛性基礎(chǔ)產(chǎn)生碰撞，假定平面光滑無摩擦。

圖13 自由下落并與基礎(chǔ)碰撞的柔性雙擺機(jī)構(gòu)Fig.13 Free-falling flexible double pendulum mechanismimpact with foundation

采用MATLAB ode45常微分方程求解器對該問題進(jìn)行數(shù)值積分，取精度控制參數(shù)εr=1.0×10-4和εa=1.0×10-6。得到雙擺結(jié)構(gòu)在不同時刻的構(gòu)型圖和碰撞力時程，分別如圖14和圖15所示。

圖14 雙擺構(gòu)型圖Fig.14 Deformed configurations of flexible double pendulum

由圖14和圖15可見，當(dāng)柔性雙擺的材料剛度較小時(E=2.1×109Pa)，第二根擺在碰撞后變形很大。由于材料較柔，碰撞后桿端并沒有彈離地面，而是進(jìn)入平順接觸狀態(tài)，沿平面滑動一段距離后脫離地面。

圖15 接觸力時程曲線Fig.15 Time lapse of normal contact force

多柔體系統(tǒng)的剛性以及接觸非線性極大的增加了數(shù)值求解的困難，數(shù)值解的高頻振蕩也往往導(dǎo)致接觸力求解失真。剛性積分器的算法阻尼往往不能有效地衰減高頻振蕩，往往通過在模型中引入少許線性阻尼或結(jié)構(gòu)阻尼的方式進(jìn)一步濾除高頻振蕩。采用文獻(xiàn)[19]所提的模型降噪方法時可采用常規(guī)的積分器(如ode45)進(jìn)行求解。此外，為了避免磨平碰撞效應(yīng)，在濾除過高頻率彈性振蕩的前提下，模型中的降噪因子應(yīng)盡可能的小，以提高數(shù)值解的精度。

4 結(jié) 論

為刻畫柔體間動態(tài)接觸力，選取當(dāng)前時刻為起點(diǎn)的短時區(qū)間[t,t+h]內(nèi)對縫隙函數(shù)進(jìn)行均勻化，建立了均勻化縫隙函數(shù)與接觸力之間的線性互補(bǔ)關(guān)系，該方法可綜合考慮碰撞和平順接觸，不必在接觸狀態(tài)發(fā)生改變時切換模型，同時保證了接觸力的非負(fù)性。為了降低高頻成分對求解接觸力的干擾，引入模型降噪方法濾除了過高頻率的高頻振蕩成分，提高了數(shù)值求解的穩(wěn)定性。

數(shù)值算例表明，所得結(jié)果在定性和定量兩個方面都是合理的。該方法不僅適用于單點(diǎn)接觸，同時也適用于柔體間的多點(diǎn)摩擦接觸問題。將所提均勻化線性互補(bǔ)模型與模型降噪方法相結(jié)合，可成為求解大型復(fù)雜柔性多體系統(tǒng)接觸碰撞問題的新途徑。