王 超，林 揚(yáng)，胡志強(qiáng)，衣瑞文

(1. 中國(guó)科學(xué)院沈陽(yáng)自動(dòng)化研究所 機(jī)器人學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 沈陽(yáng) 110016；2. 中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引 言

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

研究人員針對(duì)半球頭圓柱研究了Singhal空化模型中不同相變系數(shù)對(duì)計(jì)算結(jié)果的影響，指出不同空化數(shù)下相變系數(shù)的設(shè)置規(guī)范[6–7]。Liu等[8]在研究水泵葉輪的空化問(wèn)題后認(rèn)為，Zwart空化模型中的蒸發(fā)系數(shù)和冷凝系數(shù)的取值對(duì)葉輪上的空泡形態(tài)和尺度有重要影響。事實(shí)上，不同頭型（平頭、球頭、錐頭等）會(huì)導(dǎo)致周?chē)鲌?chǎng)產(chǎn)生不同的脈動(dòng)特性和漩渦特性，致使初生空泡形態(tài)和空化產(chǎn)生的難易程度不盡相同[9]。因此在利用空化模型模擬不同頭型的空化現(xiàn)象時(shí)，相變系數(shù)也不是一成不變的。

本文采用基于CFD的數(shù)值計(jì)算方法，研究利用Zwart空化模型模擬不同頭型的空化現(xiàn)象時(shí)相變系數(shù)對(duì)數(shù)值計(jì)算結(jié)果的影響，為進(jìn)一步提高空化問(wèn)題的數(shù)值計(jì)算精度提供指導(dǎo)性依據(jù)。

1 數(shù)值計(jì)算理論

1.1 N-S方程

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

1.2 湍流模型

耗散率ε方程為：

式（2）與式（3）中：

表 1 模型中的常數(shù)Tab. 1 Constants in model

表 1 模型中的常數(shù)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]空化模型的基礎(chǔ)上做了一定的修正，其描述液-汽質(zhì)量傳遞率的表達(dá)式如下：

Fvap和Fcond為有關(guān)相變的經(jīng)驗(yàn)系數(shù)（相變系數(shù)），分別代表蒸發(fā)系數(shù)和冷凝系數(shù)，默認(rèn)的經(jīng)驗(yàn)值為50, 0.01[3,13]。本文主要研究對(duì)于不同頭型空化問(wèn)題，這2個(gè)經(jīng)驗(yàn)系數(shù)的設(shè)置規(guī)范。

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

2 計(jì)算模型與參數(shù)

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

流場(chǎng)中無(wú)窮遠(yuǎn)處的參考?jí)毫Ω鶕?jù)流速和不同空化數(shù)確定?？栈瘮?shù)計(jì)算公式為：

圖 1 半球頭圓柱網(wǎng)格模型Fig. 1 Grid model of hemispherical cylinder

3 計(jì)算結(jié)果

表面壓力系數(shù)計(jì)算方法為：Hunter等[14]采用試驗(yàn)方法研究了不同頭型在不同空化數(shù)條件下模型的表面壓力系數(shù)分布情況，這類(lèi)頭型可以分為兩類(lèi)：流線型和非流線型。流線型頭型的圓柱表面光順，流場(chǎng)不易產(chǎn)生分離渦，如半球頭圓柱、尖頂圓角圓柱等。非流線型頭型會(huì)使流場(chǎng)產(chǎn)生強(qiáng)烈的脈動(dòng)和流動(dòng)分離，如平頭圓柱、不同錐角的錐頭圓柱等。本文研究了不同相變系數(shù)對(duì)這兩類(lèi)頭型空化計(jì)算結(jié)果的影響。

3.1 流線型頭型的空化計(jì)算

此空化區(qū)內(nèi)的壓力系數(shù)與空化數(shù)互為相反數(shù)，也是半球頭圓柱體表面的最低壓力系數(shù)。圖2（a）中未顯示出空泡表面，說(shuō)明流場(chǎng)內(nèi)部氣相體積分?jǐn)?shù)都在0.5以下，蒸發(fā)速率較慢。從表面壓力系數(shù)分布看，最低壓力低于飽和蒸汽壓，也說(shuō)明蒸發(fā)速率不夠。主要原因是半球頭圓柱表面光順，周?chē)鲌?chǎng)相對(duì)穩(wěn)定，脈動(dòng)和渦特性不明顯，不易產(chǎn)生空化現(xiàn)象。因此需要提高蒸發(fā)速率，使頭部周?chē)诘蛪涵h(huán)境下快速形成空化區(qū)，目標(biāo)是使壓力提高到飽和蒸汽壓。保持Fcond不變，增大Fvap，計(jì)算了如圖2（b）～圖2（f）所示的算例。

從圖2（b）～圖2（f）中可以看出，當(dāng)蒸發(fā)系數(shù)Fvap達(dá)到5 000時(shí)，表面壓力系數(shù)最小值達(dá)到空化區(qū)內(nèi)

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

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

在一定空化數(shù)下，計(jì)算得到的表面壓力系數(shù)出現(xiàn)一段較平穩(wěn)的最低值，這個(gè)壓力系數(shù)值與空化數(shù)互為相反數(shù)，代表這段區(qū)域壓強(qiáng)為飽和蒸汽壓，即出現(xiàn)了空泡，這段壓力系數(shù)最低值的長(zhǎng)度范圍在一定程度上表明了空泡的長(zhǎng)度尺度。由圖3可以看出，將蒸發(fā)系數(shù)設(shè)定為5 000，冷凝系數(shù)設(shè)定為0.7，1時(shí)，不同空化數(shù)下表面壓力系數(shù)的分布與試驗(yàn)值吻合的很好，間接地驗(yàn)證了數(shù)值計(jì)算得到的空泡尺度與試驗(yàn)結(jié)果一致，只有當(dāng)空化數(shù)為0.2時(shí)，在空泡潰滅區(qū)壓力分布與試驗(yàn)值相比有些異常，但總體來(lái)看壓力變化趨勢(shì)與試驗(yàn)值吻合，這樣的對(duì)比結(jié)果驗(yàn)證了計(jì)算半球頭圓柱空化問(wèn)題時(shí)相變系數(shù)設(shè)置的合理性。

為了驗(yàn)證以上結(jié)論對(duì)其他流線型頭型的空化問(wèn)題是否適用，基于以上設(shè)置規(guī)范計(jì)算了尖頂圓角圓柱在不同空化數(shù)小的壓力系數(shù)分布和空泡形態(tài)如圖4和圖5所示。

圖 2 時(shí)半球頭圓柱不同相變系數(shù)空化計(jì)算結(jié)果Fig. 2 Numerical results of hemispherical cylinder under different phase-change coefficients when

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

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

圖4中不同空化數(shù)下，數(shù)值計(jì)算得到的表面壓力系數(shù)與試驗(yàn)值變化趨勢(shì)吻合，證明計(jì)算方法的準(zhǔn)確性和相變系數(shù)設(shè)置的合理性。圖5為不同空化數(shù)下空泡的形態(tài)變化?？栈瘮?shù)為0.4時(shí)空泡初生，當(dāng)空化數(shù)達(dá)到0.24時(shí)，空泡特征為一個(gè)主空泡后夾雜即將脫落的次生空泡。

以上計(jì)算結(jié)果表明，對(duì)于流線頭型的空化問(wèn)題，這樣的相變系數(shù)設(shè)定原則合理。

3.2 非流線型圓柱的空化計(jì)算

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

從圖6可以看出，相變系數(shù)的改變對(duì)于平頭圓柱表面壓力系數(shù)分布影響不大。蒸發(fā)系數(shù)和冷凝系數(shù)的取值影響了流域內(nèi)蒸汽體積分?jǐn)?shù)的分布，并且流域內(nèi)的蒸汽體積分?jǐn)?shù)在平頭圓柱周?chē)姆植疾痪鶆?，具有?qiáng)烈的不規(guī)則性。由于在獲得表面壓力系數(shù)分布曲線時(shí)是取縱剖面與圓柱表面交線的其中一側(cè)，因此空化的不規(guī)則性在一定程度上會(huì)影響表面壓力系數(shù)的分布情況。但是綜合來(lái)看，不同相變系數(shù)下計(jì)算的壓力分布結(jié)果與試驗(yàn)值變化趨勢(shì)基本一致。

取Fvap=5000 ，F(xiàn)cond=0.7（流線頭型空化計(jì)算建議值）和Fvap=50，F(xiàn)cond=0.01（空化模型默認(rèn)值）2組相變系數(shù)計(jì)算平頭圓柱在不同空化數(shù)下的表面壓力分布如圖7所示。

圖中顯示，不同相變系數(shù)下計(jì)算得到的表面壓力系數(shù)變化趨勢(shì)符合試驗(yàn)現(xiàn)象，但是與試驗(yàn)值存在一定的偏差，在空泡潰滅區(qū)壓力上升段，計(jì)算值相對(duì)試驗(yàn)值有提前的現(xiàn)象。尤其當(dāng)空化數(shù)較大時(shí)，試驗(yàn)中平頭圓柱頭部周?chē)催_(dá)到空化壓力卻已經(jīng)出現(xiàn)空化，這主要是由于平頭圓柱頭部周?chē)a(chǎn)生嚴(yán)重的流動(dòng)分離導(dǎo)致過(guò)早的出現(xiàn)空化現(xiàn)象，當(dāng)空化數(shù)較小時(shí)（）最小壓力才接近飽和蒸汽壓[13]。而滿足較小的空化數(shù)條件時(shí)，表面壓力系數(shù)分布對(duì)相變系數(shù)并不是很敏感，保證蒸發(fā)系數(shù)和冷凝系數(shù)的數(shù)量級(jí)關(guān)系即可模擬空化現(xiàn)象。

圖 6 時(shí)平頭圓柱不同相變系數(shù)下空化計(jì)算結(jié)果Fig. 6 Numerical results of of blunt cylindrical under different phase-change coefficients when

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

以上計(jì)算說(shuō)明，對(duì)于非流線型圓柱的局部空化計(jì)算，只有當(dāng)空化數(shù)小于0.3時(shí)，數(shù)值計(jì)算結(jié)果才能更接近試驗(yàn)結(jié)果。令Fvap=5 000，F(xiàn)cond=1計(jì)算90°錐頭圓柱在不同空化數(shù)條件下的空化現(xiàn)象，如圖8所示。當(dāng)空化數(shù)時(shí)，計(jì)算結(jié)果的空化區(qū)長(zhǎng)度與試驗(yàn)結(jié)果整體變化趨勢(shì)一致，不同的是在空化潰滅區(qū)壓力變化存在一定差異。當(dāng)空化數(shù)時(shí)，同樣存在流動(dòng)分離導(dǎo)致空化過(guò)早出現(xiàn)的問(wèn)題，而數(shù)值計(jì)算中的空化模型無(wú)法捕捉這一現(xiàn)象，因此在空化初生區(qū)域數(shù)值計(jì)算得到的空化區(qū)壓力系數(shù)比試驗(yàn)結(jié)果小。

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

4 結(jié) 語(yǔ)

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

1）流線型頭型周?chē)鲌?chǎng)穩(wěn)定，不易產(chǎn)生空化，增大相變系數(shù)才能保證計(jì)算結(jié)果接近試驗(yàn)結(jié)果，建議值為Fvap=5 000，F(xiàn)cond=0.7～1.

2）非流線型頭型周?chē)鲌?chǎng)流動(dòng)分離和脈動(dòng)強(qiáng)烈，在計(jì)算空化問(wèn)題時(shí)對(duì)相變系數(shù)的變化不敏感，當(dāng)空化數(shù)較?。ǎr(shí)，保證蒸發(fā)系數(shù)和冷凝系數(shù)的數(shù)量級(jí)關(guān)系即可較準(zhǔn)確模擬空化現(xiàn)象。

后續(xù)工作將基于現(xiàn)有空化模型結(jié)構(gòu)，考慮流場(chǎng)內(nèi)流動(dòng)分離強(qiáng)度因素，對(duì)汽化壓強(qiáng)再次修正，使其能夠適應(yīng)非流線頭型在較大空化數(shù)條件下空化問(wèn)題的求解。

[ 1 ]易文俊, 熊天紅. 高速航行體的自然超空泡流阻力特性研究[J]. 艦船科學(xué)技術(shù), 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 ]王柏秋, 王聰, 黃海龍, 等. 空化模型中的相變系數(shù)影響研究[J]. 工程力學(xué), 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 ]黃彪, 王國(guó)玉, 胡常莉, 等. 繞回轉(zhuǎn)體初生空化流場(chǎng)特性的實(shí)驗(yàn)及數(shù)值研究[J]. 工程力學(xué), 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]王復(fù)峰, 王國(guó)玉, 黃彪, 等. 通氣空化多項(xiàng)流動(dòng)特性的實(shí)驗(yàn)與數(shù)值研究[J]. 工程力學(xué), 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.