MEMS导热薄膜热力耦合下的相场损伤建模

Phase-Field Modeling of Thermo-Mechanical Damage in MEMS Thermally Conductive Films

摘要: 随着MEMS技术向纳米尺度发展，MEMS器件在复杂多变的多物理场耦合载荷作用下，其内部预先存在的或在使用中萌生的微裂纹会发生演化，这种演化是导致器件力学性能发生不可逆劣化的主要物理机制之一。将相场理论与损伤力学相结合，针对热力耦合环境下的MEMS导热薄膜的裂纹演化行为开展了理论建模与分析。所构建的相场损伤模型能够准确描述导热材料在热应力作用下裂纹的萌生与扩展过程。该模型综合考虑了材料初始缺陷对断裂行为的影响，并引入了热力耦合边界条件下的多场交互效应。在ABAQUS平台上，基于UEL子程序的数值模拟，系统探究了不同热边界条件及初始损伤对裂纹扩展规律的影响。所建立的相场损伤模型为揭示导热材料在热力载荷下的断裂机制提供了可靠的理论框架，并通过数值仿真验证了其有效性，为未来MEMS器件的可靠性设计与性能评估提供重要的理论支撑与指导。

Abstract: With the advancement of MEMS technology towards the nanoscale, microcracks pre-existing or initiated during service in MEMS devices evolve under complex multiphysics coupled loads. This evolution is one of the primary physical mechanisms leading to irreversible degradation of the mechanical performance of the devices. By integrating phase-field theory with damage mechanics, this study conducts theoretical modeling and analysis of crack evolution in MEMS thermal conductive films under thermomechanical coupling environments. The developed phase-field damage model accurately captures the initiation and propagation of cracks in thermally conductive materials under thermal stress. The model incorporates the influence of initial material defects on fracture behavior and introduces multiphysics interaction effects under thermomechanically coupled boundary conditions. Numerical simulations based on the UEL subroutine on the ABAQUS platform systematically investigate the effects of different thermal boundary conditions and initial damage values on crack propagation behavior. The established phase-field damage model provides a robust theoretical framework for understanding the fracture mechanisms of thermal conductive materials under thermomechanical loads, validated through numerical simulations. This work offers important theoretical support and guidance for the reliability design and performance evaluation of future MEMS devices.

1. 引言

MEMS技术在近五十年来的持续发展和技术革新中，凭借器件小型化、低功耗、智能化和高集成度等特点使其传感、互连层和执行器等核心部件在航空航天、汽车电子和能源装备的应用前景和性能提升提供了更多可能[1]-[3]。

由于MEMS器件工作环境的特殊性与复杂性，常伴随显著的温度变化与力学载荷的共同作用，因此对材料的稳定性提出了极高要求。在此背景下，导热陶瓷因其独特性能成为关键材料，它既能通过优异的导热性快速均化温度，抑制局部热点的产生，又能凭借高刚度与耐磨性承受机械应力[4]-[6]。然而，在热‑力耦合场中，导热陶瓷虽然可以通过快速传递热量来抑制局部高温，但材料的制备和加工的工艺产生的初始缺陷(微裂纹、微孔洞、晶界等)会成为裂纹萌生与扩展的起点[7]-[9]。这些内部存在的初始缺陷成为了应力集中点，会在热循环载荷下逐渐扩展；同时，陶瓷与相邻材料的热膨胀系数失配会在界面处产生巨大的热应力，反复作用下导致界面分层或开裂。此外，在高温环境下，晶界可能发生蠕变或氧化，进一步降低材料的导热性能和机械强度[10]-[12]。导热陶瓷的失效本质上是热–力耦合作用下材料缺陷演化、界面失稳及性能退化的结果[13]。因此，从服役环境的实际工况出发，深入探究MEMS材料在多物理场耦合作用下的断裂与失效行为，对于阐明其性能退化与失效的内在物理机制、建立精准的可靠性评估模型，以及实现其服役寿命的科学预测具有至关重要的理论指导意义和工程应用价值。

近年来，在关于热诱导裂纹的行为及其对强度或其它参数的影响或已经进行了多项实验研究。Jiang等人[14]通过有限元方法，从实验和数值两方面研究了热冲击对陶瓷裂纹模式的影响。Shao等人[15]研究了陶瓷残余强度随裂纹扩展阶段的变化规律。Yousef等人[16]通过有限元模拟实验观察了大晶粒多晶氧化铝中的微裂纹演化过程，并报告了温度降低导致杨氏模量等宏观结构性能的变化。这些热疲劳裂纹扩展问题的研究中，普遍依赖预设裂纹路径或经验判据，这在复杂服役载荷条件下难以准确表征多裂纹的相互作用、分叉以及不规则扩展行为，限制了其在实际工程中的普适性与精确性[17]-[19]。相较而言，相场法所开发的损伤模型能够通过连续损伤变量描述裂纹的演化过程，不需要显式跟踪裂纹界面，能够自然处理裂纹的萌生、扩展、分问题[20]-[22]。更重要的是，相场与热传导方程和力学平衡方程可以有机耦合，建立统一的热–力–相场计算框架，从而为研究热应力驱动下的裂纹扩展机理提供了更加高效和准确的数值工具。

本文在考虑器件服役前的制备与加工阶段就产生了一定的微缺陷的前提下，依据连续损伤力学、断裂力学推导出导热陶瓷在热–力耦合下的相场损伤模型，用来描述器件内部的损伤演化过程。通过将模型应用在ABAQUS中，二次开发UEL子程序进行裂纹模拟。基于数值模拟的结果验证所建立的相场损伤模型的准确性，之后对导热陶瓷的断裂行为进行模拟，考察模型在不同工况下对裂纹扩展问题的适用性与表现。结果表明MEMS器件在应用和服役阶段裂纹的预测和模拟方面具有强力的优势和潜力。

2. 热力耦合下的相场损伤模型

2.1. 损伤的定义

由于工艺条件存在不确定性，在MEMS器件的制备过程中材料内部往往会形成的微孔洞与微裂纹等初始缺陷。这类微观缺陷难以彻底消除，会对器件的性能与可靠性造成潜在的不利影响，尤其是在长期服役过程中，这些缺陷可能会加速材料的疲劳与损伤过程。为解决这一问题，可以采用基于不可逆热力学的连续损伤理论。核心方法是假设代表性体积单元(RVE)的存在。RVE被引入作为均匀化处理材料结构损伤变量的手段，因此，损伤变量$D$ 可以通过RVE内部微孔洞的扩展来定义，进而表征整个结构的损伤程度。

如图1所示，采用代表性体积单元(RVE)模型基于连续损伤力学(CDM)理论，利用RVE来定义材料结构中的损伤，具体如下：

$D=\frac{V_D}{V_0}=\frac{d_i\times l_i}{d_i\times d_i}=\frac{l_i}{d_i}$ (1)

其中，$V_0$ 表示初始面积，$V_D$ 表示受损区域面积。$l_i$ 表示在第$i$ 个体积单元中由于损伤引起的裂纹长度，$d_i$ 则表示第$i$ 个体积单元的特征尺寸。

相场模型基于变分原理，能够很好地处理复杂的裂纹拓扑结构，包括多裂纹的相互作用和耦合。考虑材料的初始微缺陷，$l_i$ 可以通过引入相场变量${\varphi }_i$ 来定义为：

$l_i=l_c+\left(d_i-l_c\right)\cdot {\varphi }_i$ (2)

其中，$l_c$ 表示由于材料加工工艺导致的初始微缺陷的起始裂纹长度。

将方程(2)代入方程(1)中，损伤变量$D$ 可以定义为：

$D=D_0+\left(1-D_0\right){\varphi }_i$ (3)

其中，$D_0=l_c/d_i$ 表示为由于材料加工工艺导致的损伤值。

根据连续损伤力学(CMD)理论，使用有效应力取代名义应力，可以获得更加准确的解决方案，具体如下：

${\sigma }_e=\frac{1}{1-D}{\sigma }_n$ (4)

其中，${\sigma }_e$ 是有效应力，${\sigma }_n$ 是名义应力。

将上式带入，可以得到：

${\sigma }_e=\frac{1}{\left(1-D_0\right)\left(1-{\varphi }_i\right)}{\sigma }_n$ (5)

Figure 1. MEMS sensor fabrication process defects and representative volumes

图1. MEMS传感器加工工艺缺陷以及代表性体元

2.2. 含损伤变量的导热材料本构方程

考虑在2D条件下一个带有裂纹的导热材料$\Omega$ ，其边界为$\partial \Omega$ ，在热力耦合下的导热薄膜的力学行为不仅由应力–应变关系决定，还与温度梯度和热流耦合有关。因此，导热薄膜的本构方程为：

$\{\begin{array}{l}{\sigma }_n=C:\left(\epsilon -\alpha \Delta T\right)\\ q=K\cdot M\end{array}$ (6)

其中，$C$ 表示为弹性刚度张量，$\epsilon$ 为应变张量，表示为$\alpha$ 热膨胀系数张量，$M$ 是温度梯度矢量，$q$ 为热流矢量，$K$ 热导率。

上式中的应变张量$\epsilon$ 表示为：

$\{\begin{array}{l}\epsilon =\frac{\nabla u+{\left(\nabla u\right)}^{\text{T}}}{2}\\ M=-\nabla T\end{array}$ (7)

其中，$\nabla$ 是哈密顿算子，$u$ 是位移矢量，$T$ 代表温度。

综上所述，得到了考虑损伤的导热材料损伤本构方程：

${\sigma }_e=\frac{C:\left(\epsilon -\alpha \Delta T\right)}{\left(1-D_0\right)\left(1-\varphi \right)}$ (8)

2.3. 裂纹的脆性断裂

尖锐裂纹的变分问题已经通过相场模型得到了有效解决，其中引入了一个相场变量$\varphi$ 用来将尖锐裂纹近似为弥散裂纹相场变量$\varphi$ 在完全断裂区域等于1，在完好区域等于0，而$\varphi$ 在0到1之间表示从完好区域到完全断裂区域的过渡。

对于一维问题，材料的弥散裂纹可以近似表示为：

$\varphi \left(x\right)=\exp \left(-\frac{\left|x\right|}{a_c}\right)$ (9)

其中，$x=0$ 是裂纹位置，$a_c$ 表示长度尺寸参数，用于控制断裂材料与完整材料之间的过渡带宽度。

借助相场变量$\varphi$ ，尖锐裂纹面积由如下体积积分代替：

$A_s\left(\varphi \right)={\displaystyle {\int }_{\Gamma }\text{d}A}\approx {\displaystyle {\int }_{{\Gamma }_s}\gamma \left(\varphi ,\nabla \varphi \right)\text{d}V}$ (10)

其中，$A_s$ 表示扩散裂纹面积，$\gamma \left(\varphi ,\nabla \varphi \right)$ 裂纹为相场变量$\varphi$ 及其变量梯度$\nabla \varphi$ 的函数，可以定义为：

$\gamma \left(\varphi ,\nabla \varphi \right)=\frac{{\displaystyle {\int }_{\Omega }\left({\varphi }^2+a^2\nabla \varphi \cdot \nabla \varphi \right)}}{2a}$ (11)

其中$a$ 是相场正则化长度尺寸参数，即裂纹演化过程中弥散区域的带宽。

2.4. 导热材料的脆性断裂

导热材料的总能量由应变能、热自由能、断裂能和外部功组成，可以写为：

$\psi ={\psi }_1-{\psi }_2+{\psi }_3-{\psi }_4$ (12)

其中${\psi }_1$ 为弹性应变能，${\psi }_2$ 为热自由能，${\psi }_{\text{3}}$ 断裂能，${\psi }_{\text{4}}$ 为外力功。

导热材料的弹性应变能，热自由能，断裂能，外力功分别可以表示为：

$\{\begin{array}{l}{\psi }_1={\displaystyle {\int }_{\Omega }{\phi }_1\left(\epsilon ,M\right)dV}\\ {\psi }_2={\displaystyle {\int }_{\Omega }{\phi }_{\text{2}}\left(\epsilon ,M\right)dV}\\ {\psi }_{\text{3}}={\displaystyle {\int }_{\Gamma }{G}_{c}dA}\approx {\displaystyle {\int }_{\Omega }{G}_c\gamma \left(\varphi ,\nabla \varphi \right)\text{d}V}={\displaystyle {\int }_{\Omega }{G}_c\left[\frac{{\varphi }^2}{2a}+\frac{\varphi {\left(\nabla \varphi \right)}^2}{2}\right]dV}\\ {\psi }_4={\displaystyle {\int }_{\partial {\Omega }_t}h\cdot udA}-{\displaystyle {\int }_{\partial {\Omega }_m}QtdA}\end{array}$ (13)

其中，$G_c$ 为临界能量释放速率，外力功则是通过施加在边界$\partial {\Omega }_t$ 和$\partial {\Omega }_m$ 上的牵引力$h$ 和表面热流载流子$Q$ 表示。在相场模型中，弹性应变能密度函数${\phi }_1\left(\epsilon ,M\right)$ ，热自由能密度函数${\phi }_2\left(\epsilon ,M\right)$ 可以重写为：

$\{\begin{array}{l}{\phi }_1\left(\epsilon ,M\right)=\frac{1}{2}g\left(\varphi \right)\left({\sigma }_n:\epsilon \right)\\ {\phi }_{\text{2}}\left(\epsilon ,M\right)=\frac{1}{2}g\left(\varphi \right)\left(q\cdot M\right)\end{array}$ (14)

其中，$g\left(\varphi \right)={\left(1-\varphi \right)}^2$ 是相场模型中的退化函数。

将方程(5)带入到方程(13)中。有效弹性应变能密度函数可以表示为：

${\phi }_{1-2}\left(\epsilon ,M\right)=\frac{\left({\sigma }_n:\epsilon \right)f\left(\varphi \right)\cdot g\left(\varphi \right)}{2}$ (15)

其中，$P\left(\varphi \right)=f\left(\varphi \right)\cdot g\left(\varphi \right)$ 可以认为是相场损伤模型中的退化函数。代入$f\left(s\right),g\left(s\right)$ 可以得到：

$P\left(\varphi \right)=f\left(\varphi \right)\cdot g\left(\varphi \right)=k+\frac{1-\varphi }{1-D_0}$ (16)

其中，$k$ 是在模拟过程中用于避免数值奇异性的一个非常小的参数。

将上述方程带入总能量方程中，可以得到：

$\begin{array}{l}\psi =\frac{1}{2}{\displaystyle {\int }_{\Omega }p\left(\varphi \right)\left[\left({\sigma }_n:\epsilon \right)-\left(q\cdot E\right)\right]dV}+{\displaystyle {\int }_{\Omega }{G}_c\left[\frac{\left({\varphi }^2+a_c^2\nabla {\varphi }^2\right)}{2a_c}\right]\text{d}V}-{\displaystyle {\int }_{\partial {\Omega }_t}h\cdot udA}+{\displaystyle {\int }_{\partial {\Omega }_m}M\varphi dA}\\ \end{array}$ (17)

根据Griffith断裂理论，当裂纹开始扩展时，系统的总势能达到极值，此时总能量泛函的变分为：

$\begin{array}{c}\delta \psi ={\displaystyle {\int }_{\partial {\Omega }_t}\left[p\left(\varphi \right)\cdot {\sigma }_n\cdot n-h\right]\cdot \delta udA}+{\displaystyle {\int }_{\partial {\Omega }_m}\left[p\left(\varphi \right)\cdot D\cdot n+M\right]\delta \varphi dA}\\ +{\displaystyle {\int }_{\partial \Omega }{G}_c{a}_c\nabla \varphi \cdot n\delta \varphi dA}-{\displaystyle {\int }_{\Omega }\nabla \cdot p\left(\varphi \right)\cdot {\sigma }_n\cdot \delta udV}\\ -{\displaystyle {\int }_{\Omega }\nabla \cdot p\left(\varphi \right)\cdot D\cdot \delta TdV}-{\displaystyle {\int }_{\Omega }\left\{G_c{a}_c{\nabla }^2\varphi -{\psi }_5-{\psi }_6-\frac{G_c}{a_c}\varphi \right\}\delta \varphi dV}\end{array}$ (18)

其中$n$ 为单位法向量。${\psi }_5$ 和${\psi }_6$ 可以表示为：

$\{\begin{array}{l}{\psi }_5=\frac{\epsilon :\sigma }{2}\\ {\psi }_6=\frac{K\cdot M^2}{2}\end{array}$ (19)

条件中的$\delta \psi =0$ 应该对任意的$\delta u$ ，$\delta T$ 和$\delta \varphi$ 都成立，因此可以得到温度场和相场的控制方程如下：

$\{\begin{array}{l}\nabla \cdot {\sigma }_n\cdot p\left(\varphi \right)=0,\text{in}\Omega \\ \nabla \cdot D\cdot p\left(\varphi \right)=0,\text{in}\Omega \\ G_c{a}_c{\nabla }^2\varphi -\frac{{\psi }_5+{\psi }_6}{1-D_0}-\frac{G_c}{a_c}\varphi =0,\text{in}\Omega \end{array}$ (20)

相场，温度场的边界条件如下：

$\{\begin{array}{l}p\left(\varphi \right)\cdot {\sigma }_n\cdot n=h,\text{in}\partial {\Omega }_t\\ p\left(\varphi \right)\cdot D\cdot n=-Q,\text{in}\partial {\Omega }_m\\ \nabla \varphi \cdot n=0,\text{in}\Omega \end{array}$ (21)

为了确保相场变量$\varphi$ 在损伤演化过程中的不可逆性，Miehe [23]等人引入了一个历史变量场$H\left(x,t\right)$ ：

$H\left(x,t\right)=\max {\psi }_c^{+}\left(x,t\right)$ (22)

因此，可以重写相场演化方程：

$G_c\left[\frac{1}{a_c}\varphi -a_c{\nabla }^2\varphi \right]-\frac{H\left(1-\varphi \right)}{1-D_0}=0$ (23)

3. 数值结果分析

3.1. 纯弹性材料相场损伤模型的验证

本小节采用ABAQUS软件的UEL子程序对含初始损伤的纯弹性材料进行裂纹扩展分析。在ABAQUS中设定的边界条件如图2所示，方形平板底部边界沿x、y方向的位移被固定，位移步长Δu设为1 × 10⁻⁷ mm。材料弹性模量E = 210 GPa，泊松比ν = 0.3，长度尺度参数a_c = 0.0075 mm，临界能量释放率G_c = 2.7 N/m。节点力的计算结果如图3所示，与Miehe [23]和Xingxue Lu [24]等人在类似模型中获得的结果一致。考虑到初始损伤的存在，节点处的应力值较Miehe [23]和Xingxue Lu [24]等人的计算结果更高，这与损伤力学的基本原理相吻合。

Figure 2. Boundary conditions of pure elastic material

图2. 纯弹性材料的边界条件

3.2. 相场损伤模型与传统模型的对比验证

将该模型退化为导热陶瓷的相场断裂模型，并采用如图4所示的单边缺口裂纹的导热陶瓷正方形板进行验证测试。单边缺口裂纹正方形板的尺寸参数和边界条件和纯弹性材料案例相同。裂纹沿板的边缘施加，并在热流主方向施加均匀温度场。模拟材料选用高导热陶瓷材料，材料属性参数如表1所示。顶部边界温度设定为T₁，T₂沿x轴方向施加，材料发生均匀形变，材料弹性模量E = 320 GPa，泊松比ν = 0.2，临界能量释放率G_c = 3 N/m。边界条件要求顶部施加温度场T₂，初始温度边界条件保持恒定为T₁ = 300℃。

Figure 3. Support reaction-displacement curve of square plate with one-sided notch

图3. 单边缺口正方形板的支反力–位移曲线

Figure 4. Boundary conditions for a square plate with a single-sided notch (a) Finite element meshing (b)

图4. 单边缺口导热方形板的边界条件(a)有限元网格划分(b)

Table 1. Material properties of high thermal conductivity ceramics (AlN)

表1. 高导热陶瓷(AlN)的材料特性

属性	值
弹性常数(GPa)	c11 = 333	c12 = 66.6	c21 = 66.6
	C22 = 333	C33 = 133
初始温度(℃)	T1 = 300
导热系数(W/m℃)	K = 300
长度尺寸参数(mm)	ac = 0.01

仿真结果表明，初始损伤的存在对MEMS器件制造过程中导热薄膜的裂纹扩展行为具有显著影响。材料内部预先存在的缺陷会显著增强局部应力集中，从而改变系统的整体力学响应特性。初始损伤明显增加了裂纹萌生及后续扩展的敏感性。由此可见，含有预存损伤的MEMS器件在使用过程中更容易发生疲劳引发的裂纹扩展现象。为研究不同温度对导热陶瓷中裂纹扩展的影响，本实验在材料底部边界施加固定位移Δu = 1 × 10⁻⁷ mm毫米，并维持T₁ = 300℃的恒定温度场。顶部边界分别施加T₂ = 350℃，300℃，250℃三种温度。通过二次开发UEL子程序计算得到的相场云图如图5所示。

(a) ΔT = −50℃ (b) ΔT = 0℃ (c) ΔT = 50℃

Figure 5. Crack propagation cloud map under different temperatures

图5. 施加不同温度下的裂纹扩展云图

(a) t = 50 μs (b) t = 65 μs (c) t = 80 μs

Figure 6. Temperature field distribution at different time points at 300˚C

图6. 300℃下不同时刻温度场分布图

Figure 7. Thermal stress distribution map

图7. 热应力分布图

图6，图7分别为不同时刻的温度场分布图和热应力分布图，图8展示了样品在裂纹扩展与边缘载荷边界条件下的位移曲线与Ved Prakash [25]等人的结果对比。可以看出：当热流方向与y轴方向反向时，节点支撑力较小；而当热流方向与y轴方向同向时，节点支撑力较大。由此可得出结论：负温差会加速材料内部的裂纹扩展，而正温差则会延缓其裂纹形成。这一现象表明，在MEMS器件制造过程中，通过调控施加在材料上的温度，操作人员能够有效抑制材料的裂纹行为。

Figure 8. The displacement-reaction curves of a single-sided notch specimen under different temperatures

图8. 展示了单边缺口试样在不同温度作用下的位移–支反力曲线

4. 结论

本文开发了一种相场损伤模型，用于分析热力耦合条件下MEMS导热薄膜中的裂纹扩展。通过将初始微缺陷纳入模型，该方法能有效捕捉制造缺陷对断裂行为的影响，从而得到以下结论：(1) 初始微缺陷对位移约束边界处的节点支反力具有显著影响，这反映了材料承载能力的变化。当材料具有加工缺陷引起的初始损伤与未考虑初始损伤的情况对比分析，节点处的应力值更大。这种下降的情况表明，预先存在的缺陷通过局部降低有效刚度，加速了裂纹的萌生与扩展。根据相场理论，较大的初始损伤会促使应变能更早局域化，从而削弱MEMS导热陶瓷的抗断裂能力。(2) 外加温度场(T₂)对裂纹扩展具有显著影响。当施加时T₂ = 250℃，支反力达到121.2 KN，比T₂ = 300℃和T₂ = 350℃条件下的118.3 KN和117.2 KN更高。这些结果表明热流方向与材料热膨胀方向相同时能有效抑制裂纹扩展，而反向则会促进裂纹扩展。

NOTES

^*第一作者。

^#通讯作者。