基于隱式重啟Arnoldi方法的中子擴散本征值問題求解及其降階研究

向釗才 陳洽鋒 趙鵬程 張慶航

（南華大學 核科學技術學院 衡陽 421001）

多群中子擴散方程求解是核工程領域中的一個重要問題，通過求解中子擴散方程，可獲得反應堆內中子分布、功率分布以及臨界狀態(tài)等信息，對研究反應堆物理特性、優(yōu)化反應堆設計和提高反應堆的安全性具有重要意義［1］。然而，實際反應堆工況復雜多變，傳統(tǒng)源迭代方法求解收斂速度慢，計算耗時長，難以在短時間內準確預測反應堆內中子注量率分布，因此有必要發(fā)展高保真、高效率的模型以便于重構反應堆內中子注量率分布。

模型降階技術可以將高維度方程轉化為低維度方程，降低計算復雜度并提高計算效率。20世紀70年代首次出現(xiàn)在工程領域文獻中，經過幾十年的發(fā)展，模型降階已成為提高數(shù)值計算效率的有效方法［2］。隱式重啟Arnoldi方法（Implicitly Restarted Arnoldi Method，IRAM）和本征正交分解方法（Proper Orthogonal Decomposition，POD）都是一種正交投影的模型降階方法。近年來，利用隱式重啟Arnoldi方法求解中子擴散方程的高階諧波已有一定研究。黃義超等［3］基于隱式重啟Arnoldi方法對中子擴散方程高階諧波進行求解，結果顯示IRAM方法在求解速度上相較傳統(tǒng)源修正法提升了10倍，同時保持了相同的精度；王常輝等［4］利用高階諧波的完備性和正交性重構堆芯中子注量率分布，并對反應堆功率進行在線監(jiān)測，監(jiān)測誤差在±3%以內；謝金森等［5］基于重啟Arnoldi方法求解了瞬發(fā)α本征值問題，驗證了IRAM在求解瞬發(fā)α本征值問題的正確性。

POD法最早由Pearson于1901年提出［6］，后由Hotelling于1933年進一步發(fā)展［7］，1987年，Sirovich引入Snapshots方法求解POD基，極大地簡化了數(shù)據(jù)處理過程，促使POD降階方法得以快速發(fā)展［8］。關于POD法在中子擴散方程求解中的應用有了一些進展。Buchan等［9］采用基于時間依賴特征函數(shù)解的實例構造快照矩陣（Snapshot）并求解中子擴散方程的本征值問題；Sartori等［10］比較了模態(tài)法和POD法在求解時空動力學問題中的保真度，結果表明POD法具有更高的保真度；羊俊合等［11］將POD法與Galerkin投影結合求解中子輸運方程，實現(xiàn)了較高的保真度和較高的計算效率，并且可重構堆芯中子注量率分布。然而POD方法對樣本數(shù)據(jù)的模擬參數(shù)較為敏感，當樣本輸入數(shù)據(jù)質量較差時，輸出結果的準確性就會變差，因此需要高精度模型以提供樣本數(shù)據(jù)［12］。

采用隱式重啟Arnoldi方法求解本征值問題的中子擴散方程得到不同宏觀截面狀態(tài)下的高階諧波數(shù)據(jù)構建快照矩陣，基于POD與Galerkin投影相結合的方法構建POD-Galerkin低階模型，并重構穩(wěn)態(tài)TWIGL基準題中子注量率分布進行驗證。該方法將利用高保真、高效率的IRAM生成快照矩陣并與POD-Galerkin降階技術相結合，以提高計算效率和預測準確性，為堆芯中子注量率重構提供一種有效方法。

1 理論模型

1.1 中子擴散方程

中子擴散方程是中子輸運方程擴散近似的結果，描述了單位體積內中子數(shù)守恒過程。多群中子擴散方程組是對中子擴散方程能量變量的離散，即將中子能量分成若干群。非臨界狀態(tài)下穩(wěn)態(tài)表達式如下：

式中：?g是第g群中子注量率，cm-2·s-1；Dg為第g群中子擴散系數(shù)，cm；∑t，g是第g群中子移出截面，cm-1；∑s，g'→g是第g'群中子散射入g群的散射截面，cm-1；∑f，g是第g群中子裂變截面，cm-1；v為相應的裂變產額；χg為中子能量對應于能群g的能譜；keff為有效增值因子。

多群中子擴散方程用矩陣表示如下：

式中：L表示泄漏和移出矩陣；S是散射源矩陣；F是裂變源矩陣；是本征值。

通過式（2）可知中子擴散方程是一個廣義本征值問題，經過轉化可變成標準本征值問題：

式中：k是本征值。對于含增值介質的問題，k有n個正的離散解k1>k2>…>kn>0，對應諧波?1，?2，…，?n，其中n為諧波階數(shù)。最大本征值k1是基波本征值，它就是有效增值系數(shù)keff；?1是k基波，表示系統(tǒng)臨界時中子注量率分布。而當系統(tǒng)偏離臨界較遠時，?1并不表示系統(tǒng)真實的中子注量率分布。

1.2 隱式重啟Arnoldi方法（IRAM）

針對式（3）的特征值問題，Arnoldi方法選取一個初始單位向量x0構造m（m>n）維克雷洛夫（Krylov）子空間：

再通過修正的格拉姆-施密特正交化（Modified Gram-Schmidt）求得該子空間的一組正交基Vm=[v1，v2，…，vm]，其中：

利用這組正交基將矩陣A約化為m維的上海森堡矩陣（Heisenberg）即：

通過求解Hm的本征值k和本征向量y即可得到A的本征值k和本征向量（諧波）?：

Arnoldi方法將大型矩陣本征值問題投影到低維度Krylov子空間，通過求解低維度矩陣的本征值問題反演求得高維度矩陣的本征值和本征向量。然而Arnoldi方法需要儲存Krylov子空間所有基向量，而且由于計算機舍入誤差的積累可能導致基向量正交性的喪失，導致該方法計算效率和穩(wěn)定性降低。隱式重啟Arnoldi方法由Sorensen等提出［13］，該方法將Arnoldi過程與隱式位移迭代QR分解算法相結合，通過周期性地重啟迭代過程來減小Krylov子空間的規(guī)模，并不斷更新初始向量x0，從而保證基向量的正交性并降低計算成本。隱式重啟Arnoldi方法的具體過程如下：

1）進行m步Arnoldi迭代，得到AVm=VmHm+；

2）求得Hm的前q個較小的特征值作為選取位移u1，u2，…，up；

3）QR迭代是一種數(shù)值計算方法，它通過迭代地進行正交三角分解（QR分解），將原始矩陣逐漸變換為對角矩陣或塊對角矩陣的形式。對Hm進行q次隱式QR位移，得到正交矩陣Q，過程如下：

①令H1=Hm

②對i=1，…，q執(zhí)行以下循環(huán)：

③得到最終正交矩陣：Q=Q1Q2…Qp

5）利用重啟向量和截短形式回到第一步，重啟整個過程直到收斂。

由于隱式重啟Arnoldi方法在Arnoldi迭代過程中進行簡單的自動處理，因此開發(fā)相關軟件十分容易，目前已有許多程序包可以使用［14］。

1.3 POD-Galerkin法構建低階模型

POD方法可以提取數(shù)據(jù)集中的主要模態(tài)或主成分，將高維數(shù)據(jù)進行降維處理的同時保留數(shù)據(jù)主要特征，通過POD方法，可以將中子擴散方程的高維度解降低到低維度解，從而提高計算效率和減少計算資源的需求。

通過IRAM可以求得不同宏觀截面下的特征值ki和其對應諧波?i，利用這組諧波便可以構造一個快照（snapshot）矩陣M：

式中：L為快照數(shù)；?1～?L為不同宏觀截面下諧波。

利用式（9）的快照矩陣構造相關矩陣S：

式中：S為L階方陣。

為獲得快照矩陣中數(shù)據(jù)的特征信息，需求得相關矩陣S的特征值和特性向量。相關矩陣的特征值代表了數(shù)據(jù)在特定方向上的重要性，特征值越大，對數(shù)據(jù)貢獻越大；特征向量代表了數(shù)據(jù)在特定方向上的分布情況，其方向表示數(shù)據(jù)在該方向上的主要特征，是一個相互正交向量，可構成POD正交基，用于描述數(shù)據(jù)的主要特征。求解公式如下：

定向越野運動不但考驗學生的身體素質，對其綜合能力也是一種考驗。如，在開展定向越野運動時，首先，教師可多設定幾條不同路線，并在不同路線上設置多個不同的目標，路線的選擇與最終的成績相關聯(lián)，在這個思考、分析和選擇的過程中，學生的綜合能力也會得到相應的鍛煉與提高。其次，教師可在地圖上設置一些陷阱目標，從而鍛煉學生的分析能力、選擇能力和決斷能力，這對學生未來的發(fā)展具有積極意義。

式中：λi和ψi為相關矩陣S的特征值和特征向量；φi為POD基?？梢则炞CPOD基滿足正交性，即：

為提高計算效率，只需選取占比最重的前r個POD正交基，所取的POD基的個數(shù)即為模型的階數(shù)，選取準則如下：

式中：r為POD正交基的個數(shù)，一般γ≥99，W為模態(tài)矩陣。

經過上述過程求得的POD基矩陣用最少的模態(tài)數(shù)描述快照矩陣的主要特征，反應數(shù)據(jù)特征信息的同時，極大提高了計算效率。中子擴散方程諧波可用這組正交基近似表示：

式中：bi為POD基φi(i=1，2，…，r)所對應的系數(shù)。只需求系數(shù)矩陣B便可求得中子擴散方程的諧波。

將式（16）代入式（3），再運用Galerkin法便可構建出中子擴散方程低階模型：

從式（18）可知，POD-Galerkin法將高維度矩陣A本征值問題轉化為低維度矩陣WTAW的本征值問題。通過求解低階系統(tǒng)（式（18））的本征值和本征向量即可反演出全階系統(tǒng)（式（3））的本征值和本征向量。

為驗證低階模型的準確度，需要構建諧波計算精度驗證模型，采用王常輝等［4］提到的表征本征值和諧波計算精度的誤差量：

采用有限差分法離散中子擴散方程，利用Matlab實現(xiàn)隱式重啟Arnoldi方法求解本征值問題中子擴散方程獲得諧波數(shù)據(jù)，并利用所得數(shù)據(jù)編制POD-Galerkin低階模型代碼。

2 二維穩(wěn)態(tài)TWIGL基準題驗證

2.1 TWIGL基準題介紹

TWIGL基準題是一個簡化的矩形反應堆堆芯結構，包含三個區(qū)域，其中1區(qū)和2區(qū)是增值區(qū)，3區(qū)是包裹區(qū)，瞬態(tài)工況下1區(qū)的材料的熱吸收截面發(fā)生變化來驅動反應堆功率變化，幾何布置和邊界條件如圖1所示，各區(qū)域的物性參數(shù)如表1所示。

表1 二維穩(wěn)態(tài)TWIGL基準題物性參數(shù)Table 1 Physical parameters of two-dimensional steadystate TWIGL benchmark problem

圖1 TWIGL基準題區(qū)域劃分、幾何尺寸以及邊界條件Fig.1 TWIGL benchmark problem area division, geometric dimensions, and boundary conditions

2.2 利用IRAM獲取快照矩陣

由于TWIGL基準瞬態(tài)工況中，區(qū)域1中熱中子吸收截∑a，2面發(fā)生變化，且變化范圍在0.146～0.15 cm-1，因此，選擇這一范圍內的狀態(tài)變化構造快照矩陣。利用IRAM方法分別求得∑a，2取值為0.148 cm-1和0.149 cm-1的前6階本征值和前6階諧波。IRAM本征值求解結果及誤差如表2所示，前3階諧波云圖如圖2所示。

表2 不同截面下前6階本征值計算結果及誤差Table 2 Calculation results and relative deviations of the first six eigenvalues for different cross-sections

圖2 IRAM求解1區(qū)吸收截面為0.148時前3階諧波快、熱群分布（彩圖見網絡版）Fig.2 Distribution of the first three harmonics for fast and thermal groups when the absorption cross-section of Zone 1 is 0.148, solved by IRAM (color online)

從表2可知，IRAM求解εi誤差數(shù)量級在10-14，具有很高的精確度，利用IRAM所求得的12組諧波構建快照矩陣M，通過式（10）～（15）求得一組POD基，該組POD基對應特征值與能量占比如表3所示。

表3 POD基對應的特征值及能量占比Table 3 Eigenvalues and energy ratios corresponding to POD basis

從表3可知，前6個POD基總能量占比高達99.99%，且每個POD基能量占比相近，因此選擇前6個POD基構成模態(tài)矩陣W，它們用最少的模態(tài)數(shù)表示了快照矩陣的數(shù)值特征，這極大地減少了計算資源的消耗。

2.3 中子注量率分布重構

通過IRAM求解1區(qū)熱中子吸收截面分別為1.48 cm-1和1.49 cm-1時的12組諧波，利用這組諧波構建快照矩陣并提取POD基，基于POD基的總能量占比構建模態(tài)矩陣W，最后再運用IRAM求解中子擴散低階模型即可重構1區(qū)熱中子吸收截面為1.50 cm-1時的本征值和諧波。所得結果如表4所示。

表4 TWIGL基準題堆芯前6階本征值計算結果Table 4 Calculation results of the first six eigenvalues for the TWIGL benchmark problem core

從表4可知，采用POD-Galerkin方法構建的低階模型在重構TWIGL基準題本征值和諧波方面，其誤差εi隨著本征值階數(shù)而增加，但對于基波本征值和基波的預測該方法仍具有較好的精度。而實際反應堆往往在臨界附近，這時基波本征值即為有效增值系數(shù)keff數(shù)值為0.913127，與參考解0.913214的誤差僅為8.7×10-5?；礊榉磻鸭礊橹凶幼⒘柯史植?，因此可重構堆芯中子注量率分布。重構中子注量率分布云圖如圖3所示，對角線上歸一化中子注量分布與參考解［15］對比如圖4所示。

圖3 TWIGL基準題堆芯中子注量率分布 (a) 快群，(b) 熱群（彩圖見網絡版）Fig.3 TWIGL benchmark problem core neutron flux distribution for fast group (a) and thermal group (b) (color online)

圖4 對角線中子注量率分布與參考解對比 (a) 快群，(b) 熱群Fig.4 Diagonal neutron dose rate distribution and comparison with reference solution for fast group (a) and thermal group (b)

從圖3、4可以看出，降階模型重構堆芯中子注量率分布與參考解基本吻合，對角線上快群和熱群中子注量率最大相對誤差為2.56%，在計算時間上，低階模型計算一次的時間為0.039928 s，全階模型計算一次的時間0.392965 s，低階模型計算時間僅為全階模型的10.18%，因此，低階模型不僅可以準確重構堆芯中子注量率，而且可以有效減小計算用時，降低計算資源的消耗。

3 結語

采用隱式重啟Arnoldi方法求解本征值問題中子擴散方程獲得諧波數(shù)據(jù)，利用所得數(shù)據(jù)構建PODGalerkin低階模型，并重構二維穩(wěn)態(tài)TWIGL基準題中子注量率分布，所得結果如下：隱式重啟Arnoldi方法在求解中子擴散方程的高階本征值和高階諧波時表現(xiàn)出較高的精度；基于POD-Galerkin低階模型重構穩(wěn)態(tài)TWIGL基準題中子注量率分布，在計算精度上，低階模型預測對角線上快群和熱群中子注量率最大相對誤差為2.56%；在計算時間上，低階模型計算時間僅為全階模型的10.18%，因此低階模型能在保持精度的同時有效降低計算資源的消耗。

本研究為堆芯中子注量率求解提供了一種可靠且高效的方法，該方法不僅可用于重構穩(wěn)態(tài)時堆芯中子注量率分布，還具有在瞬態(tài)情況下預測中子注量率的潛力，有望在未來的應用中進一步拓展。

作者貢獻聲明向釗才負責起草文章，編寫代碼，分析/解釋數(shù)據(jù)；陳洽鋒負責采集數(shù)據(jù)，統(tǒng)計分析；趙鵬程負責對文章的知識性內容作批評性審閱，獲取研究經費，指導；張慶航負責做表畫圖。