水下滑翔機(jī)附加質(zhì)量數(shù)值計(jì)算

楊 磊，曹軍軍，姚寶恒，曾 錚，連 璉

(上海交通大學(xué) 船舶與海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

水下滑翔機(jī)附加質(zhì)量數(shù)值計(jì)算

楊 磊，曹軍軍，姚寶恒，曾 錚，連 璉

(上海交通大學(xué) 船舶與海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

水下滑翔機(jī)的附加質(zhì)量對(duì)其運(yùn)動(dòng)狀態(tài)影響較大，本文采用了 Hess-Smith（面元法）方法編制了計(jì)算任意三維物體附加質(zhì)量的程序，利用 Gambit 軟件對(duì)水下滑翔機(jī)進(jìn)行網(wǎng)格劃分，計(jì)算出實(shí)驗(yàn)室所研制的水下滑翔機(jī)附加質(zhì)量，同時(shí)利用 CFD（Computational Fluid Dynamics）技術(shù)，結(jié)合動(dòng)網(wǎng)格技術(shù)和 UDF（User-Defined Function），對(duì)水下滑翔機(jī)進(jìn)行了 PMM（Planar Motion Mechanism）試驗(yàn)?zāi)M，并與 Hess-Smith 方法得到的結(jié)果進(jìn)行對(duì)比，分析兩者之間的特點(diǎn)和各自優(yōu)勢(shì)。

附加質(zhì)量；水下滑翔機(jī)；面元法；PMM

0 引 言

水下滑翔機(jī)是通過改變自身重心和凈重力來進(jìn)行水下運(yùn)動(dòng)，本身沒有外部推進(jìn)裝置。它具有低能耗、低成本、適應(yīng)性強(qiáng)等特點(diǎn)，進(jìn)而可以完成長(zhǎng)期大范圍的調(diào)研任務(wù)。它在物理海洋、生物海洋等領(lǐng)域具有廣泛的應(yīng)用空間[1]。

隨著國(guó)內(nèi)外對(duì)深潛器研究不斷深入，對(duì)其操縱性預(yù)報(bào)要求越來越高，而潛器的慣性力水動(dòng)力可由附加質(zhì)量進(jìn)行推算，因此精確的附加質(zhì)量計(jì)算對(duì)深潛器的設(shè)計(jì)十分重要。國(guó)內(nèi)外對(duì)附加質(zhì)量進(jìn)行了大量的研究，目前獲取附加質(zhì)量有以下幾種方法，對(duì)于規(guī)則形狀物體的附加質(zhì)量可以利用圖譜進(jìn)行查找[2]，對(duì)于復(fù)雜物體的附加質(zhì)量可以使用 Hess-Smith 方法[3]，以求解出物體無(wú)升力繞流場(chǎng)的數(shù)值解。此外，也可以通過實(shí)驗(yàn)的方式得到物體的附加質(zhì)量，但其成本較高，適用范圍窄。隨著計(jì)算流體力學(xué)的飛速發(fā)展，CFD 方法也逐漸成為預(yù)報(bào)附加質(zhì)量的有效方法。在船舶領(lǐng)域、潛艇領(lǐng)域、水陸兩棲車領(lǐng)域[4-6]等均有該方法的使用。但目前對(duì)水下滑翔機(jī)附加質(zhì)量的研究較少，多數(shù)文獻(xiàn)采用簡(jiǎn)化模型近似求解甚至直接使用相似滑翔機(jī)的結(jié)果，對(duì)運(yùn)動(dòng)仿真造成一定的偏差。

為了精確求解實(shí)驗(yàn)室所研制的水下滑翔機(jī)附加質(zhì)量，本文使用了 Hess-Smith 方法設(shè)計(jì)程序，同時(shí)利用Gambit 軟件對(duì)水下滑翔機(jī)表面進(jìn)行網(wǎng)格劃分，獲取滑翔機(jī)表面節(jié)點(diǎn)信息，導(dǎo)入程序中求解矩陣從而得到結(jié)果；同時(shí)，本文利用 CFD 軟件對(duì)水下滑翔機(jī)的 PMM試驗(yàn)進(jìn)行模擬，利用動(dòng)網(wǎng)格技術(shù)實(shí)現(xiàn)網(wǎng)格時(shí)時(shí)更新，UDF 實(shí)現(xiàn)了對(duì)滑翔機(jī)純垂蕩運(yùn)動(dòng)、純橫蕩運(yùn)動(dòng)、純首搖運(yùn)動(dòng)和純俯仰運(yùn)動(dòng) 4 種運(yùn)動(dòng)控制，通過對(duì)所受到的力和力矩的分析處理，得到附加質(zhì)量。通過對(duì)比 Hess-Smith 法和 CFD 方法得到的附加質(zhì)量，驗(yàn)證了本文使用的 2 種計(jì)算方法可行性和準(zhǔn)確性，進(jìn)而對(duì)比兩者之間各自優(yōu)劣，為對(duì)以后從事相關(guān)計(jì)算的人員提供參考意見和方法。

1 基于面元法的附加質(zhì)量計(jì)算

1.1 基本方程和邊界條件

在不可壓縮勢(shì)流中，定常速度勢(shì)（Φ）在物體外部空間域中適合拉普拉斯方程，在物面上適合不可進(jìn)入條件，在無(wú)窮遠(yuǎn)處，應(yīng)與均勻來流的速度勢(shì)吻合。在速度勢(shì) Φ 中分出已知的均勻來流項(xiàng)，即可得到擾動(dòng)速度勢(shì) φ 的定解條件：

其中：式（2）在物體外部；式（3）在物體表面上；式（4）在無(wú)窮遠(yuǎn)處；V 為物體運(yùn)動(dòng)速度，n 為物面單位法線向量，指向物體內(nèi)部。

1.2 速度勢(shì)離散

假設(shè)流場(chǎng)中的場(chǎng)點(diǎn) p（x, y, z），物面上的源點(diǎn)為q，在物面上布置強(qiáng)度為 σ（q）的 Rankine 源，則流場(chǎng)內(nèi)任意點(diǎn)的速度勢(shì)為

其中：rpq為場(chǎng)點(diǎn) p 和源點(diǎn) q 之間的距離；s 為物面。

在式（5）兩端取法向?qū)?shù)，并令點(diǎn) p 沿法線方向趨于物面上的點(diǎn)。式（5）右端分布源的法向?qū)?shù)由兩部分組成，一部分是 P 點(diǎn)附近小曲面 ε 貢獻(xiàn)的，另一部分是物面其他部分 s-ε 貢獻(xiàn)的。同時(shí)代入式（3）的邊界條件，可以得到：

將積分式（6）轉(zhuǎn)換成線性代數(shù)方程組，即把物面s 分成 N 小塊，小塊形狀可為三角形或四邊形，近似代替小曲面（Δsj）。本文采用三角形劃分，利用右手法則進(jìn)行編號(hào)，得到單位向量，再通過三角形平均中心點(diǎn)進(jìn)行投影，得到新的三角形 ΔQj代替原來的小曲面。因此，物面 s 上的積分可以用 N 個(gè)三角形上的積分之和近似，即：

將式（7）代入式（6）中，可以得到 σ 的 N 階線性代數(shù)方程組：

動(dòng)點(diǎn) p 是 N 個(gè)單元的中心點(diǎn)，通過 Gambit 軟件，對(duì)物體進(jìn)行網(wǎng)格劃分，導(dǎo)出網(wǎng)格數(shù)據(jù)后，將節(jié)點(diǎn)信息導(dǎo)入程序中，即可求解出 σj的值，再通過式（5）得到速度勢(shì) φ。本文將水下滑翔機(jī)表面進(jìn)行網(wǎng)格劃分，其中網(wǎng)格為三角形，共計(jì) 907 個(gè)，利用各網(wǎng)格節(jié)點(diǎn)坐標(biāo)生成近似圖，觀察網(wǎng)格劃分質(zhì)量，如圖 1 所示。

圖 1 水下滑翔機(jī)三角形網(wǎng)格劃分Fig. 1 The generation of triangular mesh for an underwater glider

1.3 附加質(zhì)量求解

當(dāng)一個(gè)剛體在理想流體中做非定常運(yùn)動(dòng)時(shí)，會(huì)帶動(dòng)周圍流體做加速運(yùn)動(dòng)，其所受到的流體力的大小與運(yùn)動(dòng)加速度成正比，方向則與加速度方向相反，這個(gè)比例系數(shù)即為附加質(zhì)量 λij。一個(gè)任意形狀的物體運(yùn)動(dòng)時(shí)，會(huì)產(chǎn)生 36 個(gè)附加質(zhì)量。

根據(jù)勢(shì)流理論[7]，可知：

由于水下滑翔機(jī)相對(duì)于 xoz 平面左右對(duì)稱，則附加質(zhì)量中所有 i + j 為奇數(shù)的系數(shù)全部為 0。同時(shí)，本滑翔機(jī)機(jī)翼相對(duì)于 xoy 平面對(duì)稱，僅尾部上端布置了天線和平衡舵，因此相對(duì)于 xoy 平面是大致上下對(duì)稱的，則 λ13和 λ15很小，可以忽略。因此附加質(zhì)量的矩陣形式為：

按照勢(shì)流理論知物體附加質(zhì)量為：

式中，i, j=1, 2……6, φ1,2,3分別為剛體沿動(dòng)系坐標(biāo)軸 x, y, z 以單位速度平動(dòng)時(shí)引起的流體速度勢(shì)，φ4,6,6則分別為剛體以單位角速度繞 x, y, z 軸以單位角度轉(zhuǎn)動(dòng)時(shí)所引起的流體速度勢(shì)。因?yàn)閿_動(dòng)速度勢(shì) φ 已得，因此可求出物體的附加質(zhì)量。

本文對(duì)水下滑翔機(jī)進(jìn)行了網(wǎng)格劃分，同時(shí)使用Fortran 軟件編寫程序，將網(wǎng)格節(jié)點(diǎn)信息導(dǎo)入后，即可得到水下滑翔機(jī)附加質(zhì)量結(jié)果，如表 1 所示。

表 1 水下滑翔機(jī)附加質(zhì)量Tab. 1 The added mass of an underwater glider

2 基于 CFD 軟件的附加質(zhì)量計(jì)算

Fluent 軟件可提供動(dòng)網(wǎng)格技術(shù)，動(dòng)網(wǎng)格模型可以用來模擬流程形狀由于邊界運(yùn)動(dòng)而隨時(shí)間改變的問題，利用 UDF 可以定義邊界的運(yùn)動(dòng)方式。本文則是對(duì)水下滑翔機(jī)進(jìn)行 PMM 試驗(yàn)?zāi)M，從而得到附加質(zhì)量。

2.1 計(jì)算方程

水下滑翔機(jī) PMM 試驗(yàn)主要包括：純升沉運(yùn)動(dòng)、純橫蕩運(yùn)動(dòng)、純俯仰運(yùn)動(dòng)、純搖首運(yùn)動(dòng)[8]。由于純升沉運(yùn)動(dòng)與純橫蕩運(yùn)動(dòng)類似，純俯仰運(yùn)動(dòng)與純搖首運(yùn)動(dòng)類似，本文僅寫出前兩者的運(yùn)動(dòng)方程和處理方法。

對(duì)于純升沉運(yùn)動(dòng)，其運(yùn)動(dòng)方程為：

其中：ζ 為滑翔機(jī)垂向位移；a 為滑翔機(jī)純升沉運(yùn)動(dòng)的振幅；分別為水下滑翔機(jī)縱傾角度和角速度；分別為水下滑翔機(jī)垂向速度和加速度。根據(jù)運(yùn)動(dòng)方程可知，滑翔機(jī)做純升沉運(yùn)動(dòng)的軌跡如圖 2 所示。

圖 2 水下滑翔機(jī)純升沉運(yùn)動(dòng)示意圖Fig. 2 The heaving motion of an underwater glider

根據(jù)水動(dòng)力系數(shù)的含義可得到方程：

式中：Z 為運(yùn)動(dòng)過程中滑翔機(jī)受到的垂向力；M 為繞 y軸的力矩。

對(duì)于純俯仰運(yùn)動(dòng)，其運(yùn)動(dòng)方程為：

式中：θ0為滑翔機(jī)純俯仰運(yùn)動(dòng)縱傾角振幅，q 為角速度。

根據(jù)運(yùn)動(dòng)方程可知，滑翔機(jī)做純俯仰運(yùn)動(dòng)的軌跡如圖 3 所示。為了滿足純俯仰的運(yùn)動(dòng)條件，即模型縱軸時(shí)時(shí)保持與原點(diǎn)軌跡相切，可得到方程：

圖 3 水下滑翔機(jī)純俯仰運(yùn)動(dòng)示意圖Fig. 3 The pitching motion of an underwater glider

因?yàn)榛铏C(jī)的 PMM 試驗(yàn)為小振幅運(yùn)動(dòng)，因而可以近似處理，其中，U 為滑翔機(jī) x 方向上的來流速度。

2.2 CFD 計(jì)算準(zhǔn)備

精確的計(jì)算結(jié)果和高質(zhì)量的網(wǎng)格密不可分，由于水下滑翔機(jī)外形較為復(fù)雜，而四面體網(wǎng)格具有較好的貼體性，適應(yīng)各類復(fù)雜形體，同時(shí)結(jié)構(gòu)化網(wǎng)格質(zhì)量高，網(wǎng)格數(shù)量少，結(jié)果精確，節(jié)省運(yùn)算時(shí)間。因而將計(jì)算域分為內(nèi)域和外域 2 個(gè)部分，內(nèi)域采用非結(jié)構(gòu)網(wǎng)格，而外域采用結(jié)構(gòu)化網(wǎng)格。兼得兩者的優(yōu)點(diǎn)，減少計(jì)算時(shí)間，提高精度。

圖 4為計(jì)算內(nèi)、外域網(wǎng)格劃分和邊界條件說明圖，網(wǎng)格總數(shù)量約為 260 萬(wàn)，其中入口邊界給定為VELOCITY_INLET，出口邊界為 OUTFLOW，水下滑翔機(jī)表面設(shè)為 WALL，壁面設(shè)為 VELOCITY_INLET 且速度與入口邊界相同。

湍流模型選取 SST k-ω 模型，壓力速度耦合選取SIMPLEC 算法，壓力插值采取 Standard 格式，動(dòng)量、湍動(dòng)能和湍流耗散率插值采用 Second Order Upwind格式。

在 Fluent 軟件中，需要編制以描述水下滑翔機(jī)運(yùn)動(dòng)的用戶自定義函數(shù)（UDF），本文根據(jù)水下滑翔機(jī)PMM 試驗(yàn)所對(duì)應(yīng)的運(yùn)動(dòng)方程如式（12）和（14）來編制 UDF 程序，其中 DEFINE_CG_MOTION 函數(shù)控制滑翔機(jī)的運(yùn)動(dòng)，COMPUTE_FORCE_AND_MOMENT 函數(shù)用于得到滑翔機(jī)每一時(shí)刻的力和力矩，并通過FPRINTF 函數(shù)記錄和存儲(chǔ)在文件中。

2.3 結(jié)果處理與分析

以純升沉運(yùn)動(dòng)為例，設(shè)定振幅 a = 300 mm，來流速度 v = 0.3 m/s，滑翔機(jī)長(zhǎng)度為 L = 2 m，海水密度 ρ = 1 025 kg/m3。為了使迭代時(shí)間歩長(zhǎng)取整數(shù)，頻率 f 設(shè)置為：0.2, 0.4, 0.5, 0.625。以 f 為 0.4 時(shí)為例，周期是 2.5 s，因此 fluent 中設(shè)置計(jì)算步長(zhǎng)為 0.01 s，每個(gè)周期需要250 步完成，計(jì)算 3～5 個(gè)周期后結(jié)束，利用數(shù)據(jù)的結(jié)果，在 Origin 軟件中進(jìn)行處理。取時(shí)間靠后結(jié)果較為穩(wěn)定的數(shù)據(jù)進(jìn)行擬合，如圖 5 和圖 6 所示。

圖 4 計(jì)算域內(nèi)域和外域網(wǎng)格劃分Fig. 4 The mesh of internal and external flow domain

圖 5 純升沉運(yùn)動(dòng) f = 0.4 時(shí)垂向力 Z 的變化Fig. 5 The changes of force in Z-direction when the frequency of heaving motion is 0.4

圖 6 純升沉運(yùn)動(dòng) f = 0.4 時(shí)繞 Z 軸力矩的變化Fig. 6 The changes of moment around Z-direction when the frequency of heaving motion is 0.4

Origin 軟件具有自定義擬合函數(shù)的功能，利用式（13）編寫出相應(yīng)的擬合函數(shù)，從而得到高精度的擬合結(jié)果。

當(dāng) f = 0.4 時(shí)，得到擬合力和力矩的曲線函數(shù)為：

通過計(jì)算多個(gè)頻率下的純升沉運(yùn)動(dòng)，并進(jìn)行曲線擬合，最終可以得到如表 2 所示的結(jié)果。

根據(jù)表 2 的結(jié)果，通過最小二乘法擬合得到無(wú)因次的加速度系數(shù)，如圖 7 所示。

表 2 純升沉運(yùn)動(dòng)擬合結(jié)果Tab. 2 The fitting results in heaving motion

圖 7 利用最小二乘法計(jì)算加速度系數(shù)Fig. 7 The acceleration coefficients applying the least square method

類似于純升沉運(yùn)動(dòng)，可以得到其他幾組運(yùn)動(dòng)的結(jié)果，結(jié)果如表 3 和表 4 所示。

通過表 3 和表 4 可見，使用 PMM 試驗(yàn)?zāi)M結(jié)果與利用面元法得到的結(jié)果之間差別較小，最大誤差小于11%，滿足實(shí)際工程需求。

表 3 純升沉和純橫蕩運(yùn)動(dòng)結(jié)果Tab. 3 The results in heaving and swaying motions

表 4 純俯仰和純搖艏運(yùn)動(dòng)結(jié)果Tab. 4 The results in pitching and yawing motions

3 結(jié) 語(yǔ)

本文通過面元法編寫程序得到水下滑翔機(jī)附加質(zhì)量，同時(shí)結(jié)合 PMM 試驗(yàn)原理，利用 Fluent 軟件和動(dòng)網(wǎng)格技術(shù)，得到給定運(yùn)動(dòng)下的水動(dòng)力和力矩，從而得到水下滑翔機(jī)的慣性水動(dòng)力系數(shù)，通過對(duì)計(jì)算結(jié)果分析，可以得到：

1）對(duì)于水下滑翔機(jī)，Hess-Smith 面元法程序和Fluent 仿真得到的附加質(zhì)量結(jié)果均比較精確，滿足工程需求。

2）相比于 Fluent 軟件仿真需要花費(fèi)較多時(shí)間多次建模與計(jì)算，Hess-Smith 方法僅需要重新劃分物體的面網(wǎng)格，節(jié)省大量時(shí)間，但編寫程序較為復(fù)雜。

3）本文提供了計(jì)算水下滑翔機(jī)附加質(zhì)量的 2 種方法，介紹了兩者原理與優(yōu)劣，同時(shí)得到精確的附加質(zhì)量，驗(yàn)證了 2 種方法的可行性和準(zhǔn)確性，對(duì)水下滑翔機(jī)的水動(dòng)力設(shè)計(jì)和運(yùn)動(dòng)仿真具有重要指導(dǎo)意義和參考價(jià)值。

[1]ZHANG S W, YU J C. Spiraling motion of underwater gliders: modeling, analysis, and experimental results[J]. Ocean Engineering, 2013, 60: 1–13.

[2]孟凡豪. 50 kg級(jí)水下自航行器整體水動(dòng)力學(xué)性能優(yōu)化設(shè)計(jì)[D]. 杭州: 中國(guó)計(jì)量學(xué)院, 2014: 56–57. MENG Fan-hao. Shape optimization design of 50 kilograms autonomous underwater vehicle based on hydraulic characteristics[D]. Hangzhou: China Jiliang University, 2014: 56–57.

[3]林超友, 朱軍. 潛艇近海底航行附加質(zhì)量數(shù)值計(jì)算[J]. 船舶工程, 2003, 25(1): 26–29. LIN Chao-you, ZHU Jun. Numerical computation of added mass of submarine maneuvering with small clearance to sea-bottom[J]. Ship Engineering, 2003, 25(1): 26–29.

[4]楊勇, 鄒早建, 張晨曦. 深淺水中KVLCC船體橫蕩運(yùn)動(dòng)水動(dòng)力數(shù)值計(jì)算[J]. 水動(dòng)力學(xué)研究與進(jìn)展, 2011, 26(1): 85–93. YANG Yong, ZOU Zao-jian, ZHANG Chen-xi. Calculation of hydrodynamic forces on a KVLCC hull in sway motion in deep and shallow water[J]. Chinese Journal of Hydrodynamics, 2011, 26(1): 85–93.

[5]孫銘澤, 王永生, 張志宏, 等. 基于網(wǎng)格變形技術(shù)的全附體潛艇操縱性計(jì)算[J]. 武漢理工大學(xué)學(xué)報(bào)(交通科學(xué)與工程版), 2013, 37(2): 420–424. SUN Ming-ze, WANG Yong-sheng, ZHANG Zhi-hong, et al. Numerical simulation of submarine maneuverability based on mesh deformation technology[J]. Journal of Wuhan University of Technology (Transportation Science & Engineering), 2013, 37(2): 420–424.

[6]黃勁, 余志毅, 劉朝勛, 等. 基于CFD的兩棲車輛水動(dòng)力導(dǎo)數(shù)的計(jì)算方法[J]. 車輛與動(dòng)力技術(shù), 2014(2): 39–43. HUANG Jin, YU Zhi-yi, LIU Chao-xun, et al. Numerical method for hydrodynamic coefficients of an amphibious vehicle based on CFD[J]. Vehicle & Power Technology, 2014(2): 39–43.

[7]馬騁, 連璉. 水下運(yùn)載器操縱控制及模擬仿真技術(shù)[M]. 北京:國(guó)防工業(yè)出版社, 2009. MA Cheng, LIAN Lian. Maneuvering control and simulation technology of underwater vehicles[M]. Beijing: National Defense Industry Press, 2009.

[8]張曉頻. 多功能潛水器操縱性能與運(yùn)動(dòng)仿真研究[D]. 哈爾濱:哈爾濱工程大學(xué), 2008: 25–26. ZHANG Xiao-pin. Research on maneuverability and motion simulation of multifunction vehicle[D]. Harbin: Harbin Engineering University, 2008: 25–26.

Numerical method for added mass of an underwater glider

YANG Lei, CAO Jun-jun, YAO Bao-heng, ZENG Zheng, LIAN Lian
(Naval Architecture and Ocean Engineering National Laboratory, Shanghai Jiaotong University, Shanghai 200240, China)

Added mass of an underwater glider is quite important for the motions of glider. In this paper, the added mass of an arbitrary three-dimensional body is obtained through Hess-Smith method. Then an underwater glider which was designed by our laboratory is meshed by Gambit software in order to obtain its added mass. Besides, the Planar Motion Mechanism (PMM) tests of the glider are simulated by using CFD software, dynamic mesh technique and UDF. By comparing with the Hess-Smith results, the characters and advantages of Hess-Smith method and PMM are analyzed.

added Mass；underwater glider；Hess-Smith；PMM

U661.33；P715.5

A

1672–7619(2016)12–0116–05

10.3404/j.issn.1672–7619.2016.12.023

2016–05–23；

2016–06–12

國(guó)家自然科學(xué)基金資助項(xiàng)目(51279107；41527901)；上海市科委項(xiàng)目基金資助項(xiàng)目(13dz1204600)

楊磊(1993–)，男，碩士研究生，研究方向?yàn)樗潞叫衅魉畡?dòng)力分析及仿真。