王旭輝，林揚(yáng)，王定前，孟令帥，高浩，曹新星，谷海濤

(1.東北大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院，遼寧 沈陽 110819；2.中國(guó)科學(xué)院沈陽自動(dòng)化研究所機(jī)器人學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 沈陽 110016；3.中國(guó)科學(xué)院機(jī)器人與智能制造創(chuàng)新研究院，遼寧 沈陽 110169；4.中國(guó)人民解放軍32033 部隊(duì)，海南 ?？?570000)

0 引言

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

本文以AUV 的無人水面回收為背景，分析不同潛深情況下自由液面對(duì)AUV 運(yùn)動(dòng)的影響同時(shí)，結(jié)合水面回收?qǐng)鼍?，分析特定潛深下的AUV 水動(dòng)力特性，并求解水動(dòng)力系數(shù)，建立AUV 的近水面運(yùn)動(dòng)模型。最后，通過運(yùn)動(dòng)軌跡仿真對(duì)所建立的AUV 運(yùn)動(dòng)模型進(jìn)行驗(yàn)證。

1 模型

1.1 AUV 參數(shù)

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

圖1 便攜式AUV Fig.1 Portable AUV

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

圖2 簡(jiǎn)化模型Fig.2 Simplified model of AUV

1.2 數(shù)值計(jì)算方法

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

RANS 方程可寫為：

式中，τR為雷諾應(yīng)力。

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

1.3 AUV 運(yùn)動(dòng)方程

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

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

圖3 坐標(biāo)系示意圖Fig.3 Schematic of the coordinate system

在AUV 近水面運(yùn)動(dòng)過程中，假設(shè)動(dòng)系G?xyz原點(diǎn)與重心位置重合，整個(gè)運(yùn)動(dòng)過程中AUV 整體浸沒于水下，重、浮心位置不變。忽略整個(gè)運(yùn)動(dòng)過程中AUV 的橫滾，根據(jù)潛艇標(biāo)準(zhǔn)運(yùn)動(dòng)方程[13]，便攜式AUV 六自由度運(yùn)動(dòng)方程無因次形式如下式：

1.4 仿真環(huán)境

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

圖4 仿真場(chǎng)景示意Fig.4 Illustration of the simulation scene

圖5 計(jì)算域示意Fig.5 Illustration of the computed domain

1.5 網(wǎng)格無關(guān)性與精度

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

針對(duì)潛深Q=1 的直航工況，通過調(diào)整基礎(chǔ)尺寸的方式生成5 套不同密度的計(jì)算網(wǎng)格。同時(shí)，通過對(duì)比不同網(wǎng)格密度下的相對(duì)誤差e來驗(yàn)證網(wǎng)格無關(guān)性。e的定義如下式：

式中：ex為軸向力相對(duì)誤差；ez為垂向力相對(duì)誤差；X為軸向力；Z為垂向力；i為網(wǎng)格序號(hào)。

網(wǎng)格參數(shù)與計(jì)算結(jié)果如表2 所示。根據(jù)計(jì)算結(jié)果，同時(shí)考慮到實(shí)驗(yàn)室計(jì)算資源配備情況，最終選取第4 套網(wǎng)格作為后續(xù)計(jì)算標(biāo)準(zhǔn)。

表2 網(wǎng)格無關(guān)性工況設(shè)置Tab.2 Grid-independent scenario settings

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

2 數(shù)值計(jì)算過程

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

2.1 自由液面對(duì)運(yùn)動(dòng)的影響

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

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

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

2.2 直航運(yùn)動(dòng)

2.3 PMM 運(yùn)動(dòng)

PMM 屬于拘束船模試驗(yàn)的一種，通過強(qiáng)制AUV作規(guī)定運(yùn)動(dòng)，獲取AUV 受到的力與力矩，從而求得AUV的水動(dòng)力系數(shù)。通過數(shù)值模擬橫蕩、首搖、縱蕩、縱搖4 種運(yùn)動(dòng)，求解相應(yīng)的水動(dòng)力系數(shù)，運(yùn)動(dòng)規(guī)律如下式：

式中：a為運(yùn)動(dòng)振幅；ω為運(yùn)動(dòng)頻率。

通過PMM 模擬方式得到的水動(dòng)力系數(shù)會(huì)受到運(yùn)動(dòng)頻率的影響，在不同頻率下得到不同的數(shù)值，因此依靠單個(gè)工況獲得水動(dòng)力系數(shù)是不合適的。目前常用的處理方法是，求解不同頻率下的水動(dòng)力系數(shù)，之后通過過原點(diǎn)擬合或平均化處理的方式進(jìn)行近似處理，得到較為合理的最終結(jié)果[15]。

為減小仿真工況對(duì)系數(shù)的影響，擬合不同周期下的系數(shù)，通過平均化處理的方式求得平均值。同時(shí)，通過擬合優(yōu)度R2衡量擬合得到的水動(dòng)力系數(shù)的適用性。R2定義式如下式：

式中，i為不同的仿真組數(shù)標(biāo)號(hào)?？梢钥闯觯琑2越接近1，表示水動(dòng)力系數(shù)系數(shù)對(duì)原始數(shù)據(jù)的適用性越好。

3 結(jié)果與驗(yàn)證

3.1 水動(dòng)力系數(shù)

在潛深Q=1.5 工況下，不同周期下對(duì)應(yīng)的水動(dòng)力系數(shù)如表3 所示。

表3 水動(dòng)力系數(shù)表Tab.3 Summary table of hydrodynamic coefficients

擬合精度如圖7 所示。

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

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

3.2 軌跡仿真驗(yàn)證

通過CFD 模擬直航運(yùn)動(dòng)軌跡與運(yùn)動(dòng)方程模擬直航運(yùn)動(dòng)軌跡的對(duì)比，驗(yàn)證運(yùn)動(dòng)方程的準(zhǔn)確性。CFD 直航運(yùn)動(dòng)模擬場(chǎng)景如圖8 所示，潛深Q=1.5，添加一階波浪?；谶\(yùn)動(dòng)方程的直航運(yùn)動(dòng)模擬通過Simulink 實(shí)現(xiàn)，波浪力以操縱力的形式出入到運(yùn)動(dòng)方程中,仿真模型如圖9 所示[16，17]。

圖8 CFD 直航模擬場(chǎng)景Fig.8 CFD direct flight simulation scenario

圖9 運(yùn)動(dòng)方程模擬仿真模型Fig.9 Simulation model of the equation of motion simulation

基于運(yùn)動(dòng)方程的直航軌跡模擬與基于CFD 的直航軌跡模擬在3 個(gè)方向上的位移分量結(jié)果如圖10 所示。

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

數(shù)據(jù)誤差統(tǒng)計(jì)情況如表4 所示。

表4 數(shù)據(jù)誤差統(tǒng)計(jì)表Tab.4 Statistical table of data errors

可以看出，X向與Z向數(shù)據(jù)誤差較小，Y向誤差較大，其原因是Y向位移基準(zhǔn)值很小，導(dǎo)致計(jì)算出的誤差偏大。

4 結(jié)語

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

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