倪宇暉, 陽 鶯

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

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

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

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

1 預備知識

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

考慮一類在半導體領(lǐng)域中的穩(wěn)態(tài)PNP方程[9]：

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

其中ud,pd,nd是與位置坐標相關(guān)的已知函數(shù).

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

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

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

(4)

2.1 Poisson方程的離散形式

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

(5)

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

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

(6)

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

(7)

(8)

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

(9)

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

2.2 NP方程的離散形式

其中

首先定義J(p), 函數(shù)ψE和e-ψE在邊E上的調(diào)和平均α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有關(guān).利用式(16)可得

(17)

(18)

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

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

(19)

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

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

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

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

(20)

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

這里τE=qi-qj.因此

(21)

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

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

α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 數(shù)值實驗

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

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

圖1 四面體網(wǎng)格剖分Fig.1 Tetrahedral mesh division

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

進行系數(shù)規(guī)范后可得新的方程組：

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

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

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

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

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

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

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