CN104914312B - 一种计算交流阻抗谱弛豫时间分布的方法 - Google Patents
一种计算交流阻抗谱弛豫时间分布的方法 Download PDFInfo
- Publication number
- CN104914312B CN104914312B CN201510341876.7A CN201510341876A CN104914312B CN 104914312 B CN104914312 B CN 104914312B CN 201510341876 A CN201510341876 A CN 201510341876A CN 104914312 B CN104914312 B CN 104914312B
- Authority
- CN
- China
- Prior art keywords
- impedance
- relaxation time
- mrow
- msup
- frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Abstract
一种计算交流阻抗谱弛豫时间分布的方法，本发明涉及交流阻抗谱弛豫时间分布的方法。本发明的目的是为了解决现有阻抗谱分析方法频率分辨率低、无法有效解析电化学反应过程的数量、实际阻抗以及无法求解弛豫时间分布的解析方程式的问题。通过以下技术方案实现的：一、获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；二、做KK检验，使交流阻抗谱数组是稳定的并且可以解析的；三、构建弛豫时间和弛豫时间分布函数的代数方程组；四、得到弛豫时间及弛豫时间分布函数数组，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。本发明应用于电化学领域。
Description
技术领域
本发明涉及交流阻抗谱弛豫时间分布的方法。
背景技术
交流阻抗谱技术已经成为电化学研究和表征的重要工具，广泛应用于诸多领域，例如表面防腐涂层耐蚀性能的研究、固体氧化物燃料电池电化学反应过程的解析以及氧还原机理的研究、膜反应器传输过程的研究、锂离子电池，超级电容器充放电性能研究。所谓交流阻抗即为将一个特定频率的交变电压扰动施加在一个电化学电池上，当系统稳定时会产生一个相同频率的交变电流输出，交变电压与电流的比值即为这一频率下的阻抗。由于电压与电流之间通常存在相位差，阻抗通常为复数。测试一定频率范围内的阻抗即形成阻抗谱。目前，阻抗谱的表征方法通常为Nyquist图(横坐标为阻抗实部，纵坐标为负的阻抗虚部)、Bode图(横坐标为频率、纵坐标为阻抗绝对值和相位角)、和实部虚部图(横坐标为频率、纵坐标为阻抗实部和负的阻抗虚部)。电化学电池的质量传输、电荷转移、表面交换、电荷传输等过程在一定频率的扰动下形成了该频率下的阻抗，即某一频率下的阻抗包含了电化学电池的所有过程的贡献。当电化学过程的特征频率相差较大时，在阻抗谱图上可表现为不同的弧(或者峰)。Nyquist图，Bode图和实部虚部图所表现的弧的个数是一致的，即他们的频率分辨率是一致的。通常情况下，当两个电化学过程的特征频率相差一个数量级时，才能表现出不同的弧，致使阻抗谱分析方法频率分辨率低。
为了区分阻抗谱上各个电化学过程对电池阻抗的贡献，目前通常是用线性化的机理模型拟合阻抗谱数据。这种方法必须要事先假设电化学过程的数量以及实际阻抗。例如，RQ等效电路通常表示电荷转移或表面交换过程，Warburg等效电路通常表示气体扩散或质量传输过程，Gerischer等效电路通常表示表面交换和离子传输的耦合过程。等效电路模型的选取依赖于具体电化学电池的构型、材料属性、工作条件等因素。为了较为全面地表征所有可能的电化学过程，等效电路模型通常为以上三种模型的组合，甚至它们的衍生模型(如传输线模型)。实际上很难确定最佳的等效电路模型，主要有以下原因：1)很多等效电路模型都能很好地拟合阻抗谱数据；2)电化学过程的阻抗相应在频率空间内有重叠，即在特定角频率ω下的阻抗包含了特征时间在ω-1上下的所有电化学过程的贡献。所以，即便选用了某一等效电路模型，也无法检验假设的合理性。这给电化学机理的研究带来主观性和不确定性，致使无法有效解析电化学反应过程的数量以及实际阻抗的问题。所以，发展具有高频率分辨率的阻抗谱分析方法具有重要意义。
根据德拜弛豫时间公式，任何满足Kramers-Kronig关系的阻抗谱都可以表示为弛豫时间τ的积分函数：
其中，F(τ)表示阻抗Z的弛豫时间分布函数；τ是弛豫时间；Z'(∞)表示阻抗实部在角频率ω趋于无限大时的极限。以-Log102πτ为横轴，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图，阻抗实部和阻抗虚部单位为Ωcm2；
图3为该电路的模拟阻抗的虚部图，频率单位为Hz；负的阻抗虚部单位为Ωcm2；
图4为该电路的模拟阻抗的弛豫时间分布图，弛豫时间分布函数单位为Ωcm2；
图5为Ni-YSZ//YSZ//LSM-YSZ固体氧化物燃料电池在800℃下的阻抗谱虚部图，阳极燃料为1H2:4N2(3vol.％H2O)，阴极气氛为空气，频率单位为Hz；负的阻抗虚部单位为Ωcm2；
图6为该燃料电池阻抗谱的Kramers-Kronig检验结果图，纵坐标为相对误差，×为阻抗实部的相对误差，100*(Z'KK-Z'+﹤R0﹥)/RP为阻抗实部相对误差的定义，+为阻抗虚部的相对误差，100*(Z”KK-Z”)/RP为阻抗虚部相对误差的定义；
图7为该燃料电池阻抗谱的弛豫时间分布图。
具体实施方式
具体实施方式一：结合图1说明本实施方式，一种计算交流阻抗谱弛豫时间分布的方法具体是按照以下步骤进行的：
步骤一、获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；
步骤二、对阻抗实部和阻抗虚部做Kramers-Kronig检验，使步骤一得到的交流阻抗谱数组是稳定的并且可以解析的，其中，所述Kramers-Kronig检验为实部和虚部的检验；
步骤三、在确定步骤一得到的交流阻抗谱是稳定的并且可以解析的基础上，根据频率和阻抗虚部构建弛豫时间和弛豫时间分布函数的代数方程组；
步骤四、应用Tikhonov正则化方法和二次规划方法求解弛豫时间和弛豫时间分布函数的代数方程组，得到弛豫时间及弛豫时间分布函数数组{τn，F(τn)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，所述图的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗，其中，所述Tikhonov正则化方法为吉洪诺夫正则化方法。
具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部；具体过程为：
交流阻抗谱数组由电化学工作站或模拟等效电路测得；
(1)模拟等效电路测得交流阻抗谱数组的过程为：
模拟等效电路选a个RC电路的串联电路，a取值范围为任意正整数，其中R为RC电路的电阻；R单位为Ωcm2，为任意正数值；C为RC电路的电容；C单位为F/cm2，为任意正数值；其中，所述RC电路为相移电路，Ωcm2为欧姆·平方厘米，F/cm2为法/平方厘米；
交流阻抗的频率范围为107Hz～10-4Hz，每频率数量级取x个离散的频率数据和阻抗数值，x取值为10到100间的整数，取值越高阻抗谱数组越多，其中，所述每频率数量级为在频率之比为10的两个频率，频率单位为Hz，得出该等效电路的模拟阻抗谱的Nyquist图和虚部图，获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部，其中，所述Hz为赫兹；
(2)电化学工作站测得交流阻抗谱数组的过程为：
电化学电池选为固体氧化物燃料电池，测试温度为500～1000℃，阳极气氛为氢气、合成气或碳氢化合物，阴极气氛为空气或氧气；
交流阻抗的频率范围为106Hz～10-2Hz，每频率数量级取x个离散的频率数据和阻抗数值，x取值为10到100间的整数，取值越高阻抗谱数组越多，其中，所述每频率数量级为在频率之比为10的两个频率，频率单位为Hz，例如测试频率范围106～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'(ω)计算出阻抗虚部Z”(ω)cal，根据步骤一得到的阻抗虚部Z”(ω)计算出阻抗实部Z'(ω)cal；
2)对比计算出的阻抗实部Z'(ω)cal和阻抗虚部Z”(ω)cal与步骤一得到的阻抗实部Z'(ω)和阻抗虚部Z”(ω)，若计算出的阻抗实部Z'(ω)cal与步骤一得到的阻抗实部Z'(ω)的相对误差和计算出的阻抗虚部Z”(ω)cal与步骤一得到的阻抗虚部Z”(ω)的相对误差都在10％以内，则步骤一得到的交流阻抗谱是稳定的并且可以解析的。
其它步骤及参数与具体实施方式一或二相同。
具体实施方式四：本实施方式与具体实施方式一、二或三不同的是：所述步骤三中在确定步骤一得到的交流阻抗谱是稳定的并且可以解析的基础上，根据频率和阻抗虚部构建弛豫时间和弛豫时间分布函数的代数方程组；具体过程为：
弛豫时间τ与频率ω的函数关系为ωτ＝1；
弛豫时间分布函数F(τ)与步骤一得到的阻抗虚部Z”(ω)、频率ω和弛豫时间τ的函数关系为：
式中，τ为弛豫时间，单位为s；ω为频率，单位为rad/s；F(τ)为弛豫时间分布函数，单位为Ωcm2；Z”(ω)为步骤一得到的阻抗虚部，单位为Ωcm2；
由于阻抗谱数组是离散形式的，所以将离散化，求得第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时，Fn为第n个待求解的弛豫时间分布函数；
Z是N+1行列矩阵；当n＝1,2,…,N时，Z'为阻抗谱实部，Wm＝log10ωm，ωm为第m个频率；wn＝-log10τ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时，Dn,[n-1,n,n+1]＝{-1,2,-1}；D1,[1,2]＝{1,-1}；DN,[N-1,N]＝{-1,1}，其他元素为0；λ为Tikhonov正则化系数，Tikhonov正则化系数λ为10-4～0；得到弛豫时间及弛豫时间分布函数数组{τn，F(τn)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。
其它步骤及参数与具体实施方式一、二、三或四相同。
采用以下实施例验证本发明的有益效果：
实施例1
一种计算交流阻抗谱弛豫时间分布的方法具体是按照以下步骤进行的：
1)、通过模拟电路获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部。该实施方式的模拟电路选为三个RC电路的串联电路，其中R1＝R2＝R3＝1Ωcm2；C1＝1e-3F/cm2；C2＝5e-5F/cm2；C3＝1e-5F/cm2。频率范围是1e-2～1e6Hz，每频率数量级取10个阻抗数据。图2给出了该等效电路的模拟阻抗谱的Nyquist图。图3给出了相应的虚部图。由图2和图3可见，只有两个弧呈现出来。根据Nyquist和虚部图并不能确定该阻抗谱是由三个RC等效过程产生的。
2)、根据频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组。弛豫时间τ[s]与频率ω[rad/s]的函数关系为ωτ＝1，弛豫时间分布函数F(τ)[Ωcm2]与阻抗谱虚部Z”(ω)[Ωcm2]、频率ω和弛豫时间τ的函数关系为:
根据频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组为：
Γ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'为阻抗谱实部；Wn＝log10ωn；wn＝-log10τ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时，Dn,[n-1,n,n+1]＝{-1,2,-1}；D1,[1,2]＝{1,-1}；DN,[N-1,N]＝{-1,1}，其他元素为0；λ为Tikhonov正则化系数，λ取值为10-4。应用MATLAB软件的内置函数quadprog求解在弛豫时间分布函数非负条件下目标函数err最小化的二次规划问题。图4给出了该等效电路的弛豫时间分布的计算结果。由图4可见三个峰，这三个峰所围的面积分别为1.0086Ωcm2、0.9926Ωcm2和1.0072Ωcm2。与这三个RC电路的阻抗1Ωcm2基本一致。峰顶对应的弛豫时间分别为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℃，阳极气氛为1H2:4N2(3vol.％H2O)，阴极气氛为空气，频率范围是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'(ω)计算出阻抗虚部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(τ)[Ωcm2]与阻抗谱虚部Z”(ω)[Ωcm2]、频率ω和弛豫时间τ的函数关系为:
由频率和阻抗虚部数组构建弛豫时间和弛豫时间分布函数的代数方程组为：
Γ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'为阻抗谱实部；Wm＝log10ωm，ωm为第m个频率；wn＝-log10τ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时，Dn,[n-1,n,n+1]＝{-1,2,-1}；D1,[1,2]＝{1,-1}；DN,[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单位为Ωcm2，为任意正数值；C为RC电路的电容；C单位为F/cm2，为任意正数值；其中，所述RC电路为相移电路，Ωcm2为欧姆·平方厘米，F/cm2为法/平方厘米；
交流阻抗的频率范围为107Hz～10-4Hz，每频率数量级取x个离散的频率数据和阻抗数值，x取值为10到100间的整数，取值越高阻抗谱数组越多，其中，所述每频率数量级为在频率之比为10的两个频率，频率单位为Hz，得出该等效电路的模拟阻抗谱的Nyquist图和虚部图，获得交流阻抗谱数组，包括频率、阻抗实部和阻抗虚部，其中，所述Hz为赫兹；
(2)电化学工作站测得交流阻抗谱数组的过程为：
电化学电池选为固体氧化物燃料电池，测试温度为500～1000℃，阳极气氛为氢气、合成气或碳氢化合物，阴极气氛为空气或氧气；
交流阻抗的频率范围为106Hz～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关系：
<mrow>
<msup>
<mi>Z</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msup>
<msub>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mi>c</mi>
<mi>a</mi>
<mi>l</mi>
</mrow>
</msub>
<mo>=</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mn>2</mn>
<mi>&omega;</mi>
</mrow>
<mi>&pi;</mi>
</mfrac>
<mo>)</mo>
</mrow>
<msubsup>
<mo>&Integral;</mo>
<mn>0</mn>
<mi>&infin;</mi>
</msubsup>
<mfrac>
<mrow>
<msup>
<mi>Z</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msup>
<mi>Z</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
<mo>-</mo>
<msup>
<mi>&omega;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mi>d</mi>
<mi>x</mi>
<mo>;</mo>
</mrow>
<mrow>
<msup>
<mi>Z</mi>
<mo>&prime;</mo>
</msup>
<msub>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
<mrow>
<mi>c</mi>
<mi>a</mi>
<mi>l</mi>
</mrow>
</msub>
<mo>=</mo>
<msup>
<mi>Z</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>&infin;</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mfrac>
<mn>2</mn>
<mi>&pi;</mi>
</mfrac>
<msubsup>
<mo>&Integral;</mo>
<mn>0</mn>
<mi>&infin;</mi>
</msubsup>
<mfrac>
<mrow>
<msup>
<mi>xZ</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msup>
<mi>&omega;Z</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
<mo>-</mo>
<msup>
<mi>&omega;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mi>d</mi>
<mi>x</mi>
<mo>;</mo>
</mrow>
由步骤一得到的阻抗实部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”(ω)、频率ω和弛豫时间τ的函数关系为：
<mrow>
<mo>-</mo>
<msup>
<mi>Z</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msubsup>
<mo>&Integral;</mo>
<mn>0</mn>
<mi>&infin;</mi>
</msubsup>
<mfrac>
<mrow>
<mi>&omega;</mi>
<mi>F</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mo>/</mo>
<mi>l</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<mn>10</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mn>1</mn>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<mi>&omega;</mi>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mi>d</mi>
<mi>&tau;</mi>
<mo>;</mo>
</mrow>
式中，τ为弛豫时间，单位为s；ω为频率，单位为rad/s；F(τ)为弛豫时间分布函数，单位为Ωcm2；Z”(ω)为步骤一得到的阻抗虚部，单位为Ωcm2；
由于阻抗谱数组是离散形式的，所以将离散化，求得第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时，Fn为第n个待求解的弛豫时间分布函数；
Z是N+1行列矩阵；当n＝1,2,…,N时，Z'为阻抗谱实部，Wm＝log10ωm，ωm为第m个频率；wn＝-log10τ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时，Dn,[n-1,n,n+1]＝{-1,2,-1}；D1,[1,2]＝{1,-1}；DN,[N-1,N]＝{-1,1}，其他元素为0；λ为Tikhonov正则化系数，Tikhonov正则化系数λ为10-4～0；得到弛豫时间及弛豫时间分布函数数组{τn，F(τn)}，以弛豫时间的对数为横轴、弛豫时间分布函数为纵轴作图，图中的各个峰对应不同电化学过程，峰面积代表不同电化学过程的实际阻抗。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510341876.7A CN104914312B (zh) | 2015-06-18 | 2015-06-18 | 一种计算交流阻抗谱弛豫时间分布的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510341876.7A CN104914312B (zh) | 2015-06-18 | 2015-06-18 | 一种计算交流阻抗谱弛豫时间分布的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104914312A CN104914312A (zh) | 2015-09-16 |
| CN104914312B true CN104914312B (zh) | 2018-01-30 |
Family
ID=54083552
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510341876.7A Active CN104914312B (zh) | 2015-06-18 | 2015-06-18 | 一种计算交流阻抗谱弛豫时间分布的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104914312B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3848713A4 (en) * | 2018-07-23 | 2022-04-13 | National Institute for Materials Science | ANALYSIS PROCESSING METHODS USING IMPEDANCE SPECTRAL DATA, IMPEDANCE SPECTRAL DATA ANALYSIS PROCESSING SYSTEM AND IMPEDANCE SPECTRAL ANALYSIS PROCESSING PROGRAM |
| US11977124B2 (en) | 2022-03-23 | 2024-05-07 | Honda Motor Co., Ltd | Measurement device, measurement method, and storage medium |
Families Citing this family (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105137362B (zh) * | 2015-10-19 | 2017-12-29 | 清华大学 | 一种电堆的无损在线检测和故障诊断方法 |
| JP6828739B2 (ja) * | 2016-04-12 | 2021-02-10 | 株式会社村田製作所 | 解析装置、解析方法、蓄電装置、蓄電システム、電子機器、電動車両および電力システム |
| CN109307688B (zh) * | 2017-07-27 | 2021-06-01 | 通用电气公司 | 感测系统和方法 |
| CN111521645B (zh) * | 2020-03-18 | 2022-05-06 | 昆明理工大学 | 一种锌电积过程中对阴阳极实时在线测量的装置 |
| CN111781246B (zh) * | 2020-07-20 | 2023-02-07 | 哈尔滨工业大学 | 一种原位表征碳/氮/氧单向反应扩渗化学弛豫过程的直流阻抗谱方法 |
| CN112540316B (zh) * | 2020-11-03 | 2022-07-05 | 华中科技大学 | 一种复杂电池阻抗谱解析方法 |
| CN112327172B (zh) * | 2020-11-30 | 2021-09-03 | 同济大学 | 一种基于弛豫时间分布的锂离子电池建模方法 |
| CN112327171B (zh) * | 2020-11-30 | 2021-11-09 | 同济大学 | 一种基于弛豫时间分布的锂离子电池寿命估计方法 |
| CN112763545B (zh) * | 2020-12-30 | 2022-08-26 | 宁德新能源科技有限公司 | 锂离子电池eis的交流阻抗数据处理与解读的方法及电池测试设备 |
| CN112816895B (zh) * | 2020-12-31 | 2024-04-12 | 上海派能能源科技股份有限公司 | 电化学阻抗谱的分析方法、系统、设备及计算机存储介质 |
| CN113283117B (zh) * | 2021-06-17 | 2022-11-22 | 清华大学 | 一种抗干扰的燃料电池阻抗解析方法 |
| CN114137429A (zh) * | 2021-10-29 | 2022-03-04 | 国电南瑞科技股份有限公司 | 充放电过程中锂离子电池性能异常变化的参数化表征方法及装置 |
| CN115015092A (zh) * | 2022-05-07 | 2022-09-06 | 苏州科技大学 | 一种电化学阻抗谱等效模拟电路选取方法及系统 |
| WO2024090724A1 (ko) * | 2022-10-28 | 2024-05-02 | 주식회사 엘지에너지솔루션 | 전지 불량 고속 검사 방법 |
| CN117589847B (zh) * | 2024-01-04 | 2024-04-16 | 中国第一汽车股份有限公司 | 一种氧传感器起燃时间测试方法及测试电路 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0915277A (ja) * | 1995-04-24 | 1997-01-17 | Matsushita Electric Ind Co Ltd | インピーダンス演算方法及びインピーダンス演算装置 |
| CN101871974A (zh) * | 2010-06-18 | 2010-10-27 | 华南理工大学 | 一种阻抗谱的测量方法 |
| CN101952714A (zh) * | 2007-12-20 | 2011-01-19 | 3M创新有限公司 | 用于电化学阻抗谱的方法和系统 |
| CN103149441A (zh) * | 2013-03-12 | 2013-06-12 | 中国科学院半导体研究所 | 应用于电化学测量的便携式阻抗谱分析仪及阻抗谱分析方法 |
| KR101386344B1 (ko) * | 2013-05-10 | 2014-04-16 | 부경대학교 산학협력단 | 푸리에 변환 전기화학 임피던스의 실시간 모니터링장치 및 방법 |
| CN203745542U (zh) * | 2013-12-16 | 2014-07-30 | 天津科技大学 | 一种快速无损阻抗谱测量系统 |
| CN104267261A (zh) * | 2014-10-29 | 2015-01-07 | 哈尔滨工业大学 | 基于分数阶联合卡尔曼滤波的二次电池简化阻抗谱模型参数在线估计方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE102009000337A1 (de) * | 2009-01-21 | 2010-07-22 | Robert Bosch Gmbh | Verfahren zur Bestimmung eines Alterungszustandes einer Batteriezelle mittels Impedanzspektroskopie |
-
2015
- 2015-06-18 CN CN201510341876.7A patent/CN104914312B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0915277A (ja) * | 1995-04-24 | 1997-01-17 | Matsushita Electric Ind Co Ltd | インピーダンス演算方法及びインピーダンス演算装置 |
| CN101952714A (zh) * | 2007-12-20 | 2011-01-19 | 3M创新有限公司 | 用于电化学阻抗谱的方法和系统 |
| CN101871974A (zh) * | 2010-06-18 | 2010-10-27 | 华南理工大学 | 一种阻抗谱的测量方法 |
| CN103149441A (zh) * | 2013-03-12 | 2013-06-12 | 中国科学院半导体研究所 | 应用于电化学测量的便携式阻抗谱分析仪及阻抗谱分析方法 |
| KR101386344B1 (ko) * | 2013-05-10 | 2014-04-16 | 부경대학교 산학협력단 | 푸리에 변환 전기화학 임피던스의 실시간 모니터링장치 및 방법 |
| CN203745542U (zh) * | 2013-12-16 | 2014-07-30 | 天津科技大学 | 一种快速无损阻抗谱测量系统 |
| CN104267261A (zh) * | 2014-10-29 | 2015-01-07 | 哈尔滨工业大学 | 基于分数阶联合卡尔曼滤波的二次电池简化阻抗谱模型参数在线估计方法 |
Non-Patent Citations (1)
| Title |
|---|
| 铝锌系阳极型强屏蔽涂层电化学阻抗谱分析;周瑶等;《涂料工业》;20070831;第37卷(第8期);第50-52页 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3848713A4 (en) * | 2018-07-23 | 2022-04-13 | National Institute for Materials Science | ANALYSIS PROCESSING METHODS USING IMPEDANCE SPECTRAL DATA, IMPEDANCE SPECTRAL DATA ANALYSIS PROCESSING SYSTEM AND IMPEDANCE SPECTRAL ANALYSIS PROCESSING PROGRAM |
| US11977124B2 (en) | 2022-03-23 | 2024-05-07 | Honda Motor Co., Ltd | Measurement device, measurement method, and storage medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104914312A (zh) | 2015-09-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104914312B (zh) | 一种计算交流阻抗谱弛豫时间分布的方法 | |
| CN107681181B (zh) | 一种燃料电池的性能诊断方法 | |
| Forman et al. | Genetic identification and fisher identifiability analysis of the Doyle–Fuller–Newman model from experimental cycling of a LiFePO4 cell | |
| CN105116337B (zh) | 一种锂离子电池满电荷存储寿命评价方法 | |
| CN105137362B (zh) | 一种电堆的无损在线检测和故障诊断方法 | |
| CN110488194A (zh) | 一种基于电化学阻抗模型的锂电池soc估算方法及其系统 | |
| Lyu et al. | SOH estimation of lithium-ion batteries based on fast time domain impedance spectroscopy | |
| CN105606641B (zh) | 一种在线实时监测锂电池隔膜热收缩率系统及监测方法 | |
| CN108666600B (zh) | 一种基于热化学测量的全钒液流电池soc检测方法 | |
| Li et al. | Coupling multi-physics simulation and response surface methodology for the thermal optimization of ternary prismatic lithium-ion battery | |
| CN104836223A (zh) | 电网参数错误与不良数据协同辨识与估计方法 | |
| CN107145629A (zh) | 一种优化电池电极厚度的方法 | |
| CN103728570B (zh) | 一种基于电池热特性的健康状态检测方法 | |
| CN115453377B (zh) | 基于电化学-热-老化与三维降阶的电池组寿命预测方法 | |
| CN104502844A (zh) | 一种基于交流阻抗的动力锂电池劣化程度诊断方法 | |
| CN110187287A (zh) | 一种退役锂电池余能快速检测方法 | |
| CN109738806A (zh) | 模拟电池产热率的方法、装置、介质 | |
| CN113403645A (zh) | 一种电解槽工作状态的确定方法、装置及控制器 | |
| Zhang et al. | A fractional-order model of lithium-ion batteries and multi-domain parameter identification method | |
| CN114460474A (zh) | 电池分容方法及其装置、电子设备 | |
| Khatib et al. | Experimental and analytical study of open pore cellular foam material on the performance of proton exchange membrane electrolysers | |
| Kreller et al. | Modeling SOFC cathodes based on 3-D representations of electrode microstructure | |
| CN105699902B (zh) | 用于燃料电池诊断的阻抗测定装置及其方法 | |
| Xin et al. | A novel state of charge estimation method for ternary lithium batteries based on system function and extended kalman filter | |
| CN105319528B (zh) | 一种电能表的运行工况检验方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |