吳利紅，封錫盛，胡志強

（1中國科學院沈陽自動化研究所機器人學國家重點實驗室，沈陽 110016；2中國科學院研究生院，北京 100049）

三維動態(tài)混合網(wǎng)格在AUV發(fā)射過程中的應用

吳利紅1,2，封錫盛1，胡志強1,2

（1中國科學院沈陽自動化研究所機器人學國家重點實驗室，沈陽 110016；2中國科學院研究生院，北京 100049）

為研究AUV從有界流場自航發(fā)射到無界流場的運動邊界的擾流場，文章采用了三維動態(tài)混合網(wǎng)格方法進行數(shù)值模擬的策略。混合網(wǎng)格由三菱柱/四面體/六面體網(wǎng)格構成，當AUV運動時，靠近AUV的三菱柱網(wǎng)格隨AUV運動，外層是靜止的六面體網(wǎng)格,中部的四面體網(wǎng)格隨AUV運動而變形或者重構。數(shù)值仿真結果給出了不同時刻AUV表面的壓力分布、整個航程AUV的阻力系數(shù)變化，其值與理論結果吻合。同時研究了直徑比對發(fā)射管航行的附加質量和阻力系數(shù)的影響，這為水下對接AUV提供了有效手段。

AUV；自航發(fā)射；三維混合網(wǎng)格；非結構動網(wǎng)格；水下對接

1 引 言

求解包含運動邊界的非定常擾流場問題，如自由表面、多體相對運動和流體與結構耦合等問題，是計算流體力學中的一個難點。動態(tài)混合網(wǎng)格方法是近年來出現(xiàn)的解決運動邊界的有效方法，即采用混合網(wǎng)格對復雜區(qū)域進行網(wǎng)格劃分，采用非結構動網(wǎng)格實現(xiàn)運動邊界的網(wǎng)格處理。

應用混合網(wǎng)格對定常擾流下復雜區(qū)域的網(wǎng)格劃分較為成熟，已發(fā)展了多種混合網(wǎng)格：三菱柱/四面體、四面體/三菱柱/金字塔，四面體/自適應直角坐標網(wǎng)格、直角坐標網(wǎng)格/四面體/三菱柱[1]。應用非結構動網(wǎng)格對二維或者小振幅的運動邊界仿真的文獻也日益涌現(xiàn)，主要用于解決空氣動力學中機翼的振蕩擺動[2]、二維機翼—外掛物分離[3]和控制翼偏轉[4]等問題。將并行計算應用到定常流動問題也日益廣泛，代表性的有夏健[5]對DLR_F6帶攻角定常飛行流場采用1～8個分布式內存計算節(jié)點進行計算；Cavallo[6]采用并行網(wǎng)格自適應策略仿真多體分離的相對運動擾流場。

而將混合網(wǎng)格用于非定常運動中的網(wǎng)格劃分、將非結構動網(wǎng)格應用到三維大位移運動和將并行計算用于運動邊界的非定常計算中的研究較少。

本文采用動態(tài)混合網(wǎng)格結合并行計算實現(xiàn)了AUV從有界流場航行到無界流場的三維大位移擾流場的數(shù)值模擬，這為當前三維大位移的運動邊界移動較難采用數(shù)值仿真實現(xiàn)的問題提供了方法上的借鑒，同時彌補了有關AUV從有界流場中自航發(fā)射的擾流場只有理論研究而缺少數(shù)值研究的不足[7-8]。

2 動態(tài)混合網(wǎng)格方法

2.1 混合網(wǎng)格生成

假設有界流場為一端封閉的管道，簡稱發(fā)射管，AUV自航發(fā)射在初始時刻的幾何模型如圖1所示。AUV和發(fā)射管直徑比為0.625，AUV直徑333mm，長3.8m；管長為8m。初始時刻，AUV尾部離發(fā)射管端部A點0.4m，首部離發(fā)射管口B點3.8m。數(shù)值仿真AUV自航發(fā)射的三維大位移繞流場包括AUV從初始位置開始，以3Kn速度從發(fā)射管中航行到發(fā)射管外的無界流場中的繞流，AUV總航程約15m。載體坐標系Oxyz建立在AUV重心G處，O點與G點重合。

整個流場的網(wǎng)格采用混合網(wǎng)格構成，包含三菱柱、四面體、六面體三種網(wǎng)格。網(wǎng)格拓撲結構如圖2所示，包含4個網(wǎng)格區(qū)域，其中區(qū)域1為繞AUV表面的三菱柱區(qū)域；區(qū)域2為四面體網(wǎng)格；區(qū)域3是靠近發(fā)射管壁面的六面體網(wǎng)格區(qū)域；區(qū)域4為填充無界流場的四面體網(wǎng)格。

當AUV運動時，區(qū)域1隨AUV運動；區(qū)域2變形；區(qū)域3和區(qū)域4靜止不動。圖3給出了初始時刻對稱面上靠近AUV首部、AUV和發(fā)射管間隙處和靠近AUV尾部的網(wǎng)格分布，分別對應著圖3中（a）、（b）、（c）圖。

圖3 初始時刻對稱面上的網(wǎng)格分布Fig.3 Mesh in symmetry at initial time

2.2 動網(wǎng)格更新方法

動網(wǎng)格更新有三種方法：彈簧近似網(wǎng)格光順法、動態(tài)層方法和局部網(wǎng)格重構方法。彈簧近似網(wǎng)格光順法適用于小位移問題；動態(tài)層方法適應于單向運動；而局部網(wǎng)格重構方法適用于局部網(wǎng)格重構。第一種方法不改變網(wǎng)格的拓撲結構，這樣能保證網(wǎng)格質量；后兩種方法改變了網(wǎng)格的拓撲結構，新網(wǎng)格的變量需要從舊網(wǎng)格中插值得到。對于大位移和強切變問題僅靠節(jié)點松弛不能保證網(wǎng)格質量，而且有可能出現(xiàn)網(wǎng)格相交的情況。為了克服這一困難，本文采用彈簧近似網(wǎng)格光順法和局部網(wǎng)格重構法相結合，解決大位移和強切變非定常流動的網(wǎng)格移動問題。

彈簧近似模型將四面體非結構網(wǎng)格看作是一個彈簧網(wǎng)格系統(tǒng)，每條邊都認為是一根具有一定倔強系數(shù)的彈簧，初始狀態(tài)下，每個節(jié)點都處于平衡狀態(tài)，當邊界運動時，節(jié)點i所受到的合力如公式（1）所示，節(jié)點i的網(wǎng)格更新后的位置如公式（2）所示，其中，分別表示節(jié)點i和臨近節(jié)點j的位移，ni是節(jié)點i的臨近節(jié)點數(shù)目，kij是連接節(jié)點i和鄰近節(jié)點j的彈簧剛度系數(shù)。公式（2）中的， 分別表示n，n+1迭代時間步下的節(jié)點i的位置； 表示臨近節(jié)點j在內循環(huán)m步后收斂后的位移，其中上標m，c分別表示內循環(huán)迭代步數(shù)，內循環(huán)收斂的狀態(tài)。

當網(wǎng)格移動后，網(wǎng)格的品質會發(fā)生改變，表征網(wǎng)格的品質的參數(shù)如光順性、正交性、節(jié)點分布特性，高寬比（aspect ratio）和單元尺寸會發(fā)生改變。局部網(wǎng)格更新通過設置網(wǎng)格控制參數(shù)如網(wǎng)格的最小、最大長度，網(wǎng)格單元的偏斜度、尺度函數(shù)參數(shù)可以控制網(wǎng)格移動后的品質，當網(wǎng)格不滿足以上控制參數(shù)設置的值域時，網(wǎng)格被標示出，并在局部挖出一個小洞，取出洞的邊界信息，然后利用陣面推進法重新生成網(wǎng)格并進行優(yōu)化。為了使得重新生成的局部網(wǎng)格與全局網(wǎng)格保持一定的光順性，對局部網(wǎng)格采用尺度函數(shù)進行控制。

3 并行計算

并行計算就是利用多個計算節(jié)點同時進行計算。合理的網(wǎng)格分區(qū)方法是有效并行計算的前提。網(wǎng)格分區(qū)時，需要選擇生成網(wǎng)格的分割方法、設置分割數(shù)、選擇區(qū)域和記錄、以及使用的優(yōu)化方法。其目標是生成等數(shù)量單元的網(wǎng)格分塊、分割的接觸面數(shù)最小、分割的領域數(shù)最小。本文通過對不同網(wǎng)格分區(qū)方法進行試驗對比，確定出采用和初始網(wǎng)格拓撲結構一致的Cylindrical R axes的網(wǎng)格分區(qū)方法是合適的。

并行機的性能一般用加速比來評測。計算加速比的一般方法是：假設在一個處理器上運行某程序需時間T1，在P個節(jié)點（各節(jié)點的CPU性能一樣）上運行此程序需時間Tp，則加速比Sr=T1/Tp，并行效率η=Sr/P。

本文將網(wǎng)格節(jié)點數(shù)分別為449,442和592,669的計算網(wǎng)格在1～12個計算節(jié)點上并行計算，獲得其加速比如圖4。從圖4可以看出，加速比曲線與定常計算的正比曲線不一樣[5]，采用8個計算節(jié)點具有最高加速比。

圖4 并行計算加速比Fig.4 Speedup ratio for parallel computing

4 數(shù)值仿真

4.1 給定直徑比的數(shù)值仿真結果

利用上述動態(tài)混合網(wǎng)格和動態(tài)并行計算方法，數(shù)值模擬了AUV從發(fā)射管中航行到無界流場的三維大位移運動邊界的繞流，圖5為t=2s時刻對稱面上靠近首尾部的移動網(wǎng)格分布。

圖6給出了不同時刻對稱面上的壓力分布，即t=0.5s時AUV在發(fā)射管內和t=3.5s時AUV通過發(fā)射管口的壓力分布。圖7是AUV從發(fā)射管航行到無界流場中全程約15m的阻力系數(shù)圖，包括全程壓差阻力系數(shù)Cp、摩擦阻力系數(shù)Cf和總阻力系數(shù)Ct，阻力系數(shù)與力的換算關系見公式（3），其中，F(xiàn)包含摩擦阻力、壓差阻力和總阻力，C表示這三種系數(shù)，ρ為海水密度，V為AUV航速，A為AUV橫截面面積。分析總阻力系數(shù)Ct，發(fā)現(xiàn)其在發(fā)射管內的阻力系數(shù)與Richard[7]采用控制體方法計算相應載體在無限長發(fā)射管中航行的阻力系數(shù)值2.58接近，Ct在無界流場中的阻力系數(shù)段與采用CFX軟件計算定常來流獲得的系數(shù)進行對比，其中摩擦阻力與CFX軟件計算的定常值0.35接近。

圖7 AUV自航發(fā)射全程Cp，Cf和Ct Fig.7 Cp,Cf and Ct for AUV swim-out

從整個航程中還可以看出摩擦阻力系數(shù)的變化較小，而壓差阻力變化較大，尤其是在發(fā)射管內，壓差阻力占總阻力的71%，這主要是因為發(fā)射管內的有界邊界的影響，導致兩者的間隙流流速大大增加，壓力大大下降。因此非常有必要研究在發(fā)射管中運動直徑比對AUV航行阻力的影響。

4.2 改變直徑比對AUV附加質量和航行阻力的影響

發(fā)射管中運動的直徑比為β，如公式（4）所示，其中D為發(fā)射管直徑，d為AUV直徑。隨著直徑比的降低，間隙比γ（如公式（5）所示）越小，圖8給出直徑比與間隙比的關系。當AUV在不同直徑發(fā)射管中航行時，其附加質量也會變化，附加質量mh與直徑比的關系如公式（6）所示[7]，附加質量系數(shù)λ為附加質量與AUV所排開水的質量的比值，它與直徑比的關系如公式（7）所示，附加質量隨AUV航程S的關系見公式（8）所示。其中Λ為AUV排水體積；L為發(fā)射管長度，l為AUV長度,而且L＞l，l0為初始時刻AUV距離發(fā)射管底端長度。λ0為AUV完全在發(fā)射管內的附加質量系數(shù)。圖9給出附加質量系數(shù)隨航程和直徑比的關系曲線，可以看出，隨著直徑比的減少，附加質量迅速增加；隨著航程的增加，AUV漸漸靠近發(fā)射管口，其附加質量也減少。

直徑比不僅對附加質量有影響，對管道中航行的阻力也有較大的影響，其關系如圖10所示，直徑比減少，阻力迅速增加，當直徑比小于1.223時，阻力曲線陡升，影響阻力的主要原因是間隙比減小，壓差阻力大大增加，當直徑比為1.07，間隙比為0.125時，壓差阻力占總阻力的92.7%，如圖11所示。從圖10和圖11可以看出相同直徑比前提下，AUV以不同雷諾數(shù)Re航行，阻力系數(shù)變化較小。

5 結 論

本文采用三菱柱/四面體/六面體構成的動態(tài)混合網(wǎng)格數(shù)值仿真了AUV從發(fā)射管中自航發(fā)射的三維大位移運動過程，得到了如下結論：

（1）AUV從發(fā)射管航行到無界流場時，在發(fā)射管中的最大阻力約為無界流場中的7.4倍，阻力隨著航程的增加而降低，當AUV尾部通過發(fā)射管發(fā)射管口時，阻力值降低到最??；

（2）發(fā)射管中運動時，當直徑比小于2時，壓差阻力占總阻力的主要成份，直徑比為1.07，壓差阻力約為總阻力的92.7%；

（3）發(fā)射管中運動時，附加質量和阻力系數(shù)都隨直徑比減小而迅速增大，而相同直徑比下改變雷諾數(shù)對阻力系數(shù)影響較小，同時對于有限長發(fā)射管，附加質量還隨航程增加而減??；

（4）移動網(wǎng)格數(shù)值仿真非常耗時，可以采用多機并行來加速求解，本文采用8個節(jié)點的多節(jié)點共享內存并行，可以獲得最大加速比因子4.5。

本文數(shù)值仿真為AUV水下對接提供了數(shù)值仿真方法，同時數(shù)值仿真結果也為AUV水下對接工程提供了理論參考。

[1]Zhang L P,Yang Y J,Zhang H X.Numerical simulations of 3D inviscid/viscous flow fields of Cartesian/unstructured/prismatic hybrid grids[C]//Proceedings of the 4th Asian CFD Conference.Mianyang,China,2000.

[2]郭 正,劉 君,瞿章華.非結構動網(wǎng)格在三維可動邊界問題中的應用[J].力學學報,2003,35(2):140-146.

[3]張來平,王振亞,楊永健.復雜外形的動態(tài)混合網(wǎng)格生成方法[J].空氣動力學學報,2004,22(2):231-236.

[4]Mitsuhiro M Y,Kazuhiro N H.Unstructured dynamic mesh for large movement and deformation[C].AIAA,2002.

[5]夏 健,伍貽兆.基于混合網(wǎng)格的三維Navier-Stokes方程并行計算方法[J].航空學報,2005,26(3):290-293.

[6]Cavallo P A,Sinha N.Parallel unstructured mesh adaptation for transient moving body and aeropropulsive applications[C]//42nd AIAA Aerospace Sciences Meeting and Exhibit.Reno,Nevada,AIAA,2004:1057-1067.

[7]Richard F H.The near-field flow and drag on cylindrical bodies moving concentrically inside very long tubes[R].NUSC,Newport Laboratory,1991:27-33.

[8]王燕飛,張振山,張 萌.自航發(fā)射魚雷內彈道模型與仿真研究[J].系統(tǒng)仿真學報,2006,18(2):316-318,326.

Application of 3-D hybrid dynamic grids to simulate the flow in AUV swim-out

WU Li-hong1,2,FENG Xi-sheng1,HU Zhi-qiang1,2
(1 State key laboratory of Robotics,Shenyang Institute of Automation,Chinese Academy of Sciences,Shenyang 110016,China;2 Graduate School of the Chinese Academy of Sciences,Beijing 100049,China)

A 3-D hybrid dynamic grids method is presented to simulate moving boundary in Autonomous Underwater Vehicle(AUV)swim-out from bounded domain to unbounded domain.The 3-D hybrid grids are composed of prismatic/tetrahedral/hexahedral grids from inner to outer layer respecitvely.With AUV moving the prismatic grids move with AUV,while the hexahedral grids remain stationary.Meanwhile,the tetrahedral grids are deformed or re-meshed according to the motion of AUV.Results from the numerical simulation show that the pressure contour in AUV surface at different times,and the resistant coefficient of AUV at the whole sailing which agrees well with the theory are obtained.The relation between diameter ratio with AUV reistant coefficient and addmass coefficient in tube is also shown,which gives guide for AUV underwater docking to tube.

AUV;swim-out;3-D hybrid grids;dynamic unstructured grids;underwater docking

TP24

A

1007-7294（2010）07-0717-06

2009-04-29

機器人學國家重點實驗室自主課題（RLZ200810）

吳利紅（1978-），女，博士生，研究領域：水下機器人水動力分析，多體相對運動數(shù)值仿真。