电化学的相场建模学习（一）

电化学的相场建模

Phase field modeling of electrochemistry. I. Equilibrium

建立电化学系统的相场模型
通过简单的假设：质量和体积约束、泊松方程、本体溶液中理想溶液热力学及界面中竞争能量
模型捕获：平衡状态下双电层相关的电荷分离、电解质中静电势衰减。
研究结果：表面自由能、表面电荷和微分电容

1. 模型描述

1.1 组分、相和摩尔体积

电化学系统描述：

系统由两相组成：金属导电电极为$\ast \ast \alpha \ast \ast$相，离子导电电解质为$\ast \ast \beta \ast \ast$相。

包含组分：阳离子、阴离子、电子、（溶解分子）。
组分在整个界面区域平滑变化。
金属导电电极$\ast \ast \alpha \ast \ast$相包含两个组分：组分1电子$e^{-}$、组分2阳离子$M^{m+}$
离子导电电解质$\ast \ast \beta \ast \ast$相包含三个组分：组分2阳离子$M^{m+}$（主要阳离子）、组分3阳离子$N^{n+}$（次要阳离子，或者$n=0$时为溶剂分子）、组分4阴离子$A^{a-}$

主要发生电荷转移反应：

$$\begin{array}{r}{M}^{m+}(\beta )+m{e}^{-}(\alpha )\rightleftharpoons M(\alpha )\end{array}$$

假设1：所有粒子被视为具有相同偏摩尔体积$\overline{V_s}$，电子视为具有 $0$ 偏摩尔体积的间隙物质（实际各物种的偏摩尔体积至少取决于物质和相，这种变化会导致变形和流动）

定义每个组分$j$的摩尔分数$X_j$，则有：

$$\begin{array}{r}\sum _{j=1}^n{X}_j=1\end{array}$$

其中，$n=4$表示四种物种。系统摩尔体积的变化是因为组分的摩尔分数从一相通过界面进入另一相而变化。在每一点摩尔体积由下式描述：

$$\begin{array}{r}{V}_m=\sum _{j=1}^n\overline{V_j}{X}_j=\overline{V_s}\sum _{j=2}^n{X}_j\end{array}$$

其中$\overline{V_j}$表示各组分的偏摩尔体积，${\sum }_{j=1}^n$是所有组分求和，${\sum }_{j=2}^n$是除电子以为所有组分求和。每大单位体积摩尔浓度定义为$C_j\equiv X_j/V_m$，因此有：

$$\begin{array}{r}\sum _{j=1}^n\overline{V_j}{C}_j=\overline{V_s}\sum _{j=2}^n{C}_j=1\end{array}$$

这和格子气模型相似，所有物种的浓度之和等于格子位点浓度；对于任意物种$j=n$且$\overline{V_n}=\overline{V_s}\ne 0$可以由其他组分计算得到$C_n=1/\overline{V_s}-{\sum }_{j=2}^{n-1}{C}_j$

1.2 系统平衡

对整个等温带电系统用Helmholtz Free Energy描述：

$$\begin{array}{r}F(\xi ,C_1,...,C_n,\varphi )={\int }_v(f_V(\xi ,C_1,...,C_n)+\frac{{\kappa }_{\xi }}{2}|\nabla \xi {|}^2+\frac{1}{2}\rho \varphi )dV\end{array}$$

积分为整个体积域$V$，其中$f_V$是单位体积的Helmholtz Free Energy，$\xi$是相场变量，$\varphi$是静电势，${\kappa }_{\xi }$是相场的梯度能量系数，$\rho$是电荷密度由下式表示：

$$\begin{array}{r}\rho \equiv F\sum _{j=1}^n{z}_j{C}_j\end{array}$$

其中$z_j$是化合价，$F$是法拉第常数。

在式（5）中

第一项：没有梯度且没有电荷相互作用的系统能量密度；

第二项：没有静电效应的梯度能量密度；

第三项：静电能量密度；

整个系统存在三个约束条件：

约束1：每个物种$j$的总摩尔数（$n_j$）必须在体积$V$中守恒；

$$\begin{array}{r}{n}_j={\int }_v{C}_{j}dV=V\overline{C_j},j=1,...,n\end{array}$$

其中$\overline{C_j}$是物种$j$的平均摩尔浓度（在初始状态就已经确定了）。

约束2：Poisson方程必须在系统中每点都满足；

$$\begin{array}{r}\nabla \cdot [ϵ(\xi )\nabla \varphi ]+\rho =0\end{array}$$

其中$ϵ(\xi )$认为是与相相关的介电常数。

通过将Poisson方程代入式（5）中第三项可以得到，$\frac{1}{2}\rho \varphi =\frac{1}{2}D\cdot E=\frac{ϵ(\xi )}{2}|\nabla \varphi {|}^2$，由高斯定律推导，具体推导参考链接。看一看到这与梯度能量密度具有相同的形式。

约束3：满足方程（4）系统质量/体积守恒；

根据约束计算拉格朗日量：

L = F − ∫ V λ V ( x ) ( ∑ j = 1 n V j ˉ C j − 1 ) d V − ∑ j = 1 n λ j ∫ V ( C j − C j ˉ ) d V − ∫ V λ ϕ ( x ) { ∇ ⋅ [ ϵ ( ξ ) ∇ ϕ + ρ } d V \begin{align} \mathcal{L}=F-\int_V\lambda_V({\bf x})(\sum_{j=1}^n\bar{V_j}C_j-1)dV-\sum_{j=1}^n\lambda_j\int_V(C_j-\bar{C_j})dV-\int_V\lambda_\phi({\bf x})\{\nabla \cdot[\epsilon(\xi)\nabla\phi+\rho\}dV \end{align} L=F−∫V​λV​(x)(j=1∑n​Vj​ˉ​Cj​−1)dV−j=1∑n​λj​∫V​(Cj​−Cj​ˉ​)dV−∫V​λϕ​(x){∇⋅[ϵ(ξ)∇ϕ+ρ}dV​​

其中，${\lambda }_V(\mathbf{x}),{\lambda }_{\mathbf{j}},{\lambda }_{\varphi }(\mathbf{x})$分别为约束条件（4）（7）（8）引入的拉格拉日乘子。其中约束条件（4）（8）要求在整个系统中各点满足，因此${\lambda }_V(\mathbf{x}),{\lambda }_{\varphi }(\mathbf{x})$为场变量。

在系统平衡时，$\ast \ast \mathcal{L}$的每个变量变化都必须独立为零**，即：

$$\begin{array}{r}\frac{\delta \mathcal{L}}{\delta C_j}=0=\frac{\partial f_V}{\partial C_j}+\frac{1}{2}F{z}_j\varphi -{\lambda }_j-F{z}_j{\lambda }_{\varphi }-{\lambda }_V\overline{V_j}\end{array}\text{\hspace{1em}(10a)}$$

$$\begin{array}{r}j=1,...,n\end{array}$$

$$\begin{array}{r}\frac{\delta \mathcal{L}}{\delta \xi }=0=\frac{\partial f_V}{\partial \xi }-{\kappa }_{\xi }{\nabla }^2\xi +{ϵ}^{\prime }(\xi )\nabla \lambda \varphi \nabla \varphi \end{array}\text{\hspace{1em}(10b)}$$

$$\begin{array}{r}\frac{\delta \mathcal{L}}{\delta \varphi }=0=\frac{1}{2}\rho -\nabla \cdot [ϵ(\xi )\nabla {\lambda }_{\varphi }]\end{array}\text{\hspace{1em}(10c)}$$

使用Poisson方程消除等式（10c）中的$\rho$可以得到${\lambda }_{\varphi }=-\varphi /2$。

假设2：在系统域$V$的边界处保持自然边界条件$0=\partial \xi /\partial \mathbf{n}=\partial \varphi /\partial \mathbf{n}=\partial {\lambda }_{\varphi }/\partial \mathbf{n}$。

代入（10a）中得到

$$\begin{array}{r}{\lambda }_j=0=\frac{\partial f_V}{\partial C_j}+F{z}_j\varphi -{\lambda }_V(\mathbf{x})\overline{V_j},j=1,...,n\end{array}$$

由于假设1：$\overline{V_{e^{-}}}=0，\overline{V_j}=\overline{V_s}$。可以消除${\lambda }_V(\mathbf{x})$，因此可得到整个平衡控制方程**如下：

$$\begin{array}{r}{\lambda }_{jn}={\lambda }_j-\frac{\overline{V_j}}{\overline{V_s}}{\lambda }_n=[\frac{\partial f_V}{\partial C_j}+F{z}_j\varphi ]-\frac{\overline{V_j}}{\overline{V_s}}[\frac{\partial f_V}{\partial C_n}+F{z}_n\varphi ]=const\end{array}\text{\hspace{1em}(12a)}$$

$$\begin{array}{r}j=1,...,n-1\end{array}$$

$$\begin{array}{r}0=\frac{\partial f_V}{\partial \xi }-{\kappa }_{\xi }{\nabla }^2\xi -\frac{{ϵ}^{\prime }(\xi )}{2}(\nabla \varphi {)}^2\end{array}\text{\hspace{1em}(12b)}$$

$$\begin{array}{r}0=\nabla \cdot [ϵ(\xi )\nabla \varphi ]+\rho \end{array}\text{\hspace{1em}(12c)}$$

此外经典化学势表述如下：

$$\begin{array}{r}{\mu }_j=\frac{\partial f_V}{\partial C_j},j=1,...,n\end{array}$$

经典电化学势表述如下：

$$\begin{array}{r}\overline{{\mu }_j}=\frac{\partial f_V}{\partial C_j}+F{z}_j\varphi ,j=1,...,n\end{array}$$

代入方程（12a）可以得到

$$\begin{array}{r}\overline{{\mu }_j}-\frac{\overline{V_j}}{\overline{V_s}}\overline{{\mu }_n}=const\end{array}$$

Appendix A: First Integral

根据附录推导，可以得到

$$\begin{array}{r}\overline{{\mu }_j}=\overline{{\mu }_j^{\infty }}+\overline{V_j}(\frac{\kappa \xi }{2}(\nabla \xi {)}^2-\frac{ϵ(\xi )}{2}(\nabla \varphi {)}^2),j=1,...,n\end{array}$$

其中，$\overline{{\mu }_j^{\infty }}$是远离界面的电化学势能（平衡系统中）。

可以看到，在远离界面处，其中$\nabla \xi =\nabla \varphi =0$因此每个组分$j$的电化学势在每个相内是相同的，这与吉布斯热力学一直：

$$\begin{array}{r}\overline{{\mu }_j^{\alpha }}=\overline{{\mu }_j^{\beta }},j=1,...,n\end{array}$$

此外，由于电子$e^{-}$作为间隙物种，$\overline{V_{e^{-}}}=0$，因此电子的电化学势在整个系统中是均匀的。而其他物种电化学势通过界面发生变化。

因为系统平衡，在远离界面$\nabla \varphi =0$，在没有外部电荷的情况下，要求电荷在体相中为$0$，因此：

$$\begin{array}{r}\sum _{j=1}^n{z}_j{X}_j^{\alpha }=\sum _{j=1}^n{z}_j{X}_j^{\beta }=0\end{array}$$

2. 界面特性

根据电毛细管的理论，可以知道表面自由能$\gamma$、电极表面剩余电荷${\sigma }^{\alpha }$、微分电容$C_d$符合以下关系：

$$\begin{array}{r}{\sigma }^{\alpha }=-(\frac{\partial \gamma }{\partial \Delta \varphi }{)}_{{\mu }_j}\\ C_d\equiv (\frac{\partial {\sigma }^{\alpha }}{\partial \Delta \varphi }{)}_{{\mu }_j}=-(\frac{{\partial }^2\gamma }{\partial \Delta {\varphi }^2}{)}_{{\mu }_j}\end{array}$$

其中$\Delta \varphi ={\varphi }^{\alpha }-{\varphi }^{\beta }$是界面处施加的电势差。远离界面的电解质中的电势为${\varphi }^{\beta }$，电极中的电势为${\varphi }^{\alpha }$（理想导体各处电势均匀）。

Appendix B: Surface Free Energy

根据附录B推导，可以得到表面自由能：

$$\begin{array}{r}\gamma ={\int }_{-L/2}^{L/2}[{\kappa }_{\xi }(\nabla \xi {)}^2-ϵ(\xi )(\nabla \varphi {)}^2]dx\end{array}$$

其中，被积分项第一项表示：非带电界面的贡献；第二项表示：静电的贡献。因此电场的存在总是会降低表面自由能的非带电值。

Appendix C: Surface Charge and Capacitance

根据附录C推导，得到关系式：

$$\begin{array}{r}{\sigma }^{\alpha }=-(\frac{\partial \gamma }{\partial \Delta {\varphi }^0}{)}_{{\mu }_j}\end{array}$$

并定义电极侧的表面电荷密度为：

$$\begin{array}{r}{\sigma }^{\alpha }=-(\frac{\delta \gamma }{\delta \Delta {\varphi }^0}{)}_{{\mu }_j}\equiv {\int }_{-\infty }^{\infty }p(\xi )\rho dx\end{array}$$

其中$p(\xi )={\xi }^3(6{\xi }^2-15\xi +10)$为一个插值函数，在电极相内$p(\xi )=1$在电解质侧$p(\xi )=0$。符号$\Delta {\varphi }^0$不同于$\Delta \varphi$，表示的是标准状态下评估的Galvani电势差（内电势）。因此微分电容也被定义为：

$$\begin{array}{r}{C}_d\equiv (\frac{\partial {\sigma }^{\alpha }}{\partial \Delta {\varphi }^0}{)}_{{\mu }_j}\end{array}$$

3. 热力学函数与材料参数

为了让相场模型更具体并允许数值计算，需要选择特定形式的热力学函数和材料参数。

3.1 热力学函数形式的选择

假设每单位体积亥姆霍兹自由能的化学势部分通过电极和电解质组分的两个理想溶液的插值来描述（电极也假设为理想溶剂相）。

f V ( ξ , C 1 , . . . , C n ) = 1 V m f m ( ξ , C 1 , . . . , C n ) = ∑ j = 1 n C j { μ j 0 α p ( ξ ) + μ j 0 β [ 1 − p ( ξ ) ] + R T ln ⁡ C j V m + W j g ( ξ ) } \begin{align} f_V(\xi,C_1,...,C_n)=\frac{1}{V_m}f_m(\xi,C_1,...,C_n)=\sum_{j=1}^nC_j\{\mu_j^{0\alpha}p(\xi)+\mu_j^{0\beta}[1-p(\xi)]+RT\ln C_jV_m+W_jg(\xi)\} \end{align} fV​(ξ,C1​,...,Cn​)=Vm​1​fm​(ξ,C1​,...,Cn​)=j=1∑n​Cj​{μj0α​p(ξ)+μj0β​[1−p(ξ)]+RTlnCj​Vm​+Wj​g(ξ)}​​

其中$C_j{V}_m=X_j$，${\mu }_j^{0\alpha }$和${\mu }_j^{0\beta }$分别是电极或电解质中纯组分$j$的化学势。和许多凝固相场模型一样，插值函数$p(\xi )={\xi }^3(6{\xi }^2-15\xi +10)$，双阱函数$g(\xi )={\xi }^2(1-\xi {)}^2$，每个组分的势垒高度为$W_j$（$W_j$约束太宽的界面，${\kappa }_{\xi }$约束太窄的界面），也可以选择其他的$p(\xi )$和$g(\xi )$，只需要满足$p(0)=0,p(1)=1,p^{\prime }(0)=p^{\prime }(1)=0$，$g^{\prime }(0)=g^{\prime }(1)=0$。

对于方程（12b）通过对$\xi$微分，可以得到$\partial f_V/\partial \xi$：

$$\begin{array}{r}\frac{\partial f_V}{\partial \xi }=p^{\prime }(\xi )\sum _{j=1}^n{C}_j\Delta {\mu }_j^0+g^{\prime }(\xi )\sum _{j=1}^n{C}_j{W}_j\end{array}$$

其中，$\Delta {\mu }_j^0\equiv {\mu }_j^{0\alpha }-{\mu }_j^{0\beta }$

同理可以对$C_j$进行微分，得到$\partial f_V/\partial C_j$：

$$\begin{array}{r}\frac{\partial f_V}{\partial C_j}={\mu }_j={\mu }_j^{0\alpha }p(\xi )+{\mu }_j^{0\beta }[1-p(\xi )]+RT\ln C_j{V}_m+W_{j}g(\xi )\end{array}\text{\hspace{1em}(27a)}$$

从而得到电化学势：

$$\begin{array}{r}\overline{{\mu }_j}={\mu }_j^{0\beta }+\Delta {\mu }_j^{0}p(\xi )+RT\ln C_j{V}_m+W_{j}g(\xi )+z_{j}F\varphi ,j=1,...,n\end{array}$$

可以看到，电化学势$\overline{{\mu }_j}$通过摩尔体积$V_m$依赖所有物种的$C_j$。

3.2 标准化学势

在数值模型中，需要知道每个物种的$\Delta {\mu }_j^0$。在平衡时，两相的电化学势相等（17）$\overline{{\mu }_j^{\alpha }}=\overline{{\mu }_j^{\beta }}$，结合代入（27），可以得到（其中在两相体相中，插值函数在两相不同$p(0)=0,p(1)=1$，而双阱函数都为$0$）：

$$\begin{array}{rl}\overline{{\mu }_j^{\alpha }}& ={\mu }_j^{0\beta }+\Delta {\mu }_j^0+RT\ln X_j^{\alpha }+z_{j}F{\varphi }^{\alpha },j=1,...,n\\ & ={\mu }_j^{0\alpha }++RT\ln X_j^{\alpha }+z_{j}F{\varphi }^{\alpha }\\ \overline{{\mu }_j^{\beta }}& ={\mu }_j^{0\beta }+RT\ln X_j^{\beta }+z_{j}F{\varphi }^{\beta },j=1,...,n\end{array}$$

$$\begin{array}{r}\Delta \varphi =-\frac{\Delta {\mu }_j^0}{z_{j}F}+\frac{RT}{z_{j}F}\ln \frac{X_j^{\beta }}{X_j^{\alpha }},j=1,...,n\end{array}$$

其中摩尔分数受电荷中性（18）和定义（2）约束。

人为指定标准内电势$\Delta {\varphi }^0$是贯穿摩尔分数$X_j^{\alpha 0}$和$X_j^{{\beta }^0}$界面的电位差：

$$\begin{array}{r}\Delta {\mu }_j^0=RT\ln \frac{X_j^{\beta 0}}{X_j^{\alpha 0}}-z_{j}F\Delta {\varphi }^0,j=1,...,n\end{array}$$

因此可以计算四个物种的$\Delta {\mu }_j^0$，其中$\Delta {\varphi }^0$是材料属性，给定电极、电解质系统其值是固定的。

在传统电化学中，会认为将电惰性物质的参考摩尔分数在一相或另一相中设置为$0$，$X_{e^{-}}^{0\beta }=X_{A^{-}}^{0\alpha }=0$。然后将电活性物质（$M^{m+},N^{n+}$）的电化学势等同起来，最后得到：

$$\begin{array}{r}{\epsilon }_{M^{m^{+}}}^0-{\epsilon }_{N^{n^{+}}}^0=-(\frac{\Delta {\mu }_{M^{m^{+}}}^0}{mF}-\frac{\Delta {\mu }_{N^{n^{+}}}^0}{nF})\end{array}$$

但出于数值目的，需要为电惰性物质的参考浓度假设为很小但是非零值。这就相当于假设很大（正或者负）但有限的标准电位。（实际也不为零）

数值计算系统处理：

• 以$M$金属电极和$1mol/L$ $M^{2+}$和$A^{2-}$的水溶液（$N$）溶液电解质作为目标。如（2）假设，将物种$M^{2+}、N^{2-}、N$的偏摩尔体积都取为纯水的偏摩尔体积$\overline{V_s}=0.018L/mol$ $1/\overline{V_s}=55.6mol/L$。
• 电惰性物质的值限制在$X_{e^{-}}^{\beta 0}=X_{A^{2-}}^{\alpha 0}=X_N^{\alpha 0}=10^{-6}$
• 根据体相的电中性和质量守恒（2），可以获得关系$X_M\equiv \frac{3}{2}{X}_{e^{-}}$，$X_{MA}\equiv 2X_{A^{2-}}$：
  • 在电极相中，满足$X_{M^{2+}}^{\alpha }=\frac{1}{3}$，$X_{e^{-}}^{\alpha }=\frac{2}{3}$
  • 在电解质相中，满足$X_N^{\beta }=1-{\bar{V}}_s(C_{M^{2+}}+C_{A^{2-}})$，$X_{M^{2+}}^{\beta }=X_{N^{2-}}^{\beta }={\bar{V}}_s{C}_{M^{2+}}={\bar{V}}_s{C}_{A^{2-}}$

因此化学势差$\Delta {\mu }_j^0$的独立部分数值：$\ln (X_{e^{-}}^{\beta 0}/X_{e^{-}}^{\alpha 0})=-13.41$，$\ln (X_{M^{2+}}^{\beta 0}/X_{M^{2+}}^{\alpha 0})=-2.919$，$\ln (X_{A^{2-}}^{\beta 0}/X_{A^{2-}}^{\alpha 0})=9.798$，$\ln (X_N^{\beta 0}/X_N^{\alpha 0})=13.78$

3.3 其他参数

为简化，假设各物种（除电子$e^{-}$）的势垒高度$W_{j=2,...,n}=W$相同且$W_{e^{-}}=0$，这消除了（26）最后一项对$C_j$的依赖性和（12a）对$W_j$的依赖性。

在不包含静电效应的凝固相场模型中，${\gamma }_{\xi }$和${\delta }_{\xi }$可以由$W_j$和${\kappa }_{\xi }$获得：

$$\begin{array}{rl}{\gamma }_{\xi }& =\sqrt{\frac{{\kappa }_{\xi }W}{18\overline{V_s}}}\\ {\delta }_{\xi }& =\sqrt{\frac{{\kappa }_{\xi }\overline{V_s}}{2W}}\end{array}$$

其中下标$\xi$表示内含静电贡献。取$W=3.6\times 10^{5}J/mol$，${\kappa }_{\xi }=3.6\times 10^{-11}J/m$，因此${\gamma }_{\xi }=0.2J/m^2$，${\delta }_{\xi }=3\times 10^{-11}m$。

假设介电常数不依赖相（$\xi$）和浓度（$C_j$），并取$ϵ(\xi )=78.49{ϵ}_0$。（依赖的可以参考程俊的研究）

4. 数值模型

4.1 驰豫平衡

驰豫平衡的解通过将平衡方程（12）构建为时间相关的偏微分方程：

∂ C j ∂ t = ∇ ⋅ { M j ∇ δ L δ C j } = ∇ ⋅ { M j ∇ [ μ j ˉ − V j ˉ V n ˉ μ n ˉ ] } ( i f   j ≠ 1 ) = ∇ ⋅ { M j ∇ [ ( Δ μ j 0 − Δ μ n 0 ) p ( ξ ) + R T ln ⁡ C j C n + ( z j − z n ) F ϕ + ( W j − W n ) g ( ξ ) ] } ( i f   j = 1 , V e − ˉ = 0 ) = ∇ ⋅ { M e − ∇ [ Δ μ e − 0 p ( ξ ) + R T ln ⁡ V s ˉ C e − 1 + V s ˉ C e − + z e − F ϕ + W e − g ( ξ ) ] } (33a) \begin{align*} \frac{\partial C_j}{\partial t}&=\nabla\cdot\{M_j\nabla\frac{\delta \mathcal{L}}{\delta C_j} \}=\nabla\cdot\{M_j\nabla[\bar{\mu_j}-\frac{\bar{V_j}}{\bar{V_n}}\bar{\mu_n}]\}\\ (if\text{ }j\not=1)&=\nabla\cdot\{M_j\nabla[(\Delta\mu_j^0-\Delta\mu_n^0)p(\xi)+RT\ln\frac{C_j}{C_n}\\ &+(z_j-z_n)F\phi+(W_j-W_n)g(\xi)]\}\\ (if\text{ }j=1,\bar{V_{e^-}}=0)&=\nabla\cdot\{M_{e^-}\nabla[\Delta\mu_{e^-}^0p(\xi)+RT\ln\frac{\bar{V_s}C_{e^-}}{1+\bar{V_s}C_{e^-}}\\ &+z_{e^-}F\phi+W_{e^-}g(\xi)]\} \end{align*}\tag{33a} ∂t∂Cj​​(if j=1)(if j=1,Ve−​ˉ​=0)​=∇⋅{Mj​∇δCj​δL​}=∇⋅{Mj​∇[μj​ˉ​−Vn​ˉ​Vj​ˉ​​μn​ˉ​]}=∇⋅{Mj​∇[(Δμj0​−Δμn0​)p(ξ)+RTlnCn​Cj​​+(zj​−zn​)Fϕ+(Wj​−Wn​)g(ξ)]}=∇⋅{Me−​∇[Δμe−0​p(ξ)+RTln1+Vs​ˉ​Ce−​Vs​ˉ​Ce−​​+ze−​Fϕ+We−​g(ξ)]}​(33a)

$$\begin{array}{r}j=1,...,n-1\end{array}$$

$$\begin{array}{rl}\frac{\partial \xi }{\partial t}& =-M_{\xi }\frac{\delta \mathcal{L}}{\delta \xi }=-M_{\xi }[\frac{\partial f_V}{\partial \xi }-{\kappa }_{\xi }{\nabla }^2\xi -\frac{{ϵ}^{\prime }(\xi )}{2}(\nabla \varphi {)}^2]\\ & =-M_{\xi }[p^{\prime }(\xi )\sum _{j=1}^n{C}_j\Delta {\mu }_j^0+g^{\prime }(\xi )\sum _{j=1}^n{C}_j{W}_j-{\kappa }_{\xi }{\nabla }^2\xi \\ & -\frac{{ϵ}^{\prime }(\xi )}{2}(\nabla \varphi {)}^2]\end{array}\text{\hspace{1em}(33a)}$$

其中$M_j$和$M_{\xi }$分别表示组分$j$的迁移率和相的迁移率（这个在动力学里面有讨论）。除了（33）以外，还需要满足泊松方程（12c）。

（33）和（12c）组合成总的方程组，并通过有限差分、有限元求解，（12a）和（A6）满足0.1%误差即认为达到平衡。

初始条件：

• 陡峭的界面分离本体电极和电解质：${\xi }^{\alpha }=1$，${\xi }^{\beta }=0$
• 指定电解质相的浓度：$C_{M^{+}}^{\beta }$
• 由于电中性，初始整个域：$\nabla \varphi =0$

边界条件：

• 电极侧$x=0$：$\mathbf{n}\cdot \nabla \xi =0$，$\varphi =0$，指定的$C_j$（参考3.2中浓度处理）
• 电解质侧$x=L$：$\mathbf{n}\cdot \nabla \xi =0$，$\mathbf{n}\cdot \nabla \varphi =0$，指定的$C_j$（参考3.2中浓度处理）

5. comsol建模

注意：

• 代数变量$p(\xi ),g(\xi )$需要调整求解器里面的代数变量设置，向后欧拉初始步长分数调小；
• 所有两相变量都需要用$D^{eff}=D^{s}p(\xi )+D^l(1-p(\xi ))$描述
• $M_{\xi }$控制反应极化大小，越大反应极化越小，电流密度越小。同理$D_j$或者说$M_j$也控制反应极化大小，详细讨论需要看文献II（动力学）
  仿真结果：
  

6. 文献

[1] J. E. Guyer, W. J. Boettinger, J. A. Warren, and G. B. McFadden, Phys. Rev. E 69, 021603 (2004).
[2] J. E. Guyer, W. J. Boettinger, J. A. Warren, and G. B. McFadden, Phys. Rev. E 69, 021604 (2004).