第1章引言
1.1研究背景与意义
能源与环境是人类生存与发展的必要条件，能源的有效利用和环境的保护是当今社会关注的两大问题。中国目前是世界能源消费第一大国，化石能源在我国能源结构中占据了很大比重，煤炭更是占到63.7%［1］。面对我国社会主要矛盾的变化，环境保护的需求在人民对美好生活的追求中也显得更加重要。面对能源与环境这两大问题，核电将扮演不可或缺的角色。首先，核电能量密度大、运行稳定，相较于太阳能和水电等新能源，核电受季节及气候影响很小，可以保障能源应用过程的稳定和安全。同时，核电可以有效减少温室气体、粉尘等的排放，在减排环保方面能够产生重大效益。然而2016年数据显示我国核能发电占比不足4%［2］，福岛事故后我国核电发展放缓，社会对核电的经济性和安全性也有一些质疑的观点。

如何从设计上提高核电厂的安全性和经济性是核电发展的关键，核能科研工作者进行了大量研究工作。一方面，随着计算能力的提升，人们对反应堆计算程序的精度、效率和计算能力提出了更高的要求，传统的计算方法和程序已不能完全满足
对安全性和经济性方面的新需求。近年来，高保真(highfidelity)计算的理念受到越来越广泛的关注。高保真计算旨在减少计算程序在建模及模拟方法中对反应堆真实情况的近似，从而减少由计算方法导致的不确定性和保守假设，提高核电厂的安全性和经济性。
另一方面，新型反应堆设计也相继被提出，某些新型反应堆已逐步进入试验乃至工程阶段，其中具有代表性的新堆型包括： 快中子堆、球床高温气冷堆、熔盐堆、小型模块化反应堆等。另外，对压水堆新型燃料的研究，特别是
事故容错燃料(accident tolerant fuel)的研究也在如火如荼的进行当中。这些新型堆和新型燃料设计理念的提出，也对反应堆计算程序和方法提出了新的要求。

面对反应堆高保真计算的需求，以及新型堆和新型燃料的研究设计需求，蒙特卡罗方法由于其几何建模的准确性和灵活性、连续能量点截面处理的准确性以及良好的并行性，受到越来越多反应堆计算科研工作者的关注，并发挥着越来越重要的作用。世界上很多国家都研发了自主的蒙卡程序，其中美国有最经典的蒙卡程序MCNP［3］、美国海军核实验室的MC21［4］以及MIT的开源蒙卡程序OpenMC［5］。其他程序包括法国的TRIPOLI［6］、日本的MVP［7］、韩国的McCARD［8］和芬兰的Serpent［9］等。国内的一些高校及研究单位，也开发了一些自主化蒙卡程序，包括北京应用物理与计算数学研究所的JMCT［10］、中科院核能安全技术研究所的SuperMC［11］以及清华大学
工程物理系的堆用蒙卡程序RMC［12］。

随着计算能力的提高以及先进并行算法的提出，蒙卡方法和程序在全堆大规模计算方面的能力得到很大提升。然而，为了使蒙卡程序
能够被应用于实际的反应堆工程设计及研究当中，必须不断拓展蒙卡程序的分析功能，从而
使其能够应对反应堆高保真计算以及新型堆和新型燃料设计的新需求。其中，反应堆全堆全寿期高保真耦合模拟和随机介质精细输运及燃耗计算是当前蒙卡程序研究的热点和难点。

反应堆全堆全寿期高保真耦合模拟不同于传统反应堆计算中两步和三步法等均匀化的思路，其同时采用多物理(中子输运、热工水力、燃耗等)的耦合，旨在减少传统反应堆计算中的栅元/组件与堆芯分步计算及各种物理场独立处理的假设，从而减少近似，提高计算精度。在PHYSOR 2010会议上，Bill Martin和Forrest Brown在题为“Some Challenges for Largescale Reactor Calculations”的大会报告中，提出了将蒙卡方法用于大规模计算所面临的五大问题(见表1.1)，其中就包括蒙卡燃耗计算与蒙卡多物理耦合计算。在“Technical Summary of PHYSOR 2014”中，Akio Yamamoto也提出了蒙卡方法的五大未解决问题(见表1.2)，“Instability of MC based multiphysics calculation due to statistics deviation”也是其中之一。


表1.1PHYSOR 2010会议提出的将蒙卡方法用于大规模计算面临的五大问题



序号问题序号问题

1计算时间与内存4适应未来计算机
2蒙卡燃耗5多物理耦合
3源收敛


表1.2PHYSOR 2014会议提出的蒙卡方法的五大未解决问题



序号问题

1局部统计量的收敛及相关性
2基于连续能量蒙卡的高效、节省内存的敏感性和不确定性评价方法
3对多种形状/大小的随机分布颗粒的处理
4由于统计偏差造成的基于蒙卡的多物理耦合的不稳定性
5多群截面产生及其不确定度在全堆计算中的影响

“Technical Summary of PHYSOR 2014”提到蒙卡方法的五大挑战之一是“Treatment of random distributions of particles with various shapes and sizes”。这一话题之所以成为蒙卡方法研究的重要方向之一，在于其特殊的应用对象——弥散燃料。包覆颗粒弥散型燃料具有在高温及深燃耗条件下阻滞和包容裂变产物的能力，能够在保证安全性的前提下提高燃料利用的经济性，因此弥散燃料得到了越来越广泛的应用。弥散燃料及其所构成的堆芯依据材料结构的属性分类属于随机非均匀介质，精确高效的随机介质中子输运燃耗计算，是基于弥散燃料的新型燃料及反应堆设计的核心和基础。

结合目前蒙卡方法发展的热点和难点，以及蒙卡方法在反应堆设计中的需求和应用，本书基于自主反应堆蒙卡程序RMC，开展蒙卡全寿期高保真耦合模拟与随机介质精细计算研究。
1.2研 究 现 状
1.2.1蒙卡全寿期高保真耦合模拟研究现状

蒙卡方法与确定论方法并称为反应堆中子输运计算的两类重要方法。由于其在几何和能量处理方面具有灵活性和准确性，
使蒙卡方法自身就具有高保真方法的特点。由于蒙卡方法在计算效率和计算能力方面受到限制，长期以来蒙卡方法主要作为初装堆分析和屏蔽计算中重要的中子输运求解器。

然而，核反应堆是一个不同物理场相互耦合和作用的复杂系统，除了中子输运以外，还需要考虑热工水力(温度、密度)、核素燃耗、燃料性能、水化学和反应性控制(可溶硼与控制棒)等多方面的物理现象。同时，在核反应堆的多循环燃耗计算中，还需要考虑燃耗后燃料组件的倒换料问题，如图1.1所示。


图1.1数值化反应堆的基础组成部分



全寿期高保真耦合模拟对反应堆功率提升、燃耗的加深以及反应堆延寿有重要的意义，有助于提高反应堆的经济性。同时通过全寿期高保真耦合模拟，可以构建“数值化反应堆”，为反应堆工程设计和分析(如小型堆等)提供工具，也有助于高性能科学计算及相关基础科学研究的发展。
基于蒙卡方法或确定论方法为中子输运求解器的高保真多物理耦合，受到全世界研究者的关注。其中，最著名的是美国能源部的CASL计划［13］。CASL是先进轻水堆模拟联盟(Consortium for Advanced Simulation of LWRs)的简称，该联盟由美国的各大国家实验室、高校以及电力公司组成，致力于构建一个虚拟反应堆并进行与该虚拟反应堆相关的基础科学研究。
VERA是CASL计划的产物和核心［14］，是针对反应堆应用的高保真虚拟环境(
highfidelity virtual environment for reactor applications)的简称。VERA的中子输运求解器是基于特征线方法(MOC)的MPACT，热工水力程序采用子通道程序CTF。

基于确定论程序的高保真耦合程序还有韩国首尔国立大学开发的MOC程序nTRACER，其与子通道程序MATRA进行耦合［15］。基于蒙卡程序的高保真耦合包括MC21/CTF耦合［16］、MCNP/SUBCHANFLOW耦合［17］和Serpent 2/SUBCHANFLOW耦合［18］等。可以看出，目前高保真耦合程序大多数采用子通道程序作为热工水力程序。

国内研究方面，西安交通大学开发的NECPX/SUBSC是基于确定论程序的高保真耦合的代表程序［19］。NECPX也采用了2D/1D特征线方法，同时在共振处理及各向异性泄漏处理方面采用了独特的先进方法。而在蒙卡高保真耦合方面，北京应用物理与计算数学研究所也开展了JMCT与子通道程序COBRAEN的物理热工耦合研究［20］。

同时，为了验证这些高保真耦合程序的计算能力和准确性，国际上提出了一些著名的基准题。包括CASL计划的VERA系列基准题［21］和MIT提出的BEAVRS基准题［22］等。这两个基准题都是基于美国西屋公司核电站实际参数和实测数据提出的，是考验反应堆设计与分析程序的重要验证标准。这两个基准题从技术层次上可以划分为三个阶段： 热态零功率(HZP)、热态满功率首循环(HFP+Cycle 1)以及热态满功率多循环(HFP+Multicycle)。以BEAVRS基准题为例，完成第一阶段
热态零功率(HZP)计算的蒙卡程序有美国的OpenMC［22］和MC21［23］，日本的MVP［24］，韩国的McCARD［25］，中国的JMCT［26］、SuperMC［27］和RMC［28］。西安交通大学的NECPCACTI(组件程序)和NECPVIOLET(堆芯扩散程序)也完成了基于两步法的BEAVRS基准题HZP计算［29］。在
本研究初期，只有美国的MC21程序在2014年完成了1/4堆的HFP首循环计算［30］，是当时唯一完成第二阶段
热态满功率首循环(HFP+Cycle1)计算的蒙卡程序。2016年，芬兰VTT的蒙卡程序Serpent和扩散程序ARES也完成了基于两步法的BEAVRS基准题HFP首循环计算［31］。2017年，西安交通大学的陈定勇等也基于CASMO4E/SIMULATE3进行了BEAVRS首循环计算［32］。到了第三阶段
热态满功率多循环(HFP+Multicycle)，在研究初期世界范围内尚未有蒙卡程序完成BEAVRS的两循环计算，而在确定论程序方面，完成第三阶段计算的有韩国的nTRACER ［15］(2014年)和CASL的VERA程序(2017年)，这两个程序都采用了高保真的耦合计算。另外MIT的GA Gunow在2015年完成了基于CASMOSIMULATE的两步法计算［33］，国家电力投资集团中央研究院的COSINE软件包也在2016年采用确定论的两步法完成了计算。
另外一个著名的VERA系列基准题总共有10个子问题，如图1.2所示。世界范围内完成到第7个问题(全堆热态满功率问题)的确定论程序有CASL的VERA程序，蒙卡程序有美国的MC21/CTF程序。另外，西安交通大学的NECPX/SUBSC在2017年也完成了VERA系列基准题的第6个问题(组件热态满功率问题)［19］。


图1.2VERA基准题的10个子问题



这些高保真耦合基准题，是对蒙卡程序进行反应堆模拟及设计所需要的各种计算功能的考验，其中的关键算法如图1.3所示。


图1.3蒙卡全寿期耦合计算关键算法



综上所述，国内外针对高保真耦合计算的研究方兴未艾，采用高保真方法完成全寿期热态满功率计算的确定论和蒙卡程序
皆尚属少数。基于蒙卡程序的全寿期高保真耦合模拟更是世界范围内的难题。
1.2.2蒙卡随机介质精细计算研究现状
燃料设计是提高新型核能系统安全性和经济性的关键因素，弥散型燃料是其中引人关注的燃料类型之一。弥散型燃料是将核燃料弥散分布在非裂变材料(金属和石墨)中。弥散型燃料由于燃料颗粒细小，受辐照后可以达到很高的温度。同时裂变产物基本上都被包容在燃料颗粒内，基体也起到了进一步包容裂变碎片的作用，因此可以达到很深的燃耗深度。

随着工艺技术的不断成熟，包覆颗粒弥散型燃料越来越多地被应用于高温气冷堆、钍基熔盐堆和空间堆等新型先进核能系统。在福岛核事故之后，弥散型燃料作为一种事故容忍材料
(accident tolerant fuels，ATF)得到广泛关注。例如，有研究者提出在传统轻水堆中引入新型全陶瓷包覆颗粒燃料
(fully ceramic microencapsulated fuel，FCM)［34］。此外，随着降浓计划(reduced enrichment for research and test 
reactors，RERTR)的实施，弥散型燃料被广泛应用于各种实验研究堆［35］。同时弥散燃料也是核动力舰船堆芯燃料的核心，随机非均匀介质的精细输运和燃耗计算，与核动力舰船的设计、运行和安全等紧密相关。

弥散型燃料及由其构成的堆芯属于随机非均匀介质［36］。一方面，区别于传统的陶瓷型或金属型核燃料，弥散型燃料具有随机非均匀结构，即体积微小且数量巨大的燃料颗粒随机弥散在基体材料当中。另一方面，弥散型燃料元件在堆芯内还可能按照确定的(如棱柱或平板结构)或随机的(如球床结构)形式非均匀布置。以球床式高温气冷堆为例，堆芯内包含大量随机堆积的燃料球和石墨球，每个燃料球内包含大量随机分布的燃料颗粒，即所谓的双重随机非均匀性。

为了解决这种随机非均匀性，确定论程序和蒙卡程序中都开发了不同的方法。组件均匀化程序Dragon中开发了针对双重非均匀性的碰撞概率方法［37］，但是只能处理二维几何问题。西安交通大学的刘庆杰等采用统计与确定论方法相结合进行求解，通过蒙卡方法得到平均弦长，再用SN方法求解统计输运方程［38］。该方法对一些测试问题取得了较好结果，然而随着材料弦径比减小，误差将增大。球床式高温气冷堆物理设计所广泛使用的VSOP程序［39］，在建模时将全堆划分
为若干个区域，先通过碰撞概率法求解各区域的群常数，然后进行全堆扩散求解。VSOP在计算群常数时，为考虑双重非均匀性的影响，采用了近似物理模型修正群截面和共振自屏等参数。可见，确定论方法在随机介质计算的精度及普适性(多种颗粒类型以及各种体积填充率)上存在一些不足和限制。

另一方面，由于复杂几何具有精细描述、连续能量点截面处理和高效并行等能力，使蒙卡方法在解决随机介质的输运问题中
占据了独特的优势。然而传统的蒙卡程序缺乏精确高效模拟随机介质的能力。以往的基于蒙卡程序的随机介质模拟，一般采用MCNP的重复结构，假设燃料颗粒和燃料球均为规则排列［40, 41］。为达到真实的体积填充率，一些研究人员［42, 43］提出了体心立方和面心立方等排列方式，并人为调整球心间距。重复结构描述方法忽略了随机介质的随机特性，可能会引入不可忽略的系统误差。首先，重复结构会增加特定角度的中子泄漏，导致中子流效应［44］。其次，重复结构网格中的颗粒通常会被外边界截断［45］，偏离真实物理情况。此外，重复结构难以描述多尺寸、多类型随机介质的情况。譬如，在某些弥散型燃料元件/组件设计中，混合填充燃料颗粒和可燃毒物颗粒［46］，用于降低堆芯的初始反应性。对于这种燃料形式，很难通过重复结构准确描述。SCALE程序包中蒙卡程序KENO的多群模块提供了doublehet的功能［47］，采用多群截面预处理器CENTRM/PMC产生问题相关的修正多群截面，其本质也是对截面的等效和近似。

研究精确高效模拟随机介质的连续能量蒙卡中子输运方法，是蒙卡方法研究的热点问题之一。目前国际上已开展的相关工作主要包括以下三种方法。
(1) 重复结构随机栅格方法(random lattice)。该方法由Brown和Martin［48］提出，并应用在MCNP5当中。其基本思想是在几何建模时，先假设随机介质按照重复结构排列，在中子输运过程中，令重复网格内所填充的物质在其原来的位置上发生随机扰动。重复结构随机栅格方法简单易行，但也存在一定的局限： 首先，该方法无法描述紧密堆积结构(如球床)，缺少普适性； 其次，扰动幅度不能过大，否则可能造成燃料颗粒被重复结构内部网格截断； 再次，随机扰动加剧了重复结构网格被外边界截断的不确定性。

(2) 弦长抽样方法(chordlength sampling)。自20世纪90年代直至近几年，Zimmerman和Adams［49］、
Isaomurata等［50］、Ji Wei和Martin［51］以及Liang Chao和Ji Wei［52］对弦长抽样方法开展了大量研究，MVP和Serpent等程序亦采用了该方法。弦长抽样方法无须显式地描述所有随机介质，而是在中子输运过程中对随机介质的距离(即弦长)和角度进行抽样。弦长抽样方法的最大优势在于它能够简化几何建模，提高计算效率。该方法的关键难点在于准确描述弦长的概率密度函数，特别是在高体积填充率时需要对概率密度函数进行修正。

(3) 显式模拟方法(explicit modelling)。重复结构随机扰动方法和弦长抽样方法都属于在线抽样方法，即在中子输运过程中通过在线抽样来获得随机介质的空间位置。而显式模拟方法通过预先抽样获得所有随机介质的位置分布，然后采用常规的蒙卡中子输运方法进行模拟。以球床式高温气冷堆为例，需要预先产生所有燃料球和石墨球的坐标位置，以及每个燃料球内包覆颗粒的坐标位置，进而构建整个堆芯模型。显式模拟方法具有很高的保真度，因而一般用作其他方法的基准解。显式模拟方法的关键技术之一是正确且快速地抽样产生球床或包覆颗粒等随机介质的位置分布。对此国际上已开展的研究包括随机序列增加
(random sequential addition，RSA)方法［53］、准静态方法(quasidynamics method，QDM)［54］和离散单元法(discrete element method，DEM)［55］等。显式模拟方法的另一个关键技术是在蒙卡输运计算过程中提高粒子跟踪效率和减少内存占用。Monk和Serpent程序已实现可用于高温堆的显式模拟方法［56, 57］，但具体细节未见于公开文献。
西安交通大学的李志峰等采用Serpent程序的弦长抽样方法和显式模拟方法对氟盐冷却球床高温堆的燃料栅元进行了计算，并对比了规则几何、弦长抽样方法和显式模拟方法对结果的影响，发现弦长抽样方法在高体积填充率时即使对体积填充率进行修正仍会产生较大误差［58, 59］。然而基于RSA的填充方法的体积填充率上限为38%［60］，无法满足一些高体积填充率的设计(如FCM燃料)的需求。
除了对中子输运问题的模拟外，蒙卡方法研究和程序发展的另一个重要方向是蒙卡燃耗计算。所谓蒙卡燃耗计算，实质上是蒙卡中子输运与点燃耗计算过程的相互耦合。蒙卡输运过程得到中子通量和反应截面等数据，点燃耗过程求解相应的燃耗方程并更新核素密度。

对于使用包覆颗粒弥散型燃料的反应堆，蒙卡燃耗计算遇到重大挑战。原因来自多个方面。其一，如上所述，传统蒙卡程序缺乏精确高效模拟随机介质的能力。其二，传统蒙卡程序缺乏处理大规模燃耗区的能力。以高温气冷堆为例，包覆颗粒弥散型燃料元件或组件内含有的燃料颗粒数量多达104~106量级，整个堆芯的燃料颗粒数以亿计。
国际上在随机介质蒙卡燃耗计算方面已开展的工作可被概括为以下两个层面： 

(1) 弥散型燃料元件/组件燃耗计算。DeHart和Ulses发布了高温堆燃料元件燃耗基准题［61］，包括燃料颗粒无限重复栅格、燃料球栅元和柱状超栅元三种典型模型。一些参与者［47, 62, 63］使用Serpent、KENO和BGCore等程序给出了部分基准题的计算结果，这些计算大多采用重复结构描述燃料颗粒，并假设所有燃料颗粒具有相同的燃耗速率。除了上述基准题以外，Obara和Onoe［46］使用MVP程序计算了同时含燃料颗粒和可燃毒物颗粒的柱状高温堆燃料元件，Brown等［34］计算了轻水堆全陶瓷包覆颗粒燃料组件。

(2) 弥散燃料全堆燃耗计算。Kim和Venneri［64］以及Edward Read等［65］分别使用McCARD和MonteBurns程序计算了棱柱式高温堆全堆燃耗问题。为了降低问题难度，他们采用反应性等价物理转换(reactivityequivalent physical transformation，RPT)方法［66］，将包覆颗粒燃料元件近似等效为均匀材料。RPT方法是一种近似方法，其误差随燃耗增大而增大。燃耗区的划分也比较粗糙，需要人为对燃耗区进行分组。
综上所述，国际上现有的三种随机介质输运计算方法均存在缺点及局限，难以适用各种新型燃料设计的需求(如含毒物的FCM燃料)。而在燃耗计算方面，计算精细程度较低，元件/组件级别的计算只能以燃料球、燃料棒作为整个燃耗单元，无法计算颗粒的燃耗深度； 而全堆计算则采用RPT等方法，燃耗区的划分也比较粗糙。
1.3研究目标与内容
本书的研究内容包括两部分，分别是蒙卡全寿期高保真耦合模拟研究和蒙卡随机介质精细计算研究。这两部分工作都基于清华大学工程物理系REAL实验室自主开发的反应堆蒙卡程序RMC开展，依托RMC已具备的基本功能(如临界计算、燃耗计算、混合并行、数据分解与区域分解、临界搜索和平衡氙修正等)，实现RMC全寿期高保真耦合模拟和随机介质精细计算的相应核心算法及程序功能模块的研发。
蒙卡全寿期高保真耦合模拟研究的主要内容包括： 
(1) 温度相关截面全能区在线处理方法。研究蒙卡程序在可分辨共振能区、热化能区和不可分辨共振能区的在线截面处理方法。提出了可分辨共振能区的基于射线追踪法(ray tracking)的靶核运动抽样法
(target motion sampling)和基于改进高斯厄米特
（GaussHermite）求积组的在线多普勒展宽方法、热化能区热散射数据在线插值方法以及不可分辨共振能区概率表在线插值方法，并验证其效率与精度。

(2) 蒙卡物理热工耦合。在完成温度相关截面在线处理功能开发的基础上，实现了RMC和子通道热工水力程序CTF的通用耦合。并研究了统计偏差造成的不稳定性对蒙卡耦合的影响，采用了松弛因子方法进行功率