Poisson-Nernst-Planck 方程, 简称为 PNP 方程, 该方程作为电场和离子浓度分布的耦合方程是描述电扩散过程的一种理想的连续模型^[1,2], PNP 方程已被广泛用于研究离子通道、 纳米孔和燃料电池等领域^[2⇓⇓-5]. PNP 模型通常假设介电系数是常数, 而实际上, 溶液体系的介电系数与周围的电场、浓度和周围的水化层等都有关系, Booth 等^[6]的研究表明, 介电系数是随电场变化的, 这是因为在电解液体系中, 由于强电场的存在, 溶剂分子的偶极子高度定向, 进而影响了溶液体系. Wei 等^[7]提出离子的水化作用会影响离子溶液体系的介电系数. 大量实验和理论研究表明, 介电系数是随局部离子浓度的升高而降低的^[8,9].

对于介电系数依赖于浓度变化的问题, 文献 [10] 提出了非均匀溶液体系下的 Variable Dielectric Poisson-Nernst-Planck (VDPNP) 方程. 该方程引入了依赖于离子浓度的介电系数, 与经典的 PNP 方程相比, VDPNP 模型在 Nernst-Planck (NP) 方程添加了一个非线性项以及在 Poisson 方程的基础上也考虑了依赖于离子浓度的介电系数, 这个模型能更好地描述电解质溶液中复杂的动力学过程, 但这也给 VDPNP 方程的求解带来了许多困难. 据我们所知, 目前很少有人研究过有关 VDPNP 方程的相关工作, 本文针对 VDPNP 方程, 首先基于有限元方法, 对 VDPNP 模型求解. 其次, 针对纳米孔体系进行数值模拟, 考虑不同电荷密度下, PNP 方程和 VDPNP 方程在纳米孔表面浓度、电势及介电系数分布的模拟结果. 同时, 探讨在混合溶液中改变管壁电荷密度对阳离子分布的影响, 并对比分析 PNP 方程和 VDPNP 方程在纳米孔方向上的阳离子分布特征. 模拟结果表明, 改变管壁的电荷密度对阳离子的分布具有显著影响. 随着电荷密度的增加, 表面效应增强, 导致阳离子的浓度分布发生变化, VDPNP 方程能够更好地捕捉这种电荷密度对阳离子分布的影响, 并且在混合溶液中, 传统的 PNP 方程在浓度上无法有效区分钠离子和钾离子, 而 VDPNP 方程可以明显区分钠、钾离子, 且随着表面电荷越大时, 其区分效果越加明显.

本文考虑稳态 VDPNP 方程^[10]

其中 $c_i$ 是第 $i$ 种离子的浓度, $\phi$ 为静电势, $q_i$ 表示第 $i$ 种离子的电荷量, $\epsilon=\epsilon(\bar{c})$ 是一个依赖于离子浓度的非均匀电解质介电系数, 且 $\bar{c}=\sum\limits_{i=1}^k c_i$, $\epsilon^{\prime}(\bar{c})$ 为 $\epsilon$ 对浓度 $\bar{c}=\bar{c}\left(c_1, c_2, \cdots, c_K\right)$ 的导数, $\rho^f(x)$ 是电荷的奇异组合, $\Omega_s$ 表示溶剂区域, $\Omega_m$ 表示溶质区域, $D_i$ 为静电扩散系数, $\beta=\frac{1}{k_{\mathrm{B}} T}$ 为 Boltzmann 能量的倒数, 特征函数 $\lambda=\left\{\begin{array}{l}0, \text{ 在} \Omega_{\mathrm{m}} \text{上} \\ 1, \text { 在} \Omega_{\mathrm{s}} \text{上}\end{array}\right.$, 其边界条件为

其中 $\frac{\partial \phi}{\partial \vec{n}}$ 是边界处的法向导数, $\vec{n}$ 为单位外法向量, $c_i^{\infty}$ 是边界体积浓度.

设 $\Omega\subset\mathbb{R}^d$ ($d=1,2$) 是一个有界的 Lipschitz 区域. 使用标准的 Sobolev 空间记号 $W^{s, p}(\Omega)$. 其中 $\|\cdot\|_{m, \Omega}$ 和 $|\cdot| m, \Omega$ 分别表示相应范数和半范数, $(\cdot, \cdot)_{\Omega}$ 表示标准 $L^2$ 内积, 定义

$$\begin{aligned} &H_s^1(\Omega) :=\left\{\phi \in H^1(\Omega): \phi=\phi_0 \text { 在 } \Gamma_D \text{上}\right\}, \\ &H_{s, 0}^1(\Omega) :=\left\{v \in H^1(\Omega): v=0 \text { 在 } \Gamma_D \text{上}\right\}. \end{aligned}$$

令 $\rho(c)=\rho^f+\lambda \sum\limits_{i=1}^K q_i c_i$, 寻找 $\phi \in H_s^1(\Omega)$, 则

$$\int_{\Omega}-\nabla \cdot(\epsilon(\bar{c}) \nabla \phi) v{\rm d}V=\int_{\Omega} \rho(c) v{\rm d}V, \quad \forall v \in H_{s, 0}^1(\Omega).$$

令 $\Gamma=\Gamma_N+\Gamma_D$ 且 $\left.\forall v\right|_{\Gamma_D}=0$, 根据格林公式及 (2.3) 式, 可以得到 (2.2) 式的变分形式

为了得到 (2.1) 式的变分形式, 定义

$$\begin{aligned} V & :=\left\{c_i\left|c_i \in H^1\left(\Omega_s\right): c_i\right|_{\partial \Omega}=c_i^{\infty}\right\}, \\ V_0 & :=\left\{w\left|w \in H^1\left(\Omega_s\right): w\right|_{\partial \Omega}=0\right\}. \end{aligned}$$

求 $c_i \in V$, 对任意的 $w \in V_0$, 得到方程 (2.1) 的变分形式为

寻找 $\phi \in H_s^1(\Omega), c_i \in V$, 对任意的 $w \in V_0, v \in H_{s, 0}^1(\Omega)$, 采用标准的 $L^2$ 内积, (2.4)-(2.5) 式改写为

区域 $\Omega$ 的一个有限划分 $\mathcal{T}_h$ 是一个单元集合 $\{e\}$, $e$ 是剖分单元, 对每一个单元 $e \in \mathcal{T}_h$, $P(e)$ 表示线性多项式集合. 定义如下线性有限空间

$$V_h=\left\{v \in H^1(\Omega):\left.v\right|_{\partial \Omega}=0, \left.v\right|_e \in P(e), \forall e \in \mathcal{T}_h\right\},$$

则 (2.1)-(2.2) 式的标准有限元离散为: 对 $\forall v_h, w_h \in V_h$, 求 $\left(c_{i, h}, \phi_h\right) \in\left[V_h\right]^3$, $i=1,2$, 使得

本节给出数值例子验证 VDPNP 模型的可求解性. 实验环境为桂林电子科技大学超算平台 (Linux), 内存为 3.58 TB, 采用 For$\operatorname{tran} 90$ 语言编写计算程序, 所有计算结果均来自同一平台.

下面以二维为例, 考虑如下只有正、负两个离子的 VDPNP 模型, 即正离子浓度 $p$, 负离子浓度 $n$, 则 $\bar{c}=p+n$, 考虑线性介电系数模型 $\epsilon(\bar{c})=80-2 \times 10^{-4}\times(p+n) $, 求解区域为 $\Omega=[0,1] \times[0,1]$, 取 $q^1=1, q^2=-1$, 暂时不考虑其他相关物理参数, 则

$$\begin{aligned} &-\nabla \cdot\left(\nabla p+p \nabla\left(\phi-\frac{1}{2} \epsilon^{\prime}(\bar{c})|\nabla \phi|^2\right)\right) =F_1, & \text { 在 } \Omega_s \text{上},\\ &-\nabla \cdot\left(\nabla n-n \nabla\left(\phi-\frac{1}{2} \epsilon^{\prime}(\bar{c})|\nabla \phi|^2\right)\right) =F_2, & \text { 在 } \Omega_s \text{上}, \\ &-\nabla \cdot(\epsilon(\bar{c}) \nabla \phi)-(p-n) =F_3, & \text { 在 } \Omega_s \text{上}, \end{aligned}$$

右端函数依赖于解析解, 其中解析解定义为

VDPNP 方程是一类非线性耦合方程组, 常用的解耦方法是 Gummel 迭代方法, Gummel 迭代方法不仅具有较快的求解速度、良好的收敛性和稳定性, 还可以根据实际需求灵活调整迭代次数和收敛准则, 从而达到所需的数值解精度, 本文 Gummel 迭代如算法 1 所示. 令 $\epsilon(\bar{c})=\epsilon(p, n), f=\frac{1}{2} \epsilon^{\prime}(\bar{c})=f(p, n), g=|\nabla \phi|^2=g(\phi)$. 具体过程见算法 1

在计算中, 选取不同空间步长 $h$ 进行测试, 误差容限选择 $\delta=1.0 \times 10^{-6}$, 非线性迭代的最大次数设定为 500, 利用算法 1 求解方程, 精确解与有限元解的数值误差结果如表 1 所示.

表1 VDPNP 精确解和有限元解的 $L^2$ 模误差

由表 1 可知, 在 $L^2$ 范数下有限元解和精确解的误差收敛阶接近 2 阶, 尽管目前尚无 VDPNP 方程的理论推导来证明其最优收敛阶, 但数值实验结果表明, 其收敛阶达到二阶, 保证了数值计算的稳健性和可靠性.

对于生物分子周围局部高浓度的溶液, 考虑如下的非线性介电模型^[11]

其中 $\epsilon_m$ 为水化层的介电系数, $\epsilon_s$ 为纯水的的介电系数, $\alpha $ 为可调参数, $a_i$ 和 $c_i$ 分别为第 $i$ 种离子的水化层尺寸和浓度.

考虑用纳米管数值模拟 VDPNP 对体系相互作用的影响, 将要模拟的网格示意如图1 所示, 其中外层正方体盒子的边长为 $50\stackrel{\circ}{A}$, 内部纳米管的半径 $r=2 \stackrel{\circ}{A}$, 高度 $h=25 \stackrel{\circ}{A}$, 在内表面给定表面电荷密度 $\sigma$.

我们考虑 1 : 1 $ \mathrm{NaCl}$ 溶液, 在模拟中, 给定表面电荷密度 $\sigma=-0.025, -0.05$, $-0.10\mathrm{C} / \mathrm{m}^2$, 上下膜电压差为 $-0.15 V$, 模拟的 NaCl 溶液的边界体积浓度为 0.1 M, 其中 $\mathrm{Na}^{+}$和 $\mathrm{Cl}^{-}$的水化层直径分别为 $4.79 \stackrel{\circ}{A}$ 和 $6.37 \stackrel{\circ}{A}$^[12,13]. 采用非线性介电系数模型 (4.1), 其中 $\epsilon_m=2$, $\epsilon_s=78$, 可调参数 $\alpha = 2$. 下面首先给出传统 PNP 和 VDPNP 模型纳米孔带电表面附近的阳离子、电势以及介电系数的模拟结果.

由图2 的模拟结果表明, PNP 与 VDPNP 模拟在管壁表面附近的阳离子分布存在差异, 且随着表面电荷越大, 差异越显著; 当表面电荷较低时 (如图2(a) 所示), VDPNP 模拟的阳离子浓度高于 PNP; 另外图2 的 (b), (c) 表明, 当带电表面的电荷增大时, VDPNP 模拟的结果并不再始终大于 PNP, 在离管壁的一定距离处, VDPNP 模拟出的阳离子浓度低于 PNP, 且随着电荷增大, VDPNP 的阳离子浓度在靠近纳米孔表面比 PNP 的要高, 产生这种物理现象的原因可能是由于介电系数的引入使得 VDPNP 方程可能更细致地考虑了离子与带电表面之间的相互作用，如静电吸引和排斥, 这些相互作用会导致阳离子在靠近表面时更紧密的吸附, 浓度增加.

电势的模拟结果如图3 所示, 在低电荷下, VDPNP 模型模拟的电势低于 PNP 模型, 高电荷下 VDPNP 模型的电势高于 PNP. 随着表面电荷密度的增加, 两者的电势差异逐渐增大.

介电系数分布结果如图4 所示, 表面电荷密度越大, 介电系数沿管壁的径向分布变化越显著. 同时, 表面电荷密度越高, 管壁附近的阳离子浓度越大, 介电系数越小, 与介电模型相符.

下面我们进一步研究了纳米孔体系中 $\mathrm{Na}^{+}, \mathrm{K}^{+}$ 在混合溶液中的竞争行为. 溶液中含有 $\mathrm{Na}^{+}, \mathrm{K}^{+}$ 和 $\mathrm{Cl}^{-}$ 三种离子, 模拟的是 0.1 M NaCl 和 0.1 M KCl 的混合溶液, 即 $\mathrm{Na}^{+}, \mathrm{K}^{+}$ 和 $\mathrm{Cl}^{-}$ 三种离子的边界体积浓度分别为 $0.1 \mathrm{M}, 0.1 \mathrm{M}$, $0.2 \mathrm{M}$, 其水化层直径分别为 $4.79 \stackrel{\circ}{A}, 5.51 \stackrel{\circ}{A}$, 和 $6.37 \stackrel{\circ}{A}$. 施加表面电荷与电压差依然与单一离子相同, 即上下膜电压差为 $-0.15 \mathrm{ V}$, $\sigma=-0.025 C / m^2, \sigma=$$-0.05 C / m^2, \sigma=-0.10 C / m^2$, 类似的采用非线性介电系数模型 (4.1), 其中 $\epsilon_m=2$, $\epsilon_s=78$, 可调参数 $\alpha = 2$. 对比了传统 PNP 模型与 VDPNP 模型下纳米孔表面附近阳离子分布及沿孔隙方向阳离子管壁电荷变化情况.

图5 中给出了 PNP 和 VDPNP 方程模拟的纳米孔带电表面的阳离子 $Na^+$ 和 $K^+$ 浓度分布结果, 从图中可以看出, 在靠近纳米孔表面 $K^+$ 比 $Na^+$ 的浓度高, 这可能是由于 $K^+$ 的水合半径比 $Na^+$ 的水合半径更小, 在纳米孔表面受限空间中, 较小的 $K^+$ 更容易穿透水合壳层, 因此更容易接近纳米孔表面, 从而在靠近表面处的浓度更高. 并且由图5 和图6 所示, 对于混合溶液, PNP 方程在浓度上不能区分出钠、钾离子, 而 VDPNP 模型可以明显区分钠、钾离子, 且当表面电荷越大时, 其区分效果越加明显.

本文基于有限元方法对 VDPNP 模型进行了求解, 数值误差结果验证了方程的可解性和数值程序的稳健性. 通过对纳米孔数值模拟, 结果表明, 在纳米孔体系中, 传统的 PNP 方程在浓度上无法有效区分钠离子和钾离子, 而 VDPNP 模型则能够明显区分这两种离子, 并且随着表面电荷密度的增加, 离子区分效果愈加显著. 这表明 VDPNP 模型能够更精确地捕捉复杂电荷相互作用和浓度梯度效应, 适用于纳米尺度下的离子传输和选择性识别问题. 研究结果为理解电荷调控和离子选择性提供了新的理论依据, 并为纳米技术和电化学传输等领域的进一步研究提供了参考.