激光辐照下固体热输运的随机算法及应用

龚艳春，武文远，黄雁华，戴斌飞，吴成国

（解放军理工大学理学院，江苏南京 211101）

激光辐照下固体热输运的随机算法及应用

龚艳春，武文远，黄雁华，戴斌飞，吴成国

（解放军理工大学理学院，江苏南京 211101）

采用随机算法对激光辐照下固体目标温度场的分布进行了计算机模拟。该算法将热输运过程作为动力学问题来对待，不考虑具体的传输机制，从而详尽地考虑系统的细节，相比于传统的数值算法，该算法简单明了，占用计算机资源少，普适性强，尤其在处理复杂系统问题时具有独特的优点。

激光辐照；热输运；统计算法

随着激光技术的广泛应用，尤其是激光武器的出现，激光损伤机制和激光防护技术越来越受到人们的关注。研究表明激光束入射到目标体上，破坏目标体的直接因素是激光加热引起的温升，因此研究目标体的激光破坏，必须考虑目标体在激光辐照下温度场的时间和空间分布，从而为激光防护提供依据。激光辐照下目标体的温升取决于目标材料的性质，包括激光吸收系数、激光反射系数、热容率、热传导率等，而这些性质又与材料的微结构和表面状况存在密切的联系，因此研究激光的辐照效应与目标材料特征参数之间的关系，有助于人们寻求新的激光防护材料以及采取新的激光防护措施［1－2］。

关于激光辐照下目标物体温度场的计算，传统的方法是给定激光加载条件（通常把吸收的激光能量或光强作为面热源（表层吸收）或体热源（深层吸收）处理）、目标边界或界面条件、初始温度分布以及物体的热物理特性，将激光加热问题归结为不定常热传导方程的解析描述，通过有限差分（finite-difference）或有限元（finite-unite）的方法进行数值计算［3］。不过，目前大多数详细的计算方法都是针对简单系统设计的，只能处理几何形状和加热方式都比较简单的问题［4－5］。这些方法在处理复杂系统，比如复杂的材料结构、多种材料组分、材料特征参数对温度具有依赖性、空间分布不规则的扩展光源等问题时，往往需要进行各种各样的近似和简化，这些都将影响模拟计算结果的准确性。有些物理问题，比如纳米尺度上的热传输问题，由于很难进行实验测量，计算机模拟就成为基本的研究方法，而计算精度的不确定性将带来很大问题。若将一个物理问题的各种复杂因素均考虑在内，又会使数值计算变得很难处理。对于微观机制不明以致不能写出解析表达式的热过程的模拟计算，传统方法更是显得无能为力。

阐述了一种新的随机算法（stochastic methods），这种算法将热输运过程作为动力学问题来对待，不必考虑具体的传输机制，而是用统计的方法进行模拟。这种算法详尽地考虑系统的细节，而计算过程又相当简化，可以用来处理各种各样的热过程，因而具有普适性，尤其在处理复杂系统问题时具有独特的优点。笔者采用这种方法对固体硅在激光作用下的热传输过程进行了详细的计算，得出了一些有益的结论。

1 计算方法

随机方法最初是由BUNKER和GILLESPIE提出的，用以模拟复杂的化学反应，他们可以快速准确地获得一组微分方程的解，比传统的数值模拟方法要简捷得多。HOULE最先将这种算法用于固体材料在激光照射下的温度变化过程的模拟［6］，其基本思想是将热传输过程看成是一个随机的动力学过程，这一点与Monte Carlo方法类似。

人们知道激光辐照目标时，激光辐射能在目标表面处一部分被反射，一部分折射进入目标体内，使目标材料被加热、熔解以及蒸发，目标上的温度分布随着时间的变化可以用热输运方程来描述：

式中：T是目标体的温度，它是时间t和空间位置的函数；C（T），ρ（T）以及K（T）分别是目标体的热容、质量密度以及热传导率，它们与温度有关；vp是伴随蒸发导致的固体目标表面的退化速率；α（T），R（T）分别是激光吸收系数和表面反射系数，它们也与温度有关；I0（t）是入射激光强度。式（1）实际给出的是单位空间中能量的时间变化率。

按照随机算法，将时间划分成许多小的时间步长，将空间划分成许多子空间，每一子空间具有一定数量的“虚粒子”，其数目取决于该子空间的温度，或者是该子空间中的能量，一个“虚粒子”就相当于一个能量量子，当固体中发生能量的入射或能量的转移时，就用能量量子的转移来代表。

考虑具有轴对称分布的入射激光束，以入射激光束的中心轴线为主轴将目标体按照柱坐标细分，i代表径向的指标，基数用r表示，j代表厚度方向（轴向）的指标，基数用d表示，如图1所示。指标为（i，j）的子空间的体积用V i，j表示，相应的虚粒子数目为qi，j，qi，j的变化将引起该子空间能量的变化。温度随能量量子数qi，j的变化可以用式（2）来表示。

式中：Tinit是初始温度，Q0是“虚粒子”的能量单位，Ci，j是该子空间的热容。这样式（1）的热输运过程就转化成了各单元子空间上能量（或温度）的局部变化问题。

子空间的大小可以根据实际问题的需要精细划分，因而每一子空间的温度将随时间变化，但任一子空间上的温度分布是均匀一致的。同时规定一个子空间上目标体的质量密度相同，只有一

种材料组分甚至是空隙，与温度有关的特征参数如热容、热传导率、光吸收系数等在整个子空间上各向同性、均匀一致。这样统计方法即使是处理具有复杂结构和材料组分的系统，也不会带来更多复杂的计算。各子空间上相应于激光吸收、蒸发以及热传导的事件概率可表示如下（根据图1所示的柱坐标）：激光吸收概率

这里Sr和S h分别是子空间的侧面积和底面积，d r和d h分别是子空间的宽度和高度。

按照Possion统计，两个动力学事件之间的时间步长依赖于瞬时发生的所有事件的可能性，可以表示为

其中R0是介于0与1之间的随机数。

2 计算结果

采用上述算法，本文计算了激光辐照下单晶硅的热效应问题。计算中假定入射激光为高斯型连续激光，光斑半径为3.54μm，光功率可调，环境温度和介质硅初始温度为300 K，表1为常温下硅的热物理参数，各热物理参数随温度变化规律如下：

式中α，γ，β为可调参数，T0是参考温度。

图2为从介质表面离激光光斑中心不同距离处的温度分布，同一时刻介质表面温度场的径向分布与激光场能量分布相似，基本上也呈高斯状分布特征。当介质热导率较大（k＝457 W／（m·K））时，介质温升幅度较小；当介质的热导率较小（k＝157 W／m K）时，温升反而较大。这是因为随着介质热导率的增加，介质在轴向和径向方向的热扩散概率相应增大，热扩散量相对增加，介质光斑中心处从入射激光束吸收的能量能更及时地沿径向和轴向扩散出去，因此温度也相对较低。

表1 Si的热物理参数

图3给出了不同激光辐照功率下，光斑中心处介质轴向上的温度分布。由图3可知，介质轴向温升随着激光辐照功率的增加而增加，介质温度沿着轴向呈指数衰减趋势。

图4、图5给出了5 W激光功率下不同反射率及比热时，介质表面温度的径向分布和表面中心处温度随辐照时间的变化。从图中可以看出，考虑介质热物理性质随温度的变化和不考虑热物理性质随温度变化，能导致结果的较大差异。而实际情况下物体热性质一般都与温度有关，并且温度越高热性质的温度相关性越不容忽视。激光辐照固体材料，尤其是高功率激光辐照固体材料时，能短时间内把物质加热到很高温度，因此通过近似的方法解析求解激光热效应问题，能极大地影响结果的精度。

图4 不同参数下介质径向温度的时间变化

图5 不同参数下介质中心温度的时间变化

3 结 语

在热传导方程的基础上，对比分析了随机模拟方法与传统计算方法的优劣，并在柱坐标系下建立起激光辐照固体材料热效应问题随机模拟的计算模型。由于该模型中各个物理量均可调，而且也可以考虑热物理性质随温度的非线性变化，因此具有普遍的适用性。该方法的突出特点是计算过程简单，却能适用于复杂情形，并且能很容易通过计算机模拟实现最终结果。

以单晶硅材料为例进行了随机模拟计算。结果表明：温升不仅与材料热性质（热容、热导率、密度）密切相关，而且还与材料的吸收系数，激光辐照参数（功率密度、能量分布、光斑大小、作用时间）有关。

［1］刘晓凤，李 笑.SiO2保护膜对高反膜激光损伤特性的改善［J］.强激光与粒子束，2010，22（12）：2 861-2 864.

［2］徐立君，张喜和.强激光辐照硅光电探测器的损伤判别研究［J］.激光与红外，2009，39（12）：1 263-1 266.

［3］孙承伟.激光辐照效应［M］.北京：国防工业出版社，2002.

［4］郑艳丽，杜太焦，束庆邦，等.不同气流环境下激光辐照金属材料温升的数值模拟［J］.强激光与粒子束，2010，22（11）：2 531-2 534.

［5］焦路光，赵国民，陈敏孙.二维情形下组合激光辐照45＃钢靶温升效应的模拟［J］.红外与激光工程，2010，39（1）：42-46.

［6］HOULE F A，HINSBERG W D.Stochastic simulation of heat flow with application to laser-solid interactions［J］.Applied Physics A，1998，66（2）：143-151.

TN24

A

1008-1542（2011）12-0011-04

2011-06-20；责任编辑：陈书欣

龚艳春（1970-），女，江苏连云港人，副教授，主要从事激光物理与技术及电磁环境效应方面的研究工作。