宿吉強，孫中寧，張東洋

(1.哈爾濱工程大學(xué) 核安全與仿真技術(shù)國防重點學(xué)科實驗室，黑龍江 哈爾濱 150001；

2.中國船舶重工集團公司 第703研究所，黑龍江 哈爾濱 150078)

新一代核反應(yīng)堆系統(tǒng)廣泛采用非能動安全殼冷卻系統(tǒng)以保持安全殼的完整性。對于基于雙層混凝土安全殼的非能動安全殼冷卻系統(tǒng)，位于安全殼內(nèi)的冷凝器是其最關(guān)鍵的設(shè)備之一，而安全殼內(nèi)的不凝性氣體(主要為空氣)的運動與分布對冷凝器的影響至關(guān)重要。在蒸汽冷凝過程中，大量空氣的存在會導(dǎo)致傳熱的嚴(yán)重惡化，所以必須對含空氣的蒸汽冷凝開展研究[1]。

目前對于含不凝性氣體蒸汽冷凝的研究主要有實驗研究和理論分析兩大類。實驗研究得到的經(jīng)驗關(guān)聯(lián)式具有較高的實用價值[2-5]，但經(jīng)驗關(guān)聯(lián)式的適用條件苛刻，而具有廣泛適用性的經(jīng)典方法和理論尚有待形成。對于工程應(yīng)用，通過數(shù)值模擬的方法獲得評估數(shù)據(jù)是一種經(jīng)濟且快速的途徑，同時數(shù)值模擬也可獲得實驗中較難得到的流場特性及動態(tài)特性的預(yù)測結(jié)果。

含不凝性氣體蒸汽冷凝的數(shù)值模擬研究，主要利用CFD軟件和計算機語言編程相結(jié)合的方式實現(xiàn)。目前主要的數(shù)值模擬方法有兩類：一類是用基本的物理定律求解在冷凝表面的傳熱傳質(zhì)，在不凝性氣體存在條件下建立一個兩相邊界層，從能量守恒和質(zhì)量守恒原理出發(fā)提出一套數(shù)值計算模型，在數(shù)學(xué)方程的基礎(chǔ)上解釋不凝性氣體影響蒸汽冷凝的機理[6]；一類是采用已發(fā)展的傳熱傳質(zhì)關(guān)系式的積分計算，將這些關(guān)系式轉(zhuǎn)化為相應(yīng)控制方程源項，應(yīng)用于臨近冷凝表面的單元網(wǎng)格內(nèi)，從而實現(xiàn)冷凝過程的數(shù)值模擬[7]。第1類方法由于需要大量的編程及相當(dāng)細(xì)致的表格，不適用于工程應(yīng)用，而第2類方法可迅速直觀地模擬冷凝過程中流場的變化情況，具有較高的工程應(yīng)用價值，因而本文以此為基礎(chǔ)對含空氣的冷凝進行數(shù)值模擬。

為驗證數(shù)值計算模型的準(zhǔn)確性，并得到相應(yīng)的流場分布結(jié)果，本文對Su等[8]和Cornet等[9]含空氣蒸汽冷凝的實驗進行數(shù)值仿真模擬，旨在為含空氣蒸汽冷凝的傳熱傳質(zhì)研究奠定基礎(chǔ)。

1 實驗系統(tǒng)與實驗方法

文獻[8]的實驗是在壓力為0.2～0.6 MPa、空氣質(zhì)量分?jǐn)?shù)為20%～80%的條件下進行的。實驗?zāi)Ｐ褪居趫D1。含空氣蒸汽在長2 000 mm的豎直冷凝管的外表面進行冷卻，混合氣體存儲在高3.55 m的碳鋼圓柱罐體內(nèi)，其直徑為565 mm。蒸汽由罐體下側(cè)入口射入，在入口與冷凝管之間存在兩層均氣孔板，冷凝管內(nèi)通有冷卻水。

2 計算模型

2.1 冷凝模型及理論假設(shè)

含空氣蒸汽在豎直壁面冷凝的二維物理模型坐標(biāo)系如圖2所示。平行和垂直于冷凝壁面的坐標(biāo)分別為y和x，對應(yīng)的速度分量分別為v和u。主流氣體速度為u，溫度為T，空氣質(zhì)量分?jǐn)?shù)為Wa，液膜層厚度為δf，空氣積聚層厚度為δg，考慮浮生力，引入重力加速度g。

圖1 實驗?zāi)Ｐ秃唸D

圖2 坐標(biāo)系及相關(guān)物理量

為了對冷凝過程進行有效的求解計算，本文對含空氣蒸汽冷凝過程作了以下4點假設(shè)：

1) 冷凝過程中忽略液膜的影響，凝結(jié)液隨著冷凝過程被移除，以研究不凝性氣體在壁面附近的積聚及主流區(qū)的流場分布。

2) 混合氣體中的蒸汽為其分壓下飽和蒸汽，飽和蒸汽壓力ps由飽和蒸汽溫度Ts唯一確定，根據(jù)文獻[10]，飽和蒸汽壓力與飽和蒸汽溫度的關(guān)系式為：

(1)

3) 冷凝僅發(fā)生在臨近冷凝壁面的第1層網(wǎng)格內(nèi)，氣相邊界層內(nèi)飽和蒸汽受冷降溫后，還有可能出現(xiàn)液滴。而參考文獻[11]，與大量的壁面冷凝相比，氣相邊界層內(nèi)液滴的存在對蒸汽的質(zhì)量擴散僅產(chǎn)生極其微小的影響，因而對氣相邊界層中形成的液滴進行了忽略。

4) 將蒸汽空氣混合流體處理成單相雙組分流體，水蒸氣和空氣組成的雙組分氣體混合物假設(shè)為理想氣體，推導(dǎo)其密度計算式為：

(2)

式中：ρmix為混合氣體密度，kg/m3；R為摩爾氣體常數(shù)，R=8.314 5 J/(mol·K)；Wv為水蒸氣質(zhì)量分?jǐn)?shù)；Mv、Ma分別為水蒸氣、空氣摩爾質(zhì)量，Mv=18 kg/mol，Ma=28.97 kg/mol。雙組分混合物的質(zhì)擴散系數(shù)由下式計算：

D=D0(T0/T)αp0/p

(3)

式中：T0、p0分別為標(biāo)準(zhǔn)狀態(tài)下的溫度、壓力，T0=273.15 K，p0=101 325 Pa；D0為標(biāo)準(zhǔn)狀態(tài)下質(zhì)擴散系數(shù)，水蒸氣和空氣混合時其值為2.56×10-5m2/s；α取值范圍為1.5～2.0，參考文獻[12]取1.81。

2.2 含空氣蒸汽冷凝的實現(xiàn)

圖3 壁面單元示意圖

基于前文的假設(shè)，可將含空氣蒸汽冷凝過程簡化為如圖3所示的過程：在冷凝壁面第1層網(wǎng)格內(nèi)，混合氣體與壁面的對流換熱導(dǎo)致流體溫度降低，在冷凝壁面單元內(nèi)比較當(dāng)?shù)卣羝謮簆v與壁面溫度對應(yīng)的飽和蒸汽壓力ps的大小來判斷是否發(fā)生冷凝，若pv>ps，則冷凝發(fā)生，產(chǎn)生凝液量m0，同時放出熱量H0，凝液及潛熱被從系統(tǒng)中移除；反之凝液量為零，僅存在對流換熱。對流換熱量由FLUENT計算得到，冷凝的相變過程會引起混合氣相的質(zhì)量、能量、動量及組分損失，這些損失可利用編譯UDF(user defined function)在冷凝壁面的控制方程中添加源項來模擬，使用UDM(user defined memory)存儲臨近壁面各單元內(nèi)的源項值，可用于計算某一時刻蒸汽釋放的潛熱。

模擬過程中，冷凝傳熱系數(shù)的獲取采用TOSQON提出的理論[13]，將潛熱與顯熱混合物通過氣液界面向冷凝液傳遞的熱量分為對流換熱過程的顯熱和伴隨凝結(jié)傳質(zhì)過程的潛熱兩部分。依據(jù)文獻[14]定義m0的計算式為：

m0=CihoA(Tb-Tcw)/hlg

(4)

與之相對應(yīng)的冷凝焓H0為：

H0=m0(cp,mixTb-cp,airTref)

(5)

式中：Ci為調(diào)整系數(shù)，W/(m2·K)；ho為實驗得到的冷凝傳熱系數(shù)；A為近壁面網(wǎng)格單元面積，m2；hlg為飽和蒸汽的汽化潛熱，J/kg；cp,mix、cp,air分別為混合氣體和空氣的比定壓熱容，J/(kg·K)；Tref為焓為0時的參考溫度，Tref=273.15 K。

式(4)、(5)中除了Tcw，其余均為罐內(nèi)主流的物理參數(shù)。選擇Su等[12]提出的關(guān)聯(lián)式作為ho的計算公式，調(diào)試獲得的調(diào)整系數(shù)為0.98，得到m0的計算式為：

m0=0.98[10 189.3+90 416.4p-

(4 314.4+46 537p)lg(100Wa)]·

(Tb-Tcw)0.4/hlg

(6)

在模擬實驗工況時，采用瞬態(tài)逼近穩(wěn)態(tài)策略，選取基于壓力的全隱式耦合算法。根據(jù)前文的假設(shè)，壁面冷凝只發(fā)生在近壁面單元。FLUENT并不能分辨出近壁面單元，因而需設(shè)定判斷邏輯流程。為此，本文對求解近壁面單元冷凝的邏輯流程進行了設(shè)定：求解之前，通過壁面ID查找出近壁面單元的ID，并開辟存儲空間儲存這些ID，對流體域單元掃描計算控制方程時，判定該單元ID是否與存儲空間ID相符合。如果符合，說明該單元是近壁面單元，進行冷凝計算；如果不符，說明該單元不是近壁面單元，不進行冷凝。需要注意，在每個時間步長的計算收斂后，需使用UDF保存壁面單元的壓力和水蒸氣質(zhì)量分?jǐn)?shù)及壁面溫度，為下一時間步長是否發(fā)生冷凝提供依據(jù)。圖4為UDF編譯結(jié)構(gòu)的流程圖。

2.3 幾何模型與處理方法

對二維幾何模型進行網(wǎng)格劃分，簡化去除了冷卻水的進、出口段，全部采用四邊形結(jié)構(gòu)化網(wǎng)格。均氣孔板的小孔處利用尺寸函數(shù)加密，臨近冷凝壁面的第1層網(wǎng)格高度為2 cm，通過網(wǎng)格無關(guān)性解的驗證，確定二維模型的網(wǎng)格量為29 442。

圖4 冷凝程序結(jié)構(gòu)

在對模型進行計算時，對密度、動量、能量項和組分項均采用二階迎風(fēng)格式進行離散；選擇二階瞬態(tài)隱式關(guān)系式求解雙組分氣體混合過程的解；湍流模型選擇RNGk-ε模型。

冷凝壁面采用無滑移壁面邊界條件，給定壁面溫度分布，冷凝壁溫的差異會對傳熱產(chǎn)生較大影響[12]，因而模擬過程中采用實驗的壁溫參數(shù)。壁面溫度通過擬合傳熱管軸向的外壁面溫度的實驗值得到。溫度邊界由壁溫實驗值的擬合關(guān)系式得到，并通過UDF的定義邊界宏加載到冷凝壁面。除蒸汽入口外冷凝罐的其他壁面均按絕熱處理。入口邊界條件為質(zhì)量流量入口，通過UDF定義入口質(zhì)量流率等于UDM儲存的每個單元凝液量的總和，用以模擬每一穩(wěn)定工況的恒壓過程。

3 數(shù)值模型驗證及結(jié)果分析

3.1 模型驗證

Su等[12]和TOSQON的參數(shù)范圍[13]與安全殼事故工況下的參數(shù)相比，具有一定的代表性，因而本文在CFD中對該實驗進行了模擬。

為驗證本文方法的準(zhǔn)確性，對TOSQON實驗的含空氣蒸汽冷凝的工況條件進行了數(shù)值模擬。TOSQON實驗是在一高為4.8 m、直徑為1.5 m的圓柱形壓力罐內(nèi)進行的，蒸汽及不凝性氣體(空氣、氦氣)入射口位于距離罐底2.1 m處，入射口直徑為41 mm，純蒸汽入射時蒸汽流量為0.001 1 kg/s，入口蒸汽溫度為126 ℃。冷凝壁面位于距離罐底2.391 m和4.391 m之間的區(qū)域，冷凝壁面溫度為101.8 ℃，上、下兩側(cè)的非冷凝壁面溫度分別為122.0 ℃和123.5 ℃，網(wǎng)格數(shù)為6 319。實驗初始條件及模擬結(jié)果列于表1，初始條件誤差較小。同時模擬結(jié)果顯示，數(shù)值計算預(yù)測得到的溫度、速度以及水蒸氣含量同實驗結(jié)果的誤差分別為±1.89 K、±4.3%和±10%，體現(xiàn)出本工作計算方法在預(yù)測相似工況實驗方面體具有較高的準(zhǔn)確度。

表1 TOSQON初始實驗條件及模擬結(jié)果

對Su等實驗中3個壓力的9組工況進行二維數(shù)值模擬，所選實驗工況的主要物理參數(shù)列于表2。預(yù)測結(jié)果顯示，冷凝傳熱系數(shù)最大相對誤差為28.59%，主流溫度及壁面過冷度的最大誤差為3.75 ℃，空氣質(zhì)量分?jǐn)?shù)最大相對誤差為0.15%，進一步驗證了本文所用模型的準(zhǔn)確性。

表2 仿真計算結(jié)果與參考值對比

3.2 流場分析

圖5為p=0.4 MPa、Wa=34.7%工況下的穩(wěn)態(tài)模擬結(jié)果。在速度云圖中，罐內(nèi)的高速區(qū)形成于靠近冷凝罐和傳熱管的壁面附近，且冷凝罐下半部分空間的流速高于上半部分空間。在灌頂和罐底空間流場低速特點明顯，說明在這兩處流體較少參與罐內(nèi)摻混，形成滯流區(qū)。流線圖及局部矢量圖顯示冷凝罐內(nèi)形成了多處渦流，從流線圖可看出，傳熱管兩側(cè)的渦流形狀修長，渦流中心位于傳熱管下端的左、右兩側(cè)，即在傳熱管兩側(cè)各自形成了自然對流。

a——速度云圖；b——流線圖；c——局部速度矢量

圖6為p=0.4 MPa下不同空氣質(zhì)量分?jǐn)?shù)時的溫度分布云圖。從圖6可看出，隨入射蒸汽量的減小，均氣孔板以上空間的溫度場不均性逐漸減弱，在傳熱管底端20 cm以內(nèi)有低溫流體的積聚。從圖中的等值線可看出，在冷熱流體交匯處形成了較大的溫度梯度和明顯的溫度交界面，較高的溫度使蒸汽和混合氣體在混合過程中進行快速的熱量交換，讓蒸汽達到其當(dāng)?shù)胤謮合碌娘柡蜏囟?。高溫氣體沿冷凝罐內(nèi)壁上升，同時與低溫流體混合，在冷凝罐上部空間形成了均勻的溫度場。在圖6b中黑框所標(biāo)出的位置，有兩塊溫度高于周圍流體的高溫區(qū)，蒸汽沿遠(yuǎn)離入口側(cè)的冷凝罐內(nèi)壁面向上流動與罐內(nèi)低溫流體混合。結(jié)合圖6a、c發(fā)現(xiàn)，受到下降冷流體在均汽孔板上部堆積的影響，蒸汽的上升混合過程是變化的，可能在靠近入口側(cè)上升與罐內(nèi)氣體混合，也可能是在遠(yuǎn)離入口側(cè)上升與罐內(nèi)氣體混合，或在兩側(cè)同時上升與罐內(nèi)氣體混合。

a——Wa=34.7%；b——Wa=57.2%；c——Wa=83.3%

4 結(jié)論

證明了利用實驗結(jié)果實現(xiàn)含空氣蒸汽冷凝數(shù)值模擬的可行性，該方法與其他模擬方法相比具有簡單、可靠、快速等特點，具有一定的工程應(yīng)用價值。利用本文模擬方法對TOSQON實驗及Su等實驗共10組工況進行了數(shù)值模擬，模擬結(jié)果同TOSQON實驗值的相對誤差在±10%以內(nèi)，同Su等實驗結(jié)果平均溫度的最大誤差為3.75 ℃，平均壓力和空氣質(zhì)量分?jǐn)?shù)的相對誤差均在1%以內(nèi)，計算結(jié)果與實際情況吻合較好，該數(shù)值計算方法對于相似工況以及相近物理模型內(nèi)蒸汽在含空氣的環(huán)境中冷凝過程的預(yù)測有指導(dǎo)意義。通過對數(shù)值模擬結(jié)果的分析，得到了速度場和溫度場等參數(shù)細(xì)致的分布狀態(tài)。

參考文獻：

[1] ROSA J C. Review on condensation on the containment structures[J]. Progress in Nuclear Energy, 2009, 51(1): 32-66.

[2] UCHIDA H, OYAMA A, TOGO Y. Evaluation of post-incident cooling systems of light-water power reactors[C]∥Proceedings of International Conference on Peaceful Uses of Atomic Energy. [S.l.]: [s.n.], 1965: 93-102.

[3] TAGAMI T. Interim report on safety assessments and facilities establishment project[R]. Japan: Atomic Energy Research Agency, 1965.

[4] LIU H, TODREAS N E, DRISCOLL M J. An experimental investigation of a passive cooling unit for nuclear plant containment[J]. Nuclear Engineering and Design, 2000, 199(3): 243-255.

[5] DEHBI A A. Analytical and experimental investigation of the effects of non-condensable gases on steam condensation under turbulent natural convection conditions[D]. USA: Massachusetts Institute of Technology, 1990.

[6] ANDREANI M, PUTZ F, DURY T. Application of field codes to the analysis of gas mixing in large volumes[C]∥IAEA Technical Committee Meeting on the Implementation of Hydrogen Mitigation Techniques and Filtered Containment Venting. Koeln, Germany: [s.n.], 2001: 18-21.

[7] MARTIN J M, JIM M A, MARTIN F. Comparison of film condensation models in presence of non-condensable gases implemented in a CFD code[J]. Heat Mass Transfer, 2005, 41(11): 961-976.

[8] SU J Q, SUN Z N, FAN G M. Experimental study of the effect of non-condensable gases on steam condensation over a vertical tube external surface[J]. Nuclear Engineering and Design, 2013, 262(5): 201-208.

[9] CORNET P, MALET J, PORCHERON E. TOSQAN experimental results of the air-steam phase[C]∥OECD International Standard Problem on Containment Thermal-Hydraulics ISP-47. France: Institut de Radioprotection et de Sureté Nucléaire, Saclay, 2002.

[10] 劉志剛,劉咸定,趙冠春. 工質(zhì)熱物理性質(zhì)計算程序編制及應(yīng)用[M]. 北京：科學(xué)技術(shù)出版社,1992.

[11] 王朝陽,屠傳經(jīng). 霧的形成對蒸汽和不凝性氣體混合物的強迫對流膜狀凝結(jié)換熱的影響[J]. 高?；瘜W(xué)工程學(xué)報,1989,3(1):75-84.

WANG Chaoyang, TU Chuanjing. The effect of fog formation on forced convective condensation of vapor noncondensable gas mixture[J]. Journal of Chemical Engineering of Chinese Universities, 1989, 3(1): 75-84(in Chinese).

[12] MASSMAN W J. A review of the molecular diffusivities of H2O, CO2, CH4, CO, O3, SO2, NH3, N2O, NO, and NO2in air and N2near STP[J]. Atmospheric Environment, 1998, 32(6): 1 111-1 127.

[13] COLBURN A P. Design of cooler condensers for mixtures of vapors with non-condensing uses[J]. Industry Engineering Chemistry, 1934, 26(11): 1 178-1 182.

[14] IVO K, MIROSLAV B. Modeling of containment atmosphere mixing and stratification experiment using a CFD approach[J]. Nuclear Engineering and Design, 2006, 236(14): 1 682-1 692.