基于降階模型的中子擴散特征值問題的不確定性分析研究

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

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

堆芯物理計算是反應堆設計與分析的基礎,隨著堆芯建模和仿真精度提升,其計算代價和復雜度亦隨之提升[1]。此外,核數據作為堆芯物理計算的重要輸入參數,其不確定性會影響堆芯物理計算結果的精度。基于抽樣統計不確定性分析方法[2],開展精細化堆芯物理計算的不確定性量化研究,其將面臨高昂計算代價。因此有必要開展堆芯物理降階模型(ROM)在不確定性分析中的應用研究。

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

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

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

1 POD-Galerkin理論方法

1.1 本征正交分解

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

(1)

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

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

(2)

其中,r的選取規則為:

(3)

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

1.2 POD-Galerkin方法

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

(4)

(5)

(6)

(7)

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

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

(8)

(9)

(10)

(11)

根據POD方法,可設:

(12)

(13)

(14)

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

(15)

(16)

(17)

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

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

(18)

Ar=Ur,TAUr

(19)

Br=Ur,TBUr

(20)

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

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

(21)

(22)

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

1.3 降階模型構建流程

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

2 數值結果

2.1 TWIGL基準題

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

表1 TWIGL截面參數[16]

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

2.2 計算結果與分析

根據1.3節計算流程,自主編寫穩態二維雙群擴散問題的全階、降階計算模型。

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

圖3 前10階特征值占比

在原基準題條件下,構建出降階模型,全階模型本征值計算結果為0.913 23,基準解[18]為0.913 21,兩者偏差為2 pcm。降階模型本征值計算結果為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與其轉置矩陣相乘構成協方差矩陣Σ=YYT,對協方差矩陣特征值分解,即可得到特征值和特征向量,其中前幾階較大的特征值對應的特征向量即為POD基,即POD基來自于樣本協方差矩陣的特征值分解。奇異值分解被認為是最佳的分解協方差矩陣的方法之一。對于協方差矩陣Σ,存在分解使得:

Σ=USVT

(23)

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

數學上,U的各列即代表網格點數據的不同變化特征方向,以第1列兩個網格點φ1=(k1,k2)T為例,當基向量中兩點k1、k2均為正時,代表1點數據變大時,2點數據也變大,即同正向變化;當k1、k2均為負時,代表1點數據變小時,2點數據也變小,即同負向變化;當k1為正、k2為負時,代表1點數據變大時,2點數據反而變小,即相反變化。基向量中每一網格點處數字越大,代表該點處變化程度越強。

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

a——1階;b——2階

a——1階;b——2階

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

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

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

a——1階;b——2階

a——1階;b——2階

2.4 不確定性分析

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

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

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

兩種抽樣方法下不同樣本量全階與降階keff計算結果列于表2、3。可以發現,隨著樣本量的增加,keff的數學期望偏差基本穩步下降,樣本量為100時,LHS下全階與降階的數學期望偏差為1 pcm,屬于較小誤差;相同樣本量時LHS下數學期望偏差小于SRS;另外隨著樣本量增加,降階計算相對于全階計算的時間占比逐漸增加。

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

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

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

3 結論

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

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

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

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

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