把一束100飞秒的脉冲激光打到金属表面时间短到什么程度光在这段时间里只走30微米。这么短的时间里能量先被自由电子“接住”晶格还没来得及振动电子温度已经冲高到上万K——这就是飞秒激光和纳秒激光最本质的区别也是为什么经典傅里叶热传导方程在飞秒烧蚀模拟里基本失效的原因。做这类仿真的同学最终基本都会绕回同一个组合用Comsol求解双温方程再叠加热应力做热力耦合去预测烧蚀深度、残余应力和加工形貌。这篇博文不打算讲教科书式的理论我更想把“双温方程”“热力耦合”“烧蚀”这些概念放到实际建模过程里逐个拆开重点放在真正能从0搭起一个可运行模型的细节上方程怎么配置、热源怎么写、边界怎么设、网格和求解器怎么调。适合正在用Comsol模拟超快激光微加工的科研党、研究生以及做激光清洗、激光微结构制备的工程人员。1. 飞秒激光烧蚀为什么不能用常规传热方程硬算1.1 飞秒激光的“非平衡”世界先看一个时间尺度的对比。纳秒激光的脉冲宽度是10⁻⁹秒量级这个时间足以让电子和晶格完成多次能量交换激光作用期间材料基本处于局部热平衡状态经典热传导方程加一个表面热源就能算个大概。飞秒激光的脉冲宽度是10⁻¹³秒量级100飞秒意味着光在材料里只走了约30微米而电子与声子之间完成能量交换的特征时间通常在1到10皮秒量级。这意味着什么脉冲已经结束了晶格温度还基本没变化能量全部先沉积在自由电子系统里。电子通过逆韧致辐射吸收光子能量后靠电子-电子散射在几十飞秒到几百飞秒内先自己达到一个准平衡形成一个极高的电子温度T_e。随后电子再通过电子-声子耦合把能量慢慢交给晶格让晶格温度T_i在后脉冲阶段逐渐升高。这个过程中电子系统和晶格系统各自拥有完全不同的温度用一个统一的“材料温度”去描述本身就是错的。所以在飞秒烧蚀仿真里必须引入双温模型把能量在电子和晶格两个子系统之间的非平衡输运过程显式地解出来。这也是整个Comsol模型的物理基础。1.2 双温方程长什么样、每一项代表什么双温模型的基本形式是两条耦合的非线性抛物型偏微分方程。我在自己的模型里用的是经典公式C_e(T_e) * ∂T_e/∂t ∇·(k_e(T_e,T_i)·∇T_e) - G*(T_e - T_i) Q(z,t) C_i * ∂T_i/∂t ∇·(k_i·∇T_i) G*(T_e - T_i)第一条描述电子系统的能量守恒。左边是电子热容C_e乘电子温度的瞬态项右边第一项是电子热传导项第二项是电子与晶格之间的能量交换项G是电子-声子耦合系数最后Q是激光体热源。第二条描述晶格系统。晶格热容C_i通常取常数k_i是晶格热导率它和G项从电子系统接过来的能量一起决定晶格温度上升。双温方程里最常见的坑是电子热容和电子热导率都不是常数。自由电子气的电子热容近似为C_e γ·T_e其中γ是索末菲系数单位是J/(m³·K²)这对于金属在常温以上都是足够好的近似。电子热导率k_e则强烈依赖T_e和T_i工程上常用k_e k_eq·T_e/T_i来近似其中k_eq是平衡态热导率。这样处理的好处是物理上直观电子温度越高电子输运能量能力越强晶格温度越高电子-声子散射加剧导热能下降。激光热源Q(z,t)也不是随便给一个边界热流就完事。飞秒激光在金属里的光学穿透深度只有十几纳米所以通常把它当成一个沿深度指数衰减的体热源Q sqrt(4*ln(2)/pi) * F*(1-R)/(τ_p*δ) * exp(-4*ln(2)*((t-t_peak)/τ_p)^2) * exp(-z/δ)其中F是激光能量密度通量R是表面反射率τ_p是脉冲宽度δ是光学穿透深度z是距表面深度t_peak是脉冲峰值时刻。前面那个sqrt(4ln2/π)是高斯脉冲的峰值归一化因子保证这个源项对时间积分后正好等于材料吸收的总能量密度F(1-R)/δ。很多刚开始做的人会把归一化因子丢掉结果能量偏大或者偏小后面我会详细说。1.3 热力耦合为什么非加不可光算温度场其实不能满足工程需求。飞秒激光在极短时间内把巨大能量灌入材料表面极薄一层瞬间形成的温度梯度可高达10⁶-10⁷ K/m量级这个温度梯度会产生瞬态热应力应力峰值很容易超过材料的动态屈服强度或抗拉强度。实验上观察到的很多现象比如金属表面的层裂剥落、微裂纹、周期性条纹、材料以固态碎片形式被喷出都跟热力学应力直接有关单纯用“温度超过汽化点就烧蚀”这种判据根本解释不了。我在模型里加热力耦合核心思路很简单把双温方程解出来的晶格温度T_i当作温度载荷输入到固体力学模块里。固体材料产生热应变ε_th α(T_i)·(T_i - T_ref)其中α是热膨胀系数T_ref是零应力参考温度。再结合线弹性本构关系σ C·(ε - ε_th)就能得到瞬态应力场。有了应力场分布烧蚀判据就有更多选择了可以是晶格温度超过临界温度也可以是最大主应力超过材料动态拉伸强度或者两者共同作用。这里有个关键点电子温度T_e直接产生的力学效应其实很小因为力是通过原子/晶格传递的晶格还没热起来应力自然不会以电子温度为基础。所以我在做热力耦合时温度载荷必须用T_i而不能误用T_e这个细节我在审稿和帮同学看模型时碰到过好几次。2. Comsol建模方案如何把双温方程变成可算的有限元模型2.1 为什么不用内置“固体传热”模块直接算有人第一次做会想Comsol里有现成的“固体传热”接口能不能直接拿来算飞秒激光答案是不行。内置传热模块只有一条温度场方程对应的是傅里叶热传导它默认整个材料只有一个温度。双温模型需要两条相互耦合的温度方程电子系统一条晶格系统一条这超出了内置传热模块的表达范围。我的做法是绕开传热模块直接用数学模块里的“系数型PDE”接口来写双温方程。Comsol的“系数型PDE”提供了一个非常通用的模板e_a * ∂²u/∂t² d_a * ∂u/∂t ∇·(-c∇u - αu γ) β·∇u a·u f对于双温模型这种抛物型方程e_a、α、γ、β都可以设为零剩下d_a对应热容项c对应热导率f对应源项形式非常干净。我建立两个系数型PDE接口一个的因变量设为T_e另一个设为T_i然后在系数里引用对方变量和材料参数就完成了两个温度场的双向耦合。这里有个优点系数型PDE里所有系数都允许是变量表达式。比如我可以直接在“d_a”栏写gamma*T_e在“c”栏写k_eq*T_e/T_i在“f”栏写-G*(T_e-T_i)Q。Comsol在每次迭代时都会自动更新这些系数非线性双温耦合就能被直接求解不需要自己写额外的迭代循环。2.2 几何与模型域的选择几何处理是很多人第一步就卡住的地方。我通常用二维轴对称模型因为单脉冲垂直入射、光斑为高斯圆形时温度场和应力场关于光束轴线对称二维轴对称模型能把三维问题压成二维计算量下降一到两个数量级而且精度完全够用。模型域的尺寸选择很讲究。径向方向至少要大于激光光斑半径的5到10倍这样才能忽略侧面边界对热扩散的影响。深度方向要看两个尺度一是光学穿透深度δ可能只有十几纳米这是温度梯度最剧烈的区域二是热扩散和应力波传播的距离在几纳秒时间尺度内热能传播深度大约数百纳米到一微米应力波可能传播到几微米甚至更远。我一般取深度为5到10微米底部设为绝热和固定约束的组合边界这样既能覆盖热影响区又不会让应力波在边界反射过强干扰目标区域结果。薄膜结构需要单独考虑。如果做的是几十纳米厚的金属薄膜几何可以直接建成分层结构每层用不同的材料参数对多层膜在PDE接口中做材料属性映射即可。如果是半无限大块体可以在底部加一层完美匹配层或者足够厚的阻尼扩展区来吸收应力波不过这会让模型复杂不少初期不建议一上来就加。2.3 材料参数从哪来、误差有多大双温模型最麻烦的其实不是方程本身而是材料参数。不同文献给出的参数可能相差数倍而且很多参数本身具有温度依赖性和尺寸效应。我做模型时常用以下几个典型值来校准。参数符号典型值金典型值铜备注索末菲系数γ68 J/(m³·K²)96 J/(m³·K²)决定电子热容电子-声子耦合系数G2.6e16 W/(m³·K)1.0e17 W/(m³·K)文献差异较大平衡电子热导率k_eq317 W/(m·K)401 W/(m·K)用k_eq*T_e/T_i近似晶格热容C_i2.49e6 J/(m³·K)3.45e6 J/(m³·K)常取常数晶格热导率k_i约2-5 W/(m·K)约3-5 W/(m·K)通常远小于电子贡献光学穿透深度δ约15 nm约14 nm与激光波长有关反射率R0.93近红外0.97近红外随温度会变化看到没有G值不同文献里金可以从1e16到5e16差出5倍。这直接导致晶格温度峰值和烧蚀深度的预测结果明显不同。我做模型时的方法是先固定一组参数跑通流程再通过实验测量烧蚀深度或表面损伤阈值来反标定G、δ和R这三个最敏感参数。没有实验数据的纯理论预测最好在论文里明确说明参数参考来源并做敏感性分析。3. 一步步搭模型从参数定义到求解器设置3.1 全局参数和激光热源表达式在Comsol的“全局定义→参数”里我习惯把所有物理量统一写成带单位的参数这样Comsol会自动做单位检查能避免很多低级错误。tau_p 100[fs] F 0.2[J/cm^2] t_peak 5*tau_p R 0.93 delta 15[nm] r_0 20[um] // 激光光斑半径 G 2.6e16[W/(m^3*K)] gamma 68[J/(m^3*K^2)] C_l 2.49e6[J/(m^3*K)] k_eq 317[W/(m*K)] k_l 3[W/(m*K)] T_0 300[K] T_abl 8000[K] // 烧蚀判据温度按材料调节注意我在这里统一采用国际单位制时间用秒长度用米温度用开尔文。τ_p写成100[fs]没问题Comsol会换算成1e-13秒参与计算t表达式里的时间是秒这个事实不会变。激光热源我建议先定义一个“变量”Q_laser写在“定义→变量”里方便后面引用和检查Q_laser sqrt(4*ln(2)/pi)*F*(1-R)/(tau_p*delta) * exp(-4*ln(2)*(t-t_peak)^2/tau_p^2) * exp(-z/delta) * exp(-r^2/r_0^2)公式里加了径向高斯分布项exp(-r²/r₀²)r是径向坐标。有些模型里还会把这个高斯项归一化让总功率不随光斑大小变化具体看你是固定能量密度还是固定总能量。我大部分时候固定能量密度F所以不额外乘归一化系数。如果做多脉冲只需要把Q_laser写成好几个高斯峰叠加每个峰对应不同t_peak和F。3.2 用系数型PDE接口写双温方程实际操作路径是模型向导里添加“数学→PDE接口→系数型PDE”几何选二维轴对称。第一个接口因变量名改成Te第二个接口因变量名改成Ti。第一个PDE接口电子温度的系数这样填系数表达式d_agamma*Teck_eq*Te/Tia0f-G*(Te-Ti) Q_laser第二个PDE接口晶格温度的系数系数表达式d_aC_lck_la0fG*(Te-Ti)这样两个方程之间的耦合关系就建立起来了。我在第一次跑模型时犯过一个错误把Te和Ti的初始条件都设成了0 K。这会导致d_a里的gamma*Te在初始时刻为0方程出现退化求解器直接提示“初始值导致奇异矩阵”或者干脆不收敛。正确做法是全部初始化为300 K让电子热容从初始时刻就有一个合理的正值。边界条件方面所有外边界我都用默认的“零通量”条件。物理上这意味着表面没有额外的能量散失这在飞秒到亚纳秒的时间尺度内是合理的因为热辐射和空气对流在这个时间窗口内还没来得及带走多少能量。如果算的是脉冲串且时间跨度达到微秒以上就需要考虑表面热辐射和对流了不过我建议初版模型先不要加。3.3 热力耦合设置把晶格温度作为温度载荷在模型里加入“固体力学”接口并把材料模型设置为“线弹性材料”然后在“热膨胀”子节点下设置参考温度 T_ref T_0 温度场 T Ti这个设置意味着固体力学模块会读入由第二个PDE接口算出的晶格温度场Ti计算热应变。这里要特别提醒温度场必须选Ti不是Te。我自己第一次做的时候就因为图省事直接用了全局变量Te结果应力峰值大得离谱后来一查才发现热应变来自晶格温度不是电子温度。约束条件要避免刚体位移。如果模型是二维轴对称轴对称边自动约束径向位移我通常在底面或侧面选一个点做“固定约束”或者使用“弱弹簧”来抑制刚体运动。对于自由表面烧蚀的问题表面不要加任何额外约束让它能自由膨胀变形这样才能得到真实的拉压应力分布。为了控制计算量热力耦合可以做成准静态假设。飞秒激光作用后晶格温度变化时间尺度是皮秒级而弹性波在金属里传播速度快于纳米/皮秒比热扩散快得多。很多文献在处理这个问题时会把每个时间步的晶格温度场视为准静态温度载荷用静力学求解器求解固体力学方程。我在Comsol里是用“瞬态”研究但把固体力学接口的方程形式设为“准静态”物理上对应应力波瞬时传播的假设算出来结果在ps到ns时间尺度上误差不大却可以省掉极小时间步长要求的隐性约束。3.4 烧蚀边界动态后移移动网格怎么加要不要实际模拟材料去除这是模型复杂度的一个分水岭。如果只是想初步看看温度分布和应力分布完全可以不移动边界后处理阶段取晶格温度超过烧蚀阈值的等值面深度作为“有效烧蚀深度”这样最简单也最容易收敛。如果一定要看烧蚀形貌的动态演化就得用“变形几何”接口。做法是给烧蚀表面设置一个法向速度边界条件让表面在晶格温度超过阈值时向材料内部后退。我用的速度表达式是v_n v_max * flc2hs(Ti - T_abl, dT)其中v_max是表面后退的特征速度dT是阶跃平滑宽度flc2hs是Comsol内置的平滑Heaviside函数。直接写成if(TiT_abl, v_max, 0)会出问题因为阶跃跳变会导致雅可比矩阵突变求解器很容易振荡甚至发散。用flc2hs做一个10到100 K宽度的平滑过渡收敛性会好很多。v_max怎么取严格来说需要根据蒸发速率、流体动力学甚至分子动力学来确定但在简化模型里我通常取材料声速的千分之一到百分之一量级再根据实验测量的烧蚀深度反标定。这里要清醒这种移动网格模型本质上是唯象模型定量精度取决于v_max的取值但它能很好地帮助理解表面后退的时空过程。移动网格还要设置“自动重新剖分”。随着表面不断后退变形几何的网格会被逐渐拉长雅可比行列式可能出现负值求解器报“网格扭曲”错误。我在“变形几何”接口里开启“自动重新剖分网格”触发条件设为网格质量低于0.3左右这样Comsol会在计算中途自动重建网格并插值前一步的解虽然会增加计算时间但稳定性提升非常明显。3.5 网格划分与求解器调试网格是这类强瞬态多物理场问题的关键。表面附近必须捕捉光学穿透深度δ的衰减而δ只有十几纳米所以从表面开始往内部要布置非常薄的边界层网格。我的做法是表面第一层网格厚度控制在δ/10以内大约1到2纳米。用8到10层边界层网格以1.3到1.5的比例逐层增厚覆盖到约1微米深度。其余区域用自由三角形网格最大单元尺寸在1到2微米。径向方向在光斑边缘rr₀附近加密因为该处温度梯度和应力梯度也很大。全模型自由度通常在几万到十几万之间二维轴对称下这个规模在个人工作站上可以接受。求解器设置我强烈建议分两步走。第一步只求解Te和Ti两个PDE接口关闭固体力学时间范围设成0到1纳秒时间步进用BDF算法初始步长设为τ_p/10最大阶数5相对容差1e-4。等双温场跑通、结果合理了再打开固体力学接口用“分离式”求解器分别求解温度和位移避免全耦合牛顿迭代因为尺度差过大而难以收敛。时间步长的控制上一个实用技巧是设置一个“随时间变化的允许最大步长”在脉冲作用前后这段剧烈变化期强制用小步长在脉冲结束后的皮秒到纳秒阶段自动放大步长。具体方法是添加“时间步长”节点指定中间时间点比如0、t_peak、t_peak10τ_p、t_peak100τ_p、1ns等Comsol会自动在这些点附近刷新步长。我实测下来这样的时间推进效率比纯自适应快很多而且结果稳定。4. 结果后处理从云图到烧蚀深度怎么分析4.1 电子温度和晶格温度的时空演化模型算完后我最先看的是电子温度云图和晶格温度云图。在“二维绘图组”里选表面表达式分别填Te和Ti时间点选几个代表性的脉冲峰值时、1 ps后、10 ps后、100 ps后。你会看到典型的非平衡特征在1 ps时间点电子温度已经向材料内部扩散了几十纳米而晶格温度还只限于表面极薄层到10 ps后电子温度大幅下降晶格温度才逐渐升高并开始向内部传导。更好的定量分析是在轴对称轴上画一维线图。我通常在“派生值”里定义一条从表面沿z轴的截线选择“一维绘图组”x轴设为弧长即深度y轴同时画Te和Ti时间点选4到5个。这样可以清晰看到电子热扩散前沿远超光学吸收深度这个现象叫“热电子扩散”是飞秒激光与金属相互作用的一个标志性特征。如果模型里k_e取了常数你就看不到这种扩散反过来如果你发现电子温度峰高得离谱先检查是不是k_eq*T_e/T_i在低温段数值异常大。4.2 热应力分布与层裂判据应力结果方面我常画三个量最大主应力、von Mises应力和径向应力。在烧蚀初期的百飞秒到几皮秒内表面薄层受热膨胀但又被周围冷材料限制表现为很强的压应力脉冲结束之后随着晶格温度继续升高和应力波向内部传播表层会转入拉应力状态。这个压-拉转换过程是理解层裂和表面剥落的关键。我在后处理里会专门提取某个深度点比如表面下50 nm的应力随时间变化曲线看拉应力峰值出现的时刻和量级。如果最大拉应力超过材料动态拉伸强度比如金属的GPa量级就可以判断该区域可能发生机械性层裂。很多人只做热分析不做应力分析结果实验里看到的一些碎片喷溅和表面剥落完全没法解释这就是热力耦合模型的价值所在。4.3 烧蚀深度与烧蚀阈值通量的关系飞秒激光烧蚀曲线实验上常表现为烧蚀深度随通量对数线性增加可以用简化公式描述d ≈ δ_eff * ln(F / F_th)其中δ_eff是有效能量沉积深度F_th是烧蚀阈值通量。我在Comsol里验证这个关系时直接用“参数化扫描”扫描F从0.05到0.5 J/cm²对每个F分别计算晶格温度等于烧蚀阈值T_abl的最大深度画成ln(F) vs d的曲线。如果材料参数合理得到的是一条很好的直线斜率为δ_eff横截距对应F_th。这个拟合结果可以用来和实验对比误差在30%以内就算相当不错了。使用移动网格模型时烧蚀深度直接看边界位移量就行。没有移动网格时用后处理提取等值面距离本质上是一种静态判据没有考虑材料去除后对后续热传导的影响所以通量较高时误差会偏大。前期的参数扫描用静态判据足够最终精确结果再跑移动网格版本。5. 我踩过的坑收敛失败、能量不守恒、网格崩坏5.1 电子热容在低温启动时变得“奇异”我第一次把初始温度设为0 K计算直接崩了。原因很明确d_a系数里包含gamma*TeTe0时d_a0偏微分方程退化成无时间项数值上就会出现奇异矩阵。解决方案很简单把Te和Ti的初始值都设为300 K。如果你需要模拟极低温环境至少也要给一个接近环境温度的正值或者给C_e加一个小下限C_e max(gamma*Te, 10[J/(m^3*K)])这样既不会改变高温段的物理也避免了低温电子热容趋近于0的数值问题。5.2 单位制混用导致的热源能量翻车这是另一个高发问题。有人把τ_p写成100然后t表达式中直接用t单位秒结果高斯脉冲宽度等于100秒整个能量在巨长的时间里“均匀”释放表面温度完全不是飞秒脉冲该有的样子。解决办法一是参数全部带单位让Comsol自动转换二是在求解前随手在“派生值”里做一次能量积分校验。具体做法计算Q_laser对体积和时间积分结果应该等于F*(1-R)乘上光斑面积在没有归一化高斯项的情况下。如果积分值与预期差了几个数量级那一定是单位换算或者归一化因子写错了。5.3 全耦合求解器为啥老是算不下去双温PDE刚加固体力学时如果直接用全耦合牛顿迭代十有八九在脉冲初期就会因为温度梯度过大而失败。我的经验是先跑纯温度场确认Te和Ti结果合理再打开固体力学求解器类型改成“分离式”“分离式”里让Te和Ti作为一个组先迭代位移作为另一个组后迭代。分离式在强非线性问题上收敛稳健度远好于全耦合代价只是多几次迭代。另外热源Q_laser里如果写了阶跃函数导致Q瞬间从0跳到很大可以尝试把阶跃用flc2hs平滑替换也能明显改善收敛表现。5.4 移动网格“网格扭曲”的处置用变形几何做烧蚀后退时网格翻转几乎是必然遇到的。原因通常是表面法向速度过大或者平滑Heaviside参数dT太小导致局部变形梯度过于剧烈。我的排查套路是现象可能的根因调整方向表面网格严重拉长v_max设置过大把v_max降低一个数量级局部网格扭成负面积dT太小阶跃太陡把dT从1 K增大到50-100 K深层网格大面积扭曲变形几何覆盖区域过大只在表面附近1-2微米范围启用变形中途雅可比报错长时间变形累计过大开启“自动重新剖分网格”变形几何也不是必选项。如果目标只是算烧蚀深度、阈值通量和应力趋势用静态温度阈值判据完全足够而且结果可以复现、参数扫描极快。动态网格适合发论文出效果图不适合大量参数扫描。6. 这个模型还能往哪些方向扩展6.1 双脉冲和脉冲串烧蚀双温方程天然支持多脉冲只要把Q_laser改成多个高斯峰叠加即可每个峰对应不同的能量密度到达时间。做泵浦-探测研究的同学可以通过这个模型研究双脉冲时间延迟对电子温度峰值和烧蚀效率的影响实验上这对应着“预热-主加热”脉冲串工艺优化。要注意的是多脉冲计算到纳秒甚至微秒时间尺度晶格热积累效应会越来越明显表面散热边界不能再简单设零通量。6.2 多层薄膜和异质界面把几何建成多层结构后各层材料参数不同PDE接口里通过“材料”节点自动映射参数双温方程和热应力都逐层适用。界面处要特别关注电子热导率和声子耦合系数在界面处可能不连续很多文献还会在界面加界面热阻这个在Comsol里可以通过边界条件或薄层近似来实现。对于半导体上镀金属膜这类结构双温方程的适用性要谨慎评估因为半导体中的载流子动力学不止光吸收后电子弛豫这么简单。6.3 光束整形与偏振效应如果激光光斑不是简单的高斯圆斑而是平顶光束、涡旋光束或者线偏振光导致的强度分布可以通过修改Q_laser的空间分布项来模拟。二维轴对称模型只适用于圆对称分布非对称情况需要换成三维模型计算量会上升一到两个数量级。这个扩展对激光微结构加工很有价值能直接比较不同光束模式下烧蚀形貌的差异。6.4 与烧蚀产物动力学衔接双温模型止步于材料去除的条件判断至于去除后的材料是以蒸气、液滴还是固体碎片形式离开表面双温模型自己解释不了。后续如果再往上走一般有两个方向一是把去除层当作热源/动量源耦合流体接口计算羽流膨胀二是和分子动力学结合把PIC/原子尺度方法放在表面薄层宏观区域继续用双温连续介质模型。这两个方向都是研究级的大工程不适合入门阶段碰。我个人在实际操作中的体会是飞秒激光烧蚀仿真最大的不确定性永远在材料参数和烧蚀判据而不是有限元求解本身。双温方程和热力耦合的实现路径已经非常成熟Comsol只是把复杂的多物理场组合成一个“搭积木”的过程。建议所有新手先做一个一维模型验证物理认知再扩展成二维轴对称的完整模型而不是一上来就追求三维细节效果。只要把单位、热源、耦合边界、网格和时间步这五件事盯住这个模型跑出稳定、可复现的结果只是时间问题。