梁鑫源,王毅箴,郝 琛,*

(1.哈爾濱工程大學 核安全與仿真技術(shù)國防重點學科實驗室,黑龍江 哈爾濱 150001;2.清華大學 核能與新能源技術(shù)研究院,先進反應堆工程與安全教育部重點實驗室,北京 100084)

堆芯物理計算是反應堆設計與分析的基礎,隨著堆芯建模和仿真精度提升,其計算代價和復雜度亦隨之提升[1]。此外,核數(shù)據(jù)作為堆芯物理計算的重要輸入?yún)?shù),其不確定性會影響堆芯物理計算結(jié)果的精度?；诔闃咏y(tǒng)計不確定性分析方法[2],開展精細化堆芯物理計算的不確定性量化研究,其將面臨高昂計算代價。因此有必要開展堆芯物理降階模型(ROM)在不確定性分析中的應用研究。

降階模型根據(jù)構(gòu)建方法可以分為物理驅(qū)動和數(shù)據(jù)驅(qū)動兩大類。物理驅(qū)動方法包括各種投影技術(shù),如伽遼金(Galerkin)投影法,該方法需要修改控制方程,因此需要直接訪問全階模型(FOM)源代碼。某些情況下,直接剖析和修改全階模型是具有挑戰(zhàn)性的,但其優(yōu)點是機理清晰、可解釋性強。數(shù)據(jù)驅(qū)動方法包括克里金法[3-4]和動態(tài)模式分解法[5]等,僅需在“黑箱”模式下運行模擬,以生成與某些定義輸入相對應的輸出,把這些輸入和輸出數(shù)據(jù)用作訓練數(shù)據(jù),通過機器學習等方法構(gòu)建輸入?yún)?shù)與降階模型的函數(shù)關(guān)系即可,雖操作簡單,但不如物理驅(qū)動方法機理清晰。龔禾林等[1]基于本征正交分解(POD)技術(shù)和機器學習構(gòu)建了數(shù)據(jù)驅(qū)動式中子物理快速計算模型。Elzohery等[5]用一維兩群瞬態(tài)中子擴散問題對比物理驅(qū)動和數(shù)據(jù)驅(qū)動降階模型的效果,其中基于物理驅(qū)動的POD-Galerkin降階方法表現(xiàn)出最佳計算精度。

對于降階模型應用于不確定性分析的研究,Elzohery等[6]對二維兩群瞬態(tài)中子擴散問題利用貪婪采樣得到的降階模型進行不確定性分析,對比全階與降階計算的樣本均值和標準差的相對誤差,驗證了降階模型應用于不確定性分析的可行性。

本文將POD方法與經(jīng)典的求解偏微分方程的Galerkin投影法結(jié)合,針對二維兩群中子擴散問題構(gòu)建物理驅(qū)動的降階模型,探析POD基函數(shù)的物理含義,并對POD-Galerkin降階模型用于抽樣統(tǒng)計不確定性分析的可行性進行研究。

1 POD-Galerkin理論方法

1.1 本征正交分解

POD方法的核心是在核反應截面數(shù)據(jù)擾動下,從已有的中子通量場變化數(shù)據(jù)中找到1組最優(yōu)的正交基來代表其數(shù)據(jù)變化特征[7]。最常用的方法是快照法,最早由Sirovich等[8]提出?？煺帐侵覆煌朔磻孛孑斎胂?堆芯中子通量場數(shù)值解的空間分布。設Y={y1(x),y2(x),…,ym(x)}是1個由足夠多快照組成的快照矩陣,也稱樣本空間。其中:m為快照數(shù);x為空間位置向量;yi(x)為第i個核反應截面樣本輸入下堆芯中子通量場數(shù)值解的空間分布。假設核反應截面擾動下,中子通量分布樣本空間中的樣本點yj(x)函數(shù)展開如式(1)所示:

(1)

式中:φi(x)為基函數(shù)或基向量;ci為基函數(shù)對應的系數(shù)。

本征正交分解的基本步驟為:1) 對快照矩陣進行奇異值分解(SVD)[1,6,9],奇異值的平方即特征值λ,左奇異矩陣即基函數(shù)矩陣;2) 將得到的基函數(shù)按照特征值的大小降序排列,按照所需精度截取前r階基函數(shù)即POD基,用POD基對yj(x)進行低維近似,則yj(x)可表示為:

(2)

其中,r的選取規(guī)則為:

(3)

式中,ε為根據(jù)所需精度設定的值,一般情況下取99.99%即可。

1.2 POD-Galerkin方法

Galerkin投影法[10-11]的實現(xiàn)是通過將試函數(shù)本身當作權(quán)函數(shù)來構(gòu)造微分方程的積分形式,從而求得微分方程的近似解。在與POD方法結(jié)合構(gòu)建低階模型時,其試函數(shù)即POD基,即將POD基作為權(quán)函數(shù)來構(gòu)建低階模型。本文對穩(wěn)態(tài)二維兩群中子擴散問題構(gòu)建降階模型:

(4)

(5)

(6)

(7)

式中:g為群數(shù);Dg、Φg、Σa,g、νΣf,g、χg分別為g群的擴散系數(shù)、中子標通量、吸收截面、中子產(chǎn)出截面和裂變能譜;Σs,g→g′為從g群到g′群的散射截面。

將式(4)轉(zhuǎn)化為矩陣形式:

(8)

(9)

(10)

(11)

根據(jù)POD方法,可設:

(12)

(13)

(14)

根據(jù)Galerkin投影法,將式(8)投影到POD基上:

(15)

(16)

(17)

式中,c1和c2分別為1群和2群的系數(shù)向量。

將式(15)簡化為如下形式:

(18)

Ar=Ur,TAUr

(19)

Br=Ur,TBUr

(20)

式中,Ar、Br均為r×r階矩陣。

本文所研究的擴散問題為兩群,考慮分群構(gòu)建降階模型:

(21)

(22)

較多結(jié)果表明,對于大多數(shù)工程問題,r≈10～100時即能構(gòu)造出滿足精度要求的降階模型[12-15],降階模型正是因此得以提高計算速度。

1.3 降階模型構(gòu)建流程

降階模型的構(gòu)建計算流程如圖1所示。具體步驟如下:1) 對擴散系數(shù)及各種截面等不確定參數(shù)充分擾動,通過大量重復全階計算,得到不同狀態(tài)下的中子通量數(shù)據(jù),構(gòu)成足夠充分的樣本空間;2) 對樣本空間進行奇異值分解,根據(jù)式(3)進行截斷,確定基函數(shù)的階數(shù)r,從而選取最優(yōu)的r階POD基;3) 構(gòu)建全階系數(shù)矩陣A、B。再用Galerkin投影法計算出降階系數(shù)矩陣Ar、Br,完成降階模型的構(gòu)建。

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

2.1 TWIGL基準題

該基準題包含3種不同介質(zhì)區(qū)域,其1/4堆芯幾何布置如圖2所示[16],無外中子源,左邊界與下邊界為對稱邊界,右邊界與上邊界為零通量邊界,各區(qū)域截面參數(shù)列于表1[16],其中χ1=1、χ2=0。

表1 TWIGL截面參數(shù)[16]

圖2 TWIGL基準題幾何模型[16]

2.2 計算結(jié)果與分析

根據(jù)1.3節(jié)計算流程,自主編寫穩(wěn)態(tài)二維雙群擴散問題的全階、降階計算模型。

Elzohery等[6]在研究中發(fā)現(xiàn),不同的網(wǎng)格尺寸對降階模型計算精度的影響較弱,幾乎可以忽略。網(wǎng)格劃分越密,所需計算時間越長,本文采取120×120的網(wǎng)格尺寸。據(jù)文獻調(diào)研[17],核截面擾動范圍一般不超過40%,因此本文將擾動范圍定為基準題給定值的20%。擾動基準題中的14個不確定參數(shù),在擾動范圍內(nèi)隨機采樣100個點,進行全階計算生成樣本空間,對樣本空間進行POD分解,將分解所得特征值進行降序排列,計算每階特征值占比,前10階的特征值占比曲線(對數(shù)坐標)如圖3所示,可以看出,隨著階數(shù)增加,各階特征值占比迅速下降,由式(3)得精度條件ε=99.99%時,快中子群和熱中子群所需POD基的階數(shù)分別為r1=2、r2=2。

圖3 前10階特征值占比

在原基準題條件下,構(gòu)建出降階模型,全階模型本征值計算結(jié)果為0.913 23,基準解[18]為0.913 21,兩者偏差為2 pcm。降階模型本征值計算結(jié)果為0.913 22,與基準解偏差為1 pcm。通過全階模型與降階模型計算的中子通量分布分別如圖4、5所示。降階與全階的中子通量最大相對誤差為:1群0.28%,2群0.28%。單次運算全階模型的計算時間為54.50 s,降階模型的計算時間為0.56 s,僅占全階的1.02%。

a——全階模型;b——降階模型;c——相對誤差

a——全階模型;b——降階模型;c——相對誤差

2.3 POD基的物理含義

在初步驗證了降階模型的準確性后,對POD基的物理含義進行淺析。將快照矩陣Y與其轉(zhuǎn)置矩陣相乘構(gòu)成協(xié)方差矩陣Σ=YYT,對協(xié)方差矩陣特征值分解,即可得到特征值和特征向量,其中前幾階較大的特征值對應的特征向量即為POD基,即POD基來自于樣本協(xié)方差矩陣的特征值分解。奇異值分解被認為是最佳的分解協(xié)方差矩陣的方法之一。對于協(xié)方差矩陣Σ,存在分解使得:

Σ=USVT

(23)

式中:U為左奇異矩陣,U=(φ1,φ2,…,φm),其列向量為互相正交的特征向量;S為實對角矩陣,主對角元素為由大到小排列的奇異值,奇異值的平方即特征值。

數(shù)學上,U的各列即代表網(wǎng)格點數(shù)據(jù)的不同變化特征方向,以第1列兩個網(wǎng)格點φ1=(k1,k2)T為例,當基向量中兩點k1、k2均為正時,代表1點數(shù)據(jù)變大時,2點數(shù)據(jù)也變大,即同正向變化;當k1、k2均為負時,代表1點數(shù)據(jù)變小時,2點數(shù)據(jù)也變小,即同負向變化;當k1為正、k2為負時,代表1點數(shù)據(jù)變大時,2點數(shù)據(jù)反而變小,即相反變化?；蛄恐忻恳痪W(wǎng)格點處數(shù)字越大,代表該點處變化程度越強。

每一階POD基分別代表中子通量場的一種變化模式,根據(jù)特征值占比截取前幾階基就是保留中子通量場最主要的變化模式。當只擾動介質(zhì)1、2的熱群擴散系數(shù)D2時,計算得ε=99.99%時r1=2、r2=2,分別繪制快群和熱群的前兩階基向量的空間分布,如圖6、7所示,可以看出,中子通量場的變化主要是由兩種變化模式組成的。

a——1階;b——2階

a——1階;b——2階

快群第1種變化模式如圖6a所示,表明了中子通量的一種同負向變化特征,即同時減小的特征,這種變化的影響效果從內(nèi)側(cè)3介質(zhì)區(qū)向外,先由弱變強,再由強變?nèi)?越靠近反應堆邊緣,影響效果越趨于零。即隨著介質(zhì)1、2的熱群擴散系數(shù)的增大,整體中子通量均有所減小,尤其是介質(zhì)1、2所在區(qū)域通量減小更多,越靠近邊緣,減小程度越趨于零。

快群第2種變化模式如圖6b所示,有兩個區(qū)域具有明顯相反的變化特征,介質(zhì)1、2所在區(qū)域為同負向變化,而反應堆邊緣3介質(zhì)區(qū)與1、2介質(zhì)區(qū)交界處呈同正向變化特征,內(nèi)側(cè)3介質(zhì)區(qū)和邊緣3介質(zhì)區(qū)的影響效果趨于零。即隨著介質(zhì)1、2的熱群擴散系數(shù)的增大,1、2介質(zhì)區(qū)通量均有所減小,而邊緣3介質(zhì)區(qū)與1、2介質(zhì)區(qū)交界處通量均相應增加。

綜合兩種變化模式來看,隨著擴散系數(shù)增大,最主要的變化特征是整體通量減小,而更細節(jié)的變化則是1、2介質(zhì)區(qū)通量減小,交界處的通量增大。而當擾動全部14個輸入?yún)?shù)時,計算得ε=99.99%時r1=2、r2=2,繪制快群和熱群前兩階基向量空間分布如圖8、9所示,與圖6、7對比可發(fā)現(xiàn),擾動參數(shù)量不同時所得基向量分布不同,即中子通量場的變化模式不同。

a——1階;b——2階

a——1階;b——2階

2.4 不確定性分析

為研究降階模型在提高不確定性計算速度方面的潛力,進行如下測試:令所有參數(shù)服從多元正態(tài)分布,樣本均值為基準題給出的參數(shù)值,標準差是參考值的20%,在擾動范圍內(nèi)用簡單隨機抽樣(SRS)和拉丁超立方抽樣(LHS)分別采樣1 000個樣本點,選擇10～100個樣本點來構(gòu)建降階模型,其余900個樣本點作測試點,分別進行全階與降階計算。

繪制LHS樣本量為100的全階與降階keff的分位數(shù)-分位數(shù)圖(Q-Q圖),如圖10所示。由圖10可見,全階與降階計算的keff分布一致性較強。其中全階計算時間為47 631.15 s,降階計算時間(包含基向量生成時間)為5 468.10 s,降階計算時間僅占全階計算時間的11.48%。

圖10 全階與降階keff的Q-Q圖

兩種抽樣方法下不同樣本量全階與降階keff計算結(jié)果列于表2、3?？梢园l(fā)現(xiàn),隨著樣本量的增加,keff的數(shù)學期望偏差基本穩(wěn)步下降,樣本量為100時,LHS下全階與降階的數(shù)學期望偏差為1 pcm,屬于較小誤差;相同樣本量時LHS下數(shù)學期望偏差小于SRS;另外隨著樣本量增加,降階計算相對于全階計算的時間占比逐漸增加。

表2 SRS不同樣本量全階與降階keff計算結(jié)果

表3 LHS不同樣本量全階與降階keff計算結(jié)果

不同抽樣方式將得到不同的樣本空間,在相同樣本量下,LHS能夠提供更為準確的代理模型,既能夠使樣本均衡地覆蓋輸入?yún)?shù)的分布區(qū)間,而且即使在樣本數(shù)量較少的情況下,也能夠?qū)斎雲(yún)?shù)的不確定度進行準確合理的表征。在相同樣本量下,基于LHS和SRS的降階模型的不確定性分析結(jié)果都能與全階不確定性分析結(jié)果具有較小誤差,但相較而言,LHS結(jié)果的誤差更小。因此,從TWIGL基準題測試結(jié)果來看,在POD-Galerkin降階建模中,相同樣本量下,LHS方法更建議采用。

3 結(jié)論

本文利用POD-Galerkin方法構(gòu)建出一種針對穩(wěn)態(tài)二維兩群中子擴散問題的物理驅(qū)動式降階模型,并用自主編寫的程序成功實施了全階與降階模型的計算。根據(jù)TWIGL二維兩群穩(wěn)態(tài)基準題進行了數(shù)值測試,并對其進行了不確定性分析,數(shù)值結(jié)果如下。

1) 本文構(gòu)建的降階模型能夠在保證計算精度的前提下,較快地完成中子通量及本征值的求解,降階模型keff計算結(jié)果與基準解偏差為1 pcm,單次運算降階模型的計算時間僅占全階的1.02%。

2) 在不確定性分析方面,該降階模型與全階模型計算出的keff屬于同一分布,且降階與全階的數(shù)學期望、標準差和峰度具有較小誤差;將降階模型構(gòu)造所需的全階模型計算時間考慮在內(nèi),降階不確定性分析計算時間占全階的11.48%,充分展現(xiàn)了降階模型在減少不確定性分析計算時長方面的較大潛力。

3) 在相同樣本量下,基于LHS和SRS的降階模型的不確定性分析結(jié)果都能與全階不確定性分析結(jié)果具有較小誤差,但LHS的結(jié)果誤差更小。因此在POD-Galerkin降階建模中,相同樣本量下,LHS方法更建議采用。

基于本文對降階模型應用于中子擴散特征值問題的不確定性分析研究,下一步的工作可以包括:在中子輸運問題上推廣POD-Galerkin模型降階方法,以提高中子輸運計算效率;在更多輸入?yún)?shù)擾動下,考慮參數(shù)之間的相關(guān)性,采用更高效的抽樣方法獲得質(zhì)量更高的樣本空間,進一步減小降階模型與全階模型的計算誤差。