不同湍流模型在水翼前緣空化流場(chǎng)數(shù)值計(jì)算中的應(yīng)用研究

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

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

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

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

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

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

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

1 湍流計(jì)算方法

1.1 雷諾平均(RANS)法

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

(1)

(2)

(3)

(4)

公式中：ut為湍動(dòng)粘度；k為湍流動(dòng)能；δij為判定符號(hào)，當(dāng)i=j時(shí)，δij=1，當(dāng)i≠j時(shí)，δij=0.

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

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

(5)

湍動(dòng)能k方程和比耗散率ω方程分別為

(6)

(7)

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

(8)

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

(9)

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

(10)

(11)

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

1.2 大渦模擬(LES)方法

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

(12)

(13)

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

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

(14)

(15)

(16)

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

(17)

(18)

(19)

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

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

(20)

(21)

1.3 分離渦模擬(DES)方法

求解高雷諾數(shù)的邊界層流動(dòng)，在使用LES湍流計(jì)算方法成本過高，使用RANS湍流計(jì)算方法不能滿足精度要求的情況下，工程中常采用DES湍流計(jì)算方法.DES湍流計(jì)算方法將RANS和LES湍流計(jì)算方法的特點(diǎn)結(jié)合在同一個(gè)混合方程中，即在流動(dòng)邊界層內(nèi)未發(fā)生分離或流動(dòng)分離較為輕微時(shí)使用RANS求解，而在流動(dòng)分離劇烈的區(qū)域使用LES求解.在SST-DES模型中，當(dāng)RANS模型計(jì)算得到的湍流長(zhǎng)度Lt大于該區(qū)域內(nèi)網(wǎng)格尺度時(shí)，計(jì)算由SST模型轉(zhuǎn)換為L(zhǎng)ES模型[15].此時(shí)，湍動(dòng)能耗散項(xiàng)方程的長(zhǎng)度尺度變?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 計(jì)算方案

2.1 計(jì)算模型

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

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

流體計(jì)算域，如圖2所示.水翼前緣距測(cè)試段進(jìn)口2.5L；水翼尾部距測(cè)試段出口4L；水翼安放位置距測(cè)試段頂部0.7L，攻角為2.5°.在距水翼尾部中點(diǎn)延長(zhǎng)線上距水翼尾部10 mm處設(shè)置監(jiān)測(cè)點(diǎn)，并定義x方向?yàn)榱黧w流動(dòng)方向，y方向?yàn)榇怪庇诹鲃?dòng)的方向，得到如圖1的坐標(biāo)系.

圖2 二維計(jì)算域

2.2 網(wǎng)格劃分

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

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

2.3 計(jì)算設(shè)置及邊界條件

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

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

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

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

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

3 計(jì)算結(jié)果與分析

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

(26)

(27)

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

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

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

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

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

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

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

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

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

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

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

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

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

當(dāng)采樣點(diǎn)與壁面間的法向距離超過0.1，三種湍流模型計(jì)算所得x方向速度和y方向速度均與實(shí)驗(yàn)值較為吻合，即不同湍流模型在遠(yuǎn)離壁面的區(qū)域?qū)α鲌?chǎng)的模擬較為準(zhǔn)確.由于本文空化計(jì)算的假設(shè)是均相流，計(jì)算結(jié)果反映的是均相速度，與真實(shí)值存在差異，因此在近壁區(qū)，RANS、DES和LES湍流模型計(jì)算的邊界層速度與實(shí)驗(yàn)結(jié)果存在較大差異，不能準(zhǔn)確反映近壁區(qū)速度分布.

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

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

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

4 結(jié) 論

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

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

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

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