薛志恒，劉志宏，趙 晶，佘 頂

(清華大學(xué) 核能與新能源技術(shù)研究院，先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心，先進(jìn)反應(yīng)堆工程與安全教育部重點實驗室，北京 100084)

我國在役運行的第2代壓水堆，其在線監(jiān)測功能由核儀表系統(tǒng)[1]外核儀表連續(xù)監(jiān)測反應(yīng)堆功率變化及功率分布，監(jiān)測反應(yīng)堆徑向功率傾斜和軸向功率偏差等。目前，使用堆外探測數(shù)據(jù)進(jìn)行堆芯功率精細(xì)重構(gòu)的方法主要有傳輸矩陣法、諧波綜合法及神經(jīng)網(wǎng)絡(luò)法。

傳輸矩陣法認(rèn)為堆外探測器信號與功率分布之間存在線性關(guān)系，堆芯功率可由探測器示數(shù)矩陣計算得到。這種做法可重構(gòu)出堆芯軸向功率分布，但精確度尚有提升空間。清華大學(xué)核能與新能源技術(shù)研究院李富等[2]提出了諧波綜合法。諧波綜合法建立了各因素對功率分布變化影響的特征分布，但未考慮工況變化下截面的變化，可能導(dǎo)致無法準(zhǔn)確模擬真實通量。夏虹等[3]提出利用RBF神經(jīng)網(wǎng)絡(luò)實現(xiàn)對堆芯功率三維分布的實時監(jiān)測，其有統(tǒng)計數(shù)學(xué)上的意義，但忽略了擴(kuò)散方程以及有效中子增殖因數(shù)等物理約束。

本文提出一種基于堆外探測實現(xiàn)堆芯功率重構(gòu)的方法，并通過算例來驗證該方法的可用性。本方法對引起反應(yīng)堆堆芯功率波動的影響因素進(jìn)行建模，建立影響因素的擾動與截面變化、堆外探測器示數(shù)之間的關(guān)系，利用特征統(tǒng)計算法(CSA)搜索擾動模型參數(shù)并通過擴(kuò)散方程驗證參數(shù)的準(zhǔn)確性。

1 堆芯功率分布重構(gòu)方法

賈健等[4]在研究基于堆內(nèi)監(jiān)測系統(tǒng)的堆芯功率重構(gòu)方法時，提出了一種將功率偏差歸結(jié)于組件等效截面參數(shù)偏差的功率分布重構(gòu)方法。該方法的實質(zhì)是將功率分布重構(gòu)問題轉(zhuǎn)化為以截面參數(shù)為自變量、堆芯功率分布為目標(biāo)函數(shù)的優(yōu)化搜索問題。本文將上述方法作進(jìn)一步的改進(jìn)，基于某供熱堆堆型探討利用堆外探測器提供的數(shù)據(jù)來進(jìn)行堆芯功率分布的重構(gòu)方法可行性。

該方法首先分析反應(yīng)堆堆芯功率波動的原因，并對各影響因素進(jìn)行建模；建立影響因素的擾動與截面變化、堆外探測器示數(shù)之間的關(guān)系。當(dāng)堆芯功率波動時，調(diào)整影響因素模型參數(shù)，由上述關(guān)系得到堆外探測器計數(shù)，與堆外探測器的實際計數(shù)進(jìn)行比較，如果不滿足偏差要求，則重新調(diào)整影響因素的模型參數(shù)，直到滿足偏差要求。

堆外測量系統(tǒng)需由堆內(nèi)測量系統(tǒng)定期校正。本文假定堆外測量系統(tǒng)已由堆內(nèi)測量系統(tǒng)定期校正過，在分析功率分布變化原因時，可忽略富集度偏差、設(shè)備加工精度等因素，只分析有限的幾個因素，如溫度場波動、氙振蕩和控制棒移動等。本文僅分析氙毒模型和溫度場波動模型。

氙濃度的變化會導(dǎo)致堆芯局部中子注量率的變化[5]?；谀Ｐ蛥?shù)個數(shù)的要求和氙振蕩的實際現(xiàn)象，給出下式：

ΔXe(x,y,z)=(a1x+a2y+a3z)·

exp(-a4x2-a5y2-a6z2)

(1)

式中：ΔXe(x,y,z)為氙毒在空間(x,y,z)點處的變化量；a1、a2、a3、a4、a5、a6為系數(shù)，是未知量。

堆芯冷卻劑流經(jīng)堆芯帶走核燃料釋放出的能量，溫度場如果出現(xiàn)異常波動，可能的原因有燃耗變化、控制棒移動、環(huán)路流量失衡、堆芯入口溫度變化及入口流量分配、包殼破裂等原因。借助諧波綜合法的思想，使用特征分布的組合來表達(dá)溫度場波動情況。本文為說明問題，用下式表達(dá)：

ΔT(x,y,z)=b1x+b2y+

b3z+b4xy+b5xz+b6yz

(2)

式中，b1～b6為溫度擾動模型的系數(shù)。

進(jìn)行擴(kuò)散方程求解前需得到核截面數(shù)據(jù)。求解擾動變化量與截面變化之間的關(guān)系，常見的方法有線性插值法、切比雪夫多項式插值[6]等。本文采用線性插值法，公式如下：

ΔΣ1=c1ΔXe

(3)

ΔΣ2=c2ΔT

(4)

式中：ΔΣ為截面參數(shù)的實際變化量，cm-1；ΔXe為空間某位置氙毒變化百分比，100%；ΔT為空間某位置溫度變化值，℃；c1、c2為常數(shù)系數(shù)，可通過計算得到。

2 優(yōu)化問題的求解

本文選用特征統(tǒng)計算法[7]作為全局優(yōu)化算法。該算法的基本流程如下：1) 在定義域內(nèi)隨機(jī)給出若干組解N0，并計算各自的目標(biāo)函數(shù)值；2) 將得到的目標(biāo)函數(shù)值排序，選取其中最好的N1個目標(biāo)函數(shù)值，并選出其對應(yīng)的N1組解；3) 描繪出每個特征量在定義域上的分布曲線；4) 根據(jù)分布曲線繼續(xù)選取N0組解，計算出各自的目標(biāo)函數(shù)值后，在N0+N1中選取最好的N1個目標(biāo)函數(shù)值及其對應(yīng)的解；5) 重復(fù)第3步及第4步，直至循環(huán)結(jié)束。

在Qt平臺上使用C++語言對特征統(tǒng)計算法進(jìn)行編程。本文涉及到的目標(biāo)函數(shù)值是由擾動多項式的系數(shù)作為特征量組成的，目標(biāo)函數(shù)為：

(5)

式中：Ni為第i組解；n為探測器能提供的數(shù)據(jù)個數(shù)；Pk0為第k個探測器提供的數(shù)據(jù)；Pki為通過計算得到的數(shù)據(jù)。F(Ni)值越小，表示重構(gòu)的精確度越高。

3 計算測試模型

將堆芯各處的裂變中子對堆外探測器讀數(shù)的貢獻(xiàn)率作為堆外探測器響應(yīng)函數(shù)，能準(zhǔn)確描述堆內(nèi)功率分布和堆外探測器讀數(shù)的關(guān)系。實際的反應(yīng)堆響應(yīng)矩陣計算非常復(fù)雜。本文為了說明問題，進(jìn)行了簡化，忽略慢化劑溫度變化、控制棒移動等對響應(yīng)矩陣的影響，將其定為常數(shù)矩陣。程序中設(shè)定每個探測器僅受同一象限的堆芯影響，該象限對其他象限的探測器沒有貢獻(xiàn)。

本文采用某供熱堆為模型，全堆共有148個組件，結(jié)構(gòu)如圖1、2所示。在擴(kuò)散方程計算時，在堆芯x軸(數(shù)字方向)、y軸(字母方向)上按組件劃分了14個節(jié)塊；在z軸上從下到上(冷卻劑流動方向)平均劃分了16個節(jié)塊。堆外共有4組功率量程探測器，分別位于4個象限內(nèi)，每組由6個敏感段構(gòu)成。

探測器響應(yīng)矩陣為四維矩陣，每個元素表示堆芯某象限空間節(jié)塊對該象限探測器某敏感段的響應(yīng)。現(xiàn)將響應(yīng)矩陣分為兩個二維矩陣表示，第1個矩陣為堆芯軸向各層節(jié)塊對6個探測器敏感段的響應(yīng)矩陣(表1)，第2個矩陣為同一層內(nèi)各節(jié)塊對該層的貢獻(xiàn)矩陣(表2)。

圖1 某供熱堆堆外探測器布置Fig.1 Out-core detector location of nuclear heating reactor

圖2 某供熱堆堆外探測器結(jié)構(gòu)Fig.2 Out-core detector structure of nuclear heating reactor

4 計算測試結(jié)果

用C++語言在Qt平臺上將中子擴(kuò)散程序與特征統(tǒng)計算法耦合。特征統(tǒng)計算法提供特征值及計算特征量的收斂情況，中子擴(kuò)散程序計算目標(biāo)函數(shù)值，兩者耦合可重構(gòu)出堆芯功率分布。在擾動因素產(chǎn)生之前，各組件功率積分歸一化后的功率分布示于圖3。

4.1 溫度場波動

表1 各層節(jié)塊對探測器的響應(yīng)Table 1 Response of each layer of node to detector

表2 各節(jié)塊對該層的貢獻(xiàn)Table 2 Contribution of each node to this layer

圖3 某供熱堆穩(wěn)態(tài)功率分布Fig.3 Power distribution of nuclear heating reactor

在耦合程序中，限定特征量的搜索范圍為[-1.0，1.0]，同時在計算截面變化時，設(shè)定截面必須為正，同時截面相對變化絕對值不超過100%。溫度場波動模型如下：

ΔT(x,y,z)=0.1x+0.2y+0.3z+

0.1xy+0.1xz+0.1yz

(6)

經(jīng)耦合程序計算后，得到最優(yōu)的重構(gòu)結(jié)果列于表3，對應(yīng)的擾動模型多項式系數(shù)列于表4。

表3 溫度場波動重構(gòu)結(jié)果對比Table 3 Comparison of reconstruction results of temperature field fluctuation

表4 溫度場波動模型系數(shù)對比Table 4 Comparison of coefficients of temperature field fluctuation model

表3中，目標(biāo)函數(shù)值F可通過堆外探測器提供數(shù)據(jù)由式(5)計算得到，表征整體重構(gòu)結(jié)果的精確程度；ΔPremax表示各節(jié)塊中最大的相對功率偏差，其可表示搜索方法精細(xì)功率重構(gòu)的可靠程度。表4中，Set表示溫度擾動模型的設(shè)定系數(shù)值，Cal表示重構(gòu)后計算得到的擾動模型系數(shù)值，Cal系數(shù)的精確程度由系數(shù)偏差向量的2-范數(shù)‖ΔX‖2即式(7)表示，該值代表搜索方法找出功率擾動原因的可靠程度。重構(gòu)計算后，F(xiàn)為7.97×10-5，ΔPremax由6.72%降至0.067%，‖ΔX‖2為8.75×10-3，可說明重構(gòu)結(jié)果準(zhǔn)確且可靠。

(7)

4.2 氙振蕩

氙振蕩常發(fā)生在堆芯中間部位，該處的功率波動較大，探測器對此處的氙振蕩不夠敏感。為提高計算精度，在程序中合理減小特征統(tǒng)計算法的搜索范圍，核截面數(shù)據(jù)將根據(jù)實際變化率來設(shè)定核截面的變化范圍。隨機(jī)設(shè)定氙毒波動的模型如下：

ΔXe(x,y,z)=(-0.2x+0.3y+0.15z)·

exp(-0.03x2-0.05y2-0.07z2)

(8)

得到的重構(gòu)結(jié)果列于表5、6。

表5 氙振蕩重構(gòu)結(jié)果對比Table 5 Comparison of reconstruction results of xenon oscillation

表5中，F(xiàn)由4.27×10-3降至7.28×10-6，ΔPremax從1.401%降至0.053%，重構(gòu)結(jié)果精確。表6中重構(gòu)后的系數(shù)向量偏差2-范數(shù)‖ΔX‖2為4.09×10-2，重構(gòu)結(jié)果可靠。

表6 氙振蕩模型系數(shù)對比Table 6 Comparison of coefficients of xenon oscillation model

4.3 氙毒-溫度場波動

隨機(jī)設(shè)定雙影響因素的模型如下：

ΔT(x,y,z)=0.1x-0.2y+

0.3z+0.1xy+0.2xz-0.15yz

(9)

ΔXe(x,y,z)=(-0.2x+0.3y+0.15z)·

exp(-0.03x2-0.05y2-0.07z2)

(10)

確定每個參數(shù)的定義域后，進(jìn)行重構(gòu)計算。重構(gòu)結(jié)果列于表7、8。表7中，F(xiàn)由2.65×10-2降至8.42×10-5，ΔPremax由6.78%降至0.08%，重構(gòu)效果較好。表8中，重構(gòu)后的系數(shù)向量偏差2-范數(shù)‖ΔX‖2為0.118，偏差略大，這是因為溫度擾動與氙振蕩同時發(fā)生時，兩個因素對功率變化同時起作用，程序無法區(qū)分兩個變量對功率的影響，這是所有基于堆外探測數(shù)據(jù)進(jìn)行功率重構(gòu)的方法共同面對的問題。

表7 氙毒-溫度場波動重構(gòu)結(jié)果對比Table 7 Comparison of reconstruction results of xenon and temperature field fluctuation

沿著堆芯中軸線的切面、垂直于堆芯中軸線從下到上第1、5、8、12、16層節(jié)塊的重構(gòu)前、后功率相對偏差數(shù)據(jù)分別示于圖4、5。

表8 氙毒-溫度場波動模型系數(shù)對比Table 8 Comparison of coefficients of xenon and temperature field fluctuation model

圖4 堆芯中心軸向切面功率重構(gòu)前、后相對偏差分布Fig.4 z-axis distributions of power relative deviation at core center before and after reconstruction

圖5 z軸第1、5、8、12、16層節(jié)塊的功率相對偏差分布Fig.5 Distribution of power relative deviations in sections 1, 5, 8, 12 and 16 of z-axis

5 結(jié)論與展望

本文研究了基于堆外測量數(shù)據(jù)進(jìn)行功率重構(gòu)的擾動搜索方法，通過構(gòu)造氙振蕩及溫度場波動單因素算例及雙因素算例，初步驗證了方法的可行性。

單因素算例可重構(gòu)出堆芯功率分布，其最大節(jié)塊功率偏差可控制在0.1%以下，同時可精確搜索出影響因素模型的系數(shù)，該功能可用來提供堆芯功率擾動的具體原因分析；雙因素算例的堆芯功率最大節(jié)塊偏差降到0.1%以下，但擾動模型系數(shù)重構(gòu)沒有單因素算例準(zhǔn)確，是因為雙因素相互干擾，這是該類問題的固有性質(zhì)。

本文為說明問題進(jìn)行了大量簡化。該方法如在實際工程中進(jìn)行應(yīng)用，還要進(jìn)行以下幾方面的研究與開發(fā)：

1) 擾動模型的合理性與可區(qū)分性，如果擾動之間有因果關(guān)系，可進(jìn)行合并分析；

2) 截面數(shù)據(jù)與擾動之間的關(guān)系不僅是常系數(shù)線性關(guān)系，可繼續(xù)研究兩者之間的敏感性；

3) 堆內(nèi)功率分布與堆外探測器之間的響應(yīng)矩陣不是常數(shù)矩陣，該矩陣的影響因素可展開分析；

4) 在擾動模型系數(shù)較多時，堆外測量可能存在多解，特征統(tǒng)計算法對多解的甄別能力需進(jìn)行研究；

5) 目前在役壓水堆的堆外測量系統(tǒng)僅用來實現(xiàn)堆芯軸向功率分布的測量，如將該方法進(jìn)行大型壓水堆堆芯功率重構(gòu)，還需進(jìn)一步研究，結(jié)合堆芯測量系統(tǒng)的堆芯出口溫度測量系統(tǒng)進(jìn)行堆芯功率的三維立體重構(gòu)。