熔盐堆是第四代国际核能论坛(The Generation IV International Forum,GIF)提出的6种第四代反应堆之一，具有固有安全性、良好的中子经济性、放射性废物少等优势，起源于美国橡树林国家实验室(Oak Ridge National Laboratory,ORNL)^{[1, 2, 3]}。21世纪以来，俄罗斯、法国和日本的研究人员分别提出了熔盐堆MOSART(Molten Salt Actinide Recycler &Transmuter)、MSFR(Molten Salt Fast Reactor)以及小型熔盐堆FUJI系列的概念设计^{[4, 5, 6]}。中国科学院于2011年启动了“未来先进核裂变能——钍基熔盐堆核能系统”战略先导专项，致力于钍基熔盐堆(Thorium Molten Salt Reactor,TMSR)相关技术的研究，TMSR-LF是其发展目标之一^[7]。

反应堆主屏蔽是核反应堆的重要组成部分，用来有效降低反应堆运行时屏蔽体外的辐射剂量水平，以满足辐射安全的要求及反应堆部件材料对辐射限制的要求。温度是影响反应堆主屏蔽性能的重要因素之一。我国NB/T 20194-2012《压水堆核电厂辐射屏蔽设计准则》对普通硅酸盐混凝土屏蔽体温度方面作出以下要求：1)沿混凝土屏蔽体厚度方向最大温差小于100℃·m^-1；2)边界最高温度70℃^[8]。与压水堆冷却剂平均温度（秦山二期压水堆冷却剂平均温度310℃）相比，熔盐堆TMSR-LF1一回路熔盐平均温度高达610℃。熔盐堆的高温特点给屏蔽体的设计分析带来了挑战，其中特别需要对熔盐堆主屏蔽的温度场进行系统的计算分析。

针对不同反应堆的屏蔽体温度场分析已有较多研究。熔盐实验堆(Molten Salt Reactor Experiment,MSRE)的屏蔽设计中，参考反应堆屏蔽设计公式以及采用一维确定论粒子输运计算程序DTC设计了MSRE的主屏蔽，并通过热屏蔽减少辐射剂量以及向屏蔽体的传热^{[2, 9]}。快堆CEFR堆顶固定屏蔽冷却系统分析中采用了计算流体力学方法(Computational Fluid Dynamics,CFD)对系统进行三维数值模拟，并提出优化建议^[10]。针对5MW低功率堆主屏蔽设计时，假设主屏蔽内边界恒温条件下，采用程序STFO和ANISN计算其温度场^[11]。在压水堆堆内构件释热率计算方面开展了蒙特卡罗程序MCNP(Monte Carlo N Particle Transport Code)和三维离散纵标方法程序TORT的适用性研究^[12]。高温气冷堆生物屏蔽设计中，假定生物屏蔽层内边界温度恒定条件下，运用Abaqus 6.7商业有限元软件分析额定工况和事故工况下生物屏蔽层的安全性^[13]。以上研究对反应堆温度场计算模型进行了合理简化，使用了不同方法和软件成功模拟分析了反应堆主屏蔽温度分布或反应堆构件释热率分布。近年来，研究者采用计算流体力学方法对熔盐堆做了重要研究。张大林等^[14]在考虑熔盐流动对缓发中子份额的影响下（包括CFD方法计算缓发中子份额），提出了适用于液态反应堆的中子动力学模型。王成龙等^[15]采用有限元法模拟了熔盐堆事故工况下NaK热管的瞬态行为。针对氟盐冷却高温试验堆，王成龙等^[16]采用子通道一维模型获得可能出现最高温度的区域，然后采用ANSYS CFX进行三维模拟，结果显示两者符合良好。若联合计算流体力学软件和中子物理输运计算软件进行温度场计算，可以用实际计算的内边界温度分布代替假设温度，将更加逼近真实情况。本文基于现阶段2MWth液态熔盐堆TMSR-LF1的概念设计，采用MCNP和ANSYS Fluent两种程序，利用Python语言编写了MCNP后处理程序MTF(MCNP to Fluent)，联合这两个程序对熔盐堆TMSR-LF1的主屏蔽体开展热工计算分析。

熔盐堆TMSR-LF1主要由堆芯、保温层、隔热层、主屏蔽等部分组成，其中主屏蔽由普通混凝土墙构成，如图 1。堆芯燃料采用LiF-BeF₂-ThF₄-UF₄ (68-28-0.1-3.9mol%)熔盐，二回路冷却剂成分为LiF-NaF-KF(46.5-11.5-42mol%)。堆芯入口温度600℃，设计要求进出口温差为20℃。熔盐堆TMSR-LF1初始设计主要参数为：堆芯直径1.1m，热功率2MW，堆芯轴向高度1.1m，入口温度600℃，保温层0.5m，熔盐质量流量55.3kg·s^-1，隔热层0.8m，石墨密度1860kg·m^-3，主屏蔽0.7m，熔盐密度2570kg·m^-3。

图 1(Fig.1)

图 1 TMSR-LF1结构示意图 Fig.1 Structure diagram of TMSR-LF1.

堆芯物理计算方法分为确定论方法和非确定论方法，其中非确定论方法又称为蒙特卡罗方法，在描述复杂几何形状方面具有一定的优势，并且能得到足够精度的结果。MCNP是一款基于蒙特卡罗方法的粒子输运计算程序，能够进行三维中子、光子等粒子的输运计算，广泛运用于反应堆临界计算、辐射屏蔽设计等领域，计算精度得到国际的认可。它能够为热工分析提供三维内热源分布。模型是在假设控制棒完全抽出、样品通道空闲以及未设置中子毒物情况下建立^[17]。截面数据库是在ENDF/B-VII数据库的基础上，采用了针对熔盐高温特点加工的数据库ACEDATA^[18]。该数据库选取27-1927℃的30个温度点，采用NJOY加工而成。

由于TMSR-LF1几何结构基本符合1/8对称，为提高效率，堆芯物理模型采用1/8建模。建立的1/8堆芯物理模型见图 2。堆芯熔盐轴向上划分10份，径向熔盐按照熔盐通道几何形状划分为96份，总共973个栅元，采用F7计数卡对973个栅元计数供后续功率密度计算。对称轴处的半径边界设置为反射边界，反应堆容器外为真空边界条件。反应堆处于稳态工况，所以采用临界源kcode计算。裂变源放置坐标通过ksrc卡设置，迭代200次，每次迭代中子数为20000，跳过前50次迭代。计算的标准误差为0.00050，统计误差小于0.1，计算结果可以接受。

图 2(Fig.2)

图 2 1/8 TMSR-LF1堆芯纵截面(左)和横截面(右) Fig.2 Axial section (left) and cross section (right) of the 1/8 TMSR-LF1core.

在反应堆热工设计研究中，CFD方法逐渐得到应用。Fluent是一款流体力学计算的主流CFD商业软件，在反应堆热工水力计算中具有重要的价值。另外，Fluent为用户提供了二次开发的用户接口，即用户自定义函数UDF(User-Defined Function)。当Fluent标准模块不能满足用户要求时，可以通过C语言编写UDF实现预期的要求。

在Fluent计算前采用软件Gambit对反应堆进行建模划分网格，网格采用六面体网格和四面体网格混合，网格质量0.5，高于最低要求0.1，符合要求。热工模型包括堆芯、保温层、隔热层、主屏蔽，模型如图 3。1/8反应堆热工模型与物理模型在空间上完全重合。

图 3(Fig.3)

图 3 熔盐堆TMSR-LF网格模型 Fig.3 Mesh model of TMSR-LF1.

1)主屏蔽外壁面边界条件。主屏蔽和环境的换热方式主要考虑墙体与空气的自然对流传热。边界条件：混凝土外壁与环境的换热系数。空气与主屏蔽的外表面传热系数取5 W·m^-2·℃ ^-1；环境温度的设置。由于环境温度随季节变化，因此计算环境温度为5℃、18℃、25℃、30℃、35℃、40℃不同工况下的混凝土温度场。

2)内热源加载。反应堆堆芯内热源通过将MCNP计算的沉积能量转换为对应的功率密度，再通过程序MTF导入到缺失功率密度大小、空间坐标信息的UDF中，内热源根据不同坐标位置加载对应的功率密度。其中，加载功率密度的UDF采用宏DEFINE_SOURCE实现，宏DEFINE_SOURCE可以为包括能量方程在内的不同输运方程提供源项加载。旁流熔盐体积相对堆芯少，功率密度小，所以整个旁流熔盐采用均匀功率密度0.164MW·m^-3。

3)其他边界条件。对称面边界条件设置为对称边界，进口边界条件设置为速度进口，出口边界条件设置为压力出口。考虑重力影响，设置重力加速度9.8m·s^-2。采用SIMPLE算法。

一般地，在反应堆温度场计算中采用的主流软件有Fluent、CFX等，但需要与堆芯中子物理进行联合计算。MCNP可以提供热工软件计算所需的功率密度，计算结果向热工软件的导入过程一般是人工完成。因为MCNP计算结果以文本形式呈现，其结果包含信息量大、数据不够直观，所以使用传统的数据分析方法（人工提取、分析数据）容易出错，而且效率低下。同时MCNP计算结果不能直接导入到常用的热工水力计算软件（如Fluent），这给工程人员带来了不便。

本文针对熔盐堆TMSR-LF1编写了程序MTF。程序MTF抽象工作分解为5个步骤：

1)指定待处理文件：搜索MCNP计算的输出文件。

2)提取所需数据：从输出文件提取F7卡对973个栅元的平均裂变沉积能计数。联合F7卡和SD卡得到的只是每个栅元的裂变沉积能量(MeV)，需要通过式(1)转换为功率密度(W·m^-3)：

式中：p_i是第i个栅元的功率密度，W·m^-3；F_f是第i个栅元裂变沉积能量，MeV；V_i是第i个栅元的体积，m³；P是1/8反应堆热功率，W。

3)数据转换：将F7卡计算的平均裂变沉积能转换为功率密度。

4)搜索UDF样本：搜索提前准备的待完善UDF文件。

5)完善UDF：为UDF样本提供坐标信息和随空间变化的功率密度分布。

程序编写使用语言Python，编译器为Python2.7.7。

程序MTF的验证采用熔盐堆TMSR-LF1单个通道组件为模拟计算对象，对比手动输入功率密度与程序MTF计算结果来验证计算方法是否正确。单个通道组件的几何外形是一个中间流通熔盐的石墨长方体，长宽均为3.6 cm，高110 cm，如图 4。

图 4(Fig.4)

图 4 单个通道组件(a)及流道轮廓(b) Fig.4 Single tube assembly (a) and flow profile (b).

通道边界条件：1)物理边界条件设置如下：组件4个侧面设置为对称边界，组件四周设为真空环境；2)热工边界条件如下：组件侧面边界条件设置为对称边界，进口边界条件设置为速度进口，出口边界条件设置为压力出口。考虑重力影响，设置重力加速度9.8m·s^-2。将通道沿轴向划分为5份，采用MCNP计算各自的功率密度，对比手动输入和程序MTF输入功率密度到Fluent中的计算结果，熔盐温度和密度沿轴向高度分布如图 5。

图 5(Fig.5)

图 5 熔盐温度(a)、密度(b)沿轴向高度的变化 Fig.5 Variation of molten salt temperature (a) and density (b) along the axial height.

从图 5中结果来看，熔盐堆温度沿着流动方向(-0.55-0.55 m)逐渐升高。这是因为在流动过程中，燃料熔盐随着中子裂变不断加热。熔盐通道中间部分的中子通量密度高于两端，所以相对于中间部分，两端温升变缓。图 5中的密度曲线与温度曲线趋势相反，主要因为熔盐密度函数是随温度升高而减小的函数。两种导入方式的熔盐温度的最大差值为0.037，相对误差为3.8×10^-5；熔盐密度的最大差值为0.02，相对误差为7.4×10^-6。造成误差的主要原因是手动输入功率密度时，系统自动对功率密度值小数点后第一位舍入。从图 5中还可以看出，手动输入与程序MTF计算结果曲线符合较好，验证了程序MTF的正确性。

采用程序MTF处理MCNP输出文件获得TMSR-LF1堆芯功率密度分布，如图 6。从图 6中可以看出，堆芯轴向两端的功率密度低于中间部分，但不同轴向高度的功率密度分布趋势一致，其功率密度最小值位于堆芯中心，最大值位于近堆芯中心处或者堆芯边缘。这是由于堆芯中心处是真空辐照通道，靠近中心的熔盐通道中子泄露严重，所以此处的功率密度最小。此外，TMSR-LF1带有反射层，靠近反射层处，部分中子自反射层返回堆芯，使边缘的功率密度增大。

图 6(Fig.6)

图 6 不同轴向高度的功率密度分布 Fig.6 Power density distribution of the core in different axial heights.

通过上述计算，可以得到环境温度5℃、18℃、25℃、30℃、35℃、40℃下堆芯到主屏蔽各部分的温度分布，如图 7(a)。从图 7(a)中看出，随着径向距离增加，环境温度对其温度分布影响变大。主屏蔽受环境温度影响最大，这是因为主屏蔽与环境直接接触。为更清晰看到在不同环境温度下主屏蔽的温度分布曲线，给出图 7(b)。从图 7(b)中看出，不同环境温度下的主屏蔽混凝土层温度极值不同，但是温度下降斜率相同。原因是混凝土层的传热以导热为主，根据傅里叶定律可知，在热流密度和导热系数一定时，温度梯度不变（即斜率）。环境温度40℃下，反应堆的温度分布如图 7(c)。

图 7(Fig.7)

图 7 TMSR-LF1 (a)、主屏蔽(b)和环境温度40°C的温度分布(c) Fig.7 Temperature distribution of TMSR-LF1 (a), the shielding (b) and TMSR-LF1 in the condition temperature 40°C (c).

环境温度5℃、18℃、25℃、30℃、35℃、40℃下的主屏蔽混凝土墙体的温度分布，分别对应图 8(a)-(f)。其各自最高边界温度33.92-67.42℃，厚度方向最大温差72.58-78.40℃·m^-1。根据对反应堆普通混凝土屏蔽墙的两个要求：边界最高温度不能超过70℃以及沿混凝土屏蔽体厚度方向最大温差小于100℃·m^-1。从图 8中可以看出，随着环境温度的上升，主屏蔽混凝土墙温度增大。同时从表 1可知，即使环境温度达到40℃，混凝土墙的最高温度67.42℃也不会超过规定限值70℃。因此，熔盐堆在不同环境温度下都符合边界最高温度的要求。从表 1可以看出，随着环境温度增加，沿厚度方向的最大温差在减小，其最大值78.4℃·m^-1仍然小于100℃·m^-1，达到设计要求。

图 8(Fig.8)

图 8 不同环境温度下主屏蔽墙体温度分布 (a) 5°C，(b) 18°C，(c) 25°C，(d) 30°C，(e) 35°C，(f) 40°C Fig.8 Temperature distribution of the concrete wall in different temperature. (a) 5°C, (b) 18°C, (c) 25°C, (d) 30°C, (e) 35°C, (f) 40°C

表 1(Table 1)

表 1 主屏蔽混凝土墙温度分析Table 1 Temperature analysis of the concrete wall. 环境 Condition / °C 外壁 Outer wall / °C 内壁 Inner wall / °C 径向单位长度最大温差 Radial temperature contrast unit length / °C·m-1 5 18 25 13.91 26.71 33.65 33.92 46.27 53.05 78.40 76.20 75.00 30 38.58 57.85 74.20 35 43.52 62.67 73.60 40 48.45 67.42 72.58

表 1 主屏蔽混凝土墙温度分析 Table 1 Temperature analysis of the concrete wall.

本文基于现阶段2MWth液态熔盐堆TMSR-LF1的概念设计中给出的主屏蔽结构，建立了TMSR-LF1的物理和热工模型，计算分析了主屏蔽的温度分布。通过软件MCNP计算反应堆稳定工况下的功率密度空间分布，采用自行编写程序MTF将功率密度导入软件Fluent进行三维温度场计算。在混凝土屏蔽墙的温度场计算中，采用了环境温度为5℃、18℃、25℃、30℃、35℃、40℃不同边界条件分别计算。得到主屏蔽最高边界温度33.92-67.42℃，厚度方向最大温差为72.58-78.40℃·m^-1，所有计算结果的边界最高温度和沿厚度方向的最大温差均低于要求限值，2 MWth熔盐堆TMSR-LF1主屏蔽满足温度设计要求。