今年上半年我接到一个很头疼的建模任务：用 COMSOL 6.1 搭一个三维飞秒多脉冲激光烧蚀玻璃模型，并且要把双温模型、变形几何和烧蚀形貌演化全做进去。刚开始我觉得这事不难——双温模型说白了就是电子温度和晶格温度两个偏微分方程，变形几何无非给表面设个移动速度；可真把三维网格、多脉冲循环和动网格叠到一起之后，我才发现这里面的坑远比想象中多。
飞秒脉冲的时间尺度是飞秒到皮秒，脉冲间隔是微秒量级，而网格变形又要跟着晶格温度一步步推进，三个物理过程在同一个三维模型里互相撕扯，计算量和收敛难度都是指数级上升。这篇文章是我把模型从零搭到能稳定跑完五个脉冲的完整复盘，包括物理机制、COMSOL 里的方程落法、边界速度怎么设、多脉冲怎么循环，以及我实际踩过的五个问题。如果你正在做飞秒激光与透明介质相互作用仿真，或者准备从二维单脉冲模型转向三维多脉冲，这篇记录应该能帮你省掉不少试错时间。
1. 模型需求拆解：三维、多脉冲、动网格三座大山的叠加效应
1.1 飞秒脉冲烧玻璃的物理链条
飞秒激光和玻璃相互作用的物理过程，用一句话概括是：脉冲极短，电子先热，晶格后热，晶格温度超过临界值时材料被去除。拆开看可以分成四步。第一步，激光在焦点附近通过多光子电离和雪崩电离产生自由电子；第二步，自由电子通过逆轫致辐射吸收激光能量，电子系统在不到一皮秒内迅速升温；第三步，电子系统和晶格系统通过电子-声子耦合交换能量，晶格温度在几皮秒到几十皮秒内升高；第四步，晶格温度超过材料的破坏阈值，发生相变、飞溅和材料去除。
这里最麻烦的是时间尺度差异。激光脉宽只有两百飞秒左右，电子-晶格弛豫时间是几皮秒，热量向周围扩散则要持续到纳秒甚至微秒。如果把这些时间尺度全部塞进一个三维模型，求解器会在脉冲作用瞬间遇到一个极陡的热源尖峰，然后在漫长的冷却段里又面临时间步长放大的诱惑。任何一个环节处理不好，要么算出来的电子温度尖峰严重失真，要么计算步数多到不现实。
1.2 为什么许多教程跑二维，三维是另一种玩法
很多双温烧蚀模型其实是二维轴对称或者单轴一维模型，因为单脉冲、正入射、圆光斑的问题天然具备轴对称性，二维模型就可以表达温度场和烧蚀轮廓。但标题里明确要求三维，说明要关注的不是理想化对称问题，而是实际的非对称焦点、扫描路径、多脉冲热累积带来的三维形貌。三维模型每个单元体积小得多，同样的几何域，二维可能只有几十万自由度，三维轻松破百万。
飞秒烧蚀问题又必须在焦点附近加密网格，远离焦点的区域则需要足够大的体量来吸收热扩散。我最终用的是一个 40×40×60 微米的玻璃长方体，焦点区网格尺寸取 0.3 微米，远处逐步放大到 2 微米，总自由度大约七十万。如果全程都用 0.3 微米网格，自由度会冲到三四百万，求解器基本没法在可接受时间内跑完。这个取舍不是偷懒，而是三维动网格问题的必然选择。
1.3 建模目标和参数基线
| 参数 | 取值 | 说明 |
|---|---|---|
| 激光波长 | 1030 nm | 常见飞秒激光器波长 |
| 脉宽 | 220 fs | 高斯脉冲 FWHM |
| 单脉冲能量 | 10~30 μJ | 可按实验条件调整 |
| 聚焦光斑半径 | 5 μm | 焦点处 1/e² 半径 |
| 脉冲个数 | 5 | 每发间隔 1 μs |
| 玻璃域 | 40×40×60 μm | 采用四分之一模型减少计算量 |
| 初始温度 | 300 K | 电子和晶格同温 |
如果激光正入射、光斑轴对称、且不考虑偏振方向，四分之一模型是合理的；一旦焦点有偏心或者入射有角度，老老实实用全模型。我在第一次测试时先用四分之一模型跑通物理逻辑，确认源项、边界、动网格都正常之后，再切到全模型做最终验证。这个流程能省下大量调试时间。
2. 双温模型在 COMSOL 里的落地写法：两个温度场如何进 PDE
2.1 双温方程组与系数型 PDE 的对应
双温模型的核心是电子温度 Te 和晶格温度 Tl 两个变量。电子系统吸收激光能量升温，然后通过电子-声子耦合机制把能量交给晶格，所以两个方程之间存在一个能量交换项 G(Te-Tl)。在 COMSOL 中，我用了两个系数型偏微分方程（Coefficient Form PDE）接口，分别求解 Te 和 Tl。
电子温度方程：
C_e(T_e) · ∂T_e/∂t = ∇·(k_e(T_e)·∇T_e) - G·(T_e - T_l) + S
晶格温度方程：
ρ_l · C_l · ∂T_l/∂t = ∇·(k_l·∇T_l) + G·(T_e - T_l)
COMSOL 系数型 PDE 的标准形式是 d_a·∂u/∂t + ∇·(-c·∇u) + a·u = f。把双温方程逐项对照，电子方程里 d_a 取 C_e，c 取 k_e，a 取 G，右侧 f 则写成 G·T_l + S。晶格方程里 d_a 取 ρ_l·C_l，c 取 k_l，a 取 G，右侧 f 写成 G·T_e。实际操作时，G、C_e、k_e 这些系数不要写死数值，而是做成全局参数表达式，后面做参数敏感性分析时直接改一处就行。
有一点要提醒：在移动网格下求解 PDE 时，COMSOL 会自动处理 ALE 框架下网格运动带来的输运效应，不需要手动在扩散方程里加对流项。如果计算中出现数值振荡，可以在系数型 PDE 里适当开启流线扩散稳定化，但不要一开始就开，否则会引入人为扩散，让烧蚀轮廓变得模糊。
2.2 激光体热源：高斯脉冲的时间与空间表达
飞秒脉冲的激光能量沉积，我用了高斯体热源近似。这个假设的潜台词是：激光在焦点附近被非线性电离产生的自由电子吸收，而不是在玻璃表面按线性吸收定律衰减。体热源表达式写成：
S(r,z,t) = (1-R) · E_p / (π·r0²·δ) · [2·sqrt(ln2/π)/τ_p] · exp(-r²/r0²) · exp(-z/δ) · exp(-4·ln2·(t-t_p)²/τ_p²)
其中 r 表示到焦轴的距离，对四分之一模型来说 r² = x² + y²；δ 是等效吸收深度，它并不是玻璃的固有光学吸收深度，而是一个把多光子吸收、雪崩电离、自由电子吸收等复杂机制压缩在一起的经验参数。对透明玻璃，δ 的量级通常落在 2 到 10 微米，具体数值可以先用单脉冲模拟跑一版，看烧蚀区尺度是否和实验量级一致，再反过来调整。
这个源项的时间部分是一个归一化高斯脉冲，对时间积分结果等于单脉冲能量密度，所以能量标定不会出错。R 是表面反射率，正入射的玻璃大概在 0.04 左右，可以保留也可以忽略，对结果影响很小。需要特别注意的是源项中的 z 坐标，究竟用空间坐标还是参考坐标，这直接关系到多脉冲第二发之后的焦点位置，我会在第 5 章专门展开。
2.3 玻璃材料参数的取舍：哪些必须标定，哪些可以先拍脑袋
玻璃的双温模型参数不像金属那么成熟，这也是这个模型最大的不确定来源。我建立了一套参考值，能覆盖大部分常见光学玻璃的量级：
| 参数 | 参考值 | 说明 |
|---|---|---|
| 密度 ρ | 2200 kg/m³ | 熔石英玻璃量级 |
| 晶格热容 C_l | 745 J/(kg·K) | 常温附近参考值 |
| 晶格热导率 k_l | 1.4 W/(m·K) | 玻璃导热系数较低 |
| 电子热容系数 γ | 700 J/(m³·K²) | 等效线性电子热容 |
| 电子-晶格耦合系数 G | 1×10¹⁵ ~ 1×10¹⁶ W/(m³·K) | 必须实验标定 |
| 烧蚀临界温度 T_crit | 2500 ~ 3200 K | 视材料牌号调整 |
很多金属的双温参数可以直接查文献，但玻璃本质上不是典型的双温材料，Te 在这里更应该理解为一个表征激光激发态能量水平的等效场变量。正因为如此，G 和 γ 的取值带有模型性质，不能照抄金属那套。我做完敏感性分析后发现，烧蚀深度对 G 和 δ 最敏感，这两个参数务必放到标定优先级的最前面。C_l、k_l 这些常规热物性参数反而对结果影响不大，用文献参考值即可。
2.4 初始条件与边界条件：杜绝多脉冲热积累失控
模型的初始状态很简单：Te = Tl = 300 K。边界条件的设置比初始条件更讲究。四周侧面我采用了绝热或者弱对流边界，底面固定为 300 K 的恒温边界。这样做不是为了省事，而是多脉冲问题的时间尺度已经到了微秒量级，热扩散深度可能接近几十微米。如果模型底面也设成绝热，热量会不断在模型内部累积，几个脉冲之后整体温度持续抬高，烧蚀深度越算越大，和实验结果对不上。底面给一个恒温热池，相当于模拟玻璃基底下方的远场散热，是更接近真实物理的处理方式。
如果使用四分之一模型，两个对称面上的温度梯度法向分量为零，直接设成绝热边界即可。等模型切到全模型后，四个侧面和底面全部按热池或者对流边界处理，只保留顶面这个激光入射面和烧蚀表面作为主关注区。
3. 变形几何：烧蚀边界怎么动、动多快、网格怎么扛
3.1 烧蚀判据与网格法向速度模型
变形几何要解决的核心问题是：材料被烧蚀之后，几何表面如何向内退缩。工程上不需要真的把材料单元一个个删除，而是让烧蚀表面以某个速度沿法向向内运动，相当于用一个移动边界近似材料去除过程。
烧蚀判据用的是晶格温度阈值。当表面某点的晶格温度超过临界值 T_crit，该处材料就被认为进入烧蚀状态。这里最忌讳的是把速度设置成一个阶跃函数，温度一超阈值就马上全速推进，这种不连续会给瞬态求解器带来很大的麻烦。我用的写法是一个光滑过渡的形式：
v_abl = v_max · smoothstep(T_l, T_crit, T_crit + ΔT)
也就是说，晶格温度从 T_crit 升高到 T_crit+ΔT 的范围内，边界速度从 0 平滑过渡到 v_max。ΔT 取几百开尔文即可，太小近似阶跃，太大则烧蚀停止延迟，影响形貌精度。v_max 的实际取值要满足一个硬约束：单个时间步内的网格位移不能超过局部最小单元尺寸的 30%，否则网格单元会被压成负体积，求解直接崩溃。
3.2 COMSOL 6.1 里移动网格的边界设置
在 COMSOL 6.1 中，我用的接口是组件下的移动网格（Moving Mesh），原理就是 ALE 动网格。烧蚀表面设置的是法向网格速度，其余底面和侧面全部固定不动。法向速度的大小就是上一小节的 v_abl，方向取表面外法线的反方向，也就是让网格向内退缩。
几何表面在三维模型里是平面，法向速度直接沿 z 方向；但对于带曲率或者边缘倒角的表面，法向向量在拐角处会突变，速度方向容易出现锯齿。我的对策是把烧蚀区域限定在光斑半径覆盖的局部表面范围，边界外围留一圈固定区域，避免速度方向在棱边上打架。COMSOL 6.1 的求解器序列里可以直接启用自动重剖分功能，当网格单元质量跌破阈值时自动中止当前步并重新剖分、插值后续解。这个能力在 5.x 时代要靠脚本实现，6.x 版本已经集成在求解器序列中，对烧蚀类模型几乎是刚需。
3.3 网格质量与单元尺寸的实操策略
网格策略直接决定模型能不能跑完五个脉冲。焦轴附近我使用了最大值 0.3 微米的高密度网格，远离焦点区域逐步放大到 2 微米。烧蚀深度可以到十几微米，但焦点区网格只有 0.3 微米，这意味着表面每退缩一步，第一层网格单元都会被明显压缩。如果 v_max 太大，一个步长内网格被压缩超过百分之三十，马上就会出现负单元体积。
我的解决办法是先用保守的 v_max 把单脉冲算通，再在激光脉冲结束后通过冷却段让速度自然回落到零。另一个容易被忽略的细节是，烧蚀表面的初始网格需要预留足够的层数。我在焦点区下方堆了至少 20 层 0.3 微米网格，相当于给形貌演化留出 6 微米以上的退缩空间，否则形貌还没发展起来，网格就已经到了不可用的状态。
4. 多脉冲推进：时间尺度相差 12 个数量级的循环策略
4.1 三套多脉冲方案对比：一整段积分、事件接口、逐脉冲重启
多脉冲建模最核心的问题不是物理，而是时间尺度的跨度太大。激光脉冲本身的宽度只有 220 飞秒，脉冲间隔却是 1 微秒，两者相差七个数量级；如果算到更长的热扩散阶段，跨度还能拉到十个数量级以上。我对比过三种方案。
| 方案 | 原理 | 优点 | 缺点 |
|---|---|---|---|
| 长时域周期源项 | 一段瞬态直接积分到底 | 设置简单，源项用周期脉冲函数 | 时间步长被迫在飞秒和微秒之间反复横跳，步数爆炸 |
| 事件接口加瞬态 | 用事件触发脉冲加载的时刻 | 时间轴上比较灵活，不用手动拆段 | 冷却段仍然需要大量廉价积分步，总计算量没有本质降低 |
| 参数化扫描加逐脉冲重启 | 每个脉冲单独跑一段瞬态，上一段结束解作为下一段初值 | 时间尺度被拆开，每个段的步长控制都很清楚，适合动网格 | 需要处理热累积边界，源项时间偏移要按脉冲序号重写 |
我最终选择的是第三种。它最大的好处是每个脉冲对应的瞬态研究都从一个相对干净的起点开始，时间域被限制在一个脉冲周期内，求解器的步长控制压力小得多。
4.2 逐脉冲循环的具体操作：参数化扫描与初值传递
具体操作是这样：定义一个参数 n，取值范围 1 到 5，代表脉冲序号。研究步骤使用瞬态求解器，时间范围设为 0 到 1 微秒，即一个完整的脉冲周期。激光源项的时间中心设在前 200 飞秒附近，后面留足冷却段。在参数化扫描设置中选择连续求解，这样第二发脉冲开始前，第一发结束时的 Te 和 Tl 分布会自动作为初始值传入。
电子温度的处理需要额外留意。物理上，脉冲结束后的冷却段里 Te 会快速向 Tl 回落，如果冷却段足够长，第二发脉冲开始时 Te 应该已经降到接近 Tl。我在冷却段设置了 500 纳秒以上的时间窗口，确保 Te 自然衰减到 300 K 量级再进入下一发。不要试图在每个循环开始时手动重置 Te，那样反而会在初始场里制造人为的温度跳变。
激光源项的时间偏移也用脉冲序号处理：在每个瞬态段落内，时间都是从 0 开始重新计数的，所以源项的时间高斯峰固定写在 t_p = 200 fs 附近即可，不需要写成绝对累计时间。这个思路把时间域复杂度从微秒量级压回到了皮秒量级，求解器压力骤降。
4.3 求解器设置与步长控制
三维模型的求解器设置比物理建模更容易被忽略，却往往是决定成败的一环。瞬态求解器我选用了 BDF，最大阶数设为 5。初始步长给到 1×10⁻¹⁵ 秒量级，让求解器从飞秒时间域一步一步爬起来。最大步长约束为脉冲宽度的五分之一左右，也就是大约 40 飞秒，避免时间步直接跳过脉冲峰值。
非线性求解用阻尼牛顿法，相对容差取 1×10⁻³，绝对容差则根据温度量级设置。电子温度在脉冲期间可能冲到上万甚至数十万开尔文，如果绝对容差设得太松，Te 的尖峰会被平滑掉，烧蚀区域判断就会偏差很大。我后来把绝对容差的参考值设为 300 K 乘以 1×10⁻⁴ 量级，即在 0.03 K 左右，同时给 Te 设置了一个下限截断，防止数值振荡导致负温度出现。
线性求解器在 PARDISO 和 MUMPS 之间，三维模型我最终用了 PARDISO，稳定性和内存占用都比 MUMPS 更合适，我这边需要的峰值内存约 64 GB。如果你的机器显存或者内存吃紧，优先缩小远离焦点区域的网格尺寸，而不是减少焦点区加密层数。
4.4 输出与后处理：别把形貌和温度场混在一起看
输出时间点也要分密疏。脉冲作用阶段每 20 飞秒输出一步都很值得，因为要捕捉 Te 峰值和烧蚀起止时刻；冷却阶段则每 10 纳秒保存一步就够了。后处理中观察烧蚀孔形貌时，要看变形后的空间坐标下的几何表面，而不是初始参考坐标下的表面。COMSOL 的移动网格结果里，空间坐标已经包含了网格位移，所以表面凹陷路径直接就是烧蚀轮廓。这一点关系很大，后处理时选错坐标，孔形可能完全看不出来。
5. 实测中必须伺候好的五个细节：排障链路记录
5.1 电子温度尖峰：数值尖峰还是物理尖峰
第一次完整跑通时，Te 在激光作用期间冲到了几百万开尔文，直觉告诉我这个量级有问题。排查过程是这样的：先检查源项积分有没有问题，确认能量标定无误；再检查时间步长，发现求解器在脉冲中心附近自动跨过了高斯峰，导致源项的实际积分不足，然后又通过容差控制把 Te 推得极高。解决办法有两个：一是在求解器的时间步进中强制指定脉冲中心前后若干时间点，让求解器在这些点之间加密步长；二是把源项时间包络适当展宽，算完后再做能量归一化校准。调整后 Te 的峰值保持在十万开尔文量级，这个范围和同类型问题的预期比较接近。
5.2 网格翻转：法向速度与时间步长的匹配
网格翻转是我遇到次数最多的崩溃方式，报错信息基本是负单元质量或雅可比矩阵奇异。第一次我把 v_max 设成了 1000 m/s，几乎每次都会翻转。后来把约束收紧到 v_max × Δt < 0.3 × h_min，这个问题立刻缓解。还有个细节，初期我把烧蚀面附近好几个边界的位移都放开，导致侧面网格被牵动，质量迅速劣化。正确做法是把所有非烧蚀边界全部固定，只留光斑覆盖的局部表面可动，并给速度做空间上的光滑处理，避免棱角处法向不连续造成的异常速度。
5.3 第二发脉冲开始后源项偏移：坐标选择是重灾区
这个坑我印象最深。跑到第二发脉冲时，烧蚀孔剖面变得很怪，孔壁几乎是垂直向下切，而不是实验中常见的倒锥形。逐层排查后发现，激光源项里用的 z 是空间坐标，而变形几何会让空间坐标随着网格移动。第一发烧蚀之后，表面已经往下退了几微米，同一个空间坐标 z 对应的材料点已经变了，第二发脉冲的焦点位置等于整体发生了偏移，相当于激光在追着已经烧掉的空腔跑。
解决办法是把源项中的焦点位置改用参考坐标来表达，也就是 COMSOL 里基于未变形几何的坐标。烧蚀形貌的演化仍然看空间坐标，但激光焦点的空间位置必须始终固定在初始材料坐标系下的那个点。改完之后，第二发之后的源项能量才能打在变形后的真实焦点区域上，孔形才恢复正常。
5.4 三维模型的热累积非对称：截断边界造成的伪烧蚀
四分之一模型最初跑出来有个明显异常：对称边界附近等温线正常，但截断边界附近出现了不该有的二次烧蚀区。原因是侧面绝热边界把本该向外扩散的热流反射回了模型内部，多脉冲热累积之后，边界附近温度被抬高，晶格温度越过了烧蚀阈值。处理办法是把侧面绝热边界改成对流热通量边界，并把底面固定为 300 K。如果模型域够大，比如横向尺寸推到 100 微米以上，边界效应影响会小很多，但计算量会同步上升。我最终在计算量和物理真实性之间取了 40 微米横向尺寸加对流边界的组合。
5.5 后处理里烧蚀轮廓的判定争议
最后说一个后处理中很容易犯的逻辑错误。很多人会习惯性把 Tl 超过 T_crit 的等值面当成烧蚀边界，但在带变形几何的模型里，这个操作是自相矛盾的。烧蚀表面已经被移动边界处理掉了，剩余材料域内的 Tl 只代表未去除部分的温度状态，直接用阈值等值面判断烧蚀轮廓，得到的结果往往偏大或者出现碎片区域。正确做法是以变形后几何的表面形状作为烧蚀边界，用剩余材料内部的晶格温度去判断热影响区和重铸层范围。如果你要对比不同脉冲数下的孔深，记得单独导出每个脉冲结束后的空间网格表面，而不是只保留最终时刻的结果。
6. 模型复盘：几个值得继续往深处走的扩展方向
6.1 从等效体热源升级为自由电子密度耦合
我当前模型用的是等效吸收深度 δ，虽然能跑出符合量级的烧蚀形貌，但物理上毕竟是一个压缩后的经验参数。更严格的思路是引入自由电子密度 n_e 作为第三个变量，把多光子电离率和雪崩电离率显式写进去，再让吸收系数依赖于 n_e。这个扩展会增加一个瞬态变量，但对能量沉积的阈值特性、脉冲前沿的陡峭吸收过程能描述得更细腻。代价是求解规模和稳定性风险同步上升，建议先跑二维验证再往三维扩展。
6.2 把高斯本征源替换为光束传播场
当前模型的源项假设焦点处是理想高斯分布，没有考虑透镜像差、色散和非线性自聚焦效应。如果想研究焦点偏移、不对称光斑或者脉冲空间分布对烧蚀形貌的影响，就需要把源项替换为近轴光束传播方程，甚至全矢量电磁场仿真。这样做的计算代价很大，而且时间尺度没有变容易，适合作为第二阶段的专题研究。
6.3 实验标定与参数反演：先单脉冲再多脉冲
三维多脉冲模型的可信度完全取决于那几个标定参数。我的建议是标定顺序不要乱：先用单脉冲实验测出烧蚀直径和深度，反推 δ 和 G 的乘积范围；再用不同重复频率下的多脉冲孔深增量，标定 v_max 和 T_crit。如果一上来就试图用多脉冲结果同时拟合所有参数，得到的只会是一组局部最优值，换一组工艺条件就失效。
如果让我总结一句最核心的经验，我会说：不要试图让模型一开始就吃掉全部物理。先把单脉冲加双温模型跑稳，再打开变形几何，最后叠多脉冲循环，每加一层之间留出足够时间做参数敏感性分析。这样即使出现网格翻转、源项偏移这类问题，也能快速定位是哪个环节引入的。这套流程走完之后，三维多脉冲烧蚀模型才算真正从一道数学难题变成一个可复现的工程工具。
