CN104914312A - 一种计算交流阻抗谱弛豫时间分布的方法 - Google Patents

Abstract

一种计算交流阻抗谱弛豫时间分布的方法，本发明涉及交流阻抗谱弛豫时间分布的方法。本发明的目的是为了解决现有阻抗谱分析方法频率分辨率低、无法有效解析电化学反应过程的数量、实际阻抗以及无法求解弛豫时间分布的解析方程式的问题。通过以下技术方案实现的：一、获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；二、做KK检验，使交流阻抗谱数组是稳定的并且可以解析的；三、构建弛豫时间和弛豫时间分布函数的代数方程组；四、得到弛豫时间及弛豫时间分布函数数组，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。本发明应用于电化学领域。

Description

一种计算交流阻抗谱弛豫时间分布的方法

技术领域

本发明涉及交流阻抗谱弛豫时间分布的方法。

背景技术

交流阻抗谱技术已经成为电化学研究和表征的重要工具，广泛应用于诸多领域，例如表面防腐涂层耐蚀性能的研究、固体氧化物燃料电池电化学反应过程的解析以及氧还原机理的研究、膜反应器传输过程的研究、锂离子电池，超级电容器充放电性能研究。所谓交流阻抗即为将一个特定频率的交变电压扰动施加在一个电化学电池上，当系统稳定时会产生一个相同频率的交变电流输出，交变电压与电流的比值即为这一频率下的阻抗。由于电压与电流之间通常存在相位差，阻抗通常为复数。测试一定频率范围内的阻抗即形成阻抗谱。目前，阻抗谱的表征方法通常为Nyquist图(横坐标为阻抗实部，纵坐标为负的阻抗虚部)、Bode图(横坐标为频率、纵坐标为阻抗绝对值和相位角)、和实部虚部图(横坐标为频率、纵坐标为阻抗实部和负的阻抗虚部)。电化学电池的质量传输、电荷转移、表面交换、电荷传输等过程在一定频率的扰动下形成了该频率下的阻抗，即某一频率下的阻抗包含了电化学电池的所有过程的贡献。当电化学过程的特征频率相差较大时，在阻抗谱图上可表现为不同的弧(或者峰)。Nyquist图，Bode图和实部虚部图所表现的弧的个数是一致的，即他们的频率分辨率是一致的。通常情况下，当两个电化学过程的特征频率相差一个数量级时，才能表现出不同的弧，致使阻抗谱分析方法频率分辨率低。

为了区分阻抗谱上各个电化学过程对电池阻抗的贡献，目前通常是用线性化的机理模型拟合阻抗谱数据。这种方法必须要事先假设电化学过程的数量以及实际阻抗。例如，RQ等效电路通常表示电荷转移或表面交换过程，Warburg等效电路通常表示气体扩散或质量传输过程，Gerischer等效电路通常表示表面交换和离子传输的耦合过程。等效电路模型的选取依赖于具体电化学电池的构型、材料属性、工作条件等因素。为了较为全面地表征所有可能的电化学过程，等效电路模型通常为以上三种模型的组合，甚至它们的衍生模型(如传输线模型)。实际上很难确定最佳的等效电路模型，主要有以下原因：1)很多等效电路模型都能很好地拟合阻抗谱数据；2)电化学过程的阻抗相应在频率空间内有重叠，即在特定角频率ω下的阻抗包含了特征时间在ω^-1上下的所有电化学过程的贡献。所以，即便选用了某一等效电路模型，也无法检验假设的合理性。这给电化学机理的研究带来主观性和不确定性，致使无法有效解析电化学反应过程的数量以及实际阻抗的问题。所以，发展具有高频率分辨率的阻抗谱分析方法具有重要意义。

根据德拜弛豫时间公式，任何满足Kramers-Kronig关系的阻抗谱都可以表示为弛豫时间τ的积分函数：

$Z\left(\mathrm{\ω}\right)-Z^{\′}(\∞)={\∫}_0^{\∞}\frac{F\left(\mathrm{\τ}\right)}{1+\mathrm{j\ω\τ}}\mathrm{d\; l}{\mathrm{og}}_{10}\mathrm{\τ}$

其中，F(τ)表示阻抗Z的弛豫时间分布函数；τ是弛豫时间；Z'(∞)表示阻抗实部在角频率ω趋于无限大时的极限。以-Log₁₀2πτ为横轴，F(τ)为纵轴作图，就可以将电化学过程最大限度地区分开来，因为某一特定弛豫时间对应的弛豫时间分布只代表该弛豫时间过程的贡献，去除了其他弛豫时间过程的影响。所以，相比于现有的阻抗谱分析方法，弛豫时间分布具有最高的频率分辨率。1941年，Fuoss和Kirkwood在J.Am.Chem.Soc.杂志上首次发表了求解弛豫时间分布的解析方程式，但需要知道阻抗Z与角频率ω的解析关系。实际上，阻抗Z与角频率ω的解析关系是未知的，也正是电化学工作者力图知道的，致使无法求解弛豫时间分布的解析方程式。

发明内容

本发明的目的是为了解决现有阻抗谱分析方法频率分辨率低、无法有效解析电化学反应过程的数量、实际阻抗以及无法求解弛豫时间分布的解析方程式的问题，而提出了一种计算交流阻抗谱弛豫时间分布的方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；

步骤二、对阻抗实部和阻抗虚部做Kramers-Kronig检验，使步骤一得到的交流阻抗谱数组是稳定的并且可以解析的，其中，所述Kramers-Kronig检验为实部和虚部的检验；

步骤三、在确定步骤一得到的交流阻抗谱是稳定的并且可以解析的基础上，根据频率和阻抗虚部构建弛豫时间和弛豫时间分布函数的代数方程组；

步骤四、应用Tikhonov正则化方法和二次规划方法求解弛豫时间和弛豫时间分布函数的代数方程组，得到弛豫时间及弛豫时间分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，所述图的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗，其中，所述Tikhonov正则化方法为吉洪诺夫正则化方法。

发明效果

采用本发明的一种计算交流阻抗谱弛豫时间分布的方法，首先获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；然后对阻抗谱实部和虚部做Kramers-Kronig检验，保证阻抗谱数据是稳定的并且可以解析的；再根据频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组；最后应用Tikhonov正则化方法和二次规划方法求解代数方程组，得到弛豫时间及其分布函数数组，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应特定电化学过程，峰面积代表特定电化学过程的实际阻抗。相比于现有的阻抗谱分析手段，如Nyquist、Bode等方法，本发明的一种计算阻抗谱弛豫时间分布的方法提高了阻抗谱分析的频率分辨率，可以在一个频率数量级内解析出两到三个电化学过程，可区分弛豫时间严重重叠的电化学过程，并且无需经验假设就可以直接解析出电化学过程的个数和实际阻抗，求解了弛豫时间分布的解析方程式。

附图说明

图1为本发明流程图；

图2为三个串联RC电路的模拟阻抗的Nyquist图，阻抗实部和阻抗虚部单位为Ωcm²；

图3为该电路的模拟阻抗的虚部图，频率单位为Hz；负的阻抗虚部单位为Ωcm²；

图4为该电路的模拟阻抗的弛豫时间分布图，弛豫时间分布函数单位为Ωcm²；

图5为Ni-YSZ//YSZ//LSM-YSZ固体氧化物燃料电池在800℃下的阻抗谱虚部图，阳极燃料为1H₂:4N₂(3vol.％H₂O)，阴极气氛为空气，频率单位为Hz；负的阻抗虚部单位为Ωcm²；

图6为该燃料电池阻抗谱的Kramers-Kronig检验结果图，纵坐标为相对误差，×为阻抗实部的相对误差，100*(Z'_KK-Z'+﹤R₀﹥)/R_P为阻抗实部相对误差的定义，+为阻抗虚部的相对误差，100*(Z”_KK-Z”)/R_P为阻抗虚部相对误差的定义；

图7为该燃料电池阻抗谱的弛豫时间分布图。

具体实施方式

具体实施方式一：结合图1说明本实施方式，一种计算交流阻抗谱弛豫时间分布的方法具体是按照以下步骤进行的：

一种计算交流阻抗谱弛豫时间分布的方法具体是按照以下步骤进行的：

步骤一、获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；

步骤二、对阻抗实部和阻抗虚部做Kramers-Kronig检验，使步骤一得到的交流阻抗谱数组是稳定的并且可以解析的，其中，所述Kramers-Kronig检验为实部和虚部的检验；

步骤三、在确定步骤一得到的交流阻抗谱是稳定的并且可以解析的基础上，根据频率和阻抗虚部构建弛豫时间和弛豫时间分布函数的代数方程组；

步骤四、应用Tikhonov正则化方法和二次规划方法求解弛豫时间和弛豫时间分布函数的代数方程组，得到弛豫时间及弛豫时间分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，所述图的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗，其中，所述Tikhonov正则化方法为吉洪诺夫正则化方法。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；具体过程为：

交流阻抗谱数组由电化学工作站或模拟等效电路测得；

(1)模拟等效电路测得交流阻抗谱数组的过程为：

模拟等效电路选a个RC电路的串联电路，a取值范围为任意正整数，其中R为RC电路的电阻；R单位为Ωcm²，为任意正数值；C为RC电路的电容；C单位为F/cm²，为任意正数值；其中，所述RC电路为相移电路，Ωcm²为欧姆/平方厘米，F/cm²为法/平方厘米；

交流阻抗的频率范围为10⁷Hz～10^-4Hz，每频率数量级取x个离散的频率数据和阻抗数值，x取值为10到100间的整数，取值越高阻抗谱数组越多，其中，所述每频率数量级为在频率之比为10的两个频率，频率单位为Hz，得出该等效电路的模拟阻抗谱的Nyquist图和虚部图，获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部，其中，所述Hz为赫兹；

(2)电化学工作站测得交流阻抗谱数组的过程为：

电化学电池选为固体氧化物燃料电池，测试温度为500～1000℃，阳极气氛为氢气、合成气或碳氢化合物，阴极气氛为空气或氧气；

交流阻抗的频率范围为10⁶Hz～10^-2Hz，每频率数量级取x个离散的频率数据和阻抗数值，x取值为10到100间的整数，取值越高阻抗谱数组越多，其中，所述每频率数量级为在频率之比为10的两个频率，频率单位为Hz，例如测试频率范围10⁶～10^-2Hz跨越8个数量级，总共取80个离散数据点，得出电池阻抗谱的虚部图，获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中对阻抗实部和阻抗虚部做Kramers-Kronig检验，使步骤一得到的交流阻抗谱数组是稳定的并且可以解析的，其中，所述Kramers-Kronig检验为实部和虚部的检验；具体过程为：

Kramers-Kronig检验用ZSimpWin软件执行或用自编计算机程序执行；

(1)ZSimpWin软件执行Kramers-Kronig检验的具体操作步骤为：

1)把交流阻抗谱数组导入到ZSimpWin软件；

2)在工具栏窗口点击“extrapolate”(推断)下的“Apply Kramers-Kronig”(yingyong实部和虚部的检验)按钮，得到计算出的交流阻抗谱；

3)对比计算出的交流阻抗谱和步骤一得到的交流阻抗谱，若计算出的交流阻抗谱在低频区与步骤一得到的交流阻抗谱平滑连接，则步骤一得到的交流阻抗谱是稳定的并且可以解析的，其中所述低频为频率小于1Hz；

(2)自编程序执行Kramers-Kronig检验的具体操作步骤为：

1)根据阻抗实部和阻抗虚部的Kramers-Kronig关系：

$Z^{\′\′}{\left(\mathrm{\ω}\right)}_{\mathrm{cal}}=-\left(\frac{2\mathrm{\ω}}{\mathrm{\π}}\right){\∫}_0^{\∞}\frac{Z^{\′}\left(x\right)-Z^{\′}\left(\mathrm{\ω}\right)}{x^2-{\mathrm{\ω}}^2}\mathrm{dx};$

$Z^{\′}{\left(\mathrm{\ω}\right)}_{\mathrm{cal}}=Z^{\′}(\∞)+\frac{2}{\mathrm{\π}}{\∫}_0^{\∞}\frac{x{Z}^{\′\′}\left(x\right)-{\mathrm{\ωZ}}^{\′\′}\left(\mathrm{\ω}\right)}{x^2-{\mathrm{\ω}}^2}\mathrm{dx};$

由步骤一得到的阻抗实部Z'(ω)计算出阻抗虚部Z”(ω)_cal，根据步骤一得到的阻抗虚部Z”(ω)计算出阻抗实部Z'(ω)_cal；

2)对比计算出的阻抗实部Z'(ω)_cal和阻抗虚部Z”(ω)_cal与步骤一得到的阻抗实部Z'(ω)和阻抗虚部Z”(ω)，若计算出的阻抗实部Z'(ω)_cal与步骤一得到的阻抗实部Z'(ω)的相对误差和计算出的阻抗虚部Z”(ω)_cal与步骤一得到的阻抗虚部Z”(ω)的相对误差都在10％以内，则步骤一得到的交流阻抗谱是稳定的并且可以解析的。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一、二或三不同的是：所述步骤三中在确定步骤一得到的交流阻抗谱是稳定的并且可以解析的基础上，根据频率和阻抗虚部构建弛豫时间和弛豫时间分布函数的代数方程组；具体过程为：

弛豫时间τ与频率ω的函数关系为ωτ＝1；

弛豫时间分布函数F(τ)与步骤一得到的阻抗虚部Z”(ω)、频率ω和弛豫时间τ的函数关系为：

$-Z^{\′\′}\left(\mathrm{\ω}\right)={\∫}_0^{\∞}\frac{\mathrm{\ωF}\left(\mathrm{\τ}\right)/\ln \left(10\right)}{1+{\left(\mathrm{\ω\τ}\right)}^2}\mathrm{d\τ};$

式中，τ为弛豫时间，单位为s；ω为频率，单位为rad/s；F(τ)为弛豫时间分布函数，单位为Ωcm²；Z”(ω)为步骤一得到的阻抗虚部，单位为Ωcm²；

由于阻抗谱数组是离散形式的，所以将离散化，求得第n个弛豫时间τ_n下的待求解的弛豫时间分布矩阵第n个元素F(τ_n)；

离散化的弛豫时间和弛豫时间分布函数的代数方程组为：ΓF＝Z

式中，Γ和Z是已知的矩阵，F为待求解的弛豫时间分布矩阵，为列矩阵，其第n个元素为F(τ_n)；

Γ是N+1行N列矩阵，N为阻抗谱数据点个数；

当m＝1,2,…,N，n＝1,2,…,N时，m为矩阵元素行坐标，n为矩阵元素列坐标，

当n＝1,2,…,N时，Γ_N+1,n＝-δw；

F是N行列矩阵；当n＝1,2,…,N时，F_n为第n个待求解的弛豫时间分布函数；

Z是N+1行列矩阵；当n＝1,2,…,N时， $Z_n=-Z^{\′\′}\left({10}^{W_n}\right),Z_{N+1}=|Z^{\′}\left({10}^{W_N}\right)-Z^{\′}\left({10}^{W_1}\right)|;$ Z'为阻抗谱实部，W_m＝log₁₀ω_m，ω_m为第m个频率；w_n＝-log₁₀τ_n，τ_n为第n个弛豫时间。

其它步骤及参数与具体实施方式一、二或三相同。

具体实施方式五：本实施方式与具体实施方式一、二、三或四不同的是：所述步骤四中应用Tikhonov正则化方法和二次规划方法求解弛豫时间和弛豫时间分布函数的代数方程组，得到弛豫时间及弛豫时间分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，所述图的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗，其中，所述Tikhonov正则化方法为吉洪诺夫正则化方法；具体过程为：

采用Tikhonov正则化方法、二次规划方法和MATLAB软件的内置函数quadprog，通过目标函数最小化err求解弛豫时间和弛豫时间分布函数的代数方程组：

err＝(ΓF-Z)^T(ΓF-Z)+λ(DF)^T(DF)

并且满足弛豫时间分布函数非负的限制条件：

F≥0

式中，上标T表示矩阵的转置操作；D是N维Tikhonov正则化方阵；当n＝2,3,…,N-1时，D_{n,[n-1,n,n+1]}＝{-1,2,-1}；D_1,[1,2]＝{1,-1}；D_N,[N-1,N]＝{-1,1}，其他元素为0；λ为Tikhonov正则化系数，Tikhonov正则化系数λ为10^-4～0；得到弛豫时间及弛豫时间分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。

其它步骤及参数与具体实施方式一、二、三或四相同。

采用以下实施例验证本发明的有益效果：

实施例1

一种计算交流阻抗谱弛豫时间分布的方法具体是按照以下步骤进行的：

1)、通过模拟电路获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部。该实施方式的模拟电路选为三个RC电路的串联电路，其中R₁＝R₂＝R₃＝1Ωcm²；C₁＝1e-3F/cm²；C₂＝5e-5F/cm²；C₃＝1e-5F/cm²。频率范围是1e-2～1e6Hz，每频率数量级取10个阻抗数据。图2给出了该等效电路的模拟阻抗谱的Nyquist图。图3给出了相应的虚部图。由图2和图3可见，只有两个弧呈现出来。根据Nyquist和虚部图并不能确定该阻抗谱是由三个RC等效过程产生的。

2)、根据频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组。弛豫时间τ[s]与频率ω[rad/s]的函数关系为ωτ＝1，弛豫时间分布函数F(τ)[Ωcm²]与阻抗谱虚部Z”(ω)[Ωcm²]、频率ω和弛豫时间τ的函数关系为:

$-Z^{\′\′}\left(\mathrm{\ω}\right)={\∫}_0^{\∞}\frac{\mathrm{\ωF}\left(\mathrm{\τ}\right)/\ln \left(10\right)}{1+{\left(\mathrm{\ω\τ}\right)}^2}\mathrm{d\τ}.$

根据频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组为：

ΓF＝Z

其中、Γ和Z是已知的矩阵，F为待求解的弛豫时间分布矩阵；Γ是N+1行N列矩阵，N为阻抗谱数据点个数；当m＝1,2,…,N，n＝1,2,…,N时，当n＝1,2,…,N时，Γ_N+1,n＝-δw；F是N行列矩阵；当n＝1,2,…,N时，Z是N+1行列矩阵；当n＝1,2,…,N时， $Z_n=-Z^{\′\′}\left({10}^{W_n}\right);Z_{N+1}=|Z^{\′}\left({10}^{W_N}\right)-Z^{\′}\left({10}^{W_1}\right)|;$ Z'为阻抗谱实部；W_n＝log₁₀ω_n；w_n＝-log₁₀τ_n。

3)、应用Tikhonov正则化方法和二次规划方法求解代数方程组，得到弛豫时间及其分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。应用Tikhonov正则化方法和二次规划方法求解代数方程组，是通过最小化下面的目标函数求解的：

err＝(ΓF-Z)^T(ΓF-Z)+λ(DF)^T(DF)

并且满足弛豫时间分布函数非负的限制条件：

F≥0

其中、上标T表示矩阵的转置操作；D是N维Tikhonov正则化方阵；当n＝2,3,…,N-1时，D_{n,[n-1,n,n+1]}＝{-1,2,-1}；D_1,[1,2]＝{1,-1}；D_N,[N-1,N]＝{-1,1}，其他元素为0；λ为Tikhonov正则化系数，λ取值为10^-4。应用MATLAB软件的内置函数quadprog求解在弛豫时间分布函数非负条件下目标函数err最小化的二次规划问题。图4给出了该等效电路的弛豫时间分布的计算结果。由图4可见三个峰，这三个峰所围的面积分别为1.0086Ωcm²、0.9926Ωcm²和1.0072Ωcm²。与这三个RC电路的阻抗1Ωcm²基本一致。峰顶对应的弛豫时间分别为9.568e-4秒、5.27e-5秒和1.10e-5秒，与这三个RC电路的特征时间(RC)1e-3秒、5e-5秒和1e-5秒基本一致。

实施例2

一种计算交流阻抗谱弛豫时间分布的方法具体是按照以下步骤进行的：

1)、应用电化学工作站则是实际电化学电池的交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部。该实施方式的电化学电池选为Ni-YSZ//YSZ//LSM-YSZ固体氧化物燃料电池，测试温度为800℃，阳极气氛为1H₂:4N₂(3vol.％H₂O)，阴极气氛为空气，频率范围是1e-2～1e5Hz，每频率数量级取10个阻抗数据。图5给出了该电池阻抗谱的虚部图，可以看见两个比较宽泛的峰，不能确定电化学反应过程的数量。

2)、应用ZSimpWin软件，对阻抗谱实部和虚部做Kramers-Kronig检验，保证阻抗谱数据是稳定的并且可以解析的。

应用ZSimpWin软件执行Kramers-Kronig检验的具体操作步骤是：

(1)ZSimpWin软件执行Kramers-Kronig检验的具体操作步骤为：

1)把交流阻抗谱数组导入到ZSimpWin软件；

2)在工具栏窗口点击“extrapolate”(推断)下的“Apply Kramers-Kronig”按钮，得到计算出的交流阻抗谱；

3)对比计算出的交流阻抗谱和步骤一得到的交流阻抗谱，若计算出的交流阻抗谱在低频区与步骤一得到的交流阻抗谱平滑连接，则步骤一得到的交流阻抗谱是稳定的并且可以解析的，其中所述低频为频率小于1Hz；

(2)自编程序执行Kramers-Kronig检验的具体操作步骤为：

1)根据阻抗实部和阻抗虚部的Kramers-Kronig关系：

$Z^{\′\′}{\left(\mathrm{\ω}\right)}_{\mathrm{cal}}=-\left(\frac{2\mathrm{\ω}}{\mathrm{\π}}\right){\∫}_0^{\∞}\frac{Z^{\′}\left(x\right)-Z^{\′}\left(\mathrm{\ω}\right)}{x^2-{\mathrm{\ω}}^2}\mathrm{dx};$

$Z^{\′}{\left(\mathrm{\ω}\right)}_{\mathrm{cal}}=Z^{\′}(\∞)+\frac{2}{\mathrm{\π}}{\∫}_0^{\∞}\frac{x{Z}^{\′\′}\left(x\right)-\mathrm{\ω}{Z}^{\′\′}\left(\mathrm{\ω}\right)}{x^2-{\mathrm{\ω}}^2}\mathrm{dx};$

由步骤一得到的阻抗实部Z'(ω)计算出阻抗虚部Z”(ω)_cal，根据步骤一得到的阻抗虚部Z”(ω)计算出阻抗实部Z'(ω)_cal；

2)对比计算出的阻抗实部Z'(ω)_cal和阻抗虚部Z”(ω)_cal与步骤一得到的阻抗实部Z'(ω)和阻抗虚部Z”(ω)，若计算出的阻抗实部Z'(ω)_cal与步骤一得到的阻抗实部Z'(ω)的相对误差和计算出的阻抗虚部Z”(ω)_cal与步骤一得到的阻抗虚部Z”(ω)的相对误差都在10％以内，则步骤一得到的交流阻抗谱是稳定的并且可以解析的。

图6给出了该电池阻抗谱的Kramers-Kronig检验结果。可见阻抗实部和阻抗虚部的百分误差只有几个百分点，可以利用该阻抗谱数据重构阻抗谱弛豫时间分布。

3)、根据频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组。弛豫时间τ[s]与频率ω[rad/s]的函数关系为ωτ＝1，弛豫时间分布函数F(τ)[Ωcm²]与阻抗谱虚部Z”(ω)[Ωcm²]、频率ω和弛豫时间τ的函数关系为:

$-Z^{\′\′}\left(\mathrm{\ω}\right)={\∫}_0^{\∞}\frac{\mathrm{\ωF}\left(\mathrm{\τ}\right)/\ln \left(10\right)}{1+{\left(\mathrm{\ω\τ}\right)}^2}\mathrm{d\τ}.$

由频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组为：

ΓF＝Z

其中、Γ和Z是已知的矩阵，F为待求解的弛豫时间分布矩阵；Γ是N+1行N列矩阵，N为阻抗谱数据点个数；当m＝1,2,…,N，n＝1,2,…,N时，当n＝1,2,…,N时，Γ_N+1,n＝-δw；F是N行列矩阵；当n＝1,2,…,N时，Z是N+1行列矩阵；当n＝1,2,…,N时， $Z_n=-Z^{\′\′}\left({10}^{W_n}\right);Z_{N+1}=|Z^{\′}\left({10}^{W_N}\right)-Z^{\′}\left({10}^{W_1}\right)|;$ Z'为阻抗谱实部；W_m＝log₁₀ω_m，ω_m为第m个频率；w_n＝-log₁₀τ_n，τ_n为第n个弛豫时间。

4)、应用Tikhonov正则化方法和二次规划方法求解代数方程组，得到弛豫时间及其分布函数数组，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。应用Tikhonov正则化方法和二次规划方法求解代数方程组，是通过最小化下面的目标函数求解的：

err＝(ΓF-Z)^T(ΓF-Z)+λ(DF)^T(DF)

并且满足弛豫时间分布函数非负的限制条件：

F≥0

其中、上标T表示矩阵的转置操作；D是N维Tikhonov正则化方阵；当n＝2,3,…,N-1时，D_{n,[n-1,n,n+1]}＝{-1,2,-1}；D_1,[1,2]＝{1,-1}；D_N,[N-1,N]＝{-1,1}，其他元素为0；λ为Tikhonov正则化系数，λ取值为10^-4。应用MATLAB软件的内置函数quadprog求解在弛豫时间分布函数非负条件下目标函数err最小化的二次规划问题。图7给出了该电池阻抗谱的弛豫时间分布图，可以看出有6个峰，即有6个电化学过程。而这些过程在虚部图(图5)中是无法确定出来的。

Claims (5)

1.一种计算交流阻抗谱弛豫时间分布的方法，其特征在于，一种计算交流阻抗谱弛豫时间分布的方法具体是按照以下步骤进行的：

步骤一、获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；

步骤二、对阻抗实部和阻抗虚部做Kramers-Kronig检验，使步骤一得到的交流阻抗谱数组是稳定的并且可以解析的，其中，所述Kramers-Kronig检验为实部和虚部的检验；

步骤三、在确定步骤一得到的交流阻抗谱是稳定的并且可以解析的基础上，根据频率和阻抗虚部构建弛豫时间和弛豫时间分布函数的代数方程组；

步骤四、应用Tikhonov正则化方法和二次规划方法求解弛豫时间和弛豫时间分布函数的代数方程组，得到弛豫时间及弛豫时间分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，所述图的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗，其中，所述Tikhonov正则化方法为吉洪诺夫正则化方法。

2.根据权利要求1所述一种计算交流阻抗谱弛豫时间分布的方法，其特征在于，所述步骤一中获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；具体过程为：

交流阻抗谱数组由电化学工作站或模拟等效电路测得；

(1)模拟等效电路测得交流阻抗谱数组的过程为：

模拟等效电路选a个RC电路的串联电路，a取值范围为任意正整数，其中R为RC电路的电阻；R单位为Ωcm²，为任意正数值；C为RC电路的电容；C单位为F/cm²，为任意正数值；其中，所述RC电路为相移电路，Ωcm²为欧姆/平方厘米，F/cm²为法/平方厘米；

交流阻抗的频率范围为10⁷Hz～10^-4Hz，每频率数量级取x个离散的频率数据和阻抗数值，x取值为10到100间的整数，取值越高阻抗谱数组越多，其中，所述每频率数量级为在频率之比为10的两个频率，频率单位为Hz，得出该等效电路的模拟阻抗谱的Nyquist图和虚部图，获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部，其中，所述Hz为赫兹；

(2)电化学工作站测得交流阻抗谱数组的过程为：

电化学电池选为固体氧化物燃料电池，测试温度为500～1000℃，阳极气氛为氢气、合成气或碳氢化合物，阴极气氛为空气或氧气；

交流阻抗的频率范围为10⁶Hz～10^-2Hz，每频率数量级取x个离散的频率数据和阻抗数值，x取值为10到100间的整数，取值越高阻抗谱数组越多，其中，所述每频率数量级为在频率之比为10的两个频率，频率单位为Hz，得出电池阻抗谱的虚部图，获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部。

3.根据权利要求2所述一种计算交流阻抗谱弛豫时间分布的方法，其特征在于，所述步骤二中对阻抗实部和阻抗虚部做Kramers-Kronig检验，使步骤一得到的交流阻抗谱数组是稳定的并且可以解析的，其中，所述Kramers-Kronig检验为实部和虚部的检验；具体过程为：

Kramers-Kronig检验用ZSimpWin软件执行或用自编计算机程序执行；

(1)ZSimpWin软件执行Kramers-Kronig检验的具体操作步骤为：

1)把交流阻抗谱数组导入到ZSimpWin软件；

2)在工具栏窗口点击“extrapolate”下的“Apply Kramers-Kronig”按钮，得到计算出的交流阻抗谱；

3)对比计算出的交流阻抗谱和步骤一得到的交流阻抗谱，若计算出的交流阻抗谱在低频区与步骤一得到的交流阻抗谱平滑连接，则步骤一得到的交流阻抗谱是稳定的并且可以解析的，其中所述低频为频率小于1Hz；

(2)自编程序执行Kramers-Kronig检验的具体操作步骤为：

1)根据阻抗实部和阻抗虚部的Kramers-Kronig关系：

$Z^{\′\′}{\left(\mathrm{\ω}\right)}_{\mathrm{cal}}=-\left(\frac{2\mathrm{\ω}}{\mathrm{\π}}\right){\∫}_0^{\∞}\frac{Z^{\′}\left(x\right)-Z^{\′}\left(\mathrm{\ω}\right)}{x^2-{\mathrm{\ω}}^2}\mathrm{dx};$

$Z^{\′}{\left(\mathrm{\ω}\right)}_{\mathrm{cal}}=Z^{\′}(\∞)+\frac{2}{\mathrm{\π}}{\∫}_0^{\∞}\frac{x{Z}^{\′\′}\left(x\right)-\mathrm{\ω}{Z}^{\′\′}\left(\mathrm{\ω}\right)}{x^2-{\mathrm{\ω}}^2}\mathrm{dx};$

由步骤一得到的阻抗实部Z'(ω)计算出阻抗虚部Z”(ω)_cal，根据步骤一得到的阻抗虚部Z”(ω)计算出阻抗实部Z'(ω)_cal；

2)对比计算出的阻抗实部Z'(ω)_cal和阻抗虚部Z”(ω)_cal与步骤一得到的阻抗实部Z'(ω)和阻抗虚部Z”(ω)，若计算出的阻抗实部Z'(ω)_cal与步骤一得到的阻抗实部Z'(ω)的相对误差和计算出的阻抗虚部Z”(ω)_cal与步骤一得到的阻抗虚部Z”(ω)的相对误差都在10％以内，则步骤一得到的交流阻抗谱是稳定的并且可以解析的。

4.根据权利要求3所述一种计算交流阻抗谱弛豫时间分布的方法，其特征在于，所述步骤三中在确定步骤一得到的交流阻抗谱是稳定的并且可以解析的基础上，根据频率和阻抗虚部构建弛豫时间和弛豫时间分布函数的代数方程组；具体过程为：

弛豫时间τ与频率ω的函数关系为ωτ＝1；

弛豫时间分布函数F(τ)与步骤一得到的阻抗虚部Z”(ω)、频率ω和弛豫时间τ的函数关系为：

${-Z}^{\′\′}\mathrm{\ω}={\∫}_0^{\∞}\frac{\mathrm{\ωF}\left(\mathrm{\τ}\right)/\ln \left(10\right)}{1+{\left(\mathrm{\ω\τ}\right)}^2}\mathrm{d\τ};$

式中，τ为弛豫时间，单位为s；ω为频率，单位为rad/s；F(τ)为弛豫时间分布函数，单位为Ωcm²；Z”(ω)为步骤一得到的阻抗虚部，单位为Ωcm²；

由于阻抗谱数组是离散形式的，所以将离散化，求得第n个弛豫时间τ_n下的待求解的弛豫时间分布矩阵第n个元素F(τ_n)；

离散化的弛豫时间和弛豫时间分布函数的代数方程组为：ΓF＝Z

式中，Γ和Z是已知的矩阵，F为待求解的弛豫时间分布矩阵，为列矩阵，其第n个元素为F(τ_n)；

Γ是N+1行N列矩阵，N为阻抗谱数据点个数；

当m＝1,2,…,N，n＝1,2,…,N时，m为矩阵元素行坐标，n为矩阵元素列坐标，

当n＝1,2,…,N时，Γ_N+1,n＝-δw；

F是N行列矩阵；当n＝1,2,…,N时，F_n为第n个待求解的弛豫时间分布函数；

Z是N+1行列矩阵；当n＝1,2,…,N时， $Z_n=-Z^{\′\′}\left({10}^{W_n}\right),$ $Z_{N+1}=|Z^{\′}\left({10}^{W_N}\right)-Z^{\′}\left({10}^{W_1}\right)|;$ Z'为阻抗谱实部，W_m＝log₁₀ω_m，ω_m为第m个频率；w_n＝-log₁₀τ_n，τ_n为第n个弛豫时间。

5.根据权利要求4所述一种计算交流阻抗谱弛豫时间分布的方法，其特征在于，所述步骤四中应用Tikhonov正则化方法和二次规划方法求解弛豫时间和弛豫时间分布函数的代数方程组，得到弛豫时间及弛豫时间分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，所述图的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗，其中，所述Tikhonov正则化方法为吉洪诺夫正则化方法；具体过程为：

采用Tikhonov正则化方法、二次规划方法和MATLAB软件的内置函数quadprog，通过目标函数最小化err求解弛豫时间和弛豫时间分布函数的代数方程组：

err＝(ΓF-Z)^T(ΓF-Z)+λ(DF)^T(DF)

并且满足弛豫时间分布函数非负的限制条件：

F≥0

式中，上标T表示矩阵的转置操作；D是N维Tikhonov正则化方阵；当n＝2,3,…,N-1时，D_{n,[n-1,n,n+1]}＝{-1,2,-1}；D_1,[1,2]＝{1,-1}；D_N,[N-1,N]＝{-1,1}，其他元素为0；λ为Tikhonov正则化系数，Tikhonov正则化系数λ为10^-4～0；得到弛豫时间及弛豫时间分布函数数组{τ_n，F(τ_n)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。