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

(1.西安交通大學 核科學與技術(shù)學院，陜西 西安 710049；2.西安交通大學 動力工程多相流國家重點實驗室，陜西 西安 710049；3.中國核電工程有限公司，北京 100840)

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

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

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

1 數(shù)學物理模型

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

1.1 控制方程

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

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

(1)

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

動量守恒方程為：

(2)

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

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

(3)

(4)

Vgj和C0具體計算公式如下：

(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分別為氣相、液相的動力黏度，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ū)動壓頭，MPa；Δpf為回路摩擦壓降，MPa；Δploc為回路局部阻力壓降，MPa；ΔpDG為漂移流壓降，MPa；SLA為回路慣量；W為PCS自然循環(huán)質(zhì)量流量，kg·s-1。

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

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)流體與換熱器傳熱管外壁面通過含不凝結(jié)氣體的蒸汽冷凝過程進行傳熱，冷凝換熱系數(shù)采用Dehbi關(guān)系式[18-19]進行計算：

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

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

(15)

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

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

Nu=0.023Re0.8Prn

(16)

式中，加熱流體時n=0.4，冷卻流體時n=0.3。

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

(17)

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

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

α2F2(Tw-T2)

(18)

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

2) 阻力模型

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

(19)

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

(20)

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

(21)

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

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

3) 數(shù)值方法

(22)

y(t0)=y0

(23)

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

2 程序開發(fā)

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

圖1 PCS示意圖

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

圖2 程序計算流程

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

3 程序驗證

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

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

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

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

圖4 Relap5程序節(jié)點圖

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

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

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

圖6 換熱器流體溫度

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

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

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

4 PCS熱工水力特性分析

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

圖9示出不同換熱功率下PCS傳熱管內(nèi)流體的對流換熱系數(shù)及PCS自然循環(huán)流量。由圖10可見，換熱系數(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上升管路豎直段的含氣率。計算中將上升管路豎直段沿高度方向劃分為20個控制體，1號控制體在下方，20號控制體在上方。當換熱功率較小時，上升管路豎直段內(nèi)流體保持為單相水。隨換熱器換熱功率的增大，在豎直上升管道的出口附近會首先出現(xiàn)兩相的情況，這是由于隨高度的增加，管道內(nèi)壓力降低，管道內(nèi)的流體溫度達到了相應壓力下的飽和溫度，發(fā)生沸騰。當換熱功率進一步增大，出現(xiàn)兩相流的控制體數(shù)目不斷增加(即出現(xiàn)兩相流的管道長度越長)，且出現(xiàn)兩相流的控制體內(nèi)空泡份額逐漸增大。

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

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

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

圖11 換熱水箱水位

5 結(jié)論

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

本文開發(fā)的PCS瞬態(tài)熱工水力分析程序為事故后安全殼的安全分析提供了可靠工具，對PCS的設計和改善及對PCS冷卻能力的評估均有重要的現(xiàn)實意義，且為后續(xù)開發(fā)能模擬帶有PCS的安全殼內(nèi)熱工水力行為的程序打下基礎。