基于双温模型与有限元法的载流子密度与电子晶格温度模拟研究：飞秒激光源下的德鲁德模型应用

带载流子密度的双温模型matlab电子晶格温度电子密度飞秒激光源模拟有限元法解偏微分方程。 德鲁德模型带载流子密度变化。当飞秒激光哐哐砸在金属表面时电子和晶格开始上演冰火两重天的戏码。这里咱们要搞个能同时追踪电子温度、晶格温度还有载流子密度变化的双温模型这事儿用Matlab折腾起来特别带劲。先瞅瞅核心方程长啥样% 电子温度方程 C_e * dTe/dt ∇(k_e * ∇Te) - G*(Te - Tl) Q_laser % 晶格温度方程 C_l * dTl/dt G*(Te - Tl) % 载流子密度方程 dn/dt α*I(t) - β*n^3 ∇(D*∇n)这几个方程里藏着电子-晶格耦合、激光热源、载流子复合这些关键机制。特别是载流子密度那项德鲁德模型的修正项会让热导率k_e变成n的函数这就让方程组直接耦合上了。咱们用有限元法来处理空间离散时间推进用Crank-Nicolson比较稳。先定义个环形激光光斑试试laser_profile (x,y,t) exp(-((x-x0).^2 (y-y0).^2)/(2*sigma^2)) .* exp(-(t-t0).^2/(2*tau^2));这个高斯时空分布能比较好模拟飞秒激光的时空特性。注意时间项的脉宽tau要设到百飞秒量级空间分布sigma根据聚焦光斑调。网格生成方面用Matlab自带的generateMesh就行但得注意在激光作用区域加密model createpde(); geometryFromEdges(model,circleg); generateMesh(model,Hmax,0.1,Hgrad,1.5); [p,e,t] meshToPet(model.Mesh);这里Hgrad控制网格渐变避免在加密区产生突变网格影响计算稳定性。带载流子密度的双温模型matlab电子晶格温度电子密度飞秒激光源模拟有限元法解偏微分方程。 德鲁德模型带载流子密度变化。核心求解循环里藏着几个trickfor n_step 1:Nt % 更新载流子密度显式处理非线性项 n_new n_old dt*(alpha*I_now - beta*n_old.^3 D*laplacian(n_old)); % 更新电子温度隐式处理扩散项 A_Te assembleFEMatrix(C_e/dt G, k_e(n_new), ...); b_Te assembleRHS(G*Tl_old Q_laser); Te_new A_Te\b_Te; % 更新晶格温度显式足够 Tl_new Tl_old dt*(G/(C_l))*(Te_old - Tl_old); end注意载流子密度方程里的三次复合项β*n³如果全隐式处理会很难解这里用显式反而更高效。而电子温度的扩散项必须隐式否则时间步长会被限制得很小。后处理时有个好玩的现象——电子温度会在激光脉冲结束后继续上升一小会儿这是载流子密度变化导致热导率降低引起的figure; subplot(2,2,1); pdeplot(p,e,t,XYData,Te_history(:,100),Contour,on); title(电子温度(100fs));当画出三维时空演化图时会看到温度波从中心向外传播的涟漪效应而载流子密度分布则呈现火山口状——中心因为复合速率快反而密度更低。调试这种模型时最容易翻车的地方在参数量纲。有一次我把电子热容单位搞错成eV/(m³·K)结果温度瞬间飙到1e8开尔文直接蒸发了整个模型。后来加了无量纲化处理才踏实% 无量纲化处理 T_star 1e4; % 参考温度10000K t_star 1e-12; % 参考时间1ps x_star 1e-6; % 空间尺度1μm这么做不仅数值稳定还能更直观看出各物理量的相对大小。最后要提醒的是这种强耦合问题最好用Jacobian-Free Newton-Krylov方法处理非线性但对于桌面级算例用我说的分步解法已经能跑出不错的结果。毕竟咱们搞物理模拟的总不能为了绝对数值精度牺牲了探索的乐趣不是