王 寧，王任澤，楊亞鵬，馮宗洋，賈林勝，劉一寧，梁博寧

(中國輻射防護研究院，太原 030006)

福島事故表明雖然反應堆嚴重事故發(fā)生概率極低，但是它仍然可能發(fā)生。福島事故前，事故后果的預測和評價始終是核與輻射應急管理的一項重要內(nèi)容，受到各國學術界的重視[1]。福島事故后國際原子能機構(IAEA)提出，應直接根據(jù)堆芯損傷評價的結果指導防護行動[2]。在應急情況下依據(jù)幾個關鍵的反應堆工況參數(shù)實時、快速地預測未來的反應堆工況與事故進程，可以為應急決策提供支持，為場內(nèi)和場外防護行動的實施贏得更多的時間。

由于神經(jīng)網(wǎng)絡的方法具有計算速度快這一顯著優(yōu)點，所以國內(nèi)外學者使用該方法進行了一些事故工況預測的研究。

韓國學者Seung Geun Kim等使用支持向量的計算器，利用核電站狀態(tài)變量輸入在短時間內(nèi)的變化，預測嚴重事故中可能發(fā)生事件的發(fā)生時間[3]。韓國學者Soon Ho Park等使用模糊神經(jīng)網(wǎng)絡預測嚴重事故中的壓力容器水位[4]，因為在嚴重事故中壓力容器的水位無法測量。

韓國學者Dong Yeong Kim等使用人工智能的方法預測安全殼內(nèi)氫氣濃度的變化[5]，該研究將破口位置分為熱管段、冷管段以及蒸汽發(fā)生器傳熱管三種。破口尺寸分為210個步長：對于冷卻劑喪失事故(LOCA)，破口尺寸范圍為雙端斷裂尺寸的1/10000～1；對于蒸汽發(fā)生器，傳熱管破裂數(shù)量從1根變化到210根。破口尺寸的誤差為0.4%。預測安全殼氫氣濃度所用的兩個輸入值為破口尺寸、停堆時長；使用MAAP4的模擬數(shù)據(jù)進行模型的開發(fā)和驗證。

韓國學者Man Gyun Na,Sun Ho Shin等使用概率神經(jīng)網(wǎng)絡方法預測何時堆芯發(fā)生裸露、堆芯出口溫度(CET)溫度何時超過648.9 ℃(一般認為此時開始嚴重事故管理)、壓力容器何時失效[6-7]。

清華大學核能與新能源技術研究院針對高溫氣冷堆示范電站，開發(fā)了一套核應急響應支持系統(tǒng)[8]。該系統(tǒng)有兩個任務：一是基于測量儀器的異常進行故障診斷；二是基于測量儀器的讀數(shù)和操作員的行動預測事故進程。為了更好地執(zhí)行這兩項任務，開發(fā)了動態(tài)貝葉斯網(wǎng)絡，在該網(wǎng)絡中使用貝葉斯機器學習(LBP)算法進行事故診斷和預測。

SESAME由IRSN開發(fā)，用于在LOCA或蒸汽發(fā)生器傳熱管破裂(SGTR)事故下診斷、預測核電站狀態(tài)并分析源項，可以給出主要的事故特征參數(shù)，包括破口位置和尺寸、堆芯失水剩余時間、裂變產(chǎn)物釋放到環(huán)境中的活度隨時間的變化等。

日本學者Tsutomu Ishigami等為了實時預測壓水堆安全殼失效時間，開發(fā)了一項分析技術[9]。預測安全殼失效時間的思路是：熔融堆芯熔穿壓力容器后，冷卻劑經(jīng)壓力容器下腔室泄漏進入安全殼使得安全殼壓力升高，通過簡單的方程，建立未來安全殼壓力與當前安全殼壓力的關系，從而預測未來的安全殼壓力，得到安全殼失效時間。將得到的結果與THALES、STCP等成熟程序的模擬結果進行了對比，驗證了該方法的準確性。

預測的關鍵在于超實時計算，上述研究通過神經(jīng)網(wǎng)絡、簡化方程等方法實現(xiàn)了應急工況進程的預測。本文嘗試使用經(jīng)典公式預測應急工況進程，針對M310反應堆建立了大破口失水事故(LBLOCA)始發(fā)應急工況進程預測方法。根據(jù)事故特點將復雜的一回路進行合理簡化，對質(zhì)量守恒與能量守恒方程進行近似求解。將堆芯活性區(qū)劃分成4個徑向環(huán)和10個軸向?qū)?，?0個柵元。使用實驗關聯(lián)式計算堆芯換熱，得到包殼溫度，然后根據(jù)包殼溫度可判斷堆芯工況。基于本文建立的方法，可以開發(fā)大破口始發(fā)應急工況進程預測程序。

1 事故情景與假設

本文分析的事故情景為大破口失水疊加全廠斷電(SBO)始發(fā)嚴重事故。建模使用的M310反應堆部分參數(shù)列于表1。

基于事故情景做如下假設：

(1)事故前反應堆處于100%功率穩(wěn)定運行狀態(tài)；

(2)大破口和全廠斷電于零時刻同時發(fā)生；

(3)由于大破口發(fā)生后0.1 s內(nèi)一回路壓力即可降到冷卻劑飽和壓力，瞬間產(chǎn)生大量蒸汽，空泡效應引入的負反應性可使反應堆停堆，故假設0 s時反應堆停堆；

(4)高壓安注、低壓安注、主輔給水失效，非能動中壓安注能夠正常運行；

(5)主泵惰轉(zhuǎn)時間為30 s；

(6)安全殼泄露速率為0.1%/d；

(7)仿真計算以下封頭失效為終止事件，安全殼的響應不在本文的研究范圍內(nèi)；

(8)破口發(fā)生在穩(wěn)壓器所在環(huán)路。

2 應急工況預測方法

2.1 守恒方程與整體思路

流體的質(zhì)量守恒和能量守恒方程是模型的基礎，守恒方程的形式取決于所采用的兩相流模型。本文使用基于兩相分相流模型的守恒方程，如式(1)～(4)所示[10]。

表1 建模所用的部分參數(shù)Tab.1 Part of parameters of the model

①氣相質(zhì)量守恒(連續(xù)性方程)

(1)

②液相質(zhì)量守恒方程

(2)

式中，ρl為液相密度，kg/m3；Ul為液相流速，m/s。

③氣相能量守恒方程

(3)

④液相能量守恒方程

1.2.2 疼痛護理。皰疹可伴隨不同程度或不同類型的疼痛感，如刀割樣疼痛、燒灼樣等，尤其是眼周皰疹，疼痛感極易影響患者身心健康。本報告中患者為刀割樣疼痛，痛感強烈并對護理內(nèi)容實施造成影響。因此，護理人員有必要給予疼痛護理。采用通俗易懂的語言將皰疹相關信息對其進行講解，告知患者疼痛感是由于病毒侵犯神經(jīng)導致，使患者對疾病的認知增強，提高對疼痛的耐受力。如痛感較強，可給予止痛藥處理，有效減少痛感對患者造成的影響。

(4)

求解守恒方程需要初始條件、邊界條件以及結構方程，對守恒方程離散成非線性方程組，進而迭代求解。如果根據(jù)事故特點將復雜的一回路進行合理簡化，對守恒方程進行近似求解，使用經(jīng)典公式近似估算應急工況進程，計算速度將大大加快，可以實時評價嚴重事故應急工況，也可以預測事故工況的發(fā)展。

通過泄漏流量、安注流量計算堆芯水位，然后進行堆芯換熱計算，最后計算堆芯包殼溫度，通過包殼溫度判斷應急工況。

2.2 一回路的簡化

大破口失水事故具有如圖1所示的事故特點[11]，本文根據(jù)這些特點，將大破口失水事故分為噴放、再灌水/再淹沒、“長期冷卻”三個階段，針對每個階段采用不同的計算方法與模型。

圖1 大破口始發(fā)應急工況特點Fig.1 Characters of LBLCOA initiated emergency condition

在大破口失水事故期間，一回路迅速泄壓，穩(wěn)壓器很快失效，故可忽略穩(wěn)壓器的作用；在噴放階段，蒸汽發(fā)生器起到最終熱阱的作用，但是噴放階段持續(xù)時間很短，噴放階段結束后，冷卻劑流經(jīng)蒸汽發(fā)生器U型管的流量很小，可以忽略蒸汽發(fā)生器的影響。為滿足預測的計算速度要求，可將穩(wěn)壓器和蒸汽發(fā)生器忽略，反應堆一回路可等效簡化為一個圓柱體，該圓柱體的體積等于一回路冷卻劑所占體積；在該圓柱體冷管段所對應的高度上，冷卻劑通過破口發(fā)生泄漏。模型示意圖如圖2所示。

圖2 一回路簡化示意圖Fig.2 Simplification of the primary loop

2.3 堆芯換熱計算

堆芯換熱是極為重要的結構方程。有研究表明[12]，嚴重事故期間輻射換熱量對包殼溫度以及氫氣產(chǎn)生量的貢獻不大(如圖3、圖4所示)，所以在堆芯換熱計算中忽略輻射換熱。本文依據(jù)包殼溫度判斷堆芯工況，而應急工況下包殼溫度和燃料芯塊溫度基本一樣，所以可忽略燃料棒徑向溫度分布，不計算徑向?qū)崃俊?/p>

使用式5計算兩相鄰節(jié)點之間的軸向?qū)崃縖13]:

(5)

式中，Φ為沿軸向通過兩相鄰柵元之間的熱流量，W；λ為堆芯柵元的導熱系數(shù)，W/(m·K)；A為柵元的橫截面積，m2；Δt為相鄰柵元的溫度差，℃；Δl為相鄰柵元中心點距離，m。

圖3 氫氣產(chǎn)生量對比[12]Fig.3 Comparison for hydrogen production[12]

圖4 堆芯溫度對比[12]Fig.4 Comparison for core temperature[12]

在噴放、再灌水以及再淹沒階段，堆芯冷卻劑與燃料棒的對流方式屬于強迫對流，使用式(6)～式(8)計算層流強迫對流換熱系數(shù)，使用式(9)計算紊流強迫對流換熱系數(shù)，取兩者最大值作為強迫對流換熱系數(shù)[14]:

Nu=C(n)gdev

(6)

(7)

(8)

Nu=0.023Re0.8Pr0.4

(9)

式中，C(n)為常數(shù)，對于棒束和單管均取4.36；(z-z0)為計算點距離流道入口的長度，m；Dh為流道的水力學直徑，m；Re為雷諾數(shù)；Pr為普朗特數(shù)。

在事故的其他階段，對流方式為自然對流，使用式(10)～(11)計算對流換熱系數(shù)，同樣取兩式計算得到的最大值作為自然對流換熱系數(shù):

Nu=0.18Ra1/4(L/Dh)-1/9

(10)

Nu=0.065Ra1/3(L/Dh)-1/9

(11)

式中，L為流道長度，m；Ra為瑞利數(shù)。

2.4 堆芯節(jié)點劃分

將堆芯劃分為4個徑向環(huán)和14個軸向?qū)?其中10個軸向?qū)游挥诙研净钚詤^(qū))，如圖5所示。軸向?qū)?代表二次支撐組件和堆芯支撐板，軸向?qū)?代表流量分配孔板，軸向?qū)?代表堆芯下柵格板，軸向?qū)?～13表示堆芯活性區(qū)，軸向?qū)?4表示堆芯上柵格板。

圖5 堆芯節(jié)點劃分Fig.5 Node division in the core

2.5 應急工況判斷依據(jù)

根據(jù)包殼溫度可判斷應急工況，判斷依據(jù)列于表2。

表2 應急工況判斷依據(jù) [15]Tab.2 Judgment basis for emergency condition [15]

2.6 臨界流

事故初期一回路內(nèi)外壓差很大，冷卻劑泄漏出一回路管道后瞬間汽化，流動為兩相臨界流，物理現(xiàn)象比較復雜。使用Moody臨界流計算泄漏流量。兩相臨界流動的質(zhì)量流速是空泡份額的線性函數(shù)，如式(12)所示[14]:

(12)

式中，α為空泡份額；ρg為氣相密度，kg/m3；ρl為液相密度，kg/m3；ρm為兩相流密度，kg/m3，ρm=αρg+(1-α)ρl；GC(α)為兩相臨界流質(zhì)量流速，kg/(m2·s)；GC(1)為氣相臨界流質(zhì)量流速，kg/(m2·s)；GC(0)為液相臨界流質(zhì)量流速，kg/(m2·s)。

3 結論

本文針對M310反應堆，建立了大破口失水始發(fā)應急工況進程預測方法?？刂品匠袒趦上喾窒嗔髂Ｐ?。求解守恒方程需要初始條件、邊界條件以及結構方程，對守恒方程進行離散成非線性方程組，進而迭代求解。根據(jù)事故特點將復雜的一回路進行合理簡化，對守恒方程進行近似求解。使用經(jīng)典公式近似估算應急工況進程，計算速度大大加快，可以實時評價嚴重事故應急工況，也可以預測事故工況的發(fā)展。

通過泄漏流量、安注流量計算堆芯水位，然后進行堆芯換熱計算，最后計算堆芯包殼溫度，通過包殼溫度判斷應急工況。堆芯劃分為4個徑向環(huán)和14個軸向?qū)?。使用包殼溫度判斷應急工況。堆芯換熱計算是極為重要的結構方程，使用實驗關聯(lián)式進行計算。

基于本文建立的應急工況預測方法，可以開發(fā)響應的大破口始發(fā)應急工況預測程序。針對其他應急工況的預測方法正在研究中。