吳利紅，封錫盛，葉作霖，李一平

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

自主水下機(jī)器人(AUV)、潛艇通常需要快速下潛到一定深度執(zhí)行任務(wù)，精確預(yù)報(bào)這種下潛過程有助于對載體的航行軌跡和安全性進(jìn)行有效控制.模型試驗(yàn)和數(shù)值仿真是進(jìn)行載體下潛操縱性預(yù)報(bào)的主要方法.限于水池有限空間尺度，AUV下潛試驗(yàn)通常在海試中進(jìn)行，增加了風(fēng)險(xiǎn).傳統(tǒng)數(shù)值仿真方法是求解基于水動(dòng)力系數(shù)的泰勒展開的載體操縱性方程來預(yù)報(bào)垂直面操縱運(yùn)動(dòng).這種方法通過模型試驗(yàn)、數(shù)值模擬、面元法等獲得載體水動(dòng)力系數(shù)，量化載體控制面(包括舵翼)和螺旋槳的作用力和力矩，并作用于載體上，在MATLAB中的Simulink平臺中或VC平臺中實(shí)現(xiàn)載體離線的無動(dòng)力螺旋下潛或垂直下潛的操縱運(yùn)動(dòng)預(yù)報(bào)，獲得載體的運(yùn)動(dòng)軌跡、姿態(tài)、受力等量化參數(shù)[1-3].該方法無需建立載體幾何模型和求解流場，計(jì)算速度快.不足之處在于：① 較適用于以載體縱軸為主航行方向，其他方向運(yùn)動(dòng)較小的運(yùn)動(dòng)形式；② 運(yùn)動(dòng)受限于水動(dòng)力系數(shù)對應(yīng)的試驗(yàn)運(yùn)動(dòng)；③ 無法獲得載體的流場特征，無法探求物體運(yùn)動(dòng)響應(yīng)的內(nèi)在因素.

隨著計(jì)算流體力學(xué)軟件和硬件技術(shù)的發(fā)展，采用類物理數(shù)值模擬方法進(jìn)行載體自航操縱運(yùn)動(dòng)變得可能.這種操縱運(yùn)動(dòng)預(yù)報(bào)方法是直接建立載體帶舵翼和螺旋槳的全附體模型，求解流體動(dòng)力學(xué)方程，實(shí)時(shí)獲得載體的受力、運(yùn)動(dòng)，能獲得詳細(xì)的操縱運(yùn)動(dòng)流場特征.目前可以采用重疊網(wǎng)格法和動(dòng)網(wǎng)格法進(jìn)行這種類物理數(shù)值模擬，如文獻(xiàn)[4-7]采用重疊網(wǎng)格法對水面船舶和潛艇的自航操縱運(yùn)動(dòng)、z型操舵運(yùn)動(dòng)、波浪中的縱傾運(yùn)動(dòng)和水平面的回轉(zhuǎn)運(yùn)動(dòng)進(jìn)行了數(shù)值模擬.

相對于重疊網(wǎng)格，動(dòng)網(wǎng)格法是較早用于求解邊界移動(dòng)的一種方法，如對航空飛行器的武器分離仿真[8-9]，水下載體分離仿真如AUV 從魚雷管發(fā)射[10]和AUV水下對接[11].但是動(dòng)網(wǎng)格在邊界移動(dòng)時(shí)，常導(dǎo)致邊界附近流域的網(wǎng)格拉伸或壓縮比例失調(diào)，導(dǎo)致網(wǎng)格畸變，網(wǎng)格質(zhì)量變差，數(shù)值精度降低，甚至導(dǎo)致網(wǎng)格出現(xiàn)負(fù)體積而停止計(jì)算；此外，網(wǎng)格隨著邊界運(yùn)動(dòng)實(shí)時(shí)更新，網(wǎng)格數(shù)量隨著邊界移動(dòng)距離增加而急劇增加，尤其是大位移計(jì)算，網(wǎng)格總數(shù)增加導(dǎo)致計(jì)算效率降低，甚至無法繼續(xù)計(jì)算[12].這兩方面是動(dòng)網(wǎng)格發(fā)展的瓶頸，使動(dòng)網(wǎng)格在包含高速旋轉(zhuǎn)螺旋槳的類物理數(shù)值模擬方面落后于重疊網(wǎng)格.研究者開始探討采用移動(dòng)子區(qū)域代替移動(dòng)邊界法，解決動(dòng)網(wǎng)格法近場網(wǎng)格畸變問題[13-15].本文作者在已實(shí)現(xiàn)AUV自航運(yùn)動(dòng)的類物理數(shù)值模擬的基礎(chǔ)上[15-16]，將采用移動(dòng)子區(qū)域結(jié)合動(dòng)態(tài)層方法對AUV給定縱傾角變化和航速變化，而螺旋槳以恒定轉(zhuǎn)速旋轉(zhuǎn)的AUV自航下潛(以下簡稱強(qiáng)制自航[17])操縱運(yùn)動(dòng)進(jìn)行類物理數(shù)值模擬.

1 數(shù)值方法

1.1 運(yùn)動(dòng)描述

建立滿足右手系法則的大地坐標(biāo)系E-ξηζ和載體坐標(biāo)系G-xyz，如圖1所示.大地系下Eζ垂直指向地心.載體坐標(biāo)系原點(diǎn)與載體重心G重合，Gx沿AUV主對稱軸方向，向首為正.對應(yīng)大地系和載體系下的速度分量分別為v1,v2,v3和u,v,w.AUV通過調(diào)整縱傾實(shí)現(xiàn)下潛運(yùn)動(dòng)(見圖1)，分為3階段.第1階段，AUV通過打舵，使載體具有縱傾角速度，調(diào)整載體的縱傾角，準(zhǔn)備下潛；第2階段，載體在給定縱傾角下，以及螺旋槳推力作用下進(jìn)行定向直航下潛運(yùn)動(dòng)；第3階段，載體到達(dá)一定深度，回調(diào)縱傾角進(jìn)入定深直航.

圖1 AUV自航下潛示意圖Fig.1 AUV diving scheme by self-propulsion

本文將這3個(gè)階段AUV的運(yùn)動(dòng)提取出來，預(yù)定義AUV縱傾，水平和垂向速度變化，以及螺旋槳旋轉(zhuǎn)運(yùn)動(dòng)，模擬AUV強(qiáng)制自航下潛運(yùn)動(dòng).其中AUV的縱傾角變化，以縱傾角速度(ω2)給出，如圖2所示，在第1階段和第3階段，載體有縱傾角速度變化，第2階段為恒定縱傾角下潛，則縱傾角速度為0.

圖2 AUV下潛運(yùn)動(dòng)角速度變化Fig.2 Pitch rate during AUV diving

1.2 動(dòng)網(wǎng)格分區(qū)和網(wǎng)格模型

數(shù)值模擬的整個(gè)流域分區(qū)如圖3所示.流域?yàn)榉叫瘟饔?，共?個(gè)分區(qū).此方形域的中間包括核心區(qū)域C，AUV左側(cè)小流域?yàn)長，AUV右側(cè)小流域?yàn)镽.此方形域的外圍包括：上方域S1，下方域S2；左方域S3，右方域S4；前方域S5和后方域S6.網(wǎng)格模型為多塊混合網(wǎng)格，各部分的網(wǎng)格類型和運(yùn)動(dòng)形式如表1所示.C區(qū)域包括AUV載體、舵翼和螺旋槳，以及對應(yīng)的擾動(dòng)區(qū)域.此區(qū)域中，螺旋槳區(qū)域和尾跡區(qū)域采用非結(jié)構(gòu)網(wǎng)格離散.AUV和舵翼為結(jié)構(gòu)網(wǎng)格.除此外，C流域內(nèi)的體網(wǎng)格為非結(jié)構(gòu)網(wǎng)格，適用于AUV縱傾運(yùn)動(dòng)導(dǎo)致的網(wǎng)格變形和重構(gòu).其他外圍區(qū)域S1～S6為結(jié)構(gòu)網(wǎng)格.當(dāng)螺旋槳旋轉(zhuǎn)運(yùn)動(dòng)時(shí)，螺旋槳流域隨著螺旋槳旋轉(zhuǎn).AUV拖曳著螺旋槳和舵翼在C流域內(nèi)作縱傾旋轉(zhuǎn)運(yùn)動(dòng).當(dāng)AUV有縱向運(yùn)動(dòng)和垂向運(yùn)動(dòng)時(shí)，C、L和R隨AUV作縱向和垂向運(yùn)動(dòng).S1和S2作垂向運(yùn)動(dòng).外流域中運(yùn)動(dòng)流域的網(wǎng)格采用動(dòng)態(tài)層方法進(jìn)行更新.圖4給出載體自航下潛的網(wǎng)格圖，包括AUV在垂直對稱面上的網(wǎng)格圖和尾部局部放大圖.

圖3 AUV 下潛運(yùn)動(dòng)流場拓?fù)浣Y(jié)構(gòu)Fig.3 Domain topology of AUV diving

表1 網(wǎng)格類型和運(yùn)動(dòng)形式Tab.1 Mesh type and motion pattern

圖4 AUV 網(wǎng)格系統(tǒng)Fig.4 AUV grid system

1.3 用戶自定義函數(shù)(UDF)

AUV強(qiáng)制自航下潛運(yùn)動(dòng)通過編寫UDF實(shí)現(xiàn)，如圖5所示的流程圖.其中n1,n2和n3表示螺旋槳轉(zhuǎn)速繞Eξ,Eη,Eζ3軸的角速度分量；R1和R3為AUV在Gx,Gz的阻力分量；T1，T3為螺旋槳沿Gx,Gz的推力分量.此流程圖包含3個(gè)模塊：螺旋槳模塊、AUV模塊和流域模塊.在螺旋槳模塊中，根據(jù)AUV縱傾角，計(jì)算螺旋槳旋轉(zhuǎn)速度分量，計(jì)算螺旋槳推力，并將此推力傳遞給AUV；在AUV模塊，計(jì)算AUV速度在大地下的分量，并讓AUV以縱傾角速度ω2和線速度v1,v3運(yùn)動(dòng)，求解流體力學(xué)方程，獲得局部坐標(biāo)下的速度分量，并將速度分量進(jìn)行坐標(biāo)轉(zhuǎn)換；在流域模塊中，局部模塊隨著AUV作縱向和垂向平移運(yùn)動(dòng)，遠(yuǎn)場水平流域也做縱向和垂向運(yùn)動(dòng)，遠(yuǎn)場垂直流域作垂向運(yùn)動(dòng)，各界面采用非一致連接實(shí)現(xiàn)界面的變量傳遞.

圖5 AUV強(qiáng)制自航下潛運(yùn)動(dòng)流程圖Fig.5 Flowchart of AUV forced diving motion

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

圖6 敞水試驗(yàn)的試驗(yàn)結(jié)果與數(shù)值結(jié)果對比Fig.6 Comparison of experimental and numerical results of open water curve

圖7 AUV自航試驗(yàn)速度對比Fig.7 Comparison of computational and experimental data of AUV’s velocity in self-propulsion

3 數(shù)值結(jié)果

數(shù)值模擬AUV自航下潛過程，AUV以一定的速度和角速度開始下潛運(yùn)動(dòng).初始速度為接近0的一個(gè)非常小的數(shù)，待定常計(jì)算收斂后，開始以此為初始值，AUV調(diào)整縱傾下潛.計(jì)算時(shí)間步長為螺旋槳旋轉(zhuǎn)1° 所對應(yīng)的時(shí)間，垂向航行距離約1個(gè)載體長度，縱向航行距離約5個(gè)載體長度，對應(yīng)的物理時(shí)間接近7 s.圖8給出4個(gè)典型時(shí)刻(t=0.1,3.0,6.0，6.7 s) AUV下潛過程中垂直面對應(yīng)的網(wǎng)格圖.載體在垂直下潛過程中有明顯的縱傾變化，同時(shí)具有縱向和垂向運(yùn)動(dòng).在這過程中，網(wǎng)格質(zhì)量仍然保持較好.

圖8 不同時(shí)刻網(wǎng)格圖(η=0平面)Fig.8 Meshing at different times (plane of η=0)

整個(gè)下潛過程AUV受力和螺旋槳推力分別如圖9和10所示.由圖9可見，AUV調(diào)整縱傾，從靜止開始快速加速到最大速度，則AUV在Eξ和Eζ具有較大加速度，對應(yīng)的阻力R1,R3都較大.尤其是垂向阻力，在載體快速改變縱傾時(shí)，對應(yīng)的垂向阻尼變化顯著且振蕩.這在載體改變縱傾回航階段也體現(xiàn)明顯.這主要是由于載體縱傾角速度過大引起的.當(dāng)載體處于恒定縱傾角定向下潛時(shí)，阻尼趨于穩(wěn)定.推力方面，載體在強(qiáng)制下潛和回航調(diào)整縱傾時(shí)，對螺旋槳推力影響較大，縱向推力和垂向推力都變化顯著.載體無縱傾變化時(shí)，推力也趨于恒定值.因此，對于自航直航或自航下潛運(yùn)動(dòng)，航速應(yīng)該由靜止開始慢慢增加，否則載體加速度太大，流場擾動(dòng)影響較大，趨于穩(wěn)定的時(shí)間增加，則不僅需要更大的力達(dá)到對應(yīng)的航速，流場也需要更長的預(yù)行段長度才能趨于穩(wěn)定.

圖9 AUV 在縱向和垂向的阻力Fig.9 AUV resistances along longitudinal and vertical directions

圖10 螺旋槳在縱向和垂向的推力Fig.10 Propeller thrusts along longitudinal and vertical directions

圖11為AUV下潛不同時(shí)刻的AUV尾跡(v1)云圖(注：遠(yuǎn)場速度為0).由圖可見，載體縱傾變化，將導(dǎo)致AUV拖曳的螺旋槳尾跡發(fā)生變化.當(dāng)AUV繞著載體重心出現(xiàn)埋首下潛時(shí)，AUV尾部則繞重心往上運(yùn)動(dòng)，由圖11(a)可見螺旋槳尾跡往上扭曲，呈現(xiàn)上甩現(xiàn)象；當(dāng)載體無縱傾變化時(shí)(t=0.5,3.2,6.1 s)，螺旋槳尾跡沿著載體縱軸方向；當(dāng)載體達(dá)到一定深度，回調(diào)縱傾到直航時(shí)，載體抬首，對應(yīng)的載體尾部則向下運(yùn)動(dòng)，螺旋槳尾跡中可見明顯的下擺現(xiàn)象.此外，螺旋槳尾跡中可見螺旋槳旋轉(zhuǎn)運(yùn)動(dòng)有明顯的梢渦泄出，以及與梢渦運(yùn)動(dòng)方向相反的轂渦.

圖11 不同時(shí)刻AUV尾部速度放大圖Fig.11 Closed-up velocity contours at AUV’s stern at different times

4 結(jié)論

海洋載體或空中飛行器的類物理數(shù)值模擬研究是數(shù)值模擬研究的最終目標(biāo)，這種方法可以用數(shù)值模擬完全代替復(fù)雜、昂貴、風(fēng)險(xiǎn)性極高、周期長的試驗(yàn)，有利于以較小的成本設(shè)計(jì)出性能更優(yōu)的載體，預(yù)報(bào)出已設(shè)計(jì)載體復(fù)雜操縱運(yùn)動(dòng)的各種動(dòng)態(tài)特性.隨著計(jì)算流體力學(xué)的發(fā)展和超級計(jì)算機(jī)計(jì)算能力的提高，這種類物理數(shù)值模擬變得可能，已經(jīng)可以對一些載體的操縱運(yùn)動(dòng)進(jìn)行類物理數(shù)值模擬.實(shí)現(xiàn)這種類物理數(shù)值模擬的關(guān)鍵技術(shù)是預(yù)報(bào)的精度和計(jì)算效率.可以采用超級計(jì)算機(jī)提高計(jì)算效率，但是網(wǎng)格技術(shù)仍然是這種數(shù)值模擬的關(guān)鍵點(diǎn).

本文基于多塊動(dòng)態(tài)混合網(wǎng)格及動(dòng)態(tài)層區(qū)域法實(shí)現(xiàn)了AUV垂直面內(nèi)自航下潛操縱運(yùn)動(dòng)的類物理數(shù)值模擬.相對于非結(jié)構(gòu)動(dòng)網(wǎng)格方法和動(dòng)態(tài)重疊網(wǎng)格法(前者需要20～30 d的時(shí)間[16]；后者進(jìn)行自航數(shù)值模擬在高性能機(jī)上通過減少質(zhì)量進(jìn)行加速計(jì)算，仍需1～2個(gè)月[7]，且只能給出定常運(yùn)動(dòng)結(jié)果)而言，該方法計(jì)算速度快，精度高.數(shù)值模擬不需要減小質(zhì)量，能獲得載體自航下潛運(yùn)動(dòng)全程的計(jì)算結(jié)果，使數(shù)值計(jì)算在普通臺式機(jī)(i5-6400 CPU @2.70 GHz，2.70 GHz，內(nèi)存16.0 GB)上以4個(gè)處理器并行計(jì)算，歷時(shí)12 d左右完成.本文為水下載體或空中飛行器的復(fù)雜多自由度耦合的操縱運(yùn)動(dòng)的類物理數(shù)值模擬提供了方法.未來將類物理數(shù)值模擬擴(kuò)展到更多自由度耦合，更加復(fù)雜的海洋環(huán)境下的載體操縱運(yùn)動(dòng)預(yù)報(bào)中.