覃柳術(shù)， 陽(yáng) 鶯,2， 唐 鳴

(1.桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西 桂林 541004；2.桂林電子科技大學(xué) 廣西高校數(shù)據(jù)分析與計(jì)算重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541004)

Poisson-Nernst-Planck(PNP)被廣泛應(yīng)用于生物分子[1]、電化學(xué)[2]和半導(dǎo)體[3]等領(lǐng)域中。它是由描述濃度梯度和電場(chǎng)下離子通量的Nernst-Plancks(NP)方程和描述電勢(shì)的Poisson方程耦合而成。

由于PNP方程是一類(lèi)耦合的非線(xiàn)性非對(duì)稱(chēng)的偏微分方程，極少情況下能夠找到解析解，通常使用數(shù)值方法尋找其近似解。常用的數(shù)值方法包括有限差分法、有限體積法、有限元法。有限元法可以很好地處理不規(guī)則區(qū)域，且其解的收斂速度取決于解的正則性。由于NP方程存在電勢(shì)梯度這一低階量，使得有限元解近似的電勢(shì)梯度精度較低，方程存在收斂速度慢甚至發(fā)散的問(wèn)題。為了提高有限元解的導(dǎo)數(shù)的精度，超收斂方法被提出。

L2投影是實(shí)現(xiàn)超收斂的重要方法之一。文獻(xiàn)[5]給出了一致網(wǎng)格上三角線(xiàn)性元和二次元在整體區(qū)域上L2投影的超收斂理論結(jié)果。文獻(xiàn)[6]給出了任意四邊形網(wǎng)格剖分下Poisson方程的L2投影超收斂方法，并給出了超收斂的結(jié)果。文獻(xiàn)[7]提出了一種基于L2投影的梯度恢復(fù)技術(shù)，證明了它具有高階的超收斂性，并將該方法用于斯托克斯方程和彈性方程的求解。文獻(xiàn)[8]將L2投影用于求解Poisson-Boltzmann方程，達(dá)到了超收斂的效果。

鑒于此，針對(duì)一類(lèi)PNP方程，引入L2投影算子，給出了一致網(wǎng)格和非一致網(wǎng)格的誤差分析，同時(shí)通過(guò)數(shù)值實(shí)驗(yàn)進(jìn)行驗(yàn)證。

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

設(shè)Ω是R3中的有界區(qū)域，且滿(mǎn)足Lipchitz連續(xù)，?Ω為其邊界。使用標(biāo)準(zhǔn)的Sobolev空間記號(hào)Wm,p(Ω)和相應(yīng)的范數(shù)和半范數(shù)，當(dāng)p=2時(shí)，記

Ws,2(Ω)=Hs(Ω)，‖·‖s,p,Ω=‖·‖Ws,p(Ω)，

其中v|?Ω=0是在跡的意義下定義的，(·,·)表示L2內(nèi)積。

考慮如下穩(wěn)態(tài)PNP方程：

(1)

其中：pi為第i種離子的濃度；qi為第i種離子的電荷量；φ為靜電勢(shì)；ε為介電常數(shù)。方程(1)的邊界條件考慮齊次的Dirichlet邊界條件p1=p2=φ=0。

方程(1)的弱形式為：

(pi,v)+(qipiφ,v)=(Fi,v),

(2)

(3)

線(xiàn)性有限元空間定義為

Sh={v∈H1(Ω):v|?Ω=0且

v|e∈p1(e),?e∈Γh(Ω)},

其中p1(e)表示線(xiàn)性多項(xiàng)式集合。

假設(shè)方程(1)存在唯一解(φ,pi)，i=1,2，式(2)、(3)的有限元離散格式為：

(4)

(5)

2 超收斂的誤差分析

L2投影算子[5]Qh:L2(Ω)→Sh(Ω)滿(mǎn)足正交關(guān)系：

(w-Qhw,v)=0, ?v∈Sh(Ω)。

(6)

其中：C為與h無(wú)關(guān)的正常數(shù)；θ=min(s,1/2)。

證明由式(2)～(5)和Poincare不等式易證引理1成立。

引理2剖分的網(wǎng)格為一致剖分，有下列不等式成立[9-10]：

‖w-Qhw‖0,Ω≤Ch1+t‖w‖1+t,?w∈H1+t(Ω),

(7)

?v∈Sh(Ω),

(8)

其中η∈(W1,∞(Ω))3。

‖εφ-Qh(εφI)‖0,Ω≤Ch1+t,

(9)

其中0≤t≤1。

證明由式(7)得

‖εφ-Qh(εφI)‖0,Ω≤

‖εφ-Qh(εφ)‖0,Ω+

‖Qh(εφ)-Qh(εφI)‖0,Ω≤

Ch1+t+‖Qh(εφ-εφI)‖0,Ω。

(10)

對(duì)?η∈(W1,∞(Ω))3，根據(jù)式(8)可得

(Qh(εφ-εφI),η)=(φ-φI,εQhη)=

-(-φ-φI,div(εQhη))≤Ch1+t‖η‖0,Ω。

(11)

聯(lián)合式(10)、(11)可推出式(9)成立。

‖εφ-Qh(εφh)‖0,Ω≤

(12)

其中0≤t≤1。

證明由式(9)得

‖εφ-Qh(εφh)‖0,Ω≤

‖εφ-Qh(εφI)‖0,Ω+

‖Qh(εφI)-Qh(εφh)‖0,Ω≤

Ch1+t+‖Qh(εφI-εφh)‖0,Ω。

(13)

根據(jù)引理1，對(duì)任意的φ∈(L2(Ω))3，可得

(Qh(εφI-εφh),φ)=

(14)

由式(13)、(14)可推出式(12)成立。

以上結(jié)論均建立在均勻網(wǎng)格上，以下將在非均勻網(wǎng)格上對(duì)L2投影算子進(jìn)行誤差分析。

在非均勻網(wǎng)格上定義L2投影算子Qh:L2(Ω)→SH(Ω):

(w-Qhw,v)=0,?v∈SH(Ω)。

(15)

定義投影算子：

引理3[10]若

且h?1，則

(16)

其中0≤t≤1。

定理3若

則

‖εφ-Qh(εφh)‖0,Ω≤

(17)

證明由式(15)得

‖εφ-Qh(εφh)‖0,Ω≤

‖εφ-Qh(εφ)‖0,Ω+

‖Qh(εφ-εφh)‖0,Ω≤

CH1+t+‖Qh(εφ-εφh)‖0,Ω,

其中0≤t≤1。對(duì)于任意的ψ∈(L2(Ω))d，有

(Qh(εφ-εφh),ψ)=

由引理3可得

|(Qh(εφ-εφh),ψ)|≤

也即

‖Qh(εφ-εφh)‖0,Ω≤

文獻(xiàn)[4]的數(shù)值實(shí)驗(yàn)表明，

從而有

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

用數(shù)值實(shí)驗(yàn)驗(yàn)證L2投影算子構(gòu)造的超收斂格式的有效性。在Fortan4.0環(huán)境下編譯，所用計(jì)算機(jī)配置為Inter(R)Core(TM) i5-6500 CPU@2.50 GHz，內(nèi)存4 GiB。構(gòu)造如下帶有L2投影算子的Gummel迭代算法。

4)當(dāng)

成立時(shí)，停止迭代；否則，令m=m+1，返回步驟2)。

考慮一類(lèi)帶真解的方程(1)，求解區(qū)域Ω=[0,1]3，q1=1,q2=2，右端函數(shù)Fi和邊界條件依賴(lài)于解析解。解析解定義為：

用線(xiàn)性元對(duì)方程(1)進(jìn)行離散，用算法1進(jìn)行求解，得到有限元解梯度和Qh(φh)逼近真解梯度的L2誤差估計(jì)曲線(xiàn)如圖1所示。圖1中‖φ-φh‖0,Ω表示電勢(shì)真解梯度和有限元解梯度的L2模誤差；‖φ-Qh(φh)‖0,Ω表示電勢(shì)真解和用L2投影算子構(gòu)造Qh(φh)的L2模誤差。從圖1可看出，用L2投影算子構(gòu)造Qh(φh)逼近真解的梯度優(yōu)于有限元解的梯度。

圖1 有限元解梯度和Qh(φh)逼近真解梯度的L2誤差估計(jì)曲線(xiàn)

4 結(jié)束語(yǔ)

通過(guò)引入L2投影算子，對(duì)PNP方程進(jìn)行了誤差分析，對(duì)一類(lèi)帶真解的穩(wěn)態(tài)PNP方程進(jìn)行了數(shù)值實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，相比于有限元解的梯度，利用L2投影算子構(gòu)造的Qh(φh)逼近真解的梯度的精度更高。