A PHASE FIELD METHOD FOR IRRADIATION-THERMAL-MECHANICAL COUPLING FRACTURE OF THE DISPERSION NUCLEAR FUEL
-
摘要: 弥散型核燃料作为第4代核电技术的重要组成部分, 以其反应均匀、温度梯度小及高燃耗等显著优势, 在核电领域展现出重要应用前景. 随着核电站设计服役年限的延长, 为防止核裂变产物的泄露, 弥散型核燃料元件断裂和失效行为的预测更加重要. 断裂相场法是近年发展起来的计算断裂力学方法, 并在多物理场作用下的复杂介质断裂行为预测上取得了成功. 首先基于连续介质热力学建立了辐照-热-力耦合断裂相场方法用以预测弥散型核燃料在辐照、热应力和机械载荷共同作用下的断裂及传热行为. 随后, 分别对压水堆坏境下的弥散型核燃料代表性单元及整板的断裂行为进行了数值仿真计算, 获得了弥散型核燃料内部的温度场、裂纹相场和静水压力场. 结果显示: 均匀颗粒分布的弥散型核燃料具有较小的温度梯度; 基体中未发现裂纹, 结构损伤形式以燃料颗粒断裂为主, 裂纹在燃料颗粒边缘萌生并向内部扩展. 本研究可为弥散型核燃料元件的断裂行为预测提供有效的计算模拟方法和数值分析依据.Abstract: Dispersion nuclear fuels have become a pivotal component in fourth-generation nuclear power technologies, and exhibit promising applications in the nuclear energy region due to their uniform fission reaction, small temperature gradients, and high burnup capabilities. With the increase of the design service life of nuclear fuel elements, a heightened demand for predicting the fracture and failure behaviors of the fuel element is proposed to avoid the leakage of nuclear fission products. The phase field fracture method is a recently developed computational fracture mechanics approach, and has achieved considerable success in predicting the fracture behaviors of complex solid media even under multi-physics conditions. In this work, we firstly propose an irradiation-thermal-mechanical coupling phase field model of dispersion nuclear fuels based on continuum thermodynamics in order to predict the fracture and heat transfer behaviors of the nuclear fuels under irradiation, thermal stress, and mechanical loadings. Subsequently, numerical simulations for the representative volume element and the entire plate of dispersion nuclear fuel in the pressurized water reactor environment are conducted. And the temperature field, the crack phase field and the hydrostatic stress field inside the dispersion nuclear fuel are obtained. The results reveal that dispersion nuclear fuels with uniform particle distribution exhibit relatively small temperature gradients. Notably, no phase field cracks are observed in the nuclear fuel matrix, and the damage of the dispersion nuclear fuel primarily manifests as the fracture of fuel particles. Many cracks nucleate at the edges of fuel particles, and then propagate inward the particle center. This work can provide an effective simulation method and numerical analysis basis for predicting the fracture behavior of dispersion nuclear fuel element.
-
引 言
随着人类社会的快速发展, 对能源的需求日益增加, 核能作为一种高效、绿色的能源形式, 被广泛认为是未来能源发展的重要方向之一. 弥散型核燃料[1-3]通过将裂变材料颗粒(如铀和钚的化合物)均匀分散在基体材料(如金属、陶瓷或石墨)中, 形成了类似于颗粒增强复合材料的结构, 如图1所示[4]. 弥散型核燃料相比于传统核燃料棒具有良好的热传导性能和耐辐照性能. 由于燃料颗粒被基体材料有效隔离, 裂变产物对周围环境的损伤被大大降低. 这种结构使得弥散型核燃料芯块在应对核事故时具有更高的安全性. 此外, 弥散型核燃料还具有良好的可设计性, 通过调整燃料颗粒和基体材料的种类、比例以及分布方式, 可以优化燃料的性能以满足不同核反应堆的实际需求. 除核能发电领域以外, 在航天领域, 弥散型核燃料的高能量密度和长寿命等特点使其成为未来深空探测任务中一种理想的能源供应方式[5].
随着核反应堆运行时间的延长, 由于弥散型燃料的颗粒和基体材料在物理、化学和力学性能上存在差异, 它们在辐照、温度和压力等环境因素的作用下可能会产生不同的变形和损伤. 这些损伤可能引发裂纹的萌生和扩展, 最终导致燃料元件的失效[6]. 此外, 对核燃料的耐事故工况性能进行评估也颇具现实意义. 对弥散型核燃料的断裂行为进行深入研究, 揭示其断裂机理和影响因素, 对于提高燃料元件的安全性和可靠性具有重要意义.
目前, 国内外专家学者对弥散型核燃料的等效弹性性质、辐照肿胀、蠕变和热传导等方面的性能进行了较为广泛的研究[7-12]. 然而, 在断裂力学方面的研究相对较少[12], 尚待进一步深入.
断裂相场法是在连续损伤模型基础上建立的一种断裂力学研究方法, 它源自Francfort等[13]提出的断裂变分理论, 并随后由Bourdin等[14]改进为正则化裂纹的断裂理论, 用弥散损伤区域表征裂纹, 并引入标量d来描述材料的损伤程度(0代表材料完好, 1代表材料完全开裂). 相场法无需预制裂纹, 无需额外断裂准则, 可以直接模拟材料从完好到裂纹萌生、裂纹扩展和裂纹分岔/汇合等复杂的断裂过程[15-16].
断裂相场的求解在于由Griffith断裂能与应变能组成的总能量泛函最小化. 若合理地设计总能量泛函可以将多物理场因素引入到断裂相场模型中. 目前, 断裂相场法已成功应用于热-力耦合[17-18]、力-化学耦合[19-20]以及流固耦合[21-22]等多物理场情景. 在核燃料断裂力学领域, 断裂相场法也已用于UO2细观尺度的裂变气泡导致的断裂行为[23-24]、核级石墨的断裂行为[25]、核燃料棒的断裂力学仿真[26]以及TRISO颗粒的辐照收缩断裂[27].
本文旨在建立弥散型核燃料的辐照-热-力耦合断裂相场模型, 用以模拟裂纹在燃料颗粒和基体材料中的萌生、扩展过程, 从而为更好地理解弥散型核燃料的断裂行为提供数值计算方案, 为优化燃料设计、提高反应堆性能提供科学依据.
1. 数学模型
本文首先推导适用于弥散型核燃料断裂力学仿真的数值模型. 针对弥散型核燃料的辐照-热-力耦合行为, 本文建立力场、热场与相场耦合的计算断裂力学模型, 并在随后的材料参数模型中考虑辐照效应.
1.1 变量声明
在介绍控制方程前, 首先介绍偏微分方程的自变量. 控制方程包括: 热传导方程、力学平衡方程和相场演化方程. 本文的控制方程中的自变量为位移u、裂纹相场d和温度T.
1.2 基本方程
准静态平衡方程
$$ \nabla \cdot {{\boldsymbol{\sigma}} } + {{\boldsymbol{b}}} = {\boldsymbol{0}} $$ (1) 式中, ${{\boldsymbol{\sigma}} }$为Cauchy应力张量, ${{\boldsymbol{b}}}$为体力.
Gurtin认为材料相变是微观尺度力作用的结果, 而微观的力满足与常规力学量相似的平衡方程, 即对于两相材料的相场演化满足微力平衡方程[28]
$$ \nabla \cdot {{\boldsymbol{\xi}} } + h = 0 $$ (2) 式中, ${{\boldsymbol{\xi}} }$表示微应力, 为一阶张量; $h$表示体微力, 为标量.
考虑力学、相场与传热的耦合, 系统的总内能变化率可表示为
$$ \dot e = {\boldsymbol{\sigma }}:{\boldsymbol{\dot \varepsilon }} - h \cdot \dot d + {\boldsymbol{\xi }} \cdot \nabla \dot d + \dot q - \nabla \cdot {\boldsymbol{Q}} $$ (3) 式中, $e$指单位体积的内能; $\dot e$为内能对时间的一阶导数, 即内能变化率; $ {\boldsymbol{\varepsilon}} = (\nabla {\boldsymbol{u}} + {\boldsymbol{u}} \nabla) / 2 $为应变; $\dot q$为单位体积的内部热源率; ${\boldsymbol{Q}}$指热通量.
引入Legendre变换, 即
$$ e = \psi + T\eta $$ (4) 式中, $\psi $指Helmholtz自由能密度, $\eta $指单位体积的熵. 将式(4)代入式(3)可得
$$ \dot \psi + \eta \dot T + T\dot \eta = {\boldsymbol{\sigma }}:{\boldsymbol{\dot \varepsilon }} - h \cdot \dot d + {\boldsymbol{\xi }} \cdot \nabla \dot d + \dot q - \nabla \cdot {\boldsymbol{Q}} $$ (5) 进一步地, 根据热力学第二定律, 单位体积的局部耗散${\mathcal{D}}$可以表示为
$$ \begin{split} & {\mathcal{D}} = T\dot \eta + \nabla \cdot {\boldsymbol{Q}} - \dot q = \\ &\qquad {\boldsymbol{\sigma }}:{\boldsymbol{\dot \varepsilon }} - h \cdot \dot d + {\boldsymbol{\xi }} \cdot \nabla \dot d - \eta \dot T - \dot \psi \geqslant 0\end{split} $$ (6) 1.3 本构方程
假设系统的Helmholtz自由能密度函数依赖于独立变量${\boldsymbol{\varepsilon }}$, $d$, $\nabla d$和$T$, 即
$$ \psi = \psi \left( {{\boldsymbol{\varepsilon }},d,\nabla d,T} \right) $$ (7) Helmholtz自由能的变化率满足
$$ \dot \psi = \frac{{\partial \psi }}{{\partial {\boldsymbol{\varepsilon }}}}:{\boldsymbol{\dot \varepsilon }} + \frac{{\partial \psi }}{{\partial d}} \cdot \dot d + \frac{{\partial \psi }}{{\partial \nabla d}} \cdot \nabla \dot d + \frac{{\partial \psi }}{{\partial T}} \cdot \dot T $$ (8) 将式(8)代入式(6)可得
$$ \begin{split} & \left( {{\boldsymbol{\sigma }} - \frac{{\partial \psi }}{{\partial {\boldsymbol{\varepsilon }}}}} \right):{\boldsymbol{\dot \varepsilon }} - \left( {h + \frac{{\partial \psi }}{{\partial d}}} \right) \cdot \dot d + \\ &\qquad \left( {{\boldsymbol{\xi }} - \frac{{\partial \psi }}{{\partial \nabla d}}} \right) \cdot \nabla \dot d - \left( {\eta + \frac{{\partial \psi }}{{\partial T}}} \right) \dot T \geqslant 0\end{split} $$ (9) 根据Coleman-Noll过程[29], 可得如下本构方程
$$ {\boldsymbol{\sigma }} = \frac{{\partial \psi }}{{\partial {\boldsymbol{\varepsilon }}}},\;{\text{ }}h = - \frac{{\partial \psi }}{{\partial d}},\;{\text{ }}{\boldsymbol{\xi }} = \frac{{\partial \psi }}{{\partial \nabla d}},\;{\text{ }}\eta = - \frac{{\partial \psi }}{{\partial T}} $$ (10) 1.4 Helmholtz自由能密度
对于热-力-断裂相场问题, 系统的Helmholtz自由能密度可分为3部分, 即
$$ \psi = {\psi _e} + {\psi _c} + {\psi _{{\mathrm{th}}}} $$ (11) 式中, ${\psi _e}$为弹性能部分, ${\psi _c}$为断裂能部分, ${\psi _{{\mathrm{th}}}}$为热能部分.
为考虑相场裂纹处弹性能的退化, 通常在应变能部分引入能量退化函数$g\left( d \right)$. 若考虑未损伤材料满足各向同性的线弹性本构模型, 则弹性应变能部分可表示为
$$ {\psi _e} = g\left( d \right)\left[ {\lambda /2{{\left( {{\text{tr}}{{\boldsymbol{\varepsilon }}_e}} \right)}^2} + \mu {\text{tr}}\left( {{\boldsymbol{\varepsilon }}_e^2} \right)} \right] $$ (12) 式中, ${{\boldsymbol{\varepsilon }}_e}$为弹性应变, $\lambda $和$\mu $为Lamé系数. 根据Miehe等[30]的工作, 退化函数仅应该作用于拉伸部分应变能, 即
$$ {\psi _e} = g\left( d \right)\psi _e^ + + \psi _e^ - $$ (13) 式中, $\psi _e^ + $表示拉伸部分应变能密度, $\psi _e^ - $表示压缩部分应变能密度, 具体表达式为
$$ \psi _e^ \pm = \frac{\lambda }{2}\left\langle {{\varepsilon _1} + {\varepsilon _2} + {\varepsilon _3}} \right\rangle _ \pm ^2 + \mu \left( {\left\langle {{\varepsilon _1}} \right\rangle _ \pm ^2 + \left\langle {{\varepsilon _2}} \right\rangle _ \pm ^2 + \left\langle {{\varepsilon _3}} \right\rangle _ \pm ^2} \right) $$ (14) 式中, ${\varepsilon _i}{\text{ (}}i = 1,2,3{\text{)}}$为弹性应变张量${{\boldsymbol{\varepsilon }}_e}$的特征值. ${\left\langle {} \right\rangle _ \pm }$为一种算符, 满足${\left\langle x \right\rangle _ \pm } = (x \pm \left| x \right|)/2$.
对于断裂能部分的自由能密度, 采用相场近似Griffith断裂能[14, 31], 即
$$ {\psi _c} = {G_c}\gamma \left( {d,\nabla d} \right) $$ (15) 式中, ${G_c}$为临界断裂释放率; $\gamma $为裂纹密度函数, 其表达式为
$$ \gamma \left( {d,\nabla d} \right) =\frac{1}{c_0} \left[\frac{{\alpha \left( d \right)}}{{{l_0}}} + {l_0}{\left| {\nabla d} \right|^2}\right] $$ (16) 式中, $\alpha \left( d \right)$为相场裂纹的几何函数, ${c_0} = \displaystyle\int_0^d {\alpha \left( d \right){\text{d}}d} $为与几何函数相关的常数, ${l_0}$为裂纹弥散特征尺度.
热能部分的自由能密度的表达式为[17]
$$ {\psi _{{\mathrm{th}}}} = \rho c\left[ {T - {T_0} - T\ln \left( {T/{T_0}} \right)} \right] $$ (17) 式中, $\rho $为密度, $c$为比热容, ${T_0}$表示初始的温度. 将式(13)、式(15)和式(17)代入式(11)得系统总的Helmholtz自由能密度, 即
$$ \begin{split} & \psi = g\left( d \right)\psi _e^ + \left( {{{\boldsymbol{\varepsilon }}_e}} \right) + \psi _e^ - \left( {{{\boldsymbol{\varepsilon }}_e}} \right) + \\ &\qquad \frac{G_c}{c_0}\left[ {\frac{{\alpha \left( d \right)}}{{{l_0}}} + {l_0}{{\left| {\nabla d} \right|}^2}} \right] + \\ &\qquad {\rho _0}c\left[ {T - {T_0} - T\ln \left( {T/{T_0}} \right)} \right] \end{split} $$ (18) 1.5 控制方程
将式(18)和式(14)代入本构方程(10)可得
$$\begin{split} & {\boldsymbol{\sigma }} = g\left( d \right)\sum\limits_{i = 1}^3 {\left[ {\lambda {{\left\langle {{\varepsilon _1} + {\varepsilon _2} + {\varepsilon _3}} \right\rangle }_ + } + 2\mu {{\left\langle {{\varepsilon _i}} \right\rangle }_ + }} \right]{{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i}} + \\ &\qquad \sum\limits_{i = 1}^3 {\left[ {\lambda {{\left\langle {{\varepsilon _1} + {\varepsilon _2} + {\varepsilon _3}} \right\rangle }_ - } + 2\mu {{\left\langle {{\varepsilon _i}} \right\rangle }_ - }} \right]{{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i}} \end{split} $$ (19) $$ h = - g'\left( d \right)\psi _e^ + - \frac{{{G_c}}}{{{c_0}{l_0}}}\alpha '\left( d \right) $$ (20) $$ {\boldsymbol{\xi }} = 2{G_c}{l_0}\nabla d $$ (21) $$ \eta = \rho c\ln (T/{T_0}) $$ (22) 式中, ${{\boldsymbol{n}}_i}{\text{ (}}i = 1,2,3{\text{)}}$表示特征向量. 将式(20)和式(21)代入式(2)可得断裂相场演化方程, 即
$$ g'\left( d \right)\psi _e^ + + \frac{{{G_c}}}{{{c_0}{l_0}}}\alpha '\left( d \right) - \frac{2{G_c}{l_0}}{c_0}{\nabla ^2}d = 0 $$ (23) 将式(8)和式(10)代入式(5), 化简后可得余项为
$$ T\dot \eta = \dot q - \nabla \cdot {\boldsymbol{Q}} $$ (24) 将式(22)代入式(24)可得
$$ c\rho \dot T = \dot q - \nabla \cdot {\boldsymbol{Q}} $$ (25) 考虑各向同性的Fourier热传导定律$ {\boldsymbol{Q}} = - k\nabla T $, 可得热传导方程
$$ c\rho \dot T = k{\nabla ^2}T + \dot q $$ (26) 式中, k为热传导系数.
1.6 小结
综上所述, 本问题的待求方程如下:
应力平衡方程
$$ \nabla \cdot {{\boldsymbol{\sigma}} } + {{\boldsymbol{b}}} = {\boldsymbol{0}} $$ (27) 相场演化方程
$$ g'\left( d \right)\psi _e^ + + \frac{{{G_c}}}{{{c_0}{l_0}}}\alpha '\left( d \right) - \frac{{{2{G_c}{l_0}}}}{{{c_0}}} {\nabla ^2}d = 0 $$ (28) 热传导方程
$$ c\rho \dot T = k{\nabla ^2}T + \dot q $$ (29) 力学本构方程
$$ \begin{split} & {\boldsymbol{\sigma }} = g\left( d \right)\sum\limits_{i = 1}^3 {\left[ {\lambda {{\left\langle {{\varepsilon _1} + {\varepsilon _2} + {\varepsilon _3}} \right\rangle }_ + } + 2\mu {{\left\langle {{\varepsilon _i}} \right\rangle }_ + }} \right]{{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i}} + \\ &\qquad \sum\limits_{i = 1}^3 {\left[ {\lambda {{\left\langle {{\varepsilon _1} + {\varepsilon _2} + {\varepsilon _3}} \right\rangle }_ - } + 2\mu {{\left\langle {{\varepsilon _i}} \right\rangle }_ - }} \right]{{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i}}\end{split} $$ (30) 上述方程的初边值条件不再赘述.
本文采用经典AT1断裂相场模型[32], 其几何函数与能量退化函数的具体形式为
$$\qquad\qquad\qquad \alpha \left( d \right) = d $$ (31) $$\qquad\qquad\qquad g\left( d \right) = {\left( {1 - d} \right)^2} $$ (32) 相比于经典的AT2模型($\alpha = {d^2}$), AT1模型的优势为具有线弹性阶段, 从而可以明显地区分裂纹成核行为[33].
2. 材料参数模型
本节分别介绍UO2颗粒与锆合金包壳的关键材料参数模型. 由于在辐照及高温作用下, 两种材料参数模型大多与温度或辐照效应相关. 温度相关的材料属性以及热膨胀变形会使得弹性能与温度相关. 根据式(10)所得熵的表达式会出现弹性能相关的贡献, 这会使得式(22)不再准确. 但这种弹性能耦合所释放的热量远远小于核燃料裂变释放的热量, 因此本文的模型忽略了这种效应, 这使得经典的热传导方程(26)和式(29)可以近似成立, 类似的近似处理可见文献[34-35].
2.1 UO2的材料参数模型
UO2的热传导系数采用Lucuta模型[36], 该模型综合考虑了温度、燃耗、孔隙率以及裂变产物对热传导系数的影响, 具体形式为
$$ {k_{{{\mathrm{UO}}_2}}} = {k_0} \cdot {{FD}} \cdot {{FP}} \cdot {{FM}} \cdot {{FR}} $$ (33) 式中, ${k_0}$表示未经辐照的致密UO2的热传导系数, 其表达式为
$$ {k_0} = \frac{1}{{0.037\;5 + 2.165 \times {{10}^{ - 4}}T}} + \frac{{4.715 \times {{10}^9}}}{{{T^2}}}{{{\mathrm{e}}} ^{ - \frac{{16\;361}}{T}}} $$ (34) FD表示溶解的裂变产物的影响, 其表达式为
$$ \begin{split} & {{FD}} = (1.09{B^{ - 3.265}} + 0.064\;3\sqrt {{T \mathord{\left/ {\vphantom {T B}} \right. } B}} ) \cdot \\ &\qquad {{\mathrm{tan}}^{ - 1}}\left(\frac{1}{{1.09{B^{ - 3.265}} + 0.064\;3\sqrt {{T \mathord{\left/ {\vphantom {T B}} \right. } B}} }}\right) \end{split} $$ (35) 式中, B表示燃耗. FP表示析出的裂变产物的影响, 其表达式为
$$ {{FP}} = {\text{1}} + \frac{{0.019B}}{{3 - 0.019B}} \cdot \frac{1}{{1 + {{{\mathrm{e}}} ^{ - {{(T - 1200)} \mathord{\left/ {\vphantom {{(T - 1200)} {100}}} \right. } {100}}}}}} $$ (36) FM表示孔隙率的影响, 其表达式为
$$ {{FM}} = \frac{{1 - P}}{{1 + (s - 1)P}} $$ (37) 式中, P表示孔隙率; s为孔隙的形状系数, 通常对于球形孔隙取1.5. FR表示辐照效应的影响, 其表达式为
$$ {{FR}} = 1 - \frac{{0.2}}{{1 + {{{\mathrm{e}}} ^{{{(T - 900)} \mathord{\left/ {\vphantom {{(T - 900)} {80}}} \right. } {80}}}}}} $$ (38) UO2颗粒的杨氏模量E为[37]
$$ E = 2.26 \times {10^{11}}\left( {1 - 1.131 \times {{10}^{ - 4}}T} \right)\left[ {1 - 2.62\left( {1 - D} \right)} \right] $$ (39) 式中, 温度T的单位为°C; D为理论密度, 通常在92% ~ 98%之间. UO2的泊松比为$v = 0.316$. Lamé系数与杨氏模量和泊松比之间满足如下换算关系
$$ \lambda = \frac{{Ev}}{{\left( {1 + v} \right)\left( {1 - 2v} \right)}},\qquad \mu = \frac{E}{{2\left( {1 + v} \right)}} $$ (40) $$ {c_{{{\mathrm{UO}}_2}}} = \frac{{{K_1}{\theta ^2}{{{\mathrm{e}}} ^{{\theta \mathord{\left/ {\vphantom {\theta T}} \right. } T}}}}}{{{T^2}{{({{{\mathrm{e}}} ^{{\theta \mathord{\left/ {\vphantom {\theta T}} \right. } T}}} - 1)}^2}}} + {K_2}T + \frac{{Y{K_3}{E_D}}}{{2R{T^2}}}{{{\mathrm{e}}} ^{ - \frac{{{E_D}}}{{RT}}}} $$ (41) 式中, 温度T的单位为K, Y为氧与金属的比率, R为通用气体常数, 其余系数满足${K_1} = 296.7{{\;{\mathrm{J}}/({\mathrm{kg}}\cdot {\mathrm{K}})}}$, ${K_2} = 2.43 \times {10^{ - 2}}{{\;{\mathrm{J}} \cdot {\mathrm{kg}}\cdot}}{{\text{K}}^2}$, ${K_3} = 8.745 \times {10^7}{\text{ J/kg}}$, $\theta = 535.285{\text{ K}},$ ${E_D} = 1.577 \times {10^5}{\text{ J/mol}}.$
UO2颗粒的体积热生成速率与裂变反应速率呈线性关系[9], 即
$$ \dot q = {c_1} \cdot \dot f $$ (42) 式中, ${c_1} = 3.204 \times {10^{ - 11}}{\text{ J/fission}}$为常数, $\dot f$为裂变事件的反应速率.
考虑热膨胀的作用, 总应变${\boldsymbol{\varepsilon }}$可按加性分解为弹性应变${{\boldsymbol{\varepsilon }}_e}$与热膨胀应变${{\boldsymbol{\varepsilon }}_{{\mathrm{th}}}}$, 即
$$ {\boldsymbol{\varepsilon }} = {{\boldsymbol{\varepsilon }}_e} + {{\boldsymbol{\varepsilon }}_{{\mathrm{th}}}} $$ (43) 假设热膨胀应变与温度变化呈线性关系, 则
$$ {{\boldsymbol{\varepsilon}} _{{\mathrm{th}}}} = {\alpha _1}\left( {T - {T_0}} \right){\boldsymbol{I}} $$ (44) 式中, ${\alpha _1}$为线膨胀系数, ${\boldsymbol{I}}$为单位二阶张量. 本文取UO2颗粒的热膨胀系数为10.471 K−1 [40].
计算中UO2颗粒的密度取$\rho = 10.431{\text{ g/c}}{{\text{m}}^3}$, 临界断裂释放率取${G_c} = 3{\text{ J/}}{{\text{m}}^2}$[26].
2.2 锆合金的材料参数模型
锆合金的热传导系数取MATPRO模型[41], 即
$$\begin{split} & k = 7.511 + 2.088 \times {10^{ - 2}}T - \\ &\qquad 1.450 \times {10^{ - 5}}{T^2} + 7.668 \times {10^{ - 9}}{T^3} \end{split} $$ (45) 锆合金的杨氏模量[42]
$$ E = 9.9 \times {10^5} - 566.9\left( {T - 273.15} \right) \times 9.806\;7 \times {10^4} $$ (46) 锆合金的泊松比
$$ v = 0.330\;3 + 8.376 \times {10^{ - 5}}\left( {T - 273.15} \right) $$ (47) 式(45) ~ 式(47)中, 温度T的单位为K.
锆合金的比热容取MATPRO模型[41], 该参数与温度的关系满足表1中数据的线性插值.
计算中所需锆合金的其他参数有热膨胀系数${\alpha _1} = 5.5 \times {10^{ - 6}}{\text{ }}{{\text{K}}^{ - 1}}$[41], 密度$\rho = 6.56{\text{ g/c}}{{\text{m}}^3}$, 临界断裂释放率${G_c} = 1000{\text{ J/}}{{\text{m}}^2}$[26].
2.3 裂纹处的传热参数
由于UO2在发生裂变反应后会生成大量惰性气体, 如氙气(Xe), 导致裂缝处发生热流阻塞现象, 因此需要建立与相场相关的裂纹处的传热参数模型[26]. 本文假设裂缝的传热参数随相变趋于氙气的传热参数, 即
$$ k = \left( {1 - d} \right){k_{{{\mathrm{UO}}_2}}} + d{k_{{\mathrm{Xe}}}} $$ (48) $$ c = \left( {1 - d} \right){c_{{{\mathrm{UO}}_2}}} + d{c_{{\mathrm{Xe}}}} $$ (49) 式中, ${k_{{\mathrm{Xe}}}}$为氙气的热传导系数, 其表达式为[43]
$$ {k_{{\mathrm{Xe}}}} = 4.351 \times {10^{ - 5}}{T^{0.861\;6}} $$ (50) 式中, 温度T的单位为K. $ {c_{{\mathrm{Xe}}}}{\text{ }}\left( {{{{\mathrm{J}}\cdot{\mathrm{mo}}}}{{\text{l}}^{ - 1}}\cdot{{\text{K}}^{{{ - 1}}}}} \right) $为氙气的比热容, 其表达式为[44]
$$ \begin{split} & {c_{{\mathrm{Xe}}}} = 20.786 + 7.449\;32 \times {10^{ - 7}}\bar T - 2.049\;40 \times {10^{ - 7}}{{\bar T}^2} + \\ &\qquad 1.066\;66 \times {10^{ - 8}}{{\bar T}^3} + {{2.500\;26 \times {{10}^{ - 8}}} /{{{\bar T}^2}}} \\[-5pt]\end{split} $$ (51) 式中, $\bar T = T/1000$的单位为K, 氙气的摩尔质量为131.293 g/mol.
3. 数值算例
本文的数值算例使用并行有限元多物理场求解方案计算[45]. 特别地, 对于断裂相场采用变分不等式求解器进行求解以约束裂纹的不可逆性[46-47].
3.1 代表性单元
取弥散型核燃料的代表性单元及其边界条件如图2(a)所示, 其中燃料颗粒的直径为$100{\text{ μm}}$, 代表单元宽度为$250{\text{ μm}}$. 在压水堆坏境中, 核燃料由于冷却剂作用受到对流传热边界条件与压力共同作用[26]. 在代表性单元中, 对其上边界施加模拟压水堆环境的边界条件. 对流传热边界的表达式为
$$ - k\frac{{\partial T}}{{\partial n}} = {k_{{\mathrm{conv}}}}\left( {T - {T_f}} \right) $$ (52) 式中, ${k_{{\mathrm{conv}}}}$为对流热传导系数, 本文取值为$3.4 \;\times {10^4}{\text{ W}} \cdot {{\text{m}}^{ - 2}}\cdot{{\mathrm{K}}^{ - 1}}$; ${T_f} = 600{\text{ }}{\mathrm{K}}$为外围流体的温度. 外围流体压强取15.5 MPa.
本文将相场裂纹的弥散特征尺度视为模型参数, 取值为${l_0} = 1.0 \times {10^{ - 6}}{\text{ m}}$. 该取值远小于结构的尺寸(0.4%的代表性单元宽度), 从而可以确保相场裂纹良好地近似Griffith裂纹[30, 48]. 断裂相场法对断裂区域的网格密度有较高的要求以满足问题的收敛性. 由于弥散型核燃料的颗粒容易发生断裂, 颗粒处设置了更高密度的网格. 单元尺寸通常取0.25 ~ 0.5倍的${l_0}$即可满足问题的收敛性[49-50]. 本文颗粒处的单元尺寸取1/2倍的弥散裂纹特征尺度, 即$5.0 \times {10^{ - 7}}{\text{ m}}$. 图2 (b)为网格划分情况, 共划分了
56640 个线性Lagrange四边形网格.代表性单元的初始温度设置为273 K, 从$t = 0{\text{ s}}$燃料颗粒开始发生裂变反应, 裂变速率逐渐升高, 体积热释放速率也随之升高, 图3给出了UO2颗粒的裂变速率随时间的演化. 对于不同裂变增速, 当裂变速率达到约$7.76 \times {10^{21}}{\text{ fissions}} \cdot {{\text{m}}^{ - 3}} \cdot {{\text{s}}^{ - 1}}$时均发生裂纹成核及演化. 图4给出了相应的中心温度随时间的演化结果. 随着裂变速率加快, 燃料中心温度逐渐升高, 但温升速度趋于放缓.
图5给出了仿真得到代表性单元最终的温度云图、相场云图和静水压力云图. 可以发现, 在3种裂变速率增长(如图3所示)的情况下获得断裂后的云图基本一致, 这意味着弥散型燃料的断裂行为与裂变速率增长速度呈弱相关性. 此外, 代表性单元的断裂形式主要表现为燃料颗粒断裂, 燃料颗粒的裂纹由边缘萌生并向内部扩展, 与文献报道一致[51]. 在本文目前的边界条件及参数选取下, 暂未发现弥散型核燃料的基体裂纹, 表明锆合金基底的弥散型核燃料具有较高的安全性.
图6给出了沿温度提取路径(见图2)所示的中线温度曲线, 可以看到弥散型核燃料的中心与外侧的温度差约为500 K, 该数值远小于棒状核燃料的结果[26], 在燃料颗粒处存在局部温升. 此外, 在3种裂变速率增长情况下, 得到的中线温度曲线基本重合, 再次印证了裂变速率增速对最终材料内的温度、应力和相场分布影响较小.
3.2 弥散型燃料板
由于代表性单元仅能考虑单侧冷却水的作用, 而实际的核燃料各边界都可受到冷却水作用. 本文进一步对弥散型核燃料的“整板”进行仿真, 出于计算量的考虑, 选取了含45颗均匀分布燃料颗粒的燃料板进行计算, 如图7所示. 同时对整板的4个边界施加压强与对流传热边界条件. 计算模型采用了
475200 个线性Lagrange四边形网格以确保相场模拟的准确性.图8给出了发生断裂后的温度云图, 燃料板的中心温度最高并向四周呈梯度下降, 在四周冷却剂的作用下弥散型核燃料板与代表性单元相比存在更高的温度梯度. 图9(a)则给出了相应的裂纹相场云图, 可以看到高温度梯度产生了大量裂纹. 裂纹形式与代表性单元的计算结果保持一致, 大量裂纹仍在颗粒边缘成核, 并向颗粒内部扩展. 不同位置的燃料颗粒发生不同程度和形式的断裂, 中心处颗粒相比边缘处的颗粒产生了数量更多的裂纹. 图9(b)给出了颗粒断裂后的典型静水压力云图, 可以看到燃料颗粒内部存在静水压力梯度, 导致裂纹在颗粒边缘成核并向颗粒内部扩展. 不同位置处的颗粒因承受的温度梯度不同而具有不同的断裂形貌.
4. 结 论
本文提出了考虑辐照-热-力耦合作用的弥散型核燃料的断裂相场模型, 并分别对代表性单元和整板进行了仿真计算, 所得结论如下.
(1) UO2燃料颗粒的断裂形式以边缘裂纹为主, 裂纹在颗粒边缘成核并向内部扩展.
(2)裂变速率增速对弥散型核燃料的最终断裂及传热行为影响较小, 获取的温度、断裂相场与静水压力云图基本一致.
本文所建立的分析计算模型暂未考虑裂变气体压力、裂变碎片损伤区以及颗粒的非均匀分布等因素的影响, 在后续研究中将进一步完善以获取更为可靠的弥散型核燃料数值仿真评价方法. 本研究可为断裂相场法在复杂工程中的应用提供参考案例.
-
-
[1] Savchenko A, Konovalov I, Vatulin A, et al. Dispersion type zirconium matrix fuels fabricated by capillary impregnation method. Journal of Nuclear Materials, 2007, 362(2): 356-363
[2] Savchenko AM, Vatulin AV, Morozov AV, et al. Inert matrix fuel in dispersion type fuel elements. Journal of Nuclear Materials, 2006, 352(1): 372-377
[3] 舒瀚东, 周何, 马千驰等. U3O8/ZrN核壳结构微球制备及其在弥散型陶瓷核燃料中的应用研究. 中国陶瓷. 2024, 60(2): 45-52 (Shu Handong, Zhou He, Ma Qianchi, et al. Fabrication and application in dispersed ceramic nuclear fuel of U3O8/ZrN core-shell structure microspheres. China Ceramics, 2024, 60(2): 45-52 (in Chinese) Shu Handong, Zhou He, Ma Qianchi, et al. Fabrication and application in dispersed ceramic nuclear fuel of U3O8/ZrN core-shell structure microspheres. China Ceramics, 2024, 60(2): 45-52 (in Chinese)
[4] 严晓青. 弥散型核燃料板辐照力学行为的数值模拟. [硕士论文]. 上海: 复旦大学, 2009: 2-3 (Yan Xiaoqing. Numerical simulation researches on the irradiation-induced mechanical behaviors of dispersion nuclear plate. [Master Thesis]. Shanghai: Fudan University, 2009: 2-3 (in Chinese) Yan Xiaoqing. Numerical simulation researches on the irradiation-induced mechanical behaviors of dispersion nuclear plate. [Master Thesis]. Shanghai: Fudan University, 2009: 2-3 (in Chinese)
[5] 房玉良, 刘林, 孙海亮等. 核热推进反应堆燃料元件发展概述. 宇航总体技术, 2020, 4(1): 63-70 (Fang Yuliang, Liu Lin, Sun Hailiang, et al. Development of fuel elements in nuclear thermal propulsion system, Astronautical Systems Engineering Technology, 2020, 4(1): 63-70 (in Chinese) Fang Yuliang, Liu Lin, Sun Hailiang, et al. Development of fuel elements in nuclear thermal propulsion system, Astronautical Systems Engineering Technology, 2020, 4(1): 63-70 (in Chinese)
[6] 董颖璇, 吕俊男, 李群. 颗粒团聚行为对弥散型核燃料芯体失效的影响分析. 原子能科学技术, 2024, 58(4): 868-877 (Dong Yingxuan, Lyu Junnan, Li Qun. Analysis of particle agglomeration effect on failure of dispersion nuclear fuel meat. Atomic Energy Science and Technology, 2024, 58(4): 868-877 (in Chinese) Dong Yingxuan, Lyu Junnan, Li Qun. Analysis of particle agglomeration effect on failure of dispersion nuclear fuel meat. Atomic Energy Science and Technology, 2024, 58(4): 868-877 (in Chinese)
[7] 丁淑蓉, 龚辛, 赵云妹等. 弥散核燃料燃烧演化过程中的关键力学问题. 力学季刊, 2018, 39(1): 1-21 (Ding Shurong, Gong Xin, Zhao Yunmei, et al. Key mechanical problems in dispersion nuclear fuels during their burning evolution process. Chinese Quarterly of Mechanics, 2018, 39(1): 1-21 (in Chinese) Ding Shurong, Gong Xin, Zhao Yunmei, et al. Key mechanical problems in dispersion nuclear fuels during their burning evolution process. Chinese Quarterly of Mechanics, 2018, 39(1): 1-21 (in Chinese)
[8] Zhang J, Wang H, Wei H, et al. Modelling of effective irradiation swelling for inert matrix fuels. Nuclear Engineering and Technology, 2021, 53(8): 2616-2628 doi: 10.1016/j.net.2021.02.019
[9] Ding S, Huo Y, Yan X. Modeling of the heat transfer performance of plate-type dispersion nuclear fuel elements. Journal of Nuclear Materials, 2009, 392(3): 498-504 doi: 10.1016/j.jnucmat.2009.04.015
[10] Zhang J, Zhang J, Wang H, et al. Modeling of mesoscale creep behaviors and macroscale creep responses of composite fuels under irradiation conditions. Acta Mechanica Solida Sinica, 2022, 35(6): 1040-1054
[11] Zhao Y, Gong X, Ding S, et al. A numerical method for simulating the non-homogeneous irradiation effects in full-sized dispersion nuclear fuel plates. International Journal of Mechanical Sciences, 2014, 81: 174-183 doi: 10.1016/j.ijmecsci.2014.02.012
[12] Kadambi SB, Aagesen LK, Zhang Y, et al. Assessment of effective elastic constants of U-10Mo fuel: A multiscale modeling and homogenization study. Journal of Nuclear Materials, 2024, 599: 155225 doi: 10.1016/j.jnucmat.2024.155225
[13] Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, 1998, 46(8): 1319-1342 doi: 10.1016/S0022-5096(98)00034-9
[14] Bourdin B, Francfort GA, Marigo JJ. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids, 2000, 48(4): 797-826 doi: 10.1016/S0022-5096(99)00028-9
[15] 吴建营. 固体结构损伤破坏统一相场理论、算法和应用. 力学学报, 2021, 53(2): 301-329 (Wu Jianying. On the unified phase-field theory for damage and failure in solids and structures: Theoretical and numerical aspects. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 301-329 (in Chinese) doi: 10.6052/0459-1879-20-295 Wu Jianying. On the unified phase-field theory for damage and failure in solids and structures: Theoretical and numerical aspects. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 301-329 (in Chinese) doi: 10.6052/0459-1879-20-295
[16] Feng Y, Freddi F, Li J, et al. Phase-field model for 2D cohesive-frictional shear fracture: An energetic formulation. Journal of the Mechanics and Physics of Solids, 2024, 189: 105687 doi: 10.1016/j.jmps.2024.105687
[17] Miehe C, Schänzel LM, Ulmer H. Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Computer Methods in Applied Mechanics and Engineering, 2015, 294: 449-485
[18] Du C, Cui H, Zhang H. Thermal fatigue behaviors of thin-walled structures with holes: Experiments and phase field fracture modeling. International Journal of Fatigue, 2024, 185: 108338 doi: 10.1016/j.ijfatigue.2024.108338
[19] Cui C, Ma R, Martínez-Pañeda E. A phase field formulation for dissolution-driven stress corrosion cracking. Journal of the Mechanics and Physics of Solids, 2021, 147: 104254
[20] Guo Y, Na S. A reactive-transport phase-field modelling approach of chemo-assisted cracking in saturated sandstone. Computer Methods in Applied Mechanics and Engineering, 2024, 419: 116645 doi: 10.1016/j.cma.2023.116645
[21] Li Y, Peng G, Tang J, et al. Thermo-hydro-mechanical coupling simulation for fracture propagation in CO2 fracturing based on phase-field model. Energy, 2023, 284: 128629 doi: 10.1016/j.energy.2023.128629
[22] Liu SF, Wang W, Jia Y, et al. Modeling of hydro-mechanical coupled fracture propagation in quasi-brittle rocks using a variational phase-field method. Rock Mechanics and Rock Engineering, 2024, https://doi.org/10.1007/s00603-024-03896-5
[23] Zhang S, Jiang W, Gamble KA, et al. Comparing the impact of thermal stresses and bubble pressure on intergranular fracture in UO2 using 2D phase field fracture simulations. Journal of Nuclear Materials, 2023, 574: 154158 doi: 10.1016/j.jnucmat.2022.154158
[24] Jiang W, Hu T, Aagesen LK, et al. A phase-field model of quasi-brittle fracture for pressurized cracks: Application to UO2 high-burnup microstructure fragmentation. Theoretical and Applied Fracture Mechanics, 2022, 119: 103348 doi: 10.1016/j.tafmec.2022.103348
[25] Chakraborty P, Sabharwall P, Carroll MC. A phase-field approach to model multi-axial and microstructure dependent fracture in nuclear grade graphite. Journal of Nuclear Materials, 2016, 475: 200-208 doi: 10.1016/j.jnucmat.2016.04.006
[26] Li W, Shirvan K. Multiphysics phase-field modeling of quasi-static cracking in urania ceramic nuclear fuel. Ceramics International, 2021, 47(1): 793-810 doi: 10.1016/j.ceramint.2020.08.191
[27] Tan J, Wu Y, Li Q, et al. Phase field modeling of irradiation-induced shrinkage fracture in TRISO fuel particle. Journal of Nuclear Materials, 2024, 592: 154963 doi: 10.1016/j.jnucmat.2024.154963
[28] Gurtin ME. Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance. Physica D: Nonlinear Phenomena, 1996, 92(3): 178-192
[29] Coleman BD, Noll W. Foundations of linear viscoelasticity. Reviews of Modern Physics, 1961, 33(2): 239-249 doi: 10.1103/RevModPhys.33.239
[30] Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45): 2765-2778
[31] Wu JY. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. Journal of the Mechanics and Physics of Solids, 2017, 103: 72-99 doi: 10.1016/j.jmps.2017.03.015
[32] Pham K, Amor H, Marigo JJ, et al. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics, 2011, 20(4): 618-652 doi: 10.1177/1056789510386852
[33] Tanné E, Li T, Bourdin B, et al. Crack nucleation in variational phase-field models of brittle fracture. Journal of the Mechanics and Physics of Solids, 2018, 110: 80-99 doi: 10.1016/j.jmps.2017.09.006
[34] Ruan H, Rezaei S, Yang Y, et al. A thermo-mechanical phase-field fracture model: Application to hot cracking simulations in additive manufacturing. Journal of the Mechanics and Physics of Solids, 2023, 172: 105169 doi: 10.1016/j.jmps.2022.105169
[35] Naderi M, Amiri M, Khonsari MM. On the thermodynamic entropy of fatigue fracture. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2010, 466(2114): 423-438 doi: 10.1098/rspa.2009.0348
[36] Lucuta PG, Matzke H, Hastings IJ. A pragmatic approach to modelling thermal conductivity of irradiated UO2 fuel: Review and recommendations. Journal of Nuclear Materials, 1996, 232(2): 166-180
[37] Nakajima T, Ichikawa M, Iwano Y, et al. FEMAXI-III: A computer code for the analysis of thermal and mechanical behavior of fuel rods. Japan Atomic Energy Research Institute, 1985
[38] Matzke H, Lucuta PG, Verrall RA, et al. Specific heat of UO2-based SIMFUEL. Journal of Nuclear Materials, 1997, 247: 121-126 doi: 10.1016/S0022-3115(97)00069-X
[39] Luscher WG, Geelhood KJ. Material property correlations: Comparisons between FRAPCON-3.4, FRAPTRAN 1.4, and MATPRO. Pacific Northwest National Lab. (PNNL). Richland, WA, USA, 2010
[40] Kang KH, Ryu HJ, Song KC. Thermal expansion of UO2 and simulated DUPIC fuel. Journal of Nuclear Materials, 2002, 301(2): 242-244
[41] MacDonald PE, Thompson L. MATPRO: A handbook of materials properties for use in the analysis of light water reactor fuel rod behavior. SEE CODE-9502158 Aerojet Nuclear Co., Idaho Falls, Idaho, USA. 1976
[42] Fisher E, Renken CJ. Single-crystal elastic moduli and the hcp→ bcc transformation in Ti, Zr, and Hf. Physical Review, 1964, 135(2A): A482 doi: 10.1103/PhysRev.135.A482
[43] Springer GŚ, Wingeier EW. Thermal conductivity of neon, argon, and xenon at high temperatures. The Journal of Chemical Physics, 1973, 59(5): 2747-2750 doi: 10.1063/1.1680394
[44] Linstrom PJ, Mallard WG. The NIST Chemistry WebBook: A chemical data resource on the internet. Journal of Chemical & Engineering Data, 2001, 46(5): 1059-1063
[45] Permann CJ, Gaston DR, Andrš D, et al. MOOSE: Enabling massively parallel multiphysics simulation. SoftwareX, 2020, 11: 100430 doi: 10.1016/j.softx.2020.100430
[46] Gerasimov T, De Lorenzis L. On penalization in variational phase-field models of brittle fracture. Computer Methods in Applied Mechanics and Engineering, 2019, 354: 990-1026 doi: 10.1016/j.cma.2019.05.038
[47] Chao Correas A, Reinoso J, Cornetti P, et al. On the (lack of) representativeness of quasi-static variational fracture models for unstable crack propagation. Journal of the Mechanics and Physics of Solids, 2024, 186: 105573 doi: 10.1016/j.jmps.2024.105573
[48] Wu JY, Yao JR. A model scaling approach for fracture and size effect simulations in solids: Cohesive zone, smeared crack band and phase-field models. Computer Methods in Applied Mechanics and Engineering, 2022, 400: 115519
[49] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273-1311 doi: 10.1002/nme.2861
[50] Borden MJ, Verhoosel CV, Scott MA, et al. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 2012, 217-220: 77-95 doi: 10.1016/j.cma.2012.01.008
[51] 吕俊男, 杨烁, 赵毅等. 高燃耗下核燃料颗粒内部的应力状态及关键影响因素研究. 原子能科学技术, 2022, 56(12): 2646-2653 (Lyu Junnan, Yang Shuo, Zhao Yi, et al. Stress state and key influencing factor in high burnup nuclear fuel particle. Atomic Energy Science and Technology, 2022, 56(12): 2646-2653 (in Chinese) doi: 10.7538/yzk.2021.youxian.0374 Lyu Junnan, Yang Shuo, Zhao Yi, et al. Stress state and key influencing factor in high burnup nuclear fuel particle. Atomic Energy Science and Technology, 2022, 56(12): 2646-2653 (in Chinese) doi: 10.7538/yzk.2021.youxian.0374