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

(中國(guó)輻射防護(hù)研究院，太原 030006)

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

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

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

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

韓國(guó)學(xué)者M(jìn)an Gyun Na,Sun Ho Shin等使用概率神經(jīng)網(wǎng)絡(luò)方法預(yù)測(cè)何時(shí)堆芯發(fā)生裸露、堆芯出口溫度(CET)溫度何時(shí)超過(guò)648.9 ℃(一般認(rèn)為此時(shí)開(kāi)始嚴(yán)重事故管理)、壓力容器何時(shí)失效[6-7]。

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

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

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

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

1 事故情景與假設(shè)

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

基于事故情景做如下假設(shè)：

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

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

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

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

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

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

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

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

2 應(yīng)急工況預(yù)測(cè)方法

2.1 守恒方程與整體思路

流體的質(zhì)量守恒和能量守恒方程是模型的基礎(chǔ)，守恒方程的形式取決于所采用的兩相流模型。本文使用基于兩相分相流模型的守恒方程，如式(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 疼痛護(hù)理。皰疹可伴隨不同程度或不同類型的疼痛感，如刀割樣疼痛、燒灼樣等，尤其是眼周皰疹，疼痛感極易影響患者身心健康。本報(bào)告中患者為刀割樣疼痛，痛感強(qiáng)烈并對(duì)護(hù)理內(nèi)容實(shí)施造成影響。因此，護(hù)理人員有必要給予疼痛護(hù)理。采用通俗易懂的語(yǔ)言將皰疹相關(guān)信息對(duì)其進(jìn)行講解，告知患者疼痛感是由于病毒侵犯神經(jīng)導(dǎo)致，使患者對(duì)疾病的認(rèn)知增強(qiáng)，提高對(duì)疼痛的耐受力。如痛感較強(qiáng)，可給予止痛藥處理，有效減少痛感對(duì)患者造成的影響。

(4)

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

通過(guò)泄漏流量、安注流量計(jì)算堆芯水位，然后進(jìn)行堆芯換熱計(jì)算，最后計(jì)算堆芯包殼溫度，通過(guò)包殼溫度判斷應(yīng)急工況。

2.2 一回路的簡(jiǎn)化

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

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

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

圖2 一回路簡(jiǎn)化示意圖Fig.2 Simplification of the primary loop

2.3 堆芯換熱計(jì)算

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

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

(5)

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

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

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

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

Nu=C(n)gdev

(6)

(7)

(8)

Nu=0.023Re0.8Pr0.4

(9)

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

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

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

(10)

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

(11)

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

2.4 堆芯節(jié)點(diǎn)劃分

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

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

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

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

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

2.6 臨界流

事故初期一回路內(nèi)外壓差很大，冷卻劑泄漏出一回路管道后瞬間汽化，流動(dòng)為兩相臨界流，物理現(xiàn)象比較復(fù)雜。使用Moody臨界流計(jì)算泄漏流量。兩相臨界流動(dòng)的質(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 結(jié)論

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

通過(guò)泄漏流量、安注流量計(jì)算堆芯水位，然后進(jìn)行堆芯換熱計(jì)算，最后計(jì)算堆芯包殼溫度，通過(guò)包殼溫度判斷應(yīng)急工況。堆芯劃分為4個(gè)徑向環(huán)和14個(gè)軸向?qū)印Ｊ褂冒鼩囟扰袛鄳?yīng)急工況。堆芯換熱計(jì)算是極為重要的結(jié)構(gòu)方程，使用實(shí)驗(yàn)關(guān)聯(lián)式進(jìn)行計(jì)算。

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