彭良輝，湯春桃，畢光文，張宏博，楊 波

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

三維全堆芯pin-by-pin輸運計算是下一代堆芯物理分析方法(NGM)的重要研究內(nèi)容。相較于傳統(tǒng)的組件-堆芯“兩步法”，該方法可獲得更加精確的燃料棒層級的通量及功率分布。三階簡化球諧函數(shù)(SP3)方法是實現(xiàn)三維全堆芯pin-by-pin輸運計算的重要途徑之一。SP3方法的精度高于擴散方法，其計算量約為擴散方法的2倍，且遠(yuǎn)小于其他輸運算法，如特征線方法(MOC)、離散縱標(biāo)方法(SN)等[1]。選取SP3方法作為實現(xiàn)三維全堆芯pin-by-pin輸運計算的途徑，可較好地平衡計算精度與計算量之間的矛盾。此外，SP3方法的控制方程為兩個擴散類型的方程，可沿用成熟的擴散節(jié)塊方法對其進(jìn)行求解。

本文采用解析基函數(shù)展開節(jié)塊(AFEN)方法[2-3]求解SP3方法的控制方程，并對3種偏流表達(dá)形式的計算精度進(jìn)行檢驗。

1 理論模型

1.1 SP3方法的控制方程及邊界條件

多群SP3方法的控制方程[4]如下：

(1)

(2)

以上方程中：

(3)

Σr0=Σt-Σs,Σr2=Σt+0.8Σr0

(4)

(5)

(6)

基于Marshak邊界條件，可得u方向上各階偏流表達(dá)形式[5]為：

(7)

(8)

為了便于采用原有的擴散節(jié)塊展開方法(NEM)求解SP3方法的控制方程，美國普度大學(xué)研發(fā)了PARCS程序?qū)ζ鞅磉_(dá)形式進(jìn)行截斷[6]，即有：

(9)

(10)

此外，在日本NFI研發(fā)的SCOPE2程序中，為簡化節(jié)塊內(nèi)偏流及平均通量響應(yīng)矩陣，便于程序并行化，采用了與菲克定律類似的偏流表達(dá)形式[1，7]：

(11)

對于中子能群g，式(7)～(11)中的偏流表達(dá)形式可統(tǒng)一改寫為如下形式：

(12)

表1 偏流表達(dá)形式的系數(shù)矩陣

獲得偏流表達(dá)形式后，在兩相鄰節(jié)塊k和節(jié)塊k+1的交界面處，滿足如下的偏流連續(xù)條件：

(13)

(14)

1.2 通量近似解

多群SP3方法的控制方程可改寫為如下矩陣形式[3]：

(15)

設(shè)(λm,um),m=1,…,2G為矩陣[A′]的特征值與特征向量。

um=[u1,m…ui,m…u2G,m]T

m=1,…,2G

(16)

則式(15)存在如下通量解析解：

φ=[U]ψ

(17)

其中，[U]=[u1…um…u2G]。ψ中的第m分量為：

m=1,…,2G

(18)

其中：el為任意單位向量；r為位置向量；Aml和Bml為待定系數(shù)。

圖1示出單位特征向量及定解條件示意圖，其中S表示面。如圖1所示，選擇如下的3個單位特征向量，對式(18)中的輔助函數(shù)Ψm進(jìn)行3階截斷(即取l=1,2,3)：

e1=ex,e2=ey,e3=ez

(19)

圖1 單位特征向量及定解條件示意圖

此時，輔助函數(shù)Ψm的近似解表達(dá)式為：

Ψm=em·cm

(20)

其中：

(21)

cm=[Am,1Am,2Am,3Bm,1Bm,2Bm,3]T

m=1,…,2G

(22)

將式(20)代入式(17)有：

φ=[U]Ψ=[U][E]C

(23)

其中：

(24)

(25)

式(23)即為通量近似表達(dá)式。其中：矩陣[U]為式(15)系數(shù)矩陣[A′]對應(yīng)的特征向量矩陣，矩陣維數(shù)為2G×2G；矩陣[E]的矩陣元為含空間位置變量的指數(shù)函數(shù)，矩陣維數(shù)為2G×12G；C為待定系數(shù)向量，維數(shù)為12G，即每個節(jié)塊中共有12G個待定系數(shù)。

1.3 節(jié)塊響應(yīng)矩陣

在某一節(jié)塊中，為求出所有通量待定系數(shù)，需提供12G個定解條件。如圖1所示，可選取節(jié)塊6個邊界處的各階出射或入射偏流為定解條件，反解待定系數(shù)。為此首先需建立偏流與待定系數(shù)之間的關(guān)系式。令：

(26)

結(jié)合式(23)和(26)有：

(27)

(28)

其中：

(29)

其中，3種偏流表達(dá)形式下的系數(shù)矩陣[Mn]如表1所列。

jout=[P]C

(30)

同理，以節(jié)塊6個邊界處的入射偏流為定解條件，可得：

jin=[Q]C

(31)

其中：

(32)

(33)

由式(30)和(31)可得節(jié)塊內(nèi)偏流的響應(yīng)矩陣為：

jout=[P][Q]-1jin

(34)

對式(23)進(jìn)行體積平均，即對[E]的矩陣元進(jìn)行體積平均，可得：

φave=[T]C

(35)

結(jié)合式(31)可得節(jié)塊內(nèi)平均通量的響應(yīng)矩陣為：

φave=[T][Q]-1jin

(36)

式(34)及(36)即為節(jié)塊出射偏流、節(jié)塊平均通量相對于節(jié)塊入射偏流的響應(yīng)關(guān)系。

1.4 源迭代方法

基于節(jié)塊中出射偏流、平均通量與入射偏流間的響應(yīng)關(guān)系，可構(gòu)造如下的源迭代流程：1) 初始化系統(tǒng)keff和各節(jié)塊邊界處的入射偏流jin；2) 計算或更新各節(jié)塊中的響應(yīng)矩陣[P][Q]-1、[T][Q]-1；3) 根據(jù)式(34)，由入射偏流jin計算出射偏流jout；4) 根據(jù)式(36)，由入射偏流jin計算節(jié)塊平均通量φave；5) 根據(jù)節(jié)塊交界面連續(xù)條件和外邊界條件，由出射偏流jout更新入射偏流jin；6) 由節(jié)塊平均通量φave更新keff；7) 判斷keff和節(jié)塊平均通量φave是否收斂，若未收斂則返回第2步，若收斂則結(jié)束。

此外，為提高計算效率，可采用粗網(wǎng)有限差分加速方法對以上源迭代過程進(jìn)行收斂加速[4]。

2 數(shù)值檢驗

采用式(7)～(11)對應(yīng)的3種偏流表達(dá)形式，開發(fā)了用于pin-by-pin輸運計算的SP3解析基函數(shù)展開節(jié)塊程序AFEN_SP3，對該程序進(jìn)行數(shù)值檢驗。

2.1 自定義一維基準(zhǔn)例題

該例題由C5G7-MOX-2D基準(zhǔn)例題[8]演化得到，其材料及幾何布置如圖2所示，從左到右依次為H2O、UO2、4.3%MOX、7.0%MOX和8.7%MOX。各材料區(qū)的寬度均為20 cm，左邊界為真空邊界，右邊界為反射邊界。各材料的7群反應(yīng)截面均根據(jù)C5G7-MOX-2D基準(zhǔn)例題由特征線輸運程序PEACH[9]經(jīng)通量-體積均勻化得到。

圖2 自定義一維基準(zhǔn)例題示意圖

該例題的參考解由MCMG程序計算得到[4]。AFEN_SP3程序計算時，網(wǎng)格寬度為1 cm，計算結(jié)果列于表2。其中擴散程序為半解析擴散節(jié)塊程序Stella[4]，網(wǎng)格劃分方式與AFEN_SP3程序保持一致。特征值keff計算偏差為絕對偏差，柵元功率計算偏差為相對偏差，以下算例相同?？梢姡瑢τ谠撘痪S例題，3種偏流表達(dá)形式下的SP3程序計算精度基本相同，且均優(yōu)于擴散程序。

表2 自定義一維基準(zhǔn)例題計算結(jié)果

2.2 自定義二維基準(zhǔn)例題

該例題由C5G7-MOX-2D基準(zhǔn)例題[8]演化得到，其材料及幾何布置如圖3所示，且?guī)缀纬叽缂安贾梅绞脚cC5G7-MOX-2D基準(zhǔn)例題相同。各材料的7群反應(yīng)截面均根據(jù)C5G7-MOX-2D基準(zhǔn)例題[8]由特征線輸運程序PEACH[9]經(jīng)通量-體積均勻化得到。

該例題的參考解由特征線輸運程序PEACH[9]計算得到。AFEN_SP3程序計算時，網(wǎng)格大小為1.26 cm×1.26 cm，計算結(jié)果列于表3。可見，對于該二維例題，3種偏流表達(dá)形式下的SP3程序計算精度基本相同，且均優(yōu)于擴散程序。

圖3 自定義二維基準(zhǔn)例題示意圖

表3 自定義二維基準(zhǔn)例題計算結(jié)果

注：1) 各柵元功率相對偏差絕對值的平均值

2.3 三維KUCA基準(zhǔn)例題

該例題為日本京都大學(xué)建立的KUCA(Kyoto University Critical Assembly)三維輕水堆基準(zhǔn)例題，該基準(zhǔn)例題的幾何結(jié)構(gòu)如圖4所示，包含燃料區(qū)、反射層和控制棒，各區(qū)的截面參數(shù)列于表4[10]。

該例題的參考解由TEPFEM程序中的P3模型計算得到[11]。AFEN_SP3程序計算時，網(wǎng)格大小為1 cm×1 cm×1 cm，計算結(jié)果列于表5。由表5可見：SP3程序的計算精度明顯優(yōu)于擴散程序；3種偏流表達(dá)形式中，Marshak偏流表達(dá)形式計算精度最高，其余兩種偏流表達(dá)形式的計算精度在可接受范圍內(nèi)。

如圖4，取以線段L(x∈[0 cm，25 cm]，y=2.5 cm，z=0.5 cm)為中心的1排柵元，其裂變源歸一化快群平均通量分布如圖5所示。圖5中，以Marshak偏流下通量計算結(jié)果為基準(zhǔn)，進(jìn)行相對計算偏差統(tǒng)計，用于比較各偏流形式下計算結(jié)果的差別。

圖4 三維KUCA基準(zhǔn)例題示意圖

表5 三維KUCA基準(zhǔn)例題計算結(jié)果

圖5 三維KUCA基準(zhǔn)例題典型區(qū)域快群通量分布

可見，堆芯內(nèi)部(包括燃料區(qū)和控制棒區(qū))各偏流表達(dá)形式下的計算結(jié)果差別小于1%，靠近堆芯外部計算相對偏差變大，最外圍節(jié)塊相對偏差最大，約為4%。結(jié)合式(13)及各偏流表達(dá)形式可見，內(nèi)邊界偏流連續(xù)條件均可由交界面處通量及通量一階偏導(dǎo)數(shù)連續(xù)得到。根據(jù)式(14)，3種偏流表達(dá)形式間的差異主要體現(xiàn)在外邊界真空條件之上。因此，當(dāng)堆芯外圍反射層區(qū)域較大時，如上述的一維及二維基準(zhǔn)例題，3種偏流表達(dá)形式下的節(jié)塊相對功率分布計算結(jié)果差別較小。

3 結(jié)論

基于SP3解析基函數(shù)展開節(jié)塊方法，在3種不同偏流表達(dá)形式下開發(fā)了pin-by-pin輸運程序AFEN_SP3。數(shù)值檢驗結(jié)果表明：AFEN_SP3程序的開發(fā)是正確的；3種偏流表達(dá)形式下的SP3程序計算精度均優(yōu)于擴散程序；3種偏流表達(dá)形式中，Marshak偏流表達(dá)形式計算精度最高，其余兩種偏流表達(dá)形式計算精度稍差，但仍在可接受范圍內(nèi)。簡化后的偏流表達(dá)形式進(jìn)一步簡化了SP3方法，可采用成熟的擴散節(jié)塊程序稍加改動進(jìn)行求解。簡化后的偏流表達(dá)形式有利于進(jìn)一步提高SP3全堆pin-by-pin輸運程序的計算效率，有利于程序并行化算法的設(shè)計。