半球頭模型空化流場非定常特性的數(shù)值模擬

王柏秋，王聰，黃海龍，董磊，張嘉鐘

(1.哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱150001;2.北京宇航系統(tǒng)工程研究所，北京100076;3.哈爾濱工業(yè)大學(xué) 土木學(xué)院，黑龍江 哈爾濱150090;4.中國船舶重工集團公司第703 研究所，黑龍江 哈爾濱150036)

0 引言

空化現(xiàn)象具有強烈的非定常性，數(shù)值方法研究空化流場的特性時，需要針對各自特點，采用相應(yīng)理論下的空化模型。不同的數(shù)學(xué)處理方法衍生出了不同的空化模型:基于空化流場連續(xù)性方程的空化模型，例如Singhal［1］模型，Zwart［2］模型;基于經(jīng)驗得出的空化模型，例如Kunz［3］模型，Merkle［4］模型;通過空泡受力分析得出的空化模型，例如Tamura［5］模型。

以上各種空化模型在表達空化相變率時，各有特點，其中，Singhal 模型是在不可壓縮連續(xù)性方程中導(dǎo)出了氣液兩相間的質(zhì)量交換率，屬于直接方法，但是由于空泡動力學(xué)模型理論不完備，因此，在Singhal 模型的最后階段的處理過程中，引入了簡化的Rayleigh-Plesset 方程;Zwart 空化模型在相變率的推導(dǎo)上與Singhal 模型類似，但是在最后表達空泡界面運動速度時，同樣引入了簡化的Rayleigh-Plesset 方程;Merkle 模型是基于經(jīng)驗，認為空化相變率與流體密度成比例;Kunz 模型繼承了Merkle 模型中的成比例的思想，在變形后的多相流方程中，采用與流體密度成比例的方法得到了空化相變率的經(jīng)驗表達式;Tamura 空化模型是直接對空泡進行受力分析的基礎(chǔ)上得出的。

由以上分析可知，各種空化模型都有各自的數(shù)學(xué)方法上的優(yōu)點，但各種空化模型并未有效區(qū)分空化現(xiàn)象中的相變過程，使得各自模型系數(shù)帶寬差異較大，難于確定與空化數(shù)之間的關(guān)系，因而只能有效捕捉到空化現(xiàn)象中的主空泡，對一些由主空泡內(nèi)氣體泄漏或者其他原因產(chǎn)生的次生空泡無能為力，不能很好地捕捉空化流場的非定常細節(jié)。

本文基于Rayleigh-Plesset 空泡運動方程，在考慮了空泡界面上的相變作用后，導(dǎo)出了一個新的空化模型。針對半球頭模型的空化流場進行了模型系數(shù)的數(shù)值標定，進而成功模擬了空化流場中的主空泡、次生空泡和尾空泡。與實驗現(xiàn)象相比，所得空化流場更加真實，能夠反映空化流場的非定常特性。新空化模型的導(dǎo)出與應(yīng)用，為進一步研究空化機理提供了一定理論依據(jù)。

1 空化相變率

用R 表示發(fā)生相變時的空化相變率，則歐拉法描述有質(zhì)量交換的兩相流體運動的連續(xù)性方程為

式中:α 為水蒸氣體積分數(shù);ρ 表示混合物密度，ρ =αρv+(1 -α)ρl;ρv表示水蒸氣密度;ρl表示液體密度。

將(2)式進行變形，并將(1)式代入(2)式可得

不可壓縮時，混合物密度變化如下

假定空泡以球形形式在液體中均勻分布，則

式中:n 表示液體中空化核數(shù)密度;r 表示球形空泡的半徑。

對(6)式進行全微分可得

將(5)式和(7)式代入(4)式:

由于空泡的收縮與擴張速度極快，因此可認為在空泡變形過程中，泡內(nèi)氣體經(jīng)歷絕熱過程。

當(dāng)空泡擴張時，泡內(nèi)水蒸氣絕熱膨脹，溫度下降，其飽和蒸氣壓快速下降，使泡內(nèi)氣體始終處于過飽和狀態(tài)，因而不斷發(fā)生冷凝相變。此時，(8)式代表冷凝率，其中，混合物密度為相變后水的密度，且由于發(fā)生了相變，因而α= -1.于是，冷凝相變率為

同理，當(dāng)空泡在收縮時，空泡內(nèi)水蒸氣同樣經(jīng)歷絕熱過程，溫度迅速上升，泡內(nèi)飽和蒸氣壓同時快速上升，使泡內(nèi)氣體始終處于欠飽和狀態(tài)，因而，不斷產(chǎn)生蒸發(fā)相變。此時，(8)式代表蒸發(fā)率，其中，混合物密度為相變后的水蒸氣的密度，且由于發(fā)生了相變，因而α= -1.因此，由(8)式可得蒸發(fā)相變率為

以上各式中:Re表示蒸發(fā)相變率;Rc表示凝結(jié)相變率。(9)式、(10)式即為發(fā)生空化相變時的兩相之間的質(zhì)量傳遞率。

2 空化模型

為確定以上兩式中的空泡半徑和空泡泡壁運動速率，考慮不可壓縮流體中球形空泡運動方程，Rayleigh-Plesset 方程［6］:

式中:Σ 表示表面張力系數(shù);νl表示液體運動粘性系數(shù);p∞表示無窮遠壓力;pB表示空泡內(nèi)氣體壓力。(11)式中，用符號頭上一點表示該物理量對時間的一次導(dǎo)數(shù)。

忽略空泡的二階運動影響，不計表面張力和粘性力影響時，由(11)式可得空泡泡壁運動速率為

對于相變率中的空泡半徑，存在眾多不同的處理方式，其中，Zwart 等［2］取rB=1 μm;Martynov 等［7］將空泡半徑等效成空泡長度尺度，并通過空泡核數(shù)密度來表達;Hosangadi 等［8］的可壓縮空化模型和Kunz 等［3］的空化模型中，不曾出現(xiàn)空泡半徑或者空泡尺度，而是以時間尺度的形式來衡量相間質(zhì)量交換;Singhal 等［1］，張瑤等［9］通過引入經(jīng)驗表達式的方式定義了空泡半徑。Singhal 等建議用下列半徑表達式:

式中:We 表示韋伯?dāng)?shù);vrel表示空泡流場中的特征速度。

本文擬采用Singhal 方法表達空泡半徑。綜合(9)式～(13)式可得相變率如下

式中:Ce和Cc為綜合了其他各種未知因素時的模型系數(shù);Ce表示蒸發(fā)系數(shù);Cc表示冷凝系數(shù)。

至此，(16)式和(17)式表達了考慮相變作用時的相變率，即空化模型。

由于水中氣核的含量對空化過程有著重要影響，在標準狀態(tài)下，飽和水中的非凝性氣體含量約為15 ×10-6［10-12］，因此，本文取fg=1.5 ×10-5.

3 新空化模型的應(yīng)用

下面利用(16)式和(17)式對半球頭模型的空化流場進行二維軸對稱數(shù)值計算，流場的運動方程為v∞雷諾時均化N-S 方程，數(shù)值計算模型及網(wǎng)格劃分如圖1所示，其中，模型直徑為D，模型長度為L.

圖1 模型邊界條件及其網(wǎng)格劃分Fig.1 Boundary conditions and mesh

3.1 數(shù)值方法及邊界條件

根據(jù)圖1所示的數(shù)值模型的特點，計算域左側(cè)邊界和兩側(cè)邊界采用流向出口方向的速度入口，其中，左側(cè)流場邊界與模型前端距離為5D，流場兩側(cè)邊界與模型軸線距離為5D;流場右側(cè)邊界為壓力出口，與模型尾部距離為20D;計算域中間為對稱軸。

邊界條件中，由于流場的環(huán)境壓力不變，因此速度入口按下面的空化數(shù)相似準則確定，

式中:σ 表示空化數(shù);u∞表示入口來流速度;pv表示當(dāng)?shù)匾后w溫度下的飽和蒸氣壓。

采用有限體積法對流場進行空間離散。采用PISO 算法對壓力和速度場解耦。采用標準壁面函數(shù)，配合Realizable k-ε 模型模擬湍流項。

3.2 空泡形態(tài)對比分析

為與實際現(xiàn)象進行對比分析，在哈爾濱工業(yè)大學(xué)空泡水洞中進行了空泡實驗，實驗結(jié)果如圖2所示。實驗水洞為閉式循環(huán)水洞，工作段直徑0.2 m，工作壓力為常壓，采用變頻控制，電動機轉(zhuǎn)速變化在±1 r/min，保證了工作段流場穩(wěn)定。

由于空泡的不斷泄氣和再生，使得空化流場具有強烈的非定常性，因此在同一空化數(shù)下，拍攝了不同時刻的自然空泡形態(tài)進行對比。圖2中，模型直徑D=0.04 m，空化數(shù)σ=0.52，高速攝像機拍攝幀率為2 000 fps.

對比圖2中同一空化數(shù)下不同時刻的空泡實驗照片可見，由于非定常性及空泡泄氣和不斷的潰滅再生，空泡水洞實驗中產(chǎn)生的空泡分為3 種類型，分別為位于模型前段的主空泡，如圖2(a);位于模型后段與主空泡緊挨的次生空泡，如圖2(b);位于模型后面的尾空泡，如圖2(c).

圖2 空泡流實驗(σ=0.52)Fig.2 Cavitating flow (σ=0.52)

由實驗過程可知，流場中的次生空泡是由主空泡的不斷泄氣而產(chǎn)生的，因而其空間分布不規(guī)則，是空化流場具有強烈非定常性的重要來源之一，同時，次生空泡的存在給數(shù)值模擬工作帶來了難度。

圖3為半球頭模型空化流場的定常狀態(tài)數(shù)值模擬結(jié)果，其中，圖3(a)為新空化模型計算結(jié)果，圖3(b)為Singhal 模型在默認模型系數(shù)下的計算結(jié)果，圖中水蒸氣體積分數(shù)取值范圍0～1.

與圖2中的實驗現(xiàn)象進行對比可知，圖3(a)中，應(yīng)用新空化模型所得到的計算結(jié)果能夠同時得到主空泡、次生空泡和尾空泡，與實驗現(xiàn)象更為接近，更能真實反映空泡流場的細節(jié)。而圖3(b)所示的Singhal 模型的計算結(jié)果中，次生空泡幾乎不出現(xiàn)，且尾空泡流場也較弱，與實驗現(xiàn)象存在明顯差異。

另外，在進行減壓空泡實驗時，水洞擴張段內(nèi)會不斷傳出強烈的爆裂噪聲，這正是由于減壓條件增強了流場內(nèi)的次生空泡的輸運，使大量的次生空泡延遲至水洞擴張段內(nèi)才發(fā)生不斷潰滅。而常壓和增壓實驗時，次生空泡的輸運受到抑制而提前潰滅，因而，水洞擴張段內(nèi)幾乎無爆裂噪聲。

圖3 不同空化數(shù)下的空泡形態(tài)對比(從上至下，空化數(shù)σ 分別為0.5、0.4、0.3、0.2)Fig.3 Cavity shape at different cavitation number (from top to down:σ=0.5，σ=0.4，σ=0.3，σ=0.2)

3.3 非定常計算結(jié)果分析

實驗數(shù)據(jù)一般反應(yīng)的是試驗點上的某物理量在一個時間片段內(nèi)的平均值，因此，為定量驗證新空化模型的正確性和有效性，有必要取非定常數(shù)值計算結(jié)果的平均值與實驗結(jié)果進行對比分析。

下面以圖1所示數(shù)值模型為例進行非定常數(shù)值計算，模型尺寸和數(shù)值實驗點的分布情況如圖4所示，其中實驗點的布置參照了文獻［1］和文獻［8］.

圖4 數(shù)值實驗點Fig.4 Numerical experiment points

由于次生空化區(qū)內(nèi)的模型表面壓力波動明顯，具有特征意義，因此，圖5對次生空化區(qū)內(nèi)任意一點的表面壓力系數(shù)隨時間的波動情況進行了研究。圖6是圖5所示的數(shù)值計算結(jié)果的時均值與文獻［8］的對比。其中，L/D 代表用模型直徑D 無量綱化的模型沿程長度，表面壓力系數(shù)的定義為

圖5和圖6的數(shù)值計算結(jié)果表明，應(yīng)用新空化模型所得到的表面壓力系數(shù)的時均值與文獻［8］中的實驗結(jié)果非常接近，證明了新空化模型對空化流場的定量計算是正確有效的。

3.4 定常計算結(jié)果分析

圖7與圖8分別為模型表面壓力系數(shù)分布的定常計算結(jié)果與文獻［8］的對比情況(σ=0.4 和σ=0.2).

在圖7與圖8中，表面壓力系數(shù)曲線下降時，流體局部壓力降低，空泡開始生成和擴張，表面壓力系數(shù)曲線上升時，流體局部壓力升高，空泡開始收縮潰滅。

圖5 P4 表面壓力系數(shù)的非定常波動Fig.5 Surface pressure coefficient fluctuations of P4

圖6 表面壓力系數(shù)分布(σ=0.4)Fig.6 Distribution of surface pressure coefficient

由圖7中表面壓力系數(shù)曲線的對比可知，新空化模型計算得到的流場中，由于次生空泡的存在，使得模型表面壓力系數(shù)曲線沿程波動，符合實驗現(xiàn)象。而文獻［1］的Singhal 模型計算結(jié)果中，空泡的生成與潰滅只出現(xiàn)了一次，沒有捕捉到次生空泡，因而表面壓力系數(shù)曲線平滑，沒有非定常特征。

3.5 模型系數(shù)的數(shù)值標定

由(9)式與(10)式的推導(dǎo)及分析可知，在空泡的擴張過程中，空泡界面上發(fā)生的是冷凝現(xiàn)象，冷凝相變占優(yōu)，而空泡在收縮時，界面上發(fā)生了蒸發(fā)現(xiàn)象，蒸發(fā)相變占優(yōu)。因此，空化模型中的冷凝系數(shù)同時表達了空泡的擴張?zhí)匦?，蒸發(fā)系數(shù)同時又表達了空泡的收縮特性。由于模型中冷凝系數(shù)和蒸發(fā)系數(shù)互不相關(guān)，因此，在數(shù)值方法中，采用逐步逼近的方法對模型系數(shù)進行數(shù)值標定。通過數(shù)值模擬與文獻［1，8］中的表面壓力系數(shù)曲線上的主空泡位置進行比對后，針對半球頭模型的空化流場，得到了如圖9所示的模型系數(shù)與空化數(shù)之間的關(guān)系曲線。

圖7 模型表面壓力系數(shù)分布及相應(yīng)空泡流場(σ=0.4)Fig.7 Surface pressure coefficient and cavity volume fraction (σ=0.4)

由圖9(a)可見，模型系數(shù)與空化數(shù)關(guān)系密切，在當(dāng)前空化數(shù)范圍內(nèi)(0.1 ＜σ ＜0.5)，隨著空化數(shù)的降低，蒸發(fā)系數(shù)先減小后增加，而冷凝系數(shù)維持在較低水平，幾乎保持不變。模型系數(shù)隨空化數(shù)變化而變化的重要原因之一是模型自身理論不完備，例如由于缺乏空泡半徑的解析表達式，因而不得不采用與以往模型類似的處理方式，即使用近似表達式來表達模型中的空泡半徑。因此，在使用空化模型時，首先應(yīng)確定對應(yīng)空化數(shù)下的模型系數(shù)。

圖10所示為使用了不適當(dāng)?shù)哪Ｐ拖禂?shù)的計算結(jié)果。由圖可見，當(dāng)模型系數(shù)過大或過小時，空化流場中的相變作用增強或變?nèi)?，引起計算結(jié)果嚴重偏離實驗，甚至可能導(dǎo)致流場中的次生空化消失或發(fā)生嚴重變形。

3.6 超空化和弱空化

圖8 模型表面壓力系數(shù)分布及相應(yīng)的空泡流場(σ=0.2)Fig.8 Surface pressure coefficient and cavity volume fraction (σ=0.2)

以上討論的空化流場尚未超空化流場或者弱空化流場，例如當(dāng)σ≤0.1 或σ ＞0.5.

為進一步考察模型系數(shù)與空化數(shù)之間的關(guān)系，采用同樣的方法，對超空化流場和弱空化流場分別進行了數(shù)值模擬，計算結(jié)果如圖11、圖12 所示。

與圖9的非超空化流場相比，圖11 中超空化流場中的模型系數(shù)較大。由于超空泡的形成，原本附著于模型表面的次生空泡完全脫落，并與尾空泡結(jié)合形成新的尾空泡。

同時，由圖11 可見，不適當(dāng)?shù)哪Ｐ拖禂?shù)依然會導(dǎo)致對流場特征的捕捉失敗。例如，較小的模型系數(shù)將導(dǎo)致流場中的尾空泡區(qū)直接消失，如圖11(a)所示。較大的模型系數(shù)將導(dǎo)致流場發(fā)生紊亂，如圖11(c)所示。

相對于超空化流場和次生空化流場，弱空化流場的模型系數(shù)與空化數(shù)之間的關(guān)系較為簡單，如圖12所示。圖12 的弱空化流場中，沒有次生空泡，甚至沒有明顯的空化產(chǎn)生，因而，模型系數(shù)不再依賴于空化數(shù)，可以自由選取。

圖9 模型系數(shù)與空化數(shù)之間的關(guān)系Fig.9 Relationship of model coefficients and cavitation number

圖10 不適當(dāng)模型系數(shù)下的空化流場(σ=0.4)Fig.10 Cavitation flow field under improper model coefficients

至此，按照空化數(shù)的變化，空化流場可以分為3 種:超空化流，次數(shù)空化流，弱空化流。圖9、圖11和圖12 的計算結(jié)果表明，3 種空化流場中的模型系數(shù)與空化數(shù)之間的關(guān)系有較大差別，在使用空化模型之前，應(yīng)根據(jù)實驗或文獻資料首先對模型系數(shù)進行數(shù)值標定。

4 結(jié)論

圖11 不同模型系數(shù)下的超空泡流場(σ=0.1)Fig.11 Numerical results of supercavitation flow field

圖12 弱空化流場的數(shù)值模擬(σ ＞0.5)Fig.12 Numerical results of weakcavitation flow field

通過引入空泡界面上的相變作用，分析了空化流場中的相變過程。結(jié)合簡化的Rayleigh-Plesset 方程，給出了一個新的空化模型。將新空化模型應(yīng)用于半球頭水下航行體的空化問題，得到了以下結(jié)論:

1)新空化模型的成功應(yīng)用和驗證表明對相變過程的分析是有效的，即冷凝發(fā)生于空泡擴張時，而蒸發(fā)發(fā)生于空泡收縮時。

2)空化流場中存在3 種典型的空泡結(jié)構(gòu)，分別為主空泡、次生空泡、尾空泡。按照這3 種典型結(jié)構(gòu)的空泡在空化流場中出現(xiàn)的可能性，空化流場亦可分為3 種類型，分別為超空化流場、次生空化流場、弱空化流場。

3)由于考慮了空泡界面上的相變作用，應(yīng)用新的空化模型能夠模擬主空泡、次生空泡和尾空泡，得到的計算結(jié)果更加接近實驗現(xiàn)象。

4)模型中的相變系數(shù)不是常數(shù)，而是隨著空化數(shù)改變而改變。在次生空化區(qū)內(nèi)，相變系數(shù)隨著空化數(shù)的減小，先減小后增加，并且，蒸發(fā)系數(shù)遠大于冷凝系數(shù)。

5)空泡擴張過程中，發(fā)生冷凝相變，空泡的擴張系數(shù)又是相變中的冷凝系數(shù);空泡收縮過程中，發(fā)生蒸發(fā)相變，空泡的收縮系數(shù)又是相變中的蒸發(fā)系數(shù)。

新空化模型物理意義明確，適用于理論分析和數(shù)值計算，為進一步研究空化機理提供了理論基礎(chǔ)。

References)

［1］ Singhal A K，Athavale M M，LI Hui-ying，et al.Mathematical basis and validation of the full cavitation model［J］.Journal of Fluids Engineering，2002，124(3):617 -624.

［2］ Zwart P J，Gerber A G.Belamri T.A two-phase flow model for predicting cavitation dynamics［C］∥5th International Conference on Multiphase Flow.Yokohama:ICMF，2004:152.

［3］ Kunz R F，Boger D A，Stinebring D R，et al.A preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction［J］.Computers ＆ Fluids，2000，29(8):849-875.

［4］ Merkle C L，F(xiàn)eng J，Buelow P E O.Computational modeling of the dynamics of sheet cavitation［C］∥Proceedings of Third International Symposium on Cavitation.Grenoble:［s.n.］，1998:307-313.

［5］ Tamura Y，Matsumoto Y.Improvement of bubble model for cavitating flow simulations［J］.Journal of Hydrodynamics，2009，21(1):41 -46.

［6］ Brennen C E.Cavitation and Bubble Dynamics［M］.Oxford:Oxford University Press，1995:47 -50.

［7］ Martynov S B，Mason D J，Heikal M R.Numerical simulation of cavitation flows based on their hydrodynamic similarity［J］.International Journal of Engineering Research，2006，7(3):283-296.

［8］ Hosangadi A，Ahuja V，Arunajatesan S.A generalized compressible cavitation model ［C］∥Fourth International Symposium on Cavitation.Pasadena:California Institute of Technology，2001:sessionB4.003.

［9］ ZHANG Yao，LUO Xian-wu，JI Bin，et al.A thermodynamic cavitation model for cavitating flow simulation in a wide range of water temperature ［J］.Chinese Physics Letters，2010，27 (1):016401.

［10］ LI Zhi-min，PENG Xiao-feng，LEE Du-zhong.Interfacial mass transfer around a vapor bubble during nucleate boiling［J］.Heat Mass Transfer，2004，41(1):5 -11.

［11］ 褚學(xué)森，王志，顏開.自然空化流動數(shù)值模擬中參數(shù)取值影響的研究［J］.船舶力學(xué)，2007，11(1):32 -39.CHU Xue-sen，WANG Zhi，YAN Kai.Parametric study on numerical simulation of natural cavitation flow［J］.Journal of Ship Mechanics，2007，11(1):32 -39.(in Chinese)

［12］ 魏海鵬，郭鳳美，權(quán)曉波.潛射導(dǎo)彈表面空化特性研究［J］.宇航學(xué)報.2007，28(6):1506 -1523.WEI Hai-peng，GUO Feng-mei，QUAN Xiao-bo.Research on cavitation of submarine launched missile’s surface［J］.Journal of Astronautics，2007，28(6):1506 -1523.(in Chinese)