楊 超，程湯培，鄧 力，付元光，溫麗麗

(1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100094；2.中物院高性能數(shù)值模擬軟件中心，北京 100088)

基于并行自適應(yīng)結(jié)構(gòu)網(wǎng)格應(yīng)用支撐軟件框架JASMIN[1]，中物院高性能數(shù)值模擬軟件中心研制了三維并行離散縱標(biāo)(SN)屏蔽計算程序JSNT[2]，主要應(yīng)用于核電廠的屏蔽計算。JSNT程序基于區(qū)域分解實現(xiàn)了大規(guī)模并行計算和多種迭代加速算法[3]，具有可視化建模和后處理的功能，擁有對諸如反應(yīng)堆、廠房等大型裝置進行精確建模和高分辨率計算的能力[4]。H. B. Robinson 2號機組(HBR-2)是熱功率為2 300 MW的壓水堆裝置，對外公開發(fā)布了輻照監(jiān)督管處和中子劑量測量儀處的比活度測量值，是屏蔽計算程序校驗的經(jīng)典算例。本文利用JSNT對HBR-2裝置進行屏蔽計算與分析。首先利用可視化建模程序JLAMT建立HBR-2的精細(xì)模型，然后基于HBR-2基準(zhǔn)報告[5]中pin-by-pin的三維功率分布，利用源項轉(zhuǎn)換程序SourceTrans計算得到屏蔽計算所需的三維裂變中子源項，最后利用JSNT程序計算輻照監(jiān)督管處和中子劑量測量儀處的中子通量密度分布以及基準(zhǔn)報告中給出的6個核素的放射性比活度，并與實驗測量值進行對比。

1 離散縱標(biāo)方法

穩(wěn)態(tài)中子輸運方程[6]：

E′)ψ(r,E′,Ω′)dE′dΩ′+S(r,E,Ω)

(1)

式中：ψ為中子角通量密度，cm-2·s-1；r為位置向量；E為能量變量；Ω為方向向量；Σt為宏觀總截面，cm-1；Σs為從能量E′、角度Ω′散射到E和Ω的宏觀散射截面，cm-1；Σf為宏觀裂變截面，cm-1；χ為中子裂變譜；S為外源源強，cm-3·s-1。

對式(1)角度方向上采用SN方法離散，空間上采用差分，能量上采用分群近似，則有：

ωmμm(Ai+1/2,j,kψi+1/2,j,k,m,g-Ai-1/2,j,kψi-1/2,j,k,m,g)+

ωmξm(Bi,j+1/2,kψi,j+1/2,k,m,g-Bi,j-1/2,kψi,j-1/2,k,m,g)+

ωmηm(Ci,j,k+1/2ψi,j,k+1/2,m,g-Ci,j,k-1/2ψi,j,k-1/2,m,g)+

(Ai+1/2,j,k-Ai-1/2,j,k)(αm+1/2ψi,j,k,m+1/2,g-

αm-1/2ψi,j,k,m-1/2,g)+ωmΣt,i,j,k,gVi,j,kψi,j,k,m,g=

ωmVi,j,kQi,j,k,m,g

(2)

式中：ωm為方向m上的權(quán)重；g為離散后的能群號；i、j、k為離散后的空間網(wǎng)格索引；A、B、C分別為網(wǎng)格單元與R、θ、Z方向正交的網(wǎng)格面的面積，cm2；V為計算網(wǎng)格體積，cm3；Q為總源強(裂變源、散射源和外源之和)，cm-3·s-1；α為待定系數(shù)。

JSNT采用源迭代方法對式(2)進行求解，外迭代用于更新裂變源及特征值，內(nèi)迭代用于更新散射源。基于JASMIN框架、計算單元以及邊界上反射方向間數(shù)據(jù)依賴關(guān)系創(chuàng)建有向非循環(huán)圖(DAG)，根據(jù)DAG實現(xiàn)區(qū)域分解的大規(guī)模并行計算。

2 源分布計算

對于裂變源項的計算，首先根據(jù)堆芯的三維棒功率分布、每次裂變產(chǎn)生中子數(shù)和能量、裂變份額以及裂變譜等參數(shù)計算棒上的源強，然后將棒與計算網(wǎng)格進行映射。對于直角坐標(biāo)，采用體積權(quán)重方法進行映射，如圖1a所示；對于柱坐標(biāo)系，將棒離散為帶權(quán)的點，它們在棒內(nèi)均勻分布，然后根據(jù)柱計算網(wǎng)格所覆蓋的離散點的數(shù)目，計算得到該計算網(wǎng)格內(nèi)的源強，如圖1b所示。

三維空間棒網(wǎng)格源分布計算公式如下：

S(i,j,k,g)=P(i,j,k)·C′·B′(m)/K′(m)·

(3)

式中：P為棒網(wǎng)格功率密度，W/cm3；C′為單位轉(zhuǎn)換因子，即將MeV轉(zhuǎn)換為W，6.24×1012MeV·s-1·W-1；B′為組件偏倚因子；m為組件號；K為每次裂變釋放的能量，MeV；n為核素編號；f′為裂變份額；ν為每次裂變釋放粒子數(shù)。

a——直角坐標(biāo)網(wǎng)格；b—— 柱坐標(biāo)網(wǎng)格圖1 不同坐標(biāo)系下的映射Fig.1 Remapping process of different mesh ordinates

3 計算模型與計算結(jié)果分析

3.1 計算模型

HBR-2快中子注量率實驗基準(zhǔn)題由ORNL在1997年發(fā)布，其目的是驗證輻射屏蔽輸運計算程序和數(shù)據(jù)庫的正確性。HBR-2基準(zhǔn)題報告[5]給出了詳細(xì)的反應(yīng)堆幾何、材料、堆芯功率分布和反應(yīng)堆運行歷史等數(shù)據(jù)，并提供了第9循環(huán)末輻照監(jiān)督管和輻照劑量測量儀處比活度的測量數(shù)據(jù)。HBR-2熱功率為2 300 MW，堆芯裝有157盒燃料組件，每盒燃料組件內(nèi)有204個燃料元件，燃料元件布置方式為15×15。利用JLAMT軟件對HBR-2裝置進行建模，JLAMT是北京應(yīng)用物理與計算數(shù)學(xué)研究所開發(fā)的一款面向粒子輸運程序的自動、可視建模軟件[7]，可將三維模型自動輸出轉(zhuǎn)換為GDML文件與JSNT對接。HBR-2整體幾何結(jié)構(gòu)如圖2所示，堆芯中平面處的剖面如圖3所示，堆芯外圍結(jié)構(gòu)包括圍板、吊籃、熱屏、反應(yīng)堆壓力容器和生物屏蔽等部件。

3.2 中子源項計算

基于HBR-2基準(zhǔn)報告中提供的第9循環(huán)的pin-by-pin平均功率(圖4a)，利用源轉(zhuǎn)換程序SourceTrans將pin-by-pin功率轉(zhuǎn)換為pin-by-pin裂變中子源，然后映射到RθZ圓柱幾何的計算網(wǎng)格上。為便于與基準(zhǔn)報告中的計算結(jié)果比較，采用基準(zhǔn)報告中提供的47群中子裂變能譜，即采用235U和239Pu的均勻混合裂變譜[5]，其能群結(jié)構(gòu)與47群的BUGLE-96截面庫一致。圖4b示出了堆芯中平面處第18群的裂變源的分布，可看出其分布與pin-by-pin平均功率分布是一致的。

圖2 HBR-2幾何結(jié)構(gòu)圖Fig.2 Configuration of HBR-2

圖3 HBR-2剖面圖Fig.3 Schematic section of HBR-2

a——功率分布；b——源分布圖4 pin-by-pin功率分布和第18群裂變源分布Fig.4 Distributions of pin-by-pin power and fission neutron source for the 18th group

3.3 計算條件

基于47群的BUGLE-96截面庫[8]計算HBR-2裝置的輻照監(jiān)督管處和中子劑量測量儀處的中子通量密度以及放射性比活度，并分析兩種角度離散(S8和S16)和2種空間離散對計算結(jié)果的影響。2種空間離散的方式如下。

1) CASE1。201×70×96離散[9]，R方向堆芯網(wǎng)格由內(nèi)到外依次為5.0、2.0、1.0、0.8 cm，結(jié)構(gòu)材料為0.8～1.0 cm，探測位置為0.5～0.8 cm；θ方向在探測位置網(wǎng)格由0.5°逐漸遞減到0.1°，其他位置為0.5°，R-θ剖面如圖5a所示；Z方向測點位置為1.0 cm，其他位置為5.0 cm，R-Z剖面如圖6a所示。

圖5 中平面處R-θ方向網(wǎng)格Fig.5 Mesh of R-θ direction at core mid-plane

圖6 R-Z方向網(wǎng)格Fig.6 Mesh of R-Z direction

2) CASE2。274×195×104離散，R方向堆芯網(wǎng)格由內(nèi)到外依次為2.0、0.8 cm，結(jié)構(gòu)材料為0.8～1.0 cm，探測位置附近為0.5 cm；θ方向探測位置為0.125°，其他位置為0.5°，R-θ剖面如圖5b所示；Z方向測點位置為1.0 cm，其他位置為5.0 cm，R-Z剖面如圖6b所示。

3.4 計算結(jié)果與分析

1) JSNT的計算結(jié)果與分析

針對201×70×96計算網(wǎng)格，截面采用P3和方向離散采用S8的方案，JSNT在某浪潮集群上利用1個節(jié)點完成并行計算，該集群1個計算節(jié)點上共有24個處理器，其計算時間為0.8 h?？砂l(fā)現(xiàn)，JSNT通過并行較大地提高了計算效率，降低了計算時間。圖7示出了中子通量密度的三維空間分布，可看出中子通量密度沿徑向逐漸衰減。

圖7 中子通量密度分布Fig.7 Distribution of neutron flux density

2) JSNT與DORT計算結(jié)果的比較

圖8、9分別示出了兩計算點處的47群中子通量密度JSNT與DORT的計算結(jié)果的比較，其中DORT是基于BUGLE-96截面計算的結(jié)果[5]。在輻照監(jiān)督管處，最大相對偏差小于20%，大于1 MeV的能群相對偏差均小于5%；在中子劑量測量儀處，大于1 MeV的能群相對偏差小于10%，在低能區(qū)存在較大偏差。

圖8 輻照監(jiān)督管處的中子通量密度分布Fig.8 Distribution of neutron flux density at location of surveillance capsule

圖9 中子劑量測量儀處的中子通量密度分布Fig.9 Distribution of neutron flux density at location of neutron dosimeter

由于計算放射性比活度的反應(yīng)均是閾能反應(yīng)，低能區(qū)偏差大對結(jié)果影響較小。

表1、2列出了輻照監(jiān)督管處和中子劑量測量儀處的比活度的計算值與測量值的比值。隨著角度離散數(shù)和計算網(wǎng)格數(shù)的增加，計算結(jié)果趨于測量值。在輻照監(jiān)督管處，JSNT與DORT計算的比活度較一致，與測量的相對偏差均小于20%。在中子劑量測量儀處，除237Np和238U的裂變反應(yīng)外，計算的比活度與測量值之間的相對偏差也小于20%。與DORT計算結(jié)果相比，JSNT計算的237Np和238U的裂變反應(yīng)的比活度偏低，這是由于DORT的計算中子通量在高能區(qū)相對較大引起的。在中子劑量測量儀處，237Np和238U的裂變反應(yīng)比活度與測量值偏差較大，其原因可能是與堆芯距離較遠，中子通量密度相對較小，使實驗測量難度增加，且裂變反應(yīng)在工程上是通過測量它的裂變產(chǎn)物來得到反應(yīng)率的，其測量的不確定度本身比較大；另外屏蔽截面數(shù)據(jù)庫和響應(yīng)截面的精度也是導(dǎo)致偏差較大的原因之一。

表1 輻照監(jiān)督管處比活度的計算值與測量值的比值Table 1 Ratios of calculated-to-measured specific activities at location of surveillance capsule

表2 中子劑量測量儀處比活度的計算值與測量值的比值Table 2 Ratios of calculated-to-measured specific activities at location of neutron dosimeter

4 結(jié)論

基于HBR-2基準(zhǔn)，利用自主研發(fā)的并行程序JSNT完成了該裝置的屏蔽計算。計算結(jié)果表明，當(dāng)采用較多的網(wǎng)格數(shù)時，在輻照監(jiān)督管處各核素的放射性比活度的計算值和測量值的相對偏差均小于20%，滿足屏蔽計算的工程需求。但在中子劑量測量儀處，237Np和238U的裂變反應(yīng)比活度的計算值和測量值存在較大偏差，有待進一步研究。