超快激光加热薄膜的格子玻尔兹曼模拟

毛煜东吕慧丽孙浩森陈明九于明志

(山东建筑大学热能工程学院，山东济南250101)

0 引言

近年来，激光加热技术已展现出强大优势，并广泛应用于激光焊接、激光钻孔、激光熔覆、激光照明、三维印刷技术以及熔化动力学等诸多领域[1-5]。由于该技术具有高效率、高功率密度、最小的附带材料损坏、较低的烧蚀阈值等优点而受到了越来越多的关注。目前该技术在材料的热加工过程中具有很好的应用，如激光清洗、激光图案化以及超音速激光沉积等[6-8]。现代微加工技术要求微结构的高精度加工生产能力，以及材料热处理中加热时间迅速和位置的高精度控制，探索微结构的传热过程是十分重要的。1822年，傅立叶经实验研究导热过程并对其归纳总结，得到了描述宏观导热过程的基本定律—傅立叶导热定律，并逐步地应用到了机械、冶金、建筑以及电气等领域[9]。然而，傅立叶定律暗含热以扩散的方式传播且热传播速度无限大，即如果对导热体内某一点施加热扰动，则其他部分会同时感受到此热扰动带来的温度变化，即傅立叶定律忽略了热传导过程中温度梯度和热流矢量之间的弛豫时间。

玻尔兹曼输运方程BTE(Boltzmann Transport Equation)用于解决微尺度系统中的热输运问题[10]是十分有效的手段。BTE是统计力学中用以描述不平衡态分布函数演化规律的方程，但BTE中含有的碰撞算子项是一个非常复杂的非线性微分积分，极难以求解。格子玻尔兹曼方法 LBM(Lattice Boltzmann Method)是BTE常用的数值解法之一，LBM依据非平衡统计力学的运动论，在解决边界，实现编程和处理复杂的多尺度耦合问题方面具有明显的优势，是求解BTE的有效手段，受到了广泛的关注[11-12]。基于LBM而来的格子BGK模型LBGK(Lattice Bhatnagar-Gross-Krook)简化了BTE中的碰撞算子项，用弛豫时间近似表达了不同粒子间的碰撞达到的复杂机制，并且不考虑外力的影响。由于计算效率高，能精准地推导出纳维—斯托克斯方程N-S(Navier-Stokes equations)，常用于研究不可压缩流体和微尺度的导热问题，也是LBM主要的研究模型。

近年来，LBM广泛应于描述微纳米尺度传热的问题中。Escobar等[13]应用LBM研究了一维半导体材料中的瞬态多空间尺度和时间尺度的导热过程，与基于扩散-弹道输运方程、傅立叶导热方程和CV模型的结果做了对比，并对主要的导热模型的使用范围做了简要总结。华钰超等[14]用热质理论的观点分析了弹道扩散导热机理，基于声子玻尔兹曼方程推导了修正边界条件模型，数值求解了修正的普适导热定律，并与蒙特卡罗模拟进行对比。张珂等[15]基于LBM和双曲两步模型结建立了一个LBM两步方程，并利用此方法模拟研究了纳米薄膜在超快激光照射过程中的热响应特性，分析了照射过程中薄膜内温度随时间及空间的变化规律，探讨了激光强度以及薄膜厚度对金属薄膜热响应的影响。文章应用LBM方法对超快速激光加热硅薄膜的一维导热问题进行研究，来描述这种时间超短，尺寸超小的超快速纳米传热问题，展示薄膜内部能量密度分布情况，并对结果进行分析。

1 理论分析

LBGK模型的表达式由方程(1)表示为

式中:f(r，ξ，t)为分子速度分布函数；r为空间位置矢量；ξ为分子速度矢量；t为时间参数；feq为热力学平衡态时的速度分布函数；v为热载子的群速度；τ0为松弛时间(也叫弛豫时间)是粒子发生2次碰撞的时间间隔，表征粒子碰撞后达到平衡态的快慢程度。

1.1 LBM方法验证

傅立叶导热定律可以很好地描述宏观问题，文章采用LBM方法解决微观—宏观问题，应先验证LBM在解决宏观热传导问题时的有效性。建立一个简单的薄膜导热问题，即在薄膜左侧施加一个无量纲数为1的热扰动，其控制方程由式(2)表示为

式中:Ui为能量密度；x为位置参数；U0i为平衡能量密度；vi为声子群速度；τ为声子弛豫时间。

考虑常温初始条件，将式(3)～(6)所示的无量纲变换为

式中:x∗为无量纲位置参数；L为薄膜的特征长度；Cv为硅的体积热容；t∗为无量纲时间参数；U∗为无量纲能量密度参数；U1∗为向x轴正方向传输的无量纲能量密度；为向x轴负方向传输的无量纲能量密度；T0为初始时刻温度；E1为初始温度T0下的能量。

在做与式(3)～(6)相同的无量纲变换后，其边界条件由式(7)表示为

1.2 应用LBM模拟超快激光加热薄膜问题

建立基于LBM的一维热传导模型来模拟超快激光在硅薄膜中的能量分布。在硅薄膜中，声子是主要的能量载体。基于BTE，考虑超快激光加热薄膜问题的一维导热模型由式(8)表示为

式中:i=1，2；S为能量吸收率，可由式(9)[16]表示为

式中:J为激光能量发射密度；R为表面反射率；tp为激光脉冲的持续时间；δ为激光穿透深度。

使用LBM离散式求解方程(8)，忽略高阶项，并做与式(3)～(6)相同的无量纲变换后，则主导方程可由式(10)和(11)表示为

式中:W=Δt∗=Δx∗/Kn；ξ=L/δ；γ=τ/tp；W、ξ、γ均为中间参数，无特殊含义Kn=l/L为克努森数；l为声子平均自由程。

根据傅立叶定律，将上述的超快激光加热问题由经典热传导方程式进行描述，由式(12)表示为

式中:kF为硅的导热系数。

基于相同边界条件，可将函数的解展开，由式(13)～(15)表示为

1.3 应用LBM模拟双侧超快激光加热薄膜问题

利用相同的超快激光脉冲同时从硅薄膜的两侧表面进行加热，并设置与单侧加热相同的镜面反射边界条件和常温初始条件。此时描述该问题的BTE控制方程由式(16)～(18)表示为

将式(17)和(18)代入式(16)中，可将控制方程变为式(19)形式

利用LBM对式(19)进行离散求解，方法与式(10)和(11)相同，不再赘述。

2 数值结果与分析

为便于分析，超快脉冲激光加热的硅膜采用以下参数:体积热容Cv为1.66×106J/(m3·K)，硅的声子平均自由程l为41 nm，声子弛豫时间τ为6.53 ps，持续时间tp为 0.65 ps，激光通量J为312 J/m2，穿透的光学深度δ为15.3 nm，以及激光的表面反射率R为0.93。

2.1 验证LBM在微观—宏观问题的有效性

分析式(2)～(7)描述的薄膜导热问题，分别利用LBM和经典傅里叶热传导理论对该问题进行求解和模拟，观察薄膜内部的能量传递现象，如图1所示。可以看出，通过LBM获得的结果清楚地表现出扩散行为，并显示出良好的扩散性能，这与傅立叶导热理论所给出的结果一致。因此，LBM可用于模拟宏观扩散状态下的热传导。

图1 利用LBM和傅立叶定律得到的Kn=0.001的无量纲能量密度分布图

改变克努森数，将热传导过程引入到过渡区，并观察LBM和傅立叶定律得到的结果之间的差异，如图2所示。当克努森数为1时，此时薄膜内热传递处于过渡区，基于LBM和傅立叶定律的解决方案表现出截然不同的行为。可以看出，当无量纲时间为0.1时，通过LBM模拟获得的能量密度分布具有显著的特征:能量密度在受加热侧(左侧)边界有明显的跳跃；随着时间的流逝，由LBM获得的热能以热波的形式在薄膜内部传输；当t∗达到1时，能量密度又在右侧边界发生了的跃变。然而，传统的傅立叶定律却不能捕捉这一现象，由于其内含热传递速度无限大的假设，在过渡区，薄膜内任何一处都可以立即感受到受加热侧施加的热扰动。

图2 利用LBM和傅立叶定律得到的Kn=1的无量纲能量密度分布图

2.2 单侧激光加热纳米薄膜导热问题

分析式(8)和(9)所描述的利用超快激光加热厚度为82nm薄膜的导热问题。通过傅立叶定律和LBM方法得到的无量纲能量密度分布如图3所示。从图3(a)中可以看出，基于傅立叶定律的热传导是一个扩散过程。由于强大的激光脉冲加热，短时间内能量密度就在薄膜的左侧区域迅速向右侧扩展。随着时间的流逝，能量密度逐渐降低，并最终在薄膜内趋于稳定。

图3 激光加热薄膜单侧时薄膜内无量纲能量密度分布图

然而，对于相同的加热问题，LBM却得到了不同的结果，如图3(b)所示。在初始阶段，同傅立叶定律获得的结果相比，受到脉冲激光扰动的左边界附近产生了更高的能量密度。该无量纲密度峰值达到了1.18，较傅立叶结果的0.8相比，超出约47.5%。这表明，传统的傅立叶热定律在描述这此类微纳米尺度下的超快加热问题时会严重低估热载子传输的最大能量密度。

随着时间的流逝，薄膜内部能量是以波的形式从左端传递到右端，且随着能量波的传递，其能量密度峰值逐渐下降，这一现象并未出现在傅立叶定律得出的结果中。LBM的结果还发现，在镜面反射的边界条件下，能量密度会经过一个较傅立叶结果更为漫长的时间后才能在薄膜内趋于稳定。

2.3 双侧激光加热纳米薄膜导热问题

分析式(16)～(19)所述利用相同的超快激光脉冲同时从硅薄膜的两侧表面进行加热的问题。在克努森数为0.5时，利用傅里叶定律和LBM方法得到的薄膜内部的能量密度分布情况如图4所示。如同预期的那样，傅立叶的结果显示了热传导的扩散行为，也并未发现能量波状传递的行为，如图4(a)所示。初始时刻，左右两侧的无量纲能量密度峰值均为约0.8，同单侧激光加热时相同，经过无量纲时间2后，薄膜内部能量密度达到稳定，此时无量纲能量密度值为0.5。

由图4(b)可以看到，在初始时刻受激光脉冲扰动的两个边界附近的区域能量密度的快速升高(即温度快速提升)；随后，从左右边界发起的2组热波从不同方向同时向硅膜的中间区域传输，在此期间，能量密度峰值逐渐降低。当无量纲时间达到1时，2组热波峰在薄膜的中间区域汇合，相互的碰撞作用导致了薄膜中心区域的能量密度显着升高。此时，在中心点处的能量密度的峰值为1.095，约为图3(b)所示的单侧激光加热情况下所预测的峰值的两倍，这一现象不仅是傅立叶结果无法捕捉的，而且在单侧激光加热时也不会出现。最后，当无量纲时间达到16时，薄膜内能量密度趋于稳定；相比于傅立叶结果，在镜面反射的边界条件下，趋于稳定所需时间要更为漫长。

图4 激光加热薄膜两侧时薄膜内无量纲能量密度分布图

3 结论

文章应用LBM方法和傅立叶定律对超快激光加热硅薄膜的一维导热问题进行研究，分析了薄膜内部能量密度分布情况。主要结论如下:

(1)发现在克努森数较小的连续介质区，LBM得到的结果符合傅立叶定律给出的结果；然而在系统特征尺度和分子平均自由程相当的过渡区，LBM方法可以捕获到能量以波动的方式进行传输的热波特性，但传统的傅立叶定律无法得到。

(2)利用LBM数值模拟激光加热硅薄膜的一维导热问题，发现克努森数对薄膜内无量纲能量密度的分布具有重要的影响。

(3)如果利用超快激光分别在薄膜的两侧进行加热，薄膜内部的能量密度分布情况将有所改变，尤其是当沿相反方向传播的两个热波交汇在一起时，会发生碰撞并使得能量密度有显着的提升，这些现象均在傅立叶定律结果中不能呈现。此外，研究发现传统的热传导理论可能严重低估了过渡区声子输运的最大能量密度。