張彬航,畢彥釗,龔瀚源,張永紅,*,袁顯寶

(1.三峽大學 機械與動力學院,湖北 宜昌 443002;2.三峽大學 理學院,湖北 宜昌 443002)

數(shù)值反應堆技術(shù)基于高性能計算和多物理耦合實現(xiàn)對堆芯各種物理現(xiàn)象的精細刻畫,能夠準確預測反應堆在服役過程中的關(guān)鍵安全參數(shù),為先進核能系統(tǒng)的研究開發(fā)和安全分析提供了重要手段和依據(jù)[1-3]。反應堆燃耗計算是實現(xiàn)數(shù)值反應堆精細模擬的關(guān)鍵環(huán)節(jié)。燃耗計算通過求解燃耗方程獲得燃料同位素及裂變產(chǎn)物成分隨時間的演變過程,不但會通過改變?nèi)己膮^(qū)的核素種類及密度影響中子輸運過程,也為熱工水力計算提供總功率及各燃耗區(qū)的實際功率。在高保真燃耗計算中,根據(jù)目標問題的實際規(guī)模,燃耗區(qū)數(shù)量可達數(shù)十萬甚至上百萬,需要逐個計算和存儲各燃耗區(qū)的相關(guān)參數(shù),包括中子通量密度、反應截面、功率等物理量,導致計算時間和存儲開銷顯著增加[4-5]。即使借助高性能計算機的并行能力,全堆芯高保真模擬也需要數(shù)百核時的計算時間及數(shù)百GB的存儲開銷,計算成本顯著增加,同時難以獲得高保真模型實時求解所需的加速效果[6-7]。

近年來,為了減少高保真燃耗計算中的存儲開銷,提高計算效率,主要圍繞以下兩個方面開展研究工作。一是通過先進的理論模型和并行計算方法提高上游中子輸運計算和熱工水力計算的效率和精度[8-9],為燃耗計算提供精確的中子學參數(shù);二是針對反應堆燃耗計算中的物理過程和特點,開展了包括點燃耗計算方法[10]、燃耗鏈壓縮方法[11]、輸運-燃耗耦合計算方法[12]、并行計算方法[13]、模型降階方法等多方面的研究。盡管高保真堆芯物理模型對計算成本有較高要求,但模型降階方法提供了解決途徑。常見的模型降階方法包括本征正交分解(POD)、特征廣義分解(PGD)、縮減基方法(RBM)等。其中通過POD方法得到的降階模型自由度低、數(shù)據(jù)存儲量小,且為原系統(tǒng)的最佳逼近,在高維復雜系統(tǒng)的計算效率及精度等方面具有顯著優(yōu)勢[14]。該方法在中子輸運/擴散計算[15]、熱工水力計算[16]、堆芯功率在線監(jiān)測[17]等方面得到了廣泛的研究與應用,取得了良好的加速效果,但在高保真燃耗計算方面,國內(nèi)尚處于起步階段[18]。

本文基于POD方法,通過建立高保真燃耗計算的參數(shù)化降階模型,以1組正交基函數(shù)實現(xiàn)對各燃耗矩陣的快速構(gòu)造。進一步基于離線-在線數(shù)值求解策略,通過選用JAEA發(fā)布的MOX燃料柵元基準題進行計算分析,初步驗證本方法的正確性和可行性。

1 模型降階方法

1.1 POD方法

POD方法通過全階模型的仿真數(shù)據(jù)或?qū)嶒灁?shù)據(jù)構(gòu)建自相關(guān)矩陣,求解自相關(guān)矩陣的特征向量獲得降階模態(tài),并將全階模型投影到數(shù)個最大特征值對應降階模態(tài)張成的子空間上以達到降階目的。

g?n

(1)

式中:ai(t)為模態(tài)系數(shù);g為由全階模型模態(tài)能量占比截斷后的模態(tài)數(shù);φi(x)為模態(tài)正交基。

實際上,POD降階過程是將全階模型不同時刻的響應y(x,ti)與其在該正交基張成的子空間Φ上的投影誤差范數(shù)等效為一個最值問題[19]:

(2)

式中:(·,·)為Hilbert空間上的內(nèi)積;‖·‖2為L2范數(shù)的平方;I(n)為n階單位矩陣;G為響應y(x,ti)在子空間Φ上投影平方的最大值。

進一步可將式(2)的最值問題轉(zhuǎn)換為特征值求解問題:

(3)

式中:S為自相關(guān)矩陣,S∈Rm×m;λi為S的第i個特征值;Y為全階模型不同時刻的響應集合。由式(3)獲得模態(tài)正交基后利用式(1)即可完成全階模型的降階。

1.2 奇異值分解

由1.1節(jié)可知,POD方法的關(guān)鍵在于確定1組合適的正交基。本文利用奇異值分解(SVD)進行求解,該方法廣泛用于降維算法中的特征分解。奇異值分解和特征分解不同,并不要求待分解的矩陣為方陣。對1.1節(jié)中的矩陣Y進行SVD分解,可得:

(4)

由奇異值分解過程,結(jié)合式(3)和式(4)可得:正交矩陣U的前r列構(gòu)成矩陣Y中列空間的標準正交基;正交矩陣V的前r列構(gòu)成矩陣Y中行空間的標準正交基。因此,正交基Φ=Um×r=(φ1,φ2,…,φr)。

為了盡可能少地利用模態(tài)正交基完成式(1)中y(x,ti)的逼近,可以通過模態(tài)能量占比E(g)進行模態(tài)截斷。E(g)定義為前g個特征值和與總特征值和的比值,用來衡量各模態(tài)的還原能力:

(5)

式(5)中,自相關(guān)矩陣S的特征值與各奇異值存在如下關(guān)系:

(6)

E(g)越接近1,表明模態(tài)基函數(shù)所包含原全階模型的特征信息越完整。通常E(g)應大于0.95,對應所截斷的模態(tài)基函數(shù)可將原系統(tǒng)投影到更低維的子空間上,實現(xiàn)模型降階。

2 參數(shù)化燃耗降階模型構(gòu)建方法

2.1 參數(shù)化燃耗降階模型

反應堆內(nèi)燃料中核素成分隨時間的變化情況可用燃耗方程組描述,其矩陣形式及解為:

(7)

式中:N(t)為核素的密度向量;A為燃耗矩陣,由系統(tǒng)中所有核素的反應率和衰變率構(gòu)成;N(0)為初始時刻的核素密度;t為時間。矩陣A可進一步表達為如下形式:

A=Ar+Ad

(8)

式中:Ar為燃耗系數(shù)矩陣,其矩陣元由中子通量密度和微觀反應截面的乘積構(gòu)成;Ad為衰變矩陣,其矩陣元由衰變常量與衰變分支比的乘積構(gòu)成,在整個燃耗過程中該矩陣不變。

由POD降階方法可將Ar和中子通量密度φ用兩組正交基函數(shù)表示:

(9)

式中:an、bn為待定的模態(tài)系數(shù);r1、r2為模態(tài)階數(shù);ξn、φn為模態(tài)正交基。

在燃耗計算過程中需要對各燃耗區(qū)的Ar進行反復計算和存儲,計算成本較高。將φ作為輸入?yún)?shù),可得Ar的參數(shù)化形式:

Ar=f(φ)

(10)

由式(9)和式(10)可進一步得到φ與Ar的映射關(guān)系,即參數(shù)化燃耗降階模型:

(11)

2.2 模態(tài)基函數(shù)構(gòu)造方法

(12)

式中:m1為瞬態(tài)快照Ak的維數(shù);m2為中子能群數(shù);SAr∈Rm1×n;Sφ∈Rm2×n;n為快照矩陣的采樣長度,由燃耗區(qū)數(shù)、燃耗深度、功率水平等輸入?yún)?shù)確定。

對快照矩陣SAr與Sφ進行奇異值分解可得:

(13)

由奇異值分解可知,左奇異矩陣U1和U2中的列向量即為原矩陣列空間內(nèi)的1組標準正交基,即:

(14)

對于燃耗系統(tǒng)在t時刻的中子通量密度φt和燃耗系數(shù)矩陣At,根據(jù)實際計算精度要求進行模態(tài)能量截斷得到降階模態(tài)數(shù)r后,基于式(9),可將At和φt表示為:

(15)

式中:U1,cut=(ξ1,ξ2,…,ξr),U1,cut∈Rm1×r;U2,cut=(φ1,φ2,…,φr),U2,cut∈Rm2×r。

將SAr降階后的右奇異矩陣V1,cut左乘Sφ得到1組新的中子通量密度基函數(shù)Sψ:

(16)

式中:vn,r為矩陣V1,cut中第n行r列的矩陣元;Sψ∈Rm2×r。

(17)

進一步通過參數(shù)化燃耗降階模型,得到φt與對應系數(shù)矩陣At的映射關(guān)系:

U1,cut(c1σ1,c2σ2,…,crσr)T

(18)

最后由式(15)和式(18)可得:

an=cnσnn=1,2,…,r

(19)

3 離線-在線計算策略

基于參數(shù)化燃耗降階模型進行燃耗計算可分為離線和在線兩個階段:離線階段完成模態(tài)基函數(shù)的構(gòu)造;在線階段基于離線階段預先存儲的數(shù)據(jù),通過簡單的代數(shù)運算獲得中子通量密度對應的燃耗系數(shù)矩陣。

3.1 離線階段

離線階段計算流程如圖1所示:1) 根據(jù)目標問題的代表性參數(shù)確定快照矩陣的采樣長度,通過中子輸運-燃耗計算獲得燃耗系數(shù)矩陣和中子通量密度構(gòu)成的快照矩陣SAr和Sφ(式(12)),同時保存衰變矩陣Ad;2) 基于POD方法分別獲得快照矩陣SAr和Sφ的奇異矩陣和奇異值(式(13));3) 根據(jù)計算要求進行模態(tài)能量截斷確定降階模態(tài)數(shù)(式(5)),構(gòu)造燃耗系數(shù)矩陣和中子通量密度的模態(tài)基函數(shù)U1,cut和Sψ(式(15)和式(16));4) 存儲模態(tài)基函數(shù)U1,cut、Sψ及衰變矩陣Ad。

圖1 離線階段計算流程圖Fig.1 Flow chart of offline phase

3.2 在線階段

在線階段計算流程如圖2所示:1) 進行中子輸運計算獲得中子通量密度φ,無需統(tǒng)計和存儲各燃耗區(qū)各核素的微觀反應率;2) 基于離線階段預先存儲的模態(tài)基函數(shù),將φ投影至Sψ上獲得還原燃耗系數(shù)矩陣的模態(tài)系數(shù)(式(19));3) 基于離線階段預先存儲的衰變矩陣Ad,完成對燃耗矩陣A的構(gòu)造(式(8));4) 基于燃耗矩陣A求解燃耗方程,完成當前燃耗步長計算,返回第1步進行下一個燃耗步長的計算,直至步長末。

圖2 在線階段計算流程圖Fig.2 Flow chart of online phase

4 數(shù)值驗證

為驗證基于POD降階的燃耗計算方法的正確性和有效性,本文選用JAEA發(fā)布的MOX燃料柵元基準題進行計算分析,該基準題的平均卸料燃耗深度為70 GW·d/t,功率密度為36.6 W/gHM,幾何模型如圖3所示,包括3個材料填充區(qū)域,由內(nèi)至外分別為燃料、包殼和慢化劑,詳細材料及工況信息參見文獻[20]。該基準題提供了16種燃耗計算程序的計算結(jié)果,包括無限增殖因數(shù)、部分重核及裂變產(chǎn)物的核素密度等參考值,并以這些結(jié)果的平均值加減兩倍標準差作為置信區(qū)間。本文計算程序選用自主研發(fā)的燃耗計算程序AMAC[21],燃耗數(shù)據(jù)庫選用JAEA提供的簡化燃耗鏈,包括21個重核、49個裂變產(chǎn)物和1個偽核素,該簡化燃耗鏈在輕水堆燃耗計算中得到了廣泛的驗證與應用[22]。

圖3 MOX燃料柵元幾何示意圖Fig.3 Geometrical configuration of MOX fuel pin cell

4.1 無限增殖因數(shù)

該基準題提供了0、0.1、5、10、15、20、30、50和70 GW·d/t共計9個燃耗深度下的基準值。本文共計劃分35個燃耗步長:第1個步長為0.1 GW·d/t;第2個步長為1 GW·d/t,并以1 GW·d/t為步長間隔遞增至第11個步長10 GW·d/t;第12個步長為12.5 GW·d/t,并以2.5 GW·d/t為步長間隔遞增至第35個步長70 GW·d/t。徑向?qū)θ剂蠀^(qū)等體積劃分3層。圖4示出基于MOX燃料柵元建立全階燃耗模型計算的無限增殖因數(shù)kinf及與其他16種計算程序計算結(jié)果平均值的相對偏差。由圖4可知,本文計算的kinf均在平均值的置信范圍內(nèi),與平均值的相對偏差在整個燃耗過程中始終小于100 pcm。計算結(jié)果驗證了AMAC對該基準題計算的正確性。

圖4 MOX燃料柵元kinf的計算結(jié)果Fig.4 Calculation result of kinf for MOX fuel pin cell

進一步基于本文提出的參數(shù)化燃耗降階模型對全階模型進行降階。在離線階段以模態(tài)能量占比E(g)≥0.99進行模態(tài)截斷,分別獲得燃耗系數(shù)快照矩陣(SAr∈R5 041×105)和中子通量密度快照矩陣(Sφ∈R107×105)前50階的模態(tài)基函數(shù)U1,cut和Sψ。在線階段的計算過程中無需統(tǒng)計和存儲各燃耗區(qū)中核素的反應率信息,采用32階最佳一致逼近多項式方法求解燃耗方程[23],并以全階模型的計算結(jié)果作為參考值。

圖5示出采用不同模態(tài)基函數(shù)計算的kinf及與全階模型kinf相對偏差隨燃耗深度的變化,同時以基準題提供的16種計算程序計算結(jié)果平均值的正負兩倍標準差作為置信區(qū)間。由圖5可知,當采用1階模態(tài)基函數(shù)進行模態(tài)還原時,kinf與全階模型的相對偏差均在500 pcm以下,說明該模態(tài)基函數(shù)已包含了全階燃耗模型的主要特征信息,具有良好的模態(tài)還原能力。采用5階降階模型計算時,kinf與全階燃耗模型的相對偏差減少至100 pcm以內(nèi)。繼續(xù)提高模態(tài)階數(shù)至10階時,可以發(fā)現(xiàn)kinf與全階模型的相對偏差已非常小,皆在5 pcm以內(nèi)?？紤]到20、30、40和50階降階模型計算的kinf與全階模型在各燃耗步長的相對偏差皆小于5 pcm,因此在圖5中并未展示?？傮w而言,通過不同模態(tài)階數(shù)構(gòu)造的降階模型與全階模型的計算結(jié)果相比,對于kinf的計算皆具有良好的精度。

圖5 不同模態(tài)截斷下kinf的計算結(jié)果Fig.5 Calculation result of kinf with different basis functions

4.2 核素密度

基準題提供了17個重核和12個裂變產(chǎn)物的核素密度。圖6示出29個核素在燃耗過程中相對于全階模型的平均相對偏差。由圖6可知,隨階數(shù)的提高,所有核素密度的相對偏差逐漸降低。1階降階模型的計算精度較差,大部分裂變產(chǎn)物的平均相對偏差大于10%。5階降階模型得到的大部分核素密度的相對偏差低于1%,其中核素245Cm的相對偏差達到3.74%。10階降階模型的相對偏差已全部低于0.1%。

圖6 29個核素密度的平均相對偏差Fig.6 Average relative deviation of 29 nuclide densities

本文進一步選取了核素235U、239Pu、155Gd和149Sm進行詳細分析,以全階模型的計算結(jié)果作為參考值,以16種計算程序計算結(jié)果平均值的正負兩倍標準差作為置信區(qū)間。圖7示出235U、239Pu、155Gd和149Sm 4種核素的核素密度及相對偏差隨燃耗深度的變化。

圖7 不同模態(tài)截斷下核素密度的計算結(jié)果Fig.7 Calculation result of nuclide density with different basis functions

由圖7a和圖7b可知,對于重核核素235U和239Pu,基于降階模型得到的核素密度隨燃耗深度的變化趨勢與全階模型相同,且皆在基準值的置信區(qū)間內(nèi)。在燃耗深度為5 GW·d/t時,通過1階降階模型計算得到的235U 核素密度的相對偏差達到最大值,為4.175%。隨著模態(tài)階數(shù)的提高,5階和10階降階模型的計算精度遠高于1階降階模型,其相對偏差在整個燃耗過程中分別始終低于0.1%和0.005%。當模態(tài)階數(shù)進一步提高到20～50階時,通過相應降階模型計算得到的235U和239Pu的核素密度與10階降階模型的計算精度相當,約在0.005%以內(nèi)。

如圖7c和圖7d所示,基于降階模型得到的155Gd和149Sm的核素密度隨燃耗深度變化的趨勢落在基準值的置信區(qū)間內(nèi)。壽期初155Gd和149Sm的核素密度快速增長,與全階模型計算結(jié)果相比,通過1階降階模型得到155Gd和149Sm的核素密度相對偏差分別為16.645%和16.549%。隨燃耗的加深,155Gd和149Sm的相對偏差先減小后增加,壽期末的相對偏差分別為11.996%和4.519%。雖然1階降階模型對于155Gd和149Sm的計算相對偏差較大,但其核素密度的量級較小,因此對1階降階模型計算kinf的精度影響較小。表1列出不同模態(tài)截斷下155Gd和149Sm的核素密度與全階模型計算結(jié)果的相對偏差。當采用5階降階模型時,整個燃耗過程中155Gd和149Sm的核素密度相對偏差始終小于5%,處于合理的誤差范圍內(nèi)。進一步提高模態(tài)階數(shù)至10階,155Gd和149Sm的核素密度與全階模型的相對偏差在0.1%以內(nèi)?？梢妼τ?55Gd和149Sm這些具有較大中子吸收截面的核素,僅采用1階降階模型還難以精確獲得核素密度隨燃耗深度變化的計算結(jié)果,需要采用更多的模態(tài)基函數(shù)。相比于1階降階模型,5階和10階降階模型對于核素密度的計算具有更好的精度。

表1 不同模態(tài)截斷下155Gd和149Sm核素密度的相對偏差Table 1 Relative deviation of nuclide density of 155Gd and 149Sm with different basis functions

4.3 計算效率與存儲開銷

燃耗計算的耗時主要包括系數(shù)矩陣構(gòu)建、燃耗方程求解及數(shù)據(jù)更新傳遞3部分。表2列出不同階數(shù)下各降階模型和傳統(tǒng)燃耗的計算時間。由表2可見,采用降階模型構(gòu)建系數(shù)矩陣的最大加速比為82.56,整個燃耗計算耗時的最大加速比為18.23。相比于傳統(tǒng)計算方法,降階模型在在線階段無需統(tǒng)計通量歸一因子和各燃耗區(qū)內(nèi)各核素的單群反應率,減少了數(shù)據(jù)傳遞時間,提高了計算效率。由于不同階數(shù)的降階模型僅影響系數(shù)矩陣的精度和構(gòu)建時間,并不影響系數(shù)矩陣的規(guī)模,因此隨著降階模型階數(shù)提高,燃耗方程的求解耗時不變。

表2 燃耗降階模型的計算耗時Table 2 Calculation time of burnup reduced-order model

燃耗計算的內(nèi)存開銷包括定義數(shù)據(jù)和統(tǒng)計數(shù)據(jù)[7]。定義數(shù)據(jù)包括燃耗步長、功率密度、各燃耗區(qū)的通量、核素密度、截面信息等。統(tǒng)計數(shù)據(jù)主要指計算過程中更新及存儲的數(shù)據(jù)。在對本文算例的計算過程中,定義數(shù)據(jù)包括各燃耗區(qū)存儲長度為107的通量數(shù)組和長度為71的核素密度數(shù)組。此外還需存儲每個核素的4種反應截面:裂變截面、俘獲截面、(n,2n)截面、(n,3n)截面,上述數(shù)據(jù)皆為雙精度型。統(tǒng)計數(shù)據(jù)主要包括用于燃耗方程求解中的稀疏矩陣。燃耗系數(shù)矩陣采用三元組存儲,由于非零元素為1 182個,因此需要2個長度為1 182的整型數(shù)組和1個長度為1 182的雙精度型數(shù)組。同時對于各燃耗區(qū)共用的衰變矩陣,需要2個長度為84的整型數(shù)組和1個長度為84的雙精度型數(shù)組。

表3列出不同階數(shù)下各降階模型的內(nèi)存開銷及相對于傳統(tǒng)燃耗計算的內(nèi)存占比。當采用降階模型時,內(nèi)存開銷主要由奇異值、通量基函數(shù)和系數(shù)矩陣基函數(shù)所決定。由表3可知,隨著階數(shù)提高,內(nèi)存占比逐步提高。即使采用50階降階模型進行計算,其內(nèi)存開銷仍為傳統(tǒng)燃耗計算的34%。燃耗降階模型對內(nèi)存開銷的節(jié)省效果較為顯著。

表3 燃耗降階模型的內(nèi)存開銷Table 3 Memory requirement of burnup reduced-order model

由4.1節(jié)和4.2節(jié)可知,當采用5階降階模型時:在計算精度方面,kinf的相對偏差在整個燃耗過程中始終低于100 pcm,所有核素密度的相對偏差均低于5%,大部分核素密度的相對偏差低于1%;在計算成本方面,5階降階模型的加速比為17.26,內(nèi)存開銷為傳統(tǒng)燃耗計算的5%。隨著階數(shù)進一步提高,高階降階模型對計算精度和效率的提升已十分有限,但內(nèi)存占比會逐步增加。因此對于降階模型階數(shù)的選擇,在滿足計算精度的條件下,應充分考慮內(nèi)存開銷的影響,以獲得計算精度、效率和內(nèi)存開銷間的平衡。

5 結(jié)論

本文基于POD方法,建立了高保真燃耗計算的參數(shù)化降階模型,實現(xiàn)了對燃耗矩陣的快速構(gòu)造。進一步基于“離線-在線”計算策略,選用JAEA發(fā)布的MOX燃料柵元基準題對本方法進行了初步驗證與分析,結(jié)論如下。

1) 對于kinf的計算,僅需通過1階降階模型便可獲得與全階模型相對偏差小于500 pcm的計算精度;對于29個核素的計算,由5階降階模型得到的大部分核素密度的相對偏差低于1%。進一步重點分析了核素235U、239Pu、155Gd和149Sm。對于重核235U和239Pu的核素密度計算,1、5和10階降階模型的相對偏差在整個燃耗過程中分別在5%、0.1%和0.005%以內(nèi);但對于中子吸收截面較大的核素155Gd和149Sm,1階降階模型的計算精度較差,5階降階模型的相對偏差在5%以內(nèi),10階降階模型的相對偏差可達到0.1%以內(nèi)。

2) 相比于傳統(tǒng)燃耗計算方法,不同降階模型獲得的加速比最大可達18.23。相比于傳統(tǒng)燃耗計算,5階降階模型的加速比可達17.26,提高了計算效率。

3) 燃耗降階模型無需計算和存儲各燃耗區(qū)內(nèi)各核素的反應截面、反應率等參數(shù)信息,同時能夠減少統(tǒng)計數(shù)據(jù)的內(nèi)存開銷。5階降階模型能夠在獲得良好計算精度和效率的同時,內(nèi)存開銷僅為傳統(tǒng)燃耗計算的5%。

本文僅對高保真燃耗計算的參數(shù)化降階模型的正確性和可行性做了初步研究和分析。未來將進一步研究通過少量基函數(shù)完成一定規(guī)模燃耗區(qū)的矩陣構(gòu)建方法,以顯著降低全堆芯高保真燃耗的存儲開銷,提高計算效率。同時在離線階段將計算規(guī)模控制在單柵元或超柵元階段,以減少離線階段的計算成本。