鄭漢文，劉 婧，賴桂樺，曾永順，姚志峰，3

(1.中國農(nóng)業(yè)大學水利與土木工程學院，北京 100083；2.中國航發(fā)商用航空發(fā)動機有限責任公司，上海 200241；3.中國農(nóng)業(yè)大學北京市供水管網(wǎng)系統(tǒng)安全與節(jié)能工程技術研究中心，北京 100083)

伴隨能源需求量的增加，大量揚程高、尺寸大、轉速高的水力機械被投入能源生產(chǎn)，機組運行的穩(wěn)定性問題日益突出.在水力機械工作的過程中，當液壓降低至汽化壓力時，葉片前緣可能出現(xiàn)附著型空化區(qū)域[1]，且該區(qū)域內(nèi)空泡也可能脫離前緣向下游延伸，這種現(xiàn)象稱為“前緣空化”[2].前緣空化是水力機械中主要空化類型，空化可侵蝕結構表面，引起結構高幅值振動[3].前緣空化一般發(fā)生在層流邊界層分離后，該區(qū)域壓力分布與空泡長度有密切關系[4].前緣空化末端汽液混合，空泡與水流發(fā)生強烈相互作用[2]，表現(xiàn)出高度非定常特征[5].

大量學者通過實驗研究了前緣空化的發(fā)展形態(tài)、流場特性以及與水翼結構相互作用[6-8].Kawanami[6]等研究了繞水翼云狀空泡的脫落機制，認為空泡斷裂與脫落由空泡尾部反向射流引起.Wu[7-8]等探究了彈性水翼空化，發(fā)現(xiàn)空化發(fā)生會顯著增加彈性水翼的形變，同時彈性水翼動態(tài)特性也會改變空化發(fā)展的進程.

由于實驗在測量空化流場方面存在局限，數(shù)值模擬逐漸成為空化研究的重要手段[9].大量學者對數(shù)值模擬算法進行了研究，探討了空化數(shù)值模擬中湍流模型與空化模型的改進方法.Geng[10]等通過改變空化模型中凝結系數(shù)和蒸發(fā)系數(shù)參數(shù)，對比分析水翼空化流場特性后得出當凝結系數(shù)為0.08至0.1、蒸發(fā)系數(shù)為300至500時實驗與數(shù)值模擬較為吻合.

目前，求解湍流的方法分為直接求解和間接求解兩大類.直接求解方法(DNS)可以得到湍流流場的精確信息，但對計算資源的要求非常高，在實際工程中往往不采用直接求解方法.間接求解方法可分為三大類：第一類是系綜平均法，主要代表是雷諾平均法(RANS)，RANS可以計算高雷諾數(shù)的復雜運動，但計算的結果一般都是時均值，不能反應湍流的細節(jié)信息；第二類是空間篩濾法，主要代表為大渦模擬法(LES)，LES基于湍動能傳輸機制，直接計算大尺度渦流的運動，小尺度渦流的運動對大尺度渦流的影響通過建模體現(xiàn)出來，可以得到比RANS方法更多的動態(tài)信息；第三類為混合方法，主要代表為分離渦方法(DES)，它兼有RANS和LES的特點[11].

為了更加明確空化形態(tài)定量分析方法，在采用相同空化模型的條件下，本文對比了RANS、LES和DES三種湍流?；椒ㄔ诳张菪螒B(tài)預測上的差異.通過實驗數(shù)據(jù)，明確了空化模型中蒸發(fā)系數(shù)和凝結系數(shù)參數(shù)，定量分析了空泡形態(tài)與實驗結果的差異，為工程上預測前緣空化流場提供指導.

1 湍流計算方法

1.1 雷諾平均(RANS)法

雷諾平均法(RANS)廣泛使用在工程中，其原理是采用湍流統(tǒng)計理論，將非穩(wěn)態(tài)N-S方程作時間平均，求解工程中需要的時均量[12].該方法將流場中的變量分解為平均值和脈動值，并在此分解的基礎上求解均化后的運輸方程[13].雷諾平均后的N-S方程[14]：

(1)

(2)

(3)

(4)

公式中：ut為湍動粘度；k為湍流動能；δij為判定符號，當i=j時，δij=1，當i≠j時，δij=0.

根據(jù)湍動粘度ut求解方程的數(shù)量，可將渦粘模型分為零方程模型、一方程模型、兩方程模型.目前在工程應用中兩方程模型使用最為廣泛，主要有標準k-ε模型、RNGk-ε模型、k-ω模型和SSTk-ω模型等[13].

在k-ω兩方程模型中，湍動粘度ut的數(shù)值與湍動能k和比耗散率ω相關，其關系如下式所示

(5)

湍動能k方程和比耗散率ω方程分別為

(6)

(7)

公式中：β′=0.09，α=5/9，β=0.075，σk=2，σw=2為常數(shù)；密度ρ、速度矢量ui可以通過求解N-S方程得到；pkb及pωb為湍動能浮力生成項；pk為湍動能生成項[16]，其計算公式

(8)

基于k-ω的SST湍流模型將湍流切應力的運輸考慮在內(nèi)，并對逆壓梯度下流動分離的起始位置和數(shù)量有精確的預測.該湍流模型通過公式(9)定義渦粘系數(shù)，以得到湍流切應力為

(9)

公式中：νt=μt/ρ；S為應變率的不變量度；F2為混合函數(shù)，用以將該式的運用限制在邊界層內(nèi)，其表達式為

(10)

(11)

湍流模型的正確使用依賴混合函數(shù)，其表達式取決于計算點離最近的固體表面的距離以及流場流動參數(shù)，式中y即為距固體表面最近距離.

1.2 大渦模擬(LES)方法

LES湍流計算方法的基本原理是通過濾波函數(shù)將瞬時湍流運動分分解為大尺度和小尺度兩種類型.對大尺度運動采用直接求解的方法進行計算，對過濾的小尺度運動建立亞格子尺度(Subgrid-Scale，簡稱SGS)應力模型進行計算，以體現(xiàn)小尺度運動對大尺度運動的影響.LES模型未進行時均化處理，只對空間進行濾波，計算精度上相對于RANS模型具有一定優(yōu)勢，同時LES模型較直接求解可節(jié)省大量計算資源.

(12)

(13)

公式中：D為流體域；G為濾波函數(shù)；V為控制體；x′為實際流體域中的坐標值；x為濾波后大尺度空間中的坐標值.

不可壓縮流體濾波后的連續(xù)性方程及動量方程如公式(14)和公式(15)所示；其中τij為亞格子尺度應力，用以體現(xiàn)小尺度渦的運動對大尺度渦的影響，其表達公式為

(14)

(15)

(16)

在LES湍流模型求解過程中，大尺度渦流被直接求解，小尺度渦流則采用亞格子模型.CFX中采用渦流粘度方法，通過下式將亞格子尺度應力與大尺度應變率張量相關聯(lián).

(17)

(18)

(19)

公式中：Ls為亞格子尺度混合長度，計算公式為Ls=CSΔ；CS為Smagorinsky常數(shù)，該參數(shù)的取值對模型有較大影響，當取值接近0.1時能夠較好模擬大多數(shù)流動；Δ為當?shù)鼐W(wǎng)格尺度.

本文所采用的LES WALE(Wall-adapted local eddy viscosity)模型，相對于Smagorinsky模型在層流剪切層的求解更為準確，其渦粘系數(shù)μt采用公式(20)求解.

(20)

(21)

1.3 分離渦模擬(DES)方法

求解高雷諾數(shù)的邊界層流動，在使用LES湍流計算方法成本過高，使用RANS湍流計算方法不能滿足精度要求的情況下，工程中常采用DES湍流計算方法.DES湍流計算方法將RANS和LES湍流計算方法的特點結合在同一個混合方程中，即在流動邊界層內(nèi)未發(fā)生分離或流動分離較為輕微時使用RANS求解，而在流動分離劇烈的區(qū)域使用LES求解.在SST-DES模型中，當RANS模型計算得到的湍流長度Lt大于該區(qū)域內(nèi)網(wǎng)格尺度時，計算由SST模型轉換為LES模型[15].此時，湍動能耗散項方程的長度尺度變?yōu)樵搮^(qū)域的網(wǎng)格尺度Δ

ε=β*kω=k3/2/Lt→k3/2/(CDESΔ)for(CDESΔ

(22)

Δ=max(Δi)，

(23)

(24)

k方程變化為

(25)

公式中：CDES=0.61.

2 計算方案

2.1 計算模型

以二維NACA 0009鈍型尾部形狀水翼作為本文計算對象，研究水翼前緣空化流場.翼型弦長L為100 mm，尾部厚度h為3.22 mm，如圖1所示.

圖1 二維NACA 0009鈍型尾部形狀水翼模型

流體計算域，如圖2所示.水翼前緣距測試段進口2.5L；水翼尾部距測試段出口4L；水翼安放位置距測試段頂部0.7L，攻角為2.5°.在距水翼尾部中點延長線上距水翼尾部10 mm處設置監(jiān)測點，并定義x方向為流體流動方向，y方向為垂直于流動的方向，得到如圖1的坐標系.

圖2 二維計算域

2.2 網(wǎng)格劃分

空化計算中在流速20 m/s、攻角2.5°、空化數(shù)σ=0.81的條件下不同湍流模型的網(wǎng)格劃分方案和邊界條件，如表1所示.為消除計算域網(wǎng)格劃分方式對計算結果的影響，本文以水翼所受升力級及壓力為關鍵參數(shù)，采用基于理查德森外推加速法的網(wǎng)格收斂性指數(shù)[17]對上述網(wǎng)格劃分方案進行可靠性分析.結果表明，RANS型、DES湍流模型、LES湍流模型的收斂性指數(shù)分別在13%、1%、23%以內(nèi)，均滿足計算要求.

表1 不同湍流計算方法網(wǎng)格數(shù)及邊界條件(流速20 m/s、攻角2.5°、空化數(shù)0.81)

2.3 計算設置及邊界條件

采用ZGB空化模型進行空化數(shù)值計算，選擇Water at 25 ℃及Water Vapour at 25 ℃為測試段中水及水蒸氣介質，其中，水的密度ρ=997 kg/m3，水蒸氣密度為ρ=0.023 08 kg/m3.

RANS湍流計算方法瞬態(tài)項采用二階歐拉后差分格式，對流項及湍流模型相關量均采用高階精度求解；DES湍流計算方法瞬態(tài)項采用二階歐拉后差分格式，對流項采用中心差分格式，湍流模型相關量采用高階精度求解；LES湍流計算方法瞬態(tài)項采用二階歐拉后差分格式，對流項采用高階精度求解.

計算域進口邊界條件為速度進口，來流速度20 m/s；出口邊界條件為壓力出口，壓力值與Dupont[18]實驗中測試段空化數(shù)σ=0.81保持一致；計算域中上下壁面及水翼表面采用忽略粗糙度的無滑移壁面.空化計算前先進行流場的定常計算，收斂后，在定常計算的基礎上引入空化模型進行空化流場的非定常計算.非定常計算設置每個時間步內(nèi)最大迭代次數(shù)為50次，收斂殘差標準為1×10-4.

3.4 空化模型凝結系數(shù)、蒸發(fā)系數(shù)調整

結合實驗數(shù)據(jù)進行大量計算后得出，對于RANS湍流計算方法，當凝結系數(shù)為0.15、蒸發(fā)系數(shù)為140時，水翼表面壓力系數(shù)與實驗值較為吻合；對于DES湍流計算方法，當凝結系數(shù)為0.3，蒸發(fā)系數(shù)為140時，水翼表面壓力系數(shù)與實驗值較為吻合；對于LES湍流計算方法，當凝結系數(shù)為0.3，蒸發(fā)系數(shù)為30時，水翼表面壓力系數(shù)與實驗值較為吻合.后文將采用這組參數(shù)進行空化模擬.

3 計算結果與分析

3.1 前緣空化形態(tài)定量分析

(26)

(27)

Dupont[18]實驗中水翼NACA0009在2°攻角、20 m/s流速、空化數(shù)σ=0.81時的空化區(qū)域，如圖3所示.其中空泡長度lcavity=0.029 4 m，空泡高度hcavity=0.002 6 m.下文將以此數(shù)據(jù)作為衡量對象分析不同湍流模型對前緣空化形態(tài)的影響.

圖3 水翼前緣空化區(qū)域長度及高度定義

分別取三個湍流模型的一個運動周期T，截取該周期內(nèi)6個時刻對應的空化區(qū)域圖像算出空泡長度和空泡高度，并通過公式(26)、公式(27)求得該周期內(nèi)的lcavity和hcavity，結果如表2所示.計算結果顯示，與RANS和DES湍流模型相比，LES湍流模型計算所得空泡長度lcavity和空泡高度hcavity與實驗值更為吻合；對于RANS和DES湍流模型，空泡長度的最大相對誤差為16.8%，空泡高度的最大相對誤差為26.7%.

表2 不同湍流計算方法預測空泡長度及高度(流速20 m/s、攻角2.5°、空化數(shù)0.81)

在流速20 m/s、攻角2.5°、空化數(shù)0.81條件下，一個周期內(nèi)3個時刻RANS湍流模型、DES湍流模型和LES湍流模型的前緣空化流場，如圖4～圖6所示.對比分析三種湍流模型的空化區(qū)域和繞流流場發(fā)現(xiàn)：(1)RANS湍流計算方法模擬的空泡形態(tài)不隨時間改變，在空泡尾部x/L在0.22～0.3范圍內(nèi)存在明顯的回流區(qū)域；(2)DES湍流計算方法模擬的空泡在其尾部存在脫落區(qū)域，空泡脫落后不會再向流動方向移動，并且x/L在0.24～0.3這個區(qū)域內(nèi)，空泡的長度和高度會隨時間有一個較小的改變；(3)LES湍流計算方法可以模擬空泡從初生經(jīng)過生長到脫落的整個過程，在這個過程中空泡的長度和高度會隨著時間有一個較大的變化，x/L在0.08～0.31范圍內(nèi)均存在回流現(xiàn)象并且回流區(qū)域的位置會隨空泡形態(tài)的改變而改變.

圖4 RANS湍流計算方法空泡形態(tài)及繞流流場

圖5 DES湍流計算方法空泡形態(tài)及繞流流場

圖6 LES湍流計算方法空泡形態(tài)及繞流流場

3.2 近壁區(qū)速度分布

為探究不同湍流模型對前緣空化近壁區(qū)速度的影響，在水翼相對弦長x/L=0.1、0.4、0.7、0.99位置處作垂直于x軸的采樣線，每條采樣線上均布100個采樣點，測取3個流動周期內(nèi)每個采樣點對應的x方向速度和y方向速度，并將所測速度值求平均得到每點的平均速度Caverage.將速度無量綱化作為橫坐標，C=Caverage/Cinlet，其中Cinlet為計算域進口處速度；將y坐標與水翼弦長L之比作為縱坐標.將三種湍流模型對應的無量綱速度和實驗值進行比較[18]，如圖7所示.圖7(a)給出了采樣線段處的空泡的平均高度.

圖7 不同湍流模型水翼邊界層速度(流速20 m/s、攻角2.5°、空化數(shù)0.81)

對比x方向速度的實驗值和計算值發(fā)現(xiàn)：在x/L=0.4和x/L=0.7時三種湍流模型給出的預測值均與實驗值相差較大；在y/L>0.1時，三種湍流計算方法給出的預測值均與實驗值吻合的較好，但在y/L<0.1區(qū)間內(nèi)，預測值均顯著高于實驗值.

對比y方向上速度的實驗值和計算值發(fā)現(xiàn)：在x/L=0.4和x/L=0.7時，實驗值中存在沿y軸正向的速度分量，即存在速度大于0的區(qū)域，但從圖中可以看出RANS、DES、LES湍流模型均未準確的預測到這一現(xiàn)象；在y/L>0.1時，三種湍流計算方法給出的預測值與實驗值吻合得較好.

當采樣點與壁面間的法向距離超過0.1，三種湍流模型計算所得x方向速度和y方向速度均與實驗值較為吻合，即不同湍流模型在遠離壁面的區(qū)域對流場的模擬較為準確.由于本文空化計算的假設是均相流，計算結果反映的是均相速度，與真實值存在差異，因此在近壁區(qū)，RANS、DES和LES湍流模型計算的邊界層速度與實驗結果存在較大差異，不能準確反映近壁區(qū)速度分布.

三種湍流模型下同一流動周期內(nèi)水翼表面時均速度分布云圖，如圖8所示.顯示了水翼前緣空化流場中速度的分布規(guī)律.三種湍流模型對應前緣空化流場中的速度分布規(guī)律基本相同：空化區(qū)域均存在于水翼前緣且該區(qū)域的流速顯著高于該流場中的其它區(qū)域；水翼的前端和尾部均具有速度較低的的區(qū)域：前端的低速區(qū)成因是來流的撞擊，尾部的低速區(qū)成因是漩渦脫落；水翼的尾部存在因漩渦脫落形成的速度為負的區(qū)域.

圖8 不同湍流模型水翼表面時均速度分布云圖(流速20 m/s、攻角2.5°、空化數(shù)0.81)

三種湍流模型對應空化流場速度分布存在差異：在水翼的上表面，RANS和DES湍流計算方法速度為負的區(qū)域集中，而LES湍流計算方法速度為負的區(qū)域分散；RANS和DES湍流計算方法回流區(qū)域固定，LES湍流計算方法回流區(qū)域隨時間變化，說明LES湍流計算方法的出的空泡脫落位置隨時間變化，可以反映出空化的非定常特性.

4 結 論

本文以二維NACA 0009鈍型尾部形狀水翼為計算對象，在流速為20 m/s的條件下，分別用RANS、DES和LES三種湍流模型對前緣空化流場進行了數(shù)值模擬，討論了不同湍流模型在水翼前緣空化數(shù)值模擬中的應用特性.主要結論如下：

(1)RANS湍流計算方法計算所得空泡穩(wěn)定附著在水翼前緣，無空泡脫落現(xiàn)象；DES湍流計算方法在空泡尾部存在小范圍的脫落現(xiàn)象；LES湍流計算方法則可模擬得到空泡從初生到脫落并向下游移動的完整現(xiàn)象.

(2)在遠離水翼壁面的區(qū)域，RANS、DES和LES這三種湍流模型計算所得速度都與實驗吻合較好，但在近壁區(qū)這三種湍流模型計算所得速度均與實驗值存在較大差異.

(3)不同湍流模型水翼表面時均速度分布規(guī)律基本相同，僅在回流區(qū)域和負速度區(qū)域的分布上存在較小差別.