張平遜 張 斌 陳義學(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)存限制。
源強(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.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
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。
堆芯源區(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
本文針對屏蔽計(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ā)平臺。