不同湍流模型在水翼前緣空化流場數值計算中的應用研究

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

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

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

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

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

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

為了更加明確空化形態定量分析方法，在采用相同空化模型的條件下，本文對比了RANS、LES和DES三種湍流模化方法在空泡形態預測上的差異.通過實驗數據，明確了空化模型中蒸發系數和凝結系數參數，定量分析了空泡形態與實驗結果的差異，為工程上預測前緣空化流場提供指導.

1 湍流計算方法

1.1 雷諾平均(RANS)法

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

(1)

(2)

(3)

(4)

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

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

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

(5)

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

(6)

(7)

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

(8)

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

(9)

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

(10)

(11)

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

1.2 大渦模擬(LES)方法

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

(12)

(13)

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

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

(14)

(15)

(16)

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

(17)

(18)

(19)

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

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

(20)

(21)

1.3 分離渦模擬(DES)方法

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

ε=β*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處設置監測點，并定義x方向為流體流動方向，y方向為垂直于流動的方向，得到如圖1的坐標系.

圖2 二維計算域

2.2 網格劃分

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

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

2.3 計算設置及邊界條件

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

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

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

3.4 空化模型凝結系數、蒸發系數調整

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

3 計算結果與分析

3.1 前緣空化形態定量分析

(26)

(27)

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

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

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

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

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

圖4 RANS湍流計算方法空泡形態及繞流流場

圖5 DES湍流計算方法空泡形態及繞流流場

圖6 LES湍流計算方法空泡形態及繞流流場

3.2 近壁區速度分布

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

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

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

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

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

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

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

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

4 結 論

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

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

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

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