基于自然循環(huán)回路的非能動安全殼冷卻系統(tǒng)數(shù)值模擬

黃 政

（中國核電工程有限公司，北京 100840）

非能動安全系統(tǒng)對于提高電站的固有安全性、可靠性和經(jīng)濟性有重要意義，一直是研究的熱點之一。在新一代先進核電站的設計中，如AP1000、ESBWR、SMART 等，大多引入了非能動的設計理念［1］。反應堆冷卻劑系統(tǒng)發(fā)生破口時，會有大量蒸汽進入安全殼空間，使安全殼快速升溫升壓。若此時噴淋系統(tǒng)因故障或喪失電源而失效，進入安全殼的熱量無法導出，則會威脅到作為最后一道安全屏障的安全殼的完整性，造成放射性物質(zhì)釋放到環(huán)境的嚴重后果。

本文針對一種采用自然循環(huán)回路的非能動安全殼冷卻系統(tǒng)（PCCS），建立自然循環(huán)回路的一維均相流模型，研究該系統(tǒng)的穩(wěn)態(tài)運行特性，并對影響系統(tǒng)換熱能力的主要參數(shù)進行分析。

1 系統(tǒng)組成與運行方式

本文研究的PCCS原理圖如圖1所示。該系統(tǒng)由6個獨立系列組成，每個系列的主要設備包括：2臺位于安全殼內(nèi)的并聯(lián)換熱器，位于安全殼外與大氣環(huán)境連通的冷卻水箱及相關的管道和閥門。

事故發(fā)生后，大量蒸汽進入安全殼空間，氣體通過冷凝、對流和輻射作用將熱量傳遞給換熱器，依靠冷熱源密度差和高度差產(chǎn)生的驅(qū)動壓頭克服回路沿程阻力損失，建立并維持工質(zhì)的循環(huán)流動，從而將熱量導入作為最終熱阱的冷卻水箱，使安全殼內(nèi)的溫度和壓力降低。當換熱器的熱流密度較高時，流動工質(zhì)甚至會達到對應壓力的飽和溫度而汽化，此時回路上升段及部分加熱段將處于氣液兩相流動狀態(tài)，而下降段的工質(zhì)經(jīng)水箱冷卻后重新恢復為單相狀態(tài)。

圖1 PCCS原理圖Fig.1 Schematic of PCCS

2 數(shù)學模型

2.1 基本假設與節(jié)點劃分

本文主要研究對象為PCCS自然循環(huán)回路和安全殼兩部分。為建立自然循環(huán)回路均相流數(shù)學模型，提出如下假設條件：1）流體為不可壓縮；2）僅考慮流體沿垂直管道的一維流動；3）氣液兩相流速相等，且處于熱力學平衡狀態(tài)，不考慮欠熱沸騰；4）能量方程中，忽略氣液相界面的黏性耗散及流體動能和勢能；5）冷卻水箱能保證冷卻充分，使下降段入口流體為單相狀態(tài)；6）相對于安全殼的響應速度，PCCS自然循環(huán)瞬態(tài)過程可忽略；7）安全殼內(nèi)氣液相之間僅考慮質(zhì)量和能量交換，不考慮動量交換。

PCCS自然循環(huán)回路的坐標系與節(jié)點劃分如圖2所示。僅考慮沿高度z 方向的一維流動，取向上為正，原點位于最低點?；芈犯鞫味季鶆騽澐譃?00個網(wǎng)格節(jié)點，并采用交錯網(wǎng)格的形式，即壓強、焓值及流體物性等參數(shù)定義在網(wǎng)格中心，用整數(shù)編號表示；而速度則定義在網(wǎng)格邊界上，用半整數(shù)編號表示。編號順序與流動方向相同。對于換熱器傳熱管壁面，軸向相應地劃分為100段；徑向僅考慮2個溫度節(jié)點。冷卻水箱采用1個控制體模擬。

對于安全殼，僅采用1個控制體進行模擬，殼內(nèi)分別考慮氣相空間、液相空間、破口質(zhì)能釋放及PCCS帶熱作用。

2.2 控制方程

自然循環(huán)回路的質(zhì)量連續(xù)性方程為：

圖2 PCCS自然循環(huán)回路坐標系與節(jié)點劃分示意圖Fig.2 Schematic of coordinate system and nodalization of natural circulation loop of PCCS

自然循環(huán)回路的動量方程為：

自然循環(huán)回路的能量方程為：

其中：ρ為流體密度，kg／m3；u為流體流速，m／s；p 為壓強，Pa；A 為回路管道面積，m2；τw為管道壁面切應力，N／m2；K 為局部阻力系數(shù)；f 為摩擦阻力系數(shù)；Uw為壁面濕潤周長，m；UH為加熱周長，m；d為管道直徑，m；H 為流體比焓，J／kg；q為換熱器傳熱管的熱流密度，W／m2。

通過導熱傳遞的換熱器熱流量為：

其中：φ 為換熱器熱流量，W；Tcont和Tf分別為安全殼和流體溫度，℃；Tw，in和Tw，out分別為換熱器傳熱管內(nèi)、外壁面溫度，℃；hin和hout分別為傳熱管內(nèi)、外側總換熱系數(shù)，W／（m2·℃）；dw，in和dw，out分別為傳熱管內(nèi)、外徑，m；Aw，in和Aw，out分別為傳熱管內(nèi)、外側換熱面積，m2；λ 為傳熱管壁面的導熱系數(shù)，W／（m·℃）。

冷卻水箱的能量方程為：

其中：W（t）為t時刻冷卻水箱吸收的熱功率，W；cp為比定壓熱容，J／（kg·℃）；mtank為冷卻水箱內(nèi)的冷卻水質(zhì)量，kg；Ttank為冷卻水箱內(nèi)的冷卻水溫度，℃。

安全殼的質(zhì)量方程為：

安全殼的能量方程為：

其中：˙mvap和˙mliq分別為進入氣相空間和液相空間的凈質(zhì)量流量；˙mbrck，v和˙mbrck，l分別為從反應 堆冷卻劑回路破口進入安全殼氣相空間和液相空間的質(zhì)量流量；˙mevap、˙mboil和˙mcond分別為蒸發(fā)、沸騰和冷凝質(zhì)量流量；Wvap和Wliq分別為進入氣相空間和液相空間的凈能量流量；Wbrck，v和Wbrck，l分別為從反應堆冷卻劑回路破口進入安全殼氣相空間和液相空間的流體能量流量；Wevap、Wboil、Wcond和Wconv分別為蒸發(fā)、沸騰、冷凝和對流的能量流量；Wpccs為PCCS系統(tǒng)的換熱功率。

2.3 結構關系式

對于層流（Re＜2 000），f 的表達式為：

對于湍流（Re＞5 000），f 采用Colebrook-White公式［2］計算：

其中：e為管道壁面粗糙度，m；Re為流體雷諾數(shù)。

對于過渡區(qū)域（2 000≤Re≤5 000），則采用lg f-lg Re的線性插值的方法計算f。

對于換熱器傳熱管壁內(nèi)外側的換熱系數(shù)，內(nèi)側采用Churchill公式［3］為：

其中：Nuc為傳熱管外側冷卻劑流體的努賽爾數(shù)；Nufc和Nuna分別為強迫循環(huán)和自然循環(huán)對流換熱的努賽爾數(shù)；Prc、Rec和Grc分別為冷卻劑流體的普朗特數(shù)、雷諾數(shù)和格拉曉夫數(shù)。

對于傳熱管外側，主要傳熱機理為蒸汽在含有不可凝氣體條件下的冷凝，則采用Dehbi公式進行計算：

其中：L 為傳熱管長度，m；p 為混合氣體壓力，atm；w 為不可凝氣體質(zhì)量分數(shù)；T∞為混合氣體主體溫度，℃。

對于流體物性，單相情形參照IAPWS IF97計算。以下給出兩相情形汽水混合物物性計算表達式。

比容為：

動力黏度為：

氣液混合物焓值：

其中：ν為氣液混合物比容，m3／kg；x 為靜態(tài)含汽率；μ 為氣液混合物動力黏度，Pa·s；下標g和l分別表示飽和蒸汽和飽和水。

2.4 數(shù)值求解

對于自然循環(huán)回路，采用控制容積法在編號i和i-1相鄰控制體間積分，對方程進行離散。求解時，先假定一個初始質(zhì)量流量，從下降段入口開始計算，采用牛頓迭代法進行求解，即令：

其中：上標（n）表示第n步迭代；f（˙m）為構造的關于質(zhì)量流量˙m 的函數(shù)；p′tank、ptank分別為目標換熱水箱及實際水箱的壓力，MPa。

對于安全殼，可以把上述方程整理成以氣相溫度和液相溫度為自變量的能量守恒方程，再采用二階牛頓迭代法進行求解：

其中，J 為相應的二階Jacobi矩陣。

2.5 系統(tǒng)參數(shù)和邊界條件

PCCS回路幾何參數(shù)列于表1。僅考慮換熱器與安全殼的換熱，上升段和下降段均假設為絕熱邊界條件。冷卻水箱與大氣連通，壓力為0.10 MPa。

表1 自然循環(huán)回路幾何參數(shù)Table 1 Geometric parameters of PCCS loop

3 計算結果與討論

3.1 穩(wěn)態(tài)運行特性

為研究不同冷熱源溫度邊界條件下的PCCS的工作特性，取安全殼溫度變化范圍為120～150℃，冷卻水箱溫度變化范圍為40～100℃。

圖3為PCCS回路上升段出口處的靜態(tài)含汽率。結果表明，在不同的冷熱源溫度邊界條件下，PCCS可能工作在單相或兩相流動狀態(tài)下。隨冷熱源溫度的升高，工作方式逐漸由單相流動轉變?yōu)閮上嗔鲃?，且上升段出口靜態(tài)含汽率也相應逐漸升高，最大值約為0.035。

圖3 回路上升段出口靜態(tài)含汽率Fig.3 Static quality at outlet of riser for PCCS

圖4為自然循環(huán)回路穩(wěn)態(tài)質(zhì)量流量?？煽闯?，安全殼溫度越高，循環(huán)流量越大。單相自然循環(huán)流動時，流量隨冷卻水箱溫度的升高而減??；進入兩相自然循環(huán)流動時，流量隨冷卻水箱溫度的升高而增大。因此，在發(fā)生相變的過渡點附近，溫度區(qū)域?qū)氖茄h(huán)流量的最小值。

圖4 自然循環(huán)回路穩(wěn)態(tài)質(zhì)量流量Fig.4 Steady-state mass flow rate of natural circulation loop

圖5為不同冷熱源溫度下PCCS總的穩(wěn)態(tài)換熱功率?？煽闯觯簱Q熱功率隨水箱溫度降低線性減小，隨安全殼溫度升高線性增加。但在流動由單相向兩相過渡轉變時，換熱功率會有一階躍變化，使相同冷熱源溫度下兩相的功率高于單相線性外推的功率。這主要是由于進入兩相流動后質(zhì)量流量隨冷熱源溫度的升高而增大造成的（圖4）。

圖5 PCCS穩(wěn)態(tài)換熱功率Fig.5 Steady-state heat transfer power of PCCS

圖6 上升段和下降段最大流速Fig.6 Max fluid velocity of riser and downcomer

圖6為上升段和下降段的最大流速?？煽闯?，單相循環(huán)時，最大流速變化均不明顯（略微下降）；進入兩相后均有所增加，且上升段的最大流速變化更為顯著。單相循環(huán)時，下降段的最大流速大于上升段；而進入兩相循環(huán)后則相反。整個計算溫度范圍內(nèi)，上升段和下降段的最大流速均出現(xiàn)在安全殼和水箱溫度最高的工況下（Tcont=150 ℃，Ttank=100 ℃），分 別 為34.1m／s和1.51m／s。

圖7為安全殼溫度為150 ℃、水箱溫度取4個典型值下的回路節(jié)點溫度和密度分布，其中節(jié)點從下降段入口開始編號。對于溫度，下降段基本不變；換熱器段大致為線性升高；上升段內(nèi)單相情形保持不變，發(fā)生閃蒸后則線性降低。對于密度，單相階段差別不大，進入兩相區(qū)域后隨含汽率的升高下降明顯。圖8為不同安全殼和冷卻水箱溫度下的PCCS回路閃蒸起始點的高度位置?？煽闯?，安全殼和水箱溫度越高，發(fā)生閃蒸的位置越低，兩相段越長。最低位置為4.98m，位于換熱器出口位置。

3.2 瞬態(tài)運行特性

圖7 PCCS回路節(jié)點參數(shù)分布Fig.7 Fluid parameters of PCCS loop

計算瞬態(tài)時，PCCS 共有6 個系列投入運行。通過反應堆冷卻系統(tǒng)破口進入安全殼的質(zhì)能釋放數(shù)據(jù)采用冷管道雙端斷裂事故，安全殼和PCCS 冷卻水箱初始溫度為40 ℃，容積1 000m3。計算得到的安全殼和PCCS系統(tǒng)響應如圖9～11所示?？煽闯?，PCCS投入早期，由于水箱溫度低，安全殼溫度高，PCCS換熱功率較高，能使安全殼快速降溫。隨著水箱中的冷卻水溫度的升高并最終達到飽和沸騰，PCCS功率又逐漸降低，但由于此時衰變熱功率亦降低，因此，安全殼溫度變化緩慢并趨于穩(wěn)定。單相階段，PCCS循環(huán)流量較為穩(wěn)定；進入兩相階段后，流量隨含汽率的升高而上升。

圖8 PCCS回路閃蒸起始點高度Fig.8 Elevation of onset of flashing for PCCS loop

圖9 安全殼和PCCS水箱溫度Fig.9 Temperatures of containment and PCCS tank

圖10 PCCS自然循環(huán)流量和上升段出口靜態(tài)含汽率Fig.10 Mass flow rate of natural circulation and static quality at outlet of riser for PCCS

圖11 PCCS換熱功率Fig.11 Heat removal power of PCCS

4 結論

本文通過建立數(shù)學模型，研究了PCCS的穩(wěn)態(tài)運行特性和回路參數(shù)分布情況，以及事故條件下PCCS和安全殼的瞬態(tài)響應，得到如下結論。

1）PCCS在不同邊界溫度條件下處于單相或兩相自然循環(huán)。換熱功率對于安全殼和換熱水箱溫度近似呈線性分布。

2）安全殼溫度越高，循環(huán)流量越大。在發(fā)生相變的過渡點附近循環(huán)流量最小。

3）單相循環(huán)時，下降段的最大流速大于上升段；而進入兩相循環(huán)后則相反。上升段和下降段的最大流速均出現(xiàn)在安全殼和水箱溫度最高的工況下。

4）PCCS能在噴淋系統(tǒng)故障的事故條件下，在一定時間內(nèi)以非能動方式實現(xiàn)安全殼降溫。長期階段換熱能力降低，需對水箱進行冷卻和補水操作。

［1］ RAO N M，SEKHAR C C，MAITI B，et al.Steady-state performance of a two-phase natural circulation loop［J］.Int J Heat Mass Transfer，2006，33（8）：1 042-1 052.

［2］ BLEVINS R D.Conversion of the Darcy-Weisbach form to the Fanning form used in MELCOR［M］.New York：Krieger Pub Co.，1984.

［3］ DEHBI A，GUENTAY S.A model for the performance of a vertical tube condenser in the presence of noncondensable gases［J］.Nuclear Engineering and Design，1997，177（1-3）：41-52.