機器人針刺成型異形預制體的路徑規劃①

邵建娜,張雅秀,曾欣怡,蔣 云*,李嘉祿

(1.江南大學 紡織科學與工程學院,無錫 214122;2.天津工業大學,天津 300387)

0 引言

纖維預制體是采用定型劑或紡織方法將增強纖維預成型為特定形狀的纖維復合材料增強體。針刺作為一種工藝簡單、低成本的織物增強成型制備技術而被航空航天領域廣泛應用,使用該技術制備的纖維復合材料比強度、比剛度高,抗疲勞性能良好[1-2]。三維針刺成型技術的原理是使用帶倒鉤的刺針反復刺入纖維網,當刺針的倒鉤進入纖維網時,由于纖維之間摩擦力的作用,使得原本有序的纖維交叉纏繞在一起,經過大量反復的針刺,最終得到具有良好力學特性的針刺復合材料預制體。

三維針刺成型技術適用于制備平板預制體、圓環預制體、圓柱和圓錐預制體,以及具有復雜曲面形狀的預制體。在過去的幾十年時間里,國內外的研究單位及科研工作者對三維針刺預制體的成型裝備開展了廣泛的研究[3]。OLRY[4]發明了一種可以制備高厚三維針刺預制體的針刺機,該設備在厚度方向具有較大的工作范圍,可以用于制備不同厚度的平板預制體。為了進一步擴展針刺機的適用性,RENAUD等[5]公開了一種環形纖維預成型構件的制造方法,該方法可以實現環向纖維的自動鋪放和針刺,自動化程度更高。航空航天用構件通常具有較為復雜的曲面形狀。OLRY和RENAUD發明的平板和圓環預制體針刺機不能很好的滿足近凈成型的要求,限制了三維針刺預制體及其復合材料的開發和應用[3]。OLRY[6]在之前的基礎上發明了用于異形纖維預制體的針刺設備,由四自由度機械臂拖動針刺頭完成針刺,針刺頭角度可在0～90°范圍內進行調整,特別適合橢球、圓錐等特定類型的非封閉殼體制品的針刺工作[7]。后續,FRANCOIS等[8]設計了一款數控復合材料纖維增強針刺設備,設備采用龍門架式結構并配有數控系統,工作范圍更大,自動化程度更高,可進行平面、曲面的針刺任務,缺點是針刺角度不可調節[7]。由于技術封鎖的原因,國外在機器人復合材料三維針刺成型領域應用研究的資料寥寥無幾。國內有天津工業大學的陳小明[3]借助CATIA軟件開發了機器人針刺路徑規劃算法,缺點是該算法依賴其他第三方軟件。童寧[9]對PCA算法進行改進得到一種具有廣泛適應性的噴槍姿態獲取方法。周仁義[10]采用基于Zernike矩的亞像素輪廓提取方法,通過NURBS曲線擬合的方法得到機器人運動路徑。

在機器人快速發展的今天[11-15],本文在此背景下對異形預制體機器人針刺成型路徑規劃進行研究,解決異形變截面預制體變密度針刺,棱邊和大面不同植針數量末端執行器無縫銜接針刺,實現輪廓形狀可變,纖維體積含量的精確控制,推動三維針刺裝備快速發展,助力中國航天事業。

1 異形預制體機器人針刺路徑生成方法

由于自由曲面的靈活性使得異型預制體針刺路徑生成較為困難,主流的等參數法具有較大的局限性。為此,本團隊提出了基于切片法的異型預制體機器人針刺路徑生成方法[7],將預制體模型保存為STL格式并提取其中的三角面幾何信息,采用多個平行切平面與STL模型文件三角面求交得到異型預制體針刺路徑型值點集合,如圖1所示。然后利用交點冗余信息對針刺路徑點集進行優化排序。最后,使用3次NURBS方法對針刺路徑型值點進行擬合[7],從而得到機器人針刺路徑以及路徑表達式。圖中,T1、T2、…、T9為STL文件包含的三角面;i-2、i-1、…、i+2為多個切平面;a、b、…、e為i+2層切平面與三角面的交點,后文統稱為針刺路徑型值點。

圖1 異形預制體機器人針刺路徑生成示意圖Fig.1 Schematic diagram of needling path generation for special-shaped preform robot

1.1 STL模型文件特點

STL文件是計算機圖形應用系統中使用三角面集合去描述物體的一種文件格式,該格式文件可以根據不同的精度要求對模型表進行網格劃分,將模型表面離散為許多三角面的集合,每個三角面記錄著三角面的三個頂點坐標和外法向矢量[9]。不同于其他模型文件(STEP、IGES等),該格式文件只限于幾何信息數據在3D物體的表面,而沒有描述三角面之間的拓撲關系。文本格式的STL文件數據結構如下:

solid film name.stl//模型文件名稱

facet normalnxnynz//法向量的3個分量值

outer loop

vertexx1y1z1//第一個頂點坐標

vertexx2y2z2//第二個頂點坐標

vertexx3y3z3//第三個頂點坐標

endloop

endfacet //結束第一個三角面定義

……//其余三角面定義

end solid //文件結束聲明

1.2 三角面與切平面交點計算

三角面與切平面交點計算的基本算法是當前切平面與一個三角面相交,以Z軸為切片方向進行說明。經過分析,切平面與三角形面片的相交情況存在以下4種(見圖2):

(1)X-Y切平面與三角面的一個頂點相交;

(2)X-Y切平面與三角形面片的兩個頂點相交;

(3)X-Y切平面與三角形面片的一條邊和一個頂點相交;

(4)X-Y切平面與三角形面片的兩條邊相交。

設三角面的三個頂點坐標分別為v1(x1,y1,z1)、v2(x2,y2,z2)、v3(x3,y3,z3),th為切片間距,z0為Z軸方向切片的起始高度。假設當前位于第i個切片層,在Z軸方向高度為

zi=i·th+z0

(1)

當X-Y切平面與小三角形面片的一個頂點相交時,說明該頂點就是在當前切平面上的一個針刺路徑點,如圖2(a)所示,二者相交于v1點,則針刺路徑型值位置信息為

圖2 三角面相交的情況分析Fig.2 Analysis of triangular patch intersection

(2)

當X-Y切平面與小三角形面片的兩個頂點相交時,說明三角面的一條邊與切平面重合,此時該條邊的兩個頂點就是在當前切平面的兩個針刺路徑點,如圖2(b)所示,二者相交于v2和v3點,則針刺路徑型值點位置信息分別為

(3)

當X-Y切平面與小三角形面片的一條邊和一個頂點相交時,此時產生兩個針刺路徑點。其中一個為三角形的一個頂點,二者相交于v3點;另一個位于三角形邊v1v2上,可以使用相似三角形原理求解三角形邊v1v2上的針刺路徑點位置信息,如圖2(c)所示,則針刺路徑點位置信息分別為

(4)

和

(5)

當X-Y切平面與三角形面片的兩個邊相交時,說明該切片層與三角形面片相交產生兩個針刺路徑點,兩個針刺路徑點分別位于三角形面片的兩條邊上[9],如圖2(d)所示,此時兩個針刺路徑點分別為

(6)

和

(7)

1.3 基于交點冗余信息的針刺路徑型值點排序

由于STL模型丟失了原物理模型之間的拓撲關系,三角形的記錄順序雜亂無章,為了使交點坐標有序排列,需要對其重建三角形間的拓撲關系。由于STL模型三角形之間存在共邊原則,因此在切片過程中存在著冗余的交點信息,本文基于此信息利用逐點識別與剔除的思想重建拓撲關系,進而對針刺路徑型值點進行排序[10]。

根據STL文件共邊規則,三角形面片共用一條邊,T1、T2三角形與切平面的交點記錄形式[ab]、[cd][10],分別為面片T1和面片T2的針刺路徑型值點,坐標值b=c,在順時針記錄有序點坐標時可將點坐標視為點坐標的冗余交點。因此,利用三角形與切平面的交點坐標可快速查到其冗余交點坐標,得到三角形的鄰接三角形,從而得到針刺路徑點型值點順序為a、b、d或a、c、d,冗余交點如圖3所示。

圖3 交點冗余信息示意圖Fig.3 Schematic diagram of intersection redundancy information

1.4 基于NURBS曲線擬合的針刺路徑生成

經過1.2和1.3節的運算,將排序后的針刺路徑型值點擬合就形成了針刺路徑,針刺路徑型值點是后續規劃針刺點位置和姿態的基礎數據。由于NURBS方法靈活性和連續性好的優點目前已在CAD/CAM領域獲得廣泛應用,且這種方法已經成為產品數據交換的國際標準。因此,本文采用NURBS方法進行針刺路徑點擬合。一條k次的NURBS曲線可以由以下有理分式表示[7]:

(7)

其中,di(i=0,1,…,n)為控制頂點;ωi(i=0,1,…,n)為權因子,與控制頂點di一一對應;Ni,k(u)是由節點矢量U=[u0,u1,…,un+k+1]通過德布爾-考克斯遞推公式決定的k次B樣條基函數,其遞推公式為[11]

(8)

設1.3節獲取的針刺路徑型值點序列為pi(i=0,1,…,n),由于3次NURBS曲線完全滿足實際應用中的要求,所以本文應用3次NURBS曲線進行針刺路徑型值點擬合,擬合流程圖如圖4所示[7]。

圖4 針刺路徑型值點NURBS曲線擬合流程圖Fig.4 Flow chart of NURBS curve fitting for needling path type value points

(1)型值點參數化

型值點參數化實際上是為每一個型值點pi賦予參數值ui,確定節點矢量U。積累弦長參數化方法是目前公認的最佳參數化法,該方法可以如實地表現出型值點的分布情況,更符合切片法得到的數據點特點,則使用積累弦長參數化方法得到的節點矢量U[12]:

(9)

(2)控制頂點反算

將控制頂點di與權因子ωi重新組合,定義帶權控制頂點Di[ωidi,ωi](i=0,1,…,n),轉換為齊次坐標形式,它們之間的關系為

(10)

對于C2連續的3次NURBS閉曲線,式(10)中包含n+1個方程,因為首末型值點重合p0=pn,不計重復,方程數減少1個,剩余n個。又首末3個控制頂點依次相重Dn=D0,Dn+1=D1,Dn+2=D2,未知控制頂點數少了3個,剩余n個。因此,就可從n個方程構成的線性方程組求解出個n個未知不重復帶權控制頂點[13]。將式(10)改寫為矩陣形式:

(11)

其中,

(12)

對于3次NURBS開曲線,式(10)的n+1個方程不足以解決其中包含的n+3個未知控制頂點,還必須增加兩個通常由邊界條件給定的附加方程[14]。此時求解3次B樣條插值曲線未知控制頂點的線性方程組可寫成如下矩陣形式:

(13)

其中,系數矩陣中首行非零元素a1、b1、c1與右端列陣中矢量e1表示了首端點邊界條件;系數矩陣中末行非零元素an+1、bn+1、cn+1與右端列陣中矢量en+1表示了末端點邊界條件[15]。其余各行元素ai、bi、ci(i=1,2,3,…,n)與C2連續閉曲線的情況相同。

邊界條件通常采用端點切矢條件,切矢條件相當于力學中的梁的端部固支的情況,該條件具有固定的切線方向。采用端點切矢條件時,首、末端滿足下列附加方程[16]:

(14)

(15)

2 針刺點位置和姿態信息確定

在獲得針刺路徑的表達式后,根據針刺工藝對針刺路徑進行等步長插補,進而得到一系列針刺點位置和姿態信息,將生成的一系列針刺點位姿進行逆運動學求解,運動控制器對各個電機進行脈沖分配,驅動機器人TCP點沿著針刺軌跡移動到針刺點完成針刺。

2.1 針刺點位置信息確定

由于插補節點參數集合{u1,u2,…,un}與插補時間集合{t1,t2,tn}相互對應,以插補時間t為自變量,節點參數u為因變量建立函數。將參數u對時間t的函數在ti時刻進行泰勒展開,則有

(16)

記T=ti+1-ti,則在t=ti+1時刻可得

(17)

略去式(17)中高階項,只保留到二階項,記u(ti)=ui,u(ti+1)=ui+1,即可得插補參數u的遞推公式為

(18)

由式(18)可知,插補參數的計算問題轉換為了節點參數u對時間t的求導問題。若機器人末端TCP點始終沿著曲線移動,則理想插補速度Varc(ui)為

(19)

用給定進給速度V(ui)近似代替理想進給速度Vace(u)并將式(19)變形得

(20)

對式(20)求導得

(21)

因此,基于泰勒展開法的插補參數計算公式為

(22)

為了統一描述插補參數計算誤差的大小,引入速度波動率δ作為評價標準:

(23)

速度波動率δ值越小,插補參數越接近理想值。

二階泰勒展開法計算插補參數均要進行復雜的求導運算,不利于算法的實時性,同時算法截斷誤差較大。為此本文提出了基于四階Runge-Kutta的插補參數預估校正法,該方法首先運用Runge-Kutta法對插補參數進行預估計算,之后使用校正公式對得到的插補參數進行修正,使速度波動率保持在限定的范圍內,提高插補參數的計算效率和計算精度,確保針刺點間距一致。

將四階Runge-Kutta法運用到NURBS曲線插補中用以計算插補參數。插補參數計算過程中,時間t為自變量,參數u為因變量,T為時間間隔(插補周期),T=ti+1-ti,已知在ti時刻對應插補參數為ui,給定進給速度V(ui),目標是求解處時刻ti+1的插補參數ui+1,因此建立非線性微分方程:

(24)

由二階泰勒展開法可得

(25)

由四階Runge-Kutta法可得

(26)

其中,

(27)

四階Runge-Kutta法與二階泰勒展開法相比,只需要進行一階導數的運算,避免了二階導數的運算,同時提高了計算精度,但是相較于理想的插補參數還是有一定的偏差,主要原因是用給定的插補進給速度V(ui)替代理想進給速度Varc(ui)。對此以速度波動率作為優化目標設計校正算法對四階Runge-Kutta法計算出的插補參數初值進行修正,提高插補參數的計算精度,保證各個周期插補步長一致。

圖5 基于速度波動率的插補參數校正示意圖Fig.5 Schematic diagram of interpolation parameter correction based on velocity fluctuation rate

本文設計的修正策略如下:

(28)

基于四階Runge-Kutta法的預估校正法插補參數計算流程如圖6所示。

圖6 基于四階Runge-Kutta的預估校正 插補參數計算流程圖Fig.6 Calculation flow chart of predictive correction interpolation parameters based on fourth order Runge Kutta

具體計算步驟如下:

Step1:初始化i=0,ui=0;

Step3:采用本文設計的四階Runge-Kutta法計算式(26),計算初始插補參數ui+1;

Step4:計算出插補參數ui+1處的插補點坐標C(ui+1),得到當前插補參數ui+1插補時的實際進給步長ΔLi;

Step5:根據式(23)計算速度波動率,若δi=δmax,則跳轉到Step6;否則,采用式(28)對插補參數進行校正,返回Step4;

Step6:判斷是否滿足ui+1≥1,如果滿足,表明此時已經到達NURBS曲線的終點,令ui+1=1,結束插補;否則返回Step2。

將基于四階Runge-Kutta的預估校正得到的插補參數u代入式(7)中即可得到針刺點位置坐標信息。2.2 針刺點姿態信息確定

假設基于四階Runge-Kutta的預估校正法求得的某個針刺點對應的節點參數為u,從而確定節點參數u位于的節點區間,即u∈[u3+i,u4+i]。采用二分搜索法在節點矢量U的子區間{u3,u4,…,u2+n,u3+n}中搜索節點區間的下標i,即可確定針刺點在型值點pi與pi+1之間[7]。

圖7 姿態笛卡爾坐標系建立Fig.7 Establishment of attitude Cartesian coordinate system

則所有針刺點與其相鄰兩個型值點連線組成的邊的集合為

(29)

針刺點相鄰兩個型值點所在三角面的法向量集合可表示為[7]

(30)

(31)

由上述可知,機器人第i個針刺點的姿態信息可以表示為[9]

(32)

由于機器人RPY角姿態描述方式較為常用,將姿態矩陣轉換為RPY方式。機器人姿態的RPY描述方法如圖8所示,首先繞a軸旋轉α角,再繞o軸旋轉β角,最后繞n軸旋轉γ角,α、β、γ分別為其滾動角、俯仰角和偏航角[7]。

圖8 機器人姿態RPY描述方式Fig.8 RPY description mode of robot attitude

機器人RPY姿態轉換矩陣表示為

Rot(a,α)Rot(o,β)Rot(n,γ)=

(33)

式中 cα、cβ、cγ分別為代表cosα、cosβ、cosγ; sα,sβ,sγ分別代表sinα、sinβ、sinγ。

式(33)中包含α、β、γ三個耦合角,為了求解α、β、γ的具體值,需要對矩陣進行解耦。假設機器人期望姿態為

(34)

將式(34)左乘Rot(a,α)-1可得

(35)

將式(35)展開:

(36)

令式(36)中(2,1)元素對應相等:

nycα-nxsα=0

(37)

解得

α=arctan 2(ny,nx)

(38)

令式(36)中(1,1)和(3,1)項分別對應相等:

(39)

解得

β=arctan2[-nz,(nxcα+nysα)]

(40)

令式(36)中(2,2)和(2,3)項分別對應相等得

γ=arctan 2[(-aycα+axsα),(oycα-oxsα)]

(41)

3 仿真驗證

為了驗證上述算法的有效性和穩定性,選取某導彈天線罩預制體三維模型進行仿真實驗,異型預制體三維模型如圖9所示。

圖9 某導彈天線罩預制體三維模型圖Fig.9 3D model drawing of a missile radome preform

選取X軸為切片方向,切片間距設置為100 mm,得到每一個切片層的交點坐標數據并基于交點冗余信息的針刺路徑型值點排序,仿真結果如圖10所示。

圖10 針刺路徑點Fig.10 Needling route point

由于每個切片層切片形狀基本一致,為了簡化仿真流程,選取單個切片層驗證針刺路徑點NURBS曲線插補算法,仿真結果如圖11所示。圖11中,各個坐標原點為針刺點位置,綠色為針刺點法向矢量,紅色為針刺點切向矢量,藍色為針刺點軸向矢量。三個矢量構成了針刺點的姿態信息。

圖11 機器人針刺點位置和姿態生成Fig.11 Generation of robot needling position and attitude

速度波動率描述了插補誤差的大小,進而反映了針刺的均勻性。在此基礎上,運用二階泰勒展開法和基于四階Runge-Kutta的預估校正法對機器人針刺路徑進行插值仿真。仿真插補參數如下:給定進給速度V=100 mm/s,插補周期為T=0.04 s,即插補步長為4 mm,最大速度波動率δmax=10-2%,以速度波動率為判定標準比較二階泰勒展開法以及基于四階Runge-Kutta的預估校正兩種插補參數計算方法。

采用二階泰勒展開法進行仿真,結果如圖12所示。仿真結果表明,速度波動率最大為6.4%,在u=0.2、u=0.6處出現突變,此時針刺點剛好位于針刺軌跡最大曲率處,原因在于插補節點參數u計算時是用給定的弦長進給速度代替理想的切線方向進給速度,這時兩者差距較大,因此出現劇烈變化,說明此時針刺點間距與理想值偏差較大。該算法由于沒有考慮對速度波動率的限制,所以導致速度波動率遠遠超過最大允許速度波動率[17]。

圖12 二階泰勒展開法的速度波動率Fig.12 Velocity fluctuation rate of second order taylor expansion method

如圖13所示,采用基于四階Runge-Kutta的預估校正法使得最大速度波動率限制在給定10-2%內,計算精度明顯提高,原因在于四階Runge-Kutta算法精度比二階泰勒展開法高,同時對計算得到的插補參數進行校正,保證速度波動率限制在預設的范圍內。

圖13 基于四階Runge-Kutta的預估校正法速度波動率Fig.13 Velocity volatility of predictor corrector based on fourth order Runge Kutta

綜合以上分析,基于四階Runge-Kutta的預估校正法在速度波動率方面明顯優于二階泰勒展開法,插補參數的計算更加準確,保證了針刺點的均勻性,且不需要計算高階導數,計算速度也有所提升,充分體現了算法的有效性[17]。

仿真結果表明,該算法靈活性好,穩定性強。可以設置切片間距和切片方向以滿足不同針刺工藝要求的異型預制體針刺點位置和姿態生成,對于實際生產具有一定的意義。

機器人實物驗證驗證在本課題搭建的機器人針刺成形單元中進行,如圖14所示。針刺機器人采用廣州數控(RB50),針刺工具為自主研制的氣動針刺頭,如圖15所示。其余部件還包括空氣壓縮機、電磁閥、工件工裝等。

圖14 機器人針刺單元實驗平臺Fig.14 Robot needling unit experiment platform

圖15 針刺工具三維模型Fig.15 3D model of needling tool

異形預制體針刺加工驗證如圖15所示,針刺點位姿與規劃值基本一致,圖16為針刺效果的局部放大圖,經過測量,針刺間距和針刺深度滿足針刺要求。

圖16 實際加工驗證Fig.16 Actual processing verification

圖17 針刺后的異形預制體Fig.17 Special-shaped preform after needling

4 結論

(1)本文對異型預制體機器人針刺成型路徑規劃進行研究,提出了一種基于切片法的機器人針刺路徑規劃方法。將多個平行切平面與預制體STL模型文件中包含的三角面進行求交運算,從而實現了異形預制體CAD模型到機器人針刺路徑的自動生成。

(2)采用基于四階Runge-Kutta的預估校正法,用于NURBS曲線等步長插補參數的計算,得到針刺點位置;選取針刺點相鄰路徑型值點的法向量和位置坐標,結合機器人RPY姿態描述方式,得到了針刺點的姿態。

(3)仿真結果表明,該路徑規劃算法靈活有效,為異型預制體機器人針刺成型路徑規劃提供了一種簡便有效的方法。

后續,結合此路徑規劃方法進一步提高效率還可采用電動凸輪驅動提高針刺頻率和針刺力,可制作更厚、密度更高的預制體。