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

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

1 引 言

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

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

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

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

2.1 基本方程

假定汽液兩相為均相流動(dòng)，相間無(wú)速度滑移，汽液兩相的連續(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空泡生長(zhǎng)方程推導(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)過(guò)程中基于汽、液兩相流模型，由混合密度導(dǎo)出相間質(zhì)量傳輸速率，綜合考慮了空泡在相變過(guò)程中所受阻力和表面張力和實(shí)際流體中的非凝結(jié)性氣體含量的影響。

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

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ì)，考慮汽液兩相混合密度的變化對(duì)湍流粘性系數(shù)的影響，這里對(duì)RNG k-ε模型進(jìn)行了修正，應(yīng)用一個(gè)密度函數(shù)f（ρ）代替（9）式中的混合密度，湍流粘性系數(shù)采用以下兩式進(jìn)行計(jì)算：

對(duì)于（11）式中n的取值，張博等[9]通過(guò)取不同的n值，分析了對(duì)空化流動(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)條件。水翼的弦長(zhǎ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)工況對(duì)計(jì)算參數(shù)進(jìn)行相應(yīng)設(shè)置，攻角設(shè)定為α=8°，空化數(shù)設(shè)定為σ=0.8，流速 U∞設(shè)定為 10 m/s，對(duì)應(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 無(wú)量綱數(shù)的定義

計(jì)算中的主要無(wú)量綱參數(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)觀測(cè)結(jié)果[8]的對(duì)比。圖3（a）是試驗(yàn)觀測(cè)結(jié)果，圖3（b）、(c)分別為采用Kubota與Singhal空化模型計(jì)算得到的結(jié)果，計(jì)算所得的空泡形態(tài)用含汽率分布表示，藍(lán)色區(qū)域穴內(nèi)水蒸氣含量基本在80%以上，紅色部分基本上為汽相。試驗(yàn)觀測(cè)所得的空穴形態(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)觀測(cè)基本一致。在云狀空化開(kāi)始階段，如圖3，當(dāng)t=t0+3.5 ms時(shí)，翼型頭部低壓區(qū)首先產(chǎn)生厚度很小的微空泡與液滴組成的空穴，此時(shí)空穴附著在翼型表面上；隨著時(shí)間的推移，附著空穴不斷地向翼型尾部發(fā)展，其厚度沿著翼弦的方向不斷增加。在t=t0+14.7 ms時(shí)，空穴長(zhǎng)度達(dá)到最大值，空穴末端達(dá)到翼型尾部，這時(shí)空穴末端近壁面區(qū)域的水汽混合區(qū)出現(xiàn)回縮，在t=t0+32.2 ms時(shí)，翼型尾部出現(xiàn)大規(guī)模的空泡團(tuán)旋渦脫落現(xiàn)象。

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

3.2 繞翼型速度矢量場(chǎng)

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

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

圖5 典型時(shí)刻繞翼型尾部速度場(chǎng)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)過(guò)程。指出云狀空化后期，在緊貼壁面的區(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為空穴開(kāi)始回縮，反向射流逐漸向翼型前緣發(fā)展的時(shí)刻。采用不同空化模型所描述的空穴內(nèi)部的流場(chǎng)結(jié)構(gòu)是相同的：在空泡回縮過(guò)程中，采用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)的趨勢(shì)，使空泡團(tuán)逐漸脫離壁面，符合實(shí)驗(yàn)觀測(cè)現(xiàn)象。為了進(jìn)一步闡明流動(dòng)結(jié)構(gòu)的發(fā)展，圖5給出了實(shí)驗(yàn)與數(shù)值計(jì)算得到的翼型尾部局部的特寫圖像，通過(guò)對(duì)比實(shí)驗(yàn)現(xiàn)象可以發(fā)現(xiàn)：采用Singhal空化模型準(zhǔn)確地模擬出了翼型尾部的空化旋渦區(qū)，空泡的旋渦脫落現(xiàn)象非常明顯。

3.3 繞翼型旋渦特性

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

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

渦量來(lái)源于流場(chǎng)存在速度梯度，是描述有旋流動(dòng)的一個(gè)運(yùn)動(dòng)學(xué)物理量，在研究云狀空化現(xiàn)象中，涉及到空泡團(tuán)脫落等分離流動(dòng)情況，對(duì)渦量進(jìn)行分析具有很重要的意義，分別采用Kubota與Singhal空化模型對(duì)繞翼型渦量場(chǎng)進(jìn)行描述，利用相關(guān)后處理軟件得到了對(duì)應(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)在對(duì)空化旋渦區(qū)的模擬上，如圖所示，點(diǎn)劃線a與實(shí)線b與翼型壁面所圍成的區(qū)域分別為采用Singhal與Kubota空化模型計(jì)算得到的空化漩渦區(qū)，采用Singhal空化模型計(jì)算得到的空化旋渦區(qū)域明顯比采用Kubota時(shí)的大，也就是說(shuō)，采用Singhal空化模型計(jì)算結(jié)果體現(xiàn)了更為明顯的旋渦特性。

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

不同空化模型在模擬繞翼型云狀空化速度與渦量場(chǎng)上的差異反映了不同空化模型所體現(xiàn)的物理機(jī)制是不同的，Singhal空化模型在推導(dǎo)相間質(zhì)量傳輸速率過(guò)程中，添加了瞬時(shí)粘性力與湍流波動(dòng)的影響，綜合考慮了汽、液兩相間的相對(duì)速度、空泡表面張力等諸多因素，從而可更好地考慮空化現(xiàn)象中非常頻繁的質(zhì)量與速度交換，而Kubota空化模型則忽略了空化流動(dòng)相變過(guò)程中頻繁的動(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í)均速度、渦量場(chǎng)，并與試驗(yàn)結(jié)果進(jìn)行了對(duì)比分析，結(jié)論如下：

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

（2）Kubota空化模型模擬出來(lái)的空化旋渦區(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]王國(guó)玉,曹樹(shù)良,井小萩利明.剪切層區(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]蔣增輝,于開(kāi)平,張嘉鐘等.空泡尾部流場(chǎng)特性數(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]張 博,王國(guó)玉等.修正的RNG k-ε模型在云狀空化流動(dòng)計(jì)算中的應(yīng)用評(pí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]馮 光,顏 開(kāi).超空泡航行體水下彈道的數(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.