耐事故燃料能提高反应堆核燃料系统在正常运行工况以及严重事故工况下的安全性和可靠性^[1]. 自福岛第一核电站事故发生以来, 耐事故燃料的研发一直是核燃料领域的研究热点, 而筛选高铀密度和高导热系数的燃料体系是耐事故燃料开发的重要思路^[2]. 选用高铀密度的燃料可以延长反应堆的换料周期, 提高核电站运行的经济性; 采用高导热系数的燃料可以降低燃料芯块和包壳的峰值温度, 增强燃料元件的抵御严重事故能力^[3–5]. 与传统UO₂燃料相比, UN, U₃Si₂燃料都具有铀密度高和导热系数高的优点, 且它们的导热系数随着温度的升高而增加^[6]. 然而, UN在氧气、水或蒸汽环境中容易发生氧化腐蚀, 而且其粉末制备和芯块烧结的条件也较为苛刻, 在低烧结温度下制得的UN芯块很难达到满足燃料应用的密度, 而高烧结温度则会对材料的微观结构造成不利影响; U₃Si₂具有更好的抗氧化能力, 但其较低的熔点会降低安全裕度^[7,8]. 因此, UN-U₃Si₂复合燃料充分利用两种燃料相的优点, 以其高抗氧化性、高熔点等优势, 成为了一种极具研发潜力的耐事故燃料.

近年来, 国外学者开展了一系列不同温度、不同时间、不同掺杂比例、多种烧结方式下的UN-U₃Si₂复合燃料的烧结实验^[6,9–11]. Ortega等^[9]采用液相烧结工艺制备UN-U₃Si₂复合燃料样品, 并通过延长烧结时间和进行后处理去除表面孔隙, 实验最终得到了致密度达到94%的样品, 且其铀含量比UO₂高出30%以上. Johnson等^[6]的放电等离子体烧结(spark plasma sintering, SPS)实验结果表明, 在常规烧结温度下得到的复合燃料样品中, 由于UN与U₃Si₂发生相互作用, 形成了一种未知的U-Si-N三元相. White等^[10]在1873—1973 K的温度范围内, 通过液相烧结制备出了U₃Si₂体积分数在10%到40%之间的多种复合燃料, 结果表明, 每种复合燃料的热导率均随温度的升高而增加, UN和U₃Si₂之间的热膨胀差异导致了微观结构中微裂纹的产生. Lopes等^[11]采用SPS方法在1723 K烧结温度和3 min烧结时间下得到了相对密度为98%的UN-10%U₃Si₂ (质量分数)复合燃料芯块, 并从静态蒸汽高压釜测试中得出结论: U₃Si₂相提高了UN基体的稳定性, 燃料芯块的腐蚀机制从晶间开裂转变为了晶内开裂. 然而, 现有的实验技术难以对烧结过程中的微观组织演变进行实时的观测, 导致了无法更加深入地理解两相复合燃料的烧结机理. 因此, 对UN-U₃Si₂复合燃料的烧结过程开展理论模拟研究具有十分重要的意义.

目前, 相场法已被成功应用于合金腐蚀行为^[12]、核燃料裂变气体释放行为^[13]、高熵合金富Cu相析出^[14]等各方面的模拟研究中, 而在陶瓷粉末烧结过程的相场模拟方面也同样取得了一些成果^[15–19]. Liu等^[15]建立了一种新型相场模型描述两相多孔组织的烧结过程. Kumar等^[16]通过相场模拟了两种不等尺寸颗粒的烧结过程, 分析了不同子过程的热力学驱动力. Biswas等^[17]开发出一种考虑方向依赖的界面扩散各向异性和晶粒取向依赖的晶界能各向异性的相场模型, 研究了各向异性性质对固态烧结过程中微观结构演化的影响. Du等^[18]采用相场方法进一步模拟了多孔陶瓷烧结过程中的孔隙变形和晶界迁移. Hötzer等^[19]提出了一种考虑体积扩散、表面扩散和晶界扩散的相场模型描述Al₂O₃颗粒的烧结过程. 然而, 两相烧结过程不仅涉及两相之间溶质元素的扩散, 而且烧结后期的晶粒长大同时受气孔及第二相的影响. 上述模型为两相烧结相场模型的开发提供了一些思路, 但完整准确地描述两相烧结过程的形貌演化等工作需要重新构建自由能泛函数, 并准确求解相场演化方程.

本工作以Wang^[20]提出的固态烧结相场模型为基础, 建立了包含气孔和两种固相的相场模型, 采用张量形式的扩散迁移率系数以及与表面能和晶界能相关的模型参数, 模拟了UN粉末与U₃Si₂粉末的两相烧结过程, 分别对烧结颈的形成与增长过程、三叉晶界的形成过程以及不同体积分数比的两相组织演变过程等进行了研究.

本模型通过引入一系列取向场变量η₁, η₂, ···, η_p和浓度场变量ϕ_α, ϕ_β, ρ来描述不同取向的两种晶粒相和气孔相. η_i (i = 1, 2, ···, p)为非保守型相场变量, 其中η_i (i = 1, 2, ···, r)用来描述α相多晶中不同晶粒的取向, η_i (i = r + 1, r + 2, ···, p)用来描述β相多晶中不同晶粒的取向, 在第i个晶粒内部η_i的取值为1, 其余p – 1个取值为0; ϕ_α为保守型相场变量, 在α相取值为1, 在β相和气孔相取值为0; ϕ_β为保守型相场变量, 在β相取值为1, 在α相和气孔相取值为0; ρ为保守型相场变量, 在α相和β相取值为1, 在气孔相取值为0.

图1为扩散界面相场模型的示意图. 可以看出, 在晶界处η_i的取值从1连续变化为0 (或者从0连续变化为1), ϕ_α的取值从1连续变化为0, ϕ_β的取值从0连续变化为1, 在气孔边界处, ρ的取值从0连续变化为1. 已有研究^[21]表明, 当晶粒个数p > 36时, 晶粒组织演变对p取值的依赖性可以忽略.

本工作两相烧结相场模型的总自由能泛函数F采用以下形式:

式中, 第1项f为局部自由能密度函数; 第2项、第3项、第4项和第5项为梯度自由能, 分别表示气孔表面、α相的相界面、β相的相界面和晶界处的额外自由能; κ_ρ, κ_ϕ和κ_η为梯度能系数, 取决于两种材料气孔表面、晶粒相界面和晶界处的能量和宽度, 单位J/m; V为空间体积.

本工作采用的f表达形式为

式中, A和B为与两种材料气孔表面、晶粒相界面和晶界处能量和宽度有关的参数, 单位J/m³; a_GB为扩散界面参数, 取a_GB = 1.5^[22]; 方程右侧第1项为双势阱函数, 在晶粒内部与气孔内部取得极小值; 方程右侧第2项花括号中为多势阱函数, 在晶粒内部取得极小值. 本工作所采用的自由能泛函数能够保证在各个晶粒内部和气孔内部取得p + 1个极小值.

F中各参数与两种材料的物理参数之间的关系为^[23]

式中, 插值函数$ {\varphi _\alpha } = {{\displaystyle \sum\nolimits_{i = 1}^r {\eta _i^2} } \Big/ {\displaystyle \sum\nolimits_{i = 1}^p {\eta _i^2} }} $, $ {\varphi _\beta } = {{\displaystyle \sum\nolimits_{i = r + 1}^p {\eta _i^2} } \Big/ {\displaystyle \sum\nolimits_{i = 1}^p {\eta _i^2} }} $; $ \gamma _{{\text{gb}}}^\alpha $和$ \gamma _{\text{s}}^\alpha $, $ \gamma _{{\text{gb}}}^\beta $和$ \gamma _{\text{s}}^\beta $分别为两种材料的晶界能和表面能, 单位J/m²; δ为扩散界面的宽度, 单位m.

假设烧结过程中的质量守恒, 上述相场模型中ρ, $ {\phi _\alpha } $和$ {\phi _\beta } $的演化动力学方程为Cahn-Hilliard方程^[24]:

式中, t为时间; ${\boldsymbol M} $为迁移率张量, 表达形式为^[23]

其中, ν_m为摩尔体积, R为气体常数, T为温度, $ {{\boldsymbol{D}}_\alpha } $和$ {{\boldsymbol{D}}_\beta } $为两种材料的扩散系数张量.

本工作采取的体系扩散机制包括表面扩散和晶界扩散. 因此, 可以认为$ {{\boldsymbol{D}}_\alpha } = {\boldsymbol{D}}_\alpha ^{\text{s}} + {\boldsymbol{D}}_\alpha ^{{\text{gb}}} $, $ {{\boldsymbol{D}}_\beta } = {\boldsymbol{D}}_\beta ^{\text{s}} + {\boldsymbol{D}}_\beta ^{{\text{gb}}} $, 其中, $ {\boldsymbol{D}}_\alpha ^{\text{s}} $和$ {\boldsymbol{D}}_\alpha ^{{\text{gb}}} $, $ {\boldsymbol{D}}_\beta ^{\text{s}} $和$ {\boldsymbol{D}}_\beta ^{{\text{gb}}} $分别为两种材料的表面扩散系数张量和晶界扩散系数张量, 表达形式为^[23]

式中, $ D_\alpha ^{\text{s}} $和$ D_\beta ^{\text{s}} $为表面扩散系数标量, 取决于两种材料的性质; ${\boldsymbol T}^{\rm s} $为表面投影张量, 保证表面扩散在晶粒与气孔接触表面的切线方向; I为单位张量; 符号$\otimes $为张量积; ${\boldsymbol n}_{\rm s} $为垂直于表面的单位法向量.

式中, $ D_\alpha ^{{\text{gb}}} $和$ D_\beta ^{{\text{gb}}} $为晶界扩散系数标量, 取决于两种材料的性质; ${\boldsymbol{T}}^{{\text{gb}}} $为晶界投影张量, 保证晶界扩散在晶界的切线方向; ${\boldsymbol{n}}_{{\text{gb}}} $为垂直于晶界的单位法向量.

η_i的演化动力学方程为Allen-Cahn方程^[24]:

式中, L为表征晶界迁移率的系数, 与两种材料的物理参数之间的关系为^[25]

其中, M_b为晶界迁移率.

根据(7)式、(8)式、(9)式的Cahn-Hilliard方程与(19)式的Allen-Cahn方程, 对所有相场变量进行求解, 可以计算出组织演变过程中任意时刻的微观结构与空间分布. 本工作通过编写程序进行模拟, 采用有限差分法求解相场方程, 在时间上采用显式Euler算法对非线性偏微分方程进行离散化, 在空间上采用五点差分法求解Laplace项^[26].

为了直观地表征两相烧结过程, 本工作中可视化变量φ为

本工作相场模型的变量与参数采用无量纲形式. 根据(7)式、(8)式与(9)式, 保守场演化方程中含有量纲的量为t, $\boldsymbol M $, F和$ \nabla $, 通过选取参考空间长度$l^* $、参考时间$t^*$和参考能量密度$\varepsilon^*$作为参考物理量, 可以得到无量纲Hamilton算符$ \tilde \nabla = {l^ * }\nabla $、无量纲时间$\tau=t/t^*$和无量纲自由能$ \tilde F =F/\varepsilon^*$. 将无量纲量代入(7)式、(8)式与(9)式, 得到变形后的方程为

由此可得无量纲迁移率张量$ \tilde {\boldsymbol{M}} = {\boldsymbol{M}}{\varepsilon^*}{t^*}/l^{*2} $.

根据(19)式, 非保守场演化方程含有量纲的量为t, L, F和$\nabla $. 将无量纲量代入(19)式, 得到变形后的方程为

由此可得无量纲迁移系数$ \tilde L = L{\varepsilon ^*}{t^*} $.

本工作取$ {t^*} = 1/\left( {L{\varepsilon ^*}} \right) $, $ {\varepsilon ^*} = B $; $ {l^*} $的取值与δ有关, 对应无量纲化之后一个单位长度占有的界面宽度, 本工作取$ {l^*} = \delta /m $; m为扩散界面在模拟中占据的格点数, 本工作取m = 3.

根据所选取的参考物理量, 无量纲迁移率张量$\tilde {\boldsymbol{M}}$与无量纲总自由能泛函数$ \tilde F $中A, B, $ {\kappa _\rho } $, $ {\kappa _\phi } $, $ {\kappa _\eta } $等参数对应的无量纲形式分别为

本工作利用相场模型对UN粉末与U₃Si₂粉末在1823 K温度下的两相烧结过程进行了数值模拟, 以UN为α相, U₃Si₂为β相, 模拟中用到的UN和U₃Si₂的物理参数如表1和表2所列^[27–32].

将UN和U₃Si₂的物理参数代入上述方程中, 求得各无量纲参数的具体数值如表3所列.

本工作采用相场模型, 模拟了双晶粒结构中UN粉末与U₃Si₂粉末的两相烧结过程. 所选模拟区域差分网格尺寸为256 × 256, 采用周期性边界条件. 分别研究UN粉末与U₃Si₂粉末两相烧结组织的形貌演化过程、烧结颈的形成、平衡二面角以及不同晶粒尺寸对组织演变过程的影响.

图2为2个不同相等尺寸的圆形晶粒在烧结过程中形貌演化的相场模拟. 其中晶粒的半径R = 30, 黄色的部分代表UN相晶粒, 绿色的部分代表U₃Si₂相晶粒. 从图2(a)可以看出, 在烧结过程初期, U₃Si₂相晶粒具有比UN相晶粒更高的表面能, 使其在烧结颈形成过程中更趋向于减小表面积, 颈部附近的晶粒表面向晶粒内部凹陷变形的程度更明显; 在图2(a)—(d)两相晶粒烧结过程中, 先后发生了烧结颈的形成和增长, 这一现象与单相晶粒烧结模拟^[33]的结果一致; 仅考虑表面扩散和晶界扩散的模拟烧结颈增长较缓慢, 晶粒位置无明显移动.

图3为基于2个不同相等尺寸晶粒烧结模拟数据的烧结颈对数增长曲线. 由图3可以看出, 烧结颈的增长曲线斜率前后发生改变, 这表明烧结颈的增长规律并非简单的幂函数关系, 曲线可以划分为2个增长阶段: 烧结初期, 烧结过程受表面扩散主导, 颈部快速形成并增长, 斜率为0.44左右; 烧结中期, 表面扩散作用减小, 晶界扩散作用增大, 烧结颈的增长速率显著下降, 烧结颈持续缓慢增长, 斜率为0.21左右, 满足烧结颈曲线增长规律^[34]. 随着烧结过程的继续进行, 烧结颈将在后期达到稳定状态, 颈部尺寸较大, 晶粒之间紧密接触, 不再适合用烧结前中期的增长规律对烧结颈进行描述.

研究表明, 在烧结过程中, 相邻的2个晶粒在晶粒相与气孔相的相界处形成一个平衡二面角Ψ, Ψ取决于材料的表面能γ_s和晶界能γ_gb, 在各向同性条件下, 晶界处的平衡二面角方程为^[35]

图4分别为不同演化时间下2个等尺寸晶粒烧结形成的平衡二面角. 本工作中UN相晶界能与表面能的比值为$ {{\gamma _{{\text{gb}}}^\alpha } {/ } {\gamma _{\text{s}}^\alpha }} = 0.5 $, 对应计算值为Ψ = 151°, 如图4(a)和图4(d)所示; U₃Si₂相晶界能与表面能的比值为$ {{\gamma _{{\text{gb}}}^\beta } {/ } {\gamma _{\text{s}}^\beta }} = 0.65 $, 对应计算值为Ψ = 142°, 如图4(b)和图4(e)所示; 2个同相晶粒烧结形成的最终平衡二面角均与理论计算值吻合较好. 从图4可以看出, UN相晶粒与U₃Si₂相晶粒烧结形成的最终平衡二面角介于2个同相双晶粒的最终平衡二面角之间, 达到稳定状态后的平衡二面角随时间无明显变化. 这证明了2个不同相晶粒烧结形成的最终平衡二面角只由两种材料的特性所决定, 即与两种材料的晶界能和表面能的不同比值密切相关.

图5为2个不同相不同尺寸晶粒烧结过程的相场模拟. 由图5可以看出, 2个不同相不同尺寸晶粒的烧结过程与Sun等^[36]模拟的2个同相不同尺寸晶粒的烧结过程有着显著的差异: 2个同相晶粒烧结, 烧结颈向尺寸较小的晶粒一侧迁移, 大晶粒的面积不断增大, 而小晶粒的面积不断减小, 直到最终消失^[36]; 2个不同相晶粒烧结, 烧结颈向尺寸较小的晶粒一侧迁移, 迁移一段距离后稳定存在, 小晶粒的形状不断发生改变, 最终呈现为扁圆形, UN相大晶粒和U₃Si₂相小晶粒各自的面积无明显变化, 大晶粒并未吞噬小晶粒.

在2个不同相不同尺寸晶粒的烧结过程中, 晶粒面积和烧结颈尺寸随时间演化的曲线如图6所示. 从图6可以看出, 不同尺寸的2个不同相晶粒烧结过程可以划分为3个阶段: 第1阶段, 烧结颈快速形成与增长, 大小晶粒的面积几乎不变; 第2阶段, 烧结颈缓慢增长, 大小晶粒的面积基本保持稳定; 第3阶段, 烧结颈最终达到稳态, 大小晶粒的面积仍然保持不变. 在烧结过程中, 涉及表面扩散和晶界迁移两种动力学机制, 这两种机制的传质速率分别与晶粒气孔接触表面和晶界的曲率相关^[37]. 烧结颈形成之后, 大小晶粒通过表面扩散机制向烧结颈输送物质; 由于2个晶粒表面曲率的不同, 导致传质速率不同, 烧结颈向尺寸较小的晶粒一侧生长, 形成了弯曲的晶界; 在弯曲晶界曲率的驱动下, 晶界向小晶粒所在的凹侧迁移^[36]; UN与U₃Si₂之间不存在相互转化, 不同相晶粒之间不相互传递物质, 物质仅在各自的同相区域内, 向靠近烧结颈的一侧传输, 由晶界迁移引起的晶粒合并无法进行. 因此, 曲率不同的2个不同相晶粒之间不会发生“大吞小”现象.

图7为在UN相晶粒与U₃Si₂相晶粒的4种不同个数比例下, 3个晶粒烧结形成的三叉晶界, 演化的总时间步长为10⁷ $ \Delta t $. 从图7(a)和图7(d)可以看出, 在3个同相晶粒的烧结过程中, 形成 了夹角为120°的稳定三叉晶界, 这是因为3条晶界具有相同的能量, 且120°的角度使得晶界的 能量分布均匀, 从而达到了晶界能平衡的状态^[36]. 从图7(b)和图7(c)可以看出: 在2个同相晶粒与另一个不同相晶粒的烧结过程中, 形成了夹角偏离120°的三叉晶界, 这是由于两相的晶界能存在差异, 为了实现晶界能的平衡并使体系的能量最小化, 夹角需要根据不同晶界能的大小进行合理的分配; 两相晶界能的差异越大, 夹角偏离120°的程度越大; U₃Si₂相晶粒具有比UN相晶粒更高的晶界能. 因此, 在U₃Si₂相晶粒一侧的夹角趋向于大于120°, 而在UN相晶粒一侧的夹角趋向于小于120°.

在4种不同的两相晶粒个数比例下, 三叉晶界处的气孔率随时间演化的曲线如图8所示. 在模拟研究中考虑了两种扩散机制, 其中表面扩散被认为是一种非致密化过程, 主要促进烧结颈的形成和生长; 而晶界扩散被视为了一种极其重要的致密化机制^[34], 涉及到了空位的传输, 是造成三叉晶界处气孔收缩的根本原因. 从图8可以看出, 3个UN相晶粒的气孔收缩过程明显快于其他三种情况. 在晶界扩散过程中, 物质原子从晶界处迁移向与气孔接触的颈部表面, 气孔空位沿晶界扩散并最终湮灭; 由于U₃Si₂相晶粒具有比UN相晶粒更高的晶界能, 与其形成的晶界存在更高的能量势垒, 使得晶界附近的原子或空位在扩散过程中需要克服更大的能量障碍, 导致了原子与气孔空位的扩散速率减慢. 因此, 在包含U₃Si₂相晶粒的3个晶粒烧结过程中, 三叉晶界处的气孔收缩过程显著减缓. 而在包含UN相晶粒与U₃Si₂相晶粒的情况下, 两相的晶界能、表面能和扩散系数的不同都会影响三叉晶界处的气孔收缩, 使得气孔收缩速度的快慢变得更加复杂.

为了研究UN相与U₃Si₂相的不同体积分数比对烧结过程中UN-U₃Si₂两相复合燃料微观组织演变的影响, 构建了4种具有不同两相体积分数比的两相多晶组织, 并进行了两相烧结相场模拟. 其中, α相的体积分数分别为80%, 70%, 60%和50%, 而β相的体积分数则分别为20%, 30%, 40%和50%. 图9展示了这4种不同体积分数比的两相多晶组织的初始形貌.

图10为经过8 × 10⁶ $ \Delta t $的演化时间后, 4种不同体积分数比的两相多晶烧结组织的最终形貌. 通过对比图9与图10同种情况下的组织形貌可以发现, 晶粒个数和晶粒间的气孔显著减少, 这表明同相晶粒之间发生了合并, 气孔发生了收缩, 晶界扩散起到主要作用^[34]. 从图10(a)—(d)可以看出, 在α相体积分数大于β相体积分数的情况下, β相晶粒基本呈现出嵌入在α相基体中的分布, 表明在烧结过程中α相的晶粒生长占据了主导地位. 在α相晶粒与β相晶粒的接触界面处, 均存在有尺寸较小的α相晶粒或β相晶粒, 未完全与各自同相的大晶粒发生合并. 这表明α相或β相的存在会对另一相的晶界产生阻碍作用, 阻碍其晶界的迁移, 导致其小晶粒无法快速合并, 造成其晶粒长大速率减慢. 在图10(a)—(d)中, 随着β相的体积分数增加, α相的体积分数减少, 两相晶粒的接触面积会增大, 两相之间存在的阻碍作用也会增强. 在图9中, 存在只与不同相晶粒相互接触的α相小晶粒或β相小晶粒, 而在图10中, 这些晶粒消失或其面积明显减少, 这表明α相或β相晶粒也可以通过晶粒迁移方式, 向同相晶粒传输物质, 物质通过另一相晶粒的晶界与相界面进行扩散, 迁移的晶粒会逐渐变小, 最终消失^[38].

1) 建立了两相烧结相场模型, 将两相的表面能和晶界能与相场模型的参数关联起来, 模拟了UN-U₃Si₂复合燃料的烧结过程. 模拟结果表明, 烧结颈形成过程中, 具有更高表面能的U₃Si₂相晶粒在颈部附近的表面变形程度更大; 由不同相晶粒烧结形成的最终平衡二面角, 其取值由两相的晶界能和表面能的不同比值所决定; 不同相晶粒之间不存在大晶粒吞噬小晶粒的现象.

2) 3个晶粒烧结过程中的气孔收缩与三叉晶界演变的模拟结果表明, 为实现晶界能的平衡, 不同相晶粒烧结形成的三叉晶界夹角在较高晶界能的U₃Si₂相一侧大于120°, 在较低晶界能的UN相一侧小于120°; 晶界处的高能量势垒会限制气孔空位的晶界扩散, 使三叉晶界处的气孔收缩变慢.

3) 不同体积分数比值的UN-U₃Si₂两相多晶烧结组织形貌演化的模拟结果表明, 烧结组织中存在同相晶粒间的合并与气孔的收缩, 晶界扩散在两相烧结过程中起主要作用; 体积分数较大相的晶粒生长占据主导地位; 两相体系中的一相会阻碍另一相的晶界迁移; 同相晶粒之间能通过晶粒迁移方式进行物质传输.