湯春桃，畢光文，楊 波

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

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

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

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

1 時空中子動力學模型

1.1 瞬態(tài)中子輸運方程的動態(tài)頻率變換

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

(1)

(2)

考慮到在反應堆物理中通常假設裂變中子各向同性，Cm主要來自于裂變，因此式(2)中Cm假設與方向無關。

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

(3)

(4)

(5)

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

(6)

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

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

(7)

(8)

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

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

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

(9)

其中，動態(tài)參數(shù)Σ′t,g和χ′g是動態(tài)頻率的函數(shù)。迭代收斂時，kD等于1。利用式(4)和(9)可對ωg進行更新。然而式(9)求解的本征函數(shù)只是中子通量分布而非中子通量幅度。在缺乏中子通量幅度的情況下，動態(tài)頻率則無法有效更新。為克服此問題，將動態(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與位置、時間相關，而幅度頻率ωt僅與時間相關，具備全局屬性。

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

(12)

令：

(13)

利用一階泰勒展開，得到：

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

(14)

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

(15)

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

(16)

其中：

(17)

l*=

(18)

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

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

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

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

1.4 數(shù)值迭代流程

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

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

2 數(shù)值驗證

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

2.1 C5G7-TD基準題描述

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

圖2 二維C5G7-TD基準題堆芯布置

圖3 TD0～2下控制棒的移動

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

表1 TD0～2分支計算中控制棒組移動定義

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

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

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

針對TD0，控制棒采用階躍反應性引入方式，在個別拐點處(如TD0-1、TD0-4和TD0-5在2 s附近)由于控制棒價值太大引起了劇烈的功率突變幅度，PEACH-K與NECP-X程序的堆芯歸一化功率水平的最大相對偏差在2%以內(nèi)。在工程實際中，控制棒大幅度階躍提升是不太可能出現(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)中子輸運計算中的應用。采用OECD/NEA發(fā)布的C5G7-TD基準題對采用剛性限制法的PEACH-K程序進行了數(shù)值驗證，結(jié)果表明PEACH-K程序在大時間步長下仍具有很高的計算精度，剛性限制法可用于非均勻輸運瞬態(tài)計算。

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