chemheatpump.m

Matlab在化学工程中的应用

function ChemHeatPump
% 可逆气固化学热泵制冷系统的动态模拟--考虑热容随反应进度的变化
clear all

global R
global mt ms Ms Ns Mt Ff Cps1 Cps2 Cpf Cpt n2
global d do h1 h2 e1 e2 d1 d2 Vf Vs
global Rho
global koa kod Ea Ed Ma Md
global hsw hfw hfe Asw Afw Afe Wo
global dHe dHr IAD 
global Ta Pc Tfin   % Ta: 环境温度
global Toc Toh
global wgr rhos Cpgr Mg 

% 输入参数
% --------
% (1)反应器结构参数
d = 0.150;                  % 吸附器内径，m
do = 0.006;                 % 气体扩散器（中心孔）的直径，m
h1 = 0.100;                 % 反应材料的高度，m
h2 = 0.114;                 % 夹套有效高度，m
e1 = 0.0045;                % 壁厚，m
e2 = 0.012;                 % 夹套的流通宽度，m
d1 = d + 2*e1;              % 吸附器外径，m (d1=0.159 m)
d2 = d+2*e1+2*e2;           % 反应器外壁，m (d2=0.183 m)
Vf = pi/4*(d2^2-d1^2)*h2;   % 夹套中导热油的体积，m3
Vs = pi/4*(d^2-do^2)*h1;    % 反应器中盐床的体积，m3

% (2)反应器物性参数
Mt = 25;                    % 反应器质量，kg
Cpt = 50;                   % 反应器热容量，Jkg-1K-1    

% (3)反应（吸附）材料的参数
mt = 1.390;                 % kg,IMPEX的总质量,对于SrCl2IMPEX  
ms = 1037.14;               % g，对于SrCl2/NH3IMPEX
Cps1 = 697;                 % SrCl2.NH3的热容，Jm-3K-1
Cps2 = 1480;                % SrCl2.8NH3的热容，Jm-3K-1
Ms = 158.526;               % 反应盐SrCl2的摩尔质量，g/mol
Mg = 17;                    % 氨气的摩尔质量，g/mol
n2 = 7;                     % 对于SrCl2/NH3    
Ns = ms/Ms;
wgr = 0.30;                  % 可膨胀石墨的重量百分比
Cpgr = 674;                 % 可膨胀石墨的比热,J/kgK
rhos = mt/Vs;               % 反应材料的密度，kg/m3

IAD = 1;

% (4)热力学参数  
dHe = 23366;                % 氨的蒸发潜热，J/mol
dHr = 41431;                % 反应热，J/mol
K = 273.15;
R = 8.3145;

% (5)传热参数及传热面积
hsw = 500;                  % W/(Km2)
hfw = 692;                  % W/(Km2)
hfe = 2;                    % W/(Km2)  
Asw = pi*d*h1;
Afw = pi*d1*h1;
Afe = pi*d2*h2;      

% (6)动力学参数
koa = 0.0190;
Ea = 6921;
Ma = 2.96; 
kod = 0.125;
Ed = 9000;
Md = 3.02;

% (7)载热介质(导热油)物性参数：
Rho = 870;                  % 载热介质（导热油）的密度，kgm-3
Cpf = 2400;                 % 载热介质的热容，Jkg-1K-1

% (8)操作参数
Ff = 80e-6;                 % 载热介质的体积流量，m3s-1
Wo = Vf*Rho;                % 夹套中导热油的重量，kg  
Ta = 30+K;                  % 环境温度，K
Th = 140 + K;               % 绝对温度，K       
Toc = 20 + K;               % 冷油温度
Toh = 140 + K;              % 热油温度
Pa = 5;                     % 操作压力，bar
Pd = 13;
 
% 吸附过程模拟
% ------------
tspan = [0  9000];
y0 = [0 Toc Toc Toc];
[t,y] = ode23s(@Equations,tspan,y0,[],1,Pa,Toc);

% 吸附过程总反应进度的动态变化图
plot(t,y(:,1))
xlabel('Time (s)')
ylabel('X_a')
title('吸附过程总反应进度的动态变化')

% 吸附过程床层温度、壁温、载热流体温度的动态变化图
figure
plot(t,y(:,2),'k--',t,y(:,3),'b-',t,y(:,4),'r-.')
xlabel('Time (s)')
ylabel('Temperature (K)')
legend('T_s','T_w','T_f')
title('吸附过程床层温度、壁温、载热流体温度的动态变化')

% 考察吸附压力的影响
Pai = [2, 3, 4, 5];
n = length(Pai);
for i = 1:n
    [t,y] = ode23s(@Equations,tspan,y0,[],1,Pai(i),Toc);
    ti{i} = t;
    Xi{i} = y(:,1);
    Tsi{i} = y(:,2);
end

% 吸附压力对总反应进度X的影响图
figure
plot(ti{1},Xi{1},'r:',ti{2},Xi{2},'b-.',ti{3},Xi{3},'k-',ti{4},Xi{4},'b--')
xlabel('Time (s)')
ylabel('X_a')
legend('Pa = 2','Pa = 3','Pa = 4','Pa = 5')
title('吸附压力对总反应进度的影响')

% 吸附压力对反应床层温度Ts的影响图
figure
plot(ti{1},Tsi{1},'r:',ti{2},Tsi{2},'b-.',ti{3},Tsi{3},'k-',ti{4},Tsi{4},'b--')
xlabel('Time (s)')
ylabel('T_s')
legend('Pa = 2','Pa = 3','Pa = 4','Pa = 5')
title('吸附压力对反应床层温度的影响')


% 解吸过程模拟
% ------------
y0 = [0 Toh Toh Toh];
[t,y] = ode23s(@Equations,tspan,y0,[],-1,Pd,Toh);

% 解附过程图形输出
figure
plot(t,y(:,1))
xlabel('Time (s)')
ylabel('X_d')
title('解吸过程总反应进度的动态变化')
figure
% plot(t,y(:,2:4))
plot(t,y(:,2),'k--',t,y(:,3),'b-',t,y(:,4),'r-.')
xlabel('Time (s)')
ylabel('Temperature (K)')
legend('T_s','T_w','T_f')
title('解吸过程床层温度、壁温、载热流体温度的动态变化')


% ------------------------------------------------------------------
function dydt = Equations(t,y,IAD,Pc,Tfin)
% 函数功能：计算某时刻的反应盐(床层)、载热流体
% 和反应器壁的温度以及反应进度的动态变化
%
% 输入参数：
% t   --- 时间，s
% Cps --- 反应器内IMPEX的热容量，Jm-3K-1
% Cpf --- 载热流体（加热油）的比热容，Jkg-1K-1
% Cpt --- 反应器壁的热容量，Jkg-1K-1 
% Mt  --- 反应器壁的质量，kg
% Ff  --- 载热流体（加热油）的体积流率，m3s-1
% Rho --- 载热介质（导热油）的密度,kgm-3
% dHr --- 反应热
% hsw --- 床层与反应器壁之间的传热参数，W/（Km2）
% Asw --- 床层的传热面积，m2
% hfw --- 反应器壁面和载热流体之间的传热参数，W/K
% hfe --- 载热流体与环境之间的散热参数，W/K 
% Ta  --- 环境温度，K
% Pc  --- 操作压力，bar
%
% 输出变量：
%  Ts  --- 某时刻的反应盐(床层)温度
%  Tf  --- 某时刻的载热流体温度
%  Tw  --- 某时刻的反应器壁温度
%  X   --- 某时刻的总反应进度

global Ns Mt Ff Cps1 Cps2 Cpf Cpt n2
global Rho hsw hfw hfe Asw Afw Afe Wo
global dHr Ta Vf Vs
X = y(1);
Ts = y(2);
Tw = y(3);
Tf = y(4);  

dXdt = Kinetics(IAD,t,X,Pc,Ts);
Cps = CPX(IAD,X);

% 对反应盐(床层) --- 吸附时IAD=1, 解吸时IAD=-1      
dTsdt = ( hsw*Asw*(Tw-Ts)+IAD*n2*Ns*dHr*dXdt )/(Vs*Cps);

% 对反应器壁:
dTwdt = ( hsw*Asw*(Ts-Tw)-hfw*Afw*(Tw-Tf) )/(Mt*Cpt);

% 对载热流体:
dTfdt = ( hfw*Afw*(Tw-Tf)-hfe*Afe*(Tf-Ta) )/(Wo*Cpf)  ...
        - Ff*Rho*Cpf*(Tf-Tfin);

dydt = [dXdt; dTsdt; dTwdt; dTfdt];

% ------------------------------------------------------------------
function dXdt = Kinetics(IAD,t,X,Pc,Ts)
global R koa Ea Ma kod Ed Md
% koa --- 吸附系数
% kod --- 解吸系数
% Ea  --- 吸附活化能
% Ed  --- 解吸活化能
% Ma  --- 吸附幂指数
% Md  --- 解吸幂指数
if IAD == 1;
    ko = koa;
    E = Ea;
    M = Ma;
end
if IAD == -1;
    ko = kod;
    E = Ed;
    M = Md;
end
Ar = ko * exp(-E/R./Ts);
Peq = exp( -4983.16./Ts+16.00 );
Pm = IAD*( Pc-Peq )./Peq;
if (Pm<=0) | (X<0)                      % Pm <= 0表示尚未发生吸附或解吸
    X = 0;
    dXdt = 0;
else
    arg = Ar * Pm .* (M-1) .* t + 1;    
    dXdt = Ar .* ( 1-X ) .^M .* Pm;     % dXdt = f(Ts)    
end

% ------------------------------------------------------------------
function Cp = CPX(IAD,X)                % 计算给定总反应进度下的体积热容
global wgr rhos Cpgr Cps1 Cps2 Mg Ms 
% wgr       weight percentage of inert binder  
% rhos      anhydrous volumetric mass, kg/m3
% Cpgr      specific heat capacity of graphite, J/kgK
% Mg        molar weight of gas, kg/mole 
% Ms        molar weight of anhydrous salt SrCl2, kg/mol
% Cp        J/m3K

if (IAD == 1) 
    rate = X;
else
    rate = 1 - X;
end

Ms1 = rhos*(1-wgr)*(1+1*Mg/Ms)*(1-rate);    % SrCl.NH3
Ms2 = rhos*(1-wgr)*(1+8*Mg/Ms)*rate;        % SrCl.8NH3
Cp = rhos*wgr*Cpgr + Cps1*Ms1 + Cps2*Ms2;