王志浩,李洪斌,孫 磊

(1.魯東大學(xué) 交通學(xué)院，山東 煙臺 264025; 2.山東航天電子技術(shù)研究所，山東 煙臺 264670)

0 引言

數(shù)控機床等自動化裝備能夠助力智能化農(nóng)業(yè)機械的制造[1]，提高復(fù)雜零件的加工效率與精準(zhǔn)度，提升農(nóng)用機械的制造水平，實現(xiàn)裝配工藝的優(yōu)化[2]。國家在此方面進(jìn)行了大力扶持，《國家中長期科學(xué)和技術(shù)發(fā)展規(guī)劃綱要(2006—2020年)》將“高檔數(shù)控機床與基礎(chǔ)制造裝備”作為科技重大專項進(jìn)行資助[3]，取得了一定的成果。數(shù)控機床一般由末端執(zhí)行器與調(diào)姿機構(gòu)組成。調(diào)姿機構(gòu)中，基于定位器的并聯(lián)調(diào)姿機構(gòu)具有可重構(gòu)[4-6]、剛度大等優(yōu)點，越來越廣泛地應(yīng)用于智能裝配領(lǐng)域，例如大部件對接中零件的位姿調(diào)整[7]。其中，使用動平臺與多個定位器連接組成的并聯(lián)調(diào)姿機構(gòu)，可用于農(nóng)用飛機裝配時零件的位姿調(diào)整，配合龍門式結(jié)構(gòu)組成的混聯(lián)式智能裝配系統(tǒng)，能夠提供較大的作業(yè)力。

由于農(nóng)用飛機裝配所用的大型圓弧形零件面積較大，導(dǎo)致相應(yīng)的鉆鉚托架體積較大[8]，運動學(xué)標(biāo)定時[9-11]，由重力引起的變形誤差將會降低定位誤差補償精度。針對結(jié)構(gòu)變形問題，很多學(xué)者進(jìn)行了研究，為定位精度補償研究做出了巨大的貢獻(xiàn)。西北工業(yè)大學(xué)張開富等[12]針對傳統(tǒng)方法無法適應(yīng)托架變形量動態(tài)變化的問題，采用迭代方法分析了空間不同位姿下的“串聯(lián)”鉆鉚托架變形量并進(jìn)行了調(diào)平算法設(shè)計，通過仿真分析證明了算法的有效性，但未進(jìn)行試驗驗證;鄧正平等[13]采用有限元仿真的方式分析了大部件裝配過程中不同位置與姿態(tài)下各關(guān)鍵點的變形量，然后進(jìn)行了補償，在一定程度上降低了對接時的定位誤差，但是該方法工作量較大;WANG等[14]考慮重力和關(guān)節(jié)/連桿的柔度，提出一種半解析方法，可用于高速加工和強制裝配的3-SPR并聯(lián)機構(gòu)的柔度分析，但未對結(jié)構(gòu)變形進(jìn)行解析，效率較低。上述研究適用于“串聯(lián)”鉆鉚托架或仿真分析量較小的誤差補償工況，不適用于空間多姿態(tài)下的鉆鉚托架變形解析與誤差補償。

針對這一情況，本文結(jié)合積分法與有限元仿真方法，分析了空間不同位姿下鉆鉚托架的變形量，提出了一種定位點變形量計算的解析解法，建立了變形量模型，在確保變形量計算精度的同時，提高了計算效率；結(jié)合基于目標(biāo)位姿修正的定位誤差補償方法，提高了定位誤差補償精度。

1 結(jié)構(gòu)分析

本文以五自由度混聯(lián)式智能裝配系統(tǒng)中的鉆鉚托架[8]為研究對象，驗證了所提出算法的有效性。智能裝配系統(tǒng)由并聯(lián)調(diào)姿機構(gòu)、上末端執(zhí)行器、下末端執(zhí)行器組成，適用于農(nóng)用飛機壁板類零件加工。并聯(lián)調(diào)姿機構(gòu)用于調(diào)整飛機壁板類零件的位置與姿態(tài)，由4個定位器和1個鉆鉚托架組成；定位器1具有z向移動自由度；定位器2具有x向、z向移動自由度，兩移動副互相垂直；定位器3、4由3個互相垂直的移動副組成，即具有x、y、z向移動自由度；定位器與鉆鉚托架通過球鉸副連接，如圖1所示。全局坐標(biāo)系Ob-xbybzb原點Ob位于定位器1的伸縮桿軸線與地面交點，全局坐標(biāo)系各坐標(biāo)軸分別平行于定位器各方向上的移動副；鉆鉚托架坐標(biāo)系Ot-xtytzt原點位于鉆鉚托架上平面中心點，當(dāng)鉆鉚托架繞坐標(biāo)系Ob-xbybzb各軸轉(zhuǎn)角為零時，坐標(biāo)系Ot-xtytzt各坐標(biāo)軸平行于坐標(biāo)系Ob-xbybzb各坐標(biāo)軸；各坐標(biāo)系的坐標(biāo)軸方向如圖1所示。圖1b中，S1,…,S4表示球鉸副，Zp1、Xp2,…,Zp4表示移動副。

根據(jù)Kutzbach-Grübler式[15]可計算鉆鉚托架的自由度：

6×(11-13-1)+21=3。

(1)

式中：M表示機構(gòu)的自由度，nm表示構(gòu)件數(shù)目，gm表示運動副數(shù)目，fh表示第h個運動副的自由度。由式(1)可知，鉆鉚托架具有3個自由度，即繞全局坐標(biāo)系x軸、y軸的轉(zhuǎn)動自由度，沿全局坐標(biāo)系z軸的移動自由度。

鉆鉚托架長、寬、高分別為4 520 mm、2 100 mm、150 mm，跨度、自重較大，在重力的作用下將會產(chǎn)生彈性變形，最大變形量將達(dá)到0.25 mm，對于0.5 mm(0.5°)位置與姿態(tài)精度要求而言，影響較大，并且空間不同姿態(tài)下的變形量不同，因此需要進(jìn)行空間多姿態(tài)下的鉆鉚托架變形解析與誤差補償研究。

2 空間多姿態(tài)下的鉆鉚托架變形解析建模

不同位姿下，鉆鉚托架同一點處的變形量不同，因此需要進(jìn)行空間多姿態(tài)下的變形量分析。鉆鉚托架為中心對稱圖形，變形解析時，可將鉆鉚托架進(jìn)行簡化。

將鉆鉚托架簡化為簡支梁，采用積分法解析簡支梁在重力作用下各點撓度，然后與ABAQUS分析結(jié)果對比，分析二者的計算差值；采用加權(quán)整體最小二乘法對差值進(jìn)行擬合并補償，降低噪聲對變形計算誤差的影響。

鉆鉚托架長邊的長度為短邊長度的2.15倍，空間不同位姿下，長邊的撓度最大值為短邊撓度最大值的2.154倍左右。因此，在進(jìn)行鉆鉚托架簡化時，忽略短邊變形對鉆鉚托架變形量分析結(jié)果的影響；在變形解析修正時，補償因簡化短邊所造成的變形計算差值。同時，忽略鉆鉚托架制造誤差、定位器初始位置誤差對變形解析結(jié)果的影響。

2.1 基于積分法的鉆鉚托架變形解析

考慮鉆鉚托架重力對其自身變形的影響，等價于在梁上施加均布力，如圖2所示。令a1a2段、a5a6段的均布力為q=[0,0,-q]T，根據(jù)a3a4段截面的面積可計算該段均布力為1.35q?？墒褂梅e分法分析簡化后的鉆鉚托架各點的撓度。

以均布力q為例，當(dāng)鉆鉚托架轉(zhuǎn)動一定角度時，均布力q在鉆鉚托架坐標(biāo)系Ot-xtytzt各方向上的分量可以表示為式(2)：

(2)

q=mbg/l=ρVg/l=ρSg。

(3)

式中:q表示長度為l、質(zhì)量為mb、密度為ρ的梁所受均布力；S表示該段梁截面積；V表示該段梁體積。

(1)a1a2段

(4)

(2)a2a3段

圖3中，虛線表示簡化后的a2a3段梁。a2a3段截面高度可以使用式(5)計算：

(5)

hi2計算方式類似。慣性矩函數(shù)I1、I3、I5的一階導(dǎo)數(shù)與二階導(dǎo)數(shù)為零，慣性矩函數(shù)I2、I4的一階導(dǎo)數(shù)與二階導(dǎo)數(shù)約等于零。首先以均布力分量qzt為例，分析鉆鉚托架在該分量作用下的變形量vzt。由文獻(xiàn)[20]可知，式(6)成立：

d2Mzt(x)/dx2=qzt(x),

d2vzt(x)/dx2=Mzt(x)/EI(x),

(6)

(7)

AX=B。

(8)

式中，20階矩陣A為有關(guān)E、lm1、lm2、I4、I5的矩陣；20×1矩陣B為有關(guān)qzt、lm1、lm2、lm5、I1、I2、I5的矩陣，A與B的表達(dá)式已知。矩陣X為有關(guān)vzt的系數(shù)矩陣，可以使用式(8)求解，進(jìn)一步得出vzt表達(dá)式。表達(dá)式vyt可以使用類似的方法得到。qzt引起的鉆鉚托架變形量沿全局坐標(biāo)系Ob-xbybzb的xb向分量可以表示為式(9)：

(9)

需要說明的是，不同轉(zhuǎn)角條件下的撓曲線方程不同，需要分別計算撓曲線方程中的未知系數(shù)。同時，忽略定位點加工時的位置誤差引起的定位點空間位置變化。

2.2 解析精度分析

進(jìn)一步，使用ABAQUS軟件分析不同位姿下的鉆鉚托架受自重影響時定位點的坐標(biāo)變化量并與解析結(jié)果進(jìn)行對比。鉆鉚托架材料為45號鋼，彈性模量210 GPa，泊松比0.3，密度7 890 kg/m3，屈服強度355 MPa；仿真分析時，重力加速度設(shè)置為-9.8 m/s2，鉆鉚托架與定位器連接處施加位移約束，即位移量為0 mm。進(jìn)行鉆鉚托架變形量分析時，不考慮伸縮桿變形的影響。

由于鉆鉚托架上存在小尺寸凸臺、螺紋孔等復(fù)雜幾何形狀，網(wǎng)格劃分時，采用適應(yīng)性較好的四面體網(wǎng)格[21]對該類型區(qū)域進(jìn)行處理；對于其他區(qū)域，采用分析精度較高的六面體網(wǎng)格[22]進(jìn)行處理。鉆鉚托架定位點如圖4所示，圖中m表示定位點編號。

定位點m在全局坐標(biāo)系下的空間位置變化可用式(10)計算：

(10)

使用ABAQUS軟件對比分析不同位姿下鉆鉚托架的變形時，選用的位姿組數(shù)越多，定位點變形誤差的計算精度越高。本文考慮有限元分析的工作量與算法效率，選用10組位姿。位姿選用時，需要確保選取的位姿處于并聯(lián)調(diào)姿機構(gòu)的工作空間內(nèi)，并且應(yīng)確保鉆鉚托架轉(zhuǎn)角范圍由零逐漸增大至較大值，以保證定位點變形誤差的計算精度。

2.3 基于加權(quán)整體最小二乘法的變形解析修正

變形解析修正是指用解析法計算修正曲線，補償各定位點的變形計算誤差。本文擬通過五次多項式擬合變形量修正曲線。假定16×6矩陣A為系數(shù)矩陣；6×1矩陣xk為待求解的未知量矩陣；16×1矩陣b為各曲線指定點處的修正量矩陣，由各定位點處變形量差值取均值后得出。

如圖4所示，由于定位孔的影響，矩陣A中元素存在一定的誤差；由于均值計算、截面近似處理、簡化梁2(或4)的影響，矩陣b也存在一定的誤差。當(dāng)A與b皆存在噪聲(誤差)時，基于殘差平方和最小的普通最小二乘法求解出的曲線無法高精度擬合各點的變形修正量，可以使用加權(quán)整體最小二乘法[23]對曲線系數(shù)進(jìn)行求解，使修正補償結(jié)果更加貼近真實值。

假定16×6矩陣EA為矩陣A的噪聲矩陣，16×1矩陣Eb為矩陣b的噪聲矩陣，因此有：

(A+EA)xk=b+Eb。

(11)

定義：

eA=vec(EA)。

(12)

其中:vec函數(shù)表示矩陣按列拉直運算[24]；eb的定義與eA類似。將eb作為關(guān)于xk和eA的非線性函數(shù)，將式(11)展開并忽略二階及以上極小量可得：

(13)

(14)

式中:Pb為矩陣b的權(quán)矩陣[25]，為16階矩陣；PA為矩陣A的權(quán)矩陣，為96階矩陣。相應(yīng)的迭代求解算法設(shè)計為：

(2)分別定義矩陣112×102矩陣M、112×1矩陣N：

(15)

(16)

式中:?表示Kronecker積；eA為96×1矩陣；Iq表示96階單位陣；Ze表示矩陣中元素為0的矩陣。進(jìn)一步，根據(jù)拉格朗日極值條件可得：

(17)

式(17)中函數(shù)連續(xù)，λ、μ為拉格朗日乘子。對(17)中分量求偏導(dǎo)并令其值為零，經(jīng)化簡可得：

X=(MTPM)-1MTPN。

(18)

式中:

(19)

(20)

其中P為權(quán)矩陣，為112階矩陣。

使用上述方法可得出代表x向和z向的變形修正方程Cx(x)和Cz(x)。因此，修正后的鉆鉚托架變形量模型可以表示為式(21):

(21)

3 變形解析與誤差補償仿真分析

為驗證所提算法的有效性，使用MATLAB軟件與ABAQUS軟件進(jìn)行了變形解析與定位誤差補償仿真分析。使用ABAQUS軟件分析變形量時的網(wǎng)格劃分結(jié)果如圖5所示，總單元數(shù)目為525 896個；其中線性六面體(linear hexahedral)單元數(shù)目24 370個，類型為C3D8R，該類型網(wǎng)格適用于靜力與動力分析，在彎曲荷載下不易發(fā)生剪切自鎖現(xiàn)象，位移求解結(jié)果比較精確，分析的精度受網(wǎng)格扭曲影響較??；二次四面體單元(quadratic tetrahedral)數(shù)目為501 526個，類型為C3D10。

表1 變形解析所需的位置與姿態(tài)

IDxbt/mmybt/mmzbt/mmαbt/(°)βbt/(°)1-2 2051 0251 140002-2 195.86999.691 303.785.601.433-2 189.66980.701 372.438.142.084-2 208.541 016.211 214.841.96-0.785-2 196.831 039.011 304.26-3.922.216-2 200.25998.181 274.435.350.787-2 192.51989.831 338.236.731.828-2 210.551 010.641 249.603.08-1.179-2 209.621 033.591 249.73-2.24-1.5610-2 199.201 032.601 239.79-1.961.43

通過重力作用下定位點空間位置變化，對比本文所提出的鉆鉚托架變形解析方法與ABAQUS軟件分析結(jié)果的差值。鉆鉚托架被調(diào)整至表1中第1～5位姿，對比修正前本文所提出的變形量解析結(jié)果與ABAQUS分析結(jié)果。以位姿2為例，相應(yīng)的對比結(jié)果如圖7所示。

進(jìn)一步，分析了本文所提出的變形量解析方法與ABAQUS分析結(jié)果的差值，如圖8所示。

由圖8可以看出，小轉(zhuǎn)角范圍內(nèi)：

(1)兩種方法在鉆鉚托架變形量的分析上存在一定的差值，以ABAQUS軟件分析的變形量為基準(zhǔn)，本文提出的變形量解析方法需要進(jìn)行修正；

(2)本文提出的鉆鉚托架變形解析方法在全局坐標(biāo)系y向的變形差值(數(shù)量級10-3mm)相對x向和z向差值(數(shù)量級10-2mm)較小，因此僅對x向和z向差值進(jìn)行修正；

(3)在不同的空間位置與姿態(tài)下，x向差值與z向差值在一定范圍內(nèi)變化，該變化特征符合一定的數(shù)學(xué)規(guī)律，可以根據(jù)該數(shù)學(xué)規(guī)律對本文提出的鉆鉚托架變形解析方法進(jìn)行修正，建立鉆鉚托架的變形量模型，即式(21)。

為了驗證修正后的鉆鉚托架變形解析方法的精度，鉆鉚托架被調(diào)整至表1中第6～10位姿，對比變形量模型的解析結(jié)果與ABAQUS分析結(jié)果，以位姿10為例，相應(yīng)的對比結(jié)果如圖9所示。

進(jìn)一步對比修正后的解析結(jié)果與ABAQUS軟件計算結(jié)果的差值，以差值最大的位姿10為例，繪制變形差值曲線，如圖10所示；其他位姿對應(yīng)的變形差值皆小于位姿10的變形差值。

由圖8、圖10可以看出：

(1)修正后，x、z向變形量差值明顯降低；

(2)本文提出的鉆鉚托架變形量解析方法與ABAQUS分析結(jié)果間的最大差值為8.4×10-3mm，遠(yuǎn)小于各方向的變形量，滿足精度要求。

結(jié)合文獻(xiàn)[26]方法，進(jìn)行考慮結(jié)構(gòu)變形的運動學(xué)標(biāo)定與定位誤差補償，并與不考慮變形誤差方法對應(yīng)的標(biāo)定與補償結(jié)果進(jìn)行對比，分析二者的精度，如圖11所示。圖11中，針對是否考慮結(jié)構(gòu)變形與測量不確定度，定義了4種工況；橫坐標(biāo)“位姿編號”中的位姿表示鉆鉚托架處于空間中某一位姿；鉆鉚托架處于空間中某一位姿時，會存在位置誤差與姿態(tài)誤差，因此縱坐標(biāo)表示為位置誤差或姿態(tài)誤差?？紤]測量不確定度的方式為：計算仿真時的激光跟蹤儀測量的定位點坐標(biāo)時，將測量不確定度(測量誤差)施加到所計算的點坐標(biāo)中。

相對于補償前的位置與姿態(tài)誤差，誤差補償后，不同工況對應(yīng)的鉆鉚托架的定位誤差皆有所下降，但在補償精度方面存在差異，如圖11所示，在給定幾何誤差參數(shù)、測量不確定度(0.1 mm)的條件下：

(1)考慮測量不確定度的影響，誤差補償精度將會有所下降；

(2)考慮結(jié)構(gòu)變形對應(yīng)的誤差補償精度優(yōu)于忽略變形誤差的誤差補償精度。原因在于，輸入誤差模型的位姿誤差計算精度升高，確保了誤差辨識精度，在誤差補償時，各移動副的運動量計算精度也隨之提升，因此誤差補償精度較優(yōu)；

(3)本文所提出的變形解析方法對應(yīng)的變形量計算精度能夠滿足誤差補償所需的精度要求。

文獻(xiàn)[26]使用ABAQUS軟件分析了鉆鉚托架在空間不同位姿下的變形量，單次仿真運行時間約65 s。本文提出了變形量解析計算方法，并在有限元仿真的基礎(chǔ)上進(jìn)行了變形解析修正，所需的有限元仿真次數(shù)明顯減少，可有效降低空間多姿態(tài)下的鉆鉚托架變形仿真分析工作量；單次MATLAB程序運行所需時間約0.35 s。

4 誤差補償試驗驗證

為驗證所提出的變形解析算法的有效性，結(jié)合文獻(xiàn)[26]所提出的誤差補償方法進(jìn)行鉆鉚托架的定位誤差補償試驗，如圖12所示。試驗所用的激光跟蹤儀為法如X01型，上位機型號為惠普4411S，采用自行搭建的控制軟件調(diào)整鉆鉚托架位姿。試驗時的測量條件如下：工作環(huán)境10.5 ℃，空氣濕度32%，海拔高度22 m；激光跟蹤儀測量誤差(測量不確定度)為0.1 mm。

補償后的鉆鉚托架位置誤差與姿態(tài)誤差如圖13所示。圖13b中6、7、8號位姿對應(yīng)的“不考慮結(jié)構(gòu)變形”和“考慮結(jié)構(gòu)變形”的姿態(tài)誤差差別很小，可能原因是：鉆鉚托架姿態(tài)變化幅度較小，導(dǎo)致兩方法間的補償精度差別不明顯。

分析圖13中數(shù)據(jù)，可以得出誤差補償后鉆鉚托架定位誤差的最小值、最大值、平均值與均方差，如表2與表3所示。

表2 位置誤差補償精度分析

表3 姿態(tài)誤差補償精度分析

由圖13、表2與表3可以看出：

(1)誤差補償前，鉆鉚托架的位置誤差最大值為2.91 mm，姿態(tài)誤差最大值為0.57°，大于0.5 mm(0.5°)的定位精度要求，因此需要進(jìn)行定位誤差補償；

(2)相對于未考慮鉆鉚托架變形的定位誤差補償方法，考慮鉆鉚托架變形時，補償后的鉆鉚托架位置誤差最大值降低67.52%，姿態(tài)誤差最大值降低17.78%，同時滿足農(nóng)業(yè)航空裝配領(lǐng)域的定位精度要求；

(3)本文所提出的定位誤差補償方法可靠性較高、實用型較強，能夠用于農(nóng)業(yè)航空裝配領(lǐng)域大型鉆鉚托架的定位誤差補償。

5 結(jié)束語

本文提出一種空間多姿態(tài)下的鉆鉚托架變形解析方法，并進(jìn)行了誤差補償實驗：

(1)采用積分法解析了空間不同位姿下的鉆鉚托架變形量，建立了相應(yīng)的變形量模型，確保定位點變形量計算誤差小于10-2mm。

(2)相對于已有研究，本文側(cè)重于鉆鉚托架在空間不同位姿下的變形解析，推導(dǎo)了空間不同位姿下的變形量表達(dá)式，建立了變形量模型；相對于使用ABAQUS等軟件的仿真分析方法，本方法在保證計算精度的同時，計算效率得到了大幅度提高。

(3)本文使用加權(quán)整體最小二乘法，確保了鉆鉚托架的變形量擬合精度；補償后的鉆鉚托架位置誤差最大值降低67.52%，姿態(tài)誤差最大值降低17.78%，滿足了農(nóng)業(yè)航空裝配領(lǐng)域的定位精度要求。

下一步，將重點研究鉆鉚托架其余部分變形對定位精度的綜合影響。

致謝

本文在南京航空航天大學(xué)機電學(xué)院與無人機研究院進(jìn)行了試驗分析，對此表示感謝。