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

吳利紅，封錫盛，胡志強(qiáng)

（1中國科學(xué)院沈陽自動(dòng)化研究所機(jī)器人學(xué)國家重點(diǎn)實(shí)驗(yàn)室，沈陽 110016；2中國科學(xué)院研究生院，北京 100049）

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

吳利紅1,2，封錫盛1，胡志強(qiáng)1,2

（1中國科學(xué)院沈陽自動(dòng)化研究所機(jī)器人學(xué)國家重點(diǎn)實(shí)驗(yàn)室，沈陽 110016；2中國科學(xué)院研究生院，北京 100049）

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

AUV；自航發(fā)射；三維混合網(wǎng)格；非結(jié)構(gòu)動(dòng)網(wǎng)格；水下對(duì)接

1 引 言

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

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

而將混合網(wǎng)格用于非定常運(yùn)動(dòng)中的網(wǎng)格劃分、將非結(jié)構(gòu)動(dòng)網(wǎng)格應(yīng)用到三維大位移運(yùn)動(dòng)和將并行計(jì)算用于運(yùn)動(dòng)邊界的非定常計(jì)算中的研究較少。

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

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

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

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

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

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

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

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

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

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

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

3 并行計(jì)算

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

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

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

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

4 數(shù)值仿真

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

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

[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]郭 正,劉 君,瞿章華.非結(jié)構(gòu)動(dòng)網(wǎng)格在三維可動(dòng)邊界問題中的應(yīng)用[J].力學(xué)學(xué)報(bào),2003,35(2):140-146.

[3]張來平,王振亞,楊永健.復(fù)雜外形的動(dòng)態(tài)混合網(wǎng)格生成方法[J].空氣動(dòng)力學(xué)學(xué)報(bào),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ì)算方法[J].航空學(xué)報(bào),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ā)射魚雷內(nèi)彈道模型與仿真研究[J].系統(tǒng)仿真學(xué)報(bào),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

機(jī)器人學(xué)國家重點(diǎn)實(shí)驗(yàn)室自主課題（RLZ200810）

吳利紅（1978-），女，博士生，研究領(lǐng)域：水下機(jī)器人水動(dòng)力分析，多體相對(duì)運(yùn)動(dòng)數(shù)值仿真。