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

徐 倩，夏澤宇，李茂軍

(電子科技大學 數學科學學院，成都 611731)

納米光子學是結合新興納米科技與光子學的一門學科，主要研究在納米尺度上物質與光的相互作用[1]，其中，金屬納米材料由于其獨特的表面等離子體的振蕩現象[2]成為研究的重點對象之一，在生物傳感、癌癥治療、超材料、新型存儲媒介等多個領域都有著重要作用[3-4]。

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

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

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

1 線性非局部Drude模型

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

?tH+?×E=0,

(1)

?tE-?×H=-J,

(2)

?tQ-?·J=0,

(3)

(4)

其中，未知量H,E,J和Q分別表示磁場、電場、電流密度和電荷密度，常數β是金屬等離子體頻率，γ是阻尼系數,ω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做內積，可以得到：

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

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

(5)

?tE1-?yH=-J1,

(6)

?tE2+?xH=-J2,

(7)

?tQ-?xJ1-?yJ2=0,

(8)

(9)

(10)

為方便討論，考慮二維矩形區域[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}。接下來，將對系統(5)—(10)進行研究。

2 數值能量穩定性分析

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)的能量穩定性，先證下面的引理。

證明對等號右側直接進行計算，可以得到

(13)

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

2.2 BDF-DG的全離散格式

(14)

仍然滿足邊界條件：

(15)

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

3 全離散格式的誤差估計

3.1 投影算子及其估計

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

(16)

(17)

(18)

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

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

3.2 全離散格式的最優誤差估計

在對全離散格式進行最優誤差估計之前，需要對正則性做如下的假設：當k≥1 時，有

‖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是一個正常數，同時，用eu=hu-u來表示投影解與數值解的誤差。

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

其中h0和Δt0是2個較小的正常數，C是與h和Δt的選取無關的正常數，n=2，3，…，Nt。

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

(20)

(21)

由投影估計(17)與正則化假設(19) ，可以得到截斷誤差有如下估計

4.1 收斂性驗證

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

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

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

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

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

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

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

4.2 能量穩定性

為驗證能量穩定性，選取計算區域為Ω×T=[0,2π]×[0,2π]×[0,30]，初始條件為：

5 結 論

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