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

(1.大連海事大學 船舶與海洋工程學院，遼寧 大連 116026；2.中國科學院沈陽自動化研究所 機器人學國家重點實驗室，沈陽 110016)

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

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

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

1 數(shù)值方法

1.1 運動描述

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

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

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

圖2 AUV下潛運動角速度變化Fig.2 Pitch rate during AUV diving

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

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

圖3 AUV 下潛運動流場拓撲結構Fig.3 Domain topology of AUV diving

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

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

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

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

圖5 AUV強制自航下潛運動流程圖Fig.5 Flowchart of AUV forced diving motion

2 數(shù)值驗證

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

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

3 數(shù)值結果

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

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

整個下潛過程AUV受力和螺旋槳推力分別如圖9和10所示.由圖9可見，AUV調整縱傾，從靜止開始快速加速到最大速度，則AUV在Eξ和Eζ具有較大加速度，對應的阻力R1,R3都較大.尤其是垂向阻力，在載體快速改變縱傾時，對應的垂向阻尼變化顯著且振蕩.這在載體改變縱傾回航階段也體現(xiàn)明顯.這主要是由于載體縱傾角速度過大引起的.當載體處于恒定縱傾角定向下潛時，阻尼趨于穩(wěn)定.推力方面，載體在強制下潛和回航調整縱傾時，對螺旋槳推力影響較大，縱向推力和垂向推力都變化顯著.載體無縱傾變化時，推力也趨于恒定值.因此，對于自航直航或自航下潛運動，航速應該由靜止開始慢慢增加，否則載體加速度太大，流場擾動影響較大，趨于穩(wěn)定的時間增加，則不僅需要更大的力達到對應的航速，流場也需要更長的預行段長度才能趨于穩(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下潛不同時刻的AUV尾跡(v1)云圖(注：遠場速度為0).由圖可見，載體縱傾變化，將導致AUV拖曳的螺旋槳尾跡發(fā)生變化.當AUV繞著載體重心出現(xiàn)埋首下潛時，AUV尾部則繞重心往上運動，由圖11(a)可見螺旋槳尾跡往上扭曲，呈現(xiàn)上甩現(xiàn)象；當載體無縱傾變化時(t=0.5,3.2,6.1 s)，螺旋槳尾跡沿著載體縱軸方向；當載體達到一定深度，回調縱傾到直航時，載體抬首，對應的載體尾部則向下運動，螺旋槳尾跡中可見明顯的下擺現(xiàn)象.此外，螺旋槳尾跡中可見螺旋槳旋轉運動有明顯的梢渦泄出，以及與梢渦運動方向相反的轂渦.

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

4 結論

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

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