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

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

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

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

穩(wěn)態(tài)中子輸運(yùn)方程[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為外源源強(qiáng)，cm-3·s-1。

對(duì)式(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為離散后的能群號(hào)；i、j、k為離散后的空間網(wǎng)格索引；A、B、C分別為網(wǎng)格單元與R、θ、Z方向正交的網(wǎng)格面的面積，cm2；V為計(jì)算網(wǎng)格體積，cm3；Q為總源強(qiáng)(裂變?cè)础⑸⑸湓春屯庠粗?，cm-3·s-1；α為待定系數(shù)。

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

2 源分布計(jì)算

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

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

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為組件號(hào)；K為每次裂變釋放的能量，MeV；n為核素編號(hào)；f′為裂變份額；ν為每次裂變釋放粒子數(shù)。

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

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

3.1 計(jì)算模型

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

3.2 中子源項(xiàng)計(jì)算

基于HBR-2基準(zhǔn)報(bào)告中提供的第9循環(huán)的pin-by-pin平均功率(圖4a)，利用源轉(zhuǎn)換程序SourceTrans將pin-by-pin功率轉(zhuǎn)換為pin-by-pin裂變中子源，然后映射到RθZ圓柱幾何的計(jì)算網(wǎng)格上。為便于與基準(zhǔn)報(bào)告中的計(jì)算結(jié)果比較，采用基準(zhǔn)報(bào)告中提供的47群中子裂變能譜，即采用235U和239Pu的均勻混合裂變譜[5]，其能群結(jié)構(gòu)與47群的BUGLE-96截面庫一致。圖4b示出了堆芯中平面處第18群的裂變?cè)吹姆植?，可看出其分布與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群裂變?cè)捶植糉ig.4 Distributions of pin-by-pin power and fission neutron source for the 18th group

3.3 計(jì)算條件

基于47群的BUGLE-96截面庫[8]計(jì)算HBR-2裝置的輻照監(jiān)督管處和中子劑量測(cè)量儀處的中子通量密度以及放射性比活度，并分析兩種角度離散(S8和S16)和2種空間離散對(duì)計(jì)算結(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，探測(cè)位置為0.5～0.8 cm；θ方向在探測(cè)位置網(wǎng)格由0.5°逐漸遞減到0.1°，其他位置為0.5°，R-θ剖面如圖5a所示；Z方向測(cè)點(diǎn)位置為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，探測(cè)位置附近為0.5 cm；θ方向探測(cè)位置為0.125°，其他位置為0.5°，R-θ剖面如圖5b所示；Z方向測(cè)點(diǎn)位置為1.0 cm，其他位置為5.0 cm，R-Z剖面如圖6b所示。

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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