彭良輝,湯春桃,畢光文,張宏博,楊 波
(上海核工程研究設計院有限公司,上海 200233)
三維全堆芯pin-by-pin輸運計算是下一代堆芯物理分析方法(NGM)的重要研究內容。相較于傳統(tǒng)的組件-堆芯“兩步法”,該方法可獲得更加精確的燃料棒層級的通量及功率分布。三階簡化球諧函數(SP3)方法是實現三維全堆芯pin-by-pin輸運計算的重要途徑之一。SP3方法的精度高于擴散方法,其計算量約為擴散方法的2倍,且遠小于其他輸運算法,如特征線方法(MOC)、離散縱標方法(SN)等[1]。選取SP3方法作為實現三維全堆芯pin-by-pin輸運計算的途徑,可較好地平衡計算精度與計算量之間的矛盾。此外,SP3方法的控制方程為兩個擴散類型的方程,可沿用成熟的擴散節(jié)塊方法對其進行求解。
本文采用解析基函數展開節(jié)塊(AFEN)方法[2-3]求解SP3方法的控制方程,并對3種偏流表達形式的計算精度進行檢驗。
多群SP3方法的控制方程[4]如下:
(1)
(2)
以上方程中:
(3)
Σr0=Σt-Σs,Σr2=Σt+0.8Σr0
(4)
(5)
(6)
基于Marshak邊界條件,可得u方向上各階偏流表達形式[5]為:
(7)
(8)
為了便于采用原有的擴散節(jié)塊展開方法(NEM)求解SP3方法的控制方程,美國普度大學研發(fā)了PARCS程序對偏流表達形式進行截斷[6],即有:
(9)
(10)
此外,在日本NFI研發(fā)的SCOPE2程序中,為簡化節(jié)塊內偏流及平均通量響應矩陣,便于程序并行化,采用了與菲克定律類似的偏流表達形式[1,7]:
(11)
對于中子能群g,式(7)~(11)中的偏流表達形式可統(tǒng)一改寫為如下形式:
(12)
表1 偏流表達形式的系數矩陣
獲得偏流表達形式后,在兩相鄰節(jié)塊k和節(jié)塊k+1的交界面處,滿足如下的偏流連續(xù)條件:
(13)
(14)
多群SP3方法的控制方程可改寫為如下矩陣形式[3]:
(15)
設(λ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為待定系數。
圖1示出單位特征向量及定解條件示意圖,其中S表示面。如圖1所示,選擇如下的3個單位特征向量,對式(18)中的輔助函數Ψm進行3階截斷(即取l=1,2,3):
e1=ex,e2=ey,e3=ez
(19)
圖1 單位特征向量及定解條件示意圖
此時,輔助函數Ψm的近似解表達式為:
Ψ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)即為通量近似表達式。其中:矩陣[U]為式(15)系數矩陣[A′]對應的特征向量矩陣,矩陣維數為2G×2G;矩陣[E]的矩陣元為含空間位置變量的指數函數,矩陣維數為2G×12G;C為待定系數向量,維數為12G,即每個節(jié)塊中共有12G個待定系數。
在某一節(jié)塊中,為求出所有通量待定系數,需提供12G個定解條件。如圖1所示,可選取節(jié)塊6個邊界處的各階出射或入射偏流為定解條件,反解待定系數。為此首先需建立偏流與待定系數之間的關系式。令:
(26)
結合式(23)和(26)有:
(27)
(28)
其中:
(29)
其中,3種偏流表達形式下的系數矩陣[Mn]如表1所列。
jout=[P]C
(30)
同理,以節(jié)塊6個邊界處的入射偏流為定解條件,可得:
jin=[Q]C
(31)
其中:
(32)
(33)
由式(30)和(31)可得節(jié)塊內偏流的響應矩陣為:
jout=[P][Q]-1jin
(34)
對式(23)進行體積平均,即對[E]的矩陣元進行體積平均,可得:
φave=[T]C
(35)
結合式(31)可得節(jié)塊內平均通量的響應矩陣為:
φave=[T][Q]-1jin
(36)
式(34)及(36)即為節(jié)塊出射偏流、節(jié)塊平均通量相對于節(jié)塊入射偏流的響應關系。
基于節(jié)塊中出射偏流、平均通量與入射偏流間的響應關系,可構造如下的源迭代流程:1) 初始化系統(tǒng)keff和各節(jié)塊邊界處的入射偏流jin;2) 計算或更新各節(jié)塊中的響應矩陣[P][Q]-1、[T][Q]-1;3) 根據式(34),由入射偏流jin計算出射偏流jout;4) 根據式(36),由入射偏流jin計算節(jié)塊平均通量φave;5) 根據節(jié)塊交界面連續(xù)條件和外邊界條件,由出射偏流jout更新入射偏流jin;6) 由節(jié)塊平均通量φave更新keff;7) 判斷keff和節(jié)塊平均通量φave是否收斂,若未收斂則返回第2步,若收斂則結束。
此外,為提高計算效率,可采用粗網有限差分加速方法對以上源迭代過程進行收斂加速[4]。
采用式(7)~(11)對應的3種偏流表達形式,開發(fā)了用于pin-by-pin輸運計算的SP3解析基函數展開節(jié)塊程序AFEN_SP3,對該程序進行數值檢驗。
該例題由C5G7-MOX-2D基準例題[8]演化得到,其材料及幾何布置如圖2所示,從左到右依次為H2O、UO2、4.3%MOX、7.0%MOX和8.7%MOX。各材料區(qū)的寬度均為20 cm,左邊界為真空邊界,右邊界為反射邊界。各材料的7群反應截面均根據C5G7-MOX-2D基準例題由特征線輸運程序PEACH[9]經通量-體積均勻化得到。
圖2 自定義一維基準例題示意圖
該例題的參考解由MCMG程序計算得到[4]。AFEN_SP3程序計算時,網格寬度為1 cm,計算結果列于表2。其中擴散程序為半解析擴散節(jié)塊程序Stella[4],網格劃分方式與AFEN_SP3程序保持一致。特征值keff計算偏差為絕對偏差,柵元功率計算偏差為相對偏差,以下算例相同。可見,對于該一維例題,3種偏流表達形式下的SP3程序計算精度基本相同,且均優(yōu)于擴散程序。
表2 自定義一維基準例題計算結果
該例題由C5G7-MOX-2D基準例題[8]演化得到,其材料及幾何布置如圖3所示,且?guī)缀纬叽缂安贾梅绞脚cC5G7-MOX-2D基準例題相同。各材料的7群反應截面均根據C5G7-MOX-2D基準例題[8]由特征線輸運程序PEACH[9]經通量-體積均勻化得到。
該例題的參考解由特征線輸運程序PEACH[9]計算得到。AFEN_SP3程序計算時,網格大小為1.26 cm×1.26 cm,計算結果列于表3??梢?,對于該二維例題,3種偏流表達形式下的SP3程序計算精度基本相同,且均優(yōu)于擴散程序。
圖3 自定義二維基準例題示意圖
表3 自定義二維基準例題計算結果
注:1) 各柵元功率相對偏差絕對值的平均值
該例題為日本京都大學建立的KUCA(Kyoto University Critical Assembly)三維輕水堆基準例題,該基準例題的幾何結構如圖4所示,包含燃料區(qū)、反射層和控制棒,各區(qū)的截面參數列于表4[10]。
該例題的參考解由TEPFEM程序中的P3模型計算得到[11]。AFEN_SP3程序計算時,網格大小為1 cm×1 cm×1 cm,計算結果列于表5。由表5可見:SP3程序的計算精度明顯優(yōu)于擴散程序;3種偏流表達形式中,Marshak偏流表達形式計算精度最高,其余兩種偏流表達形式的計算精度在可接受范圍內。
如圖4,取以線段L(x∈[0 cm,25 cm],y=2.5 cm,z=0.5 cm)為中心的1排柵元,其裂變源歸一化快群平均通量分布如圖5所示。圖5中,以Marshak偏流下通量計算結果為基準,進行相對計算偏差統(tǒng)計,用于比較各偏流形式下計算結果的差別。
圖4 三維KUCA基準例題示意圖
表5 三維KUCA基準例題計算結果
圖5 三維KUCA基準例題典型區(qū)域快群通量分布
可見,堆芯內部(包括燃料區(qū)和控制棒區(qū))各偏流表達形式下的計算結果差別小于1%,靠近堆芯外部計算相對偏差變大,最外圍節(jié)塊相對偏差最大,約為4%。結合式(13)及各偏流表達形式可見,內邊界偏流連續(xù)條件均可由交界面處通量及通量一階偏導數連續(xù)得到。根據式(14),3種偏流表達形式間的差異主要體現在外邊界真空條件之上。因此,當堆芯外圍反射層區(qū)域較大時,如上述的一維及二維基準例題,3種偏流表達形式下的節(jié)塊相對功率分布計算結果差別較小。
基于SP3解析基函數展開節(jié)塊方法,在3種不同偏流表達形式下開發(fā)了pin-by-pin輸運程序AFEN_SP3。數值檢驗結果表明:AFEN_SP3程序的開發(fā)是正確的;3種偏流表達形式下的SP3程序計算精度均優(yōu)于擴散程序;3種偏流表達形式中,Marshak偏流表達形式計算精度最高,其余兩種偏流表達形式計算精度稍差,但仍在可接受范圍內。簡化后的偏流表達形式進一步簡化了SP3方法,可采用成熟的擴散節(jié)塊程序稍加改動進行求解。簡化后的偏流表達形式有利于進一步提高SP3全堆pin-by-pin輸運程序的計算效率,有利于程序并行化算法的設計。