GASFLOW程序液膜模型開發(fā)及初步驗證

王方年1，沈峰1，程旭1,2，黃興冠2

(1.國核科學(xué)技術(shù)研究院，北京102209；2.上海交通大學(xué)，上海200240)

摘要：本文基于三維CFD安全殼程序GASFLOW開發(fā)了熱構(gòu)件壁面上的液膜覆蓋與蒸發(fā)模型。通過選定AP1000大破口事故序列，采用耦合了液膜模型的GASFLOW程序分析了AP1000核電廠安全殼內(nèi)溫度壓力響應(yīng)及其非能動安全殼冷卻系統(tǒng)(PCS)的性能，并與相同事故序列下WGOTHIC、MELCOR、CONTAIN等程序的計算結(jié)果進(jìn)行比較。結(jié)果表明，耦合了液膜模型的GASFLOW程序可用于分析PCS的熱工水力行為，其基本功能滿足計算需要。

關(guān)鍵詞：GASFLOW程序；液膜模型；AP1000；非能動安全殼冷卻系統(tǒng)；安全殼響應(yīng)

中圖分類號：TL364.3 文獻(xiàn)標(biāo)志碼：A

收稿日期：2014-07-16;修回日期：2015-01-22

基金項目：國家科技重大專項資助項目(2013ZX06004008)

作者簡介：王方年(1985—)，男，江西九江人，工程師，碩士，核能科學(xué)與工程專業(yè)

doi：10.7538/yzk.2015.49.11.2044

Development and Preliminary Verification

of GASFLOW Code Coupled with Film Model

WANG Fang-nian1, SHEN Feng1, CHENG Xu1,2, HUANG Xing-guan2

(1.StateNuclearPowerResearchInstitute,Beijing102209,China;

2.ShanghaiJiaoTongUniversity,Shanghai200240,China)

Abstract:A model of heat structure wall surface film coverage and evaporation was developed based on 3 dimensional CFD containment code GASFLOW. The containment temperature and pressure response and the passive containment cooling system (PCS) performance of AP1000 during large break LOCA were analyzed by GASFLOW code coupled with film model. The calculation results were compared with the calculated results of other containment codes WGOTHIC, MELCOR, CONTAIN under the same accident scenario. The results show that the modified GASFLOW code coupled with film model is feasible to analyze the thermal-hydraulic behavior in PCS of PWR and the basic functions can meet the requirements for the calculation.

Key words：GASFLOW code; film model; AP1000; passive containment cooling system; containment response

GASFLOW程序是由美國LANL和德國FZK共同開發(fā)的三維計算流體力學(xué)程序[1-2]。它是一個最佳流場估計工具，可描述流場中以下三維流動的現(xiàn)象：氣體的擴散、混合分布與分層，氫氣燃燒和火焰擴散，不可凝氣體的分布對本地凝結(jié)與蒸發(fā)的影響，氣溶膠的夾帶、傳輸與沉降等。

在AP1000核電廠安全殼設(shè)計中，由于采用了重力、壓差和自然循環(huán)等非能動理念[3]，非能動安全殼冷卻系統(tǒng)(PCS)與傳統(tǒng)安全殼相比增添了許多新現(xiàn)象，尤其是鋼制安全殼外壁面液膜的蒸發(fā)冷卻現(xiàn)象。國內(nèi)外研究人員曾使用多種程序(包括WGOTHIC、MELCOR、CONTAIN、COCOSYS等)對AP1000核電廠安全殼進(jìn)行了分析[4-6]，但使用三維CFD程序分析的很少?；贏P1000核電廠安全殼的特殊性，GASFLOW能否適用于AP1000核電廠PCS分析需重新評估。GASFLOW3.2只能對所有構(gòu)筑物表面輸入初始液膜厚度，液膜不能持續(xù)不斷的補充到鋼制安全殼外壁面，因此液膜會隨著時間推移而蒸干，無法實現(xiàn)安全殼的長期冷卻。目前，最新版的GASFLOW3.5仍無法實現(xiàn)這一功能[1-2]。同時，GASFLOW用靜態(tài)液膜厚度(液膜厚度為常數(shù))來實現(xiàn)蒸發(fā)，無法表示液膜厚度隨流動方向慢慢變薄的過程，也不能區(qū)分干濕區(qū)域，無法模擬液膜在鋼制安全殼外壁面覆蓋率的變化。

本文開發(fā)熱構(gòu)件壁面上的液膜覆蓋與蒸發(fā)模型，并耦合到GASFLOW3.2中，對大破口事故下AP1000核電廠安全殼響應(yīng)和PCS排熱性能進(jìn)行分析。

1液膜模型開發(fā)

液膜模型開發(fā)主要分兩部分。首先是對GASFOW定義的鋼制安全殼壁面網(wǎng)格進(jìn)行重構(gòu)，改變鋼制安全殼表面的液膜計算順序，使其能模擬液膜流動的走向；然后通過耦合液膜覆蓋與蒸發(fā)模型對具有一定計算順序的每一個網(wǎng)格上的液膜進(jìn)行計算處理，求解液膜的質(zhì)量守恒方程、能量守恒方程，計算出液膜覆蓋率、排熱能力等。

1.1GASFLOW網(wǎng)格重構(gòu)

GASFLOW中熱構(gòu)件壁面網(wǎng)格按照式(1)從小至大的順序進(jìn)行編號，記為master[1]：

i+mboff(iblk)

(1)

其中：i、j、k為定義網(wǎng)格時的空間坐標(biāo)編號；mboff(iblk)為描述建模對象的編號。故在GASFLOW中，各墻體壁面計算順序并不是模擬AP1000液膜流動所需的先后計算順序，因此，本文對GASFLOW安全殼壁面網(wǎng)格進(jìn)行重構(gòu)，以一定的方式存儲，然后以一定的順序提供程序計算，簡單地實現(xiàn)液膜從里至外、自上而下流動的模擬。

將GASFLOW中每個壁面的網(wǎng)格進(jìn)行了重新定義，通過一4維數(shù)組master=(i，j，k，windex)將所有壁面網(wǎng)格存儲起來，然后按從里至外、自上而下進(jìn)行循環(huán)計算，該數(shù)組包含每個墻體具體位置及其網(wǎng)格面朝向等信息。其中：i、j、k定義了空間坐標(biāo)x、y、z3個方向的網(wǎng)格編號，即定義了具體位置；windex記錄著每個網(wǎng)格面的朝向，本文用1、2、3來分別代表x、y、z方向，用-1、-2、-3來分別代表-x、-y、-z方向。

1.2液膜覆蓋與蒸發(fā)模型

在完成GASFLOW壁面網(wǎng)格重組基礎(chǔ)上，通過開發(fā)液膜覆蓋與蒸發(fā)模型，計算不同位置、不同時刻的水膜厚度、水膜覆蓋率；同時也能從液膜厚度的角度模擬溪流變窄過程及蒸干現(xiàn)象；最終較為準(zhǔn)確地實現(xiàn)在安全殼壁面上由內(nèi)到外的熱量傳遞。

本文是基于Nusselt液膜理論分析液膜流動情況的，其基本假設(shè)是：液膜沿豎直方向勻速流動，根據(jù)重力與摩擦力等力的平衡求解速度。計算時認(rèn)為1個網(wǎng)格上的液膜速度是均勻的，且是液膜厚度的函數(shù)；沿流動方向上的下一個網(wǎng)格，因為蒸發(fā)作用厚度產(chǎn)生變化，其速度也跟隨變化，但網(wǎng)格上的液膜速度是同一值。

本文的液膜覆蓋與蒸發(fā)模型只是針對鋼制安全殼外壁面液膜開發(fā)的，實際情況中液膜也同時存在于鋼制安全殼內(nèi)壁面。鋼制安全殼內(nèi)壁面的液膜計算采用GASFLOW中原有模型，即靜態(tài)液膜模型，認(rèn)為液膜厚度是受限的，液膜通過冷凝達(dá)到一定厚度后將保持不變，繼續(xù)凝結(jié)的部分將被移除出內(nèi)壁面。

1) 平均液膜速度與液膜厚度

利用Nusselt理論[7]求解豎直壁面層流降膜的液膜厚度及液膜速度。動量方程為：

(2)

邊界條件為：

(3)

定義豎直壁面單位寬度上的質(zhì)量流量為：

豎直壁面上的液膜厚度為：

(4)

管外壁降膜流動的液膜厚度為：

其中：δpl為豎直平板的液膜厚度；εR為管外壁面液膜厚度與管道外徑的比值。安全殼屬于大直徑的管外降膜流動，其外壁面液膜厚度(μm級)遠(yuǎn)小于安全殼直徑(40 m或以上)，即εR≈0，可認(rèn)為δ=δpl，因此式(4)仍適用于安全殼外壁面液膜流動。式(4)在推導(dǎo)時液膜流態(tài)屬于層流，但實驗證明[8]也符合湍流計算，甚至在相當(dāng)一部分區(qū)域比層流計算更準(zhǔn)確。Nusslet理論關(guān)系式在層流及波動層流的范圍內(nèi)計算出的液膜厚度均大于其他經(jīng)驗公式，液膜越厚可導(dǎo)致覆蓋率的減小，因此本文采用Nusslet理論關(guān)系式計算是保守的。

2) 臨界液膜厚度

西屋公司的LST、SST、Flat Plate實驗[4]表明，當(dāng)液膜流速大于某一臨界值時液膜寬度保持不變，當(dāng)液膜流速降至該值后流速保持不變，液膜寬度減小直至完全蒸發(fā)。由于潤濕面積越大，蒸發(fā)流量越多，因此較大的流速限值使得液膜寬度更早減小，是保守的。本文基于文獻(xiàn)[4]的液膜蒸發(fā)實驗結(jié)果，取具有包絡(luò)性的保守值作為臨界液膜線質(zhì)量流量參考值，再代入式(4)可得到臨界液膜厚度。

3) 液膜覆蓋率及蒸發(fā)流量的計算

定義m為液膜質(zhì)量流量(即未蒸發(fā)的液膜流量)，W為濕周(液膜濕潤寬度)，Z為壁面網(wǎng)格的高度(或長度)，則有[4]：

m=WΓ

(5)

液膜流動主要分為兩個過程，即液膜寬度不變的過程(W不變，液膜厚度隨蒸發(fā)變薄)和單位長度上的質(zhì)量流量不變的過程(Γ不變，液膜已經(jīng)到達(dá)臨界厚度，并保持不變)。

總的液膜蒸發(fā)流量為液膜初始化總流量與未蒸發(fā)流量之差，即：

(6)

其中：mevp為液膜蒸發(fā)流量；minitial為液膜的初始注入流量；mn為第n個網(wǎng)格的液膜流量。

每一網(wǎng)格上液膜蒸發(fā)流量是采用傳熱傳質(zhì)類比的形式計算得到的[2]。當(dāng)液膜向下流動并蒸發(fā)，液膜厚度減薄至低于臨界液膜厚度后，鋼制安全殼外壁面網(wǎng)格上的液膜將出現(xiàn)干、濕兩部分區(qū)域，在程序模型內(nèi)認(rèn)為該網(wǎng)格的液膜寬度變窄，液膜高度不變，原理如圖1所示。

圖1　液膜覆蓋率變化原理圖 Fig.1　Schematic of variation for film coverage rate

由于液膜模型能計算出每一個網(wǎng)格的液膜寬度，因此液膜覆蓋率Rcoverage可表示為：

(7)

其中：Wn、Zn分別為第n個網(wǎng)格上液膜的寬度與高度(或長度)；W0、Z0分別為網(wǎng)格的寬度與高度初始化值。每個網(wǎng)格上的液膜蒸發(fā)流量為：

(8)

(9)

因液膜蒸發(fā)帶走的能量可表示為：

(10)

其中，hfg為液膜汽化潛熱。

氣體與熱構(gòu)件表面對流換熱量為[10]：

(11)

其中：A為網(wǎng)格面積；Tbulk,gas為大空間混合氣體的溫度；Tsurf,structure為熱構(gòu)件表面的溫度。

輻射換熱量qradiation計算采用的是壁面-氣流的簡單輻射換熱模型(壁面假設(shè)是黑體)：

(12)

其中：σ為Stefan-Boltzmann常數(shù)；εbulk,gas為氣體發(fā)射率；αbulk,gas為氣體吸收率。

2程序初步驗證

本文以分析AP1000大破口事故下安全殼響應(yīng)為例，初步驗證耦合了液膜模型的GASFLOW程序的適用性與可信度。

2.1AP1000建模

1) AP1000 GASFLOW網(wǎng)格劃分及幾何建模

根據(jù)AP1000核電廠安全殼的實際布置情況進(jìn)行了GASFLOW建模[3-4]。AP1000核電廠安全殼體積龐大且殼內(nèi)結(jié)構(gòu)復(fù)雜，若要建立一個與實際結(jié)構(gòu)完全一致的模型則需采用很精細(xì)的網(wǎng)格，這將使計算代價變得無法接受。為尋求計算代價與計算精確度的平衡，在保持自由容積、傳熱面積的參數(shù)相對準(zhǔn)確的情況下對安全殼構(gòu)筑物進(jìn)行了簡化。徑向網(wǎng)格由墻體和設(shè)備位置進(jìn)行定位，軸向網(wǎng)格由各層地板和設(shè)備位置進(jìn)行定位，盡可能保證網(wǎng)格界面與墻體、地板、設(shè)備表面重合以節(jié)省網(wǎng)格數(shù)目；圓周向采用均勻網(wǎng)格劃分。本文安全殼模型采用的是圓柱坐標(biāo)系下結(jié)構(gòu)化正交網(wǎng)格劃分方法。表1列出AP1000的網(wǎng)格劃分方式及其數(shù)目，平均網(wǎng)格體積為1.18 m3。圖2示出AP1000核電廠安全殼的GASFLOW模型。AP1000 GASLOW模型中安全殼內(nèi)自由容積(不包括外部自然循環(huán)空間)為59 200 m3，鋼制安全殼的高度為65.6 m，直徑為39.6 m。

表1　AP1000 GASFLOW網(wǎng)格劃分

2) AP1000大破口事故質(zhì)能釋放源項

AP1000大破口事故下，質(zhì)能釋放位置是在操作平臺之下的蒸汽發(fā)生器隔間，質(zhì)能釋放源項參見文獻(xiàn)[5]，包含：1) 從壓力容器側(cè)噴放出來的兩相流體；2) 從蒸汽發(fā)生器側(cè)流出的高溫蒸汽；3) 自動卸壓系統(tǒng)ADS-4釋放在蒸汽發(fā)生器隔間的蒸汽。

3) 其他初始邊界條件

安全殼內(nèi)空氣及熱構(gòu)件的初始溫度為48.89 ℃，空氣初始壓力為0.108 MPa，除內(nèi)置換料水箱隔間外，相對濕度為0。外部環(huán)境大氣初始溫度為46.11 ℃，壓力為0.101 MPa，相對濕度為22%。鋼制安全殼外壁面液膜蒸發(fā)的驅(qū)動力為界面與大氣內(nèi)水蒸氣的濃度差，濕度用于計算大氣內(nèi)水蒸氣的濃度。鋼制安全殼與混凝土屏蔽構(gòu)筑物之間的環(huán)形通道自然循環(huán)計算，本文采用定壓邊界條件。

圖2　AP1000安全殼的模型 Fig.2　Geometric model of AP1000 containment

2.2大破口事故分析

為驗證修改后GASFLOW的計算結(jié)果是否正確，與文獻(xiàn)[5，11-12]的計算結(jié)果進(jìn)行對比分析。表2列出AP1000大破口事故下安全殼內(nèi)壓力、溫度的對比。本文計算結(jié)果與WGOTHIC、MELCOR、CONTAIN等集總參數(shù)程序的計算結(jié)果相近，但略有差別。AP1000核電廠安全殼的設(shè)計壓力為0.508 MPa，設(shè)計溫度為148.89 ℃[3]，本文的計算結(jié)果再次證明AP1000當(dāng)前PCS的設(shè)計滿足大破口事故下的安全需求。

表2　大破口事故下安全殼壓力和溫度的結(jié)果對比

圖3示出相同邊界條件下AP1000大破口事故的安全殼內(nèi)壓力的對比。壓力呈現(xiàn)雙峰值形式，前一個峰值是由于壓力容器側(cè)短時間內(nèi)(1～20 s)噴放的大量高溫、高壓兩相流體造成的；在20～100 s間，安全殼內(nèi)壓力變得較為平坦，其原因是這一階段壓力容器側(cè)噴放量減少至較低的水平；在100～1 700 s期間，殼內(nèi)壓力持續(xù)上升，其原因一方面是壓力容器側(cè)釋放的流體焓值上升，另一方面是SG側(cè)開始釋放蒸汽。337 s后PCS投入，壓力曲線此時出現(xiàn)拐點。在1 700 s之后，壓力因SG側(cè)蒸汽釋放結(jié)束而緩慢下降?？梢奝CS的作用一方面是為安全殼的長期排熱提供最終熱阱，另一方面降低安全殼的峰值壓力。

圖3　不同程序計算的安全殼內(nèi)壓力對比 Fig.3　Pressure comparison inside containment by different codes

從圖3各程序計算結(jié)果對比可見，壓力變化趨勢基本一致，集總參數(shù)程序間的結(jié)果較為接近， GASFLOW計算結(jié)果較集總參數(shù)程序計算結(jié)果差別較大。這是由于各程序本身及建模過程存在差異，如GASFLOW采用傳熱傳質(zhì)類比的形式，集總參數(shù)程序一般采用冷凝經(jīng)驗關(guān)系式。

圖4示出不同程序計算的安全殼內(nèi)溫度的對比。集總參數(shù)程序一般情況下將安全殼內(nèi)穹頂當(dāng)作一個結(jié)塊或幾個結(jié)塊處理，取結(jié)塊的平均值；GASFLOW數(shù)據(jù)取自如圖2所示的A、B、C3點溫度的平均值。安全殼穹頂?shù)钠骄鶞囟染谑鹿是捌谫|(zhì)能高速噴放階段(20 s內(nèi))達(dá)到峰值，高溫流體會因噴放速度的影響迅速沖到安全殼穹頂，并隨著熱構(gòu)件對高溫蒸汽凝結(jié)的作用，使得穹頂氣流溫度達(dá)到峰值后又下降，而后因PCS的投入帶走了熱量使得溫度緩慢下降。從長期結(jié)果看，4種程序的計算結(jié)果處于同一水平且具有相同的變小趨勢。

圖4　安全殼內(nèi)穹頂處溫度對比 Fig.4　Dome temperature comparison inside containment

圖5　液膜覆蓋率的變化 Fig.5　Variation of coverage rate for film

2.3PCS排熱性能

圖5、6分別示出鋼制安全殼外壁面液膜覆蓋率的變化及不同高度處液膜厚度的變化。圖6中，a～g的位置如圖2所示。由圖5、6可見：因PCS噴淋尚未啟動，液膜覆蓋率在0～337 s時為零；PCS啟動后，液膜覆蓋率的變化與液膜厚度的變化相對應(yīng)，且大部分監(jiān)測點(a、b、c、d、e)的變化趨勢保持一致；有些區(qū)域噴放前期沒有液膜覆蓋而后出現(xiàn)了液膜覆蓋，并處于臨界液膜厚度或蒸干狀態(tài)(f、g處)?？梢?，本文程序關(guān)于液膜覆蓋與蒸發(fā)模型的耦合具有一定合理性。

圖6　鋼制安全殼外壁面液膜厚度的變化 Fig.6　Variation of film thickness on outside wall surface of steel containment

圖7示出鋼制安全殼外壁面熱流密度的變化。由圖7可見，液膜蒸發(fā)平均熱流密度約在11 kW/m2，PCS環(huán)形通道內(nèi)對流換熱、輻射換熱的平均熱流密度分別在0.5 kW/m2、0.2 kW/m2左右。因此，液膜蒸發(fā)是安全殼內(nèi)熱量導(dǎo)出的主要方式。

圖7　鋼制安全殼外壁面平均熱流密度的變化 Fig.7　Variation of average heat flux on outside wall surface of steel containment

圖8示出鋼制安全殼外壁面不同監(jiān)測點處溫度的變化。鋼制安全殼外壁面液膜溫度超過100 ℃的時間不算長，超過100 ℃的區(qū)域集中在穹頂部分，豎直壁面未超過限值。本文程序暫時不能處理液膜溫度超過100 ℃的情況，這是未來需改進(jìn)之處。

圖8　鋼制安全殼外壁面溫度的變化 Fig.8　Variation of outside wall surface temperature for steel containment

3結(jié)論

1) GASFLOW程序的二次開發(fā)滿足預(yù)設(shè)的計算功能，能夠計算液膜覆蓋率、鋼制安全殼外壁面液膜厚度的變化，能夠模擬分析AP1000核電廠PCS的排熱性能。

2) 對比分析相同事故序列下不同程序的計算結(jié)果，初步驗證程序修改后的計算結(jié)果是可信的。

3) 通過與WGOTHIC等程序計算結(jié)果的對比，再次表明目前AP1000核電廠PCS的設(shè)計在大破口事故下能夠保持安全殼的完整性，滿足設(shè)計要求。

參考文獻(xiàn)：

[1]SPORE J W, TRAVIS J R, ROYL P, et al. GASFLOW: A computational fluid dynamics code for gases aerosols, and combustion, Volume 2: User’s manual, LA-13357-MS, FZKA-5994[R]. USA: LANL; Germany: FZK, 1998.

[2]TRAVIS J R, SPORE J W, ROYL P, et al. GASFLOW: A computational fluid dynamics code for gases aerosols, and combustion, Volume 1: Theory and computational model, LA-13357-

MS, FZKA-5994[R]. USA: LANL; Germany: FZK, 1998.

[3]林誠格. 非能動安全先進(jìn)核電廠AP1000[M]. 北京：原子能出版社，2008.

[4]WOODCOCK J, ANDREYCHEK T, CONWAY L, et al. WGOTHIC application to AP600 and AP1000, Rev. 1, WCAP-15862[R]. America: Westinghouse, 2004.

[5]TILLS J, NOTAFRANCESCO A. Application of the MELCOR code to design basis PWR large dry containment analysis, SAND2009-2858[R]. America: SANDIA, 2009

[6]PHILIPP B, HANS-JOSEF A. Simulation of AP1000’s passive containment cooling with the German containment code system COCOSYS[C]∥Nuclear Energy for New Europe. Bovec, Slovenia: [s. n.], 2011.

[7]SINKUNAS S, GYLYS J. Analysis of laminar liquid film flowing down a vertical surface[C]∥Fourth International Conference on CFD in the Oil and Gas. Norway: [s. n.], 2005.

[8]YU Y Q, WEI S J, YANG Y H, et al. Experimental study of water film falling and spreading on a large vertical plate[J]. Progress in Nuclear Energy, 2012, 54: 22-28.

[9]BIRD R B, STEWART W E, LIGHTFOOT E N. Transport Phenomena[M]. New York: John Wiley and Sons, 1960.

[10]ROHSENOW W M, CHOI H. Heat, mass, and momentum transfer[M]. New Jersey: Prentice-Hall, 1961.

[11]OFSTUN R. AP1000 LOCA containment pressure analysis, Rev.2, APP-SSAR-GSC-523[R]. America: Westinghouse, 2003.

[12]Westinghouse. Shield building air inlet flow area design change, Rev. B, APP-GW-GEE-1307[R]. America: Westinghouse, 2010.