一种双变量有限元相场方法，用于铌铁素体钢中相分离的扩散-弹性耦合建模

《Mechanics of Materials》：A dual-variable FE phase-field approach for coupled diffusion-elasticity modeling of phase separation in bainitic steel

　　K. Shiva Reddy | Amirtham Rajagopal
印度海得拉巴印度理工学院土木工程系，Sangareddy 502284

**摘要**
本研究提出了一个稳健的有限元实现方法，用于Cahn–Hilliard方程在相场建模中的应用，该方法既包括与机械变形的耦合情况，也包括不包括耦合的情况。通过采用一种将化学势和碳浓度视为独立主要变量的混合公式，原始的四阶方程被分解为两个二阶方程，从而可以使用标准的C0-连续元素。该框架在Abaqus/Standard中完全实现，使用了自定义用户定义元素（UEL）子程序以及用户定义的输出变量（UVARM）子程序。引入了两层元素技术和交错算法来有效处理耦合的扩散和弹性问题。通过模拟贝氏体钢中的碳扩散和相分离过程，验证了这种方法的有效性，突显了弹性在控制扩散动力学和产生的微观结构中的重要作用。

**引言**
Cahn–Hilliard（CH）扩散方程是相场（PF）模型的一个基本组成部分，描述了由化学势梯度驱动的扩散过程（Cahn和Hilliard，1958年；Cahn和Hilliard，1959年；Cahn，1961年）。它考虑了向下和向上的扩散，使得可以模拟合金中的旋涡分解（Wheeler等人，1992年；Kim等人，1999年）。化学势来源于系统的总自由能，通常包括一个非凸的构型项，该项驱动相分离，以及一个界面/梯度项，用于规范成分梯度。CH方程已广泛应用于材料科学中，用于模拟相分离、粗化、微观结构演化、固体中的扩散，以及与力学和其他多物理现象的耦合（Travasso等人，2011年；Chen等人，2019年；Chen等人，2020年；Cao等人，2020年；Bartels等人，2021年；Wang等人，2023年；Lv等人，2024年；Wang和Lu，2025年）。它为理解相分离和界面演化的动力学提供了框架，其中晶格畸变和相变过程中的缺陷形成会影响最终的微观结构，从而影响材料的物理和机械性能（Elliott，1989年；Chen，2002年；Song等人，2015年；Tourret等人，2022年；Ma等人，2025年）。
空位的存在/演化在相变中起着重要作用，例如在马氏体转变中（Javanbakht等人，2021年）。CH方程经常与Ginzburg–Landau（GL）方程和弹性方程相结合，用于模拟空位扩散和空洞或纳米孔洞的演化。Javanbakht和Ghaedi（2021年）提出了一个马氏体相变与纳米空位之间的相互作用模型，该模型能够捕捉它们的同时演化，并考虑了关键相互作用特征——如混合（或自由）能量的相位依赖性，以及在本构关系中包含转变和失配应变张量。特别是，通过将纳米空位浓度（或相关守恒序参量）视为一个守恒场，并且其时间变化率受CH方程控制，可以模拟纳米空位纳米结构的演化（Yu和Lu，2005年）。Javanbakht和Ghaedi（2020b）引入了纳米孔洞的固气界面的热力学一致界面应力，其中CH方程控制纳米孔洞浓度的演化（从完全固态的0变化到气态的1）。
CH方程是一个四阶抛物型偏微分方程（PDE），这对标准数值方法提出了重大挑战。传统的Galerkin公式需要C1-连续近似来弱化四阶算子，这导致在自适应有限元中使用立方Hermite多项式（Stogner等人，2008年）或在等几何分析中使用B样条（Hughes等人，2005年；Gómez等人，2008年）等特殊基函数。其他的C1-连续方法包括不连续Galerkin方法（Wells等人，2006年）、谱Galerkin方法（Chen和Shen，1998年；Feng等人，2009年），以及C1-连续的无网格方法，如自然元素方法（Fischer等人，2011年；Rajagopal等人，2010年）。
为了避免C1-连续性要求的一种广泛采用的策略是通过引入辅助变量（通常是化学势（Elliott等人，1989年）或微观形态非局部浓度（Ubachs等人，2004年）将四阶CH方程重新表述为两个耦合的二阶方程系统。这种分解使得可以使用混合有限元方法，其中相场浓度和辅助变量都被视为独立的未知数。混合公式自然允许使用C0-连续的Lagrangian基函数（例如线性或二次元素），这些函数比C1-连续近似（如Hermite元素或样条）更简单易实现且在计算上成本更低，同时仍然能够提供稳健和准确的解（Fischer等人，2011年）。Kaessmair等人（2014年）和Kaessmair和Steinmann（2016年）研究了C1-连续方法和混合二阶有限元方法的比较计算分析。Kireeti等人（2025年）的近期工作中详细讨论了求解四阶CH型相场模型的各种计算方法。
在我们之前的研究中，我们已经在不同的材料系统中解决了CH方程的数值问题。在Kireeti等人（2025年）的工作中，开发了一个用于相分离与机械变形耦合的混合有限元（FE）框架。在那里，引入了一个非局部浓度场作为辅助变量，将四阶CH方程重新表述为一个更易于处理的二阶方程系统，同时自然地结合了非局部相互作用。该模型应用于锂离子电池电极材料中的锂扩散，并通过用户定义（UEL/UMAT）子程序在商业有限元软件Abaqus/Standard中实现。在另一项工作（Shiva Reddy等人，2025a）中，直接使用C1-连续的自然邻居Galerkin（无网格）方法解决了四阶CH方程。该研究关注了锆合金中的氢扩散和形态变化，这些合金广泛用于核工业中的包壳和结构部件。
有限元软件COMSOL在许多多物理问题中得到了广泛应用（Javanbakht和Ghaedi，2020a；Cui等人，2023年；Javanbakht等人，2023年；He等人，2025年）。由于COMSOL Multiphysics具有专用的PF接口，它非常适合处理四阶CH方程，可以实现快速、完全基于GUI的、无需编码的设置。它自动分割PDE，处理非线性自由能，确保质量守恒，并支持简单的多物理场耦合——使其在纯相场模拟（旋涡分解、Allen–Cahn/Ginzburg–Landau类型）或轻度耦合的相场模拟中更加优越和用户友好。相比之下，Abaqus/Standard缺乏原生的PF支持，需要自定义的用户元素UEL子程序或在某些情况下需要用户材料UMAT/UMATHT（本构耦合）子程序来实现CH方程。然而，Abaqus在复杂的力学主导的多物理问题中具有显著优势，因为它拥有行业领先的非线性求解器和有限应变塑性、超弹性、粘塑性、晶体塑性、大变形接触、显式动力学和断裂力学的元素技术。当相分离与机械场强烈耦合时（例如应力辅助扩散、大塑性变形下的相变），Abaqus提供了更好的稳健性、收敛性和与大量预验证的用户材料库的集成。用户子程序的优势在于可以定义用户定义的输出解变量，相比之下COMSOL的结构模块则不具备这一优势。有兴趣的读者可以参考Chester等人（2015年；Cui等人，2021年；Dyck等人，2025年；Navidtehrani等人，2025年；Castro等人，2026年）中关于Abaqus在某些耦合问题中的应用示例。
本研究采用了CH方程的经典双变量（混合）公式，其中化学势与浓度被视为独立的未知数。这种方法允许使用简单的C0-连续有限元素，极大地简化了实现过程。在这里，我们研究了贝氏体钢中的碳扩散和相分离——贝氏体是一种在高强度、坚韧的微观结构中形成的，在奥氏体的等温转变过程中形成。文献中的实验和理论研究表明，贝氏体形成过程中的碳分配会导致明显的相分离，形成富碳（奥氏体）和贫碳（铁素体）区域。在典型的贝氏体转变中，PF模型应捕捉贝氏体铁素体从奥氏体中的位移生长，以及驱动碳分配和碳化物沉淀的扩散过程（Düsing和Mahnken，2018年；Nagaraja和Antretter，2024年）。在低贝氏体钢中——特别是在贝氏体铁素体相内——可以通过仅模拟碳扩散（向上/非标准）和相应的相分离来近似内部碳化物沉淀，而不引入碳化物相本身的显式非守恒序参量（参见Düsing和Mahnken，2016a；Düsing和Mahnken，2019年）。这种方法是有效的，因为贝氏体铁素体是通过位移机制快速形成的，一旦过饱和铁素体建立，内部的碳化物沉淀主要是由碳重新分布驱动的扩散相分离过程。关于贝氏体从奥氏体转变的更多细节以及低贝氏体转变和高贝氏体转变的定义在Düsing和Mahnken（2016b）中进行了讨论。当前的模拟成功捕捉了碳扩散和相分离，提供了对贝氏体转变过程中碳的动力学和空间分布的更深入理解。本研究的新颖方面可以总结如下：
•开发了一个稳健且灵活的FE框架，能够处理未耦合的相场CH型扩散和耦合的相分离力学相互作用。
•引入了双层元素方法，使得用户元素（UEL）与辅助用户定义输出变量（UVARM）子程序的有效集成成为可能。
•实现了交错求解算法，以便在每个时间增量内依次解决机械平衡和相场演化方程。

**论文的其余部分组织如下：** 第2节概述了耦合的PF模型，第3节讨论了数值实现，第4节展示了用于验证的数值模拟，最后第5节总结了研究并指出了潜在的未来研究方向。

**附录**
**耦合多场模型**
图1展示了二维弹性固体连续体Ω的示意图，其中耦合问题涉及三个主要场变量：浓度c、化学势μ和位移u。对于机械问题，在Γu上规定位移Dirichlet边界条件u?，而在Γt上施加牵引Neumann边界条件t?，且Γu∪Γt=Γ以及Γu∩Γt=0?。对于浓度，在Γξ上施加通量Neumann边界条件ξ?=?c?n?。

**数值实现**
在本节中提供了数值实现方法。如上所述，双变量方法将原始的四阶微分方程转换为两个耦合的二阶微分方程系统。这种重新表述提供了显著的计算优势，因为它允许使用标准的Lagrangian有限元近似。一个关键的好处是放宽了对测试（或权重）函数的可微性要求。

**数值示例**
在本节中，通过数值示例说明了所讨论的耦合扩散-变形模型。第3节描述的框架应用于解决贝氏体钢中的碳扩散问题。模拟中使用的材料和数值参数列在表3中。由于其中大部分参数还是暂定的，因此需要进一步研究。图2显示了局部化学自由能随浓度的变化情况。

**结论**
通过使用将浓度和化学势视为独立主要变量的混合有限元公式，成功解决了四阶Cahn–Hilliard方程。这种策略将问题简化为两个二阶方程，并允许使用标准的C0-连续Lagrangian元素，从而大大简化了实现过程。整个框架已在Abaqus/Standard中通过UEL和UVARM子程序实现，得到了一个稳健的结果。

**作者貢献声明**
K. Shiva Reddy：撰写——审阅与编辑、撰写——原始草稿、可视化、验证、软件、方法论、概念化。
Amirtham Rajagopal：撰写——审阅与编辑、可视化、监督、方法论、概念化。

**利益冲突声明**
作者声明他们没有已知的竞争利益或个人关系可能影响本文报告的工作。