董洪輝，王寶壽，陳瑋琪，張珂

（中國(guó)船舶科學(xué)研究中心，江蘇無(wú)錫214082）

泡狀流中空化與激波相互作用的數(shù)值模擬

董洪輝，王寶壽，陳瑋琪，張珂

（中國(guó)船舶科學(xué)研究中心，江蘇無(wú)錫214082）

文章針對(duì)定?？蓧嚎s泡狀流中回轉(zhuǎn)體的空化繞流進(jìn)行了數(shù)值模擬，在給定來(lái)流速度和環(huán)境壓力的條件下，對(duì)不同含氣率下的空化與激波的相互作用進(jìn)行了研究。首先，對(duì)含氣率為0的情形進(jìn)行計(jì)算并同試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，證明所采用的計(jì)算模型是可信的。其次，改變流場(chǎng)含氣率從0至0.5進(jìn)行計(jì)算，結(jié)果表明，隨著含氣率的增大，流場(chǎng)的可壓縮性隨之增強(qiáng)，體現(xiàn)為數(shù)值模擬得到了回轉(zhuǎn)體繞流流場(chǎng)的3種波系，包括回轉(zhuǎn)體頭部處的脫體激波、分離面處的膨脹波和空泡尾端的斜激波。而空化與激波的相互作用則體現(xiàn)為：頭激波削弱了空化效應(yīng)，使得空化區(qū)域減小；由于空泡外形的影響，使得空泡尾端出現(xiàn)了斜激波；含氣率超過(guò)一定值，空泡已經(jīng)無(wú)法閉合在物面上，而是閉合在空泡尾端的激波面上，體現(xiàn)為空泡尾端壁面逆壓梯度趨于平緩。計(jì)算結(jié)果揭示了泡狀流中空化與激波相互作用的新的物理現(xiàn)象。

泡狀流；激波；空化

0 引言

均勻混合的氣液兩相流具有可壓縮性，Karplus等人[1]的試驗(yàn)研究表明，當(dāng)偏離單相極點(diǎn)時(shí)，氣液混合物的聲速迅速降低，特別是當(dāng)氣相體積分?jǐn)?shù)為0.5時(shí)，氣液兩相流的聲速接近最小，只有幾十米每秒，如圖1為不同壓力下聲速與含氣率的關(guān)系。因此，在等溫條件下，泡狀流中就可以產(chǎn)生激波，此類激波稱為動(dòng)力激波。

Eddington[3]對(duì)氣液兩相流中的超音速現(xiàn)象進(jìn)行了研究，得到并證實(shí)了當(dāng)氣相體積分?jǐn)?shù)為0.5時(shí)泡狀流中的正激波基本關(guān)系。Jhon等人[4]在等溫條件下，忽略相間速度滑移，提出了一種勢(shì)流線性分析方法，對(duì)在泡狀流中高速移動(dòng)的擾動(dòng)點(diǎn)所產(chǎn)生激波的波型進(jìn)行了研究。Saito等人[2]使用單相流模型通過(guò)求解可壓縮N-S方程組，對(duì)NACA0015水翼的三維非定常云空化繞流進(jìn)行了數(shù)值模擬。陳紅兵[5]使用速度入口條件計(jì)算了水翼在氣相體積分?jǐn)?shù)為0.5的來(lái)流中的跨聲速流動(dòng)，在超聲速來(lái)流條件下得到了水翼前端的弓形脫體激波，雖然計(jì)算中使用了空化模型但是沒(méi)有觀察到空化現(xiàn)象。

圖1 聲速與含氣率的關(guān)系[2]Fig.1 Relation between the speed of sound and the void fraction[2]

可見(jiàn)，對(duì)泡狀流動(dòng)的可壓縮性的理論和數(shù)值研究，多集中于氣相體積分?jǐn)?shù)為0.5的情形，尚沒(méi)有對(duì)不同含氣率情況下泡狀流中空化現(xiàn)象開(kāi)展研究，因而也沒(méi)有對(duì)空化與激波的相互影響進(jìn)行研究。本文基于CFD并在給定來(lái)流速度的條件下，對(duì)回轉(zhuǎn)體在不同氣相體積分?jǐn)?shù)下的泡狀繞流場(chǎng)中空化與激波的相互作用進(jìn)行了數(shù)值模擬。

1 多相流基本控制方程

本文采用商用軟件Fluent 12中Mixture多相流模型進(jìn)行數(shù)值模擬，忽略相間速度滑移[4]，其主要控制方程為[6]：

混合物連續(xù)性方程

混合物動(dòng)量方程

其中：μt為湍動(dòng)粘度，混合物粘度

混合物能量方程

氣相體積分?jǐn)?shù)方程

理想氣體狀態(tài)方程

標(biāo)準(zhǔn)k-ε湍流模型

其中：常數(shù)取值為C1ε=1.44，C2ε=1.92，Cμ=0.09，σk=1.0，σε=1.3。

Zwart-Gerber-Belamri空化模型的汽相輸運(yùn)方程

其中：Re和Rc分別表示蒸汽泡生成和潰滅的質(zhì)量源項(xiàng)：

其中：蒸汽泡半徑RB=10-6m，氣核體積分?jǐn)?shù)αnuc=5×10-4，汽化系數(shù)Fvap=50，凝結(jié)系數(shù)Fvap=0.001。

2 計(jì)算網(wǎng)格、求解器設(shè)置及驗(yàn)證

2.1 計(jì)算網(wǎng)格與求解器設(shè)置

研究對(duì)象為回轉(zhuǎn)體，其頭部接近圓錐，后體為半徑R=100 mm、長(zhǎng)度L=1 000 mm的圓柱，不計(jì)及尾部影響，將圓柱尾端延伸至流場(chǎng)出口邊界。由于流動(dòng)具有對(duì)稱性，因此計(jì)算域?yàn)槎S軸對(duì)稱網(wǎng)格，其中入口邊界距回轉(zhuǎn)體頭部50R，圓柱長(zhǎng)度L=100R，出口邊界距回轉(zhuǎn)體頭部100R，如圖2。使用ICEM CFD 12.0進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分，共劃分71 631個(gè)四邊形單元，回轉(zhuǎn)體附近網(wǎng)格采用C型網(wǎng)格，如圖2所示。

圖2 計(jì)算域與回轉(zhuǎn)體附近網(wǎng)格Fig.2 Computational domain and mesh around the body

陳紅兵[5]的數(shù)值試驗(yàn)表明，當(dāng)泡狀流中氣相本身為低亞聲速流動(dòng)或氣相馬赫數(shù)小于0.1時(shí)，速度入口條件可用于泡狀流的無(wú)窮遠(yuǎn)或入口邊界條件。因此，入口和外部邊界均采用速度入口條件，給定來(lái)流速度為V∞＝30 m/s；回轉(zhuǎn)體表面采用壁面無(wú)滑移無(wú)穿透條件；出口采用壓力出口條件；計(jì)算域的中軸線設(shè)置為軸對(duì)稱邊界條件；流場(chǎng)初始溫度T=300K，工作壓力P=138 540 Pa，空化數(shù)σ=0.3。

采用基于壓力—速度耦合的SIMPLE算法進(jìn)行求解，計(jì)算過(guò)程是非定常的，采用Mixture多相流模型，忽略各相間的速度滑移，第一相為水，第二相分別為空氣和水蒸氣，其中空氣為理想氣體。以圓柱直徑為特征長(zhǎng)度的雷諾數(shù)Re=3×106，采用標(biāo)準(zhǔn)k-ε湍流模型，近壁處理采用標(biāo)準(zhǔn)壁面函數(shù)。邊界湍流參數(shù)設(shè)置同褚學(xué)森等人[7]所采用的一致，為T(mén)urbulent Intensity=0.5%，Turbulent Viscosity Ratio=5。

空化模型采用Zwart-Gerber-Belamri模型，飽和蒸氣壓為3 540 Pa。

2.2 計(jì)算模型驗(yàn)證

來(lái)流速度V∞＝30 m/s、工作壓力P=138 540 Pa（即空化數(shù)σ=0.3），對(duì)含氣率αg=0的情形進(jìn)行了計(jì)算，并同權(quán)曉波等人[8]的試驗(yàn)結(jié)果進(jìn)行了對(duì)比，如圖3所示。

由圖3的回轉(zhuǎn)體繞流流場(chǎng)的密度云圖可以清晰地看見(jiàn)空泡外形以及尾部的回射流結(jié)構(gòu)，而壁面壓力系數(shù)曲線與試驗(yàn)結(jié)果吻合得較好，這表明所采用的計(jì)算模型是可信的。

圖3 回轉(zhuǎn)體密度云圖和壁面壓力系數(shù)曲線Fig.3 Axis-symmetrical body’s density color contour and pressure coefficient curve on wall

3 計(jì)算結(jié)果及討論

基于以上經(jīng)過(guò)驗(yàn)證的計(jì)算模型，本文保持來(lái)流速度V∞＝30 m/s、工作壓力P=138 540 Pa，計(jì)算了氣液均勻混合介質(zhì)含氣率0～0.5的流動(dòng)狀況，以期得到含氣率導(dǎo)致的可壓縮性對(duì)回轉(zhuǎn)體繞流流場(chǎng)的影響。

3.1 含氣率對(duì)繞流流場(chǎng)的影響

圖4至圖6分別給出了各含氣率下的水相體積分?jǐn)?shù)云圖、絕對(duì)壓力云圖和各含氣率下的來(lái)流馬赫數(shù)和無(wú)量綱空泡長(zhǎng)度。由圖4可見(jiàn)，不同含氣率下，水相在流場(chǎng)中的體積分布差異很大。而造成這種差異的原因正是含氣率導(dǎo)致的流場(chǎng)可壓縮性的不同：從圖5可以看到回轉(zhuǎn)體繞流流場(chǎng)的壓力分布，由于流場(chǎng)含氣，氣相在壓力高的區(qū)域受壓縮，在壓力低的區(qū)域膨脹，從而導(dǎo)致回轉(zhuǎn)體頭部駐點(diǎn)附近的高壓區(qū)水相體積分?jǐn)?shù)最高，而分離面后的低壓區(qū)水相體積分?jǐn)?shù)最低；含氣率的升高也導(dǎo)致了流場(chǎng)可壓縮性的增強(qiáng)，由圖6可見(jiàn)，來(lái)流速度不變的情況下，隨著含氣率的增大，來(lái)流馬赫數(shù)逐漸增大，當(dāng)含氣率增至0.3附近，來(lái)流馬赫數(shù)開(kāi)始大于1，故在圖4（e）、（f）和圖5（e）、（f）中可以清晰地看到回轉(zhuǎn)體頭部前端存在脫體激波。

圖4 各含氣率下水相體積分?jǐn)?shù)云圖Fig.4 Water volume color contours of each gas fraction

圖5 各含氣率下絕對(duì)壓力云圖Fig.5 Absolute pressure color contours of each gas fraction

圖4給出的是水相體積分?jǐn)?shù)云圖，分離面后水相體積分?jǐn)?shù)低的區(qū)域即是空泡，若將空泡界面定義為ρ=0.9ρm，則圖6給出了不同含氣率下無(wú)量綱空泡長(zhǎng)度。由圖6可見(jiàn)，隨著含氣率的增大，空泡長(zhǎng)度先增大后減小，其原因可以由空泡的氣體成分來(lái)解釋：空泡內(nèi)氣體由水蒸汽和空氣兩部分組成，隨著含氣率的增大，由圖5可見(jiàn)回轉(zhuǎn)體分離面后的低壓區(qū)逐漸減小，則空化效應(yīng)被減弱，因而空泡內(nèi)的水蒸汽部分逐漸變少；但是流場(chǎng)中均勻混合的空氣對(duì)空泡可起到通氣作用，隨著含氣率的增大其作用愈發(fā)明顯；含氣率較低時(shí)，通氣作用對(duì)空泡長(zhǎng)度的增益占優(yōu)，含氣率較高時(shí)，空化效應(yīng)減弱對(duì)空泡長(zhǎng)度的影響占優(yōu)，兩種因素綜合影響下，導(dǎo)致了空泡長(zhǎng)度先增大后減小。

圖6 不同含氣率下來(lái)流馬赫數(shù)和空泡無(wú)量綱長(zhǎng)度Fig.6 Incoming flow Mach number and dimensionless cavitylength of each gas fraction

3.2 空化與激波的相互作用

上一節(jié)中，我們已經(jīng)知道隨著含氣率的升高，流場(chǎng)可壓縮性增強(qiáng)，流場(chǎng)中將出現(xiàn)激波。現(xiàn)在，讓我們來(lái)研究當(dāng)流場(chǎng)可壓縮性增強(qiáng)至出現(xiàn)激波時(shí)，激波與空化的相互作用。

圖7 各含氣率下Ma數(shù)云圖Fig.7 Mach number color contours of each gas fraction

圖7給出了各含氣率下的Ma數(shù)云圖，從中可以對(duì)比流場(chǎng)的可壓縮性（其中為了便于進(jìn)行對(duì)比，云圖標(biāo)尺是相同的）。圖7（a）表明，含氣率為0時(shí)，流場(chǎng)基本是不可壓縮的；圖7（a）至（f）表明，隨著含氣率的增大，來(lái)流馬赫數(shù)逐漸增大，當(dāng)αg=0.5時(shí)，可以清楚地看到回轉(zhuǎn)體頭部出現(xiàn)脫體激波，這說(shuō)明隨著含氣率的升高，來(lái)流流場(chǎng)的可壓縮性逐漸增強(qiáng)；由圖7（b）-（f），可以看到頭部分離面處都出現(xiàn)了膨脹波，而空泡尾部出現(xiàn)了激波，這也從另一個(gè)側(cè)面說(shuō)明氣液摻混的混合介質(zhì)具有了可壓縮性。對(duì)比圖7（b）和（f），發(fā)現(xiàn)含氣率為0.1時(shí)的Mamax比含氣率為0.5時(shí)要大，這是由于空化效應(yīng)和含氣率的相互作用導(dǎo)致的：如圖1所示，氣液混合物的聲速是壓力和含氣率的函數(shù)，壓力越低、含氣率越接近于0.5，聲速越低；當(dāng)含氣率為0.1時(shí)，回轉(zhuǎn)體頭部還未出現(xiàn)頭激波，分離面后的空化效應(yīng)比較明顯，該處壓力僅為水的飽和蒸汽壓量級(jí)并且局部含氣率接近0.5；當(dāng)含氣率為0.5時(shí)，回轉(zhuǎn)體頭部存在頭激波，波后壓力較高；因此，含氣率為0.1時(shí)流場(chǎng)中的Mamax比含氣率為0.5時(shí)要大。另外，由圖7（b）也驗(yàn)證了氣液混合物的聲速與含氣率的關(guān)系—在回轉(zhuǎn)體頭部駐點(diǎn)處水相占支配地位，而在回轉(zhuǎn)體分離面后空化區(qū)域蒸氣相占支配地位，因而兩處的Ma數(shù)均較小；而在空泡界面附近，各相摻混導(dǎo)致局部含氣率接近0.5，故Ma數(shù)最大。

由前面的分析知，回轉(zhuǎn)體繞流流場(chǎng)可能存在頭部脫體激波、分離面處的膨脹波和空泡尾部的壓縮波，圖8給出了典型情況下回轉(zhuǎn)體空泡繞流流場(chǎng)的波系示意圖，其中空泡界面定義為αl=0.1。下面對(duì)回轉(zhuǎn)體泡狀繞流流場(chǎng)的波系產(chǎn)生的原因進(jìn)行分析。

圖8 回轉(zhuǎn)體空泡繞流流場(chǎng)的波系Fig.8 Wave system of flow field around axis-symmetrical body

圖9為含氣率αg分別為0～0.5時(shí)的壓力等值線圖。對(duì)比可以發(fā)現(xiàn)，當(dāng)含氣率增大至一定程度時(shí)（由圖6，αg≈0.3時(shí)，來(lái)流Ma≥1），在回轉(zhuǎn)體頭部會(huì)出現(xiàn)脫體激波，而產(chǎn)生脫體激波后，波后壓力升高，這也是圖6中當(dāng)αg＞0.1后空化長(zhǎng)度減小的原因之一；對(duì)比有無(wú)含氣的情形可以發(fā)現(xiàn)，由于含氣率導(dǎo)致的流場(chǎng)大范圍的可壓縮性，使得空泡尾端無(wú)法閉合在物面上，而是閉合在激波面上，如圖9（b）-（f）所示。

圖9 各含氣率下壓力等值線圖Fig.9 Absolute pressure contour line of each gas fraction

圖10為αg=0.1時(shí)蒸氣相體積分?jǐn)?shù)云圖和Ma數(shù)云圖，其中圖10（b）為便于顯示Ma≥1的區(qū)域，云圖標(biāo)尺下限定為Ma=1?？梢园l(fā)現(xiàn)，雖然含氣率αg=0.1時(shí)，來(lái)流Ma數(shù)小于1，但是由于在蒸汽泡界面附近，空隙率接近0.5，由圖1可知，則空泡界面處馬赫數(shù)大于1為超聲速流，故如圖10（b）所示在流經(jīng)凸鈍角時(shí)產(chǎn)生了扇形膨脹波；空泡界面附近的流場(chǎng)膨脹加速后依舊為超聲速流，而由于空泡外形的影響，空泡尾端相當(dāng)于形成了凹鈍角，超聲速流動(dòng)流經(jīng)該處形成了斜激波。此外，由于繞流流場(chǎng)的壓力分布和含氣率分布不均勻，導(dǎo)致各地聲速的分布也不均勻，最終導(dǎo)致的結(jié)果就是空泡尾端的激波面為曲面，如圖10（b）所示。

圖10 αg=0.1時(shí)蒸氣相體積分?jǐn)?shù)云圖和Ma數(shù)云圖Fig.10 Vapor volume color contour and Mach number color contour of gas fraction 0.1

3.3 不同含氣率對(duì)空化水動(dòng)力的影響

為研究不同含氣率對(duì)回轉(zhuǎn)體繞流流場(chǎng)空化水動(dòng)力的影響，對(duì)含氣率αg=0～0.5工況的壁面壓力系數(shù)進(jìn)行了對(duì)比，如圖11所示，其中壓力系數(shù)的定義為

圖11 各含氣率下壁面壓力系數(shù)曲線Fig.11 Pressure coefficient curve on wall of each gas fraction

圖11表明，隨著含氣率的增大，壓力系數(shù)最小的區(qū)域（即可能產(chǎn)生空化的區(qū)域）逐漸變短，其原因?yàn)殡S著含氣率增大，混合介質(zhì)的可壓縮性增強(qiáng)（特別是頭部脫體激波的產(chǎn)生），導(dǎo)致分離面后的低壓也隨之升高，從而導(dǎo)致低壓力的空化區(qū)域減小。因此，隨著含氣率的增大，流場(chǎng)的空化效應(yīng)被削弱。另外，由圖11也可以看出，隨著流場(chǎng)含氣率的增大，空泡尾端的逆壓梯度逐漸減??；當(dāng)含氣率αg≥0.2時(shí)，逆壓梯度已變得相當(dāng)平緩；這印證了前面的結(jié)論，由于流場(chǎng)含氣帶來(lái)的大范圍可壓縮性的影響，空泡已經(jīng)無(wú)法閉合在物面上，而是閉合在空泡尾端的激波面上，從而體現(xiàn)為空泡尾端壁面逆壓梯度減小。

4 結(jié)論

本文對(duì)回轉(zhuǎn)體的可壓縮泡狀繞流進(jìn)行了數(shù)值計(jì)算，其中含氣率αg=0的計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行了對(duì)比，驗(yàn)證了本文計(jì)算模型的合理性。各含氣率下的計(jì)算結(jié)果表明，隨著含氣率的增大，流場(chǎng)的可壓縮性也相應(yīng)增強(qiáng)，體現(xiàn)為數(shù)值模擬得到了回轉(zhuǎn)體繞流流場(chǎng)的3種波系，包括回轉(zhuǎn)體頭部處的脫體激波、分離面處的膨脹波和空泡尾端的斜激波。而空化與激波的相互作用則體現(xiàn)為：頭激波削弱了空化效應(yīng)，使得空化區(qū)域減??；由于空泡外形的影響，使得空泡尾端出現(xiàn)了斜激波；當(dāng)含氣率超過(guò)一定值時(shí)，空泡已經(jīng)無(wú)法閉合在物面上，而是閉合在空泡尾端的激波面上，體現(xiàn)為空泡尾端壁面逆壓梯度趨于平緩。本文的計(jì)算結(jié)果揭示了可壓縮泡狀流動(dòng)的新的物理現(xiàn)象，有助于實(shí)驗(yàn)結(jié)果的驗(yàn)證。

[1]Karplus H B.The velocity of sound in a liquid containing gas bubbles[J].Armour Research Foundation of Illinois Institute of Technology,1958:C00-248.

[2]Saito Yoshinori,Takami Rieko,Nakamori Ichiro,Ikohagi Toshiaki.Numerical analysis of unsteady behavior of cloud cavitation around a NACA0015 foil[J].Computational Mechanics,2007,40(1):85-96.

[3]Eddington.Investigation of supersonic shock phenomena in a two-phase(liquid-gas)tunnel[J].AIAA,1970,08(01):65-74.

[4]Grant J R,Kirschner I N.High-speed motion in bubbly flows[C]//Fifth International Symposium on Cavitation(CAV2003). Osaka,Japan,2003:1-10.

[5]陳紅兵.定?？蓧嚎s泡狀流動(dòng)激波的數(shù)值仿真[J].船舶力學(xué),2012,12(10):1107-1114. Chen Hongbing.Numerical study on the shocks in thesteady compressible bubbly flows[J].Journal of Ship Mechanics, 2012,12(10):1107-1114.

[6]ANSYS Inc.ANSYS FLUENT 12.0 Theory Guide[K].2009.

[7]褚學(xué)森,王志,顏開(kāi).自然空化流動(dòng)數(shù)值模擬中參數(shù)取值影響的研究[J].船舶力學(xué),2007,11(1):32-39. Chu Xuesen,Wang Zhi,Yan Kai.Parametric study on numerical simulation of natural cavitation flow[J].Journal of Ship Mechanics,2007,11(1):32-39.

[8]權(quán)曉波,李巖,魏海鵬,王寶壽,孔德才.大攻角下軸對(duì)稱航行體空化流動(dòng)特性試驗(yàn)研究[J].水動(dòng)力學(xué)研究與進(jìn)展A輯,2008,23(6):662-667. Quan Xiaobo,Li Yan,Wei Haipeng,Wang Baoshou,Kong Decai.An experiment study on cavitation of underwater vehicle’s surface at large angles of attack[J].Journal of Hydrodynamics,Ser.A,2008,23(6):662-667.

Numerical simulation on interaction between cavitation and shock wave in compressible bubbly flows

DONG Hong-hui,WANG Bao-shou,CHEN Wei-qi,ZHANG Ke
(China Ship Scientific Research Center,Wuxi 214082,China)

The numerical simulation on the cavitating flow around an axis-symmetrical body in the steady compressible bubbly flow was simulated,and the interaction between cavitation and shock wave under different air volume fractions with the constant flow velocity and operating pressure was investigated.Firstly,to verify the used calculation model,the case which air volume fraction takes 0 in the incoming flow was calculated and shows a good agreement with the experimental result.Secondly,more cases with air volume fractions varying from 0 to 0.5 were calculated.The results show that with the increase of air volume fractions, flow field becomes more compressible.This result was confirmed by the appearance of three kinds of wave systems,including detached shock wave before the head,expansion wave at separation surface and oblique shock wave after cavity tail.The interaction between cavitation and shock wave is reflected as follows:bow shock before caviting region weakens cavitation effects,making caviting region decrease;cavity generated by cavitation results in occurrence of oblique shock;cavity closure point cannot locate on the body surface but on the shock surface,which can be confirmed by the adverse pressure gradient at the end of cavity.These results revealed a new physical phenomenon which existed in the interaction process between the cavity and the shock wave in compressible bubbly flows.

bubbly flow;shock wave;cavitation

U661.3

A

10.3969/j.issn.1007-7294.2015.11.001

1007-7294（2015）11-1295-09

2015-08-15

董洪輝（1989-），男，碩士研究生，E-mail：donghonghui@163.com；

王寶壽（1963-），男，研究員。