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’);

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;

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;

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

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;

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_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);

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

% 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

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

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_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)

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;

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);

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)

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;

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);

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_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