一、理论基础
热传导主要基于傅里叶热传导定律。在燃料电池中,除了各组件内部的热传导外,还有冷却流体与双极板的对流换热。公式略。
燃料电池内部稳态导热:
d
2
T
d
x
2
+
q
i
n
t
k
=
0
\frac{d^{2}T}{dx^{2}}+\frac{q_{int}}{k}=0
dx2d2T+kqint=0
注意
1.如果界面间采用第三类边界条件,则在两组间交界面处温度连续;如果采用第二类边界条件,考虑接触热阻,则在界面处存在温差,温度不连续。
2、催化层中产热可以假设为除了输出电压、离子阻抗损失、电阻抗损失之外的损耗。
二、编程实践
0 ~ 0.1为冷却流道,0.1 ~ 0.42为双极板, 0.42 ~ 0.5为GDL区域。
% --------- 燃料电池中一维导热问题 ---------- %
% | | | | |
% | | | | |
% | | | | |
% | | | | |
% cooling BP GDL Cl
% Tf T1 T2/T3 T4
% Q1 Q2/Q3 Q4
clc;clear;
%% ------------------ Parameters ------------------
j = 1; % 电流密度 A/cm2
E = 0.6; % 输出电压 V
R_ionic = 0.1; % 离子阻抗 Ω·cm2
Rou_GDL = 0.08; % GDL的电阻率 Ω·cm
Rou_BP = 0.006; % BP的电阻率 Ω·cm
lamda_GDL = 0.017; % GDL的有效热导率 W·cm-1·K-1
lamda_BP = 0.19; % BP的有效热导率 W·cm-1·K-1
R_tc = 0.005; % GDL与Cl的接触电阻 Ω·cm2
R_tc1 = 1; % GDL与Cl的接触热阻 K/W
L_GDL = 0.04; % GDL的厚度 cm
L_BP = 0.32; % BP的厚度 cm
L_chanel = 0.1; % 气体流道的宽度
Er = 1.254; % 最大可逆电压
Tf = 60; % 阴极气体流道温度
h = 0.16; % 热传递系数 W·cm-2·K-1
%% ------------------ 各界面热通量 ------------------
% 计算各部件的欧姆损耗
V_GDL = j * Rou_GDL * L_GDL;
V_BP = j * Rou_BP * L_BP;
V_tc = R_tc * j;
% 假设阴阳极两侧的欧姆损耗相同
V_elec = 2 * (V_GDL + V_BP + V_tc);
V_ionic = R_ionic * j;
V_ohmic = V_elec + V_ionic;
% GDL/CL界面的热通量为除了离子阻抗和电阻抗之外所有损耗
V_h = Er - E - V_ohmic;
Q4 = V_h * j^2 + 0.5 * R_ionic * j^2;
Q_GDL = j * V_GDL;
Q_BP = j * V_BP;
Q_tc = j * V_tc;
Q1 = Q4 + Q_GDL + Q_BP + Q_tc;
Q2 = Q1 - Q_BP;
Q3 = Q2 - Q_tc;
%% ------------------ 求解温度分布 ------------------
dx = 0.01; % spatial grid size
N = 50; % grid number
T_dis = zeros(N+1,2);
T1 = Tf + Q1/h;
D = L_chanel+L_BP;
for i = 1 : N+1
x = dx * (i-1);
if i >= 11 && i < 43
q_intBP = Q_BP / L_BP;
T2 = T1 + (Q2/lamda_BP) * (x-L_chanel) + (q_intBP/lamda_BP) * (L_BP*(x-L_chanel)- ((x-L_chanel)^2)/2);
T_dis(i,2) = T2;
T_dis(i,1) = x;
elseif i == 43
T2 = T1 + (Q2/lamda_BP) * (x-L_chanel) + (q_intBP/lamda_BP) * (L_BP*(x-L_chanel)- ((x-L_chanel)^2)/2);
T_dis(i,2) = T2;
T_dis(i,1) = x;
T3 = T1 + (Q2/lamda_BP) * (x-L_chanel) + (q_intBP/lamda_BP) * (L_BP*(x-L_chanel)- ((x-L_chanel)^2)/2) + R_tc1 * Q3;
T_dis(i+1,2) = T3;
T_dis(i+1,1) = x;
elseif (i >= 44 && i < 52)
q_intGDL = Q_GDL / L_GDL;
T4 = T3 + (Q4/lamda_GDL) * (x-L_chanel-L_BP) + (q_intGDL/lamda_GDL) * (L_GDL*(x-L_chanel-L_BP)- ((x-L_chanel-L_BP)^2)/2);
T_dis(i+1,2) = T4;
T_dis(i+1,1) = dx * (i-1);
end
end
plot(T_dis(11:52,1),T_dis(11:52,2),'r*');
xlabel('x [cm]','FontSize',12,'FontWeight','Bold')
ylabel('Temperature [℃]','FontSize',12,'FontWeight','Bold')
ylim([63 68])
xlim([0 0.55])
line([0.42 0.42],[63 68],'Color','Blue','LineStyle','-','LineWidth',2);
line([0.1 0.1],[63 68],'Color','Blue','LineStyle','-','LineWidth',2);
line([0.55 0.55],[63 68],'Color','Blue','LineStyle','-','LineWidth',2);
line([0.5 0.5],[63 68],'Color','Blue','LineStyle','-','LineWidth',2);
line([0 0],[63 68],'Color','Blue','LineStyle','-','LineWidth',2);