CN111414722B - 一种核反应堆堆芯物理与热工耦合的模拟方法 - Google Patents

Abstract

本发明公开了一种核反应堆堆芯物理与热工耦合的模拟方法，首先采用Fluent前处理软件建立计算域非结构化网格模型，随后将网格模型导入Fluent中建立算例；采用Fluent软件的用户自定义函数功能提取Fluent算例中的网格信息及流场数据，自动生成MCNP5计算所需的输入文件；借助Fluent的用户自定义函数功能控制MCNP5完成中子物理场的计算，解析其输出文件，提取核燃料区的裂变功率数据；将裂变功率传入Fluent软件中并执行计算；当Fluent计算初步收敛后，使用其用户自定义函数重新生成MCNP5输入文件,再次调用MCNP5计算裂变功率数据并传递给Fluent；重复以上步骤直至Fluent最终收敛。本发明可以方便的实现物理热工耦合，为数值反应堆的研究奠定了基础。

Description

一种核反应堆堆芯物理与热工耦合的模拟方法

技术领域

本发明涉及核反应堆堆芯设计与安全分析领域，具体涉及一种核反应堆堆芯物理与热工耦合的模拟方法。

背景技术

反应堆中子物理场与冷却剂流场是相互依赖、相互制约的。在反应堆实际运行过程中，当反应堆功率变化时，堆芯温度及其分布也要发生变化，从而引起燃料温度、慢化剂密度、中子截面和可溶解硼溶解度的变化。这些参数的变化均将导致堆芯有效增值因子的变化，进而作用于堆芯功率。由于中子物理场和冷却剂热工水力流场存在显著的时间尺度、空间维度上的差异，且这种反馈具有显著的非线性特征，在早期的分析中被迫引入人为假设，将两类物理场的解耦单独进行分析。

随着核反应堆理论研究的不断深入和计算机性能的迅速提升，进行多物理场耦合仿真、消除早期反应堆设计分析中引入的近似与误差逐渐变得现实可行。自2010年美国能源部提出数值反应堆计划以来，研究人员不断尝试用各种方法对反应堆内中子物理场与流场的物理热工耦合现象进行分析。

其中基于有限体积理论与随机输运理论的耦合被视为精度最高的方法，有望彻底解决反应堆内的物理热工耦合问题。由于随机输运方法通常采用构造实体方法刻画计算域，这与有限理论采用的非结构网格有着本质差异，这一差异导致两类程序间的控制体很难一致、程序间数据传递困难。具体而言表现在两方面：控制体分割困难，在对计算域划分控制体时需要兼顾到另一程序的控制划分方法，需要引入大量附加操作；数据映射困难，需要针对计算域的具体情况，设计将要交换的数据的结构。

发明内容

为了克服上述现有技术存在的问题，本发明的目的在于提供一种核反应堆堆芯物理与热工耦合的模拟方法，降低耦合建模难度，为实现高精度的物理热工耦合分析提供了基础。

为了达到上述目的，本发明采用如下技术方案：

一种核反应堆堆芯物理与热工耦合的模拟方法，基于蒙特卡罗中子输运程序与计算流体动力学程序耦合，包括如下步骤：

步骤1：针对核反应堆燃料组件结构，依次采用计算流体动力学程序中的前处理软件即几何建模软件和网格划分软件分别建立燃料组件的几何模型和网格模型，其中网格模型需为全六面体；

步骤2：将步骤1中的网格模型导入到计算流体动力学程序中，根据软件使用说明书对计算流体动力学程序进行设置；在设置监视器时，设置监视器监测计算域最高温度，用作后续计算收敛的判定依据；采用计算流体动力学程序用户自定函数对燃料区能量方程增加源项；

步骤3：利用计算流体动力学程序的用户自定函数，提取计算流体动力学程序模型中的网格及流场参数，生成蒙特卡罗中子输运程序计算所需输入文件，具体分为以下步骤：

步骤3-1：利用计算流体动力学程序的用户自定义函数遍历计算流体动力学程序算例中的网格模型，获取各网格单元的编号、几何中心坐标、体积、材料密度、温度、各侧面中心点坐标及其面法向量；

步骤3-2：利用步骤3-1中获得的网格单元各侧面中心点坐标及面法向量构建平面方程，利用这些平面方程在蒙特卡罗中子输运程序的输入文件中重构对应网格单元；

步骤3-3：根据步骤3-1中获得的网格单元内材料密度判断材料类别，即判断材料为裂变材料、包壳、冷却剂还是慢化剂，由此确定蒙特卡罗中子输运程序对应网格单元中材料的核素构成；当判定材料为裂变材料时，在蒙特卡罗中子输运程序的输入文件中对该网格单元设置裂变功率监测器；

步骤3-4：根据步骤3-1获得的网格单元内温度对蒙特卡罗中子输运程序输入文件中对应网格单元的材料和截面库进行平方根插值，插值公式如下：

∑(T)＝f_L∑(T_L)+(1-f_L)∑(T_H)

其中T为网格单元内温度，T_L为低温核截面库温度，T_H为高温核截面库，以上这些温度单位均为K；∑(T)为插值后核截面库，∑(T_L)为低温核截面库，∑(T_H)为高温核截面库，f_L为插值后截面库中低温核截面库所占百分数，

步骤4：使用计算流体动力学程序的用户自定义函数控制蒙特卡罗中子输运程序读入步骤3中生成的输入文件，完成中子物理场的计算，得到裂变能；

步骤5：使用计算流体动力学程序的用户自定函数解析蒙特卡罗中子输运程序的输出文件，提取裂变数据并传递给计算流体动力学程序，具体步骤如下：

步骤5-1：将步骤4中蒙特卡罗中子输运程序中计算得到的裂变能转换成计算流体动力学程序计算所需的体积释热率，转换公式如下：

式中：Q_f为计算流体动力学程序计算所需的体积释热率，单位为W/m³，H_f为蒙卡核物理计算软件统计获得了每个燃料控制体的裂变能，单位为MeV/g，P为核燃料组件设计热功率，v为每次裂变释放的平均中子数；Q每次裂变释放的能量，单位为MeV；ρ_f燃料密度；k_eff有效反应性系数；

步骤5-2：通过匹配计算流体动力学程序控制体编号与蒙特卡罗中子输运程序控制体编号，将步骤5-1获得的体积释热率传递给计算流体动力学程序；

步骤6：在步骤5的基础上，采用计算流体动力学程序进行迭代计算，当步骤2中设置的监视器监测计算域最高温度变化小于1K时，认为计算收敛，完成计算；若设置的监视器监测计算域最高温度大于等于1K时，进入步骤7；

步骤7：使用计算流体动力学程序的用户自定义函数获取计算流体动力学程序软件的连续性残差，当连续性残差大于等于10^-3，继续计算流体动力学程序的迭代；当连续性残差小于10^-3，重复步骤3-步骤6。

有益效果：

相比常见的物理与热工耦合处理，本发明带来以下有益效果：

1)本发明方法提供了一种反应堆堆芯物理与热工的高精度耦合分析方法，能提供高精度的反应堆局部信息，指导反应堆燃料元件设计；

2)本发明方法利用计算流体动力学程序的用户自定义函数自动将计算流体动力学程序算例模型映射成MCNP5求解模型，无需用户自行建立蒙特卡罗中子输运程序的求解模型，减少了用户的工作量；

3)本发明方法在耦合过程中，实现了计算流体动力学程序与蒙特卡罗中子输运程序的一一映射，避免了早期耦合程序中的均匀化近似处理，可以保证耦合计算精度；

4)本发明方法仅采用了一种脚本语言——计算流体动力学程序的自带用户自定义函数，降低了对用户的编程语言要求。

附图说明

图1本发明核反应堆堆芯物理与热工耦合的模拟方法流程图。

图2板型燃料组件几何模型，其中，图a为全局图，图b为局部图。

图3板型燃料组件网格模型。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。

以下结合图1所示流程图，以蒙特卡罗中子输运程序与计算流体动力学程序耦合计算板型燃料为例为例，对本发明作进一步的详细描述：以下将蒙特卡罗中子输运程序称作MCNP5程序，将计算流体动力学程序称作Fluent程序。

步骤1：针对板型燃料组件结构，采用Fluent程序中的前处理软件即几何建模软件ANSYS DesignModeler建立燃料组件的几何模型，对位于燃料与冷却剂间的薄包壳如实刻画，如图2中图a和图b所示；采用Fluent程序中的网格划分软件ANSYS Meshing程序对其进行网格划分，建立网格模型时，指定模型中各边线上网格节点数，生成纯六面体网格，如图3所示；

步骤2：将步骤1中的网格模型导入到Fluent程序中，根据软件使用说明书对Fluent程序进行设置；在设置监视器时，设置监视器监测计算域最高温度，用作后续计算收敛的判定依据；采用Fluent程序用户自定函数对燃料区能量方程增加源项；

步骤3：利用Fluent程序的用户自定函数，提取计算流体动力学程序模型中的网格及流场参数，生成MCNP5程序计算所需输入文件，具体分为以下步骤：

步骤3-1：利用Fluent程序的用户自定义函数遍历Fluent程序算例中的网格模型，获取各网格单元的编号、几何中心坐标、体积、密度、温度、各侧面中心点坐标及其面法向量；

步骤3-2：利用步骤3-1中获得的网格单元侧面中心点坐标及面法向量构建平面方程，利用这些平面方程在MCNP5程序的输入文件中重构对应网格单元；

步骤3-3：根据步骤3-1中获得的网格单元内材料的密度数值判断材料类别，即判断材料为裂变材料、包壳、冷却剂还是慢化剂，由此自动确定MCNP5程序对应网格单元中材料的核素构成；当判定材料为裂变材料时，在MCNP5程序的输入文件中对该网格单元设置裂变功率监测器；

步骤3-4：根据步骤3-1获得的网格单元内温度对MCNP5程序输入文件中对应网格单元的材料和截面库进行平方根插值；

步骤4：使用Fluent程序的用户自定义函数控制MCNP5程序读入步骤3中生成的输入文件，完成中子物理场的计算，得到裂变能；

步骤5：使用Fluent程序的用户自定函数解析MCNP5程序的输出文件，提取裂变数据并传递给Fluent程序，具体步骤如下：

步骤5-1：将步骤4中MCNP5程序中计算得到的裂变数据转换成Fluent程序计算所需的体积释热率，转换公式如下：

式中：Q_f为Fluent程序计算所需的体积释热率，单位为W/m³，H_f为蒙卡核物理计算软件统计获得了每个燃料控制体的裂变能，单位为MeV/g，P为核燃料组件设计热功率，v为每次裂变释放的平均中子数；Q每次裂变释放的能量，单位为MeV；ρ_f燃料密度；k_eff有效反应性系数；

步骤5-2：通过匹配Fluent程序控制体编号与MCNP5程序控制体编号，将步骤5-1获得的体积释热率传递给Fluent程序；

步骤6：在步骤5的基础上，采用Fluent程序进行迭代计算，当步骤1中设置的监视器监测计算域最高温度变化小于1K时，认为计算收敛，完成计算；若设置的监视器监测计算域最高温度大于等于1K时，进入步骤7；

步骤7：使用Fluent程序的用户自定义函数获取Fluent程序的连续性残差，当连续性残差大于等于10^-3，继续Fluent程序的迭代；当连续性残差小于10^-3，重复步骤3-步骤6。

本发明未详细说明部分属本领域技术人员公知常识。在以上步骤中涉及若干软件，亦可采用相似软件替换。

Claims (1)

1.一种核反应堆堆芯物理与热工耦合的模拟方法，其特征在于，基于蒙特卡罗中子输运程序与计算流体动力学程序耦合，包括如下步骤：

步骤1：针对核反应堆燃料组件结构，依次采用计算流体动力学程序中的前处理软件即几何建模软件和网格划分软件分别建立燃料组件的几何模型和网格模型，其中网格模型需为全六面体；

步骤2：将步骤1中的网格模型导入到计算流体动力学程序中，根据软件使用说明书对计算流体动力学程序进行设置；在设置监视器时，设置监视器监测计算域最高温度，用作后续计算收敛的判定依据；采用计算流体动力学程序用户自定函数对燃料区能量方程增加源项；

步骤3：利用计算流体动力学程序的用户自定函数，提取计算流体动力学程序模型中的网格及流场参数，生成蒙特卡罗中子输运程序计算所需输入文件，具体分为以下步骤：

步骤3-1：利用计算流体动力学程序的用户自定义函数遍历计算流体动力学程序算例中的网格模型，获取各网格单元的编号、几何中心坐标、体积、材料密度、温度、各侧面中心点坐标及其面法向量；

步骤3-2：利用步骤3-1中获得的网格单元各侧面中心点坐标及面法向量构建平面方程，利用这些平面方程在蒙特卡罗中子输运程序的输入文件中重构对应网格单元；

步骤3-3：根据步骤3-1中获得的网格单元内材料密度判断材料类别，即判断材料为裂变材料、包壳、冷却剂还是慢化剂，由此确定蒙特卡罗中子输运程序对应网格单元中材料的核素构成；当判定材料为裂变材料时，在蒙特卡罗中子输运程序的输入文件中对该网格单元设置裂变功率监测器；

步骤3-4：根据步骤3-1获得的网格单元内温度对蒙特卡罗中子输运程序输入文件中对应网格单元的材料和截面库进行平方根插值，插值公式如下：

∑(T)＝f_L∑(T_L)+(1-f_L)∑(T_H)

其中T为网格单元内温度，T_L为低温核截面库温度，T_H为高温核截面库温度，以上这些温度单位均为K；∑(T)为插值后核截面库，∑(T_L)为低温核截面库，∑(T_H)为高温核截面库，f_L为插值后截面库中低温核截面库所占百分数，

步骤4：使用计算流体动力学程序的用户自定义函数控制蒙特卡罗中子输运程序读入步骤3中生成的输入文件，完成中子物理场的计算，得到裂变能；

步骤5：使用计算流体动力学程序的用户自定函数解析蒙特卡罗中子输运程序的输出文件，提取裂变数据并传递给计算流体动力学程序，具体步骤如下：

步骤5-1：将步骤4中蒙特卡罗中子输运程序中计算得到的裂变能转换成计算流体动力学程序计算所需的体积释热率，转换公式如下：

式中：Q_f为计算流体动力学程序计算所需的体积释热率，单位为W/m³，H_f为蒙卡核物理计算软件统计获得了每个燃料控制体的裂变能，单位为MeV/g，P为核燃料组件设计热功率，v为每次裂变释放的平均中子数；Q每次裂变释放的能量，单位为MeV；ρ_f燃料密度；k_eff有效反应性系数；

步骤5-2：通过匹配计算流体动力学程序控制体编号与蒙特卡罗中子输运程序控制体编号，将步骤5-1获得的体积释热率传递给计算流体动力学程序；

步骤6：在步骤5的基础上，采用计算流体动力学程序进行迭代计算，当步骤2中设置的监视器监测计算域最高温度变化小于1K时，认为计算收敛，完成计算；若设置的监视器监测计算域最高温度大于等于1K时，进入步骤7；

步骤7：使用计算流体动力学程序的用户自定义函数获取计算流体动力学程序软件的连续性残差，当连续性残差大于等于10^-3，继续计算流体动力学程序的迭代；当连续性残差小于10^-3，重复步骤3-步骤6。

Images