离子通道一维泊松-能斯特-普朗克（Poisson-Nernst-Planck）模型
参考文献：Rectification in synthetic conical nanopores: A one-dimensional Poisson-Nernst-Planck model
理论推导
如上图所示离子通道，体系具有轴对称性，孔半径
\begin{equation} R(z)=R(0)+[R(L)-R(0)]z/L \label{radius} \end{equation}
满足如下方程\begin{equation} \begin{split} \epsilon_0\nabla\cdot[\epsilon(r,z)\nabla\Phi(r,z)]=&-\sum_i\rho_i(r,z)-\rho_f(r,z)\\ =&-10^3eN_A\sum_i\nu_ic_i(r,z)-\delta[r-R(z)]\sigma \end{split} \label{Poisson} \end{equation}
\begin{equation} \begin{split} \frac{\partial c_i(r,z)}{\partial t}=&-\nabla\cdot \vec{j}_i(r,z)\\ =&D_i\nabla\cdot [\nabla c_i(r,z)+e\beta \nu_ic_i(r,z)\nabla \Phi(r,z)]\\ =&0 \end{split} \label{Nernst-Planck} \end{equation}
其中 $\epsilon_0$ 为真空电容率，介电常数 $\epsilon(r,z)=\epsilon_{\mathrm p}\Theta[r-R(z)]+\epsilon_{\mathrm w}\Theta[-r+R(z)]$，这里 $\epsilon_{\mathrm p}$ 和 $\epsilon_{\mathrm w}$ 分别为膜和水的介电常数。$\Phi$为电势。$i$ 标记组分，$\rho_i$ 为电荷密度，单位为$\mathrm{C/m^3}$，$c_i$ 为浓度，单位$\mathrm M$，即 $\mathrm{mol/L}$。$e=1.602\times 10^{-19}\mathrm C$ 为基本电量，$N_A=6.023\times 10^{23}$ 为阿伏加德罗常数，$D_i$为扩散系数，$\beta=1/(k_BT)$。由于孔道很长，可以认为离子和电势在孔道截面上均匀分布，即
\begin{equation} c_i(r,z)\approx c_i(z), \quad \Phi(r,z)\approx \Phi(z) \label{r-uniform} \end{equation}
又散度满足\begin{equation} \begin{split} \nabla\cdot \vec{u}=&\lim_{\Delta_z \to 0}\frac{1}{A(z)\Delta_z}\oint \vec{u}\cdot d\vec{S}\\ =&\lim\limits_{\Delta_z \to 0}\frac{A(z+\Delta_z)u(z+\Delta_z)-A(z)u(z)}{A(z)\Delta_z}=\frac{1}{A(z)}\frac{\partial Au}{\partial z} \end{split} \label{div} \end{equation}
其中$A(z)=\pi R^2(z)$ 为截面面积。将 \eqref{div}式代入\eqref{Poisson}式，得
\begin{equation} \epsilon_0\epsilon_w\frac{d^2\Phi(z)}{dz^2}+\epsilon_0\epsilon_w\frac{d\Phi(z)}{dz}\frac{d\ln A(z)}{dz}=-F\sum_i\nu_i\frac{\overline{c}_i(z)}{A(z)}-\frac{2\sigma}{R(z)} \label{1DPoisson} \end{equation}
其中 $F=eN_A$ 为法拉第常数，$\overline{c}_i(z)=10^3c_i(z)A(z)$ 为一维浓度，单位 $\mathrm{mol/m}$。将 \eqref{div}式代入\eqref{Nernst-Planck}式，得
\begin{equation} \frac{dJ_i}{dz}=-D_i\frac{d}{dz}\left [A(z)\frac{d}{dz}\frac{\overline{c}_i(z)}{A(z)}+e\beta \nu_i\overline{c}_i(r)\frac{d\Phi(z)}{dz} \right ]=0 \label{1DNernst-Planck} \end{equation}
电流 $I_i=F\nu_iJ_i$，则 \eqref{1DNernst-Planck} 可化为\begin{equation} \frac{I_i}{F\nu_iD_i}+\frac{d\overline{c}_i(z)}{dz}-\frac{\overline{c}_i(z)}{A(z)}\frac{dA(z)}{dz}+e\beta \nu_i\overline{c}_i(r)\frac{d\Phi(z)}{dz}=0 \label{1DNernst-PlanckI} \end{equation}
令 $e\beta \Phi(z)=\Phi^*(z)$，$z=z^*L$，$R(z)=R^*(z)R(0)$，$\overline{c}_i=\overline{c}_i^*c_0$，$I_iL/(FD_ic_0)=I_i^*$，分别代入 \eqref{1DPoisson}式和\eqref{1DNernst-PlanckI}式，并考虑到 $A(z)=\pi R^2(z)$，分别得\begin{equation} \frac{1}{\lambda^2}\left(\frac{d^2\Phi^*(z)}{dz^{*2}}+\frac{2}{R^*}\frac{dR^*}{dz^*}\frac{d\Phi^*}{dz^*}\right )+\frac{\varepsilon}{R^*(z)}+\frac{1}{\pi R^{*2}(z)}\sum_i\nu_i\overline{c}_i^*(z)=0 \label{1DPoissonN} \end{equation}
其中 $\frac{1}{\lambda^2}=\frac{\epsilon_0\epsilon_w R^2(0)}{c_0F\beta eL^2}=\frac{R_0^2}{4\pi l_BN_Ac_0L^2}$，$\varepsilon=\frac{2\sigma R(0)}{Fc_0}$\begin{equation} \frac{I_i^*}{\nu_i}+\frac{d\overline{c}_i^*(z)}{dz^*}-\frac{2\overline{c}_i^*(z)}{R^*}\frac{dR^*}{dz}+\nu_i\overline{c}_i^*(r,z)\frac{d\Phi^*(z)}{dz^*}=0 \label{1DNernst-PlanckIN} \end{equation}
边界条件为
\begin{equation} \begin{split} \overline{c}_i^*(0)=&10^3\pi R(0)^2c_{\mathrm{bulk}}/c_0 \\ \overline{c}_i^*(1)=&10^3\pi R(L)^2c_{\mathrm{bulk}}/c_0 \\ \Phi^*(0)=&0\\ \Phi^*(0)=&-e\beta U \end{split} \label{BC} \end{equation}
其中 $c_{\mathrm{bulk}}$ 为本体盐溶液离子浓度，单位为$\mathrm M$，$U$ 为电源电势差，单位为$\mathrm V$。
方程中的参数，$c_0=10^{-14}/6.023 \quad\mathrm{mol/m}=10^{-23}/6.023\quad \mathrm{mol/nm}$，$\lambda^2=4\pi\times 0.71[L/R(0)]^2$，$\varepsilon=2\sigma R(0)$，这里$\sigma$ 单位为$\mathrm{e/nm^2}$，$R(0)$单位为$\mathrm {nm}$，$e\beta=38.91/\mathrm V$。
解方程\eqref{1DPoissonN} 和 \eqref{1DNernst-PlanckIN}即得结果。可以用Matlab bvp4c求解，程序如下。
数值求解
以下程序适用于求解$L$比较大的情形，用continuation得初值，对于比较短的孔道，不必这么麻烦。
|
|