弥散核燃料等效辐照肿胀计算模拟

蔡 维，赵云妹，龚 辛，丁淑蓉，霍永忠

（复旦大学力学与工程科学系，上海 200433）

弥散核燃料等效辐照肿胀计算模拟

蔡 维，赵云妹，龚 辛，丁淑蓉＊，霍永忠

（复旦大学力学与工程科学系，上海 200433）

本文将弥散核燃料芯体看作一种特殊的颗粒复合材料，利用细观计算力学的方法，假设燃料颗粒在芯体中周期性分布，建立了对芯体等效辐照肿胀进行计算模拟的有限元模型。考虑颗粒的辐照肿胀和基体材料的辐照硬化效应，分别建立了燃料颗粒和基体材料的应力更新算法，编制了用户材料子程序，在Abaqus软件中实现了芯体等效辐照肿胀的有限元模拟。计算分析了颗粒大小和体积含量对芯体等效辐照肿胀的影响，并得到了等效辐照肿胀的拟合公式。研究结果表明，影响芯体等效辐照肿胀的主要因素是颗粒的辐照肿胀和体积含量。

弥散核燃料；有限元法；均匀化方法；等效辐照肿胀

相对于传统的棒状核燃料元件，弥散核燃料元件［1］具有高热导，可达较高的燃耗，已广泛应用于研究试验堆中，在先进核反应堆和核废料处理方面也具有良好的应用前景［1－2］。

弥散核燃料元件由弥散核燃料芯体和包壳组成，分棒状和板型两种形式，是组成弥散核燃料组件的基本构件。弥散核燃料芯体在结构上类似于颗粒复合材料，由一定体积含量的燃料颗粒弥散分布于金属基体中构成。为达到一定的功率密度，芯体中需装载大量的燃料颗粒。

为对弥散核燃料组件进行优化设计，需研究其在堆内的辐照－热力耦合行为。辐照试验研究耗时、花费极高且无法在线观测，建立元件或组件堆内辐照热力耦合行为的理论模型和数值模拟手段，对其堆内行为的演化进行预测成为一种重要的辅助手段［3－5］。对弥散核燃料元件或组件的堆内热力耦合行为进行数值模拟，需将弥散核燃料芯体均匀化［6］，等效成一种均匀的介质，不再考虑颗粒和基体的相互作用。

文献［7－8］运用均匀化理论，通过计算模拟对弥散核燃料的等效热传导率和弹性性能进行了研究，而对其等效辐照肿胀的计算模拟研究尚未见报道。在核反应堆苛刻的环境中，核燃料颗粒所产生的裂变产物导致其发生辐照肿胀，并随燃耗增长而增大［9］，成为影响燃料元件或组件堆内安全的一重要的载荷。

本文假设燃料颗粒在芯体中周期性分布，根据燃料芯体的微观结构，选取代表性体元（RVE）和有限元模型，考虑金属基体的辐照硬化效应，建立计算模拟燃料芯体等效辐照肿胀的数值方法。同时，计算分析颗粒体积含量和大小对芯体等效辐照肿胀的影响，并对解析公式进行拟合。

1 有限元模型

1．1 有限元几何模型

考虑无限大的弥散核燃料，如图1a所示，假设燃料颗粒在基体中呈简单立方分布。此弥散核燃料芯体可看作图1b所示的代表性体元在空间周期性排列而成。由于代表性体元包含了燃料芯体所有的微观结构信息，可用代表性体元代替整个燃料芯体，对其等效辐照肿胀性能进行研究。在代表性体元中，燃料颗粒随燃耗的发展产生辐照肿胀，这导致体元随燃耗的发展而产生相对的体积变化，称为等效辐照肿胀。体元位于无限大复合燃料芯体中，根据对称性，随燃耗的发展体元边界保持为平面。考虑到体元结构和载荷的对称性，选取其1／8作为有限元模型进行研究，如图1c所示。将其进行空间离散，如图1d所示。

为计算分析燃料颗粒体积含量变化对芯体等效肿胀的影响，本文分别选取颗粒体积含量为10．657 4%、15．514 0%、19．081 6%、23．832 4%、30．300 9%的几种情况进行建模。假设颗粒的半径r为0．05mm，有限元几何模型的边长列于表1。对于颗粒尺寸不同的燃料芯体，可类似地建立有限元模型。

1．2 边界条件

根据前面的分析可知，弥散核燃料芯体的等效辐照肿胀需通过计算有限元模型的体积变化得到。有限元模型中的主要载荷来自于颗粒的辐照肿胀，其使得颗粒和基体之间产生强烈的力学相互作用，导致颗粒和基体内部产生大应变、有限元模型的棱边产生变化。针对图1c所施加的位移边界条件如下：

图1 燃料颗粒在芯体中的分布形式（a）、代表性体元（b）、有限元几何模型（c）及有限元网格（d）Fig．1 Distribution pattern of fuel particles in fuel meat（a），RVE（b），finite element geometric model（c）and finite element grid（d）

表1 颗粒体积含量与有限元几何模型边长的关系Table 1 Relationship between particle volume fraction and size of finite element geometric model

1）分别对x＝L、y＝L和z＝L的这3个面，施加节点位移约束方程，使这3个面始终保持为平面；

2）对x＝0、y＝0和z＝0的这3个面，根据结构和载荷对称性，施加对称边界条件。

1．3 材料参数

1）燃料颗粒的材料参数［10－12］

（1）弹性模量和泊松比

式中：E为弹性模量，Pa；T为温度，K；D为颗粒的相对理论密度，取95%。

泊松比ν取0．316。

（2）塑性模型

（3）肿胀经验公式

燃料颗粒的体积肿胀SWfuel包括裂变气体肿胀（gs）、裂变固体肿胀（ss）和裂变收缩（ds）3部分。

式中，V0为燃料颗粒的原始体积。裂变气体肿胀、裂变固体肿胀导致的每产生1020cm－3的裂变密度下的体积变化分数［13－14］分别为：

裂变密实所导致的体积收缩［15］为：

式中：BU为燃耗，MW·d／kgUO2；ADST为调整因子，取0．6。

2）金属基体的材料参数

受快中子的冲击，基体材料将产生辐照硬化效应，下面给出考虑辐照硬化效应的基体材料性能，与文献［11］中所用模型相同。

（1）弹性模量和泊松比［16］

式中：φ为快中子注量率，m－2·s－1；t为时间，s。

（2）辐照硬化塑性模型［17］

式中：σ为真应力，Pa；ε为真应变；˙ε为真应变率，若˙ε＜10－5s－1，设˙ε＝10－5s－1；参数K、n和m分别为强度系数、应变硬化指数和应变率敏感指数。本文中快中子注量率φ取2× 1018m－2·s－1，T取300～730K。

2 数值模拟方法

为对弥散核燃料芯体的等效辐照肿胀进行计算，需将所考察的燃耗区间分为多个时间步。针对前面给出的特殊的材料性能，考虑大变形效应，需在旋转坐标系下［12］分别建立燃料颗粒和基体材料的三维应力更新算法。

2．1 燃料颗粒的应力更新算法

对于一典型的时间步［t，t＋Δt］，t时刻的柯西应力σtij和弹性对数应变εiej，t满足如下关系：

式中，λ和G为拉梅常数，由材料的弹性模量和泊松比计算得到。则t＋Δt时刻的柯西应力σt＋Δtij为：

式中：Δσij为柯西应力增量；Δεeij为弹性对数应变的增量。对于燃料颗粒内部的每个积分点，有：

其中：Δεtotalij为总的应变增量；Δεswij和Δεpij分别为辐照肿胀和塑性应变增量。

先假设无塑性应变产生，即Δεpij＝0，计算得到试验应力σprij，并由此得到Mises等效应力，其中，sprij为试验应力的偏斜部分，sprij＝σprij－σprkkδij／3。若为根据式（2）计算所得的t时刻的等效塑性应变），则说明此增量步有新的塑性应变产生，否则，试验应力σprij为t＋Δt时刻的真实应力σt＋Δtij，应力更新完成。

对于有新塑性应变产生的情况，由塑性流动法则［12］：

对式（17）进行后欧拉积分，得到：

则

2．2 金属基体的应力更新算法

考虑辐照硬化效应，基体材料的弹性参数和塑性模型随时间而变化。对于时间步［t，t＋Δt］，积分点的弹性应变增量为：

t＋Δt时刻的柯西应力σt＋Δtij为：

式中：G（t＋Δt）和λ（t＋Δt）分别为积分点t＋Δt时刻的拉梅常数；ΔG和Δλ为因时间变化所导致的拉梅常数的变化。

针对2．1和2．2节的应力更新算法，分别编制了定义材料力学本构关系的子程序，并引入ABAQUS软件的有限元计算中，实现弥散核燃料等效辐照肿胀的计算模拟。

3 数值模拟方法有效性的验证

针对体积含量为10．657 4%、温度为600K的弥散核燃料，计算其随燃耗增长的微观的位移、应力和应变场，并根据有限元计算的结果来验证所建立的数值模拟方法的有效性。

3．1 辐照硬化塑性的验证

3．2 颗粒辐照肿胀的验证

图2 基体中的节点A（a）及节点A计算结果的验证（b）Fig．2 Node A in matrix（a）and validation of computation results of node A（b）

图3 燃料颗粒体积变化的有限元计算结果与辐照肿胀理论值的比较Fig．3 Comparison of volumetric variation result obtained by FEM with theoretical irradiation swelling result

一方面，通过有限元计算可得到燃料颗粒边缘处节点在不同燃耗下的位移ΔR，进而可计算燃料颗粒在不同燃耗下的实际体积变化ΔV／V0＝（（R0＋ΔR）3－R30）／R30，其中，R0为燃料颗粒的原始半径。另一方面，根据式（3）可计算出燃料颗粒不同燃耗下的理论辐照肿胀结果（需对不同的燃耗单位进行换算，参见文献［10］）。通过上述两种途径得到的结果示于图3。由图3a可明显看出，在很低的燃耗下，燃料颗粒的体积变化主要是辐照所致裂变收缩引起的，实际体积变化为负，在FIMA为0．3%左右出现了最低点。且在很低的燃耗下，有限元计算所得燃料颗粒的实际体积缩小量（绝对值）低于辐照肿胀的理论计算结果（绝对值）。由于燃料颗粒的实际体积变化包括辐照肿胀和弹性变化两部分。弹性变形所导致的体积变化主要与燃料颗粒内各点的静水压力－σkk／3有关（图4）。

由图4a所示的FIMA为0．3%的静水压力云图可看出，燃料颗粒内各点的静水压力均为负，即静水拉力有使燃料颗粒体积增大的贡献，抵消了一部分辐照所致裂变收缩的贡献，所以，燃料颗粒在很低燃耗下的实际体积缩小量低于理论的辐照肿胀结果。与图3a不同，由图3b可看出，在较高燃耗下，燃料颗粒的实际体积变化略低于理论辐照肿胀结果。同样，从图4b所示的FIMA为30%的静水压力云图看，FIMA为30%的燃料颗粒各点具有很大的静水压力，其驱使燃料颗粒的体积缩小，抵消了一部分燃料颗粒的辐照肿胀，所以，有限元计算所得实际体积变化低于辐照肿胀理论结果是合理的。通过上述分析，可看出燃料颗粒的辐照肿胀得到了正确的模拟。

图4 不同燃耗下的静水压力云图Fig．4 Contour plots of hydrostatic pressure with different burnups

4 结果与讨论

本节针对温度为600K的弥散核燃料，计算分析颗粒体积含量和大小对其等效辐照肿胀的影响，并得到等效辐照肿胀的拟合公式。

4．1 网格收敛性考察

为保证有限元计算结果的精度和计算的经济性，对不同颗粒体积含量情况下的有限元模型，由粗到细逐步加密网格，并考察结果的收敛性，力求得到一种经济高效的有限元网格。网格的逐步加密图如图5所示。针对颗粒体积含量为20%的情况，给出了FIMA为30%下有限元模型Path 1上的位移和Mises应力随网格加密而变化的情况，如图6所示。从图6可看出，位移很易收敛；而对于Mises应力，随着网格的进一步加密才得到收敛。考虑到经济性，可采用具有3 024个单元的网格进行计算。对于其他体积含量的有限元模型，也采用类似的方法获得经济性和计算精度兼具的网格。

图5 网格逐步加密图Fig．5 Refinement of finite element mesh

图6 不同网格下Path 1上各点的位移和Mises应力Fig．6 Displacements and Mises stresses on Path 1for different cases of mesh grid

4．2 颗粒体积含量对芯体等效辐照肿胀的影响

对不同体积含量的有限元模型进行计算，得到的芯体的等效辐照肿胀如图7所示。从图7可看到，芯体的等效辐照肿胀随颗粒体积含量的增加而增加。在较低燃耗下，高体积含量的情况体积收缩也较大；对于不同体积含量的情况，由收缩转为膨胀的转变燃耗几乎相同，与燃料颗粒的转变燃耗相近；转变燃耗后，等效辐照肿胀与燃耗呈线性关系；对于颗粒体积含量大的芯体，其等效辐照肿胀的变化率也较高。由于受到金属基体的约束，在FIMA为30%以下时，等效辐照肿胀均小于6%。

4．3 颗粒大小对芯体等效辐照肿胀的影响

对于体积含量分别为10．657 4%和30．300 9%两种情况，分别针对0．02、0．05、0．08mm的颗粒建立有限元模型，对芯体的等效肿胀进行计算。计算的不同颗粒大小下的等效辐照肿胀如图8所示。可看出颗粒大小对芯体等效肿胀几乎无影响。在一定燃耗下，在不同颗粒大小的有限元模型中，颗粒和基体之间的力学相互作用存在差异，而这种差异对芯体等效辐照肿胀的总体贡献几乎为零，故造成图4所示的结果。

图7 不同体积含量下的等效辐照肿胀Fig．7 Equivalent irradiation swelling for different particle volume fractions

图8 颗粒大小对等效辐照肿胀的影响Fig．8 Effect of particle size on equivalent irradiation swelling

4．4 芯体等效辐照肿胀的拟合公式

通过以上计算分析可得出，芯体的等效辐照肿胀主要决定于颗粒的辐照肿胀和体积含量以及燃料颗粒和基体中弹性变形的贡献，而弹性变形的贡献相对较小。对于所建立的有限元模型，将其体积的绝对变化首先近似为：

其中：ΔVfuel为核燃料颗粒部分的绝对体积变化；Vfuel0为核燃料颗粒部分的原始体积。

据此，可将芯体的等效辐照肿胀初步近似为：

其中，Vf为颗粒的原始体积含量。SW为有限元计算结果，Vf×SWfuel为初步近似值。根据有限元计算所得的位移结果可得到不同燃耗下的实际棱边长度a＋Δa（a为原始的棱边长度，Δa为棱边的绝对变形量），进而可得到SW＝（（a＋Δa）3－a3）／a3。同时，由式（3）所得燃料颗粒的SWfuel可计算得到不同燃耗下的Vf×SWfuel。

将芯体等效辐照肿胀的有限元计算值与通过式（25）得到的初步近似结果进行比较，结果如图9所示。从图9可看出，二者在较低燃耗下仍存在较大误差，还需进一步拟合。

为达到更精确的近似值，进行分段考虑。

当FIMA≤3%时，有：

当FIMA＞3%时，有：

有限元计算结果与利用式（26）、（27）得到的拟合结果的比较示于图10。可看出，利用式（26）、（27）可得到更好的拟合结果，可将式（26）、（27）作为弥散核燃料等效辐照肿胀的拟合公式。

图9 芯体等效肿胀有限元计算结果与近似值的比较Fig．9 Comparison of finite element calculation result with approximate result

图10 芯体等效肿胀有限元计算结果与分段考虑近似值的比较Fig．10 Comparison of finite element calculation result with segment approximate result

5 结论

本研究建立了对弥散核燃料等效辐照肿胀进行计算模拟的有限元模型，并考虑基体材料的辐照硬化效应，建立了对等效辐照肿胀进行计算模拟的数值方法，分别考察了燃料颗粒大小、体积含量对芯体等效辐照肿胀的影响，得到了等效辐照肿胀的拟合公式。在所考察的模型和参数范围内，得到的主要结论如下：

1）芯体等效辐照肿胀主要受颗粒辐照肿胀和体积含量的影响。在较低燃耗下，燃料颗粒的裂变收缩起主要作用，等效辐照肿胀为负值；在较高燃耗下，等效辐照肿胀主要来自于颗粒固体裂变产物的贡献，与燃耗呈线性关系。

2）本文所建立的有限元模型和数值模拟方法是有效的。

［1］ XU Z．Design strategies for optimizing high burnup fuel in pressurized water reactors［D］．US：Massachu－setts Institute of Technology，2003．

［2］ DUYN L．Evaluation of the mechanical behavior of a metal－matrix dispersion fuel for plutonium burning［D］．US：Georgia Institute of Technology，2003．

［3］ WANG Q，YAN X，DING S，et al．Research on the interfacial behaviors of plate－type dispersion nuclear fuel elements［J］．Journal of Nuclear Materials，2010，399（1）：41－54．

［4］ WANG Q，YAN X，DING S，et al．Mechanical behaviors of the dispersion nuclear fuel plates induced by fuel particle swelling and thermal effectⅠ：Effects of variations of the fuel particle volume fractions［J］．Journal of Nuclear Materials，2010，400（2）：157－174．

［5］ DING S，WANG Q，HUO Y．Mechanical behaviors of the dispersion nuclear fuel plates induced by fuel particle swelling and thermal effectⅡ：Effects of variations of the fuel particle diameters［J］．Journal of Nuclear Materials，2010，397（1－3）：80－91．

［6］ 乔生儒．复合材料细观力学性能［M］．西安：西北工业大学出版社，1997．

［7］ 姜馨，丁淑蓉，霍永忠．弥散型燃料等效弹性性质的有限元模拟［J］．核动力工程，2011，32（3）：43－47．

JIANG Xin，DING Shurong，HUO Yongzhong．Finite element simulation for equivalent elastic properties of dispersion fuel elements［J］．Nuclear Power Engineering，2011，32（3）：43－47（in Chinese）．

［8］ 姜馨，丁淑蓉，霍永忠．弥散型燃料元件等效热传导系数的有限元模拟［J］．核动力工程，2010，31（3）：45－49．

JIANG Xin，DING Shurong，HUO Yongzhong．Finite element simulation for effective thermal conductivity coefficient of dispersion fuel element［J］．Nuclear Power Engineering，2010，31（3）：45－49（in Chinese）．

［9］ MATPRO－09：A handbook of materials properties for use in the analysis of light water reactor fuel rod behavior［M］．US：NRC，1976．

［10］SUZUKI M，SAITOU H．Light water reactor fuel analysis code FEMAⅪ－6［R］．Japan：Atomic Energy Agency，2005．

［11］GONG X，DING S，ZHAO Y，et al．Effects of irradiation hardening and creep on the thermomechanical behaviors in inert matrix fuel elements［J］．Mechanics of Materials，2013，65：110－123．

［12］BELYTSCHKO T，LIU W K，MORAN B．Nonlinear finite elements for continua and structures［M］．UK：John Wiley ＆Sons Ltd．，2000．

［13］CHUBB W，STORHOK V W，KELLER D L．Factors affecting the swelling of nuclear fuel at high temperatures［J］．Nuclear Technology，1973，18：231－255．

［14］NAKAJIMA T，ICHIKAWA M，IWANO Y，et al．FEMAⅪ－Ⅲ：A computer code for the analysis of thermal and mechanical behavior of fuel rods［R］．Japan：JAERI，1985．

［15］WIESENACK W，VANKEERGHEN M，THANKAPPAN R．Assessment of UO2conductivity degradation based on in－pile temperature data［R］．Norway：HWR，1996．

［16］FISHER E F，RENKEN C J．Single－crystal elastic moduli and the hcp－transformation in Ti，Zr and Hf［J］．Physical Review A，1954，135（2）：482－494．

［17］HAGRMAN D L，REYMAN G A．MATPRO－11：A handbook of materials properties for use in the analysis of light water reactor fuel rod behavior［M］．US：NRC，1979．

Calculation Simulation of Equivalent Irradiation Swelling for Dispersion Nuclear Fuel

CAI Wei，ZHAO Yun－mei，GONG Xin，DING Shu－rong＊，HUO Yong－zhong
（Department of Mechanics and Engineering Science，Fudan University，Shanghai 200433，China）

The dispersion nuclear fuel was regarded as a kind of special particle composites．Assuming that the fuel particles are periodically distributed in the dispersion nuclear fuel meat，the finite element model to calculate its equivalent irradiation swelling was developed with the method of computational micro－mechanics．Considering irradiation swelling in the fuel particles and the irradiation hardening effect in the metal matrix，the stress update algorithms were established respectively for the fuel particles and metal matrix．The corresponding user subroutines were programmed，and the finite element simulation of equivalent irradiation swelling for the fuel meat was performed in Abaqus．The effects of the particle size and volume fraction on the equivalent irradiation swelling were investigated，and the fitting formula of equivalent irradiation swelling was obtained．The results indicate that the main factors to influence equivalent irradiation swelling of the fuel meat are the irradiation swelling and volume fraction of fuel particles．

dispersion nuclear fuel；finite element method；homogenization method；equivalent irradiation swelling

TL341

：A

：1000－6931（2015）03－0502－09

10．7538／yzk．2015．49．03．0502

2013－12－20；

2014－03－08

国家自然科学基金资助项目（11172068，91226101，11272092）；教育部博士点基金资助项目（20110071110013）

蔡 维（1987—），男，安徽六安人，硕士研究生，从事弥散燃料的等效辐照性能研究

＊通信作者：丁淑蓉，E－mail：dsr1971＠163．com