CN102661983B - 一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法 - Google Patents

Abstract

本发明公开了一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法，本发明属于电化学阻抗谱的拟合方法，特别涉及一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法。可根据不同的需要选择矩阵法或串-并联法，均能取得相同的结果。以均匀设计法对拟合参数的初始值进行选择，大大提高了运算的速度与时间。对应于反应的可逆与否，扩散层是半无限、有限还是块层扩散类型，Warburg阻抗应用不同的表达公式。

Description

一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法

技术领域

本发明属于电化学阻抗谱的拟合方法，特别涉及一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法。

背景技术

电化学阻抗谱The electrochemical impedance impedancespectroscopy(EIS)或者交流阻抗谱是得到广泛应用的研究电极表面电化学机理的技术。由于这种技术是测量稳态状况下的电化学反应，应用较宽的频率域，在每个周期内测量信号的随时间变化的平均值。

多孔型材料是各类型电池和燃料电池的主要电极材料。这是由于其为电化学反应提供了其巨大的内表面积。然而，随着多孔型材料颗粒的变小，相应的使用效率并未随之提高。造成这种现象的原因部分来自于电化学有效使用面积受界面上的累积的溶液电阻控制。而这种电阻是由离子孔内的扩散决定的。另一方面，也是由各个颗粒之间的累计接触电阻，传导电阻和法拉第电阻等影响。这些都是多孔型电极的内在的电极动力学障碍。

通常情况下，一个多孔型电极的所有内表面不能全部被利用到。由于各种孔的形状，孔径，孔深不同，这些结构对外加的激励信号的响应时间不同。因此在小孔内离子的运动比在大孔中要迟缓，这就造成离子在多孔型电极中扩散速度与在本体溶液中不同。这种不同随孔类型和尺寸大小而变化。

传输线模型，The Transmission Line model(TML)是由DeLevie提出的用于解释多孔型电极的电极过程。根据这一模型，在孔的不同的深度对外加的激励信号有不同的响应。单位距离(从孔口)的电解质电阻和双电层电容均匀沿孔壁分布。一个理想孔的电极行为可以用均匀的电阻-电容传输线模型来描述。这个模型很好地描述了从孔口到孔底的一个连续增长的电阻过程。这个模型可用于各种不同类型的电极：如半导体电极，薄层电极，燃料电池，多孔型电池电极，聚合物与聚合物包裹电极等。改进型传输线模型，例如孔尺寸分布传输线模型Pore Size Distribution(PSD)TML被用于模拟只涉及非法拉第过程的多孔型电极(如石墨，金等)的阻抗谱。

本法中使用传输线模型对应的等效电路如图1所示。

其中R_s表示单位长度电解质溶液电阻(Ω/cm)，R_ct是电荷转移电阻(Ω·cm²)，Z_w是Warburg阻抗(Ω·cm²)，C_d1双电层比电容(F/cm²)。图中Z单元沿孔壁分布，表示电极表面与电解质溶液界面的电解过程，其单位是单位面积的比阻抗(Ω·cm²)。一个孔的阻抗，Z_pore，可由公式(1)计算：

$Z_{\mathrm{pore}}=\sqrt{R_{s}Z}\·\coth \left(l_p\sqrt{\frac{R_s}{Z}}\right)---\left(1\right)$

其中1_p在等价园柱体结构中孔的深度。

在De Levie的传输线模型中，Z单元的阻抗在孔壁的各处都相同。由此而导出的公式(1)用于计算一个孔的阻抗。因而可用复数非线性最小二乘法(CNLLS)来进行阻抗谱拟合。这样的拟合方法是常用的对多孔型电极的阻抗谱的拟合方法。这样的方法本身包含着不合理的假设，因而运用CNLLS拟合法不能达到较好的拟合精度，而且往往不能完成全频率域的拟合。本法中，为了应用CNLLS方法作为对比，用公式(2)表明C_s与C_b的关系：

$C_s=\frac{C_b}{N_c}---\left(2\right)$

然而，在交流阻抗的测量中，如果激励信号是正玄波电势，相应的等价体是多孔圆柱体，电势沿孔壁的分布，根据De Levie的推导，按如下规律(3)：

$V=V_0{e}^{(-x\·\sqrt{\frac{R_s}{Z}})}\·\sin \left(\mathrm{\ωt}\right)---\left(3\right)$

其中ω角频率，V₀波幅，R_s电解质溶液电阻，Z等效电路中Z单位的阻抗，x从孔口到孔底的某深度。

公式(3)表明电势沿孔壁呈现指数下降分布。这种分布使得活性电解质离子在某一频率下在孔壁各处有不同的电化学反应速率：离孔口越远，速度越慢，相应的电流也是越深入孔底越小。因此阻抗也随着孔深而变化。本法中，电解质溶液的电导率被认为不随孔深而变化，i、e、R_s是常数。溶液相的扩散是造成Z是孔深x的函数的原因。本法从理论推导出的C_s与C_b的关系式而发展出数值迭代法计算Z阻抗，进而发展出数值迭代拟合法。

发明内容

本发明的目的在于克服现有技术存在的以上问题，提供一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法。

为实现上述技术目的，达到上述技术效果，本发明通过以下技术方案实现：

一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法包括以下步骤：

步骤1)选择几何等价结构；

步骤2)选择Z_w模式，所述Z_w表示Warburg阻抗；

步骤3)为参数分配值域；

步骤4)均匀设计表；

步骤5)为参数赋初值；

步骤6)设i为角频率的序号，i初值为1，所述ω_i表示第i个角频率；

步骤7)在某一角频率ω_i下，使用数值迭代法计算Z₁到Z₁₀；

步骤8)通过Z₁，Z₂...Z₁₀，使用Mathematica中的内插值命令计算Z(x)；

步骤9)将Z(x)变换为每个Z单元的值，使用矩阵法或串并联法得出Z_proe，所述Z_proe表示等价体中一个孔的阻抗；

步骤10)Z_total＝Z_proe/N，所述Z_total表示等价体中总阻抗，所述N表示等价体中包含孔数；

步骤11)判断相位角是否接近其实验值，是，进入下一步，否，调整第一参数(C_d1)并返回步骤7，总是不，根据均匀设计表选择下一套参数初始值赋值方案，返回步骤5；

步骤12)判断实部是否接近其实验值，是，进入下一步，否，调整第二参数(D)并返回步骤7，总是不，根据均匀设计表选择下一套参数初始值赋值方案，返回步骤5；

步骤13)i＝i+1，返回步骤6。

进一步的，步骤7中所述循环迭代数值法计算Z_i包括以下步骤：

步骤a)将C_b作为C_s的起始值；

步骤b)通过

$I_F=\mathrm{nF}\·k\·C_s^{\mathrm{\γ}}$

${\mathrm{\σ}}_w=\frac{\mathrm{\γ}{R}_{\mathrm{ct}}\left|I_F\right|}{n{\mathrm{FC}}_s\sqrt{2\mathrm{\ωD}}}$

计算得I_F及σ_w，所述C_b表示电解质溶液本体浓度，所述C_s表示电解质离子在电极表面的浓度，所述I_F表示法拉第电流密度，所述F表示法拉第常数，所述k表示电化学反应常数，所述γ表示电化学反应级数，所述σ_w表示Warburg阻抗表达式前置因子，所述R_ct表示电荷转移电阻，所述ω为ω_i表示第i个角频率，所述D表示扩散系数；

步骤c)计算得Z_w，所述Z_w表示Warburg阻抗；

步骤d)通过

$Z=\frac{1}{\frac{1}{R_{\mathrm{ct}}+Z_w}+\mathrm{j\ω}{C}_{d1}}$

计算得Z，式中Z单元的阻抗，Z_w表示Warburg阻抗，R_ct表示电荷转移电阻，C_d1表示双电层电容，j表示虚数单位，ω_i为第i个角频率；

步骤e)通过

$V=V_0{e}^{(-x\·\sqrt{\frac{R_s}{Z}})}\·\sin \left(\mathrm{\ωt}\right)$

计算得V，式中V表示所加正弦波电压，V₀表示V的波幅，R_s表示对应于在传输线模型等价电路中的电解质溶液电阻，Z为Z单元的阻抗，ω_i为第i个角频率；

步骤f)通过

$k=\frac{i_0[e^{-\mathrm{\αf\η}}-e^{(1-\mathrm{\α})\mathrm{f\η}}]}{\mathrm{nF}\·C_s^{\mathrm{\γ}}}$

计算得k，式中k表示电化学反应常数，i0表示交换电流密度，α表示对称因子，F表示法拉第常数，η表示过电势，C_s表示电解质离子在电极表面的浓度，γ表示电化学反应级数，f＝F/RT，F表示法拉第常数，R表示气体常数，T表示绝对温度；

步骤g)通过

$C_s=C_b\·\mathrm{Csch}(\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}}\·l_p)\·\mathrm{Sinh}[\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}}\·(l_p-x)]$

计算得C_s，与之前所得C_s值进行比较，若前后两次C_s值之差小于10^-3则认为C_s已经达到稳定并进入步骤h，否则进入步骤b；

步骤h)计算得Z_w，所述Z_w表示Warburg阻抗；

步骤i)通过

$Z=\frac{1}{\frac{1}{R_{\mathrm{ct}}+Z_w}+\mathrm{j\ω}{C}_{d1}}$

计算得Z，式中Z为Z单元的阻抗，Z_w表示Warburg阻抗，R_ct表示电荷转移电阻，C_d1表示双电层电容，j表示虚数单位，ω_i为第i个角频率。

进一步的，所述步骤c与步骤h中计算Z_w对于半无限扩散类型公式为

Z_w＝σ_w·(1-j)

对于有线扩散类型公式为

$Z_w={\mathrm{\σ}}_w\·(\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}-j\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)})$

对于块层扩散类型公式为

$Z_w={\mathrm{\σ}}_w\·(\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}-j\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)})$

式中Z_w表示Warburg阻抗，σ_w表示Warburg阻抗表达式前置因子，ω_i为第i个角频率，D表示扩散系数，δ表示扩散层厚度，j表示虚数单位。

进一步的，步骤9中对于包含由较大单元数目的传输线模型，通过电流在每个单元的分布遵守基耳霍夫原理，使用所述矩阵法计算出一个孔的电化学阻抗；对于包含由较大单元数目的传输线模型，通过等效电路中各个部件的串并联关系，使用所述串并联法来计算出一个孔的电化学阻抗。

进一步的，步骤11中所述第一参数为C_d1，调整方式为当相位角大于其实验值，C_d1＝C_d1+C_d1/100；当相位角小于其实验值，C_d1＝C_d1-C_d1/100，式中C_d1表示双电层电容。

进一步的，步骤12中所述第二参数为D，调整方式为当实部大于其实验值，D＝D+D/100；当实部小于其实验值，D＝D-D/100，式中D表示扩散系数。

本发明的有益效果是：

可根据不同的需要选择矩阵法或串-并联法，均能取得相同的结果。以均匀设计法对拟合参数的初始值进行选择，大大提高了运算的速度与时间。对应于反应的可逆与否，扩散层是半无限、有限还是块层扩散类型，Warburg阻抗应用不同的表达公式。

上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，并可依照说明书的内容予以实施，以下以本发明的较佳实施例并配合附图详细说明如后。本发明的具体实施方式由以下实施例及其附图详细给出。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1传输线模型等效电路用于数值迭代拟合法；

图2计算V，k，C_s，Z_w，和Z的循环图；

图3循环数值法计算Z1；

图4|C_s|的值随孔深度的变化；

图5图示了Z(x)与孔深度的关系和计算Z_i的方法；

图6电流在传输线模型等价电路中的分布；

图7数值迭代拟合法拟合多孔型电极的电化学阻抗谱的流程图；

图8实施例1中使用CNLLS法拟合值与实验值的Nyquist图；

图9实施例1中使用数值迭代法拟合值与实验值的Nyquist图；

图10实施例2中使用CNLLS法拟合值与实验值的Nyquist图；

图11实施例2中使用数值迭代法拟合值与实验值的Nyquist图；

图12实施例3中使用CNLLS法拟合值与实验值的Nyquist图；

图13实施例3中使用数值迭代法拟合值与实验值的Nyquist图。

具体实施方式

下面将参考附图并结合实施例，来详细说明本发明。

1.数值迭代拟合法有以下假定：

(1)所有的等价多孔圆体都有相同尺寸大小的孔(半径r₀，深度l_p)，这些孔都相互平行排列，总孔数目是N。

(2)电极材料没有电阻。

(3)等效电路中每个单元在等效圆柱体中有相同的高度。

(4)电极过程包括法拉第和非法拉第过程。

(5)正玄波电势作为扰动信号。

(6)扩散采用平面扩散模式。

2.C_s与C_b数学关系式的推导和对Warburg阻抗Z_w表达式的影响

因为电势沿孔壁从孔口到孔底的变化是孔深度x的函数，电活性离子从本体溶液向孔壁扩散按照费克第二定律：

$\frac{{\∂C}_s}{\∂t}=D\frac{{\∂}^2{C}_s}{{\∂x}^2}-\frac{2I_f}{\mathrm{nF}{r}_0}{C}_s^{\mathrm{\γ}}---\left(4\right)$

解这个方程，边界条件是C_s(0)＝C_b and C_s(1_p)＝0，解出的C_s表达式是：

$C_s=C_b\·\mathrm{Csch}(\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}}\·l_p)\·\mathrm{Sinh}[\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}}\·(l_p-x)]---\left(5\right)$

公式中k电化学反应常数，V激励电势或过电势，D是扩散系数，n电化学反应转移电子数，R_ct是电荷转移电阻，r₀和l_p分别是等价体孔半径和孔深。方程(4)先进行Laplace变换使其从时间域变成频率域然后解得公式(5)。

对于一个不可逆反应，Z_w的前置系数是：

${\mathrm{\σ}}_w=\frac{\mathrm{\γ}{R}_{\mathrm{ct}}\left|I_F\right|}{n{\mathrm{FC}}_s\sqrt{2\mathrm{\ωD}}}---\left(6\right)$

$I_F=i_0[e^{-\mathrm{\αf\η}}-e^{(1-\mathrm{\α})\mathrm{f\η}}]=n{\mathrm{Fv}}_r=\mathrm{nF}\·k\·C_s^{\mathrm{\γ}}---\left(7\right)$

式中i₀为交换电流密度，其中R是气体常数，T是Kelvin温度，n是电子转移数，F是法拉第常数。

式中γ电化学反应级数，|I_F|是在图1中通过R_ct-Z_w分支的法拉第电流密度的大小，C_sd电活性电解质离子在孔表面的浓度，C_b电活性电解质离子在本体溶液中的浓度。

继而根据对不同扩散类型，对不可逆电化学反应的Warburg阻抗(Z_w)可计算获得：

对半无限扩散类型：

Z_w＝σ_w·(1-j)             (8)

对有限扩散类型：

$Z_w={\mathrm{\σ}}_w\·(\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}-j\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)})---\left(9\right)$

对块层扩散类型：

$Z_w={\mathrm{\σ}}_w\·(\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}-j\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)})---\left(10\right)$

其中ω是角频率，D是扩散系数，δ是扩散距离，σ_w是Z_w的前置系数。其根据反应的可逆或不可逆性质有不同的表达式。

当把公式(5)代入到公式(6)，σ_w就成了孔深，x的函数。如此Z_w和Z也成了x-的函数。

3.Z的数值计算法

在图1中，Z单元的阻抗可由公式(11)计算：

$Z=\frac{1}{\frac{1}{R_{\mathrm{ct}}+Z_w}+\mathrm{j\ω}{C}_{d1}}---\left(11\right)$

公式(11)中要计算Z，C_d1，R_ct，和Z_w的值必须已知。对一个不可逆电化学反应，它的Z_w可从公式(8)到(10)计算而来。但是公式(5)不能直接用于C_s的计算。

在公式(5)中，C_s由η(V)和k，决定。k是电化学反应的反应常数，其通过公式(12)计算获得。公式(12)从Butler-Volmer公式推导而来：

$k=\frac{i_0[e^{-\mathrm{\αf\η}}-e^{(1-\mathrm{\α})\mathrm{f\η}}]}{\mathrm{nF}\·C_s^{\mathrm{\γ}}}---\left(12\right)$

式中α对称因子，f＝F/RT。

从公式(3)可知，η(V)由Z决定，而Z由Z_w通过公式(11)决定。Z_w又是由C_s从公式(6)至(10)决定，C_s则由k和η(V)通过公式(5)决定，同时k由η(V)和C_s通过公式(12)决定，再次，η(V)由Z通过公式(3)决定...，这样的计算导致了循环嵌套，如图2所示。

对于一个孔深x，在某一角频率下ω，当达到稳态后，电活性电解质离子有一个稳定的值。从此原理出发，可发展出一个数值迭代法用于计算在某一孔深，在某一角频率下的值。作为起始，C_b被当成C_s的起始值。然后通过公式(6)，公式(7)(8)/(9)和公式(10)计算Z_w。然后运用公式(11)计算Z。获得Z后通过公式(3)计算V，再从V通过公式(12)计算k。当k值已知则通过公式(5)计算而得第一次循环所得C_s值。用这个C_s值作为初始值再进行下一次循环得到第二次循环C_s值。如果前后两次C_s值之差小于10^-3，则认为C_s已经达到稳定值。否则再进行一次循环，直到C_s稳定值的要求达到。这一过程的流程图如图3。

4.C_s收敛性的证明

在这个数值迭代法计算最终的C_s值不是通过一个公式而得来。需要证明在某一角频率下，在孔壁的某一深度处，其C_s值有只有一个稳定的最终值。本法中，当前后两个C_s值相差不超过10^-3时则认为C_s达到了稳定值。大多数情况下需要循环至少三次C_s才能达到稳定。因此需要证明在孔壁的更多的深度C_s都能达到稳定值。

图4表明了|C_s|值在十个孔深处的变化。对于某一点，其值经过循环计算后逐渐趋向一个最终值。箭头指向循环计算所得值的方向。从孔口到孔底|C_s|逐渐变小。这是因为随着深度的增加，电势减小导致法拉第电流也减小。电化学反应速率也相应减小，因此|C_s|变小。在图4中，箭头指向每次循环计算后|C_s|的变化方向。对所有的十个点从0.11_p到0.951_p，显示了相同的变化趋势。这就说明每个点都有其最终的稳定值。从0.11_p到0.951_p，每个点|C_s|的最终稳定值逐渐变小。这符合前面的分析。

除了图4中所显示的十个点，其余的点，包括0.151_p，0.251_p，0.351_p，0.451_p，0.551_p，0.651_p，0.751_p，0.851_p，0.131_p，0.231_p，0.331_p，0.431_p，0.531_p，0.631_p，0.731_p，0.831_p，0.931_p，0.181_p，0.281_p，0.381_p，0.481_p，0.581_p，0.681_p，0.781_p，0.881_p，和0.981_p，它们的|C_s|最终的稳定值也进行了计算它们也都符合相同的规律。从以上的各点的|C_s|的计算表明，我们可以认为通过数值迭代计算后，所有的孔深度处在某一角频率下都有一个最终稳定的C_s。

5.数值Z(x)的构建和Z_i的计算

Z是等价电路中的Z单元的阻抗。如前所述其是孔深度，x的函数。在某一ω下，选择沿孔壁的十个深度点：0.11_p，0.21_p，0.31_p，0.41_p，0.51_p，0.61_p，0.71_p，0.81_p，0.91_p，和0.951_p来计算各自的Z阻抗值。运用Mathematica软件包的“Interpolation”(插值)命令，Z的随x变化的函数：

Z(x)＝Interpolation[{Z1，0.11_p}，{Z2，0.21_p}，{Z 3，0.31_p}，{Z4，0.41_p}，{Z5，0.51_p}，{Z6，0.61_p}，{Z7，0.71_p}，{Z8，0.81_p}，{Z9，0.91_p}，{Z10，0.951_p}]

内插值命令的工作原理是在连续的点间进行多项式曲线拟合。Wolfram公司(Mathematica的开发商)的技术文件详细介绍了内插值命令的工作原理。

Z(x)是一个随x连续变化的函数。对每个Z单元，其值Z_i，由Z_i+1和Z_i的差值获得。角标i是Z单元的序列号，从孔口到孔底。在从Z(x)计算Z_i，x要首先变成Z单元的序列号。如图5所示，在Z(x)通过运行Mathematica的Interpolation命令获得后，第i个Z的值从Z_i+1-Z_i获得。i是Z单元的序列号。

6.Z_pore和Z_total的计算

图1中，一个孔的阻抗(Z_pore)是由R_s和Z决定的。因为Z是孔深x的函数，公式(1)不能用来计算Z_pore。当每个Z单元的阻抗值Z_i和溶液电阻R_s已知，一个孔的阻抗可籍矩阵法或串-并联法计算获得。本法中，多孔型电极的等价体是多孔的圆柱体。所有的孔都有相同的大小，它们之间平行排列。当获得一个孔的阻抗后，整个圆柱体或者整个多孔型电极的阻抗由公式(13)计算：

$Z_{\mathrm{total}}=\frac{Z_{\mathrm{pore}}}{N}---\left(13\right)$

公式中N等价圆柱体中孔结构数目。因为它们之间是并联排列，所以总阻抗是单个孔的N分之一。

6.1.矩阵法

如图6显示了电流在等价电路各单元的分布。

通过每个单元的电流(I₂，I₃，...，I_n)符合基耳霍夫原理(Kirchoff’s law)。公式(14)每个单元中R_s，Z_i和I_i的关系：

$\left[\begin{array}{c}V=(n{R}_s+Z_n)I_1-R_s(I_2+I_3+...+I_n)\\ 0=(Z_1+Z_2+R_s)I_2-R_s{I}_1-Z_2{I}_3\\ 0=(Z_2+Z_3+R_s)I_3-R_s{I}_1-Z_2{I}_2-Z_3{I}_4\\ 0=(Z_3+Z_4+R_s)I_4-R_s{I}_1-Z_3{I}_3-Z_4{I}_5\\ ....\\ 0=(Z_{n-3}+Z_{n-2}+R_s)I_{n-2}-R_s{I}_1-Z_{n-3}{I}_{n-3}-Z_{n-2}{I}_{n-1}\\ 0=(Z_{n-2}+Z_{n-1}+R_s)I_{n-1}-R_s{I}_1-Z_{n-2}{I}_{n-2}-Z_{n-1}{I}_n\\ 0=(Z_{n-1}+Z_n+R_s)I_n-R_s{I}_1-Z_{n-1}{I}_{n-1}-(Z_n+R_s)I_1\end{array}\right.---\left(14\right)$

公式(15)是对公式(14)的线性矩阵描述：

$\underset{n\×n}{\left(I_i\right)}\×\underset{n\×1}{(R_s/Z_i)}=\underset{n\×1}{\left(\begin{array}{c}{V}_1\\ 0\\ 0\\ 0\\ 0\end{array}\right)}---\left(15\right)$

解此矩阵，I₁到I_n可得，一个孔的阻抗Z_pore，可由此计算：

$Z_{\mathrm{pore}}=\frac{V_1}{I_1}---\left(16\right)$

整个电极的阻抗可由公式(13)计算获得。

这种方法需要较长的运算时间。但可以获得电流在各个单元的分布，这可用于通透百分数和通透深度的计算。

6.2.串-并联方法

这一方法使用等效电路中各个部件的串并联关系来计算阻抗。图6中，Z_n与R_s处于串联状态。这二者的总阻抗S_n＝Z_n+R_s，又与Z_n-1，处于并联状态。以上的总阻抗S_n-1＝1/(1/S_n+Z_n-1)，又与R_s处于串联状态......通过这样的过程，等价电路的总阻抗可计算获得。这个方法的优点在于计算时间短。

如果拟合的阻抗谱与实验值差别较大，可供选择的方案有：或者改变参数初始值的范围，或者采用不同的Z_w模式，或者使用不同形状的等价体等直到达到拟合精度。这些步骤可通过流程图表示，如图7就是数值迭代拟合法的流程图。根据此图编写了Mathematica6.0程序进行阻抗谱拟合。

这里对比了应用CNLLS和数值迭代法对三个电解二氧化锰电极，一个充满电的硫酸-铅电池的负极板的电化学阻抗谱的拟合。在两种方法中，不可逆-有限扩散类型作为Warburg阻抗Zw的类型。

实施例1

对三个电解二氧化锰电极的电化学阻抗谱的拟合

A.使用CNLLS法

如图8，三个EMD的Nyquist图，全频率域拟合值vs。试验值，应用不可逆-有限扩散Z_w。实心图标代表实验值；空心图标代表拟合值。CNLLS作为拟合计算方法；EMD1：■实验值，口拟合值；EMD2：▲实验值，△拟合值；EMD3：●实验值，○拟合值。

B.使用数值迭代法

如图9，三个EMD的Nyquist图，全频率域拟合值vs。试验值，应用不可逆-有限扩散Z_w。实心图标代表实验值；空心图标代表拟合值；数值迭代法作为拟合方法；EMD1：■实验值，口拟合值；EMD2：▲实验值，△拟合值；EMD3：●实验值，○拟合值。

实施例2

对一个一个充满电的硫酸-铅电池的负极板的电化学阻抗谱的拟合

A.使用CNLLS法

如图10，一个充满电的硫酸-铅电池负极的Nyquist图，全频率域拟合值vs。试验值，应用不可逆-有限扩散Z_w。CNLLS作为拟合计算方法；电池容量：7Ah；I＝0；◆实验值，■拟合值。

B。使用数值迭代法

如图11，一个充满电的硫酸-铅电池负极的Nyquist图，全频率域拟合值vs。试验值，应用不可逆-有限扩散Z_w。数值迭代法作为拟合计算方法；◆实验值，■拟合值。

实施例3

对一个充满电的硫酸-铅电池的电化学阻抗谱的拟合

A.使用CNLLS法

如图12，一个充满电的硫酸-铅电池的Nyquist图，全频率域拟合值vs。试验值，应用不可逆-有限扩散Z_w。CNLLS作为拟合计算方法；电池容量：7Ah；I＝0；◆实验值，■拟合值。

B.使用数值迭代法

如图13，一个充满电的硫酸-铅电池的Nyquist图，全频率域拟合值vs。试验值，应用不可逆-有限扩散Z_w。数值迭代法作为拟合计算方法；电池容量：7Ah；I＝0；◆实验值，■拟合值。

当对比这两种拟合方法，使用数值迭代法可以获得比使用CNLLS方法更好的拟合精度和拟合频率范围.

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种对于多孔性电极的电化学阻抗谱的数值迭代拟合方法，其特征在于，包括以下步骤：

步骤1)选择几何等价结构；

步骤2)选择Z_w模式，所述Z_w表示Warburg阻抗；

步骤3)为参数分配值域；

步骤4)均匀设计表；

步骤5)为参数赋初值；

步骤6)设i为角频率的序号，i初值为1，所述ω_i表示第i个角频率；

步骤7)在某一角频率ω_i下，使用数值迭代法计算Z₁到Z₁₀；

步骤8)通过Z₁，Z₂...Z₁₀，使用Mathematica中的内插值命令计算Z(x)；

步骤9)将Z(x)变换为每个Z单元的值，再使用矩阵法或串并联法得出Z_pore，所述Z_pore表示等价体中一个孔的阻抗；

步骤10)Z_total＝Z_pore/N，所述Z_total表示等价体中总阻抗，所述N表示等价体中包含的孔数；

步骤11)判断相位角是否接近其实验值，是，进入下一步，否，调整第一参数(C_d1)并返回步骤7，总是不，根据均匀设计表选择下一套参数初始值赋值方案，返回步骤5；

步骤12)判断实部是否接近其实验值，是，进入下一步，否，调整第二参数(D)并返回步骤7，总是不，根据均匀设计表选择下一套参数初始值赋值方案，返回步骤5；

步骤13)i＝i+1，返回步骤6。

2.根据权利要求1所述的对于多孔性电极的电化学阻抗谱的循环数值迭代拟合方法，其特征在于，由于激励电势沿孔壁从孔口到孔底的变化是孔深度x的函数，活性离子从本体溶液向孔壁扩散按照费克第二定律：解这个方程，边界条件是C_s(0)＝C_b和C_s(1_p)＝0，解出的C_s表达式是： $C_s=C_b\·\mathrm{Cosh}\left(\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}\·l_p}\right)\·\mathrm{Sinh}[\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}}\·(l_p-x)];$

公式中k电化学反应常数，V激励电势或过电势，D是扩散系数，n电化学反应转移电子数，R_ct是电荷转移电阻，r₀和l_p分别是等价体孔半径和孔深；费克第二定律方程先进行Laplace变换使其从时间域变成频率域然后解得C_s公式。

3.根据权利要求1所述的对于多孔性电极的电化学阻抗谱的循环数值迭代拟合方法，其特征在于，步骤7中所述数值迭代法计算Z_i包括以下步骤：

步骤a)将C_b作为C_s的起始值；

步骤b)通过

$I_F=\mathrm{nF}\·k\·C_s^{\mathrm{\γ}}$

${\mathrm{\σ}}_w=\frac{\mathrm{\γ}{R}_{\mathrm{ct}}\left|I_F\right|}{{\mathrm{nFC}}_s\sqrt{2{\mathrm{\ω}}_{i}D}};$

，计算得I_F及σ_w，所述C_b表示电解质溶液本体浓度，所述C_s表示电解质离子在电极表面的浓度，所述γ表示电化学反应级数，所述表示C_s的γ次方，所述I_F表示法拉第电流密度，所述F表示法拉第常数，所述k表示电化学反应常数，所述σ_w表示Warburg阻抗表达式前置因子，所述R_ct表示电荷转移电阻，所述ω_i表示第i个角频率，所述D表示扩散系数；

步骤c)计算得Z_w，所述Z_w表示Warburg阻抗；

步骤d)通过

$Z=\frac{1}{\frac{1}{R_{\mathrm{ct}}+Z_w}+j{\mathrm{\ω}}_i{C}_{d1}}$

计算得Z，式中Z为Z单元的阻抗，Z_w表示Warburg阻抗，R_ct表示电荷转移电阻，C_d1表示双电层电容，j表示虚数单位，ω_i为第i个角频率；

步骤e)通过

$V=V_0{e}^{(-x\·\sqrt{\frac{R_s}{Z}})}\·\sin \left({\mathrm{\ω}}_{i}t\right)$

计算得V，式中V表示所加正弦波电压，V₀表示V的波幅，R_s表示对应于在传输线模型等价电路中的电解质溶液电阻，Z为Z单元的阻抗，ω_i为第i个角频率；

步骤f)通过

$k=\frac{i_0[e^{-\mathrm{\αf\η}}-e^{(1-\mathrm{\α})\mathrm{f\η}}]}{\mathrm{nF}\·C_s^{\mathrm{\γ}}}$

计算得k，式中k表示电化学反应常数，i₀表示交换电流密度，α表示对称因子，F表示法拉第常数，η表示过电势，C_s表示电解质离子在电极表面的浓度，所述γ表示电化学反应级数，所述表示C_s的γ次方，f＝F/RT，F表示法拉第常数，R表示气体常数，T表示绝对温度；

步骤g)通过

$C_s=C_b\·\mathrm{Csch}(\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}}\·l_p)\·\mathrm{Sinh}[\sqrt{\frac{2k\·V}{D\·n\·r_0\·R_{\mathrm{ct}}}}\·(l_p-x)]$

计算得C_s，与之前所得C_s值进行比较，若前后两次C_s值之差小于10^-3则认为C_s已经达到稳定并进入步骤h，否则以所得C_s进入步骤b；

步骤h)计算得Z_w，所述Z_w表示Warburg阻抗；

步骤i)通过

$Z=\frac{1}{\frac{1}{R_{\mathrm{ct}}+Z_w}+j{\mathrm{\ω}}_i{C}_{d1}}$

计算得Z，式中Z为Z单元的阻抗，Z_w表示Warburg阻抗，R_ct表示电荷转移电阻，C_d1表示双电层电容，j表示虚数单位，ω_i为第i个角频率。

4.根据权利要求3所述的对于多孔性电极的电化学阻抗谱的数值迭代拟合方法，其特征在于：所述步骤c与步骤h中计算Z_w对于半无限扩散类型公式为：

Z_w＝σ_w·(1-j)

对于有线扩散类型公式为

$Z_w={\mathrm{\σ}}_w\·(\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}-j\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)})$

对于块层扩散类型公式为

$Z_w={\mathrm{\σ}}_w\·(\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}-j\frac{\sinh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)+\sin \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)}{\cosh \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)-\cos \left(\sqrt{\frac{2\mathrm{\ω}}{D}}\mathrm{\δ}\right)})$

式中Z_w表示Warburg阻抗，σ_w表示Warburg阻抗表达式前置因子，ω为ω_i，表示第i个角频率，D表示扩散系数，δ表示扩散层厚度，j表示虚数单位。

5.根据权利要求1所述的对于多孔性电极的电化学阻抗谱的数值迭代拟合方法，其特征在于：步骤9中对于包含由较大单元数目的传输线模型，通过电流在每个单元的分布遵守基耳霍夫原理，使用所述矩阵法计算出一个孔的电化学阻抗；对于包含由较大单元数目的传输线模型，通过等效电路中各个部件的串并联关系，使用所述串并联法来计算出一个孔的电化学阻抗。

6.根据权利要求1所述的对于多孔性电极的电化学阻抗谱的数值迭代拟合方法，其特征在于：步骤11中所述第一参数为C_d1，调整方式为当相位角大于其实验值，C_d1＝C_d1+C_d1/100；当相位角小于其实验值，C_d1＝C_d1-C_d1/100，式中C_d1表示双电层电容。

7.根据权利要求1所述的对于多孔性电极的电化学阻抗谱的数值迭代拟合方法，其特征在于：步骤12中所述第二参数为D，调整方式为当实部大于其实验值，D＝D+D/100；当实部小于其实验值，D＝D-D/100，式中D表示扩散系数。