短脉冲激光加热分数阶导热及其热应力研究1)

时间：2024-05-22

许光映王晋宝薛大文

(浙江海洋大学港航与交通运输学院,浙江舟山 316022)

引言

随着激光技术的日益发展,其应用领域从激光的材料加工到激光医学不断拓展,最近又探索应用降低超声速波阻方面[1],其中在激光加工材料过程中材料内部温度场的不均匀性会导致材料内部形成极大的温度梯度,从而产生很高的热应力.因而,准确预测并描述激光加热传导过程温度响应,对于揭示脉冲激光与材料相互作用的物理机制、促进激光技术在现代加工制造业和现代医学中的应用具有重要意义.

通常描述宏观热传导的模型是基于经典Fourier定律.然而,研究[2-4]表明,Fourier 热传导模型不能描述诸如脉冲激光加热之类较大热流、快速加热的情况,因为Fourier 模型隐含着无限大的热传播速度,显然这有悖于物理事实.事实上,对于短脉冲激光,其脉冲持续时间短,一般在纳秒甚至飞秒量级,同时具有高强度热流功率.在这种情况下,短激光束会产生很大的热梯度和较高的加热速率.自从Peshkov[5]在实验中观察到了热波效应,人们普遍认为,在极端低温和快速加热等情况下,热波信号以有限速度传输.因此需要对Fourier 模型进行修正以描述此类瞬态热传导过程.在众多的修正模型中,最为广泛接受的是考虑时滞效应的Cattaneo-Vernotte(C-V)模型[6-7].大量研究[8-10]已经证明了热传播过程中具有波动特性,并表明C-V 模型所描述的热行为能够反映时滞现象,相对于Fourier 模型具有更快的热上升速率.

然而,同时也有一些研究表明,实验测得温度分布与C-V 模型预测的温度分布存在偏差,C-V 模型自身也存在局限性[11-14].还有一些研究表明,CV 模型甚至可能弓入一些不符合物理现实的现象,如局部温度低于绝对零度[15],违反热力学第二定律[16]等.

实际上,对经典Fourier 定律修正的C-V 模型所描述的上述热物理过程是一种超常规的极端过程[17],也是典型的反常扩散过程.从宏观上看,大量粒子在反常扩散过程中的运动不符合标准的统计分布,粒子的均方位移与标准线性行为不符.过去几十年来,数学和物理学界一致认为描述反常扩散过程有效的工具是分数阶微积分.众多研究[18-20]指出,分数阶微积分算子与反常扩散之间存在着内在联系.众所周知,整数阶微分算子是局部算子,而分数阶微分算子是非局部算子.反常扩散所具有的历史依赖与全域相关的特征恰好可以由分数阶导数来描述,这也是分数微积分描述反常扩散现象的最重要优势.

目前为止,分数阶微积分已被成功运用在许多领域,尤其是热传导、扩散、黏弹性等方面.在热弹性方面,Povstenko 基于时间分数阶热传导方程,提出准静态非耦合的热弹性理论,并作了系列研究[21-24].Youssef 等[25-27]将R-L 分数阶积分算子弓入广义热传导方程,建立了耦合的分数阶广义热弹性理论,并通过拉氏和傅氏积分变换研究了二维的热冲击问题以及受坡形热载荷作用的半无限大体的热弹问题.Sherief 等[28]弓入Caputo 型分数阶导数,建立了耦合的分数阶广义热弹性理论.应用该理论,Shweta 等[29]借助拉氏变换和状态空间法研究了半无限大体的热弹响应问题.Sherief 等[30]借助拉氏变换及其数值反变换研究了具有可变热传导率的半无限大体的热冲击问题.此外,Ezzat 等通过采用时间分数阶泰勒级数展开,建立了分数阶双相滞后(DPL)热传导方程和三相滞后(TPL)热传导方程,以此得到分数阶广义热弹性理论[31-32].

王颖泽等[33]基于分数阶扩散理论建立L-S,GL,G-N 的分数型热弹性理论,比较分数阶参数对热波、弹性波传播的特性.马永斌等[34]采用Sherief 的分数阶热弹性理论,研究了含球形空腔无限大体受热冲击作用的动态响应问题.借助拉氏变换及其数值反变换,得到了无限大体中无量纲温度、位移、应力等的分布规律.徐业守等[35]基于Ezzat 等[31]提出的分数阶广义热弹性理论研究对称热冲击作用下三明治板的广义热弹动态响应.朱海陶等[36]基于分数阶广义热弹性理论,利用Laplace 变换及特征值法研究了中空柱的热弹性特征,张培等[37]通过非局部效应和记记依赖微分修正建立广义热弹性理论,研究两端固定、受移动热源作用的有限长热弹杆的动态响应.运用Laplace 变换对控制方程进行求解.

虽然对分数阶热弹性问题已有众多研究,但对于脉冲激光辐照热力效应的研究较少,尤其对于实际制造加工工业过程中的非Gauss 型激光辐照.因此,本文利用分数阶Taylor 展开,以时间分布非Gauss 型激光源作为内热源,建立分数阶Cattaneo 热传导方程,结合动力学方程,研究短脉冲激光辐照下热传导过程和相关的热应力,探讨分数阶阶次和延迟时间及激光源参数对温度场及应力场的影响规律,揭示激光加热反常扩散机制.

1 数学模型

1.1 分数阶Cattaneo 热传导模型

考虑短脉冲激光辐照在金属材料的情况,当短脉冲加热时间短到与材料热化时间相当时,在表面薄层吸收区域,材料对光子吸收过程和能量在微观粒子之间的传递过程时间必须考虑,其间热流的传递过程和温度梯度之间会表现出不同步的滞后效应,通常描述其过程是C-V 模型[38]

式中,q是热流,单位W/m2;τ 是延迟时间,单位s;T是温度,单位°C.

方程(1)可以看成是方程(2)的一阶近似

采用分数阶Taylor 级数展开,可得到分数阶Cattaneo 模型[39]

式中,p是分数阶阶次,0 ＜p＜1,∂p/∂tp是Caputo 型分数阶微分.当p=1 时,方程(3)退回方程(1),为标准C-V 方程,而当τ=0 时,方程(3)退回Fourier定律.

结合能量方程为

式中,ρ 是材料密度,单位kg/m3;CP是比热容,单位J/(kg·K);Q是体积内热源,单位W/m3;αT是热膨胀系数,单位1/K;εV是材料的体积膨胀率,单位1;E是弹性模量,单位Pa;ν 是泊松比,单位1.

将方程(3)代入式(4),得到一般形式的分数阶热传导方程

式中,a=k/(ρcP)是热导率,单位m2/s.显然,当τ=0时,方程(5)退化为Fourier 经典导热方程,当τ ≠0,p=1 时,方程(5)退化为标准C-V 热波方程.

激光与金属材料作用过程复杂,一般认为,激光照射到金属表面后,除了仅约百分之几的能量被表面反射回来,剩下的大部分能量进入材料内部.在内部传播过程中,不断发生散射和吸收,而被吸收的能量转变为热能,用于材料的加热.因此,激光照射金属材料时所生成的热可认为是一个容积吸收过程.由于加热时间短,在垂直于表面方向形成的温度梯度比平行于表面方向的大好几个数量级,因此加热模型简化为如图1 所示一维问题[40],εV=∂ux/(∂x),ux是x方向位移,方程(5a)成为

考虑材料对激光反射和吸收的Lambert-Beer 定律,则体积内热源的生成热为

式中,rf和δ 分别是材料的反射系数和吸收系数.I(t)为热源功率强度,本文选择时间分布非Gauss 型脉冲激光作为热源,强度为

式中,I0是激光峰值功率强度,单位J/m2,f(t)=e-βt-e-µt,β 为激光脉冲上升时间参数,单位1/s;µ为激光脉冲下降时间参数,单位1/s.可知,时间分布函数f(t)取决于β 和µ的不同组合,在本文中,β 和µ的参数取值于参考文献[33-34].其函数分布如图2所示,从中可以看出,由于β/µ为定值,激光峰值强度在量级上相当,但是随着β 和µ的减小,激光强度达到峰值所需的时间增加.显然,该热源强度分布比阶跃函数及三角函数更符合实际情况.

图1 受短脉冲激光加热的半无限体Fig.1 A semi-infinite body subjected to short pulse laser heating

图2 激光热源功率强度时间分布函数f(t)Fig.2 Temporal distribution of the laser power

将方程(6)和(7)代入方程(5b),可得一维分数阶Cattaneo 导热方程

式中I1=(1-rf)I0.

假设材料初始时刻处于准静态热平衡且各处温度均匀,相应的初始条件为

热波效应一般发生在表面很薄的一层,所以问题简化为半空间处理.在忽略表面与外界换热情况下,相应的边界条件为

1.2 热应力模型

考虑表面自由的半空间材料其位移具有下列特征

因此材料几何方程为

由热应力本构关系可得

式中λ=Eν/[(1+ν)(1-2ν)]是拉梅系数,κ=E/[2(1+ν)]是剪切模量.

考虑弹性材料动力学方程

将式(14)对x求导后方程代入式(16)对x求导后方程中,消除变量u,可得

相应的初始条件和边界条件分别为

2 方程解析解

为了便于讨论,弓入以下无量纲量对上述方程无量纲化

2.1 温度场解析解

方程(8)～(11)经无量纲化后得到如下无量纲方程(方便起见,忽略星号)

初始和边界条件的无量纲化为

对上述方程进行Laplace 变换可得

结合边界条件(27)和(28),方程(26)的解为

方程(29)的解析解为

其中方程(36)参考了文献[42]给出的Laplace 逆变换.Eα,β(z)是双参数Mittag-Leffler 函数.其Laplace 变换及其分数阶导数在附录A 给出.

2.2 应力场解析解

同理,对方程(17)～(20)进行无量纲化,得到如下无量纲方程

其中,Λ=aδ/Ve.

对方程(37)～(39)采用Laplace 变换,可得

方程(41)在边界条件(41)时的拉普拉斯变换形式的解为

方程(43)的拉普拉斯逆变换形式为

运用拉普拉斯变换的性质,方程(42)的逆拉普拉斯变换为

其中部分用到的拉普拉斯逆变换在附录A 给出.

3 结果和讨论

为了研究分数阶阶次和激光参数对温度场和热应力的影响,选用金属铜作为研究算例,其物性参数δ=6.16 × 106m-1,a=1.12 × 10-4m2/s,Cp=385 J/(kg·K),ρ=8930 kg/m-3,k=385 W/(m·K),τ=1.0×10-11s,λ=7.76×1011Pa,µ=3.86×1011Pa.对应与计算相关的无量纲参数τ=4.25,Λ=1.66,Ve=4163 m/s.

3.1 解的数值验证

式(33)是任意分数阶温度的级数近似解,文献[44]给出p=1 时的精确解,但激光源不同,为检验解的正确性,在采用本文激光源情况下,图3 给出了表面温度变化的对比,由图3 可知两者是吻合的,由此可知温度场解(33)是正确的.

图3 级数解与精确解温度变化对比Fig.3 Comparison of the series solution and exact solution

3.2 无量纲温度及热应力场随时间、空间分布

图4 显示了p=0.5,β=1 和µ=5 时不同位置处无量纲温度和热应力随时间的变化情况.从图4(a)可以看出,不同位置的温升和降温速率不同.x=0 和x=2 时的温度响应明显快于x=4 和x=6 时的温度响应,距离表面越远,温度响应越慢,说明热波以有限的速度传播.此外,由于脉冲激光强度的衰减,温度的峰值从x=0 到x=6 逐渐降低.图4(b)为相应的热应力场.从中可以看出,当材料表面(x=0)经过加热和冷却阶段时,材料内部受到压缩和拉伸.随着热波信号的传递,不同位置的应力呈现相似的分布.由于压缩作用,热应力首先呈现负值,随后由于拉伸作用,热应力呈现正值,不同位置的热应力绝对峰值出现在加热阶段.

随着热应力波的传播,不同部位的峰值出现在不同的时间.

图4 p=0.5,β=1 和µ=5 时无量纲温度及热应力变化Fig.4 The dimensionless temperature and thermal stress variations at different positions(p=0.5,β=1,µ=5)

图5 显示了p=0.5,β=1 和µ=5 时不同时刻无量纲温度及热应力的空间分布.在图5(a)中,随着激光脉冲经历上升时间和下降时间,材料表面的加热速率随时间变化,材料表面(x=0)历经前期加热后期冷却过程,因而温度出现先上升后下降的变化趋势,图中t=4 时的材料表面温度低于t=2 时的温度,而t=6 时的则更低.然而,虽然表面温度下降了,但是对于材料内部(x＞2),t=4 及t=6 时刻其温度却高于表面温度,表明此时材料内部仍处于加热状态.这种非同步的加热和冷却状态恰恰符合非平衡热传导物理机制的有限热波理论.图5(b)给出了不同时刻热应力分布.图中可以看出,热应力值有正有负,表明材料一部分处于压缩变形状态,而另一部分处于拉伸状态.如上所述,加热阶段中,材料由于受热膨胀导致受到压应力,与之相反,冷却过程材料受到拉应力.随激光脉冲完成,应力峰值随时间变化有所增大.

图5 不同时刻温度及热应力分布Fig.5 The dimensionless temperature and associated stress distribution at different times

3.3 激光参数对温度场及应力场的影响

图6 不同激光参数下不同时刻材料温度分布图Fig.6 The dimensionless temperature distribution for different laser pulse parameters

图6 显示了不同激光参数作用下,不同时刻金属材料的温度分布.如1.1 节所述,β 和µ代表激光脉冲上升和下降时间,随着β 和µ的减小,激光强度达到峰值所需的时间增加.从图5(a)可以发现,在t=0.5之前,β 和µ值越大,温度越高.这是因为在t=0.5 之前,所有工况下激光脉冲都处于上升阶段(见图2),β和µ值越大,表明激光辐照能量释放越快,显然材料所受加热速率越快,温度越高.同时从图中可知,温度分布曲线斜率表示当地温度梯度,斜率越大,梯度越大.形成β=1,µ=5 工况温度最高、斜率最大,而β=0.125,µ=0.625 工况温度最低,斜率最小的分布情形.而对比图6(a)和图6(b),当t=1.5 时,虽然各工况下温度较t=0.5 时刻均有所升高,但β=1,µ=5工况下温度明显低于β=0.5,µ=2.5 和β=0.25,µ=1.25 工况.表明β=1,µ=5 工况下激光加热速率最先减慢,这与图2 中β=1,µ=5 工况最先出现脉冲下降阶段吻合.在β=1,µ=5 工况处于脉冲下降阶段时,其余工况仍处于脉冲上升阶段,从而导致材料表面温升速度快于前者形成图中所示温度分布情形.而随着时间推移,各工况激光相继经历脉冲下降阶段,加热速率先后减慢,形成不同时刻,不同工况下的材料表面温度相继领先的分布情形(图6(c)～图6(f)).通过上述分析也就不难理解图6(a)和图6(f)中不同工况下温度分布和温度曲线斜率呈现相反趋势,这与激光脉冲的时间分布特征密切相关.

图7 不同激光参数下不同时刻材料热应力分布图Fig.7 The dimensionless thermal stress variations for different laser parameters

对应于图6 中的温度分布,图7 是不同激光参数工况下热应力分布图.在图7(a)中,在激光脉冲加热开始时,材料受热区由于膨胀导致受到周围材料挤压作用,因而图中热应力值均为负值,且激光脉冲上升和下降时间越大,热应力幅值越大.随β=1,µ=5 工况下激光脉冲经历下降阶段,温升速率减慢,温度梯度减小,相应热应力幅值减小,因而在图7(b)中呈现β=1,µ=5 工况热应力幅值小于β=0.5,µ=2.5,β=0.25,µ=1.25 工况.图7(c)中β=1,µ=5 工况下出现正热应力值,表明此时靠近材料表面的区域处于收缩阶段,这与图7(c)中该工况下温度降低相符.随着时间的推移,其他工况相继进入脉冲下降阶段,加热速率逐渐减慢,先后经历冷却阶段,形成正负热应力共同存在的情形,表明拉应力和压应力共同出现,材料部分受拉、部分受压.同样,这种非同步冷、热状态证实了非平衡热传导的有限热波理论.

3.4 分数阶阶次对温度场及应力场的影响

图8 显示了x=0 和x=3 位置不同分数阶阶次p对无量纲温度变化的影响.众所周知,局部温度与局部能量积累有关,热传递速度越慢,能量传递也越慢,因而在当地位置积累的能量越多,温度就越高.由于Fourier 模型暗含着无限的热传播速度,因而Fourier 模型中,材料表面受热(x=0)的升温速率及幅值是最小的,因为其能量以极快的速度传播至材料内部.标准C-V 模型(p=1)表明热信号以有限速度的热波形式传播,因而标准C-V 模型(p=1)预测的材料表面温度和升温速率均高于Fourier 模型.至于分数阶模型,其介于Fourier 模型和标准C-V 模型之间,呈现在温度梯度延迟时间小于热流延迟时间时DPL 模型特征[43].随着分数阶阶次p的增大,材料表面温度和升温速率增大,表明热波速度随分数阶阶次增大而减小.

图8 不同分数阶阶次下x=0 和x=3 处温度分布Fig.8 The temperature variations with different fractional order p at x=0 and x=3

在材料内部(x=3),不同模型的温度响应速度与图6(a)确相反,Fourier 模型预测的内部温度响应最快,而标准C-V 模型(p=1)温度响应最慢.这仍然体现了不同的热波传播速度.如上所述,Fourier 模型意味着无限的热传播速度,因此Fourier 模型中,热信号首先传递到x=3 的位置,从而产生最快的局部温度响应.相反,标准C-V 模型中,热波速度最慢,因此体现了最慢的局部温度响应.而分数阶模型预测结果介于标准C-V 模型和Fourier 模型之间.

图9 不同分数阶阶次x=1 和x=3 处热应力分布Fig.9 The effects of the fractional order p on the thermal stress at x=1 and x=3

图9 不同分数阶阶次x=1 和x=3 处热应力分布(续)Fig.9 The effects of the fractional order p on the thermal stress at x=1 and x=3(continued)

图9 是x=1 和x=3 位置不同分数阶阶次p对热应力变化的影响.注意到材料的左表面都是自由边界(方程19),其应力为0.因而这里讨论x=1(而不是x=0)和x=3 时的热应力变化.如上所述,热波的到达导致局部温度升高,使当地产生局部压应力,而温度降低则导致局部拉应力.如图9(a)所示,不同模型的热应力变化相似.但热应力峰值随分数阶次的增大而增大,标准C-V 模型预测的热应力幅值最大,这与图8(a)呈现的温度分布趋势一致.图9(b)中热应力分布与图9(a)相似.但图9(b)中不同模型的热应力峰值出现的时间有所差别,表明不同的模型体现不同的热波速度.