数值求解纳米尺度热传导分数阶抛物两步模型

沈淑君

(1.华侨大学 数学科学学院，福建 泉州 362021；2.华侨大学 计算科学福建省高校重点实验室，福建 泉州 362021)

纳米传热现象的模拟已经引起人们的广泛关注，特别是超短脉冲激光加热引起的纳米尺度热传导.由于效率高、功率密度高、附带材料损坏小、加热高精度控制等优点，超短脉冲激光加热技术已广泛应用于生物、化学、医学、物理和其他的热加工领域，如薄金属薄膜的结构监测，激光微加工，激光微薄膜的制图、结构裁剪和薄膜沉积中的激光加工等[1].分数阶微分算子具有非局部性，非常适用于描述现实世界中具有记忆及遗传性质的材料，已成为描述各类复杂力学和物理行为的重要工具之一[2].分数阶微积分已成功地应用于多种热传导模型的模拟和研究[3-5]，分数阶模型已被证实是一个很好的纳米级传热的候选者[6].继续研究文献[1]中的分数阶抛物型两步超短脉冲激光加热引起的纳米级热传导模型，即

(1)

(2)

(3)

1892年，Hadamard提出了Hadamard分数阶导数[10]，这种分数阶导数与Riemann-Liouville和Caputo分数阶导数的区别在于积分核含有任意指数的对数函数.研究者发现，Hadamard分数阶微积分在描述材料疲劳裂纹扩展方面具有潜在的应用价值[11-14].但是，Riemann-Liouville导数和Hadamard导数的其中一个缺点就是对常数的导数一般不等于零.2012年，Jarad等[7]修改了Hadamard分数阶导数，使其对一个常数的导数为零，并且使其在物理上具有类似于Caputo分数阶导数那样的可解释的初始条件，由此产生了Caputo-Hadamard分数阶导数.比起Hadamard分数阶导数，Caputo-Hadamard分数阶导数更适合应用在工程和科学中[7，10，14].

与文献[1]做法一样，引入无量纲参数，便得到无量纲的分数阶抛物型两步热传导模型为

(4)

(5)

式(4)中：Kn表示Knudsen数.

此外，需要指出的是，在考虑纳米级传热时[15]，非跳温边界条件忽略了边界声子散射的影响，导致在边界附近的结果令人不满意.通过观察发现，碳纳米管的热流是阻塞的，温度在管的两端可以看到跳跃[16-17].为了捕捉纳米几何结构内部的边界声子散射效应，提出温度跳跃(Robin′s)边界条件[1]，即

(6)

为使记号简单，采用u(x，t)和v(x，t)分别代替模型中的Te(x，t)和Tl(x，t).

1 模型方程

考虑无量纲分数阶两步模型，即

(7)

(8)

带有Robin边界和初值条件为

(9)

u(x，a)=T1(x)，v(x，a)=T2(x)，x∈[0，1].

(10)

(11)

同理，用类似的过程可以得到v(x，t)的第2个边界条件为

(12)

2 紧有限差分格式

引理1[1]如果f(x)∈C6[x0，xM] ，则有

引理2[9]设0<α<1，g(t)∈C2[a，tk]，则有

(13)

(14)

(15)

(16)

把边界条件(9)和方程(16)一起代入方程(14)，可得

(17)

(18)

于是，方程(17)变为

(19)

同理，用类似上面的过程可以推导出当i=M时方程(7)的紧格式为

(20)

当1≤i≤M-1时，利用引理1，2可得

(22)

(23)

(24)

(25)

(26)

(27)

式(23)～(26)中:1≤k≤N.

备注2用类似文献[1]中的能量分析的方法可以证明差分格式(23)～(27)是稳定的，并且收敛阶为O(τ2-α+h4).

3 数值算例

考虑带Robin边界的Caputo-Hadamard分数阶抛物型两步模型，即

(28)

(29)

(30)

u(x，1)=0，v(x，1)=0，x∈[0，1].

(31)

式(28)，(29)中:0<α<1.

当时间tF=2，空间步长h=1/50时，随着时间步长τ的减少，取不同时间分数阶α时，例1数值解的最大误差和时间上的收敛阶的情况，如表1～4所示.表1～4中：EU，max，EV，max分别为U，V的最大误差；RateU，RateV分别为U，V的收敛阶.

表1 α=0.2 时的最大误差和时间上的收敛阶(例1)Tab.1 Maximum errors and convergence rates in time for α=0.2 (example 1)

表2 α=0.5 时的最大误差和时间上的收敛阶(例1)Tab.2 Maximum errors and convergence rates in time for α=0.5 (example 1)

表3 α=0.8 时的最大误差和时间上的收敛阶(例1)Tab.3 Maximum errors and convergence rates in time for α=0.8 (example 1)

表4 α=0.999 时的最大误差和时间上的收敛阶(例1)Tab.4 Maximum errors and convergence rates in time for α=0.999 (example 1)

由表1～4可知：时间上的收敛阶大概是2-α.这与引理2的结论吻合.

当时间tF=2，时间分数阶α=0.5，τ=1/20 000时，随着空间步长的减少，例1数值解的最大误差和空间上的收敛阶，如表5所示.由表5可知：空间上的收敛阶大概是4阶.

表5 α=0.5时的最大误差和空间上的收敛阶(例1)Tab.5 Maximum errors and convergence rates in space for α=0.5 (example 1)

源项

参数取为B=G=c0=Kn=1.可以算出精确解为

当时间tF=2，h=1/50时，随着时间步长的减少，取不同的时间分数阶α时，例2数值解的最大误差和时间上的收敛阶，如表6～9所示.

表6 α=0.2 时的最大误差和时间上的收敛阶(例2)Tab.6 Maximum errors and convergence rates in time for α=0.2 (example 2)

表7 α=0.5 时的最大误差和时间上的收敛阶(例2)Tab.7 Maximum errors and convergence rates in time for α=0.5 (example 2)

表8 α=0.8 时的最大误差和时间上的收敛阶(例2)Tab.8 Maximum errors and convergence rates in time for α=0.8 (example 2)

表9 α=0.999 时的最大误差和时间上的收敛阶(例2)Tab.9 Maximum errors and convergence rates in time for α=0.999 (example 2)

由表6～9可知：时间上的收敛阶大概也是2-α，这与引理2的结论吻合.

当时间tF=2，α=0.2，τ=1/20 000时，随着空间步长的减少，例2数值解的最大误差和空间上的收敛阶，如表10所示.由表10可知：空间上的收敛阶也大概是4阶.

表10 α=0.2 时的最大误差和空间上的收敛阶(例2)Tab.10 Maximum errors and convergence rates in space for α=0.2 (example 2)

4 结束语

对无量纲的分数阶抛物型两步热传导模型提出了高阶有效的数值算法.方程中采用的是Caputo-Hadamard分数阶导数，并且考虑的是Robin边界条件，目前国内外文献还鲜有研究.利用空间四阶紧格式和Caputo-Hadamard时间分数阶导数的L1逼近格式建立了数值格式.文中的两个数值例子验证了提出的数值算法的有效性.通过改变Knudsen数和分数阶导数的数值以及边界条件中的参数，模拟可以作为分析超短脉冲激光作用下多孔介质(如多孔金属薄膜)纳米尺度热传导的工具.后续可进一步深入探索算法的理论研究.