一類Poisson-Nernst-Planck方程的邊平均有限元計算

倪宇暉, 陽 鶯

(桂林電子科技大學 數學與計算科學學院, 廣西高校數據分析與計算重點實驗室, 廣西密碼學與信息安全重點實驗室， 廣西 桂林 541004)

本文考慮一類Poisson-Nernst-Planck(PNP)數學模型, 它是由Nernst-Planck(NP)方程和Poisson方程耦合而成. PNP方程組在半導體[1]和生物分子[2]等領域應用廣泛. 例如, 在半導體領域中, 半導體器件中的電流通常由電子和空穴兩類載流子移動產生, PNP模型用于描述半導體載流子的運動規律, 通過計算方程可得到半導體內的電勢分布和載流子濃度分布.

由于PNP模型是一個強耦合、 強非線性的偏微分方程組, 所以僅在極少數情形下能找到它的解析解. Gummel[3]首次用數值方法代替解析方法討論了一維雙極晶體管的電學性質, 從而開啟了計算機數值模擬半導體器件的時代. 基于半導體器件復雜的空間結構和邊界條件, 善于構造高精度格式的有限元方法[4]和滿足流守恒的有限體積法[5]在實際應用中被廣泛應用.

基于半導體的電學性質, PNP方程常有對流占優的特點, 而利用標準有限元方法對對流占優問題進行求解時易出現偽振蕩. 對于一般的對流占優問題, 目前已有一些特殊的有限元方法, 例如Petrov-Galerkin有限元法[6]、 流線擴散有限元法[7]以及邊平均有限元方法[8]等. 邊平均有限元方法是一類迎風有限元方法, 其能自動選擇迎風方向, 用該方法可使離散所得的剛度矩陣是M-矩陣. 本文針對半導體中的一類PNP方程構造一種邊平均有限元離散形式. 在合適的網格假設下, 離散所得的總剛度矩陣是M-矩陣. 將邊平均有限元的數值結果和標準有限元進行對比, 得出邊平均有限元方法在求解這一類PNP方程所需時間僅為標準有限元的40%, 且計算出的解誤差略低.

1 預備知識

設Ω是3中的有界區域, 滿足Lipchitz連續, ?Ω表示邊界.本文用H1(Ω)表示一階導數屬于L2的函數構成的Soblev空間,C0(Ω)表示區域上的連續函數空間.Γh表示區域Ω分解的單元集合,T∈Γh表示四面體單元.Nh表示節點總數, 用xi(i=1,2,…,Nh)表示區域Ω內的所有節點, 用qi(i=1,2,3,4)表示單元T上的節點,E表示邊,Eij表示以qi為起點、 以qj為終點的邊.

考慮一類在半導體領域中的穩態PNP方程[9]：

方程(1)為Poisson方程, 描述了靜電勢與空穴、 電子濃度分布的關系, 在半導體領域稱為電勢方程. 方程(2),(3)為NP方程, 描述了空穴和電子兩種載流子的運動規律, 又稱為載流子連續性方程. 方程的邊界為Dirichlet邊界條件：

其中ud,pd,nd是與位置坐標相關的已知函數.

2 邊平均有限元的離散形式

首先在線性有限元空間Vh上定義一組基函數φi(i=1,2,…,Nh), 它們均在區域Ω內連續, 在每個單元內是線性函數, 且

φi(xi)=1,φi(xj)=0,j≠i.

(4)

2.1 Poisson方程的離散形式

迭代計算中,p和n會用上一步求解得出的值或為給定初值, 所以方程(1)會變成Poisson方程的形式.下面給出邊平均有限元離散步驟[8].

(5)

針對任意有限元空間的函數Φ定義記號δE為

δEΦ=Φ(qi)-Φ(qj).

(6)

基于式(5)的結果, 可得Poisson方程的邊平均有限元的離散形式:

(7)

(8)

引理1[8]Poisson方程邊平均有限元的離散形式對應的總剛度矩陣是M-矩陣, 當且僅當對所有邊E, 下列等式均成立:

(9)

這里,T?E表示包含邊E的所有單元T.

2.2 NP方程的離散形式

其中

首先定義J(p), 函數ψE和e-ψE在邊E上的調和平均αE如下:

J(p)=p+pu,

(10)

(11)

(12)

(13)

這里,δE由式(6)定義.

如果將J(p)在單元T上視為一個常向量JT(p), 則由式(12)和式(13)可推出

JT(p)·τE≈αEδE(eψEp).

(14)

利用式(7)和式(8)可知, 對任意的vh∈Vh, 均有

(15)

將式(10)和式(14)代入式(15), 可得在單元上滿足

(16)

式(16)右端的ψE與u有關.利用式(16)可得

(17)

(18)

因此可得方程(2)的邊平均有限元離散形式: 尋找ph∈Vh, 滿足

ah(ph,vh)=f(vh),vh∈Vh.

(19)

引理2[8]NP方程的有限元離散形式(19)對應的總剛度矩陣是M-矩陣的充分必要條件是式(9)成立.

由引理1和引理2可見, 當式(9)成立時, PNP方程的總剛度矩陣均為M-矩陣, 使求解過程變得更穩定.

2.3 邊平均有限元的實際計算

當Vh是線性有限元空間時,uh在每個單元上是常數, 此時式(18)可進一步簡化.下面討論ah(ph,φi), 首先將φi代入式(18)并利用基函數的特性(4), 可得

(20)

在實際計算中, 式(20)中αEeψE(xi)和αEeψE(xj)的計算較復雜, 所以本文針對這兩項做進一步處理.如果x是邊E上的一點, 則該點必能表示為x=s(xi-xj)+xj(0≤s≤1), 根據式(11)的定義有

這里τE=qi-qj.因此

(21)

對式(21)沿邊E積分, 并利用微分學基本定理, 可得

這里B(s)是一個Bernoulli函數, 定義為

αEeψE(xj)=B(uh·τE).

(22)

由式(22)可推出

αEeψE(xi)=B(-u·τE).

(23)

將式(22)和式(23)代入式(19), 可得方便計算的邊平均有限元離散形式: 尋找ph∈Vh, 使得對?φi(i=1,2,…,Nh), 滿足

ah(ph,φi)=f(φi),

這里,

2.4 數值實驗

為驗證邊平均有限元方法在PNP方程計算中的有效性, 本文選取一個有真解的數值算例, 并與標準有限元方法進行比較. 實驗環境為個人電腦, CPU為2.10 GHz, 內存為32 GB, 在電腦上運行MATLAB 2018b軟件.

計算區域Ω=[-1/2,1/2]3是邊長為1 nm的正方體區域, 邊界?Ω為正方體的6個面. 網格內所有單元的二面角均小于等于90°, 所以采取的網格剖分滿足總剛度矩陣為M-矩陣的充分必要條件(9), 如圖1所示.

圖1 四面體網格剖分Fig.1 Tetrahedral mesh division

這里,cbulk是為模型真解設定的濃度, 此處取為6.022×1026m-3, 即10-3mol/L.將各物理參數值和單位代入并無量綱化后可見, 真解pr和nr的值非常大, 達到1025數量級, 由于計算機的精度限制, 該值可能對結果有影響, 因此本文做如下系數規范化：

進行系數規范后可得新的方程組：

這里cλ=0.179 07, 邊界條件和真解為

顯然, 求出該方程即可得原問題的解.在不同的自由度下, 本文分別使用邊平均有限元方法和標準有限元方法對該方程進行求解, 結果分別列于表1和表2.

表1 精確解與邊平均有限元解的L2模誤差

表2 精確解與有限元解的L2模誤差

由表1和表2可見, 在自由度相同的情形下, 邊平均有限元方法所用的CPU時間約是有限元方法的40%.

綜上, 本文給出了半導體器件中一類PNP方程的邊平均有限元離散形式, 該離散形式下的總剛度矩陣為M-矩陣. 數值實驗結果表明, 邊平均有限元方法得出的解在誤差上略小于標準有限元法, 且CPU時間更短.

感謝中國科學院數學與系統科學研究院張倩茹在數值實驗方面提供的幫助.