基于大渦模擬的波狀前緣水翼空化抑制研究

馬楷東，曹留帥，萬德成

（上海交通大學 船海計算水動力學研究中心 船舶海洋與建筑工程學院，上海 200240）

空化是一種復雜的相變現象，一直是水動力學領域的經典和基礎性問題。當液體介質中的局部壓力低于環境溫度下的飽和蒸汽壓時，液體介質中就會出現空泡。這種空泡的形成、發展和潰滅過程以及由此產生的一系列物理變化稱為空化?？栈哂忻黠@的空間流動性和劇烈的非定常性，往往會造成嚴重的水動力性能損失，并造成侵蝕、噪聲、振動等不利影響［1］。Wang 等［2］利用Clark-Y 翼型進行了空化試驗，通過改變空化數和水翼攻角，對不同條件下空化流動的結構進行了識別與分析，總結了水翼周圍空化流的定常和非定常特性。Ji 等［3］采用均勻空化模型和大渦模擬相結合的方法，對NACA66 水翼周圍的非定?？栈鬟M行了數值模擬，對空化周期性演化、空化—渦流相互作用和空化引起的壓力脈動等復雜流動行為進行了機理研究和總結。

由于空化現象的產生大多伴隨負面效應，因此對空化的抑制便成為了新的研究熱點。根據當前研究成果，可以將空化控制手段分為兩類：主動控制方法和被動控制方法。主動控制方法是指從系統外部向內部注入能量控制流動結構，進而控制空化，例如注入聚合物、水、空氣或利用超聲波等。Lu等［4］在NACA66水翼吸力面進行開孔射流，實現了對水翼片狀及云空化的高效抑制。王子豪［5］以NACA0015 水翼模型作為研究對象對局部自然空化流動的演變機理進行深入研究，通過空泡形態、壓力脈動、速度場等特性，重點分析沖擊波機制下空化流場的演化規律。主動控制方法雖然在空化抑制方面取得了一定的效果，但需要諸如射流、通氣等外部條件進行干預，在工程應用中實現難度大且成本較高。因此，作用力來自于系統本身且不借助于外部能量的被動控制方法被廣泛研究。車邦祥［6］通過微渦發生器方法對水翼空化進行了試驗和模擬兩方面的研究，對附著型空化的被動控制進行了詳實的研究，文章還總結了常見的被動控制方法，如加裝橫向障礙物、挖設溝槽或改變表面粗糙度等。其中，有一種運用仿生學原理，對水翼進行表面構型改變，以控制近壁面流體流動特性進而控制空化的方法引起了學者們的關注。研究發現，海洋生物座頭鯨的鰭肢前部有波狀的凸結，這種凸結結構可以幫助座頭鯨鰭肢在很大攻角下擺動，而不會產生失速現象，并且可以有效地控制流動分離，減小游動阻力。于是，開始有學者將這種波狀凸結構型應用到了控制空化的研究中。Custodio等［7］以NACA634-021翼型為基準，試驗測量了不同波長和波幅組合的波狀水翼在不同攻角下的初生空化數、升阻力系數和升阻比等數據，結果發現，波幅較大的波狀水翼空化區域范圍明顯減小，而波長似乎對空化現象的影響很小。陳柳等［8］采用SSTk-ω湍流模型對多種波狀水翼進行了數值模擬分析，發現凸結構型導致了一系列規律的對轉渦帶，限制了空泡的展向發展，對空化具有抑制效果。但由于雷諾時均方法（Reynolds-averaged Navier-Stokes，簡稱RANS）無法顯示空化流場的細節信息，大渦模擬方法（large eddy simulation，簡稱LES）開始得到應用。Pendar 等［9］對波狀水翼非定?？栈匦赃M行了大渦模擬研究，發現凸結的存在使來流匯聚到波谷處，從而加速波谷處的流體流動，使其壓力變低，空化易于從此處產生。文章重點研究了不同構型波狀水翼表面的流動分離現象對空化的影響，但并沒有計算諸如空泡體積等數據。Li等［10］基于LES方法對大攻角下的波狀水翼進行了數值模擬研究，結果表明減小波長和增加振幅不僅可以使空化體積減小約30%，還可以控制空化脫落的不穩定性，降低壓力幅值。但該規律特性對小攻角工況下的波狀水翼是否適用仍有待考量。因此，Zhao 等［11］運用等升力原理提出了不同的觀點，在詳細比較了小攻角波狀水翼空化的周期特性、壓力脈動以及空化對流向渦的影響等關鍵流動特性后發現，在該工況下波浪型前緣的空泡體積并沒有明顯減小，但是有效抑制了壓力脈動，文章分析可能是等升力原理導致了不同結論的產生。

通過以上分析發現，目前關于波狀水翼空化抑制的試驗和數值模擬研究還不夠充分，對其內在機理和影響規律的認識還不全面，因此，采用大渦模擬方法對三維波狀水翼的空化現象進行精細化仿真，通過空化周期、升阻力系數、壓力脈動以及流向渦結構等流場特性分析，探究前緣波浪曲線的波幅和波長對水翼空化的影響規律和內在機理，為波狀水翼的空化抑制研究提供參考和經驗。

1 數學模型

1.1 控制方程和LES方法

在模擬所采用的混合介質模型中，假設蒸汽和液體的速度和壓力相同，基本控制方程由連續性方程和動量守恒方程組成，表示形式如下：

式中：ui為i方向上的速度分量，p是混合物壓力。假設蒸汽和液體的壓力值相同，μ是層流黏度，ρ是混合介質密度。μ和ρ的表達式為：

大渦模擬方法通過濾波函數將大尺度的渦和小尺度的渦分離開，大尺度渦直接模擬，小尺度渦用亞格子模型來封閉，提供細節化的流場信息。將式（1）和式（2）進行濾波運算后，得到的LES方程為：

式中：下標l和v分別代表液相和氣相，αv是蒸汽體積分數。τij被稱為亞格子（SGS）應力，其定義如下：

為了構造亞格子應力的封閉模型，需要對未知的亞格子應力進行建模處理。假設亞格子應力與可求解尺度的應變張量系數-Sij成比例，將小渦對大渦的影響關聯起來：

式中：τkk為亞格子應力中各向同性部分；μt為亞格子模型的湍流黏度；δij為克羅內克符號，i=j時，δij=1，i≠j時，δij=0。采用的亞格子模型為Nicoud 和Ducros［12］提出的WALE（wall-adapting local eddy-viscosity）亞格子模型，該模型在求解空化問題上具有很好的適用性。其湍流黏度μt表達式為：

其中，Ls是亞格子尺度的混合長度，其定義如下：

其中，k為von Karman 常數，d為到最近壁面的距離，V為計算單元的體積，Cw為默認的WALE 常數，取值為0.5。

1.2 空化模型

采用Schnerr 和Sauer 基于Rayleigh-Plesset 方程所提出的Schnerr-Sauer 空化模型來描述氣、液傳質過程［13］。該模型的優點是源項表達式不包含經驗常數，具有較高的普適性?？栈Ｐ偷姆匠虨椋?/p>

其中，代表蒸發率m?+和冷凝率m?-的質量源項定義為：

2 數值模擬

2.1 幾何模型

將座頭鯨鰭肢前緣凸結引入到水翼設計中，采用波浪狀構型做仿生處理。由于座頭鯨鰭肢前緣比較肥大，因此選取NACA634-021翼型截面進行三維翼型建模，水翼弦長c=102 mm，厚翼特點明顯。模型截面的橫坐標以最大厚度0.35c處為分界點，前緣重構，尾緣不變，如圖1所示。波狀水翼前緣曲線由余弦曲線構造：

圖1 模型參數示意Fig.1 Model parameters and definition

式中：z為水翼的展向坐標，A為波幅，λ為波長。

為了研究基準水翼與波狀水翼空化的區別，并探究不同前緣波幅和波長對波狀水翼的影響和規律，進行了參數化建模，各水翼幾何模型如圖2所示，模型參數及命名見表1。

表1 前緣波幅與波長參數Tab.1 Leading edge amplitude and wavelength parameters

圖2 水翼幾何模型Fig.2 Geometric model of hydrofoils

2.2 計算域和邊界條件

計算域和邊界條件設置如圖3所示，水翼攻角αAoA= 6°，計算域入口位于翼前緣4倍弦長處，出口位于翼尾緣7倍弦長處，可以保證穩定的來流，并使尾流充分發展。

圖3 計算域和邊界條件Fig.3 Computational domain and boundary conditions

為了防止壁面的影響，兩側邊界做對稱平面處理；上下邊界設置為自由滑移壁面；入口和出口邊界分別設置為速度入口和壓力出口，出口壓力值p∞根據空化數計算公式得出：

式中：σ為空化數，pv為試驗環境下的飽和蒸汽壓（模擬選取的環境溫度為20 ℃），ρl為液體介質的密度，U為來流速度。具體的模擬工況在表2中給出。

表2 模擬工況參數Tab.2 Simulated operating parameters

在進行非定?？栈瘮抵登蠼庥嬎銜r，時間步長的選取尤為關鍵。一般通過庫朗數確定時間步長，這種方法結合時間參數與空間參數，通常選取庫朗數為1或小于1的值進行計算：

式中：Δt為時間步長，Δx為網格特征尺寸。模擬的時間步長設置為5×10-5s。

2.3 網格驗證與對比

2.3.1 網格無關性驗證

網格采用結構化網格劃分形式。為了確保模擬的準確性和精確度，設置了3 套網格進行網格無關性驗證。3套網格按網格量分為粗網格、中等網格、細網格，對無空化下的升阻力系數進行了對比，公式為：

式中：FL和FD分別為升力和阻力，S為水翼的正投影面積。網格無關性對比結果在表3中給出。由于大渦模擬對網格質量要求嚴格，3 套網格通過調整第一層網格高度保證水翼表面處的無量綱距離y+值均小于1，如圖4所示。

表3 網格無關性驗證Tab.3 Grid independence verification

圖4 網格y+值分布Fig.4 Distributions of y+ under different grids

在表3中，中等網格和細網格的模擬結果比較接近。在保證模擬精度和平衡計算量的前提下，選取中等網格進行后續的數值模擬計算分析。圖5給出了基準水翼和不同的波狀水翼表面網格和側面邊界網格劃分情況。

圖5 不同水翼網格劃分Fig.5 The meshes around different hydrofoils

2.3.2 試驗結果對比

為了驗證所建立數值模型的適用性和正確性，選取基準水翼和2M 水翼進行升阻力系數的試驗驗證。參考Johari［14］得出的雷諾數Re=4.5×105的試驗數據，對比結果在表4 中給出。結果表明一致性較好，說明所建立的數值模型具有良好的適用性。

表4 計算與試驗對比Tab.4 Simulation and experimental comparison

3 結果分析

3.1 空化周期

為了對比研究基準水翼和波狀水翼空化現象的區別，對不同水翼的典型空化周期進行可視化處理，如圖6所示，將空化蒸汽體積分數為0.49的等值面和水翼表面壓力系數Cp疊加顯示，以便直觀地揭示波狀水翼前緣凸結對空化產生和發展的直接影響。

圖6 典型空化周期（等值面αv=0.49）Fig.6 Typical cavitation cycle (The iso-surface of αv=0.49)

對于基準水翼，在圖6（a）的1/8T～4/8T過程中，空泡從前緣初生，附著在水翼表面并在整個展向上連續分布，不斷向后發展，當逆壓力梯度足以克服水翼表面附近的湍流強度時，造成空泡斷裂的回射流產生。如圖6（a）的5/8T所示，空化由片狀空化發展成云狀空化，并出現空泡破裂和脫落現象。前端斷裂的片空泡不斷向前收縮，直至全部消失，圖6（a）的6/8T～7/8T可以清晰地反映這一現象。而如圖6（a）的T所示，尾端脫落的云空泡隨流場向后運動，當局部壓力恢復到飽和蒸汽壓以上時，云空泡消失，并在水翼前緣開始新一輪的空化現象。觀察空化周期云圖發現，波狀水翼的空化過程和基準水翼具有相似性，但由于前緣波浪構型使得空泡只在波谷處產生，如圖6（b）、（c）、（d）所示，初生的片空泡在整個展向上是不連續的。對2M 水翼進行分析，發現波峰處的壓力相對較高不利于空泡的產生，空泡起始于波谷處的低壓區并向后延伸，而波峰似乎起到了分隔作用，限制了空泡的展向發展。當波谷處產生的片空泡延伸到一定長度后，受到來自前緣波峰的影響越來越小，空泡在展向上開始融合，并發展成與基準水翼相同的脫落模式，往復循環，周而復始。結合相關文獻，將這種波峰波谷對空化發展的影響定義為分區效應。其內在機理是通過前緣波浪構型將來流進行整合，使流經波谷處的流速加快，壓力降低，利于空化的產生，而波峰處由于流體介質在波谷處的聚集，相對壓力較高，則沒有空化的產生。

為了更好地探究前緣參數對空化的控制效果和作用規律，通過改變波狀水翼前緣構型的波幅和波長進行對比分析。如圖6（c）的2L 水翼所示，增大了前緣波幅，其整體的空化流動模式與2M 水翼相同，但仔細對比后發現，經過發展后的片空泡并沒有在展向進行很好地融合，其每個波谷產生的空泡在展向上幾乎相互獨立，在回射流的作用下片空泡斷裂，脫落形成的云空泡也是各自獨立的，只在遠場會產生干擾，波幅較大的波狀水翼其空化分區效應更加顯著。觀察圖6（c）的T發現，新一輪的空化周期開始時，兩側波谷的空泡要比中間的空泡發展更快，說明各個波谷產生的空化自成體系，循環周期不同步。對于波長的影響，在圖6（d）所示的4M 水翼空化過程中，不同波長下的波狀水翼空化發展模式相同，但是在展向空泡融合上，小波長水翼要弱于大波長水翼，波谷處產生的空泡雖然獨立性明顯，卻有著同步的循環周期。因此，波幅和波長的改變都會對空化流場產生影響，增大波幅和減小波長都增強波狀水翼的分區效應，控制空泡的產生和發展。

為了進一步研究前緣波浪構型對空化控制的影響，圖7對水翼的空泡體積進行了時均處理。對比發現，波狀水翼的時均空泡體積都低于基準水翼，說明波浪型前緣的引入可以對空化進行抑制。進行定量分析后發現，2M 水翼的空化抑制率為15.7%，2L 水翼的空化抑制率比2M 水翼稍高為18.6%，說明波幅較大的波狀水翼對空化的控制效果更明顯。而4M 水翼的空化抑制率為27.9%，相比于2M 水翼對空化的控制水平大大提高，說明減小波長可以大幅提高波狀水翼對空化的抑制效果。這種大波幅和小波長控制效果良好的結果應該與分區效應導致空泡間相互獨立，限制了展向抬升有關。

圖7 時均空泡體積Fig.7 Time-averaged vapor volumes

3.2 升阻力性能

圖8顯示了不同水翼的時均升力系數和阻力系數。在空化過程中，由于波狀水翼的空化受到抑制效果，其升力系數相比于基準水翼都得到了不同程度的提升。對比2M和2L水翼，升力系數隨波幅的增大而減小，與波幅大小成反比，阻力系數則與波幅大小成正比關系，升阻比隨波幅的增大而減??；對比2M 水翼和4M 水翼，升力系數與波長大小成正比關系，前緣波長越大，升力系數越大，阻力系數與波長成反比關系，升阻比隨波長的減小而增大。結果表明，波浪型前緣可以一定程度上提高水翼的升力系數，但阻力系數也會隨之變大，無法保證升阻比的變化規律，故合適的幅長比有利于提高水翼的水動力性能，超過一定限值后性能則會下降。

圖8 時均升阻力系數Fig.8 Time-averaged lift coefficients and drag coefficients

下面進一步對升力系數進行時域統計。從圖6中可以看出水翼空化是具有周期性的，故將物理時間t無量綱化，處理后的特征周期參數T定義如下［11］：

如圖9 所示，截取了4 種水翼空化過程中一個典型周期內的升力系數和空泡體積分布?？梢郧逦匕l現基準水翼升力系數的數值振蕩比較劇烈，非定常特性明顯，存在較大的升力波動。升力的產生源于水翼上下表面的壓力差，故也可從側面分析出基準水翼的壓力脈動是比較明顯的。而波狀水翼升力系數的振蕩處于較低水平，沒有明顯的峰值，表明在空化狀態下，波狀水翼獲得的升力相對穩定，受到空泡發展和潰滅的影響較小，水動力性能優越。

圖9 一個典型周期內空泡體積和升力系數的時域曲線Fig.9 Time-domain curve of the vapor volume and lift coefficient in one typical cycle

3.3 壓力脈動

將表面壓力系數做時均化處理，可以更好地分析波狀水翼對空化控制的內在機理，幫助理解空化的非定常特性及規律。如圖10所示，深色附近為低壓區，即空化初生的區域。在基準水翼表面低壓區連續均勻，所以空泡也是在整個水翼展向均勻分布的。而波狀水翼的低壓區只在波谷處存在且不連續，波峰處顯示為高壓區，故空化不從此處產生。

圖10 吸力側時均壓力分布Fig.10 Distribution of time-averaged pressure on suction side

對比圖10（b）和圖10（c），不同波幅下波狀水翼波谷處低壓區在展向和流向并沒有太大區別，但是由于波峰和波谷的間距不同，波幅越長的水翼，前端無空化區域（淺色區域）面積就越大，空泡的分區效應越明顯，展向發展受到的限制越強，在尾部也越容易脫落，空化現象減弱。對比圖10（b）和圖10（d），由于波長的改變，波谷區域的低壓區也沿展向隨之縮小或增大。波狀水翼的波長越小，波谷低壓區面積越小，產生的空泡在展向上越難以發展和融合成大型附著型空泡。當受到由于逆壓梯度產生的回射流作用時，空泡很容易發生斷裂和脫落，故空泡長度也會隨著前緣波長的縮小而變短，總體空化現象減弱。

為了更加直觀地反映波狀水翼對壓力脈動的控制效果，在基準水翼和波狀水翼中縱截面的吸力側從空泡初生位置到空泡潰滅位置等比例依次布置了5 個壓力測點，對表面壓力進行監測。測點的具體位置在圖11中給出。

圖11 壓力測點位置分布Fig.11 Position distribution of pressure measuring points

對監測得到的壓力進行快速傅里葉變換（FFT）后得到頻域圖，如圖12 所示，每個水翼的測點位置都在相應的頻域圖中給出。在圖12（a）中，壓力譜曲線上出現一個峰值，說明空化發展具有較強的周期性。從P1測點到P5 測點壓力系數幅值逐漸升高，對應位置會發生片空泡的斷裂和云空泡的潰滅兩個重要現象，說明這兩種現象是大幅度壓力脈動的主要來源。

圖12 吸力側壓力頻域圖Fig.12 Frequency-domain diagram of suction side pressure

同樣，波狀水翼也具有周期性，對比圖12（a）和圖12（b）～（d），結果表明波狀水翼的壓力脈動得到了明顯抑制，壓力脈動幅值顯著降低。但規律性與基準水翼類似，都受到了片空泡斷裂和云空化潰滅的影響，從而產生大幅度的壓力脈動。值得注意的是，4M 水翼靠近尾緣的P5 測點存在一個突出的幅值，這是由于4M 水翼受到前緣構型的影響，空化泡長度較短，只有少部分的空泡能夠到達P5 測點位置，故壓力較高幅值較大。在不考慮4M 水翼P5測點的前提下，相比于基準水翼，2M 水翼、2L水翼和4M 水翼的壓力脈動幅值分別減小了55.3%、67.3%和74.6%，說明前緣波浪構型可以顯著降低水翼表面的壓力脈動。繼而分析不同前緣參數的波狀水翼對壓力脈動的影響，發現小波幅水翼對壓力脈動的控制不如大波幅水翼，2L水翼比2M水翼的壓力脈動幅值減小了43.3%；而前緣波長對壓力脈動的影響相對較弱，4M 水翼相比2M 水翼的壓力脈動幅值只減小了26.9%。綜上所述，前緣波浪構型對水翼表面的壓力脈動都有一定的控制效果，這種控制效果的規律與對空泡體積的控制相似，隨著波幅的增大和波長的減小，控制效果越增強。

3.4 時均流向渦

大量研究表明，非定?？栈Y構與渦旋運動密切相關?？栈莸漠a生和發展會產生復雜的渦結構，進而改變流場中的渦量分布。因此，有必要探究空化過程中由于波浪型前緣誘導產生的流向渦的發展規律和流場特征。如圖13 所示，對時均流向渦進行0.2c、0.6c和1.2c這3 個截面的可視化分析，并與時均蒸汽體積分數為0.49的等值面疊加顯示。結果顯示，基準水翼的渦結構碎片化嚴重，沒有明顯分布規律，而波狀水翼表面則存在明顯的對轉渦結構。

圖13 時均蒸汽體積分數為0.49的等值面和時均流向渦Fig.13 Time-averaged volume fraction iso-surface of 0.49 and time-averaged streamwise vortices

流向渦由于前緣波浪構型的影響，導致相鄰的渦旋相向移動，形成一個收斂模式，或者相反運動，形成一個發散模式，分區效應明顯。相鄰渦的合并導致空泡從水翼表面抬升并限制了渦量在展向的拉伸。在0.2c平面處，可以清楚地觀察到這種對轉渦結構，且沿流向向上抬升，這種抬升應該與波谷產生的空泡有關，形成的空腔越厚，抬升也越高?？紤]不同前緣參數的影響，發現波幅越大或者波長越小，對轉渦結構越明顯，抬升也越高。在0.6c平面處，由于云空泡潰滅的影響，流向渦的分布開始變得無序，這在圖13（a）和13（d）中尤為明顯，而圖13（b）和13（c）的時均空泡長度較長，還顯示為規則的對轉渦結構。在1.2c平面處，波狀水翼的流向渦也表現出與基準水翼相同的不規則性，這主要是由于水翼尾流的流向渦結構受云空泡潰滅的影響較大，幾乎不再受波浪前緣分區效應的影響。

而分析時均空化等值面發現，波狀前緣構型的分區效應明顯，每個波谷處產生的時均空泡獨立存在，沿流向貫通前后，互不干擾。對比分析2M、2L 和2M、4M 兩組水翼后發現，大波幅構型限制了空泡在波谷處的展向發展，小波長構型則限制了沿流向的空泡長度，也從另一方面證實了3.3節中對時均表面壓力分析得到的結論。

4 結 語

基于數值模擬方法對具有波浪型前緣的水翼進行了非定?？栈鲃犹匦苑治?，采用LES 方法對流場進行精細化模擬，對比分析了不同波幅和波長參數下波狀水翼的瞬時和時均流場特性，主要結論如下：

1）相比于基準水翼，波浪型前緣的引入可以對空化進行有效地抑制，空化從波谷處產生，分區效應明顯。隨著波幅的增大和波長的減小，波狀水翼對空化的控制效果變得明顯，流場中的空泡體積逐漸降低。

2）前緣波浪構型可以提升水翼的升力系數，并顯著降低水翼表面的壓力脈動幅值，波幅越大，波長越小，壓力脈動幅值越低，流動越趨于平穩。合適的波幅和波長可以使升阻力性能得到最優解。

3）不同前緣參數下波狀水翼渦結構的演化是相似的，即由于前緣的分區效應導致沿流向呈現出對轉渦結構。進一步分析表明，空泡發展與潰滅的整個過程對渦結構的發展具有顯著影響，進而改變水翼的水動力性能。