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

（南華大學(xué) 核科學(xué)技術(shù)學(xué)院 衡陽 421001）

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

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

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

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

1 理論模型

1.1 中子擴(kuò)散方程

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

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

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

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

通過式（2）可知中子擴(kuò)散方程是一個廣義本征值問題，經(jīng)過轉(zhuǎn)化可變成標(biāo)準(zhǔn)本征值問題：

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

1.2 隱式重啟Arnoldi方法（IRAM）

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

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

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

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

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

1）進(jìn)行m步Arnoldi迭代，得到AVm=VmHm+；

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

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

①令H1=Hm

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

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

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

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

1.3 POD-Galerkin法構(gòu)建低階模型

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

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

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

利用式（9）的快照矩陣構(gòu)造相關(guān)矩陣S：

式中：S為L階方陣。

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

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

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

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

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

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

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

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

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

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

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

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

2.1 TWIGL基準(zhǔn)題介紹

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

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

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

2.2 利用IRAM獲取快照矩陣

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

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

圖2 IRAM求解1區(qū)吸收截面為0.148時前3階諧波快、熱群分布（彩圖見網(wǎng)絡(luò)版）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組諧波構(gòu)建快照矩陣M，通過式（10）～（15）求得一組POD基，該組POD基對應(yīng)特征值與能量占比如表3所示。

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

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

2.3 中子注量率分布重構(gòu)

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

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

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

圖3 TWIGL基準(zhǔn)題堆芯中子注量率分布 (a) 快群，(b) 熱群（彩圖見網(wǎng)絡(luò)版）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可以看出，降階模型重構(gòu)堆芯中子注量率分布與參考解基本吻合，對角線上快群和熱群中子注量率最大相對誤差為2.56%，在計算時間上，低階模型計算一次的時間為0.039928 s，全階模型計算一次的時間0.392965 s，低階模型計算時間僅為全階模型的10.18%，因此，低階模型不僅可以準(zhǔn)確重構(gòu)堆芯中子注量率，而且可以有效減小計算用時，降低計算資源的消耗。

3 結(jié)語

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

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

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