湯春桃，畢光文，楊 波

(上海核工程研究設(shè)計(jì)院有限公司，上海 200233)

中子動(dòng)力學(xué)方程是描述核反應(yīng)堆系統(tǒng)瞬態(tài)過(guò)程中堆芯動(dòng)態(tài)特性的基本控制方程之一。目前工業(yè)級(jí)設(shè)計(jì)程序主要采用兩群瞬態(tài)中子擴(kuò)散模型求解點(diǎn)堆動(dòng)力學(xué)或三維時(shí)空動(dòng)力學(xué)方程，該方法的成功運(yùn)用在核電廠低功率物理啟動(dòng)試驗(yàn)動(dòng)態(tài)控制棒價(jià)值測(cè)量和反應(yīng)堆事故分析中發(fā)揮了關(guān)鍵作用。隨著新型反應(yīng)堆堆芯(如小型堆、超臨界水堆、鉛基堆)和新型燃料組件(如MOX燃料)的研發(fā)與設(shè)計(jì)，現(xiàn)行兩群瞬態(tài)中子擴(kuò)散模型正面臨越來(lái)越多的技術(shù)挑戰(zhàn)。伴隨數(shù)值反應(yīng)堆技術(shù)的蓬勃發(fā)展與硬件計(jì)算能力的飛速提升，有必要開(kāi)展瞬態(tài)多群中子輸運(yùn)方程求解方法的研究。

剛性限制法(SCM)利用動(dòng)態(tài)頻率變換技術(shù)，實(shí)現(xiàn)中子通量密度平衡方程與緩發(fā)中子先驅(qū)核密度平衡方程的解耦，將動(dòng)力學(xué)方程中的剛性問(wèn)題淹沒(méi)在穩(wěn)態(tài)方程的求解過(guò)程中[1-2]，確保在較大時(shí)間離散步長(zhǎng)情況下獲得較高的計(jì)算精度，大幅提升了計(jì)算效率。SCM已在兩群瞬態(tài)中子擴(kuò)散方程求解中獲得了廣泛應(yīng)用[3-5]。

本文將SCM拓展應(yīng)用于求解多群瞬態(tài)特征線中子輸運(yùn)方程，在原有中子輸運(yùn)方程特征線方法(MOC)求解程序PEACH[6]的基礎(chǔ)上，增添瞬態(tài)求解功能，開(kāi)發(fā)PEACH-K程序，并通過(guò)數(shù)值模擬驗(yàn)證其精度與穩(wěn)定性。

1 時(shí)空中子動(dòng)力學(xué)模型

1.1 瞬態(tài)中子輸運(yùn)方程的動(dòng)態(tài)頻率變換

時(shí)間相關(guān)的多群中子輸運(yùn)方程和緩發(fā)中子先驅(qū)核密度方程為：

(1)

(2)

考慮到在反應(yīng)堆物理中通常假設(shè)裂變中子各向同性，Cm主要來(lái)自于裂變，因此式(2)中Cm假設(shè)與方向無(wú)關(guān)。

引入以下中子通量動(dòng)態(tài)頻率ωg和緩發(fā)中子先驅(qū)核密度動(dòng)態(tài)頻率μm：

(3)

(4)

(5)

將式(3)和(5)代入式(1)、(2)，再引入ωg(r,Ω,t)各向同性假設(shè)，用ωg(r,t)近似替代ωg(r,Ω,t)，動(dòng)力學(xué)方程變換為如下形式：

(6)

其中：Σ′t,g為動(dòng)態(tài)總截面；χ′g為動(dòng)態(tài)瞬發(fā)中子裂變能譜。

式(6)的中子通量仍是各向異性的。Σ′t,g和χ′g的計(jì)算表達(dá)式如下：

(7)

(8)

式(6)與穩(wěn)態(tài)中子輸運(yùn)方程具有相同的形式。在動(dòng)力學(xué)數(shù)值迭代過(guò)程中，一旦獲得動(dòng)態(tài)頻率ωg和μm后，可采用穩(wěn)態(tài)中子輸運(yùn)方程MOC對(duì)式(6)進(jìn)行求解。需要說(shuō)明的是，式(4)中的φg(r，t)為MOC最小平源近似區(qū)的標(biāo)量中子通量。

1.2 動(dòng)態(tài)頻率的迭代更新

式(6)是一本征值問(wèn)題。與求解穩(wěn)態(tài)中子輸運(yùn)方程類似，可引入一動(dòng)態(tài)本征值kD，這樣式(6)可轉(zhuǎn)換為：

(9)

其中，動(dòng)態(tài)參數(shù)Σ′t,g和χ′g是動(dòng)態(tài)頻率的函數(shù)。迭代收斂時(shí)，kD等于1。利用式(4)和(9)可對(duì)ωg進(jìn)行更新。然而式(9)求解的本征函數(shù)只是中子通量分布而非中子通量幅度。在缺乏中子通量幅度的情況下，動(dòng)態(tài)頻率則無(wú)法有效更新。為克服此問(wèn)題，將動(dòng)態(tài)頻率分解為形狀頻率和幅度頻率兩部分：

ωg(r,t)=ωs,g(r,t)+ωt(t)

(10)

μm(r,t)=μs,m(r,t)+ωt(t)

(11)

其中，形狀頻率ωs，g、μs，m與位置、時(shí)間相關(guān)，而幅度頻率ωt僅與時(shí)間相關(guān)，具備全局屬性。

形狀頻率根據(jù)本征函數(shù)進(jìn)行更新，幅度頻率根據(jù)本征值進(jìn)行更新。在迭代過(guò)程中，如果式(9)中的kD不等于1，那么可在動(dòng)態(tài)頻率的基礎(chǔ)上增加一定的修正量Δωt，使得kD接近于1，此時(shí)要求以下關(guān)系成立：

(12)

令：

(13)

利用一階泰勒展開(kāi)，得到：

γg(ωt)+γ′g(ωt)·Δωt

(14)

其中，γ′g(ωt)=dγg(ωt)/dωt。將式(13)和(14)代入式(12)可得：

(15)

選取1組權(quán)重函數(shù)向量W，采用剩余權(quán)重法對(duì)式(15)進(jìn)行能量和空間積分，可得：

(16)

其中：

(17)

l*=

(18)

至此，中子通量的動(dòng)態(tài)頻率可根據(jù)式(4)、(10)和(16)進(jìn)行更新，緩發(fā)中子先驅(qū)核密度的動(dòng)態(tài)頻率可根據(jù)式(5)和(11)進(jìn)行更新。針對(duì)權(quán)重函數(shù)向量W，在兩群瞬態(tài)中子擴(kuò)散方程求解時(shí)，建議選用共軛通量分布函數(shù)；在多群瞬態(tài)中子輸運(yùn)方程求解時(shí)，由于共軛通量分布函數(shù)需通過(guò)求解共軛方程才能獲得，影響計(jì)算效率，建議選用裂變率分布函數(shù)。

1.3 緩發(fā)中子先驅(qū)核密度的求解

緩發(fā)中子先驅(qū)核密度采用解析法求解。在獲得tn時(shí)刻的先驅(qū)核密度Cm(r,tn)后，將式(2)在tn～tn+1時(shí)間范圍內(nèi)積分，可得tn+1時(shí)刻的先驅(qū)核密度Cm(r,tn+1)：

Cm(r,tn+1)=Cm(r,tn)·e-λmΔtn+e-λmΔtn·

(19)

其中，Δtn=tn+1-tn。

F(r,t)=F(r,tn)+

(20)

根據(jù)式(19)和(20)可解得：

(21)

這樣，通過(guò)數(shù)值迭代求解獲得新的中子通量密度后，即可根據(jù)式(21)有效地更新緩發(fā)中子先驅(qū)核密度。

1.4 數(shù)值迭代流程

在穩(wěn)態(tài)中子輸運(yùn)方程MOC求解程序PEACH的基礎(chǔ)上，根據(jù)上述動(dòng)力學(xué)模型，開(kāi)發(fā)了動(dòng)力學(xué)求解程序PEACH-K，數(shù)值迭代流程如圖1所示。

圖1 PEACH-K程序動(dòng)力學(xué)模型的數(shù)值迭代流程

2 數(shù)值驗(yàn)證

國(guó)際上現(xiàn)有瞬態(tài)基準(zhǔn)問(wèn)題多用于驗(yàn)證基于組件均勻化和少群擴(kuò)散近似的動(dòng)力學(xué)模型，鮮有針對(duì)瞬態(tài)中子輸運(yùn)方程的非均勻多群?jiǎn)栴}。為此，OECD/NEA于2016年發(fā)布了基準(zhǔn)題C5G7-TD[7-8]，含二維和三維，本文基于該基準(zhǔn)題(二維)對(duì)PEACH-K程序進(jìn)行數(shù)值驗(yàn)證。

2.1 C5G7-TD基準(zhǔn)題描述

C5G7-TD基于OECD/NEA于2001年發(fā)布的基準(zhǔn)題C5G7-MOX改造而來(lái)，用于瞬態(tài)非均勻輸運(yùn)模型驗(yàn)證。該問(wèn)題堆芯由UO2燃料組件和MOX燃料組件混合裝載，共計(jì)16盒燃料組件，呈1/8對(duì)稱，具有強(qiáng)泄漏、組件間能譜差異大、非均勻性強(qiáng)等特點(diǎn)。圖2示出該問(wèn)題的1/4堆芯裝載圖，其中每個(gè)組件的幾何尺寸為21.42 cm×21.42 cm，水反射層的寬度也為21.42 cm。該問(wèn)題所使用的UO2和MOX組件均為17×17的元件排列方式，其中UO2組件內(nèi)只有一種富集度的燃料元件棒，而每個(gè)MOX組件內(nèi)則都有3種含量不同的MOX燃料柵元，詳見(jiàn)文獻(xiàn)[7-8]。

圖2 二維C5G7-TD基準(zhǔn)題堆芯布置

圖3 TD0～2下控制棒的移動(dòng)

C5G7-TD基準(zhǔn)題基于不同的反應(yīng)性引入方式，共設(shè)計(jì)了4種算例，針對(duì)每種算例又設(shè)計(jì)了多項(xiàng)分支工況。其中，前3種算例(TD0～2)采用控制棒引入反應(yīng)性，第4種算例(TD3)采用慢化劑密度變化引入反應(yīng)性。針對(duì)控制棒引入反應(yīng)性的方式，圖3示出了TD0～2下控制棒插入份額隨時(shí)間的變化曲線。算例TD0和TD1分別設(shè)計(jì)了5項(xiàng)分支工況，算例TD2設(shè)計(jì)了3項(xiàng)分支工況，對(duì)應(yīng)不同的控制棒組定義(表1)。針對(duì)慢化劑密度變化引入反應(yīng)性的方式，圖4示出了TD3各分支工況堆芯平均慢化劑密度隨時(shí)間的變化曲線，圖4中ω表示慢化劑密度隨時(shí)間變化與初始工況(零時(shí)刻)的比例關(guān)系。

表1 TD0～2分支計(jì)算中控制棒組移動(dòng)定義

圖4 TD3堆芯平均慢化劑密度的變化

2.2 基準(zhǔn)題數(shù)值結(jié)果

C5G7-TD基準(zhǔn)題的參考解目前還在搜集和整理階段，本文采用NECP-X程序[9]作為比較基準(zhǔn)，該程序采用預(yù)測(cè)-修正準(zhǔn)靜態(tài)的動(dòng)力學(xué)求解方法。針對(duì)本基準(zhǔn)題在前2 s反應(yīng)性變化較劇烈的特點(diǎn)，NECP-X程序0～2 s采用的時(shí)間步長(zhǎng)為0.025 s，2～2.5 s采用的時(shí)間步長(zhǎng)為0.05 s，2.5～3 s采用的時(shí)間步長(zhǎng)為0.1 s，3～10 s采用的時(shí)間步長(zhǎng)為0.5 s，而PEACH-K程序在全時(shí)間區(qū)域內(nèi)時(shí)間步長(zhǎng)均采用0.25 s。圖5～8示出各算例堆芯歸一化功率水平的變化?？傮w來(lái)看，PEACH-K與NECP-X程序結(jié)果符合非常好，TD1～3的堆芯歸一化功率水平的相對(duì)偏差均在1%以內(nèi)，表明PEACH-K程序在大時(shí)間步長(zhǎng)下仍具有很高的計(jì)算精度。

針對(duì)TD0，控制棒采用階躍反應(yīng)性引入方式，在個(gè)別拐點(diǎn)處(如TD0-1、TD0-4和TD0-5在2 s附近)由于控制棒價(jià)值太大引起了劇烈的功率突變幅度，PEACH-K與NECP-X程序的堆芯歸一化功率水平的最大相對(duì)偏差在2%以內(nèi)。在工程實(shí)際中，控制棒大幅度階躍提升是不太可能出現(xiàn)的。

a——TD0-1；b——TD0-2；c——TD0-3；d——TD0-4；e——TD0-5圖5 TD0堆芯歸一化功率水平的變化

a——TD1-1；b——TD1-2；c——TD1-3；d——TD1-4；e——TD1-5圖6 TD1堆芯歸一化功率水平的變化

a——TD2-1；b——TD2-2；c——TD2-3圖7 TD2堆芯歸一化功率水平的變化

a——TD3-1；b——TD3-2；c——TD3-3；d——TD3-4圖8 TD3堆芯歸一化功率水平的變化

3 結(jié)論

本文介紹了剛性限制法在瞬態(tài)中子輸運(yùn)計(jì)算中的應(yīng)用。采用OECD/NEA發(fā)布的C5G7-TD基準(zhǔn)題對(duì)采用剛性限制法的PEACH-K程序進(jìn)行了數(shù)值驗(yàn)證，結(jié)果表明PEACH-K程序在大時(shí)間步長(zhǎng)下仍具有很高的計(jì)算精度，剛性限制法可用于非均勻輸運(yùn)瞬態(tài)計(jì)算。

感謝美國(guó)西屋公司退休專家趙榮安博士在剛性限制法研究方面的鼓勵(lì)。感謝西安交通大學(xué)劉宙宇博士提供C5G7-TD基準(zhǔn)題NECP-X程序的對(duì)比結(jié)果。