多腔體氣動軟體致動器的建模與仿真

徐齊平， 劉錦陽

(上海交通大學 船舶海洋與建筑工程學院， 上海 200240)

與通常用于特定環(huán)境下重復作業(yè)的傳統(tǒng)剛性機器人不同，軟體機器人因其變形量大、安全性高、柔順性好、自適應性強、與環(huán)境相容性好等優(yōu)點而逐漸涌現(xiàn)并受到很高的關注度.軟體機器人是一種以自然界的柔軟體生物為原型的新型連續(xù)體仿生機器人，是機器人研究領域的拓展，其利用諸如硅膠、水凝膠、橡膠、電活性聚合物等柔軟性材料制成[1-3]，可以實現(xiàn)抓取、爬行、游動、彈跳等復雜運動[4-6].軟體致動器可以廣泛用于實現(xiàn)各種復雜精細動作，其由氣壓、加壓液體、熱、電場、化學反應等驅動.其中，氣動軟體致動器因輸入簡單、質輕、材料成本低和制造效率高等優(yōu)勢而具有廣闊的應用前景[7-8].

在大量氣動軟體致動器的設計中，纖維增強型致動器[9]和氣動網(wǎng)絡型致動器[10]是兩種十分受青睞的類型.纖維增強型致動器主要由具有整體腔室的彈性體、螺旋排列的纖維和不可延伸層構成.通過改變纖維角度，可以實現(xiàn)伸縮、擴展和彎曲等動作.氣動網(wǎng)絡型致動器由多個相互連接的、網(wǎng)絡狀的腔室構成[11]，通過在腔室內部施加氣壓，使致動器充氣膨脹，實現(xiàn)預定的彎曲變形，執(zhí)行各種基本動作.

氣動軟體致動器除了能夠完成平面內的基本運動形式之外，還可以實現(xiàn)空間內的彎曲、扭轉、螺旋等耦合運動以便進行靈活操作[12-13].哈佛大學的Whitesides及其團隊[4, 11]設計制作了大量氣動軟體致動器并進行一系列實驗，這些致動器可以執(zhí)行許多動作，從而完成移動、爬行、游動、彈跳等復雜運動.Hirai等[14]提出一種分析模型并使用可變形圓柱型致動器以展現(xiàn)這些運動形式，并利用由多個單彈性管組成的一組氣動致動器實現(xiàn)復合運動.Moser等[15]應用一種分析模型，通過平行組合纖維和采用額外的纖維創(chuàng)建了蛇狀致動器.Connolly等[3, 9]詳細分析了纖維增強型致動器的運動如何受到纖維角度的影響，并提出了快速設計的分析模型.Sun等[16]提出了一種通過周向改變氣動軟體致動器的應變約束層的每個步驟，以實現(xiàn)螺旋運動的方法.

雖然國內外學者在氣動軟體致動器的設計方面開展了卓有成效的工作，但是如何對其建立準確的力學模型，并且采用何種方法是行之有效的，這仍是個具有挑戰(zhàn)性的難題.有限元思想提供了一種切實可行的解決途徑，考慮到軟體致動器的材料和幾何非線性，為了降低計算成本、揭示物理現(xiàn)象、滿足精度并提高效率，需要選擇合適的單元類型和非線性本構關系.針對經(jīng)典位移場型的實體單元，高階單元通常比低階單元更精確有效，然而從計算效率來看，高階單元難以應用于實際場合，因為其難以構造且計算更耗時，尤其對于自由度龐大和幾何構型復雜的三維軟體結構更是如此，因此，在保證計算精度的前提下，采用低階單元更加高效.此外，在計算單元內力陣和剛度陣時，如果采用完全積分法，會出現(xiàn)體積鎖定現(xiàn)象.為了避免體積鎖定問題且保證計算效率，采用選擇性縮減積分技術以解決該問題.綜上所述，采用縮減積分法的低階實體單元對氣動軟體致動器進行建模是一種合理的解決方案.

上述軟體致動器的共同之處在于它們具有氣動網(wǎng)絡構造，不可壓縮超彈性體嵌入到一系列相互連接的氣腔中.之前大部分研究主要集中在設計制作、實驗及采用桿、梁等簡單模型進行整體變形分析上[10, 17-18].相對而言，在考慮到相鄰兩氣腔互相接觸的情況下，對整個氣動軟體致動器進行建模與分析則很少涉及，然而這正是我們感興趣和關注的地方.盡管以前的實驗研究已經(jīng)明顯觀察到軟體致動器相鄰兩氣腔彼此接觸的現(xiàn)象，但仍缺乏關于相鄰兩氣腔之間如何接觸的研究.迄今為止，在考慮接觸的情況下如何使致動器實現(xiàn)足夠彎曲變形、建立準確的力學模型并進行整體構型和應力分布分析尚是一個未得到有效解決的難題.

本文針對多腔體氣動軟體致動器相鄰兩氣腔之間的互相接觸問題，考慮到結構的復雜性、材料和幾何非線性，采用多點接觸的面-面接觸方法，對其進行建模與分析，解決了不考慮接觸時相鄰兩氣腔之間的穿透問題.仿真結果表明，接觸對致動器的大變形影響顯著.最后，對軟體四爪機構抓取圓柱體的整個過程進行了準靜態(tài)研究，詳細闡述了整體構型和Mises應力變化，這有助于指導軟體致動器的優(yōu)化設計及力學特性研究.

1 氣動軟體致動器

圖1所示為常見的內嵌4個氣腔的氣動軟體致動器，制作過程詳見文獻[19].為了清楚地顯示致動器的內部結構，考慮到結構的對稱性及簡化計算，首先以其一半作為研究對象，內部構造如圖2所示.圖1和圖2中：L=44 mm，W=20 mm，W1=10 mm，H=13 mm，t=4 mm，t1=2 mm，t2=10 mm，t3=6 mm，t4=8 mm.

圖1 氣動軟體致動器

圖2 氣動軟體致動器的一半

2 建模及非線性力陣

采用能夠消除體積鎖定的、具有選擇性縮減積分的三維實體單元進行建模，該單元所有場變量依賴于總體坐標x、y、z.實體單元在空間中具有任意的幾何構型、材料屬性和邊界條件，其形狀可以為平面或曲面的四面體或六面體.實體單元的每個節(jié)點具有3個移動自由度，從而可以在三維空間中實現(xiàn)任意變形.

對于三維實體單元，容易構造出形函數(shù)并進行數(shù)值積分，實體單元上任意一點的位移矢量在總體坐標系下的坐標陣為

u(ξ,η,ζ,t)=N(ξ,η,ζ)qe(t)

(1)

(2)

形函數(shù)對自然坐標與對笛卡爾坐標導數(shù)的轉換關系為

(3)

式中：J為雅克比矩陣，可表示為

(4)

2.1 單元彈性力陣及其導數(shù)陣

設r為單元上任意一點的位置矢量，可表示為r=[xyz]T+u，定義G=[?r/?x?r/?y?r/?z]為變形梯度張量，G為變形梯度張量的行列式.右柯西-格林變形張量C=GTG，tr(C)為C的跡.

引入體積能量罰函數(shù)，超彈性不可壓縮材料的非線性Yeoh模型的單元應變能為[20]

(5)

采用選擇性縮減積分技術以消除體積鎖定問題，則

(6)

式中：GC=G(ξC,ηC,ζC)，C為單元的中心點；κ為不可壓縮常數(shù)；C10,C20,C30為材料常數(shù)，通過拉伸試驗確定[20].I1為右柯西-格林變形張量的第一不變量，可表示為

(7)

將式(6)對單元坐標陣qe求導，可得單元彈性力陣為

(8)

將式(6)代入式(8)，則單元彈性力陣可表示為

(9)

單元彈性力陣的導數(shù)陣為Qek對qe的導數(shù)，即

(10)

2.2 單元壓力陣及其導數(shù)陣

為了反映相鄰兩氣腔之間的接觸情況，首先需要求出作用于氣動軟體致動器的單元壓力陣.氣壓所做虛功可表示為

(11)

(12)

(13)

2.3 單元接觸力陣及其導數(shù)陣

之前對接觸的研究大多集中于硬質材料，而對軟體致動器的研究也只專注于整體變形，沒有考慮到軟體致動器相鄰兩氣腔之間的互相接觸，且整體構型和應力分析也很少涉及.本文不計相鄰兩氣腔之間的摩擦，假設給定的接觸點總是唯一地投射到任意表面，針對軟材料的軟接觸問題，引入多點接觸的面-面接觸方法研究相鄰兩氣腔之間的接觸情況.多點接觸的面-面接觸方法是指一個從單元的多個高斯點投影到相應主單元上并與主單元的多個高斯點接觸的方法.為了避免穿透并提高精度，采用4×4個高斯點，當接觸剛度系數(shù)逐漸增大到一定值時，結果收斂.同時為了避免在每個時間步對任意兩個單元進行窮舉幾何搜索，接觸檢測可分為全局檢測和局部檢測兩個階段，全局檢測階段消除了大部分組合，并且選擇一系列單元對進行局部檢測，文中采用的全局搜索法是層次包圍盒算法[21].

接觸分析中的接觸對由主面和從面構成，如圖3所示.設v為從面上的從節(jié)點，m為從節(jié)點v到主面上的最近投影點[22].R是從節(jié)點v的全局位置矢量，ρ為最近投影點m的全局位置矢量，相對位置矢量為rmv=R-ρ，nm為接觸主面的法向量，

ρ=ρ(ξm,ηm)

(14)

圖3 最近投影點

首先，構造如下函數(shù)：

(15)

然后，求最近投影點主要是求解式(16)的極小值問題

‖R-ρ‖=(R-ρ)·(R-ρ)→min

(16)

通過求解式(17)可得式(16)的解

(17)

rmv在接觸力的單位法線方向nm上投影為

(18)

顯然，當d>0時，從面上任意一點和主面上投影點受到的接觸應力pv=pm=0；而當d<0時，穿透量δ=|d|.由罰函數(shù)法，從面上任意一點和主面上的投影點受到的接觸應力為

pv=kδnm,pm=-kδnm

(19)

式中：k為接觸剛度系數(shù).從面和主面的單元接觸力陣為

(20)

(21)

2.4 靜平衡方程

針對圖2所示一端夾持的多腔體氣動軟體致動器，可得具有約束的靜平衡方程

(22)

式中：λ為約束方程對應的拉格朗日乘子列陣；Qk、Qp和Qc分別為系統(tǒng)總體廣義彈性力陣、廣義壓力陣和廣義接觸力陣；Φq為雅克比矩陣；Φ=0為氣動軟體致動器受到固定端約束的約束方程.

(23)

式中：nk為實體單元數(shù)；np為受氣壓作用的單元數(shù)；nc為接觸單元對數(shù).

綜合使用Newton-Raphson迭代法和高斯積分法，對式(22)進行求解，在迭代計算中，第i步和第i+1步具有以下關系:

(24)

(25)

式中：Kk、Kp和Kc分別為Qk、Qp和Qk對總體廣義坐標的導數(shù)陣，

(26)

3 仿真分析、實驗驗證與模型比較

3.1 單個軟體致動器的仿真分析

在考慮到相鄰兩氣腔之間的互相接觸情況下，通過求解式(22)，采用多點接觸的面-面接觸方法，對圖2所示一端夾持的多腔體氣動軟體致動器在受氣壓作用下的變形進行仿真分析，這里忽略重力的影響.

本數(shù)值算例的氣動軟體致動器的外形尺寸如圖2所示，硅膠材料的各項參數(shù)如下：密度ρ=1 100 kg/m3，泊松比ν=0.499 5，Yeoh模型的材料系數(shù)C10=0.271 2 MPa，C20=0.030 53 MPa，C30=-4.013×10-4MPa，κ=1 000 MPa，接觸剛度系數(shù)k=4 N/mm3.

為了清楚地觀察到軟體致動器氣腔的內部變形和整體構型變化，按照1 mm的間隔對該致動器劃分網(wǎng)格，當氣壓p=35 kPa時，在考慮接觸與不考慮接觸兩種情況下，借助于MATLAB數(shù)值軟件進行仿真對比分析，半個軟體致動器的變形后構型分別如圖4和5所示.圖4為考慮接觸時半個氣動軟體致動器的變形后構型，圖5為不考慮接觸時半個氣動軟體致動器的變形后構型.從圖5可以看出在不考慮接觸的情況下，相鄰兩氣腔相互穿透，即圖中黑色虛線包圍的區(qū)域，這顯然不符合實際情況；而在考慮接觸的情況下，相鄰兩氣腔互相接觸，無穿透現(xiàn)象，這可以從圖4明顯地反映出.此外，為了清晰地顯示軟體致動器的相鄰兩氣腔之間的應力分布，在考慮接觸時對該致動器的Mises應力分布進行計算，Mises應力為

σs={[(σx-σy)2+(σy-σz)2+(σz-σx)2+

(27)

式中：σx、σy、σz分別為第一、第二、第三主應力；τxy、τyz、τzx為切應力.從而得出變形之后構型的 Mises 應力分布情況，如圖6所示.

根據(jù)結構的對稱性且考慮到相鄰兩氣腔之間的互相接觸作用，采用MATLAB軟件對圖1所示的整個多腔體氣動軟體致動器進行仿真，可得在氣壓p=35 kPa時變形后構型和Mises應力分布，分別如圖7和圖8所示.圖7為考慮接觸時整個氣動軟體致動器的變形后構型，圖8為考慮接觸時整個氣動軟體致動器的Mises應力分布情況.從圖6和圖8可以清晰地發(fā)現(xiàn)應力最大值位于相鄰兩氣腔之間的接觸區(qū)域與連接相鄰兩氣腔的根部區(qū)域，這對氣腔壁承受氣壓是個極大的考驗，計算和分析這些區(qū)域的應力變化有助于軟體致動器的優(yōu)化設計和力學性能分析.

圖4 考慮接觸時半個軟體致動器的變形后構型

圖5 不考慮接觸時半個軟體致動器的變形后構型

圖6 考慮接觸時半個軟體致動器的Mises應力分布

圖7 考慮接觸時整個軟體致動器的變形后構型

圖8 考慮接觸時整個軟體致動器的Mises應力分布

3.2 實驗驗證與模型比較

為了驗證所提出的多點接觸的面-面接觸模型對多腔體氣動軟體致動器建模的有效性和準確性，并測試該致動器受不同氣壓作用下的彎曲性能，對受不同氣壓載荷作用下的軟體致動器進行相應的實驗研究，得到其變形狀態(tài)，如圖9所示.從圖9中可以清晰地看出，隨著氣壓逐漸增大，軟體致動器的變形逐漸變大，即固有曲率輪廓逐漸增大.

文獻[10]采用桿模型對氣動軟體致動器進行研究，該模型基于彈性體的歐拉理論，即假定致動器的底層為彈性桿，是一種未強調相鄰兩氣腔互相接觸的簡化模型.

圖10給出了利用兩種不同模型得到的加壓氣動軟體致動器的變形形狀.圖中氣壓p依次取值為0，5，15，25，35 kPa.圖中實線、虛線及離散點分別對應于利用多點接觸的面-面接觸模型獲得的結果、桿模型獲得的結果及實驗數(shù)據(jù)，其中，實驗數(shù)據(jù)可以通過測量圖9中致動器上的各個黑色標記點的位置而獲得.從圖10中可以看出，軟體致動器的變形隨著氣壓增加而變大，當氣壓小于p=15 kPa時，即在小變形范圍內，采用上述兩種模型得到的仿真結果與實驗數(shù)據(jù)基本一致.當氣壓大于p=25 kPa時，隨著彈性變形逐漸增大，使用多點接觸的面-面接觸模型得到的結果與實驗數(shù)據(jù)吻合，而使用桿模型得到的結果逐漸偏離實驗數(shù)據(jù)，誤差逐漸變大.這表明與桿模型相比，在大變形階段，利用多點接觸的面-面接觸模型能夠有效提高模擬大變形彎曲的精度.

圖9 受不同氣壓載荷作用的軟體致動器的變形狀態(tài)

圖10 利用兩種不同模型得到的軟體致動器的變形形狀

3.3 軟體四爪機構的仿真分析

利用3D打印技術制造一個裝載氣動軟體致動器的平臺，如圖11(a)所示，具體尺寸為Lg=96 mm，l=34.5 mm，h=3 mm，h1=10 mm，h2=8 mm，w=27 mm，l1=23 mm，w1=20 mm，l2=20 mm，w2=17 mm，m=5 mm，其余各項參數(shù)與單個致動器的參數(shù)相同.

將圖1所示的4個相同的氣動軟體致動器裝載到圖11(a)所示的平臺中，組裝成圖11(b)所示的軟體四爪機構.利用圖11(b)所示的軟體四爪機構對抓取圓柱體的整個過程進行準靜態(tài)分析.軟體四爪機構與圓柱體的位置關系如圖12所示.

圖11 平臺和軟體四爪機構

圖12 抓取圓柱體的軟體四爪機構

在分析了上述單個軟體致動器之后，為了驗證該致動器的性能，在考慮接觸的情況下，采用多點接觸的面-面接觸模型對該軟體四爪機構抓取圓柱體的整個過程進行準靜態(tài)仿真分析.圓柱體的密度為 3 388 kg/m3，直徑為32 mm，高為36 mm.初始時，軟體四爪機構與圓柱體之間互不接觸，圓柱體在支撐平臺上處于靜止狀態(tài).隨后，該四爪機構的4個致動器開始均勻充氣直到在t=1 s時達到35 kPa為止，之后保持該氣壓值不變，氣壓的加載曲線如圖13所示.在t=2 s時，將支撐平臺撤掉，在摩擦力作用下，圓柱體保持平衡.

考慮到結構的對稱性，取四分之一圓柱體進行研究，軟體四爪機構與圓柱體之間為點接觸，如圖14所示，這是機器人抓爪或人類手指抓取類型的最常見方式.

作用于四分之一圓柱體的力由帶箭頭的紅色實線表示，法向接觸力FN和摩擦力Ff必須滿足以下靜力學準則：

Ff≤μFN

(28)

式中：μ為靜摩擦系數(shù)，取值為0.25.摩擦力Ff的大小為四分之一圓柱體的重力.

考慮摩擦力作用和位移約束，可得單個氣動軟體致動器的靜平衡方程：

(29)

式中：Qg和Qf分別為系統(tǒng)總體廣義重力陣和廣義摩擦力陣；ΦN=0為軟體致動器與圓柱體接觸點的位移約束方程，λN為位移約束方程ΦN=0對應的拉格朗日乘子列陣.根據(jù)式(29)，可以計算出幾個時刻的法向接觸力FN和摩擦力Ff的變化情況，如表1所示.

為了更具體和清晰地反映出軟體四爪機構抓取圓柱體的整個準靜態(tài)過程，繪制出其不同時刻的構型，如圖15所示.

結合圖13、圖15和表1可知，在初始時刻t=0 s時，整個機構處于靜止狀態(tài)；在t=0.5 s時，軟體四爪機構的4個致動器均勻充氣膨脹產(chǎn)生彎曲變形而與圓柱體接觸，法向接觸力FN大于0，此時圓柱體還處在支撐平臺上保持不動；在t=1 s時，充氣氣壓達到最大值而使4個致動器與圓柱體緊緊接觸，F(xiàn)N增大；直至t=2 s為止，圓柱體的重力與支撐力平衡而仍保持不動，摩擦力Ff的大小為0，F(xiàn)N保持不變.在t=3 s時，達到新平衡狀態(tài)，F(xiàn)N和Ff滿足平衡條件.

表1 法向接觸力FN和摩擦力Ff

圖13 氣壓載荷-時間曲線

圖14 四分之一圓柱體

圖15 在不同時刻軟體四爪機構抓取圓柱體的構型

此外，為了顯示軟體四爪機構抓取圓柱體時的應力分布，根據(jù)式(27)，可得未變形構型和變形后構型的Mises應力分布情況，如圖16所示.由圖16可知，在初始時刻t=0 s時，整個機構處于靜止狀態(tài)，其在未變形構型時的Mises應力分布為0；在t=3 s時，圓柱體在摩擦力作用下保持平衡，圓柱體與軟體四爪機構之間的壓力必須足夠大，以滿足摩擦力小于極限摩擦力.從顏色條可以看出相鄰兩氣腔之間的接觸區(qū)域Mises應力達到最大值，這個區(qū)域的應力變化情況特別值得注意，計算和分析這個區(qū)域的應力，有助于在不破壞氣腔壁的前提下進一步優(yōu)化設計軟體致動器并分析力學特性.

圖16 軟體四爪機構抓取圓柱體的Mises應力分布

4 結論

針對現(xiàn)有梁、桿等簡單模型對氣動軟體致動器進行整體變形研究的局限性，且與剛性材料的硬接觸不同，本文對多腔體氣動軟體致動器的軟接觸問題，建立了準確的力學模型，進行了整體構型和應力分布分析.綜合考慮相鄰兩氣腔之間的互相接觸作用以及氣腔結構的復雜性、材料和幾何非線性，采用多點接觸的面-面接觸方法，解決了不考慮接觸時相鄰兩氣腔之間的相互穿透問題.通過相應的實驗，驗證了本模型的有效性，與桿模型相比，有效地提高了模擬大變形彎曲的精度.

本力學模型能夠模擬在不同氣壓作用下軟體致動器的彎曲變形，適用于研究各類致動器的大變形問題，不僅可以分析軟體致動器的整體變形，而且能夠描述氣腔結構內部的應力分布情況.在考慮到相鄰兩氣腔之間的互相接觸時，能夠準確地模擬其變形后構型；而在不考慮互相接觸時，相鄰兩氣腔相互穿透，這明顯不符合實際情況.

軟體四爪機構在變形過程中的整體構型和Mises應力變化揭示，相鄰兩氣腔之間的接觸對致動器的變形影響顯著.相鄰兩氣腔之間的接觸區(qū)域為應力集中區(qū)域，如何在保證材料不屈服的條件下增大氣壓以實現(xiàn)成功抓取，是未來軟體致動器設計需要考慮的一個重要問題.本文提出的力學模型對于氣動軟體致動器的力學特性分析、優(yōu)化設計以及驅動性能研究具有一定的理論參考價值.