水下機(jī)器人試航速度的類物理數(shù)值方法預(yù)報(bào)

吳利紅， 張愛鋒， 李一平， 封錫盛， 王詩文

(1. 大連海事大學(xué) 船舶與海洋工程學(xué)院，遼寧 大連 116026；2. 中國科學(xué)院沈陽自動(dòng)化研究所 機(jī)器人學(xué)國家重點(diǎn)實(shí)驗(yàn)室，遼寧 沈陽 110016)

預(yù)報(bào)水下機(jī)器人(autonomous underwater vehicle，AUV)定常航行的巡航速度和非定常運(yùn)動(dòng)的操縱性能不僅是評估其推進(jìn)系統(tǒng)性能的重要指標(biāo)，還有助于提高AUV航程估計(jì)，能量估算和減少定位誤差和AUV丟失的風(fēng)險(xiǎn)。采用類物理的數(shù)值方法進(jìn)行海洋載體操縱預(yù)報(bào)是數(shù)值模擬的目標(biāo)。此方法應(yīng)用真實(shí)模型模擬螺旋槳高速旋轉(zhuǎn)，能反映船槳舵內(nèi)在作用機(jī)理，再現(xiàn)流場特征。但它具有模型復(fù)雜，網(wǎng)格尺度相差大，載體各部件相對運(yùn)動(dòng)復(fù)雜，高頻與低頻運(yùn)動(dòng)混合的模擬難點(diǎn)[1]。自20世紀(jì)90年代美國海軍提出這個(gè)概念以來，直至2010年左右陸續(xù)有研究成果。如Chase等[1]、Pankajakshan等[2]、Poremba等[3]對潛艇自航運(yùn)動(dòng)進(jìn)行了數(shù)值模擬；Mofidi等[4]對集裝箱船KCS帶全附體槳舵的ZigZag操縱運(yùn)動(dòng)的模擬；Carrica等[5]對水面戰(zhàn)艦的隨浪橫甩運(yùn)動(dòng)進(jìn)行自航模擬；沈志榮[6]對水面船舶操縱運(yùn)動(dòng)進(jìn)行數(shù)值模擬。以上研究均采用動(dòng)態(tài)重疊網(wǎng)格法，在邊界移動(dòng)時(shí)無需網(wǎng)格更新，但須實(shí)時(shí)確定移動(dòng)邊界與背景網(wǎng)格的重疊區(qū)域，這是數(shù)值誤差的源頭[7]。而移動(dòng)網(wǎng)格法在邊界移動(dòng)時(shí)網(wǎng)格也隨著運(yùn)動(dòng)，運(yùn)動(dòng)更符合拉格朗日描述法，但邊界移動(dòng)常導(dǎo)致復(fù)雜的網(wǎng)格運(yùn)動(dòng)，從而引起網(wǎng)格畸變和網(wǎng)格質(zhì)量降低，影響數(shù)值計(jì)算精度。受此限制，移動(dòng)網(wǎng)格法主要應(yīng)用在二維物體的運(yùn)動(dòng)：圓球的振動(dòng)，魚類的擺鰭運(yùn)動(dòng)[8]和三維物體的有限區(qū)間或簡單的運(yùn)動(dòng)，如翼的偏轉(zhuǎn)[9]，預(yù)定義AUV直航對接運(yùn)動(dòng)[10]。而本文則采用多塊動(dòng)態(tài)混合網(wǎng)格方法[11]解決移動(dòng)網(wǎng)格的網(wǎng)格畸變和數(shù)值精度差的問題[12]，實(shí)現(xiàn)AUV自航試驗(yàn)的數(shù)值模擬，獲得了AUV巡航速度預(yù)報(bào)，再現(xiàn)了AUV自航過程非定常運(yùn)動(dòng)的航行特性。

1 水下機(jī)器人自航模型建立

AUV為改進(jìn)型REMUS AUV，L=1.7 m，D=0.191m[13]。具有橢球型艏、圓柱形平行中體和圓錐型艉。艉部有“十字”型舵翼以及一個(gè)MAU4-40，直徑d=0.156 m的螺旋槳。AUV艏、舯和艉3段以及舵翼采用結(jié)構(gòu)網(wǎng)格繪制，螺旋槳采用非結(jié)構(gòu)網(wǎng)格離散。圖1為AUV帶舵翼和螺旋槳的網(wǎng)格圖。大地坐標(biāo)系為oxyz,ox為AUV對稱軸，向艏為正(向前)；oy垂直向上為正；oz指向右舷。

圖1 帶槳和舵的AUV網(wǎng)格模型Fig.1 Grid model of AUV appended propeller and rudders

圖2為AUV帶槳舵流域在對稱面上網(wǎng)格分區(qū)，分區(qū)劃分是根據(jù)AUV的運(yùn)動(dòng)趨勢確定的。流域包含 6個(gè)分區(qū)。分區(qū)Ⅰ包含AUV和舵翼表面的四邊形網(wǎng)格以及近場邊界層區(qū)域的六面體結(jié)構(gòu)網(wǎng)格；分區(qū)Ⅱ包括螺旋槳表面的三角形網(wǎng)格以及螺旋槳擾動(dòng)子區(qū)域的四面體非結(jié)構(gòu)網(wǎng)格；分區(qū)Ⅲ包含AUV、螺旋槳和舵翼擾動(dòng)區(qū)域的六面體結(jié)構(gòu)網(wǎng)格；分區(qū)Ⅳ為AUV從艏部延伸到遠(yuǎn)場的航行區(qū)域網(wǎng)格，靠近艏部加密，遠(yuǎn)場較為稀疏，由六面體結(jié)構(gòu)網(wǎng)格構(gòu)成；分區(qū)Ⅴ為從AUV艉部延伸到遠(yuǎn)場的結(jié)構(gòu)網(wǎng)格區(qū)域，靠近AUV艉部為加密區(qū)，遠(yuǎn)場為稀疏區(qū)，也由六面體結(jié)構(gòu)網(wǎng)格構(gòu)成；外圍非結(jié)構(gòu)網(wǎng)格區(qū)域?yàn)榉謪^(qū)Ⅵ。

圖2 動(dòng)網(wǎng)格分區(qū)Fig.2 Dynamic meshing zones

在AUV自航推進(jìn)數(shù)值模擬過程中，網(wǎng)格運(yùn)動(dòng)區(qū)域限制在圖2的AUV直航窄條帶區(qū)域，這有利于局部擾動(dòng)區(qū)域梯度的捕捉。其中分區(qū)Ⅱ既有隨螺旋槳的旋轉(zhuǎn)運(yùn)動(dòng)，又有隨AUV的直航運(yùn)動(dòng)；分區(qū)Ⅰ和Ⅲ僅隨AUV作直航運(yùn)動(dòng)；分區(qū)Ⅳ和Ⅴ分別為網(wǎng)格壓縮區(qū)域和拉伸區(qū)域，在分區(qū)Ⅳ的最前端為靜止壁面，采用動(dòng)態(tài)層方法對區(qū)域Ⅳ的最前端網(wǎng)格進(jìn)行壓縮合并，因此最前端為網(wǎng)格消亡區(qū)域；分區(qū)Ⅴ的最后端為動(dòng)態(tài)層網(wǎng)格的拉伸區(qū)域，為網(wǎng)格增長區(qū)域。外圍區(qū)域Ⅵ為網(wǎng)格變形區(qū)。各個(gè)分區(qū)之間的連接采用非一致連接。這種多塊結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格混合的網(wǎng)格構(gòu)成方法，以及采用動(dòng)態(tài)層網(wǎng)格更新方法，是借鑒了動(dòng)態(tài)重疊網(wǎng)格的子區(qū)域網(wǎng)格運(yùn)動(dòng)的思想，從網(wǎng)格構(gòu)建質(zhì)量、網(wǎng)格更新質(zhì)量和網(wǎng)格更新速度等方面著手，提高數(shù)值求解的速度和計(jì)算精度。

AUV自航推進(jìn)過程數(shù)值模擬中，AUV和螺旋槳之間力和速度的傳遞可以采用用戶自定義函數(shù)UDF(user defined function)進(jìn)行編寫，嵌入到程序中[14]：首先，AUV靜止，螺旋槳以恒定轉(zhuǎn)速旋轉(zhuǎn)，產(chǎn)生推力，將此推力保存同時(shí)將推力傳遞給AUV；AUV讀取螺旋槳推力，同時(shí)實(shí)時(shí)求解URANS(unsteady reynolds averaged navier-stokes)方程，獲得AUV的阻力。AUV在推力和阻力的聯(lián)合作用下，通過在6-DOF方程加載推力和實(shí)時(shí)阻力，求解載體的加速度，并采用龍格庫塔積分獲得AUV實(shí)時(shí)速度V；此速度傳遞給螺旋槳，使AUV拖曳螺旋槳以V航速運(yùn)動(dòng)；隨著AUV航行，AUV的尾流場更新，影響螺旋槳的進(jìn)流場，螺旋槳在新的進(jìn)速和轉(zhuǎn)速下產(chǎn)生新的推力，重新記錄此推力并傳遞給AUV；AUV在阻力作用下，獲得新的航速，如此重復(fù)，直到AUV的推力T與其航行阻力R接近平衡，速度曲線趨于穩(wěn)定，加速度接近為零?；?qū)Ρ群剿倥c多參考系坐標(biāo)[15](multiple frames of reference, MFR)的定常自航點(diǎn)航速，兩者誤差在5%內(nèi)，也認(rèn)為達(dá)到自航點(diǎn)。記錄此時(shí)的推力、阻力、速度歷史數(shù)據(jù)，保存和分析。

2 水下機(jī)器人自航數(shù)值結(jié)果與分析

2.1 受力和速度

AUV 以不同轉(zhuǎn)速從靜止自航推進(jìn)，最終會(huì)趨于勻速狀態(tài)。分別假定AUV以300、450、600 r/min 3種轉(zhuǎn)速航行，每種轉(zhuǎn)速自航一次，迭代步長為螺旋槳旋轉(zhuǎn)1°對應(yīng)的時(shí)間，每個(gè)外循環(huán)包含內(nèi)循環(huán)20次，收斂精度為10-4。如圖3分別為轉(zhuǎn)速為300、450和600 r/min的自航螺旋槳推力和AUV阻力的歷史變化曲線。推力和阻力曲線在數(shù)值求解中有點(diǎn)振蕩，這是動(dòng)網(wǎng)格數(shù)值求解引起的不穩(wěn)定現(xiàn)象。阻力部分振蕩主要是由于載體低速巡航采用了和螺旋槳高速旋轉(zhuǎn)的小時(shí)間步引起的。隨著速度增加，推力減小，阻力增大，最終推力和阻力平衡，速度會(huì)趨于一個(gè)恒定的速度，這個(gè)速度即為該轉(zhuǎn)速下的自航點(diǎn)。3種轉(zhuǎn)速可得3次航速歷史變化，如圖4所示，可得3個(gè)不同的自航點(diǎn)。轉(zhuǎn)速越大，達(dá)到勻速的時(shí)間越短。600 r/min，需要約9 s, 可達(dá)到勻速1.55 m/s；450 r/min，需要約10.53 s, 可達(dá)到勻速1.15 m/s； 300 r/min，需要約14.39 s, 可達(dá)到勻速0.74 m/s。同時(shí)，轉(zhuǎn)速越大，起始加速度越大。

圖3 不同轉(zhuǎn)速下AUV自航阻力和推力變化Fig.3 The resistance (R) and Thrust (T) of AUV self-propulsion at different rotation speeds

圖4 不同轉(zhuǎn)速下AUV自航過程的航速變化Fig.4 Velocity history (V) in AUV self-propulsion at different rotation speeds (n)

2.2 數(shù)值驗(yàn)證

為了驗(yàn)證數(shù)值模擬結(jié)果的準(zhǔn)確性，將AUV自航的數(shù)值結(jié)果與水池試驗(yàn)(試驗(yàn)在中國科學(xué)院沈陽自動(dòng)化所的深水水池：L×B×D=25 m×15 m×9 m)結(jié)果進(jìn)行對比。圖5記錄了螺旋槳轉(zhuǎn)速為300 r/min的AUV速度歷史數(shù)據(jù)。包括5次試驗(yàn)結(jié)果，以離散數(shù)據(jù)點(diǎn)表示：Exp1、Exp2、Exp3、Exp4 和Exp5；數(shù)值結(jié)果以實(shí)線表示。數(shù)值模擬不受試驗(yàn)水池邊界的限制，因此數(shù)值模擬的時(shí)間較長。數(shù)值結(jié)果和試驗(yàn)結(jié)果吻合良好。此外，為了驗(yàn)證不同轉(zhuǎn)速下AUV自航點(diǎn)的可靠性，將本文的結(jié)果與多參考坐標(biāo)系MFR(multi-frame of reference)的結(jié)果進(jìn)行對比驗(yàn)證，如圖6所示，獲得的自航點(diǎn)速度略低于MFR，誤差在5%以內(nèi)。在一段轉(zhuǎn)速范圍內(nèi)，航速和轉(zhuǎn)速成正比。這是因?yàn)橄嘟囊欢无D(zhuǎn)速范圍內(nèi)，推力系數(shù)接近常數(shù)，阻力系數(shù)接近常數(shù)。螺旋槳推力與轉(zhuǎn)速成平方，勻速后推力與阻力相等，阻力與速度成平方，則自航點(diǎn)速度和轉(zhuǎn)速成正比。這個(gè)結(jié)論有助于更快找到不同轉(zhuǎn)速下的自航點(diǎn)。由轉(zhuǎn)速和航速的關(guān)系，根據(jù)設(shè)計(jì)航速1.5 m/s的要求，對圖6進(jìn)行插值，可得對應(yīng)自航點(diǎn)的轉(zhuǎn)速為570 r/min。

圖5 AUV自航速度變化的試驗(yàn)和數(shù)值結(jié)果 (n=300 r/min)Fig.5 The velocity history from experimental and numerical results in AUV self-propulsion at n=300 r/min

圖6 自航點(diǎn)對應(yīng)的轉(zhuǎn)速Fig.6 Approaching velocities vs. rotating speeds

2.3 流場分析

圖7為轉(zhuǎn)速n為300、 450和600 r/min時(shí)達(dá)到勻速直航的AUV對稱面上的速度云圖(注：遠(yuǎn)場速度為零)?？梢?，轉(zhuǎn)速越大，螺旋槳的艉跡越明顯，梢渦強(qiáng)度越大，向后推水的水速度越大，反作用力越大，則推力越大。且轉(zhuǎn)速越大，AUV艏部的速度擾動(dòng)流場范圍越大，表明AUV達(dá)到自航點(diǎn)對應(yīng)的速度越大。3種勻速狀態(tài)時(shí)，AUV艏部最大截面處都產(chǎn)生了艏渦，且隨著速度增加，艏渦范圍加大，表明AUV艏部粘壓阻力增加。隨著AUV速度增加，AUV艉部截面變化區(qū)域開始產(chǎn)生艉渦(n=600 r/min)，表明艉部速降增加，艉部粘壓阻力增加。同時(shí)，3種航速時(shí)，AUV表面都有較明顯的邊界層，從艏至艉逐步增厚。

圖7 AUV不同轉(zhuǎn)速對應(yīng)的自航點(diǎn)速度云圖Fig.7 Velocity contours of AUV at different self-propulsion points with varying rotation speeds

圖8為設(shè)計(jì)航速1.5 m/s對應(yīng)轉(zhuǎn)速570 r/min 的AUV自航過程4個(gè)典型時(shí)刻(t為0.1，3,4.9和6.9 s)對稱面上AUV的速度分布(右圖)和螺旋槳艉跡局部放大圖(左圖)。相應(yīng)的壓力云圖如圖9所示?？梢?，當(dāng)轉(zhuǎn)速一定時(shí)，AUV速度較低時(shí)，螺旋槳艉跡強(qiáng)，表明螺旋槳此時(shí)的推力較大(這與螺旋槳的敞水性能曲線一致，和圖3的推力曲線一致)。螺旋槳葉梢曳出的梢渦出現(xiàn)重疊；隨著AUV航速增加，梢渦強(qiáng)度降低，表明螺旋槳推力減少。但是梢渦分離清晰可見，這是由AUV航行時(shí)的艉跡反向作用引起的。

圖8 AUV自航對稱面上不同時(shí)刻的速度云圖(n=570 r/min)Fig.8 Velocity contours of AUV self-propulsion on the symmetry plane at different times with n=570 r/min

圖9 AUV自航對稱面上不同時(shí)刻的壓力云圖(n=570 r/min)Fig.9 Pressure contours of AUV self-propulsion at the symmetry plane at different times with n=570 r/min

隨著AUV航速的增加， 螺旋槳產(chǎn)生的推力在變小，這點(diǎn)在圖9的左側(cè)壓力云圖中清晰可見。低速時(shí)，螺旋槳葉面的高壓值很大，葉背的低壓值也較大，導(dǎo)致兩者的壓差較大，因此推力較大；隨著AUV速度增加，葉面高壓值減小，葉背低壓值也減小，兩者的壓差變小，則推力也在減少。但是，隨著螺旋槳持續(xù)旋轉(zhuǎn)，只要AUV推力大于阻力，AUV就獲得加速，隨著AUV速度增加，AUV艏部速度繞流場范圍加大，艏部駐點(diǎn)的壓強(qiáng)增加，導(dǎo)致AUV艏艉壓差阻力增加；同時(shí)艏肩部速度增加，導(dǎo)致摩擦阻力增加，AUV總阻力也在增加，直至與逐漸減小的推力平衡，達(dá)到勻速直航，即達(dá)到自航點(diǎn)。隨著螺旋槳推進(jìn)AUV自航，AUV艉跡延長。螺旋槳艉跡梢渦強(qiáng)度隨AUV速度增加弱化。同時(shí)，在螺旋槳轂帽附近有一股根渦，和AUV航行方向一致，流速大小隨著AUV速度增加而增加。

3 結(jié)論

1) 水下機(jī)器人不同轉(zhuǎn)速對應(yīng)的最終巡航速度不同。300、450和600 r/min分別對應(yīng)的巡航速度為0.74、1.15和1.55 m/s。AUV的巡航速度1.5 m/s對應(yīng)的轉(zhuǎn)速是570 r/min。

2) AUV自航穩(wěn)定所需時(shí)間與AUV轉(zhuǎn)速大小有關(guān)。轉(zhuǎn)速越大，加速度越大，穩(wěn)定需要時(shí)間越短。300、450和600 r/min穩(wěn)定分別需要14.39、10.53和 9 s。

3) 在螺旋槳轉(zhuǎn)速變化較小的范圍內(nèi)，AUV的自航點(diǎn)速度和轉(zhuǎn)速成正比關(guān)系。

4) AUV自航時(shí)，螺旋槳旋轉(zhuǎn)運(yùn)動(dòng)曳出梢渦和根渦。梢渦強(qiáng)度隨著AUV航速增加而降低，運(yùn)動(dòng)方向與AUV航行方向相反。根渦大小隨AUV航速增加而增加，方向與AUV航行方向一致。

5) 螺旋槳推力的變化起源于其葉面和葉背壓差變化。AUV以恒定轉(zhuǎn)速自航，推力隨速度增加而減小。