使用comsol仿真软件 利用双温方程模拟飞秒激光二维移动烧蚀材料 可看观察温度与应力分布 周期为10us,变形几何部分本人还在完善学习中 三维的也有 还有翻阅的论文文献一起打包
最近折腾飞秒激光加工仿真搞得头大,特别是那个材料烧蚀过程中电子和晶格的热传递过程。COMSOL里双温方程建模说难不难,但二维动态移动热源配置起来真能让人薅秃头发。
先丢个核心方程镇楼:
% 双温方程COMSOL弱形式 test(Te_t) = (G*(Tl-Te)) + ∇⋅(k_e*∇Te) + Q_laser; test(Tl_t) = (G*(Te-Tl)) + ∇⋅(k_l*∇Tl);这俩方程看着简单,实操时电子热导率ke和晶格热导率kl的温度依赖性能让计算疯狂发散。之前用铜做测试,参数表直接写成这样:
k_e = 300*(Te/300)^(-0.5); % 电子热导率随温度变化 G = 1e17*(1 + 0.2*exp(-Te/1e4)); % 电子-声子耦合系数结果迭代到第三个脉冲就崩了,后来发现是晶格温度梯度太大导致网格畸变。现在改成移动网格配合变形几何模块,边界条件这么设:
physics.create('defgeom', 'DeformedGeometry', 'geom'); physics.feature('defgeom').feature('dmm1').set('d', {'u' 'v'}); physics.feature('defgeom').feature('dmm1').set('shape', 'free');移动热源这块有个骚操作——用解析函数实现扫描路径:
% 飞秒激光光斑移动函数 v = 0.5; % 扫描速度m/s Q_laser = @(x,y,t) P0*exp(-((x-v*t)^2+y^2)/(2*r0^2)).*heaviside(t-n*T).*heaviside((n+1)*T-t);参数扫掠时发现周期10μs的脉冲积累效应明显。某次作死把脉冲重叠率调到80%,瞬间得到火山口状的熔池轮廓——跟去年Nature子刊里那个钛合金烧蚀形貌神似。
目前卡在应力场的动态显示上,特别是热应力与相变应力的耦合。翻到篇2019年的APL论文提到用增量应力更新算法,正在尝试移植到COMSOL的固体力学模块:
model.component('comp1').physics('solid').feature('i1').set('s', {'sX' 'sY' 'sZ' 'sXY'}); model.component('comp1').physics('solid').feature('dwe1').set('weak', 'dw_lin_elasticity_iso_planestress');三维模型更刺激,光边界条件就设了六个不同热流条件。有个坑爹现象:当Z方向尺寸超过50μm时,边缘会出现鬼影应力集中,估计是网格映射时雅可比矩阵出问题。
最近在憋大招想把热-力耦合跟流体模块联动,模拟熔融金属的飞溅过程。不过变形几何模块的ALE算法还没玩透,每次计算到7μs左右就报错"Jacobian matrix has non-positive determinant",得手动重置网格——这酸爽,谁做谁知道。
(文献包里的那篇《Ultrafast laser ablation of copper: a molecular dynamics and experimental study》简直是救命稻草,里面提供的电子弛豫时间参数让收敛速度快了三倍不止)