王曉坤，齊少璞，楊 軍，葉尚尚，王利霞，馮宗蕊，種道彤，賈鴻玉，楊曉燕，劉一哲，楊紅義,*

(1.中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413；2.西安交通大學(xué) 動(dòng)力工程多相流國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049)

鈉冷快堆是第4代核能系統(tǒng)中最成熟的堆型，我國已建成中國實(shí)驗(yàn)快堆(CEFR)，正在建設(shè)600 MW示范快堆(CFR600)。系統(tǒng)程序是重要的反應(yīng)堆設(shè)計(jì)和安全分析軟件，用于工況設(shè)計(jì)、瞬態(tài)分析、確定論事故分析、概率安全分析的熱工水力計(jì)算。

鈉冷快堆系統(tǒng)程序有兩種來源，一是由水堆系統(tǒng)程序改造而來，RELAP5[1]、ATHLET[2]、TRACE[3]、CATHARE[4]等水堆系統(tǒng)程序都擴(kuò)展了鈉冷快堆計(jì)算能力。水堆系統(tǒng)程序的一個(gè)重要功能是分析失水事故(LOCA)，其流動(dòng)換熱計(jì)算通?；诹匠虄上嗔髂Ｐ汀ｂc冷快堆冷卻劑流動(dòng)換熱的主要特點(diǎn)是單相、低壓、高溫、低普朗特?cái)?shù)，基于這些特點(diǎn)開發(fā)專用系統(tǒng)程序就成為鈉冷快堆系統(tǒng)程序的另一來源。美國阿貢國家實(shí)驗(yàn)室(ANL)開發(fā)的SAS4A/SASSYS-1[5]和SAM[6]，法國CEA開發(fā)的OASIS[7]，美國布魯克海文國家實(shí)驗(yàn)室(BNL)開發(fā)的SSC，俄羅斯開發(fā)的GRIF、DINROS、RUBIN，中國原子能科學(xué)研究院開發(fā)的FASYS[8]，華北電力大學(xué)開發(fā)的SAC-CFR[9]，西安交通大學(xué)開發(fā)的THACS[10]等都是典型的鈉冷快堆專用系統(tǒng)程序。這些程序主要用于事故分析，計(jì)算范圍包括堆芯和一回路系統(tǒng)，一部分軟件可計(jì)算二回路和事故余熱排出系統(tǒng)(DHRS)。

CFR600是我國首座自主設(shè)計(jì)的電廠規(guī)模鈉冷快堆，其工況設(shè)計(jì)較CEFR更為復(fù)雜，如多模塊蒸汽發(fā)生器(SG)的協(xié)調(diào)配合，為解決這類實(shí)際的工程問題，中國原子能科學(xué)研究院在CFR600工程科研項(xiàng)目的支持下，開發(fā)了鈉冷快堆系統(tǒng)程序FR-Sdaso，并對(duì)其進(jìn)行了初步驗(yàn)證。相比國內(nèi)外其他程序，F(xiàn)R-Sdaso計(jì)算范圍更廣，包括堆芯、一回路、二回路、三回路、四回路、DHRS，以及必要的控制邏輯。FR-Sdaso的主要模型考慮了我國鈉冷快堆技術(shù)路線和特點(diǎn)，如其鈉池模型能很好地模擬我國快堆鈉池流動(dòng)換熱特性。FR-Sdaso已用于CFR600設(shè)計(jì)[11-12]和安全分析。

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

1.1 中子學(xué)模型

中子學(xué)模型用于計(jì)算裂變功率、衰變功率和反應(yīng)性反饋。FR-Sdaso采用點(diǎn)堆中子學(xué)模型，并以功率代替經(jīng)典方程中的中子通量密度、緩發(fā)中子濃度、裂變產(chǎn)物或錒系元素總量，變形后的裂變功率方程和衰變功率方程分別為：

i=1,…,M

(1)

(2)

其中：N為裂變功率；t為時(shí)間；ρX為反應(yīng)性；β為有效緩發(fā)中子份額；Λ為中子代時(shí)間；λ為衰變常量；C為緩發(fā)中子引發(fā)的裂變功率；M為緩發(fā)中子組數(shù)；H為裂變產(chǎn)物或錒系元素總量與每次衰變釋放能量的乘積；E為各組核素對(duì)應(yīng)的功率份額；J為裂變產(chǎn)物或錒系元素組數(shù)。

反應(yīng)性反饋模型以反饋系數(shù)和特征溫度計(jì)算各項(xiàng)反應(yīng)性，考慮Doppler效應(yīng)、冷卻劑溫度效應(yīng)、燃料軸向膨脹效應(yīng)、堆芯徑向膨脹效應(yīng)、燃料彎曲效應(yīng)，總反應(yīng)性ρX(t)表示為：

ρX(t)=ρe(t)+αfa(Tf(t)-Tf(0))+

αc(Tc(t)-Tc(0))+αb(Tout(t)-Tin(0))

(3)

其中：ρe為外加反應(yīng)性；αfa為燃料軸向膨脹反饋系數(shù)；αfr為堆芯徑向膨脹反饋系數(shù)；αfD為Doppler反饋系數(shù)；αc為冷卻劑溫度反饋系數(shù)；αb為燃料彎曲反饋系數(shù)；Tf為燃料平均溫度；Tin為堆芯入口冷卻劑溫度；Tc為冷卻劑平均溫度；Tout為堆芯出口冷卻劑平均溫度。

1.2 堆芯熱工模型

堆芯熱工模型用于計(jì)算燃料、包殼、冷卻劑的溫度分布，忽略燃料和包殼材料的軸向?qū)岷铜h(huán)向?qū)幔剂虾桶鼩さ臏囟确匠毯?jiǎn)化為極坐標(biāo)下一維擴(kuò)散方程：

(4)

邊界條件為燃料表面和包殼內(nèi)表面的氣隙換熱(忽略氣體的比熱)、包殼外表面與冷卻劑的對(duì)流換熱，即：

(5)

其中：ρ為密度；cp為比定壓熱容；r為徑向坐標(biāo)；k為熱導(dǎo)率；qv為體積功率，忽略包殼釋熱時(shí)，包殼控制方程中qv=0；Tu為燃料表面溫度；Tci為包殼內(nèi)壁溫度；Tco為包殼外壁溫度；Tm為冷卻劑溫度；ru為芯塊外徑；rci為包殼內(nèi)徑；rco為包殼外徑；hc為冷卻劑與包殼對(duì)流換熱系數(shù)；hg為芯包間隙換熱系數(shù)；ku為燃料熱導(dǎo)率；kc為包殼熱導(dǎo)率。

冷卻劑溫度計(jì)算使用多管單通道模型，忽略冷卻劑導(dǎo)熱，其能量守恒方程一維對(duì)流方程如下：

2πrcohc(Tco-Tm)

(6)

其中：A為流道截面積；w為流量；z為軸向坐標(biāo)。

1.3 控制體模型和流網(wǎng)模型

控制體模型是熱工計(jì)算的基本模型，是構(gòu)成流網(wǎng)和過程設(shè)備模型的基礎(chǔ)，F(xiàn)R-Sdaso的控制體默認(rèn)為不可壓縮流體，其模型只有能量守恒方程：

(7)

其中：V為控制體體積；φ為源項(xiàng)。

由控制體模型可派生出管道、閥門、泵等流體部件模型。管道模型包括管道的直徑、長度、阻力系數(shù)等屬性。閥門模型包括阻力系數(shù)和開度等屬性和表示阻力特性的函數(shù)。泵模型包括泵的轉(zhuǎn)速、壓頭、惰轉(zhuǎn)半時(shí)間、惰轉(zhuǎn)到0的時(shí)間等屬性，以及特性曲線、惰轉(zhuǎn)曲線等函數(shù)。FR-Sdaso的泵模型特別考慮了適應(yīng)設(shè)計(jì)初期尚無設(shè)備詳細(xì)參數(shù)的條件下開展工作的需要，輸入條件僅有額定轉(zhuǎn)速下特性曲線和惰轉(zhuǎn)半時(shí)間、惰轉(zhuǎn)到0的時(shí)間，實(shí)際應(yīng)用中可通過安全評(píng)價(jià)確定惰轉(zhuǎn)半時(shí)間和惰轉(zhuǎn)到0的時(shí)間，并以設(shè)計(jì)要求的形式提出，在設(shè)備設(shè)計(jì)中實(shí)現(xiàn)。第一象限(正轉(zhuǎn)正流量)的揚(yáng)程-轉(zhuǎn)速-流量曲線用二次曲線近似表示：

h=aw2+bwn+cn2

(8)

其中：h為揚(yáng)程；n為相對(duì)轉(zhuǎn)速；a、b、c為常數(shù)，可通過泵在額定轉(zhuǎn)速下的特性曲線擬合得到。

采用正切模型計(jì)算泵的惰轉(zhuǎn)曲線[13]：

(9)

其中：tv為惰轉(zhuǎn)半時(shí)間；tw為惰轉(zhuǎn)到0的時(shí)間；s為常數(shù)。

流網(wǎng)模型描述系統(tǒng)構(gòu)成和系統(tǒng)內(nèi)設(shè)備的連接關(guān)系，用于求解各支路的流量。流網(wǎng)由節(jié)點(diǎn)和支路組成，首先基于伯努利方程建立支路方程：

(10)

其中：p1、p2分別為支路進(jìn)出口壓力，即支路兩側(cè)的節(jié)點(diǎn)壓力；Y為支路流通能力，是與支路流通截面、流道長度、流道形狀等有關(guān)的擬合參數(shù)，用于宏觀表示支路的阻力特性；f(w)為源項(xiàng)，如泵的壓頭或自然循環(huán)壓頭等流動(dòng)的驅(qū)動(dòng)力。

將支路方程在t時(shí)刻Taylor展開，并整理得到支路流量的表達(dá)式：

(11)

節(jié)點(diǎn)模型由控制體模型增加節(jié)點(diǎn)連續(xù)方程∑w=0得到，將支路流量方程和節(jié)點(diǎn)連續(xù)方程聯(lián)立得到流網(wǎng)壓力矩陣，通過特征線法[14]求解可得流網(wǎng)中各節(jié)點(diǎn)的壓力和各支路的流量。

鈉池的流動(dòng)換熱呈現(xiàn)明顯的三維特征，但在更宏觀的層面，將熱池、冷池、泵出口、柵板聯(lián)箱抽象為節(jié)點(diǎn)，則鈉池仍可抽象為一維流網(wǎng)，從而可解出堆芯、中間熱交換器、壓力管、泵支撐冷卻、主容器冷卻等主流道和主要輔助流道流量。考慮鈉池三維特征，將熱池和冷池分別劃分為多個(gè)控制體，可模擬鈉池出口溫度響應(yīng)入口溫度變化的時(shí)間效應(yīng)，并在自然循環(huán)工況下計(jì)算密度分布。采用一維流網(wǎng)與多控制體鈉池相結(jié)合的模型，可較好地模擬鈉池的流動(dòng)換熱特性，滿足系統(tǒng)程序分析的需要。鈉池流網(wǎng)和控制體劃分如圖1、2所示。圖2中空心箭頭表示主流道的流動(dòng)方向，紅色表示熱鈉的主流方向，藍(lán)色表示冷鈉的主流方向。

圖1 鈉池流網(wǎng)圖Fig.1 Flow net of sodium pool

圖2 鈉池控制體劃分Fig.2 Control volume division of sodium pool

熱交換器兩側(cè)均可表示為若干個(gè)串聯(lián)控制體，源項(xiàng)為換熱量。SG水側(cè)模型基于一維均相流假設(shè)[15]，采用滑移網(wǎng)格模型[16]，考慮過冷、核態(tài)沸騰、膜態(tài)沸騰、過熱4個(gè)換熱區(qū)，在各區(qū)根據(jù)流量選用適當(dāng)?shù)膶?duì)流換熱關(guān)系式[16]。

1.4 常規(guī)島

FR-Sdaso將快堆系統(tǒng)程序的計(jì)算范圍擴(kuò)展到常規(guī)島的主要設(shè)備，以滿足開展全廠工況設(shè)計(jì)和堆機(jī)匹配設(shè)計(jì)的需要，可用于分析緊急停堆后常規(guī)島排余熱、甩負(fù)荷至廠用電等涉及全廠響應(yīng)的工況。三回路和四回路(循環(huán)水)分別以流網(wǎng)模擬，其中汽輪機(jī)熱慣性小，采用準(zhǔn)穩(wěn)態(tài)模型模擬，級(jí)組的支路模型基于Flugel公式建立[17]：

(12)

其中，下標(biāo)i、o、s分別表示入口、出口、額定值。

常規(guī)島其他模型采用集總參數(shù)模型。汽水分離器以分離效率表示其去除濕蒸汽中飽和水的比例。凝汽器分為殼側(cè)蒸汽區(qū)、殼側(cè)水區(qū)(熱井)、管側(cè)3部分，對(duì)各部分分別建立質(zhì)量守恒和能量守恒模型。給水加熱器兩側(cè)分別為焓值表示的控制體模型，源項(xiàng)為換熱量。除氧器內(nèi)氣相和液相處于相平衡。電功率由汽輪機(jī)功率乘以發(fā)電機(jī)效率得到，汽輪機(jī)的功率由各級(jí)理想焓降乘以級(jí)組效率得到。

1.5 控制模型

FR-Sdaso控制模型包括PI調(diào)節(jié)模型和邏輯模型兩部分。PI模型用于根據(jù)目標(biāo)值調(diào)節(jié)主要控制變量，如根據(jù)功率調(diào)節(jié)控制棒棒位，根據(jù)流量調(diào)節(jié)一回路主泵、二回路主泵，根據(jù)蒸汽發(fā)生器出口鈉溫調(diào)節(jié)給水泵轉(zhuǎn)速和給水調(diào)節(jié)閥開度，根據(jù)蒸汽壓力調(diào)節(jié)主蒸汽閥開度，根據(jù)凝汽器和除氧器液位調(diào)節(jié)凝結(jié)水泵轉(zhuǎn)速等，從而維持給定功率運(yùn)行或?qū)崿F(xiàn)啟動(dòng)、停堆等變工況過程的控制。邏輯模型用于觸發(fā)保護(hù)，SG模塊隔離等根據(jù)觸發(fā)條件自動(dòng)觸發(fā)的動(dòng)作，或在某一指定的時(shí)刻或達(dá)到某一指定條件時(shí)引入故障，如控制棒失控提升、主泵停運(yùn)、失去廠外電源、失去給水等。

2 程序架構(gòu)

FR-Sdaso程序包括前處理、計(jì)算、輸出、數(shù)據(jù)管理4個(gè)模塊。

2.1 前處理模塊

前處理模塊的功能是讀入輸入數(shù)據(jù)并對(duì)數(shù)據(jù)庫進(jìn)行初始化，數(shù)據(jù)庫中的部分變量通過輸入卡讀入，另一部分變量需經(jīng)過初始化計(jì)算得到，如換熱器內(nèi)部的溫度分布等。同時(shí)，F(xiàn)R-Sdaso有續(xù)算功能，此時(shí)數(shù)據(jù)庫中所有的參數(shù)均從輸入卡讀入。

2.2 計(jì)算模塊

計(jì)算模塊是FR-Sdaso的核心功能模塊，根據(jù)輸入卡設(shè)定的條件，逐時(shí)層計(jì)算瞬態(tài)過程的參數(shù)。每一時(shí)層的計(jì)算流程如圖3所示。

2.3 輸出模塊

輸出模塊根據(jù)輸入卡的設(shè)置，以規(guī)定的時(shí)層間隔輸出規(guī)定參數(shù)的值，也可在特定的時(shí)刻輸出特定設(shè)備的詳細(xì)分布參數(shù)。

2.4 數(shù)據(jù)管理模塊

FR-Sdaso通過共享內(nèi)存和輸出文件保存數(shù)據(jù)。共享內(nèi)存中保存3類參數(shù)：第1類是描述電廠結(jié)構(gòu)的參數(shù)，如換熱器的換熱管根數(shù)、換熱管直徑、壁厚等，運(yùn)行過程中一般不允許修改；第2類是運(yùn)行狀態(tài)控制參數(shù)，如引入故障的時(shí)間點(diǎn)或觸發(fā)條件、保護(hù)參數(shù)的觸發(fā)狀態(tài)等，這些參數(shù)將根據(jù)運(yùn)行狀態(tài)變化；第3類是時(shí)層的工藝參數(shù)，如堆芯流量、堆芯出口溫度、泵轉(zhuǎn)速等，也包括換熱器內(nèi)部的溫度分布等參數(shù)，這類參數(shù)在內(nèi)存中僅保存當(dāng)前時(shí)層和上一時(shí)層的值，完成1個(gè)時(shí)層的計(jì)算后，以當(dāng)前時(shí)層值代替上一時(shí)層值，進(jìn)入下一時(shí)層的計(jì)算。輸出文件根據(jù)輸入卡的設(shè)置，保存第1類和第2類參數(shù)中的部分或全部數(shù)據(jù)，并按給定的時(shí)層間隔保存第3類參數(shù)中的指定數(shù)據(jù)。

圖3 每一時(shí)層的計(jì)算流程圖Fig.3 Calculation flowchart of each time layer

3 驗(yàn)證和確認(rèn)(V&V)

FR-Sdaso的V&V過程包括需求評(píng)價(jià)、應(yīng)用范圍分析、物理模型適用性評(píng)價(jià)、關(guān)鍵現(xiàn)象識(shí)別和排序表(PIRT)，測(cè)試矩陣和驗(yàn)證用例分析。需求分析、應(yīng)用范圍評(píng)價(jià)、物理模型評(píng)價(jià)均表明FR-Sdaso程序開發(fā)過程符合預(yù)期。

3.1 PIRT

FR-Sdaso程序的PIRT分析了其用于設(shè)計(jì)和安全分析的包絡(luò)算例，能覆蓋其使用范圍的所有瞬態(tài)。其中一回路主泵卡軸事故的PIRT中堆芯和一回路部分結(jié)果列于表1。根據(jù)事故特征將事故進(jìn)程分為2個(gè)階段：事故初期由事故初因引起的熱工參數(shù)劇烈變化過程(階段1)和自緊急停堆開始電廠向安全狀態(tài)過渡的過程(階段2)。對(duì)階段1和階段2的重要度評(píng)價(jià)和知識(shí)水平評(píng)價(jià)用高(H)、中(M)、低(L)表示，重要度越高表示現(xiàn)象對(duì)計(jì)算結(jié)果的影響越大，知識(shí)水平越高表示對(duì)現(xiàn)象的理論分析越完備。重要度高而知識(shí)水平低的現(xiàn)象需經(jīng)過實(shí)驗(yàn)算例驗(yàn)證。

表1 FR-Sdaso的PIRT示例Table 1 PIRT of FR-Sdaso

3.2 驗(yàn)證算例

FR-Sdaso程序測(cè)試矩陣可覆蓋所有重要度高且知識(shí)水平低或中的關(guān)鍵現(xiàn)象。本文以鳳凰堆壽期末自然循環(huán)實(shí)驗(yàn)基準(zhǔn)例題[14]和CEFR 1臺(tái)一回路主泵停運(yùn)實(shí)驗(yàn)2個(gè)算例進(jìn)行分析。

圖4 堆芯入口/出口溫度變化曲線Fig.4 Change curve of core inlet/outlet temperature

鳳凰堆壽期末自然循環(huán)實(shí)驗(yàn)基準(zhǔn)例題堆芯入口/出口溫度計(jì)算值與實(shí)驗(yàn)值對(duì)比示于圖4。由圖4可知，堆芯入口溫度升高后趨于穩(wěn)定，實(shí)驗(yàn)值受測(cè)點(diǎn)布置和測(cè)量誤差的影響波動(dòng)明顯，計(jì)算值無波動(dòng)，兩者開始升高的時(shí)間、升高的幅度、趨于穩(wěn)定的時(shí)間均符合，說明計(jì)算值反映了失熱阱后堆內(nèi)溫度升高的特點(diǎn)。實(shí)驗(yàn)開始后，受負(fù)反饋影響堆芯出口溫度降低，緊急停堆后因堆芯流量迅速降低，堆芯出口溫度升高，建立自然循環(huán)后堆芯出口溫度趨穩(wěn)。計(jì)算結(jié)果反映了這一過程，但變化幅度大于測(cè)量值，原因是測(cè)點(diǎn)布置在組件出口上方，受攪渾影響實(shí)驗(yàn)值未反映堆芯出口溫度快速變化。對(duì)實(shí)驗(yàn)過程和實(shí)驗(yàn)結(jié)果的分析表明FR-Sdaso能清晰反映實(shí)驗(yàn)過程中失熱阱、二回路建立有效熱阱等特點(diǎn)，可較為可靠地模擬該實(shí)驗(yàn)過程。

CEFR 1臺(tái)一回路主泵切除后，40%功率下，0時(shí)刻切除一回路2環(huán)路主泵，觸發(fā)反應(yīng)堆緊急停堆，一回路和二回路的2環(huán)路主泵轉(zhuǎn)速惰轉(zhuǎn)至0 r/min，一回路1環(huán)路主泵轉(zhuǎn)速線性下降至150 r/min，二回路1環(huán)路主泵轉(zhuǎn)速線性下降至300 r/min。以泵轉(zhuǎn)速為輸入，計(jì)算一回路各流道流量和主要節(jié)點(diǎn)溫度，以驗(yàn)證一回路流網(wǎng)模型的正確性。圖5為泵旁通流道流量計(jì)算值與實(shí)驗(yàn)值的對(duì)比(實(shí)驗(yàn)中僅能測(cè)量旁通流道的流量)，圖6為堆芯出口溫度計(jì)算值與實(shí)驗(yàn)值的對(duì)比。

圖5 泵旁通流道流量變化曲線Fig.5 Change curve of flow in bypass channel of pump

圖6 堆芯出口溫度變化曲線Fig.6 Change curve of core outlet temperature

本算例覆蓋堆芯和一回路系統(tǒng)主要流道的強(qiáng)迫流動(dòng)、堆芯換熱等重要現(xiàn)象。由圖5可見，泵旁通流道流量的FR-Sdaso計(jì)算值與實(shí)驗(yàn)值符合良好，證明FR-Sdaso程序?qū)ε酝鞯懒髁坑?jì)算準(zhǔn)確；堆芯出口溫度計(jì)算值與測(cè)量值變化趨勢(shì)、峰值出現(xiàn)的時(shí)間以及峰值的大小均符合良好，考慮到堆芯出口溫度與堆芯流量的強(qiáng)相關(guān)關(guān)系，證明FR-Sdaso程序?qū)Χ研玖髁坑?jì)算準(zhǔn)確。由此推斷FR-Sdaso程序?qū)σ换芈妨鲃?dòng)計(jì)算準(zhǔn)確。實(shí)驗(yàn)從正常運(yùn)行開始，經(jīng)歷過渡過程，最后穩(wěn)定在不對(duì)稱工況，全過程中計(jì)算值與測(cè)量值都能較好地符合，說明FR-Sdaso程序一回路流動(dòng)模型體現(xiàn)了鈉池流動(dòng)特點(diǎn)，適用于從正常運(yùn)行過渡到不對(duì)稱運(yùn)行的各種工況。本算例也表明FR-Sdaso程序?qū)崿F(xiàn)了基于我國鈉冷快堆技術(shù)特點(diǎn)開發(fā)的初衷。

4 結(jié)論

FR-Sdaso程序是在CFR600工程科研項(xiàng)目支持下開發(fā)的，計(jì)算范圍涵蓋鈉冷快堆一、二、三、四回路，DHRS，控制邏輯，通過開展以PIRT和測(cè)試矩陣為核心的V&V工作，初步證明了其正確性和適用性。目前FR-Sdaso已用于CFR600設(shè)計(jì)和安全分析，滿足CFR600工程需要，實(shí)現(xiàn)了預(yù)期的目標(biāo)。

未來，隨著不斷將CEFR、CFR600和其他快堆的運(yùn)行經(jīng)驗(yàn)，以及理論和實(shí)驗(yàn)研究的成果集成于系統(tǒng)程序，我國鈉冷快堆系統(tǒng)程序水平將不斷提高，逐步建立起有自主知識(shí)產(chǎn)權(quán)和國際影響力的鈉冷快堆系統(tǒng)程序體系，推動(dòng)我國快堆研究和設(shè)計(jì)能力的進(jìn)步。