豎直管外含空氣蒸汽冷凝-層流液膜耦合分析

周書航, 邊浩志, 吳桐雨, 李文濤, 丁銘

(1.哈爾濱工程大學(xué) 核動(dòng)力裝置性能與設(shè)備黑龍江省重點(diǎn)實(shí)驗(yàn)室,黑龍江 哈爾濱 150001; 2.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院,黑龍江 哈爾濱 150001)

當(dāng)反應(yīng)堆發(fā)生破口事故時(shí),會(huì)向安全殼內(nèi)噴放大量的高溫高壓蒸汽,致使安全殼內(nèi)壓力和溫度急劇上升乃至威脅其完整性[1-2]。為了確保安全殼邊界完整性,非能動(dòng)安全殼熱量導(dǎo)出系統(tǒng)普遍應(yīng)用在第三代核電廠中,其主要利用流體密度差產(chǎn)生驅(qū)動(dòng)力,通過蒸汽冷凝的方式持續(xù)帶走安全殼內(nèi)熱量[3-5]。然而在蒸汽冷凝過程中,安全殼內(nèi)不可避免地會(huì)存在不凝性氣體[6-9]。

文獻(xiàn)[10]研究表明,不凝性氣體的存在會(huì)極大削弱蒸汽冷凝能力,即使只存在1%的空氣,冷凝傳熱系數(shù)也會(huì)下降60%。隨后國內(nèi)外學(xué)者對(duì)含不凝性氣體蒸汽冷凝現(xiàn)象開展了大量的實(shí)驗(yàn)研究[11-13]以及數(shù)值分析[14-17],其中,數(shù)值分析中普遍采用的模型是擴(kuò)散邊界層模型。該模型可以在較廣的參數(shù)范圍預(yù)測(cè)蒸汽冷凝局部流動(dòng)和傳熱現(xiàn)象,但在高壓、高蒸汽濃度下的適用性仍然存在問題。

針對(duì)這一問題,Bird[18]、Dehbi[19]、Bian[20]等分別提出了修正因子進(jìn)一步提高了擴(kuò)散邊界層模型的準(zhǔn)確性和適用性。但目前國內(nèi)外在對(duì)含不凝性氣體蒸汽冷凝過程進(jìn)行數(shù)值分析時(shí),由于不凝性氣體的存在,幾乎忽略了液膜熱阻。然而,在高壓、高蒸汽濃度下液膜熱阻是否可忽略,以及液膜將如何影響蒸汽冷凝特性尚未揭示。因此,本文采用CFD 計(jì)算軟件STAR-CCM+(17.02)開展含空氣蒸汽冷凝和層流液膜耦合模型研究,評(píng)估液膜效應(yīng)如何影響蒸汽冷凝特性以及在不同參數(shù)范圍內(nèi)的變化機(jī)制。

1 含空氣蒸汽冷凝-層流液膜耦合模型建立

含空氣蒸汽冷凝-層流液膜耦合模型的建立需要針對(duì)多組分氣體側(cè)、液膜側(cè)以及氣-液交界面分別建立控制方程描述其冷凝的完整過程。

1.1 氣體側(cè)控制方程

氣體側(cè)控制方程主要描述多組分氣體的流動(dòng)與傳熱過程,其中控制方程主要包括質(zhì)量守恒方程、動(dòng)量守恒方程、能量守恒方程和組分輸運(yùn)方程。

質(zhì)量守恒方程：

(1)

動(dòng)量守恒方程：

ρgfg+Sg,pw

(2)

能量守恒方程：

(3)

組分輸運(yùn)方程：

(4)

式中：ρ為流體密度,kg/m3;t為時(shí)間,s;w為速度矢量,m/s;Sm為質(zhì)量源項(xiàng),kg/(m3·s);P表示表面力,N/m2;τ為剪切力,m2/s;f為體積力,N/m3;Spw為動(dòng)量源項(xiàng),N/m3;E為流動(dòng)流體所具有的能量,J;ω為質(zhì)量分?jǐn)?shù);keff為流體的等效傳熱系數(shù),W/(m·K);Sh為能量源項(xiàng),J/(m3·s);D為擴(kuò)散系數(shù),m2/s;下標(biāo)i表示組分,下標(biāo)g表示氣體。

此外,氣體側(cè)流速較高時(shí),其流態(tài)會(huì)由層流過渡到湍流態(tài),需在此基礎(chǔ)上添加湍流模型。本文采用可實(shí)現(xiàn)的k-ε兩層湍流模型。該模型適用于各種流動(dòng)的過程,相比于標(biāo)準(zhǔn)k-ε模型有更好的預(yù)測(cè)精度,具有適用性廣、精度高的優(yōu)點(diǎn)。

1.2 層流液膜控制方程

層流液膜控制方程主要描述液膜的層流流動(dòng)與傳熱過程,其中控制方程主要包括質(zhì)量守恒方程、動(dòng)量守恒方程、能量守恒方程和組分輸運(yùn)方程。

質(zhì)量守恒方程：

(5)

動(dòng)量守恒方程：

(6)

能量守恒方程：

(7)

組分輸運(yùn)方程：

(8)

式中：δ是液膜厚度,m;下標(biāo)l表示液體。

1.3 氣-液交界面控制方程

氣-液交界面控制方程主要描述蒸汽冷凝傳熱過程,其冷凝傳熱過程采用基于擴(kuò)散邊界層理論的冷凝模型,該模型可以很好地預(yù)測(cè)局部流動(dòng)和傳熱現(xiàn)象。當(dāng)蒸汽在氣-液交界面冷凝時(shí),氣相和液相之間將發(fā)生質(zhì)量、動(dòng)量和能量交換。

質(zhì)量源項(xiàng)：

(9)

式中Δ為近壁面第1層網(wǎng)格厚度,m。其中,混合氣體的擴(kuò)散系數(shù)D計(jì)算公式為[21]：

(10)

式中：下標(biāo)0表示標(biāo)準(zhǔn)狀態(tài)。

動(dòng)量源項(xiàng)：

Spw=Sm·w

(11)

能量源項(xiàng)：

Sh=Smhv

(12)

1.4 模型驗(yàn)證

為評(píng)估基于擴(kuò)散邊界層建立的含空氣蒸汽冷凝-層流液膜耦合模型,需開展模型驗(yàn)證工作。COAST實(shí)驗(yàn)[6]的參數(shù)范圍(壓力：0.2～1.6 MPa,空氣質(zhì)量分?jǐn)?shù)：0.16～0.71,過冷度：45～117 ℃)較廣,相比于安全殼事故工況,其具有一定的代表性。為在較廣的參數(shù)范圍內(nèi)評(píng)估模型,本文針對(duì)COAST實(shí)驗(yàn)工況進(jìn)行了選取并建立了數(shù)值模型。該模型為直徑1.5 m、高度2.5 m的豎直殼體,殼體中心設(shè)置有長(zhǎng)度1 m,直徑19 mm的單管,如圖1所示。冷凝段設(shè)置為恒壁溫壁面,入口采用質(zhì)量流量入口。

圖1 幾何與網(wǎng)格模型示意Fig.1 Geometry and mesh model

為確保該模型網(wǎng)格條件的計(jì)算精度,需開展網(wǎng)格無關(guān)性驗(yàn)證工作,由于本文所采用幾何模型與作者前序研究所采用的模型相同,因此網(wǎng)格條件的選取(主流網(wǎng)格尺寸m= 0.04 m 和Y+= 1,網(wǎng)格數(shù)量為30萬)與前序研究保持一致[20]。

通過對(duì)COAST實(shí)驗(yàn)工況數(shù)值模擬發(fā)現(xiàn),如圖2所示,隨著壓力、蒸汽質(zhì)量分?jǐn)?shù)的提高,擴(kuò)散邊界層模型預(yù)測(cè)的蒸汽冷凝傳熱系數(shù)與實(shí)驗(yàn)值偏差逐漸變大,尤其在高壓、高蒸汽質(zhì)量分?jǐn)?shù)條件下,偏差可達(dá)到50%以上。主要原因考慮是擴(kuò)散邊界層模型無法預(yù)測(cè)冷凝率增強(qiáng)引發(fā)的抽吸效應(yīng)以及液膜流動(dòng)對(duì)氣側(cè)影響的液膜效應(yīng),因此需對(duì)模型進(jìn)行修正。

圖2 擴(kuò)散邊界層模型驗(yàn)證Fig.2 Validation of diffusion boundary layer model

2 修正的含空氣蒸汽冷凝-層流液膜耦合模型

2.1 抽吸效應(yīng)

由于擴(kuò)散邊界層模型本身的缺陷,Bird認(rèn)為當(dāng)蒸汽在換熱面發(fā)生冷凝后,換熱面附近和主流區(qū)域之間會(huì)產(chǎn)生較大的氣體濃度梯度,從而增強(qiáng)氣體之間的對(duì)流擴(kuò)散以強(qiáng)化換熱,這種現(xiàn)象稱為抽吸效應(yīng)[19-20]。

為了在模型中考慮這一額外對(duì)流/擴(kuò)散傳質(zhì)項(xiàng),本文通過自定義場(chǎng)函數(shù)的形式在冷凝模型中添加了Bird修正因子B和θB,并將擴(kuò)散系數(shù)DBird添加到模型中：

(13)

(14)

(15)

DBird=D0·θB

(16)

此外,Dehbi認(rèn)為Bird修正因子會(huì)高估其冷凝率,因此基于Bird修正因子提出了新的修正因子θD,并與Dehbi實(shí)驗(yàn)進(jìn)行了驗(yàn)證,關(guān)系式為：

θD=0.346 64+0.495 24θB+0.165 02θB2

(17)

本文在模型中分別添加Bird和Dehbi修正因子與實(shí)驗(yàn)值進(jìn)行對(duì)比,如圖3所示。在模型中添加Bird和Dehbi修正因子后,其預(yù)測(cè)的傳熱系數(shù)偏差雖相較于擴(kuò)散邊界層模型來說已大幅減小,但整體偏差大于20%。

圖3 各種修正因子比較Fig.3 Comparisons among various correction factors

其主要原因考慮是Bird和Dehbi修正因子均不是在自然對(duì)流條件下提出,且都是基于忽略液膜的假設(shè),而為使模型適用于自然對(duì)流條件下分析且考慮液膜效應(yīng),需對(duì)其進(jìn)一步修正。

2.2 液膜效應(yīng)

由于該模型考慮了液膜的影響,當(dāng)處在高壓、高蒸汽濃度工況下時(shí),蒸汽冷凝率會(huì)顯著提升,液膜厚度也隨之增加,導(dǎo)致更大的液膜熱阻以及液膜軸向流速,如圖4所示。

圖4 液膜厚度及流速變化規(guī)律Fig.4 Variation of liquid film thickness and velocity

然而,液膜熱阻增大會(huì)導(dǎo)致?lián)Q熱系數(shù)降低,液膜軸向流速增大會(huì)提升冷凝換熱系數(shù),因此液膜效應(yīng)對(duì)冷凝換熱既有強(qiáng)化作用,也有抑制作用。為更好地預(yù)測(cè)自然對(duì)流條件下的液膜效應(yīng),本文基于COAST實(shí)驗(yàn)提出新的修正因子θZ：

(18)

為評(píng)估新修正因子θZ,本文基于COAST實(shí)驗(yàn)、Dehbi實(shí)驗(yàn)工況[22]進(jìn)行了計(jì)算對(duì)比驗(yàn)證,如圖5、6所示。結(jié)果表明：與Bird和Dehbi修正因子相比,新修正因子θZ與實(shí)驗(yàn)結(jié)果吻合較好,偏差普遍在10%以內(nèi)。

圖5 修正后模型驗(yàn)證Fig.5 Validation of modified model

圖6 耦合模型預(yù)測(cè)值與實(shí)驗(yàn)值對(duì)比Fig.6 Comparison between predicted results of the coupled model and experimental results

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

3.1 液膜效應(yīng)的影響

由于目前CFD研究中通常在含不凝性氣體蒸汽冷凝條件下忽略液膜熱阻,本文為評(píng)估忽略液膜產(chǎn)生的偏差以及液膜效應(yīng)的作用機(jī)制,在P=1.6 MPa,ΔT=20 ℃工況,對(duì)不同空氣質(zhì)量分?jǐn)?shù)下忽略液膜的單管外含空氣蒸汽冷凝進(jìn)行了數(shù)值計(jì)算,如圖7所示。結(jié)果表明：忽略液膜后,當(dāng)空氣質(zhì)量分?jǐn)?shù)高于0.5時(shí),偏差小于10%,而當(dāng)空氣質(zhì)量分?jǐn)?shù)小于0.5,偏差將迅速上升,在空氣質(zhì)量分?jǐn)?shù)ωa=0.1時(shí)偏差達(dá)到20%以上。由圖8也可以看出,隨著空氣質(zhì)量分?jǐn)?shù)的減小,總熱阻和氣側(cè)熱阻呈現(xiàn)下降的趨勢(shì),而液膜熱阻呈現(xiàn)上升的趨勢(shì),并且在空氣質(zhì)量分?jǐn)?shù)小于0.5時(shí),液膜的熱阻占比會(huì)快速增加,最高可達(dá)60%以上。基于上述結(jié)果得出,當(dāng)處在高壓、高蒸汽濃度工況時(shí),液膜熱阻的影響是不可忽略的,有必要在模型中考慮液膜效應(yīng)的影響。

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

圖8 不同空氣質(zhì)量分?jǐn)?shù)下熱阻分布Fig.8 Thermal resistance distribution under different air concentration

此外,本文評(píng)估了液膜效應(yīng)在不同工況參數(shù)下的作用機(jī)制,如圖9所示。結(jié)果表明：忽略液膜模型將明顯低估實(shí)際傳熱過程中的換熱面氣膜軸向流動(dòng)速度,并且隨著壓力、蒸汽質(zhì)量分?jǐn)?shù)的增大,其模型低估導(dǎo)致的偏差將越來越大?；谝陨戏治?液膜對(duì)換熱既有強(qiáng)化作用也有抑制作用,強(qiáng)化作用主要是液膜會(huì)帶動(dòng)近壁面氣膜起到軸向加速的作用,從而增強(qiáng)冷凝換熱。抑制作用主要是液膜本身存在熱阻導(dǎo)致。這一結(jié)果進(jìn)一步證實(shí)了液膜的影響是不可忽略的,尤其是在高壓、高蒸汽濃度條件下。

圖9 有液膜與無液膜下氣體流速對(duì)比Fig.9 Comparison of gas velocity with and without liquid film

3.2 流態(tài)區(qū)域劃分

在高壓高蒸汽濃度條件下,冷凝量會(huì)大幅提升,從而使得冷凝液膜急劇變厚,除了液膜熱阻變得不可忽略外,其液膜的流態(tài)可能不限于層流狀態(tài)。為獲取實(shí)驗(yàn)中液膜流態(tài)區(qū)域(層流區(qū)域、湍流區(qū)域),本文針對(duì)實(shí)驗(yàn)參數(shù)范圍進(jìn)行了數(shù)值計(jì)算,提取對(duì)應(yīng)工況下的液膜最大雷諾數(shù),如圖10所示。結(jié)果表明：COAST實(shí)驗(yàn)中壓力P≤ 1.3 MPa的實(shí)驗(yàn)工況,其冷凝液膜流態(tài)均處在層流區(qū)域,模型可較好地預(yù)測(cè)其局部流動(dòng)與換熱特性。

圖10 不同工況下最大雷諾數(shù)Fig.10 Maximum Reynolds number under various conditions

當(dāng)處于壓力P=1.6 MPa、空氣質(zhì)量分?jǐn)?shù)ωa=0.16工況,其冷凝液膜最大雷諾數(shù)達(dá)到3 000,表征其冷凝液膜在部分高度已處于湍流區(qū)域。為了獲取液膜湍流區(qū)域長(zhǎng)度,圖11示出了壓力P=1.6 MPa、空氣質(zhì)量分?jǐn)?shù)ωa=0.16時(shí)不同過冷度下的局部雷諾數(shù)分布。結(jié)果表明：在P=1.6 MPa、ωa=0.16工況下,其液膜層流區(qū)域長(zhǎng)度約為0.6～0.7 m,處在換熱管頭部;湍流區(qū)域長(zhǎng)度約為0.3～0.4 m,處在換熱管尾部。

圖11 P = 1.6 MPa、ωa=0.16下局部雷諾數(shù)分布Fig.11 Local Reynolds number distribution at P=1.6 MPa, ωa=0.16

4 結(jié)論

1)提出的含空氣蒸汽冷凝-層流液膜耦合模型,可在廣的參數(shù)范圍內(nèi)較好地預(yù)測(cè)自然對(duì)流條件下的含空氣蒸汽冷凝過程,其偏差普遍在10%以內(nèi)。

2)冷凝液膜對(duì)換熱既有強(qiáng)化作用,也有抑制作用。抑制作用是由于其本身存在熱阻,強(qiáng)化作用主要是液膜對(duì)近壁面的氣膜具有軸向加速作用。當(dāng)空氣質(zhì)量分?jǐn)?shù)低于0.5時(shí),忽略液膜會(huì)產(chǎn)生較大偏差,最大可達(dá)到20%以上。

3)總結(jié)得到了COAST實(shí)驗(yàn)的液膜流態(tài)區(qū)域。在高壓、高蒸汽濃度條件下,管外冷凝液膜會(huì)達(dá)到湍流態(tài)。當(dāng)處在壓力P≤1.3 MPa的實(shí)驗(yàn)工況,其冷凝液膜流態(tài)均處在層流區(qū)域,當(dāng)處于壓力P=1.6 MPa工況,其尾部液膜已處于湍流區(qū)域。