張德勝 馮 健 萬福來 李俊峰 施衛(wèi)東

(1.江蘇大學流體機械工程技術研究中心，鎮(zhèn)江212013;2.南通大學機械工程學院，南通226019)

0 引言

空化是一種復雜的兩相流動現(xiàn)象［1-3］，普遍存在于多種流體機械中，例如水輪機、透平機、螺旋槳等。當流體機械內部的局部流體壓力低于相應溫度下的飽和蒸汽壓時，就會產(chǎn)生空化?？栈瘯斐奢^大的振動和噪聲，損壞過流部件，加劇軸承的磨損，影響流體機械的使用壽命。

湍流模型在空化問題的研究中起著重要作用。目前，基于雷諾時均方程的模擬方法被廣泛應用于空化數(shù)值模擬中。文獻［4］采用SST k-ω 湍流模型進行數(shù)值模擬，并結合實驗分析了高速誘導輪離心泵內的空化特性。文獻［5］采用RNG k-ε 湍流模型分析了不同空化數(shù)下翼型的云狀空化脫落特性。近年來，隨著計算流體力學的快速發(fā)展，一些學者采用大渦模擬(Large eddy simulation，LES)研究瞬態(tài)流動特性［6-8］，大渦模擬可以直接模擬大尺度湍流流動，利用次網(wǎng)格尺度模型模擬小尺度湍流流動對大尺度湍流流動的影響，在實際工程中得到一定應用。目前，有關湍流模型對空化流模擬影響的研究主要以常溫水為工作介質，常溫水的飽和蒸汽壓強隨溫度變化梯度較小，其空化過程可視為等溫過程，而高溫水、液化天然氣、液氫等介質的飽和蒸汽壓強對溫度變化敏感，屬于熱敏介質，其空化過程不可視為等溫過程。一些學者著重分析了空化模型對熱敏介質空化流的影響［9-10］，而針對湍流模型對熱敏介質空化流模擬影響的研究較少。

本文結合FBM 模型和DCM 模型，對現(xiàn)有3 種湍流模型(k-ε、RNG k-ε 和SST k-ω)進行修正，并對不同溫度水繞NACA0015 翼型進行空化模擬，結合實驗對比驗證修正模型的數(shù)值模擬精度。

1 控制方程與數(shù)值方法

1.1 基本方程

本文采用均質平衡流模型模擬氣液兩相空化流動，該模型將氣液混合物的密度看成均勻單一的密度，并且混合物具有相同的流速和壓力?；究刂品匠贪ㄟB續(xù)性方程、動量方程和能量方程，依次表示為

其中

式中 ρm——混合介質密度，kg/m3

ρl——液相密度，kg/m3

ρv——氣相密度，kg/m3

u——速度，m/s

αv——氣相體積分數(shù)

p——壓力，Pa

μ——混合介質層流粘度，Pa·s

μt——混合介質湍流粘度，Pa·s

h——焓，J

L——液體的氣化潛熱，kJ/kg

fv——氣相質量分數(shù)

PrL——層流普朗特數(shù)

Prt——湍流普朗特數(shù)

t——時間，s

δij——克羅內克函數(shù)

下角標i、j、k 表示坐標方向。

1.2 湍流模型

湍流模型可以分成兩方程模型、一方程模型和零方程模型。本文主要針對兩方程模型中的3 種湍流模型(k-ε、RNG k-ε 和SST k-ω)進行修正和對比分析。

1.2.1 k-ε 湍流模型

Launder 和Spalding 在1972 年提出了k-ε 湍流模型，需要求解兩個附加方程，通過求解湍動能和湍動能耗散率來求解湍流粘度，以此得到湍流應力。

湍動能k 方程為

式中 σk——湍動能的湍流普朗特數(shù)Gk——湍流動能

湍動能耗散率ε 方程為

式中 σε——耗散率的湍流普朗特數(shù)經(jīng)驗系數(shù)Cε1、Cε2分別取1.44 和1.92。

湍流粘度μt方程為

式中 Cμ——經(jīng)驗系數(shù)

1.2.2 RNG k-ε 湍流模型

RNG k-ε 湍流模型［11］考慮了平均流動中的旋轉流動情況，在湍動能耗散率ε 方程中添加一項，反映了時均應變率Eij。

湍動能k 方程為

式中 αk——經(jīng)驗系數(shù)，取1.393

μeff——有效粘度，Pa·s

湍動能耗散率ε 方程為

其中

式中 αε——經(jīng)驗系數(shù)

經(jīng)驗系數(shù)分別為:αε=1.393，Cε1=1.42，Cε2=1.68，Cμ=0.084 5，η0=4.38，β=0.012。

1.2.3 SST k-ω 湍流模型

SST k-ω 湍流模型［12］結合k-ω 模型和k-ε 模型，近壁面處采用k-ω 模型計算低雷諾數(shù)問題，在遠壁處采用k-ε 模型計算。

湍動能k 方程為

式中 Pk——湍流生成速率

β′——經(jīng)驗系數(shù)，取0.09

湍流頻率ω 方程為

渦流粘度νt方程為

式中 S——固定應變率估計值

F1、F2——混合函數(shù)

經(jīng)驗系數(shù)分別為:α1= 0.3，α = 5/9，σω2=1/0.856。

1.3 湍流模型修正

文獻［13］考慮氣液兩相混合密度的可壓縮性對湍流粘度的影響，提出采用DCM 模型修正湍流粘度，公式為

其中

根據(jù)文獻［13］，n 取10。

文獻［14］提出一種FBM 模型，保證運輸方程不變的情況下，在湍流粘度中引入一個濾波函數(shù)fFBM，公式為

其中

式中 C3——經(jīng)驗系數(shù)，取1.0

Δ——濾波尺寸，其值取決于網(wǎng)格尺寸

結合濾波器模型與密度修正模型對湍流粘度進行修正［15］，公式為

其中

式中，χ 指函數(shù)符號，根據(jù)文獻［15］，經(jīng)驗系數(shù)C1、C2、C3、C4分別取4、0.2、0.6、0.2。

1.4 Merkle 空化模型

Merkle 空化模型［16］推導了氣液混合物中的相間傳質方程。氣相體積分數(shù)輸運方程和蒸發(fā)凝結源相表達式為

其中

式

中 Fe、Fc——蒸發(fā)、凝結系數(shù)，取1、80 u

in——進口速度，m/s

t∞——特征時間，s

pv——飽和壓力，Pa

1.5 數(shù)值設置與驗證

文獻［17］對不同溫度水繞NACA0015 翼型進行了系統(tǒng)的空化實驗研究，分析了熱力學效應對空化性能的影響。實驗用的NACA0015 翼型結構如圖1 所示。實驗段進出口各設有3 個測壓孔，翼型吸力面設有10 個測壓孔，壓力面設有2 個測壓孔。

考慮到流動的充分發(fā)展以及穩(wěn)定性，對計算域適當延伸，進口與出口分別延伸至2 倍與4 倍翼型弦長，翼型攻角為5°，三維模型如圖2 所示。

圖1 結構示意圖Fig.1 Schematic of test section

圖2 NACA0015 翼型三維模型Fig.2 3D model of NACA0015 hydrofoil

邊界條件設置與實驗對應:進口設置為速度入口，出口設置為壓力出口，壁面設置為無滑移邊界。不同溫度水下的進出口參數(shù)值如表1 所示。

表1 邊界條件Tab.1 Boundary conditions

本文采用ICEMCFD 軟件對NACA0015 翼型劃分六面體結構網(wǎng)格，計算域網(wǎng)格如圖3 所示。翼型邊界層網(wǎng)格質量對數(shù)值計算精度有較大影響，通常采用y+來驗證網(wǎng)格質量是否達到計算精度要求。大多數(shù)學者采用的y+范圍為60 以下，此范圍內得到的結果與實驗值較為吻合［18-19］。為提高翼型網(wǎng)格質量，對翼型前緣采用C 型網(wǎng)格劃分，并對翼型前緣和尾緣部分進行網(wǎng)格加密，如圖4 所示。調整后的翼型上、下表面y+分布如圖5 所示，其中Lc為翼型長度百分比。采用3 組不同加密程度的網(wǎng)格，在定常無空化工況下，以25℃水為介質，進行繞NACA0015 翼型數(shù)值計算。運用無量綱壓力系數(shù)Cp與實驗結果［17］對比，公式為

圖3 NACA0015 翼型網(wǎng)格Fig.3 Mesh of NACA0015 hydrofoil

圖4 前緣和尾緣網(wǎng)格分布Fig.4 Mesh distribution of leading edge and trailing edge

圖5 NACA0015 翼型y +分布Fig.5 Distribution of NACA0015 hydrofoil

計算結果如圖6 所示，3 組網(wǎng)格節(jié)點數(shù)分別為114 萬(方案1)、285 萬(方案2)和456 萬(方案3)。計算結果表明3 組方案的模擬壓力系數(shù)與實驗值吻合度較高，結合網(wǎng)格實用性與計算成本，選取方案2作為本文的研究方案。

圖6 網(wǎng)格無關性驗證Fig.6 Verification of mesh independence

2 結果與討論

2.1 25℃時湍流模型

采用3 種湍流模型以及修正后的湍流模型，進行常溫下繞NACA0015 翼型的非空化單相數(shù)值模擬，得到翼型吸力面的壓力系數(shù)分布，如圖7(圖中FBDCM 表示修正模型)所示。結果表明，采用修正前后的湍流模型得到的計算結果相差不大，均與實驗值比較接近，驗證了網(wǎng)格和湍流模型的適用性。

圖7 25℃時非空化工況下NACA0015 翼型表面壓力系數(shù)分布Fig.7 Pressure coefficient distribution of NACA0015 hydrofoil on non-cavitation at 25℃

采用3 種湍流模型以及修正后的湍流模型，在常溫下，繞NACA0015 翼型進行空化氣液兩相數(shù)值模擬，圖8 給出了翼型吸力面的壓力系數(shù)分布。結果表明，修正前3 種湍流模型的計算結果與實驗值在翼型長度前30%的低壓區(qū)存在一定誤差，壓力開始恢復的起點均早于實驗結果，在空化核心區(qū)域，模擬壓力系數(shù)與實驗壓力系數(shù)存在誤差。在翼型長度后70%段，RNG k-ε 與SST k-ω 湍流模型的計算結果更接近實驗值，k-ε 模型的模擬結果與實驗結果相差較大;修正后的RNG k-ε 模型與修正后的SST k-ω 模型的計算結果得到明顯的改善，壓力恢復起始點與實驗值較為接近，而修正后的k-ε 模型模擬結果與修正前相比誤差明顯減小，修正效果較為明顯。這是因為兩方程湍流模型模擬空化流動時，計算的氣液混合區(qū)的湍流粘度過大，造成過高的能量耗散，影響空化模擬的準確性，而FBDCM 模型通過分域函數(shù)結合FBM 模型和DCM 模型，在不同的區(qū)域采用不同的方式修正湍流粘度，提高了空化模擬的準確性。

圖8 25℃時空化工況下NACA0015 翼型表面壓力系數(shù)分布Fig.8 Pressure coefficient distribution of NACA0015 hydrofoil on cavitation at 25℃

圖9 給出了常溫下3 種湍流模型修正前后模擬得到的翼型表面蒸汽體積分數(shù)分布?？梢钥闯觯琸-ε模型在空穴尾部的閉合區(qū)域出現(xiàn)較大尺度渦，空化核心區(qū)域發(fā)展過程較長，而RNG k-ε 與SST k-ω 模型在空化核心區(qū)域發(fā)展過程較短。3 種湍流模型在修正前與實驗結果均存在一定誤差，其中k-ε 模型的模擬結果誤差最大。修正后的k-ε 模型消除了kε 模型空穴尾部閉合區(qū)域的渦，但空化核心區(qū)域長度較修正前明顯減小。修正的RNG k-ε 模型與修正的SST k-ω 模型計算得到的空化核心區(qū)域較修正前明顯擴大，與實驗結果更接近，但空化核心區(qū)域蒸汽體積分數(shù)增大，空化程度更嚴重。原因為k-ε 模型中引入的耗散率ε 方程是由經(jīng)驗方法得到的，不是由推導所得，并且采用的3 種湍流模型中只有k-ε模型沒有考慮應變率的影響，不能反映時均應變率，導致在空穴尾部區(qū)域具有較大的壓力梯度與大曲率流動的情況下，k-ε 模型的模擬精度較差。采用FBDCM 模型修正各湍流模型后，降低了空穴區(qū)域的湍流粘度，減少了能量耗散，空化流場得到發(fā)展。

2.2 50℃時湍流模型

根據(jù)上述討論，修正的RNG k-ε 模型與修正的SST k-ω 模型具有較好的適用性，故本節(jié)主要采用修正后的RNG k-ε 模型與修正后的SST k-ω 模型對NACA0015 翼型以50℃水為介質進行空化數(shù)值模擬分析。

圖10 給出了修正的RNG k-ε 模型與修正的SST k-ω 模型模擬所得的NACA0015 翼型吸力面的壓力分布。結果表明，在空化核心區(qū)域，兩種模型計算得到的壓力均低于實驗值，空化現(xiàn)象較實驗結果更為嚴重。修正的SST k-ω 模型的計算結果與實驗結果較符合，但壓力恢復階段的梯度與實驗結果相比較大。

圖9 25℃時NACA0015 翼型表面蒸汽體積分數(shù)分布Fig.9 Vapor volume fraction distributions on surface of NACA0015 hydrofoil

圖10 50℃時翼型表面壓力分布Fig.10 Pressure distribution of hydrofoil at 50℃

圖11 為50℃時修正的RNG k-ε 模型與修正的SST k-ω 模型得到的翼型表面相間質量傳輸速率分布，其中蒸發(fā)過程為正值，凝結過程為負值。結果表明，修正的SST k-ω 凝結速率明顯大于修正的RNG k-ε 凝結速率，但修正的SST k-ω 凝結區(qū)域明顯小于修正的RNG k-ε 凝結區(qū)域。此外，與修正的SST k-ω模型相比，修正的RNG k-ε 模型凝結起始點位置更靠近翼型頭部。在50℃時，兩種模型的模擬結果與實驗結果均有一定誤差。

2.3 70℃時湍流模型

圖12 為70℃時修正的RNG k-ε 模型與修正的SST k-ω 模型模擬得到的NACA0015 翼型吸力面壓力系數(shù)分布圖。結果表明，兩種湍流模型計算得到的空化核心區(qū)域范圍均比實驗結果大，但修正的RNG k-ε 模型模擬得到的空泡生成和潰滅趨勢與實驗結果較為一致，修正的SST k-ω 模型空化發(fā)展過程較長，在壓力恢復階段梯度較大。

圖11 NACA0015 翼型表面相間質量傳輸速率分布Fig.11 Bulk interphase mass transfer rate distributions on surface of NACA0015 hydrofoil

圖13 給出了70℃時修正的RNG k-ε 模型與修正的SST k-ω 模型模擬得到的NACA0015 翼型表面湍動能分布情況。結果表明，修正的SST k-ω 模型計算所得湍動能明顯大于修正的RNG k-ε 模型計算所得湍動能，近壁處流動較為紊亂，并且空化更加嚴重，與實驗所得結果相差較大。

2.4 修正的RNG k-ε 模型

為了分析不同水溫下的翼型空化特性，定義了無量綱空化數(shù)

圖12 70℃時翼型表面壓力系數(shù)分布Fig.12 Pressure coefficient distribution of hydrofoil at 70℃

式中 pin——進口壓力

圖14 為在相同空化數(shù)(σ=1.5)，不同水溫下，由修正的RNG k-ε 模型計算得到的翼型表面蒸汽體積分數(shù)分布圖。結果表明:隨著溫度的升高，蒸汽體積分數(shù)減小，空化強度下降，空穴面積減小，空穴尾端閉合區(qū)域的逆壓梯度減小，氣液界面變得模糊［20-23］。由于溫度的升高，水蒸氣含量降低，氣液混合區(qū)密度增大，氣液界面密度梯度增大，導致氣液界面模糊［24-26］，說明隨著溫度的升高，水的熱力學效應抑制空化的效果更加明顯，空化強度相對減弱，這與實驗觀察到的規(guī)律相一致。

圖13 NACA0015 翼型表面湍動能分布Fig.13 Turbulence kinetic energy distribution on surface of NACA0015 hydrofoil

圖14 NACA0015 翼型表面蒸汽體積分數(shù)分布Fig.14 Vapor volume fraction distributions on surface of NACA0015 hydrofoil

3 結論

(1)3 種方案的計算結果均與實驗結果吻合較好，同時驗證了網(wǎng)格的無關性。在常溫非空化條件下，修正前后3 種湍流模型的計算結果與實驗結果一致。

(2)25℃時，修正的k-ε 模型對湍流尺度具有明顯的修正效果，空泡尾部閉合區(qū)域的渦被消除，修正的RNG k-ε 模型與修正的SST k-ω 模型的模擬結果在空化區(qū)域與實驗結果較為接近。50℃時，在蒸發(fā)過程中由修正的RNG k-ε 模型計算得到的蒸發(fā)區(qū)域較小，而在凝結過程中由修正的SST k-ω 模型計算得到的凝結速率較大。70℃時，修正的SST k-ω模型在近壁處的湍動能明顯大于修正的RNG k-ε模型的結果，空化程度更嚴重，與實驗結果相差較大。

(3)在不同溫度、相同空化數(shù)下，修正的RNG k-ε模型揭示了NACA0015 翼型空化情況隨溫度的變化規(guī)律，并且與實驗觀察結果相一致，驗證了修正的RNG k-ε 模型準確性。

(4)對3 種溫度下湍流模型修正前后的結果進行了綜合分析，并與實驗結果進行了比較。結果表明，修正的RNG k-ε 模型對于熱敏介質的氣液兩相流模擬具有較好的適用性。