平板噴氣粘性流場數值計算方法研究

吳　浩，歐勇鵬（海軍工程大學 艦船工程系，湖北　武漢 430033）

平板噴氣粘性流場數值計算方法研究

吳浩，歐勇鵬
（海軍工程大學 艦船工程系，湖北武漢 430033）

針對平板底部直接噴氣形成氣液混合流的復雜情況，采用 Mixture 模型與 RANS 方程相結合的方法建立了氣液混合流粘性流場數值計算模型。通過對 4 種湍流模型、4 種網格、3 種壁面處理方法進行組合，形成了 8種不同的數值計算方法，分析了壁面函數、壁面第 1 層網格、湍流模型等對數值計算結果的影響，并與試驗結果進行對比，獲得了可有效模擬氣液混合流的數值模型。研究結果表明：采用 RNG k-ε 湍流模型、標準壁面函數、第 1層網格 1 mm、y + 為 31～35 的計算方案，所得結果可用于平底船底部氣液混合流分析。

氣液混合流；Mixture 模型；壁面函數；湍流模型

0　引　言

氣層減阻技術可大幅降低船舶阻力，是實現我國航運業節能減排新目標的重要途徑。相對于模型試驗，目前氣泡船的數值研究相對較少，以下學者開展了有意義的研究：Jin-Keun Choi 等［1］聯合采用面元法與求解 UN-RANS（Unsteady Reynolds Averaged Navier-Stokes equation）的方法分別計算了噴氣下瘦長型船的興波、船底氣層形態、粘性阻力；D.Kim 等［2］采用直接數值模擬（DNS）的方法計算了斷階后氣層的非定常流動形態，并與試驗結果進行對比，數值結果較為精確的與試驗結果吻合；Hoang Cong Liem 等［3］結合邊界層積分方程及經驗公式，計算了噴氣對船舶阻力的影響，該方法可用于計算實尺度氣泡船的減阻效果；李云波等［4 - 6］基于勢流理論，采用邊界元方法計算了船底凹槽內的氣層形態；王家楣等［7］基于 PHOENICS 軟件中的 IPSA（相間滑移模型）計算了二維氣泡船的噴氣減阻效果，計算結果與試驗結果在趨勢上一致；蔡紅玲等［8 - 10］開展了高速氣泡船流場數值模擬，在Fluent平臺上采用氣液兩相流模型進行計算，所得結果在規律上合理，但是由于在船體表明僅采用非結構網格，因而不能較好地模擬邊界層流動。

上述研究均為湍流模型、網格、壁面函數等進行系列探討，所得結果尚需提高。本論文針對氣液混合流，系統研究湍流模型、網格、壁面函數對計算結果的影響，最終形成可較好模擬氣液混合流的數值方案，對肥大型氣泡船氣層形態及減阻率預報具有指導意義。

1　Mixture 兩相流模型

Mixture 模型的控制方程由連續性方程和動量方程構成，其具體表達形式如下：

連續性方程為：

混合相的動量方程由各分相的動量方程相加得到，可表示為：

該模型同時考慮了相間的滑移速度，即第 p 相相對于第 q 相的滑移速度，根據 Manninen等的研究結果，滑移速度可表示為：

2　湍流模型及壁面處理方法

采用光滑平板為對象研究氣液混合流數值模型的構建方法。平板的外形及尺寸如圖 1 所示。平板總長1 200 mm，寬 380 mm，厚度 10 mm，在距離頭部 410 mm處安裝噴氣板。噴氣板寬 120 mm，長 70 mm。

圖 1　計算平板主尺度示意圖Fig. 1　The main parameters of plate

計算流域的網格布局如圖 2 所示。流域總長為 6倍板長，總寬為 4 倍板寬。流域入口距離平板頭部 1倍板長，設置為速度入口；出口距離平板尾部 4 倍板長，為靜壓力出口；噴氣入口設置為質量流量入口。由于氣體不可壓縮，則根據質量流量可換算得到相應的體積流量。

圖 2　平板流場區域及網格布局Fig. 2　The flow field of plate and its mesh distribution

眾多研究表明［11 - 15］：只有當氣泡分布在湍流邊界層內時才具備較好的減阻效果，因此采用 Mixture 模型計算混合流邊界層時，壁面附近網格的處理方法較為關鍵。

目前，在船舶粘性流場計算中，壁面的處理方法主要有標準壁面函數法、增強壁面函數法及無壁面函數的直接計算方法。

根據湍流模型、壁面函數的自身特征以及數值計算經驗，壁面附近網格的處理辦法一般基于如下原則：

1）標準壁面函數法。y+應大于 30～60，最好接近30，湍流邊界層內應布置一定的網格數量。標準壁面函數不適用于層流底層，因此第 1 層網格應布置在對數律層內。

2）增強壁面函數法。該方法可用于計算層流底層的流動，當用于計算層流底層時，y+最好取值為 1 左右，但如果第 1 層網格布置在了層流層內，略高的 y+值也可行。計算過程中，層流底層內應布置至少 10 個網格節點。

3）對于 S-A 模型（Spalart-Allmaras），若采用增強壁面函數，則 y+= 1；若采用標準壁面函數，則 y+＞ 30。

4）對于低雷諾數 k-ω 模型，網格要求與增強壁面函數相同；采用高雷諾數 k-ε 模型時，網格要求與標準壁面函數相同。

5）對于 LES 模型，壁面附近的網格要求與增強壁面函數的要求相同。

本文計算了 8 種壁面函數及湍流模型的組合方案，如表 1 所示。表中 y+值所對應的來流速度為 1.241 m/s。

表 1　湍流模型及壁面函數Tab. 1　The turbulence models and wall functions

3　氣層形態分析

3.1氣層宏觀形態

采用上述湍流模型及壁面函數計算了氣流量 Q = 0.5 m3/h，速度 V = 1.241 m/s 時平板底部的氣層形態，所得結果如圖 3～圖 10 所示。圖 11 為相同氣流量及來流速度下氣層的模型試驗觀測結果。

從氣層的分布上看：計算方案 1 與計算方案 3 所得結果與試驗結果在變化規律上相差甚遠，說明使用S-A 模型與 LES 模型模擬氣液混合流時具有較大局限性，在本文建立的網格方案下，尚需進行優化。計算方案 2、方案 4～方案 8 在氣層的分布規律上與試驗結果較一致，但各方案之間存在差別。

從氣泡濃度與氣層流態上看：計算方案 2、方案 3與方案 5 所得氣體在壁面上的分布濃度大于 0.95，平板表面上幾乎形成了分層流；計算方案 7 與計算方案8 所得氣體在壁面上的濃度為 0.5～0.8，主要為氣體與液體的混合流動。從圖 11 的試驗圖片可以看出：氣體主要氣泡群的形式存在，未形成分層流。說明計算方案 7 或計算方案 8 所得結果與試驗結果較為接近。

圖 3　計算方案 1（S-A 模型，y+= 3）Fig. 3　The first calculation case by S-A turbulence model at y+= 3

圖 4　計算方案 2（S-A 模型，y+= 10）Fig. 4　The second calculation case by turbulence S-A model at y+= 10

圖 5　計算方案 3（k-ω 模型）Fig. 5　The third calculation case by k-ω model

圖 6　計算方案 4（LES 模型）Fig. 6　The fourth calculation case by LES model

從圖 7～圖 9（計算方案 5～方案 7）還可以看出：在湍流模型及壁面函數不變的情況下，壁面法向的第 1 層網格厚度（y+值）對平板表面氣泡濃度的計算結果有較大影響。當氣流量及來流速度相同時，第1 層網格厚度增加，壁面上的氣泡濃度降低。

圖 7　計算方案 5（RNG k-ε 模型及增強壁面函數，y+= 2）Fig. 7　The fifth calculation case by RNG k-ε model and enhance wall function at y+= 2

圖 8　計算方案 6（RNG k-ε 模型及增強壁面函數，y+= 10）Fig. 8　The sixth calculation case by RNG k-ε model and enhance wall function at y+= 10

圖 9　計算方案 7（RNG k-ε 模型及增強壁面函數，y+= 20）Fig. 9　The seventh calculation case by RNG k-ε model and enhance wall function at y+= 20

圖 10　計算方案 8（RNG k-ε 模型及標準壁面函數，y+= 30）Fig. 10　The eighth calculation case by RNG k-ε model and standard wall function at y+= 30

圖 11　平板底部氣層模型試驗照片

3.2氣層厚度及寬度

采用 3.1 節的計算方案 8 計算不同流動工況下平板底部的氣層，并將計算所得氣層厚度與實驗結果進行對比（試驗中氣層厚度的測量見文獻［6］），如圖12～圖 14 所示。圖中橫坐標 x 表示與平板頭部的距離，其中 x = 415 mm 處為噴氣入口；縱坐標表示沿流動方向不同位置處的氣層厚度。

從圖 12～圖 14 可看出，采用 Mixture 兩相流模型與 RNG k-ε 湍流模型相結合的數值方法，可有效模擬氣層厚度沿流動方向的變化規律。

圖 12　V = 1.287 m/s，Q = 1.08 m3/h 氣層厚度沿平板長度方向變化的計算值與實驗結果Fig. 12　The air layer thickness of numerical results and test result at V = 1.287 m/s and Q = 1.08 m3/h

圖 13　V = 0.838 m/s，Q = 1.11 m3/h 氣層厚度沿平板長度方向變化的計算值與實驗結果Fig. 13　The air layer thickness of numerical result and test result at V = 0.838 m/s and Q = 1.11 m3/h

圖 14　V = 1.287 m/s，Q = 2.061 m3/h 氣層厚度沿平板長度方向變化的計算值與實驗結果Fig. 14　The air layer thickness of numerical result and test result at V = 1.287 m/s and Q = 2.06 m3/h

圖 15 和圖 16 給出了不同流動工況下氣層寬度數值計算結果與實驗結果的對比。圖中 x = 0 mm 處為噴氣入口；縱坐標表示在流動方向上的氣層寬度。

從圖 15 和圖 16 可看出：采用 Mixture 模型與 RNG k-ε 湍流模型相結合的數值方法可有效模擬氣層沿平板寬度方向的擴散邊界。

圖 15　V = 1.287 m/s，Q = 2.061 m3/h 沿平板長度方向氣層厚度的計算值與實驗結果Fig. 15　The air layer thickness of numerical result and test result at V = 1.287 m/s and Q = 2.06 m3/h

圖 16　V = 0.837 m/s，Q = 1.04 m3/h 沿平板長度方向氣層厚度的計算值與實驗結果Fig. 16　The air layer thickness of numerical result and test result at V = 0.837 m/s and Q = 1.04 m3/h

4　結　語

1）采用 RNG k-ε 湍流模型、標準壁面函數、第 1層網格 1 mm、y + 為 31～35 的計算方案可較好地模擬平板底部氣液混合層的宏觀形態。

2）采用 Mixture 模型與 RNG k-ε 湍流模型相結合的方法可有效模擬氣層沿平板寬度方向的擴散邊界和氣層厚度沿流動方向的變化規律，說明搞方法對模擬氣液混合流氣層形態有效。

CHOI J K，HSIAO C T，CHAHINE G L. Numerical studies on

［1］ the hydrodynamic performance and the startup stability of high speed ship hulls with air plenums and air tunnels［C］//Proceedings of the Ninth International Conference on Fast Sea Transportation FAST2007. Shanghai，China，2007：476-484.

［2］KIM D，MOIN P. Direct numerical Simulation of air layer drag reduction over a backward-facing step［C］//Proceedings of the 63rd annual Meeting of the APS Division of Fluid Dynamics. Long Beach，California：American Physical Society，2010：351-363.

［3］LIEM H C，TODA Y，SANADA Y. A consideration on drag reduction by air lubrication using integral type boundary layer computation［J］. Journal of the Japan Society of Naval Architects and Ocean Engineers，2011，13：59-65.

［4］LI Y B，WU X Y，MA Y，et al. A method based on potential theory for calculating air cavity formation of an air cavity resistance reduction ship［J］. Journal of Marine Science and Application，2008，7（2）：98-101.

［5］吳曉宇. 氣泡船大尺寸多凹槽穩定氣穴形成理論研究［D］. 哈爾濱：哈爾濱工程大學，2008.

［6］董文才，郭日修，朱鳳榮，等. 平板湍流邊界層內氣泡流流動實驗研究［J］. 海軍工程大學學報，2001，13（3）：34-37.

［7］鄭曉偉，王家楣，曹春燕. 二維船舶微氣泡減阻數值模擬［J］.船舶工程，2005，27（6）：15-18.

［8］蔡紅玲. 高速氣泡船流場數值模擬［D］. 武漢：武漢理工大學，2008.

［9］楊鵬. 氣泡船三維粘性繞流的數值模擬［D］. 武漢：武漢理工大學，2008.

［10］曹春燕. 船舶微氣泡減阻數值模擬［D］. 武漢：武漢理工大學，2003.

［11］荒賀浩一，松井良輔，脇本辰郎，等. 界面活性剤水溶液の水平円管內流れに及ぼす微細気泡の影響［J］. 実験力學，2010，10（3）：304-311.

［12］梁志勇，陳池. 微氣泡對平板摩擦阻力影響的分析［J］. 上海大學學報（自然科學版），2002，8（3）：267-272.

［13］FELTON K，LOTH E. Diffusion of spherical bubbles in a turbulent boundary layer［J］. International Journal of Multiphase Flow，2002，28（1）：69-92.

［14］郭峰，歐勇鵬，董文才，等. 平板微氣泡減阻預報及影響因素研究［J］. 中國造船，2008，49（S）：66-74.

［15］MORIGUCHI Y，KATO H. Influence of microbubble diameter and distribution on frictional resistance reduction［J］. Journal of Marine Science and Technology，2002，7（2）：79-85.

Numerical study of method of flat plate viscous flow field with bubble

WU Hao，OU Yong-peng
（Department of Naval Architecture Engineering，Naval University of Engineering，Wuhan 430033，China）

In order to investigate the complexity of gas-liquid mixing under different condition of flat plate，a method with the combination of RANS equations and Mixture model is proposed for the viscous-flow calculation of a large flat bottom ship. There are eight different numerical calculation method formed by a combination 4 kinds of turbulence models，grid，and 3 kinds of wall treatment. The influence of wall function，the first layer of mesh wall and turbulence models for numerical results was analyzed. Experimental results were compared with the numerical study. The results show that：RNG k-ε turbulence model，standard wall function，1mm first layer of mesh，y + for the calculation of 31 to 35，the results can be used for mixed-flow analysis.

gas-liquid mixing；Mixture model；wall function；turbulence model

U631.1

A

1672 - 7619（2016）08 - 0047 - 05

10.3404/j.issn.1672 - 7619.2016.08.010

2015 - 09 - 11；

2015 - 10 - 13

工信部高技術船舶科研資助項目（［2011］530）；高性能船舶技術教育部重點實驗室開放基金資助項目（2013033102）

吳浩（1987 - ），男，博士研究生，研究方向為船舶水動力學。