一類含時Poisson-Nernst-Planck方程的時間自適應計算
2024-01-01 00:00:00韋鵬沈瑞剛
吉林大學學報(理學版) 2024年6期
摘要： 基于相鄰時間步的時間自適應有限元方法求解一類含時Poisson-Nernst-Planck（PNP）方程， 以加速求解含有多種離子的二維PNP方程的長時間數值模擬效率. 數值實驗結果表明， 該方法在兩種數值格式下都可以有效加速計算， 對于長時間PNP方程的數值計算有效.
關鍵詞： Poisson-Nernst-Planck方程; 有限元方法; 時間自適應; 數值計算
中圖分類號： O241.82""文獻標志碼： 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）方程是一類用來描述靜電擴散反應過程和帶電粒子傳輸強耦合的偏微分方程組. PNP方程在生物化學領域［1］應用廣泛， 在半導體器件［2-3］、 生物離子通道［4-5］和電化學系統［6-8］等領域也具有重要作用. 由于強耦合非線性因素的制約， 一般情況下很難獲得該方程的解析解， 為求該方程的近似解， 目前已提出了許多數值方法， 如有限體積法［9］、 有限差分法［10］、 有限元法［11］和虛單元法［12］等. 有限差分法因其形式簡單、 編程易實現， 因而在一些簡單模型的求解中得到廣泛應用［13］， 但對于復雜的大分子生物系統， 如細胞膜、 離子通道等復雜問題， 有限元法具有較好的邊界適應性和網格靈
敏度， 在生物大分子模擬研究中得到廣泛關注［14］.
由于含時PNP方程隨著時間演化在某一時間段通常會出現數值解的劇烈變化， 因此， 為觀察并刻畫這一階段的物理擴散過程， 通常需要較小的時間步長和精細的數值計算.
針對上述問題， Fu等［15］提出了一種新的高階時空有限元格式， 并使用經典的PI步長控制算法指導PNP方程時間步變化; Yan等［16］提出了一種半隱式（VSSBDF）的時間離散格式， 并設計基于兩個時間步長下控制算法指導步長變化. 盡管這些方法都取得了很好的效果， 但沒有同時考慮時間離散格式和相鄰時間層之間的影響， 因此這些方法還可以進一步優化. 本文針對一類含時PNP方程， 空間上采取有限元離散格式， 時間上采用向后Euler格式和向
后微分格式， 設計一種基于相鄰時間層的時間自適應算法， 進行有限元計算. 數值實驗結果表明， 該方法對于含時PNP方程長時間演化模擬有效.
1"預備知識
設Ω是2中的有界多邊形區域， Ω是Ω的Lipchitz連續邊界. 使用標準的Sobolev空間記號Ws，p（Ω）. 特別地， 當S=2時， 可記H2（Ω）=W2，p（Ω）， 其中， ‖·‖m ，Ω
和·m，Ω分別表示相應的范數和半范數， （·，·）Ω表示標準L2內積. 假設Th（Ω）為三角形單元τ的并， 且是大小為h形狀正則的網格. 令h=maxτ∈Th hτ， 其中hτ表示單元τ的直徑. 定義線性有限元空間Sh（Ω）={v∈H1（Ω）： vτ∈P1（τ）， τ∈Th（Ω）}，
Sh0（Ω）=Sh（Ω）∩H10（Ω），其中P1（τ）為τ上的線性多項式集合.
考慮如下的含時PNP模型：
pit-·（pi+qipiφ）=Fi，在Ω上，"t∈（0，T），"i=1，2，
-·（φ） =∑2i=1qipi+F3，在Ω上，"t∈（0，T），
pi（0）=pi，0，"φ（0）=φ0，"Fi（0）=F0i，在Ω上，（1）
及齊次Dirichlet邊界條件
φ=0，在Ω上，"t∈（0，T），pi=0，在Ω上，"t∈（0，T），（2）
其中： pi是第i種離子的濃度; φ表示靜電勢; qi 表示第i種離子的電荷量; Fi為反應源項; pi，0，φ0，Fi0表示初始值.
方程組（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"有限元離散格式
假設方程組（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）的全離散格式， 需要對時間區間［0，T］進行劃分， 劃分為N個區間， 分點為0=S0lt;S1lt;…lt;Sn=T. 每個時間層的時間步長為tn，
且Sn=Sn-1+tn-1=S0+t1+t2+…+tn-1. 于是有方程組（5）-（6）的全離散向后微分格式為： 給定（φn-2h，p
1，n-2h，p2，n-2h）∈［Sh0（Ω）］3和（φn-1h，p1，n-1h，p2，n-1h）∈［Sh0（Ω）］3， 求
（φnh，p1，nh，p2，nh）∈［Sh0（Ω）］3滿足：
3pi，nh-4pi，n-1h+pi，n-2h2tn，vh+（pi，nh，vh）+（pi，nhφnh，vh）=（Fi，vh），"vh∈Sh0（Ω），（7）
（φnh，wh）=∑2i=1d（pi，nh，wh）+（F3，wh），"wh∈Sh0（Ω），（8）
其中i=1，2. 將式（7）-（8）變形可得
（3pi，nh，vh）+2tn（pi，nh+qipi，nh
φnh，vh）=2tn（Fi，vh）+（4pi，n-1h-pi，n-2h， vh），"vh∈Sh0（Ω），
（φnh，wh）=∑2i=1qj（pi，nh，wh）+（F3，wh），"wh∈Sh0（Ω）.
求解方程組（7）-（8）得到的解（φnh，p1，nh，p2，nh）稱為全離散向后微分格式下的有限元解.
同理， 對于方程組（5）-（6）的全離散向后Euler格式為： 給定（φn-1h，p1，n-1h，p2，n-1h）∈
［Sh0（Ω）］3， 求（φnh，p1，nh，p2，nh）∈［Sh0（Ω）］3滿足：
pi，nh-pi，n-1htn，vh+（pi，nh，vh）+（qipi，nhφnh，vh）
=（Fi，vh），"vh∈Sh0（Ω），（9）
（φnh，wh）=∑2i=1qi（pi，nh，wh）+（F3，wh），"wh∈Sh0（Ω），（10）
其中i=1，2. 將式（9）-（10）變形可得
（pi，nh，vh）+tn（pi，nh+qipi，nhφnh，vh）=tn（Fi，vh）+（pi，n-1h，vh），"vh∈Sh0（Ω），
（φnh，wh）=∑2i=1qi（pi，nh，wh）+（F3，wh），"wh∈Sh0（Ω）.
求解方程組（9）-（10）得到的解（φnh，p1，nh，p2，nh）稱為全離散向后Euler格式下的有限元解.
文獻［17］給出了方程組（9）-（10）在二維線性有限元解的誤差估計結果.
引理1［17］"設（φn，pi，n）和（φnh，pi，nh）（i=1，2）分別是方程組（5）-（6）和方程組（9）-（10）
的解. 則存在兩個不依賴于網格剖分h的常數τ0和h0， 使得對于n=0，1，…，N， 有
max0≤n≤N‖φn-φnh‖+∑2i=1‖pi，n-pi，nh ‖≤C（t+h2），
其中t≤τ0， h≤h0.
2.2"自適應時間步算法
為解決方程組（9）-（10）的強耦合性和非線性， 目前常用的方法是Gummel迭代法［18］和兩網格方法［17，19］， 由于兩網格方法
對算法實現和調試方面比單網格方法復雜， 因此本文采用相對簡單的帶有Gummel迭代的有限元算法.
算法1"Gummel有限元算法.
輸入： 迭代的初值（?0h，p1，0h，p2，0h）， 時間層時間n及其步長tn， 容錯誤差δ=10-6;
輸出： （?T，p1，T，p2，T）;
步驟1） 初始化非線性迭代： 當n≥0， 且l=0時， （?nh，p1，nh，p2，nh）=（?n+t，0h，pn+t，1，0h，pn+t，2，0h）;
步驟2） 在每個時間層上進行有限元計算： 當（φn+t，lh，p1，n+t，lh，p2，n+t，lh）∈［Sh0（Ω）］3時， 對l≥0且vh，wh∈Sh0（Ω）， 有
（p1，n+tn，l+1h，vh）+tn（p1，n+tn，l+1h+q1p1，n+tn，l+1hφn+tn，lh，vh）=tn（F1，vh）+（p1，nh，vh），
（p2，n+tn，l+1h，vh）+tn（p2，n+tn，l+1h+q1p2，n+tn，l+1hφn+tn，lh，vh）=tn（F2，vh）+（p2，nh，vh），
（φn+tn，l+1h，wh）=（q1p1，n+tn，l+1h+q2p2，n+tn，l+1h，wh）+（F3，wh） ;
步驟3） 中止非線性迭代： 對給定的容錯誤差δ， 當
‖φn+tn，l+1h-φn+tn，lh‖≤δ
時即停止迭代， 此時令
（?n+tn，p1，n+tn，p2，n+tn）=（?n+tn，l+1，p1，n+tn，l+1，p2，n+tn，l+1），
否則， 令l=l+1， 返回步驟2）繼續迭代;
步驟4） 終止循環： 當n+tn=T時停止循環， 否則， 令n=n+tn， 返回步驟1）.
由于一致時間剖分下PNP方程的數值效率較低， 且隨著長時間演化， 誤差累積的現象急劇增強， 給長時間離子通道問題的模擬研究帶來了較大困難. 自適應時間步算法能
提升數值模擬的效率， 特別是在處理含時方程時. 目前， PNP方程的自適應時間步算法有經典的PI控制算法［15］和雙時間步長算法［16］等. 這些算法的主要思想
是根據不同的準則調整時間步長， 通常用戶需指定最大（為準確性）和最小時間步長（為效率性）以及容許誤差. 基于這些數值方法， 本文提出一種基于相鄰時間層的時間自適應算法，
以優化時間步長并加快模擬速度.
算法2"自適應時間步算法.
輸入： （?n-1，p1，n-1，p2，n-1，tn）;
輸出： tn+1;
步驟1） 求有限元解： 在給定的網格上計算電勢和濃度的有限元解（?n，p1，n，p2，n）;
步驟2） 計算單個自適應指示值： ei，n=‖pi，n-pi，n-1‖‖pi，n‖;
步驟3） 計算總的自適應指示值： en=max{e1，n，e2，n}.
if engt;tol "then
更新tn+1=tolen0.5tn;
else
更新tn+1=minmaxtolen0.5tn，tmin，tmax.
end if
算法2給出了自適應時間步長的算法， 該算法以最大時間步長tmax、 最小時間步長 tmin和容限誤差tol為參數來調整自適應水平.
3"數值實驗
下面給出數值實例驗證向后微分和向后Euler數值格式的正確性， 并驗證所提出時間自適應方法的有效性. 實驗環境為桂林電子科技大學超算平臺（Linux）， 內存
為3.58 TB， 采用Fortran90語言編寫計算程序， 所有計算結果均來自同一平臺. 在數值實驗中， 比較了時間自適應方法和一致時間步方法的迭代步數（用IT表示， 單位為次）和迭代時間（用tCPU表示， 單位為s）.
下面以二維為例， 求解的對象為方程組（1）-（2）組成的PNP方程組， 求解區域為Ω=［0，1］×［0，1］， 取q1=1， q2=-1， T=0.5， 右端函數依賴于解析解， 其中解析解定義為
φ=（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范數為uL2=‖u-uh‖22，其中， u表示式（11）中p1，p2，φ的精確解， uh是有限元解. 選取如下收斂階計算公式：
order=log（uL2，i+1/uL2，i）log（hi+1/hi），""i=1，2，….
為驗證向后Euler和向后微分數值格式的有效性及數值結果的正確可靠性， 選取不同空間步長h， T=0.5， 時間步長t1=t2=…=tn=h2， 初值（φ0h，p1，0h，p2，0h）=（0，0，0）進行測試.
表1列出了p1在L2范數下的誤差. 由表1可見， 在L2范數下有限元解和精確解的誤差收斂階接近2階， 表明兩種格式下有限元方法的數值結果收斂階與定理1理論結果一致，
因此可以繼續進行長時間的數值模擬. 于是， 選擇tmax=0.1， tmin=0.000 1和tol=0.01作為自適應時間步算法的參數， 對上述給出兩種數值格式下的PNP方程組（7）-（10）進行計算.
表２列出了自適應算法和恒定步長算法的實驗結果. 表2結果表明， 在固定空間步長h=1/64， 初始時間步長t1=0.001， 并采用兩種不同的格式時， 使用自適應步長算法比恒定步長需要更少的時間和迭代
次數， 從而有效減少了大規模計算的迭代次數和計算成本， 從tCPU和IT的角度出發向后微分數值格式比向后Euler格式迭代次數少， 說明自適應算法對于向后微分格式
效果更好. 此外， 兩種格式都避免了額外的計算開銷而導致的計算時間過長問題.
圖1為兩種格式下濃度p2的數值誤差對比結果. 由圖1可見， 在兩種格式下自適應算法和一致步長算法對于p2的L2誤差
最大的差距是3.1×10-4， 雖然微分格式在L2范數下的誤差波動情況比Euler格式下大， 但微分格式出現了比一致步長算法誤差還小的情況， 因此這是在可接受
范圍內的. 此外， 兩種格式都避免了算法過度調整而導致的不穩定性.
圖2為兩種格式下電勢φ的數值誤差對比結果. 由圖2可見， 對于電勢φ的L2誤差， 兩種格式下自適應算法和一致步長算法下的結果不存在差別， 說明自適應算法對于兩種數值格式下的PNP方程有效.
圖3為兩種格式下的時間層步長對比結果. 由圖3可見， 兩種格式在每個時間層的時間步變換是在最大最小時間步長的限制下， 同時當解趨近于穩定時， 時間步長開始變粗， 趨近最大時間步長， 當解開始震蕩時， 時間步長開始變細. 但微分格式下的時間步長變換比Euler格式下的更快， 表明自適應時間步算法對于向后微分格式更精確.
綜上所述， 本文給出了一種求解含時PNP方程的時間自適應算法， 對一類帶真解的含時PNP方程進行了數值實驗， 數值實驗結果表明了時間自適應方法的有效性. 該方法可進一步推
廣到更復雜的PNP模型中.
參考文獻
［1］"彭波， 李翰林， 盧本卓. 生物分子模擬中的靜電計算 ［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］. 吉林大學學報（理學版）， 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方程的兩網格有限元離散方法［J］. 吉林大學學報（理學版）， 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.）
（責任編輯： 李"琦）
