周書航, 邊浩志, 吳桐雨, 李文濤, 丁銘
(1.哈爾濱工程大學 核動力裝置性能與設備黑龍江省重點實驗室,黑龍江 哈爾濱 150001; 2.哈爾濱工程大學 核科學與技術學院,黑龍江 哈爾濱 150001)
當反應堆發生破口事故時,會向安全殼內噴放大量的高溫高壓蒸汽,致使安全殼內壓力和溫度急劇上升乃至威脅其完整性[1-2]。為了確保安全殼邊界完整性,非能動安全殼熱量導出系統普遍應用在第三代核電廠中,其主要利用流體密度差產生驅動力,通過蒸汽冷凝的方式持續帶走安全殼內熱量[3-5]。然而在蒸汽冷凝過程中,安全殼內不可避免地會存在不凝性氣體[6-9]。
文獻[10]研究表明,不凝性氣體的存在會極大削弱蒸汽冷凝能力,即使只存在1%的空氣,冷凝傳熱系數也會下降60%。隨后國內外學者對含不凝性氣體蒸汽冷凝現象開展了大量的實驗研究[11-13]以及數值分析[14-17],其中,數值分析中普遍采用的模型是擴散邊界層模型。該模型可以在較廣的參數范圍預測蒸汽冷凝局部流動和傳熱現象,但在高壓、高蒸汽濃度下的適用性仍然存在問題。
針對這一問題,Bird[18]、Dehbi[19]、Bian[20]等分別提出了修正因子進一步提高了擴散邊界層模型的準確性和適用性。但目前國內外在對含不凝性氣體蒸汽冷凝過程進行數值分析時,由于不凝性氣體的存在,幾乎忽略了液膜熱阻。然而,在高壓、高蒸汽濃度下液膜熱阻是否可忽略,以及液膜將如何影響蒸汽冷凝特性尚未揭示。因此,本文采用CFD 計算軟件STAR-CCM+(17.02)開展含空氣蒸汽冷凝和層流液膜耦合模型研究,評估液膜效應如何影響蒸汽冷凝特性以及在不同參數范圍內的變化機制。
含空氣蒸汽冷凝-層流液膜耦合模型的建立需要針對多組分氣體側、液膜側以及氣-液交界面分別建立控制方程描述其冷凝的完整過程。
氣體側控制方程主要描述多組分氣體的流動與傳熱過程,其中控制方程主要包括質量守恒方程、動量守恒方程、能量守恒方程和組分輸運方程。
質量守恒方程:
(1)
動量守恒方程:
ρgfg+Sg,pw
(2)
能量守恒方程:
(3)
組分輸運方程:
(4)
式中:ρ為流體密度,kg/m3;t為時間,s;w為速度矢量,m/s;Sm為質量源項,kg/(m3·s);P表示表面力,N/m2;τ為剪切力,m2/s;f為體積力,N/m3;Spw為動量源項,N/m3;E為流動流體所具有的能量,J;ω為質量分數;keff為流體的等效傳熱系數,W/(m·K);Sh為能量源項,J/(m3·s);D為擴散系數,m2/s;下標i表示組分,下標g表示氣體。
此外,氣體側流速較高時,其流態會由層流過渡到湍流態,需在此基礎上添加湍流模型。本文采用可實現的k-ε兩層湍流模型。該模型適用于各種流動的過程,相比于標準k-ε模型有更好的預測精度,具有適用性廣、精度高的優點。
層流液膜控制方程主要描述液膜的層流流動與傳熱過程,其中控制方程主要包括質量守恒方程、動量守恒方程、能量守恒方程和組分輸運方程。
質量守恒方程:
(5)
動量守恒方程:
(6)
能量守恒方程:
(7)
組分輸運方程:
(8)
式中:δ是液膜厚度,m;下標l表示液體。
氣-液交界面控制方程主要描述蒸汽冷凝傳熱過程,其冷凝傳熱過程采用基于擴散邊界層理論的冷凝模型,該模型可以很好地預測局部流動和傳熱現象。當蒸汽在氣-液交界面冷凝時,氣相和液相之間將發生質量、動量和能量交換。
質量源項:
(9)
式中Δ為近壁面第1層網格厚度,m。其中,混合氣體的擴散系數D計算公式為[21]:
(10)
式中:下標0表示標準狀態。
動量源項:
Spw=Sm·w
(11)
能量源項:
Sh=Smhv
(12)
為評估基于擴散邊界層建立的含空氣蒸汽冷凝-層流液膜耦合模型,需開展模型驗證工作。COAST實驗[6]的參數范圍(壓力:0.2~1.6 MPa,空氣質量分數:0.16~0.71,過冷度:45~117 ℃)較廣,相比于安全殼事故工況,其具有一定的代表性。為在較廣的參數范圍內評估模型,本文針對COAST實驗工況進行了選取并建立了數值模型。該模型為直徑1.5 m、高度2.5 m的豎直殼體,殼體中心設置有長度1 m,直徑19 mm的單管,如圖1所示。冷凝段設置為恒壁溫壁面,入口采用質量流量入口。

圖1 幾何與網格模型示意Fig.1 Geometry and mesh model
為確保該模型網格條件的計算精度,需開展網格無關性驗證工作,由于本文所采用幾何模型與作者前序研究所采用的模型相同,因此網格條件的選取(主流網格尺寸m= 0.04 m 和Y+= 1,網格數量為30萬)與前序研究保持一致[20]。
通過對COAST實驗工況數值模擬發現,如圖2所示,隨著壓力、蒸汽質量分數的提高,擴散邊界層模型預測的蒸汽冷凝傳熱系數與實驗值偏差逐漸變大,尤其在高壓、高蒸汽質量分數條件下,偏差可達到50%以上。主要原因考慮是擴散邊界層模型無法預測冷凝率增強引發的抽吸效應以及液膜流動對氣側影響的液膜效應,因此需對模型進行修正。

圖2 擴散邊界層模型驗證Fig.2 Validation of diffusion boundary layer model
由于擴散邊界層模型本身的缺陷,Bird認為當蒸汽在換熱面發生冷凝后,換熱面附近和主流區域之間會產生較大的氣體濃度梯度,從而增強氣體之間的對流擴散以強化換熱,這種現象稱為抽吸效應[19-20]。
為了在模型中考慮這一額外對流/擴散傳質項,本文通過自定義場函數的形式在冷凝模型中添加了Bird修正因子B和θB,并將擴散系數DBird添加到模型中:
(13)
(14)
(15)
DBird=D0·θB
(16)
此外,Dehbi認為Bird修正因子會高估其冷凝率,因此基于Bird修正因子提出了新的修正因子θD,并與Dehbi實驗進行了驗證,關系式為:
θD=0.346 64+0.495 24θB+0.165 02θB2
(17)
本文在模型中分別添加Bird和Dehbi修正因子與實驗值進行對比,如圖3所示。在模型中添加Bird和Dehbi修正因子后,其預測的傳熱系數偏差雖相較于擴散邊界層模型來說已大幅減小,但整體偏差大于20%。

圖3 各種修正因子比較Fig.3 Comparisons among various correction factors
其主要原因考慮是Bird和Dehbi修正因子均不是在自然對流條件下提出,且都是基于忽略液膜的假設,而為使模型適用于自然對流條件下分析且考慮液膜效應,需對其進一步修正。
由于該模型考慮了液膜的影響,當處在高壓、高蒸汽濃度工況下時,蒸汽冷凝率會顯著提升,液膜厚度也隨之增加,導致更大的液膜熱阻以及液膜軸向流速,如圖4所示。

圖4 液膜厚度及流速變化規律Fig.4 Variation of liquid film thickness and velocity
然而,液膜熱阻增大會導致換熱系數降低,液膜軸向流速增大會提升冷凝換熱系數,因此液膜效應對冷凝換熱既有強化作用,也有抑制作用。為更好地預測自然對流條件下的液膜效應,本文基于COAST實驗提出新的修正因子θZ:
(18)
為評估新修正因子θZ,本文基于COAST實驗、Dehbi實驗工況[22]進行了計算對比驗證,如圖5、6所示。結果表明:與Bird和Dehbi修正因子相比,新修正因子θZ與實驗結果吻合較好,偏差普遍在10%以內。

圖5 修正后模型驗證Fig.5 Validation of modified model

圖6 耦合模型預測值與實驗值對比Fig.6 Comparison between predicted results of the coupled model and experimental results
由于目前CFD研究中通常在含不凝性氣體蒸汽冷凝條件下忽略液膜熱阻,本文為評估忽略液膜產生的偏差以及液膜效應的作用機制,在P=1.6 MPa,ΔT=20 ℃工況,對不同空氣質量分數下忽略液膜的單管外含空氣蒸汽冷凝進行了數值計算,如圖7所示。結果表明:忽略液膜后,當空氣質量分數高于0.5時,偏差小于10%,而當空氣質量分數小于0.5,偏差將迅速上升,在空氣質量分數ωa=0.1時偏差達到20%以上。由圖8也可以看出,隨著空氣質量分數的減小,總熱阻和氣側熱阻呈現下降的趨勢,而液膜熱阻呈現上升的趨勢,并且在空氣質量分數小于0.5時,液膜的熱阻占比會快速增加,最高可達60%以上?;谏鲜鼋Y果得出,當處在高壓、高蒸汽濃度工況時,液膜熱阻的影響是不可忽略的,有必要在模型中考慮液膜效應的影響。

圖7 忽略液膜偏差分布Fig.7 Deviation distribution under ignoring liquid film

圖8 不同空氣質量分數下熱阻分布Fig.8 Thermal resistance distribution under different air concentration
此外,本文評估了液膜效應在不同工況參數下的作用機制,如圖9所示。結果表明:忽略液膜模型將明顯低估實際傳熱過程中的換熱面氣膜軸向流動速度,并且隨著壓力、蒸汽質量分數的增大,其模型低估導致的偏差將越來越大。基于以上分析,液膜對換熱既有強化作用也有抑制作用,強化作用主要是液膜會帶動近壁面氣膜起到軸向加速的作用,從而增強冷凝換熱。抑制作用主要是液膜本身存在熱阻導致。這一結果進一步證實了液膜的影響是不可忽略的,尤其是在高壓、高蒸汽濃度條件下。

圖9 有液膜與無液膜下氣體流速對比Fig.9 Comparison of gas velocity with and without liquid film
在高壓高蒸汽濃度條件下,冷凝量會大幅提升,從而使得冷凝液膜急劇變厚,除了液膜熱阻變得不可忽略外,其液膜的流態可能不限于層流狀態。為獲取實驗中液膜流態區域(層流區域、湍流區域),本文針對實驗參數范圍進行了數值計算,提取對應工況下的液膜最大雷諾數,如圖10所示。結果表明:COAST實驗中壓力P≤ 1.3 MPa的實驗工況,其冷凝液膜流態均處在層流區域,模型可較好地預測其局部流動與換熱特性。

圖10 不同工況下最大雷諾數Fig.10 Maximum Reynolds number under various conditions
當處于壓力P=1.6 MPa、空氣質量分數ωa=0.16工況,其冷凝液膜最大雷諾數達到3 000,表征其冷凝液膜在部分高度已處于湍流區域。為了獲取液膜湍流區域長度,圖11示出了壓力P=1.6 MPa、空氣質量分數ωa=0.16時不同過冷度下的局部雷諾數分布。結果表明:在P=1.6 MPa、ωa=0.16工況下,其液膜層流區域長度約為0.6~0.7 m,處在換熱管頭部;湍流區域長度約為0.3~0.4 m,處在換熱管尾部。

圖11 P = 1.6 MPa、ωa=0.16下局部雷諾數分布Fig.11 Local Reynolds number distribution at P=1.6 MPa, ωa=0.16
1)提出的含空氣蒸汽冷凝-層流液膜耦合模型,可在廣的參數范圍內較好地預測自然對流條件下的含空氣蒸汽冷凝過程,其偏差普遍在10%以內。
2)冷凝液膜對換熱既有強化作用,也有抑制作用。抑制作用是由于其本身存在熱阻,強化作用主要是液膜對近壁面的氣膜具有軸向加速作用。當空氣質量分數低于0.5時,忽略液膜會產生較大偏差,最大可達到20%以上。
3)總結得到了COAST實驗的液膜流態區域。在高壓、高蒸汽濃度條件下,管外冷凝液膜會達到湍流態。當處在壓力P≤1.3 MPa的實驗工況,其冷凝液膜流態均處在層流區域,當處于壓力P=1.6 MPa工況,其尾部液膜已處于湍流區域。