飞秒激光打到材料表面后纳米尺度内发生了什么如果你的第一反应是打开MATLAB直接用普通热传导方程算算完大概率会发现电子温度和晶格温度根本对不上实验。原因在于飞秒尺度下电子和晶格处于完全不同的温度状态必须用双温模型把电子子系统和晶格子系统分开描述如果材料是半导体还要再把载流子密度作为独立变量拉进来。这篇文章我来完整拆解带载流子密度的双温模型在MATLAB里如何用有限元法求解从方程推导、参数标定、代码框架到调试经验一步不落。适合正在做超快光学、飞秒激光加工或者计算材料模拟的读者参考。1. 飞秒激光模拟为什么绕不开双温模型从物理图像说起1.1 双温模型的基本假设与适用场景飞秒激光脉冲的能量首先被材料中的自由电子吸收而不是直接均匀地变成整个材料的热量。光子激发电子跃迁到高能态电子子系统内部通过电子-电子散射在极短时间几十飞秒到几百飞秒内建立准平衡此时电子温度T_e可以升到几千甚至上万开尔文。但晶格原子核因为质量大通过电子-声子耦合接受能量是一个相对缓慢的过程典型时间常数在几百飞秒到几皮秒。于是在脉冲后的一小段时间内出现了电子热、晶格冷的非平衡状态。双温模型就是描述这种非平衡态的经典框架。它的核心是把材料划分为电子子系统和晶格子系统用两个耦合的热传导方程分别求解T_e和T_l。两个方程之间通过一个能量交换项G(T_e−T_l)绑定G就是电子-声子耦合系数。这个模型从Anisimov在上世纪70年代提出至今依然是飞秒激光与金属、半导体相互作用模拟的默认起点因为它的形式足够简洁物理图像足够清晰又能在很宽的激光参数范围内复现实验规律。什么时候必须用双温模型而不是普通热传导判据很简单激光脉宽和电子-晶格能量弛豫时间可比或者激发的能量密度足够让电子温度与晶格温度出现显著差异。脉宽在几十飞秒到几皮秒的场景基本都满足这一条。反过来如果是纳秒激光或者连续激光电子和晶格来得及在脉冲持续时间内保持热平衡双温模型退化为经典热传导模型普通求解器就够用了。1.2 什么时候必须引入载流子密度经典双温模型里载流子密度不是独立变量。对金属来说自由电子密度近似常数电子热容C_e常用γT_e近似热导率k_e用随温度变化的经验表达式这样已经能描述大部分飞秒激光-金属相互作用。但在半导体、绝缘体、相变材料甚至石墨烯等体系中飞秒激光产生的非平衡载流子密度变化剧烈这个密度n不再是一个常数它会显著影响光学吸收、电导率、电子热容以及复合放热。一个典型的例子是飞秒激光加工硅。激光脉冲在几百飞秒内把表层电子密度从本底的10^16 cm^−3量级激发到10^20~10^21 cm^−3量级材料的光学性质、热输运性质随之突变电子-空穴对通过俄歇复合把大量能量释放给晶格又进一步改变烧蚀行为。此时如果不把n作为独立变量求解电子热容、吸收深度、复合热源全都定不下来模拟结果会明显失真。因此完整的“带载流子密度的双温模型”实际上是一个三方程系统载流子密度方程、电子能量/温度方程、晶格温度方程。金属情况下n固定为平衡密度方程退化为标准双温半导体情况下三方程耦合求解才是严谨做法。下面所有讨论都基于这个框架展开。1.3 常见误区把双温模型当成普通热传导我见过不少刚入门的同学拿着双温模型的方程却用普通热传导的思路去处理最后得到一堆看起来合理、实际没有物理意义的结果。最常见的三个误区把激光源当恒定体热源飞秒激光的功率密度在脉冲内剧烈变化脉宽150 fs的脉冲在1 ps时间尺度上已经结束。如果按恒定功率算峰值温度和时间演化会完全错误。忽略电子热导率的温度依赖性电子热导率k_e(T_e, T_l)在飞秒尺度不是常数。电子温度越高热扩散越快峰值温度往往低于常数k_e的预测。把表面温度直接当损伤判据实验上观测的损伤阈值更多与晶格温度达到某个临界值相关而表面电子温度峰值出现得早、数值高但电子温度高并不直接等同于损伤。需要用晶格温度、载流子密度或者结合相变判据来判断烧蚀深度和阈值。理解了这些误区后面方程离散和参数设定就有了方向。2. 带载流子密度的双温方程组从系数标定到方程离散2.1 三个方程如何耦合电子、晶格、载流子以一维厚度方向z为例表面在z0材料向z0方向延伸。三个方程的强耦合形式如下载流子密度方程 ∂n/∂t ∂/∂z[D_n(n, T_e) ∂n/∂z] S_n(z,t) − n/τ_r − Γ_n n^3右边第一项是载流子扩散第二项是激光光生载流子源项第三项是线性复合缺陷辅助复合等第四项是俄歇复合。电子能量/温度方程 C_e(T_e, n) ∂T_e/∂t ∂/∂z[k_e(T_e, T_l) ∂T_e/∂z] − G(T_e, T_l)(T_e − T_l) S_e(z,t) 化学势项化学势项在飞秒时间尺度内通常贡献小很多文献直接省略。保留它会让方程变成带对流特征的方程数值上更难处理。晶格温度方程 C_l ∂T_l/∂t ∂/∂z[k_l(T_l) ∂T_l/∂z] G(T_e, T_l)(T_e − T_l) 载流子复合放热项这里最关键的耦合体现在四处G项双向交换能量k_e和C_e依赖T_e、T_l甚至n载流子复合放热直接给晶格加热光生载流子源项S_n和电子热源项S_e都与激光强度和吸收深度绑定。金属情况下n保持常数S_n和复合项消失方程就回到大家熟悉的双温模型。2.2 材料参数的物理含义与数值选择参数是整个模拟的灵魂。模型写错了可以改参数给错了结果基本没法用。下表给出金和单晶硅在飞秒激光模拟中常用的一组参数注意这些是室温附近的常用值实际使用时必须以目标材料的实验表征数据为准。参数符号金Au常用值单晶硅Si常用值单位电子比热系数γ21—J/(m^3·K^2)电子/载流子热容C_eγ·T_e约等于 1.5n k_BJ/(m^3·K)晶格热容C_l2.5e60.7e6~1e6J/(m^3·K)电子-晶格耦合系数G2.1e160.1e16~1e16W/(m^3·K)电子热导率k_e0315硅中较小载流子导流W/(m·K)晶格热导率k_l315130W/(m·K)光学穿透深度δ15500~700近红外nm反射率R0.93近红外0.3~0.5无量纲电子热容CeγTe是金属的自由电子气近似在Te远低于费米温度时成立。硅这类半导体光生载流子浓度极高时热容表达式更接近简并电子气形式有时用Ce≈(3/2)nk_B估计。G的取值对结果影响最大而且文献值分散度很高后面我会专门说这个坑。k_e的温度依赖在经典双温模型里常用k_e k_e0·(T_e/T_l)作为一阶近似。物理含义是电子热导率正比于电子热容乘电子扩散系数而电子-声子散射限制了平均自由程所以T_e升高会增强热扩散。但在极端高电子温度下电子-电子散射开始贡献电阻k_e实际增长速度会放缓有人用k_e k_e0·(T_e/T_l)·(A/(AB·T_e^2))这类修正表达式实际模拟时建议先跑通一阶形式再逐步加复杂修正。2.3 激光热源项的数学表达式激光源项在飞秒双温模拟里是最容易被写错的地方因为涉及到脉冲时间分布、空间吸收分布和能量归一化三件事。假设激光脉冲时间强度分布为高斯型 I(t) F·sqrt(4 ln2 / (π t_p^2))·exp(−4 ln2·(t−t_0)^2 / t_p^2)其中F是激光能量密度J/m^2t_p是脉冲半高全宽FWHMt_0是脉冲中心时刻。这个表达式严格保证∫I(t)dt F不会因为脉宽不同而出现能量缩放错误。光在材料内部按指数衰减吸收单位体积的吸收功率密度为 S_total(z,t) (1−R)·I(t)·exp(−z/δ)/δ这个表达式里除以δ是关键因为∫0^∞ exp(−z/δ)/δ dz 1这样才能保证单位面积上吸收的总能量等于(1−R)·F。很多初学者漏掉这个δ导致结果能量偏差几十倍。对于半导体和载流子激发问题S_total需要拆成两部分一部分用于产生电子空穴对进入载流子源项S_n剩余能量转化为电子热能进入电子温度方程S_e。常用做法是按光子能量hν和带隙E_g拆分把E_g/(hν)份额分给载流子产生剩下的分给电子热。2.4 边界条件的合理简化边界条件是另一个容易翻车的地方。表面对空气/真空的热对流在飞秒时间尺度内可以完全忽略所以边界上通常取绝热条件 −k_e ∂T_e/∂z|(z0)0−k_l ∂T_l/∂z|(z0)0∂n/∂z|_(z0)0远侧边界上如果模拟时间足够短几十纳秒内热量还没传到远处直接取固定温度和平衡载流子密度 T_e(L)T_l(L)T_0n(L)n_eq一维模型厚度L要足够大一般取几微米到几十微米确保模拟时间窗口内远侧温度不发生变化。如果模拟时间拉到微秒以上远边界条件建议换成对流换热或者继续加厚求解域不要让固定温度边界影响结果。3. MATLAB有限元求解三条路径工具箱、符号推导还是手写组装3.1 PDE Toolbox系数型方程映射技巧MATLAB自带的PDE Toolbox支持2D/3D的有限元求解1D问题也能用createPDEModel(1)建立然后通过pdeCoefficients定义系数型方程。系数型PDE的标准形式是 d·∂u/∂t − ∇·(c∇u) a·u f双温模型的三个方程可以写成这个形式dC对应热容ck对应热导率a0f包含G项和激光源项。PDE Toolbox支持的u是分量列向量所以三个方程可以分别映射到u1Te, u2Tl, u3n。实际用下来PDE Toolbox在处理简单线性问题上很舒服但双温模型是非线性强耦合的系数c、d、f都依赖u。每次时间步需要重新调用assembleFEMatrices更新矩阵这会引入不少额外开销和配置复杂度。加上PDE Toolbox在1D问题上的网格控制不如手写灵活所以我个人只在做2D轴对称快速验证时用它。1D厚度方向模拟手写有限元完全够用而且代码透明、可控、容易排查问题。3.2 一维线性单元有限元刚度矩阵推导一维线性单元是最基础的有限元单元。设节点坐标z_i单元e连接节点i和i1单元长度为h_ez_{i1}−z_i。单元自由度上线性形函数为N_1(z_{i1}−z)/h_eN_2(z−z_i)/h_e。以电子温度方程为例弱形式是∫v·C_e ∂T_e/∂t dz ∫k_e ∂v/∂z·∂T_e/∂z dz ∫v·f_e dz其中f_e包含−G(T_e−T_l)S_e。在单元上质量矩阵元素M_ij∫N_i·N_j·C_e dz刚度矩阵元素K_ij∫∂N_i/∂z·∂N_j/∂z·k_e dz载荷向量f_i∫N_i·f_e dz。对线性单元这些积分有解析值单元质量矩阵是(h_e·C_e/6)·[2 1; 1 2]单元刚度矩阵是(k_e/h_e)·[1 −1; −1 1]。单元载荷向量如果f_e在单元内变化不大可以用(h_e/2)f_e近似如果激光源项在表面附近变化剧烈建议把该区域的网格加密到源项空间变化能被很好分辨的程度必要时用两点高斯积分提高精度。这里要注意C_e和k_e是温度的函数每个单元上取单元中心温度计算才能保持有限元离散的合理性。这个“系数取值位置”是手写代码和教科书简化之间最容易出偏差的地方。3.3 时间离散与刚性问题隐式格式和ode15s的选择三个方程的特征时间尺度差异很大。激光脉冲期间电子温度变化极快载流子密度复合寿命可能在纳秒量级晶格温度响应则在皮秒到纳秒。整个系统是典型的刚性方程用显式时间推进会严重受限于稳定性条件步长小到无法接受。我推荐两种策略一是手写隐式欧拉或者Crank-Nicolson。隐式欧拉简单稳定时间步长可以比显式格式大几个量级Crank-Nicolson精度高但有振荡的风险。对强非线性问题隐式欧拉配合足够小的时间步是优先选择。二是直接用MATLAB的ode15s或ode23t工具箱函数把空间离散后的半离散方程组也就是一个大尺度常微分方程组直接丢给求解器。优点是变步长自动处理缺点是耦合方程的Jacobian计算有时比较脆弱。不过我实测下来对于节点数量在几百到一两千的一维问题ode15s很可靠适合作为交叉验证手段。4. 从方程到代码核心MATLAB实现框架4.1 主程序结构和数据流我习惯把整个模拟拆成四个层级参数定义、网格生成、矩阵组装、时间推进。这样改参数、改材料、加物理项都只动局部不用把整个脚本推倒重来。参数定义集中在脚本开头包括激光参数、材料参数、网格参数和时间参数。网格生成用线性分布加加密表面前几百纳米用均匀细网格纵深方向网格逐步拉大。矩阵组装把质量矩阵、刚度矩阵、载荷向量都做成稀疏矩阵避免全矩阵存储。时间推进用循环控制步长每步更新温度相关系数后再组装求解。数据流就是参数→网格→初始条件→时间循环→后处理。下面给出核心代码片段以金属金为例n固定为常数退化为双温模型扩展示例里我会说明怎么把半导体载流子方程加回来。4.2 单元组装和稀疏矩阵构建% 1D 有限元组装示例金双温退化版n固定 L 2e-6; % 材料厚度 2 um N 800; % 节点数 z linspace(0, L, N); dz z(2) - z(1); % 均匀网格后续可替换为非均匀 % 单元参数 gamma 21; % 电子比热系数 J/(m^3 K^2) Cl 2.5e6; % 晶格热容 G 2.1e16; % 电子-晶格耦合系数 k0 315; % 电子平衡热导率 kl 315; % 晶格热导率 F 50; % 激光能量密度 J/m^2 tp 150e-15; % 脉宽 FWHM t0 3 * tp; % 脉冲中心 R 0.93; % 反射率 delta 15e-9; % 光学穿透深度 % 初始条件 Te T0 * ones(N,1); Tl T0 * ones(N,1); % 稀疏矩阵预分配三对角结构 M sparse(N,N); K sparse(N,N); Ml sparse(N,N); Kl sparse(N,N); % 线性单元组装 for e 1:N-1 n1 e; n2 e1; he z(n2) - z(n1); % 电子质量矩阵/刚度矩阵系数基于当前Te, Tl计算 Te_mid 0.5*(Te(n1)Te(n2)); Tl_mid 0.5*(Tl(n1)Tl(n2)); Ce_mid gamma * Te_mid; ke_mid k0 * Te_mid / max(Tl_mid, 1); % 防止除零 Me_loc he * Ce_mid / 6 * [2 1; 1 2]; Ke_loc ke_mid / he * [1 -1; -1 1]; M(n1,n1) M(n1,n1) Me_loc(1,1); M(n1,n2) M(n1,n2) Me_loc(1,2); M(n2,n1) M(n2,n1) Me_loc(2,1); M(n2,n2) M(n2,n2) Me_loc(2,2); K(n1,n1) K(n1,n1) Ke_loc(1,1); K(n1,n2) K(n1,n2) Ke_loc(1,2); K(n2,n1) K(n2,n1) Ke_loc(2,1); K(n2,n2) K(n2,n2) Ke_loc(2,2); end这段代码放在时间循环外作为“冷组装”框架时间循环里只需要用更新后的系数重复组装局部矩阵。注意稀疏矩阵的更新不是零成本但对800个节点的一维问题每次组装都很快不要为了优化而过早把代码搞复杂。4.3 时间推进循环的实现细节时间循环的核心是求解耦合方程组。以隐式欧拉为例电子方程离散为 (M_e/Δt K_e) T_e_new (M_e/Δt) T_e_old − G .* (T_e_old − T_l_old) S_e晶格方程离散为 (M_l/Δt K_l) T_l_new (M_l/Δt) T_l_old G .* (T_e_new − T_l_new)这里的处理有一个常见简化G项在晶格方程里用了半隐式也就是把T_e_new代入晶格方程同时T_l_new出现在右边这样可以提高稳定性。但如果在同一大步里T_e_new更新过于激进T_l_new可能出现非物理振荡。实际更稳妥的做法是再内迭代一两次。dt 0.02e-15; % 初始时间步试算后调整 Nt 5000; for it 1:Nt t it * dt; % 激光源项高斯脉冲时间分布 It F * sqrt(4*log(2)/(pi*tp^2)) * exp(-4*log(2)*(t-t0)^2/tp^2); Se (1-R) * It * exp(-z/delta) / delta; Ae M / dt K; be M * Te / dt Se - G * (Te - Tl); Te_new Ae \ be; Al M / dt Kl; bl M * Tl / dt G * (Te_new - Tl); Tl_new Al \ bl; Te Te_new; Tl Tl_new; % 后半段自动放长时间步 if t 5e-12 dt 0.1e-12; elseif t 1e-9 dt 1e-12; end end上面这段代码里Se加到了电子方程右边但在激光脉冲早期电子温度可能瞬间跳变。如果时间步太大Te_new会有明显振荡。实战中我的习惯是在脉冲期间强制用小于1 fs的步长脉冲结束后逐步放宽。条件切换步长的方法放在循环里虽然简单但要注意切换时MATLAB的三对角求解不会因为步长突变出问题代数方程本身还是良态的。4.4 后处理电子温度、晶格温度和载流子密度的时空演化算完之后最常用的画图有两类一类是某个时刻的温度空间分布曲线另一类是温度随深度和时间的二维伪彩色图。伪彩色图用contourf或者pcolor都可以但时间轴要取对数因为飞秒到纳秒跨越了7个数量级。figure; contourf(z*1e9, t*1e12, Te_matrix, 50, LineStyle, none); set(gca, YScale, log); xlabel(深度 z (nm)); ylabel(时间 t (ps)); colorbar; title(电子温度演化);晶格温度、载流子密度同样处理。如果想提取表面峰值电子温度随能量密度的变化可以用一个循环扫描不同F值然后把峰值T_e或峰值T_l画成F的函数这样就能和损伤阈值实验对照。半导体情况加载流子方程的做法是在矩阵组装和时间推进里增加第三个变量。组装方式和电子方程类似C_n1D_n作为扩散系数源项S_n(1−R)I(t)exp(−z/δ)/(hν·δ)再把复合项−n/τ_r和俄歇项−Γn^3放进f_n。晶格方程里增加复合放热项E_g·n/τ_r。代码结构完全一致只是矩阵从2×2分块变成3×3分块。5. 实测结果与调试经验网格、时间步、参数敏感性的坑5.1 网格布置在吸收深度内加密单元光学穿透深度δ决定了激光能量集中沉积的深度范围。以金为例δ只有15 nm这个区域内源项强度是深层区域的几十倍网格如果太粗有限元离散会把尖峰源项平均掉导致表层温升偏低、能量沉积分布被抹平。我的经验是表层0~3δ范围内单元尺寸取δ/10甚至更小3δ~10δ范围内逐步放大到δ量级更深处按几何级数加密到整体L。800节点的网格里前200 nm分配给300~400个节点其余节点均匀铺开。这样做既能保证源项分辨率又不至于让节点数爆炸。网格加密后要注意检查结果是否收敛。把节点数从500提到1000如果峰值温度变化超过5%就继续加密如果变化在1%以内基本可以认为网格足够。5.2 时间步长选择从飞秒到皮秒/纳秒的跨越时间步长策略是整个模拟最容易让新手崩溃的地方。用固定步长跑全程要么慢到无法接受要么因为步长太大出现温度振荡甚至负值。我的推荐策略分三个阶段时间范围物理过程建议步长0~2 ps脉冲吸收与电子热化0.01~0.1 fs2 ps~1 ns电子-晶格能量交换与热扩散0.1~5 fs 逐步放宽1 ns以上晶格主导热传导10 fs~1 ps为什么脉冲期间要那么小的步长因为高斯脉冲在峰值附近的时间导数极大150 fs脉宽中心附近每飞秒能量变化相对幅度达到百分之一量级。如果时间步超过脉宽的百分之一源项时间采样误差就会明显影响峰值电子温度。时间步放大别直接跳。我常用对数渐变让步长每隔几步均匀放大20%~50%而不是从0.01 fs直接跳到1 ps。你会发现系统对步长的突然变化很敏感本质上是因为温度演化速率在不同阶段切换但半隐式格式对系数突变会引入额外误差。5.3 结果合理性检查能量守恒和温升特征任何仿真跑完第一件事不是看曲线好不好看而是做能量守恒检查。把整个求解域上积分的总吸收能量算出来 E_abs (1−R)·F·A_surface然后在每个时间步积分T_e和T_l的能量增量 E_sim ∫C_e(T_e)dT_e体积积分 ∫C_l T_l体积积分 − 初始内能两个值偏差超过5%基本可以断定某个环节出了问题。我遇到过的违规原因有反射率R用错波长、δ漏除导致源项能量翻倍、边界条件在长时间模拟中导致能量泄漏。另一个常见的物理特征检查是温度弛豫时间。以金为例电子温度和晶格温度在脉冲结束后经历一个快速拉平过程典型时间在几皮秒量级。如果模拟中观察到T_e和T_l在皮秒内就完全重合可能是G设得偏大如果过了几十皮秒还差几千K则G可能偏小。这个“弛豫时间标尺”比单纯看温度峰值更能暴露参数错误。5.4 参数敏感性G值对结果的影响在所有材料参数里G值对结果的影响几乎是一票否决级的。文献中金的G值从1.1e16到4e16 W/(m^3·K)都有报道不同实验手段和薄膜质量会让这个值出现明显差异。我自己做过一个粗糙的敏感性测试固定其他参数把G从1.5e16改到3.5e16结果表面峰值晶格温度相差近20%电子-晶格平衡时间从1.5 ps变到0.7 ps。这种差距远超数值误差。所以不要迷信某个固定G值。正确做法是先用中值参数跑一次然后做G的敏感性扫描把T_g峰值随G变化画出来。再对照你实验室的损伤阈值或者反射率变化数据反推当前材料状态下的有效G值。这个过程比盲调其他参数有效得多。k_e的温度依赖表达式也会带来显著差异。用k_ek0·Te/Tl和常数k0这两种选择峰值电子温度差可以达到30%以上。建议至少跑两种表达式做上下界估计心里有数之后再确定最终模型形式。5.5 半导体载流子密度方程的稳定性问题最后单独提醒一下半导体情况。载流子密度方程里的复合项−n/τ_r和俄歇项−Γn^3是强非线性项当n到达10^27 m^−3量级时俄歇项可能引发数值刚性问题。我遇到过最典型的现象是载流子密度正常但晶格温度在复合放热项加入后出现负值振荡。根因是俄歇复合放热速率极高电子方程和晶格方程的时间尺度严重失配用大步长时把局部热源瞬间注入晶格温度就过冲了。解决思路有两个一是用隐式格式处理复合项让放热速率与当前温度同步更新二是把复合放热项拆分一部分在电子方程中扣除一部分进入晶格方程避免能量凭空产生。无论哪种方案都要确保总能量守恒校验通过后再跑正式参数。我对这套模型的实际体会是物理上最脆弱的地方永远是参数而不是方程和算法。双温模型及其载流子扩展形式本身已经足够成熟但每个材料的G、k_e、C_e表达式都带着强烈的样品依赖性。无论你是在复现文献还是做新材料的预测建议把参数敏感性分析作为标准流程而不是最后补救的手段。这样跑出来的结果你在审稿人面前也能站得住脚。