數(shù)值求解納米尺度熱傳導分數(shù)階拋物兩步模型

沈淑君

(1.華僑大學 數(shù)學科學學院，福建 泉州 362021；2.華僑大學 計算科學福建省高校重點實驗室，福建 泉州 362021)

納米傳熱現(xiàn)象的模擬已經(jīng)引起人們的廣泛關注，特別是超短脈沖激光加熱引起的納米尺度熱傳導.由于效率高、功率密度高、附帶材料損壞小、加熱高精度控制等優(yōu)點，超短脈沖激光加熱技術已廣泛應用于生物、化學、醫(yī)學、物理和其他的熱加工領域，如薄金屬薄膜的結構監(jiān)測，激光微加工，激光微薄膜的制圖、結構裁剪和薄膜沉積中的激光加工等[1].分數(shù)階微分算子具有非局部性，非常適用于描述現(xiàn)實世界中具有記憶及遺傳性質的材料，已成為描述各類復雜力學和物理行為的重要工具之一[2].分數(shù)階微積分已成功地應用于多種熱傳導模型的模擬和研究[3-5]，分數(shù)階模型已被證實是一個很好的納米級傳熱的候選者[6].繼續(xù)研究文獻[1]中的分數(shù)階拋物型兩步超短脈沖激光加熱引起的納米級熱傳導模型，即

(1)

(2)

(3)

1892年，Hadamard提出了Hadamard分數(shù)階導數(shù)[10]，這種分數(shù)階導數(shù)與Riemann-Liouville和Caputo分數(shù)階導數(shù)的區(qū)別在于積分核含有任意指數(shù)的對數(shù)函數(shù).研究者發(fā)現(xiàn)，Hadamard分數(shù)階微積分在描述材料疲勞裂紋擴展方面具有潛在的應用價值[11-14].但是，Riemann-Liouville導數(shù)和Hadamard導數(shù)的其中一個缺點就是對常數(shù)的導數(shù)一般不等于零.2012年，Jarad等[7]修改了Hadamard分數(shù)階導數(shù)，使其對一個常數(shù)的導數(shù)為零，并且使其在物理上具有類似于Caputo分數(shù)階導數(shù)那樣的可解釋的初始條件，由此產(chǎn)生了Caputo-Hadamard分數(shù)階導數(shù).比起Hadamard分數(shù)階導數(shù)，Caputo-Hadamard分數(shù)階導數(shù)更適合應用在工程和科學中[7，10，14].

與文獻[1]做法一樣，引入無量綱參數(shù)，便得到無量綱的分數(shù)階拋物型兩步熱傳導模型為

(4)

(5)

式(4)中：Kn表示Knudsen數(shù).

此外，需要指出的是，在考慮納米級傳熱時[15]，非跳溫邊界條件忽略了邊界聲子散射的影響，導致在邊界附近的結果令人不滿意.通過觀察發(fā)現(xiàn)，碳納米管的熱流是阻塞的，溫度在管的兩端可以看到跳躍[16-17].為了捕捉納米幾何結構內(nèi)部的邊界聲子散射效應，提出溫度跳躍(Robin′s)邊界條件[1]，即

(6)

為使記號簡單，采用u(x，t)和v(x，t)分別代替模型中的Te(x，t)和Tl(x，t).

1 模型方程

考慮無量綱分數(shù)階兩步模型，即

(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)變?yōu)?/p>

(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)是穩(wěn)定的，并且收斂階為O(τ2-α+h4).

3 數(shù)值算例

考慮帶Robin邊界的Caputo-Hadamard分數(shù)階拋物型兩步模型，即

(28)

(29)

(30)

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

(31)

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

當時間tF=2，空間步長h=1/50時，隨著時間步長τ的減少，取不同時間分數(shù)階α時，例1數(shù)值解的最大誤差和時間上的收斂階的情況，如表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，時間分數(shù)階α=0.5，τ=1/20 000時，隨著空間步長的減少，例1數(shù)值解的最大誤差和空間上的收斂階，如表5所示.由表5可知：空間上的收斂階大概是4階.

表5 α=0.5時的最大誤差和空間上的收斂階(例1)Tab.5 Maximum errors and convergence rates in space for α=0.5 (example 1)

源項

參數(shù)取為B=G=c0=Kn=1.可以算出精確解為

當時間tF=2，h=1/50時，隨著時間步長的減少，取不同的時間分數(shù)階α時，例2數(shù)值解的最大誤差和時間上的收斂階，如表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數(shù)值解的最大誤差和空間上的收斂階，如表10所示.由表10可知：空間上的收斂階也大概是4階.

表10 α=0.2 時的最大誤差和空間上的收斂階(例2)Tab.10 Maximum errors and convergence rates in space for α=0.2 (example 2)

4 結束語

對無量綱的分數(shù)階拋物型兩步熱傳導模型提出了高階有效的數(shù)值算法.方程中采用的是Caputo-Hadamard分數(shù)階導數(shù)，并且考慮的是Robin邊界條件，目前國內(nèi)外文獻還鮮有研究.利用空間四階緊格式和Caputo-Hadamard時間分數(shù)階導數(shù)的L1逼近格式建立了數(shù)值格式.文中的兩個數(shù)值例子驗證了提出的數(shù)值算法的有效性.通過改變Knudsen數(shù)和分數(shù)階導數(shù)的數(shù)值以及邊界條件中的參數(shù)，模擬可以作為分析超短脈沖激光作用下多孔介質(如多孔金屬薄膜)納米尺度熱傳導的工具.后續(xù)可進一步深入探索算法的理論研究.