劉若晴,馬 濤,張立國(guó)
(清華大學(xué) 核能與新能源技術(shù)研究院,先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心,先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)
作為第4代先進(jìn)核能系統(tǒng)的典型代表,球床式高溫氣冷堆(簡(jiǎn)稱球床堆)具有國(guó)際公認(rèn)的良好安全特性[1-2]。隨著球床式高溫氣冷堆的商業(yè)開發(fā),核安保與核不擴(kuò)散問(wèn)題備受關(guān)注[3-4]。對(duì)球床式高溫氣冷堆的乏燃料進(jìn)行超鈾含量的準(zhǔn)確計(jì)算,是核安全保障的必然要求[5-6]。
10 MW高溫氣冷實(shí)驗(yàn)堆(HTR-10)是由清華大學(xué)核能與新能源技術(shù)研究院研發(fā)的、主要用于實(shí)驗(yàn)的球床式高溫氣冷堆[7-8]。與傳統(tǒng)反應(yīng)堆相比,該實(shí)驗(yàn)堆的核燃料以大量燃料球的形式依靠重力作用在反應(yīng)堆內(nèi)流動(dòng),具有不停堆換料的特點(diǎn)[9]。燃料球穿越堆底后,由在線燃耗測(cè)量系統(tǒng)對(duì)其燃耗、部分活度容易測(cè)量的核素活度(簡(jiǎn)稱易測(cè)核素活度)進(jìn)行測(cè)量[10]。燃料球的燃耗達(dá)到卸料燃耗整定值將會(huì)被卸出堆芯并成為乏燃料球,否則返回堆芯再次循環(huán)[11]。
對(duì)于輕水反應(yīng)堆,核燃料在反應(yīng)堆內(nèi)的位置可跟蹤標(biāo)記,其進(jìn)堆時(shí)間、出堆時(shí)間可準(zhǔn)確掌握,卸出堆芯成為乏燃料后,結(jié)合反應(yīng)堆運(yùn)行歷史,可用核素生成與耗減計(jì)算軟件算出乏燃料的燃耗和各種核素的活度,即可根據(jù)初態(tài)和過(guò)程信息,利用軟件正向計(jì)算末態(tài)信息[12-14]。對(duì)于球床式高溫氣冷堆,針對(duì)特定的乏燃料球,可用在線燃耗測(cè)量系統(tǒng)測(cè)量其末態(tài)的易測(cè)核素活度,而無(wú)法掌握其進(jìn)堆時(shí)間、在反應(yīng)堆內(nèi)的流動(dòng)情況和循環(huán)次數(shù)等信息,故無(wú)法直接應(yīng)用核素生成與耗減計(jì)算軟件計(jì)算其燃耗和各種核素活度。本工作的難點(diǎn)在于,需根據(jù)一部分已知末態(tài)信息(易測(cè)核素活度)計(jì)算另一部分未知末態(tài)信息(不易測(cè)核素活度),而根據(jù)已知末態(tài)信息難以反推出初態(tài)和過(guò)程信息,故依賴于初態(tài)和過(guò)程信息進(jìn)行正向計(jì)算來(lái)得到未知末態(tài)信息的路徑不可行。因而直接通過(guò)計(jì)算并結(jié)合少量試驗(yàn)分析的方法確定輻照后燃料中超鈾核素含量的方法并不適用。
考慮到核反應(yīng)過(guò)程中各種核素的生成和耗減彼此關(guān)聯(lián),本文使用神經(jīng)網(wǎng)絡(luò)方法,針對(duì)在線燃耗測(cè)量系統(tǒng)中易測(cè)核素活度與不易測(cè)核素活度間的對(duì)應(yīng)關(guān)系進(jìn)行建模,通過(guò)易測(cè)核素活度估算輻照后球形燃料元件中超鈾元素含量。
在反應(yīng)堆的平衡堆芯中,各種核素的生成和耗減滿足以下方程[15]:
(λp+φσp)Xpp=1,…,N
(1)
其中:λq為第q種核素的衰變常量,s-1;σp為第p種核素譜平均中子吸收截面,cm2;lpq為第q種核素衰變成第p種核素的份額;fpr為第r種核素吸收中子后轉(zhuǎn)變?yōu)榈趐種核素的份額;φ為按位置、能量平均的中子注量率,cm-2·s-1;Xq、Xr和Xp分別為第q、r和p種核素的含量;t為時(shí)間。
式(1)中,λq、σp、lpq、fpr與反應(yīng)堆運(yùn)行歷史無(wú)關(guān),φ與反應(yīng)堆運(yùn)行歷史相關(guān)。文獻(xiàn)[16]采用VSOP程序結(jié)合核測(cè)數(shù)據(jù)和校正的反應(yīng)堆功率,追蹤評(píng)價(jià)HTR-10的運(yùn)行過(guò)程,給出不同時(shí)段堆芯中子注量率分布。
一方面,對(duì)于HTR-10,堆芯連續(xù)地裝載新元件,同時(shí)卸出乏燃料元件,最終長(zhǎng)期處于穩(wěn)定的平衡堆芯狀態(tài)[17],不同核素的生成和耗減滿足該方程,具有彼此關(guān)聯(lián)的特性。另一方面,該微分方程的復(fù)雜性和耦合性決定了以乏燃料球中的易測(cè)核素活度為已知量,難以通過(guò)直接求逆問(wèn)題計(jì)算不易測(cè)核素的含量。復(fù)雜運(yùn)行歷史條件下,特別是含有長(zhǎng)時(shí)間停堆時(shí)段,直接針對(duì)式(1)求逆將面臨不適定問(wèn)題[18]。
由于未知量與已知量的關(guān)聯(lián)性,以及問(wèn)題直接求解的復(fù)雜性,本文采用深度神經(jīng)網(wǎng)絡(luò)對(duì)易測(cè)核素活度與不易測(cè)核素活度之間的關(guān)系進(jìn)行正向建模[19]。深度神經(jīng)網(wǎng)絡(luò)作為一種多層前饋神經(jīng)網(wǎng)絡(luò)[20],采用誤差反向傳播算法[21],能在反復(fù)迭代過(guò)程中調(diào)整網(wǎng)絡(luò)權(quán)重和偏置,對(duì)于存在確定內(nèi)部機(jī)制的問(wèn)題具有強(qiáng)大的擬合能力,在求解復(fù)雜問(wèn)題上具有明顯優(yōu)勢(shì)。
文獻(xiàn)[22]給出HTR-10在線燃耗測(cè)量系統(tǒng)中較易測(cè)量的核素,共有11種,包括95Zr、106Rh、137Cs、103Ru、140Ba、140La、134Cs、95Nbm、133Xem、125Sb和131Tem,其計(jì)數(shù)率分別為101.6、30.27、21.94、20.83、10.2、10.04、9.2、8.87、8.87、5.31和5.26 s-1。包括超鈾元素同位素在內(nèi)的其余核素活度不易測(cè)量。因此,在深度神經(jīng)網(wǎng)絡(luò)訓(xùn)練過(guò)程中無(wú)法直接采用不易測(cè)核素活度的實(shí)際測(cè)量值。受此技術(shù)條件所限,本文神經(jīng)網(wǎng)絡(luò)訓(xùn)練數(shù)據(jù)采用軟件模擬值。
文獻(xiàn)[23]應(yīng)用球流模擬軟件LAMMPS模擬了HTR-10大量燃料球在堆芯內(nèi)的流動(dòng),結(jié)合HTR-10實(shí)際運(yùn)行功率歷史,應(yīng)用軟件KORIGEN模擬了燃料球進(jìn)入堆芯后的核素生成與耗減過(guò)程,計(jì)算出乏燃料球的燃耗及75種核素的活度,其中包括11種易測(cè)核素活度。
多種參數(shù)計(jì)算值與實(shí)驗(yàn)值的對(duì)照表明,軟件模擬結(jié)果能較好地反映燃料球在HTR-10堆芯內(nèi)實(shí)際流動(dòng)和燃耗情況,軟件模擬結(jié)果與實(shí)際物理過(guò)程的吻合程度較高[23]。因此,本文以LAMMPS和KORIGEN模擬的大量燃料球的燃耗及各種核素活度為數(shù)據(jù)基礎(chǔ)建立神經(jīng)網(wǎng)絡(luò)。其中11種易測(cè)核素活度作為神經(jīng)網(wǎng)絡(luò)的輸入端數(shù)據(jù),燃耗和其余64種不易測(cè)核素活度作為神經(jīng)網(wǎng)絡(luò)的輸出端數(shù)據(jù)。每個(gè)燃料球的燃耗和75種核素活度為1組數(shù)據(jù),每組數(shù)據(jù)包含11個(gè)輸入數(shù)據(jù)和65個(gè)輸出數(shù)據(jù)。
在LAMMPS和KORIGEN中,導(dǎo)入HTR-10 2003年初至2007年5月中旬的實(shí)際運(yùn)行功率歷史數(shù)據(jù)以為KORIGEN進(jìn)行核反應(yīng)模擬提供實(shí)時(shí)功率數(shù)據(jù),令20 000個(gè)模擬的燃料球先后進(jìn)入堆芯,并對(duì)球流和核反應(yīng)過(guò)程進(jìn)行模擬,得出了20 000組模擬數(shù)據(jù)。其中19 500組模擬數(shù)據(jù)作為訓(xùn)練組,用于深度神經(jīng)網(wǎng)絡(luò)的訓(xùn)練;其余500組模擬數(shù)據(jù)作為校驗(yàn)組,用于驗(yàn)證該深度神經(jīng)網(wǎng)絡(luò)的性能。
由于不同核素活度數(shù)量級(jí)相差較大,為避免數(shù)值計(jì)算中病態(tài)方程組的發(fā)散問(wèn)題[24],用式(2)對(duì)20 000個(gè)燃料球的燃耗和75種核素活度進(jìn)行歸一化處理:
i=1,…,20 000,j=1,…,76
(2)
經(jīng)歸一化處理,燃耗和所有核素的活度變?yōu)榫禐?、方差為1的分布,有利于深度神經(jīng)網(wǎng)絡(luò)初值賦值和最終結(jié)果收斂。
深度神經(jīng)網(wǎng)絡(luò)包括輸入層、隱含層和輸出層[25]。輸入層節(jié)點(diǎn)數(shù)與輸入數(shù)據(jù)個(gè)數(shù)一致,為11;輸出層節(jié)點(diǎn)數(shù)與輸出數(shù)據(jù)個(gè)數(shù)一致,為65。經(jīng)反復(fù)摸索調(diào)整,將隱含層的層數(shù)設(shè)置為8,每層節(jié)點(diǎn)數(shù)分別為100、300、500、500、500、300、100和80。
在梯度下降方法的選擇上,采用批量隨機(jī)梯度下降方法[26]。經(jīng)探索優(yōu)化,每批量的數(shù)據(jù)組數(shù)設(shè)置為20。若每批量的數(shù)據(jù)組數(shù)較大,有利于早期快速收斂,但其缺點(diǎn)在于精度不高,且收斂接近完成時(shí)波動(dòng)較大。
學(xué)習(xí)效率設(shè)置為0.001。研究發(fā)現(xiàn),學(xué)習(xí)效率設(shè)置過(guò)大,將導(dǎo)致深度神經(jīng)網(wǎng)絡(luò)訓(xùn)練時(shí)無(wú)法收斂,設(shè)置過(guò)小則不利于提高收斂速度。
模型訓(xùn)練中的懲罰函數(shù)[26]采用下式:
(3)
深度神經(jīng)網(wǎng)絡(luò)的權(quán)重和偏置,通過(guò)誤差反向傳播算法進(jìn)行更新。深度神經(jīng)網(wǎng)絡(luò)基于目前業(yè)界流行的、支持GPU加速的Tensorflow框架,使用Python編程語(yǔ)言進(jìn)行構(gòu)建[27-29]。
將19 500組模擬數(shù)據(jù)導(dǎo)入深度神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練,發(fā)現(xiàn)收斂曲線(圖1)在初期下降十分迅速,后期下降緩慢,這一結(jié)果符合批量隨機(jī)梯度下降方法的特征。
圖1 深度神經(jīng)網(wǎng)絡(luò)的收斂曲線Fig.1 Convergence curve of deep neural network
訓(xùn)練完成后,將校驗(yàn)組數(shù)據(jù)導(dǎo)入所生成的深度神經(jīng)網(wǎng)絡(luò)模型,每組數(shù)據(jù)中的11種易測(cè)核素活度作為模型的輸入數(shù)據(jù),經(jīng)計(jì)算,模型會(huì)輸出燃耗和不易測(cè)超鈾元素同位素的活度預(yù)測(cè)值。
圖2 239Pu活度的深度神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值與校驗(yàn)組原始值的對(duì)比Fig.2 239Pu activity comparison between prediction value by deep neural network and original value
將深度神經(jīng)網(wǎng)絡(luò)模型的預(yù)測(cè)值與校驗(yàn)組數(shù)據(jù)中相應(yīng)的原始值進(jìn)行比較,發(fā)現(xiàn)二者吻合程度較好。以239Pu核素為例,如圖2所示,模型預(yù)測(cè)數(shù)據(jù)與其對(duì)應(yīng)的測(cè)試數(shù)據(jù)基本分布在y=x直線附近,模型對(duì)于239Pu的活度預(yù)測(cè)比較準(zhǔn)確。
燃耗預(yù)測(cè)如圖3所示,模型預(yù)測(cè)的燃耗與對(duì)應(yīng)的原始數(shù)據(jù)中燃耗組成的點(diǎn)也基本分布在y=x直線附近。因此,模型對(duì)于燃耗的預(yù)測(cè)也比較準(zhǔn)確。
圖3 燃耗的深度神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值與校驗(yàn)組原始值的對(duì)比Fig.3 Burnup comparison between prediction value by deep neural network and original value
本文后續(xù)計(jì)算中的誤差均指神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值與LAMMPS和KORIGEN軟件模擬值之間的誤差,而非神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值與該數(shù)據(jù)實(shí)測(cè)值間的誤差。
使用下式求各超鈾元素同位素活度的神經(jīng)網(wǎng)絡(luò)平均相對(duì)誤差:
(4)
使用下式求超鈾元素同位素活度的相關(guān)系數(shù)R2:
(5)
(6)
燃料球燃耗和239Pu活度的相對(duì)誤差如圖4所示。
對(duì)燃料球燃耗和超鈾元素活度的部分統(tǒng)計(jì)學(xué)指標(biāo)進(jìn)行計(jì)算,計(jì)算結(jié)果列于表1。其中,相關(guān)系數(shù)R2基本在1.0附近,燃耗、236U活度、239Pu活度相對(duì)誤差在±5%以內(nèi)的燃料球份額超過(guò)99%,相對(duì)誤差在±2.5%以內(nèi)的燃料球份額超過(guò)81.6%??煽闯?,深度神經(jīng)網(wǎng)絡(luò)模型對(duì)軟件模擬數(shù)據(jù)的預(yù)測(cè)精度總體上較好。
圖4 燃料球燃耗和239Pu活度的相對(duì)誤差Fig.4 Relative error of burnup and 239Pu activity
表1 超鈾元素同位素活度預(yù)測(cè)準(zhǔn)確性的統(tǒng)計(jì)Table 1 Statistics of model predicted transuranic isotope activity
本文基于深度神經(jīng)網(wǎng)絡(luò),使用HTR-10 2003年初至2007年5月中旬的運(yùn)行功率歷史數(shù)據(jù)作為KORIGEN核反應(yīng)模擬的功率數(shù)據(jù)進(jìn)行超鈾元素同位素含量計(jì)算。用11種易測(cè)核素的活度推算出的各超鈾元素同位素活度與校驗(yàn)組中對(duì)應(yīng)值符合程度很好,如校驗(yàn)組數(shù)據(jù)中239Pu預(yù)測(cè)結(jié)果相對(duì)誤差在±2.5%以內(nèi)的份額為81.6%,預(yù)測(cè)精度較高。該結(jié)果說(shuō)明本文提出的方法處理球床堆燃料中核素含量估計(jì)問(wèn)題時(shí)具有較好的發(fā)展?jié)摿ΑT摲椒ㄒ灿芯窒扌?。首先,不易測(cè)核素活度缺乏實(shí)測(cè)數(shù)據(jù),該方法采用軟件模擬的數(shù)據(jù)作為基礎(chǔ),除深度神經(jīng)網(wǎng)絡(luò)模型對(duì)軟件模擬數(shù)據(jù)的擬合不到位造成的誤差外,還可能疊加軟件對(duì)實(shí)際物理過(guò)程模擬不到位造成的誤差;其次,易測(cè)核素中不同核素活度的測(cè)量難易程度和精度存在差別,使用較多的易測(cè)核素活度作為輸入端數(shù)據(jù),增加了測(cè)量工作量,降低了測(cè)量精度,不利于模型的實(shí)際應(yīng)用。因此,改進(jìn)LAMMPS和KORIGEN等軟件對(duì)實(shí)際物理過(guò)程的模擬精度、改進(jìn)神經(jīng)網(wǎng)絡(luò)模型對(duì)軟件模擬數(shù)據(jù)的擬合精度、以盡可能少的易測(cè)量核素活度實(shí)現(xiàn)對(duì)超鈾元素同位素活度的估算,將是下一步的研究目標(biāo)。