带载流子密度的双温模型：MATLAB有限元实现飞秒激光仿真

飞秒激光打金属或半导体薄膜很多刚上手的人会被一张温度分布图问住脉冲只有几百飞秒怎么材料表面温度就上万K了其实这里的温度并不是我们平时说的“材料温度”而是电子子系统的温度。几百飞秒的激光脉冲能量集中打在电子气上电子-电子散射时间极短电子气迅速形成高温分布而晶格因为原子质量大、电声耦合响应慢温度几乎还留在原地。这正是双温模型TTM要描述的核心场景。项目标题里的“带载流子密度的双温模型”是在基础双温方程之上引入光生载流子密度N的瞬态方程让模型能同时覆盖金属自由电子和半导体激发载流子的情形。这篇文章我围绕一个实际问题展开如何在MATLAB里用有限元法把电子温度场、晶格温度场和载流子密度场三个耦合的偏微分方程解出来并把飞秒激光源项正确加载进去。1. 双温模型到底在模拟什么物理图景与方程推导1.1 飞秒加热的“双温”本质先讲清楚物理时间尺度不然方程会看得一头雾水。飞秒脉冲激光的典型脉宽在几十到几百飞秒这个时间尺度比电子-电子散射时间大约10~100飞秒略长但远短于电子-声子耦合时间大约几百飞秒到几皮秒。激光能量首先被电子通过逆韧致辐射、带间吸收等机制吸收于是电子系统被加热到几千甚至上万K而晶格系统还没来得及响应。随后电子通过电子-声子耦合把能量慢慢交给晶格最终在几皮秒到几十皮秒内两者达到热平衡。所以这里的“双温”不是数值上的花架子而是超快激光与材料相互作用过程中必然出现的热力学中间状态。传统热传导方程假设电子和晶格始终处于同一温度在纳秒脉冲甚至连续激光下基本成立但在飞秒脉冲下完全不适用。做模拟时必须让电子子系统有自己的温度场T_e(x,t)晶格子系统有自己的温度场T_l(x,t)两者通过电声耦合项交换能量。这就是标题里“电子晶格温度”的来历——不是某个固定温度而是两条温度曲线的时空演化。1.2 载流子密度方程为什么必须加进来基础双温模型只写T_e和T_l两个方程适用于金属因为金属内部自由电子密度基本固定能量吸收和输运行为可以完全划入电子热容和热导率。但一旦把目光转到大带隙半导体、或者需要研究薄膜中光生载流子的输运时就不能默认电子密度不变了。飞秒脉冲激发出高密度电子-空穴对这部分载流子会通过扩散迁移、通过辐射和非辐射复合再释放能量改变后续能量沉积的空间分布。载流子密度场N(x,t)本身也要满足一个带源项和复合项的扩散方程并且通过复合项把能量回馈给电子或晶格系统。很多做激光损伤研究的同行容易忽略这一点直接用纯双温模型算半导体结果电子温度峰值明显偏高尤其是表面温度。原因就在于没有考虑非平衡载流子的扩散和Auger复合等耗能过程。把载流子密度方程加进来之后模型会真实很多代价是方程数量从两个变成三个耦合关系更加紧密数值求解难度也上一个台阶。这也正好对应了这个项目标题的前半句“带载流子密度”。1.3 不同材料的方程形式差异双温模型的具体方程形式需要根据材料类别调整。对于金属一个比较经典的写法是电子系统C_e(T_e) ∂T_e/∂t ∇·(k_e(T_e,T_l)∇T_e) - G(T_e - T_l) S(x,t)晶格系统C_l ∂T_l/∂t G(T_e - T_l)这里C_e是电子热容k_e是电子热导率G是电声耦合系数S是激光体热源。晶格热传导项通常被省略因为在皮秒量级晶格热扩散对温度场影响很小但如果你关注几纳秒以后的热分布可以加回k_l∇T_l项。对于半导体需要增加载流子密度方程常见形式是∂N/∂t ∇·(D∇N) - N/τ_r - γ_A N³ S_N(x,t)其中D是载流子扩散系数τ_r是载流子复合寿命γ_A是Auger复合系数S_N是光生载流子源项。同时电子能量方程里要加一项载流子复合释放的热量比如用能带带隙乘以复合率或者用一个等效热源项耦合进去。实际项目里具体采用哪种形式取决于你手里的材料参数和实验要对比的物理量。如果目标材料是金、铜、铝这类金属载流子方程可以直接关闭节省大量计算量如果目标是硅、砷化镓或者氧化铟锡这类宽禁带/透明导电材料就必须把N方程打开。2. 参数准备与边界条件仿真的地基2.1 核心热物性参数的表达式与取值双温模型里最让人头疼的就是一堆温度相关的物性参数而且不同文献给的数值能差出好几倍。先说最常见的处理方式电子热容用自由电子气模型写成C_e γ T_e其中γ是电子热容系数单位是J/m³K²。电子热导率不能直接用室温常数因为飞秒加热后电子温度和晶格温度差别巨大。工程上常用低激发近似k_e k_0·T_e/T_l其中k_0是室温附近电子热导率。这个表达式在T_e远大于T_l时会让电子热导率变得很大仿真时容易出现表面温度被快速抹平的情况属于模型本身的局限。也有一些文献用平滑过渡公式k_e k_0·T_e^5/(T_e^5T_l^5)来避免低温奇异性我自己在调试中更推荐后者数值上更稳定。下表给出几种常见材料的参考参数注意这些数值在不同文献里差异较大用于初步调试可以真要对比实验需要自己拟合或标定。材料γ (J/m³K²)k_0 (W/mK)G (W/m³K)C_l (J/m³K)备注金 Au约68约315约2.2×10^16约2.5×10^6电子热导率大电声耦合弱铜 Cu约97约386约1×10^17约3.4×10^6电声耦合比金强铝 Al约135约235约2~3×10^17约2.4×10^6轻金属升温更快硅 Si查表或用Drude模型约150约1×10^16量级约1.66×10^6半导体需开N方程吸收系数波长相关我个人的经验是G的取值对最终温度峰值影响最敏感尤其是金属薄膜的损伤阈值模拟。文献里给出的G都是静态值但其实它随电子温度和材料微结构变化很大。如果实验能测到电子-晶格弛豫时间τ用G≈C_e/τ反推最靠谱。2.2 飞秒激光源项时空分布与吸收激光源项是驱动整个模型的关键也是最容易出错的环节。先写时间分布飞秒脉冲通常用高斯近似I(t) I_0 · exp(-4ln2·(t-t_0)²/τ_FWHM²)式中τ_FWHM是脉冲半高全宽t_0是脉冲中心时刻。如果给定脉冲能量密度F单位J/m²峰值强度I_0可以由能量守恒算出来F ∫_{-∞}^{∞} I_0·exp(-4ln2·(t-t_0)²/τ_FWHM²) dt I_0·τ_FWHM·√(π/(4ln2))所以I_0 ≈ 0.94·F/τ_FWHM。这个换算很多人第一次会漏直接用F除以τ_FWHM算I_0差别不大但不够严谨。空间上激光在材料内部的吸收按Beer-Lambert定律放在方程源项里而不是放在边界条件里S(x,t) α·(1-R)·I(t)·exp(-αx)其中α是材料对激光波长的吸收系数或有效吸收深度倒数R是表面反射率x是深度。这里的单位是W/m³注意和表面热流密度W/m²区分。对金属而言光吸收主要是趋肤效应有效吸收深度通常是十几纳米对半导体则直接用带间吸收系数。如果在二维里模拟聚焦光斑还要在横向乘一个高斯空间因子exp(-2r²/w_0²)。2.3 边界条件与初始温度设置边界条件的设置比想象中重要。薄膜表面X0处激光能量已经被放进体源项里所以表面边界通常设为绝热即∂T/∂x0避免重复计算能量。薄膜背面和衬底接触那侧可以按实际情况设为恒温300K或者绝热。如果你模拟的是自由悬挂薄膜且关注时间在几十皮秒内背面绝热一般够用。载流子密度方程在表面处通常取零通量边界背面同样零通量因为载流子跨界面输运涉及表面复合模型基础版本里不处理。初始条件相对简单电子温度、晶格温度都设为环境温度300K初始载流子密度设为平衡值或者0。但要注意在时间起点之前需要留一段零激光阶段让求解器起步否则瞬间加载强源项会导致初始时间步过小。我习惯把脉冲中心放在t_05倍脉宽左右时间窗口从0开始前几百飞秒没有激光数值温度曲线会更平滑。3. MATLAB有限元实现从方程到可运行代码3.1 网格划分与PDE工具箱配置很多人看到“有限元法”就以为一定要上PDE工具箱。工具箱确实能用但说实话面对这种强非线性、三方程耦合、源项时间尺度跨越好几个数量级的问题PDE工具箱的parabolic函数并不好用它处理线性抛物型PDE比较顺手但T_e相关热容和热导率是随T_e变化的系数矩阵每个时间步都要重新组装工具箱内部实现起来很别扭。我的建议是先用MATLAB内置的几何和网格函数做一个结构网格然后自己组装有限元刚度矩阵和质量矩阵。这样计算流程透明调试起来也容易定位问题。一维问题是最佳起点沿着薄膜厚度方向建网格节点数通常取100到400。网格步长要小于电子热扩散长度我一般控制在5纳米以内对于200纳米厚的薄膜大约划分200个节点。二维轴对称情况下可以用PDE工具箱的geometryFromEdges和generateMesh生成网格再把节点坐标和单元连接信息导出后续矩阵组装逻辑和一维完全一致。3.2 非线性系数的迭代处理模型里的非线性体现在三处电子热容C_e(Te)随Te变化、电子热导率k_e(Te,Tl)随温度变化、载流子Auger复合项γ_A·N³。如果直接交给MATLAB的ode15s去处理需要在右端函数里实时组装修正后的稀疏矩阵。基本流程是每个时间步开始时用当前解更新温度依赖系数组装质量矩阵M和刚度矩阵K然后解线性方程组得到时间导数。这里有一个关键经验不要在每个时间步里面手动做Newton迭代收敛。对大多数参数范围Picard迭代就够了也就是用上一步温度场计算系数然后推进一步。如果Te变化太剧烈导致振荡再考虑在同一时间层内循环两到三次Picard。我自己实测下来ode15s配合Picard迭代比完全Newton迭代稳定很多也快很多。3.3 时间积分的刚性问题与ode15s双温模型的刚性非常强激光脉冲期间电子系统特征时间在飞秒量级脉冲结束后晶格弛豫在皮秒到纳秒量级两者差了上千倍。如果自己写显式Euler或RK4时间步长会被一个很小的CFL条件锁死根本算不下去。我在代码里求解随边界条件的隐式方案但更方便的做法是直接调MATLAB的ode15s它对刚性问题有内置自适应算法只要把空间离散后的半离散ODEs写出来就可以获得自动控制的时间步长。需要注意ode15s要求右端函数dudt f(t,u)而我们的质量矩阵依赖Te所以最好把质量矩阵形式传给ode15s或者直接把质量矩阵移到右边使用M的逆乘以其他项。前者的好处是ode15s可以利用稀疏矩阵结构后者则简单直观但可能引入额外数值误差。我推荐使用Mass选项即便配置麻烦一点对长时间积分稳定性更有利。3.4 核心代码一维薄靶的建模仿真下面给出一维金属双温模型的核心代码骨架。这个版本暂时关闭载流子方程先把Te和Tl两个通道跑通再加N方程。% 参数定义 L 200e-9; % 膜厚 200 nm nx 200; % 节点数 x linspace(0, L, nx); dx x(2) - x(1); gamma 68; % 电子热容系数 J/(m^3*K^2) k0 315; % 室温电子热导率 W/(m*K) G 2.2e16; % 电声耦合 W/(m^3*K) Cl 2.5e6; % 晶格热容 J/(m^3*K) T0 300; % 初始温度 % 激光参数 F 0.5; % 能量密度 J/cm^2 - 换算 F F * 1e4; % 转成 J/m^2 tFWHM 100e-15; % 脉宽 100 fs t0 5 * tFWHM; alpha 1e7; % 吸收系数 m^-1 R 0.5; % 表面反射率 I0 0.94 * F / tFWHM; % 峰值强度 % 时间跨度 tspan [0 20e-12]; % 0 ~ 20 ps u0 [T0*ones(nx,1); T0*ones(nx,1)]; % 求解 Moptions odeset(RelTol, 1e-5, AbsTol, 1); [tout, uout] ode15s((t,u) rhsModel(t,u,x,... gamma,k0,G,Cl,I0,tFWHM,t0,alpha,R), tspan, u0, Moptions); function dudt rhsModel(t,u,x,gamma,k0,G,Cl,I0,tFWHM,t0,alpha,R) nx length(x); Te u(1:nx); Tl u(nx1:2*nx); % 电子热导率 ke k0 * Te ./ max(Tl, 300); % 一维有限元近似这里用中心差分矩阵 % 实际应组装稀疏三对角矩阵下式示意 K stiffness1D(x, ke); M mass1D(x, gamma * max(Te, 300)); % 激光体热源 It I0 * exp(-4*log(2)*(t-t0)^2/tFWHM^2); S alpha * (1-R) * It * exp(-alpha * x); % 电子方程空间项和耦合项 dTedt M \ (-K*Te - G*(Te-Tl) S); dTldt (G/Cl) * (Te - Tl); dudt [dTedt; dTldt]; end这段代码里的stiffness1D和mass1D需要根据一维线性有限元组装。不过多解释了实际中每个单元刚度矩阵是[k/h -k/h; -k/h k/h]质量矩阵是h/6乘以[2 1; 1 2]。组装好以后用sparse函数生成稀疏矩阵会大幅提升速度。加上载流子方程时只需要把N作为第三段解向量在rhsModel里添加相应的N扩散项、复合项和光生源项并把复合热源回加到电子方程即可。对半导体Si吸收系数α和载流子产生项S_N α(1-R)It·exp(-αx)/(hν)要同时出现在电子能量方程和N方程里注意单位一个是W/m³一个是m⁻³s⁻¹。4. 结果怎么算、怎么看、怎么验4.1 电子温度与晶格温度的分离与弛豫算出来的第一件事是画Te-Tl随时间演化曲线。典型结果长这样激光脉冲上升沿Te快速攀升在脉冲末端附近达到峰值峰值温度可以到几千甚至上万KTl几乎没有动。随后Te开始下降Tl缓慢上升两者在若干皮秒后逐渐靠近。电子-声子弛豫时间大致可由τ ≈ C_e/G估算金在低激发下大约零点几皮秒到一皮秒铜更快一些。如果Te峰值出现在t_0以前或者Te曲线在激光脉冲里出现锯齿状跳变先怀疑源项时间坐标或时间分辨率问题。把输出时间点加密到脉冲附近基本上每个飞秒一个点曲线就会光滑。理论峰值时间对脉宽不太敏感但对能量密度比较敏感因为C_e随Te增加高注量下电子热容变大温升反而更钝。4.2 载流子密度演化与能量沉积深度打开N方程之后载流子密度在表面附近最高因为吸收源项随深度指数衰减。早期载流子分布基本跟着Beer-Lambert分布走后续扩散会让分布前端逐渐往前推进推进距离由√(D·t)给出。对硅来说百飞秒脉冲激发的N峰值可高达10²⁰甚至10²¹ cm⁻³量级如此大的载流子密度下Auger复合极快本质上改变了早期能量耗散路径。如果把Auger项关掉只留线性复合载流子寿命会远长于实际值仿真的电子温度下降曲线会偏慢。这里建议同时画N沿深度的一维剖面演化云图或者几个时刻的曲线叠加能够很直观地看到载流子峰值位置和扩散前沿。如果表面N出现负值多半是时间步长过大导致扩散项振荡需要缩小时间精度或者给N方程加一点人工扩散稳定性处理。4.3 网格收敛性和时间步敏感性有限元仿真最怕网格无关性没验证就挖结果。我的做法是固定一组物理参数分别用100、200、400个节点跑一遍对比表面电子温度峰值T_e,max。如果200节点和400节点结果差异小于1%基本可以认为网格收敛。对时间步长ode15s自带自适应控制主要检查输出时间点是否覆盖脉冲形状有时输出点间隔过大峰值会被漏掉。网格节点数表面Te峰值 (K)备注1008210偏粗峰有一定过冲2008175可接受4008170基本收敛网格步长的物理约束是热扩散长度。飞秒加热后电子热扩散长度大约为√(D_e·τ_pulse)对金来说大约是几十纳米因此网格步长必须远小于这个值。我一般把膜厚方向网格控制在2~5纳米二维光斑方向反而可以粗一些因为横向温度梯度相对平缓。4.4 用实验数据对标阈值与烧蚀深度仿真结果不能只在笔记本里自嗨最好和实验损伤阈值对照。比如对金薄膜百飞秒脉冲的单发损伤阈值普遍在百mJ/cm²到1 J/cm²量级不同实验室差异很大原因是表面粗糙度、氧化层、厚度和衬底散热条件都不一样。我的经验是把仿真里晶格温度接近熔点的注量作为理论损伤阈值通常落在实验值偏低一侧因为实际损伤还涉及去除、飞溅等复杂机制。如果发现仿真阈值和实验对不上先调G和k_0这两个参数。G决定电子到晶格的能量传输效率直接影响晶格升温速度k_0决定电子热量在内部扩散的速度。损伤阈值和晶格温度强相关所以G的影响极明显。另一个常被忽略的是表面反射率R不同激光波长下R差异很大而且激发过程中R可能变化基础模型里直接用静态R会引入误差。5. 常见数值问题与排查技巧实录5.1 负温度和数值振荡电子热容的低温截断双温模型里最常见的错误是算出一帧Te负值。原因主要是电子热容C_eγTe在Te趋近0时趋近0导致方程退化和强振荡。解决方案简单粗暴给Te一个下限比如300K即C_eγ·max(Te,300)。物理上虽然不完全严谨但工程上能有效避免非物理振荡。更讲究的做法是引入Debye温度附近电子热容的平滑曲线但多数场合用截断就够了。T_l出现负值更少见如果出现检查电声耦合项数值是否有不对称处理。我习惯在整个方程中保持Te-Tl作为一个整体计算不要分开算两部分做减法否则浮点误差在小温差时容易引发振荡。5.2 刚性时间积分与初始爆速ode15s在初始阶段经常走到极端小步长如果发现计算卡住不动先看初始条件附近是否加载了过强的源项。把激光源项改为从零开始平滑爬升或者提前用一段零激光预热时间通常能缓解。另一个办法是设置MaxStep参数比如tFWHM的十分之一避免求解器在时间跨度后端盲目大步长跳过脉冲峰值。如果自己写隐式Euler迭代建议每个时间步的初始猜测用上一步解做线性外推而不是直接用上一步解能减少一到两次迭代。这个技巧在大温度梯度情况下尤其管用。5.3 源项单位与入射边界处理单位问题绝对是重灾区。激光能量密度F如果用J/cm²必须乘1e4换算成J/m²不然所有温度都会偏小几个量级。吸收项Sα(1-R)I(t)exp(-αx)出来是W/m³参与电子能量方程时单位也要一致。At t0时刻虽然激光还没开始但表面反射率和吸收深度没变模型会自然给出初始300K均匀场。如果误把激光能量写成边界热流温度会集中在最表面一两个网格里比体源项结果高出很多而且网格依赖极其严重。我判断边界/体源项是否搞错的方法很简单缩网格一倍看结果是否基本不变。如果表面Te峰值明显变化说明问题处理方式不对。5.4 参数敏感性分析哪些参数不能乱拍给一张排查速查表方便对照现象可能原因检查方法Te峰值过高R偏低、源项单位换算错打印峰值处S单位检查F换算Te峰值过低G偏大、k0偏大缩短弛豫时间τC_e/G验证Tl几乎不动G太小或时间窗口太短延长tspan到50 ps再看N出现负值时间步长过大设MaxStep或缩小AbsTol结果网格依赖网格不够细做两倍网格加密对比曲线锯齿输出点太疏加密tspan尤其在脉冲附近参数敏感性方面G、R、α对表面Te峰值影响最大C_l只影响晶格上升速率k_e对内部温度分布影响明显但对表面峰值影响相对小。建议把参数扫描写成循环跑一个网格图看一下趋势不要只凭一组参数下结论。6. 一点个人体会和后续扩展建议我自己的经验是双温模型在MATLAB里的实现难度不在物理方程而在数值细节。物理方程列出来也就三行真正调试时间几乎全花在单位换算、参数拟合和消除数值振荡上。建议新上手的同行不要直接奔着二维去先弄一个一维薄靶模型把Te和Tl跑顺再加N方程然后才是聚焦光斑的二维模型。每一步都验证好再往上搭否则出了问题你根本不知道是网格的锅还是方程的锅。后续可以考虑加入相变项比如Pennes式等焓法处理熔化和蒸发潜热让模型能预测烧蚀深度也可以把载流子能量项和迁移率随温度变化考虑进去用于半导体激光烧蚀模拟。这个项目标题的框架本身很标准真正体现水平的地方在参数标定和模型扩展。有朝一日你需要做飞秒激光微纳加工工艺优化时这套“双温载流子FEM”的组合能给你省下大量试错成本希望对你的实验和论文都有帮助。