葛 魁，王 輝，王明軍,*，田文喜，蘇光輝，秋穗正

(1.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049；2.西安交通大學(xué) 動(dòng)力工程多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049；3.中國(guó)核電工程有限公司，北京 100840)

安全殼是事故工況下防止反應(yīng)堆內(nèi)的放射性物質(zhì)向外泄漏的最后一道屏障，其完整性對(duì)核反應(yīng)堆安全至關(guān)重要。第三代核電技術(shù)[1-2]因具有很高的安全性和經(jīng)濟(jì)性逐漸受到全世界核電行業(yè)的認(rèn)可。第三代核電技術(shù)廣泛采取了非能動(dòng)與能動(dòng)相結(jié)合的安全設(shè)計(jì)理念[3]，滿足了“國(guó)際最高核電安全標(biāo)準(zhǔn)”的要求[4-5]，可以靠核電廠自身的安全系統(tǒng)來(lái)保證安全性，減少發(fā)生事故后人為操作不當(dāng)對(duì)核電站安全性的影響。除普通壓水堆[6]外，非能動(dòng)安全系統(tǒng)還被應(yīng)用于一體化輕水堆[7]中。

黃政[8]基于均相流模型建立了一維非能動(dòng)安全殼熱量導(dǎo)出系統(tǒng)分析程序，利用牛頓迭代法求解，模擬了穩(wěn)態(tài)運(yùn)行及在事故工況下安全殼和非能動(dòng)安全殼熱量導(dǎo)出系統(tǒng)的瞬態(tài)響應(yīng)過(guò)程，得到自然循環(huán)系統(tǒng)的流動(dòng)換熱特性。李軍等[9]基于計(jì)算流體力學(xué)(CFD)軟件對(duì)華龍一號(hào)循環(huán)水箱升溫過(guò)程進(jìn)行了三維流動(dòng)傳熱的數(shù)值模擬。Guo等[10]基于兩相均勻流模型開(kāi)發(fā)了模擬PCS自然循環(huán)的瞬態(tài)計(jì)算程序，并對(duì)PCS的循環(huán)流量、壓降、溫度及傳熱系數(shù)等熱工水力參數(shù)進(jìn)行分析研究。白晉華等[11]基于Relap5程序從啟動(dòng)時(shí)間、運(yùn)行工況的穩(wěn)定性等方面對(duì)多種PCS設(shè)計(jì)方案進(jìn)行了評(píng)價(jià)。

綜上，目前針對(duì)PCS熱工水力分析大多基于簡(jiǎn)單的均相流模型，或基于商用Relap5程序進(jìn)行分析計(jì)算。由于均相流模型本身限制，某些工況下的計(jì)算結(jié)果不能真實(shí)反映系統(tǒng)特性及狀態(tài)參量的變化，而漂移流模型能更好地反映系統(tǒng)特性[12]。因此，本文針對(duì)華龍一號(hào)PCS，基于兩相漂移流模型自主可控開(kāi)發(fā)PCS一維自然循環(huán)瞬態(tài)計(jì)算程序，并利用該程序?qū)CS的熱工水力特性進(jìn)行分析研究。

1 數(shù)學(xué)物理模型

華龍一號(hào)換熱水箱中溫度較低流體由于重力作用沿下降管路流入換熱器，流體流經(jīng)換熱器時(shí)被加熱、膨脹，之后流入上升管路最終流回至換熱水箱。隨著流入換熱水箱的流體溫度不斷升高，當(dāng)水箱中的流體溫度達(dá)到飽和時(shí)，將有部分水蒸氣排出至外界環(huán)境，換熱水箱內(nèi)水位降低。當(dāng)水位降低到一定程度時(shí)，水箱補(bǔ)水功能啟動(dòng)，補(bǔ)水管路開(kāi)始以設(shè)定的流速向水箱中注入低溫的水，直至水箱水位恢復(fù)至原來(lái)高度。程序中假設(shè)工質(zhì)在系統(tǒng)內(nèi)的流動(dòng)是一維的，即工質(zhì)的熱工水力參數(shù)僅沿軸向變化，且忽略了換熱器傳熱管壁沿軸向的導(dǎo)熱。

1.1 控制方程

PCS內(nèi)流體熱工水力參數(shù)的變化遵循質(zhì)量、動(dòng)量和能量守恒[13]的基本規(guī)律，再輔以流體的狀態(tài)方程及相應(yīng)的輔助方程可構(gòu)成封閉的方程組，采用Gear算法求解，就可得到上述的狀態(tài)參數(shù)。為能更為準(zhǔn)確地模擬兩相過(guò)程中的瞬態(tài)特性，本研究采用漂移流模型[13]。相較于其他模型，漂移流模型采用代表兩相介質(zhì)橫向分布的量C0和代表兩相之間局部相對(duì)速度的量Vgj來(lái)描述兩相流動(dòng)特性，且動(dòng)量方程中的壓降除考慮摩擦阻力壓降、重力壓降和加速壓降外，還要考慮漂移流壓降。

質(zhì)量守恒方程為：

(1)

式中：t為時(shí)間，s；α為空泡份額；ρf、ρg分別為液相和氣相的密度，kg·m-3；A為流通面積，m2；W為流量，kg·s-1；Z為流動(dòng)方向長(zhǎng)度，m。

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

(2)

式中：p為壓力，MPa；τ為剪切力，N·m-2；Uh為加熱周長(zhǎng)，m；ρ為流體密度，kg·m-3；g為重力加速度，m·s-2；G為質(zhì)量流速，kg·m-2·s-1；x為含氣率。

在漂移流模型中，動(dòng)量方程表示的壓降除摩擦阻力壓降、重力壓降和加速壓降外，還要考慮漂移流壓降梯度：

(3)

(4)

Vgj和C0具體計(jì)算公式如下：

(5)

(6)

(7)

Ed=1-exp(0.138ρf/ρg-0.23α/μg)

(8)

(gDLPL(ρf/ρg-1))0.5

(9)

(10)

(11)

式中：σ為表面張力系數(shù)，N·m-1；μg和μf分別為氣相、液相的動(dòng)力黏度，N·s·m-2；De為管道的水力直徑，m。

能量守恒方程為：

(12)

式中：hf、hg分別為液相和氣相的焓值，kJ·kg-1；Vf為液相的流速，m·s-1；Vg為氣相的流速，m·s-1；q為熱流密度，W·m-2。

自然循環(huán)質(zhì)量流量方程為：

(13)

式中：B為自然循環(huán)驅(qū)動(dòng)壓頭，MPa；Δpf為回路摩擦壓降，MPa；Δploc為回路局部阻力壓降，MPa；ΔpDG為漂移流壓降，MPa；SLA為回路慣量；W為PCS自然循環(huán)質(zhì)量流量，kg·s-1。

換熱器由多根傳熱管組成，程序中將自然循環(huán)流量平均分配到各傳熱管中，且自然循環(huán)流量計(jì)算中考慮下降管路至多根傳熱管之間的形阻壓降。

1.2 輔助模型

1) 傳熱模型

在PCS中，安全殼內(nèi)流體與換熱器傳熱管外壁面、換熱器傳熱管內(nèi)的流體與傳熱管內(nèi)壁面的換熱，均采用牛頓冷卻定律[14-17]：

Q=h·A·ΔT

(14)

式中：Q為熱流量，W；ΔT為傳熱管內(nèi)壁面與傳熱管內(nèi)流體的溫差，K。

安全殼內(nèi)流體與換熱器傳熱管外壁面通過(guò)含不凝結(jié)氣體的蒸汽冷凝過(guò)程進(jìn)行傳熱，冷凝換熱系數(shù)采用Dehbi關(guān)系式[18-19]進(jìn)行計(jì)算：

hDehbi=1.25L0.05((3.7+28.7pb)-

(2 443+458.3pb)lgWa)(T1-Tw)-0.25

(15)

式中：hDehbi為冷凝換熱系數(shù)，W/(m2·K)；L為傳熱管管長(zhǎng)，m；pb為安全殼壓力，atm；Wa為安全殼內(nèi)混合氣體中空氣質(zhì)量分?jǐn)?shù)；Tl為安全殼流體溫度，K；Tw為傳熱管壁面溫度，K。公式適用范圍為：壓力0.15～0.45 MPa，管長(zhǎng)0.3～3.5 m，壁面過(guò)冷度10～50 ℃。

換熱器傳熱管內(nèi)的流體與傳熱管內(nèi)壁面的換熱根據(jù)不同情況選取不同的換熱系數(shù)關(guān)系式。對(duì)于單相流動(dòng)換熱系數(shù)采用D-B公式計(jì)算：

Nu=0.023Re0.8Prn

(16)

式中，加熱流體時(shí)n=0.4，冷卻流體時(shí)n=0.3。

對(duì)于兩相流動(dòng)換熱采用詹斯-洛特斯(Jens-Lottes)公式計(jì)算：

(17)

式中：p為壓力，MPa；q為熱流密度，W/m2。

高溫流體把熱量傳給換熱器傳熱管，同時(shí)傳熱管內(nèi)的冷流體在管內(nèi)流動(dòng)，帶走一部分能量，使換熱器傳熱管冷卻。換熱器傳熱管管壁導(dǎo)熱方程為：

α2F2(Tw-T2)

(18)

PCS中冷卻劑在換熱器吸收安全殼內(nèi)熱量后流進(jìn)換熱水箱，與換熱水箱中流體混合，導(dǎo)致?lián)Q熱水箱內(nèi)水溫升高。當(dāng)換熱水箱的水溫達(dá)到飽和時(shí)，將有部分水蒸氣排出至外界環(huán)境，使換熱水箱水位降低。

2) 阻力模型

單相流動(dòng)的摩擦阻力系數(shù)f[10]為：

(19)

兩相流動(dòng)摩擦壓降的計(jì)算要比單相的復(fù)雜得多，直接計(jì)算有一定難度。通常采用先計(jì)算全液相壓降梯度，然后再乘以一個(gè)因子[13]：

(20)

兩相摩擦壓降倍增因子選用McAdams關(guān)系式[20]：

(21)

式中，vg、vf分別為氣相和液相的比體積，m3/kg。

程序針對(duì)管排式換熱器進(jìn)行計(jì)算，并假設(shè)工質(zhì)的流動(dòng)是一維的，自然循環(huán)總流量平均分給每根換熱器管道。每根換熱器管道的重力壓降、摩擦壓降、漂移流壓降的計(jì)算方法與其他部件管道類似，同時(shí)換熱器的進(jìn)口控制體、出口控制體要考慮形阻壓降。

3) 數(shù)值方法

(22)

y(t0)=y0

(23)

由于核動(dòng)力系統(tǒng)的復(fù)雜性，非線性常微分方程組通常是剛性方程組。本文在求解這類方程組時(shí)選用Gear算法[21-22]，它采用向后差分的隱式方法，是求解剛性問(wèn)題的通用方法，該方法中配備了Adams方法和Gear剛性方法，可根據(jù)情況進(jìn)行變階或變步長(zhǎng)積分，具有較好的穩(wěn)定性。

2 程序開(kāi)發(fā)

華龍一號(hào)PCS示意圖如圖1所示，主要結(jié)構(gòu)包括上升管路、下降管路、換熱器、換熱水箱及水箱補(bǔ)水管路。其中PCS換熱器通過(guò)自然循環(huán)將事故后安全殼內(nèi)的熱量導(dǎo)出至換熱水箱，水箱補(bǔ)水管路在換熱水箱的水位過(guò)低時(shí)投入運(yùn)行，為換熱水箱補(bǔ)水。

圖1 PCS示意圖

圖2為程序計(jì)算流程，首先讀取PCS各管道幾何參數(shù)、安全殼隔間參數(shù)、初始熱工水力狀態(tài)等重要參數(shù)，之后對(duì)程序中各管道的溫度、焓值等進(jìn)行初始化。在程序計(jì)算中，首先計(jì)算流量、焓值的動(dòng)態(tài)導(dǎo)數(shù)，之后調(diào)用Gear算法進(jìn)行求解，并對(duì)參數(shù)進(jìn)行更新。當(dāng)運(yùn)行到設(shè)定的停止時(shí)間時(shí)，程序結(jié)束。

圖2 程序計(jì)算流程

圖3為程序的結(jié)構(gòu)層次圖，主程序調(diào)用數(shù)據(jù)輸入模塊、初始化模塊、動(dòng)態(tài)計(jì)算模塊和系統(tǒng)變量模塊。系統(tǒng)變量模塊用于COMMON模塊的定義，方便各模塊間的變量傳遞。動(dòng)態(tài)計(jì)算模塊分為微分方程數(shù)值求解模塊及在每一步長(zhǎng)結(jié)束后輸出數(shù)據(jù)的數(shù)據(jù)輸出模塊。通過(guò)計(jì)算動(dòng)態(tài)導(dǎo)數(shù)，形成雅各比矩陣，并調(diào)用Gear算法來(lái)求解微分方程組。動(dòng)態(tài)導(dǎo)數(shù)計(jì)算模塊主要計(jì)算各控制體的壓降導(dǎo)數(shù)、焓值導(dǎo)數(shù)兩部分。此外，動(dòng)態(tài)導(dǎo)數(shù)計(jì)算模塊還調(diào)用輔助模塊進(jìn)行物性、換熱系數(shù)、摩擦阻力等的計(jì)算。

3 程序驗(yàn)證

將一維自然循環(huán)瞬態(tài)計(jì)算程序計(jì)算結(jié)果與Relap5程序計(jì)算結(jié)果進(jìn)行對(duì)比，以驗(yàn)證程序計(jì)算結(jié)果的合理性。Relap5程序節(jié)點(diǎn)圖如圖4所示，整個(gè)回路的建模包括了PCS換熱器、換熱水箱、上升管路水平段、上升管路豎直段、下降管路水平段及下降管路豎直段。將回路中下降管路豎直段用部件125模擬，下降管路水平段用部件130模擬；上升管路豎直段用部件140模擬，上升管路水平段用部件145模擬，各部件中包含了所含彎頭的阻力系數(shù)。Relap5中部件135為換熱器，帶有熱構(gòu)件用于提供功率，一維自然循環(huán)瞬態(tài)計(jì)算程序也使用恒定功率進(jìn)行計(jì)算。

程序計(jì)算中，管道控制體數(shù)量劃分會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生一定影響。由于上升管路豎直段可能會(huì)出現(xiàn)兩相流動(dòng)的情形，因此需劃分更多的節(jié)點(diǎn)。本文測(cè)試了4個(gè)算例，分別將上升管路豎直段劃分為5、10、15、20控制體管道進(jìn)行計(jì)算。4個(gè)算例的換熱器流體溫度如圖5所示。

研究結(jié)果表明，在前期的單相階段，控制體劃分影響很小，各算例的計(jì)算結(jié)果基本重合。各算例的差別主要在后期。5控制體管道算例在后期直接發(fā)散，無(wú)法繼續(xù)計(jì)算；10控制體管道算例的計(jì)算結(jié)果與15、20控制體管道算例有較為明顯差別；15控制體管道與20控制體管道算例計(jì)算結(jié)果基本重合。

圖3 程序結(jié)構(gòu)層次圖

圖4 Relap5程序節(jié)點(diǎn)圖

圖5 不同上升管路豎直段控制體數(shù)目下的換熱器流體溫度

在本次驗(yàn)證中，換熱器模塊沿流動(dòng)方向被劃分為5個(gè)熱工水力控制體；換熱水箱模塊沿流動(dòng)方向被劃分為5個(gè)水力學(xué)控制體；由于上升管路豎直段可能會(huì)出現(xiàn)兩相流動(dòng)的情形，因此沿流動(dòng)方向被劃分為20個(gè)熱工水力控制體，上升管路水平段被劃分為5個(gè)熱工水力控制體；下降管路水平段和豎直段均沿流動(dòng)方向被劃分為5個(gè)熱工水力控制體。利用Relap5程序和本文開(kāi)發(fā)的程序計(jì)算了1個(gè)換熱器功率為332 kW恒定值的工況，環(huán)境壓力為0.1 MPa，回路初始水溫為50 ℃，計(jì)算持續(xù)時(shí)間為50 000 s。

圖6示出PCS換熱器流體溫度。由圖6可見(jiàn)，本文程序與Relap5程序的溫度變化趨勢(shì)相同，且換熱器溫度在計(jì)算后期均基本達(dá)到穩(wěn)定值，本文程序計(jì)算得到的換熱器底部(控制體1)和頂部(控制體5)的流體溫度分別為378 K和397 K，而Relap5程序計(jì)算得到的相應(yīng)值為373 K和395 K。兩個(gè)程序計(jì)算結(jié)果差異較小，產(chǎn)生差異的原因可能是由于Relap5中使用的某些輔助模型與本文程序不同，以及兩流體模型與漂移流模型的差異所導(dǎo)致。由于Relap5計(jì)算中重要的輔助模型，如兩相摩擦壓降模型、換熱模型等都無(wú)法自由修改，而本文程序針對(duì)PCS自然循環(huán)這一物理現(xiàn)象選取了更加合適的模型，因此也會(huì)導(dǎo)致計(jì)算結(jié)果產(chǎn)生了一定差異。

圖6 換熱器流體溫度

圖7 自然循環(huán)流量

圖7示出兩個(gè)程序計(jì)算得到的自然循環(huán)流量。在計(jì)算后期，隨換熱器流體溫度的上升，系統(tǒng)內(nèi)會(huì)出現(xiàn)兩相流，導(dǎo)致流量出現(xiàn)震蕩。由于Relap5程序使用兩流體模型，而本文程序采用漂移流模型，因此流量變化有所差別。在計(jì)算后期，兩個(gè)程序計(jì)算得到的流量均基本達(dá)到穩(wěn)定值，本文程序計(jì)算得到穩(wěn)定流量為2.34 kg/s，Relap5程序計(jì)算得到的穩(wěn)定流量為2.43 kg/s。

驗(yàn)證結(jié)果表明，本文開(kāi)發(fā)的一維PCS自然循環(huán)回路程序能用于分析核反應(yīng)堆事故下PCS投入運(yùn)行過(guò)程中重要參數(shù)的變化，評(píng)估PCS帶走安全殼內(nèi)熱量的能力。

4 PCS熱工水力特性分析

利用本文開(kāi)發(fā)的程序?qū)Σ煌瑩Q熱功率下PCS長(zhǎng)期運(yùn)行條件下的熱工水力特性進(jìn)行分析研究。圖8示出不同的換熱功率下PCS換熱水箱出口流體的溫度及換熱器進(jìn)出口控制體內(nèi)的流體溫度。由于未考慮換熱水箱與環(huán)境空氣的對(duì)流換熱，換熱水箱內(nèi)水溫持續(xù)上升，最終達(dá)到飽和狀態(tài)，所以換熱水箱出口的流體溫度保持不變，為環(huán)境壓力下的飽和溫度。不同加熱功率下，換熱器進(jìn)口溫度變化不大，由于流體在換熱器內(nèi)被加熱，換熱器出口流體的溫度大于進(jìn)口溫度，且隨換熱功率呈近似線性增長(zhǎng)的趨勢(shì)。換熱功率為292.5 kW時(shí)換熱器進(jìn)出口溫差約為7 K，而當(dāng)換熱功率為2 044 kW時(shí)換熱器進(jìn)出口溫差增長(zhǎng)到了約12 K。

圖9示出不同換熱功率下PCS傳熱管內(nèi)流體的對(duì)流換熱系數(shù)及PCS自然循環(huán)流量。由圖10可見(jiàn)，換熱系數(shù)和自然循環(huán)流量均隨換熱功率的增大而增加。在292.5～2 044 kW范圍內(nèi)，換熱系數(shù)在1 000 kW·m-2·K-1左右，而流量變化較大，由292.5 kW下的約9 kg·s-1上升至2 044 kW下的約36 kg·s-1。

圖8 不同位置的溫度

圖10示出不同換熱功率下PCS上升管路豎直段的含氣率。計(jì)算中將上升管路豎直段沿高度方向劃分為20個(gè)控制體，1號(hào)控制體在下方，20號(hào)控制體在上方。當(dāng)換熱功率較小時(shí)，上升管路豎直段內(nèi)流體保持為單相水。隨換熱器換熱功率的增大，在豎直上升管道的出口附近會(huì)首先出現(xiàn)兩相的情況，這是由于隨高度的增加，管道內(nèi)壓力降低，管道內(nèi)的流體溫度達(dá)到了相應(yīng)壓力下的飽和溫度，發(fā)生沸騰。當(dāng)換熱功率進(jìn)一步增大，出現(xiàn)兩相流的控制體數(shù)目不斷增加(即出現(xiàn)兩相流的管道長(zhǎng)度越長(zhǎng))，且出現(xiàn)兩相流的控制體內(nèi)空泡份額逐漸增大。

圖11示出瞬態(tài)計(jì)算時(shí)不同換熱功率下?lián)Q熱水箱內(nèi)的水位。水箱補(bǔ)水管路是本文程序針對(duì)華龍一號(hào)PCS開(kāi)發(fā)的，它可在水箱水位較低時(shí)為水箱注水，保證PCS的正常運(yùn)行。即在瞬態(tài)計(jì)算中，水箱溫度達(dá)到飽和后，由于水的蒸發(fā)和蒸汽的排放導(dǎo)致水箱水位不斷降低，當(dāng)打開(kāi)水箱補(bǔ)水開(kāi)關(guān)后，水位降至某一設(shè)定值(2.5 m)時(shí)，水箱補(bǔ)水功能啟動(dòng)，補(bǔ)水管路開(kāi)始以設(shè)定的流速向水箱中注入低溫的水，直至水箱水位恢復(fù)至原來(lái)高度。由圖11可見(jiàn)，換熱功率越大，水箱水位下降越快。水箱補(bǔ)滿時(shí)，水位在一段時(shí)間內(nèi)維持不變的原因是低溫水注入水箱后，水箱內(nèi)的水并未達(dá)到飽和溫度，未發(fā)生蒸發(fā)。

圖9 換熱系數(shù)及自然循環(huán)流量

圖10 上升管路豎直段的含氣率

圖11 換熱水箱水位

5 結(jié)論

本文針對(duì)華龍一號(hào)PCS的結(jié)構(gòu)和特點(diǎn)，基于兩相漂移流模型建立了一套適用于計(jì)算PCS熱工水力現(xiàn)象的數(shù)學(xué)物理模型，并輔之以傳熱模型、阻力模型等，開(kāi)發(fā)了專門適用于華龍一號(hào)PCS的瞬態(tài)熱工水力分析程序。經(jīng)過(guò)與Relap5程序結(jié)果的對(duì)比，二者誤差較小，驗(yàn)證了本文程序的正確性和可靠性。采用該程序?qū)CS內(nèi)的關(guān)鍵熱工水力參數(shù)進(jìn)行分析計(jì)算，得到了自然循環(huán)流量、換熱功率、溫度及含氣率等的分布規(guī)律。

本文開(kāi)發(fā)的PCS瞬態(tài)熱工水力分析程序?yàn)槭鹿屎蟀踩珰さ陌踩治鎏峁┝丝煽抗ぞ?，?duì)PCS的設(shè)計(jì)和改善及對(duì)PCS冷卻能力的評(píng)估均有重要的現(xiàn)實(shí)意義，且為后續(xù)開(kāi)發(fā)能模擬帶有PCS的安全殼內(nèi)熱工水力行為的程序打下基礎(chǔ)。