三維動態混合網格在AUV發射過程中的應用

吳利紅，封錫盛，胡志強

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

三維動態混合網格在AUV發射過程中的應用

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

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

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

AUV；自航發射；三維混合網格；非結構動網格；水下對接

1 引 言

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

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

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

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

2 動態混合網格方法

2.1 混合網格生成

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

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

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

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

2.2 動網格更新方法

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

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

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

3 并行計算

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

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

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

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

4 數值仿真

4.1 給定直徑比的數值仿真結果

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

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

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

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

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

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

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

5 結 論

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

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

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

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

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

本文數值仿真為AUV水下對接提供了數值仿真方法，同時數值仿真結果也為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]郭 正,劉 君,瞿章華.非結構動網格在三維可動邊界問題中的應用[J].力學學報,2003,35(2):140-146.

[3]張來平,王振亞,楊永健.復雜外形的動態混合網格生成方法[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]夏 健,伍貽兆.基于混合網格的三維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]王燕飛,張振山,張 萌.自航發射魚雷內彈道模型與仿真研究[J].系統仿真學報,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-），女，博士生，研究領域：水下機器人水動力分析，多體相對運動數值仿真。