黃 彪，王國玉，張 博，時(shí)素果

（北京理工大學(xué) 機(jī)械與車輛工程學(xué)院，北京 100081）

1 引 言

空化現(xiàn)象不僅和汽液相變過程相關(guān)，還涉及到汽液兩相的大規(guī)模旋渦運(yùn)動(dòng)，是一種復(fù)雜的非定常多相湍流流動(dòng)。Arakeri和Acosta采用全息攝影的方法研究了繞軸對稱物體的水流中的空化現(xiàn)象，發(fā)現(xiàn)空化發(fā)生與水流的分離區(qū)域的旋渦運(yùn)動(dòng)有關(guān)。Katz[1]發(fā)現(xiàn)分離區(qū)域的軸向剪切渦結(jié)構(gòu)影響空化的發(fā)展。Laberteraux等[2]利用高速攝像觀察發(fā)現(xiàn)附著型空化閉合區(qū)域存在著空化渦形態(tài)結(jié)構(gòu)，王國玉等[3]用高速攝影和錄像的方法觀察了發(fā)生在射流放水閥內(nèi)部的粘性非定常水流運(yùn)動(dòng)，指出空化的形成過程和旋渦運(yùn)動(dòng)有關(guān)?？栈郎u流動(dòng)對非定常空化流動(dòng)的數(shù)值模擬提出了巨大的挑戰(zhàn)，不僅要求空化模型能精確描述空化的形態(tài)發(fā)展過程，還要求能夠?qū)Π殡S的大規(guī)模旋渦運(yùn)動(dòng)進(jìn)行預(yù)測。

近年來，人們?yōu)榱诉M(jìn)行非定?？栈鲃?dòng)計(jì)算，主要采用了兩類空化模型，即正壓狀態(tài)方程的單流體模型和基于相間傳輸?shù)木嗔髂Ｐ?。由于基于正壓流體法則，密度和壓力項(xiàng)具有相同的變化梯度，斜壓矩為零，使得使用正壓流體方程不能很好地捕捉實(shí)驗(yàn)觀測到的空化旋渦流動(dòng)特性。Kubata[4]、Kunz[5]、Singhal[6]、Senocak和Shyy[8]分別基于相間的傳輸模型，通過添加不同源項(xiàng)來調(diào)節(jié)汽液兩相間的傳輸。其中Kubota空化模型和Singhal空化模型在目前得到了比較廣泛的應(yīng)用。蔣增輝等[7]使用商業(yè)軟件Fluent中的Singhal空化模型，對不同的空泡尾部流場特性進(jìn)行數(shù)值模擬研究。張博等[9]應(yīng)用商業(yè)軟件CFX的Kubota空化模型，描述了非定常云狀空化形態(tài)的發(fā)展過程。雖然目前已有較多的使用這些模型來研究非定?？栈奈墨I(xiàn)，但是，由于在一種軟件中僅包含一種模型，缺乏對兩種模型的應(yīng)用對比，特別是對空化流場的旋渦流場結(jié)構(gòu)模擬的對比研究還不是很充分。

為了評價(jià)空化模型對非定?？栈鲃?dòng)中的應(yīng)用，作者通過商業(yè)軟件CFX的二次開發(fā)技術(shù)將Singhal空化模型以及湍流模型鑲?cè)肓擞?jì)算軟件，分別采用Kubota和Singhal模型計(jì)算了繞Clark-y型水翼云狀空化流動(dòng)，分析了不同空化模型對計(jì)算結(jié)果的影響，特別是對空化流動(dòng)中的時(shí)均速度及渦量場的特性的影響，通過這兩種空化模型計(jì)算結(jié)果和實(shí)驗(yàn)現(xiàn)象的對比，對空化模型進(jìn)行了討論和評價(jià)。

2 數(shù)學(xué)模型和數(shù)值計(jì)算方法

2.1 基本方程

假定汽液兩相為均相流動(dòng)，相間無速度滑移，汽液兩相的連續(xù)方程和動(dòng)量方程如下所示：

式中，ρ為混合密度，μ為混相介質(zhì)的動(dòng)力粘性系數(shù)，為速度，P為壓強(qiáng)，μt為湍流粘性系數(shù)。

2.2 空化模型

2.2.1 Kubota空化模型[4]

Kubota空化模型基于單個(gè)空泡的生成和發(fā)展時(shí)空泡體積變化，基于Rayleigh-Plesset空泡生長方程推導(dǎo)出了如下的蒸發(fā)和凝結(jié)相的表達(dá)式：

式中，RB為平均空泡半徑，F(xiàn)e和Fc分別是蒸發(fā)和凝結(jié)常數(shù)項(xiàng)，隨蒸發(fā)和凝結(jié)而異，設(shè)計(jì)時(shí)應(yīng)考慮它們發(fā)生的速率不同。（凝結(jié)通常比蒸發(fā)慢得多）

2.2.2 Singhal空化模型[6]

Singhal空化模型在推導(dǎo)過程中基于汽、液兩相流模型，由混合密度導(dǎo)出相間質(zhì)量傳輸速率，綜合考慮了空泡在相變過程中所受阻力和表面張力和實(shí)際流體中的非凝結(jié)性氣體含量的影響。

式中，Ce和Cc是隨蒸發(fā)與凝結(jié)程度不同而變化的經(jīng)驗(yàn)系數(shù)，fv、fg分別為汽相、不溶性氣核的質(zhì)量分?jǐn)?shù)，vch是特征速度，反映氣、液兩相之間的相對速度。

2.3 RNG k-ε湍流模型及其修正

計(jì)算采用了由Yakhot和Orzag提出的RNG k-ε湍流模型：

式中，k、ε 分別為湍動(dòng)能和湍流耗散率，Gk為湍動(dòng)能生成項(xiàng)，C1ε和 C2ε為經(jīng)驗(yàn)常數(shù)，σk=σε=1.39，Cμ=0.09。

由于空化區(qū)內(nèi)含有大量的水蒸汽，是一種水汽混合介質(zhì)，考慮汽液兩相混合密度的變化對湍流粘性系數(shù)的影響，這里對RNG k-ε模型進(jìn)行了修正，應(yīng)用一個(gè)密度函數(shù)f（ρ）代替（9）式中的混合密度，湍流粘性系數(shù)采用以下兩式進(jìn)行計(jì)算：

對于（11）式中n的取值，張博等[9]通過取不同的n值，分析了對空化流動(dòng)數(shù)值計(jì)算結(jié)果的影響，認(rèn)為在云狀空化階段n取10是合理的。

2.4 數(shù)值計(jì)算方法

2.4.1 計(jì)算網(wǎng)格和邊界條件

計(jì)算采用了和試驗(yàn)[10]相同的Clark-y型水翼和流動(dòng)條件。水翼的弦長C=0.07 m。圖1給出了計(jì)算區(qū)域的網(wǎng)格分區(qū)及其邊界條件。翼型前端的區(qū)域采用c型結(jié)構(gòu)化網(wǎng)格劃分，這樣可以較好地匹配翼型頭部的形狀。計(jì)算區(qū)域的入口距翼型前緣為2.5C，出口距翼型尾緣的距離為5C。采用了文獻(xiàn)[9]中的網(wǎng)格設(shè)置，總網(wǎng)格數(shù)為43 000，并如圖2所示在翼型周圍近壁區(qū)域進(jìn)行了網(wǎng)格加密，近壁面y+值為20～80之間，滿足壁面函數(shù)要求。

本文采用速度入口、壓力出口的邊界條件，流動(dòng)區(qū)域上下邊界為自由滑移壁面條件，翼型表面采用絕熱、自由滑移固壁條件。根據(jù)試驗(yàn)工況對計(jì)算參數(shù)進(jìn)行相應(yīng)設(shè)置，攻角設(shè)定為α=8°，空化數(shù)設(shè)定為σ=0.8，流速 U∞設(shè)定為 10 m/s，對應(yīng)的雷諾數(shù)為 7×105。

圖1 計(jì)算區(qū)域網(wǎng)格塊和邊界條件Fig.1 Outline of the computational domain with boundary condition

圖2 翼型周圍網(wǎng)格Fig.2 Computational grid for flow analysis

2.4.2 無量綱數(shù)的定義

計(jì)算中的主要無量綱參數(shù)為空化數(shù)σ和壓力系數(shù)cp，分別定義為：

這里p∞、U∞和pv分別為距試驗(yàn)段上游入口210 mm處參考斷面上的平均靜壓強(qiáng)、斷面平均速度和汽化壓強(qiáng)，p為所取試驗(yàn)點(diǎn)的當(dāng)?shù)貕簭?qiáng)。

3 結(jié)果與討論

3.1 非定常空穴形態(tài)

圖3給出了分別采用兩種空化模型計(jì)算所得空泡形態(tài)隨時(shí)間的變化及其與試驗(yàn)觀測結(jié)果[8]的對比。圖3（a）是試驗(yàn)觀測結(jié)果，圖3（b）、(c)分別為采用Kubota與Singhal空化模型計(jì)算得到的結(jié)果，計(jì)算所得的空泡形態(tài)用含汽率分布表示，藍(lán)色區(qū)域穴內(nèi)水蒸氣含量基本在80%以上，紅色部分基本上為汽相。試驗(yàn)觀測所得的空穴形態(tài)由圖中白色汽相區(qū)域表示。

圖3 空穴形態(tài)隨時(shí)間的變化Fig.3 Time evolution of cavity shape in the experiment and the calculation

采用Kubota和Singhal空化模型所得的計(jì)算結(jié)果都清楚地描述了云狀空化的產(chǎn)生—發(fā)展—脫落的準(zhǔn)周期性變化，這與試驗(yàn)觀測基本一致。在云狀空化開始階段，如圖3，當(dāng)t=t0+3.5 ms時(shí)，翼型頭部低壓區(qū)首先產(chǎn)生厚度很小的微空泡與液滴組成的空穴，此時(shí)空穴附著在翼型表面上；隨著時(shí)間的推移，附著空穴不斷地向翼型尾部發(fā)展，其厚度沿著翼弦的方向不斷增加。在t=t0+14.7 ms時(shí)，空穴長度達(dá)到最大值，空穴末端達(dá)到翼型尾部，這時(shí)空穴末端近壁面區(qū)域的水汽混合區(qū)出現(xiàn)回縮，在t=t0+32.2 ms時(shí)，翼型尾部出現(xiàn)大規(guī)模的空泡團(tuán)旋渦脫落現(xiàn)象。

兩種空化模型對非定?？栈螒B(tài)的計(jì)算結(jié)果的主要差異表現(xiàn)在描述空穴回縮、空泡團(tuán)脫落過程中。采用Singhal空化模型計(jì)算結(jié)果空泡團(tuán)是以大尺度旋渦形式脫落，和試驗(yàn)觀測的結(jié)果基本一致，而采用Kubota空化模型計(jì)算的空穴變化主體為空穴整體壓縮，隨后，在t=t0+32.2 ms時(shí)刻空穴斷裂成兩個(gè)部分，空穴前端附著在翼型吸力面上，空穴尾段形成了斷裂的空泡團(tuán)，并向下游運(yùn)動(dòng)，和試驗(yàn)觀測的流動(dòng)現(xiàn)象有較大的差異。兩種空化模型在模擬空穴形態(tài)變化過程中的差異，是由于兩種空化模型建立過程對空化現(xiàn)象所做的假設(shè)不同造成的。盡管Kubota與Singhal空化模型都是用來描述汽液兩相間質(zhì)量交換的，但Kubota空化模型在推導(dǎo)過程中，將空泡體積變化率作為重要參量，并且假定蒸發(fā)與凝結(jié)源項(xiàng)與平均空泡半徑成反比，空穴的增長和潰滅可以視為無數(shù)空泡體積變化的整體行為，空穴內(nèi)部的含汽量幾乎恒定，而Singhal空化模型基于空化流動(dòng)中敏感的密度變化，空穴內(nèi)部的密度處于非定常變化的狀態(tài)，因此采用Kubota描述的云狀空化的脫離過程更近視于單個(gè)空泡的潰滅過程，而不能體現(xiàn)空化的復(fù)雜多相流動(dòng)特性，而Singhal模型則更多地反映了汽液多相流動(dòng)過程的大規(guī)模漩渦運(yùn)動(dòng)，能更精確地描述云狀空化的非定常過程。

3.2 繞翼型速度矢量場

為了對比兩種空化模型所計(jì)算的空化區(qū)域的流動(dòng)結(jié)構(gòu)，分別提取了空泡脫落過程中幾個(gè)典型時(shí)刻繞翼型瞬時(shí)速度矢量圖，圖4為分別采用Kubota、Singhal空化模型計(jì)算得到繞翼型瞬時(shí)速度場，圖5給出了對應(yīng)時(shí)刻翼型尾部局部放大圖和試驗(yàn)觀測圖[10]。

圖4 典型時(shí)刻繞翼型速度場Fig.4 Instantaneous velocity of cloud cavitation around of the foil

圖5 典型時(shí)刻繞翼型尾部速度場Fig.5 Instantaneous velocity of cloud cavitation in the rear of the foil

文獻(xiàn)[10]中給出了試驗(yàn)中空穴區(qū)域內(nèi)部大規(guī)模的旋渦運(yùn)動(dòng)過程。指出云狀空化后期，在緊貼壁面的區(qū)域內(nèi)誘導(dǎo)了一股指向上游的反向流動(dòng)，在反向射流的作用下，翼型尾部的大尺度空泡團(tuán)是以強(qiáng)烈的旋渦分離的形式脫離壁面的。

采用Kubota、Singhal空化模型數(shù)值計(jì)算的結(jié)果均體現(xiàn)了反向射流的產(chǎn)生。t0+17.2 ms、t0+19.6 ms、t0+26.6 ms為空穴開始回縮，反向射流逐漸向翼型前緣發(fā)展的時(shí)刻。采用不同空化模型所描述的空穴內(nèi)部的流場結(jié)構(gòu)是相同的：在空泡回縮過程中，采用Kubota空化模型計(jì)算得到的結(jié)果（圖4（a)）表明：在翼型尾段近壁面，水汽兩相混合區(qū)沿壁面反向流動(dòng)，空穴區(qū)域內(nèi)遠(yuǎn)離壁面的流體與主流方向一致，在空穴內(nèi)部，不存在明顯的旋渦，而采用Singhal空化模型計(jì)算結(jié)果（圖4（b）)會(huì)有明顯的改善：反向流動(dòng)向壁面發(fā)展，并且在逐漸遠(yuǎn)離壁面的空化區(qū)域內(nèi)，會(huì)產(chǎn)生非常明顯的順時(shí)針的旋渦，旋渦有向翼型后上方移動(dòng)的趨勢，使空泡團(tuán)逐漸脫離壁面，符合實(shí)驗(yàn)觀測現(xiàn)象。為了進(jìn)一步闡明流動(dòng)結(jié)構(gòu)的發(fā)展，圖5給出了實(shí)驗(yàn)與數(shù)值計(jì)算得到的翼型尾部局部的特寫圖像，通過對比實(shí)驗(yàn)現(xiàn)象可以發(fā)現(xiàn)：采用Singhal空化模型準(zhǔn)確地模擬出了翼型尾部的空化旋渦區(qū)，空泡的旋渦脫落現(xiàn)象非常明顯。

3.3 繞翼型旋渦特性

圖6為試驗(yàn)[10]和采用不同空化模型獲得的云狀空化階段沿來流方向上的速度時(shí)均值u沿翼弦的分布，從圖中可以看出：在云狀空化階段，主流與空化流中的速度分布存在明顯的界面，空穴內(nèi)部的混合物流動(dòng)速度由壁面附近呈現(xiàn)出較大的速度梯度變化，并逐漸與主流區(qū)的速度相同，點(diǎn)劃線a與實(shí)線b與翼型表面之間的區(qū)域?yàn)榭栈闪髅}動(dòng)區(qū)域，對應(yīng)于空化旋渦區(qū)。對比試驗(yàn)與數(shù)值計(jì)算的結(jié)果，采用Singhal空化模型得到的時(shí)均空化紊流脈動(dòng)區(qū)域大小與試驗(yàn)比較接近，Kubota空化模型過高地預(yù)計(jì)了剪切層的速度梯度變化。

圖6 云狀空化階段速度u時(shí)均分布圖Fig.6 Average streamwise velocity distribution at the cloud cavitation stage

渦量來源于流場存在速度梯度，是描述有旋流動(dòng)的一個(gè)運(yùn)動(dòng)學(xué)物理量，在研究云狀空化現(xiàn)象中，涉及到空泡團(tuán)脫落等分離流動(dòng)情況，對渦量進(jìn)行分析具有很重要的意義，分別采用Kubota與Singhal空化模型對繞翼型渦量場進(jìn)行描述，利用相關(guān)后處理軟件得到了對應(yīng)空化條件下瞬時(shí)z向渦量分布圖，這里，z向渦量被定義為：

圖7為數(shù)值計(jì)算得到的z向時(shí)均渦量沿翼弦的分布。從圖中可以看出，由于反向射流的存在，采用不同空化模型計(jì)算均可得到：在近壁面處，都會(huì)產(chǎn)生旋渦空化區(qū)，隨著離翼型壁面垂直距離y的不同而變化，同一截面上的渦量均會(huì)發(fā)生變化，在遠(yuǎn)離壁面區(qū)域，渦量趨近于0。采用不同空化模型計(jì)算結(jié)果的差異主要體現(xiàn)在對空化旋渦區(qū)的模擬上，如圖所示，點(diǎn)劃線a與實(shí)線b與翼型壁面所圍成的區(qū)域分別為采用Singhal與Kubota空化模型計(jì)算得到的空化漩渦區(qū)，采用Singhal空化模型計(jì)算得到的空化旋渦區(qū)域明顯比采用Kubota時(shí)的大，也就是說，采用Singhal空化模型計(jì)算結(jié)果體現(xiàn)了更為明顯的旋渦特性。

圖7 云狀空化階段z向渦量時(shí)均分布圖Fig.7 Average vorticity distribution at the cloud cavitation stage

不同空化模型在模擬繞翼型云狀空化速度與渦量場上的差異反映了不同空化模型所體現(xiàn)的物理機(jī)制是不同的，Singhal空化模型在推導(dǎo)相間質(zhì)量傳輸速率過程中，添加了瞬時(shí)粘性力與湍流波動(dòng)的影響，綜合考慮了汽、液兩相間的相對速度、空泡表面張力等諸多因素，從而可更好地考慮空化現(xiàn)象中非常頻繁的質(zhì)量與速度交換，而Kubota空化模型則忽略了空化流動(dòng)相變過程中頻繁的動(dòng)量交換與湍流波動(dòng)等因素，所以，在空化區(qū)域內(nèi)，采用Singhal空化模型計(jì)算得到的速度以及動(dòng)量交換要比Kubota空化模型劇烈，體現(xiàn)出更為明顯的空化旋渦特性。

4 結(jié) 論

分別應(yīng)用Kubota、Singhal空化模型計(jì)算了繞Clark-y型水翼云狀空化流動(dòng)，分析了繞翼型云狀空化的空穴形態(tài)、空泡團(tuán)脫落形式、時(shí)均速度、渦量場，并與試驗(yàn)結(jié)果進(jìn)行了對比分析，結(jié)論如下：

（1）兩種空化模型數(shù)值計(jì)算結(jié)果與試驗(yàn)現(xiàn)象均清楚地描述了云狀非定?？栈^程中的產(chǎn)生—發(fā)展—脫落的準(zhǔn)周期性變化，但兩種模型對空泡團(tuán)的脫落機(jī)理描述是不同的。Kubota模型將脫落過程描述為空泡團(tuán)的壓縮，這和試驗(yàn)結(jié)果是不一致的，而Singhal空化模型可以獲得和試驗(yàn)觀測一致的計(jì)算結(jié)果。

（2）Kubota空化模型模擬出來的空化旋渦區(qū)明顯小于試驗(yàn)，而采用Singhal空化模型計(jì)算得到的沿翼弦方向的時(shí)均速度分布與試驗(yàn)結(jié)果一致，可較好地模擬出空化區(qū)域內(nèi)的旋渦空化區(qū)域，與實(shí)驗(yàn)現(xiàn)象相符。

[1]Katz J.Cavitation phenomena within region of flow separation[J].Journal of Fluid Mechanics,1984,140(4):2397-2436.

[2]Laberteraux K R,Ceccio S L,Mast rocola V J,et al.High speed digital imaging of cavitating vortices[J].Experiments in Fluids,1998,24(6):489-498.

[3]王國玉,曹樹良,井小萩利明.剪切層區(qū)域旋渦空化的發(fā)生機(jī)理[J].清華大學(xué)學(xué)報(bào):自然科學(xué)版,2001,41(10):62-64,89.Wang Guoyu,Cao Shuliang,Ikohagi Toshiaki.Vortex cavitation mechanisms in shear layer flow[J].Tsinghua Univ:Sci&Tech,2001,41(10):62-64,89.

[4]Kubota A,Kato H.Unsteady structure measurement of cloud cavitation on a foil section using conditional sampling techniques[J].Journal of Fluids Engineering,1989,111(3):204-210.

[5]Kunz R,Boger D,Stinebring D,Chyczewski T.A preconditioned implicit method for two-phase flows with application to cavitation prediction[J].Comput.Fluids,2000,29(8):849-875.

[6]Singhal A K,Athavale M M.Mathematical basis and validation of the full cavitation model[J].Journal of Fluids Engineering,2002,124:617-624.

[7]蔣增輝,于開平,張嘉鐘等.空泡尾部流場特性數(shù)值模擬研究[J].船舶力學(xué),2008,12(2):225-230.Jiang Zenghui,Yu Kaiping,Zhang Jiazhong.Numerical simulation on character of the cavity wake region[J].Journal of Ship Mechanics,2008,12(2):225-230.

[8]Senocak I,Shyy W.A pressure-based method for turbulent cavitating flow computations[J].Journal of Computional Physics,2002,176(2):363-383.

[9]張 博,王國玉等.修正的RNG k-ε模型在云狀空化流動(dòng)計(jì)算中的應(yīng)用評價(jià)[J].北京理工大學(xué)學(xué)報(bào),2008,28(12):1065-1069.Zhang B,Wang G,et al.Evaluation of a modified RNG k-ε model for computations of cloud cavitating flow[J].Transac-tions of Beijing Institute of Technology,2008,28(12):1065-1069.

[10]Wang G,Senocak I,Shyy W,et al.Dynamics of attached turbulent cavitating flows[J].Progress in Aerospace Sciences,2001,37:551-581.

[11]王 志,陳九錫.回轉(zhuǎn)體空泡繞流現(xiàn)象的數(shù)值模擬[J].船舶力學(xué),2006,8(4):38-43.Wang Zhi,Chen jiuxi.Numerical simulation of cavitation flow on revolution body[J].Journal of Ship Mechanics,2006,8(4):38-43.

[12]Zhou L J,Wang Z W.Numerical simulation of cavitation around a hydrofoil and evaluation of a RNG k-ε model[J].Journal of Fluids Engineering,2008,130(1):1-7.

[13]Coutier-Delgosha O.Numerical prediction of cavitation flow on a two-dimensional symmetrical hydrofoil and comparison to experiments[J].Journal of Fluids Engineering,2007,129(3):279-291.

[14]馮 光,顏 開.超空泡航行體水下彈道的數(shù)值計(jì)算[J].船舶力學(xué),2005,9(2):1-8.Feng guang,Yan kai.Numerical calculation of underwater trajectory of supercavitating bodies[J].Journal of Ship Mechanics,2005,9(2):1-8.