基于Rayleigh-Plesset方程的空化模型改進與應用

洪 鋒 高振軍 袁建平

(1.三峽大學機械與動力學院， 宜昌 443002； 2.三峽大學水電機械設備設計與維護湖北省重點實驗室， 宜昌 443002；3.江蘇大學國家水泵工程技術研究中心， 鎮江 212013)

0 引言

空化是一種包含相變、介質壓縮性、粘性及湍動效應的復雜多相流。在水力機械領域，空化的發生總會造成不同程度的負面作用。迄今為止，空化和空泡動力學特性一直是國際水動力學界的熱點和前沿科學問題[1]。

當前，應用于工程問題中的空化數值模型大部分是由只考慮空泡半徑一階導數項及壓力驅動項的Rayleigh-Plesset(R-P)方程推導而來[2-8]，該類模型沒有考慮二階導數項、粘性項及表面張力項對空泡半徑增長的影響。文獻[9]指出，完整R-P方程中的二階導數項對超聲誘導空泡半徑的增長有著重要的影響，且不可被忽略；基于簡化R-P方程的空化模型在模擬繞水翼及螺旋槳空化流場時均存在不足之處[10-11]。上述3種空化模型在推導過程中只考慮了兩種組分(空泡相和液相)，而實際流體中還包含了不可凝結的氣體(NCG)，一個標準大氣壓下純水所能溶解NCG的飽和質量濃度為1.5×10-5mg/L。該部分氣體在空化發生過程中充當空化核，空化核對空化初生有著重要的作用[12]。

為了改善這種基于簡化R-P方程的空化模型模擬空化流動的能力，本文基于一種改進的R-P球形空泡動力學模型、液相-空泡相-不可凝結氣體相3種組分的連續性方程和均相流假設推導一種修正的R-P模型；然后運用UDF(User-define function)技術，把該模型嵌入到ANSYS FLUENT 14.5平臺，并聯立一種改進的濾波器湍流模型(MFBM)對繞二維Clark-Y水翼的非定常空化流場進行數值計算。通過與Schnerr-Sauer(S-S)空化模型的數值計算結果及文獻中公開的實驗結果作對比，以驗證修正的R-P模型在預測非定常空化流動特征的可行性及準確性。

1 基于R-P方程的改進空化模型

廣義的空泡動力學方程為[12]

(1)

式中R——空泡半徑νl——液相運動粘度

S——表面張力系數

pB(t) ——空泡內壓力

p∞(t)——無窮遠處壓力

ρl——液相密度t——時間

等號左邊依次為二階導數項、一階導數項、粘性項和表面張力項，等號右邊為壓力驅動項。式(1)中二階導數項在空泡體積急劇變化時發揮重要作用，而粘性項與表面張力項對空泡半徑變化的影響尚不明確，因此研究兩種簡化形式的空泡動力學模型

(2)

(3)

利用Matlab軟件，采用4～5階RUNGE-KUTTA算法分別對式(2)和式(3)進行求解。假定初始半徑為R0的空化核在時間為0時進入低壓區，可以得到該空化核流過低壓區時R(t)的一個典型解，如圖1(圖中T表示球形空泡生長周期)所示。從圖中可以看到，在t=10T之后，粘性項對空泡半徑變化的影響較小，且相比于表面張力項，粘性的影響是可以忽略不計的，這主要是因為液態水及因空化形成的空泡均屬于低粘度工質，其運動粘度系數為10-6數量級。

圖1 不同空泡動力學模型計算得到的空泡增長過程Fig.1 Variation of bubble radium based on different bubble dynamic models

因此，考慮由壓力驅動項、空泡一階導數項、二階導數項及表面張力項所組成的動力學模型，如式(3)所示。將式(3)左右兩邊同時乘以2R2dR/dt，應用初始條件DR/Dt(0)=0[12]積分可以得到

(4)

因此，改進后的球形空泡動力學方程為

(5)

當i=1時表示凝結過程；當i=2時表示蒸發過程。

為了真實反映空化流動的本質，將介質處理成空泡相、液相及不可凝結氣體的混合物，但是空化引起的質量傳輸過程只在液相和空泡相之間產生，且液相、空泡相及不可凝結氣體相3種組分的連續性方程分別為

(6)

(7)

(8)

其中

αl+αv+αg=1

(9)

式中ρ——混合密度V——速度

Se——蒸發率Sc——凝結率

α——體積分數

下標l、v和g分別代表液相、空泡相和不可凝結氣體相。

均相流思想是把介質假想成一種均勻的液相-空泡相-不可凝結氣體的混合物，這種處理方式可以避免因求解每一組分獨立的連續性方程及其動量方程所需要的巨大數值耗散時間。同時，假定各組分之間無相對滑移速度，即

Vl=Vv=Vg=V

(10)

式(6)～(8)相加，可以得到混合相的連續性方程，即

(11)

混合密度ρ及混合動力粘度系數μ分別表示為

ρ=αlρl+αvρv+αgρg

(12)

μ=αlμl+αvμv+αgμg

(13)

在單元控制體內，為了解決液相與空泡相之間由于發生相變而產生劇烈的相對密度變化引起的數值計算困難，故使用一種非守恒型連續性方程來描述混合流體的運動。該方程為[13]

(14)

因為

(15)

由式(10)及式(15)可得

(16)

由式(9)可得

(17)

把式(16)代入到式(17)可得

(18)

聯立式(12)、(14)及(18)可得

(19)

把式(19)代入到式(7)中可得

(20)

空泡體積分數αv由空泡半徑R及單位體積內空泡數量N0表示[7]，即

(21)

由式(21)可推導出

(22)

將式(22)代入到式(20)可得

(23)

把式(5)代入到式(23)中可得最終形式的蒸發率表達式為

(24)

同理，最終形式的凝結率表達式為

(25)

N0=1×1013m-3、R0=1 μm，NCG的體積分數計算公式為

αg=fgρ/ρg

(26)

ρg按照密度形式理想氣體狀態方程計算，質量分數fg=1.5×10-6。

計算時為了考慮湍流脈動對空化初生的影響，將空化壓力pv修正為

pv=psat+0.195ρk

(27)

式中psat——飽和蒸汽壓，取3 169 Pa

k——湍動能

式(24)～(27)所組成的數學模型即為修正的R-P模型。該模型與Schnerr-Sauer(S-S)模型均是基于均相流理論建立，也都采用式(20)來描述空泡體積分數與空泡半徑及單位體積內空化核數量之間的關系，但是修正的R-P模型具有以下優點：①修正的R-P模型把空泡流真實地還原成液相、空泡相及不可凝結氣體的混合物，而S-S模型則把介質處理成由液相和空泡相組成的混合物。②修正的R-P模型推導是基于改進后的Rayleigh-Plesset空泡動力學方程，該方程考慮了空泡半徑的二階變化及表面張力的影響。③修正的R-P模型考慮了湍流脈動對空化初生的影響。

2 修正的濾波器湍流模型(MFBM)

湍流模型的合理選取對空化流動數值計算尤其重要，傳統的雙方程湍流模型會過度預測閉合空穴尾部的湍流粘度，從而會阻止回射流向水翼前緣運動，使得空穴脫落不會被預測到，而直接模擬或者大渦模擬會消耗大量的計算機資源，尤其在實際的水力機械計算中這一問題會顯得更加突出[14]。因此，本文采用文獻[15]中提出的一種修正的濾波器湍流模型(MFBM)。MFBM實際上是一種修正的RNGk-ε的混合模型，模型的思想是通過一個函數fMFBM聯立濾波器模型(FBM)和密度修正模型(DCM)來共同修正傳統RANS模型的湍流粘度，從而使得修正后的模型兼備FBM和DCM的優點。基于RNGk-ε模型修正后的湍流粘度系數為

(28)

其中

fMFBM=min(fDCM,fFBM)

(29)

(30)

(31)

式中λ——濾波尺寸

n=10，Cμ=0.09，C3=1.0。

3 計算網格及邊界條件

以二維Clark-Y型水翼作為研究對象，翼型弦長c=70 mm，攻角A=8°。計算域為矩形，其寬度和長度分別為2.7c和10c，水翼前緣到計算域進口距離為3c，水翼尾緣到計算域出口距離為6c，該水翼的空化水洞實驗由王國玉等[16-19]完成。邊界條件設定均與文獻[16]保持一致，即入口采用速度入口，計算域進口速度Uin=10 m/s，出口采用壓力出口，通過調整出口壓力，從而獲取相應的空化數條件。計算域上下邊界均為自由無滑移壁面條件，為獲得更精確的計算結果，翼型壁面函數選擇增強型壁面條件。

對計算域劃分了4種不同尺寸的網格，并作網格無關性分析來消除網格數量造成的計算誤差。不同網格方案對應的網格節點數分別為N1=59 576，N2=106 523，N3=157 296及N4=198 776。通過計算空化數σ=0.8條件下水翼上下表面壓力系數Cp隨網格節點數的變化規律來確定最終的計算網格。計算結果如圖2所示，從圖中可以看出，不同網格方案計算得到的翼型上表面壓力分布隨著網格節點數的增加逐漸趨于一致。細網格與極細網格之間的壓力分布幾乎相等，表明繼續增大網格數量對計算結果影響可以忽略不計，因此最終選擇細網格作為本文的計算網格。

圖2 網格無關性分析Fig.2 Mesh independence analysis

4 計算結果分析

為了研究修正的R-P模型對水翼空化流場的模擬效果，采用數值計算的方法得到了Schnerr-Sauer模型及修正的R-P模型對片狀空化(σ=1.4)及云狀空化(σ=0.8)兩種典型空化條件下翼型所受的水動力系數及其吸力面上的空穴形態分布規律，并與文獻中的實驗結果作對比。數值計算均采用非定常求解器，時間步長Δt依據CFL條件(Courant-Friedrichs-Lewy condition)確定[1]，結合計算成本，最終取Δt=1.0×10-6s。定義空化數、升力系數、阻力系數及空泡總體積分別為

(32)

(33)

(34)

(35)

式中pout——出口壓力N——網格單元總數

Fl——升力Fd——阻力

L——翼型展向長度

αi——網格單元內的空泡體積分數

Vi——網格單元的體積

4.1 準靜態片狀空化

張博等[20]對Clark-Y型水翼片狀空化進行了大量實驗研究。根據實驗結果，當空化數降低至σ=1.4時，片狀空化逐漸形成，此時翼型前緣附著透明的空泡，而翼型尾緣處存在波動，且有較小的空泡團脫落。

圖3為兩種模型模擬得到的4個不同時刻空泡體積分數分布圖，從圖中可以看到，在翼型前緣處，S-S模型得到的高質量分數的空化區逐漸減小，表明該處結構穩定的片狀空穴長度發生了波動，這與實驗結果不符。相比較而言，采用修正的R-P模型計算得到的片狀空穴長度基本上維持不變，而且在空穴尾部出現了明顯的相間破碎界面并伴有體積分數較低的空泡團脫落現象，而這與實驗結果的吻合度較好。同時，不同空化模型計算得到的平均升力、阻力系數對比如表1所示，從表中可以看到，相比于S-S模型，修正的R-P模型模擬得到的平均升阻力系數與實驗值吻合度更好。綜上分析，修正的R-P模型對片狀空穴長度、空穴閉合末端流動細節的捕捉能力以及對水動力系數預測的能力均優于S-S模型。

圖3 兩種不同空化模型計算得到的片狀空穴形態對比Fig.3 Comparisons of evolution of sheet cavitation predicted by Schnerr-Sauer and proposed models

模型類型 Cl Cd計算值實驗值[18]相對誤差/%計算值實驗值[18]相對誤差/%修正的R-P模型1.0511.169.390.0390.0414.88Schnerr-Sauer模型1.0341.1610.860.0370.0419.76

4.2 周期性云狀空化

依據文獻[21]，當空化數降低至σ=0.8時，翼型前緣處附著的片狀空穴所表現出的不穩定性加劇，其長度和厚度都明顯增加，且尾緣處大體積空泡團脫落現象也更加顯著。

圖4為兩種模型計算得到的瞬時升力系數與實驗值[21]的對比。從圖中可以看到，修正的R-P模型模擬得到的瞬時升力系數與實驗值的變化趨勢吻合程度更好，而S-S模型模擬得到的升力系數出現了較多的波動。表2所示為不同空化模型模擬得到的平均升力、阻力系數與實驗值對比，從表中可以看到，修正的R-P模型模擬得到的平均升力系數與實驗值吻合度更好，其相對誤差僅為6.93%，且修正的R-P模型預測到的平均阻力系數相比于Schnerr-Sauer模型預測結果也與實驗值更接近。Schnerr-Sauer模型模擬得到的平均升、阻力系數與實驗值的相對誤差均超過10%。

圖5所示為兩種空化模型計算得到的空泡總體積隨時間變化圖，從圖中可以看到，二者計算得到的空泡體積隨時間呈較為明顯的周期性變化，且修正的R-P模型和S-S模型計算得到的云狀空化平均演變周期分別為T1=45.5 ms和T2=38.4 ms，實驗觀測到的云狀空化演變周期為T3=50 ms[18]。由此可見，修正的R-P模型在模擬云狀空化周期性變化過程時相比于S-S模型更具有優越性；另一方面，修正的R-P模型模擬得到的空泡體積明顯大于S-S模型，表明非線性模型可以計算得到面積更大的空化區域，這主要是由于修正的R-P模型考慮了NCG、湍流脈動對空化初生的影響。

圖6(圖中t0表示翼型尾緣處空泡團體積增長至最大的時刻)所示為模擬得到的一個周期內云狀空化演變與實驗結果[21]的對比，圖6a～6e為實驗結果，圖6f～6j為修正R-P模型的預測結果，圖6k～6o為S-S模型的預測結果。對比可以發現，不同空化模型模擬得到的空泡演變過程都基本上反映出了云狀空化中位于翼型前緣處片狀空穴長度和厚度的增長，以及隨后云狀空化脫落和向下游潰滅的周期性變化過程。但是，二者的不同之處主要有以下4點：總體上來看，不同時刻修正的R-P模型模擬得到的空化區域均大于S-S模型的模擬結果；在t=t0時刻，片狀空穴回縮至翼型前緣幾乎消失，尾緣處空泡形成團狀結構并開始逐漸脫離壁面，對比發現，修正的R-P模型預測到的空泡團尺寸與實驗拍攝到的結果更加吻合，且空泡團位置仍靠近翼型尾緣，而S-S模型預測到的空泡團已經離開翼型吸力面；在t=t0+58%Ti(i=1、2、3)時刻，空穴內部表現出強烈的不穩定性且兩相界面出現明顯的破碎，對比可以發現S-S模型預測的界面破碎程度要高于修正的R-P模型的模擬結果，而空穴內強烈的不穩定性與回射流結構有著緊密聯系[1]，表明S-S模型預測到的回射流強度大于修正的R-P模型的預測結果；在t=t0+84%Ti(i=1、2、3)時刻，實驗結果中的空穴長度近似為0.5倍弦長，S-S模型模擬的空穴長度明顯偏小。綜上分析，盡管兩個空化模型均可以模擬出云狀空化周期性演變過程，但在一些空化流動細節的捕捉上，修正的R-P模型更具有優越性。

圖4 兩種模型模擬得到的升力系數與實驗值對比Fig.4 Comparisons of transient lift coefficients of simulation and experiment

模型類型ClCd計算值實驗值[18]相對誤差/%計算值實驗值[18]相對誤差/%修正的R-P模型0.7520.8086.930.1280.11511.30Schnerr-Sauer模型0.7230.80810.520.1310.11513.91

圖5 兩種模型模擬得到的空泡總體積對比Fig.5 Comparisons of total vapor volume of simulations by using two models

圖6 兩種模型模擬得到的云空化演變與實驗值對比Fig.6 Comparisons of predicted evolution of cloud cavitation of simulations and experiment by two models

5 結論

(1)在準靜態片狀空化階段，修正的R-P模型計算得到的片狀空穴長度基本上維持不變，而且在空穴尾部出現了明顯的相間破碎界面并伴有體積分數較低的空泡團脫落現象，與實驗描述較為一致。

(2)修正的R-P模型與Schnerr-Sauer模型預測云空穴演變平均周期分別為45.5 ms和38.4 ms，且修正的R-P模型能較清晰地預測云狀空泡的初生、發展、潰滅和脫落的全過程。

(3)云空化階段，修正的R-P模型模擬得到的瞬時升力系數與實驗值的變化趨勢吻合程度更高，盡管兩個空化模型均可以模擬出云狀空化周期性演變過程，但在一些流動細節的捕捉上，修正的R-P模型對云空化周期性演變的預測能力更好。

(4)修正的R-P模型聯立改進的濾波器湍流模型，聯合求解繞水翼非定常空化流動具有較高的準確性和可行性。

1 洪鋒，袁建平，周幫倫.空泡半徑非線性變化的空化模型及應用[J]．華中科技大學學報，2015，43(10)：15-20.

HONG Feng，YUAN Jianping，ZHOU Banglun. Cavitation model considering non-linear variation of bubble radius and its application [J]. Journal of Huazhong University of Science and Technology，2015，43(10)：15-20. (in Chinese)

2 SINGHAL A K，ATHAVALE M M，LI H，et al．Mathematical basis and validation of the full cavitation model[J]. ASME Journal of Fluids Engineering，2002，124(3)：617-624．

3 RHEE S H, KAWAMURA T, LI H．Propeller cavitation study using an unstructured grid based Navier-Stoker solver[J]．ASME Journal of Fluids Engineering，2005，127(5)：986-994．

4 ATHAVALE M M, LI H Y, JIANG Y U, et al. Application of the full cavitation model to pumps and inducers[J]. International Journal of Rotating Machinery, 2002, 8(1):45-56．

5 ZWART P J，GERBER A G，BELAMRI T．A two-phase flow model for predicting cavitation dynamics[C]∥ICMF Fifth International Conference on Multiphase Flow，2004．

6 LIU H, WANG J, WANG Y, et al．Influence of the empirical coefficients of cavitation model on predicting cavitating flow in the centrifugal pump[J]．International Journal of Naval Architecture and Ocean Engineering, 2014, 6(1): 119-131．

7 SCHNERR G H，SAUER J．Physical and numerical modeling of unsteady cavitation dynamics[C]∥ICMF Fourth International Conference on Multiphase Flow，2001．

8 JI B，LUO X W，ARNDT R E A，et al．Large eddy simulation and theoretical investigations of the transient cavitating vortical flow structure around a NACA66 hydrofoil[J]．International Journal of Multiphase Flow，2015，68：121-134．

10 MORGUT M，NOBILE E，BILUI．Comparison of mass transfer models for the numerical prediction of sheet cavitation around a hydrofoil[J]．International Journal of Multiphase Flow，2011，37(6)：620-626．

11 楊瓊方，王永生，張志宏. 改進空化模型和修正湍流模型在螺旋槳空化模擬中的評估分析[J]. 機械工程學報，2012，48(9)：178-185.

YANG Qiongfang, WANG Yongsheng, ZHANG Zhihong. Assessment of the improved cavitation model and modified turbulence model for ship propeller cavitation simulation[J]. Journal of Mechanical Engineering, 2012, 48(9): 178-185. (in Chinese)

12 BRENNEN C E. Cavitation and bubble dynamics [M]. Cambridge: Cambridge University Press, 2013.

13 YUAN W，SCHNERR G H．Numerical simulation of two-phase flow in injection nozzles： interaction of cavitation and external jet formation[J]．ASME Journal of Fluids Engineering，2003，125(6)：963-969．

14 COUTIER-DELGOSHA O，FORTES-PATELLA R，REBOUD J L，et al．Numerical simulation of cavitating flow in 2D and 3D inducer geometries[J]．International Journal for Numerical Methods in Fluids，2005，48(2)：135-167．

15 LI Z, ZHENG D, HONG F, et al. Numerical simulation of the sheet/cloud cavitation around a two-dimensional hydrofoil using a modified URANS approach [J]. Journal of Mechanical Science & Technology, 2017, 31(1):215-224.

16 HUANG B，YOUNG Y L，WANG G Y，et al．Combined experimental and computational investigation of unsteady structure of sheet/cloud cavitation[J]．ASME Journal of Fluids Engineering，2013，135(7)：071301．

17 黃彪，王國玉，王復峰，等．繞水翼非定常空化流場的實驗研究[J]．兵工學報，2012，33(4)：401-407．

HUANG Biao, WANG Guoyu, WANG Fufeng, et al. Experrimental study on unsteady cavitating flow around hydrofoil with particle imaging velocimetry system[J]. Acta Armamentarii, 2012,33(4):401-407. (in Chinese)

18 WANG G Y，SENOCAK I，SHYY W，et al．Dynamics of attached turbulent cavitating flows[J]．Progress in Aerospace Sciences，2001，37(6)：551-581．

19 HUANG B，WANG G Y，YU Z，et al．Detached-eddy simulation for time-dependent turbulent cavitating flows[J]．Chinese Journal of Mechanical Engineering，2012，25(3)：484-490．

20 張博，王國玉，李向賓，等.繞水翼片狀空化流動結構的數值與實驗研究[J]．工程熱物理學報，2008，29(11)：1847-1851．

ZHANG Bo, WANG Guoyu, LI Xiangbin, et al. Analysis of sheet cavitation flow structure by numerical and experimental studies [J]. Journal of Engineering Thermophsics, 2008, 29 (11):1847-1851. (in Chinese)

21 張光建．水翼及軸流泵非定常空化特性數值研究[D]．鎮江：江蘇大學，2014．

ZHANG Guangjian. Numerical simulation of the unsteady cavitation characteristics in hydrofoil and axial flow pump [D].Zhenjiang: Jiangsu University, 2014. (in Chinese)