二、PEMFC基础之电化学与反应动力学

1.电流、电流密度
2.反应速率常数
3.交换电流密度
4.电化学动力学奠基石B-V方程
5.活化损失计算Tafel公式
6.计算案例

1.电流、电流密度

由法拉第定律
i = d Q d t = n F d N d t i=\frac{dQ}{dt}=\frac{nFdN}{dt} i=dtdQ=dtnFdN j = i A j=\frac{i}{A} j=Ai
单位面积反应速率（mol·s-1·cm-2）：
υ = 1 A d N d t = i n F A = j n F \upsilon=\frac{1}{A}\frac{dN}{dt}=\frac{i}{nFA}=\frac{j}{nF} υ=A1dtdN=nFAi=nFj

2.反应速率常数

反应速率常数k（s-1）由过渡态理论获取：
k = k B T h e x p ( − △ G R T ) k=\frac{k_{B}T}{h}exp\left ( \frac{-\bigtriangleup G}{RT} \right ) k=hkBTexp(RT−△G)
其中，吉布斯自由能∆G可以分为化学能项（下标ac）和电能项两部分。
△ G = △ G a c + α R d F E r \bigtriangleup G=\bigtriangleup G_{ac}+\alpha _{Rd}FE_{r} △G=△Gac+αRdFEr △ G = △ G a c − α O x F E r \bigtriangleup G=\bigtriangleup G_{ac}-\alpha _{Ox}FE_{r} △G=△Gac−αOxFEr
注意：
a.α为转移系数，部分文献将对称因子β与α混淆使用，严格意义上来说是不对的，区别在于：
对称因子β严格用于涉及单个电子的单步反应，且阴阳极对称因子之和为1；
转移系数α描述多步过程，阴阳极的转移因子和不一定为1
b.阳极的传输系数一般取0.5，阴极在0.1~0.5之间比较合适。

3.交换电流密度

正向（下标f）和逆向（下标b）反应的单位面积反应速率又可以表示为：
υ f = k f C O x \upsilon_{f}=k_{f}C_{Ox} υf=kfCOx υ b = k b C R d \upsilon_{b}=k_{b}C_{Rd} υb=kbCRd
正向电流密度和逆向电流密度分别为：
j f = n F k 0 , f C O x e x p [ − α R d F E r R T ] j_{f}=nFk_{0,f}C_{Ox}exp\left [ -\frac{\alpha _{RdFE_{r}}}{RT} \right ] jf=nFk0,fCOxexp[−RTαRdFEr]
j b = n F k 0 , b C R d e x p [ − α O x F E r R T ] j_{b}=nFk_{0,b}C_{Rd}exp\left [ -\frac{\alpha _{OxFE_{r}}}{RT} \right ] jb=nFk0,bCRdexp[−RTαOxFEr]
当jf=jb时，即没有净电流输出时，此时电流密度为交换电流密度，jf=jb=j0
注意：
a.交换电流密度越大，活化过电势越小，净电流密度越大，总之该项越大电池性能越好。
b.阳极交换电流密度比阴极交换电流密度高几个量级。
c.有效交换电流密度与交换电流密度的区别：
交换电流密度代表电化学反应的内在动力学，由电极材料的性质、电化学反应和电极表面的反应物浓度决定；
有效交换电流密度不仅考虑了电化学反应的内在动力学，而且还考虑了燃料电池电极结构的影响，如催化剂层的孔隙率、催化剂负载和催化剂表面积的利用。在计算燃料电池的活化损失时，通常使用有效交换电流密度，因为它能更真实地反映实际电极性能。
有效交换电流密度的计算公式如下：
j = j 0 r e f a c L c ( P r P r r e f ) γ e x p [ − E c R T ( 1 − T T r e f ) ] j=j_{0}^{ref}a_{c}L_{c}\left ( \frac{P_{r}}{P_{r}^{ref}} \right )^{\gamma } exp\left [ -\frac{E_{c}}{RT}\left ( 1-\frac{T}{T_{ref}} \right ) \right ] j=j0refacLc(PrrefPr)γexp[−RTEc(1−TrefT)]

4.电化学动力学奠基石B-V方程

Bulter-Volumer Equation的两种形式如下：
j = j 0 [ e x p ( α R d F ( E − E r ) R T ) − e x p ( α O x F ( E − E r ) R T ) ] j=j_{0}\left [ exp\left (\frac{\alpha _{Rd}F\left ( E-E_{r} \right )}{RT} \right ) -exp\left (\frac{\alpha _{Ox}F\left ( E-E_{r} \right )}{RT} \right )\right ] j=j0[exp(RTαRdF(E−Er))−exp(RTαOxF(E−Er))]

j = j 0 [ e x p ( α n F η a c t R T ) − e x p ( − ( 1 − α ) n F η a c t R T ) ] j=j_{0}\left [ exp\left (\frac{\alpha nF\eta_{act}}{RT} \right ) -exp\left (\frac{-\left ( 1-\alpha \right ) nF\eta_{act} }{RT} \right )\right ] j=j0[exp(RTαnFηact)−exp(RT−(1−α)nFηact)]
注意：
a.可以看出，活化过电势越大，电流密度越大。
b.仔细观察，两种形式右侧分子项存在区别（是否有n），这是因为
α R d = α ∗ n \alpha _{Rd} = \alpha * n αRd=α∗n
因此，对于阳极来说n=2，阴极来说n=4。此外，在仿真时尤其是使用商业软件仿真时，需要注意软件中的BV公式是如何描述的。ps.FLUENT中燃料电池模块用的第一个公式。

5.活化损失计算Tafel公式

η a c t = a + b l n i \eta _{act}=a+blni ηact=a+blni a = − R T n F l n ( i 0 ) a=-\frac{RT}{nF}ln\left ( i_{0} \right ) a=−nFRTln(i0) b = − R T n F b=-\frac{RT}{nF} b=−nFRT
也可以表达为：
η a c t a n o d e + η a c t c a t h = R T n F α l n ( i i 0 ) a n o d e + R T n F α l n ( i i 0 ) c a t h o d e \eta _{actanode}+\eta _{actcath}=\frac{RT}{nF\alpha }ln\left ( \frac{i}{i_{0}}\right )_{anode}+\frac{RT}{nF\alpha }ln\left ( \frac{i}{i_{0}}\right )_{cathode} ηactanode+ηactcath=nFαRTln(i0i)anode+nFαRTln(i0i)cathode

考虑到阳极的交换电流密度大的多，因此一般可以忽略阳极的活化损失。

6.计算案例

%  0维模型计算极化曲线
clc;clear;
% 参数设定
R = 8.314;                      % 理想气体常数 j/mol*K
n = 4;                          % 每mol的O2转移的电子摩尔数
Alpha = 0.25;                   % 传输系数
i0 = 10^(-6.912);               % 交换电流密度
iL = 1.41;                      % 极限电流密度
F = 96485;                      % 法拉第常数
r = 0.19;                       % 内阻 Ω/cm2
Tk = 333;                       % K
Tc = 60;                        % ℃
P_H2    = 3;                    % 氢气压力 atm
P_air   = 3;                    % 空气压力 atm
Et = 1.19;
% step.1 压力计算及能斯特电压计算
P_H2O = 0.9869*10^(-2.1794+0.02953*Tc-9.1837e-5*Tc^2+1.4454e-7*Tc^3);       % 水饱和蒸气压 loop = 1;
for fid = 1 : 1400i = 0.001 * fid;             % 电流密度 A/cm2PP_H2(loop) = 0.5*P_H2./exp(1.653*i/(Tk^1.334))-P_H2O;            PP_O2(loop) = P_air./exp(4.192*i/(Tk^1.334))-P_H2O;                                 % 计算三部分电压损失B = R*Tk/(n*F*Alpha);V_act(loop) = -B * log(i/i0);     % 极化损失VV_ohmic(loop) = -(i*r);           % 欧姆损失Vterm = 1 - i/iL;if term > 0V_conc(loop) = (R*Tk/n/F)*(1+1/Alpha)*log(1-(i/iL));elseV_conc(loop) = 0;end% 能斯特电压计算 包括温度修正V_nerst(loop) =  Et - R*Tk*log(P_H2O/(PP_H2(loop)*PP_O2(loop)^0.5))/(2*F);V_out(loop) = V_nerst(loop) + V_act(loop) + V_ohmic(loop) + V_conc(loop);loop = loop + 1;
end
current_density = 0.001:0.001:1.40;figure1 = figure('color',[1 1 1]);
hdlp = plot(current_density,V_out);
title('Fuel cell polarization curve','FontSize',14,'FontWeight','Bold');
xlabel('Current density /A·cm-2','FontSize',12,'FontWeight','Bold')
ylabel('Activation losses /V','FontSize',12,'FontWeight','Bold')
set(hdlp,'LineWidth',1.5);
grid on;