李琛坤，王 帥，黃 俊

（1中南大學(xué)化學(xué)化工學(xué)院，湖南 長(zhǎng)沙 410083；2烏爾姆大學(xué)理論化學(xué)研究所，德國(guó) 烏爾姆 89069）

電化學(xué)阻抗譜(electrochemical impedance spectroscopy，EIS)具有原位、無(wú)損的優(yōu)點(diǎn)，被廣泛應(yīng)用于電化學(xué)儲(chǔ)能(二次電池、燃料電池)和轉(zhuǎn)化器件(電解池等)的性能監(jiān)測(cè)、故障診斷和參數(shù)辨識(shí)[1-3]。在EIS實(shí)驗(yàn)中，我們用電化學(xué)工作站先將測(cè)試體系控制在一個(gè)穩(wěn)定狀態(tài)，然后施加微弱的電壓或電流激勵(lì)，測(cè)試相應(yīng)的電流或電壓響應(yīng)。傅里葉變換后的電壓和電流的比值即定義為阻抗。

其中，F(xiàn)表示傅里葉變換算子。測(cè)量不同頻率下的阻抗便構(gòu)成了EIS。

對(duì)測(cè)量得到的EIS 進(jìn)行解析，需要相應(yīng)的模型。常見(jiàn)的模型分為等效電路模型[4-5]和物理模型[1,6-7]。等效電路模型是指用電路元件(電阻、電容和電感等)的串并聯(lián)來(lái)模擬研究體系。如圖1 所示，物理模型是指將研究體系內(nèi)發(fā)生的物理過(guò)程用相應(yīng)的控制方程來(lái)描述。求解控制方程，即得到EIS的解析解或數(shù)值解。李響等人[7]詳細(xì)推導(dǎo)了多孔電極EIS 的表達(dá)式，并細(xì)致分析了不同過(guò)程對(duì)EIS 響應(yīng)的影響。等效電路模型的優(yōu)點(diǎn)在于通用性好、計(jì)算量小和操作簡(jiǎn)單，缺點(diǎn)是物理機(jī)理不明確。程蕾等[8]推導(dǎo)了玻碳電極界面EIS的等效電路模型。通過(guò)擬合實(shí)驗(yàn)數(shù)據(jù)，定量對(duì)比分析了修飾前后電極表面的特征變化。由等效電路模型衍生而來(lái)的弛豫時(shí)間分布(distribution of relaxation times，DRT)[9]可以有效辨識(shí)不同過(guò)程的時(shí)間常數(shù)和阻抗，一定程度上辨析內(nèi)部機(jī)理，因此得到了廣泛的應(yīng)用。物理模型的優(yōu)點(diǎn)在于對(duì)體系內(nèi)的物理過(guò)程描述清晰，適用于研究?jī)?nèi)部機(jī)理?；诳刂品匠探⒌奈锢砟Ｐ蛣t有助于解析實(shí)驗(yàn)數(shù)據(jù)，辨析機(jī)理。因此，本文主要圍繞物理模型展開介紹。

圖1 (a)EIS測(cè)試原理；(b)典型的EIS實(shí)驗(yàn)結(jié)果；(c)物理模型構(gòu)建原理(其中i0為交換電流密度，c為濃度，D為擴(kuò)散系數(shù))；(d)2階RC等效電路模型Fig.1 (a)Principle of EIS measurement;(b)Typical experimental results of EIS;(c)Principle of physical model construction,(i0 is the exchange current density,c the concentration,D the diffusion coefficient);(d)Equivalent circuit model of second-order RC

本文首先介紹了EIS物理模型三種常見(jiàn)的求解方法：解析求解，時(shí)域數(shù)值方法和頻域數(shù)值方法。其中，解析求解指通過(guò)合理假設(shè)，推導(dǎo)出EIS的解析表達(dá)式。時(shí)域數(shù)值方法指利用有限元等方法在時(shí)間-空間維度求解微分方程，而后通過(guò)傅里葉變化得到EIS 的數(shù)值解。頻域數(shù)值方法指在頻率-空間維度求解微分方程，而后計(jì)算EIS數(shù)值解。介紹完三種方法后，本文以電池中常見(jiàn)的金屬離子沉積反應(yīng)為例，對(duì)比了不同方法的優(yōu)缺點(diǎn)，并對(duì)每種方法的適用場(chǎng)景進(jìn)行了分析。

1 三種計(jì)算方法介紹

下文將以一個(gè)經(jīng)典的帶有法拉第反應(yīng)的雙電層界面為例，具體介紹上述三種方法。筆者[10]此前對(duì)雙電層的平衡、非平衡、交流狀態(tài)下的建模進(jìn)行了詳細(xì)的推導(dǎo)和總結(jié)，可以作為初學(xué)者的指南，感興趣的讀者可以進(jìn)行查閱。筆者[11]同時(shí)梳理了EIS 在電催化界面和反應(yīng)中的研究和最新的進(jìn)展。相比之下，本文重點(diǎn)關(guān)注EIS物理模型的求解方法。

圖2 (a)帶氧化還原反應(yīng)的電極/電解質(zhì)溶液界面示意圖；(b)EIS解析解的Nyquist圖；(c)基于時(shí)域數(shù)值方法的EIS數(shù)值解與解析解的對(duì)比；(d)基于頻域數(shù)值方法的EIS數(shù)值解與解析解的對(duì)比。頻率范圍為105～10-4 Hz。其余參數(shù)如下：c0 = 0.1 mol/L，D = 1 × 10-10 m2/s，j0 = 1 A/m2，α = 0.5Fig.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 105～10-4 Hz.Other parameters are as follows:c0=0.1 mol/L,D=1×10-10 m2/s,j0=1 A/m2,α=0.5

溶液中離子傳輸僅由濃度梯度驅(qū)動(dòng)，控制方程為菲克第二定律。

其中，Di為離子的擴(kuò)散系數(shù)，i= Ox表示氧化物，i= Re表示還原物。

假設(shè)本體電解質(zhì)位于無(wú)窮遠(yuǎn)處，x= ∞，離子濃度為本體濃度

在電極表面，x= 0，反應(yīng)物被消耗，產(chǎn)物被生成。消耗/生成的速率和反應(yīng)速率的關(guān)系如下

其中負(fù)號(hào)對(duì)應(yīng)氧化物，正號(hào)對(duì)應(yīng)還原物，氧化反應(yīng)的電流定義為正電流，F(xiàn)為法拉第常數(shù)，n為轉(zhuǎn)移的電子數(shù)，jct為法拉第反應(yīng)的電流，一般通過(guò)Butler-Volmer方程計(jì)算[12-13]

其中，j0為交換電流密度，α為電荷轉(zhuǎn)移系數(shù)，F(xiàn)為法拉第常數(shù)，R為理想氣體常數(shù)，T為溫度，φM為電極電勢(shì)，φeq為平衡電勢(shì)，利用能斯特方程計(jì)算[12-13]

其中φ為標(biāo)準(zhǔn)平衡電勢(shì)，γi為活度，c為表面濃度。

總電流除了反應(yīng)電流，還有雙電層電流jdl

其中，Cdl為雙電層電容，一般依據(jù)Gouy-Chapman-Stern模型計(jì)算[12-13]

其中，CGC=?s/λD為Gouy-Chapman電容，?s為溶液本體的介電常數(shù)，λD為德拜長(zhǎng)度。CH=?HP/δHP為Helmholtz 電容，?HP和δHP分別為Helmholtz 層的介電常數(shù)和厚度。

1.1 解析求解

對(duì)方程式(2)施加傅里葉變換可以得到頻域的控制方程

傅里葉變換后的邊界條件如下，本體x= ∞處

電極表面，x= 0處

結(jié)合邊界條件式(10)和式(11)，可以得到方程式(9)的解析解為

對(duì)Bulter-Volmer 方程做線性化展開，利用傅里葉變換可以得到

在稀溶液中，γOx=γRe= 1，把方程式(14)代入方程式(13)得到

對(duì)方程式(7)做傅里葉變換，可以得到頻域的雙電層電流

1.2 時(shí)域數(shù)值方法

其中，f(t)為系統(tǒng)的時(shí)域響應(yīng)信號(hào)，F(xiàn)為傅里葉變換算子。

對(duì)信號(hào)f(t)做線性插值處理可以得到

其中，σ(t)為歸一化的階躍函數(shù)。

把方程式(19)代入方程式(18)可以得到AFT 的表達(dá)式

基于TSM 求解，再通過(guò)AFT 變換得到的EIS數(shù)值解和解析解的對(duì)比如圖2(c)所示。在高頻區(qū)域，數(shù)值解和解析解基本完全重合。低頻區(qū)域的偏差略有增加，主要原因是有限元等方法求解微分方程會(huì)累計(jì)誤差。

1.3 頻域數(shù)值方法

基于頻域數(shù)值方法得到的EIS 數(shù)值解和解析解的對(duì)比如圖2(d)所示(具體的MATLAB 程序放在支撐文件)?？梢悦黠@看到，數(shù)值解在整個(gè)頻率范圍內(nèi)都和解析解很吻合，表明FSM 數(shù)值解精度較高。

2 求解案例

2.1 金屬離子沉積反應(yīng)體系

2.1.1 控制方程

溶液中離子的傳輸采用Poisson-Nernst-Planck(PNP)方程[14,20-22]描述。

其中，c+和c-分別為陽(yáng)離子和陰離子的濃度，D為擴(kuò)散系數(shù)，?s為本體溶液的介電常數(shù)。

2.1.2 邊界條件和初始條件

研究體系的邊界分為左邊界和右邊界，左邊界為Helmholtz 平面，x= 0 處。邊界條件設(shè)置如下，對(duì)于陽(yáng)離子，發(fā)生電子轉(zhuǎn)移反應(yīng)，陽(yáng)離子通量和反應(yīng)電流的關(guān)系見(jiàn)式(24)

其中，jr為法拉第反應(yīng)的電流，這里采用Frumkin-Butler-Volmer(FBV)方程[13]計(jì)算

其中，c0為陽(yáng)離子的本體濃度，k0為反應(yīng)速率常數(shù)，η為反應(yīng)的過(guò)電勢(shì)。η=φM-φHP-φ，φHP為Helmholtz平面處的電勢(shì)。

對(duì)于陰離子，不參與反應(yīng)，所以離子的通量為0

對(duì)于電勢(shì)，從電極表面到Helmholtz 平面，電勢(shì)線性下降[22]，所以電勢(shì)的邊界條件為

由于所有的電荷均為體電荷，Helmholtz 平面體積為0，不存在凈電荷。因此，根據(jù)高斯定律，Helmholtz平面兩側(cè)的電勢(shì)梯度滿足如下關(guān)系

方程式(28)代入方程式(27)可得電勢(shì)左邊界的最終形式

右邊界為溶液本體，x=xb，電勢(shì)為0，離子濃度均為本體濃度

結(jié)合邊界條件對(duì)控制方程進(jìn)行求解，我們可以得到零電荷電勢(shì)(Potential of zero charge，PZC)下的解析表達(dá)式[10]

其 中，rc=CGC/CH為Gouy-Chapman 電 容 和Helmholtz電容的比值，ωnd=ω/D為無(wú)量綱的角頻率，Xb=xb/λD為無(wú)量綱的電解質(zhì)溶液厚度，λ1= 1 +jωnd，λ2=jωnd，和為沉積反應(yīng)相關(guān)的變量。

其中，ν1=k0{αexp(αΠ0)+(1-α)×exp[-(1 -α)Π]0]}，ν2=k0exp[ -(1-α)Π0]，Π0=(--φeq)F/RT為穩(wěn)態(tài)的過(guò)電勢(shì)，上標(biāo)0 表示該變量為穩(wěn)態(tài)變量。

2.1.3 數(shù)值求解

需要注意的是，在金屬沉積反應(yīng)中，邊界條件中嚴(yán)格考慮了雙電層的存在。因此，雙電層電流不能再以Cdl為常數(shù)進(jìn)行計(jì)算，而應(yīng)該從電荷量求導(dǎo)去計(jì)算[22-23]

其中，qdl表示雙電層中積累的凈電荷量，負(fù)號(hào)表示積累的電荷和電極所帶的電荷相反。

總電流的表達(dá)式如下

基于頻域方程的數(shù)值求解則先要對(duì)控制方程和邊界條件做傅里葉變換，具體的計(jì)算詳情這里不再贅述[22-23]。頻域下的控制方程(這里只考慮PZC下的情況)如下

在x= 0處，頻域下的左邊界條件為

在x=xb處，頻域下的右邊界條件為

2.1.4 三種方法結(jié)果對(duì)比

圖3(b)～(d)依次為金屬沉積反應(yīng)EIS 解析解，時(shí)域數(shù)值方法的EIS數(shù)值解與解析解的對(duì)比，和頻域數(shù)值方法的EIS 數(shù)值解和解析解的對(duì)比(具體的MATLAB 程序放在支撐文件)。需要注意的是，這里的解析解和數(shù)值解都只適用于PZC[22-23]。對(duì)于非PZC下的求解，初始條件不再滿足均勻假設(shè)，離子濃度和電勢(shì)均存在分布。特別是在雙電層中，電勢(shì)梯度和濃度梯度極大。筆者尚未看到有相關(guān)研究在考慮界面反應(yīng)和雙電層的條件下，自洽求解了非PZC下的EIS。

圖3 (a)金屬離子沉積反應(yīng)電極/電解質(zhì)溶液界面示意圖；(b)EIS解析解的Nyquist圖；(c)時(shí)域數(shù)值方法的EIS數(shù)值解與解析解的對(duì)比；(d)頻域數(shù)值方法的EIS數(shù)值解與解析解的對(duì)比。頻率范圍為105～10-1 Hz。其余參數(shù)如下：c0 = 0.1 mol/L，D = 1 × 10-12 m2/s，k0 = 4 × 10-4 mol/(m2·s)，α = 0.5Fig.3 (a)Schematic diagram of electrode/electrolyte solution interface with metal ion deposition reaction;(b)Nyquist plot of EIS analytical solution;(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 105～10-1 Hz.Other parameters are as follows:c0=0.1 mol/L,D=1×10-12 m2/s,k0=4×10-4mol/(m2·s),α=0.5

從圖3(b)解析解的結(jié)果可以看出，EIS 由高頻反應(yīng)半圓和低頻擴(kuò)散半圓構(gòu)成。解析解作為精確解，用來(lái)評(píng)價(jià)不同數(shù)值方法求解EIS 的精度。圖3(c)為時(shí)域方法的數(shù)值解和解析解。二者在中低頻區(qū)域吻合較好，高頻區(qū)域數(shù)值解出現(xiàn)感抗行為。從筆者的經(jīng)驗(yàn)來(lái)看，主要是因?yàn)橛?jì)算雙電層電流采用的數(shù)值積分及微分存在較大的誤差，進(jìn)而導(dǎo)致高頻感抗的行為。圖3(d)對(duì)比了頻域方法的數(shù)值解和解析解。不管是高頻還是低頻，數(shù)值解和解析解均很好吻合，表明數(shù)值解精度較高。

從三種方法對(duì)比結(jié)果可以得出，解析解避免了數(shù)值求解微分方程，因而精度最高。但解析解往往需要大量數(shù)學(xué)運(yùn)算，才有可能推導(dǎo)出復(fù)雜體系的解析EIS。相比之下，數(shù)值解借助有限元等方法，可以直接得到體系的EIS數(shù)值解。對(duì)比不同數(shù)值方法的結(jié)果，可以得出，頻域數(shù)值方法的EIS數(shù)值解精度較高，而時(shí)域數(shù)值方法的EIS數(shù)值解，低頻區(qū)域精度較高，高頻區(qū)域精度較低。因此，筆者建議優(yōu)先選取頻域數(shù)值方法計(jì)算EIS。

2.2 單價(jià)離子和多價(jià)離子阻抗的對(duì)比

金屬離子沉積反應(yīng)中，常見(jiàn)的離子價(jià)態(tài)為一價(jià)、二價(jià)和三價(jià)。圖4 對(duì)比了不同價(jià)態(tài)離子的EIS(具體的MATLAB程序放在支撐文件)。其中，圖4(a)和(b)分別為TSM 數(shù)值方法和FSM 數(shù)值方法的結(jié)果。在TSM 數(shù)值方法的EIS 結(jié)果可以看出，高頻仍然存在感抗特征，其余部分和FSM 方法的結(jié)果一致，這也說(shuō)明了FSM 方法的求解精度更高。對(duì)于不同價(jià)態(tài)的離子，在保持其余參數(shù)不變的情況下，離子價(jià)態(tài)越高，阻抗越小。主要原因在于，一方面FBV方程中考慮離子價(jià)態(tài)后，反應(yīng)電流極大地增加。另一方面，通過(guò)積分電荷量計(jì)算的雙電層電流也增加。因此，總電流極大的增加使得高價(jià)態(tài)離子的EIS減小。但需要注意的是，對(duì)于真實(shí)的研究體系，不同價(jià)態(tài)離子的動(dòng)力學(xué)參數(shù)往往差別較大，讀者需要結(jié)合實(shí)驗(yàn)標(biāo)定相關(guān)參數(shù)，然后利用模型去擬合和解析實(shí)驗(yàn)數(shù)據(jù)。

圖4 (a)基于TSM數(shù)值方法的EIS結(jié)果對(duì)比；(b)基于FSM數(shù)值方法的EIS結(jié)果對(duì)比(除考慮不同離子價(jià)態(tài)外，其余參數(shù)均與圖3 保持相同)Fig.4 (a)comparison of EIS results based on TSM numerical method;(b)comparison of EIS results based on FSM numerical method(Except for different ion valence states,other parameters remain the same as Fig.3 )

3 三種方法優(yōu)缺點(diǎn)對(duì)比

圖5對(duì)比了三種求解方法的優(yōu)缺點(diǎn)。其中，解析求解優(yōu)點(diǎn)在于物理含義清晰，易于分析高低頻的漸進(jìn)行為。缺點(diǎn)在于推導(dǎo)復(fù)雜，復(fù)雜體系難以解析求解。數(shù)值方法優(yōu)點(diǎn)在于求解簡(jiǎn)單，無(wú)復(fù)雜推導(dǎo)。缺點(diǎn)在于計(jì)算耗時(shí)。其中，時(shí)域數(shù)值方法求解精度低于頻域數(shù)值方法。但頻域數(shù)值方法無(wú)法獲得時(shí)域信號(hào)，難以有效判斷中間變量的正確性。

圖5 三種求解方法優(yōu)缺點(diǎn)對(duì)比Fig.5 Comparison of advantages and disadvantages of three solving methods

4 結(jié)論和展望

本文簡(jiǎn)要介紹了EIS 物理模型常見(jiàn)的三種求解方法：解析求解，時(shí)域數(shù)值方法和頻域數(shù)值方法。隨后結(jié)合具體的例子對(duì)比了不同方法的優(yōu)缺點(diǎn)，得到的主要結(jié)論如下：解析解精度最高，但一般只適用于PZC 條件；對(duì)于非PZC，考慮雙電層的條件下，無(wú)法求出解析解。對(duì)于數(shù)值解，優(yōu)先推薦頻域數(shù)值方法，精度較高且沒(méi)有高頻感抗行為。

為拓寬EIS在電化學(xué)中的應(yīng)用，尚有很多工作值得開展。

（1）本文的解析解和數(shù)值解都只適用于PZC條件，對(duì)于非PZC下的EIS求解，筆者尚未看到有相關(guān)研究。如何自洽求解非PZC下的EIS是一個(gè)值得研究的課題。

（2）本文只考慮了單電子轉(zhuǎn)移的反應(yīng)，而燃料電池、電解水等電化學(xué)裝置中的反應(yīng)均為4電子和2電子反應(yīng)。分步考慮電子轉(zhuǎn)移反應(yīng)，建立微觀動(dòng)力學(xué)模型，再計(jì)算EIS對(duì)于理解反應(yīng)速控步驟，辨析反應(yīng)機(jī)理等至關(guān)重要[24-25]。

（3）本文僅考慮了一維平板電極情況，而真實(shí)的研究體系多為粗糙界面、團(tuán)聚體顆粒和多孔電極。Song等人[26]和Kant等人[27-29]的研究表明，幾何效應(yīng)對(duì)EIS響應(yīng)有很大的影響。