基于數值方法的AUV 近水面運動特性研究

王旭輝，林揚，王定前，孟令帥，高浩，曹新星，谷海濤

(1.東北大學 機械工程與自動化學院，遼寧 沈陽 110819；2.中國科學院沈陽自動化研究所機器人學國家重點實驗室，遼寧 沈陽 110016；3.中國科學院機器人與智能制造創新研究院，遼寧 沈陽 110169；4.中國人民解放軍32033 部隊，海南 海口 570000)

0 引言

自主式水下航行器（autonomous underwater vehicle，AUV）可自身攜帶能源，并配備多種傳感器，依靠自制能力在水下完成被賦予的各種任務[1]。AUV 自主性高，機動性好，在軍事和民用領域得到了廣泛應用。但受到目前能源技術水平的限制，AUV 不能長時間在水下作業，為及時補充能源和傳輸數據，需對其進行回收。目前AUV 回收方式可分為有人回收和無人回收2 種，有人回收方式主要依靠操作人員在母船上通過起吊裝置將AUV 起吊至甲板，而無人回收又可細分為水面回收和水下回收2 種方式[2]。無論是有人回收還是無人水面回收，回收過程的最終階段都需要AUV 在近水面航行，以便將AUV 回收至母船甲板上。為了提高回收成功率，需對AUV 近水面動力學特性進行研究。學者們對各種不同類型AUV 的動力學特性進行了大量研究工作，但大部分的研究是以深水場景為前提的[3– 6]。在近水面情況下，AUV 的運動受到自由液面的影響，會有與深水狀態不同的運動特性[7– 10]。

本文以AUV 的無人水面回收為背景，分析不同潛深情況下自由液面對AUV 運動的影響同時，結合水面回收場景，分析特定潛深下的AUV 水動力特性，并求解水動力系數，建立AUV 的近水面運動模型。最后，通過運動軌跡仿真對所建立的AUV 運動模型進行驗證。

1 模型

1.1 AUV 參數

研究對象為50 千克級AUV，按照美軍2000～2004 年間發布的《美海軍UUV 總體規劃》中給定的潛航器分級標準，將該小型AUV 稱為便攜式AUV[11]。便攜式AUV 的原始模型如圖1 所示。

圖1 便攜式AUV Fig.1 Portable AUV

便攜式AUV 主體上有天線、舵板以及推進器導流罩等，主體以外的裝置稱為附體。相對于AUV 主體來說，附體在CFD 計算中，附體會極大增加劃分網格的工作量，增加計算時間，降低計算精度，同時使計算結果難以收斂。為了解決上述問題，將AUV 原始模型進行適當簡化，在CFD 計算模型中，去除對計算結果影響較小但對計算難度影響較大的附體。簡化的AUV 模型如圖2 所示，習慣上稱之為光體。其體長為2 593 mm，直徑為250 mm。

圖2 簡化模型Fig.2 Simplified model of AUV

1.2 數值計算方法

粘性流動數值求解的本質是求解納維-斯托科斯（Navier–Stokes，N–S）方程[12]。根據雷諾數計算公式可以看出，AUV 近水面運動過程是一個湍流過程,目前湍流的數值模擬方法主要有直接數值模擬（direct numerical simulation,DNS）、大渦模擬（large eddy simulation,LES）、雷諾平均方程（reynolds-averaged navier-stokes equations，RANS）等。DNS 方法沒有對湍流流動做任何簡化，直接求解N–S 方程，理論上可以得到精確的計算結果。但在實際應用中，DNS 方法耗費巨大，目前僅能應用在低雷諾數的簡單湍流中；LES 方法對N–S 方程做了一定簡化，該方法認為，大尺度渦對流動的影響較大。通過N–S 方程直接求解，而小尺度渦對流動的影響較小，可通過模型對其進行刻畫。LES 方法在計算量上小于DNS 方法，但對于實際應用來說，該方法的計算量仍然偏大，需要高性能的計算設備才能實現；RANS 方法顧名思義，是對N–S 方程進行平均化處理，忽略了流動的細節，放棄模擬非定常流動，計算平均意義下的流動結果。N–S 中的非線性項平均化后會產生雷諾應力項，需引入湍流模型對其進行處理。RANS 方法計算量小，在工程實踐上得到了廣泛應用。考慮到時間成本與現有計算資源，選用RANS 方法進行便攜式AUV 的近水面運動模擬。

RANS 方程可寫為：

式中，τR為雷諾應力。

要求解RANS 方程，必須引入湍流模型，應用k?ε兩方程模型。湍流動能k和湍流耗散率 ε的傳輸方程可寫成如下形式：

1.3 AUV 運動方程

建立圖3 所示的慣性坐標系和隨體坐標系，并定義相關符號。E?ξηζ為慣性坐標系，固定于地球，不隨AUV 運動，可簡稱為定系；G?xyz為隨體坐標系，隨AUV 一同運動，可簡稱為動系，x，y，z為動系相對于定系的坐標值；φ，θ，ψ分別為動系相對于定系的姿態角（橫傾角、縱傾角、首向角）；V為AUV 在定系下的速度；?為AUV 在定系下的角速度；F為作用在AUV 上的力；M為作用在AUV 上的力矩。V，?，F，M又可表示為動系G?xyz上的分量，如表1 所示。

表1 矢量的各軸分量Tab.1 The components of each axis of the vector

圖3 坐標系示意圖Fig.3 Schematic of the coordinate system

在AUV 近水面運動過程中，假設動系G?xyz原點與重心位置重合，整個運動過程中AUV 整體浸沒于水下，重、浮心位置不變。忽略整個運動過程中AUV 的橫滾，根據潛艇標準運動方程[13]，便攜式AUV 六自由度運動方程無因次形式如下式：

1.4 仿真環境

采用CFD 軟件進行仿真計算，利用Q=h/D表示AUV 的下潛深度，h為AUV 軸線距自由液面距離，D為AUV 直徑如圖4 所示。L為AUV 軸向長度，計算域幾何尺寸如圖5 所示。計算域外形為立方體，長6L，寬2.4L，高3.5L，AUV 尾部指向面為壓力出口，其余5 個面設置為速度入口。

圖4 仿真場景示意Fig.4 Illustration of the simulation scene

圖5 計算域示意Fig.5 Illustration of the computed domain

1.5 網格無關性與精度

CFD 方法求解水動力問題需要消耗大量計算資源。本文采用的求解方法是RANS方法，對于RANS 方法，由于湍流模型的引入，當網格達到一定密度之后，繼續增加網格密度對計算結果影響不大。同時，由于離散誤差的存在，網格密度過大時誤差反而可能增大。為了合理利用計算資源，首先進行網格無關性驗證，即在誤差合理的情況下，用最粗糙的網格進行計算。為了保證更改網格密度后網格質量不變，最優方法是成倍數加密網格，但是對于三維問題來說，每次網格量需要至少增加8 倍，這種方法在實際操作中不可行。本文驗證網格無關性的方法如下：首先按照經驗大致劃分一套初始網格；再根據網格情況通過調整基礎尺寸的方式調整網格密度，將不同網格密度下的計算結果進行比較；選取最合理的一套網格為最終的計算網格。

針對潛深Q=1 的直航工況，通過調整基礎尺寸的方式生成5 套不同密度的計算網格。同時，通過對比不同網格密度下的相對誤差e來驗證網格無關性。e的定義如下式：

式中：ex為軸向力相對誤差；ez為垂向力相對誤差；X為軸向力；Z為垂向力；i為網格序號。

網格參數與計算結果如表2 所示。根據計算結果，同時考慮到實驗室計算資源配備情況，最終選取第4 套網格作為后續計算標準。

表2 網格無關性工況設置Tab.2 Grid-independent scenario settings

為了確定網格劃分方法的計算精度，采用REMUS 100 AUV 進行驗證。REMUS 100 AUV 外形尺寸與本文研究的AUV 相近，其在試驗中測得的阻力系數為cd=0.118347[14]。在上述網格參數設置的情況下進行REMUS 100 AUV 的阻力系數數值計算，得到結果為cd=0.124262，與試驗數據誤差為4.998%，滿足工程要求。

2 數值計算過程

x，y，z分別為AUV 隨體坐標系在慣性坐標系下的坐標；u，v，w分別為AUV 速度V在隨體坐標系3 個坐標軸上的投影；p，q，r分別為AUV 角速度 ?在隨體坐標系3 個軸上的投影；φ，θ，ψ分別為橫傾角、縱傾角和首向角。

2.1 自由液面對運動的影響

為了分析自由液面對AUV 運動的影響，分別設置潛深Q=1，1.5，2，2.5，3，3.5，4，航速0.5 kn，1 kn，1.5 kn，通過分析阻力變化探究自由液面對AUV 運動的影響。給定來流速度為1.0288 m/s，不同情況下的阻力變化如圖6 所示。

圖6 軸向力與垂向力變化情況Fig.6 Changes in axial and vertical forces

可以看出，自由液面會增大AUV 航行過程中的阻力，同時產生垂直向上的吸力。隨著潛深逐漸增加，自由液面所帶來的影響逐漸減小，在潛深達到4 倍直徑時，影響可忽略不計。

2.2 直航運動

2.3 PMM 運動

PMM 屬于拘束船模試驗的一種，通過強制AUV作規定運動，獲取AUV 受到的力與力矩，從而求得AUV的水動力系數。通過數值模擬橫蕩、首搖、縱蕩、縱搖4 種運動，求解相應的水動力系數，運動規律如下式：

式中：a為運動振幅；ω為運動頻率。

通過PMM 模擬方式得到的水動力系數會受到運動頻率的影響，在不同頻率下得到不同的數值，因此依靠單個工況獲得水動力系數是不合適的。目前常用的處理方法是，求解不同頻率下的水動力系數，之后通過過原點擬合或平均化處理的方式進行近似處理，得到較為合理的最終結果[15]。

為減小仿真工況對系數的影響，擬合不同周期下的系數，通過平均化處理的方式求得平均值。同時，通過擬合優度R2衡量擬合得到的水動力系數的適用性。R2定義式如下式：

式中，i為不同的仿真組數標號。可以看出，R2越接近1，表示水動力系數系數對原始數據的適用性越好。

3 結果與驗證

3.1 水動力系數

在潛深Q=1.5 工況下，不同周期下對應的水動力系數如表3 所示。

表3 水動力系數表Tab.3 Summary table of hydrodynamic coefficients

擬合精度如圖7 所示。

圖7 擬合精度匯總圖Fig.7 Fit the accuracy summary graph

可以看出，擬合精度R2絕大部分都在1 附近，證明水動力系數擬合精度良好。

3.2 軌跡仿真驗證

通過CFD 模擬直航運動軌跡與運動方程模擬直航運動軌跡的對比，驗證運動方程的準確性。CFD 直航運動模擬場景如圖8 所示，潛深Q=1.5，添加一階波浪。基于運動方程的直航運動模擬通過Simulink 實現，波浪力以操縱力的形式出入到運動方程中,仿真模型如圖9 所示[16，17]。

圖8 CFD 直航模擬場景Fig.8 CFD direct flight simulation scenario

圖9 運動方程模擬仿真模型Fig.9 Simulation model of the equation of motion simulation

基于運動方程的直航軌跡模擬與基于CFD 的直航軌跡模擬在3 個方向上的位移分量結果如圖10 所示。

圖10 直航軌跡仿真結果圖Fig.10 Trajectory simulation result graph

數據誤差統計情況如表4 所示。

表4 數據誤差統計表Tab.4 Statistical table of data errors

可以看出，X向與Z向數據誤差較小，Y向誤差較大，其原因是Y向位移基準值很小，導致計算出的誤差偏大。

4 結語

本文以水面自主回收AUV 為背景，對便攜式AUV 的近水面運動特性進行研究。在分析自由液面對AUV 運動影響的基礎上，針對其潛深Q=1.5 時的近水面水動力系數進行求解，并建立運動方程。最后，通過仿真運動軌跡方式對所得到的運動方程進行驗證。

結果表明，受到自由液面的影響，AUV 航行過程中會受到X軸負向的阻力和Z軸正向的升力，且影響程度隨著深度的增加逐漸減小。在深度超過4 倍AUV 直徑后，影響可以忽略不計。CFD 方法模擬直航軌跡與運動方程模擬直航軌跡得到的結果誤差較小，證明水動力系數及運動方程的準確性。