文献来源：

[1] Chen L, Zhang HW, Liang LY, et al. Modulation of dendritic patterns during electrodeposition: A nonlinear phase-field model[J]. Journal of Power Sources, 2015, 300: 376～38

[2] Hong Z, Viswanathan V. Phase-field simulations of lithium dendrite growth with open-source software[J]. ACS Energy Letters, 2018, 3(7): 1737～1743.

写在前面：说是笔记，其实也就是流水帐，本人相场学得也一般，博文单纯记录分享用，反正我自己的博客，我自己的服务器，想怎么写都没关系，不会像写论文，一个标点符号错了都不行。随便写写，也许会有很多错误，欢迎批评指正。如果有人看就就搞一搞排版，优化一下内容，没人看就这么放着，直接服务器到期吧。OK，下面正式开始。

论文[1]是相场法模拟锂枝晶的核心文献，很多做相关领域的仿真都会参考该文献，论文[2]中附带一个锂枝晶仿真的源码，该源码采用开源有限元仿真框架MOOSE运行。现在的问题是论文[1]我能看懂，但是百分百按照论文复现根本算不出来，论文[2]需要的巨势理论我不太熟悉，统计力学学完早忘了。不过源码能跑动。所以下文将主要讨论我对论文[1]的理解以及论文[2]的源码中遇到的一些问题，论文[2]的理论我直接放弃吧，懒得再把统计力学捡起来了。。。

1. 公式的推导

论文[1]中的模型整体来看是这三个方程耦合得到

$$
\frac{\partial \xi}{\partial t} = -L_0 \left( \xi'(\xi) – \kappa \nabla^2 \xi \right) – L_\eta h'(\xi) \left\{ \exp \left[ \frac{(1 – \alpha) n F \eta_a}{RT} \right] \right.\left. – \tilde{c}_+ \exp \left[ \frac{-\alpha n F \eta_a}{RT} \right] \right\} \tag{1}
$$

$$
\frac{\partial \tilde{c}_+}{\partial t} = \nabla \cdot \left[ D^{\text{eff}} \nabla \tilde{c}_+ + \frac{D^{\text{eff}} \tilde{c}_+}{RT} nF \nabla \phi \right] – \frac{c_s}{c_0} \frac{\partial \tilde{c}}{\partial t} \tag{2}
$$

$$
\nabla \cdot \left[ \sigma^{\text{eff}} \nabla \left( \phi(r, t) \right) \right] = I_R\tag{3}
$$

三个公式分别代表模型中的三个变量：序变量，阳离子浓度，电势的分布，这里重点需要提及的是公式(1)的推导过程，也就是论文[1]的附录部分。

使用相场法推导模型，最重要的部分就是自由能的推导，相场法的核心就是对自由能使用变分。一般来说相场法的论文只在开头提一下自由能，后面都是使用自由能密度这个词。这里简单讲一下，不要把这两个概念搞混。这是由于变分法是用于求解泛函的极值的，在相场法中泛函就是指自由能。形如$F[c(x,y)] = \int f(x, y, c,c_x,c_y) \, dV $这种带有积分形式的泛函被称为最简泛函，欧拉-拉格朗日是专门对付这种形式泛函求解极值的神兵利器，书上是这么解释的：本定理的重要作用在于，它把求泛函的极值问题转化为求解欧拉方程在满足边界条件下的定解问题。

欧拉-拉格朗日公式的形式如下：

$$
\frac{\delta}{\delta c} F = \left\{ – \nabla \frac{\partial}{\partial \nabla c} + \frac{\partial}{\partial c} \right\} f\tag{4}
$$

可以看到，变分的目标是针对泛函F求解极值，但是在求解过程中使用的是被积函数f，变分法中称其为核。也就是说F是目标，f是手段，而相场法中F是自由能，f是的自由能密度，这也导致在公式推导过程中，大量在说自由能密度而非自由能。

本文中自由能密度的来源分为两个，一个是电化学自由能密度，一个是界面（文中用的是gradient）能密度。

首先简单讲一下界面能，这里界面能一般是相场法用来打补丁的一个东西，在相场法的推导过程中，如果纯用扩散形式动力学公式分析，会发生从相分离系统出发，生成均相系统的矛盾，因此加入界面能这样一个纯数学形式的补丁（可能界面能有自己的物理意义吧，不过我没有查过）。大部分的相场论文中，界面能都是直接放上去的，不讲任何道理。

文中在求解电化学自由能密度时，使用的是$\text{粒子浓度}\times\text{化学势}$的形式文中的粒子有三种，阳离子、阴离子和不带电的原子，三种粒子的浓度全部作归一化处理。

先讲解一下化学势，一般来说不带电溶液中的化学势使用公式：
$$
\mu = \mu^\Theta + RT \text{ln}a \tag{5}
$$来计算，其中$\mu^{\Theta}$是标准化学势，$a$是活度。根据化学势的定义，溶液中物种B的化学势应理解为：在保持温度、压强和其他物种的物质的量不变且不做非体积功的情况下，单独向巨大的溶液系统中加入1mol物种B时，溶液Gibbs函数的增加。而在电解质溶液中，无法保持其他离子数量均不变的情况下单独加入1mol离子B，因为需要保持溶液的电中性，或者无法在不做非体积功的情况下单独加入1mol离子B，因为必须做电功。所以电解质的化学势的定义中去除掉不做非体积功的条件，电解质的电化学势计算公式为
$$
\tilde{\mu}_B = \mu_B + z_B F \Phi\tag{6}
$$式中第一项是非带电物质的化学势，也就是式(5)，第二项是电功。这里有一个假想过程，物种B是一种离子，先假想物种B不带电放入溶液中，此时增加的Gibbs函数就是$\mu_B$，在溶液中人为令其带电，此时增加的Gibbs函数就是电功$z_B F \Phi$。为了令非理想溶液服从Raoult定律，对浓度乘一个活度系数$\gamma$，得到式(5)中的活度，将考虑活度系数的影响下的化学势与标准化学势相减，就得到了论文[1]中的$\mu_i^{ex}$，超额化学势，化学势的计算公式中又恢复到了浓度的形式。

在了解化学势后，化学势乘上浓度就是自由能密度，上文中将化学势区分为普通溶液化学势和电功，原文中进一步将普通溶液化学势分为两类：中性原子的贡献与离子的贡献，文中假设不带电的原子全部沉积到电极表面，溶液中不存在中性粒子，即电极中$\tilde c=1$，溶液中$\tilde c=0$，$\tilde{c}$只有离散的0和1，不存在其他值。原文式(A7)中性原子对吉布斯自由能贡献公式为$ c_s RT\tilde{c}ln\tilde{c}$，该式在自变量为1或0时都取值为0。原文将离散的中性离子浓度$\tilde{c}$改成连续相场自变量$\xi$，它的吉布斯自由能贡献公式$ c_s RT\tilde{c}ln\tilde{c}$被改为双势阱函数$W\xi^{2}(1 – \xi)^{2}$。

至此，自由能中相场序变量的改造完成，自由能中其他部分没有改动。电能，界面能保持原样，Helmholtz自由能中原子的贡献被改为双势阱函数。

1.2 动力学分析（附录B）

这一节的动力学分析不是特别理解，在电化学课本中，动力学分析会使用Arrhenius公式，求出正向和逆向的反应速度并作差，得到的是电流-电势的特征关系式。在原论文中使用的是类似的思路，但是利用的是化学势与反应速度的关系式。

在代入有三种物种的化学势公式和Nernst公式后，将前面的化学势与反应速度的关系式转换成类似butlver-Volmer的形式，并且得到了$R_0$、$I_0$、$\bar{\mu}^{ex}_t$的求解公式，这些求解公式为下一节相场方程的推导做准备。

不过这一节确实有些粗糙，原文(A15)中$\bar{\mu}^{ex}_t$表示标准活化超额化学势，一般电化学表示活化状态使用符号是$\ddagger$，原论文使用的是$t$。我翻了参考文献和我的教材，都没有用$t$的，这就让我有点难受了。此外式(A23)到式(A24)很明显跳过了对活化能势垒的介绍，只给了一幅示意图。当然，仅仅是论文的附录，没有必要说得太详细，可以理解。

1.3 相场方程（附录C、附录D）

这是公式推导最核心的一节，也是写得最混乱的一节，我只能是猜测，可能是这节是作者改得太多了，所以公式的标号都是乱的，接下来我先把错误的公式标号改正，表格的备注中掺杂一些评论：

原论文错误位置和内容	修改内容	备注
公式(A30)下方：Subtracting Eq. (A28) from Eq. (A22)	从公式(A23)中减去公式(A30)
公式(A32)上方：In Eq. (A27), the activity for M-atom	在公式(A27)中，M原子的活度	公式(A32)中多出来的cs不知道是从哪里来的
公式(A33)	偏导应该改成变分
公式(A34)中的$\eta_\alpha$	$\eta_a$	$\eta_a$表示活化过电位，后文所有的活化过电位下标都用的是$\alpha$，式中已有$\alpha$用于表示对称系数，下标应该用$a$表示活化过电位
公式(A37)上方一段：Performing Taylor expansion on Eq. (A34)	对式(A36)进行泰勒展开
公式(A36)上方一段：we further write Eq. (A32) as	我们进一步将公式(A34)表示为
公式(A37)	第二行严格来说是用$\approx$	这里把$x$当成常数，以$y$为自变量，在$y=0$处泰勒展开
公式(A41)上面一段：Eq. (A30) into	公式(A33)
公式(A47)上面一段：Combining Eqs. (A43e44) yields	结合等式(A45-46)得到

相场公式的推导思路还是很清晰的，核心就是将动力学公式中改成关于粒子浓度的函数，将复杂的动力学公式简化为x和y，然后做泰勒展开，这样将动力学公式拆分成两类，一个是$R_\sigma$对应界面能的驱动力，一个是$R_\eta$对应于电化学反应的驱动力，这两部分分别转换为与迁移率相乘的形式后相加，就得到的序变量的时间导数。

这里$L_{\sigma}$是界面能项的迁移率，公式为：

$$
\begin{align*}
L_\sigma = R_0 \{(1 – \alpha) \exp[(1 – \alpha) x] + \alpha \exp(-\alpha x)\} / c_s RT\tag{7}
\end{align*}
$$

其中$x=\frac{n e \eta_\alpha}{R T} – \ln \tilde{c}_+$，$L_{\sigma}$是一个关于电势与离子浓度的函数。

$L_\eta$是电化学反应项的迁移率，公式为：

$$
\begin{align*}
L_\eta = k_0a_M^\alpha/\gamma_t\tag{8}
\end{align*}
$$

$L_\eta$是一个关于原子浓度的函数。

在原论文中，将这两者当作常数，那么当模型的浓度或电势与原论文不一致时，是否需要重新计算呢。

在附录D中，扩散方程中等号右侧第一项由C-H方程改写而来，第二项是Nernst-Planck方程中的迁移项改写而来，具体推导过程没有。原子项扩散系数视作0，针对阳离子与原子分别列出扩散方程，两个方程联立可得到公式(2)。

2. 求解

求解这一节主要还是问题，如果有了解的大佬麻烦指点一下。

2.1 针对论文[1]

首先是针对界面能项，界面自由能密度公式为$1/2\nabla \overrightarrow{\tilde{c}} \cdot \kappa \nabla \overrightarrow{\tilde{c}}$，但是这种形式的界面能完全无法长出来枝晶，只能长出柱状晶，因此，界面自由能采用$1/2\nabla \overrightarrow{\tilde{c}} \cdot \kappa^2 \nabla \overrightarrow{\tilde{c}}$这种形式是否更好，能否长出枝晶呢，这里我仅仅是借鉴小林雪花枝晶的模型，也没有什么理论基础。

在求解时，电化学反应项中，指数参数项中活化过电位数量级应在1e-1，因为$F\alpha/RT$的结果近似20，在这个数量级的exp结果在1e13，如果活化过电位比1大，该项将完全吞没相场方程其他项，方程就会不收敛，我想如果施加电压更大更过1，是否需要改其他常数，或者重新做归一化，在论文中的电势可是达到了-3V或-1.5V，且也有提到标准电势设为0的。

归一化参数方面的问题：界面迁移率、梯度能系数，电极中的扩散系数和溶液中的扩散系数在归一化后计算结果都与文章给出的数值不符。

2.2 针对论文[2]的代码

论文[2]的公式中，序变量公式完全借用的是论文[1]的序变量公式，而在求解浓度时，没有直接求解浓度，而是先求解化学势，再求解浓度，电场方程两者也是一样的。至于为什么做这种替换，论文是这样写的：

当ci→0时，由于熵，自由能曲线的斜率变得陡峭，ci中非常小的波动导致能量的大变化，导致数值不稳定。这种现象似乎限制了其他相场模型中可行的电解质成分的范围。然而，以μ为自然变量，能量变化对波动的敏感性要低得多，并且在低ci下更加稳健。

论文[2]的原理我看得不是特别懂，所以接下来就直接说源码了。

论文提供的源码更像是一demo，不像是一个完全体。在解读源码时，也有一些问题，放在这里：

在Mobility_coefficient和Migration_coefficient这两个块中，分别用于求解迁移和传导项的系数，对应原文公式(3)中的$Dc_{Li^+}$，代码中乘了两次1-h以限制该项只在液项中生效，乘一次1-h可能就足够了吧。
Bultervolmer块中对应电化学驱动力的那一项，源码中是这样写函数的：Ls*(exp(pot*AA/2.)+14.89*cl*(1-h)*exp(-pot*AA/2.))*dh，这里是将两项相加的，但是实际论文中的公式是两项相减的，不知为何代码与论文公式是不一致的。
材料参数中有一项A也不知为何物，他取值为1也不知有何作用。

以上，由于电脑太渣，跑一个模型要好久，调试很麻烦，很久才能看效果，所以有一些调试的想法也很难实现，就把想法放在这里了。另外我尽力对源作了一些微调，提高了收敛性，放在我的github中，github博客首页有放我的github链接，github中也有一些调试代码的建议，欢迎查看。