基于線性互補(bǔ)問(wèn)題的多體系統(tǒng)接觸/碰撞動(dòng)力學(xué)研究

張欣剛，齊朝暉，王 剛，國(guó)樹(shù)東，吳志剛

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

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

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

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

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

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

1.1 連續(xù)化模型

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

0≤gi⊥τi≥0

(1)

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

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

(2)

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

(3)

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

(a)

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

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

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

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

(4)

(5)

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

(6)

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

(7)

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

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

(8)

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

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

(9)

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

接觸力的虛功率可表示為

(10)

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

(11)

(12)

所有接觸力虛功率之和為

(13)

式中

Kf=[Tf1e1,Tf2e2,…,Tfnen]

(14)

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

(15)

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

(16)

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

(17)

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

(18)

(19)

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

(20)

并通過(guò)求解其等價(jià)的非線性方程組(7)來(lái)求解接觸力f，其中

(21)

(22)

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

用均勻化的縫隙函數(shù)與法向接觸力建立互補(bǔ)方程，實(shí)際上是將突變的接觸力進(jìn)行了光滑化處理，當(dāng)物體臨近接觸時(shí)接觸力已經(jīng)開(kāi)始起作用，而當(dāng)縫隙函數(shù)及其變化率均為零時(shí)，方程又可退化為平順接觸時(shí)的接觸力方程。一方面能夠解出接觸力隨時(shí)間的演化，另一方面嚴(yán)格滿足互補(bǔ)條件，因此不會(huì)解出負(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號(hào)質(zhì)點(diǎn)在外力作用下向左移動(dòng)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ù)接觸力模型對(duì)柔體本構(gòu)關(guān)系的影響，我們以Flores模型為例，分別計(jì)算不同恢復(fù)系數(shù)條件下的系統(tǒng)響應(yīng)。圖3分別給出了不同恢復(fù)系數(shù)、相同時(shí)刻不同恢復(fù)系數(shù)兩種情況下的系統(tǒng)機(jī)械能時(shí)程曲線，圖中cr表示恢復(fù)系數(shù)。

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

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

(a) 時(shí)程曲線

(b) 相同時(shí)刻不同恢復(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

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

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

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

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

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

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

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

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

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

圖10給出了圖9相應(yīng)的接觸力與實(shí)際縫隙函數(shù)的關(guān)系。與連續(xù)接觸力模型相似，兩者均可視為局部碰撞的連續(xù)化。但由于采用均勻化縫隙與接觸力建立互補(bǔ)關(guān)系，因此物體在發(fā)生碰撞時(shí)也會(huì)有一定的嵌入量。另外，與沖量形式的摩擦碰撞線性互補(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ū)間長(zhǎng)度h對(duì)系統(tǒng)動(dòng)力學(xué)響應(yīng)以及接觸力有著很大的影響，圖11給出了不同光滑參數(shù)h下滑塊3的位移時(shí)程曲線。

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

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

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

(a) 位移

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

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

(23)

(24)

由于單元應(yīng)變?cè)诠残鴺?biāo)系中為小量，因此式(8)可用單元?jiǎng)偠染仃嚤硎?，并按照?24)修改為

(25)

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

(26)

式中：Ke為單元?jiǎng)偠染仃?；ue為單元局部參數(shù)，描述截面在單元坐標(biāo)系下的位移和轉(zhuǎn)動(dòng)。

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

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

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

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

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

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

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

4 結(jié) 論

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

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