本文根据对官方案例的初步学习，将自己理解记录在该文档，欢迎大家指出问题，共同讨论学习。
本文共3250字，预计阅读时间：18min。

本案例（如图所示）的主要是围绕求解以下方程组来进行求解的。案例中考虑五种离子的扩散传输问题，即$F{e}^{2+}、N{a}^{+}、C{l}^{-}、H^{+}和O{H}^{-}$。
$$\{\begin{array}{l}{J}_i=-D_i\nabla c_i-z_i{u}_{m,i}F{c}_i\nabla {\varphi }_l\\ i_l=F\sum z_i{N}_i=F\sum z_i{J}_i\\ i_{loc,expr}=i_0(exp(\frac{{\alpha }_{a}F\eta }{RT})-exp(\frac{-{\alpha }_{c}F\eta }{RT})\\ {\eta }_m={\varphi }_s-{\varphi }_l-E_{eq,m}\end{array}\text{\hspace{1em}(1)}$$

点蚀三次电流模型

反应原理

溶解的铁然后在靠近电极表面的电解质中形成氢氧化亚铁，使得靠近电极表面的$O{H}^{-}$被消耗，电极表面附近的pH降低，在本案例中，将假设pH值范围在10到8之间，并且铁氧混合电位相当低，这是氢氧化铁最有可能成为氧化产物的条件¹²。

氧还原反应：$$4H^{+}+4e^{-}+O_2\to 2H_{2}O\text{\hspace{1em}(2)}$$
亚铁氧化反应：$$Fe(s)\to F{e}^{2+}+2e^{-}\text{\hspace{1em}(3)}$$
$Fe(OH{)}_2$的生成反应：$$F{e}^{2+}+2O{H}^{-}\to Fe(OH{)}_2(s)\text{\hspace{1em}(4)}$$
水的自水解：$$H_{2}O\to H^{+}+O{H}^{-}\text{\hspace{1em}(5)}$$
$C{l}^{-}$和$N{a}^{-}$通过在参数中设置$$cC{l}_0=cN{a}_0+c{H}_0+2\ast cF{e}_0-cO{H}_0$$用来保持溶液中的电荷稳定，不参与到动力学反应当中。

反应动力学

Nernst-Planck方程

电解质中稀物质的一般质量平衡由以下每种物质 i 的方程描述³：

$$\frac{\partial C_i}{\partial t}+\nabla \cdot N_i=R_{i,tot}\text{\hspace{1em}(6)}$$
式中$N_i$是物质i的总通量$mol/(m^2\cdot s)$，$R_{i,tot}$表示物质i的随时间的变化率，应该就是反应速率的含义。该通量由Nernst-Planck方程来描述。
$${\mathbf{N}}_i=-D_i\nabla c_i-z_i{u}_{m,i}F{c}_i\nabla {\varphi }_l+c_i\mathbf{u}={\mathbf{J}}_i+c_i\mathbf{u}\text{\hspace{1em}(7)}$$
式中：$u_{m,i}$迁移率($m^2/s$)，${\Phi }_l$表示电解质电位，$J_i$表示相对于对流传输的摩尔通量，即，
$$J_i=-D_i\nabla c_i-z_i{u}_{m,i}F{c}_i\nabla {\varphi }_l\text{\hspace{1em}(8)}$$
电解质中的净电流密度：
$$i_l=F\sum z_i{N}_i=F\sum z_i{J}_i\text{\hspace{1em}(9)}$$
一般是符合电中性理论（多孔电极除外，需要用电流密度和比表面积相乘即电流源 $Q_l$来表示）这个在该模型中，已经用cCl_0的浓度来表示了，即，
$$\nabla \cdot i_l=0\text{\hspace{1em}(10)}$$
回到Nernst-Plack方程中，该方程用来求解距电极表面x处的物质i的流量，在没有搅拌或没有密度梯度的静止溶液中，本模型应该只是考虑了迁移和扩散。给定物质i的初始浓度+电位的边界条件即可求解利用COMSOL求解器求解该方程⁴。根据方程（6）通过得到物质的传输速率。在反应1模块可以看出，本模型重点就是求解R_FeOH2的问题。

Butler–Volmer (BV)方程

电极反应的控制方程就是$$i_{loc,expr}=i_0(exp(\frac{{\alpha }_{a}F\eta }{RT})-exp(\frac{-{\alpha }_{c}F\eta }{RT})\text{\hspace{1em}(11)}$$

式中，${\alpha }_a=0.25$（对于单电子反应，大多数近似为0.5，这个为什么为0.25我也存疑，有兴趣的可以参考这篇回答⁵），${\alpha }_a+{\alpha }_c=1$。

参考交换电流密度$i_0$，在电极反应1方程式图和参数中是用以下公式表达：

$$i_0=i_{0,\mathrm{ref}}(T)\prod _{i:{\nu }_i>0}{\left(\frac{c_i}{c_{i,\text{ ref }}}\right)}^{{\alpha }_c{\nu }_i/n}\prod _{i:{\nu }_i<0}{\left(\frac{c_i}{c_{i,\text{ ref }}}\right)}^{-{\alpha }_{\mathrm{a}}{v}_i/n}\text{\hspace{1em}(12)}$$

模型用该参数i_0_ref=i0_ref_Fe*(tcd.cH/cH_ref)来表示$i_{0,\mathrm{ref}}(T)$，可以看出该值与域中的浓度是相关的，会随时间变化。

i0_ref_Fe  1e0[A/m^2]  1 A/m²  Exchange current density at reference conditions

cH_ref  1e-7[M]  1E-4 mol/m³  Reference concentration of protons for Fe dissolution kinetics 铁溶解动力学的质子参考浓度

tcd.cH  if(tcd.ctemp<0,tcd.cHexpr,tcd.Kw/tcd.cOHexpr)  mol/m³  质子浓度  域 1

tcd.cHexpr  -0.5*tcd.ctemp+sqrt(0.25*tcd.ctemp^2+tcd.Kw)  mol/m³  质子浓度  域 1

tcd.cOHexpr  0.5*tcd.ctemp+sqrt(0.25*tcd.ctemp^2+tcd.Kw)  mol/m³  氢氧化物浓度  域 1

tcd.ctemp  cNa*tcd.z_cNa+cCl*tcd.z_cCl+cFe*tcd.z_cFe  mol/m³  像是除了H+和OH-的其他离子浓度和  域 1，式中z都是代表计量数

边界条件

浓度

在案例中设置了离子的初始浓度，并通过pH来设定了$H^{+}和O{H}^{-}$。初始反应离子浓度参数如下，cCl_0可以根据[1][3][4][5]求得，式中初始浓度都是给定的常数。在给定温度下，水的Water dissociation是固定值，如[7]。因此根据体溶液的pH值可以求得初始的cH_0和cOH_0。再根据电荷守恒可以得出cCl_0。已知溶质溶解的溶度积平衡常数表达式为$$Fe(OH{)}_2(s)⟷F{e}^{2+}+2O{H}^{-}\text{\hspace{1em}(13)}$$，因此可以根据式[3]求的cFe_0。

cNa_0  1[mM]  1 mol/m³   Bulk Na+ concentration [1]

cCl_0  cNa_0+cH_0+2*cFe_0-cOH_0  0.9321 mol/m³  Bulk Cl- concentration [2]

cFe_0  Keq_FeOH2/cOH_0^2  9.6945E-5 mol/m³   Bulk Fe2+ concentration [3]

cH_0  10^(-pKw-pH)[M]  6.8089E-22 mol/m³  Bulk H+ concentration [4]

cOH_0  10^(-pKw+pH)[M]  0.068089 mol/m³  Bulk OH- concentration 体溶液氢氧根 [5]

Keq_FeOH2  (7.66e-6[M])^3  4.4946E-7 mol³/m⁹  Equilibrium constant FeOH2平衡常数 [6]

pKw  14.94-0.04209*(-273.15+T)/1[K]+1.718E-4*(-273.15+T)^2/1[K^2]  14.167   Water dissociation constant 水的离解常数和温度的关系 [7]

pH  10  10   Bulk pH [8]

电位

在模型电极表面模块设置外部电位E_metal，为金属混合电位。应该是根据极化曲线测试求得拟合的自腐蚀电位来进行设置的。Fe的溶解平衡电位应该就是查表⁶所得。电化学反应的速率可以通过反应速率和活化过电位来关联描述，对于反应m，

$${\eta }_m={\varphi }_s-{\varphi }_l-E_{eq,m}\text{\hspace{1em}(14)}$$
其中$E_{eq,m}$可以将Eeq_Fe作为$E_{eq,ref}$带入如下的能斯特方程（8 ）进行求解
$$E_{eq}=E_{eq,ref}(T)-\frac{RT}{nF}ln\prod _i(\frac{c_i}{c_{i,ref}}{)}^{v^i}\text{\hspace{1em}(15)}$$
据文档⁷描述，电解质电位等于边界电解质电位，即${\varphi }_l={\varphi }_{l,bnd}=0V$，已在电解质电位模块设置。电势等于外部电位${\varphi }_s={\varphi }_{s,ext}$，本模型中就是下面的E_metal。因此可以根据方程（14）可以求得电极表面的过电位。模型方程视图中对应的是tcd.eta_er1=tcd.phis_es1-phil-tcd.Eeq_er1，其中tcd.phis_es1是恒定的，等于金属混合电位。后面两项可能是与物质i的浓度相关的。其中电解质电位为正值，如图1所示，但是并没有在方程视图中找到相关的计算公式，可能是根据边界电解质电位来进行随着物质浓度变化的梯度求解。其中平衡电位为正值，如图2所示。

E_metal  -0.31[V]  -0.31 V  Metal mixed potential vs SHE 金属混合电位

Eeq_Fe  -0.44[V]  -0.44 V  Iron dissolution equilibrium potential vs SHE Fe溶解平衡电位

图1 不同时间步下边界的电解质电位

图2 不同时间步下边界的平衡电位

扩散系数

本案例中的扩散系数固定，不随时间和空间位置变化，在参数中设置如下。

D_Cl  2.0e-9[m^2/s]  2E-9 m²/s  Diffusion coefficient, Na+ 扩散系数

D_Fe  0.72e-9[m^2/s]  7.2E-10 m²/s  Diffusion coefficient, Fe2+

D_H  9.3e-9[m^2/s]  9.3E-9 m²/s  Diffusion coefficient, H+

D_Na  1.3e-9[m^2/s]  1.3E-9 m²/s  Diffusion coefficient, Cl-

D_OH  5.3e-9[m^2/s]  5.3E-9 m²/s  Diffusion coefficient, OH-

总结

通过联立方程（8、9、11、14）即可得到方程组（1）。本文是对该模型的总体思路和主体研究进行归纳，还存在模型具体参数来源不清和方程的设置原因不明等问题，后面继续学习研究。对于文中存在的任何问题，欢迎大家在评论区讨论指正。

---

1. 点蚀模型案例 ↩︎

2. COMSOL博客-点蚀讲解 ↩︎

3. COMSOL文档-用户指南 > 电化学界面 > 电流分布接口理论 >Nernst-Planck 方程 ↩︎

4. 李娴娟,许传炬.分数阶Nernst-Planck方程的有限差分/谱元法求解[J].工程数学学报,2010,27(02):207-218. ↩︎

5. Transfer Coefficient究竟是什么？ - 知乎 ↩︎

6. Peter Atkins (1997). Physical Chemistry, 6th edition (W.H. Freeman and Company, New York). ↩︎

7. COMSOL文档-电化学模块 > 用户指南 > 电化学界面 > 电流分布接口理论 > 电化学反应以及一次电流和二次电流分布之间的差异 ↩︎