剛性限制法在瞬態中子輸運計算中的應用

湯春桃，畢光文，楊 波

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

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

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

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

1 時空中子動力學模型

1.1 瞬態中子輸運方程的動態頻率變換

時間相關的多群中子輸運方程和緩發中子先驅核密度方程為：

(1)

(2)

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

引入以下中子通量動態頻率ωg和緩發中子先驅核密度動態頻率μm：

(3)

(4)

(5)

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

(6)

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

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

(7)

(8)

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

1.2 動態頻率的迭代更新

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

(9)

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

ω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僅與時間相關，具備全局屬性。

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

(12)

令：

(13)

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

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

(14)

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

(15)

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

(16)

其中：

(17)

l*=

(18)

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

1.3 緩發中子先驅核密度的求解

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

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

(21)

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

1.4 數值迭代流程

在穩態中子輸運方程MOC求解程序PEACH的基礎上，根據上述動力學模型，開發了動力學求解程序PEACH-K，數值迭代流程如圖1所示。

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

2 數值驗證

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

2.1 C5G7-TD基準題描述

C5G7-TD基于OECD/NEA于2001年發布的基準題C5G7-MOX改造而來，用于瞬態非均勻輸運模型驗證。該問題堆芯由UO2燃料組件和MOX燃料組件混合裝載，共計16盒燃料組件，呈1/8對稱，具有強泄漏、組件間能譜差異大、非均勻性強等特點。圖2示出該問題的1/4堆芯裝載圖，其中每個組件的幾何尺寸為21.42 cm×21.42 cm，水反射層的寬度也為21.42 cm。該問題所使用的UO2和MOX組件均為17×17的元件排列方式，其中UO2組件內只有一種富集度的燃料元件棒，而每個MOX組件內則都有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 基準題數值結果

C5G7-TD基準題的參考解目前還在搜集和整理階段，本文采用NECP-X程序[9]作為比較基準，該程序采用預測-修正準靜態的動力學求解方法。針對本基準題在前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程序在全時間區域內時間步長均采用0.25 s。圖5～8示出各算例堆芯歸一化功率水平的變化?？傮w來看，PEACH-K與NECP-X程序結果符合非常好，TD1～3的堆芯歸一化功率水平的相對偏差均在1%以內，表明PEACH-K程序在大時間步長下仍具有很高的計算精度。

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

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 結論

本文介紹了剛性限制法在瞬態中子輸運計算中的應用。采用OECD/NEA發布的C5G7-TD基準題對采用剛性限制法的PEACH-K程序進行了數值驗證，結果表明PEACH-K程序在大時間步長下仍具有很高的計算精度，剛性限制法可用于非均勻輸運瞬態計算。

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