基于傳質(zhì)源項(xiàng)的水翼空蝕風(fēng)險(xiǎn)數(shù)值預(yù)報(bào)

鄭恩慧，曹彥濤，程懷玉，季 斌，彭曉星

（1.中國船舶科學(xué)研究中心船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室，江蘇無錫 214082；2.武漢大學(xué)水利水電學(xué)院，武漢 430072）

0 引 言

隨著現(xiàn)代艦船噸位和航速的不斷提升，空化現(xiàn)象也愈發(fā)常見?？栈鶗?huì)顯著增大艦船推進(jìn)器的振動(dòng)、噪聲，引起性能下降，嚴(yán)重時(shí)甚至?xí)苯釉斐赏七M(jìn)器葉片表面材料的剝蝕，引起空蝕，嚴(yán)重危害艦船的安全。在推進(jìn)器的設(shè)計(jì)過程中，需要在保證推進(jìn)效率的同時(shí)，盡可能減小或避免空蝕的影響。因此需要一個(gè)準(zhǔn)確客觀的空蝕風(fēng)險(xiǎn)預(yù)報(bào)方法，以獲取空蝕可能出現(xiàn)的區(qū)域和強(qiáng)度，從而為設(shè)計(jì)方案的優(yōu)化提供依據(jù)。

傳統(tǒng)推進(jìn)器設(shè)計(jì)中主要依靠模型涂層試驗(yàn)來判斷空蝕風(fēng)險(xiǎn)的潛在位置及強(qiáng)度，即通過在模型表面噴涂軟面涂層開展設(shè)計(jì)工況下的空化試驗(yàn)，利用涂層剝離區(qū)的狀況判斷空蝕風(fēng)險(xiǎn)。Cao 等[1]利用涂層試驗(yàn)得到了兩種工況下NACA 0009 三維扭曲水翼的空蝕風(fēng)險(xiǎn)區(qū)域；Li 等[2]分別對(duì)8°攻角NACA 0015 水翼及6.5°攻角NACA 0018-45 水翼進(jìn)行設(shè)計(jì)工況下的涂層試驗(yàn)。但采用試驗(yàn)方法研究空蝕風(fēng)險(xiǎn)問題通常面臨著試驗(yàn)周期長、費(fèi)用高、數(shù)據(jù)有限等問題。近年來隨著CFD 的不斷發(fā)展，空化流動(dòng)數(shù)值模擬技術(shù)日趨成熟，其在推進(jìn)器空化性能的設(shè)計(jì)評(píng)估中得到了廣泛的應(yīng)用[3-7]。在可靠的空化流場數(shù)值模擬結(jié)果基礎(chǔ)上，對(duì)空蝕風(fēng)險(xiǎn)進(jìn)行數(shù)值預(yù)報(bào)也成為了可能。目前，研究者們通常以空化流場壓力、蒸汽體積分?jǐn)?shù)等特征參數(shù)為依托，通過一些表征參數(shù)構(gòu)建空蝕風(fēng)險(xiǎn)預(yù)報(bào)的方法。這些方法近年來獲得了較大發(fā)展，在空化流場空蝕風(fēng)險(xiǎn)的評(píng)估中得到了一定的應(yīng)用。Dular[8]將蒸汽體積分?jǐn)?shù)的標(biāo)準(zhǔn)差與空蝕風(fēng)險(xiǎn)聯(lián)系起來；Nohmi[9]通過觀察數(shù)值模擬中的汽泡行為，簡單定義了固壁壓力及蒸汽體積分?jǐn)?shù)的函數(shù)來衡量空蝕風(fēng)險(xiǎn)；Ochiai[10]提出了離散泡法，其空蝕風(fēng)險(xiǎn)參數(shù)是通過R-P 方程計(jì)算微觀汽泡內(nèi)部的壓力來估計(jì)的。

準(zhǔn)確可靠的數(shù)值模擬結(jié)果是開展空蝕風(fēng)險(xiǎn)數(shù)值預(yù)報(bào)的基礎(chǔ)。對(duì)于空化流動(dòng)的數(shù)值模擬而言，當(dāng)前基于均質(zhì)平衡流模型框架的質(zhì)量輸運(yùn)方程空化模型得到了廣泛應(yīng)用[11-15]。結(jié)合改進(jìn)型RANS、DES或LES等湍流模擬方法[16-20]，可以較好地模擬非定常空化流動(dòng)的主要特征。

而在空蝕風(fēng)險(xiǎn)表征參數(shù)的構(gòu)建方面，當(dāng)前存在兩種主要研究思路：一是精細(xì)化捕捉空化泡群潰滅特征，利用歐拉-拉格朗日方法精細(xì)化捕捉空化潰滅的瞬態(tài)載荷特征，但該類方法所需計(jì)算資源巨大，當(dāng)前尚未達(dá)到工程實(shí)用化的程度[21-25]；另一種思路是拋開復(fù)雜的流動(dòng)細(xì)節(jié)，基于宏觀流場的參數(shù)構(gòu)建能夠反映空蝕風(fēng)險(xiǎn)特性的表征方法。其中，以勢(shì)能假設(shè)為基礎(chǔ)的空蝕風(fēng)險(xiǎn)表征參數(shù)方法已經(jīng)獲得較大發(fā)展。Patella[26]提出了空泡潰滅后沖擊波的能量傳遞機(jī)制；Li[2]提出了基于壓力時(shí)間導(dǎo)數(shù)的空蝕風(fēng)險(xiǎn)表征參數(shù)；Leclercq[27]利用局部勢(shì)能變化率結(jié)合立體角模型計(jì)算了壁面吸收到的輻射能量；Melissar?is[28]假定空泡結(jié)構(gòu)的初始勢(shì)能在潰滅前首先轉(zhuǎn)化為空泡界面處的動(dòng)能，隨著空泡體積的減小，動(dòng)能逐漸向空泡中心聚焦，一旦達(dá)到潰滅條件，全部能量以壓力波的形式向外釋放，壓力波以投影的方式輻射到壁面上，從而獲得壁面上的空蝕風(fēng)險(xiǎn)區(qū)域；Mohammad[29]則直接計(jì)算空泡潰滅過程中周圍液體的動(dòng)能，避免了勢(shì)能計(jì)算過程中驅(qū)動(dòng)壓強(qiáng)的不確定性問題，并且同時(shí)考慮了沖擊波與微射流兩種空蝕作用機(jī)制。

上述研究加深了人們對(duì)于空蝕風(fēng)險(xiǎn)預(yù)報(bào)理論與方法的認(rèn)識(shí)，為本文的研究提供了非常有意義的參考。本文基于單泡潰滅模型，分析影響空泡潰滅的因素，在當(dāng)?shù)貕簭?qiáng)大于飽和蒸汽壓及壓強(qiáng)變化率大于某個(gè)閾值的前提下，利用傳質(zhì)源項(xiàng)計(jì)算一段時(shí)間內(nèi)由于空泡潰滅引起的壁面上的能量累積，通過為壁面累積能量劃定閾值確定空蝕風(fēng)險(xiǎn)區(qū)域的范圍，通過壁面累積能量的大小判定空蝕風(fēng)險(xiǎn)的相對(duì)嚴(yán)重程度。通過該方法預(yù)報(bào)NACA 0009三維扭曲水翼的空蝕風(fēng)險(xiǎn)結(jié)果，并與試驗(yàn)結(jié)果進(jìn)行比較，驗(yàn)證該方法的準(zhǔn)確性。

1 水翼空化演化過程數(shù)值模擬

1.1 控制方程與大渦模擬方法

本文空化流動(dòng)計(jì)算采用的兩相流模型為mixture 模型。控制方程為連續(xù)性方程和動(dòng)量方程，以張量形式表述如下：

式中，ui是混合物在i方向的速度，p是汽液混合物的壓力?；旌衔锏恼承驭膛c密度ρ視為蒸汽體積分?jǐn)?shù)αv的函數(shù)：

其中，下標(biāo)v和l分別代表蒸汽相及液體相。

為了能較好地模擬各種湍流尺度及更好地捕捉渦旋結(jié)構(gòu)，本文采用LES 方法模擬湍流。大渦模擬的控制方程為

式中，應(yīng)力項(xiàng)τij用亞格子模型求解。通過WALE 模型模擬亞格子尺度湍流特性。該模型在此類空化流動(dòng)中已經(jīng)得到了較為廣泛的應(yīng)用，取得了令人滿意的預(yù)報(bào)結(jié)果[16,18]

1.2 Schnerr-Sauer空化模型

本文采用的空化模型為Schnerr-Sauer(S-S)空化模型[15]，汽、液兩相間的質(zhì)量傳輸由蒸汽體積分?jǐn)?shù)輸運(yùn)方程來控制：

式中：αv表示蒸汽體積分?jǐn)?shù)；ρl為液體密度；ρv為蒸汽密度；m?為汽液兩相間的凈質(zhì)量輸運(yùn)源項(xiàng)，又分為蒸汽生成項(xiàng)m?+及凝結(jié)項(xiàng)m?-，表征如下：

其中：pd為泡外驅(qū)動(dòng)壓力；pv為飽和蒸汽壓；n為汽泡數(shù)密度，本文取n=1×1013。

1.3 計(jì)算域及網(wǎng)格劃分

本文采用的水翼模型為三維扭曲水翼，其截面翼型為NACA 0009。弦長為112.5 mm，展長為225 mm，左右對(duì)稱，水翼中部攻角為11°，兩端攻角為0°，水翼模型的外形如圖1 所示。圖2 給出了數(shù)值模擬中使用的計(jì)算域及相應(yīng)的邊界條件的設(shè)置，C為水翼弦長，邊界條件采用速度進(jìn)口與壓力出口，進(jìn)口速度為14 m/s，通過調(diào)節(jié)出口壓力值使空化數(shù)等于1.2，空化數(shù)定義為

圖1 水翼模型Fig.1 Hydrofoil model

圖2 計(jì)算域及邊界條件Fig.2 Computational domain and boundary conditions

需要注意的是，為了節(jié)省計(jì)算資源，本文在數(shù)值模擬中僅對(duì)該水翼的一半進(jìn)行了模擬計(jì)算，中間切面選取對(duì)稱面邊界條件，其余壁面均選取無滑移壁面條件。

圖3 給出了計(jì)算域內(nèi)的網(wǎng)格劃分情況。為保證網(wǎng)格質(zhì)量，全域均采用結(jié)構(gòu)化網(wǎng)格。為精細(xì)地捕捉水翼吸力面的空泡結(jié)構(gòu)演變特性，在水翼吸力面附近進(jìn)行網(wǎng)格細(xì)化，在水翼壁面附近添加邊界層，使得壁面函數(shù)y+在1 左右。為兼顧計(jì)算效率與計(jì)算精度，參考Long 等[30]的工作，本文的總網(wǎng)格數(shù)取1100萬，非定常計(jì)算時(shí)間步長選取為Δt=1×10-5s。

圖3 網(wǎng)格劃分情況Fig.3 Mesh generation

1.4 計(jì)算結(jié)果與對(duì)比分析

以數(shù)值模擬手段對(duì)空蝕風(fēng)險(xiǎn)預(yù)報(bào)的前提是該模擬獲得的空化流場形態(tài)應(yīng)當(dāng)與實(shí)際一致。為了驗(yàn)證空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果的準(zhǔn)確性，本文在涂層試驗(yàn)過程中使用高速攝影拍攝水翼吸力面的空化形態(tài)。試驗(yàn)工況與數(shù)值模擬一致。

圖4 為一個(gè)脫落周期（T0）內(nèi)水翼吸力面試驗(yàn)空化形態(tài)與數(shù)值模擬蒸汽體積分?jǐn)?shù)αv=0.1 的等值面圖像的對(duì)比。結(jié)果顯示，試驗(yàn)和數(shù)值結(jié)果較一致地展現(xiàn)了片空化的產(chǎn)生、發(fā)展、脫落及潰滅的準(zhǔn)周期過程。該工況下空化形態(tài)主要特點(diǎn)為：水翼兩側(cè)流動(dòng)狀態(tài)相對(duì)穩(wěn)定，片空化的脫落、再生過程集中在水翼中央；片空化充分發(fā)展時(shí)被向上游發(fā)展的回射流從導(dǎo)邊處切斷形成脫落云，見圖4（a）～（b）；隨后脫落云中間部分向上卷起，兩端附著在水翼表面，逐漸形成U型渦，見圖4（c）～（d）；脫落云與U 型渦隨主流向下游運(yùn)動(dòng)過程中，伴隨著大量空泡潰滅，直到U 型渦離開隨邊時(shí)，空泡幾乎全部潰滅，見圖4（e）～（f）。

圖4 水翼吸力面試驗(yàn)中的空化形態(tài)與數(shù)值模擬αv=0.1 空化形態(tài)Fig.4 Cavitation morphology captured by test and iso-surface of αv=0.1 captured by numerical simulation on hydrofoil suction surface

2 空蝕風(fēng)險(xiǎn)涂層試驗(yàn)

2.1 試驗(yàn)設(shè)備及試驗(yàn)方法

本次試驗(yàn)主要設(shè)備包括循環(huán)式高速空泡水筒、高速攝影相機(jī)等。圖5為高速空泡水筒的示意圖。高速空泡水筒試驗(yàn)段的長度為1600 mm，寬和高均為225 mm，最高水速為25 m/s，試驗(yàn)段壓力可調(diào)整為5 kPa至500 kPa，邊壁為透明有機(jī)玻璃，以便高速攝影拍攝試驗(yàn)情況。高速攝影采用3200 fps的拍攝頻率，對(duì)水翼的吸力面進(jìn)行拍攝。

圖5 高速空泡水筒示意圖Fig.5 Schematic diagram of high speed cavitation tank

水翼模型參數(shù)與數(shù)值模擬模型完全相同，鋁合金材質(zhì)，在水翼表面噴涂軟面材料，便于更快地觀測空蝕風(fēng)險(xiǎn)情況，試驗(yàn)前水翼吸力面如圖6所示。本次水翼安裝采用0°攻角，吸力面朝下，在試驗(yàn)段的安裝位置如圖7所示。

圖6 試驗(yàn)前水翼吸力面Fig.6 Hydrofoil suction surface before test

圖7 水翼安裝位置Fig.7 Installation position of hydrofoil

試驗(yàn)水速為14 m/s，試驗(yàn)段空化數(shù)為1.2。涂層試驗(yàn)始終在該工況下進(jìn)行，試驗(yàn)開始后，空化形態(tài)很快趨于穩(wěn)定。試驗(yàn)過程中，通過高速攝影判斷水翼空蝕風(fēng)險(xiǎn)的大致情況，當(dāng)發(fā)現(xiàn)高速攝影中出現(xiàn)明顯的空蝕風(fēng)險(xiǎn)區(qū)域時(shí)，停止運(yùn)行空泡水筒，再在無空化情況下確認(rèn)空蝕風(fēng)險(xiǎn)區(qū)域是否形成。本次試驗(yàn)進(jìn)行1 h后，發(fā)現(xiàn)水翼吸力面出現(xiàn)了大量的空蝕風(fēng)險(xiǎn)區(qū)域，隨即停止試驗(yàn)，將水翼從高速水洞中取出，觀察并記錄水翼表面空蝕風(fēng)險(xiǎn)情況。

2.2 涂層試驗(yàn)空蝕風(fēng)險(xiǎn)結(jié)果及分析

空蝕風(fēng)險(xiǎn)試驗(yàn)結(jié)果如圖8 所示。其中，黑色區(qū)域?yàn)橥繉硬牧鲜艿絼兾g的區(qū)域。圖8（a）中標(biāo)出了空蝕風(fēng)險(xiǎn)區(qū)域，圖9則給出了相應(yīng)區(qū)域的高速攝影典型空化形態(tài)。

圖8 涂層試驗(yàn)空蝕風(fēng)險(xiǎn)結(jié)果Fig.8 Results of cavitation erosion risk

圖9 試驗(yàn)空蝕風(fēng)險(xiǎn)結(jié)果與空化形態(tài)的對(duì)應(yīng)Fig.9 Correspondence between cavitation erosion risk results and cavitation morphology

涂層試驗(yàn)結(jié)果顯示，空蝕風(fēng)險(xiǎn)區(qū)域主要分為兩部分：主空蝕風(fēng)險(xiǎn)區(qū)域及條帶狀空蝕風(fēng)險(xiǎn)區(qū)域。其中主空蝕風(fēng)險(xiǎn)區(qū)域大約處于水翼弦向2/11弦長到8/11弦長之間，其空蝕風(fēng)險(xiǎn)最為嚴(yán)重，對(duì)應(yīng)于片空化脫落區(qū)域。與主空蝕風(fēng)險(xiǎn)區(qū)域相連接的兩個(gè)條帶狀的空蝕風(fēng)險(xiǎn)區(qū)域一直延伸到隨邊，存在明顯的點(diǎn)狀空蝕風(fēng)險(xiǎn)區(qū)域，對(duì)應(yīng)于U 型渦兩腿經(jīng)過的區(qū)域。此外，U 型渦脫離開隨邊的位置處，空蝕風(fēng)險(xiǎn)明顯加重，對(duì)應(yīng)于條帶狀空蝕風(fēng)險(xiǎn)區(qū)域與隨邊相連的位置，如圖9（b）中紅色箭頭指向的位置。

3 空蝕風(fēng)險(xiǎn)預(yù)報(bào)方法

3.1 空蝕機(jī)理回顧

空泡向水翼下游移動(dòng)進(jìn)入壓力恢復(fù)區(qū)域時(shí)，由于泡外壓力不斷增大，空泡開始潰滅。人們普遍認(rèn)為，空蝕是蒸汽空泡潰滅引起的。對(duì)于初始半徑為a、初始泡壁速度為u0=0 的蒸汽空泡，不考慮表面張力、粘性、不可凝氣體在內(nèi)的潰滅動(dòng)力學(xué)方程，即R-P方程為

式中，r0為汽泡半徑，pd(t)為泡外壓強(qiáng)，pv為飽和蒸汽壓，ρl為液體密度。式（10）的解，即潰滅過程泡壁速度為

空泡收縮至潰滅的時(shí)間為

顯然，當(dāng)空泡半徑r0→0，泡壁面流體沿徑向速度u0→∞。實(shí)際流動(dòng)中，由于粘性泡內(nèi)存在不可凝結(jié)氣體，使得空泡不能潰滅至半徑為0，泡壁面速度也不能達(dá)到無窮大。盡管如此，空泡潰滅至最小半徑時(shí)，其泡壁速度也將是很大的數(shù)值，這實(shí)際上就是造成空蝕的微沖擊波及微射流的來源。

當(dāng)然，只有非?？拷诿娴目张轁绮艑?duì)空蝕有貢獻(xiàn)，否則微沖擊波與微射流在尚未到達(dá)壁面時(shí)將衰減耗盡。此外，無論是球形繞流體近壁面的順流而下，還是球形泡的收縮都會(huì)產(chǎn)生趨壁效應(yīng)[31]。定性而言，趨壁效應(yīng)使得近壁面?zhèn)葔毫M(jìn)一步降低，原本中心對(duì)稱的沖擊波在靠近壁面方向的強(qiáng)度會(huì)大于其他方向，而微射流的方向也會(huì)指向壁面，使得向壁面方向傳遞的能量更加集中。

所以，決定空蝕風(fēng)險(xiǎn)的主要是近壁面的流場，本文提出的空蝕風(fēng)險(xiǎn)參數(shù)，只需在貼壁流場中去判別與尋找。Phillip 和Lauterborn[32]也通過試驗(yàn)證明了最大的沖擊載荷是由與壁面直接接觸的空泡潰滅造成的。

3.2 蒸汽泡潰滅的能量轉(zhuǎn)化過程

根據(jù)前述，本文兩相流的數(shù)值模擬方法與三維扭曲水翼空化流試驗(yàn)結(jié)果吻合較好。但該兩相流數(shù)值模擬無法直接顯示蒸汽泡潰滅引起的徑向流場速度u0的時(shí)空分布，更無法表征微沖擊波與微射流。為此需要尋找其它表征潰滅的物理參數(shù)。

3.2.1 蒸汽泡勢(shì)能

將半徑為r0的單個(gè)蒸汽球泡與周圍半徑為R（R→∞）的球狀液體空間組成一個(gè)系統(tǒng)，如圖10所示。由于系統(tǒng)蒸汽泡內(nèi)部氣壓為飽和蒸汽壓Pv，而系統(tǒng)外則為大于飽和蒸汽壓的當(dāng)?shù)貕毫d(t)，因此，在t時(shí)刻系統(tǒng)內(nèi)外存在壓差Pd(t)-Pv。

圖10 單空泡潰滅模型Fig.10 Single bubble collapse model

收縮過程中，蒸汽球泡泡徑為r0時(shí)，若泡徑變化dr0，pv對(duì)液體作負(fù)功為

泡外充滿水的充分大半徑為R的同心球體因連續(xù)性方程同時(shí)需縮小dR，有

pd(t)對(duì)液體作正功為

故泡收縮時(shí)，外力對(duì)液體實(shí)際作功為

式中，Vv為蒸汽泡的體積。

因此，當(dāng)蒸汽泡收縮dVv時(shí)，外力對(duì)泡外液體作功dW；而蒸汽泡損失的勢(shì)能為

蒸汽泡潰滅釋放出全部勢(shì)能

3.2.2 潰滅過程中蒸汽泡的勢(shì)能轉(zhuǎn)化為流場的動(dòng)能

初始泡壁速度為0的空泡，當(dāng)半徑由a縮減至任意值r*≤a時(shí)，空泡的勢(shì)能減小為

式中，u*為泡壁r*處的速度。

當(dāng)空泡半徑由a縮減至r*時(shí)，流場動(dòng)能增加為

則ΔEp=-ΔEk。可知蒸汽泡減小的勢(shì)能全部轉(zhuǎn)化為周圍流場的動(dòng)能。所以，尋求蒸汽泡潰滅所引起的流場動(dòng)能的時(shí)空分布，可以等價(jià)轉(zhuǎn)換為尋求流場蒸汽泡勢(shì)能減少的時(shí)空分布。

3.3 空蝕風(fēng)險(xiǎn)表征

本文采用的S-S 空化模型中選用單位體積液體氣泡數(shù)密度為n=1×1013，空蝕風(fēng)險(xiǎn)是群泡的物理行為，就是對(duì)單位體積中存在的1×1013個(gè)蒸汽泡的時(shí)空特性的描述。

3.3.1 空蝕風(fēng)險(xiǎn)表征參數(shù)的選取

（1）環(huán)境壓強(qiáng)

式（10）可改寫為

（2）環(huán)境壓強(qiáng)變化率

（3）壁面累積能量

由于目前數(shù)值模擬尚不能給出空泡群潰滅過程準(zhǔn)確定量的描述，且空蝕風(fēng)險(xiǎn)預(yù)報(bào)無需定量確定壁面用于剝蝕能量的大小，只需要一個(gè)相對(duì)值。但有一點(diǎn)可以肯定的是，近壁面流場釋放的空泡潰滅的能量越大越多，壁面累積的用于空蝕的能量便越多，空蝕就越容易發(fā)生，空蝕程度就越嚴(yán)重。

式（18）給出了單泡潰滅時(shí)勢(shì)能釋放的表達(dá)式，則對(duì)于群泡，在單位體積中釋放的勢(shì)能變化率為

式中，αv為單位體積的蒸汽體積分?jǐn)?shù)，ep為單位體積泡群釋放的勢(shì)能。

因?yàn)榭瘴g是微沖擊波與微射流長期作用在壁面上的結(jié)果，參數(shù)一、二的滿足可視為單次沖擊具備了引起空蝕風(fēng)險(xiǎn)的潛力，但需要多次類似的沖擊疊加才可能導(dǎo)致空蝕?？紤]到空化流的非定常特性，故選取一充分長的時(shí)間T，將壁面累積的空泡潰滅過程中釋放的勢(shì)能值超過閾值C2作為參數(shù)三。

3.3.2 空蝕風(fēng)險(xiǎn)表征方法

（3）壁面空蝕風(fēng)險(xiǎn)表征參數(shù)三為Ic.e＞C2。

根據(jù)前人的研究，本文在計(jì)算壁面累積能量Ic.e時(shí)采用了時(shí)間平均壓力場pˉd(t)作為驅(qū)動(dòng)空泡潰滅的環(huán)境壓力場[28]。

需要說明的是，式（26）中蒸汽體積分?jǐn)?shù)變化率項(xiàng)若使用物質(zhì)導(dǎo)數(shù)的展開式進(jìn)行計(jì)算，其時(shí)間偏導(dǎo)數(shù)項(xiàng)及對(duì)流項(xiàng)均需采用差分格式，若無法保障數(shù)值計(jì)算結(jié)果的時(shí)間步長足夠小及網(wǎng)格分辨率足夠低，將會(huì)對(duì)計(jì)算精度產(chǎn)生較大影響。本文采用傳質(zhì)源項(xiàng)對(duì)蒸汽體積分?jǐn)?shù)變化率進(jìn)行替代，式（6）將蒸汽體積變化與相變過程中的傳質(zhì)源項(xiàng)聯(lián)系起來，傳質(zhì)源項(xiàng)的表達(dá)式見式（7），可發(fā)現(xiàn)其全部變量數(shù)據(jù)取自當(dāng)前時(shí)間步，且無需求解對(duì)流項(xiàng)，因此解決了差分格式帶來的誤差問題，相對(duì)而言，精度更高。根據(jù)式（6）～（7），潰滅過程中的蒸汽體積分?jǐn)?shù)變化率可表示為

在滿足當(dāng)?shù)貕簭?qiáng)大于飽和蒸汽壓，壓強(qiáng)隨時(shí)間變化率大于閾值C1的條件后，計(jì)算壁面累積能量Ic.e，壁面累積能量大于C2可視為該區(qū)域存在空蝕風(fēng)險(xiǎn)。

3.4 空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果

3.4.1 閾值C1、C2的選定

閾值C1、C2作為壓強(qiáng)隨時(shí)間變化率及壁面累積能量的下限，共同決定了空蝕風(fēng)險(xiǎn)區(qū)域的范圍。在此基礎(chǔ)上，壁面累積能量Ic.e的相對(duì)大小則可表征空蝕風(fēng)險(xiǎn)的嚴(yán)重程度。

分別選取閾值C1=2×105、1×106、5×106及C2=1×1010、2×1010、3×1010，將閾值C1、C2兩兩配對(duì)計(jì)算水翼表面的空蝕風(fēng)險(xiǎn)。由此得到9個(gè)空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果，并將其與試驗(yàn)空蝕風(fēng)險(xiǎn)區(qū)域進(jìn)行比較，如圖11所示。可以看出，隨著閾值C1的增加，水翼中前部的空蝕風(fēng)險(xiǎn)區(qū)域逐漸縮小，中后部空蝕風(fēng)險(xiǎn)越來越明顯；隨著閾值C2的增加，空蝕風(fēng)險(xiǎn)區(qū)域的外圍輪廓，尤其是水翼導(dǎo)邊附近區(qū)域隨之縮減。通過比較，與試驗(yàn)結(jié)果吻合最好的為圖11（e）。因此選定閾值C1=1×106，C2=2×1010。

圖11 9組閾值計(jì)算的空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果與試驗(yàn)結(jié)果比較Fig.11 Prediction results of cavitation erosion risk calculated by 9 groups of threshold values compared with the experimental results

3.4.2 空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果

圖12 為采用基于傳質(zhì)源項(xiàng)的空蝕風(fēng)險(xiǎn)預(yù)報(bào)方法得到的空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果與試驗(yàn)結(jié)果的對(duì)比。可以看出，空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果在片空化脫落的區(qū)域預(yù)報(bào)出了空蝕風(fēng)險(xiǎn)，同時(shí)在U型渦經(jīng)過的位置也預(yù)報(bào)出了兩個(gè)條帶狀的空蝕風(fēng)險(xiǎn)區(qū)域，與試驗(yàn)結(jié)果顯示的空蝕風(fēng)險(xiǎn)區(qū)域基本吻合。此外，預(yù)報(bào)結(jié)果顯示片空化脫落區(qū)域的空蝕風(fēng)險(xiǎn)程度更為嚴(yán)重，與試驗(yàn)結(jié)果亦較為吻合。盡管預(yù)報(bào)結(jié)果與試驗(yàn)結(jié)果在大部分區(qū)域吻合較好，但預(yù)報(bào)結(jié)果顯示在水翼兩端即紅框范圍內(nèi)存在空蝕風(fēng)險(xiǎn)區(qū)域，與試驗(yàn)結(jié)果不同。數(shù)值模擬和試驗(yàn)顯示在水翼兩端的空化形態(tài)為片空化的周期性波動(dòng)，與水翼中央出現(xiàn)的云空化不同，片空化周期性變化引起的空蝕風(fēng)險(xiǎn)預(yù)報(bào)方法需進(jìn)一步研究。

圖12 空蝕風(fēng)險(xiǎn)預(yù)報(bào)結(jié)果與試驗(yàn)結(jié)果的對(duì)比Fig.12 Comparison between cavitation erosion risk prediction results and painting test results

4 結(jié) 論

本文以NACA 0009 三維扭曲水翼為數(shù)值模擬及涂層試驗(yàn)對(duì)象，將數(shù)值模擬空化形態(tài)與試驗(yàn)高速攝影結(jié)果進(jìn)行了對(duì)比。基于單個(gè)蒸汽球泡的潰滅模型，分析了潰滅過程中能量轉(zhuǎn)化過程及相質(zhì)量的轉(zhuǎn)化過程，提出了空蝕風(fēng)險(xiǎn)的三個(gè)表征參數(shù)，并合理選定了表征參數(shù)的閾值。基于該方法，本文對(duì)繞NACA 0009 三維扭曲水翼表面的空蝕風(fēng)險(xiǎn)進(jìn)行了預(yù)報(bào)，并與涂層試驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。主要結(jié)論如下：

（1）采用大渦模擬方法可較好地捕捉繞NACA 0009 三維扭曲水翼空化流動(dòng)的初生、發(fā)展、脫落、潰滅等過程，與試驗(yàn)結(jié)果吻合較好；

（2）提出空蝕風(fēng)險(xiǎn)預(yù)報(bào)的三個(gè)表征參數(shù)：當(dāng)?shù)貕簭?qiáng)、當(dāng)?shù)貕簭?qiáng)變化率及壁面累積能量Ic.e。在當(dāng)?shù)貕簭?qiáng)大于飽和蒸汽壓及當(dāng)?shù)貕簭?qiáng)變化率大于閾值C1的前提下，通過壁面累積能量Ic.e大于閾值C2可以劃定空蝕風(fēng)險(xiǎn)的區(qū)域，通過Ic.e的大小可以判定區(qū)域內(nèi)空蝕風(fēng)險(xiǎn)的相對(duì)嚴(yán)重程度。其中，Ic.e的求解采用了基于傳質(zhì)源項(xiàng)的表達(dá)，避免了數(shù)值誤差。

（3）通過與涂層試驗(yàn)比較，本文建立的空蝕風(fēng)險(xiǎn)數(shù)值預(yù)報(bào)方法獲得的空蝕風(fēng)險(xiǎn)結(jié)果與試驗(yàn)空蝕風(fēng)險(xiǎn)結(jié)果在大部分區(qū)域吻合良好。兩者顯示的空蝕風(fēng)險(xiǎn)區(qū)域主要分為兩部分，一是位于片空化脫落的區(qū)域，二是位于U型渦兩腿經(jīng)過的區(qū)域，其中片空化脫落區(qū)域空蝕風(fēng)險(xiǎn)更為嚴(yán)重。

雖然本文空蝕風(fēng)險(xiǎn)預(yù)報(bào)僅對(duì)三維扭曲水翼進(jìn)行了驗(yàn)證，而且閾值C1、C2的確定帶有經(jīng)驗(yàn)性質(zhì)，但因?yàn)槿齻€(gè)空蝕風(fēng)險(xiǎn)參數(shù)均建立于空蝕機(jī)理之上，所以具有一定的普適性。下一步將研究如何將閾值C1、C2與影響空蝕的物理量聯(lián)系起來，從理論上確定閾值C1、C2。