考慮熱力學效應的空化模型修正及其適用性研究

張德勝 萬福來 許 彬 王超超 施衛東

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

0 引言

空化是一種包含湍流、多相流、相變、噪聲、可壓縮和非定常特性的復雜流動現象[1-3]。近年來，低溫介質的運用促進了低溫空化的研究。液化天然氣等低溫介質已經成為我國能源領域迅猛發展的新型產業[4-5]，而其溫度為-160℃左右，極容易導致泵內發生空化現象，這對低溫介質的空化性能提出了較高的研究要求。由于大部分空化都發生在常溫水中，所以通常假設空化過程為等溫過程，而低溫流體介質物理屬性對溫度變化敏感，故熱力學效應決定低溫流體空化流數值模擬的準確性顯著[6-7]。

在低溫流體空化的數值計算和實驗研究方面，STEPANOFF等[8]最早提出了一種B因子理論，用于預測泵中熱力學效應對揚程的影響；HORD[9]以液氮和液氫為介質，做了較為全面系統的低溫流體的空化實驗，進一步加深了對低溫空化特性的認識。CERVONE等[10]對不同溫度水為介質的NACA0015翼形水翼空化開展了一系列實驗研究，分析了熱力學效應對其空化特性的影響。王巍等[11]對不同溫度水為介質的NACA0066(MOD)翼型開展了翼型表面空化場的流動分析。由于低溫流體工作環境的限制，極大增加了實驗的難度。近年來，數值計算成為研究低溫流體空化的重要手段，HOSANGADI等[12]采用Merkle空化模型對液氫和液氮繞翼型空化流動進行了模擬計算，但模擬結果與實驗對比誤差較大。為了提高數值模擬預測精度，TANI等[13]在B因子理論基礎上提出了一種新的空化模型，通過與實驗數據相比較確定了空化模型的蒸發凝結系數，使得計算結果較為吻合。國內學者也對低溫流體空化做了一些研究[14-16]，評價了不同空化模型對低溫空化的數值計算情況，并改進了個別空化模型，使得數值計算結果與實驗數據更為接近，但這些研究都是二維翼型的數值模擬計算。因此，本文重點開展三維翼型考慮熱力學效應的空化模型改進和驗證研究。

本文通過CFX前處理編輯，將3種典型空化模型嵌入CFX軟件中，并將水、液氮熱傳導系數以及飽和蒸汽壓強等隨溫度變化的物性參數引入到求解代碼程序。 同時考慮空化過程汽化潛熱的影響，并在B因子理論的基礎上，提出一種空化模型的修正方法。采用原始的和修正的空化模型，計算水及液氮繞翼型的空化流動，并與實驗結果進行對比，評估不同空化模型對空化流動的熱力學效應的敏感性和修正的空化模型的適用性，為開展低溫空化研究提供理論支撐。

1 控制方程與數值方法

1.1 基本方程

對于考慮熱力學效應的空化流動問題數值計算的控制方程，除了連續性方程和動量方程外，還包括包含能量源項的能量方程，依次為

(1)

(2)

(3)

式中ρm——混合相密度p——壓強

μ——混合介質的動力黏性系數

keff——熱傳導系數

L——汽化潛熱

Cp——定壓比熱容

fv——氣相的質量分數

T——溫度t——時間

xi——x在i方向分量

xj——x在j方向分量

δij——克羅內克函數

ui、uj——速度在i、j方向分量

μt——湍流粘度

1.2 湍流模型

計算采用了由MENTER[17]提出的SSTk-ω兩方程湍流模型，公式為

(4)

(5)

渦黏系數計算公式為

(6)

式中ρ——密度k——湍動能

ω——湍流頻率

σk3、σω2、σω3、α3、β3、β′——湍流模型常數

F1、F2——混合函數

S——剪切力張量的常數項

Pkb——浮力湍動能生成項

Pk——粘性力湍動能生成項

Pωb——湍流模型自定義項

a1——湍流模型常數

湍流模型常數取值參照文獻[18]。

1.3 空化模型

(7)

式中u——速度矢量

αv——氣相體積分數

ρv——氣相密度

1.3.1Zwart空化模型

Zwart空化模型[19]是基于單個空泡的生成和發展時空泡體積的變化，基于Rayleigh-Plesse空泡生長方程推導出的蒸發源項和凝結源項表達式

(8)

(9)

式中RB——平均空泡半徑，取1×10-6m

ρl——液相密度pv——氣相壓力

αnuc——氣核體積分數，取5×10-4

Fvap、Fcond——蒸發、凝結系數，取50、0.01

1.3.2Merkle空化模型

Merkle空化模型[20]是基于汽液兩相流模型，由混合密度導出相間質量傳輸率，公式為

(10)

(11)

式中U∞——參考速度

t∞——參考時間

經驗系數Cvap、Ccond分別取1、80。

1.3.3Singhal空化模型

Singhal空化模型[21]是基于汽液兩相流模型，綜合考慮了空泡在相變過程中表面張力和非凝結性氣體的影響，公式為

(12)

(13)

式中fg——非凝結氣體的質量分數，取10-8

σ——液體表面張力系數

經驗系數Cvap、Ccond分別取0.02、0.01。

1.4 空化模型的修正方法

在上述空化模型中，通常認為質量傳輸只由液相和氣相的壓力差來驅動，但在高溫水和低溫液體中，發生空化時，液體汽化吸收汽化潛熱，導致空泡附近液體溫度降低，從而形成空泡外的薄液體邊界層，假設泡內的氣體溫度是均勻的，其壓力等于飽和壓力，溫度等于飽和壓力對應下的飽和溫度[13]，由于熱力學邊界層的存在，泡內和泡外形成了一溫度差ΔT，則相變的熱平衡公式為

ρvυvL=ρlυlCplΔT

(14)

式中υv、υl——氣相、液相的體積質量流量

Cpl——液相的定壓比熱容

B因子方法[22]及理想氣液混合時B因子的表達形式為[23]

(15)

式中B——無量綱因子

國內外化工行業常用的純液體飽和蒸汽壓的三參數方程Antoine方程[24]為

lgpsat=a-b/(T+c)

(16)

式中a、b、c——Antoine常數，取值參照文獻[23]

psat——液相飽和蒸汽壓力

假定空泡周圍熱力學狀態保持為飽和狀態，則可得由于潛熱傳遞而形成的局部飽和蒸汽壓的變化ΔpL，其公式為

(17)

除此之外，研究表明湍動能對空化產生重要的影響[24],采用文獻[24]中所提出的方法來計算湍動能k對當地汽化壓強的影響，其公式為

Δpturb=0.195ρmk

(18)

式中 Δpturb——湍動能引起的壓力變化量

基于以上分析，可對1.3節中3個空化模型中蒸發源項和凝結源項就考慮熱力學效應的影響進行修正,從而可得

pv=psat+ΔpL+Δpturb

(19)

2 不同溫度水的空化流數值模擬

2.1 計算區域及邊界條件

為了評估各空化模型的熱力學修正效果，首先對以不同溫度水為介質的NACA0015翼形的空化流場進行了模擬計算，并與CERVONE等[10]的實驗值進行比較。翼型尺寸與計算域選取與實驗條件保持一致，如圖1所示。設置邊界條件：進口給定速度與溫度，出口為平均靜壓；壁面為無滑移絕緣壁面條件。為更準確地計算空化流場，采用結構網格，并在翼型近壁面區域進行網格加密，使分子到最近壁面的無量綱距離y+為10～50之間，大多數學者認為y+值不超過60[25]，則計算y+值滿足壁面函數要求，如圖2所示。為了對計算所用網格進行驗證，本文選取了3套網格，分別模擬計算了無空化流場，如圖3a所示。所得壓力系數pc與實驗值[10]進行對比驗證，其中Lc為翼型長度百分數。可知計算所得壓力系數與實驗較為接近,如圖3b所示，并考慮計算經濟性，選取方案2進行后續的計算。不同方案的網格信息如表1所示，計算工況見表2。

圖1 三維水體模型及邊界條件Fig.1 Schematic of 3D model with water and boundary conditions of cases

圖2 NACA0015翼型網格分布Fig.2 Grids distribution near NACA0015 hydrofoil

圖3 NACA0015網格無關性驗證Fig.3 Verification of NACA0015 mesh independence

計算中采用的無量綱參數為無窮遠處空化數σ∞、當地空化數σc、壓力系數pc、當地空化數與無窮遠處空化數之差σ，分別定義為

(20)

(21)

(22)

σ=σc-σ∞

(23)

式中pout——無窮遠處壓力

T∞——無窮遠處溫度

uin——進口速度Tc——當地溫度

pv(Tc) ——溫度為Tc時的飽和蒸氣壓力

pv(T∞) ——溫度為T∞時的飽和蒸氣壓力

表1 網格信息Tab.1 Information of different grids

表2 不同溫度水空化數值模擬邊界條件Tab.2 Boundary conditions of numerical simulation about cavitation in water at different temperatures

2.2 計算結果與討論

為了分析修正空化模型在不同溫度水中空化數值模擬中的應用情況，分別運用1.4節中考慮熱力學效應修正的Zwart、Merkle及Singhal 3種空化模型，并對298、323、343 K溫度水繞NACA0015翼型空化進行了數值模擬研究，并與實驗[8]所得翼型吸力面的壓力系數分布進行對比分析，如圖4所示。結果表明，在3種溫度下，修正前后Singhal模型及Zwart模型模擬所得翼型吸力面壓力系數分布無明顯差別。溫度為298 K和323 K時，水的熱力學效應比較穩定，壓力隨溫度變化梯度較小，3種空化模型修正前后的模擬結果差異較小。而修正的Merkle模型在水溫為343 K所模擬的結果表現出較強的修正效果，修正前后有較明顯差別。修正Merkle模型數值模擬得到的翼型吸力面低壓區范圍較之未修正前的模擬結果均有所增加，對應壓力有所降低，并隨溫度的增加，降低得越多，體現出不同溫度水中空化吸收潛熱，引起飽和蒸汽壓強不同程度降低的熱力學效應作用，且越靠近臨界溫度，其敏感度越高。而同時對應空穴閉合區壓力變化梯度均明顯增大，且隨溫度的增加，其壓力變化梯度也隨之增大，從而更加接近于實驗值。

圖4 不同溫度下不同修正空化模型模擬得到的NACA0015翼型吸力面壓力系數分布Fig.4 Pressure coefficient distribution on suction surface of NACA0015 hydrofoil

3 不同溫度液氮的空化流數值模擬

3.1 計算區域及邊界條件

HORD[9]針對不同溫度液氮繞水翼空化進行了較為全面系統的低溫流體空化實驗研究，其實驗段及水翼結構如圖5和圖6所示。

通過三維造型生成與實驗所對應三維水翼，并根據實驗中實驗段的尺寸進行了實驗流道的建模。

圖5 實驗段結構Fig.5 Test section structure

圖6 HORD水翼結構圖Fig.6 Diagram of hydrofoil structure

同時由于其對稱性結構，在數值模擬過程中可進行對稱邊界條件的設置，故只構建了一半的對稱水體模型。三維對稱水體模型結構如圖7所示。考慮與實驗的一致性，進口設置為壓力進口，其靜壓與實驗對應進口壓力相一致，來流溫度根據實驗進口所測的溫度來設置；出口設置為速度出口，其速度與實驗保持一致，邊界條件如表3所示。為精準地模擬計算空化流場，在ICEM 軟件中進行六面體結構網格的劃分，考慮到水翼前緣的圓頭形狀的影響，采用C型結構化網格，在水翼近壁面區域進行網格加密如圖8所示。選取3套網格計算液氮溫度83.06 K時翼型表面壓力系數分布情況，如圖9a所示,各方案網格信息如表4所示。模擬所得壓力系數與實驗值[9]進行比較，如圖9b所示。考慮數值計算的準確性與經濟性，選擇方案2作為后續模擬計算網格。

圖7 三維對稱水體模型及邊界條件Fig.7 Schematic of symmetrical 3D model with water and boundary conditions of cases

3.2 計算結果與討論

為了分析不同空化模型在低溫流體空化數值模擬中的應用情況，分別運用Zwart、Merkle及Singhal空化模型，采用等溫isothermal的流場熱傳輸模式，分別對83.06、77.64 K兩種不同溫度液氮繞水翼的空化進行數值計算模擬研究，并與HORD的實驗數據[9]進行對比分析，所得水翼表面的壓力分布如圖10所示。在兩種液氮中，3種模型模擬所得低壓區壓強保持為對應遠場溫度下的飽和蒸汽壓強不變，且均大于實驗值，說明液氮中熱力學效應明顯，不容忽視，其中空化與常溫水中不同，不能假設為等溫過程。而模擬所得低壓區范圍差異較大，Merkle模型模擬結果較為接近實驗值，而Zwart模型與Singhal模型兩者模擬結果相近均偏離實驗值較多，過多地預測了低壓區即空化區域的范圍。

表3 不同液氮中空化數值模擬邊界條件Tab.3 Boundary conditions of numerical simulation about cavitation in nitrogen at different temperatures

圖8 HORD翼型網格分布Fig.8 Grids distribution near HORD hydrofoil

圖9 HORD翼型網格變化Fig.9 Changes of HORD hydrofoil mesh

方案網格節點數最小角度最大寬高比最小質量1325409653.462920.8042388265653.462940.8043486105653.462920.804

圖10 空化模型模擬所得HORD水翼表面的壓力分布Fig.10 Pressure distribution on surface of HORD hydrofoil with different cavitation models

綜合上述分析可知，Merkle模型在低溫流體空化數值模擬中有較好的可行性，其與實驗吻合度較為接近。為了分析修正空化模型在不同溫度液氮空化數值模擬中的應用情況，運用考慮熱力學效應修正的Merkle模型，分別對83.06、77.64 K兩種不同溫度液氮繞水翼的空化進行數值計算模擬，并與實驗數據進行對比分析，所得水翼表面的壓力分布如圖11所示。結果表明，修正前后Merkle模型所模擬的結果有較明顯差別，修正的Merkle模型模擬所得水翼表面低壓區范圍較之未修正前的模擬結果均有所減小，低壓區壓力有所降低，且83.06 K中降低量較77.64 K下更多，使得兩者模擬結果更加吻合實驗低壓區壓力，這體現出低溫流體中空化吸收潛熱，引起飽和蒸汽壓強不同程度降低的熱力學效應作用，也進一步說明了越靠近其臨界溫度，其敏感度越高。

圖11 修正空化模型數值模擬所得HORD水翼表面的壓力分布Fig.11 Pressure distribution on surface of HORD hydrofoil with modified cavitation model

采用考慮熱力學效應修正及等溫模式的Merkle模型模擬所得蒸汽體積分數分布情況，如圖12所示，圖中y軸表示翼型絕對尺寸，x軸表示翼型相對尺寸。結果表明，等溫模式下，無熱力學效應作用，空穴區域內的壓強恒定為對應遠場溫度下的飽和蒸汽壓強，并且83.06 K中對應遠場空化數為1.706，而77.64 K中為1.765，故在83.06 K中的空化區域大于77.64 K的，其空穴長度明顯較長。而考慮熱力學效應時，局部的溫降，導致飽和蒸汽壓強的降低，從而降低了空化強度，使得空化區域的蒸汽體積分數減小，空穴長度明顯縮短。

根據當地空化數σc的定義，為了更好地表征熱力學效應在低溫流體空化中的影響，計算了兩種液氮中，修正Merkle模型模擬所得空化流場的當地空化數與無窮遠處空化數之差σ分布情況，如圖13所示。考慮熱力學效應時，當地空化數最大的區域與最大壓降區、最大溫降區相對應，都在空化核心區域，其當地空化數均大于相應的遠場空化數，由此與等溫模式下相比，考慮熱力學效應時的空化強度降低。

圖12 不同液氮中模擬所得HORD水翼表面蒸汽體積分數云圖Fig.12 Vapour volume fraction on surface of HORD hydrofoil simulated in different liquid nitrogen

圖13 不同液氮中模擬所得HOAD水翼表面σ分布云圖Fig.13 Distribution nephogram of σ on surface of HORD hydrofoil simulated in different liquid nitrogen

4 結論

(1)采用Zwart、Merkle及Singhal共3種空化模型及對應修正后的空化模型,分別對3種溫度水進行數值模擬，并與實驗進行對比，結果表明修正的Merkle模型能更好地體現熱力學效應對空化的影響。

(2)采用3種空化模型分別數值模擬83.06、77.64 K兩種液氮繞水翼空化流場，并與實驗數據進行了對比，發現Merkle模型數值模擬結果與實驗值較為接近，適用性相對較強。考慮熱力學效應修正的Merkle空化模型，對兩種不同溫度液氮中空化進行了數值模擬研究，并與實驗所得水翼表面壓力分布數據進行對比分析，再次驗證了修正空化模型的適用性及空化模型修正方法在低溫流體空化數值模擬中的可行性。

(3)在修正Merkle模型數值模擬結果的基礎上，與等溫模式的數值模擬結果相對比，其空化區域由于溫降的產生，使得對應飽和蒸汽壓強有所降低，致使空化強度減弱，空穴長度縮短，當地空化數均大于遠場空化數。且相比77.64 K，83.06 K的溫度較高，對于熱力學效應的敏感性更高，其空化區域較多的壓降使其熱力學效應更顯著，空化受到抑制作用更明顯，空穴縮短更多，空化區域蒸汽體積分數降低更多。