軟體尺蠖爬行機(jī)器人建模與仿真分析

張力文， 徐齊平， 劉錦陽

(上海交通大學(xué) 工程力學(xué)系， 上海 200240)

不同于具有高剛度和復(fù)雜結(jié)構(gòu)的傳統(tǒng)剛性機(jī)器人，軟體機(jī)器人具有結(jié)構(gòu)簡單、人機(jī)交互安全性強(qiáng)、靈活性高、與復(fù)雜工作環(huán)境兼容性好的優(yōu)點，是目前機(jī)器人行業(yè)中最具吸引力的領(lǐng)域之一.軟體機(jī)器人主要利用形狀記憶合金、橡膠材料、硅膠、介電彈性體及芳綸纖維等柔性材料制成[1-5]，能夠完成蠕動、游動、越障、跳躍及抓握等復(fù)雜的動作[6-10].其驅(qū)動方式主要包括流體驅(qū)動、形狀記憶合金驅(qū)動及電活性聚合物驅(qū)動等[11-13].

目前，軟體機(jī)器人的研究還處于起步階段，研究者們通常從自然界生物的外形結(jié)構(gòu)和運動模式獲得靈感，設(shè)計適用于特殊環(huán)境作業(yè)的仿生機(jī)器人.Rafsanjani等[5]研制出了一種由硅橡膠管組成的蛇形機(jī)器人，當(dāng)機(jī)器人充氣膨脹時，附著于機(jī)器人表面的可拉伸的塑料薄片會被迫彈出，從而錨定于地面上，整個機(jī)器人即可利用這些錨定的塑料片的拉動前進(jìn)，這種蠕動爬行的運動方式為眾多仿蠕蟲和仿蛇的軟體爬行機(jī)器人的運動研究提供了新思路.王江北等[14]設(shè)計了一種多氣囊軟體爬行仿生機(jī)器人，通過對其硅膠材料的本構(gòu)關(guān)系的研究，得出了機(jī)器人的運動步幅和內(nèi)部氣壓的關(guān)系，并利用實驗驗證了該理論關(guān)系的可行性，為此類氣動驅(qū)動的軟體爬行機(jī)器人爬行步幅的實驗研究提供了理論指導(dǎo).

在所有的仿生軟體爬行機(jī)器人中，尺蠖形四足軟體爬行機(jī)器人靈活性好、實用性強(qiáng)，是很多科研人員的重點研究對象.哈佛大學(xué)的Shepherd及其團(tuán)隊對該類四足軟體爬行機(jī)器人進(jìn)行了實驗研究，通過對機(jī)器人四肢上簡單的閥門系統(tǒng)進(jìn)行充放氣控制，可實現(xiàn)四足之間的協(xié)調(diào)運動，從而能完成較為復(fù)雜的爬行動作[15].隨后，他們還利用高彈性硅膠材料制造出了抗干擾性更強(qiáng)的軟體機(jī)器人，該機(jī)器人由電動空氣壓縮機(jī)驅(qū)動，可以在嚴(yán)寒、炎熱等惡劣條件下長時間自主行進(jìn)，并能抵抗高強(qiáng)度沖擊、碾壓等不利因素的干擾[16].雖然Shepherd及其團(tuán)隊對四足軟體爬行機(jī)器人的結(jié)構(gòu)、材料、運動步態(tài)以及適應(yīng)能力進(jìn)行了系統(tǒng)的研究，但是未對軟體機(jī)器人的建模進(jìn)行深入研究.

氣動驅(qū)動的軟體機(jī)器人由多個氣腔構(gòu)成，對其進(jìn)行力學(xué)建模是一個復(fù)雜的過程.目前，對氣動驅(qū)動的軟體機(jī)器人的建模方法主要包括兩種：基于歐拉伯努利理論的曲梁建模方法和非線性有限元建模方法.加州大學(xué)伯克利分校的O’Reilly等[17-18]用曲梁建模方法對四足軟體爬行機(jī)器人運動時和地面的黏附情況進(jìn)行了建模和理論分析，得出了與地面線接觸情況下軟體機(jī)器人運動的穩(wěn)定性條件.隨后，他們還對于該類機(jī)器人和地面的摩擦接觸情況進(jìn)行了分析，提出了一種能夠預(yù)測其運動方向的基于摩擦定律的彈性桿理論模型[19].文獻(xiàn)[20]還基于已經(jīng)得出的穩(wěn)定性條件和摩擦模型，將四足軟體爬行機(jī)器人的力學(xué)模型簡化為曲梁模型，對該類機(jī)器人的運動情況進(jìn)行了理論上的建模和分析.然后，他們繼續(xù)通過實驗和理論推導(dǎo)，得出了四足軟體爬行機(jī)器人單個氣動致動器的曲率、彎矩以及內(nèi)部壓強(qiáng)之間的關(guān)系，實驗上驗證了理論模型的正確性[21].此外，Majidi等[19]基于歐拉伯努利曲梁模型和庫倫摩擦理論建立了機(jī)器人運動方向和接觸面上摩擦力之間的理論關(guān)系.Matia等[22]采用基于歐拉伯努利梁模型估算了驅(qū)動的氣壓和變形之間的函數(shù)關(guān)系.Suzumori等[7,23]采用非線性有限元方法進(jìn)行仿真計算，能夠得到相對準(zhǔn)確的結(jié)果，但是計算時間較長，仿真效率較低.

以上學(xué)者對尺蠖四足軟體爬行機(jī)器人的運動提出了較為完善的理論模型，但是他們在研究過程中未考慮軟體機(jī)器人在不同運動階段之間的聯(lián)系，因此也未能反映機(jī)器人在整個運動周期內(nèi)的連續(xù)變化情況.

本文針對軟體尺蠖爬行機(jī)器人建立了滑塊-曲梁的簡化模型，并對運動控制方程進(jìn)行離散和無量綱化，利用牛頓迭代法求解準(zhǔn)靜態(tài)條件的邊值問題.首先對軟體機(jī)器人在氣壓作用下的彎曲變形進(jìn)行仿真，通過靜力學(xué)實驗驗證了理論模型的可行性，在此基礎(chǔ)上分析了軟體機(jī)器人在整個運動周期全過程中的姿態(tài)變化、行進(jìn)運動規(guī)律情況.通過對所得結(jié)果進(jìn)行詳細(xì)的討論和分析，得到軟體尺蠖機(jī)器人爬行運動的一般性規(guī)律，對于仿生爬行機(jī)器人的設(shè)計和運動分析具有一定的指導(dǎo)意義和價值.

1 軟體尺蠖爬行機(jī)器人的簡化力學(xué)模型

1.1 運動過程的描述

在初始狀態(tài)下對軟體尺蠖機(jī)器人充氣，可使其獲得一定大小的初始曲率，如圖1所示.本文利用具有初始曲率的曲梁表示機(jī)器人的彎曲構(gòu)型，從而把整個軟體機(jī)器人簡化為由剛性滑塊和曲梁構(gòu)成的簡化模型，稱之為軟體尺蠖機(jī)器人的滑塊-曲梁模型[8]，如圖2所示.圖中曲梁的線密度為ρ，滑塊的質(zhì)量為ms，重力加速度為g，假設(shè)水平向右為E1方向，豎直向上為E2方向，建立整體坐標(biāo)系OE1E2，原點O位于初始條件下滑塊與曲梁的連接處，且假設(shè)曲梁中線上任意一點沿E1方向和E2方向的坐標(biāo)分別為x、y.

圖1 軟體尺蠖機(jī)器人Fig.1 Soft inchworm-like robot

圖2 滑塊-曲梁模型Fig.2 Slider-curved beam model

圖3 切向量r′(s)Fig.3 Tangential vector r′(s)

設(shè)曲梁中軸線上某一點距離其最左端處的弧長為s，曲梁總長度為l，該點處曲梁中軸線的切向量為r′(s)，如圖3所示.假設(shè)曲梁上處于0和s之間的任意一點的弧長為ξ，則切向量相對E1方向的姿態(tài)角θ為s的函數(shù)，其中，

(1)

設(shè)曲梁的初始曲率為κ0.根據(jù) O’Reilly 的研究發(fā)現(xiàn)，初始曲率與氣腔內(nèi)的氣壓成正比[21]，因此可通過控制κ0的變化對軟體機(jī)器人的運動進(jìn)行控制.本文采用的κ0具有如下形式：

(2)

由于未知量的個數(shù)以及邊界條件的不同，將軟體尺蠖機(jī)器人在整個運動周期的運動狀況劃分為3個不同的階段，如圖4所示，圖中Δ為位移大小.

根據(jù)軟體機(jī)器人的實際運動狀況，其在一個周期內(nèi)的3個不同階段的大致特征為：

(1) 階段C：機(jī)器人部分黏附于地面，隨著滑塊曲梁的初始曲率幅值從0逐漸開始增大，滑塊逐漸向右滑移，曲梁和地面的黏附部分長度也逐漸減少，即曲梁與地面逐漸剝離直到完全脫離，從而過渡到階段A.

(2) 階段A：機(jī)器人右端和地面點接觸，且該接觸點和滑塊都黏滯不動，隨著曲梁的初始曲率幅值減小，曲梁右端點處所需的摩擦力逐漸增大，該點逐漸產(chǎn)生向右滑移的趨勢.當(dāng)曲梁的初始曲率幅值減小到一定值時，曲梁右端點開始滑移，從而過渡到階段B.

(3) 階段B：機(jī)器人右端仍然和地面點接觸，滑塊仍然處于黏滯狀態(tài)，但是曲梁的右端點會向右滑移，隨著曲梁的初始曲率幅值減小，該點逐漸向右移動直到初始曲率幅值減小到0，此時曲梁右端會開始貼附于地面，回到階段C.

可見，軟體尺蠖機(jī)器人完全可以通過上述3個不同階段之間的過渡和轉(zhuǎn)換，完成整個周期的運動，并向前滑移一段位移，實現(xiàn)連續(xù)向前爬行運動的目的.

1.2 應(yīng)變能和重力勢能

不考慮曲梁模型的軸向變形，假設(shè)曲梁上任意一點的軸向應(yīng)變?yōu)棣牛S向應(yīng)力為σ，設(shè)曲梁截面彎矩為M，相應(yīng)的抗彎剛度為EI(E為彈性模量，I為截面慣性矩).根據(jù)O’Reilly由實驗得出的結(jié)論，曲梁截面彎矩M和曲率的改變量存在線性關(guān)系[21]，此處可表示為M=EI(θ′-κ0).假設(shè)曲梁截面上任意一點相對中軸線上投影點的位置矢量在中軸線法向上的坐標(biāo)分量為yu，以曲梁的體積V、截面積A以及弧長l′為積分自變量，由線彈性理論可以得到應(yīng)變能為

(3)

當(dāng)曲梁右端和地面點接觸時，若以曲梁左端點所在水平線為重力勢能的零點位置，在距離曲梁左端點長度為s處的豎直高度為Y(s)，重力加速度的大小為g，則整個曲梁的重力勢能WG為

(4)

當(dāng)曲梁和地面線接觸時，設(shè)γ為曲梁未黏附段的長度，h為慣性基原點相對于地面的高度，則曲梁上未黏附段的重力勢能WG1和黏附段的重力勢能WG2可分別表示為

(5)

1.3 外力勢能

當(dāng)曲梁右端和地面點接觸時，設(shè)接觸點處支持力的大小為N1，摩擦力的大小為Ff1，并以沿著曲梁在左端點處為勢能零點，設(shè)FL=[Ff1N1]T，r(l)=[x(l)y(l)]T，此時由曲梁右端力產(chǎn)生的外力勢能為

(6)

WF1=-n(γ)·r(γ).

圖5 曲梁和地面線接觸時受外載荷的情況Fig.5 External force on curved beam when in line contact with the ground

1.4 控制方程的建立

當(dāng)曲梁和地面點接觸時，結(jié)合曲梁的應(yīng)變能、重力勢能以及外力勢能，可得到點接觸階段中曲梁的總勢能V1為

FL·r(l)=

(7)

當(dāng)曲梁和地面線接觸時，曲梁的未黏附段的勢能V2為

(8)

(9)

得到該情況下的控制方程的具體形式：

(Ff1sinθ-N1cosθ)=0,s∈[0,l]

(10)

(n1(γ)sinθ-n2(γ)cosθ)=0,s∈[0,γ]

(11)

1.5 邊界條件的建立

對于階段A，曲梁和地面點接觸，由于曲梁左端和剛性滑塊固接，即θ(0)=0.由于曲梁右端為自由端，并且曲梁右端截面處也無外力矩，則此處的彎矩應(yīng)該為0，可知M(l)=EI(θ′(l)-κ0(l))=0，即θ′(l)=κ0(l).

曲梁的左端點與滑塊固連，以該點為勢能零點，該點相對地面的高度為-h，且h>0，從而得到系統(tǒng)在豎直方向的幾何邊界條件：

在曲梁右端點黏滯的階段A中，假設(shè)了曲梁右端點相對于滑塊黏滯，可設(shè)兩者間距離的大小始終保持為d，則有水平方向的幾何邊界條件：

綜上，可聯(lián)立曲梁右端和地面點接觸時的控制方程以及邊界條件，階段A用于求解的滿足給定邊界條件的常微分方程組為

(12)

對于階段B，曲梁右端點滑移.設(shè)機(jī)器人與地面接觸點的靜摩擦因數(shù)和動摩擦因數(shù)分別為μs和μk，則可設(shè)曲梁右端處的水平速度為ν(l)>0，階段B在水平方向的邊界條件可寫為

此外，在準(zhǔn)靜態(tài)條件下，設(shè)滑塊和曲梁的相互作用力沿E1和E2方向分量為n1和n2，對階段B中整個滑塊-曲梁系統(tǒng)進(jìn)行靜平衡分析，如圖6所示，可以得到相應(yīng)的靜力平衡方程N(yùn)2=msg+ρgl-N1,Ff1=Ff2.階段B用于求解的滿足給定邊界條件的常微分方程組為：

(13)

如果|Ff2|>μsN2，曲梁左端點也發(fā)生滑移，F(xiàn)f2=-μkN2.

圖6 階段B中滑塊-曲梁系統(tǒng)的靜平衡分析Fig.6 Static equilibrium analysis of slider-curved beam system in phase B

對于階段C，在曲梁和地面線接觸的情況下，由于其在s=γ處的未黏附段右端點即為曲梁和地面接觸的起始點，則θ(γ)=0.此外，設(shè)滑塊和曲梁的相互作用力沿E1和E2方向分量為n1(0)和n2(0)，根據(jù)靜平衡條件先對滑塊進(jìn)行力平衡分析，如圖7所示，得到如下關(guān)系式：

(14)

再對曲梁上未黏附段進(jìn)行力平衡分析，如圖7所示：

(15)

圖7 階段C中滑塊-曲梁系統(tǒng)的靜平衡分析Fig.7 Static equilibrium analysis of slider-curved beam system in phase C

由于滑塊在階段C情況下滑動，設(shè)滑塊處的水平速度為νh，動摩擦因數(shù)為μk，結(jié)合靜平衡方程得出的表達(dá)式，從而得到水平方向的邊界條件：

n1(γ)=μkN2sign(νh)=

μk(msg-n2(0))sign(νh)=

μk(msg+ρgγ-n2(γ))sign(νh)

(16)

式中：sign(νh)為符號函數(shù)，當(dāng)νh≥0，sign(νh)=1；當(dāng)νh=0，sign(νh)=0；當(dāng)νh<0，sign(νh)=1.

當(dāng)曲梁和地面線接觸時，還需要考慮由黏附作用產(chǎn)生的黏附勢能.設(shè)曲梁黏附段上單位長度的黏附勢能為-ω，曲梁未黏附段的黏附邊界條件為[17]

(17)

結(jié)合階段C的曲梁未黏附段的邊界條件和控制方程，則得到如下滿足給定邊界條件的常微分方程組：

(18)

2 軟體尺蠖爬行機(jī)器人控制方程的仿真分析

2.1 方程組的離散化和無量綱化

在曲梁和地面點接觸時，將連續(xù)的曲梁劃分為n段微段，如圖8所示，設(shè)其中每一微段的長度為ds，則nds=l.設(shè)第i段微段的左端點轉(zhuǎn)角為θi，則si=s(θi)=(i-1)ds.

圖8 曲梁和地面點接觸時曲梁的離散構(gòu)型Fig.8 Discrete configuration of curved beam when in point contact with the ground

圖9 曲梁和地面線接觸時曲梁的離散構(gòu)型Fig.9 Discrete configuration of curved beam when in line contact with the ground

假設(shè)曲梁右端末尾還存在第n+1段微段，該段長度仍為ds，且不納入曲梁的總長之中.則此時方程組的未知量共為n+3個，分別為θ1,θ2,…,θn,θn+1,N1,Ff1.

在曲梁和地面線接觸時，由于曲梁黏附段的構(gòu)型已知，通過離散分析曲梁的未黏附段，完成在微段數(shù)目不變，微段長度改變的條件下曲梁構(gòu)型的求解.將連續(xù)的曲梁未黏附段劃分為m段微段，每一微段的長度為ds，使得mds=γ，并且其他假設(shè)和點接觸階段相同.此時方程組的未知量個數(shù)共為m+4個，分別為θ1,θ2, …,θm,θm+1,n2(γ),n1(γ),γ.

本文中對所有的物理量都采取了無量綱化的處理，并將無量綱化的物理量做了新的符號規(guī)定，如表1所示.表中L′、F和T分別表示長度、力和時間的量綱.

從而可以寫出階段A最終的離散化、無量綱化的非線性代數(shù)方程組：

fn+2=θ1=0

同理，要寫出階段B的離散化、無量綱化的非線性代數(shù)方程組，只需采用如下fn+1替換上式中的水平方向邊界條件，即fn+1=Ff1+μkN1=0.

還可以寫出階段C最終的離散化、無量綱化的非線性代數(shù)方程組：

表1 各物理量的無量綱化Tab.1 Dimensionalization of physical quantities

fm+1=θ1=0，fm+2=θm+1=0

2.2 方程組的求解方法

此處將對牛頓迭代法的求解思路以曲梁和地面點接觸的情況為例進(jìn)行說明，而曲梁和地面線接觸情況的求解過程與其一致.設(shè)方程組在第k次迭代時為列向量F(k)，此時由上一次迭代求得的未知量列向量為X(k)，求導(dǎo)結(jié)果構(gòu)成維數(shù)為n+3的方陣即雅可比矩陣F′(k)，并設(shè)未知量迭代矩陣為H(k)：

H(k)=-F′(k)-1F(k)
X(k+1)=X(k)+H(k)

2.3 充氣變形的靜力學(xué)實驗

為了驗證p和κ0之間正比關(guān)系的有效性，本文進(jìn)行了四氣腔軟體機(jī)器人充氣變形的靜力學(xué)實驗.將機(jī)器人的一端固定，研究其受不同氣壓載荷作用下的變形狀態(tài)，如圖10所示.

圖10 受不同氣壓載荷作用的軟體制動器的變形狀態(tài)Fig.10 Deformed states of soft actuator subject to different air pressure loads

圖11 實驗和仿真得到的軟體制動器變形狀態(tài)的對比Fig.11 Comparisons of deformation state of soft actuator obtained by experiment and simulation

2.4 機(jī)器人運動姿態(tài)的準(zhǔn)靜態(tài)實驗

為了驗證所建立的方程組求解結(jié)果的準(zhǔn)確性，本文還通過準(zhǔn)靜態(tài)實驗研究了長度為44 mm的單臂軟體機(jī)器人在爬行過程中某一個時刻的運動姿態(tài)，如圖12所示.單臂軟體機(jī)器人的特點是左端無姿態(tài)約束，滑塊的等效質(zhì)量為0.

圖12 實驗過程中的軟體機(jī)器人Fig.12 Soft robot in experiment

在實驗過程中，隨著充放氣不斷交替變化，軟體爬行機(jī)器人逐漸向前運動，在不同的時刻產(chǎn)生了相應(yīng)的變形形狀曲線，其中在時間t=2 s處軟體機(jī)器人的形狀如圖13所示.

圖13 t=2 s時軟體機(jī)器人形狀的實驗和仿真結(jié)果Fig.13 Experimental and simulation results of soft robot shape at t=2 s

2.5 點接觸階段的計算結(jié)果分析

圖14 階段A在上的構(gòu)型圖Fig.14 Configuration diagram in phase A at

圖15 階段A中曲梁右端外力和的關(guān)系曲線Fig.15 External force of right end of curved beam versus in phase A

圖16 階段B在上的構(gòu)型圖Fig.16 Configuration diagram in phase B at

圖17 階段B中滑塊-曲梁系統(tǒng)外力絕對值和的關(guān)系曲線Fig.17 External force of slider-curved beam system versus in phase B

圖18 階段B中和的關(guān)系曲線 versus in phase B

2.6 線接觸階段的計算結(jié)果分析

圖19 階段C在上的構(gòu)型圖Fig.19 Configuration diagram in phase C at

圖20 滑塊的位移和的關(guān)系曲線Fig.20 Slider displacement versus in phase C

圖21 階段C中滑塊處外力絕對值和的關(guān)系曲線Fig.21 External force at the slider versus in phase C

2.7 軟體尺蠖機(jī)器人的運動周期分析

圖22 曲梁在整個運動周期內(nèi)的構(gòu)型變化Fig.22 Configuration changes of curved beams during the entire motion cycle

此外，若以階段C為起始階段，則也可以獲得初始時曲梁左右兩端點相對的水平位置距離.根據(jù)階段C的計算結(jié)果，曲梁的右端點在初始時的水平投影長度為 0.980 3l.而由階段B的計算結(jié)果，曲梁的水平投影長度在整個運動周期結(jié)束時為 0.990 5l，從而得知曲梁左右端點的水平距離在運動周期前后近似相等，意味著整個軟體機(jī)器人在完成了1個周期的運動后能大致回到運動前的構(gòu)型和狀態(tài)，使得每個周期的運動過程近似一致，從而保證了本文的所有計算結(jié)果和分析結(jié)論對機(jī)器人在任意周期運動內(nèi)的通用性、一致性.

3 結(jié)語

針對現(xiàn)有軟體爬行機(jī)器人的模型在模擬整體運動時的不連貫性、計算效率較低以及在整個周期內(nèi)未能反映連續(xù)運動變化的情況，本文對軟體尺蠖爬行機(jī)器人建立了滑塊-曲梁的簡化力學(xué)模型，推導(dǎo)了用于數(shù)值求解的無量綱離散化的非線性代數(shù)方程組，對該機(jī)器人的整體運動狀況進(jìn)行了準(zhǔn)靜態(tài)建模和仿真分析.