Nernst Plank 方程有限差分求解

Nernst-Planck方程是描述电解质在电场中运动的方程。其表达式为：

$$\frac{\partial c_i}{\partial t} = -\nabla\cdot J_i + R_i$$

其中，$c_i$表示第$i$种离子的浓度，$J_i$表示第$i$种离子的电迁移流量，$R_i$表示第$i$种离子的反应速率。

为了数值求解Nernst-Planck方程，可以采用有限差分方法。将空间分成网格，时间分成若干个步长，可以将方程离散化。采用中心差分法，可以得到第$i$种离子在网格点$(x,y,z)$处的浓度变化率：

$$\frac{\partial c_i}{\partial t}\Big|_{(x,y,z)} = \frac{c_i^{n+1}(x,y,z)-c_i^n(x,y,z)}{\Delta t}$$

其中，$n$表示时间步数，$\Delta t$表示时间步长。

电迁移流量可以表示为：

$$J_i = -D_i\nabla c_i z_i F\psi$$

其中，$D_i$表示第$i$种离子的扩散系数，$z_i$表示离子的电荷数，$F$表示法拉第常数，$\psi$表示电势。采用中心差分法，可以得到第$i$种离子在网格点$(x,y,z)$处的电迁移流量：

$$J_i\Big|{(x,y,z)} = -\frac{D_i}{\Delta x}\left(\frac{c{i+1/2}(x,y,z)-c_{i-1/2}(x,y,z)}{\Delta x}\right)z_i F\psi(x,y,z)$$

其中，$c_{i+1/2}$表示第$i$种离子在网格点$(x+\Delta x/2,y,z)$处的浓度，$c_{i-1/2}$表示第$i$种离子在网格点$(x-\Delta x/2,y,z)$处的浓度，$\Delta x$表示网格大小。

反应速率可以表示为：

$$R_i = k_i c_ic_j$$

其中，$k_i$表示第$i$种离子的反应速率常数，$c_j$表示其他离子的浓度。采用中心差分法，可以得到第$i$种离子在网格点$(x,y,z)$处的反应速率：

$$R_i\Big|_{(x,y,z)} = k_i c_i(x,y,z)c_j(x,y,z)$$

将上述三个式子代入Nernst-Planck方程中，可以得到：

$$\frac{c_i^{n+1}(x,y,z)-c_i^n(x,y,z)}{\Delta t} = \frac{1}{\Delta x^2}\left[D_i\left(\frac{c_{i+1/2}(x,y,z)-c_{i-1/2}(x,y,z)}{\Delta x}\right)z_i F\psi(x,y,z) - D_i\left(\frac{c_{i+1/2}(x,y,z)-c_{i-1/2}(x,y,z)}{\Delta x}\right)z_i F\psi(x,y,z) + k_i c_i(x,y,z)c_j(x,y,z)\right]$$

这是一个显式差分方程，可以用迭代法求解。初始条件为$c_i^0(x,y,z)$，边界条件为$c_i(x,y,z)$在边界处的值