cosRMC复杂曲面开发及在聚变中子学的应用

刘仕倡，申靖文，严伊蔓，杜 华，孙业帅，秦 瑶，郑 俞， 卢 棚，李 夏，张 显，刘松林，陈义学

(1.华北电力大学 核科学与工程学院，北京 102206； 2.国家电投集团科学技术研究院有限公司 核电软件技术中心，北京 102209； 3.中国科学院 等离子体物理研究所，安徽 合肥 230031)

中国聚变工程实验堆(CFETR)[1]是我国于2011年开始设计研究的超导托卡马克反应堆，在DEMO级别聚变电站与ITER间起到了过渡衔接的作用。托卡马克型的聚变堆堆芯主要包括：第一壁、包层、偏滤器、真空室、辐射屏蔽、生物屏蔽等。包层是聚变堆堆芯的一个重要组成部分，在强中子辐照下，包层主要有氚增殖、能量转换及辐射屏蔽等作用，因此，包层中子学的设计与聚变堆中的大量问题均相关。此外，托卡马克装置具有非常复杂的几何结构，包含各种结合体并且结构十分精细，而蒙特卡罗方法几何描述精细的特点无疑非常适宜于解决这样的问题[2]。但如果用于模拟计算的蒙特卡罗程序不具备高次曲面的描述能力，仅使用一阶或二阶曲面并不能满足几何建模精细化的需求。那么对于这类复杂几何模型的建立就仅能采取几何近似的方式，如可将托卡马克装置的环形真空室的圆环面近似为很多段直线，不仅耗时耗力，还会导致精度损失，不利于聚变堆的准确模拟。国际上具有处理聚变堆复杂几何能力的蒙特卡罗程序有美国Los Alamos实验室的MCNP[3]、日本原子能机构的MVP[4]、法国原子能委员会的TRIPOLI[5]等。国内方面，中国科学院核能安全技术研究所的SuperMC[6]在磁约束核聚变装置的蒙特卡罗模拟方面进行了大量研究，北京应用物理与计算数学研究所的JMCT[7]也进行了一些研究。其中，MCNP是应用最为广泛的经典蒙特卡罗程序，其复杂曲面处理方法被广泛应用于国内外蒙特卡罗程序，然而对该方法的数学模型和计算原理进行系统推导和阐述的研究报道较少。

本文基于自主化蒙特卡罗程序cosRMC，研究及阐述蒙特卡罗复杂曲面建模的数学模型和计算方法，并对原有的几何模块进行拓展，使其具备对复杂几何模型的描述能力，最后将其应用于CFETR三维模型的聚变中子学分析中，对聚变堆包层设计时的3个关键的中子学参数进行计算分析。

1 cosRMC复杂几何建模功能

1.1 cosRMC及其几何建模功能介绍

cosRMC是由清华大学工程物理系反应堆工程计算分析实验室(REAL团队)与国家电投集团科学技术研究院有限公司核电软件技术中心联合研发的以蒙特卡罗粒子输运方法为理论基础的自主化蒙特卡罗程序[8]。主要用于核工程领域中固定源计算、临界计算、燃耗计算、瞬态计算和均匀化并群计算等。

cosRMC采用文本建模方式，基于构造几何实体(CSG)对系统几何结构进行描述。即通过对交界面不同指向的区域进行布尔运算，如交、并和非等运算，构成对体的描述。目前，cosRMC程序可进行描述的面包括一般平面、任意球面、轴心位于坐标轴上或平行于坐标轴的圆柱面或圆锥面等标准形式的曲面，所有曲面的曲面方程均使用解析方程表示方式。在本文工作之前，cosRMC无法实现对一些复杂曲面的描述，如任意位置的圆柱面、圆锥面、椭圆面和圆环面等。

1.2 复杂几何建模功能实现

粒子在蒙特卡罗程序几何中的输运需确定粒子沿飞行方向与曲面的交点，以圆环面为例，通过解析方法简化为求解一个四次多项式方程解析解的问题。

1) 中子轨迹与复杂曲面的联立求解

(1)

进一步在粒子输运问题中，在射线簇{x=ξ+aX,y=η+BX,z=ζ+γX;-∞

{[(1-β2)+ρβ2]X2+[2(aξ+γζ)+2ρβη-

2ρβy0]+(ξ2+ζ2+ρη2-2ρηy0+B0)}2=

A0[(1-β2)X2+2(aξ+γζ)X+ξ2+ζ2]

(2)

式(2)可转化为：

(GX2+MX+U)2=A0(FX2+LX+T)

(3)

其中，F=1-β2，G=F+ρβ2，L=2(aξ+γζ)，M=L+2ρβ(η-y0)，T=ξ2+ζ2，U=T+ρη(η-2y0)+B0。

所以，射线簇{x=ξ+aX,y=η+BX,z=ζ+γX;X>0}与圆环面求交点(x,y,z)的问题可转化成一个求解式(4)所示四次方程的正根X的问题，仅需判断射线簇是与圆环面的内表面还是外表面存在交点即可。

X4+BX3+CX2+DX+E=0

(4)

2) 复杂曲面多项式方程的求解

由于舍入误差的存在，使得利用计算机程序求解形如式(4)的四次多项式方程均非常困难。Maisonnier和Everett提出了一种求解“射线簇与三维或四维曲面与交点”的单精度四次算法[9]，该算法可分为两步，第1步是将粒子穿过圆环面时四次多项式方程的求解转化为三次方程的求解。

首先讨论三阶曲面的解析方程，可利用泰勒展开求解：

f(x)=d+cx+bx2+x3=

(5)

对式(5)求导代换后可得到新的三次解析方程：

g(x)=y3+py+q

(6)

Δ3=-4(27)W

(7)

同样对四阶曲面的解析方程也进行泰勒展开：

F(X)=E+DX+CX2+BX3+X4=

(8)

通过推导可得出，四阶曲面方程的判别式为Δ4=Δ3。根据判别式条件就变成了一个类似于三阶曲面方程求解的过程。

第2步是根据多项式系数判定来确定是否存在任意正实根，主要包括3种情况：当Δ3>0(W<0)时，有3个不同的实根；当Δ3=0(W=0)时，有3个含有重复根的实根；当Δ3<0(W>0)时，有单个实根和1对共轭复根。

根据以上理论方法，本文在cosRMC原有的几何模块中新添加了椭圆面、双曲面、抛物面、圆环面等复杂曲面，依旧采用助记名的方式存储，如表1所列。

1.3 功能测试

聚变堆概念设计堆PPCS(power plant conceptual study)是2000年于欧洲启动的最为系统的聚变堆研究计划，目的是帮助评估聚变能形式和确立欧洲聚变计划的可持续性和优先权，属于世界上目前较具规模和代表性的磁约束聚变动力反应堆概念设计系列之一，如图1所示。作为典型的磁约束聚变堆概念设计系列之一，PPCS具有非常复杂的几何机构，在使用蒙特卡罗方法对其进行聚变中子学分析时，需构建诸多包含复杂曲面的几何体，以及包层模块重复结构的描述，如图1c中标注所示。

表1 新增cosRMC复杂曲面卡Table 1 New cosRMC complex surface card

a——PPCS的cosRMC模型；b——PPCS的MCNP模型；c——PPCS模型中选取的栅元 图1 PPCS的cosRMC及MCNP的模型示意图Fig.1 Model diagram of PPCS with cosRMC and MCNP

利用cosRMC程序对PPCS进行建模(图1)后进行模拟计算，中子源采用14 MeV的各向同性点源。MCNP和cosRMC模拟的总粒子数均为1×108。选取了18个有高次曲面的栅元，图1c中标注的数字即为选取栅元的栅元号。统计其体通量，与MCNP的计算结果进行比较，令相对统计误差为：

(9)

σMCNP和σcosRMC为MCNP和cosRMC的计数器相对统计误差。栅元通量的对比如图2所示。可看出，cosRMC程序对所有几何体上包含高次曲面的体通量的计算结果，与MCNP的计算结果吻合良好，相对偏差几乎均小于10%。且cosRMC单核计算时间为250.980 6 min，MCNP单核计算时间为258.74 min，计算效率相当。说明了cosRMC程序中实现的复杂几何建模功能的正确性。

图2 PPCS模型cosRMC与MCNP模拟计算结果比较Fig.2 Comparison of simulation results of PPCS between cosRMC and MCNP

2 CFETR建模与聚变中子学分析

2.1 CFETR建模

CFETR的三维模型体积大、几何结构复杂，托卡马克装置中包含环形真空室、环向线圈、极向线圈等结构部件使得在使用蒙特卡罗方法对其进行聚变中子学分析时，需构建诸多包含四阶曲面的几何体。使用拓展了复杂几何建模功能的cosRMC程序对CFETR三维模型进行精细化建模，如图3所示，与MCNP程序构建的模型图进行对比，如图4所示。从模型图可看出cosRMC程序描述的CFETR三维模型的正确性。

图3 CFETR的cosRMC模型径向与纵向剖面图Fig.3 Radial and longitudinal profiles of cosRMC model of CFETR

图4 CFETR的MCNP模型纵向剖面图Fig.4 Longitudinal profile of MCNP model of CFETR

进一步对CFETR的包层中子学中的3个重要参数进行计算分析。CFETR的包层模块如图5所示，沿环向上和极向上进行分割，分为内包层模块和外包层模块。中子源采用氘氚聚变中子源，分布于5个栅元中，中子源强分别为1%、6%、13%、30%和50%。源中子能量为14 MeV的高斯聚变谱，方向各向同性。CFETR的氚增殖比(TBR)、中子壁载荷(NWL)和核热沉积计算过程中，聚变功率采用第1阶段的运行功率200 MW。MCNP和cosRMC模拟的总粒子数均为6.9×106。MCNP和cosRMC均采用FENDL-2.1核截面数据库进行计算。

图5 CFETR包层结构编号及位置示意图Fig.5 Diagram of CFETR blanket structure number and location

2.2 TBR

氘氚聚变反应堆能否投入商业运行的首要条件是能实现氚的自持。聚变中子学分析计算时，增殖材料中6Li、7Li的产氚率为：

T6=NLi∬σ6Li(n,α)(r,E)φ(r,E)dEdr

(10)

T7=NLi∬σ7Li(n,n′α)(r,E)φ(r,E)dEdr

(11)

其中，σ6Li(n,α)和σ7Li(n,n′α)分别为6Li的(n,α)反应截面和7Li的(n,n′α)反应截面。那么，堆芯等离子体平均1次氘氚聚变反应在包层内增殖的TBR为：

(12)

其中，Sn为聚变中子源强。利用cosRMC对其增殖包层各模块的氚增殖比进行计算，与MCNP的计算结果进行比较(表2)。可看出，两程序对各增殖包层氚增殖比的计算结果吻合良好，相对偏差均小于1%。并且，通过模拟计算可发现，增殖包层总的TBR为1.2，能满足氘氚聚变反应堆运行时的氚自持需求。

2.3 NWL

NWL即加载到第一壁表面(W-PFC)单位面积上的平均聚变中子功率。之所以要计算NWL，是因在托卡马克装置中第一壁是面向等离子体的部件，在包层中起着重要作用，它的安全直接影响包层模块的稳定性，其计算方法如式(13)所示，其中A为表面积。

(13)

表3列出利用cosRMC计算的CFETR第一壁表面的NWL分布，各包层模块的结果与MCNP符合良好，相对偏差均在1%以内。进一步计算可得到第一壁表面的平均NWL为0.305 1 MW/m2，NWL峰值出现在位于赤道平面的BLK3包层。

表2 cosRMC和MCNP计算的CFETR增殖包层各模块TBR的比较Table 2 Comparison of TBR in each module of CFETR blanket

表3 CFETR第一壁表面的NWL分布计算结果比较Table 3 Comparison of NWL in each module of CFETR first wall surface

2.4 核热沉积

托卡马克装置中用于约束等离子体的磁铁结构通常是由带有不锈钢包层的铜制成的，以提供抵抗变形的强度。铜线圈具有极好的电导率，允许通过的电流强度高达数十MA。但同时会在铜线圈上产生较大的热沉积。限制等离子体放电的因素之一即是当高强度的电流通过磁铁后，磁铁温度的升高，应使核热沉积最小化。因此，研究磁铁线圈系统上产生的核热沉积是十分必要的，其计算方法如式(14)所示。

(14)

其中：H(E)为热数；σt(E)为总截面；ρ为核数密度；m为栅元质量。

表4列出利用cosRMC程序计算的CFETR的BLK1包层内各部件的核热沉积，两程序模拟计算结果吻合良好，相对偏差均小于5%。同一包层上不同部件的核热沉积不同是因为各部件所处位置不同，中子通量密度不同，受到的中子辐照强度也不同。

3 结论

本文研究及推导了蒙特卡罗复杂曲面建模的数学模型和计算方法，并在自主化蒙特卡罗程序cosRMC原有的几何模块上拓展了复杂曲面建模功能。通过构建聚变堆概念设计PPCS模型并进行计算，栅元通量与MCNP的相对偏差小于10%，说明了该功能的正确性。

表4 CFETR的BLK1包层各部件核热沉积计算结果比较Table 4 Comparison of nuclear thermal deposition in different components of BLK1 blanket

进一步利用cosRMC对CFETR进行三维精细化建模，并对包层模块设计中的3个重要中子学参数进行计算分析。结果表明，cosRMC与MCNP计算的各增殖包层TBR的相对偏差均小于1%。且增殖包层总的TBR为1.2，能满足氘氚聚变堆运行时的氚自持需求。cosRMC与MCNP计算的NWL相对偏差均在1%以内，NWL最大值出现在位于赤道平面的BLK3包层。cosRMC与MCNP计算的BLK1包层内各部件的核热沉积的相对偏差均小于5%。研究结果证明了cosRMC应用于聚变堆包层中子学分析的正确性和有效性。CFETR中子学参数的计算分析，也为其设计和优化提供了参考。

感谢清华大学工程物理系核能科学与工程管理研究所反应堆工程计算分析实验室(REAL团队)在程序使用过程中提供的技术支持和帮助。