一种基于反距离加权插值的核热强耦合方法与流程

1.本发明涉及核反应堆堆芯设计与安全分析领域，具体涉及一种基于反距离加权插值的核热强耦合方法。


背景技术：

2.反应堆多物理场耦合是对反应堆进行综合分析的一项关键技术，在传统方式下，将中子物理计算与热工流体计算分割开来，得到的结果通常是偏保守的，无法精确反映反应堆或堆内组件的真实物理情况与现象。因此将中子学计算程序与热工流体计算程序耦合起来进行核热耦合的稳态分析将为反应堆设计提供更加精细的物理参数与相关设计准则。
3.主流的蒙特卡罗方法输运计算程序主要有mcnp、openmc、serpent以及scale；主流的计算流体力学(computational fluid dynamics,cfd)软件主要有：ansys fluent、ansys cfx以及star
‑
ccm 。目前核热耦合程序大多采用mcnp与fluent、mcnp与star
‑
ccm 、mcnp与openfoam进行耦合，其中二者的网格映射是一大难点，一般使用的方法为网格一一映射：即将mcnp与fluent使用同一套网格，二者的网格数量、网格栅元的中心坐标与各个面的坐标均完全相同；但该方法缺点是在目标结构非常复杂时，会造成mcnp生成网格的困难，并且mcnp在输入卡上有最大栅元数目与最大曲面数目的限制，所以在流体力学计算需要较大数目网格量时，为了保证mcnp的成功运行，使用网格一一映射方法不得不将流体力学计算的网格数量进行下调，这样可能会造成流体力学计算网格质量较差，误差较大的情况。
4.还有一种较为常用的方法为体积权重方法，该方法较上面的方法有较大改进，中子学计算与流体力学计算不在使用同一套网格，而是在mcnp上使用粗网格，在流体力学计算上采用细网格，这样可以在一定程度上克服上述中受mcnp网格最大数量限制的因素，但该方法在求解体积与权重因子时会造成计算量大，求解困难的情况；并且该方法容易受到几何模型与网格划分方式的影响，若几何模型较为复杂或者mcnp与流体力学计算网格中有任意一方采用非结构化网格，则在精细耦合求解时会在求解体积上进行大量的运算，由此而造成一定的计算困难。
5.反距离加权插值算法也可称为距离倒数乘算法。其计算结果不依赖于相关的数学条件。目标节点的值是其与插值节点之间的距离和插值节点的值进行加权平均后计算得到的。当插值节点与目标节点相距较远时，证明该插值节点对于目标节点的影响较小，故该权值较小；反之，当插值节点与目标节点相距较近时，则证明该插值节点对于目标节点的影响较大，故该权值较大。
6.在使用反距离加权插值的三维核热耦合方法时，关键点在于网格的中心点坐标，所以反距离加权插值可以使用粗细两种网格进行计算。其中ansys fluent计算采用精细的细网格，mcnp采用较为粗糙的粗网格进行计算。这样可以有效提高mcnp的运行速度并且提高ansys fluent计算结果的精度，并且突破mcnp的最大栅元数目的限制。在反距离加权插值算法进行的核热耦合的过程中，几何体的相对形状则是次要的，只需要保证粗细网格的原点与相对位置是相同的即可，该方法在处理复杂几何时会非常有效。
7.综上所述，反距离加权插值算法的核热耦合具有以下几点优势：
8.1.突破mcnp最大栅元数与最大曲面数限制；
9.2.可以划分更加精细的ansys fluent网格模型来进行更加精确的流体力学计算；
10.3.提高mcnp运行速度；
11.4.反距离加权插值算法与几何体的形状无关，更加高效地处理复杂几何；
12.中国先进实验快堆的板型燃料组件的流道为闭式流道，组件内冷却剂流道间几乎不存在冷却剂质量和动量的交换，组件间的水隙也较小，板型燃料组件中冷却剂流速很高，传热性能也随之较大；其与相同体积下的棒束型燃料组件相比，中子通量密度更高，工作温度较低，燃耗较深，为了进行更加精细的求解板型燃料元件与组件中的中子物理与热工水力情况，开发了基于反距离加权插值的精细三维核热强耦合方法。


技术实现要素：

13.为了克服上述现有技术存在的问题，本发明的目的在于提供一种基于反距离加权插值的核热强耦合方法，应用本发明的方法所得到的数据可以为反应堆堆内燃料元件、组件、堆芯等相关设计研究工作提供参考。
14.为了达到上述目的，本发明采用如下技术方案：
15.一种基于反距离加权插值的核热强耦合方法，包括如下步骤：
16.步骤一：使用前处理建模软件、计算流体力学软件以及蒙特卡洛方法输运计算程序进行计算对象的模型构造与网格划分；使用前处理建模软件构造核燃料组件几何模型，并将该几何模型导入计算流体力学软件中，再使用计算流体力学软件划分核燃料组件的精细结构化网格并进行流体力学计算；使用蒙特卡罗方法输运计算程序对核燃料组件进行建模，模型中包含燃料元件、冷却剂与包壳；通过计算流体力学软件中的用户自定义函数对应用于蒙特卡罗方法输运计算程序中的核燃料组件模型完成粗网格自动划分；其中，通过计算流体力学软件中的用户自定义函数对应用于蒙特卡罗方法输运计算程序中核燃料组件模型完成粗网格自动划分的具体步骤如下：
17.步骤1：向计算流体力学软件的用户自定义函数中指定核燃料组件几何模型的尺寸以及材料；
18.步骤2：指定在使用计算流体力学软件的用户自定义函数对核燃料组件几何模型进行网格划分时坐标原点位置；
19.步骤3：指定x方向、y方向以及z方向的网格划分数目，通过计算流体力学软件的用户自定义函数获取核燃料组件的网格中心点坐标、曲面位置坐标；
20.步骤4：依据步骤3中得到的网格中心点坐标以及曲面位置坐标，生成用于中子物理计算的蒙特卡洛输运计算程序的曲面输入卡；
21.步骤5：依据步骤4中生成的曲面输入卡，指定相对应的布尔运算规则，生成用于中子物理计算的蒙特卡洛方法输运计算程序的栅元输入卡；
22.步骤6：依据步骤5中生成的栅元输入卡，通过指定不同区域材料标号生成用于中子物理计算的蒙特卡洛方法输运计算程序的材料输入卡；
23.步骤7：指定本征值计算中模拟粒子数、迭代次数以及裂变源项分布生成本征值计算输入卡；
24.步骤8：指定裂变能量统计计数，生成裂变能量计数输入卡；
25.步骤9：对步骤4中生成的曲面输入卡、步骤5中生成的栅元输入卡、步骤6中生成的材料输入卡、步骤7中生成的本征值计算输入卡以及步骤8中生成的裂变能量计数输入卡进行合并，得到最终的的蒙特卡洛方法输运计算程序的输入卡；
26.步骤二：分别设定计算流体力学软件、蒙特卡洛方法输运计算程序中对于核燃料组件的边界条件；其中计算流体力学软件设置组件为速度入口边界条件以及压力为零的出口边界条件，组件四周为对称边界条件；蒙特卡洛方法输运计算程序设置组件四周为全反射边界条件，上下两端为真空边界条件；
27.步骤三：使用计算流体力学软件进行流体力学计算，待连续性残差降至1e
‑
3时即认为流体力学计算已经收敛，然后通过用户自定义函数提取流体力学计算结果中目标几何体各个网格的温度与密度，再利用反距离加权插值算法进行计算，将流体力学计算模型中的精细网格上的数据映射进蒙特卡洛方法输运计算程序的粗糙网格中；
28.反距离加权插值算法：
29.如果d(x,x
i
)≥err，则：
[0030][0031]
否则：
[0032]
u(x)＝u
i
[0033]
其中：
[0034][0035]
式中：
[0036]
i为插值节点的索引；
[0037]
d(x,x
i
)为目标节点与插值节点的距离；
[0038]
err为一个极小正数，用于计算机中的浮点判断；
[0039]
w
i
(x)为位置x处且索引为i的节点的权值；
[0040]
n为插值节点的个数；
[0041]
u
i
为索引为i的插值节点值；
[0042]
u(x)为目标节点值；
[0043]
p为幂值；
[0044]
步骤四：依据步骤三中得到的温度数据与密度数据，通过温度均方根插值函数计算所选温度区间内核素的截面数据；
[0045]
温度均方根插值公式为：
[0046][0047]
σ(t)＝f
l
σ(t
l
) (1
‑
f
l
)σ(t
h
)
[0048]
式中：
[0049]
t为核燃料组件网格温度；
[0050]
t
h
为所选温度区间的最高温度/k；
[0051]
t
l
为所选温度区间的最低温度/k；
[0052]
σ(t
h
)为所选温度区间最高温度的截面数据；
[0053]
σ(t
l
)为所选温度区间最低温度的截面数据；
[0054]
σ(t)为宏观截面/m
‑1；
[0055]
f
l
为该温度区间中最低温度的和数据库中所占份额；
[0056]
通过温度均方根插值公式进行蒙特卡洛方法输运计算程序的截面数据选取与生成，进而生成蒙特卡洛方法输运计算程序的输入卡，调用蒙特卡洛方法输运计算程序进行计算；
[0057]
步骤五：蒙特卡洛方法输运计算程序计算完成后，读取输出文件中的有效增殖因子k
eff
与归一化功率分布，随后通过功率密度分配函数进行功率分配，来进行核燃料组件的真实功率的计算：
[0058]
功率密度分配函数为：
[0059][0060]
式中：
[0061]
g
i
为蒙特卡洛方法输运计算程序输出文件中每个栅元计算归一化结果；
[0062]
∑g
i
为所有归一化结果的求和；
[0063]
p为元件或组件的总功率密度；
[0064]
p
i
为每一个栅元加权计算后的真实功率密度；
[0065]
将利用上式计算完毕后的每一个栅元的真实功率密度进行反距离加权插值，即实现了从蒙特卡洛方法输运计算程序的粗网格中向流体力学计算模型的细网格进行网格映射与数据传输；插值过程中利用用户自定义函数中的内存结构进行数据存储，然后通过用户自定义函数对流体力学计算模型中指定区域赋予能量源项，用以进行下一步流体力学的计算工作；
[0066]
步骤六：提取几次耦合计算中蒙特卡洛方法输运计算程序计算完成后的有效增殖因子k
eff
，若几次迭代的k
eff
波动范围小于1e
‑
3，则认为该耦合计算已经收敛，计算完毕，反之则返回步骤一开始下一次迭代计算。
[0067]
和现有技术相比较，本发明具备如下优点：
[0068]
1.突破蒙特卡罗方法输运计算程序中最大栅元数与最大曲面数限制；
[0069]
2.可以划分更加精细的计算流体力学软件网格模型来进行更加精确的流体力学计算；
[0070]
3.提高蒙特卡罗方法输运计算程序的运行速度；
[0071]
4.反距离加权插值算法与几何体的形状无关，可以高效地处理复杂几何；
附图说明
[0072]
图1为本发明的方法流程图；
[0073]
图2为中国先进实验快堆板型燃料组件几何模型；
[0074]
图3为中国先进实验快堆板型燃料组件网格模型；
具体实施方式
[0075]
下面结合附图与具体实施案例中国先进实验快堆板型燃料组件对本发明做进一步详细说明，案例中使用的三维建模软件为ansys design modeler，计算流体力学软件为ansys fluent，蒙特卡洛方法输运计算程序为mcnp。
[0076]
如图1所示，本实施例基于反距离加权插值的精细三维核热强耦合方法，包括如下步骤：
[0077]
步骤一：使用ansys design modeler、ansys fluent以及mcnp进行板型燃料组件的模型构造与网格划分；ansys design modeler构造用于流体力学计算的板型燃料组件几何模型，如图2所示；ansys fluent构造用于流体力学计算的板型燃料组件网格模型，如图3所示；mcnp构造用于中子物理计算的几何模型与网格模型，其中中子物理计算的网格模型使用ansys fluent中的用户自定义函数完成构造，网格模型的构造具体步骤如下：
[0078]
步骤1：向用户自定义函数中指定板型燃料组件的几何模型尺寸以及材料；
[0079]
步骤2：指定使用ansys fluent中用户自定义函数进行网格划分时板型燃料组件的坐标原点位置；
[0080]
步骤3：指定x方向、y方向以及z方向的网格划分数目，通过ansysfluent中的用户自定义函数获取板型燃料组件的网格中心点坐标与曲面位置坐标；
[0081]
步骤4：依据步骤3中得到的网格中心点坐标以及曲面位置坐标，生成mcnp的曲面输入卡；
[0082]
步骤5：依据步骤4中生成的曲面输入卡，指定相对应的布尔运算规则，生成mcnp的栅元输入卡；
[0083]
步骤6：依据步骤5中生成的栅元输入卡，通过指定不同区域材料标号生成mcnp的材料输入卡；
[0084]
步骤7：指定本征值计算中模拟粒子数、迭代次数以及裂变源项分布生成本征值计算输入卡；
[0085]
步骤8：指定裂变能量统计计数，生成裂变能量计数卡；
[0086]
步骤9：对步骤4中生成的曲面输入卡、步骤5中生成的栅元输入卡、步骤6中生成的材料输入卡、步骤7中生成的本征值计算输入卡以及步骤8中生成的裂变能量计数卡进行合并，得到mcnp的整体输入卡；
[0087]
步骤二：设定ansys fluent与mcnp对于板型燃料组件相对应的物理边界条件；其中ansys fluent设置板型燃料组件为速度入口边界条件、压力为零的出口边界条件，四周为对称边界条件；mcnp设置组件四周为全反射边界条件，上下两端为真空边界条件；
[0088]
步骤三：首先使用ansys fluent进行流体力学计算，使用ansys fluent的用户自定义函数判断流体力学计算中的连续性残差是否小于1e
‑
3，若小于1e
‑
3，则认为此时流体力学计算已经收敛，结束ansys fluent的流体力学计算；ansys fluent计算完毕后，通过用户自定义函数提取ansys fluent中目标几何体各个网格的温度与密度，再利用反距离加权插值算法将ansys fluent中板型燃料组件的精细网格上的数据映射进mcnp中板型燃料组件的粗糙网格中；
[0089]
反距离加权插值算法：
[0090]
如果d(x,x
i
)≥err，则：
[0091][0092]
否则：
[0093]
u(x)＝u
i
[0094]
其中：
[0095][0096]
式中：
[0097]
i为插值节点的索引；
[0098]
d(x,x
i
)为目标节点与插值节点的距离；
[0099]
err为一个极小正数，用于计算机中的浮点判断；
[0100]
w
i
(x)为位置x处且索引为i的节点的权值；
[0101]
n为插值节点的个数；
[0102]
u
i
为索引为i的插值节点值；
[0103]
u(x)为目标节点值；
[0104]
p为幂值；
[0105]
步骤四：依据步骤三中得到的温度数据与密度数据，通过温度均方根插值函数计算所选温度区间内的核素截面数据；
[0106]
温度均方根插值公式为：
[0107][0108]
σ(t)＝f
l
σ(t
l
) (1
‑
f
l
)σ(t
h
)
[0109]
式中：
[0110]
t为核燃料组件网格温度；
[0111]
t
h
为所选温度区间的最高温度/k；
[0112]
t
l
为所选温度区间的最低温度/k；
[0113]
σ(t
h
)为所选温度区间最高温度的截面数据；
[0114]
σ(t
l
)为所选温度区间最低温度的截面数据；
[0115]
σ(t)为宏观截面/m
‑1；
[0116]
f
l
为该温度区间中最低温度的和数据库中所占份额；
[0117]
通过温度均方根插值公式进行核素截面数据选取与生成，进而生成mcnp输入卡，调用mcnp进行中子物理计算；
[0118]
步骤五：mcnp计算完成后，读取输出文件中的有效增殖因子k
eff
与功率分布，mcnp输出中为归一化的功率密度，需要进行真实功率的计算，故通过功率密度分配函数进行分配功率：
[0119]
功率密度分配函数为：
[0120][0121]
式中：
[0122]
g
i
为mcnp输出文件中每个栅元计算归一化结果；
[0123]
∑g
i
为所有归一化结果的求和；
[0124]
p为元件或组件的总功率密度；
[0125]
p
i
为每一个栅元加权计算后的真实功率密度；
[0126]
将利用上式计算完毕后的每一个栅元的真实功率密度进行反距离加权插值，即从mcnp的板型燃料组件的粗网格向ansys fluent的板型燃料组件的细网格进行网格映射与数据传输并赋予能量源项，用以进行下一步ansys fluent的计算工作。
[0127]
步骤六：提取几次耦合计算中mcnp计算完成后的k
eff
，若几次的k
eff
波动范围小于1e
‑
3，则认为该耦合计算已经收敛，计算完毕，反之则返回步骤一开始下一次迭代计算。