一、理论基础

热传导主要基于傅里叶热传导定律。在燃料电池中，除了各组件内部的热传导外，还有冷却流体与双极板的对流换热。公式略。
燃料电池内部稳态导热：
$$\frac{d^{2}T}{d{x}^2}+\frac{q_{int}}{k}=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);