倪宇暉, 陽 鶯

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

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

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

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

1 預(yù)備知識(shí)

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

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

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

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

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

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

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

(4)

2.1 Poisson方程的離散形式

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

(5)

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

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

(6)

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

(7)

(8)

引理1[8]Poisson方程邊平均有限元的離散形式對(duì)應(yīng)的總剛度矩陣是M-矩陣, 當(dāng)且僅當(dāng)對(duì)所有邊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上視為一個(gè)常向量JT(p), 則由式(12)和式(13)可推出

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

(14)

利用式(7)和式(8)可知, 對(duì)任意的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)對(duì)應(yīng)的總剛度矩陣是M-矩陣的充分必要條件是式(9)成立.

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

2.3 邊平均有限元的實(shí)際計(jì)算

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

(20)

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

這里τE=qi-qj.因此

(21)

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

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

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

(22)

由式(22)可推出

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

(23)

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

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

這里,

2.4 數(shù)值實(shí)驗(yàn)

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

計(jì)算區(qū)域Ω=[-1/2,1/2]3是邊長為1 nm的正方體區(qū)域, 邊界?Ω為正方體的6個(gè)面. 網(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的值非常大, 達(dá)到1025數(shù)量級(jí), 由于計(jì)算機(jī)的精度限制, 該值可能對(duì)結(jié)果有影響, 因此本文做如下系數(shù)規(guī)范化：

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

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

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

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

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

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

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

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