殷 晗 張 斌 劉曉晶 張滕飛

1（上海交通大學(xué)核科學(xué)與工程學(xué)院 上海 200240）

2（中國(guó)核動(dòng)力研究設(shè)計(jì)院核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室 成都 610213）

在新型反應(yīng)堆概念設(shè)計(jì)中，快中子反應(yīng)堆因具有裂變材料增殖和高放廢物嬗變等優(yōu)勢(shì)成為各國(guó)的研究重點(diǎn)。然而，適用于熱堆的中子物理計(jì)算方法大多基于矩形幾何組件，無(wú)法處理六角形快堆組件；另一方面，快中子更長(zhǎng)的平均自由程導(dǎo)致快堆計(jì)算需要采用高階輸運(yùn)方法以處理快堆的強(qiáng)角度各向異性和強(qiáng)泄漏問(wèn)題，因此，精確高效的中子學(xué)計(jì)算方法成為了快堆概念設(shè)計(jì)和方案優(yōu)化的必備條件［1］。

變 分 節(jié) 塊 法［2］（Variational Nodal Method，VNM）作為一種典型的中子輸運(yùn)算法，最早在VARIANT程序中實(shí)現(xiàn)，并廣泛用于快堆物理計(jì)算［3］。VNM從二階偶宇稱形式的中子輸運(yùn)方程出發(fā)，在問(wèn)題求解域建立包含二階中子輸運(yùn)方程和自然邊界條件的泛函，并利用變分原理對(duì)空間項(xiàng)和角度項(xiàng)進(jìn)行離散，通過(guò)分別求解各個(gè)節(jié)塊內(nèi)的響應(yīng)矩陣方程，最終得到問(wèn)題域內(nèi)的中子通量密度分布。近 年 來(lái) ，VNM 已 成 功 應(yīng) 用 于 ERANOS［4］、VIOLET［5］、PANX［6］等多種中子計(jì)算程序中。

然而，傳統(tǒng)的變分節(jié)塊法在角度離散上通常采用球諧函數(shù)展開(kāi)法（PN）［7］，隨著角度展開(kāi)階數(shù)的增加，總的自由度數(shù)目和響應(yīng)矩陣的規(guī)模將迅速增大，帶來(lái)巨大的內(nèi)存占用和時(shí)間消耗。因此，針對(duì)VNM的加速技術(shù)得到廣泛研究，發(fā)展出了矩陣分離算法［2］、積 分 方 法［8］和 準(zhǔn) 反 射 邊 界 條 件（Quasi-Reflected Interface Condition，QRIC）［9］等一系列加速方法。

在之前的工作中［10］，積分方法被成功應(yīng)用于三維均勻矩形節(jié)塊問(wèn)題中，并被證明可有效提升響應(yīng)矩陣的構(gòu)造效率。然而，其僅能處理矩形節(jié)塊問(wèn)題，而不具備處理六角形節(jié)塊的能力，因此不滿足六角形組件幾何快堆模擬的需要；另外，快堆的高階輸運(yùn)模擬帶來(lái)大量的表面角度自由度，這對(duì)程序的計(jì)算效率提出了嚴(yán)峻挑戰(zhàn)，必須采用有效的加速方法以高效求解矩陣方程。在眾多加速方法中，QRIC方法對(duì)高階角度項(xiàng)采用反射邊界條件，從而可以實(shí)現(xiàn)減小節(jié)塊表面角度項(xiàng)的目的。因此，QRIC方法在求解包含大量表面角度項(xiàng)的矩陣方程上有天然的優(yōu)勢(shì)。目前QRIC方法在求解壓水堆問(wèn)題上有效實(shí)現(xiàn)了加速，但是其處理具有更強(qiáng)角度各向異性的快堆問(wèn)題的可行性還有待進(jìn)一步驗(yàn)證。

基于上述兩點(diǎn)，本文對(duì)積分方法進(jìn)行了進(jìn)一步研究，將其應(yīng)用于六角形組件反應(yīng)堆中，并采用QRIC方法對(duì)求解過(guò)程進(jìn)行加速。通過(guò)對(duì)TAKEDA-4基準(zhǔn)題的計(jì)算，驗(yàn)證了該方法高效求解六角形組件幾何快堆問(wèn)題的可行性。

1 理論模型

通過(guò)變分方法，建立包含二階中子輸運(yùn)方程和自然邊界條件的節(jié)塊泛函［11］：

式中：Fv為單個(gè)節(jié)塊內(nèi)部的泛函；ψ+為偶階中子角通量密度；φ為中子標(biāo)通量密度；q為中子源項(xiàng)；Σt為中子總截面；Σs為宏觀散射截面；Ω為方位角向量。值得注意的是，與原有積分形式變分節(jié)塊法［10］的泛函形式不同，此處引入ω=nu?Ωψ-這一偶宇稱變量，避免了處理周期性邊界條件時(shí)復(fù)雜的表面坐標(biāo)系旋轉(zhuǎn)操作。

節(jié)塊內(nèi)部的偶宇稱中子角通量密度的離散形式［10］為：

式 中 ：f(x，y，z)為 正 交 多 項(xiàng) 式 向 量 ，滿 足(x，y，z)dV=IV，其中IV為對(duì)角矩陣，對(duì)角項(xiàng)為節(jié)塊體積，需要注意的是，對(duì)于六角形反應(yīng)堆，正交多項(xiàng)式構(gòu)造的Gramm-Schimdt過(guò)程是在六角形節(jié)塊上進(jìn)行，而非矩形節(jié)塊［10］。ψ(Ω)—節(jié)塊內(nèi)部角度相關(guān)的中子角通量密度矩向量，滿足關(guān)系

節(jié)塊表面偶宇稱中子角通量密度離散形式為：

式中：Y+γ(Ω)為γ面上的偶宇稱球諧函數(shù)向量；ωγ為γ面上的偶宇稱中子角通量密度矩向量。

將式（2）、（3）代入式（1）中得到泛函的離散形式：

相關(guān)的系數(shù)矩陣的具體表達(dá)式為：

為了提高矩陣方程的求解速度，本文采用了準(zhǔn)反射邊界條件（Quasi-Reflected Interface Condition，QRIC）加速技術(shù)［9］。QRIC方法的思想在于：通過(guò)將反射邊界條件作用于節(jié)塊表面的高階角度項(xiàng)，以減少節(jié)塊表面的角度相關(guān)自由度數(shù)目，將響應(yīng)矩陣的規(guī)模減小至與低階角度近似的水平，從而提高求解速度。

為推導(dǎo)QRIC方法的加速原理，將式（3）重寫(xiě)為：

式中：Y+γl(Ω)為低階偶宇稱球諧函數(shù)向量；Y+γh(Ω)為高階偶宇稱球諧函數(shù)向量；ωl為低階矩向量；ωh為高階矩向量。

對(duì)于高階項(xiàng)，采用反射邊界條件處理，即令高階球諧函數(shù)向量中的奇數(shù)m項(xiàng)為零，只保留偶數(shù)m項(xiàng)，處理后的高階角度項(xiàng)用下標(biāo)e表示。將式（2）、式（5）代入式（1）中，得到離散泛函：

施用鈣肥切莫做了無(wú)用功。近些年，部分蔬菜因缺鈣導(dǎo)致品質(zhì)降低甚至壞死，所以很多菜農(nóng)都加大鈣肥施用量。加上很多經(jīng)銷商也一再宣傳鈣肥要在施基肥時(shí)施足。這讓菜農(nóng)認(rèn)為，施基肥時(shí)應(yīng)大量施用鈣肥。但事實(shí)不是鈣肥施得越多，蔬菜就不會(huì)缺鈣。雖然鈣肥在土壤中和植株體內(nèi)移動(dòng)性差，以基施為主，但蔬菜大量表現(xiàn)出缺鈣癥狀卻不是鈣肥施用量不足所致。而是近年來(lái)，由于土壤酸化、追求高產(chǎn)等因素共同作用下，土壤中鹽離子濃度過(guò)高，影響了蔬菜對(duì)鈣的吸收。因此，問(wèn)題的關(guān)鍵是讓蔬菜充分吸收鈣肥，而不是一味地增加鈣肥施用量。而且過(guò)多施用鈣肥，還會(huì)改變土壤的酸堿性，對(duì)蔬菜生長(zhǎng)非常不利。

令方程（6）關(guān)于ωγl和ωγe的一階變分為零，可得到下面兩個(gè)量在節(jié)塊表面的連續(xù)性條件：

準(zhǔn)反射邊界條件對(duì)應(yīng)的系數(shù)矩陣為：

2 數(shù)值結(jié)果

基于上述理論模型，本研究在中子輸運(yùn)計(jì)算程序VITAS［10］的基礎(chǔ)上開(kāi)發(fā)了具備六角形節(jié)塊處理能力的新版VITAS程序。

新版VITAS在原有程序基礎(chǔ)上增加了六角形節(jié)塊輸運(yùn)模型，使其不僅能處理矩形節(jié)塊，同時(shí)能夠處理六角形節(jié)塊，以滿足六角形幾何快堆模擬的要求；另一方面，新版VITAS在原有廣義矩陣分離（Generalized Partitioned Matrix）加速算法［10］的基礎(chǔ)上引入了QRIC加速算法，減少了表面角度項(xiàng)的數(shù)目，從而提高了程序的求解效率。

本文采用TAKEDA-4六角形組件基準(zhǔn)題［12］對(duì)新版VITAS程序進(jìn)行了驗(yàn)證。該基準(zhǔn)題模型是基于KNK-II反應(yīng)堆堆芯建立，其堆芯徑向分布和軸向布置分別如圖1、2所示。整個(gè)堆芯由8圈六角形組件構(gòu)成，堆芯總高度為190 cm（包含軸向反射層），每個(gè)六角形組件的邊長(zhǎng)為7.5 cm。計(jì)算考慮三種不同的控制棒棒位情況：控制棒拔出（CR Out）；控制棒半插（CR Half）；控制棒全插（CR In）。

圖1 TAKEDA-4基準(zhǔn)題的組件徑向分布Fig.1 The assembly radial layout of the TAKEDA-4 benchmark

圖2 TAKEDA-4基準(zhǔn)題的組件軸向分布Fig.2 The assembly axial layout of the TAKEDA-4 benchmark

通過(guò)對(duì)空間展開(kāi)階數(shù)的敏感性分析，最終確定中子角通量密度的空間離散在節(jié)塊內(nèi)部采用7階正交多項(xiàng)式，在節(jié)塊表面采用2階正交多項(xiàng)式，繼續(xù)增大階數(shù)，對(duì)計(jì)算精度不會(huì)有提升。

表1給出了標(biāo)準(zhǔn)PN法和積分方法在不同控制棒插入情況、不同PN階數(shù)的特征值計(jì)算偏差，并由此繪制特征值計(jì)算偏差隨PN階數(shù)的變化圖，如圖3所示，參考解采用Monte Carlo方法計(jì)算得到［12］。計(jì)算結(jié)果表明：1）兩種方法的特征值計(jì)算偏差隨PN階數(shù)的增大呈現(xiàn)相反方向的漸進(jìn)收斂勢(shì)，這主要是由于二者在節(jié)塊內(nèi)部的角度處理方式不同導(dǎo)致的；2）低階角度近似情況下，積分方法具有更高的計(jì)算精度。P1近似下，積分方法可將特征值的計(jì)算偏差降低5倍左右，P3近似下，可將特征值的計(jì)算偏差降低兩倍左右；3）相比標(biāo)準(zhǔn)PN法，積分方法隨著角度階數(shù)的增加可以更快地達(dá)到漸進(jìn)收斂。研究表明：不同控制棒棒位下的計(jì)算精度和計(jì)算效率呈現(xiàn)相同的趨勢(shì)，因此為簡(jiǎn)潔起見(jiàn)，后文中的分析僅以控制棒拔出情況為例。

圖3 特征值計(jì)算偏差隨角度階數(shù)的變化Fig.3 Variation of calculation errors of eigenvalue with angular order

表1 標(biāo)準(zhǔn)PN法和積分方法的特征值計(jì)算偏差對(duì)比Table 1 Eigenvalue error comparisons of the standard method and the integral method

TAKEDA-4控制棒拔出情況下積分方法和標(biāo)準(zhǔn)PN方法的CPU時(shí)間隨角度展開(kāi)階數(shù)的變化情況如圖4所示，其中Total、RM、Solve分別表示總時(shí)間、響應(yīng)矩陣構(gòu)造時(shí)間和矩陣方程求解時(shí)間。此處將標(biāo)準(zhǔn)P7計(jì)算的CPU時(shí)間作為基準(zhǔn)值并歸一。圖4表明：1）標(biāo)準(zhǔn)PN法在階數(shù)增大時(shí)要付出更大的時(shí)間代價(jià)，當(dāng)PN階數(shù)由P5增加到P7，積分方法的計(jì)算時(shí)間由0.04增加到0.20，而標(biāo)準(zhǔn)PN法的總時(shí)間由0.12增加到1.00；2）標(biāo)準(zhǔn)PN法的求解時(shí)間更短，這得益于其較少的表面角度自由度數(shù)目。然而，這并不足以抵消其矩陣構(gòu)造過(guò)程中巨大的時(shí)間開(kāi)銷，因此從總CPU時(shí)間來(lái)看，積分方法的計(jì)算效率更高，在高階角度近似如P7近似下，積分方法相比標(biāo)準(zhǔn)PN法具有5.0的加速比。為驗(yàn)證準(zhǔn)反射邊界條件方法的加速效果，對(duì)積分方法表面全階數(shù)PN處理（‘Full’表示）和表面準(zhǔn)反射邊界條件處理（‘QRIC’表示）在不同PN階數(shù)下的計(jì)算時(shí)間進(jìn)行了比較，如圖5所示。圖中QRIC-P3處理是指把角度項(xiàng)保留至P3，高于P3的角度項(xiàng)采用反射邊界條件處理。另外，為便于比較，此處將表面全階數(shù)P9處理的構(gòu)造時(shí)間作為基準(zhǔn)值并歸一，圖5結(jié)果表明，在P9近似下，QRIC加速算法能夠提升75%的構(gòu)造速度，75%的求解速度，總加速比為6.7。

同時(shí)，本研究針對(duì)QRIC加速算法對(duì)計(jì)算精度的影響進(jìn)行了分析，包括特征值的計(jì)算和中子通量密度的計(jì)算，計(jì)算結(jié)果見(jiàn)圖5和表3。計(jì)算結(jié)果顯示：采用QRIC加速算法所造成的特征值計(jì)算精度的損失小于0.000 12，中子通量密度計(jì)算精度的損失小于0.1%。因此，認(rèn)為采用QRIC加速算法對(duì)計(jì)算精度的影響可以忽略。

表2 控制棒拔出情況下積分方法和標(biāo)準(zhǔn)法的中子通量密度計(jì)算偏差對(duì)比Table 2 Comparisons of flux errors between the integral method and the standard method for the CR out case

圖4TAKEDA-4基準(zhǔn)題控制棒拔出情況下CPU時(shí)間隨角度展開(kāi)階數(shù)的變化Fig.4 Variation of CPU time with degree of angular orders for the TAKEDA-4 benchmark CR Out case

圖5 全角度處理和準(zhǔn)反射邊界條件處理效率-精度對(duì)比Fig.5 Cost-accuracy trade-off comparisons between the full angular treatment and the quasi-reflected interface condition

表3 控制棒全插情況下全角度處理和準(zhǔn)反射邊界條件處理的中子通量密度百分偏差對(duì)比Table 3 Comparisons of neutron flux percentage errors between the full angular treatment and the quasi-reflected interface condition for the CR in case

3 結(jié)語(yǔ)

本文對(duì)積分形式的變分節(jié)塊法［10］進(jìn)行了進(jìn)一步研究和改進(jìn)，提出一種可用于六角形組件幾何快堆模擬的積分輸運(yùn)變分節(jié)塊法。在中子角通量密度的角度處理上，采用積分方法處理節(jié)塊內(nèi)部的中子角通量密度，采用偶宇稱球諧函數(shù)處理節(jié)塊表面的中子角通量密度，同時(shí)采用QRIC加速算法以提高計(jì)算效率。針對(duì)TAKEDA-4基準(zhǔn)題的驗(yàn)證結(jié)果表明：

1）在低階角度近似下，相比于傳統(tǒng)基于球諧函數(shù)離散的變分節(jié)塊法，本方法所采用的積分方法可提高特征值和中子通量密度的計(jì)算精度，使特征值的計(jì)算偏差降低2～5倍，中子通量密度計(jì)算的百分偏差減少約2%。

2）在高階角度近似下，積分方法和QRIC加速算法的共同作用能夠顯著提高計(jì)算效率。在P7近似下，積分方法和QRIC_P3加速算法分別具有5.0和6.6的加速比，相比于傳統(tǒng)變分節(jié)塊法，本方法可實(shí)現(xiàn)33倍的加速比。