空化模型中的相變系數對不同頭型的適用性研究

王 超，林 揚，胡志強，衣瑞文

(1. 中國科學院沈陽自動化研究所 機器人學國家重點實驗室，遼寧 沈陽 110016；2. 中國科學院大學，北京 100049)

0 引 言

采用超空泡減阻技術的水下航行體，其減阻率可達95%以上，因此該項技術獲得廣泛應用研究[1]。為了描述自然空化過程中汽-液之間的相變過程，目前發展出了多種空化模型，這其中包括基于Rayleigh-Plesset氣泡動力學方程的Singhal模型[2]、Zwart模型[3]，基于經驗的Kunz模型[4]、Hosangadi模型[5]等。不同的空化模型基本都包含可人工調節的相變系數，對于不同的空化模型、不同類型的空化問題，相變系數對數值計算結果影響很大。

研究人員針對半球頭圓柱研究了Singhal空化模型中不同相變系數對計算結果的影響，指出不同空化數下相變系數的設置規范[6–7]。Liu等[8]在研究水泵葉輪的空化問題后認為，Zwart空化模型中的蒸發系數和冷凝系數的取值對葉輪上的空泡形態和尺度有重要影響。事實上，不同頭型（平頭、球頭、錐頭等）會導致周圍流場產生不同的脈動特性和漩渦特性，致使初生空泡形態和空化產生的難易程度不盡相同[9]。因此在利用空化模型模擬不同頭型的空化現象時，相變系數也不是一成不變的。

本文采用基于CFD的數值計算方法，研究利用Zwart空化模型模擬不同頭型的空化現象時相變系數對數值計算結果的影響，為進一步提高空化問題的數值計算精度提供指導性依據。

1 數值計算理論

1.1 N-S方程

針對汽-液兩相流，混合介質的N-S方程可以表示為：

1.2 湍流模型

耗散率ε方程為：

式（2）與式（3）中：

表 1 模型中的常數Tab. 1 Constants in model

表 1 模型中的常數Tab. 1 Constants in model

Cμ σk σε Cε1 Cε2 0.09 1.0 1.3 1.44 1.92

式中：

式中：n一般取值范圍為7～15，本文中n=10。

1.3 空化模型

Zwart空化模型是在Kubota[12]空化模型的基礎上做了一定的修正，其描述液-汽質量傳遞率的表達式如下：

Fvap和Fcond為有關相變的經驗系數（相變系數），分別代表蒸發系數和冷凝系數，默認的經驗值為50, 0.01[3,13]。本文主要研究對于不同頭型空化問題，這2個經驗系數的設置規范。

另外，研究表明湍動能對空化也有重要影響，目前普遍采用以下方式對汽化壓強修正[2]：

2 計算模型與參數

采用三維結構網格劃分不同頭型的圓柱模型所在流場，以半球頭圓柱為例，其網格模型如圖1所示。其中，邊界層首層厚度保證y+=50左右。計算中采用速度入口，流速U=10 m/s；壓力出口，相對壓力0 Pa。

流場中無窮遠處的參考壓力根據流速和不同空化數確定。空化數計算公式為：

圖 1 半球頭圓柱網格模型Fig. 1 Grid model of hemispherical cylinder

3 計算結果

表面壓力系數計算方法為：Hunter等[14]采用試驗方法研究了不同頭型在不同空化數條件下模型的表面壓力系數分布情況，這類頭型可以分為兩類：流線型和非流線型。流線型頭型的圓柱表面光順，流場不易產生分離渦，如半球頭圓柱、尖頂圓角圓柱等。非流線型頭型會使流場產生強烈的脈動和流動分離，如平頭圓柱、不同錐角的錐頭圓柱等。本文研究了不同相變系數對這兩類頭型空化計算結果的影響。

3.1 流線型頭型的空化計算

此空化區內的壓力系數與空化數互為相反數，也是半球頭圓柱體表面的最低壓力系數。圖2（a）中未顯示出空泡表面，說明流場內部氣相體積分數都在0.5以下，蒸發速率較慢。從表面壓力系數分布看，最低壓力低于飽和蒸汽壓，也說明蒸發速率不夠。主要原因是半球頭圓柱表面光順，周圍流場相對穩定，脈動和渦特性不明顯，不易產生空化現象。因此需要提高蒸發速率，使頭部周圍在低壓環境下快速形成空化區，目標是使壓力提高到飽和蒸汽壓。保持Fcond不變，增大Fvap，計算了如圖2（b）～圖2（f）所示的算例。

從圖2（b）～圖2（f）中可以看出，當蒸發系數Fvap達到5 000時，表面壓力系數最小值達到空化區內

當發生空化時，空泡內壓力為飽和蒸汽壓psat，因飽和蒸汽壓產生的壓力系數值，半球頭圓柱頭部圓周出現穩定空泡。然而圓柱表面壓力達到最低值后呈緩慢趨勢增加至零壓，未出現正壓力的尖峰突變。這種現象可以理解為冷凝速率較低，回射流不及時，空泡呈較慢的速度潰滅消失。因此需提高冷凝系數，驗證冷凝系數對空化現象的影響，計算圖2（g）～圖2（k）所示的算例。

從圖2（g）～圖2（k）可以看出，當冷凝系數Fcond達到0.4后，表面壓力系數分布已經很接近試驗值，出現正壓力峰值，圓柱頭部有一個主空泡，主空泡后出現次生空泡脫落現象，與局部空化的試驗現象相符。當冷凝系數Fcond達到0.7～1時，表面壓力系數與試驗值基本吻合。因此，對于半球頭圓柱，當空化數時，在數值計算過程中Fvap和Fcond宜分別設定為5 000，0.7～1。為驗證這樣的參數設置方案是否滿足其他空化數條件下的計算精度要求，計算如圖3所示的算例。

在一定空化數下，計算得到的表面壓力系數出現一段較平穩的最低值，這個壓力系數值與空化數互為相反數，代表這段區域壓強為飽和蒸汽壓，即出現了空泡，這段壓力系數最低值的長度范圍在一定程度上表明了空泡的長度尺度。由圖3可以看出，將蒸發系數設定為5 000，冷凝系數設定為0.7，1時，不同空化數下表面壓力系數的分布與試驗值吻合的很好，間接地驗證了數值計算得到的空泡尺度與試驗結果一致，只有當空化數為0.2時，在空泡潰滅區壓力分布與試驗值相比有些異常，但總體來看壓力變化趨勢與試驗值吻合，這樣的對比結果驗證了計算半球頭圓柱空化問題時相變系數設置的合理性。

為了驗證以上結論對其他流線型頭型的空化問題是否適用，基于以上設置規范計算了尖頂圓角圓柱在不同空化數小的壓力系數分布和空泡形態如圖4和圖5所示。

圖 2 時半球頭圓柱不同相變系數空化計算結果Fig. 2 Numerical results of hemispherical cylinder under different phase-change coefficients when

圖 3 半球頭圓柱不同空化數下表面壓力系數分布Fig. 3 Surface pressure coefficient distribution of hemispherical cylinder under different cavitation numbers

圖 4 尖頂圓角圓柱不同空化數下表面壓力系數分布Fig. 4 The surface pressure coefficient distribution of ogival rounded cylinder under different cavitation numbers

圖4中不同空化數下，數值計算得到的表面壓力系數與試驗值變化趨勢吻合，證明計算方法的準確性和相變系數設置的合理性。圖5為不同空化數下空泡的形態變化?？栈瘮禐?.4時空泡初生，當空化數達到0.24時，空泡特征為一個主空泡后夾雜即將脫落的次生空泡。

以上計算結果表明，對于流線頭型的空化問題，這樣的相變系數設定原則合理。

3.2 非流線型圓柱的空化計算

以平頭圓柱作為計算對象，計算其在空化數為0.4時不同相變系數下的壓力分布情況和縱剖面蒸汽體積分數分布云圖，如圖6所示。

從圖6可以看出，相變系數的改變對于平頭圓柱表面壓力系數分布影響不大。蒸發系數和冷凝系數的取值影響了流域內蒸汽體積分數的分布，并且流域內的蒸汽體積分數在平頭圓柱周圍的分布不均勻，具有強烈的不規則性。由于在獲得表面壓力系數分布曲線時是取縱剖面與圓柱表面交線的其中一側，因此空化的不規則性在一定程度上會影響表面壓力系數的分布情況。但是綜合來看，不同相變系數下計算的壓力分布結果與試驗值變化趨勢基本一致。

取Fvap=5000 ，Fcond=0.7（流線頭型空化計算建議值）和Fvap=50，Fcond=0.01（空化模型默認值）2組相變系數計算平頭圓柱在不同空化數下的表面壓力分布如圖7所示。

圖中顯示，不同相變系數下計算得到的表面壓力系數變化趨勢符合試驗現象，但是與試驗值存在一定的偏差，在空泡潰滅區壓力上升段，計算值相對試驗值有提前的現象。尤其當空化數較大時，試驗中平頭圓柱頭部周圍未達到空化壓力卻已經出現空化，這主要是由于平頭圓柱頭部周圍產生嚴重的流動分離導致過早的出現空化現象，當空化數較小時（）最小壓力才接近飽和蒸汽壓[13]。而滿足較小的空化數條件時，表面壓力系數分布對相變系數并不是很敏感，保證蒸發系數和冷凝系數的數量級關系即可模擬空化現象。

圖 6 時平頭圓柱不同相變系數下空化計算結果Fig. 6 Numerical results of of blunt cylindrical under different phase-change coefficients when

圖 7 平頭圓柱不同空化數表面壓力系數分布Fig. 7 Surface pressure coefficient distribution of blunt cylinder under different cavitation numbers

以上計算說明，對于非流線型圓柱的局部空化計算，只有當空化數小于0.3時，數值計算結果才能更接近試驗結果。令Fvap=5 000，Fcond=1計算90°錐頭圓柱在不同空化數條件下的空化現象，如圖8所示。當空化數時，計算結果的空化區長度與試驗結果整體變化趨勢一致，不同的是在空化潰滅區壓力變化存在一定差異。當空化數時，同樣存在流動分離導致空化過早出現的問題，而數值計算中的空化模型無法捕捉這一現象，因此在空化初生區域數值計算得到的空化區壓力系數比試驗結果小。

圖 8 90°錐頭圓柱表面壓力系數分布Fig. 8 Surface pressure coefficient distribution of 90°conical cylinder

4 結 語

采用基于CFD的數值計算方法結合Zwart空化模型求解不同頭型圓柱在不同空化數條件下的表面壓力系數，并與試驗數據對比，可以得到以下結論：

1）流線型頭型周圍流場穩定，不易產生空化，增大相變系數才能保證計算結果接近試驗結果，建議值為Fvap=5 000，Fcond=0.7～1.

2）非流線型頭型周圍流場流動分離和脈動強烈，在計算空化問題時對相變系數的變化不敏感，當空化數較?。ǎr，保證蒸發系數和冷凝系數的數量級關系即可較準確模擬空化現象。

后續工作將基于現有空化模型結構，考慮流場內流動分離強度因素，對汽化壓強再次修正，使其能夠適應非流線頭型在較大空化數條件下空化問題的求解。

[ 1 ]易文俊, 熊天紅. 高速航行體的自然超空泡流阻力特性研究[J]. 艦船科學技術, 2009, 31(1): 38–42.YI Wen-jun, XIONG Tian-hong. Research on drag characteristics of natural supercavitation profile for high speed bodies[J]. Ship Science and Technology, 2009, 31(1): 38–42.

[ 2 ]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.

[ 3 ]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 Japan: ICMF,2004.

[ 4 ]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.

[ 5 ]MERKLE C L, FENG 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.

[ 6 ]王柏秋, 王聰, 黃海龍, 等. 空化模型中的相變系數影響研究[J]. 工程力學, 2012, 29(8): 378–384.WANG Bai-qiu, WANG Cong, HUANG Hai-long, et al. Study of the influence of phase-chage coefficient in the cavitation model[J]. Engineering Mechanics, 2012, 29(8): 378–384.

[ 7 ]WANG Bai-qiu, WANG Cong, HUANG Hai-long, et al. Study of the phase-change coefficients of the cavitation model in cavitation flow fields generated from cone cavitator[J]. Journal of Harbin Institute of Technology, 2012, 19(4): 1–8.

[ 8 ]LIU Hou-lin, WANG Jian, WANG Yong, et al. Influence of the empirical coefficients of cavitation model on predicting cavitating flow in the centrifugal pump[J]. Int. J. Nav. Archit.Ocean Eng. 2014, 6: 119–131.

[ 9 ]黃彪, 王國玉, 胡常莉, 等. 繞回轉體初生空化流場特性的實驗及數值研究[J]. 工程力學, 2012, 29(6): 320–325.HUANG Biao, WANG Guo-yu, HU Chang-li, et al.Experimental study on fluctuating hydrodynamics around axisymmetric bodies[J]. Engineering Mechanics, 2012, 29(6):320–325.

[10]DULAR M, BACHERT R, STOFFEL B, et al. Experimental evaluation of numerical simulation of cavitating flow around hydrofoil[J]. European Journal of Mechanics-B/Fluids, 2005,24(4): 522–538.

[11]COUTIER-DELGOSHA O, FORTES-PATELLA R,REBOUD J L. Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation[J]. Journal of Fluids Engineering, 2003, 125(1): 38–45.

[12]KUBOTA A, KATO H, YAMAGUCHI H. A new modeling of cavitation flows: a numerical study of unsteady cavitation on a hydrofoil section[J]. Journal of Fluid Mechanics, 1992, 240(1):59–96.

[13]王復峰, 王國玉, 黃彪, 等. 通氣空化多項流動特性的實驗與數值研究[J]. 工程力學, 2016, 33(9): 220–226.WANG Fu-feng, WANG Guo-yu, HUANG Biao, et al.Experimental and numerical study on characteristic of ventilated cavitation multiphase flow[J]. Engineering Mechanics, 2016, 33(9): 220–226.

[14]HUNTER R, JOHN S M. Cavitation and pressure distribution,head forms at zero angle of yaw[R]. State University of Iowa,Studies in Engineering, Bulletin 32.