2. 中国科学院等离子体物理研究所 合肥 230031
2. Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, China
中国聚变工程实验堆(China Fusion Engineering Test Reactor, CFETR)是正处于概念设计中的超导托卡马克装置，其目标是达到50-200MW的聚变功率，年聚变燃烧时间(Duty Cycle)达到30%-50%，氚增值比(Tritium Breeding Ratio, TBR)大于1，来展示聚变能的可行性并验证堆系统的氚自持[1]。
CFETR由多个部件构成，设计过程中需要考虑复杂的物理和工程间的约束关系。正在研发的CFETR集成设计平台是为了提高聚变堆的设计效率，实现聚变堆的自洽设计[2]。CFETR集成设计平台包括物理设计平台和工程设计平台，工程设计平台采用模块化方式，由磁体、真空室、偏滤器、中子学等模块构成[3-5]。
中子学模块的功能之一是提供中子学计算结果的后处理，为其他各工程设计模块提供必要的中子学参数评估和载荷输入。热工水力分析是聚变堆的堆内部件设计的基础，主要使用CFD (Computational Fluid Dynamics)软件如ANSYS CFX、FLUENT等进行。堆内部件产生的中子核热是热工水力分析的重要热源项，中子核热数据来源于中子学计算。而对于中子学软件MCNP (Monte Carlo N particle transport code)输出记载核热的网格计数文件，存在数据量大、人工分析效率低、CFD软件无法直接读取作为热源项加载等问题，为此需要开发中子学与CFD核热耦合模块，作为中子学模块的一个子模块，以提高CFETR中子学计算后处理的效率。
在中子学与CFD核热传递过程中，目前的聚变堆部件热工设计中一般采用一维或者二维核热数据[6-8]，精确的三维核热耦合有利于聚变堆的工程分析，在本核热耦合模块中采用三维核热耦合方式。
由于中子学与CFD计算网格不匹配，需要通过网格插值将MCNP得到的核热数据映射到CFD网格上以实现三维核热耦合。目前网格插值方法主要有网格-网格插值法和点-点插值法，两种方法在核裂变和核聚变领域均有采用。在核裂变领域，Thomas等[9]基于点-点插值实现中子学软件DeCART和Star-CD之间的核热耦合；何杰等[10]基于点-点插值实现MCNP与Fluent之间的核热耦合；陈军等[11]基于网格-网格插值实现MCNP和STAR-CCM+之间的核热耦合。在核聚变领域，Ying等[12]基于网格-网格插值实现MCNP和SC/Tetra软件之间的核热耦合；Wilson等[13]基于点-点插值实现MCNP和Star-CD软件之间的核热耦合；Qiu等[14]基于网格-网格插值实现MCNP和CFX、FLUENT之间的核热耦合。在耦合模块中，分别提供上述两种插值方法的程序进行核热传递，以方便工程设计和分析。
本文首先介绍了核热耦合模块的流程，为实现核热耦合功能，进行了耦合程序的开发。然后，为了验证核热耦合模块的有效性，使用MCNP与CFX软件，以CFETR的一种氦冷陶瓷包层(Helium Cooled Ceramic Breeder blanket, HCCB)的增殖单元模块为分析对象，采用三维核热耦合方式进行了热工水力分析，也单独进行了常用的一维核热数据传递下的计算用于对比。最后对计算结果中的核热分布和温度场分布进行了分析讨论。
1 核热耦合模块的开发 1.1 模块流程核热耦合模块的流程图如图 1所示。
|
图 1 核热耦合模块流程图 Figure 1 Workflow of neutronics-thermal hydraulics coupling module. |
1) 获取三维核热数据：对于给定的分析模型，经MCNP三维中子学计算后，产生网格计数文件(meshtal file)，该文件记录了重叠网格计数卡(Superimposed Mesh tally, FMesh)划分的每个网格上的核热密度。
2) 获取CFD网格文件：用CFX/Fluent计算时，首先要生成网格，这里采用CGNS (CFD General Notation System)网格格式进行网格数据信息的传递。CGNS是美国波音公司和美国国家航空航天局于1994年联合提出的一项旨在标准化CFD输入与输出的系统，包括网格、区域连接、边界条件、流场解及其它辅助信息的数据结构等[15]。
3) 耦合程序插值：MCNP的FMesh划分的网格是正交结构化网格，而CFD计算网格为非结构化网格，MCNP和CFD软件所采用网格类型及尺寸不一致产生了网格不匹配问题，因此进行核热数据传递时，需要进行两者网格之间的插值映射。
4) 核热输出至CFD软件：耦合程序插值计算出CFD需要的每个网格上的核热密度后，以C语言文件格式的用户自定义函数(User-Defined Function, UDF)形式输出给Fluent，Fluent计算时编译该文件即可；以Fortran语言文件格式的用户函数形式输出给CFX，CFX计算时进行编译生成.dll格式的动态库，在前处理设置中调用该动态库即可。
5) 进行热工水力分析，给出热工水力分析结果。
1.2 耦合程序 1.2.1 耦合原理在CFD计算中，每种材料的计算域都需要单独地设置热源，热源只能是这一种材料的产热。MCNP能按两种方式对核热结果进行统计：一种是按实际材料统计的混合材料方式；另一种是假设所有材料均为一种材料统计的单一材料方式[12]。混合材料方式统计的每个网格体热值是这个网格内的均值，当网格包含多种材料时无法分离出每种材料的产热输出给CFD，所以需采用单一材料方式统计。即假设所有材料均为一种材料，单独地统计每种材料每种粒子（如中子、光子等）的产热，将每种材料的所有粒子产热累加得到该材料产热，最后映射到该材料的CFD计算网格上。
耦合模块中采用的网格-网格插值法即体积权重法，是根据MCNP网格和CFD网格的重叠体积为权重来计算CFD网格的核热密度值；耦合模块中采用的点-点插值法是根据CFD网格几何中心点坐标确定其所在的MCNP网格，将该MCNP网格的核热密度值赋予CFD网格。通常，CFD的网格要比MCNP网格精细，一个MCNP网格包含多个CFD网格。当CFD网格被MCNP网格包含时，不会产生核热传递误差。而当CFD网格横跨多个MCNP网格时，由于MCNP相邻网格中子通量密度一般不会出现突变，相应地核热密度也不会出现突变，忽略横跨的其他MCNP网格一般不会造成大的核热传递误差。此外，在Fluent中区域离散化采用单元中心法(Cell Centered Scheme)[16]，而在CFX中区域离散化采用单元顶点法(Cell Vertex Scheme)[16]，故在MCNP与CFX耦合时，求出单元中心点核热密度值后，采取最临近点方法求出单元顶点的核热密度值。
1.2.2 耦合程序实现网格-网格插值法是通过基于开源平台SALOME开发的McMeshTran[14]程序实现。
点-点插值法是通过开发基于文献[17]预留的CFD接口程序实现。程序以C/C++作为开发语言，在MCNP计数文件的数据逻辑操作中加入坐标变换功能，以解决中子学模型和热工模型的坐标系不一致问题。该程序命名为MT2X (Mesh tally to X)，X可以是Ensight/Paraview/ANSYS/CFX。
2 核热耦合模块功能验证为了验证耦合模块功能的可靠性，选取CFETR的一种S型弯氦冷陶瓷包层的增殖单元模块为测试模型。CFD软件采用CFX，应用MT2X、McMeshTran进行三维核热耦合计算，同时也进行了一维核热耦合方式下的计算。对三种情况下的耦合计算结果进行相互对比分析，重点分析了核热的分布和温度场的分布。一维核热耦合(1D)方式是通过使用MCNP中的分段计数卡(Tally Segment card, FS)，得出各个材料沿单一维度的中子核热分布，最终加载到CFX软件中进行后续热工计算。
2.1 验证模型目前S型弯氦冷陶瓷包层还在设计优化中，其主要优点是减少了氦气联箱数目，为增殖区留出了更大的空间，结构示意图如图 2所示。考虑到模型的几何对称性，最终采用包含两个S型氦冷流道的切片模型作为CFX计算模型（图 3），进行分析以简便计算，整个模型尺度为630mm×1002mm× 11mm (x×y×z)，切片模型的氦冷流道采取了相应的倒角措施。该包层管道结构材料采用低活化铁素体/马氏体(Reduced Activation Ferritic/Martensitic, RAFM)钢，氚增殖剂采用正硅酸锂(Li2SiO4)，中子倍增剂采用铍(Be)。各种材料的物性参数取自CFETR集成平台的材料数据库[3]。
|
图 2 S型弯氦冷陶瓷包层模型 Figure 2 S-shape HCCB blanket model. |
|
图 3 CFX计算模型 Figure 3 Simulation model of CFX. |
中子学计算采用的网格尺度在x、y方向最小为1.5mm，z方向为2.5mm，并且根据文献[11]的热量传递误差分析，在各材料边界处尽可能地进行了边界区分以避免大的核热传递误差。CFX计算网格单元数分别为：RAFM区为1.2×105，Li2SiO4区为1.7×105，Be区为1.5×105，流体He区为3.3×105。
CFX计算模型的边界条件为：氦气采取质量流率进口和压力出口的组合，进口质量流率为0.004kg·s-1，进口温度为573K，出口压力为8MPa。靠近第一壁边界取恒温668K（由于S型弯包层第一壁氦气流向与切片方向垂直，切片处对应的第一壁各条管道位置相似，因此近似假设成一个恒温边界条件），其他边界均为绝热。参考文献[6]，CFX计算采用k-ε模型和标准壁面函数。
2.2 结果分析与讨论 2.2.1 总核热对比从MCNP传递到CFD的核热理想的结果是传递前后每种材料的总核热值均保持守恒，但实际核热传递中会出现插值误差，为此需要对传递后的核热进行统计来分析误差。统计三种耦合方式传递后各材料区域的总热量，将MCNP使用栅元计数法统计得到的总热量作为误差分析的基准值，得出的数据统计表如表 1、2所示。
|
表 1 各材料的总核热 Table 1 Total power deposition of each material. |
|
|
表 2 各材料的总核热误差 Table 2 Relative errors of total power deposition of each material. |
注意到表 1中各个材料区都出现了总核热数值1D < McMeshTran < MT2X的情况，需要特别说明的是，这是插值过程中由于插值误差造成的偶然情况。从表 2中可以看出，1D、McMeshTran、MT2X三种情况下的各个材料核热误差最大值均小于0.5%。而MCNP统计误差在5%以下认为结果可信，这说明三种核热耦合方式均能使传递后的总核热误差控制在允许范围内，误差值说明McMeshTran和MT2X实现了高精度的核热耦合。
2.2.2 核热空间分布对比进一步地，对比三者核热密度空间分布情况。由于MCNP和CFX网格不同，虽然用MCNP统计出的网格核热密度是混合材料的核热值，数据存储于网格中心点，而插值后CFX中显示的每个网格核热值是单一材料的，数据存储于网格中心，但是仍能作直观对比。核热密度分布情况如图 4，其中图 4(a)是利用Paraview进行了核热可视化的结果。从图 4中可以看到，两种三维核热耦合方式相比一维更准确地体现了MCNP计算结果的分布情况。
|
图 4 核热密度分布 (a) MCNP混合材料，(b) 1D，(c) McMeshTran，(d) MT2X Figure 4 Nuclear heating distributions. (a) MCNP material-mixed results, (b)1D, (c) McMeshTran, (d) MT2X |
为更详细地对比核热局部分布情况，选取核热密度变化较为明显的z=0m平面和y=0.054m平面的交界线L（如图 3中的虚线部分所示），得到直线L沿x轴的核热密度分布。考虑对称性，取L/2作出核热密度分布图（图 5）。从图 5可以看出，1D、McMeshTran、MT2X和MCNP混合材料的核热密度在RAFM区和Be区基本一致，而在Li2SiO4区差别较大。在Li2SiO4的边界处，McMeshTran、MT2X和MCNP混合材料的峰值不同。这个现象一方面是因为MCNP混合材料在边界处统计的核热密度是材料混合后的统计值，而非单一材料的实际值，另一方面是因为不同插值方法会产生误差。例如，考虑到Li2SiO4区产生的核热高于其他材料区，混合材料方式统计的硅酸锂边界处的核热密度值低于单一Li2SiO4材料统计方式边界处的核热密度值。McMeshTran方法所得的核热密度值是通过单一材料统计方式得到，故而会出现图 5中在Li2SiO4边界处McMeshTran比MCNP混合材料统计结果偏高的情况。
|
图 5 L/2的核热密度分布 Figure 5 Nuclear heating distribution on line L/2. |
1D与McMeshTran、MT2X相比，在氦气进出口处（x=-0.45m）和Li2SiO4边界处差别明显，相对误差最大值为24.4%（x=-0.452m处）。在氦气进出口处相差较大是因为MCNP计算进出口处的实际核热分布较中间低，1D将核热密度平均后造成核热密度在Li2SiO4区平整。
2.2.3 温度场对比核热耦合计算的最终目的是获取各部件的温度场分布，以便进行下一步的结构热应力分析。三种核热耦合方法得到的温度分布如图 6所示，同样可以提取L/2处的局部温度分布情况（图 7）。从图 6、7可以看出，McMeshTran和MT2X计算出的温度分布大致相同。1D与McMeshTran、MT2X相比，在RAFM区和Be区差别不大，而在Li2SiO4区差异较明显，尤其体现在靠近氦气进出口处（x=-0.45m），最高温差可达30K，相对误差最大值为3%（x=-0.45m处）。从图 5的核热密度分布情况可以看出，进出口处温度差异较大是该处1D核热密度高于McMeshTran、MT2X造成的。
|
图 6 1D (a)、McMeshTran (b)、MT2X (c)的温度分布 Figure 6 Temperature distributions of 1D (a), McMeshTran (b) and MT2X (c). |
|
图 7 L/2温度分布曲线 Figure 7 Temperature distribution on line L/2. |
各材料计算域的最高温度如表 3所示，Be、Li2SiO4、RAFM的允许最高温度分别为923 K、1193K和823 K[3]。从表 3可以看出各个材料温度均在允许范围内，三种方式在最高温度上体现出了较好的一致性，整体相差较小。McMeshTran和MT2X最大差异在1K左右，而1D与McMeshTran和MT2X最大差异在3K左右。
|
|
表 3 各固体域最高温度 Table 3 Highest temperature of solid domains. |
综合核热和温度结果的对比，可以看出三维耦程序实现了高精度的核热传递，和一维核热传递下的结果对比显示出了三维核热耦合在局部的精确性，验证了核热耦合模块功能的可靠性。
3 结语本文基于CFETR集成设计平台开发了核热耦合模块，采用网格-网格插值法和点-点插值法，通过MT2X和McMeshTran程序，提供两种途径解决了MCNP核热输出文件不能直接导入CFD软件进行分析等问题。核热耦合模块实现了高精度的三维核热耦合，有利于提高聚变堆设计效率。以CFETR的一种S型弯氦冷陶瓷包层的增殖单元为对象进行了热工水力分析，分析结果表明了核热耦合模块功能的可靠性。
| [1] |
Wan Y X. Design and strategy for CFETR[C]. Proceedings of 25th Symposium on Fusion Engineering (SOFE), San Francisco, CA, USA, 2013.
|
| [2] |
Ye M Y, Wang S J, Wang Z W, et al. Development of system code of CFETR design[C]. Proceedings of 26th Symposium on Fusion Engineering (SOFE), Austin, TX, USA, 2015.
|
| [3] |
Wang S J, Ye M Y, Zhu C, et al. The Integration platform development of system code for CFETR[J]. IEEE Transactions on Plasma Science, 2016, 44(9): 1745-1750. DOI:10.1109/TPS.2016.2598599 |
| [4] |
Zhang J W, Ye M Y, Peng X B, et al. Preliminary design and verification of divertor module of CFETR system code[J]. IEEE Transactions on Plasma Science, 2016, 44(9): 1763-1767. DOI:10.1109/TPS.2016.2599264 |
| [5] |
Zhu C, Ye M Y, Liu X F, et al. Development of vacuum vessel design and analysis module for CFETR integration design platform[J]. Science and Technology of Nuclear Installations, 2016, 5321057. DOI:10.1155/2016/5321057 |
| [6] |
Zhou G M, Li M, Liu Q W, et al. Thermal analysis of breeder unit for helium cooled solid breeder blanket of Chinese fusion engineering test reactor[J]. Journal of Fusion Energy, 2015, 34(2): 339-345. DOI:10.1007/s10894-014-9798-y |
| [7] |
Zanino R, Richard L S, Subba F, et al. CFD analysis of a regular sector of the ITER vacuum vessel Part Ⅱ: Thermal-hydraulic effects of the nuclear heat load[J]. Fusion Engineering and Design, 2013, 88(12): 3248-3262. DOI:10.1016/j.fusengdes.2013.10.005 |
| [8] |
Chen X M, Ma X B, Jiang K C, et al. Thermal hydraulic design and analysis of a water-cooled ceramic breeder blanket with superheated steam for CFETR[J]. Plasma Science and Technology, 2015, 17(9): 787-791. DOI:10.1088/1009-0630/17/9/11 |
| [9] |
Thomas J W, Zhong Z, Sofu T, et al. Methodology for coupling computational fluid dynamics and integral transport neutronics[C]. Proceedings of PHYSOR 2004-The Physics of Fuel Cycles and Advanced Nuclear Systems: Global Developments, Chicago, IL, USA, 2004.
|
| [10] |
何杰, 夏晓彬, 蔡军, 等. 2 MW液态钍基熔盐实验堆主屏蔽温度场分析[J]. 核技术, 2016, 39(4): 040601. HE Jie, XIA Xiaobin, CAI Jun, et al. Temperature field analysis for the main shielding of the 2-MW thorium-based molten salt experimental reactor[J]. Nuclear Techniques, 2016, 39(4): 040601. DOI:10.11889/j.0253-3219.2016.hjs.39.040601 |
| [11] |
陈军, 曹良志, 郑友琦, 等. 基于蒙特卡罗方法和CFD方法的物理-热工耦合计算[J]. 原子能科学技术, 2016, 50(2): 301-305. CHEN Jun, CAO Liangzhi, ZHENG Youqi, et al. Neutronics and thermal-hydraulics coupling calculation based on Monte Carlo and CFD method[J]. Atomic Energy Science and Technology, 2016, 50(2): 301-305. DOI:10.7538/yzk.2016.50.02.0301 |
| [12] |
Ying A, Narula M, Abdou M, et al. Toward an integrated simulation predictive capability for fusion plasma chamber systems[J]. Fusion Science and Technology, 2009, 56(2): 918-924. DOI:10.13182/FST09-A9027 |
| [13] |
Wilson P, Tautges T, Kraftcheck J, et al. Acceleration techniques for the direct use of CAD-based geometry in fusion neutronics analysis[J]. Fusion Engineering and Design, 2010, 85(10): 1759-1765. DOI:10.1016/j.fusengdes.2010.05.030 |
| [14] |
Qiu Y F, Lu P, Ulrich F, et al. A generic data translation scheme for the coupling of high-fidelity fusion neutronics and CFD calculations[J]. Fusion Engineering and Design, 2014, 89(7): 1330-1335. DOI:10.1016/j.fusengdes.2014.02.044 |
| [15] |
Legensky S, Edwards D, Bush R, et al. CFD general notation system (CGNS)-status and future directions[C]. 40th AIAA Aerospace Sciences Meeting & Exhibit, 2002: 752.
|
| [16] |
陶文铨. 数值传热学[M]. 西安: 西安交通大学出版社, 2001. TAO Wenquan. Numerical heat transfer[M]. Xi'an: Xi'an Jiaotong University Press, 2001. |
| [17] |
何泽鸿, 叶民友, 王忠伟, 等. MCNP与后处理软件耦合应用程序开发[J]. 核技术, 2015, 38(6): 060601. HE Zehong, YE Minyou, WANG Zhongwei, et al. MCNP coupled with post-processing software application development[J]. Nuclear Techniques, 2015, 38(6): 060601. DOI:10.11889/j.0253-3219.2015.hjs.38.060601 |