余澤爽 毛筱菲

(武漢理工大學(xué)交通學(xué)院1) 武漢 430063) (武漢理工大學(xué)教育部高性能船重點實驗室2) 武漢 430063)

基于重疊網(wǎng)格M型船水氣二相流數(shù)值模擬研究

余澤爽1,2)毛筱菲1,2)

(武漢理工大學(xué)交通學(xué)院1)武漢 430063) (武漢理工大學(xué)教育部高性能船重點實驗室2)武漢 430063)

M型船是一種新型滑行艇，同時具備普通滑行艇的阻力特性、高速多體船的耐波性能，以及微氣泡減阻船優(yōu)異的阻力性能.為了實現(xiàn)對M型船不同航行階段中復(fù)雜變航態(tài)過程的動態(tài)捕捉、對阻力的精確預(yù)報，以及對非排水航行階段復(fù)雜的水氣二相流場的模擬，基于粘性CFD軟件STAR-CCM+的重疊網(wǎng)格技術(shù)，對M型船在靜水中從排水、過渡到滑行航行階段進(jìn)行數(shù)值模擬，從船體表面網(wǎng)格大小討論了網(wǎng)格獨立性，結(jié)合計算資源以確定最佳的網(wǎng)格方案；將阻力和縱傾角的計算值與試驗值進(jìn)行對比，從湍流模型與y+值以及時間步長獨立性討論了不同航行階段的模擬需求；從水氣二相流的模擬討論了M型船在不同航行階段的減阻原理.

M型船；重疊網(wǎng)格；水氣二相流；湍流模型；時間步長適應(yīng)性

0 引 言

M型船是為避免游船在高速航行時產(chǎn)生的興波對河岸的沖刷破壞而發(fā)展的新型消波船型，因其橫剖面形狀為M形而得名，由此還衍生出W型船、雙M型船等船型.M型船因其獨特的幾何形狀和優(yōu)秀的流體動力性能，綜合了普通滑行艇、高速多體船和微氣泡減阻船等高性能船的優(yōu)點，能夠充分利用滑行效應(yīng)、浮力、表面效應(yīng)，以及首波與空氣捕捉等原理，產(chǎn)生高效的浮力、水氣動升力來支持艇重，從而降低船體的總阻力、提高航速.M型船由中心排水體(CDS)、滑行涵道(PT)和剛性側(cè)裙(RS)組成.

雙M槽道滑行艇(短劍Stiletto)在美國進(jìn)行了試航和驗收，并交付近海特殊作戰(zhàn)部隊，此外還是首次采用了碳纖維濕法纏繞環(huán)氧復(fù)合材料的整體成型技術(shù)的新型高速艇[1].遲玉鵬等[2]通過系列船模試驗，對槽道艇和單體V型艇阻力和耐波性能進(jìn)行了對比分析，發(fā)現(xiàn)槽道艇具有橫搖小、噴濺少和垂向加速度小的優(yōu)點，并提出了阻力最優(yōu)的重心縱向位置.蘇永昌等[3]通過系列試驗，對槽道滑行艇的基本原理、運動特性及其變化規(guī)律進(jìn)行了研究，并重點討論了船體型線、槽道尺寸和形狀，以及重心縱向位置等影響因素.陳輝[4]基于Fluent對M型船靜水阻力進(jìn)行數(shù)值模擬，并與試驗結(jié)果進(jìn)行對比，在非滑行階段(Fr▽<3)較為精確，但當(dāng)Fr▽>3時數(shù)值計算誤差較大.唐建飛等[5-6]通過系列船模試驗，討論了排水量和重心縱向位置對M船型阻力和航態(tài)的影響，提出了該船型存在兩個阻力峰值的特性，同時還比較了M船型與普通槽道滑行艇的阻力與耐波性能，提出M船型優(yōu)良的水動力性能.Ghassabzadeh等[7]采用動網(wǎng)格技術(shù)(dynamic mesh)，考慮了自由面水氣二相流場，對槽道多體船不同航行階段進(jìn)行了模擬，實現(xiàn)了對高速運動的模擬和流場信息的捕捉.鄒勁等[8]基于STAR-CCM+，討論了船表網(wǎng)格尺寸、近船面網(wǎng)格節(jié)點分布、船體周圍加密區(qū)網(wǎng)格尺寸以及變形網(wǎng)格技術(shù)等對三體滑行艇的阻力預(yù)報、收斂速度以及運動穩(wěn)定性的影響.

綜上，M型船數(shù)值模擬研究的難點在于：①航行過程中大幅度的航態(tài)變化，對網(wǎng)格要求極高，普通網(wǎng)格已無法滿足模擬需求；②M型船航速變化范圍很大，對不同航行階段，有不同的模擬需求；③滑行涵道中水氣兩相流場的模擬，對網(wǎng)格、壁面處理和湍流模型提出更嚴(yán)格的需求.

1 理論背景

1.1 控制方程

流體的流動需滿足質(zhì)量守恒定律，對于不可壓縮流體，即有連續(xù)性方程：

(1)

式中：u，v，w分別為x，y，z方向上的分量.

流體流動還應(yīng)滿足動量守恒定律，即由Naviver-Stokes提出的N-S方程：

(2)

Su=Fx+sx=0

Sv=Fy+sy=0

Sw=Fz+sz=ρg

(3)

直接求解N-S方程過于困難，而雷諾平均的RANS方程則提供了便利.RANS方程是平均變量的控制方程，基于湍流模式假定，將高頻脈動的流場變量分解為時均量與脈動量之和.

1.2 自由面捕捉

自由面捕捉采用的VOF方法，根據(jù)網(wǎng)格單元中不同流體的體積分?jǐn)?shù)的變化來模擬自由面的位置和形狀，流體的體積分?jǐn)?shù)即為流體體積與網(wǎng)格單元體積之比：

(4)

式中：f1，f2分別為空氣相和水相的體積分?jǐn)?shù)，其控制方程為

(5)

以體積分?jǐn)?shù)f表示時，f=1即為水相，f=0即為空氣相，f=0.5即為自由面.

1.3 湍流模型

由于湍流現(xiàn)象的高度復(fù)雜性，所以至今未有一種通用的方法能夠精確模擬所有流動問題中的湍流現(xiàn)象.雷諾平均的N-S方程雖然解決了方程的變量問題，但是計算中涉及的湍流問題，都是以湍流模型的形式添加到方程中的[9].

Star-CCM+中提供了k-ε模型、SSTk-ε模型、Spalart-Allmaras模型和雷諾應(yīng)力湍流模型，本文使用和討論了前兩個模型.

1.4 邊界層的處理

粘性流體繞流靜止的壁面時，由于物面為不可滑移面，故在緊靠物面的薄層內(nèi)，流速由物面上的零增加到與來流類同，這種在高雷諾數(shù)下，在緊靠物面處形成的薄層即為邊界層.在邊界層中，速度梯度很大，湍流增強，因此壁面處網(wǎng)格的處理對捕捉該處流場的物理特性格外重要.

在使用湍流模型時，一般采用壁面函數(shù)法對壁面處網(wǎng)格做特殊處理，.通常采用無因次距離y+來定義船體表面第一層網(wǎng)格厚度，可采用經(jīng)驗公式進(jìn)行估算

(6)

式中：Δy為船體表面第一層網(wǎng)格厚度；L為濕體長度.

2 數(shù)值模擬與試驗驗證

2.1 模型試驗

本文采用的模型和實驗數(shù)據(jù)均來自文獻(xiàn)[5-6]，圖1為該M型船橫剖面圖、使用CATIA三維建模軟件建立的M型船三維模型、實驗用船模圖.M型船主要參數(shù)及實驗數(shù)據(jù)見表1～2.

圖1 M型船

表1 M型船主要參數(shù)

表2 M型船模型阻力試驗數(shù)據(jù)

2.2 計算設(shè)置

圖2為計算域和網(wǎng)格區(qū)域的劃分情況，分別表示了重疊網(wǎng)格區(qū)域、船表周圍加密區(qū)、自由面加密區(qū)以及為了捕捉船后興波的興波加密區(qū).

圖2 區(qū)域劃分情況

2.3 網(wǎng)格獨立性研究

由重疊網(wǎng)格特點，重疊網(wǎng)格區(qū)域網(wǎng)格尺寸與船體表面網(wǎng)格一致，背景網(wǎng)格為重疊網(wǎng)格的四倍，由于重疊網(wǎng)格與背景網(wǎng)格之間采用線性插值，因此船表周圍加密區(qū)網(wǎng)格為重疊網(wǎng)格兩倍，自由面加密區(qū)網(wǎng)格(Z向)為重疊網(wǎng)格的80%，船后興波加密區(qū)(XY向)為重疊網(wǎng)格的2倍.表3為網(wǎng)格數(shù)量隨船表網(wǎng)格尺寸變化得到的不同網(wǎng)格方案.

表3 網(wǎng)格數(shù)量隨船表網(wǎng)格尺寸變化

選取k-ε湍流模型，邊界層以固定壁厚給定，計算時，保持y+值為30～200，時間步長取為0.001 s，時間離散為二階形式.圖3為不同網(wǎng)格方案計算的阻力與縱傾角隨網(wǎng)格數(shù)量變化的趨勢與試驗值的對比.

圖3 不同網(wǎng)格方案計算值隨網(wǎng)格因子的變化趨勢

網(wǎng)格因子RGi為單元網(wǎng)格體積之比，即網(wǎng)格數(shù)量反比的三次方根，其表達(dá)式為

RGi=(Ni/N1)1/3

(7)

式中：Ni為方案i的網(wǎng)格數(shù)量.

由圖3可知，計算結(jié)果隨網(wǎng)格數(shù)量增加呈現(xiàn)單調(diào)趨勢，說明網(wǎng)格收斂性得到驗證.通過比較數(shù)值變化趨勢，綜合考慮計算資源，認(rèn)定方案3作為后續(xù)研究的網(wǎng)格方案，因為在此方案基礎(chǔ)上，增加網(wǎng)格數(shù)量對計算結(jié)果的影響并不顯著，即網(wǎng)格獨立性得到滿足.

2.4 湍流模型與y+值適應(yīng)性研究

2.4.1k-ε湍流模型

選用k-ε湍流模型時，分別取y+為50，100和200，根據(jù)式(6)估算得到船表第一層網(wǎng)格高度Δy，取邊界層數(shù)為5，厚度增長率為1.3，則可得到邊界層厚度δ.表4為采用k-ε湍流模型時，不同航速下阻力和縱傾角隨y+值變化的計算值與試驗值的誤差.

表4 y+值對阻力和縱傾角的影響

由表4可知，在排水航行階段(Fr▽<1)，y+值對計算結(jié)果的影響并不大，與實驗結(jié)果的誤差均保持在5%以內(nèi)，橫向比較可知選用y+=200為最佳；在過渡航行階段(13)，誤差則存在極大不確定性，可見在該航行階段，k-ε湍流模型已無法滿足模擬需求，因此需要考慮SSTk-ε湍流模型.

2.4.2SSTk-ε湍流模型

船舶在流體中航行時，隨著航速的提高，雷諾數(shù)Re也隨之提高，邊界層則隨之變薄，由壁面函數(shù)公式(6)可知，此時的y+值必然降低，而由文獻(xiàn)[9]的相關(guān)研究，在選用SSTk-ε湍流模型時，y+應(yīng)不大于1，故這里選取y+為0.5和1.表5為選用SSTk-ε湍流模型時，不同航速下阻力與縱傾角隨y+值變化的計算值與試驗值的誤差.

表5 y+值對阻力和縱傾角的影響

由表5可知，在非滑行航行階段(Fr▽<3)，計算結(jié)果整體的誤差均較大，不如采用k-ε湍流模型時的結(jié)果；而在滑行航行階段(Fr▽>3)，當(dāng)y+=1時的結(jié)果誤差較小低于8%，故在采用k-ε湍流模型模擬滑行航行階段時，y+應(yīng)接近1為最佳.

2.5 時間步長適應(yīng)性研究

在數(shù)值模擬中，時間步長的選取對計算結(jié)果會造成一定影響，一個合適的時間步長能使非定常的數(shù)值模擬得到同時滿足準(zhǔn)確性、穩(wěn)定性和計算成本的需求[10-11].

這里主要討論滑行航行階段采用k-ω湍流模型的時間步長適應(yīng)性，這里選取時間步長為0.01,0.005,0.0025,0.001和0.000 5 s，時間離散采用二階形式.表6為采用k-ε湍流模型時，滑行航行階段阻力和縱傾角隨時間步長變化的計算值與試驗值的比較.

表6 時間步長對阻力和縱傾角的影響

由表6可知，隨著時間步長的減小，計算結(jié)果誤差并不一定越小，但是結(jié)果基本上呈現(xiàn)單調(diào)的收斂關(guān)系.同時，對于高速情況，時間步長不宜過大，否則會出現(xiàn)發(fā)散的情況，這也一定程度說明k-ε湍流模型穩(wěn)定性較差的特點，此外，綜合分析阻力與縱傾角角的計算值與試驗值的誤差為標(biāo)準(zhǔn)，再考慮到時間步長對計算收斂過程的影響，兼顧計算資源的經(jīng)濟(jì)性，選擇最佳時間步長為0.005 s，時間離散為二階形式.

2.6 試驗驗證

由上述網(wǎng)格獨立性、湍流模型、y+值以及時間步長適應(yīng)性的研究，對M型船從排水、過渡到滑行航行階段的靜水直航運動進(jìn)行數(shù)值模擬，并將計算的阻力與縱傾角與試驗結(jié)果對比見圖4.在全航行階段，采用上述數(shù)值模擬方法計算的阻力與縱傾角，與試驗結(jié)果基本吻合，誤差均在12%以內(nèi).

圖4 數(shù)值模擬與試驗驗證

3 M型船特性研究

3.1 水氣兩相流場特性分析

由于滑行涵道特殊的喇叭形結(jié)構(gòu)，M型船在航行過程中，中心排水體艏部產(chǎn)生的興波，隨空氣一同進(jìn)入涵道，而涵道出口一般僅僅略高于靜水面，此外高速航行時，滑行面與液面的劇烈摩擦，產(chǎn)生大量水氣，此時，涵道中形成復(fù)雜的水氣混合物，并且壓力逐漸上升，最終在艉部以噴濺的形式從出口溢出，形成集中的尾興波.

圖5為不同航速時滑行涵道中水氣分布圖.由圖5可知，隨著航速的提高，滑行涵道中水氣混合物逐漸減少，同時對比不同位置處的水氣分布情況，可見接近船尾的水氣分布更密集更復(fù)雜，這是因為高速情況下，由于涵道中急劇升高的壓強，水氣混合物從出口壓出.

圖5 不同航速時滑行涵道中水氣分布

圖6為不同航速時船體表面水氣分布圖.由圖6可知，在排水和過渡航行階段，隨著航速的提高，船體表面水氣分布呈現(xiàn)有規(guī)則的水氣分層現(xiàn)象，一方面是由于航態(tài)縱傾，另一方面是由于水氣混合物此時主要集中于涵道中，并未出現(xiàn)劇烈的溢出現(xiàn)象；而在滑行航行階段，船表水氣分布顯得更為紊亂，這是由于涵道中的水氣混合物以噴濺的形式從出口溢出，同時滑行面與液面之間高速的相對摩擦運動，引起近壁面處出現(xiàn)局部的氣化現(xiàn)象.

圖6 不同航速時船體表面水氣分布

圖7為不同航速時的船后興波圖.由圖7可知，在排水航行階段，船后興波的橫波系和散波系與普通滑行艇類似；在過渡航行階段，滑行涵道中的水氣混合物開始少量溢出，由于涵道收縮的結(jié)構(gòu)，使液面被壓縮，流體以集束的形式從出口流出，由于兩股束流之間的相互干擾，同時側(cè)壁剛?cè)箤χ行呐潘w興波的阻隔作用，船后興波相比排水航行階段要更為收縮；在滑行階段，涵道中的水氣混合物因高壓從出口急劇溢出，使船后興波進(jìn)一步收縮，而不會出現(xiàn)明顯的“空穴”和“雞尾流”現(xiàn)象.

圖7 不同航速時的船后興波

3.2 減阻原理分析

3.2.1滑行升力

與普通滑行艇類似，特殊的滑行面形狀在航行中提供一定的水動升力，使船體抬離水面，達(dá)到減阻的目的.圖8為不同航行階段船表壓力.由圖8可知，在不同航行階段，隨著航速的提高，船表壓力逐漸提升.

圖8 不同航行階段船表壓力

3.2.2表面效應(yīng)升力

由于涵道的特殊形狀，中心排水體產(chǎn)生的艏興波隨同空氣，以螺旋狀的運動方式進(jìn)入涵道，并在涵道內(nèi)被加壓、分離，形成外旋的水氣渦流，可視為縱向的二元渦，帶自由面的二元渦產(chǎn)生表面效應(yīng)升力.表面效應(yīng)升力L表達(dá)式為

L∝kρω2r3l

(8)

式中：k為系數(shù)；ω為二元渦角速度；r為二元渦半徑；l為二元渦的長度.

圖9為不同航行階段涵道中的渦流，由圖9可知，隨航速增加，渦流由船的艏部逐漸向后延伸，由式(8)可知，該升力也隨之增加，說明航速越高，槽道渦表面效應(yīng)升力越大.

圖9 不同航行階段涵道中的渦流

3.2.3水氣的潤濕作用

由圖5～6涵道和船表水氣分布情況可知，尤其在滑行階段滑行面上產(chǎn)生的氣化現(xiàn)象，能顯著降低摩擦阻力系數(shù)，從而達(dá)到減阻效果.

3.2.4抑制船行波系

由于側(cè)壁剛?cè)蛊街钡耐獗谀芙档痛w壓力梯度，減小形狀阻力，而內(nèi)壁可以屏蔽中心排水體產(chǎn)生的興波.由圖7可知，隨著航速增加，船后興波呈現(xiàn)收縮的趨勢.圖10為船后不同距離處，船行波波面高度分布圖，由圖10可知，隨航速增加，船行波波面高度分布更集中.

圖10 不同航速的船行波

4 結(jié) 論

1) 通過網(wǎng)格收斂性和獨立性研究，建議M型船數(shù)值模擬時，船表網(wǎng)格尺寸宜取為0.012 5m，即0.625%L.

2) 在排水航行階段(Fr▽<1)，宜采用k-ε湍流模型，且y+值應(yīng)取200為宜；在過渡航行階段,宜采用k-ε湍流模型,并且其前半段(13)，宜采用SSTk-ε湍流模型，且y+值應(yīng)接近1為佳.

3) 選用SSTk-ε湍流模型時，通過時間步長適應(yīng)性研究，推薦選用時間步長為0.005 s，時間離散為二階.

4) 通過對M型船水氣二相流的數(shù)值模擬研究，發(fā)現(xiàn)滑行涵道中的水氣混合的現(xiàn)象，主要發(fā)生在排水和過渡航行階段，而滑行階段主要是在滑行面上發(fā)生了氣化現(xiàn)象.

[1] 趙連恩,王慶.高性能船(艦)原理與技術(shù)[M].哈爾濱：哈爾濱工程大學(xué)出版社,2013.

[2] 遲云鵬,孟憲欽.高速槽道艇阻力及耐波性能試驗研究[J].船舶工程,1995(3):27-31.

[3] 蘇永昌,趙連恩.高性能槽道滑行艇的運動特性[J].中國造船,1996(1):11-16.

[4] 陳輝.M船型水氣兩相流場特性研究[D].武漢：武漢理工大學(xué),2011.

[5] 唐建飛,黃武剛.M船型阻力模型試驗研究[J].中國艦船研究,2014(5):49-52.

[6] 黃武剛.M型艇與槽道型艇的阻力和耐波性比較[J].船海工程,2015,44(1):56-59.

[7] GHASSABZADEH M, GHASSEMI H. Determining of the hydrodynamic forces on the multi-hull tunnel vessel in steady motion[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering,2014,36(4):697-708.

[8] 鄒勁,姬朋輝,孫寒冰,等.網(wǎng)格因素對三體滑行艇阻力計算影響探究[J].船舶,2016,27(3):8-14.

[9] 章麗麗,孫寒冰,蔣一,等.三體滑行艇槽道的水氣動力特性研究[J].哈爾濱工程大學(xué)學(xué)報,2017,38(1):31-36.

[10] 王福軍.計算流體動力學(xué)分析：CFD軟件原理與應(yīng)用[J].熱能動力工程,2005(5):468-468.

[11] 孫華偉.滑行面形狀對滑行艇阻力與航態(tài)影響數(shù)值分析[D].哈爾濱：哈爾濱工程大學(xué),2012.

Research of the RANS Simulation of Two-phase Flow Around M-hull Ship Based on Overset Grid

YUZeshuang1,2)MAOXiaofei1,2)

(SchoolofTransportation,WuhanUniversityofTechnology,Wuhan430063,China)1)(KeyLaboratoryofHighPerformanceShipTechnologyofMinistryofEducation,WuhanUniversityofTechnology,Wuhan430063,China)2)

M-hull ship is a new kind of planning craft with the resistance characteristics of the normal planning craft, the sea-keeping performance of the high-speed multi-hull ship and the excellent resistance performance of the micro-bubble drag reduction ship. In order to achieve the goal to dynamically capture the complex changing navigation state during different navigation phases, accurately predict the resistance performance and simulate the two-phase flow during non-navigation phase of displacement, the RANS simulation from the navigation phase of displacement to transition and planning in the static water is carried out based on the overset gird technology of Star-CCM+. The mesh independence study is performed according to the grid size closed to the ship surface and an optimal mesh scheme is determined combined with the computing resource. Comparing the compute values of the resistance and trim angle and the experiment values, the influence between turbulence model and they+value, and the time step adaptive is discussed during different navigation phases. The drag reduction principle during different navigation phase is also discussed according to the simulation of the two-phase flow around the M-hull ship.

M-hull ship; overset grid; two-phase flow; turbulence model; time step adaptive

U661.32

10.3963/j.issn.2095-3844.2017.05.030

2017-08-31

余澤爽(1992—)：男，碩士生，主要研究領(lǐng)域為船舶性能