線性非局部Drude模型的一種解耦格式的穩(wěn)定性和收斂性分析

徐 倩，夏澤宇，李茂軍

(電子科技大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，成都 611731)

納米光子學(xué)是結(jié)合新興納米科技與光子學(xué)的一門學(xué)科，主要研究在納米尺度上物質(zhì)與光的相互作用[1]，其中，金屬納米材料由于其獨(dú)特的表面等離子體的振蕩現(xiàn)象[2]成為研究的重點(diǎn)對(duì)象之一，在生物傳感、癌癥治療、超材料、新型存儲(chǔ)媒介等多個(gè)領(lǐng)域都有著重要作用[3-4]。

在納米光子學(xué)中，常采用Drude模型來描述電子在金屬中的輸運(yùn)性質(zhì)[5]，近幾十年來，許多學(xué)者致力于研究該模型的空間離散化[6]。有限差分法(Finite Difference Method，F(xiàn)DM)由于其簡(jiǎn)單、易實(shí)現(xiàn)的特點(diǎn)被初步應(yīng)用于求解Drude模型[7]，在此基礎(chǔ)上，Shibayama等[8]利用時(shí)域有限差分(Finite Difference Time Domain，F(xiàn)DTD)方法求解非線性Drude模型，通過對(duì)金屬超材料中光脈沖傳播和二次諧波的模擬，證明了該算法的可行性。此外，Hiremath等[9]采用有限元方法(Finite Element Method，F(xiàn)EM)基于弱形式進(jìn)行了求解。Prokopeva等[10]利用有限體積法(Finite Volume Method，F(xiàn)VM)對(duì)Drude模型進(jìn)行求解，以及對(duì)周期性超材料進(jìn)行二維模擬。

對(duì)于Drude模型的數(shù)值計(jì)算方法，間斷伽遼金(Discontinuous Galerkin，DG)方法是另一種可行的方法。Hesthaven[11]提出一種求解Drude模型的低存儲(chǔ)高精度DG方法。Stannigel等[12]推動(dòng)了Drude模型DG方法的發(fā)展。在此基礎(chǔ)上，文獻(xiàn)[13]提出利用時(shí)域間斷伽遼金(Discontinuous Galerkin Time Domain，DGTD)方法來求解Drude模型。Fezoui等[14]設(shè)計(jì)了一種能夠處理金屬結(jié)構(gòu)局部色散模型的DGTD方法，并給出了穩(wěn)定性和收斂性分析。Li等[15]設(shè)計(jì)了一種用于求解超材料的時(shí)域麥克斯韋方程的節(jié)點(diǎn)DGTD方法。文獻(xiàn)[16]與文獻(xiàn)[17]針對(duì)非局部Drude模型分別提出了采用中心數(shù)值通量和迎風(fēng)數(shù)值通量的DGTD方法。

本文針對(duì)線性非局部Drude模型，基于二階向后差分(BDF)間斷伽遼金(DG)方法提出了一種解耦的能量穩(wěn)定格式，以期證明半離散格式和全離散格式的能量穩(wěn)定性，以及全離散格式的最優(yōu)收斂速率。

1 線性非局部Drude模型

考慮在凸域Ω×[0,T](Ω?3)上的一個(gè)線性非局部Drude模型，表示如下：

?tH+?×E=0,

(1)

?tE-?×H=-J,

(2)

?tQ-?·J=0,

(3)

(4)

其中，未知量H,E,J和Q分別表示磁場(chǎng)、電場(chǎng)、電流密度和電荷密度，常數(shù)β是金屬等離子體頻率，γ是阻尼系數(shù),ωp是金屬等離子體頻率。為了求解模型，給出初始條件和邊界條件[18]：

H|t=0=H0，E|t=0=E0，Q|t=0=Q0,J|t=0=J0，n·J|?Ω=0。

將公式(1)—(4)分別與H,E,Q,J做內(nèi)積，可以得到：

在TE模式下，E=E1e1+E2e2,H=H3e3，其中em,m={1,2,3}為笛卡爾基向量。記H=H3，從而在TE模式下式(1)—(4)可以改寫為如下系統(tǒng)：

?tH+(?xE2-?yE1)=0,

(5)

?tE1-?yH=-J1,

(6)

?tE2+?xH=-J2,

(7)

?tQ-?xJ1-?yJ2=0,

(8)

(9)

(10)

為方便討論，考慮二維矩形區(qū)域[a0,a1]×[b0,b1]以及周期邊界條件：

u(x,y,t)=u(x+Δa,y,t)，u(x,y,t)=u(x,y+Δb,t)，

(11)

其中，Δa=a1-a0，Δb=b1-b0，u∈{H,E1,E2,Q,J1,J2}。接下來，將對(duì)系統(tǒng)(5)—(10)進(jìn)行研究。

2 數(shù)值能量穩(wěn)定性分析

2.1 半離散的DG格式

(12)

且滿足邊界條件：uh(x,y,t)=uh(x+Δa,y,t)，uh(x,y,t)=uh(x,y+Δb,t)，其中Δa=a1-a0，Δb=b1-b0，uh∈{Hh,E1h,E2h,Qh,J1h,J2h}。為證明(12)的能量穩(wěn)定性，先證下面的引理。

證明對(duì)等號(hào)右側(cè)直接進(jìn)行計(jì)算，可以得到

(13)

在所有單元Ii,j，i=1,2,…,Nx，j=1,2,…,Ny上把上面的6個(gè)式子累加，可得

2.2 BDF-DG的全離散格式

(14)

仍然滿足邊界條件：

(15)

注2 為了求解二維系統(tǒng)(5)—(10)，本文設(shè)計(jì)了一種解耦的BDF-DG格式，將(5)—(10)的6個(gè)方程簡(jiǎn)化為2個(gè)小系統(tǒng)。第一部分為(5)—(7)，其未知量為磁場(chǎng)H和電場(chǎng)E，第二部分為(8)—(10)，其未知量為電荷密度Q和電流密度J。即在(5)—(7)中，顯式處理電流密度J，而在(8)—(10)中，顯式處理電場(chǎng)E，因此這兩個(gè)子系統(tǒng)是獨(dú)立求解的。

3 全離散格式的誤差估計(jì)

3.1 投影算子及其估計(jì)

在本節(jié)將給出全離散格式(14)的誤差估計(jì)。在此之前，介紹2個(gè)重要的投影算子。

(16)

(17)

(18)

由于上式中所有項(xiàng)都是線性的，所以利用投影誤差估計(jì)(16)和(17)，可以得到

其中，時(shí)間誤差是利用泰勒展開和柯西不等式得到的。

3.2 全離散格式的最優(yōu)誤差估計(jì)

在對(duì)全離散格式進(jìn)行最優(yōu)誤差估計(jì)之前，需要對(duì)正則性做如下的假設(shè)：當(dāng)k≥1 時(shí)，有

‖u‖l∞([0,T]；Hk+1(Ω))+‖ut‖l∞([0,T]；L2(Ω))+‖utt‖l∞([0,T]；L2(Ω))+‖uttt‖l∞([0,T]；L2(Ω))≤M,

(19)

其中u∈{H,E1,E2,Q,J1,J2},M是一個(gè)正常數(shù)，同時(shí)，用eu=hu-u來表示投影解與數(shù)值解的誤差。

定理3 在滿足式(19)的前提下，由全離散格式(14)得到的數(shù)值解(Hn,En,Qn,Jn)是唯一的，且當(dāng)h

其中h0和Δt0是2個(gè)較小的正常數(shù)，C是與h和Δt的選取無關(guān)的正常數(shù)，n=2，3，…，Nt。

證明對(duì)于全離散格式(14)解的存在唯一性，直接通過方程的線性性，即線性方程組對(duì)應(yīng)的齊次方程只有零解就能夠證明。接下來證明誤差估計(jì)，將投影格式(18)與全離散格式(14)作差，可以得到誤差方程為：

(20)

(21)

由投影估計(jì)(17)與正則化假設(shè)(19) ，可以得到截?cái)嗾`差有如下估計(jì)

4.1 收斂性驗(yàn)證

這并不滿足PDE系統(tǒng)(5)—(10)，因此必須加上源項(xiàng)f=7t6sin2(x)sin2(y)，h=7t6sin2(x)sin2(y),

固定Nx=Ny=128，使空間誤差可以忽略，分別取時(shí)間步長為Δt=1/16,1/32,1/64。此外為了簡(jiǎn)便起見，將原方程中的所有常數(shù)均設(shè)為1，得到解耦的BDF-DG法在h=π/64時(shí)的時(shí)間誤差與收斂階(表1)。由表1的數(shù)值結(jié)果可以看出，當(dāng)Δt變小時(shí)，誤差也會(huì)變小，數(shù)值格式在時(shí)間上是收斂的，并且可以得到數(shù)值格式的時(shí)間精度近似為2階，與誤差的理論分析一致。

表1 解耦的BDF-DG法在h=π/64時(shí)的時(shí)間誤差與收斂階

同時(shí)，取Nx=Ny=150，使空間誤差可以忽略，分別取時(shí)間步長為Δt=1/16,1/32,1/64，得到解耦的BDF-DG法在h=π/75時(shí)的時(shí)間誤差與收斂階(表2)。由表2的數(shù)值結(jié)果可以看出，對(duì)于不同的h，當(dāng)Δt變小時(shí)，誤差也會(huì)變小，在時(shí)間上數(shù)值方法是收斂的，并且同樣可以得到數(shù)值方法的時(shí)間精度近似為2階，與誤差的理論分析一致。

表2 解耦的BDF-DG法在h=π/75時(shí)的時(shí)間誤差與收斂階

固定Nt=1 000，分別取空間步長為h=π/5,π/10,π/20,π/40，得到解耦的BDF-DG法在Δt=1/1 000時(shí)的空間誤差與收斂階(表3)。由表3的數(shù)值結(jié)果可以看出，當(dāng)h變小時(shí)，誤差也會(huì)變小，在空間上數(shù)值方法是收斂的，并且可以得到數(shù)值方法的空間精度近似為2階，與誤差的理論分析一致。

表3 解耦的BDF-DG法在Δt=1/1000時(shí)的空間誤差與收斂階

4.2 能量穩(wěn)定性

為驗(yàn)證能量穩(wěn)定性，選取計(jì)算區(qū)域?yàn)棣浮罷=[0,2π]×[0,2π]×[0,30]，初始條件為：

5 結(jié) 論

為了求解線性非局部Drude模型(5)—(10)，本文采用了一種解耦的二階BDF-DG方法對(duì)其進(jìn)行數(shù)值離散，即用解耦的二階BDF方法進(jìn)行時(shí)間離散，結(jié)合DG方法進(jìn)行空間離散，構(gòu)造相應(yīng)的解耦BDF-DG數(shù)值格式。從理論角度分別證明了DG空間半離散格式和全離散格式的能量穩(wěn)定性。在全離散格式中，將整個(gè)系統(tǒng)解耦成2個(gè)較小的部分，即采用電流密度函數(shù)在電場(chǎng)方程中顯式，且電場(chǎng)函數(shù)在電流密度方程中顯式。這種“解耦”技術(shù)在系統(tǒng)較大的時(shí)候可以提高計(jì)算效率。此外還給出了全離散格式的最優(yōu)收斂速率的分析。最后，通過一些數(shù)值算例進(jìn)一步驗(yàn)證了該理論分析，即解耦的BDF-DG格式具有二階時(shí)間和空間精度，且能量函數(shù)是遞減的。