李 勇,劉 成,萬德成
(上海交通大學 船海計算水動力學研究中心 (CMHL) 船舶海洋與建筑工程學院 海洋工程國家重點實驗室, 上海 200240)
空化流動是一種復雜的氣液兩相流動,它涉及到流體力學中諸多的流動問題,包括相變、質量傳輸、可壓縮性和非定常性。當液體中的局部壓力低于當前溫度下的飽和蒸汽壓力時,空化現象便發生了,所以空化常常發生在速度和壓力梯度大的水利設備中,例如泵和螺旋槳[1]??栈鲃邮欠嵌ǔ5模驗榭栈F象伴隨著空泡周期性的產生與潰滅,這會給水利設備造成不利的壓力震蕩與疲勞破壞。不管是試驗手段還是理論手段,都很難直接研究空化現象的機理機制,現在對空化現象的研究都是從結構簡單的水翼著手,逐步往結構復雜的螺旋槳或者泵發展。
根據空化不同的形態和形式,空化可以被分為片空化,云空化,超空化和渦空化等。片空化本身就是非定常的,在生成并發展到一定的程度后就會坍塌破碎成一束束細碎的空泡,并從附著在水利設備上的片空化脫離開來,形成云空化[2]。大量的觀察和研究被用來分析片空化發展成云空化的機理機制,最先被發現的片空化脫落機制是回射流機制。Knapp等[3]發現在空化流動中存在與主流方向相反的回射流。Furness等[4]對一個二維的縮放噴管進行理論和試驗研究,發現在片空化發展的過程中,一股與來流方向相反的回射流,從片空化末端的底部,緊靠著管壁向上游移動,并在噴管喉部將片空化切斷。Kawanami等[5]通過高速攝像技術和壓力檢測技術對云空化進行一系列的試驗觀察,發現片空化的脫落和云空化的產生是由于片空化的底部存在著由尾部向導邊移動的回射流。同年,Kawanami等為了驗證回射流在空化流動中的作用,在水翼上表面的不同位置布置了長短不同的擋板,來觀察擋板的存在對回射流以及片空化脫落效果的影響,研究發現,片空化底下的擋板能阻止回射流向上游移動,從而一定程度上緩解片空化的脫落。以上研究都定性地說明回射流的存在會使片空化脫落并產生云空化。Callenaere等[6]通過觀察回射流相對于片空化厚度的大小,定量地分析了回射流的相對厚度對片空化穩定性的影響。
近年來,片空化斷裂和脫落的另外一種機制被發現,也就是大尺度空化瞬間潰滅的激波機制。Arndt等[7]通過綜合試驗和數值模擬的手段,對一個二維NACA0015水翼的空化發展過程進行研究,研究表明,根據σ/(2α)取值的大小(式中的σ為空化數,α為水翼的攻角),可以把片空化過渡到云空化的機制分為兩類,當σ/(2α)的值比較大時,回射流機制是片空化發展成云空化的主要原因;當σ/(2α)的值比較小時,激波機制是片空化脫落生成云空化的主要原因,同時Arndt等[7]指出σ/(2α)=4是取值大小的過渡階段,σ/(2α)≈4時兩種控制機制都會對片空化的脫落產生影響。也就是說,水翼在相同的攻角下,空化數σ越小,激波機制越明顯;空化數σ越大,回射流機制則成為空化脫落的主要原因。Ganesh等[8]通過高速攝像技術和X射線來觀察片空化向云空化的過渡過程,通過系統地改變空化數σ的大小,在楔形頂點的流動分離區觀察空化的水動力性能,發現在空化數較小時,激波機制成為了片空化脫落的主要因素。
對空化的研究一般采用試驗與數值模擬相結合的方式。在試驗中激光多譜樂測速(LDV)技術和粒子圖像測速(PIV)技術是觀察空化最常用的方法。Kubota等[9]通過高速攝像機和LDV技術觀察到云空化一般處在最大渦強中間,并且發現云空化是由很多細小的空泡構成的,而片空泡則要純凈得多,幾乎是塊狀的單質空泡。Foeth等[10-11]則通過時間分辨粒子圖像測速(TR-PIV)技術觀察回射流導致片空化的脫落和云空化的產生?,F在采用數值模擬對空化現象進行研究的方法也流行起來,新型湍流模型和空化模型的更新和選取,讓數值模擬的結果與試驗結果更好的吻合。Hsiao等[12]提出了一種基于歐拉-拉格朗日耦合方法的多尺度空化模型,通過小尺度模型捕捉微觀空泡結構,大尺度模型追蹤大空泡團的動力學特性,從而清晰地描述小空化核生長成為宏觀空泡的過程。
針對回射流對片空化脫落的作用,研究阻斷回射流對片空化發展的影響。為此,分別對平板和水翼的吸力面添加擋板,用來阻擋逆流而上的回射流,從而研究不同布置不同長度的擋板對回射流的影響,以及比較分析添加擋板后的水動力性能。數值模擬采用開源CFD平臺OpenFOAM的interPhaseChangeFoam求解器,選取Schnerr-Sauer空泡模型和修正后的SST k-ω湍流模型對帶擋板的水翼進行數值求解,并將模擬結果與試驗結果相互比較。
在給出控制方程之前,先定義一個體積分數α,在網格中代表液體的相分數。在控制方程中,把氣液兩相看作是單質混合模型,共享速度和壓力,而速度場需要相分數α來確定??紤]到氣液兩相系統需要在同一個網格里面表示,網格內全是液體時,則α=1;網格內全是氣體時α=0;當相分數α介于0和1之間時,網格內是氣液混合物,同時也看作是水氣界面。由于在控制方程中只能用到一個體積分數,因此把液相的體積分數αl寫作α,而氣相的體積分數αv由αv=1-α得出。求解空化現象的控制方程由三部分組成:液體的相方程,連續性方程和動量方程。

(1)
(2)

(3)

U=αlUl+αvUv
(4)
p_rgh=p-ρgh
(5)
采用的空化模型是基于輸運方程模型(TEM)的Schnerr-Sauer空化模型[13],該模型將水氣混合物看做是包含大量球形蒸汽泡的混合物。根據Rayleigh-Plesset方程得到的不同源項用來模擬空化的汽化與液化過程,當局部壓力低于飽和蒸汽壓力pv時,發生劇烈汽化;當局部壓力高于飽和蒸汽壓力pv時,空泡液化為液態??栈P涂梢杂墒?6)表示:

(6)
其中,Cv和Cc是經驗常數,一般取1;ρl,ρv和ρ分別是液體、氣體和氣液混合物的密度;p是局部的壓強,pv是液體在當前溫度下的飽和蒸汽壓力;Rb是氣核半徑,可由式(7)求出,式中的n是液相中氣泡的數目。
(7)
Larsen和Fuhrman[14]指出,常用的標準二方程閉合模型(k-ω和k-ε)都是無條件不穩定的,基于標準k-ω湍流模型變形來的SST k-ω湍流模型也有著同樣的問題,即無條件不穩定,且有著漸近的不穩定增長率,使得部分區域的黏度總是大得不符合物理現象。為了避免過度預測湍流閉合區域的湍動能和渦黏度,采用一種新的穩定化的SST k-ω湍流模型來閉合控制方程。在新的SST k-ω湍流模型中,Larsen和Fuhrman推薦采用一種簡單的辦法來將湍流模型穩定化處理:通過對湍流黏性系數進行簡單的修正,從而達到穩定湍流模型的目的。具體的思路就是將渦黏度νT由式(8)重新定義為式(9)。

(8)

(9)
在修正后的渦黏度公式中,分母中的max函數,除了原來就存在的兩個判斷大小的依據,新添加了第三個判斷依據,在接近勢流的流動條件下才會被激活。式中的a1,λ2和β*都是經驗常數,其值分別取0.31,0.44和0.09,p0為湍流產生量,k為湍流動能密度,ω為湍流耗散率,p0為湍流產生量,pΩ為修正后的湍流產生量,F2是與k和ω有關的正切函數。根據解析分析,并通過這樣的修正,湍動能和渦黏度不再隨著計算的進行增長,而是漸進指數衰減[13]。
為了研究擋板對回射流的影響以及回射流與片空化脫落之間的關系,采用了兩個物理模型來進行分析。第一個是根據Zhang等[15]的試驗設計一個平板水翼,其形狀如圖1所示,是一塊尺寸為150 mm×200 mm×12 mm(長×寬×高)的平板水翼,其中150 mm是平板水翼的弦長,攻角固定為0°,一條全展長的擋板設置在平板水翼的上表面,擋板高2 mm,寬2 mm,距離導邊0.37倍弦長。
在這個工況下,來流速度為10 m/s,飽和蒸汽壓力取3 540 Pa。出口壓力根據空化數求得,空化數的定義由式(10)給出,其中P取出口壓強,Pv是飽和蒸汽壓。采取的空化數和文獻[14]中采用的一致,取值0.7。在整個計算域中,進口距離平板水翼的導邊為2倍弦長,出口距離平板水翼的隨邊為4倍弦長,計算域的高度為2倍弦長,展長選取0.3倍弦長。具體的邊界條件由圖1給出。

圖1 平板水翼的尺寸及其計算域Fig. 1 Sizes of the flat hydrofoil and its computation domain

(10)
由OpenFOAM自帶的blockMesh工具來畫整個計算域的背景網格,然后用snappyHexMesh工具,對平板水翼附近的一個矩形的加密區進行網格細化。分別采用10×105,15×105,18×105和20×105的網格進行計算,18×105網格的計算結果和20×105網格的計算結果已經比較相近,為了保證水翼附近的網格足夠細,最終選取了20×105的網格來進行數值計算,平板水翼中剖面處的網格分布如圖2所示。

圖2 中剖面處的網格分布Fig. 2 Mesh distribution in mid-plane
第二個物理模型采用的是Kawanami[5]所做的一系列試驗中用到的對稱型橢圓鼻水翼(elliptic nose foil),橢圓鼻水翼的形狀和理論坐標由圖3給出。橢圓鼻水翼弦長150 mm,并在一個高和寬為600 mm×150 mm的螺旋槳空泡水筒中進行試驗。

圖3 橢圓鼻水翼的理論坐標系統Fig. 3 E.N. foil and coordinate system
在第二種工況下,計算域的尺寸與進行試驗的空泡水筒尺寸一致。計算域的入口距離橢圓鼻水翼的導邊為2倍弦長,計算域的出口距離橢圓鼻水翼的導邊6倍弦長,計算域的寬度和高度與空泡水筒的寬度和高度一致,也就是1倍弦長和4倍弦長。計算域的布局和邊界條件由圖4給出。

圖4 橢圓鼻水翼的計算域和邊界條件Fig. 4 Computational domain and boundary conditions of the E.N. hydrofoil
為了研究不同長度不同位置的擋板對回射流的影響,同時考慮到試驗中的最大空化長度為0.5倍弦長,分別選擇在橢圓鼻水翼吸力面距離導邊0.37倍弦長處和0.60倍弦長處設置全展長的擋板。同時為了研究不同長度擋板的抑制效果,一個只有十分之一展長的擋板同樣布置在0.37倍弦長處。擋板布置示意如圖5所示。橢圓鼻水翼的計算工況由表1給出,表中的k和ω分別為湍動能和比耗散率。不同于平板水翼的第一個工況,橢圓鼻水翼的計算工況為20°C下的飽和蒸汽壓力2 338.8 Pa,水翼的攻角為6°,來流速度為5 m/s。

圖5 擋板長度與位置分布Fig. 5 Arrangement of the obstacles used in simulations
網格的生成方法與工況一相同,都是通過OpenFOAM自帶的blockMesh工具生成背景網格,然后用snappyHexMesh生成加密區。計算采用的網格有33×105,35×105,37×105和40×105這4套,為了使得計算準確同時提高計算效率,最終選取了一套37×105的網格,添加擋板與不添加擋板的網格量都近似在37×105。圖6給出了橢圓鼻水翼在中剖面處的網格,其中,圖的左右下角,分別是橢圓鼻水翼導邊和隨邊的網格分布。圖7為在0.37倍弦長處,十分之一倍展長的擋板局部網格分布。

圖6 中剖面以及水翼兩邊的網格分布Fig. 6 Mesh distribution in mid-plane and near the leading edge and the trailing edge

圖7 布置(d)中的擋板附近網格分布Fig. 7 Mesh distribution around the obstacle in configuration (d)
首先分析平板水翼的物理模型。圖8給出了平板水翼在空化數為0.7時的升力系數歷時曲線。從圖中可以看出,擋板的存在會改變升力系數曲線的相位,尤其是在空化發展到0.06 s之后,添加擋板的水翼與不添加擋板的水翼呈現出完全相反的相位。

圖8 空化數為0.7時的升力系數歷時曲線Fig. 8 Time history of lift coefficient at cavitation number 0.7
值得注意的是,在給平板水翼0.37倍弦長處添加擋板的前后,數值模擬計算得到的平均阻力系數幾乎沒變,不加擋板前的平均阻力系數為0.076 6,添加擋板后的阻力系數為0.079 67,變化不大。而計算得到的升力系數卻發生了很大的變化,盡管添加擋板前后,升力系數的值都比較小,但對平板水翼添加擋板使得水翼的平均升力系數由0.003 565上升到了0.010 733,升力系數為無擋板時的3倍??梢缘贸?,在給水翼添加適當的擋板后,升力系數的增幅會遠大于阻力系數的增幅,也就是說,擋板的存在,可以以相對較小的阻力代價,換取巨大的升力提高,擋板的合理使用,在工程應用中有著巨大的潛能。在后文對平板水翼表面壓力分布和空泡形態的分析中,解釋了擋板的存在會使得升力系數得到大幅度提高的原因。
圖9和圖10分別為平板水翼上表面的平均壓力分布和瞬態壓力分布。從圖中可以看出,擋板的存在改變了平板水翼上表面的壓力分布,無論是從瞬態上觀察還是直接分析平均壓力分布,都能發現擋板的存在降低了擋板后面區域的壓力。

圖9 平板水翼上表面的平均壓力分布Fig. 9 Average pressure distribution on the upper surface
圖9的平均壓力分布更加能說明問題,在圖9中橢圓標記出的區域,擋板的存在會在擋板后的一片區域制造低壓區,在橢圓外,添加擋板與不添加擋板的壓力區別不大。添加擋板會在平板水翼上表面造成低壓區,使得上下表面的壓差增大,這很好地解釋了為什么擋板的存在會使平板水翼的升力系數大幅度的提高。后文通過對上表面空化形態的觀察會發現擋板的存在會抑制片空化向后發展,也因為空化在擋板后面難以發展,壓力得以降低,低壓區得以形成。圖11為平板水翼在0.25 s時的瞬態壓力分布云圖。

圖10 平板水翼上表面的瞬態壓力分布Fig. 10 Transient pressure distribution on the upper surface

圖11 0.25 s瞬時壓力分布Fig. 11 Contour of pressure at 0.25 s
圖12展示了平板水翼空化界面與回射流的相互作用。從圖中可以看出,在片空化發展到一定程度時,片空化底部的回射流開始逆著來流方向向上游移動,擋板的存在會阻礙回射流繼續向上發展,從而有效抑制片空化的脫落。雖然回射流的速度會比主流速度要小,但回射流速度和主流的速度是同一量級的。

圖12 平板水翼周圍的空化形態和速度矢量Fig. 12 Cavitation shape and velocity vectors around the hydrofoil
橢圓鼻水翼的抑制效應則更加明顯。如前文中提到的空化脫落的兩種機制,根據σ/(2α)取值的大小,當σ/(2α)的值比較大時,回射流機制是片空化發展成云空化的主要原因;當σ/(2α)的值比較小時,激波機制是片空化脫落生成云空化的主要原因,兩種機制的分水嶺在于σ/(2α)=4,基于文中的數值模擬條件計算得出該工況下σ/(2α)=5.11,略高于這個分水嶺,可認為回射流為該工況下空化脫落的主要機制。圖13給出了橢圓鼻水翼在0.37倍弦長處的擋板對回射流的阻擋現象。在片空化底部生成的回射流受擋板的阻礙后改變了流向,不再逆流而上。與擋板相撞后的回射流不再有繼續向上的沖勁,會順著主流,隨脫落的空泡一起翻卷流向下游,并在擋板的后方形成一個渦。此后,片空化不再發展到擋板之后,擋板的存在抑制了空化區的形成以及空化的脫落。

圖13 布置(b)工況下擋板對回射流的阻擋過程Fig. 13 Block phenomenon of placed obstacles on the re-entrant jet in configuration (b)
當沒有產生空化現象時,在一個光滑的水翼表面安裝一小塊擋板無疑會增加水翼的阻力,從而導致水翼的升阻比減小,效率降低。而當空化現象發生后,情況就出現了反轉。在水翼表面布置合適的擋板可以以較小的阻力代價來獲得巨額的升力提升,從而提高升阻比,提高水翼的效率。這是因為在水翼吸力面低壓區的局部壓力下降到飽和蒸汽壓后,壓力就因為空化的產生而無法繼續下降。而添加合適的擋板可以改變這種情況,擋板的存在可以改變吸力面的壓力分布,減小空化區,并在擋板后的一小塊區域形成新的低壓區,從而有效地提高升力,提高升阻比,提高水翼的效率。圖14是圖5給出的不同擋板布置下,橢圓鼻水翼在上表面的瞬態壓力分布圖。
從圖14可以看出,不同布置的擋板會不同程度地改變橢圓鼻水翼吸力面的壓力分布,且趨勢都是類似的,都會在擋板的前端讓壓力突增,而在擋板后的一小塊區域內形成低壓區。為了比較不同布置擋板的升阻力變化,將計算得到的升阻力系數在表2中給出。

表2 不同擋板布置下的平均升阻力系數Tab. 2 The calculated average lift and drag coefficients under different configurations
布置a就是不加擋板的工況,布置b相對于它來說,升阻力系數都有提高。遺憾的是,因為橢圓鼻水翼有一定的攻角,升力系數的增幅并沒有遠高于阻力系數的增幅,但可以肯定的是在空化區設置擋板的確有效地提高了升力系數。在試驗中,橢圓鼻水翼的最大空化長度為0.5倍弦長,布置c的擋板設置在了空化充分發展的區域以外,這種情況就如同沒有空化產生時在水翼表面添加了多余的附件,不僅阻力系數大幅增加,升力系數也因為擋板未對空化區域產生影響而減小,這種布置的擋板在實踐應用中應該盡可能的避免。布置d和布置b一樣都將擋板設置在了0.37倍弦長處,但其長度只有展長的十分之一,因為擋板設置在空化區域內,有效地提高了升力系數,同時又不會過多的增加阻力系數,根據升阻比來看,這種布置的擋板是最有希望提高水翼效率的一種布置方式。
通過對添加了擋板的平板水翼和橢圓鼻水翼進行數值模擬,研究了不同布置的擋板對空化回射流和空化脫落的影響。在數值模擬中,采用一種新修正的穩定化SST k-ω湍流模型,將黏度公式重新定義,從而不至于過度預測湍動能和渦黏度。研究結果表明:
1) 在空化區域內適當地布置擋板,可以改變水翼吸力面的壓力分布,擋板下游形成的低壓區可以使水翼上下表面的壓差增大,從而提高升力系數。當水翼的攻角較小時,合適的擋板布置可以以較小的阻力代價換取巨額的升力提升,從而提高升阻比,提高水翼的效率,這在工程應用中有著巨大的潛能。
2) 擋板的存在可以有效地減小片空泡的脫落和云空化的產生,在片空化底部生成的回射流受擋板的阻礙后將改變流向,不再逆流而上。與擋板相撞后的回射流不再有繼續向上的沖勁,會順著主流,隨脫落的空泡一起翻卷流向下游,并在擋板的后方形成一個渦。
3) 擋板的布置需要依據不同的工作條件來定,把擋板設置在充分發展的空化區域以外,就如同沒有空化現象產生時在水翼表面添加了多余的附件,不僅阻力系數大幅增加,升力系數也會減小。即使擋板設置在充分發展的空化區域當中,也要考慮擋板的形狀和擋板的長度,一個工程中合適的擋板需要綜合多方面因素來選取。