论文部分内容阅读
摘 要：发展新一代反应堆物理数值计算理论和方法的目标是实现三维全堆芯Pin-by-Pin输运计算，MOC、SN、SPN是其中三种最具有代表意义的全堆芯中子输运计算方法。该研究介绍了这三种全堆芯中子输运计算方法的概念及特点，通过调研国内外研究现状比较分析了三种计算方法的计算精度、计算效率及几何适应性，总结了三种全堆芯输运方法的优缺点。最终归纳了两种具备工程应用可行性的全堆芯中子输运计算方案：（1）二维MOC耦合一维输运/节块法直接进行全堆芯输运计算；（2）采用二维MOC进行栅格计算，预制参数截面库，再运用三维SPN方法进行全堆芯输运计算。
关键词：全堆芯中子输运 MOC SN SPN
中图分类号：TL329+.2 文献标识码：A 文章编号：1672-3791（2015）07（c）-0001-04
近年来，我国大力发展第四代反应堆，多样的堆芯环境给核设计带来了更复杂的中子学特性和非均匀性，并对中子输运方法的堆芯分析能力提出了更高的要求，主要体现在计算精度、计算效率、处理几何复杂度等方面[1]。
新一代反应堆物理数值计算理论和方法的目标是实现三维全堆芯Pin-by-Pin 输运计算，精确求解任意几何结构的中子输运问题，避免任何均匀化近似；具备强大的几何处理能力，能适应复杂中子输运问题，可以并行计算[2]。目前，国内外开展了大量研究工作，并开发了多个全堆芯输运计算程序。主要形成了以下三种最具代表意义的全堆芯Pin-by-Pin 输运计算方法：①以非均匀栅元为基础的三维特征线法（Method of Characteristics，MOC）；②以均匀栅元为基础的三维离散纵标法（Discrete Ordinate Method，SN）。③以均匀栅元为基础的三维简化球谐函数法（Simplified Spherical Harmonics Method，SPN）。
1 特征线法（MOC）
特征线法（MOC）从中子输运方程的微分形式出发，沿着中子运行轨迹（特征线）进行积分求解。从理论上说，其整个求解过程不受几何形状的限制，无需做均匀化处理，能准确处理强各向异性问题，是全堆芯Pin-by-Pin 输运计算的理想方法[3]。目前，MOC已逐渐成为许多程序包中的标准求解器。表1列出了部分基于MOC开发的程序。
MOC方法采用特征线扫描方式求解中子输运方程，通常采用平源近似，在求解过程中无需保存角中子通量和有关系数矩阵，比较节省内存[4]。几何处理大量采用模块化特征线追踪方法，重复结构共用一套特征线信息，从而节省大量存储空间，提高计算效率。因空间、能群、角度离散，各条特征线的求解相对独立，所以MOC具有天生的并行性。
但是MOC为取得高计算精度，必需采用细致的特征线和平源近似区划分，导致求解复杂问题的速度缓慢，在三维计算中更为突出，需要研究加速并行。其次，MOC在处理三维模型的反射边界条件时，存在入射线和反射线空间或角度错位的问题[5]。再次，模块化特征线追踪方法仅适用于直线规则边界问题，限制了求解问题的任意性。
直接采用MOC方法进行三维复杂中子输运问题计算，速度缓慢且占用大量存储空间。考虑堆芯的非均匀性主要来源于堆芯径向布置，轴向布置相对均匀，国际上提出了一种“退化”的三维MOC方法，即径向二维模块化MOC耦合轴向一维扩散或输运方法。此方法大大提升了计算效率，同时可以保持较高的计算精度。需要提到的是，若轴向一维方向采用扩散方法，在计算控制棒移动等工况时会出现一定误差。基于此方法开发的程序主要有DeCART[6]、Tiger-3D[7]、MOCHA_2D1D[8]等。
2 离散纵标法（SN）
离散纵标法（SN）将方向自变量离散，每个方向对应一定的立体角，所有立体角之和为4π，在特定方向上求解中子输运方程，并假定立体角内的角通量是常数，采用求积代替积分，得到中子标通量。SN在反应堆屏蔽计算的中子及光子深穿透问题、压力壳辐照损伤分析中得到了广泛应用[3]。表2列出了一些基于SN开发的程序。
SN方法每个离散方向的方程相对独立，应用迭代法求解时，数值过程比较简单，便于编写通用程序。相比于扩散近似、球谐函数等近似方法，SN对边界条件的离散精确自洽，便于处理真空边界。
但是为得到足够的精度，需要足够细密的离散点，这导致计算耗用大量内存和时间，尤其对于离散方向数较大的三维问题，挑战更为严峻。其次，SN通常采用简化方法来处理几何和源模型，使其在处理复杂几何及精确源模拟上存在严重缺陷[3]。再次，SN在进行多维问题计算时，由于角度离散本身固有的缺陷，容易出现“射线效应”并引起假散射，导致计算误差。SN通常采用菱形差分计算模型，这往往导致计算出现与物理不符的负中子通量密度，需要采用置零修正，这也进一步影响了计算精度。
3 简化球谐函数法（SPN）
球谐函数方法（PN）针对输运方程中的角度自变量做近似，将中子通量密度用球谐函数展开成N阶级数，将中子输运方程化成相互耦合微分方程组，并确定级数中的系数，再采用迭代策略进行求解。但是当阶数N增大时，尤其对于多维度复杂几何，PN方法将变得非常复杂和困难[9]。简化球谐函数法（SPN）在一维平板几何下导出球谐函数方程组，将一维空间微分算子直接换成三维空间微分算子得到三维直角坐标系下的SPN方程组。此方法由PN方法适当化简而来，减小了方程的复杂性，保持了较高的精度[10]。表3列出了一些SPN程序。
SPN方法仅对空间位置离散，对角度变量的处理是连续的，具有角度旋转不变性，不存在射线效应；可以采用有限元方法，特别适用于不规则几何的非结构网格计算[11]。相比高阶PN和SN方法，SPN更容易在满足精度的前提下节省计算时间和资源。
但是，SPN方法在栅元均匀化的基础上实施，需要进行计算网格离散，通过SPH因子修正均匀化截面，预制参数截面库，这些均给全堆芯pin-by-pin计算引入了误差。其次，SPN方法作为低阶输运近似本身具有一定误差。 4 MOC、SN、SPN对比分析
三种全堆芯输运计算方法各有优势，但又都存在不同程度的缺点。国际上针对三种方法开展了一些对比研究，主要分析了计算精度、计算效率、几何适应性等。
日本名古屋大学[12]通过BWR组件计算比较了MOC、SPN、扩散方法的计算精度和计算时间。结果表明：在采用细致网格进行均匀栅元的组件计算时，经SPH因子修正后的SPN在计算精度方面与MOC计相当，比扩散方法高2倍；在计算时间方面，比MOC快很多，仅为扩散方法的1.5倍。SPN与基于非均匀栅元的MOC相比，Keff的误差约为0.02%，反应率的误差约为0.5%。另一方面，考虑到第四代反应堆中包含不同形式的组件，SPN方法必须具备处理交错网格的能力，才有工程上的实用性。但目前尚未有人开展此方面的工作。
加拿大蒙特利尔大学以OECD/NEA发布的NEA3D-TAB-2007中子输运基准题为例，分析比较了DRAGON程序中三维MOC方法和SN方法对空间参数剧烈变化问题的计算性能，结果表明：MOC和SN方法均具备较好的计算精度，但SN方法在非均匀性较强的介质中具有明显的射线效应；而MOC方法在平源近似情况下，角度离散和网格划分数量的不足将导致计算误差急剧增加，且在同等的计算条件下，MOC方法的计算时间远超于SN方法。
5 MOC、SN、SPN联合研究
考虑到三种方法各有优劣，因此可以取长补短，采用两种方法联合进行全堆芯pin-by-pin计算。目前国际上开展的相关主要研究有：
日本核工程公司、名古屋大学和核燃料工业公司合作开发了PWR堆芯燃料管理系统程序包AEGIS/SCOPE2[12]，其中AEGIS为栅格输运计算程序，采用MOC方法计算得到栅元均匀化截面并进行SPH因子修正；SCOPE2为堆芯程序，采用近似输运的节块SPN方法。AEGIS/SCOPE2采用MOC加SPN的“两步法”进行全堆芯输运计算。计算的Keff与实验值相差152 pcm，棒功率相差1%以内。
2005年，欧盟提出了NURESIM计划[8]，旨在研究更先进、精确的堆芯计算方法，以满足四代反应堆的设计计算要求。该计划提出了两条技术路线：先进蒙卡方法和确定论方法研究。其中先进确定论方法的研究路线为：采用二维MOC进行组件栅格计算，开发了APOLLO2、HELIOS-2等程序；采用三维节块SPN进行堆芯计算，研制了DYN3D程序等；并以蒙卡程序Serpent为基准，验证了HELIOS-2/DYN3D系统计算分析钠冷快堆的准确性[10]。
从以上研究情况和国际发展趋势来看，采用MOC加SPN的“两步法”进行全堆芯输运计算可以得取得足够的计算精度，计算代价也在可接受范围内。因此，该方案可作为工程应用的选择之一。
6 结语
该研究介绍了MOC、SN、SPN三种全堆芯中子输运计算方法的概念及特点，通过国内外研究现状分析对比了三种方法的全堆芯输运计算能力，得到以下结论：①计算精度：MOC不存在均匀化，精度最高。SN存在射线效应以及计算出负，计算精度略低。SPN需栅元均匀化，在细致网格条件下精度接近输运计算，但在非均匀性强的区域误差偏大。②计算效率：MOC计算时间最长，占用大量内存，采用模块化MOC可提高效率，节省存储空间。SN同样需耗用较大的计算内存和时间。SPN计算时间最短，计算效率最高。③几何适应性：MOC理论上适用于任意几何，若采用模块化MOC则限制了几何适应能力。SN和有限元SPN方法同样具备很强的几何适应性。但由于计算耗时、耗内存，MOC、SN还多限于二维问题计算。
三种方法取长补短，已形成两套可行的全堆芯pin-by-pin输运计算方案，可为今后全堆芯输运计算方法的工程应用提供参考依据：①二维MOC耦合一维输运/节块法直接进行全堆芯计算，如韩国DeCART程序。②二维MOC进行栅格计算，预制参数截面库，再采用三维SPN进行全堆芯计算，如日本AEGIS/SCOPE2、欧盟APOLLP2/DYN3D程序。
参考文献
[1] 陈珍平，王电喜，何桃，等.基于 CAD 技术的特征线中子输运计算程序开发[J].核科学与工程，2013，32（4）：354-359.
[2] 汤春桃.中子输运方程特征线解法及嵌入式组件均匀化方法的研究[D].上海：上海交通大学，2009.
[3] 韩静茹.三维蒙特卡罗—离散纵标双向耦合屏蔽计算方法研究[D].北京：华北电力大学，2012.
[4] 汤春桃.线性源近似的中子输运方程特征线解法[J].原子能科学技术，2011，45（12）：1414 -1420.
[5] 柴晓明，姚栋，王侃.基于特征线方法的三维中子输运程序（Ⅰ）—边界条件的插值处理[J].核动力工程，2010（2）：11-15.
[6] Cho J Y，Kim K S，Lee C C，et al.Axial SPN and radial MOC coupled whole core transport calculation[J].Journal of Nuclear Science and Technology，2007，44（9）：1156-1171.
[7] 吴文斌，李庆，王侃.基于模块化射线追踪的矩阵 MOC方法 （1）—理论研究[J].核动力工程，2014，35（3）：129-133.
[8] 梁亮，吴宏春，郑友琦，等.基于模块化特征线方法的二维/一维耦合输运程序开发[J].核动力工程，2014（S2）：131-134.
[9] 马文娟.基于辐射传输方程高阶球谐近似模型的时域 DOT/FDOT 成像方法研究[D].天津：天津大学，2012.
[10] 曹良志，吴宏春.非结构网格中子输运方程的球谐函数解法研究[J].核动力工程，2004，25（5）： 395-398.
[11] 黄细珍.辐射传递的球谐函数法模拟[D].哈尔滨：哈尔滨工业大学，2010.
[12] Tatsumi M，Ohoka Y，Endo T，et al.sVerification and Validation of AEGIS/SCOPE2：The State-of-the-Art In-Core Fuel Management System for PWRs[C].Xi’an：18th International Conference on Nuclear Engineering，2010：9-15.
关键词：全堆芯中子输运 MOC SN SPN
中图分类号：TL329+.2 文献标识码：A 文章编号：1672-3791（2015）07（c）-0001-04
近年来，我国大力发展第四代反应堆，多样的堆芯环境给核设计带来了更复杂的中子学特性和非均匀性，并对中子输运方法的堆芯分析能力提出了更高的要求，主要体现在计算精度、计算效率、处理几何复杂度等方面[1]。
新一代反应堆物理数值计算理论和方法的目标是实现三维全堆芯Pin-by-Pin 输运计算，精确求解任意几何结构的中子输运问题，避免任何均匀化近似；具备强大的几何处理能力，能适应复杂中子输运问题，可以并行计算[2]。目前，国内外开展了大量研究工作，并开发了多个全堆芯输运计算程序。主要形成了以下三种最具代表意义的全堆芯Pin-by-Pin 输运计算方法：①以非均匀栅元为基础的三维特征线法（Method of Characteristics，MOC）；②以均匀栅元为基础的三维离散纵标法（Discrete Ordinate Method，SN）。③以均匀栅元为基础的三维简化球谐函数法（Simplified Spherical Harmonics Method，SPN）。
1 特征线法（MOC）
特征线法（MOC）从中子输运方程的微分形式出发，沿着中子运行轨迹（特征线）进行积分求解。从理论上说，其整个求解过程不受几何形状的限制，无需做均匀化处理，能准确处理强各向异性问题，是全堆芯Pin-by-Pin 输运计算的理想方法[3]。目前，MOC已逐渐成为许多程序包中的标准求解器。表1列出了部分基于MOC开发的程序。
MOC方法采用特征线扫描方式求解中子输运方程，通常采用平源近似，在求解过程中无需保存角中子通量和有关系数矩阵，比较节省内存[4]。几何处理大量采用模块化特征线追踪方法，重复结构共用一套特征线信息，从而节省大量存储空间，提高计算效率。因空间、能群、角度离散，各条特征线的求解相对独立，所以MOC具有天生的并行性。
但是MOC为取得高计算精度，必需采用细致的特征线和平源近似区划分，导致求解复杂问题的速度缓慢，在三维计算中更为突出，需要研究加速并行。其次，MOC在处理三维模型的反射边界条件时，存在入射线和反射线空间或角度错位的问题[5]。再次，模块化特征线追踪方法仅适用于直线规则边界问题，限制了求解问题的任意性。
直接采用MOC方法进行三维复杂中子输运问题计算，速度缓慢且占用大量存储空间。考虑堆芯的非均匀性主要来源于堆芯径向布置，轴向布置相对均匀，国际上提出了一种“退化”的三维MOC方法，即径向二维模块化MOC耦合轴向一维扩散或输运方法。此方法大大提升了计算效率，同时可以保持较高的计算精度。需要提到的是，若轴向一维方向采用扩散方法，在计算控制棒移动等工况时会出现一定误差。基于此方法开发的程序主要有DeCART[6]、Tiger-3D[7]、MOCHA_2D1D[8]等。
2 离散纵标法（SN）
离散纵标法（SN）将方向自变量离散，每个方向对应一定的立体角，所有立体角之和为4π，在特定方向上求解中子输运方程，并假定立体角内的角通量是常数，采用求积代替积分，得到中子标通量。SN在反应堆屏蔽计算的中子及光子深穿透问题、压力壳辐照损伤分析中得到了广泛应用[3]。表2列出了一些基于SN开发的程序。
SN方法每个离散方向的方程相对独立，应用迭代法求解时，数值过程比较简单，便于编写通用程序。相比于扩散近似、球谐函数等近似方法，SN对边界条件的离散精确自洽，便于处理真空边界。
但是为得到足够的精度，需要足够细密的离散点，这导致计算耗用大量内存和时间，尤其对于离散方向数较大的三维问题，挑战更为严峻。其次，SN通常采用简化方法来处理几何和源模型，使其在处理复杂几何及精确源模拟上存在严重缺陷[3]。再次，SN在进行多维问题计算时，由于角度离散本身固有的缺陷，容易出现“射线效应”并引起假散射，导致计算误差。SN通常采用菱形差分计算模型，这往往导致计算出现与物理不符的负中子通量密度，需要采用置零修正，这也进一步影响了计算精度。
3 简化球谐函数法（SPN）
球谐函数方法（PN）针对输运方程中的角度自变量做近似，将中子通量密度用球谐函数展开成N阶级数，将中子输运方程化成相互耦合微分方程组，并确定级数中的系数，再采用迭代策略进行求解。但是当阶数N增大时，尤其对于多维度复杂几何，PN方法将变得非常复杂和困难[9]。简化球谐函数法（SPN）在一维平板几何下导出球谐函数方程组，将一维空间微分算子直接换成三维空间微分算子得到三维直角坐标系下的SPN方程组。此方法由PN方法适当化简而来，减小了方程的复杂性，保持了较高的精度[10]。表3列出了一些SPN程序。
SPN方法仅对空间位置离散，对角度变量的处理是连续的，具有角度旋转不变性，不存在射线效应；可以采用有限元方法，特别适用于不规则几何的非结构网格计算[11]。相比高阶PN和SN方法，SPN更容易在满足精度的前提下节省计算时间和资源。
但是，SPN方法在栅元均匀化的基础上实施，需要进行计算网格离散，通过SPH因子修正均匀化截面，预制参数截面库，这些均给全堆芯pin-by-pin计算引入了误差。其次，SPN方法作为低阶输运近似本身具有一定误差。 4 MOC、SN、SPN对比分析
三种全堆芯输运计算方法各有优势，但又都存在不同程度的缺点。国际上针对三种方法开展了一些对比研究，主要分析了计算精度、计算效率、几何适应性等。
日本名古屋大学[12]通过BWR组件计算比较了MOC、SPN、扩散方法的计算精度和计算时间。结果表明：在采用细致网格进行均匀栅元的组件计算时，经SPH因子修正后的SPN在计算精度方面与MOC计相当，比扩散方法高2倍；在计算时间方面，比MOC快很多，仅为扩散方法的1.5倍。SPN与基于非均匀栅元的MOC相比，Keff的误差约为0.02%，反应率的误差约为0.5%。另一方面，考虑到第四代反应堆中包含不同形式的组件，SPN方法必须具备处理交错网格的能力，才有工程上的实用性。但目前尚未有人开展此方面的工作。
加拿大蒙特利尔大学以OECD/NEA发布的NEA3D-TAB-2007中子输运基准题为例，分析比较了DRAGON程序中三维MOC方法和SN方法对空间参数剧烈变化问题的计算性能，结果表明：MOC和SN方法均具备较好的计算精度，但SN方法在非均匀性较强的介质中具有明显的射线效应；而MOC方法在平源近似情况下，角度离散和网格划分数量的不足将导致计算误差急剧增加，且在同等的计算条件下，MOC方法的计算时间远超于SN方法。
5 MOC、SN、SPN联合研究
考虑到三种方法各有优劣，因此可以取长补短，采用两种方法联合进行全堆芯pin-by-pin计算。目前国际上开展的相关主要研究有：
日本核工程公司、名古屋大学和核燃料工业公司合作开发了PWR堆芯燃料管理系统程序包AEGIS/SCOPE2[12]，其中AEGIS为栅格输运计算程序，采用MOC方法计算得到栅元均匀化截面并进行SPH因子修正；SCOPE2为堆芯程序，采用近似输运的节块SPN方法。AEGIS/SCOPE2采用MOC加SPN的“两步法”进行全堆芯输运计算。计算的Keff与实验值相差152 pcm，棒功率相差1%以内。
2005年，欧盟提出了NURESIM计划[8]，旨在研究更先进、精确的堆芯计算方法，以满足四代反应堆的设计计算要求。该计划提出了两条技术路线：先进蒙卡方法和确定论方法研究。其中先进确定论方法的研究路线为：采用二维MOC进行组件栅格计算，开发了APOLLO2、HELIOS-2等程序；采用三维节块SPN进行堆芯计算，研制了DYN3D程序等；并以蒙卡程序Serpent为基准，验证了HELIOS-2/DYN3D系统计算分析钠冷快堆的准确性[10]。
从以上研究情况和国际发展趋势来看，采用MOC加SPN的“两步法”进行全堆芯输运计算可以得取得足够的计算精度，计算代价也在可接受范围内。因此，该方案可作为工程应用的选择之一。
6 结语
该研究介绍了MOC、SN、SPN三种全堆芯中子输运计算方法的概念及特点，通过国内外研究现状分析对比了三种方法的全堆芯输运计算能力，得到以下结论：①计算精度：MOC不存在均匀化，精度最高。SN存在射线效应以及计算出负，计算精度略低。SPN需栅元均匀化，在细致网格条件下精度接近输运计算，但在非均匀性强的区域误差偏大。②计算效率：MOC计算时间最长，占用大量内存，采用模块化MOC可提高效率，节省存储空间。SN同样需耗用较大的计算内存和时间。SPN计算时间最短，计算效率最高。③几何适应性：MOC理论上适用于任意几何，若采用模块化MOC则限制了几何适应能力。SN和有限元SPN方法同样具备很强的几何适应性。但由于计算耗时、耗内存，MOC、SN还多限于二维问题计算。
三种方法取长补短，已形成两套可行的全堆芯pin-by-pin输运计算方案，可为今后全堆芯输运计算方法的工程应用提供参考依据：①二维MOC耦合一维输运/节块法直接进行全堆芯计算，如韩国DeCART程序。②二维MOC进行栅格计算，预制参数截面库，再采用三维SPN进行全堆芯计算，如日本AEGIS/SCOPE2、欧盟APOLLP2/DYN3D程序。
参考文献
[1] 陈珍平，王电喜，何桃，等.基于 CAD 技术的特征线中子输运计算程序开发[J].核科学与工程，2013，32（4）：354-359.
[2] 汤春桃.中子输运方程特征线解法及嵌入式组件均匀化方法的研究[D].上海：上海交通大学，2009.
[3] 韩静茹.三维蒙特卡罗—离散纵标双向耦合屏蔽计算方法研究[D].北京：华北电力大学，2012.
[4] 汤春桃.线性源近似的中子输运方程特征线解法[J].原子能科学技术，2011，45（12）：1414 -1420.
[5] 柴晓明，姚栋，王侃.基于特征线方法的三维中子输运程序（Ⅰ）—边界条件的插值处理[J].核动力工程，2010（2）：11-15.
[6] Cho J Y，Kim K S，Lee C C，et al.Axial SPN and radial MOC coupled whole core transport calculation[J].Journal of Nuclear Science and Technology，2007，44（9）：1156-1171.
[7] 吴文斌，李庆，王侃.基于模块化射线追踪的矩阵 MOC方法 （1）—理论研究[J].核动力工程，2014，35（3）：129-133.
[8] 梁亮，吴宏春，郑友琦，等.基于模块化特征线方法的二维/一维耦合输运程序开发[J].核动力工程，2014（S2）：131-134.
[9] 马文娟.基于辐射传输方程高阶球谐近似模型的时域 DOT/FDOT 成像方法研究[D].天津：天津大学，2012.
[10] 曹良志，吴宏春.非结构网格中子输运方程的球谐函数解法研究[J].核动力工程，2004，25（5）： 395-398.
[11] 黄细珍.辐射传递的球谐函数法模拟[D].哈尔滨：哈尔滨工业大学，2010.
[12] Tatsumi M，Ohoka Y，Endo T，et al.sVerification and Validation of AEGIS/SCOPE2：The State-of-the-Art In-Core Fuel Management System for PWRs[C].Xi’an：18th International Conference on Nuclear Engineering，2010：9-15.