深度學(xué)習(xí)方法求解中子輸運方程的微分變階理論

劉 東,王雪強(qiáng),張 斌,俞蔡陽,宮兆虎,陳奇隆

(1.中國核動力研究設(shè)計院,四川 成都 610213;2.中核核能軟件與數(shù)字化反應(yīng)堆工程技術(shù)研究中心,四川 成都 610213;3.四川大學(xué) 計算機(jī)學(xué)院,四川 成都 610065)

在反應(yīng)堆中子物理學(xué)領(lǐng)域,求解中子輸運方程是反應(yīng)堆物理計算的基本問題,對反應(yīng)堆工程的堆芯設(shè)計研發(fā)與運行支持具有至關(guān)重要的作用。求解輸運方程的傳統(tǒng)方法是以離散縱標(biāo)法(SN)、特征線法(MOC)為代表的確定論方法[1-2],和以蒙特卡羅方法為代表的統(tǒng)計學(xué)方法[3]。這些方法總體思路是從方程能量離散、空間離散、角度離散、并行加速等方面開展研究,從而獲得輸運方程的數(shù)值解。傳統(tǒng)方法各有特點,互為補充,在實際工程中已得到大量成功應(yīng)用。

近年來,隨著以深度學(xué)習(xí)為代表的新一代人工智能技術(shù)的發(fā)展,利用深度學(xué)習(xí)技術(shù)求解復(fù)雜偏微分方程的方法正成為數(shù)值計算領(lǐng)域前沿研究熱點[4]。不同于傳統(tǒng)的數(shù)據(jù)驅(qū)動機(jī)器學(xué)習(xí)方法,基于物理信息的神經(jīng)網(wǎng)絡(luò)(PINN)等深度學(xué)習(xí)框架[5]不需要預(yù)先收集大量數(shù)據(jù)作為機(jī)器學(xué)習(xí)的訓(xùn)練集,可直接正向求解復(fù)雜高階多維方程,并具有空間幾何限制少、邊界條件處理靈活、計算流程規(guī)范性好、實驗數(shù)據(jù)同化能力強(qiáng)等特點,已在求解熱工水力N-S方程等多個領(lǐng)域取得了重要突破[5-6]。在求解反應(yīng)堆中子擴(kuò)散方程領(lǐng)域也取得了系列進(jìn)展[7-10],這為中子輸運方程的求解開辟了一條新的潛在技術(shù)途徑,越來越受到業(yè)界的重視。

當(dāng)前深度學(xué)習(xí)技術(shù)成功求解的方程類型是微分方程[5-10],而中子輸運方程具有微分-積分方程的特點,其存在能量、角通量密度表示的中子散射與裂變多重定積分項,這造成了直接將現(xiàn)有方法應(yīng)用到中子輸運方程時存在很大困難。針對該問題,文獻(xiàn)[11]參考傳統(tǒng)SN方法[1-2],將定積分項的被積角度變量進(jìn)行離散,通過高斯求積組[1-2]構(gòu)造定積分項的近似離散和,形成損失函數(shù)進(jìn)行機(jī)器學(xué)習(xí),并給出了瞬態(tài)情況下一維幾何的離散時間、角度標(biāo)通量密度計算結(jié)果。

本文基于深度學(xué)習(xí)技術(shù)求解微分方程的基本方法,針對中子輸運方程存在積分項的特點,提出求解中子輸運方程的微分變階理論,將角通量密度定積分變換為其對應(yīng)的高階微分形式原函數(shù)進(jìn)行深度學(xué)習(xí),成功后再降階為原輸運方程,從而為解決輸運方程多重積分項帶來的深度學(xué)習(xí)方法困難探索新技術(shù)途徑,并可獲得完整的幾何-角通量密度連續(xù)性計算結(jié)果。

1 中子輸運方程的各種形式

中子輸運方程的完整形式參考文獻(xiàn)[1],其離散能群、裂變源/散射源各向同性的穩(wěn)態(tài)輸運方程在不考慮外源下的形式為:

(1)

式中:φ(r,Ω)為在r位置處角度為Ω的中子角通量密度,m-2·s-1;Σt為總反應(yīng)截面,m-1;Σs為散射截面,m-1;ν為每次裂變放出的平均中子數(shù);xg為裂變中子的能譜分布;Σf為裂變截面,m-1;keff為有效增殖因數(shù);g為分段能群編號,g=1,2,…,G。

當(dāng)考慮中子能量為單能時,式(1)簡化為:

(2)

特別地,式(2)的平板幾何形式[2]為:

(3)

式中:x為平板厚度;μ為直角坐標(biāo)下的方向余弦,其詳細(xì)定義參考文獻(xiàn)[1-2]。

式(2)的球幾何形式[2]為:

(4)

式中:r為球半徑;μ為球坐標(biāo)下的方向余弦,其詳細(xì)定義參考文獻(xiàn)[1-2]。

式(2)的無限圓柱幾何形式[2]為:

(5)

式(2)～(5)可描述為統(tǒng)一的偏微分-積分方程形式[12]:

(6)

2 微分變階理論基本原理

深度學(xué)習(xí)技術(shù)求解輸運方程的微分變階理論基本原理是:利用微積分?jǐn)?shù)理方法對中子輸運方程進(jìn)行微分升階,將其變換為關(guān)于角通量密度積分原函數(shù)表示的高階微分方程,進(jìn)而用原函數(shù)線性組合形式消除輸運方程角通量密度定積分項;同時,給出微分升階后的原函數(shù)特定解形成方法,確定不同角通量密度積分形式對應(yīng)的原函數(shù)定解約束,并將輸運方程的邊界條件映射為原函數(shù)表示的邊界條件;以微分升階變換后的原函數(shù)為目標(biāo),將控制方程殘差、邊界條件、原函數(shù)定解條件以及特征值約束加權(quán)形成統(tǒng)一損失函數(shù)進(jìn)行深度學(xué)習(xí),待神經(jīng)網(wǎng)絡(luò)函數(shù)逼近原函數(shù),可獲得原函數(shù)的數(shù)值解;此后,將原函數(shù)關(guān)于通量密度角度變量進(jìn)行求導(dǎo)實現(xiàn)微分降階,即得到原中子輸運方程的角通量密度函數(shù)分布,從而實現(xiàn)利用深度學(xué)習(xí)技術(shù)求解中子輸運方程。

為簡化起見,以單能臨界、穩(wěn)態(tài)具有一維/二維角通量密度形式的中子輸運方程為對象進(jìn)行詳細(xì)討論,其他形式輸運方程包含的角通量密度定積分項形式在數(shù)理方法上與其一致,可直接類推。

式(1)角通量密度φ(r,Ω)關(guān)于角度變量Ω為連續(xù)函數(shù)[1-3],按照微積分學(xué)原理[13],其存在關(guān)于角度變量Ω積分原函數(shù)F(r,Ω),由牛頓-萊布尼茲公式可表示如下。

1)Ω為一元變量

用方向余弦μ表示Ω[13],有:

F(r,μ1)-F(r,μ0),

(7)

式中,μ1、μ0分別為定積分項中Ω′積分上下限。

2)Ω為二元變量

用方向余弦ξ、幅角φ表示Ω[14],有:

F(r,ξ1,φ1)-F(r,ξ0,φ1)-

F(r,ξ1,φ0)+F(r,ξ0,φ0),

(8)

式中,ξ1、ξ0和φ1、φ0分別為方向余弦ξ、幅角φ表示的定積分項中Ω′積分上下限。

(9)

據(jù)此給出平板幾何、球幾何、圓柱幾何微分升階后的由原函數(shù)表示的輸運方程形式,其他幾何形式的方程可據(jù)此類推。

1) 平板幾何微分升階形式

將式(7)代入式(3),并令μ1=1,μ0=-1,則有:

(10)

2) 球幾何微分升階形式

將式(7)代入式(4),并令μ1=1,μ0=-1,則有:

(F(r,1)-F(r,-1))

(11)

3) 無限圓柱幾何微分升階形式

將(8)式代入(5)式,并令ξ1=1,ξ0=0,φ1=π,φ0=0,則有:

F(r,0,π)-F(r,1,0)+F(r,0,0),

(F(r,1,π)-F(r,0,π)-F(r,1,0)+F(r,0,0))

(12)

3 原函數(shù)確定解生成方法

根據(jù)微積分的性質(zhì),同一連續(xù)函數(shù)的積分原函數(shù)為一簇函數(shù)[13],若不確定特定的原函數(shù)進(jìn)行深度學(xué)習(xí),則深度學(xué)習(xí)方法難以收斂,以下分別就中子通量密度的角度Ω為一維與二維情況進(jìn)行討論,給出原函數(shù)確定解的方法。

3.1 Ω為一維變量

設(shè)μ0為任意固定角度變量,令Cn=-F1(r,μ0),則有特定的確定原函數(shù):

F0(r,μ)=F1(r,μ)-F1(r,μ0)

(13)

由此,原函數(shù)簇Fn(r,μ)在μ滿足定解約束條件下,找到了一個確定的原函數(shù)F0(r,μ)。

針對式(10),令μ0=-1,則Cn=-F1(x,-1)。根據(jù)式(13)有F0(x,-1)=0,將式(13)確定的原函數(shù)F0(x,μ)及相關(guān)關(guān)系代入式(10),得到升階形式平板幾何方程為:

(14)

此時,原函數(shù)定解約束條件是μ0=-1,F0(x,-1)=0。

同理,式(13)代入式(11),升階形式球幾何方程變?yōu)?

(15)

此時,原函數(shù)定解約束條件是μ0=-1,F0(r,-1)=0。

3.2 Ω為二維變量

設(shè)ξ0、φ0為任意固定角度變量,令Kn(r,ξ)=-F1(r,ξ,φ0),Jn(r,φ)=-F1(r,ξ0,φ),Mn(r)+Cn=F1(r,ξ0,φ0),則有特定的確定原函數(shù):

F0(r,ξ,φ)=F1(r,ξ,φ)-F1(r,ξ0,φ)-

F1(r,ξ,φ0)+F1(r,ξ0,φ0)

(16)

由此,就為原函數(shù)簇Fn(r,ξ,φ)在ξ、φ滿足定解約束條件下找到了一個確定的原函數(shù)F0(r,ξ,φ)。

針對式(12),令ξ0=0,φ0=0,則Kn(r,ξ)=-F1(r,ξ,0),Jn(r,φ)=-F1(r,0,φ),Mn(r)+Cn=F1(r,0,0)。根據(jù)式(16)有F0(r,0,φ)=0及其特例F0(r,0,π)=0,F0(r,ξ,0)=0及其特例F0(r,1,0)=0,F0(r,0,0)=0。將式(16)中確定的原函數(shù)F0(r,ξ,φ)及相關(guān)關(guān)系代入式(12),升階形式無限圓柱方程變?yōu)?

F0(r,1,π)

(17)

此時,原函數(shù)的定解約束項為:ξ0=0,F0(r,0,φ)=0;φ0=0,F0(r,ξ,0)=0;ξ0=0,φ0=0,F0(r,0,0)=0。

同理,根據(jù)圓柱坐標(biāo)下的二維有限圓柱幾何形式[1]給出相應(yīng)的升階表達(dá)方式:

(18)

式中,z為圓柱高度變量。原函數(shù)的定解約束項為:ξ0=-1,F0(r,z,-1,φ)=0;φ0=0,F0(r,z,ξ,0)=0;ξ0=-1,φ0=0,F0(r,z,-1,0)=0。

其他多維幾何方程可據(jù)此類推,定積分項的處理方法相同。

在式(9)基礎(chǔ)上,通過確定原函數(shù)表達(dá)的升階后輸運方程形式如下:

(19)

4 微分升階后的方程邊界條件表達(dá)方法

常見的中子輸運方程邊界條件包括:中子通量密度為有限非負(fù)實數(shù)、連續(xù)邊界條件、自由邊界條件、全反射邊界條件、反照邊界條件等[1-3]。原中子輸運方程微分升階后,由于中子輸運方程是由原函數(shù)表達(dá)的,需參照式(7)、(8)將原函數(shù)關(guān)于通量密度角度進(jìn)行求導(dǎo),變換為F(r,μ)′、Fξ,φ(r,ξ,φ)表示角通量密度函數(shù)φ(r,μ)、φ(r,ξ,φ),然后進(jìn)行原中子輸運方程的邊界表達(dá),部分邊界的形式列于表1。其他形式邊界可據(jù)此根據(jù)邊界類型與方程的幾何形式類推。

表1 邊界表達(dá)形式Table 1 Boundary representation

5 深度學(xué)習(xí)數(shù)值計算方法

基于PINN基本網(wǎng)絡(luò)架構(gòu)[5,8],結(jié)合升階后中子輸運方程原函數(shù)特點,進(jìn)行深度學(xué)習(xí)方法設(shè)計與實現(xiàn),其原理與流程是:將人工神經(jīng)網(wǎng)絡(luò)函數(shù)作為原函數(shù)F0(r,Ω)的試函數(shù),代入方程中得到的殘差形成控制方程損失函數(shù),與邊界條件約束、原函數(shù)定解約束、特征值約束形成的損失函數(shù),加權(quán)后得到統(tǒng)一的損失函數(shù);將神經(jīng)網(wǎng)絡(luò)神經(jīng)元連接權(quán)重及有效增殖因數(shù)keff作為可以調(diào)節(jié)的優(yōu)化參數(shù),以降低加權(quán)損失函數(shù)為目標(biāo)進(jìn)行深度學(xué)習(xí);當(dāng)損失函數(shù)趨近于極小值時,神經(jīng)網(wǎng)絡(luò)函數(shù)趨近于原函數(shù)的解。最后,對原函數(shù)關(guān)于角度變量求導(dǎo)降階,獲得原中子輸運方程角通量密度的數(shù)值解。

l層神經(jīng)網(wǎng)絡(luò)函數(shù)形式[15]為:

Nl(x)=fl(wlfl-1(wl-1fl-2(…f1(w1x+

b1)…+bl-2)+bl-1)+bl)

(20)

式中:w為神經(jīng)網(wǎng)絡(luò)連接權(quán)重;b為神經(jīng)網(wǎng)絡(luò)偏置項;Nl(x)為神經(jīng)網(wǎng)絡(luò)輸出;x為神經(jīng)網(wǎng)絡(luò)輸入向量;f為激活函數(shù)。

將式(20)作為試函數(shù),代入到控制方程、邊界條件約束、原函數(shù)定解約束、特征值約束中,生成相應(yīng)的損失函數(shù)。

1) 控制方程損失函數(shù)

令N(xr,xΩ)=F0(r,Ω),其中神經(jīng)網(wǎng)絡(luò)向量xr為方程空間幾何變量,xΩ為方程角通量密度角度變量,xr、xΩ可根據(jù)r、Ω維度做適應(yīng)性調(diào)整。

將式(20)及上述等價關(guān)系代入式(19),得到控制方程殘差表示的方程損失函數(shù)為:

(21)

式中,i為機(jī)器學(xué)習(xí)所需要的樣本點,可按照一定的概率密度分布生成,或根據(jù)方程的形式、邊界條件采用不同密度分布策略生成[8]。

2) 邊界條件約束損失函數(shù)

(22)

3) 原函數(shù)定解約束損失函數(shù)

設(shè)φ(r,Ω)定義域中,rd,Ωd∈Гd為原函數(shù)的定解約束取值域,F0(rd,Ωd)為第3節(jié)所示的定解約束函數(shù)項,令式(20)為N(xrd,xΩd)=F0(rd,Ωd),則得到原函數(shù)定解約束項損失函數(shù)為:

(23)

式中,m為定解約束項取值域上的機(jī)器學(xué)習(xí)樣本點。

4) 特征值約束損失函數(shù)

如文獻(xiàn)[1-3]所述,由于穩(wěn)態(tài)中子輸運特征值方程相同特征值的解有無窮多個,直接用深度學(xué)習(xí)方法求解在收斂速度上存在一定困難。解決的方法是通過增加特定固定的幾何空間點,設(shè)定歸一化的通量密度預(yù)設(shè)值,構(gòu)建相應(yīng)的特征值約束損失函數(shù)來提高收斂速度[8]。固定的幾何空間點一般選擇通量密度的對稱點、中心點等,預(yù)設(shè)固定值理論上可以是任意實數(shù)值,取值的大小會影響最終通量密度分布的幅度,但不會影響其分布形狀。實際應(yīng)用中,可根據(jù)堆芯功率設(shè)定等應(yīng)用需要,類似傳統(tǒng)求解輸運方程的方法,對最終結(jié)果做歸一化處理,保證計算結(jié)果的唯一性。

(24)

式中,k為預(yù)設(shè)值特定樣本點。

5) 統(tǒng)一加權(quán)損失函數(shù)

將上述損失函數(shù)加權(quán),得到統(tǒng)一的機(jī)器學(xué)習(xí)損失函數(shù):

Fall-Loss=PfFf-loss+PbFb-loss+

PdFd-loss+PcFc-loss

(25)

式中,Pf、Pb、Pd、Pc分別為控制方程、邊界條件約束、原函數(shù)定解約束、特征值約束損失函數(shù)權(quán)重。

深度學(xué)習(xí)過程是調(diào)整式(20)中的wl、bl,若為非臨界問題同時調(diào)節(jié)keff,使得Fall-Loss趨近于0,N(xr,xΩ)趨近于F0(r,Ω)。此后,將N(xr,xΩ)關(guān)于角度Ω進(jìn)行求導(dǎo)降階,得到的N′(xr,xΩ)即原中子輸運方程的角通量密度數(shù)值解φ(r,Ω),此時的keff也趨近于系統(tǒng)有效增殖因數(shù)的解。實踐中一般為Fall-Loss設(shè)定一定的收斂準(zhǔn)則,機(jī)器學(xué)習(xí)滿足收斂準(zhǔn)則后停止學(xué)習(xí)迭代,N′(xr,xΩ)求導(dǎo)采用機(jī)器學(xué)習(xí)自動微分技術(shù)實現(xiàn)[16]。

6 數(shù)值驗證

6.1 臨界單群平板幾何輸運問題

針對式(3)描述的平板幾何單群輸運問題,設(shè)平板材料特性Σt=0.050 cm-1,Σs=0.030 cm-1,νΣf=0.022 5 cm-1,平板兩側(cè)為真空邊界條件。按文獻(xiàn)[2]理論,本題的截面數(shù)據(jù)相當(dāng)于C=(Σs+νΣf)/Σt=1.05(C值定義參考文獻(xiàn)[2]),根據(jù)文獻(xiàn)[2]給出的C=1.05時的幾何平板臨界半厚度,本問題臨界半厚度b=3.300 263 772λ=66.005 275 44 cm,其中λ=1/Σt為中子自由程,此時keff=1。利用升階形式方程(式(14))進(jìn)行深度機(jī)器學(xué)習(xí)求解。

表2 平板幾何機(jī)器學(xué)習(xí)損失函數(shù)及樣本生成方式Table 2 Machine learning loss function and sample generation method of slab geometry

a——升階原函數(shù)分布散點圖;b——原函數(shù)降階后的角通量密度分布散點圖圖1 平板幾何單群輸運問題數(shù)值計算結(jié)果Fig.1 Numerical calculation result for single group transport problem of slab geometry

驗證實驗中,樣本數(shù)量大一般具有相對較高精度,但學(xué)習(xí)時間較長,反之樣本數(shù)量少,學(xué)習(xí)效率高,但精度下降。為了采用較為少量樣本點獲得較高的機(jī)器學(xué)習(xí)精度,在控制方程邊界附近采用了增加樣本點的技術(shù)[8],實踐中可根據(jù)性能要求及硬件環(huán)境綜合考慮樣本數(shù)量與方式;同時,權(quán)重系數(shù)、Pf、Pb、Pd、Pc為經(jīng)驗參數(shù),本文的策略是以Pf=1為基準(zhǔn),在訓(xùn)練初期設(shè)定不同的Pb、Pd、Pc,選取使得Fall-Loss下降較快的值進(jìn)行后續(xù)訓(xùn)練。由于神經(jīng)網(wǎng)絡(luò)具有一定的不可解釋性,對于不同的問題一般有不同的優(yōu)化值,主要基于實驗過程進(jìn)行選擇。

圖1示出平板幾何單群輸運問題數(shù)值計算結(jié)果,圖2示出平板幾何下中心歸一化定積分項標(biāo)通量密度散點圖。平板幾何標(biāo)通量計算結(jié)果與理論值對比列于表3。

如圖1、2和表3所示,深度學(xué)習(xí)技術(shù)求解出標(biāo)通量密度與理論值對比具有良好的精度。同時,由于機(jī)器學(xué)習(xí)獲得通量密度分布是通過式(19)形式表達(dá)的,而式(19)及其導(dǎo)數(shù)形式均為連續(xù)函數(shù),通過神經(jīng)網(wǎng)絡(luò)的泛化計算可給出角通量密度的連續(xù)分布,這相對于SN、MOC等傳統(tǒng)方法具有顯著的特點。

圖2 平板幾何下中心歸一化定積分項標(biāo)通量密度散點圖Fig.2 Scatter plot of normalized definite integral term scalar flux density under slab geometry

6.2 臨界單群球幾何輸運問題

針對式(4)描述的球幾何單群單材料區(qū)域問題,其材料特性、C與邊界條件同6.1節(jié),由文獻(xiàn)[2]理論可知C=1.05時的歸一化球幾何半徑R=7.277 181 79λ=145.543 635 8 cm,此時系統(tǒng)為臨界時狀態(tài),keff=1。利用升階形式方程(式(15))進(jìn)行深度機(jī)器學(xué)習(xí)求解。

表3 平板幾何標(biāo)通量密度計算結(jié)果與理論值對比Table 3 Comparison between calculated result and theoretical value of slab geometric scalar flux density

圖3示出球幾何單群輸運問題數(shù)值計算結(jié)果。圖4示出球幾何下中心歸一化定積分項標(biāo)通量密度散點圖。臨界球幾何標(biāo)通量計算結(jié)果與理論值對比列于表5。可見,數(shù)值計算結(jié)果驗證了機(jī)器學(xué)習(xí)方法具有良好的精度,同時也給出了角通量密度的連續(xù)分布。

a——升階原函數(shù)分布散點圖;b——原函數(shù)降階后的角通量密度分布散點圖圖3 球幾何單群輸運問題數(shù)值計算結(jié)果Fig.3 Numerical calculation result for single group transport problem of spherical geometry

圖4 球幾何下中心歸一化定積分項標(biāo)通量密度散點圖Fig.4 Scatter plot of normalized definite integral term scalar flux density under spherical geometry

6.3 非臨界單群球幾何輸運問題

針對6.2節(jié)中球幾何單群單材料區(qū)域問題,將式(4)中的Σf引入誤差Δ,使得Σ′f=Σf·Δ,此時反應(yīng)堆為非臨界系統(tǒng)。當(dāng)keff=Δ時式(4)成立,將此作為數(shù)值計算參考標(biāo)準(zhǔn)理論值,表5中的歸一化理論值仍然適用。機(jī)器學(xué)習(xí)方法中交替更新神經(jīng)網(wǎng)絡(luò)權(quán)重與keff,其余網(wǎng)絡(luò)參數(shù)、損失函數(shù)及樣本生成方式同6.2節(jié)。

非臨界球幾何標(biāo)通量密度計算結(jié)果與理論值對比列于表6,其升階原函數(shù)與降階后的通量密度形狀與圖3、4相似。

本示例可作為評估深度學(xué)習(xí)求解非臨界系統(tǒng)keff誤差的標(biāo)定方法。

6.4 非臨界雙群平板幾何輸運問題

針對平板幾何雙群單材料區(qū)域問題,雙群平板理論及輸運方程形式參考文獻(xiàn)[1-3],定積分項處理方法如式(14),幾何參數(shù)同6.1節(jié),材料參數(shù)列于表7。

表5 臨界球幾何標(biāo)通量密度計算結(jié)果與理論值對比Table 5 Comparison between calculated result and theoretical value of critical scalar flux density for spherical geometry

機(jī)器學(xué)習(xí)方法與參數(shù)選擇為:用兩個神經(jīng)網(wǎng)絡(luò)分別表示熱群與快群,機(jī)器學(xué)習(xí)過程中,交替更新兩個神經(jīng)網(wǎng)絡(luò)權(quán)重、keff。特征值約束針對快群為:F′0(0, - 1) =F′0(0,1)=0.1,不加入熱群特征值約束,即熱群Pc=0,其余網(wǎng)絡(luò)參數(shù)、損失函數(shù)及樣本生成方式同6.1節(jié),keff更新方法同6.3節(jié)。

由于此問題本身沒有解析解或理論估計值,將深度學(xué)習(xí)結(jié)果與輸運計算程序OpenMC[17]獲得的標(biāo)通量密度計算結(jié)果進(jìn)行比較(均以快群的x=0處標(biāo)通量密度進(jìn)行歸一化),結(jié)果如圖5、表8所示。由圖5可見,數(shù)值計算結(jié)果也可獲得雙群系統(tǒng)連續(xù)性角通量密度分布。

a——快群通量密度計算結(jié)果;b——熱群通量密度計算結(jié)果;c——快群標(biāo)通量密度與OpenMC軟件對比圖;d——熱群標(biāo)通量密度與OpenMC軟件對比圖圖5 非臨界雙群平板幾何輸運問題數(shù)值計算結(jié)果Fig.5 Numerical result of transport problem of non-critical two-group slab geometry

表8 非臨界雙群平板幾何標(biāo)通量密度計算結(jié)果Table 8 Calculation result for non-critical two-group scalar flux density of slab geometry

6.5 非臨界單群多材料區(qū)域平板幾何輸運問題

非臨界單群多材料區(qū)域平板幾何如圖6所示。

平板幾何由兩種材料3個區(qū)域組成,b定義與6.1節(jié)相同,材料參數(shù)列于表9。

機(jī)器學(xué)習(xí)方法、網(wǎng)絡(luò)參數(shù)、損失函數(shù)及樣本生成方式同6.1節(jié),keff更新方法同6.3節(jié)。將計算結(jié)果與OpenMC[17]結(jié)果進(jìn)行比較(按照x=0處的標(biāo)通量密度進(jìn)行歸一化處理),結(jié)果如圖7、表10所示。

與前述例題類似,圖7a顯示出角通量密度連續(xù)分布的特點。對照單材料臨界系統(tǒng)(圖2),機(jī)器學(xué)習(xí)結(jié)果(圖7b)的標(biāo)通量密度分布兩側(cè)存在明顯的“內(nèi)陷”,與OpenMC結(jié)果呈現(xiàn)出的形狀相比具有良好的一致性。

圖6 非臨界單群多材料問題幾何描述Fig.6 Geometric diagram for non-critical single group multi-material problem

表9 非臨界單群多材料區(qū)域平板幾何區(qū)域材料特性Table 9 Material property of non-critical single group multi-material region of slab geometry

a——角通量密度計算結(jié)果;b——標(biāo)通量密度計算結(jié)果與OpenMC結(jié)果對比圖圖7 非臨界單群多材料區(qū)域平板標(biāo)通量密度數(shù)值計算結(jié)果Fig.7 Numerical calculation result of non-critical single group multi-material region scalar flux density for slab geometry

表10 非臨界單群平板幾何標(biāo)通量密度計算結(jié)果Table 10 Calculation result of non-critical single group scalar flux density of slab geometry

綜上所述,上述驗證表明微分升階理論數(shù)理方法基本原理在臨界與非臨界系統(tǒng)下均具有良好的適用性。由于非臨界系統(tǒng)keff計算值存在一定誤差,盡管這種誤差范圍與PINN基礎(chǔ)框架[5]中的示例相當(dāng),但這是導(dǎo)致非臨界系統(tǒng)的角通量密度精度低于臨界系統(tǒng)的重要因素。

7 結(jié)論

本文針對當(dāng)前深度學(xué)習(xí)技術(shù)求解微分方程的方法不能直接適用于具有積分項的中子輸運方程面臨的困難,提出了將輸運方程通過數(shù)理方法升階為高階微分形式的原函數(shù)進(jìn)行機(jī)器學(xué)習(xí),成功后再將其降階為原方程數(shù)值解的微分變階理論。同時,給出了相應(yīng)的原函數(shù)定解條件約束、邊界條件約束、特征值約束形式等損失函數(shù)構(gòu)造方法。最后,分別通過臨界與非臨界平板、球幾何例題,驗證了該理論及相關(guān)方法的正確性。研究表明機(jī)器學(xué)習(xí)方法求解得到數(shù)值結(jié)果具有幾何-角通量密度分布連續(xù)性的特點。研究工作為中子輸運方程的數(shù)值求解方法探索出了新的技術(shù)途徑,并為利用深度學(xué)習(xí)方法求解幾何、能群更加復(fù)雜的中子輸運方程提供了良好的技術(shù)支撐。