拋棄式探頭下沉運動數值計算的網格劃分方法研究

陳振濤，劉 鳳，王曉蕾，葉 松

（國防科技大學 氣象海洋學院，江蘇 南京 211101）

拋棄式探頭下沉運動數值計算的網格劃分方法研究

陳振濤，劉 鳳*，王曉蕾，葉 松

（國防科技大學 氣象海洋學院，江蘇 南京 211101）

各類拋棄式探頭下沉運動的數值計算研究都需要對探頭表面和整個計算域進行網格劃分的前處理過程，選擇不同的網格劃分方法和劃分精度，會影響數值計算速度和計算精度，最終生成的網格質量決定了計算結果的收斂性和準確性。針對拋棄式探頭的復雜結構，對比了兩種網格劃分方法、三種網格劃分精度，完成探頭計算區域的網格劃分，采用k-ε湍流模型，進行拋棄式探頭下沉運動的數值計算，并將計算結果與水箱和水庫實驗結果進行對比，驗證了混合網格劃分方法和普通精度網格對拋棄式探頭下沉運動數值計算的適用性，研究結果對類似較為復雜結構的水下運動體的數值計算前處理過程具有一定的參考和借鑒意義。

拋棄式探頭；數值計算；網格劃分；混合網格

目前針對各類拋棄式探頭下沉運動的數值計算，主要采用 CFD（Computational Fluid Dynamics，計算流體動力學）技術[1-3]，其處理過程主要包括前處理、有限元計算以及后處理三部分，其中網格劃分是前處理過程的主要工作，較為繁瑣和耗時，但影響到有限元計算成功的關鍵[4]。

不同的流體問題需要采用不同的數值求解方法，相應的網格形式也有所不同，網格生成的類型和數量會影響計算速度和計算精度，網格質量的好壞則直接影響到計算結果的收斂性和準確性[5]。

針對拋棄式探頭（本文特指NMOHEMS[6]探頭，以下同）的復雜結構，對比不同網格劃分方案，選擇最佳的網格劃分方法和劃分精度，完成探頭表面和計算區域的網格劃分，以開展探頭下沉運動數值計算。

1 網格劃分方法

1.1 網格劃分技術

利用GAMBIT前處理軟件[7]，其可劃分的網格類型主要包括結構網格和非結構網格。其中結構網格[8]結構清晰，網格節點排列有序，相鄰節點之間的關系明確，所有節點周圍網格的拓撲結構都相同。2D中表現為四邊形網格，3D中為六面體網格。結構網格處理簡單的幾何體網格時十分有效，但是不適用于復雜外形的幾何體網格生成。

非結構網格[9]消除了網格節點的結構性約束，可以任意分布其網格節點和單元，方便靈活地處理一些復雜的邊界情況。但需要控制好網格的大小和節點密度，以保證網格的質量。非結構網格的主要缺點是計算效率較低，特別是對計算機內存的要求較高，相應計算時間較長，所以其普適性受到制約。

針對復雜外形的網格問題，近年來發展起來的還有一種自適應笛卡爾網格[10]，根據幾何體表面和流場特點，在邊界等局部區域不斷細化加密網格，以實現較高精度的一種非均勻網格。自適應笛卡爾網格的生成過程省時方便且數據結構較為簡單，但是外形描述的網格精度較低，無法做到貼體網格，對于一些邊界層流動問題，無法滿足計算要求。三種網格方法生成的2D網格對比，如圖1所示。

圖1 三種網格方法生成的2D網格對比

在解決拋棄式探頭類復雜外形網格的生成問題時，可以采用將結構網格和非結構網格相結合的混合網格方法[11]，即在幾何體的復雜表面區域采用非結構網格，而在遠離幾何體的區域采用結構網格，這樣生成的網格質量較高，同時保證較高的計算效率。

1.2 計算域及邊界條件

探頭的整個計算域由進流、出流、壁面和控制域邊界組成，圖2為Z=0時計算區域模型的截面。圖中，L為探頭的長度，R為探頭最大半徑，X為水平方向，Y為重力加速度的反方向。

圖2 探頭模型計算域示意圖

各邊界條件如下：采用速度進口邊界條件，入流斷面為均勻來流；采用壓力出流邊界條件，認為湍流已經達到平衡，以加快收斂速度；探頭模型表面和計算域的其余外邊界均設為無滑移壁面；模型近壁采用標準壁面函數進行處理。

其中，速度進口的湍流參數項選擇為湍流強度I和特征長度l（Intensity and Length Scale）方法。充分發展的管內流動中，湍流強度I可以利用經驗公式[12]給出，湍流強度是一個很小的量，通常在1%～10%范圍內，對于邊界條件為外部入流，其計算公式為：

式中：u'為速度脈動；uˉ為平均速度；，ReDH為入口邊界處的雷諾數。根據上式，當雷諾數為104時，湍流強度I約為5.06%。

特征長度l是與湍流中大渦尺度相關的物理量，圓管內流動時與管長相關，其計算公式[12]為：

式中：L為圓管長度，本文中探頭運動為外部繞流問題；L即為探頭沿水流方向的長度。

1.3 計算域的網格劃分方法

因為探頭外形較為復雜，無法采用結構網格的劃分方法，采用完全非結構網格和混合網格兩種方法進行比較。

完全非結構網格的劃分方法，在探頭表面劃分邊界層，其余區域均劃分為四面體形的非結構網格，探頭附近區域網格較密集，遠離探頭區域網格較稀疏，模型總網格數為2 546 718個。

采用混合網格劃分計算域，在探頭表面附近劃分邊界層，劃分為三棱柱形的半結構網格；邊界層外的探頭附近區域，網格類型采用四面體形的非結構網格；遠離探頭區域則劃分為六面體形的結構網格。進一步加密導流腔和尾翼附近網格，模型總網格數為1 964 509個。

兩種方法的探頭表面網格均采用三角形非結構網格進行劃分，如圖3所示；計算域的縱向對稱面網格分布如圖4所示。

圖3 探頭表面的網格劃分

圖4 兩種網格的縱向對稱面網格分布

選擇Realizable k-ε湍流模型，分別采用上述兩種網格對流場進行計算，網格數、迭代次數n和計算的阻力系數Cd如表1所示。其中實驗結果采用實驗室玻璃水箱，高精度調整探頭質量，利用高速攝像和圖像處理技術，經過詳細的誤差分析和修正后，得到多個較小質量探頭的下沉運動的下沉曲線和極限速度[13]。

從表中可以看到，混合網格的網格數比完全非結構網格減少22.86%，迭代次數也相應地減少，并且隨著速度的增加，迭代次數減少得更加明顯，而兩者的計算結果非常接近。證明了混合網格的計算速度和計算效率比完全非結構網格有相當大的提高，同時保證了計算結果的穩定性。

1.4 網格劃分精度

網格劃分精度的不同對占用內存、計算時間和計算結果有較大的影響。采用混合網格的劃分方法，設定來流速度分別為0.1 m/s，0.2 m/s，0.3 m/s，選擇以下3種網格精度進行比較：（a）較稀疏（Coarse mesh）、（b）普通（Nominal mesh）、（c）較精密（Fine mesh）。

三種網格精度的最小、最大網格尺寸、網格數、迭代次數和計算的阻力系數比較如表2所示。

表2 三種網格精度及阻力系數的比較

表中可以看到，網格較稀疏時的網格數最小，迭代次數最少，但是計算的阻力系數結果誤差偏大；網格較精密時的網格數最大，迭代次數最多，且成指數倍增加，計算耗時較長，且容易產生內存溢出；普通網格精度的計算耗時較短，與較精密網格的計算結果相差不大，所以最終選定網格精度為普通精度。

1.5 求解控制參數的設置

采用有限體積法離散控制方程和湍流模型，壓力方程的離散格式包括：Standard、PRESTO、Second Order和Body Force四種，其中PRESTO主要應用于高速流動，特別是含有旋轉及高曲率的情況下；Second Order用于可壓流動；Body Force用于含有大體力的流動；所以采用Standard離散格式進行離散壓力方程。

動量方程、湍流方程和雷諾應力方程均采用二階迎風格式進行離散，壓力速度耦合迭代采用SIMPLEC算法，欠松弛因子取為1.0。

此外，求解過程中需要監視各變量的殘差，殘差收斂的精度越高，湍流發展越充分，計算的結果也越精確。殘差收斂精度的選取對計算時間和計算精度的影響也較為明顯，所以在選擇時應當盡量滿足湍流充分發展。選擇以下3種殘差收斂精度進行比較：（a）10-3，（b）10-4，（c）10-5。

分別將來流速度設為 0.1 m/s，0.5 m/s，1.0 m/s，1.5 m/s和2.0 m/s進行計算，迭代次數和計算得到的阻力大小如表3所示。

表3 三種殘差收斂精度時的迭代次數和阻力系數比較

表中可以看到，殘差收斂精度為10-3時的迭代次數較少，但是阻力系數的計算結果偏大；殘差收斂精度為10-4時的迭代次數都非常大，其中“＞10 000”是指迭代10 000次時殘差仍然沒有達到指定的收斂精度，湍流發展較為充分，所以阻力系數的計算結果較為精確；殘差收斂精度為10-4時的迭代次數較少，而阻力系數的計算結果與10-5時較為接近，所以最終選定殘差收斂精度為10-4時。

2 計算結果及分析

采用混合網格劃分方法、普通網格精度，對探頭表面和計算域進行網格劃分，選擇Realizable k-ε湍流模型，采用Standard離散格式離散壓力方程，采用二階迎風格式離散動量方程、湍流方程和雷諾應力方程，SIMPLEC算法進行壓力速度耦合迭代處理，殘差收斂精度為10-4，對探頭及其計算域進行數值計算。

圖5給出了水箱實驗得到的不同質量探頭，極限速度對應的雷諾數與阻力系數的關系；*為實驗結果，□為數值計算結果，雷諾數范圍在104左右。

圖5 雷諾數與阻力系數的關系

可以看到兩個結果的變化趨勢較為一致，阻力系數均隨著雷諾數的增大而減小，模擬結果整體均小于實驗結果，并且隨著雷諾數的增加，兩者的差別逐漸增大。分析可能有兩部分原因，一是建立的數值計算模型較為理想和光滑，造成了計算得到的阻力偏小；二是沒有考慮探頭下沉時水箱中存在的波動干擾等因素的影響。

表4為六種不同質量探頭對應阻力系數的數值計算結果與水箱實驗結果對比關系。

表4 數值計算結果與水箱實驗結果對比

可以看到，不同質量時兩個結果都非常接近，誤差δ控制在2%以內，驗證了雷諾數范圍在104左右時，網格劃分方法和數值計算模型的正確性和可行性。

此外，還進行了水庫實驗，對較大質量探頭的下沉運動進行測量，得到0.48 kg探頭的極限速度為1.92 m/s，雷諾數在105范圍，實驗得到的測量結果經過水面風速和裝置摩擦等誤差處理后，得到測量修正值，三者的對比情況見表5。

表5 水庫測量結果、修正值與數值計算結果對比

表中可見，經過誤差處理后的測量修正值與數值計算結果較為接近，阻力系數相對誤差為1.55%，極限速度的相對誤差只有-0.31%，最大誤差均控制在2%以內，驗證了雷諾數范圍在105左右時，網格劃分和數值計算模型的正確性和可行性。

圖6～圖7分別給出了v為0.001 m/s，0.01 m/s，0.1 m/s和2.0 m/s時探頭縱向對稱面的壓力分布和y+分布情況。

圖6 不同速度時縱向對稱面的壓力分布曲線

從圖中可以看到，不同速度時縱向對稱面上的高低壓分布情況基本一致，只是壓力的強度和范圍有所區別。流場在靠近探頭頭部（-200 mm處）時壓力值迅速增大，之后跟隨一個低壓區，壓力又迅速減小；尾翼前端（-20 mm處）因為尾翼的阻礙作用，也形成一個明顯的高壓區，并且跟隨一個小的低壓區；通過尾翼后在尾翼后端（0 mm處）形成一個較小的高壓區。v=0.001 m/s時探頭頭部的壓力強度較小，跟隨的低壓區也不明顯，但是尾翼前端的高壓區和低壓區較為明顯，強度甚至超過探頭頭部；v=0.01 m/s時探頭頭部的壓力強度增大，跟隨低壓區的強度和范圍都有所增大，尾翼后端的高壓區也稍有增強；v=0.1 m/s時探頭頭部的高壓區和跟隨低壓區的壓力強度和范圍繼續增大，尾翼前后的高壓區和跟隨低壓區的增長速度相對減弱；v=2.0 m/s時探頭頭部和跟隨低壓區的壓力強度和范圍繼續增大，尾翼前后的高低壓區相對減弱得較為明顯。

圖7 不同速度時縱向對稱面上的y+分布曲線

從不同速度時縱向對稱面上的y+分布可以看到，橫坐標在-200 mm處，即探頭頭部位置的y+值最大，之后由于導流腔和收縮段的影響y+值迅速減小，到-20 mm處，即尾翼前部位置時y+值又有所增加；v=0.001 m/s時探頭壁面的y+最大值只有0.05左右，v=0.01 m/s時探頭壁面的y+最大值約為0.2，v=0.1 m/s時探頭壁面的y+最大值約為1.2，v=2.0 m/s時探頭壁面的y+值達到8；并且隨著速度的增大，探頭尾部的y+值的增幅相對較小。0～2.0 m/s范圍內探頭表面的y+值均小于8，即探頭表面一直處于粘性底層，證明了邊界層網格劃分和近壁處理能夠真實反應探頭表面邊界層的流動情況。

不同速度時探頭的周圍流場分布情況有所差異，圖8～圖9分別給出探頭在不同速度時周圍的壓力和速度分布情況。可以看到，混合網格能夠較為完整地顯示不同速度時，探頭周圍流場的壓力、速度分布情況。

圖8 不同速度時探頭周圍流場的壓力分布圖

圖9 不同速度時探頭周圍流場的速度分布圖

均勻來流均受到了探頭的阻礙，在靠近探頭頭部附近形成了一個局部靜壓高于來流靜壓的高壓區，探頭近表面附近的流場速度小于外圍流場速度；速度為0.001 m/s時，壓力和速度較小，尚未出現分離現象，也未產生湍流渦旋，所以此時探頭周圍仍然是層流狀態；速度為0.01 m/s時，各流場分布有所加強，壓力和速度開始出現分離現象，探頭尾部出現較為明顯的尾渦區，但是強度和范圍較小，說明此時探頭周圍流場已由層流進入轉捩區；速度為0.1 m/s時，壓力和速度的流場分離現象更加明顯，尾渦區的強度和范圍明顯增強，說明此時流場的湍流得到較強的發展；速度為2.0 m/s時，流場分布進一步增強，湍流已得到充分發展，速度場經過上表面的阻礙后反射，與原流場進行疊加，使得整個計算域的流場分布更加復雜。

3 結束語

本文根據拋棄式探頭特點，對網格劃分方法和劃分精度進行了討論，計算其不同下沉速度時對應的阻力系數；最后對計算結果進行簡單分析，并分別與水箱和水庫實驗結果進行對比，得到主要結論如下：

（1）在解決拋棄式探頭類的復雜外形網格生成問題時，可以采用結構網格和非結構網格相結合的混合網格方法，這樣生成的網格質量較高、網格數較少、運算速度較快且計算結果穩定。

（2）網格劃分精度對占用內存、計算時間和計算結果有較大的影響，為保證計算結果更高的準確性，應根據計算機內存和運算性能選擇較高的網格精度。拋棄式探頭選用普通網格精度的網格數適中，計算耗時較短，計算的阻力系數結果與實驗結果較為接近。

（3）CFD仿真計算結果除了與網格劃分相關，還與選擇的湍流模型、離散求解方案，壓力速度耦合迭代算法，以及殘差收斂精度等有關。

有限元網格生成技術還有曲面網格生成、并行網格生成等很多方法，針對不同問題和需求，選擇的網格劃分方法也不盡相同，需要在實際應用中不斷總結經驗才能快速地得到最佳的網格劃分模型。

[1]肖鴻,劉長根,陶建華.拋棄式溫鹽探頭阻力系數的數值模擬及其實驗驗證[J].海洋技術,2006,25(3):35-37.

[2]孫濤,黃銀水,陶建華.拋棄式溫鹽探頭運動狀態的數值模擬及其實驗驗證[J].海洋通報,2002,21(2):69-76.

[3]劉寧,張瑞,陳文義,等.XCP投棄式海洋探頭阻力系數數值計算及試驗研究[J].海洋技術,2010,29(4):12-14.

[4]杜平安.有限元網格劃分的基本原則[J].機械設計與制造,2000(1):34-36.

[5]邢靜忠,李軍.ANSYS的建模方法和網格劃分[J].中國水運(學術版),2006,9(6):116-117.

[6]葉松,王曉蕾,焦冰，等.NMOHEMS的概念與設計[J].海洋技術,2010,29(1):28-31.

[7]Bhutta M M A,Hayat N,Bashir M H,et al.CFD applications in various heat exchangers design:A review[J].Applied Thermal Engineering,2012,32:1-12.

[8]常煜,張志榮,趙峰.多塊結構化網格在含附體水面船模粘性流場數值計算中的應用[J].船舶力學,2004,8(1):19-25.

[9]何曉聰,曹紅松,趙捍東，等.彈箭外部繞流數值仿真中網格的選擇[J].彈箭與制導學報,2009,29(2):191-194.

[10]韓玉琪,張常賢.基于自適應笛卡爾網格的可壓縮黏性流動數值模擬[J].科學技術與工程,2015(24):115-120.

[11]吳明.艦船不同浮態下阻力的CFD方法研究[D].大連：海軍大連艦艇學院,2009.

[12]王福軍.計算流體動力學分析-CFD軟件原理與應用[M].北京:清華大學出版社,2004.

[13]陳振濤,鐘中,葉松，等.NMOHEMS探頭模型的制作及實驗方法研究[J].實驗流體力學,2013(2):100-105.

Research on the Mesh Dividing Method for Expendable Probe Sinking Simulation

CHEN Zhen-tao,LIU Feng,WANG Xiao-lei,YE Song
Institute of Meteorology and Oceanography,National University of Defense Technology,Nanjing 211101,Jiangsu Province,China

The pre-treatment process with mesh dividing is required for the surface of probe and the whole computational field when conducting expendable probe sinking numerical simulation.Different meshing methods and meshing densities will affect the calculation speed and accuracy,and the quality of the resulting mesh determines the convergence and accuracy of the calculation results.In view of the complex structure of the expendable probe,two meshing methods and three mesh densities are compared,with the mesh division of the probe computational field completed.Numerical simulation of the probe sinking motion is carried out by selected turbulence model,whose results are compared with the experimental results of the laboratory pool and lake.This paper verifies the applicability of the hybrid mesh dividing method and the general precision of mesh density in expendable probe sinking numerical simulation.The results of the study can be used as a reference for the numerical simulation of underwater moving bodies with complex structures.

expendable probe;numerical simulation;mesh dividing;hybrid mesh

P716+.1

A

1003-2029（2017）05-0032-06

10.3969/j.issn.1003-2029.2017.05.006

2017-05-31

國家自然科學基金資助項目（41406107）

陳振濤（1983-），男，博士，講師，研究方向為海洋探測技術。E-mail:czt1212@126.com

劉鳳（1983-），女，碩士，講師，研究方向為海洋探測裝備。E-mail:li_ufeng@126.com