使用并行3-D辐射输运代码的反应堆剂量测定应用

具体实施方式

本发明涉及系统的辐射场分布的计算方法。辐射场可以包括中子和伽马辐射。在本发明的一方面中，可以计算核反应堆腔中的容器外中子剂量测定响应。核反应堆的种类不受限制，且可以包括本领域中已知的各种商业设计。合适的反应堆可以包括但不限于：压水反应堆(PWR)和沸水反应堆(BWR)。为了说明简便，将参考2-环路商业PWR来说明本发明的该方面。在PWR的设计和操作中使用产生剂量测定响应的PWR的计算机建模。

本发明的方法包括3-D并行辐射输运代码的应用，3-D并行辐射输运代码在这里称为RAPTOR-M3G(辐射的快速并行输运一多3D几何构造)。输运代码提供用于求解S_N方程的并行算法的集合。该方法基于域分解算法，其中，空间域、角度域和/或能量域被划分为能够在多处理器体系结构上被独立地分配及处理的子集。已知的合适的3-D并行确定性输运代码的例子包括PENTRAN^TM和PARTISN(Sjoden G.E.和Haghighat A.，“PENTRAN-Parallel EnvironmentNeutralparticle TRANport in 3-D Cartesian Geometry，”Proceedingsof the Joint International Conference on Mathematical Methods andSupercomputing for Nuclear Applications，Vol.1，pp.232-234，Saratoga Springs，NY(1997))。与传统的单处理器应用相比，本发明的方法减少了每个处理器的计算负载以及存储器需求，为大的3-D问题提供了高效的求解方法。

使用消息传递接口(MPI)并行库在Fortran 90中开发了RAPTOR-M3G计算机代码。(Gropp W.，Lusk E.，and Skjellum A.，Using MPI Portable Parallel Programming with the Message PassingInterface，The MIT Press，Cambridge，Massachusetts(1999))。RAPTOR-M3G的若干特征包括下面的内容：

·非均匀正交构造的网格上的3-D笛卡尔(RAPTOR-XYZ)和圆柱几何构造(RAPTOR-RTZ)上的多群(multi-group)S_N方程的解(M.A.Hunter，G.Longoni，and S.L.Anderson，“Extension of RAPTOR-M3Gto r-θ-z geometry for use in reactor dosimetry applicatons，”Proceedingsof the 13^th Internatioanl Symposium on Reactor Dosimetry，TheNetherlands(2008))；

·空间域、角度域以及耦合的空间/角度域分解算法；

·正定加权差分方案(positive definite weighted differencingscheme)：0/θ加权的以及定向θ加权的；

·直到级数20的能级对称求积组的自动生成(Longoni G.et al.，“Investigation of New Quadrature Sets for Discrete Ordinates Methodwith Application to Non-Conventional Problems，”Transactions of theAmerican Nuclear Society，Vol.84，pp.224-226(2001))；

·并行存储器：允许空间子域和角度子域的本地(local)分配，从而减少每个处理器的存储器需求；

·并行分配任务：在多处理器上对S_N方程并行求解，与单处理器技术相比减少了计算时间；

·并行I/O：每个处理器本地访问其存储设备，以减少I/O时间；以及

·与BOT3P(R.Orsi，“Potential Enhanced Performances inRadiation Transport Analysis on Structured Mesh Grids Made Availableby BOT3P，”Nuclear Science and Engineering，Vol.157，pp.110-116(2007))、自动网格生成器、以及GIP多群截面预处理器的兼容性和集成性。

图1a示出本发明实施例的3-D输运模型的几何构造，例如2-环路PWR。PWR可以包括：12英尺核堆芯(nuclear core)、热屏蔽设计以及3英寸反应堆腔空隙。模型几何构造包括：堆芯-水混合物(core-water mixture)、堆芯罩(core shroud)、堆芯筒(core barrel)、热屏蔽、包括不锈钢衬里(liner)的反应堆压力容器(RPV)以及反射绝缘体。PWR的RPV典型地具有一般圆柱形状，并且例如通过底封头(bottom head)和可移除的顶封头(top head)在两端处封闭。使用钢-水混合物对反应堆堆芯上面和下面的上层内部(upperinternal)区域和下层内部区域进行建模。RPV的下层内部包括堆芯筒(即，堆芯支撑结构)。该堆芯筒被堆芯筒和RPV的内壁之间的热屏蔽围住。在某些情况下，使用中子垫代替热屏蔽。堆芯罩设置在堆芯筒内部。环形下水管围绕反应堆堆芯筒。使冷却液(一般是水)循环到下水管中。

图1b示出2-环路PWR的在z＝0.0cm处的x-y平面上的该模型的2-D截面。此外还示出PWR中的材料分布。使用BOT3P代码的5.2版本来产生模型几何构造和网格离散化。该模型沿着x轴和y轴从0.0cm延伸到245.0cm，并沿着z轴从-200.0cm延伸到200.0cm。在整个模型中应用均匀的网格；沿着x轴、y轴和z轴分别指定2.0×2.0×4.0cm的网格大小，得到总共1,464,100个网格。

使用BUGLE-96截面库(RSICC Data Library CollectionBUGLE-96，“Coupled 47Neutron，20Gamma-Ray Group CrossSection Library Derived from ENDF/B-VI for LWR Shielding andPressure Vessel Dosimetry Applications，”Oak Ridge NationalLaboratory，Oak Ridge，TN(1999))和GIP计算机代码(即DOORS包的一部分)(RSICC Computer Code Collection DOORS 3.2a，“One-，Two-and Three-Dimensional Discrete Ordinates Neutron/PhotonTransport Code System，”Oak Ridge National Laboratory，OakRidge，TN(2003))来处理输运模型中材料混合物的截面。S₈能级对称求积组和散射核的P₃球谐展开被用于输运计算。无源中子检测器的系统可以被安装在反射绝缘体和压力容器之间的反应堆腔空隙中。剂量测定系统可以提供关于反应堆容器的束带区之上的快中子照射的准确信息。纯金属箔可以被安装在反应堆腔中，装入铝壳内，其使得快中子谱的失真最小，有效地产生自由场测量。输运模型中没有明确地定义安装在反应堆腔空隙中的中子剂量测定器。

本发明的一方面包括用于S_N方程的离散化的域分解算法，该算法使用这样的方法：其中，空间能量域和/或角度能量域被划分为能够在多处理器体系结构上被独立地分配及处理的各子集。

S_N方程的空间和角度离散化、以及这里描述的角度域分解算法是代码的3D笛卡尔XYZ版专用。由于散射重新分布期间的存在，为RAPTOR-RTZ开发的S_N方程的方程式与RAPTOR-XYZ不同。

S_N方程的相空间被离散化，即，角度、空间和能量。因此，所得到的线性代数方程的集合适于在数字计算机上求解。使用多群方法将能量域离散化为很多离散间隔，即，g＝1...G，以最高能量粒子开始(g＝1)，以最低能量粒子结束(g＝G)。多群近似中的输运方程(即，LBE)以方程(1)表示。

$\hat{\Omega }·\overrightarrow{▿}{\psi }_g(\overrightarrow{r},\hat{\Omega })+{\sigma }_g\left(\overrightarrow{r}\right){\psi }_g(\overrightarrow{r},\hat{\Omega })=\underset{g^{\prime }=1}{\overset{G}{\Sigma }}\underset{4x}{\int }d{\Omega }^{\prime }{\sigma }_{{\mathrm{gg}}^{\prime }}(\overrightarrow{r},{\hat{\Omega }}^{\prime }·\hat{\Omega }){\psi }_{g^{\prime }}(\overrightarrow{r},{\hat{\Omega }}^{\prime })$

$+\frac{1}{k}{\chi }_g\underset{g^{\prime }=1}{\overset{G}{\Sigma }}v{\sigma }_{f,g^{\prime }}\left(\overrightarrow{r}\right){\phi }_{g^{\prime }}\left(\overrightarrow{r}\right)+q_g^e(\overrightarrow{r},\hat{\Omega }).---\left(1\right)$

通过考虑方向的有限集并通过应用适当的求积积分方案(quadrature integration scheme)对角度域进行离散化。每个离散方向可以被形象化为具有相关表面面积的单位球的表面上的点，所述相关表面面积在数学上与求积方案的权重相对应。离散方向和对应权重的组合被称为求积组。通常，求积组应该满足很多条件，以便准确并在数学上被确定；可以采用若干方法，例如，能级对称求积组(LQn)和基于Legendre多项式的求积组(Longoni G.andHaghighat A.，“Development of New Quadrature Sets with theOrdinate Splitting Technique，”Proceedings of the ANS InternationalMeeting on Mathematical Methods for Nuclear Applications(M&C2001)，Salt Lake City，UT，September 9-13，2001，American NuclearSociety，Inc.，La Grange Park，IL(2001))。在RAPTOR-M3G中开发和实现的求积组基于LQn方法。

可以用若干技术对空间变量进行离散化，例如，有限差分方法和有限元方法。在RAPTOR-M3G中开发的方程式基于有限差分方法，该有限差分方法包括将空间域划分为例如细网格的计算单元，其中，假设在每个单元内的截面是恒定的。在3D笛卡尔几何构造中，使用方程(2)估计在单元中心位置处的角通量。

${\psi }_{i,j,k,m,g}=\frac{q_{i,j,k}+\frac{\left|{\mu }_m\right|}{a_{i,j,k,m,g}\mathrm{\Delta x}}{\psi }_{\mathrm{in},x}+\frac{\left|{\eta }_m\right|}{b_{i,j,k,m,g}\mathrm{\Delta y}}{\psi }_{\mathrm{in},y}+\frac{\left|{\xi }_m\right|}{c_{i,j,k,m,g}\mathrm{\Delta z}}{\psi }_{\mathrm{in},z}}{\frac{\left|{\mu }_m\right|}{a_{i,j,k,m,g}\mathrm{\Delta x}}+\frac{\left|{\eta }_m\right|}{b_{i,j,k,m,g}\mathrm{\Delta y}}+\frac{\left|{\xi }_m\right|}{c_{i,j,k,m,g}\mathrm{\Delta z}}+{\sigma }_{i,j,k}}---\left(2\right)$

在方程(2)中，角度和能量相关性分别由下标m和g来表示。项q_i，j，k表示位于单元中心的散射、裂变以及外部源的和。下标i、j、k表示单元中心值，且权重a_{i，j，k，m，g}、b_{i，j，k，m，g}和c_{i，j，k，m，g}被限制到0.5和1.0之间的范围；RAPTOR-M3G使用θ加权(TW)、0加权(ZW)或的自适应定向θ加权(DTW)差分方案来计算输运扫描(transport sweep)期间的各权重(B.Petrovic and A.Haghighat，“New DirectionalTheta-Weighted S_N Differencing Scheme and Its Application to PressureVessel Fluence Calculations”Proceedings of the 1996RadiationProtection and Shielding Topical Meeting，Falmouth，MA，Vol.1，pp.3-10(1996))。

通过经历从问题域的边界开始的每个方向来解S_N方程；该求解过程也被称为输运扫描。从边界条件或从先前在相邻单元中计算的边界角通量开始估计在中心单元位置定义的角通量。使用方程(2)计算单元中心角通量。使用称为“差分方案”的附加关系来计算离开(exit)计算单元的角通量。

在称为源叠代的叠代过程内执行输运扫描，该源叠代也被公知为定点叠代，或理查森叠代(Richardson iteration)。继续该过程，直到满足适当的收敛性判别标准，即，在两个叠代之间的以任何标准的标量通量上的相对误差在某个截止值之下(Adams M.L.and Larsen E.W.，“Fast Iterative Methods for Discrete-Ordinates Particle TransportCalculations，”progress in Nuclear Energy，Vol.40，n.1(2002))。对于辐射屏蔽计算，该截止值一般被设定为1.0e^-3或1.0e^-4。

在RAPTOR-M3G中开发的并行算法基于处理器网络上角度域和/或空间域的分解。RAPTOR-M3G创建虚拟拓扑，相应地，多个处理器被分配到角度域和空间域，被分配到角度域和空间域的处理器的数目分别被指定为P_a和P_s。任何分解所需的处理器的总数是P_n＝P_a·P_s。基于该信息，处理器网络被映射到空间域和角度域上，创建虚拟拓扑，该虚拟拓扑将每个处理器关联到该处理器的本地子域。

角度域被按八分区(octant)划分，其中，在角度域上指定的处理器被按顺序分配到本地的八分区。每个处理器所分配的八分区的本地数目由方程(3)给定。

$N_{\mathrm{loct}}=\frac{8}{P_a}---\left(3\right)$

在P_a个处理器上的N_loct个八分区上本地执行输运扫描；用于角度域的MPI通信器被用于对各处理器之间的角通量进行同步，并导致反射边界条件。

通过按顺序将P_s个处理器分配给若干x-y平面，沿z轴划分空间域。在P_s个处理器上划分沿z轴的细网格的总数(即，km)；映射阵列(即kmloc)用于将x-y平面分配给P_s个处理器。分配给P_s个处理器的x-y平面的数目是任意的；然而，需要满足方程(4)中的条件，以便限定与问题几何构造在拓扑上一致的空间分解。

$\underset{i=1}{\overset{P_s}{\Sigma }}\mathrm{kmloc}\left(i\right)=\mathrm{km}---\left(4\right)$

将空间域上的处理器映射到任意数目的x-y平面的灵活性可以取决于这样的事实：z平面的数目可能不能被空间域上的处理器数目整除。P_s个处理器上的x-y平面的不均匀划分可能导致处理器负载不均衡，从而导致性能的损失。在本发明中，可以应用混合角度/空间分解策略以克服该困难。混合分解包括角度域和空间域的结合，以包括这些域的并发划分。在下面的例1中进一步说明混合域分解。

图5示出本发明实施例的流程图，其中，实现3-D辐射输运计算机程序，以产生辐射场分布。本实施例包括获得关于要建模的系统的几何和材料信息。该信息可以从各种源获得，诸如但不限于核反应堆图。然后可以选择计算所需的充足的S_N级数。此外，可以选择角度域(即，散射核和角通量)的P_n展开级数。还可以选择适当的差分方案(即，TW、ZW或DTW)。可以使用笛卡尔XYZ或RTZ几何构造的适当网格产生器来产生和离散化系统的3-D模型。通过将每种材料与适当的数据集(诸如BUGLE-96)混合，为每种材料产生截面表。针对要解决的问题选择处理器的数目和相应的分解策略。然后执行计算。在计算后可以进行对所产生的结果的后处理以及分析。

用于本发明的RAPTOR-M3G计算机代码提供辐射输运问题的准确和高效(例如，减少的计算时间)的解决方案。在本发明的一方面中，针对核反应堆容器中的腔计算容器外中子剂量测定响应。与实际测量相比较，通过RAPTOR-M3G计算的2-环路PWR的反应堆腔空隙中的快中子反应的平均准确度为96％。此外，在20个处理器计算机集群上使用混合角度/空间域分解策略在大约106分钟的时钟时间内获得输运问题的解。

虽然已经详细说明了本发明的具体实施例，本领域技术人员应该了解根据本公开的全部教导可以开发对那些细节的各种修改和替代。例如，这里说明的本发明的一个方面涉及核工业，特别是核反应堆。然而，本发明还可以用于诸如医学领域的广泛范围的其它应用。例如，本发明可以被用于确定为治疗和/或治愈癌症对患者施送(deliver)的辐射剂量。因此，公开的特定方案对于本发明的范围只是说明性的，而不是限制性的，本发明的范围由所附权利要求书及其所有等同物的全部范围给出。

例子

例1——RAPTOR-M3G并行性能分析

在“例子”中讨论的输运计算用运行于20个处理器计算机集群(即，EAGLE-1)上的RAPTOR-M3G执行。该集群由具有2双核双处理器AMD Opteron 64位体系结构的5个节点组成。可用的集群总存储器(即，RAM)是40G字节；网络互联的特征在于1GBit/s(比特/s)带宽。利用该硬件配置，RAPTOR-M3G在20个处理器上在大约106分钟内完成对2-环路PWR的全部3D输运计算。使用DTW、TW或ZW差分方案观察不到显著的性能差异。

此外，设立简单的测试问题以分析代码的并行性能。该测试问题由具有均匀分布的固定源的、以1cm的均匀网格离散化的50×50×50cm的箱(box)组成。与一能量群截面集一起使用S₈求积组和P0各向同性散射。使用壁钟(wall-clock)时间、加速和并行效率来评估RAPTOR-M3G的并行性能。在方程5和6中分别定义了加速和并行效率。

S_p＝T_s/T_p    (5)

η_p＝S_p/N_p    (6)

其中，T_s和T_p分别是单处理器计算和多处理器计算所需的壁钟时间。N_p是用于实现壁钟时间T_p的处理器的数目。图4示出对于多达20个处理器使用不同的分解策略获得的加速的比较。

利用空间分解获得的加速随着处理器数目的增加而逐渐减小。据信，该行为是由于每个处理器的更细的计算粒度引起的；随着空间域被分解为更小的子域，每个处理器的操作数目减少，而处理器之间的通信时间增加；因此导致性能下降。各节点之间的网络数据传输通常是分布式存储器体系结构的限制因素。使问题收敛所需的更大数量的叠代进一步使得空间分解策略的性能降低。然而，角度域和空间域被并发划分的混合分解产生了更好的结果。据信，该行为是由于由该分解引入的更粗糙的计算粒度引起的；而且对于混合分解，使问题收敛所需的叠代数目不像空间分解增加的那么多。

例子2——测量的剂量测定响应与RAPTOR-M3G计算的响应的比较

在测量的剂量测定响应和用RAPTOR-M3G获得的相应预测之间进行比较。使用IRDF-2002剂量测定库(I.Kodeli and A.Trkov，“Validation of the IRDF-2002Dosimetry Library”，NuclearInstruments and Methods in Physics Research Section A：Accelerators，Spectrometers，Detectors and Associated Equipment，Vol.57，Issue 3，pp.664-681(2007))，以产生针对表1中所列中子反应的计算出的剂量测定响应。

将表1中所列反应的测量的剂量测定响应与用RAPTOR-M3G计算的响应相比较。

表1.用剂量测定系统测量的中子反应

材料      反应

铜        ⁶³Cu(n，α)⁶⁰Co

铁        ⁵⁴Fe(n，p)⁵⁴Mn

镍        ⁵⁸Ni(n，p)⁵⁸Co

铀        ²³⁸U(n，f)¹³⁷Cs

镎        ²³⁷Np(n，f)¹³⁷Cs

使用镉屏蔽的金属箔来测量表1中列出的反应，由此，中子谱的热分量被抑制。

在反应堆腔空隙中的堆芯中平面处的4个方位角位置(即0°、15°、30°和45°)获得测量的响应。因为2-环路PWR反应堆通常特征在于0°位置处的峰值快中子影响，由于核燃料与RPV极接近，因此在该位置处获得附加的测量。具体来说，在0°处，在活动堆芯的顶部和底部处在轴向上获得测量。最初，发现所计算的剂量测定响应一贯地过预测(over-predict)了所测量的数据。

进一步的研究揭示：在输运模型中使用的RPV厚度小于在反应堆压力容器(RPV)的使用中检测(ISI)期间所测量的RPV厚度，证实了最初的发现。将新的RPV厚度引入输运模型，计算的剂量测定数据的准确度提高了平均8％。图2中示出了使用DTW自适应差分方案计算的剂量测定数据的测量值对计算值(M/C)的比率。如图2中所示，在每个位置并且针对每种剂量测定材料，M/C比率一直在10％范围内变化。可以通过使用在30°和45°方位角这些位置处的非均匀网格细化来减少在30°和45°方位角位置处的过预测，其中，系统的曲率变得更相关。在所有剂量测定位置处的平均M/C比率是0.96。

图3a到图3c示出在有ISI校正的RPV厚度的情况下获得的M/C值与在没有厚度校正的情况下获得的值的比较。对位于0°方位角位置处的所有剂量测定样本进行比较并呈现该比较。使用ISI测量的校正的RPV厚度改进了在每个剂量测定位置处所计算的响应的准确度。在15°、30°和45°方位角位置也获得相似的结果。

在表2中示出了在每个剂量测定位置处的测量的和计算的反应率，以及使用利用ISI测量所校正的RPV厚度的M/C比率。在表2中列出的在所有剂量测定位置上的反应的平均M/C比率是0.96。

表2.利用DTW差分方案获得测量的和计算的反应率