Thermodynamics Homework
function main = ADD_2bed_Sg()
clc; clear all; close all; clf;
global Flag1
global Rw E D_so Rp Ea R Cpads wm_d n_d Qst_d A
global Tcw_in_d T_chw_d_in Thw_in_d m_dot_cw_cond_d m_dot_cw_bed_d m_dot_hw_bed_d m_dot_chw_d
global M_ref_L_cond_d m_dot_sw_in_i m_dot_br_out_i
global T_w_bedA_d_out T_w_bedD_d_out T_w_cond_d_out T_chw_d_out
global MCphx_bed_d MCphx_cond_d MCphx_evp_d
global U_bedA_d U_bedD_d A_bed_d U_cond_d A_cond_d U_evp_d A_evp_d
global Mads_d P_bedA_d P_bedD_d hfg_bedD_d
t_normal = 560; %Time for normal operation
t_sw = 40; %Time for switching
t_cycle = (t_normal+t_sw); %Total cycle time
step = 5; %Evaluate every 5 seconds
N_cyc = 12; %No. of cycles evaluated
N_cyc_day = (24*60*60)/t_cycle;
Tcw_in_d = 20+273; %Kelvin
Thw_in_d = 70+273; %Kelvin
T_chw_d_in = 14.8+273;
%Constant parameters
D_so = 2.54E-4; % m2/s
Rp = 1E-4; %m
Ea = 45.5; %Kj/mol
R = 8.3144621e-3; %kJ/mol.K
Cpads = 0.805; %Cp of silica gel (kJ/kg.K)
Qst_d = 3110; %kJ/kg
%A = 9E-7;
Rw = 8.3144621e-3; %kj/kg-k
E = 4; %(1/Pa)
wm_d = 0.21;
n_d = 5;
Mads_d = 20; %kg
%=========================================================================%
%Heat Exchanger Parameters
%==========================================================================
U_bedA_d = 0.175; % kW/m2.K
U_bedD_d = 0.175; % kW/m2.K
A_bed_d = 58.259; % m2
MCphx_bed_d = 140; % Kg*Kj/Kg-k
m_dot_cw_bed_d = 1.19; %kg/s
m_dot_hw_bed_d = 1.20; %Hot water mass flow rate (kg/s)
U_cond_d = 4.115;
A_cond_d = 3.71;
MCphx_cond_d = 18.61; % Kg*Kj/Kg-k
m_dot_cw_cond_d = 1.0;
M_ref_L_cond_d = 10;
U_evp_d = 1.715;
A_evp_d = 3.5;
m_dot_chw_d = 1.0;
MCphx_evp_d = 25.44; % Kg*Kj/Kg-k
M_sw_evp_d_i = 250;
m_dot_sw_in_i = 0.065; %kg/s
m_dot_br_out_i = 0.03; %kg/s
%==========================================================================
%Initial Values
%==========================================================================
Time_i = 0;
T_cond_d_i = 30 + 273;
T_evp_d_i = 30 + 273;
T_bedA_d_i = 30+273; %Kelvin
Ps_bedA_d_i = 10^3*(-5.007261711E-15*T_bedA_d_i^6 + 1.249032418E-11*T_bedA_d_i^5 – 1.135454979E-08*T_bedA_d_i^4 + 5.164337223E-06*T_bedA_d_i^3 – 1.275119405E-03*T_bedA_d_i^2 + 1.643475641E-01*T_bedA_d_i – 8.706432224E+00);
P_bedA_d_i = 10^3*(-5.007261711E-15*T_evp_d_i^6 + 1.249032418E-11*T_evp_d_i^5 – 1.135454979E-08*T_evp_d_i^4 + 5.164337223E-06*T_evp_d_i^3 – 1.275119405E-03*T_evp_d_i^2 + 1.643475641E-01*T_evp_d_i – 8.706432224E+00);
if(Ps_bedA_d_i <= P_bedA_d_i)
%w_bedA_d_i = wm_d;
w_bedA_d_i = 0;
else
w_bedA_d_i = wm_d*exp(-(((Rw*T_bedA_d_i)/E)*log(Ps_bedA_d_i/P_bedA_d_i))^n_d);
end
T_bedD_d_i = 30+273;
Ps_bedD_d_i = 10^3*(-5.007261711E-15*T_bedD_d_i^6 + 1.249032418E-11*T_bedD_d_i^5 – 1.135454979E-08*T_bedD_d_i^4 + 5.164337223E-06*T_bedD_d_i^3 – 1.275119405E-03*T_bedD_d_i^2 + 1.643475641E-01*T_bedD_d_i – 8.706432224E+00);
P_bedD_d_i = 10^3*(-5.007261711E-15*T_cond_d_i^6 + 1.249032418E-11*T_cond_d_i^5 – 1.135454979E-08*T_cond_d_i^4 + 5.164337223E-06*T_cond_d_i^3 – 1.275119405E-03*T_cond_d_i^2 + 1.643475641E-01*T_cond_d_i – 8.706432224E+00);
if(Ps_bedD_d_i <= P_bedD_d_i)
w_bedD_d_i = wm_d;
else
w_bedD_d_i = wm_d*exp(-(((Rw*T_bedD_d_i)/E)*log(Ps_bedD_d_i/P_bedD_d_i))^n_d);
end
%==========================================================================
t = Time_i;
s = 0;
for M =1:N_cyc %Main loop for evaluating cycles
t_round = 0;
Q_A_d = 0;
Q_D_d = 0;
Q_C_d = 0;
Q_E_d = 0;
SDWP_A_i = 0;
PR_A = 0;
for L = 1:step:t_cycle
s = s+1; %s is a counter. add counter
t = t+step; %count elapsed time
Time(s) = t; %Time is an array containing t where the process is evaluated at
t_round = t_round + step; %count elapsed time in one for-next loop
if(t_round <= t_normal) %check whether the elapsed time has passed the normal operation time
Flag1 = 1 ; %adsorption/desorption
elseif((t_round > t_normal)&&(t_round <= t_cycle))
Flag1 = 0 ; %switching mode
end;
Timerange1 = [Time_i t]; %set time range to be used to solve differential eqns.
Initialvalues = [w_bedA_d_i T_bedA_d_i w_bedD_d_i T_bedD_d_i T_cond_d_i M_sw_evp_d_i T_evp_d_i]; %Initial values being used to solve diff. eqns. Being stored first to Y before it is being used.
option1 = odeset(‘RelTol’, 1E-4, ‘AbsTol’, 1E-4); %setting error tolerance
Y = ode45(@ddydwdt, Timerange1, Initialvalues, option1); %solving all the diff eqns. in ddydwdt.m file –> the step is being decided by matlab
Y_t = deval(Y,t); %Y_t’s values are only Y’s value at t; not all the values of Y evaluated
%Storing the newly calculated value into array of variables
w_bedA_d(s) = Y_t(1);
T_bedA_d(s) = Y_t(2);
w_bedD_d(s) = Y_t(3);
T_bedD_d(s) = Y_t(4);
T_cond_d(s) = Y_t(5);
M_sw_evp_d(s) = Y_t(6);
T_evp_d(s) = Y_t(7);
Tcw_bed_d_in = Tcw_in_d;
Thw_bed_d_in = Thw_in_d;
if(Flag1 == 1)
Tcw_bed_d_out = T_w_bedA_d_out;
Thw_bed_d_out = T_w_bedD_d_out;
elseif(Flag1 == 0)
Tcw_bed_d_out = T_w_bedD_d_out;
Thw_bed_d_out = T_w_bedA_d_out;
end
T_cw_cond_d_in = Tcw_in_d;
T_cw_cond_d_in_plot(s) = Tcw_in_d;
T_cw_cond_d_out = T_w_cond_d_out;
T_cw_cond_d_out_plot(s) = T_w_cond_d_out;
T_chw_evp_d_in = T_chw_d_in;
T_chw_evp_d_out = T_chw_d_out;
%==================================================================
%Adsorption/Desorption Switching %==================================================================
if ((mod(M,2)==1)||(M == 1))
T_bed2_d(s) = T_bedD_d(s);
W_bed2_d(s) = w_bedD_d(s);
T_bed1_d(s) = T_bedA_d(s);
W_bed1_d(s) = w_bedA_d(s);
P_bed1_d(s) = P_bedA_d;
P_bed2_d(s) = P_bedD_d;
else
T_bed1_d(s) = T_bedD_d(s);
W_bed1_d(s) = w_bedD_d(s);
T_bed2_d(s) = T_bedA_d(s);
W_bed2_d(s) = w_bedA_d(s);
P_bed2_d(s) = P_bedA_d;
P_bed1_d(s) = P_bedD_d;
end
Hcw_bed_d_out = 1E-3*10^3*(1.95664E-10*Tcw_bed_d_out^5 – 3.36772E-07*Tcw_bed_d_out^4 + 2.34408E-04*Tcw_bed_d_out^3 – 8.21171E-02*Tcw_bed_d_out^2 + 1.86089E+01*Tcw_bed_d_out – 2.15617E+03);
Hcw_bed_d_in = 1E-3*10^3*(1.95664E-10*Tcw_bed_d_in^5 – 3.36772E-07*Tcw_bed_d_in^4 + 2.34408E-04*Tcw_bed_d_in^3 – 8.21171E-02*Tcw_bed_d_in^2 + 1.86089E+01*Tcw_bed_d_in – 2.15617E+03);
Hhw_bed_d_in = 1E-3*10^3*(1.95664E-10*Thw_bed_d_in^5 – 3.36772E-07*Thw_bed_d_in^4 + 2.34408E-04*Thw_bed_d_in^3 – 8.21171E-02*Thw_bed_d_in^2 + 1.86089E+01*Thw_bed_d_in – 2.15617E+03);
Hhw_bed_d_out = 1E-3*10^3*(1.95664E-10*Thw_bed_d_out^5 – 3.36772E-07*Thw_bed_d_out^4 + 2.34408E-04*Thw_bed_d_out^3 – 8.21171E-02*Thw_bed_d_out^2 + 1.86089E+01*Thw_bed_d_out – 2.15617E+03);
H_cw_cond_d_out = 1E-3*10^3*(1.95664E-10*T_cw_cond_d_out^5 – 3.36772E-07*T_cw_cond_d_out^4 + 2.34408E-04*T_cw_cond_d_out^3 – 8.21171E-02*T_cw_cond_d_out^2 + 1.86089E+01*T_cw_cond_d_out – 2.15617E+03);
H_cw_cond_d_in = 1E-3*10^3*(1.95664E-10*T_cw_cond_d_in^5 – 3.36772E-07*T_cw_cond_d_in^4 + 2.34408E-04*T_cw_cond_d_in^3 – 8.21171E-02*T_cw_cond_d_in^2 + 1.86089E+01*T_cw_cond_d_in – 2.15617E+03);
Hfg_cond_d =(10^3*(6.40224E-12*Y_t(5)^5 – 3.93212E-08*Y_t(5)^4 + 3.63907E-05*Y_t(5)^3 – 1.39601E-02*Y_t(5)^2 + 4.34707E+00*Y_t(5) + 1.82254E+03)- 10^3*(1.95664E-10*Y_t(5)^5 – 3.36772E-07*Y_t(5)^4 + 2.34408E-04*Y_t(5)^3 – 8.21171E-02*Y_t(5)^2 + 1.86089E+01*Y_t(5) – 2.15617E+03))*1E-3;
H_chw_d_in = 1E-3*10^3*(1.95664E-10*T_chw_evp_d_in^5 – 3.36772E-07*T_chw_evp_d_in^4 + 2.34408E-04*T_chw_evp_d_in^3 – 8.21171E-02*T_chw_evp_d_in^2 + 1.86089E+01*T_chw_evp_d_in – 2.15617E+03);
H_chw_d_out = 1E-3*10^3*(1.95664E-10*T_chw_evp_d_out^5 – 3.36772E-07*T_chw_evp_d_out^4 + 2.34408E-04*T_chw_evp_d_out^3 – 8.21171E-02*T_chw_evp_d_out^2 + 1.86089E+01*T_chw_evp_d_out – 2.15617E+03);
dq_A_d = (m_dot_cw_bed_d*step)*(Hcw_bed_d_out – Hcw_bed_d_in)/t_cycle;
Q_A_d = Q_A_d + dq_A_d;
dq_D_d = (m_dot_hw_bed_d*step)*(Hhw_bed_d_in – Hhw_bed_d_out)/t_cycle;
Q_D_d = Q_D_d + dq_D_d;
dq_C_d = (m_dot_cw_cond_d*step)*(H_cw_cond_d_out – H_cw_cond_d_in)/t_cycle;
Q_C_d = Q_C_d + dq_C_d;
dq_E_d = (m_dot_chw_d*step)*(H_chw_d_in – H_chw_d_out)/t_cycle;
Q_E_d = Q_E_d + dq_E_d;
W_Pr_d(s) = (m_dot_cw_cond_d*(H_cw_cond_d_out – H_cw_cond_d_in)/Hfg_cond_d)*60;
SDWP_A_sum = (SDWP_A_i + (m_dot_cw_cond_d*step*(H_cw_cond_d_out – H_cw_cond_d_in)*N_cyc_day/(2*Mads_d*Hfg_cond_d)));
SDWP_A_i = SDWP_A_sum;
dPR_A = ((((m_dot_cw_cond_d*(H_cw_cond_d_out – H_cw_cond_d_in)/Hfg_cond_d)*Hfg_cond_d)/(m_dot_hw_bed_d*(Hhw_bed_d_in – Hhw_bed_d_out)))*step)/t_cycle;
PR_A = PR_A + dPR_A;
t_xls = [Time’, (T_bed1_d-273)’, (T_bed2_d-273)’, (T_cond_d-273)’, (T_evp_d-273)’, W_bed1_d’, W_bed2_d’, P_bed1_d’, P_bed2_d’, W_Pr_d’];
%======================================================================
%Storing the newly calculated value as the initial values for the next round of calculation %======================================================================
Time_i = t;
w_bedA_d_i = Y_t(1);
T_bedA_d_i = Y_t(2);
w_bedD_d_i = Y_t(3);
T_bedD_d_i = Y_t(4);
T_cond_d_i = Y_t(5);
M_sw_evp_d_i = Y_t(6);
T_evp_d_i = Y_t(7);
end
Q_ads_d(M) = Q_A_d;
Q_des_d(M) = Q_D_d;
Q_cond_d(M)= Q_C_d;
Q_evp_d(M) = Q_E_d;
Q_add_d(M) = (Q_evp_d(M) + Q_des_d(M));
Q_rej_d(M) = (Q_ads_d(M) + Q_cond_d(M));
Error_d(M) = (Q_add_d(M) – Q_rej_d(M))*100/Q_rej_d(M);
COP_d(M) = Q_evp_d(M)/Q_des_d(M)
%SCP_d(M) = Q_evp_d(M)/(2*Mads_d)
SCP_d(M) = Q_evp_d(M)*1000/(3.517*2*Mads_d)
%SCP_d(M) = (Q_evp_d(M)/(2*Mads_d))*(1000/3.517)
SDWP_A_d(M) = SDWP_A_sum
PR_A_d(M) = PR_A
OCR(M)=(Q_evp_d(M)+ Q_cond_d(M))/Q_des_d(M)
cycle(M) = M;
C_xls = [cycle’, COP_d’, SCP_d’, Q_evp_d’, Q_des_d’, SDWP_A_d’, PR_A_d’,OCR’];
w_bedA_d_i = Y_t(3);
T_bedA_d_i = Y_t(4);
w_bedD_d_i = Y_t(1);
T_bedD_d_i = Y_t(2);
end
%xlswrite(‘C:\Users\Administrator\Dropbox\Shetu\NTU\Paper Plan\Multi bed multi stage adsorption cooling cum desalination system\Submission-2\Write up\SDWP_Changed\Simulation\2bed Desalination\Data_2B1E’,C_xls,’section I’, ‘A2’);
%xlswrite(‘C:\Users\Administrator\Dropbox\Shetu\NTU\Paper Plan\Multi bed multi stage adsorption cooling cum desalination system\Submission-2\Write up\SDWP_Changed\Simulation\2bed Desalination\Data_2B1E’,t_xls,’section II’, ‘A2’);
figure (1)
plot(Time, (T_bed1_d-273), Time, (T_bed2_d-273), Time, (T_cond_d-273), Time, (T_evp_d -273));
hold all;
figure (2)
plot(Time, (P_bed1_d), Time, (P_bed2_d));
hold all;
figure(3)
plot(Time, W_bed1_d, Time, W_bed2_d);
hold all;
figure(4)
plot(Time, M_sw_evp_d);
hold all;
figure(5)
plot(Time, W_Pr_d);
hold all;
figure(6)
plot(Time, T_cw_cond_d_in_plot);
hold all;
figure(7)
plot(Time, T_cw_cond_d_out_plot);
hold all;
figure(8)
plot(cycle, Q_add_d, cycle,Q_rej_d);
hold all;
figure(9)
plot(cycle, COP_d);
hold all;
figure(10)
plot(cycle, SCP_d);
hold all;
figure(11)
plot(cycle, SDWP_A_d);
hold all;
figure(12)
plot(cycle, PR_A_d);
hold all;
figure(13)
plot(cycle, Q_evp_d);
hold all;
figure(14)
plot(cycle, Q_cond_d);
hold all;
figure(15)
plot(cycle, Q_des_d);
hold all;
end
function dy = ddydwdt(t, Y)
global Flag1
global Rw E D_so Rp Ea R Cpads wm_d n_d Qst_d A
global m_dot_sw_in_i m_dot_br_out_i
global Tcw_in_d T_chw_d_in Thw_in_d m_dot_cw_bed_d m_dot_hw_bed_d m_dot_cw_cond_d m_dot_chw_d
global T_w_bedA_d_out T_w_bedD_d_out T_w_cond_d_out T_chw_d_out
global MCphx_bed_d MCphx_cond_d MCphx_evp_d M_ref_L_cond_d
global U_bedA_d U_bedD_d A_bed_d U_cond_d A_cond_d U_evp_d A_evp_d
global Mads_d P_bedA_d P_bedD_d hfg_bedD_d
dy = zeros(7,1); %dy is a 6×1 matrix; syntax = zeros(row, column)
%Storing initial values from Main function to local variables
w_bedA_d = Y(1); T_bedA_d = Y(2); w_bedD_d = Y(3); T_bedD_d = Y(4);
T_evp_d = Y(7); T_cond_d = Y(5); M_sw_evp_d = Y(6);
if (M_sw_evp_d >= 252)
m_dot_sw_in = 0;
m_dot_br_out = m_dot_br_out_i;
elseif (M_sw_evp_d <= 248)
m_dot_sw_in = m_dot_sw_in_i;
m_dot_br_out = 0;
else
m_dot_sw_in = m_dot_sw_in_i;
m_dot_br_out = m_dot_br_out_i;
end
%==========================================================================
%Desalination Cycle
%==============================Evaporator==================================
T_chw_d_out = EvpHex(T_evp_d, T_chw_d_in, m_dot_chw_d, U_evp_d, A_evp_d);
h_fg_evp_d =(10^3*(6.40224E-12*T_evp_d^5 – 3.93212E-08*T_evp_d^4 + 3.63907E-05*T_evp_d^3 – 1.39601E-02*T_evp_d^2 + 4.34707E+00*T_evp_d + 1.82254E+03)- 10^3*(1.95664E-10*T_evp_d^5 – 3.36772E-07*T_evp_d^4 + 2.34408E-04*T_evp_d^3 – 8.21171E-02*T_evp_d^2 + 1.86089E+01*T_evp_d – 2.15617E+03))*1E-3;
h_f_chw_d_in = 1E-3*10^3*(1.95664E-10*T_chw_d_in^5 – 3.36772E-07*T_chw_d_in^4 + 2.34408E-04*T_chw_d_in^3 – 8.21171E-02*T_chw_d_in^2 + 1.86089E+01*T_chw_d_in – 2.15617E+03);
h_f_chw_d_out = 1E-3*10^3*(1.95664E-10*T_chw_d_out^5 – 3.36772E-07*T_chw_d_out^4 + 2.34408E-04*T_chw_d_out^3 – 8.21171E-02*T_chw_d_out^2 + 1.86089E+01*T_chw_d_out- 2.15617E+03);
Cp_ref_L_evp_d = 10^3*(-1.35090E-11*T_evp_d^5 + 2.40993E-08*T_evp_d^4 – 1.71323E-05*T_evp_d^3 + 6.07609E-03*T_evp_d^2 – 1.07591*T_evp_d + 8.02982E01)*1E-3;
h_ref_L_evp_d = 1E-3*10^3*(1.95664E-10*T_evp_d^5 – 3.36772E-07*T_evp_d^4 + 2.34408E-04*T_evp_d^3 – 8.21171E-02*T_evp_d^2 + 1.86089E+01*T_evp_d – 2.15617E+03);
%==============================Bed A======================================%
%checking the current mode; normal mode or switching mode
if(Flag1 == 1) %Bed1 is under adsorption mode
T_w_bedA_d_in = Tcw_in_d;
m_dot_w_bedA_d = m_dot_cw_bed_d;
elseif(Flag1 == 0) %Bed1 is under switching mode
T_w_bedA_d_in = Thw_in_d; %Pre-heating for the adsorption process
m_dot_w_bedA_d = m_dot_hw_bed_d;
end
if(Flag1 == 1)
P_bedA_d = 10^3*(-5.007261711E-15*T_evp_d^6 + 1.249032418E-11*T_evp_d^5 – 1.135454979E-08*T_evp_d^4 + 5.164337223E-06*T_evp_d^3 – 1.275119405E-03*T_evp_d^2 + 1.643475641E-01*T_evp_d – 8.706432224E+00);
Ps_bedA_d = 10^3*(-5.007261711E-15*T_bedA_d^6 + 1.249032418E-11*T_bedA_d^5 – 1.135454979E-08*T_bedA_d^4 + 5.164337223E-06*T_bedA_d^3 – 1.275119405E-03*T_bedA_d^2 + 1.643475641E-01*T_bedA_d – 8.706432224E+00);
hfg_bedA_d = (10^3*(6.40224E-12*T_bedA_d^5 – 3.93212E-08*T_bedA_d^4 + 3.63907E-05*T_bedA_d^3 – 1.39601E-02*T_bedA_d^2 + 4.34707E+00*T_bedA_d + 1.82254E+03)- 10^3*(1.95664E-10*T_bedA_d^5 – 3.36772E-07*T_bedA_d^4 + 2.34408E-04*T_bedA_d^3 – 8.21171E-02*T_bedA_d^2 + 1.86089E+01*T_bedA_d – 2.15617E+03))*1E-3;
if(Ps_bedA_d <= P_bedA_d)
w_star_bedA_d = wm_d;
else
w_star_bedA_d = wm_d*exp(-(((Rw*T_bedA_d)/E)*log(Ps_bedA_d/P_bedA_d))^n_d);
end;
w_bedA_const_d= 15*D_so/Rp^2;
elseif(Flag1 == 0)
w_star_bedA_d = w_bedA_d;
theta_bedA_d = w_star_bedA_d/wm_d;
Ps_bedA_d = 10^3*(-5.007261711E-15*T_bedA_d^6 + 1.249032418E-11*T_bedA_d^5 – 1.135454979E-08*T_bedA_d^4 + 5.164337223E-06*T_bedA_d^3 – 1.275119405E-03*T_bedA_d^2 + 1.643475641E-01*T_bedA_d – 8.706432224E+00);
P_bedA_d = Ps_bedA_d/exp((E*(-log(theta_bedA_d))^(1/n_d))/(R*T_bedA_d));
w_bedA_const_d = 0;
end
T_w_bedA_d_out = BedHex(T_bedA_d, T_w_bedA_d_in, m_dot_w_bedA_d, U_bedA_d, A_bed_d);
h_w_bedA_d_in = 1E-3*10^3*(1.95664E-10*T_w_bedA_d_in^5 – 3.36772E-07*T_w_bedA_d_in^4 + 2.34408E-04*T_w_bedA_d_in^3 – 8.21171E-02*T_w_bedA_d_in^2 + 1.86089E+01*T_w_bedA_d_in – 2.15617E+03);
h_w_bedA_d_out = 1E-3*10^3*(1.95664E-10*T_w_bedA_d_out^5 – 3.36772E-07*T_w_bedA_d_out^4 + 2.34408E-04*T_w_bedA_d_out^3 – 8.21171E-02*T_w_bedA_d_out^2 + 1.86089E+01*T_w_bedA_d_out – 2.15617E+03);
CpabeA_d = 10^3*(-1.35090E-11*T_bedA_d^5 + 2.40993E-08*T_bedA_d^4 – 1.71323E-05*T_bedA_d^3 + 6.07609E-03*T_bedA_d^2 – 1.07591*T_bedA_d + 8.02982E01)*1E-3;
%==============================Bed D======================================%
%checking the current mode; normal mode or switching mode
if(Flag1 == 1) %Bed2 is under desorption mode
T_w_bedD_d_in = Thw_in_d;
m_dot_w_bedD_d = m_dot_hw_bed_d;
elseif(Flag1 == 0) %Bed2 is under switching mode
T_w_bedD_d_in = Tcw_in_d; %Pre-cooling for the desorption process
m_dot_w_bedD_d = m_dot_cw_bed_d;
end
if(Flag1 == 1)
P_bedD_d = 10^3*(-5.007261711E-15*T_cond_d^6 + 1.249032418E-11*T_cond_d^5 – 1.135454979E-08*T_cond_d^4 + 5.164337223E-06*T_cond_d^3 – 1.275119405E-03*T_cond_d^2 + 1.643475641E-01*T_cond_d – 8.706432224E+00);
Ps_bedD_d = 10^3*(-5.007261711E-15*T_bedD_d^6 + 1.249032418E-11*T_bedD_d^5 – 1.135454979E-08*T_bedD_d^4 + 5.164337223E-06*T_bedD_d^3 – 1.275119405E-03*T_bedD_d^2 + 1.643475641E-01*T_bedD_d – 8.706432224E+00);
hfg_bedD_d = (10^3*(6.40224E-12*T_bedD_d^5 – 3.93212E-08*T_bedD_d^4 + 3.63907E-05*T_bedD_d^3 – 1.39601E-02*T_bedD_d^2 + 4.34707E+00*T_bedD_d + 1.82254E+03)- 10^3*(1.95664E-10*T_bedD_d^5 – 3.36772E-07*T_bedD_d^4 + 2.34408E-04*T_bedD_d^3 – 8.21171E-02*T_bedD_d^2 + 1.86089E+01*T_bedD_d – 2.15617E+03))*1E-3;
if(Ps_bedD_d <= P_bedD_d)
w_star_bedD_d = wm_d;
else
w_star_bedD_d = wm_d*exp(-(((Rw*T_bedD_d)/E)*log(Ps_bedD_d/P_bedD_d))^n_d);
end
w_bedD_const_d= 15*D_so/Rp^2;
elseif(Flag1 == 0)
w_star_bedD_d = w_bedD_d;
theta_bedD_d = w_star_bedD_d/wm_d;
Ps_bedD_d = 10^3*(-5.007261711E-15*T_bedD_d^6 + 1.249032418E-11*T_bedD_d^5 – 1.135454979E-08*T_bedD_d^4 + 5.164337223E-06*T_bedD_d^3 – 1.275119405E-03*T_bedD_d^2 + 1.643475641E-01*T_bedD_d – 8.706432224E+00);
P_bedD_d = Ps_bedD_d/exp((E*(-log(theta_bedD_d))^(1/n_d))/(R*T_bedD_d));
w_bedD_const_d = 0;
end
T_w_bedD_d_out = BedHex(T_bedD_d, T_w_bedD_d_in, m_dot_w_bedD_d, U_bedD_d, A_bed_d);
h_w_bedD_d_in = 1E-3*10^3*(1.95664E-10*T_w_bedD_d_in^5 – 3.36772E-07*T_w_bedD_d_in^4 + 2.34408E-04*T_w_bedD_d_in^3 – 8.21171E-02*T_w_bedD_d_in^2 + 1.86089E+01*T_w_bedD_d_in – 2.15617E+03);
h_w_bedD_d_out = 1E-3*10^3*(1.95664E-10*T_w_bedD_d_out^5 – 3.36772E-07*T_w_bedD_d_out^4 + 2.34408E-04*T_w_bedD_d_out^3 – 8.21171E-02*T_w_bedD_d_out^2 + 1.86089E+01*T_w_bedD_d_out – 2.15617E+03);
CpabeD_d = 10^3*(-1.35090E-11*T_bedD_d^5 + 2.40993E-08*T_bedD_d^4 – 1.71323E-05*T_bedD_d^3 + 6.07609E-03*T_bedD_d^2 – 1.07591*T_bedD_d + 8.02982E01)*1E-3;
%==============================Condenser===================================
T_w_cond_d_in = Tcw_in_d;
m_dot_w_cond_d = m_dot_cw_cond_d;
T_w_cond_d_out = CondHex(T_cond_d, T_w_cond_d_in, m_dot_w_cond_d, U_cond_d, A_cond_d);
h_fg_cond_d =(10^3*(6.40224E-12*T_cond_d^5 – 3.93212E-08*T_cond_d^4 + 3.63907E-05*T_cond_d^3 – 1.39601E-02*T_cond_d^2 + 4.34707E+00*T_cond_d + 1.82254E+03)- 10^3*(1.95664E-10*T_cond_d^5 – 3.36772E-07*T_cond_d^4 + 2.34408E-04*T_cond_d^3 – 8.21171E-02*T_cond_d^2 + 1.86089E+01*T_cond_d – 2.15617E+03))*1E-3;
h_f_cond_d = 1E-3*10^3*(1.95664E-10*T_cond_d^5 – 3.36772E-07*T_cond_d^4 + 2.34408E-04*T_cond_d^3 – 8.21171E-02*T_cond_d^2 + 1.86089E+01*T_cond_d – 2.15617E+03);
h_f_w_cond_d_in = 1E-3*10^3*(1.95664E-10*T_w_cond_d_in^5 – 3.36772E-07*T_w_cond_d_in^4 + 2.34408E-04*T_w_cond_d_in^3 – 8.21171E-02*T_w_cond_d_in^2 + 1.86089E+01*T_w_cond_d_in- 2.15617E+03);
h_f_w_cond_d_out = 1E-3*10^3*(1.95664E-10*T_w_cond_d_out^5 – 3.36772E-07*T_w_cond_d_out^4 + 2.34408E-04*T_w_cond_d_out^3 – 8.21171E-02*T_w_cond_d_out^2 + 1.86089E+01*T_w_cond_d_out- 2.15617E+03);
Cp_ref_L_cond_d = 10^3*(-1.35090E-11*T_cond_d^5 + 2.40993E-08*T_cond_d^4 – 1.71323E-05*T_cond_d^3 + 6.07609E-03*T_cond_d^2 – 1.07591*T_cond_d + 8.02982E01)*1E-3;
%==========================================================================
%Desalination
%==========================================================================%
%Linear Driving Force
dy(1) = (w_bedA_const_d*exp(-Ea/(R*Y(2)))*(w_star_bedA_d – Y(1)));
%dT/dt from Energy Balance equation for BedA
dy(2) = ((Qst_d* Mads_d * dy(1))-(m_dot_w_bedA_d*(h_w_bedA_d_out-h_w_bedA_d_in)))/(Mads_d*Cpads + MCphx_bed_d + Mads_d*Y(1)*CpabeA_d);
%Linear Driving Force
dy(3) = (w_bedD_const_d*exp(-Ea/(R*Y(4)))*(w_star_bedD_d – Y(3)));
%dT/dt from Energy Balance equation for BedD
dy(4) = (((m_dot_w_bedD_d*(h_w_bedD_d_in – h_w_bedD_d_out) + Qst_d* Mads_d * dy(3))))/(Mads_d*Cpads + MCphx_bed_d + Mads_d*Y(3)*CpabeD_d);
%Energy balance for condenser
dy(5) = (-Mads_d*dy(3)*h_fg_cond_d – m_dot_w_cond_d*(h_f_w_cond_d_out – h_f_w_cond_d_in) + Mads_d*dy(3)*h_f_cond_d)/(MCphx_cond_d + M_ref_L_cond_d*Cp_ref_L_cond_d);
%Mass Balance of Evaporator
dy(6) = m_dot_sw_in – m_dot_br_out – Mads_d*dy(1);
%Energy balance for evaporator
dy(7) = (m_dot_chw_d*(h_f_chw_d_in – h_f_chw_d_out) + h_ref_L_evp_d*(m_dot_sw_in – m_dot_br_out) – h_fg_evp_d*Mads_d*dy(1))/(MCphx_evp_d + Y(6)*Cp_ref_L_evp_d);
end
function [ T_w_bed_out] = BedHex(T_bed, T_w_bed_in, m_dot_w_bed, U_bed, A_bed)
Cp_w_bed = 10^3*(-1.35090E-11*T_bed^5 + 2.40993E-08*T_bed^4 – 1.71323E-05*T_bed^3 + 6.07609E-03*T_bed^2 – 1.07591*T_bed + 8.02982E01)*1E-3;
T_w_bed_out = T_bed + (T_w_bed_in – T_bed)*exp((-U_bed*A_bed)/(m_dot_w_bed*Cp_w_bed));
end
function T_w_cond_out = CondHex(T_cond, T_w_cond_in, m_dot_w_cond, U_cond, A_cond)
Cp_w_cond = 10^3*(-1.35090E-11*T_cond^5 + 2.40993E-08*T_cond^4 – 1.71323E-05*T_cond^3 + 6.07609E-03*T_cond^2 – 1.07591*T_cond + 8.02982E01)*1E-3;
T_w_cond_out = T_cond + (T_w_cond_in – T_cond)*exp((-U_cond*A_cond)/(m_dot_w_cond*Cp_w_cond));
end
function T_chw_out = EvpHex(T_evp, T_chw_in, m_dot_chw, U_evp, A_evp)
Cp_w_evp = 10^3*(-1.35090E-11*T_evp^5 + 2.40993E-08*T_evp^4 – 1.71323E-05*T_evp^3 + 6.07609E-03*T_evp^2 – 1.07591*T_evp + 8.02982E01)*1E-3;
T_chw_out = T_evp + (T_chw_in – T_evp)*exp((-U_evp*A_evp)/(m_dot_chw*Cp_w_evp));
end
Solution
function main = ADD_2bed_Sg()
clc; clear all; close all; clf;
global Flag1
global Rw E D_so Rp Ea R Cpads wm_d n_d Qst_d A
global Tcw_in_d T_chw_d_in Thw_in_d m_dot_cw_cond_d m_dot_cw_bed_d m_dot_hw_bed_d m_dot_chw_d
global M_ref_L_cond_d m_dot_sw_in_i m_dot_br_out_i
global T_w_bedA_d_out T_w_bedD_d_out T_w_cond_d_out T_chw_d_out
global MCphx_bed_d MCphx_cond_d MCphx_evp_d
global U_bedA_d U_bedD_d A_bed_d U_cond_d A_cond_d U_evp_d A_evp_d
global Mads_d P_bedA_d P_bedD_d hfg_bedD_d
t_normal_range = 960; %60:960
Thw_in_d_range = (55+273):(90+273);
COPmax = zeros(1,size(Thw_in_d_range,2));
COPmaxtime = zeros(1,size(Thw_in_d_range,2));
SCPmax = zeros(1,size(Thw_in_d_range,2));
SCPmaxtime = zeros(1,size(Thw_in_d_range,2));
SDWPmax = zeros(1,size(Thw_in_d_range,2));
SDWPmaxtime = zeros(1,size(Thw_in_d_range,2));
for i=1:size(Thw_in_d_range,2)
for j=1:size(t_normal_range,2)
t_normal = t_normal_range(j); %Time for normal operation
t_sw = 40; %Time for switching
t_cycle = (t_normal+t_sw); %Total cycle time
step = 5; %Evaluate every 5 seconds
N_cyc = 12; %No. of cycles evaluated
N_cyc_day = (24*60*60)/t_cycle;
Tcw_in_d = 20+273; %Kelvin
Thw_in_d = Thw_in_d_range(i); %Kelvin
T_chw_d_in = 14.8+273;
%Constant parameters
D_so = 2.54E-4; % m2/s
Rp = 1E-4; %m
Ea = 45.5; %Kj/mol
R = 8.3144621e-3; %kJ/mol.K
Cpads = 0.805; %Cp of silica gel (kJ/kg.K)
Qst_d = 3110; %kJ/kg
%A = 9E-7;
Rw = 8.3144621e-3; %kj/kg-k
E = 4; %(1/Pa)
wm_d = 0.21;
n_d = 5;
Mads_d = 20; %kg
%=========================================================================%
%Heat Exchanger Parameters
%==========================================================================
U_bedA_d = 0.175; % kW/m2.K
U_bedD_d = 0.175; % kW/m2.K
A_bed_d = 58.259; % m2
MCphx_bed_d = 140; % Kg*Kj/Kg-k
m_dot_cw_bed_d = 1.19; %kg/s
m_dot_hw_bed_d = 1.20; %Hot water mass flow rate (kg/s)
U_cond_d = 4.115;
A_cond_d = 3.71;
MCphx_cond_d = 18.61; % Kg*Kj/Kg-k
m_dot_cw_cond_d = 1.0;
M_ref_L_cond_d = 10;
U_evp_d = 1.715;
A_evp_d = 3.5;
m_dot_chw_d = 1.0;
MCphx_evp_d = 25.44; % Kg*Kj/Kg-k
M_sw_evp_d_i = 250;
m_dot_sw_in_i = 0.065; %kg/s
m_dot_br_out_i = 0.03; %kg/s
%==========================================================================
%Initial Values
%==========================================================================
Time_i = 0;
T_cond_d_i = 30 + 273;
T_evp_d_i = 30 + 273;
T_bedA_d_i = 30+273; %Kelvin
Ps_bedA_d_i = 10^3*(-5.007261711E-15*T_bedA_d_i^6 + 1.249032418E-11*T_bedA_d_i^5 – 1.135454979E-08*T_bedA_d_i^4 + 5.164337223E-06*T_bedA_d_i^3 – 1.275119405E-03*T_bedA_d_i^2 + 1.643475641E-01*T_bedA_d_i – 8.706432224E+00);
P_bedA_d_i = 10^3*(-5.007261711E-15*T_evp_d_i^6 + 1.249032418E-11*T_evp_d_i^5 – 1.135454979E-08*T_evp_d_i^4 + 5.164337223E-06*T_evp_d_i^3 – 1.275119405E-03*T_evp_d_i^2 + 1.643475641E-01*T_evp_d_i – 8.706432224E+00);
if(Ps_bedA_d_i <= P_bedA_d_i)
%w_bedA_d_i = wm_d;
w_bedA_d_i = 0;
else
w_bedA_d_i = wm_d*exp(-(((Rw*T_bedA_d_i)/E)*log(Ps_bedA_d_i/P_bedA_d_i))^n_d);
end
T_bedD_d_i = 30+273;
Ps_bedD_d_i = 10^3*(-5.007261711E-15*T_bedD_d_i^6 + 1.249032418E-11*T_bedD_d_i^5 – 1.135454979E-08*T_bedD_d_i^4 + 5.164337223E-06*T_bedD_d_i^3 – 1.275119405E-03*T_bedD_d_i^2 + 1.643475641E-01*T_bedD_d_i – 8.706432224E+00);
P_bedD_d_i = 10^3*(-5.007261711E-15*T_cond_d_i^6 + 1.249032418E-11*T_cond_d_i^5 – 1.135454979E-08*T_cond_d_i^4 + 5.164337223E-06*T_cond_d_i^3 – 1.275119405E-03*T_cond_d_i^2 + 1.643475641E-01*T_cond_d_i – 8.706432224E+00);
if(Ps_bedD_d_i <= P_bedD_d_i)
w_bedD_d_i = wm_d;
else
w_bedD_d_i = wm_d*exp(-(((Rw*T_bedD_d_i)/E)*log(Ps_bedD_d_i/P_bedD_d_i))^n_d);
end
%==========================================================================
t = Time_i;
s = 0;
for M =1:N_cyc %Main loop for evaluating cycles
t_round = 0;
Q_A_d = 0;
Q_D_d = 0;
Q_C_d = 0;
Q_E_d = 0;
SDWP_A_i = 0;
PR_A = 0;
for L = 1:step:t_cycle
s = s+1; %s is a counter. add counter
t = t+step; %count elapsed time
Time(s) = t; %Time is an array containing t where the process is evaluated at
t_round = t_round + step; %count elapsed time in one for-next loop
if(t_round <= t_normal) %check whether the elapsed time has passed the normal operation time
Flag1 = 1 ; %adsorption/desorption
elseif((t_round > t_normal)&&(t_round <= t_cycle))
Flag1 = 0 ; %switching mode
end;
Timerange1 = [Time_i t]; %set time range to be used to solve differential eqns.
Initialvalues = [w_bedA_d_i T_bedA_d_i w_bedD_d_i T_bedD_d_i T_cond_d_i M_sw_evp_d_i T_evp_d_i]; %Initial values being used to solve diff. eqns. Being stored first to Y before it is being used.
option1 = odeset(‘RelTol’, 1E-4, ‘AbsTol’, 1E-4); %setting error tolerance
Y = ode45(@ddydwdt, Timerange1, Initialvalues, option1); %solving all the diff eqns. in ddydwdt.m file –> the step is being decided by matlab
Y_t = deval(Y,t); %Y_t’s values are only Y’s value at t; not all the values of Y evaluated
%Storing the newly calculated value into array of variables
w_bedA_d(s) = Y_t(1);
T_bedA_d(s) = Y_t(2);
w_bedD_d(s) = Y_t(3);
T_bedD_d(s) = Y_t(4);
T_cond_d(s) = Y_t(5);
M_sw_evp_d(s) = Y_t(6);
T_evp_d(s) = Y_t(7);
Tcw_bed_d_in = Tcw_in_d;
Thw_bed_d_in = Thw_in_d;
if(Flag1 == 1)
Tcw_bed_d_out = T_w_bedA_d_out;
Thw_bed_d_out = T_w_bedD_d_out;
elseif(Flag1 == 0)
Tcw_bed_d_out = T_w_bedD_d_out;
Thw_bed_d_out = T_w_bedA_d_out;
end
T_cw_cond_d_in = Tcw_in_d;
T_cw_cond_d_in_plot(s) = Tcw_in_d;
T_cw_cond_d_out = T_w_cond_d_out;
T_cw_cond_d_out_plot(s) = T_w_cond_d_out;
T_chw_evp_d_in = T_chw_d_in;
T_chw_evp_d_out = T_chw_d_out;
%==================================================================
%Adsorption/Desorption Switching %==================================================================
if ((mod(M,2)==1)||(M == 1))
T_bed2_d(s) = T_bedD_d(s);
W_bed2_d(s) = w_bedD_d(s);
T_bed1_d(s) = T_bedA_d(s);
W_bed1_d(s) = w_bedA_d(s);
P_bed1_d(s) = P_bedA_d;
P_bed2_d(s) = P_bedD_d;
else
T_bed1_d(s) = T_bedD_d(s);
W_bed1_d(s) = w_bedD_d(s);
T_bed2_d(s) = T_bedA_d(s);
W_bed2_d(s) = w_bedA_d(s);
P_bed2_d(s) = P_bedA_d;
P_bed1_d(s) = P_bedD_d;
end
Hcw_bed_d_out = 1E-3*10^3*(1.95664E-10*Tcw_bed_d_out^5 – 3.36772E-07*Tcw_bed_d_out^4 + 2.34408E-04*Tcw_bed_d_out^3 – 8.21171E-02*Tcw_bed_d_out^2 + 1.86089E+01*Tcw_bed_d_out – 2.15617E+03);
Hcw_bed_d_in = 1E-3*10^3*(1.95664E-10*Tcw_bed_d_in^5 – 3.36772E-07*Tcw_bed_d_in^4 + 2.34408E-04*Tcw_bed_d_in^3 – 8.21171E-02*Tcw_bed_d_in^2 + 1.86089E+01*Tcw_bed_d_in – 2.15617E+03);
Hhw_bed_d_in = 1E-3*10^3*(1.95664E-10*Thw_bed_d_in^5 – 3.36772E-07*Thw_bed_d_in^4 + 2.34408E-04*Thw_bed_d_in^3 – 8.21171E-02*Thw_bed_d_in^2 + 1.86089E+01*Thw_bed_d_in – 2.15617E+03);
Hhw_bed_d_out = 1E-3*10^3*(1.95664E-10*Thw_bed_d_out^5 – 3.36772E-07*Thw_bed_d_out^4 + 2.34408E-04*Thw_bed_d_out^3 – 8.21171E-02*Thw_bed_d_out^2 + 1.86089E+01*Thw_bed_d_out – 2.15617E+03);
H_cw_cond_d_out = 1E-3*10^3*(1.95664E-10*T_cw_cond_d_out^5 – 3.36772E-07*T_cw_cond_d_out^4 + 2.34408E-04*T_cw_cond_d_out^3 – 8.21171E-02*T_cw_cond_d_out^2 + 1.86089E+01*T_cw_cond_d_out – 2.15617E+03);
H_cw_cond_d_in = 1E-3*10^3*(1.95664E-10*T_cw_cond_d_in^5 – 3.36772E-07*T_cw_cond_d_in^4 + 2.34408E-04*T_cw_cond_d_in^3 – 8.21171E-02*T_cw_cond_d_in^2 + 1.86089E+01*T_cw_cond_d_in – 2.15617E+03);
Hfg_cond_d =(10^3*(6.40224E-12*Y_t(5)^5 – 3.93212E-08*Y_t(5)^4 + 3.63907E-05*Y_t(5)^3 – 1.39601E-02*Y_t(5)^2 + 4.34707E+00*Y_t(5) + 1.82254E+03)- 10^3*(1.95664E-10*Y_t(5)^5 – 3.36772E-07*Y_t(5)^4 + 2.34408E-04*Y_t(5)^3 – 8.21171E-02*Y_t(5)^2 + 1.86089E+01*Y_t(5) – 2.15617E+03))*1E-3;
H_chw_d_in = 1E-3*10^3*(1.95664E-10*T_chw_evp_d_in^5 – 3.36772E-07*T_chw_evp_d_in^4 + 2.34408E-04*T_chw_evp_d_in^3 – 8.21171E-02*T_chw_evp_d_in^2 + 1.86089E+01*T_chw_evp_d_in – 2.15617E+03);
H_chw_d_out = 1E-3*10^3*(1.95664E-10*T_chw_evp_d_out^5 – 3.36772E-07*T_chw_evp_d_out^4 + 2.34408E-04*T_chw_evp_d_out^3 – 8.21171E-02*T_chw_evp_d_out^2 + 1.86089E+01*T_chw_evp_d_out – 2.15617E+03);
dq_A_d = (m_dot_cw_bed_d*step)*(Hcw_bed_d_out – Hcw_bed_d_in)/t_cycle;
Q_A_d = Q_A_d + dq_A_d;
dq_D_d = (m_dot_hw_bed_d*step)*(Hhw_bed_d_in – Hhw_bed_d_out)/t_cycle;
Q_D_d = Q_D_d + dq_D_d;
dq_C_d = (m_dot_cw_cond_d*step)*(H_cw_cond_d_out – H_cw_cond_d_in)/t_cycle;
Q_C_d = Q_C_d + dq_C_d;
dq_E_d = (m_dot_chw_d*step)*(H_chw_d_in – H_chw_d_out)/t_cycle;
Q_E_d = Q_E_d + dq_E_d;
W_Pr_d(s) = (m_dot_cw_cond_d*(H_cw_cond_d_out – H_cw_cond_d_in)/Hfg_cond_d)*60;
SDWP_A_sum = (SDWP_A_i + (m_dot_cw_cond_d*step*(H_cw_cond_d_out – H_cw_cond_d_in)*N_cyc_day/(2*Mads_d*Hfg_cond_d)));
SDWP_A_i = SDWP_A_sum;
dPR_A = ((((m_dot_cw_cond_d*(H_cw_cond_d_out – H_cw_cond_d_in)/Hfg_cond_d)*Hfg_cond_d)/(m_dot_hw_bed_d*(Hhw_bed_d_in – Hhw_bed_d_out)))*step)/t_cycle;
PR_A = PR_A + dPR_A;
t_xls = [Time’, (T_bed1_d-273)’, (T_bed2_d-273)’, (T_cond_d-273)’, (T_evp_d-273)’, W_bed1_d’, W_bed2_d’, P_bed1_d’, P_bed2_d’, W_Pr_d’];
%======================================================================
%Storing the newly calculated value as the initial values for the next round of calculation %======================================================================
Time_i = t;
w_bedA_d_i = Y_t(1);
T_bedA_d_i = Y_t(2);
w_bedD_d_i = Y_t(3);
T_bedD_d_i = Y_t(4);
T_cond_d_i = Y_t(5);
M_sw_evp_d_i = Y_t(6);
T_evp_d_i = Y_t(7);
end
Q_ads_d(M) = Q_A_d;
Q_des_d(M) = Q_D_d;
Q_cond_d(M)= Q_C_d;
Q_evp_d(M) = Q_E_d;
Q_add_d(M) = (Q_evp_d(M) + Q_des_d(M));
Q_rej_d(M) = (Q_ads_d(M) + Q_cond_d(M));
Error_d(M) = (Q_add_d(M) – Q_rej_d(M))*100/Q_rej_d(M);
COP_d(M) = Q_evp_d(M)/Q_des_d(M);
%SCP_d(M) = Q_evp_d(M)/(2*Mads_d)
SCP_d(M) = Q_evp_d(M)*1000/(3.517*2*Mads_d);
%SCP_d(M) = (Q_evp_d(M)/(2*Mads_d))*(1000/3.517)
SDWP_A_d(M) = SDWP_A_sum ;
PR_A_d(M) = PR_A ;
OCR(M)=(Q_evp_d(M)+ Q_cond_d(M))/Q_des_d(M);
cycle(M) = M;
C_xls = [cycle’, COP_d’, SCP_d’, Q_evp_d’, Q_des_d’, SDWP_A_d’, PR_A_d’,OCR’];
w_bedA_d_i = Y_t(3);
T_bedA_d_i = Y_t(4);
w_bedD_d_i = Y_t(1);
T_bedD_d_i = Y_t(2);
end
%xlswrite(‘C:\Users\Administrator\Dropbox\Shetu\NTU\Paper Plan\Multi bed multi stage adsorption cooling cum desalination system\Submission-2\Write up\SDWP_Changed\Simulation\2bed Desalination\Data_2B1E’,C_xls,’section I’, ‘A2’);
%xlswrite(‘C:\Users\Administrator\Dropbox\Shetu\NTU\Paper Plan\Multi bed multi stage adsorption cooling cum desalination system\Submission-2\Write up\SDWP_Changed\Simulation\2bed Desalination\Data_2B1E’,t_xls,’section II’, ‘A2’);
% figure (1)
% plot(Time, (T_bed1_d-273), Time, (T_bed2_d-273), Time, (T_cond_d-273), Time, (T_evp_d -273));
% hold all;
%
% figure (2)
% plot(Time, (P_bed1_d), Time, (P_bed2_d));
% hold all;
%
% figure(3)
% plot(Time, W_bed1_d, Time, W_bed2_d);
% hold all;
%
% figure(4)
% plot(Time, M_sw_evp_d);
% hold all;
%
% figure(5)
% plot(Time, W_Pr_d);
% hold all;
%
% figure(6)
% plot(Time, T_cw_cond_d_in_plot);
% hold all;
%
% figure(7)
% plot(Time, T_cw_cond_d_out_plot);
% hold all;
%
% figure(8)
% plot(cycle, Q_add_d, cycle,Q_rej_d);
% hold all;
%
% figure(9)
% plot(cycle, COP_d);
% hold all;
%
% figure(10)
% plot(cycle, SCP_d);
% hold all;
%
% figure(11)
% plot(cycle, SDWP_A_d);
% hold all;
%
% figure(12)
% plot(cycle, PR_A_d);
% hold all;
%
% figure(13)
% plot(cycle, Q_evp_d);
% hold all;
%
% figure(14)
% plot(cycle, Q_cond_d);
% hold all;
%
% figure(15)
% plot(cycle, Q_des_d);
% hold all;
if COP_d(1,8)>COPmax(1,i)
COPmax(1,i) = COP_d(1,8);
COPmaxtime(1,i) = t_cycle;
end
if SCP_d(1,8)>SCPmax(1,i)
SCPmax(1,i) = SCP_d(1,8);
SCPmaxtime(1,i) = t_cycle;
end
if SDWP_A_d(1,8)>SDWPmax(1,i)
SDWPmax(1,i) = SDWP_A_d(1,8);
SDWPmaxtime(1,i) = t_cycle;
end
end
disp([‘Temp: ‘,num2str(Thw_in_d-273) ,’degree celsius’,’ – Max value of COP is ‘,num2str(COPmax(1,i)),’ obtained at ‘,num2str(COPmaxtime(1,i))])
disp([‘Temp: ‘,num2str(Thw_in_d-273) ,’degree celsius’,’ – Max value of SCP is ‘,num2str(SCPmax(1,i)),’ obtained at ‘,num2str(SCPmaxtime(1,i))])
disp([‘Temp: ‘,num2str(Thw_in_d-273) ,’degree celsius’,’ – Max value of SDWP is ‘,num2str(SDWPmax(1,i)),’ obtained at ‘,num2str(SDWPmaxtime(1,i))])
end
end
function dy = ddydwdt(t, Y)
global Flag1
global Rw E D_so Rp Ea R Cpads wm_d n_d Qst_d A
global m_dot_sw_in_i m_dot_br_out_i
global Tcw_in_d T_chw_d_in Thw_in_d m_dot_cw_bed_d m_dot_hw_bed_d m_dot_cw_cond_d m_dot_chw_d
global T_w_bedA_d_out T_w_bedD_d_out T_w_cond_d_out T_chw_d_out
global MCphx_bed_d MCphx_cond_d MCphx_evp_d M_ref_L_cond_d
global U_bedA_d U_bedD_d A_bed_d U_cond_d A_cond_d U_evp_d A_evp_d
global Mads_d P_bedA_d P_bedD_d hfg_bedD_d
dy = zeros(7,1); %dy is a 6×1 matrix; syntax = zeros(row, column)
%Storing initial values from Main function to local variables
w_bedA_d = Y(1); T_bedA_d = Y(2); w_bedD_d = Y(3); T_bedD_d = Y(4);
T_evp_d = Y(7); T_cond_d = Y(5); M_sw_evp_d = Y(6);
if (M_sw_evp_d >= 252)
m_dot_sw_in = 0;
m_dot_br_out = m_dot_br_out_i;
elseif (M_sw_evp_d <= 248)
m_dot_sw_in = m_dot_sw_in_i;
m_dot_br_out = 0;
else
m_dot_sw_in = m_dot_sw_in_i;
m_dot_br_out = m_dot_br_out_i;
end
%==========================================================================
%Desalination Cycle
%==============================Evaporator==================================
T_chw_d_out = EvpHex(T_evp_d, T_chw_d_in, m_dot_chw_d, U_evp_d, A_evp_d);
h_fg_evp_d =(10^3*(6.40224E-12*T_evp_d^5 – 3.93212E-08*T_evp_d^4 + 3.63907E-05*T_evp_d^3 – 1.39601E-02*T_evp_d^2 + 4.34707E+00*T_evp_d + 1.82254E+03)- 10^3*(1.95664E-10*T_evp_d^5 – 3.36772E-07*T_evp_d^4 + 2.34408E-04*T_evp_d^3 – 8.21171E-02*T_evp_d^2 + 1.86089E+01*T_evp_d – 2.15617E+03))*1E-3;
h_f_chw_d_in = 1E-3*10^3*(1.95664E-10*T_chw_d_in^5 – 3.36772E-07*T_chw_d_in^4 + 2.34408E-04*T_chw_d_in^3 – 8.21171E-02*T_chw_d_in^2 + 1.86089E+01*T_chw_d_in – 2.15617E+03);
h_f_chw_d_out = 1E-3*10^3*(1.95664E-10*T_chw_d_out^5 – 3.36772E-07*T_chw_d_out^4 + 2.34408E-04*T_chw_d_out^3 – 8.21171E-02*T_chw_d_out^2 + 1.86089E+01*T_chw_d_out- 2.15617E+03);
Cp_ref_L_evp_d = 10^3*(-1.35090E-11*T_evp_d^5 + 2.40993E-08*T_evp_d^4 – 1.71323E-05*T_evp_d^3 + 6.07609E-03*T_evp_d^2 – 1.07591*T_evp_d + 8.02982E01)*1E-3;
h_ref_L_evp_d = 1E-3*10^3*(1.95664E-10*T_evp_d^5 – 3.36772E-07*T_evp_d^4 + 2.34408E-04*T_evp_d^3 – 8.21171E-02*T_evp_d^2 + 1.86089E+01*T_evp_d – 2.15617E+03);
%==============================Bed A======================================%
%checking the current mode; normal mode or switching mode
if(Flag1 == 1) %Bed1 is under adsorption mode
T_w_bedA_d_in = Tcw_in_d;
m_dot_w_bedA_d = m_dot_cw_bed_d;
elseif(Flag1 == 0) %Bed1 is under switching mode
T_w_bedA_d_in = Thw_in_d; %Pre-heating for the adsorption process
m_dot_w_bedA_d = m_dot_hw_bed_d;
end
if(Flag1 == 1)
P_bedA_d = 10^3*(-5.007261711E-15*T_evp_d^6 + 1.249032418E-11*T_evp_d^5 – 1.135454979E-08*T_evp_d^4 + 5.164337223E-06*T_evp_d^3 – 1.275119405E-03*T_evp_d^2 + 1.643475641E-01*T_evp_d – 8.706432224E+00);
Ps_bedA_d = 10^3*(-5.007261711E-15*T_bedA_d^6 + 1.249032418E-11*T_bedA_d^5 – 1.135454979E-08*T_bedA_d^4 + 5.164337223E-06*T_bedA_d^3 – 1.275119405E-03*T_bedA_d^2 + 1.643475641E-01*T_bedA_d – 8.706432224E+00);
hfg_bedA_d = (10^3*(6.40224E-12*T_bedA_d^5 – 3.93212E-08*T_bedA_d^4 + 3.63907E-05*T_bedA_d^3 – 1.39601E-02*T_bedA_d^2 + 4.34707E+00*T_bedA_d + 1.82254E+03)- 10^3*(1.95664E-10*T_bedA_d^5 – 3.36772E-07*T_bedA_d^4 + 2.34408E-04*T_bedA_d^3 – 8.21171E-02*T_bedA_d^2 + 1.86089E+01*T_bedA_d – 2.15617E+03))*1E-3;
if(Ps_bedA_d <= P_bedA_d)
w_star_bedA_d = wm_d;
else
w_star_bedA_d = wm_d*exp(-(((Rw*T_bedA_d)/E)*log(Ps_bedA_d/P_bedA_d))^n_d);
end;
w_bedA_const_d= 15*D_so/Rp^2;
elseif(Flag1 == 0)
w_star_bedA_d = w_bedA_d;
theta_bedA_d = w_star_bedA_d/wm_d;
Ps_bedA_d = 10^3*(-5.007261711E-15*T_bedA_d^6 + 1.249032418E-11*T_bedA_d^5 – 1.135454979E-08*T_bedA_d^4 + 5.164337223E-06*T_bedA_d^3 – 1.275119405E-03*T_bedA_d^2 + 1.643475641E-01*T_bedA_d – 8.706432224E+00);
P_bedA_d = Ps_bedA_d/exp((E*(-log(theta_bedA_d))^(1/n_d))/(R*T_bedA_d));
w_bedA_const_d = 0;
end
T_w_bedA_d_out = BedHex(T_bedA_d, T_w_bedA_d_in, m_dot_w_bedA_d, U_bedA_d, A_bed_d);
h_w_bedA_d_in = 1E-3*10^3*(1.95664E-10*T_w_bedA_d_in^5 – 3.36772E-07*T_w_bedA_d_in^4 + 2.34408E-04*T_w_bedA_d_in^3 – 8.21171E-02*T_w_bedA_d_in^2 + 1.86089E+01*T_w_bedA_d_in – 2.15617E+03);
h_w_bedA_d_out = 1E-3*10^3*(1.95664E-10*T_w_bedA_d_out^5 – 3.36772E-07*T_w_bedA_d_out^4 + 2.34408E-04*T_w_bedA_d_out^3 – 8.21171E-02*T_w_bedA_d_out^2 + 1.86089E+01*T_w_bedA_d_out – 2.15617E+03);
CpabeA_d = 10^3*(-1.35090E-11*T_bedA_d^5 + 2.40993E-08*T_bedA_d^4 – 1.71323E-05*T_bedA_d^3 + 6.07609E-03*T_bedA_d^2 – 1.07591*T_bedA_d + 8.02982E01)*1E-3;
%==============================Bed D======================================%
%checking the current mode; normal mode or switching mode
if(Flag1 == 1) %Bed2 is under desorption mode
T_w_bedD_d_in = Thw_in_d;
m_dot_w_bedD_d = m_dot_hw_bed_d;
elseif(Flag1 == 0) %Bed2 is under switching mode
T_w_bedD_d_in = Tcw_in_d; %Pre-cooling for the desorption process
m_dot_w_bedD_d = m_dot_cw_bed_d;
end
if(Flag1 == 1)
P_bedD_d = 10^3*(-5.007261711E-15*T_cond_d^6 + 1.249032418E-11*T_cond_d^5 – 1.135454979E-08*T_cond_d^4 + 5.164337223E-06*T_cond_d^3 – 1.275119405E-03*T_cond_d^2 + 1.643475641E-01*T_cond_d – 8.706432224E+00);
Ps_bedD_d = 10^3*(-5.007261711E-15*T_bedD_d^6 + 1.249032418E-11*T_bedD_d^5 – 1.135454979E-08*T_bedD_d^4 + 5.164337223E-06*T_bedD_d^3 – 1.275119405E-03*T_bedD_d^2 + 1.643475641E-01*T_bedD_d – 8.706432224E+00);
hfg_bedD_d = (10^3*(6.40224E-12*T_bedD_d^5 – 3.93212E-08*T_bedD_d^4 + 3.63907E-05*T_bedD_d^3 – 1.39601E-02*T_bedD_d^2 + 4.34707E+00*T_bedD_d + 1.82254E+03)- 10^3*(1.95664E-10*T_bedD_d^5 – 3.36772E-07*T_bedD_d^4 + 2.34408E-04*T_bedD_d^3 – 8.21171E-02*T_bedD_d^2 + 1.86089E+01*T_bedD_d – 2.15617E+03))*1E-3;
if(Ps_bedD_d <= P_bedD_d)
w_star_bedD_d = wm_d;
else
w_star_bedD_d = wm_d*exp(-(((Rw*T_bedD_d)/E)*log(Ps_bedD_d/P_bedD_d))^n_d);
end
w_bedD_const_d= 15*D_so/Rp^2;
elseif(Flag1 == 0)
w_star_bedD_d = w_bedD_d;
theta_bedD_d = w_star_bedD_d/wm_d;
Ps_bedD_d = 10^3*(-5.007261711E-15*T_bedD_d^6 + 1.249032418E-11*T_bedD_d^5 – 1.135454979E-08*T_bedD_d^4 + 5.164337223E-06*T_bedD_d^3 – 1.275119405E-03*T_bedD_d^2 + 1.643475641E-01*T_bedD_d – 8.706432224E+00);
P_bedD_d = Ps_bedD_d/exp((E*(-log(theta_bedD_d))^(1/n_d))/(R*T_bedD_d));
w_bedD_const_d = 0;
end
T_w_bedD_d_out = BedHex(T_bedD_d, T_w_bedD_d_in, m_dot_w_bedD_d, U_bedD_d, A_bed_d);
h_w_bedD_d_in = 1E-3*10^3*(1.95664E-10*T_w_bedD_d_in^5 – 3.36772E-07*T_w_bedD_d_in^4 + 2.34408E-04*T_w_bedD_d_in^3 – 8.21171E-02*T_w_bedD_d_in^2 + 1.86089E+01*T_w_bedD_d_in – 2.15617E+03);
h_w_bedD_d_out = 1E-3*10^3*(1.95664E-10*T_w_bedD_d_out^5 – 3.36772E-07*T_w_bedD_d_out^4 + 2.34408E-04*T_w_bedD_d_out^3 – 8.21171E-02*T_w_bedD_d_out^2 + 1.86089E+01*T_w_bedD_d_out – 2.15617E+03);
CpabeD_d = 10^3*(-1.35090E-11*T_bedD_d^5 + 2.40993E-08*T_bedD_d^4 – 1.71323E-05*T_bedD_d^3 + 6.07609E-03*T_bedD_d^2 – 1.07591*T_bedD_d + 8.02982E01)*1E-3;
%==============================Condenser===================================
T_w_cond_d_in = Tcw_in_d;
m_dot_w_cond_d = m_dot_cw_cond_d;
T_w_cond_d_out = CondHex(T_cond_d, T_w_cond_d_in, m_dot_w_cond_d, U_cond_d, A_cond_d);
h_fg_cond_d =(10^3*(6.40224E-12*T_cond_d^5 – 3.93212E-08*T_cond_d^4 + 3.63907E-05*T_cond_d^3 – 1.39601E-02*T_cond_d^2 + 4.34707E+00*T_cond_d + 1.82254E+03)- 10^3*(1.95664E-10*T_cond_d^5 – 3.36772E-07*T_cond_d^4 + 2.34408E-04*T_cond_d^3 – 8.21171E-02*T_cond_d^2 + 1.86089E+01*T_cond_d – 2.15617E+03))*1E-3;
h_f_cond_d = 1E-3*10^3*(1.95664E-10*T_cond_d^5 – 3.36772E-07*T_cond_d^4 + 2.34408E-04*T_cond_d^3 – 8.21171E-02*T_cond_d^2 + 1.86089E+01*T_cond_d – 2.15617E+03);
h_f_w_cond_d_in = 1E-3*10^3*(1.95664E-10*T_w_cond_d_in^5 – 3.36772E-07*T_w_cond_d_in^4 + 2.34408E-04*T_w_cond_d_in^3 – 8.21171E-02*T_w_cond_d_in^2 + 1.86089E+01*T_w_cond_d_in- 2.15617E+03);
h_f_w_cond_d_out = 1E-3*10^3*(1.95664E-10*T_w_cond_d_out^5 – 3.36772E-07*T_w_cond_d_out^4 + 2.34408E-04*T_w_cond_d_out^3 – 8.21171E-02*T_w_cond_d_out^2 + 1.86089E+01*T_w_cond_d_out- 2.15617E+03);
Cp_ref_L_cond_d = 10^3*(-1.35090E-11*T_cond_d^5 + 2.40993E-08*T_cond_d^4 – 1.71323E-05*T_cond_d^3 + 6.07609E-03*T_cond_d^2 – 1.07591*T_cond_d + 8.02982E01)*1E-3;
%==========================================================================
%Desalination
%==========================================================================%
%Linear Driving Force
dy(1) = (w_bedA_const_d*exp(-Ea/(R*Y(2)))*(w_star_bedA_d – Y(1)));
%dT/dt from Energy Balance equation for BedA
dy(2) = ((Qst_d* Mads_d * dy(1))-(m_dot_w_bedA_d*(h_w_bedA_d_out-h_w_bedA_d_in)))/(Mads_d*Cpads + MCphx_bed_d + Mads_d*Y(1)*CpabeA_d);
%Linear Driving Force
dy(3) = (w_bedD_const_d*exp(-Ea/(R*Y(4)))*(w_star_bedD_d – Y(3)));
%dT/dt from Energy Balance equation for BedD
dy(4) = (((m_dot_w_bedD_d*(h_w_bedD_d_in – h_w_bedD_d_out) + Qst_d* Mads_d * dy(3))))/(Mads_d*Cpads + MCphx_bed_d + Mads_d*Y(3)*CpabeD_d);
%Energy balance for condenser
dy(5) = (-Mads_d*dy(3)*h_fg_cond_d – m_dot_w_cond_d*(h_f_w_cond_d_out – h_f_w_cond_d_in) + Mads_d*dy(3)*h_f_cond_d)/(MCphx_cond_d + M_ref_L_cond_d*Cp_ref_L_cond_d);
%Mass Balance of Evaporator
dy(6) = m_dot_sw_in – m_dot_br_out – Mads_d*dy(1);
%Energy balance for evaporator
dy(7) = (m_dot_chw_d*(h_f_chw_d_in – h_f_chw_d_out) + h_ref_L_evp_d*(m_dot_sw_in – m_dot_br_out) – h_fg_evp_d*Mads_d*dy(1))/(MCphx_evp_d + Y(6)*Cp_ref_L_evp_d);
end
function [ T_w_bed_out] = BedHex(T_bed, T_w_bed_in, m_dot_w_bed, U_bed, A_bed)
Cp_w_bed = 10^3*(-1.35090E-11*T_bed^5 + 2.40993E-08*T_bed^4 – 1.71323E-05*T_bed^3 + 6.07609E-03*T_bed^2 – 1.07591*T_bed + 8.02982E01)*1E-3;
T_w_bed_out = T_bed + (T_w_bed_in – T_bed)*exp((-U_bed*A_bed)/(m_dot_w_bed*Cp_w_bed));
end
function T_w_cond_out = CondHex(T_cond, T_w_cond_in, m_dot_w_cond, U_cond, A_cond)
Cp_w_cond = 10^3*(-1.35090E-11*T_cond^5 + 2.40993E-08*T_cond^4 – 1.71323E-05*T_cond^3 + 6.07609E-03*T_cond^2 – 1.07591*T_cond + 8.02982E01)*1E-3;
T_w_cond_out = T_cond + (T_w_cond_in – T_cond)*exp((-U_cond*A_cond)/(m_dot_w_cond*Cp_w_cond));
end
function T_chw_out = EvpHex(T_evp, T_chw_in, m_dot_chw, U_evp, A_evp)
Cp_w_evp = 10^3*(-1.35090E-11*T_evp^5 + 2.40993E-08*T_evp^4 – 1.71323E-05*T_evp^3 + 6.07609E-03*T_evp^2 – 1.07591*T_evp + 8.02982E01)*1E-3;
T_chw_out = T_evp + (T_chw_in – T_evp)*exp((-U_evp*A_evp)/(m_dot_chw*Cp_w_evp));
end