求解PNP方程的虛單元兩水平算法

關鍵詞：Poisson-Nernst-Planck方程；虛單元方法；兩水平算法；Gummel算法中圖分類號：O241.82 文獻標志碼：A 文章編號：1671-5489（2025）04-1051-08

Virtual Element Two-Level Algorithm for Solving PNP Equations

MAO Wantao，YANG Ying

（Guangxi Applied Mathematics Center （GUET），Guangxi Colleges and Universities Key Laboratory of Data Analysis and Com putation， School of Mathematics amp; Computing Science ，

Guilin University of Electronic Technology， Guilin 541O04，Guangxi Zhuang Autonomous Region， China）

Abstract：We proposed a two-level algorithm based on the virtual element discretization for both steady-state and time-dependent PNP （Poisson-Nernst-Planck） equations. This algorithm first decoupled and linearized the PNP equations using the linear virtual element solution，and then solved them in the quadratic virtual element space. Compared with the commonly used Gummel algorithm for solving the PNP equations， this algorithm could accelerate the solving speed. Numerical experimental results，including both steady-state and time-dependent PNP equations， show that compared with the Gummel algorithm with the linear virtual element，the two-level algorithm has higher accuracy， and consumes less CPU time and is more efficient at comparable accuracy.

Keywords： Poisson-Nernst-Planck equations； virtual element method;two-level algorithm； Gummelalgorithm

考慮一類由 Poisson 方程和 Nernst-Planck方程耦合而成的 PNP方程，該方程通常用于描述離子質量守恒和靜電勢擴散反應過程，在生物分子系統、半導體器件和電化學系統等領域應用廣泛.

由于PNP方程自身有強耦合性和非線性性，通常很難尋求解析解，因此，有限差分法、有限元法和有限體積法等許多數值方法被應用于PNP方程的數值求解中．虛單元法是在有限元法的基礎上發展的一種新型偏微分方程數值離散技術，與傳統有限元法相比，虛單元法有良好的網格適應性，已廣泛應用于偏微分方程數值求解的多個領域．文獻[1]和文獻[2]分別針對穩態和含時PNP方程給出了虛單元方法及相應的誤差分析.

PNP方程是一類非線性耦合方程組，一種常用的解耦和線性化方法是Gummel迭代法[3-5].文獻[6]和文獻[7]分別針對穩態和含時PNP方程提出了有限元兩網格算法．該算法先利用粗網格解對方程進行解耦和線性化，再在細網格上求解方程．相較于傳統的Gummel有限元算法，能加快求解過程．然而，該算法依賴于粗細嵌套網格的構造，而在處理復雜實際問題時，嵌套網格的設計通常很困難，不利于應用．文獻[8]提出了一種有限元兩水平算法，并將其應用到非對稱不定的橢圓問題中.該算法的基本思想是先利用低階有限元解對問題進行對稱化和正定化，然后在高階有限元空間上求解對稱正定問題．由于兩水平算法僅使用一套網格，因此與兩網格算法相比更適用于實際計算.

本文針對PNP問題，采用虛單元法進行離散并結合兩水平算法進行計算．該算法的主要過程是：首先在線性虛單元空間上求解原問題獲得線性元解，利用線性元解對方程進行解耦及線性化，然后在高階虛單元空間上求解解耦后的方程，從而獲得兩水平解．數值實驗結果表明，對于穩態和含時 PNP方程，與常用的線性虛單元的Gummel算法相比，兩水平算法能達到更高的精度．此外，在相當精度下，兩水平算法所需的CPU時間顯著低于Gummel算法，證明了兩水平算法的有效性.

1預備知識

設 是有界多邊形區域．本文采用Sobolev空間 W^k，p（Ω） 中的標準記號及相關的范數和半范數的定義，分別記為 和 |?|_k，p，Ω ．特別地，當 p=2 時，記 ，H₀¹（Ω）={v∣v∈H¹（Ω）}， ， （???，??）_Ω 表示標準的 L² 內積.

一類典型的含時PNP方程為

其初邊值條件為

其中： Ω_T=Ω×（0，T] ， ?Ω_T=?Ω×（0，T] ; ? 為靜電勢， ?¹?p² 表示兩種離子濃度； q¹=1，q²=-1

表示離子的帶電量； c_λ 表示無量綱化后的常數； F₁，F₂，f 為反應源項.

方程（1）所對應的穩態PNP方程[9]為

相應的Dirichlet邊界條件為

i程（1）-（2）的變分形式為：求 pⁱ∈L²（0，T;H₀¹（Ω））（i=1，2） 及 ?∈L²（0，T;H₀¹（Ω）） ，滿足

方程（3）-（4）的變分形式為：求 pⁱ∈H_o¹（Ω）（i=1，2） 及 ，滿足

2 虛單元離散及兩水平算法

下面先給出虛單元空間及相關投影算子的定義，再以問題（1）-（2）為例，給出虛單元離散格式及相應的兩水平算法．問題（3）-（4）的離散格式及兩水平算法可類似給出.

2.1 虛單元空間

假設 T^h={E} 為 的網格剖分，滿足如下假設[10]：對每個單元 E∈T^h ，存在常數 γgt;0 和 Φ_cgt;0 ，使得：

1）單元 E 對半徑大于等于 γh_E 的圓盤上的點都是星形區域；

2）單元 E 任意兩個頂點之間的距離大于等于 ch_E

對 ?E∈T^h ，首先引入空間

定義橢圓投影算子 ： H¹（E）P_k（E） ，對 ?v∈H¹（E） ，滿足

其中

定義 L² 投影算子 Π_k^o：L²（E）P_k（E） ，對 ，滿足

（v-Π_k^?0v，q）_?0，E=0，?q∈P_k（E）.

定義局部虛單元空間為

其中 P_k/P_k-2（E） 表示 P_k（E） 中在 L²（E） 內積下與 P_k-2（E） 正交的多項式．利用上述局部虛單元空間給出整體虛單元空間的定義為

2.2 虛單元離散形式

對 ?u，v，?，u¹，u²∈H_o¹（Ω） ，定義如下形式：

設 S_m^E（?，?） 和 S_a^E（?，?） 為定義在 Q_h^k（E）×Q_h^k（E） 上的對稱雙線性形式[2]，滿足：

α₁（v_h，v_h）_E?S_m^E（v_h，v_h）?α₂（v_h，v_h）_E，?v_h∈Q_h^k（E），II_k⁰v_h=0，

其中 α_i，β_i（i=1，2） 為正常數.

對 ?u_h，v_h，?_h，u_h¹，u_h²∈Q_h^k ，定義

于是方程（5）-（6）對應的虛單元半離散格式為：求 p_hⁱ∈L²（0，T，Q_h^k）（i=1，2） 及 ?_h∈L²（0，T，Q_h^k） ，滿足：

），（11）

下面給出方程（5）-（6）的虛單元全離散格式．先將時間區間等距剖分為 0=t^°1lt;…N=T ，其中時間步長 τ=T/N ， tⁿ=nτ ， n∈Z ．對 ?u_h∈Q_h^k ，記

u_hⁿ=u_h（?，tⁿ），

則方程（5）-（6）對應的向后Euler全離散格式為：尋找 ?_h^i，n（i=1，2） 及 ?_hⁿ∈Q_h^k ，使得

2.3 虛單元兩水平算法

含時 PNP方程的全離散格式（12）是一個耦合的非線性系統，通常采用Gummel迭代法對其進行解耦和線性化．參考文獻[2]，可給出式（12）的基于線性虛單元（即 k=1 ）的Gummel算法.

算法1 Gummel虛單元算法（ k=1AA ）：

輸入：迭代初始值 （p_h^1，0，p_h^2，0，?_h⁰） ，時間步長 τ ，時間 tⁿ ，停機標準（;

輸出： （p_h^1，N，p_h^2，N，?_h^N） ：

步驟1）當 n?1 時，令 p_h^i，n，0=p_h^i，n-1，i=1，2，?_h^n，0=?_h^n-1

步驟2）對 l?1 ，求 （?_h^i，n，l，?_h^n，l）∈[Q_h¹]³ ， i=1，2 ，使得對任意的 v_h?ω_h∈Q_h¹ ，有i=1，2

步驟3）對給定的誤差容限 ? ，當

時即停止迭代，并令 （p_h^1，n，p_h^2，n，?_hⁿ）=（p_h^1，n，l，p_h^2，n，l，?_h^n，l） ；否則，令 ι←ι+1 并返回步驟 2）繼續非線性迭代；

步驟4）當 n+1=N 時停止，否則令 n^?n⁺¹ 并返回步驟1）.

對于穩態PNP方程離散系統的Gummel算法可參考文獻[1]．下面給出本文的主要算法，即虛單元兩水平算法．該算法的主要步驟是：首先在線性虛單元空間上求解原問題，然后利用線性元解對方程進行解耦及線性化，進而在二次虛單元空間上求解解耦后的方程得到兩水平解，方程（12）對應的虛單元兩水平算法如下.

算法2 虛單元兩水平算法.

輸入：迭代初始值 （p_h^1，0，p_h^2，0，?_h⁰） ，時間步長 τ ，時間 tⁿ ：

輸出： （ph*1，N，ph*2，N，?h*N）

步驟1）對 $＼begin{array}{c} n { ＼stackrel { ＼textstyle gt; 1 } { ＼end{array} } }$ ，已知 ?_h^i，n-1∈Q_h¹ ， i=1，2 ，求 ，使得對 ?v_h ，滿足（號 2，

步驟2）已知 p_h^i，n-1∈Q_h² ， i=1，2 ，利用步驟1）中所得的 ?_hⁿ ，求 ，使得對 ?v_h ， ω_h∈Q_h² ，滿足

步驟3）當 n+1=N 時停止，否則令 并返回步驟1）.

對穩態PNP方程（3）的離散系統對應的虛單元兩水平算法也可仿照算法2給出．此外，算法2的誤差估計可參考文獻[7]中的兩網格算法類似給出，其收斂階通常比算法1高一階.

3 數值實驗

下面將虛單元兩水平算法分別應用于二維穩態和含時PNP方程．采用Fortran9O語言編寫計算程序，數值實驗操作環境為個人電腦，配置為CPU-1.80GHz（AMD（R）Ryzen（TM）7-4800U），RAM-16GB．分別給出算法1的線性元解以及算法2的兩水平解的誤差結果，通過比較在達到相當精度時線性元解和兩水平解所需的CPU時間，驗證虛單元兩水平算法的有效性.

數值結果中 H¹ 和 L² 模誤差分別定義為

其中 u 表示 ?，p¹ 和 p² 的精確解， u_h 表示虛單元解．選取如下的收斂階計算公式：

其中 h_i 表示網格尺寸， e_H¹，i 表示在尺寸 h_i 的網格剖分下的 e_H¹，L² 模誤差的收斂階可類似計算．記∣E∣ 為單元 E 的面積，網格尺寸定義為 ：

3.1 穩態PNP模型

考慮如下二維穩態 PNP方程[9]：

相應的邊界條件為

右端 F_i（i=1，2），f 及邊界條件由如下定義的真解給出：

注意到方程（13）是半導體器件中的PNP方程經過無量綱化和規范化后的二維形式[9]．原始求解域為 Ω=[-L/2，L/2]² ，其中 L 表示半導體器件的尺寸．經過無量綱化和規范化處理后，方程（13）的計算域變為 Ω=[-1/2，1/2]² ， c_λ=0.179L² .若 L=1nm ，則 c_λ=0.179 ．因此，對方程（13）的計算，選取 Ω=[-1/2，1/2]² ， c_λ=0.179 ：

在計算中，采用一致三角形網格剖分．表 1～ 表3分別列出了利用算法1和算法2求解方程（13）-（14）的誤差結果．表1和表2的結果表明：算法1的線性元解的 H¹ 模誤差達到1階收斂階；算法2的兩水平解的 H¹ 模誤差達到2階收斂階．說明兩水平算法獲得了更高的精度（這是因為兩水平算法利用二次虛單元求解方程，因而精度更高）．此外，在達到相當精度的情況下，算法2所需的CPU 時間顯著低于算法1．例如，由表1和表2可見：當電勢的線性元解達到精度 ?_h=5.46×10^-2 （20（ h=1/64） 時，所需的 CPU 時間為 15.66s ；當電勢的兩水平解達到更高的精度 ?_h?=3.26×10^-2 中 Δh=1/8） 時，所需的CPU時間僅為 1.20s ，可見，算法2能顯著提高求解效率，證明了虛單元兩水平算法的有效性.

表3展示了穩態PNP方程兩水平解的 L² 模誤差，其收斂階達到三階．考慮到在文獻[1]中已給出算法1的 L² 模誤差的二階收斂結果，且對比分析結果與 H¹ 模誤差相似，因此本文不再列出．圖1為求解穩態PNP方程得到的線性元解和兩水平解總CPU時間和 H¹ 模誤差的關系曲線．由圖1可見，算法2的計算效率遠高于算法1.

3.2 含時PNP模型

考慮問題（1）及初邊值條件：

計算域 的取法與3.1節相同．取時間 T=0，25 ，設置網格比為 τ=h² ，右端 F_i（i=1，2），f 及邊界

條件可由下列定義的真解獲得：

分別用算法1和算法2求解含時PNP方程（1）．表 4～ 表6分別列出了其在最后一個時間層t=0.25 上的線性元解和兩水平解的誤差結果．圖2為含時PNP方程的總CPU時間和 H¹ 模誤差關系曲線.

由表4和表5及圖2可見，算法2比算法1的精度更高，且在相當精度下，兩水平解所需CPU時間遠小于線性元解所需的CPU時間．如當電勢的線性元解 ?_h=1.21×10^-2 （ h=1/64） 時，所需CPU時間為 4446.38s ，而當電勢的兩水平解在達到更高的精度 ?_h*=7.33×10^-3 （ h=1/8） 時，所需CPU時間僅為 4.05s ，表明了兩水平算法求解含時PNP方程（1）的有效性.

表6結果表明，兩水平解的 L² 模誤差收斂階為二階（若將網格比調整為 τ=h³ ， L² 模誤差能達到更高的精度，其收斂階達到三階）．與線性元解的 L² 模誤差相比，算法2所需的CPU時間更少，效率更高．該問題的 L² 模誤差需要更高的網格比才能達到三階，這可能與該問題的特殊性有關（這是半導體器件中的一個BenchMark模型），具體原因需進一步探討.

綜上所述，針對穩態和含時兩種PNP方程，本文提出了一種虛單元兩水平算法．該算法先利用線性虛單元解對方程進行解耦和線性化，然后在二次虛單元空間上求解解耦后的問題．數值實驗結果表明，在達到相當的數值精度下，兩水平算法所需的CPU時間顯著低于常用的線性虛單元的Gummel算法，表明了兩水平算法求解 PNP方程的有效性．由于兩水平算法只需要生成一套網格，因此對實際問題有很強的適用性.

參考文獻

[1]劉楊．二維穩態 Poisson-Nernst-Planck方程的一種虛單元方法［D]．湘潭：湘潭大學，2020.（LIUY.A Virtual Element Method for Two-Dimensional Steady-State Poisson-Nernst-Planck Equations [D]. Xiangtan： Xiangtan University，2020.）

[2]YANG Y，LIU Y，LIU Y，et al. Eror Analysis of Virtual Element Methods for the Time-Dependent PoissonNernst-Planck Equations [J]. Journal of Computational Mathematics，2025，43（3）： 731-770.

[3]BURGER M，PINNAUR. A Globally Convergent Gummel Map for Optimal Dopant Profing [J]. Mathematical Models and Methods in Applied Sciences，2009，19（5）：769-786.

[4]JEROME JW， BROSOWSKI B. Evolution Systems in Semiconductor Device Modeling：A Cyclic Uncoupled Line Analysis for the Gummel Map [J]. Mathematical Methods in the Applied Sciences，1987，9（1）： 455-492.

[5]LU B Z， ZHOU YC. Poisson-Nernst-Planck Equations for Simulating Biomolecular Diffusion-Reaction Proceses ⅡI： Size Efects on Ionic Distributions and Diffusion-Reaction Rates [J]. Biophysical Journal，20l1，100（10）： 2475-2485.

[6]YANG Y，LU B Z， XIE Y. A Decoupling Two-Grid Method for the Steady-State Poison-Nernst-Planck Equations [J]. Journal of Computational Mathematics，2019，37（4）：556-578.

[7]SHEN RG，SHU S，YANG Y，et al. A Decoupling Two-Grid Method for the Time-Dependent Poisson-NernstPlanck Equations [J]. Numerical Algorithms，2020，83（4）：1613-1651.

[8]TANG M， XING X Q，YANG Y，et al. Iterative Two-Level Algorithm for Nonsymmetric or Indefinite Elliptic Problems [J]. Applied MathematicsLetters，2023，140：108594-1-108594-5.

[9]倪宇暉，陽鶯．一類 Poisson-Nernst-Planck 方程的邊平均有限元計算［J]．吉林大學學報（理學版），2022， 60（1）：73-78. （NI Y H，YANG Y. Edge-Averaged Finite Element Calculation for a Class of Poisson-NernstPlanck Equations [J]. Journal of Jilin University （Science Edition），2O22，6O（1）：73-78.）

[10]BEIRAO DA VEIGA L，BREZZI F，MARINI L D，et al. Virtual Element Method for General Second-Order Elliptic Problems on Polygonal Meshes [J].Mathematical Models and Methods in Applied Sciences，2016，26（4）： 729-750.

[11]BEIRAO DA VEIGA L， BREZZIF， MARINI L D，et al. The Hitchiker’s Guide to the Virtual Element Method [J]. Mathematical Models and Methods in Applied Sciences，2014，24（8）：1541-1573. （責任編輯，趙立芹）