摘要： 基于相鄰時(shí)間步的時(shí)間自適應(yīng)有限元方法求解一類含時(shí)Poisson-Nernst-Planck（PNP）方程， 以加速求解含有多種離子的二維PNP方程的長時(shí)間數(shù)值模擬效率. 數(shù)值實(shí)驗(yàn)結(jié)果表明， 該方法在兩種數(shù)值格式下都可以有效加速計(jì)算， 對于長時(shí)間PNP方程的數(shù)值計(jì)算有效.

關(guān)鍵詞： Poisson-Nernst-Planck方程; 有限元方法; 時(shí)間自適應(yīng); 數(shù)值計(jì)算

中圖分類號： O241.82""文獻(xiàn)標(biāo)志碼： A""文章編號： 1671-5489（2024）06-1352-07

Time-Adaptive Computation for a Class of Time-DependentPoisson-Nernst-Planck Equations

WEI Peng， SHEN Ruigang

（Guangxi Colleges and Universities Key Laboratory of Data Analysis and Computation，Guangxi Mathematics "Applied Center （GUET），

School of Mathematics and Computing Science， Guilin University of Electronic Technology，Guilin 541004， Guangxi Zhuang Autonomous Region， China）

Abstract： The time-adaptive finite element method based on adjacent time steps was used to solve a class of time-dependent Poisson-Nernst-Planck （PNP）

equations， which accelerated the efficiency of long-term numerical simulations for solving two-dimensional PNP equations containing multiple ions. Numerical experimental results show

that the method can effectively accelerate computations in both numerical formats and is effective for long-term numerical computation of PNP equations.

Keywords： Poisson-Nernst-Planck equation; finite element method; time-adaptive; numerical computation

Poisson-Nernst-Planck（PNP）方程是一類用來描述靜電擴(kuò)散反應(yīng)過程和帶電粒子傳輸強(qiáng)耦合的偏微分方程組. PNP方程在生物化學(xué)領(lǐng)域^［1］應(yīng)用廣泛， 在半導(dǎo)體器件^［2-3］、 生物離子通道^［4-5］和電化學(xué)系統(tǒng)^［6-8］等領(lǐng)域也具有重要作用. 由于強(qiáng)耦合非線性因素的制約， 一般情況下很難獲得該方程的解析解， 為求該方程的近似解， 目前已提出了許多數(shù)值方法， 如有限體積法^［9］、 有限差分法^［10］、 有限元法^［11］和虛單元法^［12］等. 有限差分法因其形式簡單、 編程易實(shí)現(xiàn)， 因而在一些簡單模型的求解中得到廣泛應(yīng)用^［13］， 但對于復(fù)雜的大分子生物系統(tǒng)， 如細(xì)胞膜、 離子通道等復(fù)雜問題， 有限元法具有較好的邊界適應(yīng)性和網(wǎng)格靈

敏度， 在生物大分子模擬研究中得到廣泛關(guān)注^［14］.

由于含時(shí)PNP方程隨著時(shí)間演化在某一時(shí)間段通常會(huì)出現(xiàn)數(shù)值解的劇烈變化， 因此， 為觀察并刻畫這一階段的物理擴(kuò)散過程， 通常需要較小的時(shí)間步長和精細(xì)的數(shù)值計(jì)算.

針對上述問題， Fu等^［15］提出了一種新的高階時(shí)空有限元格式， 并使用經(jīng)典的PI步長控制算法指導(dǎo)PNP方程時(shí)間步變化; Yan等^［16］提出了一種半隱式（VSSBDF）的時(shí)間離散格式， 并設(shè)計(jì)基于兩個(gè)時(shí)間步長下控制算法指導(dǎo)步長變化. 盡管這些方法都取得了很好的效果， 但沒有同時(shí)考慮時(shí)間離散格式和相鄰時(shí)間層之間的影響， 因此這些方法還可以進(jìn)一步優(yōu)化. 本文針對一類含時(shí)PNP方程， 空間上采取有限元離散格式， 時(shí)間上采用向后Euler格式和向

后微分格式， 設(shè)計(jì)一種基于相鄰時(shí)間層的時(shí)間自適應(yīng)算法， 進(jìn)行有限元計(jì)算. 數(shù)值實(shí)驗(yàn)結(jié)果表明， 該方法對于含時(shí)PNP方程長時(shí)間演化模擬有效.

1"預(yù)備知識

設(shè)Ω是2中的有界多邊形區(qū)域， Ω是Ω的Lipchitz連續(xù)邊界. 使用標(biāo)準(zhǔn)的Sobolev空間記號W^s，p（Ω）. 特別地， 當(dāng)S=2時(shí)， 可記H2（Ω）=W^2，p（Ω）， 其中， ‖·‖m ，Ω

和·m，Ω分別表示相應(yīng)的范數(shù)和半范數(shù)， （·，·）Ω表示標(biāo)準(zhǔn)L2內(nèi)積. 假設(shè)Th（Ω）為三角形單元τ的并， 且是大小為h形狀正則的網(wǎng)格. 令h=maxτ∈Th hτ， 其中hτ表示單元τ的直徑. 定義線性有限元空間Sh（Ω）={v∈H1（Ω）： vτ∈P1（τ）， τ∈Th（Ω）}，

Sh0（Ω）=Sh（Ω）∩H10（Ω），其中P1（τ）為τ上的線性多項(xiàng)式集合.

考慮如下的含時(shí)PNP模型：

pit-·（pi+qipiφ）=Fi，在Ω上，"t∈（0，T），"i=1，2，

-·（φ） =∑2i=1qipi+F3，在Ω上，"t∈（0，T），

pi（0）=p^i，0，"φ（0）=φ0，"Fi（0）=F0i，在Ω上，（1）

及齊次Dirichlet邊界條件

φ=0，在Ω上，"t∈（0，T），pi=0，在Ω上，"t∈（0，T），（2）

其中： pi是第i種離子的濃度; φ表示靜電勢; qi 表示第i種離子的電荷量; Fi為反應(yīng)源項(xiàng); p^i，0，φ0，F(xiàn)i0表示初始值.

方程組（1）-（2）的變分形式為： 求piH10（Ω）（i=1，2）和φH10（Ω）， 使得

pit，v+（pi，v）+（qipiφ，v）=（Fi，v），"v∈H10（Ω），（3）

（φ，w）=∑2i=1qi（pi，w）+（F3，w），"w∈H10（Ω）.（4）

2"PNP方程全離散格式

2.1"有限元離散格式

假設(shè)方程組（3）-（4）存在唯一解（φ0，p1，p2）， 則可給出方程組（3）-（4）的半離散格式： 求（φh，p1h，p2h）∈［Sh0（Ω）］3滿足：

piht，vh+（pih，vh）+（qipihφh，vh）=（Fi，vh），"vh∈Sh0（Ω），"i=1，2，（5）

（φh，wh）=∑2i=1qi（pih，wh）+（F3，wh），"wh∈Sh0（Ω）.（6）

為給出方程組（5）的全離散格式， 需要對時(shí)間區(qū)間［0，T］進(jìn)行劃分， 劃分為N個(gè)區(qū)間， 分點(diǎn)為0=S0lt;S1lt;…lt;Sn=T. 每個(gè)時(shí)間層的時(shí)間步長為tn，

且Sn=S^n-1+t^n-1=S0+t1+t2+…+t^n-1. 于是有方程組（5）-（6）的全離散向后微分格式為： 給定（φ^n-2h，p

^1，n-2h，p^2，n-2h）∈［Sh0（Ω）］3和（φ^n-1h，p^1，n-1h，p^2，n-1h）∈［Sh0（Ω）］3， 求

（φnh，p^1，nh，p^2，nh）∈［Sh0（Ω）］3滿足：

3p^i，nh-4p^i，n-1h+p^i，n-2h2tn，vh+（p^i，nh，vh）+（p^i，nhφnh，vh）=（Fi，vh），"vh∈Sh0（Ω），（7）

（φnh，wh）=∑2i=1d（p^i，nh，wh）+（F3，wh），"wh∈Sh0（Ω），（8）

其中i=1，2. 將式（7）-（8）變形可得

（3p^i，nh，vh）+2tn（p^i，nh+qip^i，nh

φnh，vh）=2tn（Fi，vh）+（4p^i，n-1h-p^i，n-2h， vh），"vh∈Sh0（Ω），

（φnh，wh）=∑2i=1qj（p^i，nh，wh）+（F3，wh），"wh∈Sh0（Ω）.

求解方程組（7）-（8）得到的解（φnh，p^1，nh，p^2，nh）稱為全離散向后微分格式下的有限元解.

同理， 對于方程組（5）-（6）的全離散向后Euler格式為： 給定（φ^n-1h，p^1，n-1h，p^2，n-1h）∈

［Sh0（Ω）］3， 求（φnh，p^1，nh，p^2，nh）∈［Sh0（Ω）］3滿足：

p^i，nh-p^i，n-1htn，vh+（p^i，nh，vh）+（qip^i，nhφnh，vh）

=（Fi，vh），"vh∈Sh0（Ω），（9）

（φnh，wh）=∑2i=1qi（p^i，nh，wh）+（F3，wh），"wh∈Sh0（Ω），（10）

其中i=1，2. 將式（9）-（10）變形可得

（p^i，nh，vh）+tn（p^i，nh+qip^i，nhφnh，vh）=tn（Fi，vh）+（p^i，n-1h，vh），"vh∈Sh0（Ω），

（φnh，wh）=∑2i=1qi（p^i，nh，wh）+（F3，wh），"wh∈Sh0（Ω）.

求解方程組（9）-（10）得到的解（φnh，p^1，nh，p^2，nh）稱為全離散向后Euler格式下的有限元解.

文獻(xiàn)［17］給出了方程組（9）-（10）在二維線性有限元解的誤差估計(jì)結(jié)果.

引理1^［17］"設(shè)（φn，p^i，n）和（φnh，p^i，nh）（i=1，2）分別是方程組（5）-（6）和方程組（9）-（10）

的解. 則存在兩個(gè)不依賴于網(wǎng)格剖分h的常數(shù)τ0和h0， 使得對于n=0，1，…，N， 有

max0≤n≤N‖φn-φnh‖+∑2i=1‖p^i，n-p^i，nh ‖≤C（t+h2），

其中t≤τ0， h≤h0.

2.2"自適應(yīng)時(shí)間步算法

為解決方程組（9）-（10）的強(qiáng)耦合性和非線性， 目前常用的方法是Gummel迭代法^［18］和兩網(wǎng)格方法^{［17，19］}， 由于兩網(wǎng)格方法

對算法實(shí)現(xiàn)和調(diào)試方面比單網(wǎng)格方法復(fù)雜， 因此本文采用相對簡單的帶有Gummel迭代的有限元算法.

算法1"Gummel有限元算法.

輸入： 迭代的初值（?0h，p^1，0h，p^2，0h）， 時(shí)間層時(shí)間n及其步長tn， 容錯(cuò)誤差δ=10^-6;

輸出： （?T，p^1，T，p^2，T）;

步驟1） 初始化非線性迭代： 當(dāng)n≥0， 且l=0時(shí)， （?nh，p^1，nh，p^2，nh）=（?^n+t，0h，p^n+t，1，0h，p^{n+t，2，0}h）;

步驟2） 在每個(gè)時(shí)間層上進(jìn)行有限元計(jì)算： 當(dāng)（φ^n+t，lh，p^{1，n+t，l}h，p^{2，n+t，l}h）∈［Sh0（Ω）］3時(shí)， 對l≥0且vh，wh∈Sh0（Ω）， 有

（p^{1，n+tn，l+1}h，vh）+tn（p^{1，n+tn，l+1}h+q1p^{1，n+tn，l+1}hφ^n+tn，lh，vh）=tn（F1，vh）+（p^1，nh，vh），

（p^{2，n+tn，l+1}h，vh）+tn（p^{2，n+tn，l+1}h+q1p^{2，n+tn，l+1}hφ^n+tn，lh，vh）=tn（F2，vh）+（p^2，nh，vh），

（φ^n+tn，l+1h，wh）=（q1p^{1，n+tn，l+1}h+q2p^{2，n+tn，l+1}h，wh）+（F3，wh） ;

步驟3） 中止非線性迭代： 對給定的容錯(cuò)誤差δ， 當(dāng)

‖φ^n+tn，l+1h-φ^n+tn，lh‖≤δ

時(shí)即停止迭代， 此時(shí)令

（?^n+tn，p^1，n+tn，p^2，n+tn）=（?^{n+tn，l+1}，p^{1，n+tn，l+1}，p^{2，n+tn，l+1}），

否則， 令l=l+1， 返回步驟2）繼續(xù)迭代;

步驟4） 終止循環(huán)： 當(dāng)n+tn=T時(shí)停止循環(huán)， 否則， 令n=n+tn， 返回步驟1）.

由于一致時(shí)間剖分下PNP方程的數(shù)值效率較低， 且隨著長時(shí)間演化， 誤差累積的現(xiàn)象急劇增強(qiáng)， 給長時(shí)間離子通道問題的模擬研究帶來了較大困難. 自適應(yīng)時(shí)間步算法能

提升數(shù)值模擬的效率， 特別是在處理含時(shí)方程時(shí). 目前， PNP方程的自適應(yīng)時(shí)間步算法有經(jīng)典的PI控制算法^［15］和雙時(shí)間步長算法^［16］等. 這些算法的主要思想

是根據(jù)不同的準(zhǔn)則調(diào)整時(shí)間步長， 通常用戶需指定最大（為準(zhǔn)確性）和最小時(shí)間步長（為效率性）以及容許誤差. 基于這些數(shù)值方法， 本文提出一種基于相鄰時(shí)間層的時(shí)間自適應(yīng)算法，

以優(yōu)化時(shí)間步長并加快模擬速度.

算法2"自適應(yīng)時(shí)間步算法.

輸入： （?^n-1，p^1，n-1，p^2，n-1，tn）;

輸出： tⁿ⁺¹;

步驟1） 求有限元解： 在給定的網(wǎng)格上計(jì)算電勢和濃度的有限元解（?n，p^1，n，p^2，n）;

步驟2） 計(jì)算單個(gè)自適應(yīng)指示值： e^i，n=‖p^i，n-p^i，n-1‖‖p^i，n‖;

步驟3） 計(jì)算總的自適應(yīng)指示值： en=max{e^1，n，e^2，n}.

if engt;tol "then

更新tⁿ⁺¹=tolen^0.5tn;

else

更新tⁿ⁺¹=minmaxtolen^0.5tn，tmin，tmax.

end if

算法2給出了自適應(yīng)時(shí)間步長的算法， 該算法以最大時(shí)間步長tmax、 最小時(shí)間步長 tmin和容限誤差tol為參數(shù)來調(diào)整自適應(yīng)水平.

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

下面給出數(shù)值實(shí)例驗(yàn)證向后微分和向后Euler數(shù)值格式的正確性， 并驗(yàn)證所提出時(shí)間自適應(yīng)方法的有效性. 實(shí)驗(yàn)環(huán)境為桂林電子科技大學(xué)超算平臺(tái)（Linux）， 內(nèi)存

為3.58 TB， 采用Fortran90語言編寫計(jì)算程序， 所有計(jì)算結(jié)果均來自同一平臺(tái). 在數(shù)值實(shí)驗(yàn)中， 比較了時(shí)間自適應(yīng)方法和一致時(shí)間步方法的迭代步數(shù)（用IT表示， 單位為次）和迭代時(shí)間（用tCPU表示， 單位為s）.

下面以二維為例， 求解的對象為方程組（1）-（2）組成的PNP方程組， 求解區(qū)域?yàn)棣?［0，1］×［0，1］， 取q1=1， q2=-1， T=0.5， 右端函數(shù)依賴于解析解， 其中解析解定義為

φ=（1-exp{-t}）sin（πx）sin（πy），p1=sin（2πx）sin（2πy）sin（t），

p2=sin（3πx）sin（3πy）sin（2t）.（11）

定義L2范數(shù)為uL2=‖u-uh‖22，其中， u表示式（11）中p1，p2，φ的精確解， uh是有限元解. 選取如下收斂階計(jì)算公式：

order=log（uL2，i+1/uL2，i）log（hi+1/hi），""i=1，2，….

為驗(yàn)證向后Euler和向后微分?jǐn)?shù)值格式的有效性及數(shù)值結(jié)果的正確可靠性， 選取不同空間步長h， T=0.5， 時(shí)間步長t1=t2=…=tn=h2， 初值（φ0h，p^1，0h，p^2，0h）=（0，0，0）進(jìn)行測試.

表1列出了p1在L2范數(shù)下的誤差. 由表1可見， 在L2范數(shù)下有限元解和精確解的誤差收斂階接近2階， 表明兩種格式下有限元方法的數(shù)值結(jié)果收斂階與定理1理論結(jié)果一致，

因此可以繼續(xù)進(jìn)行長時(shí)間的數(shù)值模擬. 于是， 選擇tmax=0.1， tmin=0.000 1和tol=0.01作為自適應(yīng)時(shí)間步算法的參數(shù)， 對上述給出兩種數(shù)值格式下的PNP方程組（7）-（10）進(jìn)行計(jì)算.

表２列出了自適應(yīng)算法和恒定步長算法的實(shí)驗(yàn)結(jié)果. 表2結(jié)果表明， 在固定空間步長h=1/64， 初始時(shí)間步長t1=0.001， 并采用兩種不同的格式時(shí)， 使用自適應(yīng)步長算法比恒定步長需要更少的時(shí)間和迭代

次數(shù)， 從而有效減少了大規(guī)模計(jì)算的迭代次數(shù)和計(jì)算成本， 從tCPU和IT的角度出發(fā)向后微分?jǐn)?shù)值格式比向后Euler格式迭代次數(shù)少， 說明自適應(yīng)算法對于向后微分格式

效果更好. 此外， 兩種格式都避免了額外的計(jì)算開銷而導(dǎo)致的計(jì)算時(shí)間過長問題.

圖1為兩種格式下濃度p2的數(shù)值誤差對比結(jié)果. 由圖1可見， 在兩種格式下自適應(yīng)算法和一致步長算法對于p2的L2誤差

最大的差距是3.1×10^-4， 雖然微分格式在L2范數(shù)下的誤差波動(dòng)情況比Euler格式下大， 但微分格式出現(xiàn)了比一致步長算法誤差還小的情況， 因此這是在可接受

范圍內(nèi)的. 此外， 兩種格式都避免了算法過度調(diào)整而導(dǎo)致的不穩(wěn)定性.

圖2為兩種格式下電勢φ的數(shù)值誤差對比結(jié)果. 由圖2可見， 對于電勢φ的L2誤差， 兩種格式下自適應(yīng)算法和一致步長算法下的結(jié)果不存在差別， 說明自適應(yīng)算法對于兩種數(shù)值格式下的PNP方程有效.

圖3為兩種格式下的時(shí)間層步長對比結(jié)果. 由圖3可見， 兩種格式在每個(gè)時(shí)間層的時(shí)間步變換是在最大最小時(shí)間步長的限制下， 同時(shí)當(dāng)解趨近于穩(wěn)定時(shí)， 時(shí)間步長開始變粗， 趨近最大時(shí)間步長， 當(dāng)解開始震蕩時(shí)， 時(shí)間步長開始變細(xì). 但微分格式下的時(shí)間步長變換比Euler格式下的更快， 表明自適應(yīng)時(shí)間步算法對于向后微分格式更精確.

綜上所述， 本文給出了一種求解含時(shí)PNP方程的時(shí)間自適應(yīng)算法， 對一類帶真解的含時(shí)PNP方程進(jìn)行了數(shù)值實(shí)驗(yàn)， 數(shù)值實(shí)驗(yàn)結(jié)果表明了時(shí)間自適應(yīng)方法的有效性. 該方法可進(jìn)一步推

廣到更復(fù)雜的PNP模型中.

參考文獻(xiàn)

［1］"彭波， 李翰林， 盧本卓. 生物分子模擬中的靜電計(jì)算 ［J］. 計(jì)算物理， 2015，

32（2）： 127-159. （PENG B， LI H L， LU B Z. Electrostatic Calculation in Biomolecular Modeling ［J］. Chinese Journal Computational Physics， 2015， 32（2）： 127-159.）

［2］"VAN ROOSBROECK W. Theory of the Flow of Electrons an

d Holes in Germanium and Other Semiconductors ［J］. The Bell System Technical Journal， 1950， 29（4）： 560-607.

［3］"JEROME J W. Analysis of Charge Transport： A Mathematical Study of Semiconductor Devices［M］. New York： Springer-Verlag， 1996： 1-167.

［4］"SINGER A， NORBURY J. A Poisson-Nernst-Planck Model

for Biological Ion Channels—An Asymptotic Analysis in a Three-Dimensional Narrow Funnel ［J］. SIAM Journal on Applied Mathematics， 2009， 70（3）： 949-968.

［5］"EISENBERG B. Ionic Channels in Biological Membranes-Electrostatic

Analysis of a Natural Nanotube ［J］. Contemporary Physics， 1998， 39（6）： 447-466.

［6］"VAN SOESTBERGEN M， BIESHEUVEL P M， BAZANT M Z. Diffuse-Charge Effe

cts on the Transient Response of Electrochemical Cells ［J］. Physical Review E， 2010， 81（2）： 021503-1-021503-13.

［7］"RICHARDSON G， KING J R. Time-Dependent Modelling and Asymptoti

c Analysis of Electrochemical Cells ［J］. Journal of Engineering Mathematics， 2007， 59（3）： 239-275.

［8］"MARCICKI J， CONLISK A T， RIZZONI G. Comparison of Limiting Descr

iptions of the Electrical Double Layer Using a Simplied Lithium-Ion Battery Model ［J］. ECS Transactions， 2012， 41（14）： 9-21.

［9］"BANGK R E， "COUGHRAN W M， Jr， COWSAR L C. The Finite Volume Scharfetter-Gummel Method for Steady Convection Diffusion Equations ［

J］. Computing and Visualization in Science， 1998， 1（3）： 123-136.

［10］"CRDENAS A E， COALSON R D， KURNIKOVA M G. Three-Dimensional Po

isson-Nernst-Planck Theory Studies： Influence of Membrane Electrostatics on Gramicidin a Channel Conductance ［J］. Biophysical Journal， 2000， 79（1）： 80-93.

［11］"LU B Z， HOLST M J， MCCAMMON J A， et al. Poisson-Ne

rnst-Planck Equations for Simulating Biomolecular Diffusion-Reaction Processes Ⅰ： Finite Element So

lutions ［J］. Journal of Computational Physics， 2010， 229（19）： 6979-6994.

［12］"丁聰， 劉楊， 陽鶯， 等. 一種三維Poisson-Nernst-Planck方程的虛單元計(jì)算［J］. 吉林大學(xué)學(xué)報(bào)（理學(xué)版）， 2024， 62（2）： 293-301.

（DING C， LIU Y， YANG Y， et al. Virtual Element Computation for a Three-Dimensional Poisson-Nernst-Planck Equations［J］. Journal of Jilin University （Science Edition）， 2024， 62（2）： 293-301.）

［13］"JEROME J W， KERKHOVEN T. A Finite Element Approxima

tion Theory for the Drift Diffusion Semi-conductor Model ［J］. SIAM Journal on Numerical Analysis， 1991， 28（2）： 403-422.

［14］"SONG Y H， ZHANG Y J， SHEN T Y， et al. Finite Elemen

t Solution of the Steady-State Smoluchowski Equation for Rate Constant Calculations ［J］. Biophysical Journal， 2004， 86（4）： 2017-2029.

［15］"FU G S， XU Z L. High-Order Space-Time Finite Element Methods for the Poisson-Nernst-Planck Equations： Positivity and Uncondi

tional Energy Stability ［J］. Computer Methods in Applied Mechanics and Engineering， 2022， 395： 115031-1-115031-15.

［16］"YAN D， PUGH M C， DAWSON F P. Adaptive Time-Stepping Schemes for the Solution

of the Poisson-Nernst-Planck Equations ［J］. Applied Numerical Mathematics， 2021， 163： 254-269.

［17］"SHEN R G， SHU S， YANG Y， et al. A Decoupling Two-Grid Method

for the Time-Dependent Poisson-Nernst-Planck Equations ［J］. Numerical Algorithms， 2020， 83（4）： 1613-1651.

［18］"LU B Z， ZHOU Y C. Poisson-Nernst-Planck Equations

for Simulating Biomolecular Diffusion-Reaction Processes Ⅱ： Size Effects on I

onic Distributions and Diffusion-Reaction Rates ［J］. Biophysical Journal， 2011， 100（10）： 2475-2485.

［19］"唐鳴， 陽鶯， 李雪芳. 一類Poisson-Nernst-Planck方程的兩網(wǎng)格有限元離散方法［J］. 吉林大學(xué)學(xué)報(bào)（理學(xué)版）， 2019， 57（3）： 523-529.

（TANG M， YANG Y， LI X F. Two-Grid Finite Element Discretization Methods for a Class of Poisson-Nernst-Planck Equations［J］. Journal of Jilin University （Science Edition）， 2019， 57（3）： 523-529.）

（責(zé)任編輯： 李"琦）