《储能科学与技术》推荐|李琛坤等：电化学阻抗谱物理模型求解方法

作者：李琛坤 1  王帅1 黄俊 2

单位：1. 中南大学化学化工学院；2. 乌尔姆大学理论化学研究所

引用： 李琛坤,王帅,黄俊.电化学阻抗谱物理模型求解方法[J].储能科学与技术,2022,11(03):912-920

DOI：10.19799/j.cnki.2095-4239.2022.0032

摘 要 电化学阻抗谱(EIS)是一种常用的电化学测试方法。通过EIS可以获取研究体系的各种动力学信息，比如电荷转移，双电层充电和物质传输等。EIS数据的解析依赖于合理的模型，常见的模型分为等效电路模型和物理模型。本文简要介绍了EIS物理模型三种常见的求解方法：解析求解，时域数值方法和频域数值方法。我们首先以一个经典的带有法拉第反应的双电层界面为例，仅考虑扩散控制传质的情况下，推导了相应的解析解，介绍了数值求解EIS的时域方法和频域方法，并对比了不同方法的结果。然后，我们以电池中常见的金属离子沉积反应为例，给出了零电荷电势下的EIS的解析解，并对比了两种数值方法的求解精度，同时我们也对比了不同价态的金属离子的EIS结果。计算结果表明，频域数值方法的精度高于时域数值方法，可以作为EIS数值求解的首选。保持其余参数不变的情况下，高价离子的EIS小于低价离子。最后，我们总结对比了三种方法的优缺点。本工作研究结果一方面可以作为EIS物理模型的入门指南，另一方面可以用于分析金属离子沉积反应的EIS结果。关键词 电化学阻抗谱；物理模型；数值方法

     电化学阻抗谱(electrochemical impedance spectroscopy，EIS)具有原位、无损的优点，被广泛应用于电化学储能(二次电池、燃料电池)和转化器件(电解池等)的性能监测、故障诊断和参数辨识。在EIS实验中，我们用电化学工作站先将测试体系控制在一个稳定状态，然后施加微弱的电压或电流激励，测试相应的电流或电压响应。傅里叶变换后的电压和电流的比值即定义为阻抗。

其中，表示傅里叶变换算子。测量不同频率下的阻抗便构成了EIS。对测量得到的EIS进行解析，需要相应的模型。常见的模型分为等效电路模型和物理模型^[。等效电路模型是指用电路元件(电阻、电容和电感等)的串并联来模拟研究体系。如图1所示，物理模型是指将研究体系内发生的物理过程用相应的控制方程来描述。求解控制方程，即得到EIS的解析解或数值解。李响等人详细推导了多孔电极EIS的表达式，并细致分析了不同过程对EIS响应的影响。等效电路模型的优点在于通用性好、计算量小和操作简单，缺点是物理机理不明确。程蕾等推导了玻碳电极界面EIS的等效电路模型。通过拟合实验数据，定量对比分析了修饰前后电极表面的特征变化。由等效电路模型衍生而来的弛豫时间分布(distribution of relaxation times，DRT)^[9]可以有效辨识不同过程的时间常数和阻抗，一定程度上辨析内部机理，因此得到了广泛的应用。物理模型的优点在于对体系内的物理过程描述清晰，适用于研究内部机理。基于控制方程建立的物理模型则有助于解析实验数据，辨析机理。因此，本文主要围绕物理模型展开介绍。图 1   (a) EIS测试原理；(b) 典型的EIS实验结果；(c) 物理模型构建原理(其中为交换电流密度，为浓度，为扩散系数)；(d) 2阶RC等效电路模型本文首先介绍了EIS物理模型三种常见的求解方法：解析求解，时域数值方法和频域数值方法。其中，解析求解指通过合理假设，推导出EIS的解析表达式。时域数值方法指利用有限元等方法在时间-空间维度求解微分方程，而后通过傅里叶变化得到EIS的数值解。频域数值方法指在频率-空间维度求解微分方程，而后计算EIS数值解。介绍完三种方法后，本文以电池中常见的金属离子沉积反应为例，对比了不同方法的优缺点，并对每种方法的适用场景进行了分析。

1 三种计算方法介绍

下文将以一个经典的带有法拉第反应的双电层界面为例，具体介绍上述三种方法。笔者^[10]此前对双电层的平衡、非平衡、交流状态下的建模进行了详细的推导和总结，可以作为初学者的指南，感兴趣的读者可以进行查阅。笔者^[11]同时梳理了EIS在电催化界面和反应中的研究和最新的进展。相比之下，本文重点关注EIS物理模型的求解方法。如图2(a)所示，在电极表面，处，发生氧化还原反应：。图2   (a) 带氧化还原反应的电极/电解质溶液界面示意图；(b) EIS解析解的Nyquist图；(c) 基于时域数值方法的EIS数值解与解析解的对比；(d) 基于频域数值方法的EIS数值解与解析解的对比。频率范围为～。其余参数如下：，，

，

Fig. 2   (a) Schematic diagram of electrode/electrolyte solution interface with a redox reaction; (b) Nyquist plot of analytical solution of EIS; (c) Comparison of numerical and analytical solutions of EIS based on TSM numerical method; (d) Comparison of numerical and analytical solutions of EIS based on FSM numerical method; Frequency range is 10⁵～10^-4 Hz. Other parameters are as follows: =0.1 mol/L, D=1×10^-10 m²/s, =1 A/m², α=0.5溶液中离子传输仅由浓度梯度驱动，控制方程为菲克第二定律。其中，为离子的扩散系数，表示氧化物，表示还原物。假设本体电解质位于无穷远处，，离子浓度为本体浓度

在电极表面，，反应物被消耗，产物被生成。消耗/生成的速率和反应速率的关系如下

     其中负号对应氧化物，正号对应还原物，氧化反应的电流定义为正电流，为法拉第常数，为转移的电子数，为法拉第反应的电流，一般通过Butler-Volmer方程计算

     其中，为交换电流密度，为电荷转移系数，为法拉第常数，为理想气体常数，为温度，为电极电势，为平衡电势，利用能斯特方程计算

其中为标准平衡电势，为活度，为表面浓度。总电流除了反应电流，还有双电层电流

     其中，为双电层电容，一般依据Gouy-Chapman-Stern模型计算

其中，为Gouy-Chapman电容，为溶液本体的介电常数，为德拜长度。为Helmholtz电容，和分别为Helmholtz层的介电常数和厚度。

1.1　解析求解

对方程式(2)施加傅里叶变换可以得到频域的控制方程其中，表示为频域下的扰动浓度。本文“～”表示扰动量。傅里叶变换后的边界条件如下，本体处

电极表面，处

结合边界条件式(10)和式(11)，可以得到方程式(9)的解析解为

其中正号对应氧化物，负号对应还原物，为频率依赖的扩散长度。对Bulter-Volmer方程做线性化展开，利用傅里叶变换可以得到

其中可以通过对能斯特方程做微分得到

^{在稀溶液中，，把方程式(14)代入方程式(13)得到}

^{对方程式(7)做傅里叶变换，可以得到频域的双电层电流}

^{基于总电流，和阻抗定义，我们可以得到}

^{其中为电荷转移电阻，为Warburg扩散阻抗，}

^。

至此，我们得到了该体系EIS的解析表达式。为明确其在各个频率区间的特征，我们可以做如下渐进分析。在高频区域，，方程式(17)简化为，即电阻和电容并联。因此，高频区域的阻抗为半圆。在低频区域，，，方程式(17)简化为，即低频区域为45°斜直线。

1.2　时域数值方法

明确了控制方程和边界条件后，我们可以在数值计算软件(MATLAB、Maple等)中利用有限元等方法对方程式(2)在时间-空间维度(time-space method，TSM)进行求解^[14-15]，把求解后的浓度带入Butler-Volmer方程即可得到不同时间下的法拉第反应电流。双电层电流的计算独立于浓度，可以直接进行计算。得到总电流后，对其做解析傅里叶变换(analytical fourier transform，AFT)可以得到频域下的总电流^[16]。对施加的交流电压做AFT变换得到频域下的电压。AFT变换后的电压和电流的比值即为阻抗，。这里简要介绍一下AFT。傅里叶变换的定义如下
其中，为系统的时域响应信号，为傅里叶变换算子。对信号做线性插值处理可以得到
其中，为归一化的阶跃函数。把方程式(19)代入方程式(18)可以得到AFT的表达式基于TSM求解，再通过AFT变换得到的EIS数值解和解析解的对比如图2(c)所示。在高频区域，数值解和解析解基本完全重合。低频区域的偏差略有增加，主要原因是有限元等方法求解微分方程会累计误差。

1.3　频域数值方法

相比TSM求解，再做AFT，频域数值方法则是先对控制方程和边界条件做傅里叶变换，然后通过有限元等方法在频域求解(frequency-space method，FSM)，省去了后续对时域电流和电压的AFT。在解析求解部分，已经对控制方程和边界条件做了傅里叶变换。因此，方程式(9)即为频域的控制方程，方程式(10)和式(11)为对应的边界条件。数值求解方程式(9)得到频域下的浓度，代入方程式(13)便得到了频域的法拉第反应电流，加上频域下的双电层电流，即为频域下的总电流。阻抗定义为交流电压幅值和频域下的总电流的比值，/J。其中，上标amp表示幅值。基于频域数值方法得到的EIS数值解和解析解的对比如图2(d)所示(具体的MATLAB程序放在支撑文件)。可以明显看到，数值解在整个频率范围内都和解析解很吻合，表明FSM数值解精度较高。

2 求解案例

2.1　金属离子沉积反应体系

如图3(a)所示，考虑一维平板电极，氧化物为金属离子M⁺，还原物为金属原子M。在Helmholtz平面发生法拉第反应：。电解质溶液为稀溶液，阴阳离子的化合价均为单价，且具有相同的扩散系数。图3   (a) 金属离子沉积反应电极/电解质溶液界面示意图；(b) EIS解析解的Nyquist图；(c) 时域数值方法的EIS数值解与解析解的对比；(d) 频域数值方法的EIS数值解与解析解的对比。频率范围为～。其余参数如下：，，，2.1.1　控制方程溶液中离子的传输采用Poisson-Nernst-Planck(PNP)方程描述。其中，和分别为阳离子和阴离子的浓度，为扩散系数，为本体溶液的介电常数。2.1.2　边界条件和初始条件研究体系的边界分为左边界和右边界，左边界为Helmholtz平面，处。边界条件设置如下，对于阳离子，发生电子转移反应，阳离子通量和反应电流的关系见式(24)

其中，为法拉第反应的电流，这里采用Frumkin-Butler-Volmer(FBV)方程计算

其中，为阳离子的本体浓度，为反应速率常数，为反应的过电势。，为Helmholtz平面处的电势。对于阴离子，不参与反应，所以离子的通量为0

对于电势，从电极表面到Helmholtz平面，电势线性下降，所以电势的边界条件为

     由于所有的电荷均为体电荷，Helmholtz平面体积为0，不存在净电荷。因此，根据高斯定律，Helmholtz平面两侧的电势梯度满足如下关系

方程式(28)代入方程式(27)可得电势左边界的最终形式

右边界为溶液本体，，电势为0，离子浓度均为本体浓度

结合边界条件对控制方程进行求解，我们可以得到零电荷电势(Potential of zero charge，PZC)下的解析表达式

     其中，为Gouy-Chapman电容和Helmholtz电容的比值，为无量纲的角频率，为无量纲的电解质溶液厚度，，，和为沉积反应相关的变量。

其‍中‍，exp，，为稳态的过电势，上标0表示该变量为稳态变量。2.1.3　数值求解需要注意的是，在金属沉积反应中，边界条件中严格考虑了双电层的存在。因此，双电层电流不能再以为常数进行计算，而应该从电荷量求导去计算其中，表示双电层中积累的净电荷量，负号表示积累的电荷和电极所带的电荷相反。总电流的表达式如下

     基于频域方程的数值求解则先要对控制方程和边界条件做傅里叶变换，具体的计算详情这里不再赘述。频域下的控制方程(这里只考虑PZC下的情况)如下

在处，频域下的左边界条件为其中，为频域的反应电流，。在处，频域下的右边界条件为2.1.4　三种方法结果对比图3(b)～(d)依次为金属沉积反应EIS解析解，时域数值方法的EIS数值解与解析解的对比，和频域数值方法的EIS数值解和解析解的对比(具体的MATLAB程序放在支撑文件)。需要注意的是，这里的解析解和数值解都只适用于PZC。对于非PZC下的求解，初始条件不再满足均匀假设，离子浓度和电势均存在分布。特别是在双电层中，电势梯度和浓度梯度极大。笔者尚未看到有相关研究在考虑界面反应和双电层的条件下，自洽求解了非PZC下的EIS。从图3(b)解析解的结果可以看出，EIS由高频反应半圆和低频扩散半圆构成。解析解作为精确解，用来评价不同数值方法求解EIS的精度。图3(c)为时域方法的数值解和解析解。二者在中低频区域吻合较好，高频区域数值解出现感抗行为。从笔者的经验来看，主要是因为计算双电层电流采用的数值积分及微分存在较大的误差，进而导致高频感抗的行为。图3(d)对比了频域方法的数值解和解析解。不管是高频还是低频，数值解和解析解均很好吻合，表明数值解精度较高。从三种方法对比结果可以得出，解析解避免了数值求解微分方程，因而精度最高。但解析解往往需要大量数学运算，才有可能推导出复杂体系的解析EIS。相比之下，数值解借助有限元等方法，可以直接得到体系的EIS数值解。对比不同数值方法的结果，可以得出，频域数值方法的EIS数值解精度较高，而时域数值方法的EIS数值解，低频区域精度较高，高频区域精度较低。因此，笔者建议优先选取频域数值方法计算EIS。

2.2　单价离子和多价离子阻抗的对比

金属离子沉积反应中，常见的离子价态为一价、二价和三价。图4对比了不同价态离子的EIS(具体的MATLAB程序放在支撑文件)。其中，图4(a)和(b)分别为TSM数值方法和FSM数值方法的结果。在TSM数值方法的EIS结果可以看出，高频仍然存在感抗特征，其余部分和FSM方法的结果一致，这也说明了FSM方法的求解精度更高。对于不同价态的离子，在保持其余参数不变的情况下，离子价态越高，阻抗越小。主要原因在于，一方面FBV方程中考虑离子价态后，反应电流极大地增加。另一方面，通过积分电荷量计算的双电层电流也增加。因此，总电流极大的增加使得高价态离子的EIS减小。但需要注意的是，对于真实的研究体系，不同价态离子的动力学参数往往差别较大，读者需要结合实验标定相关参数，然后利用模型去拟合和解析实验数据。图4   (a) 基于TSM数值方法的EIS结果对比；(b) 基于FSM数值方法的EIS结果对比(除考虑不同离子价态外，其余参数均与图3保持相同)

3 三种方法优缺点对比

图5对比了三种求解方法的优缺点。其中，解析求解优点在于物理含义清晰，易于分析高低频的渐进行为。缺点在于推导复杂，复杂体系难以解析求解。数值方法优点在于求解简单，无复杂推导。缺点在于计算耗时。其中，时域数值方法求解精度低于频域数值方法。但频域数值方法无法获得时域信号，难以有效判断中间变量的正确性。图5   三种求解方法优缺点对比

4 结论和展望

本文简要介绍了EIS物理模型常见的三种求解方法：解析求解，时域数值方法和频域数值方法。随后结合具体的例子对比了不同方法的优缺点，得到的主要结论如下：解析解精度最高，但一般只适用于PZC条件；对于非PZC，考虑双电层的条件下，无法求出解析解。对于数值解，优先推荐频域数值方法，精度较高且没有高频感抗行为。为拓宽EIS在电化学中的应用，尚有很多工作值得开展。（1）本文的解析解和数值解都只适用于PZC条件，对于非PZC下的EIS求解，笔者尚未看到有相关研究。如何自洽求解非PZC下的EIS是一个值得研究的课题。（2）本文只考虑了单电子转移的反应，而燃料电池、电解水等电化学装置中的反应均为4电子和2电子反应。分步考虑电子转移反应，建立微观动力学模型，再计算EIS对于理解反应速控步骤，辨析反应机理等至关重要^[24-25]。（3）本文仅考虑了一维平板电极情况，而真实的研究体系多为粗糙界面、团聚体颗粒和多孔电极。Song等人和Kant等人的研究表明，几何效应对EIS响应有很大的影响。

通讯作者：李琛坤（1997—），男，硕士研究生，研究方向为理论电化学，E-mail：[e-mail]。

相关文章

邮发代号：80-732

联系热线：010-64519601/9602/9643

投稿网址：http://esst.cip.com.cn/CN/2095-4239/home.shtml