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

（1.中國船舶科學研究中心船舶振動噪聲重點實驗室，江蘇無錫 214082；2.武漢大學水利水電學院，武漢 430072）

0 引 言

隨著現(xiàn)代艦船噸位和航速的不斷提升，空化現(xiàn)象也愈發(fā)常見?？栈鶗@著增大艦船推進器的振動、噪聲，引起性能下降，嚴重時甚至會直接造成推進器葉片表面材料的剝蝕，引起空蝕，嚴重危害艦船的安全。在推進器的設計過程中，需要在保證推進效率的同時，盡可能減小或避免空蝕的影響。因此需要一個準確客觀的空蝕風險預報方法，以獲取空蝕可能出現(xiàn)的區(qū)域和強度，從而為設計方案的優(yōu)化提供依據(jù)。

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

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

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

上述研究加深了人們對于空蝕風險預報理論與方法的認識，為本文的研究提供了非常有意義的參考。本文基于單泡潰滅模型，分析影響空泡潰滅的因素，在當?shù)貕簭姶笥陲柡驼羝麎杭皦簭娮兓蚀笥谀硞€閾值的前提下，利用傳質(zhì)源項計算一段時間內(nèi)由于空泡潰滅引起的壁面上的能量累積，通過為壁面累積能量劃定閾值確定空蝕風險區(qū)域的范圍，通過壁面累積能量的大小判定空蝕風險的相對嚴重程度。通過該方法預報NACA 0009三維扭曲水翼的空蝕風險結(jié)果，并與試驗結(jié)果進行比較，驗證該方法的準確性。

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

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

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

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

其中，下標v和l分別代表蒸汽相及液體相。

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

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

1.2 Schnerr-Sauer空化模型

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

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

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

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

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

圖1 水翼模型Fig.1 Hydrofoil model

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

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

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

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

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

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

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

圖4 水翼吸力面試驗中的空化形態(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 空蝕風險涂層試驗

2.1 試驗設備及試驗方法

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

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

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

圖6 試驗前水翼吸力面Fig.6 Hydrofoil suction surface before test

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

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

2.2 涂層試驗空蝕風險結(jié)果及分析

空蝕風險試驗結(jié)果如圖8 所示。其中，黑色區(qū)域為涂層材料受到剝蝕的區(qū)域。圖8（a）中標出了空蝕風險區(qū)域，圖9則給出了相應區(qū)域的高速攝影典型空化形態(tài)。

圖8 涂層試驗空蝕風險結(jié)果Fig.8 Results of cavitation erosion risk

圖9 試驗空蝕風險結(jié)果與空化形態(tài)的對應Fig.9 Correspondence between cavitation erosion risk results and cavitation morphology

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

3 空蝕風險預報方法

3.1 空蝕機理回顧

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

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

空泡收縮至潰滅的時間為

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

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

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

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

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

3.2.1 蒸汽泡勢能

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

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

收縮過程中，蒸汽球泡泡徑為r0時，若泡徑變化dr0，pv對液體作負功為

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

pd(t)對液體作正功為

故泡收縮時，外力對液體實際作功為

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

因此，當蒸汽泡收縮dVv時，外力對泡外液體作功dW；而蒸汽泡損失的勢能為

蒸汽泡潰滅釋放出全部勢能

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

初始泡壁速度為0的空泡，當半徑由a縮減至任意值r*≤a時，空泡的勢能減小為

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

當空泡半徑由a縮減至r*時，流場動能增加為

則ΔEp=-ΔEk?？芍羝轀p小的勢能全部轉(zhuǎn)化為周圍流場的動能。所以，尋求蒸汽泡潰滅所引起的流場動能的時空分布，可以等價轉(zhuǎn)換為尋求流場蒸汽泡勢能減少的時空分布。

3.3 空蝕風險表征

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

3.3.1 空蝕風險表征參數(shù)的選取

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

式（10）可改寫為

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

（3）壁面累積能量

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

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

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

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

3.3.2 空蝕風險表征方法

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

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

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

在滿足當?shù)貕簭姶笥陲柡驼羝麎?，壓強隨時間變化率大于閾值C1的條件后，計算壁面累積能量Ic.e，壁面累積能量大于C2可視為該區(qū)域存在空蝕風險。

3.4 空蝕風險預報結(jié)果

3.4.1 閾值C1、C2的選定

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

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

圖11 9組閾值計算的空蝕風險預報結(jié)果與試驗結(jié)果比較Fig.11 Prediction results of cavitation erosion risk calculated by 9 groups of threshold values compared with the experimental results

3.4.2 空蝕風險預報結(jié)果

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

圖12 空蝕風險預報結(jié)果與試驗結(jié)果的對比Fig.12 Comparison between cavitation erosion risk prediction results and painting test results

4 結(jié) 論

本文以NACA 0009 三維扭曲水翼為數(shù)值模擬及涂層試驗對象，將數(shù)值模擬空化形態(tài)與試驗高速攝影結(jié)果進行了對比?；趩蝹€蒸汽球泡的潰滅模型，分析了潰滅過程中能量轉(zhuǎn)化過程及相質(zhì)量的轉(zhuǎn)化過程，提出了空蝕風險的三個表征參數(shù)，并合理選定了表征參數(shù)的閾值?；谠摲椒ǎ疚膶@NACA 0009 三維扭曲水翼表面的空蝕風險進行了預報，并與涂層試驗結(jié)果進行了對比分析。主要結(jié)論如下：

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

（2）提出空蝕風險預報的三個表征參數(shù)：當?shù)貕簭?、當?shù)貕簭娮兓始氨诿胬鄯e能量Ic.e。在當?shù)貕簭姶笥陲柡驼羝麎杭爱數(shù)貕簭娮兓蚀笥陂撝礐1的前提下，通過壁面累積能量Ic.e大于閾值C2可以劃定空蝕風險的區(qū)域，通過Ic.e的大小可以判定區(qū)域內(nèi)空蝕風險的相對嚴重程度。其中，Ic.e的求解采用了基于傳質(zhì)源項的表達，避免了數(shù)值誤差。

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

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