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;