超快激光烧蚀成孔带有热应力的COMSOL模型，采用双PDE方程模拟双温以及热应力模块，动态图所...

超快激光烧蚀成孔带有热应力的COMSOL模型，采用双PDE方程模拟双温以及热应力模块，动态图所示是10ps下的结果，另外附有参考文献

超快激光加工中温度场和应力场的耦合效应一直是仿真难点。最近折腾了个基于COMSOL的双温模型热应力耦合方案，针对10皮秒量级脉冲的烧蚀成孔过程，这里分享些关键实现细节。

模型核心是双温方程与热弹性方程的联立求解。电子温度Te和晶格温度Tl的动态平衡过程用两个耦合的PDE描述：

% 电子温度方程 d(Te,t) - ∇·(k_e*∇Te) = -G*(Te - Tl) + Q_laser; % 晶格温度方程 d(Tl,t) - ∇·(k_l*∇Tl) = G*(Te - Tl);

这里G是电子-声子耦合系数，ke和kl分别为电子和晶格的导热系数。有趣的是在COMSOL中需要把这两个方程拆成两个独立的PDE模块——别直接使用内置传热模块，否则会丢失电子温度特有的非线性项。

材料参数设置要注意温度依赖特性。比如金属的电子热容C_e随温度变化显著，这里用了分段函数处理：

if Te < T_Fermi C_e = gamma*Te; else C_e = gamma*T_Fermi*(1 + (Te-T_Fermi)/(2*T_Fermi)); end

这种处理方式比简单线性假设更符合实验数据，但需要小心处理数值稳定性。建议先用稳态求解器确定初始温度场，再切换到瞬态求解。

超快激光烧蚀成孔带有热应力的COMSOL模型，采用双PDE方程模拟双温以及热应力模块，动态图所示是10ps下的结果，另外附有参考文献

热应力模块的难点在于应变率的匹配。当激光脉冲宽度在皮秒量级时，必须开启几何非线性选项，同时材料模型建议采用各向异性弹塑性模型。这里有个容易被忽略的参数——热膨胀系数的温度相关性：

alpha = alpha_0 + alpha_1*Tl + alpha_2*Tl^2;

二次项的存在会导致应力波传播方向发生偏转，这在后处理时需要特别关注。

动态图显示在10ps时，材料表面出现了明显的热弹力波（见附图红色波动区域）。有趣的是应力集中区并不与最高温度区完全重合，这与传统连续波激光的热应力分布特征明显不同。这种非平衡态现象验证了双温模型的必要性。

求解器配置方面，推荐使用分离式求解器并启用自动牛顿迭代。时间步长设置为0.1ps能满足精度要求，但要注意当电子温度梯度超过1e15 K/m时需要动态缩小步长。内存消耗大概在32GB左右，建议使用分布式求解加速计算。

参考文献可重点看Rethfeld的电子弛豫模型（Phys. Rev. B 2004）和孙教授的飞秒激光烧蚀理论（IJHMT 2019）。这些理论框架为参数设置提供了重要依据，特别是电子-声子耦合时间常数的取值逻辑。