王 巍 張慶典 唐滔 安昭陽 佟天浩 王曉放
(大連理工大學能源與動力學院,海洋能源利用與節能教育部重點實驗室,遼寧大連 116024)
空化存在于水力機械過流部件的低壓區域,是一種復雜的水動力學現象,涉及湍流、相變、可壓縮、非定常等復雜的流體力學問題[1-3].空化的發生常伴隨振動和噪聲[4],尤其是云空化發生時的非定常準周期性潰滅過程還會對船舶推進器[5]以及水利水電設備[6]、液壓機械[7]等過流部件造成持續性的沖擊和疲勞損傷.空化是造成流場流動不穩定的主要來源之一,為了確保水力機械高效、可靠和安全地運行,就必須要對流場流動的不穩定性進行主動控制和管理,至少避免和部分消除與系統過流部件發生共振現象的可能性.這就需要研究者們深入研究流場不穩定的原因,系統產生低頻壓力脈動的原因,進而分析空化發生的機制和演變過程,揭示空化的抑制策略和抑制機理,尋找更高效和更有針對性的抑制方法.
近年來,國內外學者對流場流動不穩定的原因進行深入研究并將其分為兩類:水力系統的不穩定性和固有不穩定性[8].系統的不穩定性受系統不同過流部件的結構設計所影響而較少的被研究.固有不穩定性由空泡本身產生,在壓力脈動頻譜上表現出比前者更高的脈動頻率,因而產生的危害更大,受到廣大學者的重視.例如研究較為成熟的回射流機制[9-10]和近年來通過試驗發現的激波機制[11-14]就屬于固有不穩定性.在不同的工況下,回射流機制和激波機制各占主導地位,造成片空化失穩和片空化向云空化的演變.
水力機械云空化的控制目的在于減小因空化非定常特性造成的壓力脈動,從而減小水力機械受到的應力.Tsujimoto 等[15-17]對水翼流動模型進行一維簡化和分析,發現空化流動中的低頻壓力脈動與空穴積對時間的二階導數成正比,揭示了空化流場中低頻壓力脈動產生的根源,將空化激振力與空化狀態的關系進行了定量描述.在空化狀態和發展的控制策略上,一般采取被動控制和主動控制的方法.被動控制是通過某種方式改變壁面特性來實現的,不需要向流場提供能量,因而容易實現,可操作性強.但很難實現對不同工況的交互精準調節[18-20].主動控制通過采取注入氣體、聚合物和水等的方法,實現對流場的控制.Timoshevskiy 等[21-22]和王巍等[23]分別采用在水翼吸力面布置切向射流水槽和不同角度的射流水孔進行射流的方法,都實現了對水翼云空化的抑制.
先前的研究者通過實驗和模擬等手段對云空化的發展過程進行了詳細的研究,揭示了空化的演變過程和回射流與空穴間的相互作用.但先前的研究更側重于對云空化準周期過程選取若干個特征時刻,對這些特征時刻的流場進行研究從而來概括流場的整個非定常演變過程.本文設想既然人們更關注于空化對水翼葉型吸力面或其他過流部件造成的影響,則吸力面附近(吸力線)是研究的重點.另一方面,空化流動的渦動力特性逐漸成為空化研究的熱點.Ji 等[24]研究者綜合Gopalan[25]的實驗結果,得出相比于氣體體積分數云圖,Q準則渦量圖更能準確地反應流場中復雜的流場結構.因此,本文采用對數值模擬的瞬態結果一維簡化的的方法,再增加時間維度重構二維的時空分布云圖,用于研究云空化的非定常周期性演變過程和射流抑制空化的抑制機制.對流場的非定常空化形態與渦結構的相互關系進行研究,旨在研究射流抑制空化的深層次抑制機理.
基于均相流模型的控制方程.假定汽液兩相為均質平衡流,兩相間無滑移速度,則汽液兩相的連續性方程和動量方程為[26]

混合相密度和黏度為

式中,下標i和j代表坐標方向,u和p分別為混合相速度和壓力,ρm,μm和μt分別是混合相密度和混合相層流/湍流黏性系數,α 代表體積分數,下標l 和v 分別表示液相和汽相.由于引入了新的未知量ρl和ρv,在求解空化流場時,除了需要對湍流進行封閉外,還需要求解汽液相密度,因此需要尋求混合相密度與其他量的關系,或引入空化模型進行求解.本文的研究過程采用引入質量輸運模型即空化模型進行方程封閉求解.
空化模型是描述汽液兩相相間轉化關系的數學模型.本文采用Schnerr-Sauer 空化模型對汽液相密度進行求解.Schnerr-Sauer 模型基于汽液均相流的假設,由Rayleigh-Plesset 方程推導而來[27].蒸汽體積分數具有如下的一般表達形式

這里,靜質量源項的表達式為

式中,Re為氣泡的蒸發速率,Rc為氣泡的冷凝速率.
與其他主流空化模型不同,Schnerr 和Sauer 將蒸汽體積分數和單位體積的液體中的氣泡數量進行聯系,有

式中,RB為氣泡直徑.因此,經過整理后的Schnerr-Sauer 空化模型為

本文采用密度修正的的RNGk-ε 湍流模型來對雷諾方程進行封閉.該模型在傳統的RNGk-ε 模型上考慮了汽液兩相密度變化對湍流黏度的影響[28].在標準RNGk-ε 模型中,湍流黏度為

修正的RNGk-ε 湍流模型考慮了汽液兩相密度變化對湍流黏度的影響,修正后的湍流模型采用以下兩式計算湍流黏度

式中,ε 為湍動能耗散率;cμ為模型常數,取0.085;通過密度修正指數n來對湍流黏度進行修正,本文中取n=3.
本文基于數值計算方法開展繞NACA66(mod)水翼的空化流動分析.計算的幾何模型,計算域的尺寸和工況條件的設置均與實驗條件保持一致.水洞試驗段的來流速度U∞為7.832 m/s,來流攻角為8?,流場溫度為290 K,飽和蒸汽壓1940 Pa,動力黏度μ=1.08 g/(m·s),空化數σ=0.99,雷諾數Re=5.1×105.計算區域劃分如圖1,水翼表面進行加密處理,最終網格數約為73 000.文獻[29]已對本文選用的網格和時間步長的選擇進行了無關性驗證,此處不再贅述.非定常計算的時間步長選擇為?t=5.0×10?4.
為了驗證所采用數值方法的可行性,將數值模擬結果與實驗結果進行對比.實驗設備及條件如文獻[30]所述.在實驗中采用高速全流場相機對瞬時空化流場形態進行捕捉.瞬時空化形態如圖2 所示,對比原始水翼空化流動的實驗測試結果,發現采用數值分析方法較好地預測出空泡準周期的生長、發展、脫落和潰滅的過程.周期為Tcycle=66 ms,頻率為f=15.15 Hz,則基于水翼弦長的斯特勞哈爾數為S t=fc/U=0.135 0,相比較該工況下實驗所得的斯特勞哈爾數S t=0.137 3,誤差僅為1.7%,說明所選用的數值分析模型和模擬結果能夠較好地預測云狀空化流場的周期性.但是,在對水翼尾緣脫落空穴的捕捉相比于實驗有一定的縮短或者厚度變薄,但周期性和變化趨勢基本一致,如圖3 所示.

圖1 計算區域與近壁網格劃分Fig.1 Computational domain and near foil meshing

圖2 一個周期內特征時刻空穴形態對比圖(σ=0.99,Re=5.1×105)Fig.2 Instantaneous outline of cavities of experiment and simulation in a time period

圖3 二個周期內空穴無量綱長度對比Fig.3 Instantaneous length of cavities of experiment and simulation in two time periods
本文首先對原始水翼模型,即表面無射流孔的NACA66(mod)水翼模型進行數值模擬,模擬工況為σ=0.99,研究其發生云空化時的非定常周期特性,其次對含有射流孔水翼的數值分析結果進行對比研究.原始水翼和射流水翼結構如圖4 所示.先前的研究表明對應安裝攻角下的水翼吸力面最高點附近的流場速度最高,且往往存在流動分離和損失.在最高點處布置射流水孔,有利于阻擋回射流,降低損失[29].因此,本研究中將射流水翼的射流位置布置于該安裝攻角下的吸力面最高點(即距水翼前緣0.19 倍水翼弦長位置).

圖4 水翼結構示意圖Fig.4 Schematic diagram of hydrofoil
圖5 給出了繞原始水翼空穴在特征時刻的瞬時形態變化.為了研究回射流區和空化區的相互作用,除了空穴的瞬時形態變化,圖上還分別顯示了壓力瞬時分布(等值線,單位Pa)和回射流區的瞬時形態變化.從圖中可以看出,模擬的結果較好地顯示了云空化發生時空穴的生長、回縮、發展和潰滅的準周期過程.在t1時刻,游離型空穴在水翼尾緣的下游處潰滅,其潰滅的瞬間出現局部的高壓,在從水翼尾緣到水翼前緣的壓力梯度的作用下著向水翼前緣傳播,并導致了水翼前緣附著型空穴的一次回縮.從t1到t2為附著型空穴的生長過程,t2時刻在空穴的尾緣第一次出現了回射流.回射流出現后空穴繼續發展,在t3時刻空穴充分發展,而回射流也在這個過程中推進到最靠近水翼前緣的位置,并幾乎遍布滿整個水翼的吸力面,回射流和空穴存在強烈的相互作用.t4時刻水翼表面的附著型空穴幾乎完全消失,出現了附著型空穴的二次回縮過程.脫落的空穴在回射流的作用下被抬起,脫離水翼表面.t5時刻幾乎完全消失的附著型空穴開始再次出現,而在水翼的尾緣附近,游離型空穴逐漸遠離水翼壁面.值得一提的是,游離型空穴在此過程體積有所增大,呈現膨脹的形態.t5和t6水翼前緣的附著型空穴繼續生長,而在尾緣附近的游離型空穴體積最大,相應于t6時刻,在水翼下游潰滅,局部高壓區再次出現,云空化的非定常準周期性發展進入下一周期.


圖5 繞NACA66 水翼非定常云狀空穴結構演化過程(空穴含汽率αv>0.15,周期Tcycle=66 ms)Fig.5 The numerically predicted vapor fraction(αv>0.15)of cavitation pattern around a NACA66 hydrofoil(Tcycle=66 ms)
為進一步深入地研究繞流水翼云空化發展過程和影響因素,本文在水翼吸力面近壁區設立一條監測線(monitoring line A-A),對水翼吸力面近壁區的流場特征參數進行監測,如圖6(a)所示.同時,為了研究邊界層空化流場特性,在時均空化流場的邊界層(UBL,mean<0.99U∞)的外部輪廓線(monitoring line B-B)上設立同步監測線,如圖6(b)所示.開展上述兩條監測線上空化流動特征參數的對比分析.由于云空化非定常準周期性,本文對一維的瞬時空化流場增加時間維度,獲得水翼吸力面近壁區和邊界層處的汽相體積分數、回射流速度、壓力系數等特征參數的時空分布云圖.進而分析繞流水翼云空化的發展和射流抑制空化機理,如圖6 所示.
定義當地壓力系數和無量綱壓力梯度

式中,p為當地壓力,而p∞為來流壓力.而壓力梯度值的正負代表著方向,正值代表著逆壓梯度,是造成回射流的主要原因之一;負值代表著從前緣到尾緣的順壓梯度.為了便于分析,壓力梯度云圖只顯示了壓力梯度的較高值(即|gradCp| >4).圖6 中的t1~t6時刻線均與圖5 相對應.
汽相體積分數的時空分布如圖6(c)和圖6(d)所示.圖6(c)清晰地反映了云空化非定常周期性過程,吸力面附近的空穴在整個周期內可以被分為3 個部分:當地含汽率大于0.7 的附著型空化區、小于0.7 的附著型空化Ⅱ區(汽液混合區)和位于翼型尾緣附近的游離型空穴區(空化Ⅲ區);圖6(d)的邊界層外部輪廓線監測結果顯示了空穴的發展厚度與邊界層的關系:在水翼前部邊界層監測線之外仍然存在汽相體積分數較大的Ⅰ區域,水翼中部的Ⅱ空穴區域基本被包含在邊界層內部,而脫落的游離型空穴在t4~t6時刻內已經被抬起到邊界層附近,不斷的向下游發展并最終潰滅.
回射流的時空分布云圖6(e)和圖6(f)所示,在t2時刻初現的回射流位于該瞬時附著型空穴的尾部(t[s]=t2,x/C=0.4).之后回射流的存在區域不斷擴大,并向水翼前緣推進造成空化Ⅰ區的回縮;隨著空穴的發展,回射流存在區域向水翼尾緣擴張并與空化Ⅱ區相摻混,造成當地含汽率的下降,說明回射流不斷被空穴汽化[31].圖6(f)可以發現在時均邊界層輪廓線上,幾乎沒有回射流的存在,主要原因是回射流起源于空穴尾部,并在水翼吸力面和空穴之間流動.另一方面說明了沿壁面法線方向,邊界層內部速度由負到正,存在很大的速度梯度.
壓力系數和壓力梯度的云圖如圖6(g)~圖6(i)所示,可以發現在游離型空穴潰滅的瞬時(t1時刻),產生了局部高壓[32],壓力系數可以達0.8 左右,且維持時間較短,產生了向水翼前緣的壓力梯度,造成了附著型空穴的一次回縮.此外,由圖6(g)和圖6(c)的對比可以發現,空化區位于流場的低壓區域中(Cp0.8).由圖6(i)和圖6(e)的對比可以發現,回射流區的發展受限于較高的壓力梯度(|gradCp|>4).

圖6 繞NACA66 水翼非定常云狀空化流場的時空分布云圖(空穴含汽率αv>0.15)Fig.6 Spatiotemporal distribution of unsteady cloud cavitation fl w around NACA66 hydrofoil

圖6 繞NACA66 水翼非定常云狀空化流場的時空分布云圖(空穴含汽率αv >0.15)(續)Fig.6 Spatiotemporal distribution of unsteady cloud cavitation fl w around NACA66 hydrofoil(continued)
圖6(l)顯示整個周期空穴無量綱面積變化,其中Scav為空穴面積,Sc表示水翼端部截面積,而Sre?cav為空穴與回射流的相互摻混部分的面積.可以發現空穴對回射流的汽化作用和回射流對空穴體的沖擊作用在回射流初現的時刻(t2時刻)就開始發生,直至周期結束(t6時刻).其中t5~t6時間內面積的變化顯示了空穴的充分發展,然后迅速潰滅的過程.
綜合汽相體積分數、回射流區域、壓力系數分布、壓力梯度分布的時空云圖可以得出以下結論:游離型空穴潰滅時產生的局部高壓導致附著型空穴的一次回縮,同時升力系數也有所下降(如圖6(k)所示).回射流出現于附著型空穴的尾部,并與空穴相互作用,其向水翼前緣推進,對附著型空穴產生沖擊作用;一部分隨空穴向下游發展被汽化,造成當地含汽率的下降.空化區存在于流場的低壓區域中(Cp0.8),回射流的區域也受限于較高的壓力梯度.高的壓力梯度一直存在,但回射流的首次出現需要時間積累.游離型空穴在回射流的作用下脫落和抬升、遠離固壁,導致升力系數的第二次急劇下降.在游離型空穴抬升和潰滅的同時(t3~t5),水翼尾緣處出現正負交接的壓力梯度,推測有渦的出現.
基于上述原始水翼的空化條件,開展了表面開孔射流這一主動控制策略對繞流水翼空化流場的影響研究,射流孔位于距水翼前緣0.19 倍弦長位置.圖7 顯示了射流水翼云空化發展的非定常周期性過程.t1時刻依然為游離型空穴潰滅產生局部高壓的瞬時.此后前緣的附著型空穴開始慢慢發展,t2時刻為附著型空穴發展到射流孔位置處,由于存在主動射流,此時空穴和射流開始相互接觸并摻混.而回射流在時刻在空穴尾端才首次出現,與原始水翼的回射流初現時刻相比,推遲了15.8%.回射流初現后不斷的發展(t3~t4時間內),之后回射流、射流和空穴三者間相互作用.在t4時刻回射流推進到最靠近水翼前緣點位置處,同時水翼尾緣處的空穴被托起,開始遠離壁面,脫落型空穴的面積不斷縮小,并最終在t6時刻潰滅.相比原始水翼的空穴發展周期66 ms,射流水翼明顯增加,達到82.5 ms,空穴脫落頻率降低了20%.


圖7 繞NACA66 射流水翼非定常云狀空穴結構演化過程(空穴含汽率αv>0.15,Fig.7 The numerically predicted vapor fraction(αv>0.15)of cavitation pattern around a NACA66 hydrofoil with jet fl w
射流水翼的時空分布云圖則更能清晰地反映射流、回射流和空穴三者間的相互作用過程,進而分析射流抑制空化的機理,如圖8 所示.圖中的t1~t6時刻線均與圖7 相對應.圖8(c)汽相體積分數的時空分布云圖反映了射流水翼云空化的周期性發展過程.與原始水翼不同,射流水翼前緣的附著型空穴發展較晚,在t1時刻產生的局部高壓也較小(Cp=0.2),并沒有監測到附著型空穴的回縮現象;在t2時刻附著型空穴發展到射流孔位置處;在t2~t4過程中,附著型空化區被射流切分為Ⅰ和Ⅱ兩個部分,射流孔及其附近含汽率較小.此外,圖7 的瞬時云圖顯示了該工況下依舊存在空穴的脫落,而在圖8(d)水翼尾緣邊界層外輪廓線附近并沒有監測到脫落的的游離型空穴.這說明射流水翼的空穴脫落和潰滅在剛離開固壁就已經發生,整個過程位于時均邊界層內部.相比于原始水翼,圖8(e)和圖8(f)顯示采用射流使得回射流強度減弱,回射流區域減小.這也可從圖8(i)和圖8(j)看出,在射流孔位置處出現較高的順壓梯度,該順壓梯度與造成回射流的逆壓梯度的分布相隔斷,對回射流強度和范圍起到了消減的作用.射流減弱了回射流對空穴的沖擊作用,空化發展受到抑制.
相比于圖6(l)無量綱空穴面積分析結果,采用主動射流后,繞流射流水翼的無量綱空穴面積有大幅度的減少,如圖8(l)所示,證明了該工況下射流抑制空化的有效性.

圖8 繞NACA66 射流水翼非定常云狀空化流場的時空分布云圖(空穴含汽率αv>0.15)Fig.8 Spatiotemporal distribution of unsteady cloud cavitation fl w around NACA66 hydrofoil with jet fl w

圖8 繞NACA66 射流水翼非定常云狀空化流場的時空分布云圖(空穴含汽率αv>0.15)(續)Fig.8 Spatiotemporal distribution of unsteady cloud cavitation fl w around NACA66 hydrofoil with jet fl w(continued)
Q是在伽利略變換下的速度梯度張量第二不變量,近年來,研究者認為以Q顯示的渦流能夠更清晰的反應流場內的渦旋結構,且空化與旋渦間存在交互作用[33-35].為了深入研究空化及射流抑制空化機理,本文利用Q判據對空化流場進行分析.對于三維流場,Q定義如下

對于二維的流場研究,上式可以簡化成

圖9 顯示了原始水翼的固壁監測線(monitoring line A-A)和射流水翼的固壁監測線(monitoring line的Q值時空分布云圖.Q的局部最大正值可以用來識別渦核,而負值則表示可能存在剪切但沒有旋渦運動的流動區域.圖9(a)原始水翼的Q值時空分布云圖顯示,在某一瞬時沿弦長方向壁面監測線上的Q值正負交錯,在水翼非定常發展過程中,Q值也發生變化.對比分析Q值和圖6 中的汽相體積分數、回射流的時空分布云圖,可以發現,Q值分布云圖中的正值與附著型空化I 區和游離型空化區(III)相契合,如圖中黃色顯示區域;負值與回射流的區域外輪廓線相契合,如圖中黑色虛線包圍的區域,而在正值與負值相互圍繞的中間區域則與汽液混合區(附著型空化II 區)相契合.這說明了在水翼前緣的含汽率較高的附著型空化I 區和水翼尾緣的游離型空穴內,存在著旋渦運動;回射流對周圍的空穴存在著剪切作用;而在水翼弦長靠中間位置,即存在著旋渦又存在著回射流的剪切作用,為含汽率相對較低的汽液混合物.圖(b)為在0.19 倍弦長位置上開設射流后的Q時空分布云圖,可以看出Q值的強度和范圍都得到了減小,且在射流孔的位置左右出現了正負值相互對立的隔斷面,這說明射流也存在著剪切作用,但該剪切作用能對前緣的附著型空穴向后發展和回射流向水翼前緣的推進形成阻擋,從而對云空化起到抑制作用.

圖9 水翼壁面監測線的Q 值時空分布云圖Fig.9 Spatiotemporal distribution contours of Q-criterion around NACA66 hydrofoil
本文采用數值分析方法,對NACA66(mod)原始水翼和射流水翼的云空化非定常過程開展研究,通過采用設置水翼固壁監測線的方法對流場一維簡化并加入時間維度,得到時空分布云圖,引入Q判據對時空流場中的渦旋特性進行研究,分析空化和空化抑制的機理.得到的主要結論如下:
(1)水翼尾緣游離型空穴潰滅時存在瞬間高壓,并向水翼前緣傳播;而采用主動射流后,該瞬間高壓得到顯著消減.
(2)回射流區域的發展受限于較高的壓力梯度.采用射流后,水翼吸力面壓力梯度減小,射流孔附近壓力增高,回射流初生的時間延遲,回射流的強度降低,有效抑制了空化的脫落發展.
(3)水翼空化的發展和流場中的渦結構有著緊密的聯系.Q正值分布與水翼前緣附著型空穴和尾緣的游離型空穴形狀相契合;Q負值分布與回射流的外部輪廓相契合.說明前緣和尾緣空穴內部存在旋渦作用,而回射流對空穴也存在剪切作用.
(4)射流對前緣附著型空穴和回射流具有剪切作用,抑制了回射流向前緣的推進和前緣附著型空穴的向后發展.