水翼空化流場的調控和優化研究

趙偉,劉真威,陳石陽,丁耀東,李平

(西安交通大學熱流科學與工程教育部重點實驗室,710049,西安)

在水力機械運行過程中,當局部壓力降低為飽和蒸汽壓時,液相變為氣相產生空化。空化會在葉輪旋轉過程中周期性生成和潰滅,引起水力機械的振動和噪聲[1-2]。空化潰滅時形成的沖擊波會對設備表面產生高度集中的表面應力,附著型空化的周期性脫落與不穩定流動會導致表面應力的反復作用,造成表面材料的侵蝕,影響設備的安全運行[3-4]。因此,改善空化流場,抑制和削弱空化對水力機械產生的破壞,對于提高效率、保持水力機械安全穩定的運行至關重要。

NACA翼型在水力機械中獲得大量應用,但NACA翼型的開發主要以飛機機翼為應用目標,沒有考慮空化問題。因此,在保持水動力性能的基礎上,有必要設計符合水力機械抗汽蝕性能的NACA翼型。常欣等[5]對三維水翼型線進行了優化設計,優化后翼型前緣厚度減小并呈流線型,水翼升阻比提高。黃勝等[6]利用多目標優化算法,以提高升阻比和改善水翼表面壓力分布為優化目標,對不同翼型的型線進行了改進,改進后水翼的厚度減小,升力效率和空泡性能提高。黃斌等[7]利用粒子群優化算法對NACA66(MOD)水翼型線進行了改進,獲得的最優方案中,最大拱度位置向尾緣移動,最大厚度略微增加。李靖璐等[8]采用自由變形法對NACA0012水翼進行了型線調整,研究發現,水翼厚度減小可使阻力系數降低,后緣扭轉可使升力系數增加。Luo等[9]利用NSGA-Ⅱ算法對水輪機葉片翼型進行了優化設計,改進后水翼前緣附近的曲率半徑增大,翼型厚度減小,升阻比和最小壓力系數增加,運行效率提高,減小了空化發生的可能性。上述研究結果表明,改進水翼型線可以實現空化流場的整體調控,并改善水翼的水動力性能。

以流場中空化發生的局部位置及后續發展趨勢為著眼點,其調控方法主要有主動式控制和被動式控制。被動控制方法中,布設流動控制結構是一種易于實施和維護的控制方式。Kadivar等[10]在水翼表面分別布置了圓柱形和半球形渦發生器,研究發現在合適位置布置渦發生器可以抑制空泡脫落,減小水翼尾跡區的壓力脈動。顧魏等[11]分別在35%弦長和75%弦長位置處布置單個展向擋流條,發現35%弦長處的擋流條在一定的空化數范圍內增強了空化的穩定性。Qiu等[12]在NACA0015水翼表面布置微渦發生器,發現微渦發生器對流體進行擾動,促進了不同能量流體的混合,形成穩定的渦流區,使空泡破滅時劇烈程度降低,減小了對水翼表面的空蝕。以上研究表明,流動控制結構可以實現局部空化流場的有效控制。在空化范圍較大時,單個結構的控制作用有限,因此連續流動控制結構成為了優選方法[13-15]。

在布設多個流動控制結構時,其幾何參數和相對位置會顯著影響空化控制效果,而將優化算法耦合到數值求解和參數化分析過程中,則可以快速得到理想的控制結構參數。Lin等[16]將簡化的共軛梯度法與數值模擬方法相結合,對微通道中肋的寬度和厚度等結構參數進行了優化設計,獲得不同工況下最優設計變量組合。Mazaheri等[17]利用差分進化優化算法對翼型表面的凸起結構進行了優化設計,優化后的翼型阻力和流動分離區域減小,增強了流場的穩定性。Xue等[18]利用遺傳算法對凹槽結構的數量和深度等進行了全局優化,獲得的最優凹槽結構改善了凹槽區域的壓力分布,抑制了空化的產生。以上研究表明,耦合優化算法能夠為空化性能和水動力性能的優化提升提供有效的快速實現方法。

綜上分析:本文通過改進水翼型線,進行空化流場的整體優化;通過布設連續流動控制結構,并基于廣義模式搜索算法進行結構優化,快速實現空化流場的局部精準調控。

1 數值計算方法

1.1 物理模型和邊界條件

在水翼空化流場研究中,本文選擇NACA0009、NACA0012和Clark Y 3種典型水翼型線進行對比分析。NACA0009和NACA0012為低速對稱翼型,彎度均為0,具有低阻力系數,是水力機械中的基礎翼型;Clark Y水翼為非對稱翼型,具有高升力系數。圖1展示了NACA0009、NACA0012和Clark Y 3種水翼型線的對比,3種水翼弦長相等,均為0.1 m。其中,NACA0009水翼最大厚度為9.90 mm,位于49.5%弦長處;NACA0012水翼最大厚度為12.0 mm,位于30%弦長處;Clark Y水翼最大厚度為11.7 mm,位于28%弦長處,吸力面最大凸起高度為9.16 mm,位于36%弦長處,最大彎度位于42%弦長處。

圖1 3種水翼型線對比

計算域如圖2所示,采用速度入口和靜壓出口,水翼表面和上下壁面為無滑移壁面。流體溫度為20℃,在此條件下飽和蒸汽壓為2 300 Pa,水的密度為998 kg·m-3,蒸汽密度為0.017 kg·m-3,水的動力黏度為0.001 kg·m-1·s-1,蒸汽動力黏度為9.60×10-6kg·m-1·s-1。在入口速度為10 m·s-1條件下,針對攻角為6°和8°,空化數為1.5、1.3、1.0和0.8等工況時,對比分析不同水翼的流場特性。

圖2 水翼空化模型計算域

1.2 數值計算方法

1.2.1 基本控制方程

本文研究采用均質混合流模型[19],其中液相和氣相具有相同的速度和壓力,空化流動非穩態求解控制方程如下。

混合相的連續性方程為

(1)

(2)

(3)

其中混合相的密度和黏度分別為

ρ=αvρv+ρl(1-αv)

(4)

μ=αvμv+μl(1-αv)

(5)

式中:ρ為混合相密度;ui(i=1,2,3)表示與坐標軸xi平行的速度分量;αv為氣相體積分數;下標l代表液相,下標v代表氣相;m為兩相之間的傳質速率。

1.2.2 湍流模型

本文采用剪切應力輸運(SSTk-ω)湍流模型[20]對空化流動進行計算,該模型如下

(6)

(7)

式中:k為湍動能;Pk為湍流產生率;ω為湍流頻率;σk=0.85;σω1=0.5;σω2=0.856;F1為混合函數。

由于原始的SSTk-ω模型過度預測了空化區域的湍流黏度,因此Ducoin等[21]對SSTk-ω模型進行了修正,考慮了混合相局部可壓縮性的影響,提出了基于局部氣相體積分數的混合湍流黏度求解方法,將原公式中的湍流黏度μt修正為

μt_mod=μtf(ρ)

(8)

(9)

(10)

式中:a1=0.31;S為液體表面張力系數;F2為混合函數。通過調節n降低空化區域的密度,從而減小湍流黏度,參照文獻[21-22],當n取3時,對空化區域湍流黏度的預測滿足本研究的精確度要求。

1.2.3 空化模型

Zwart空化模型在非穩態空化模擬中具有較好的效果,自提出至今,進行了大量驗證[23]。考慮Rayleigh-Plesset方程的二階項對空化數值模擬的影響,Geng等[24]通過數學推導對Zwart空化模型進行了修正,得到更接近實驗的計算結果,修正后的空化模型為

(11)

式中:me為蒸發傳質速率;mc為冷凝傳質速率;αnuc為汽核體積分數,取5×10-4;R為1×10-6m;Fcond為凝結系數,取0.01;Fvap為汽化系數,取50。

1.2.4 參數定義

空化數σ是描述空化狀態的無量綱參數,其定義為

(12)

式中:p∞為出口壓力;V∞為來流速度。

壓力系數定義為

(13)

反映水翼水力特性的升力系數、阻力系數分別定義為

(14)

(15)

式中:Fl、Fd分別為水翼所受的升力和阻力;A為翼型投影面積,二維翼型時為翼型弦長。

1.3 模型驗證

1.3.1 網格無關性驗證

本文中對流體域進行結構化網格劃分,并對水翼近壁區進行網格加密,調整第一層網格高度,確保離壁面最近的第一個點位于黏性底層內,基于下式保證y+<1,網格劃分如圖3所示。

圖3 水翼周圍及兩端網格示意圖

(16)

式中:y為單元中心到壁面的距離;τw為壁面切應力。

為確定合理的網格數,對網格進行無關性檢驗,保持第一層網格高度和增長率不變,不斷對網格進行加密,最終網格數分別為66 048、94 748和119 608。以升力系數、阻力系數作為評價標準,結果見表1。網格3與網格2相比,參數相對誤差均低于0.1%,因此,以網格3為基準,進行后續水翼網格的劃分。

表1 網格無關性驗證

1.3.2 數值模型驗證

以穩態空化流場計算結果作為初始值,對非穩態空化流場進行計算。將NACA0009水翼在入口速度為20 m/s、攻角為2.5°、空化數為0.9時,吸力面壓力系數的計算值與Dupont[25]的實驗數據進行對比,驗證結果如圖4所示,可以看出,時均壓力系數與實驗值吻合較好。

圖5給出了Clark Y水翼在入口速度為10 m/s、攻角為8°、空化數為0.8時,計算所得空泡形態隨時間的變化及其與實驗結果[26]的對比,數值計算結果清楚地描述了云狀空化的產生、發展和脫落的準周期性變化,模擬空泡形態與實驗結果吻合度較高。

圖5 Clark Y水翼云空化空泡演變

綜上可見,本文計算模型可以有效模擬水翼流場及其空化特性,模型精度足夠可靠,滿足下一步研究需求。

2 水翼空化流場整體調控

2.1 流線分布和水動力性能

圖6和圖7為不同空化數下,一個周期內空化流場隨時間的變化。水翼云空化的發展過程可歸納為:空泡從翼型前緣開始產生,并逐步發展成附著在水翼表面的片狀空泡;隨著回射流向前緣方向運動,空泡尾端被拉伸,并逐漸脫落;脫落的空泡隨流體向翼型尾緣方向移動,隨著壓力增大,空泡逐漸變小,最后在翼型尾緣附近消失。由于空泡的大量脫落和潰滅,造成壓力脈動和沖擊,對水翼壁面產生破壞。隨著空化數的降低,流場壓力降低,附著空泡的長度和厚度增加。最終,空泡發展的長度超過水翼尾緣,變成超空化狀態。在相同空化數和不同攻角下,同一時刻空泡形態相似,但較大攻角下的空泡體積略大。

圖6 α=6°和σ=1.5時氣體體積分數分布

圖7 α=6°和σ=0.8時氣體體積分數分布

對不同水翼的不同工況進行了計算,圖8為不同空化數下各水翼周圍的流線分布。空泡區域屬于低壓區,當水翼表面的附著空化發展到一定程度,由于逆壓梯度的存在,流體沿水翼吸力面反方向流動,形成回射流。回射流所在區域速度與主流速度方向相反,對空泡產生反作用,同時與空泡體接觸部分不斷汽化,在相互作用的過程中,區域整體的運動方向與主流方向相同,形成回射流“倒退”現象。空化數較小時,水翼表面附著空化發展到尾緣附近,導致尾緣吸力面壓力降低,與壓力面的壓差增大,流體從壓力面流向吸力面尾緣并形成回射流,因此回射流產生的位置在尾緣附近。空化數較大時,附著空化的長度較小,離尾緣有一定的距離,尾緣附近吸力面和壓力面的壓差較小,提供給尾緣壓力面流體流向吸力面的動力較小,因此尾緣產生的回射流較小。對比不同水翼的回射流分布可以看出,在較小空化數下,NACA0009和NACA0012水翼吸力面尾緣產生的回射流運動的更加深入,而相同空化數下Clark Y水翼吸力面尾緣產生的回射流運動的距離較短。比較3種水翼的型線可以看出,Clark Y水翼吸力面較大的凸起高度,改變了回射流運動的角度,使附著型空化更容易被剪斷,從而抑制了回射流向前緣的發展,抑制了空化的發展。因此,在保持吸力面較大凸起高度的條件下,將最大凸起高度位置沿來流方向適當后移,可以抑制回射流的發展,減小空泡的脫落。

圖8 α=6°時不同空化數下水翼流線分布

升阻比是評價水力機械水動力性能的重要參數。升阻比越高,意味著效率越高,同樣的工作條件下可以節省更多的能源。圖9給出了不同工況下3種水翼的升阻比大小。如圖9所示,隨著空化數的升高,各水翼的升阻比也逐漸升高。僅在攻角為6°、空化數為1.0時,NACA0009水翼升阻比大于Clark Y水翼升阻比,其余工況下Clark Y水翼的水動力性能均優于NACA0009水翼和NACA0012水翼。對于不同水翼,攻角為6°時的升阻比均大于攻角為8°時的升阻比。盡管攻角增大,升力系數也隨之增大,但阻力中的誘導阻力增加比重更大,因而最終的結果是升阻比減小。Clark Y水翼的升阻力特性與NACA0012水翼、NACA0012水翼的升阻力特性不同。如表2所示,對于NACA0012水翼,升力系數和阻力系數均小于Clark Y水翼,而阻力系數對升阻比的影響較大。因此,NACA0012水翼和具有高升力系數的Clark Y水翼在升阻比上差別較小。NACA0012水翼具有最小阻力系數,NACA0009水翼的阻力系數介于NACA0012水翼和Clark Y水翼之間。結合圖1中3種水翼的型線對比可以得出,翼型吸力面凸起高度過小和過大都會使阻力系數增大。因此,為保持高升力系數的同時擁有較小阻力系數,型線的凸起高度需要介于Clark Y水翼和NACA0012水翼的凸起高度之間。

表2 α=6°的時水翼升力系數、阻力系數

圖9 不同工況下的升阻比

2.2 新水翼結構與流場分析

為保證水動力性能的同時,實現對空化的整體調控,本研究依據多項式形式的型線方程設計新水翼型線[27]。在對比分析的3種水翼中,Clark Y水翼綜合性能最佳,通過Clark Y水翼的數據點擬合型線方程,確定多項式中的各項系數和冪指數,得到待優化的初始型線方程

y=pxa1(c-x)b1±qxc1(c-x)d1±

rxs1(c-x)t1

(17)

式中:p、q、r均為大于0的常數。對于上下型線,指數a1、b1、c1、d1、s1、t1可以取不同的值,但均為大于0的常數。

基于前文翼型對比分析得出的型線優化思路,以NACA0012水翼和Clark Y水翼型線作為新水翼型線優化的參數變化范圍。以改善水翼的空化性能和水動力性能為目標,通過對比不同的常數值,調整型線方程對型線進行優化,進而確立新水翼型線,方程式如下

(18)

式中:ys、yp分別為吸力面、壓力面型線方程。由型線方程建立新水翼模型,命名為NF水翼模型。

新水翼型線和現有3種水翼的對比如圖10所示。新水翼的最大厚度為11.3 mm,在33%弦長處,吸力面最大凸起高度為8.34 mm,在38.4%弦長處。新水翼吸力面凸起高度介于NACA0012水翼和Clark Y水翼之間,并且與Clark Y翼型相比,吸力面最大凸起高度的位置沿來流方向后移。

圖10 新水翼型線示意圖

通過翼型優化方案,建立了新的翼型結構。圖11展示了不同空化數下新水翼流線分布,可以看出,空化數較小時,回射流從尾緣產生并向前緣流動;空化數較大時,回射流產生的位置提前。與圖8各水翼流線分布相比,在相同的空化數下,新水翼吸力面回射流向前緣運動的距離比NACA0009水翼和NACA0012水翼小。雖然新水翼吸力面最大凸起高度小于Clark Y水翼,但新水翼吸力面最大凸起高度向后移動了一定的距離,所以回射流發展的距離相比Clark Y水翼增大不明顯。因此,利用優化方案建立的新水翼抑制了回射流的發展,減小了回射流對空泡脫落的影響。

圖11 α=6°時不同空化數下新水翼流線分布

新水翼與Clark Y水翼在攻角為6°和8°時的水動力性能如圖12所示,可以看出,在相同攻角和相同空化數條件下,新水翼的水動力性能優于Clark Y水翼,升阻比最大提升48.8%。這意味著,根據不同流場特性對比分析得出的型線優化思路,為提高水翼水動力性能提供了方向。新水翼升阻比變化規律與其他3種水翼相同,即升阻比隨空化數的升高而升高,且攻角為6°時的升阻比大于攻角為8°時的升阻比。另外,由于新水翼吸力面具有較大的凸起高度,所以新水翼阻力系數較大的同時具有較大的升力系數,升阻比特性與Clark Y水翼類似。

圖12 新水翼升阻比特性

3 新水翼流場局部調控

3.1 TR水翼及優化算法

為進一步對新翼型的空化流場進行局部調控,在新水翼上添加流動控制結構。參考空化流場中流動控制結構布置的原則[28-30]可知,控制結構離前緣過近或過遠都難以對空化起到有效的控制作用,而控制結構高度過高會對流場產生強烈擾動,對空化造成不良影響。因此,控制結構需選取合適的高度和間距,并且布置在水翼空化覆蓋范圍內。針對空化區域較大的σ為0.8和1.0兩種工況,本文在新水翼吸力面均勻添加3個三角形凸起控制結構,研究多個控制結構對空化流場的影響。新水翼表面添加多個控制結構的模型如圖13所示,命名為TR水翼模型。

圖13 TR水翼模型

利用給定的流動控制結構進行水翼空化調控研究時,存在遺漏流動控制結構最優設計參數的可能。為了實現流動控制結構設計的快速尋優,以及水翼性能的最大化提升,針對TR水翼,利用廣義模式搜索算法進行優化設計。在優化過程中,以流動控制結構起始位置C0、控制結構高度h以及相鄰控制結構的間距d作為設計變量,以代表水翼水動力性能的升阻比作為目標參數。為保證結構的正常添加,避免遠離水翼吸力面前緣,以及控制結構由于間距過小而互相干擾等問題,給出各設計變量[C0,h,d]變化的上下限,對控制結構的添加進行限制,下限為[0.1,0.03,0.003],上限為[0.7,0.3,0.01]。控制結構優化流程如圖14所示,利用Matlab調用相應的程序,自動完成幾何建模、網格劃分以及數值計算,并對結果進行處理。最后,利用Matlab中的優化算法完成目標參數的尋優和對比,并將調整的設計變量參數重新輸入,繼續下一次迭代。

圖14 添加控制結構的水翼優化流程圖

利用模式搜索算法對函數式(17)進行求解,驗證該優化方法的可靠性,得到

f=-x2-2y2+0.4cos(3πx)+

0.6cos(4πy)

(19)

式中:x,y∈[-10,10]。該函數為多峰函數,在[0,0]處具有最大值1。通過模式搜索方法進行最大值搜索,經過多次迭代之后,在x=-0.534 3×10-15,y=-0.104 1×10-15處得到最大值1,兩坐標值接近理論值0,由此驗證了優化程序的可靠性。

3.2 優化結果與分析

在優化設計中,將攻角為6°和8°時優化前的模型分別命名為TR6和TR8,同時將各工況優化后的模型按照以下規則進行命名:TRα-σ。基于新水翼流場分析結果,綜合考慮空泡和回射流產生的位置及大小等,選定本次優化的初始點為x0=[0.025 m,0.000 6 m,0.005 m],即控制結構起始位置C0在弦長25%位置處,控制結構高度h為0.6 mm,相鄰控制結構間距d為5 mm。表3為優化前后控制結構參數對比,從優化結果可以看出,優化后流動控制結構初始位置后移,同時控制結構的高度以及間距都增加。

表3 不同工況下流動控制結構優化結果

優化后水翼特性與NF水翼、TR水翼特性的對比如圖15和圖16所示,可以看出,相對于NF水翼,TR水翼升阻比減小,空化厚度Cy與長度Cx增加。通過優化設計后,水翼的升阻比增加,與NF水翼相比,TR6-0.8水翼升阻比提升了9.5%,TR6-1.0水翼升阻比提升了13.6%,TR8-0.8水翼升阻比提升了4.8%。同時,優化后水翼空化的發展得到了抑制,空化長度Cx與厚度Cy均有所減小。這意味著,使用優化設計程序實現了流動控制結構的優化,空化流場得到局部調控。

圖15 優化前后水翼空泡脫落頻率和升阻比

圖16 優化前后水翼空化厚度和空化長度

圖17和圖18分別給出了攻角為6°、空化數為0.8條件下,控制結構優化前后,即TR6水翼和TR6-0.8水翼的氣體體積分數與流線分布在一個周期內的變化,其中圖17(a)和圖18(a)為水翼周圍流場圖,圖17(b)和圖18(b)為控制結構附近局部放大圖。優化后的TR6-0.8水翼對空化的抑制作用更大,空化的長度和厚度明顯減小。與TR6水翼相比,TR6-0.8水翼的氣體體積分數減小,空化流場得到改善。從流線分布可以看出,在空泡初生時刻,水翼尾緣處回射流還未產生,控制結構處由于自身形狀原因產生渦結構。由于TR6-0.8水翼控制結構之間的間距大于TR6水翼,TR6-0.8水翼控制結構處的渦結構較大,因此TR6-0.8水翼控制結構對空化流場的改善作用更大。在空泡脫落時,由于回射流的切割作用,使空泡變的不穩定。從空泡脫落時控制結構處的流線分布可以看出,回射流受到3個連續控制結構的阻礙作用,向前緣的深入運動得到連續抑制。TR6水翼控制結構初始位置在弦長25%處,回射流繞過后面兩個控制結構,雖然控制結構對回射流起到了阻礙作用,但是TR6水翼的控制結構位置離前緣較近,不僅沒有起到抑制空化的作用,還導致空化長度和厚度增加。對于TR6-0.8水翼,由于控制結構的位置在弦長32%處,并且高度更大,回射流在繞過第二個控制結構到達第一個控制結構后,沒有再繞過第一個控制結構,充分抑制了回射流向前的深入運動。同樣的,優化后得到的TR6-1.0水翼和TR8-0.8水翼都對回射流產生了較大的抑制作用。通過布設連續流動控制結構,能夠實現空化流場的局部調控,但需要合適的控制結構參數及布局,通過優化算法,快速實現了這一過程。

(a)水翼周圍 (b)控制結構附近

(a)水翼周圍 (b)控制結構附近

4 結 論

本文在不同攻角和空化數的工況下,對NACA0009、NACA0012和Clark Y 3種水翼的流場進行了對比分析。通過改變型線,設計了新水翼,對空化流場進行整體控制。進一步地,結合廣義模式搜索算法在新水翼吸力面布設多個控制結構進行了優化設計,對空化流場進行局部調控。主要研究結論如下。

(1)在相同工況下,3種水翼中Clark Y水翼的水動力性能最優,同一時刻不同水翼產生的空泡形態相似。對于同一水翼,空泡體積會隨攻角的減小和空化數的增加而減小,相應的升阻比也會升高,水動力性能提升。

(2)為保持高升力系數的同時擁有較小阻力系數,新水翼型線的凸起高度需要介于Clark Y水翼和NACA0012水翼之間;在保持吸力面較大的凸起高度條件下,將最大凸起高度位置沿來流方向適當后移,可以抑制回射流的發展。

(3)新水翼的水動力性能優于NACA0009、NACA0012和Clark Y 3種水翼,升阻比明顯提高,最大提高48.8%;與NACA0009和NACA0012水翼相比,新水翼抑制了回射流的發展,減小了空泡脫落頻率,實現了對空化流場的整體調控。

(4)在新水翼吸力面布設連續控制結構,通過優化得到的TR6-0.8水翼升阻比進一步提高了9.5%,TR6-1.0水翼提高了13.6%,TR8-0.8水翼提高了4.8%,優化后水翼的水動力性能進一步提升,空化長度與厚度均有所減小,實現了對空化流場的局部調控。

(5)控制結構位置需要與水翼前緣具有一定的距離,如果控制結構離前緣較近,不僅不會起到抑制空化的作用,還會導致空化惡化。均勻布置的多個控制結構能夠對回射流起到連續抑制作用,減小回射流對空化區域的影響。