CN110598324A - 一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法 - Google Patents

Abstract

本发明公开了一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法，步骤如下：1、建立核反应堆弥散型板型燃料组件流体域计算模型；2、建立核反应堆弥散型板型燃料组件固体域计算模型；3、将堆芯燃耗过程中辐照效应以子程序形式编写进固体有限元分析程序中；4、对第一步中的模型开展稳态计算；5、利用中间数据交换接口程序调取第四步结果，调用第一、二步中的计算模型，流体动力学计算程序和第三步中二次开发后的固体有限元分析程序之间不断交换数据，迭代计算；6、迭代计算到设定好的计算终止时间，计算停止。通过计算获得核反应堆板型燃料元件堆芯辐照条件下，流固两侧详尽的瞬时热工水力现象，对核反应堆的设计与安全分析具有重要意义。

Description

一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法

技术领域

本发明属于方法发明技术领域，具体涉及一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法。

背景技术

弥散型板型燃料元件优点诸多，在研究型实验堆、舰艇反应堆和小型动力堆上具有广泛的应用。在反应堆运行后期，由于快中子辐照及裂变产物的积累将产生明显的辐照肿胀，加上燃料芯块及包壳的热膨胀，燃料元件厚度和流道的宽度将发生变化，导致燃料的流量分配、燃料板的温度场及冷却剂出口温度偏离反应堆设计值。非均匀辐照条件下，燃料组件受到热工、机械和辐照的作用，发生热膨胀、机械变形、辐照损伤等行为，使得燃料组件发生宏观变形，进而引起燃料组件内部流道间隙尺寸的变化，使冷却剂流动状态发生改变，反过来会燃料组件的热工、机械和材料性能产生影响，这是一个具有反馈作用的流固耦合问题。过去计算资源相对匮乏，流固耦合一般采用单向耦合的方式，即忽略了固体的几何变形对流体速度场及温度场的影响。单向流固耦合是对双向流固耦合的一种简化，在单相耦合计算过程中，固体域计算得出的变形位移不会作用到流体域上，仅支持流体域参数向固体域单向传递，然而在真实的物理过程中，数据应该在流体和固体之间迭代传输，双向流固耦合契合真实的物理过程，在计算的过程中，流体域计算得到的压力、温度会通过耦合面传递给固体，而固体得到的壁面位移、壁面温度也通过耦合面传递到流体域。

在过去几十年，研究人员主要关注于燃料元件的热力耦合行为，一些1.5D和2D的燃料力学分析程序被开发出来，例如美国的PRAPCON,日本的FEMAX,西安交通大学的PROBA。随着对燃料力学行为精细化模拟的需求，研究工作者的关注点开始向3D程序过渡，例如基于MOOSE开发出来的BISON以及BEEs程序，基于COMSOL开发出来的CAMPUS程序。然而，这些程序仅仅将流体简化成零维流动或一维流动，并没有考虑六体侧复杂的热工水力行为。对于非均匀辐照的燃料组件，考虑三维流体与燃料组件辐照变形行为之间的相互作用，即双向流固耦合可以使辐照变形计算结果更加符合实际。

双向流固耦合是随着计算机能力的不断提升而逐渐发展起来的一种更精确的求解技术，目前正处于蓬勃发展期，然而，对于目前国内外许多有限元分析软件或程序代码而言，例如ABAQUS、ANSYS、COMSOL等，这些软件和程序都没有考虑辐照效应带来的辐照肿胀、密实化即蠕变等效应。为了可以完整的实现考虑辐照影响的弥散型板型燃料元件流固耦合计算，需要考虑堆芯的辐照效应。本文这种计算方法可以为核反应堆系统燃料组件的设计优化、热工水力安全分析提供参考。

发明内容

为了解决上述问题，本发明提供了一种考虑辐照效应的核反应堆弥散型板型燃料元件堆芯流固耦合计算方法，该方法将流体动力学分析程序、固体有限元分析程序(经过二次开发，编写考虑辐照效应的子程序)通过数据接口程序耦合起来，完整而准确地模拟辐照效应下核反应堆弥散型板型燃料元件堆芯的热工水力现象及流致振动情况。

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

一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法，包括如下步骤：

步骤一：获取核反应堆弥散型板型燃料元件堆芯参数，根据所获堆芯参数对堆芯流体部分进行几何建模，得到流体域几何模型，并在几何模型中标识出耦合面；对流体域几何模型划分网格，得到流体域网格模型，利用流体动力学计算程序对流体域网格模型设置边界条件和初始条件，得到流体域计算模型；

步骤二：根据步骤一中所获得的堆芯参数，对堆芯固体部分进行几何建模，得到固体域几何模型，在几何模型中标识出耦合面，固体域几何模型中耦合面要和流体域几何模型耦合面位置坐标一致；对固体域几何模型划分网格并设置边界条件，确定其分析类型为显性动力学分析，得到固体域计算模型；

步骤三：在固体有限元分析程序中，材料的物性设定为常数，没有考虑到其他因素对材料物性造成的影响，此外，外界因素还可能引起材料的本征变形或体积变化；对于弥散型板型燃料元件的堆芯，真实燃耗过程中存在辐照效应，辐照效应会导致弥散型芯体的导热系数、热膨胀系数、弹性模量及泊松比发生变化，还会导致芯体的肿胀及密实化现象；因此需要考虑辐照效应对固体域计算模型的影响；弥散型核燃料由UO₂分布在锆基体中构成，通过用户材料子程序UMAT和定义模型本征体积变化子程序对固体有限元分析程序进行如下二次开发；

用户材料子程序UMAT用来定义材料的热本构行为，包括随温度变化的等效导热系数、等效热膨胀系数、弹性模量及泊松比等材料物性；

基于麦克斯韦二相复合材料的导热模型进行弥散型芯体的导热系数计算，关系式如下：

式中，k_eff——弥散型芯体的等效导热系数，k₁——UO₂颗粒辐照条件下的导热系数/W·K^-1·m^-1，k₂——锆合金基体辐照条件下的导热系数/W·K^-1·m^-1，V——UO₂颗粒的体积百分比；k₁和k₂也是随温度变化的物理量，k₁基于UO₂四因子导热关系式如下：

k₁＝k₀·FD·FP·FM·FR

式中：k₀——UO₂颗粒无辐照条件下的导热系数/W·K-1·m-1；k₁——UO₂颗粒辐照条件下的导热系数/W·K-1·m-1；T₁——UO₂颗粒的温度/K；Bu——反应堆燃耗深度/％atom；p——燃料芯块孔隙率；s——形状因子；k₂用下面关系式来描述：

k₂＝7.51+2.09×10^-2T₂-1.45×10^-5T₂ ²+7.67×10^-9T₂ ³

式中：T₂——锆合金基体的温度/K；

等效热膨胀系数基于线性法则提出的混合模型，热膨胀系数计算关系式如下：

α_c＝α_mV_m+α_pV_p

式中：V_m、V_p分别为UO₂颗粒和锆基体合金的体积百分比；α_c为弥散型芯体的等效热膨胀系数；α_m为锆合金基体的热膨胀系数，温度变化对锆合金基体热膨胀系数影响较小，取值为α_m＝5.58e^-6/K；α_p为UO₂颗粒的热膨胀系数；UO₂颗粒的热膨胀系数计算关系式如下：

α_p＝-3.0289×10^-4+8.2417×10^-6(T₁-273)+2.1481×10^-9(T₁-273)²

式中，T₁是UO₂颗粒的温度/K；

UO₂颗粒的弹性模量和泊松比是反映燃料芯体力学性能的重要参数，采用的计算关系式如下：

E₁＝2.26×10¹¹[1-1.131×10^-4(T₁-273.15)][1-2.62(1-D_E)

v＝0.316

其中，E₁是UO₂颗粒的弹性模量/Pa；v是UO₂颗粒的泊松比；T₁是UO₂颗粒的温度/K；D_E是燃料芯块理论密度百分比，取95％；锆合金基体的弹性模量和泊松比计算关系式如下：

E₂＝[9.9×10⁵-566.9×(T₂-273.15)]×9.8067×10⁴

v＝0.3303+8.376×10^-5(T₂-273.15)

E₂是锆合金基体的弹性模量/Pa；v是锆合金的泊松比；T₂是锆合金基体的温度/K；

二氧化铀在裂变过程中会产生固体裂变产物及气体裂变产物，这些产物会导致二氧化铀的肿胀行为，称为辐照肿胀现象；固体裂变产物在反应堆运行过程中，其体积肿胀率与燃耗基本呈线性关系；裂变反应产生的气体会从晶粒扩散到晶界中，气体会聚集成核，引起二氧化铀发生肿胀；在低燃耗下，二氧化铀的塑性流动引起芯体内部的空隙被填充，造成芯体的体积减小，这就是燃料芯块的密实化现象；

在定义模型本征体积变化子程序中定义下面材料的本征体积变化关系式，对应弥散型核燃料辐照肿胀和密实化导致的燃料本征体积变化；

辐照肿胀的计算关系式参考帕拉斯基夫提出的经验关系式：

其中，SW(Bu)为肿胀程度；Bu是反应堆燃耗/％FIMA，Bu为用户预设值，不同Bu值代表不同燃耗程度下的效应；β₁和β₂分别表示由于气体裂变产物和固体裂变产物引起的肿胀率；T₁是UO₂颗粒的温度/K；ΔV为UO₂颗粒体积变化，V₀为UO₂颗粒初始体积；密实化的计算关系式如下：

β₃＝-[0.51exp(-59.9Bu)+4.76×10^-2×exp(-10.07Bu)

其中，Bu是反应堆燃耗/％FIMA，与辐照关系式中定义相同，均为用户预设值，不同Bu值代表不同燃耗程度下的效应；β₃表示由于密实化引起的体积变化率；

通过在材料子程序UMAT和定义材料本征体积变化子程序中定义以上公式，实现固体有限元分析程序在计算分析固体域计算模型过程中，固体域计算模型的物性变化、本征体积变化和辐照效应下的弥散型核燃料物性变化、本征体积变化情况一致；

步骤四：流体动力学计算程序对步骤一中流体域计算模型进行稳态计算，计算得到流体域的流量分布、压力、温度及换热系数结果参数；

步骤五：中间数据交换接口程序调取步骤四中所获得的结果参数，作为耦合瞬态计算的T₀时刻初始边界条件，同时调用步骤一中的流体域计算模型和步骤二中的固体域计算模型；两模型调用之后分别在流体动力学计算程序中和经过步骤三二次开发后的固体有限元分析程序中求解计算；求解过程如图3所示，首先，二次开发后固体有限元分析程序对固体域计算模型开展迭代计算，计算得到T₁时刻固体域的温度分布、位移及应力应变，耦合面的温度和节点位移通过中间数据交换接口程序反馈给流体域计算模型，作为其T₀时刻计算的初始边界条件，流体动力学分析程序计算得到T₁时刻流体域流量大小、换热系数及温度和压力分布，耦合面的温度、压力及换热系数通过中间数据交换接口传递给固体域计算模型，作为其T₁时刻计算的初始边界条件，二次开发固体有限元分析程序再次计算得到T₂时刻固体域的温度分布、位移及应力应变，耦合面的温度和节点位移再反馈给流体域计算模型，作为其T₁时刻计算的初始边界条件，流体动力学分析程序计算得到T₂时刻流体域流量大小、换热系数及温度和压力分布，耦合面的温度、压力及换热系数再传递给固体域计算模型；

步骤六：通过步骤五流体域和固体域每个时间步内的数据传递与反馈，计算得到核反应堆弥散型板型燃料元件堆芯在辐照效应的影响下，流体域即堆芯冷却剂的瞬时流量大小及温度、压力分布，以及固体域即燃料板的瞬时温度、位移及应力应变分布情况，当固体有限元分析程序计算到T_n时刻时，n为预设时间值，将计算所得耦合面固体域温度和节点位移反馈给流体域模型，作为其T_n-1时刻的初始边界条件，流体动力学计算模型计算得到T_n时刻流体域流量分布、温度、压力和换热系数，计算停止。

本发明具有以下优点和有益效果：

1.该方法提供了一种核反应堆弥散型板型燃料元件堆芯流固耦合方案，该方法适用于现有的大多数流体动力学计算程序，如FLUENT、ANSYS、FLOWMASTER等，以及固体有限元分析程序ABAQUS、COMSOL等。

2.该方法基于双向流固耦合方法和相应的数值模拟手段，通过压力场、温度场、以及位移场之间的交互耦合，在模拟流体侧三维流动换热状态的同时，可以模拟固体侧的热量传输和位移变化过程，实现了辐照条件下弥散型板型燃料元件堆芯流固耦合分析。

3.考虑辐照效应对固体域计算模型的影响，使得模型更加符合实际核反应堆堆芯的运行情况，从而可以更准确地观察核反应堆在辐照效应下，弥散型板型燃料元件堆芯中复杂的三维热工水力现象。

4.该方法不改变原程序程序的建模特性，可根据各程序的特性对燃料元件流体和固体部分进行分别建模，简单方便。

5.该方法能够对核反应堆不同燃耗程度下堆芯的三维水里热工现象进行模拟，从而研究不同燃耗深度对堆芯的影响。

本发明提出的核反应堆弥散型板型燃料元件堆芯流固耦合计算方法适用于弥散型板型燃料元件，可以采取不同的结构参数，还可以通过修改边界条件和初始条件进行反应性引入事故、失流事故等瞬态事故工况的计算。

附图说明

图1为本发明方法的流程图。

图2为核反应堆弥散型板型燃料组件简化模型，其中图2a为简化模型主视图，图2b为简化模型侧视图。

图3为步骤五中耦合数据交换流程图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步详细描述：

本发明提供了一种如图1所示核反应堆板型燃料元件堆芯流固耦合计算方法，具体方法如下：

步骤一：获取核反应堆弥散型板型燃料元件堆芯参数，根据所获得的堆芯参数，对堆芯流体部分进行几何建模，得到流体域几何模型，并在几何模型中标识处耦合面；考虑到弥散型板型燃料元件结构的对称性及计算资源的利用率，整个模型简化为三个流道两块板结构，如图2中图2a和图2b所示，则流体域几何模型由三个冷却剂流道，加上入口段和出口段部分组成；将建立好的几何模型导入网格划分工具中进行网格划分，得到流体域网格模型；然后利用流体动力学分析程序读取网格文件并设置初始条件和边界条件，选择适配的湍流模型，得到流体域计算模型；

步骤二：根据步骤一中所获得的堆芯参数，对图2简化模型中对应的固体部分进行几何建模，得到固体域几何模型，芯块和包壳之间无缝连接，同样在几何模型中标识出耦合面，固体域几何模型中耦合面要和流体域几何模型耦合面位置坐标一致；对固体域几何模型划分网格并设置边界条件，确定其分析类型为显性动力学分析，得到固体域计算模型；

步骤三：在固体有限元分析程序中，材料的物性设定为常数，没有考虑到其他因素对材料物性造成的影响，此外，外界因素还可能引起材料的本征变形或体积变化；对于弥散型板型燃料元件的堆芯，真实燃耗过程中存在辐照效应，辐照效应会导致弥散型芯体的导热系数、热膨胀系数、弹性模量及泊松比发生变化，还会导致芯体的肿胀及密实化现象；因此需要考虑辐照效应对固体域计算模型的影响；弥散型核燃料由UO₂分布在锆基体中构成，通过用户材料子程序UMAT和定义模型本征体积变化子程序对固体有限元分析程序进行如下二次开发；

用户材料子程序UMAT用来定义材料的热本构行为，包括随温度变化的等效导热系数、等效热膨胀系数、弹性模量及泊松比等材料物性；

基于麦克斯韦二相复合材料的导热模型进行弥散型芯体的导热系数计算，关系式如下：

式中，k_eff——弥散型芯体的等效导热系数，k₁——UO₂颗粒辐照条件下的导热系数/W·K^-1·m^-1，k₂——锆合金基体辐照条件下的导热系数/W·K^-1·m^-1，V——UO₂颗粒的体积百分比；k₁和k₂也是随温度变化的物理量，k₁基于UO₂四因子导热关系式如下：

k₁＝k₀·FD·FP·FM·FR

式中：k₀——UO₂颗粒无辐照条件下的导热系数/W·K-1·m-1；k₁——UO₂颗粒辐照条件下的导热系数/W·K-1·m-1；T₁——UO₂颗粒的温度/K；Bu——反应堆燃耗深度/％atom；p——燃料芯块孔隙率；s——形状因子；k₂用下面关系式来描述：

k₂＝7.51+2.09×10^-2T₂-1.45×10^-5T₂ ²+7.67×10^-9T₂ ³

式中：T₂——锆合金基体的温度/K；

等效热膨胀系数基于线性法则提出的混合模型，热膨胀系数计算关系式如下：

α_c＝α_mV_m+α_pV_p

式中：V_m、V_p分别为UO₂颗粒和锆基体合金的体积百分比；α_c为弥散型芯体的等效热膨胀系数；α_m为锆合金基体的热膨胀系数，温度变化对锆合金基体热膨胀系数影响较小，取值为α_m＝5.58e^-6/K；α_p为UO₂颗粒的热膨胀系数；UO₂颗粒的热膨胀系数计算关系式如下：

α_p＝-3.0289×10^-4+8.2417×10^-6(T₁-273)+2.1481×10^-9(T₁-273)²

式中，T₁是UO₂颗粒的温度/K；

UO₂颗粒的弹性模量和泊松比是反映燃料芯体力学性能的重要参数，采用的计算关系式如下：

E₁＝2.26×10¹¹[1-1.131×10^-4(T₁-273.15)][1-2.62(1-D_E)

v＝0.316

其中，E₁是UO₂颗粒的弹性模量/Pa；v是UO₂颗粒的泊松比；T₁是UO₂颗粒的温度/K；D_E是燃料芯块理论密度百分比，取95％；锆合金基体的弹性模量和泊松比计算关系式如下：

E₂＝[9.9×10⁵-566.9×(T₂-273.15)]×9.8067×10⁴

v＝0.3303+8.376×10^-5(T₂-273.15)

E₂是锆合金基体的弹性模量/Pa；v是锆合金的泊松比；T₂是锆合金基体的温度/K；

二氧化铀在裂变过程中会产生固体裂变产物及气体裂变产物，这些产物会导致二氧化铀的肿胀行为，称为辐照肿胀现象；固体裂变产物在反应堆运行过程中，其体积肿胀率与燃耗基本呈线性关系；裂变反应产生的气体会从晶粒扩散到晶界中，气体会聚集成核，引起二氧化铀发生肿胀；在低燃耗下，二氧化铀的塑性流动引起芯体内部的空隙被填充，造成芯体的体积减小，这就是燃料芯块的密实化现象；

在定义模型本征体积变化子程序中定义下面材料的本征体积变化关系式，对应弥散型核燃料辐照肿胀和密实化导致的燃料本征体积变化；

辐照肿胀的计算关系式参考帕拉斯基夫提出的经验关系式：

其中，SW(Bu)为肿胀程度；Bu是反应堆燃耗/％FIMA，Bu为用户预设值，不同Bu值代表不同燃耗程度下的效应；β₁和β₂分别表示由于气体裂变产物和固体裂变产物引起的肿胀率；T₁是UO₂颗粒的温度/K；ΔV为UO₂颗粒体积变化，V₀为UO₂颗粒初始体积；密实化的计算关系式如下：

β₃＝-[0.51exp(-59.9Bu)+4.76×10^-2×exp(-10.07Bu)

其中，Bu是反应堆燃耗/％FIMA，与辐照关系式中定义相同，均为用户预设值，不同Bu值代表不同燃耗程度下的效应；β₃表示由于密实化引起的体积变化率；

通过在材料子程序UMAT和定义材料本征体积变化子程序中定义以上公式，实现固体有限元分析程序在计算分析固体域计算模型过程中，固体域计算模型的物性变化、本征体积变化和辐照效应下的弥散型核燃料物性变化、本征体积变化情况一致；

步骤四：流体动力学计算程序对步骤一中流体域计算模型进行稳态计算，计算得到流体域的流量分布、压力、温度及换热系数结果参数；

步骤五：中间数据交换接口程序调取步骤四中所获得的结果参数，作为耦合瞬态计算的T₀时刻初始边界条件，同时调用步骤一中的流体域计算模型和步骤二中的固体域计算模型；两模型调用之后分别在流体动力学计算程序中和经过步骤三二次开发后的固体有限元分析程序中求解计算；求解过程如图3所示，首先，二次开发后固体有限元分析程序对固体域计算模型开展迭代计算，计算得到T₁时刻固体域的温度分布、位移及应力应变，耦合面的温度和节点位移通过中间数据交换接口程序反馈给流体域计算模型，作为其T₀时刻计算的初始边界条件，流体动力学分析程序计算得到T₁时刻流体域流量大小、换热系数及温度和压力分布，耦合面的温度、压力及换热系数通过中间数据交换接口传递给固体域计算模型，作为其T₁时刻计算的初始边界条件，二次开发固体有限元分析程序再次计算得到T₂时刻固体域的温度分布、位移及应力应变，耦合面的温度和节点位移再反馈给流体域计算模型，作为其T₁时刻计算的初始边界条件，流体动力学分析程序计算得到T₂时刻流体域流量大小、换热系数及温度和压力分布，耦合面的温度、压力及换热系数再传递给固体域计算模型；

步骤六：通过步骤五流体域和固体域每个时间步内的数据传递与反馈，计算得到核反应堆弥散型板型燃料元件堆芯在辐照效应的影响下，流体域即堆芯冷却剂的瞬时流量大小及温度、压力分布，以及固体域即燃料板的瞬时温度、位移及应力应变分布情况，当固体有限元分析程序计算到T_n时刻时，n为预设时间值，将计算所得耦合面固体域温度和节点位移反馈给流体域模型，作为其T_n-1时刻的初始边界条件，流体动力学计算模型计算得到T_n时刻流体域流量分布、温度、压力和换热系数，计算停止。

以上内容是结合具体的优选实施方式对本发明所做的进一步详细说明，不能确定本发明的具体实施方式仅限于此，对本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单的推演或替换，都应当视为属于本发明由所提交的权利要求书确定专利保护范围。

Claims (1)

1.一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法，其特征在于：针对核反应堆弥散型板型燃料元件堆芯中的热工水力现象，由流体动力学计算程序、固体有限元分析程序和中间数据交换接口程序，通过流体和固体之间实时双向的数据传递，实现程序间的耦合；同时，在固体有限元分析程序中考虑弥散型板型燃料元件在燃耗过程中辐照效应引起的变物性、肿胀及密实化现象，最终实现考虑辐照影响的堆芯三维热工水力现象双向流固耦合计算；

该方法包括以下步骤：

步骤一：获取核反应堆弥散型板型燃料元件堆芯参数，根据所获堆芯参数对堆芯流体部分进行几何建模，得到流体域几何模型，并在几何模型中标识出耦合面；对流体域几何模型划分网格，得到流体域网格模型，利用流体动力学计算程序对流体域网格模型设置边界条件和初始条件，得到流体域计算模型；

步骤二：根据步骤一中所获得的堆芯参数，对堆芯固体部分进行几何建模，得到固体域几何模型，在几何模型中标识出耦合面，固体域几何模型中耦合面要和流体域几何模型耦合面位置坐标一致；对固体域几何模型划分网格并设置边界条件，确定其分析类型为显性动力学分析，得到固体域计算模型；

步骤三：在固体有限元分析程序中，材料的物性设定为常数，没有考虑到其他因素对材料物性造成的影响，此外，外界因素还可能引起材料的本征变形或体积变化；对于弥散型板型燃料元件的堆芯，真实燃耗过程中存在辐照效应，辐照效应会导致弥散型芯体的导热系数、热膨胀系数、弹性模量及泊松比发生变化，还会导致芯体的肿胀及密实化现象；因此需要考虑辐照效应对固体域计算模型的影响；弥散型核燃料由UO₂分布在锆基体中构成，通过用户材料子程序UMAT和定义模型本征体积变化子程序对固体有限元分析程序进行如下二次开发；

用户材料子程序UMAT用来定义材料的热本构行为，包括随温度变化的等效导热系数、等效热膨胀系数、弹性模量及泊松比等材料物性；

基于麦克斯韦二相复合材料的导热模型进行弥散型芯体的导热系数计算，关系式如下：

式中，k_eff——弥散型芯体的等效导热系数，k₁——UO₂颗粒辐照条件下的导热系数/W·K^-1·m^-1，k₂——锆合金基体辐照条件下的导热系数/W·K^-1·m^-1，V——UO₂颗粒的体积百分比；k₁和k₂也是随温度变化的物理量，k₁基于UO₂四因子导热关系式如下：

k₁＝k₀·FD·FP·FM·FR

式中：k₀——UO₂颗粒无辐照条件下的导热系数/W·K-1·m-1；k₁——UO₂颗粒辐照条件下的导热系数/W·K-1·m-1；T₁——UO₂颗粒的温度/K；Bu——反应堆燃耗深度/％atom；p——燃料芯块孔隙率；s——形状因子；k₂用下面关系式来描述：

k₂＝7.51+2.09×10^-2T₂-1.45×10^-5T₂ ²+7.67×10^-9T₂ ³

式中：T₂——锆合金基体的温度/K；

等效热膨胀系数基于线性法则提出的混合模型，热膨胀系数计算关系式如下：

α_c＝α_mV_m+α_pV_p

式中：V_m、V_p分别为UO₂颗粒和锆基体合金的体积百分比；α_c为弥散型芯体的等效热膨胀系数；α_m为锆合金基体的热膨胀系数，温度变化对锆合金基体热膨胀系数影响较小，取值为α_m＝5.58e^-6/K；α_p为UO₂颗粒的热膨胀系数；UO₂颗粒的热膨胀系数计算关系式如下：

α_p＝-3.0289×10^-4+8.2417×10^-6(T₁-273)+2.1481×10^-9(T₁-273)²

式中，T₁是UO₂颗粒的温度/K；

UO₂颗粒的弹性模量和泊松比是反映燃料芯体力学性能的重要参数，采用的计算关系式如下：

E₁＝2.26×10¹¹[1-1.131×10^-4(T₁-273.15)][1-2.62(1-D_E)

v＝0.316

其中，E₁是UO₂颗粒的弹性模量/Pa；v是UO₂颗粒的泊松比；T₁是UO₂颗粒的温度/K；D_E是燃料芯块理论密度百分比，取95％；锆合金基体的弹性模量和泊松比计算关系式如下：

E₂＝[9.9×10⁵-566.9×(T₂-273.15)]×9.8067×10⁴

v＝0.3303+8.376×10^-5(T₂-273.15)

E₂是锆合金基体的弹性模量/Pa；v是锆合金的泊松比；T₂是锆合金基体的温度/K；

二氧化铀在裂变过程中会产生固体裂变产物及气体裂变产物，这些产物会导致二氧化铀的肿胀行为，称为辐照肿胀现象；固体裂变产物在反应堆运行过程中，其体积肿胀率与燃耗基本呈线性关系；裂变反应产生的气体会从晶粒扩散到晶界中，气体会聚集成核，引起二氧化铀发生肿胀；在低燃耗下，二氧化铀的塑性流动引起芯体内部的空隙被填充，造成芯体的体积减小，这就是燃料芯块的密实化现象；

在定义模型本征体积变化子程序中定义下面材料的本征体积变化关系式，对应弥散型核燃料辐照肿胀和密实化导致的燃料本征体积变化；

辐照肿胀的计算关系式参考帕拉斯基夫提出的经验关系式：

其中，SW(Bu)为肿胀程度；Bu是反应堆燃耗/％FIMA，Bu为用户预设值，不同Bu值代表不同燃耗程度下的效应；β₁和β₂分别表示由于气体裂变产物和固体裂变产物引起的肿胀率；T₁是UO₂颗粒的温度/K；ΔV为UO₂颗粒体积变化，V₀为UO₂颗粒初始体积；密实化的计算关系式如下：

β₃＝-[0.51exp(-59.9Bu)+4.76×10^-2×exp(-10.07Bu)

其中，Bu是反应堆燃耗/％FIMA，与辐照关系式中定义相同，均为用户预设值，不同Bu值代表不同燃耗程度下的效应；β₃表示由于密实化引起的体积变化率；

通过在材料子程序UMAT和定义材料本征体积变化子程序中定义以上公式，实现固体有限元分析程序在计算分析固体域计算模型过程中，固体域计算模型的物性变化、本征体积变化和辐照效应下的弥散型核燃料物性变化、本征体积变化情况一致；

步骤四：流体动力学计算程序对步骤一中流体域计算模型进行稳态计算，计算得到流体域的流量分布、压力、温度及换热系数结果参数；

步骤五：中间数据交换接口程序调取步骤四中所获得的结果参数，作为耦合瞬态计算的T₀时刻初始边界条件，同时调用步骤一中的流体域计算模型和步骤二中的固体域计算模型；两模型调用之后分别在流体动力学计算程序中和经过步骤三二次开发后的固体有限元分析程序中求解计算；求解过程如图3所示，首先，二次开发后固体有限元分析程序对固体域计算模型开展迭代计算，计算得到T₁时刻固体域的温度分布、位移及应力应变，耦合面的温度和节点位移通过中间数据交换接口程序反馈给流体域计算模型，作为其T₀时刻计算的初始边界条件，流体动力学分析程序计算得到T₁时刻流体域流量大小、换热系数及温度和压力分布，耦合面的温度、压力及换热系数通过中间数据交换接口传递给固体域计算模型，作为其T₁时刻计算的初始边界条件，二次开发固体有限元分析程序再次计算得到T₂时刻固体域的温度分布、位移及应力应变，耦合面的温度和节点位移再反馈给流体域计算模型，作为其T₁时刻计算的初始边界条件，流体动力学分析程序计算得到T₂时刻流体域流量大小、换热系数及温度和压力分布，耦合面的温度、压力及换热系数再传递给固体域计算模型；

步骤六：通过步骤五流体域和固体域每个时间步内的数据传递与反馈，计算得到核反应堆弥散型板型燃料元件堆芯在辐照效应的影响下，流体域即堆芯冷却剂的瞬时流量大小及温度、压力分布，以及固体域即燃料板的瞬时温度、位移及应力应变分布情况，当固体有限元分析程序计算到T_n时刻时，n为预设时间值，将计算所得耦合面固体域温度和节点位移反馈给流体域模型，作为其T_n-1时刻的初始边界条件，流体动力学计算模型计算得到T_n时刻流体域流量分布、温度、压力和换热系数，计算停止。