核反應(yīng)堆屏蔽計(jì)算堆芯中子源強(qiáng)生成方法研究

張平遜 張 斌 陳義學(xué)

（華北電力大學(xué) 核科學(xué)與工程學(xué)院 北京 102206）

離散縱標(biāo)法（SN）是粒子輸運(yùn)方程求解中重要的確定論方法之一，在解決深穿透屏蔽問題中具有顯著優(yōu)勢［1］。源強(qiáng)作為重要的輸入?yún)?shù)，計(jì)算精度直接影響屏蔽計(jì)算的可靠性。源強(qiáng)受幾何模型、燃耗與功率分布的影響［2］，計(jì)算結(jié)果不確定性［3-4］較大。全堆芯Pin-by-Pin源強(qiáng)處理可較準(zhǔn)確地描述堆芯源強(qiáng)分布，但幾何建模復(fù)雜繁瑣，實(shí)際應(yīng)用并不理想。針對堆芯中子源強(qiáng)計(jì)算方法，國內(nèi)外開展了大量研究。Haghighat［5］研究發(fā)現(xiàn)，當(dāng)源強(qiáng)計(jì)算考慮堆芯不同裂變核素與采用不同核素處理方法時，將對結(jié)果最高引入10%以上的偏差。西屋公司綜合考慮源強(qiáng)計(jì)算相關(guān)影響因素，為DORT［6］等離散縱標(biāo)程序開發(fā)了SORCERY前端源項(xiàng)處理程序，但受計(jì)算機(jī)內(nèi)存限制，無法處理大規(guī)模三維問題［7］。上海核工程研究設(shè)計(jì)院完善了TORT程序的源項(xiàng)生成功能，提高程序的應(yīng)用范圍［8］?？傮w上，國內(nèi)外源強(qiáng)計(jì)算方法研究已較為成熟，為獲取精確的堆芯源強(qiáng)分布，針對大規(guī)模屏蔽問題采用Pin-by-Pin源強(qiáng)計(jì)算時，需對燃料柵元進(jìn)行大量冗余幾何建模，計(jì)算效率存在較大挑戰(zhàn)。

本文研究基于多維離散縱標(biāo)輸運(yùn)計(jì)算程序ARES［9］進(jìn)行展開，分析源強(qiáng)計(jì)算方法對屏蔽計(jì)算精度的影響。以平均源強(qiáng)計(jì)算方法為基礎(chǔ)，構(gòu)造多權(quán)重源強(qiáng)網(wǎng)格映射算法，避免復(fù)雜Pin-by-Pin幾何建模，提高源強(qiáng)計(jì)算效率與精度，降低程序內(nèi)存限制。

1 理論模型

1.1 中子源強(qiáng)計(jì)算方法

源強(qiáng)計(jì)算作為SN方法屏蔽計(jì)算的前處理過程，輸運(yùn)計(jì)算前需進(jìn)行幾何網(wǎng)格的源強(qiáng)賦值，基本公式為：

式中：χig為固定源能譜；Ci為功率源強(qiáng)轉(zhuǎn)換因子；Pi為燃料組件功率；i為堆芯組件編號；g為能群編號。針對單一裂變核素進(jìn)行源強(qiáng)計(jì)算，功率源強(qiáng)轉(zhuǎn)換因子定義為：

式中：υ為核素每次裂變平均釋放的中子數(shù)；Er為每次裂變平均釋放的能量。由于組件燃耗變化，如裂變核素钚的產(chǎn)生與沉積將對堆芯源強(qiáng)產(chǎn)生重要影響［10-11］?？紤]各裂變核素對源強(qiáng)的貢獻(xiàn)，采用截面平均方法對固定源能譜與功率源強(qiáng)轉(zhuǎn)換因子進(jìn)行加權(quán)平均［5］：

1.2 多權(quán)重源強(qiáng)網(wǎng)格映射算法

1.2.1 中子價(jià)值理論

根據(jù)中子守恒原理，屏蔽問題中不含裂變源的穩(wěn)態(tài)固定源中子輸運(yùn)方程表示為［12］：

式中：從左往右依次為泄漏項(xiàng)、移出項(xiàng)、散射源項(xiàng)和外中子源項(xiàng)為中子角通量密度為宏觀總截面為宏觀散射截面分別為空間、能量和角度變量。

式中：H*為H的共軛算子。根據(jù)價(jià)值守恒理論推導(dǎo)共軛中子輸運(yùn)方程，詳細(xì)過程見參考文獻(xiàn)［1］：

對中子輸運(yùn)方程與共軛輸運(yùn)方程進(jìn)行簡寫：

分別將式（8）和式（9）與ψ*和ψ作內(nèi)積，并將兩式相減，邊界條件為真空邊界，則有：

若函數(shù)Σd(E)為中子與某探測器的響應(yīng)值，則共軛函數(shù)表示探測器對此中子的響應(yīng)值，即中子價(jià)值。屏蔽計(jì)算中，通過一次共軛計(jì)算即可獲得堆內(nèi)各燃料組件對壓力容器處探測器響應(yīng)的貢獻(xiàn)。

1.2.2 多權(quán)重源強(qiáng)網(wǎng)格映射算法

反應(yīng)堆各組件源強(qiáng)對壓力容器處中子注量率的貢獻(xiàn)不同［13］，采用共軛計(jì)算可獲得堆芯共軛通量密度分布，共軛通量密度大的區(qū)域其中子價(jià)值高，產(chǎn)生的中子對壓力容器處中子注量率的影響較大，幾何模型需以燃料柵元為單元進(jìn)行Pin-by-Pin精細(xì)建模。

在壓力容器處構(gòu)建共軛源，經(jīng)共軛輸運(yùn)可得堆芯源強(qiáng)區(qū)域中子價(jià)值分布。如圖1所示，燃料組件由外向內(nèi)中子價(jià)值迅速降低，外圍組件產(chǎn)生的中子對壓力容器處中子注量率起主要貢獻(xiàn)。多權(quán)重源強(qiáng)網(wǎng)格映射算法充分考慮堆芯組件中子價(jià)值分布，對中子價(jià)值大的外圍組件進(jìn)行精細(xì)源強(qiáng)計(jì)算。源強(qiáng)網(wǎng)格與幾何網(wǎng)格在相同笛卡爾直角坐標(biāo)系下處理，與幾何網(wǎng)格劃分不同，源強(qiáng)網(wǎng)格為程序自動處理計(jì)算的區(qū)域坐標(biāo)。幾何網(wǎng)格源強(qiáng)賦值流程如圖2所示，主要分為以下三步：

圖2 幾何網(wǎng)格源強(qiáng)賦值流程Fig.2 Flow chart of geometric grid source assignment

1）由中子價(jià)值對功率分布進(jìn)行近似：外圍兩層組件中子價(jià)值分布為0.1~1.0，對輸運(yùn)計(jì)算產(chǎn)生顯著影響，組件使用Pin-by-Pin計(jì)算，內(nèi)層組件的中子價(jià)值小于0.1，對輸運(yùn)結(jié)果影響較小，可采用均勻近似，經(jīng)用戶定義的幾何輸入卡獲得不同層級的源強(qiáng)區(qū)域邊界與功率分布；

2）根據(jù)各組件燃耗信息，使用截面平均方法對源強(qiáng)裂變譜與功率源強(qiáng)轉(zhuǎn)換因子進(jìn)行加權(quán)平均，結(jié)合源區(qū)邊界與功率獲得區(qū)域源強(qiáng)；

3）遍歷源區(qū)各幾何網(wǎng)格，采用映射算法實(shí)現(xiàn)區(qū)域源強(qiáng)向幾何網(wǎng)格的源強(qiáng)賦值，其中非源區(qū)網(wǎng)格源強(qiáng)直接賦值為零，節(jié)約時間成本。

源強(qiáng)區(qū)域劃分為計(jì)算機(jī)基于功率分布自動處理的虛擬區(qū)域，實(shí)際模型并未進(jìn)行劃分，且源區(qū)邊界與網(wǎng)格邊界位置不同，無法直接進(jìn)行源強(qiáng)賦值。為實(shí)現(xiàn)源強(qiáng)賦值，保證總源強(qiáng)守恒，采用體積權(quán)重法對幾何網(wǎng)格進(jìn)行分割，并對劃分的子網(wǎng)格分別進(jìn)行源強(qiáng)計(jì)算。其中有：

式中：Smesh為幾何網(wǎng)格源強(qiáng)；Vmesh為幾何網(wǎng)格體積；Si為按源強(qiáng)區(qū)域邊界劃分的子網(wǎng)格源強(qiáng)；Vi為子網(wǎng)格體積；M為劃分的子網(wǎng)格數(shù)目。如圖3所示，堆芯某幾何網(wǎng)格跨多個源強(qiáng)區(qū)域，按源區(qū)邊界將網(wǎng)格進(jìn)行分割。分別對分割產(chǎn)生的8個子網(wǎng)格區(qū)域進(jìn)行源強(qiáng)計(jì)算，采用式（13）進(jìn)行幾何網(wǎng)格的源強(qiáng)賦值求和。依次對各燃料組件進(jìn)行網(wǎng)格源強(qiáng)映射計(jì)算，實(shí)現(xiàn)堆芯幾何網(wǎng)格源強(qiáng)賦值。

圖3 幾何網(wǎng)格源強(qiáng)映射方法 (a) 幾何網(wǎng)格分割，(b) 子網(wǎng)格源強(qiáng)計(jì)算，(c) 幾何網(wǎng)格源強(qiáng)賦值Fig.3 Geometric grid source mapping method (a) Geometric grid segmentation, (b) Sub-grid source calculation, (c) Geometric grid source assignment

2 基準(zhǔn)驗(yàn)證

2.1 基準(zhǔn)題描述

NUREG/CR-6115基準(zhǔn)題［14］中標(biāo)準(zhǔn)堆芯燃料裝載方案為典型壓水堆模型，包含204個燃料組件，總功率為2 527.73 MW，該基準(zhǔn)題可用于屏蔽程序的評估與分析?；鶞?zhǔn)題模型堆芯外圍組件徑向功率分布梯度較大，且軸向功率峰因子對源強(qiáng)結(jié)果影響明顯，使用多權(quán)重源強(qiáng)網(wǎng)格映射算法與平均源強(qiáng)算法分別進(jìn)行計(jì)算，分析多權(quán)重源強(qiáng)網(wǎng)格映射算法的計(jì)算精度與效率。圖4為NUREG/CR-6115基準(zhǔn)題的幾何模型，反應(yīng)堆由堆芯向外包含吊籃、下降段、熱屏、水隙、壓力容器、空腔和混凝土，壓力容器內(nèi)半徑為219.075 cm，厚21.59 cm。反應(yīng)堆內(nèi)其他結(jié)構(gòu)尺寸和材料成分在NUREG/CR-6115基準(zhǔn)報(bào)告中詳細(xì)給出。

圖4 NUREG/CR-6115基準(zhǔn)題幾何模型1-圍板，2-吊籃，3-內(nèi)部水隙，4-熱屏，5-壓力容器內(nèi)襯，6-空隙，7-壓力容器隔熱層，8-混凝土內(nèi)襯Fig.4 Geometric diagram of the NUREG/CR-6115 benchmark model 1-Core baffle, 2-Core barrel, 3-Inner inlet water gap,4-Thermal shield, 5-Reactor pressure vessel liner, 6-Air gap, 7-Pressure vessel insulation, 8-Biological shield liner

基于KASHIL-E70多群截面數(shù)據(jù)庫［15］進(jìn)行截面處理，采用離散縱標(biāo)屏蔽程序ARES進(jìn)行三維輸運(yùn)計(jì)算。網(wǎng)格劃分為284×284×161，空間離散格式選擇指數(shù)定向θ權(quán)重差分方法（Exponential Directional Weighted difference method，EDW），各向異性散射截面使用P3階勒讓德多項(xiàng)式進(jìn)行展開近似，求積組采用勒讓德-切比雪夫求積組（PNTN），求積階數(shù)設(shè)置為S8，迭代收斂精度為1×10-3。

2.2 結(jié)果分析

堆芯源區(qū)以燃料組件為基本單元進(jìn)行幾何建模，使用多權(quán)重源強(qiáng)網(wǎng)格映射算法與組件平均源強(qiáng)計(jì)算方法分別進(jìn)行源強(qiáng)計(jì)算，輸運(yùn)參數(shù)設(shè)置相同。圖5為兩種源強(qiáng)計(jì)算方法下三維堆芯源強(qiáng)分布，多權(quán)重源強(qiáng)網(wǎng)格算法軸向源強(qiáng)分布由中間位置向兩端遞減，徑向外圍組件源強(qiáng)呈梯形分布，相比于平均源強(qiáng)計(jì)算方法能更準(zhǔn)確描述堆芯源強(qiáng)。選取壓力容器內(nèi)表面軸向峰值處與焊縫處探測點(diǎn)位置的中子能譜分布進(jìn)行分析，能譜分布結(jié)果如圖6所示。壓力容器軸向峰值處中子能譜分布趨勢總體一致，吻合程度相對較好，但焊縫處中子能譜分布吻合相對較差，平均源強(qiáng)計(jì)算方法輸運(yùn)結(jié)果偏大。偏差產(chǎn)生的主要原因?yàn)檩S向功率峰因子對不同軸向區(qū)域快中子注量率的影響不同，造成輸出結(jié)果穩(wěn)定性較差，圖7為外圍組件軸向功率分布。

圖5 NUREG/CR-6115模型源強(qiáng)分布 (a) 平均源強(qiáng)計(jì)算方法，(b) 多權(quán)重源強(qiáng)網(wǎng)格映射算法Fig.5 Source distribution of the NUREG/CR-6115 model(a) Average source calculation algorithm, (b) Multi-weight source mesh mapping algorithm

圖6 壓力容器內(nèi)表面中子能譜分布 (a) 壓力容器軸向峰值處，(b) 壓力容器焊縫處Fig.6 Neutron spectrum distribution on internal surface of pressure vessel (a) Axial peak location, (b) RPV lower weld location

圖7 NUREG/CR-6115模型軸向歸一化相對功率分布Fig.7 Axial-normalized relative power distribution of the NUREG/CR-6115 model

為比較兩種源強(qiáng)計(jì)算方法的優(yōu)劣，對總體輸運(yùn)計(jì)算結(jié)果進(jìn)行分析。選取壓力容器不同軸向位置隨方位角變化（0°~45°）的探測點(diǎn)統(tǒng)計(jì)快中子注量率（E＞1.0 MeV）輸運(yùn)結(jié)果，按式（14）計(jì)算與基準(zhǔn)題報(bào)告中MCNP和DORT參考值的相對誤差均方根：

式中：φcalc，i為輸運(yùn)計(jì)算值；φref，i為基準(zhǔn)參考值；N為選取的探測點(diǎn)數(shù)目。

表1為NUREG/CR-6115模型壓力容器處不同軸向位置快中子注量率ERMS結(jié)果，總體上多權(quán)重源強(qiáng)網(wǎng)格映射算法與平均源強(qiáng)計(jì)算方法相比ERMS偏小。其中多權(quán)重源強(qiáng)網(wǎng)格映射算法輸運(yùn)結(jié)果與DORT參考值相比ERMS最大約為14.4%，與MCNP參考值相比最大約為5.7%；平均源強(qiáng)計(jì)算方法輸運(yùn)結(jié)果與DORT參考值相比ERMS最大約為29.7%，與MCNP參考值相比最大約為21.2%。由于平均源強(qiáng)計(jì)算方法將各組件源強(qiáng)進(jìn)行平均，造成壓力容器不同軸向高度位置快中子注量率ERMS不同，焊縫處ERMS比軸向峰值處增大近一倍。相比于平均源強(qiáng)計(jì)算方法，多權(quán)重源強(qiáng)網(wǎng)格映射算法輸運(yùn)結(jié)果更精確，其中壓力容器焊縫處計(jì)算結(jié)果與MCNP參考值相比ERMS降低約18.46%。

表1 反應(yīng)堆壓力容器快中子注量率（E＞1.0 MeV）相對誤差均方根Table 1 Root mean square (RMS) of the relative error for fast neutron fluence in reactor pressure vessel (E＞1.0 MeV)

多權(quán)重源強(qiáng)映射方法建模方便，使用靈活的網(wǎng)格映射技術(shù)，可實(shí)現(xiàn)任意位置源強(qiáng)賦值。以外圍三層組件Pin-by-Pin輸運(yùn)結(jié)果為參考值，選擇壓力容器處探測點(diǎn)位置，分別對比不同分層下快中子注量率分布結(jié)果。表2為源強(qiáng)映射方法下不同層級數(shù)源強(qiáng)相對誤差分布，能量E＞1.0 MeV時，平均源強(qiáng)計(jì)算方法與外圍三層組件Pin-by-Pin快中子注量率相對偏差的平均值為20.68%，能量E＞0.1 MeV時相對偏差為19.7%。外圍一層組件與外圍兩層組件平均相對功率偏差小于0.5%，其中外圍兩層組件輸運(yùn)計(jì)算結(jié)果與參考值基本吻合，滿足計(jì)算精度要求。

表2 中子源強(qiáng)分布影響Table 2 Influence of neutron source distribution

3 結(jié)語

本文針對屏蔽計(jì)算中復(fù)雜源強(qiáng)處理問題，提出了多權(quán)重源強(qiáng)網(wǎng)格映射算法。通過對幾何權(quán)重較大的組件進(jìn)行精細(xì)計(jì)算，并采用體積權(quán)重法將區(qū)域源強(qiáng)向幾何網(wǎng)格進(jìn)行映射，保證總源強(qiáng)守恒。避免復(fù)雜幾何建模，提高計(jì)算效率，并有效降低計(jì)算機(jī)內(nèi)存限制。數(shù)值計(jì)算表明，多權(quán)重源強(qiáng)網(wǎng)格映射算法相比于平均源強(qiáng)計(jì)算方法，提高了輸運(yùn)計(jì)算的精度，且保障軸向位置計(jì)算結(jié)果的穩(wěn)定性。另外，光子輸運(yùn)在屏蔽計(jì)算中占有重要作用，光子源強(qiáng)將是未來研究的工作重點(diǎn)。

作者貢獻(xiàn)聲明張平遜負(fù)責(zé)文獻(xiàn)調(diào)研、程序功能開發(fā)與測試，以及文章撰寫和修訂；張斌負(fù)責(zé)相關(guān)文獻(xiàn)資料的搜集與整理并對文章進(jìn)行修訂；陳義學(xué)提供程序研發(fā)平臺。