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

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

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

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

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

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

1 POD-Galerkin理論方法

1.1 本征正交分解

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

(1)

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

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

(2)

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

(3)

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

1.2 POD-Galerkin方法

Galerkin投影法[10-11]的實現(xiàn)是通過將試函數(shù)本身當(dāng)作權(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ù)、中子標(biāo)通量、吸收截面、中子產(chǎn)出截面和裂變能譜;Σs,g→g′為從g群到g′群的散射截面。

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

(8)

(9)

(10)

(11)

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

(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ù)充分?jǐn)_動,通過大量重復(fù)全階計算,得到不同狀態(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基準(zhǔn)題

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

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

圖2 TWIGL基準(zhǔn)題幾何模型[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ù)文獻(xiàn)調(diào)研[17],核截面擾動范圍一般不超過40%,因此本文將擾動范圍定為基準(zhǔn)題給定值的20%。擾動基準(zhǔn)題中的14個不確定參數(shù),在擾動范圍內(nèi)隨機采樣100個點,進行全階計算生成樣本空間,對樣本空間進行POD分解,將分解所得特征值進行降序排列,計算每階特征值占比,前10階的特征值占比曲線(對數(shù)坐標(biāo))如圖3所示,可以看出,隨著階數(shù)增加,各階特征值占比迅速下降,由式(3)得精度條件ε=99.99%時,快中子群和熱中子群所需POD基的階數(shù)分別為r1=2、r2=2。

圖3 前10階特征值占比

在原基準(zhǔn)題條件下,構(gòu)建出降階模型,全階模型本征值計算結(jié)果為0.913 23,基準(zhǔn)解[18]為0.913 21,兩者偏差為2 pcm。降階模型本征值計算結(jié)果為0.913 22,與基準(zhǔn)解偏差為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基的物理含義

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

Σ=USVT

(23)

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

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

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

a——1階;b——2階

a——1階;b——2階

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

快群第2種變化模式如圖6b所示,有兩個區(qū)域具有明顯相反的變化特征,介質(zhì)1、2所在區(qū)域為同負(fù)向變化,而反應(yīng)堆邊緣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ū)交界處通量均相應(yīng)增加。

綜合兩種變化模式來看,隨著擴散系數(shù)增大,最主要的變化特征是整體通量減小,而更細(xì)節(jié)的變化則是1、2介質(zhì)區(qū)通量減小,交界處的通量增大。而當(dāng)擾動全部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)分布,樣本均值為基準(zhǔn)題給出的參數(shù)值,標(biāo)準(zhǔn)差是參考值的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。可以發(fā)現(xiàn),隨著樣本量的增加,keff的數(shù)學(xué)期望偏差基本穩(wěn)步下降,樣本量為100時,LHS下全階與降階的數(shù)學(xué)期望偏差為1 pcm,屬于較小誤差;相同樣本量時LHS下數(shù)學(xué)期望偏差小于SRS;另外隨著樣本量增加,降階計算相對于全階計算的時間占比逐漸增加。

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

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

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

3 結(jié)論

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

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

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

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

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