Code doesn't give any output.

1 visualización (últimos 30 días)
TM Abir Ahsan
TM Abir Ahsan el 1 de Feb. de 2022
Comentada: Star Strider el 1 de Feb. de 2022
Ran in:
I kind of coded a energy simulation, but I think I have some syntax issue for which the program is not showing any output. The data is also not getting loaded to workspace. Could any body provide some support?
close all; clear all; clc;
global Tcw_in Thw_in T_chw_in T_chw_out T_w_bedA_out T_w_bedD_out T_w_cond_out P_bedA P_bedD w_bedA w_bedD
global Ea Q R Rp Dso A Ao A1 A2 A3 B Bo B1 B2 B3
global m_dot_cw_bed m_dot_cw_cond m_dot_hw_bed m_dot_chw Cpw Cpv Lv
global MhexCphex_bed Mhex_cond Cphex_cond Mhex_evap Cphex_evap Mads Mw_cond Mw_evap_max Cpads A_bedA U_bedA A_bedD U_bedD A_cond U_cond A_evap U_evap UA_bedA UA_bedD
global t_normal t_sw t_cycle N_cycle step
global FLAG1
%%%%%%%%%%%%%%%%%%%%%%%%%%% Cycle Duration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_normal=500; %Time for normal operation [s]
t_sw=30; %Time for switching operation [s]
t_cycle=t_normal+t_sw; %Total cycle time [s]
step=5; %Evaluate data every 5 sec
N_cycle=20; %No of cycles to be evaluated
%%%%%%%%%%%%%%%%%%%%%%%%%%% Operating conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tcw_in=30; %Temperature of cold water in [C]
Thw_in=85; %Temperature of hot water in [C]
T_chw_in=14; %Temperature of chiled water in [C]
m_dot_cw_bed=1.7; %Temperature of bed cooling water [C]
m_dot_cw_cond=1.3; %Temperature of condenser cooling water [C]
m_dot_hw_bed=1.7; %Temperature of bed heating water [C]
m_dot_chw=0.7; %Temperature of chilled water flow in evaporator [C]
%%%%%%%%%%%%%% Physical and thermal properties of adsorbate-water %%%%%%%%%%%%%
Cpw=4180; %Specific heat of water [J/kgK]
Cpv=1890; %Specific heat of water vapor [J/kgK]
Lv=2500*10^3; %Latent heat of vaporization
%%%%%%%%%%%%%% Physical and thermal properties of adsorbent %%%%%%%%%%%%%%%%
Cpads=924; %Specific heat capacity of silica gel, [J/kg.K]
Mads=47; %Mass of silica, [kg]
Rp=1.6*10^-4; %Radious of silica gel particle, [m]
Dso=2.54*10^-4; %Pre-exponential constant, (m^2 S^-1)
Ea=4.2*10^4; %Activation energy of surface diffusion, (J/mol)
R=8.314; %Ideal gas constant, (j/mol.K)
Q=2.8*10^6; %Heat of adsorption, (J/kg)
%%%%%%%%%%%%%%%%%%%%%%% SBK mode constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ao=-6.5314;
A1=0.072452;
A2=-0.23951*10^-3;
A3=0.25493*10^-6;
Bo=-15.587;
B1=0.15915;
B2=-0.50612*10^-3;
B3=0.5329*10^-6;
%%%%%%%%%%Physical and thermal properties of adsorber/desorber%%%%%%%
MhexCphex_bed=77719.4; %Mass capacitance of heat exchanger (Adsorber/desorber), [J/K]
A_bedA=2.46; %Surface area of adsorber, [m^2]
U_bedA=1602.56; %Heat transfer coefficient of desorber, [W/m^2.K]
UA_bedA=U_bedA*A_bedA;
A_bedD=2.46;
U_bedD=1724.14;
UA_bedD=A_bedD*U_bedD;
%%%%%%%%%Physical and thermal Properties of Condenser%%%%%%%%%%%%%%%%
A_cond=3.73; %Heat transfer area of condenser, [m^2]
Mhex_cond=24.28; %Mass of heat exchanger (Condenser), [kg]
U_cond=4115.23; %Heat transfer coefficient of condenser, [W/m^2.K]
Mw_cond=20; %Mass of water in Condenser [kg]
Cphex_cond=386; %Specific heat of condenser heat exchanger tubes [J/kgK]
%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%%
A_evap=1.91; %Surface area of evaporator, (m^2)
Mhex_evap=12.45; %Mass of heat exchanger (Evaporator), (kg)
U_evap=2557.54; %Heat transfer coefficient of evaporator, (W/m^2.K)
Cphex_evap=386; %Specific heat of heat exchanger tube of condenser and evaporator [J/kgK]
Mw_evap_max=50; %Initial amount of water in evaporator [kg]
%%%%%%%%%%%Initial values%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Time_i=0;
t=Time_i;
T_bedA_i=27;
T_bedD_i=27;
T_evap_i=27;
T_cond_i=27;
w_bedA_i=0.05;
w_bedD_i=0.05;
Mw_evap_i=Mw_evap_max;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=0;
for M=1:N_cycle %Loop for cycle iteration
t_round=0; %flag determinator
Q_D=0;
Q_E=0;
for L=1:step:t_cycle %Loop for inside cycle iteration
S=S+1;
t=t+step;
Time(S)=t;
t_round=t_round+step;
if (t_round<=t_normal)
FLAG1=1; %Normal cycle operatoon
else
FLAG1=0; %Switching operation (preheating/precooling)
end
Timerange=[Time_i t];
Initial=[w_bedA_i T_bedA_i w_bedD_i T_bedD_i T_cond_i T_evap_i Mw_evap_i];
options=odeset('RelTol',1e-4,'AbsTol',1e-4);
Y=ode45(@dYdt,Timerange,Initial,options); %solving all the differential equation
Y_t=deval(Y,t); %Y_t = deval(Y,t) evaluates the solutions of the differential equations at the points contained in t.
w_bedA(S)=Y_t(1);
T_bedA(S)=Y_t(2);
w_bedD(S)=Y_t(3);
T_bedD(S)=Y_t(4);
T_cond(S)=Y_t(5);
T_evap(S)=Y_t(6);
Mw_evap(S)=Y_t(7);
if(FLAG1==1)
Thw_bed_in=Thw_in;
Thw_bed_out=T_w_bedD_out;
Tcw_bed_in=Tcw_in;
Tcw_bed_out=T_w_bedA_out;
elseif(FLAG1==0)
Thw_bed_in=Thw_in;
Thw_bed_out=T_w_bedA_out;
Tcw_bed_in=Tcw_in;
Tcw_bed_out=T_w_bedD_out;
end
T_cw_cond_in=Tcw_in;
T_cw_cond_out=T_w_cond_out;
%%%%%%%%%%%Adsorption/Desorption switching%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if((mod(M,2)==1)||(M==1))
T_bed2(S)=T_bedD(S);
W_bed2(S)=w_bedD(S);
T_bed1(S)=T_bedA(S);
W_bed1(S)=w_bedA(S);
else
T_bed1(S)=T_bedD(S);
W_bed1(S)=w_bedD(S);
T_bed2(S)=T_bedA(S);
W_bed2(S)=w_bedA(S);
end
Tcond_in(S)=T_cw_cond_in;
Tcond_out(S)=T_cw_cond_out;
T_cw_in(S)=Tcw_bed_in;
T_cw_out(S)=Tcw_bed_out;
T_hw_in(S)=Thw_bed_in;
T_hw_out(S)=Thw_bed_out;
T_chill_in(S)=T_chw_in;
T_chill_out(S)=T_chw_out;
dq_D=(m_dot_hw_bed*step*Cpw*(Thw_bed_in-Thw_bed_out))/t_cycle;
Q_D=Q_D+dq_D;
dq_E=(m_dot_chw*Cpw*step*(T_chw_in-T_chw_out))/t_cycle;
Q_E=Q_E+dq_E;
%Storing the newly calculated values as the new initial
%values for the next round of caluculation (ncycle value)
w_bedA_i=Y_t(1);
T_bedA_i=Y_t(2);
w_bedD_i=Y_t(3);
T_bedD_i=Y_t(4);
T_cond_i=Y_t(5);
T_evap_i=Y_t(6);
Mw_evap_i=Y_t(7);
end
Q_evap(M)=Q_E;
Q_des(M)=Q_D;
COP(M)=Q_evap(M)/Q_des(M);
SCP(M)=Q_evap(M)/(2*Mads);
cycle(M)=M;
w_bedA_i=Y_t(3);
T_bedA_i=Y_t(4);
w_bedD_i=Y_t(1);
T_bedD_i=Y_t(2);
end
figure(1)
plot(Time,T_bed1,Time,T_bed2,Time,T_cond,Time,T_evap);
hold all;
figure(2)
plot(Time,W_bed1,Time,W_bed2);
hold all;
figure(3)
plot(Time,T_cw_out,Time,T_hw_out,Time,T_chill_out,Time,Tcond_out);
hold all;
figure(4)
plot(cycle,SCP);
hold all;
figure(5)
plot(cycle,COP);
hold all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Solver %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=dYdt(t,Y)
global Tcw_in Thw_in T_chw_in T_chw_out T_w_bedA_out T_w_bedD_out T_w_cond_out P_bedA P_bedD w_bedA w_bedD
global Ea Q R Rp Dso A Ao A1 A2 A3 B Bo B1 B2 B3
global m_dot_cw_bed m_dot_cw_cond m_dot_hw_bed m_dot_chw Cpw Cpv Lv
global MhexCphex_bed Mhex_cond Cphex_cond Mhex_evap Cphex_evap Mads Mw_cond Mw_evap_max Cpads A_bedA U_bedA A_bedD U_bedD A_cond U_cond A_evap U_evap UA_bedA UA_bedD
global t_normal t_sw t_cycle N_cycle step
global FLAG1
dy=zeros(7,1); %dy is a 7*1 matrix
%Storing initial values from main function to local variables
w_bedA=Y(1);
T_bedA=Y(2);
w_bedD=Y(3);
T_bedD=Y(4);
T_cond=Y(5);
T_evap=Y(6);
Mw_evap=Y(7);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bed 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(FLAG1==1)
T_w_bedA_in=Tcw_in;
m_dot_w_bedA=m_dot_cw_bed;
elseif(FLAG1==0)
T_w_bedA_in=Thw_in;
m_dot_w_bedA=m_dot_hw_bed;
end
T_w_bedA_out=BedHex(T_bedA,T_w_bedA_in,m_dot_w_bedA);
if(FLAG1==1)
P_bedA=10^(8.07131-(1730.63/(T_evap+233.426)));
Ps_bedA=10^(8.07131-(1730.63/(T_bedA+233.426)));
w_star_bedA=(Ao+(A1*(T_bedA+273.16))+(A2*((T_bedA+273.16)^2))+(A3*((T_bedA+273.16)^3)))*(P_bedA/Ps_bedA)^(Bo+(B1*(T_bedA+273.16))+(B2*((T_bedA+273.16)^2))+(B3*((T_bedA+273.16)^3)));
w_bed_const1=(15*Dso)/(Rp^2);
elseif(FLAG1==0)
P_bedA=10^(8.07131-(1730.63/(T_bedA+233.426)));
Ps_bedA=10^(8.07131-(1730.63/(T_bedA+233.426)));
w_star_bedA=(Ao+(A1*(T_bedA+273.16))+(A2*((T_bedA+273.16)^2))+(A3*((T_bedA+273.16)^3)))*(P_bedA/Ps_bedA)^(Bo+(B1*(T_bedA+273.16))+(B2*((T_bedA+273.16)^2))+(B3*((T_bedA+273.16)^3)));
w_bed_const1=0;
end
T_bed_const1=FLAG1*Mads*Q*(15*Dso/(Rp^2));
T_bed_const2=FLAG1*Mads*Cpv*(T_evap-T_bedA)*(15*Dso/(Rp^2));
T_bed_const3=m_dot_w_bedA*Cpw*(T_w_bedA_in-T_w_bedA_out);
T_bed_const4=(Mads*Cpads)+(MhexCphex_bed);
T_bed_const5=Mads*Cpw;
if(FLAG1==1)
T_w_bedD_in=Thw_in;
m_dot_w_bedD=m_dot_hw_bed;
elseif(FLAG1==0)
T_w_bedD_in=Tcw_in;
m_dot_w_bedD=m_dot_cw_bed;
end
T_w_bedD_out=BedHex(T_bedD,T_w_bedD_in,m_dot_w_bedD);
if(FLAG1==1)
P_bedD=10^(8.07131-(1730.63/(T_cond+233.426)));
Ps_bedD=10^(8.07131-(1730.63/(T_bedD+233.426)));
w_star_bedD=(Ao+(A1*(T_bedD+273.16))+(A2*((T_bedD+273.16)^2))+(A3*((T_bedD+273.16)^3)))*(P_bedD/Ps_bedD)^(Bo+(B1*(T_bedD+273.16))+(B2*((T_bedD+273.16)^2))+(B3*((T_bedD+273.16)^3)));
w_bed_const2=(15*Dso)/(Rp^2);
elseif(FLAG1==0)
P_bedD=10^(8.07131-(1730.63/(T_bedD+233.426)));
Ps_bedD=10^(8.07131-(1730.63/(T_bedD+233.426)));
w_star_bedD=(Ao+(A1*(T_bedD+273.16))+(A2*((T_bedD+273.16)^2))+(A3*((T_bedD+273.16)^3)))*(P_bedD/Ps_bedD)^(Bo+(B1*(T_bedD+273.16))+(B2*((T_bedD+273.16)^2))+(B3*((T_bedD+273.16)^3)));
w_bed_const2=0;
end
T_bed_const6=FLAG1*Mads*Q*(15*Dso/(Rp^2));
T_bed_const7=m_dot_w_bedD*Cpw*(T_w_bedD_in-T_w_bedD_out);
T_bed_const8=(Mads*Cpads)+(MhexCphex_bed);
T_bed_const9=Mads*Cpw;
T_w_cond_in=Tcw_in;
m_dot_w_cond=m_dot_cw_cond;
T_w_cond_out=CondHex(T_cond,T_w_cond_in,m_dot_w_cond);
T_cond_const1=-FLAG1*Mads*(15*Dso/(Rp^2));
T_cond_const2=-FLAG1*Mads*Cpv*(T_bedD-T_cond)*(15*Dso/(Rp^2));
T_cond_const3=m_dot_w_cond*Cpw*(T_w_cond_in-T_w_cond_out);
T_cond_const4=Mw_cond*Cpw+Mhex_cond*Cphex_cond;
T_chw_out=EvapHex(T_evap,T_chw_in,m_dot_chw);
T_evap_const1=-FLAG1*Mads*Lv*(15*Dso/(Rp^2));
T_evap_const2=-FLAG1*Mads*Cpw*(T_cond-T_evap)*(15*Dso/(Rp^2));
T_evap_const3=m_dot_chw*Cpw*(T_chw_in-T_chw_out);
T_evap_const4=Mw_evap*Cpw+Mhex_evap*Cphex_evap;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dy(1)=w_bed_const1*(-Ea/(R*(Y(2)+273.16)))*(w_star_bedA-Y(1));
dy(2)=((T_bed_const1+T_bed_const2)*exp(-Ea/(R*(Y(2)+273.16)))*(w_star_bedA-Y(1))+T_bed_const3)/(T_bed_const4+(T_bed_const5*Y(1)));
dy(3)=w_bed_const2*(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3));
dy(4)=(T_bed_const6*exp(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3))+T_bed_const7)/(T_bed_const8+(T_bed_const9*Y(3)));
dy(5)=((T_cond_const1+T_cond_const2)*exp(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3))+T_cond_const3)/T_cond_const4;
dy(6)=((T_evap_const1+T_evap_const2)*exp(-Ea/(R*(Y(2)+273.16)))*(w_star_bedA-Y(1))+T_evap_const3)/T_evap_const4;
dy(7)=-Mads*FLAG1*((w_bed_const1*(w_star_bedA-Y(1))*exp(-Ea/(R*(Y(2)+273.16))))+(w_bed_const2*exp(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3))));
end
function [T_w_bed_out]=BedHex(T_bed,T_w_bed_in,m_dot_w_bed)
global UA_bedA UA_bedD Thw_in Tcw_in Cpw
if(T_w_bed_in==Tcw_in)
T_w_bed_out=T_bed+(T_w_bed_in-T_bed)*exp((-UA_bedA)/(m_dot_w_bed*Cpw));
elseif(T_w_bed_in==Thw_in)
T_w_bed_out=T_bed+(T_w_bed_in-T_bed)*exp((-UA_bedD)/(m_dot_w_bed*Cpw));
end
end
function T_w_cond_out=CondHex(T_cond,T_w_cond_in,m_dot_w_cond)
global U_cond A_cond Cpw
T_w_cond_out=T_cond+(T_w_cond_in-T_cond)*exp((-U_cond*A_cond)/(m_dot_w_cond*Cpw));
end
function T_chw_out=EvapHex(T_evap,T_chw_in,m_dot_chw)
global U_evap A_evap Cpw
T_chw_out=T_evap+(T_chw_in-T_evap)*exp((-U_evap*A_evap)/(m_dot_chw*Cpw));
end
  1 comentario
Stephen23
Stephen23 el 1 de Feb. de 2022
Get rid of those anti-pattern GLOBAL variables by parameterizing the functions:
https://www.mathworks.com/help/matlab/math/parameterizing-functions.html
Either:
  • use nested functions (probably best for your case)
  • put the parameters into one/several structures and pass those via anonymous functions.

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Cris LaPierre
Cris LaPierre el 1 de Feb. de 2022
Editada: Cris LaPierre el 1 de Feb. de 2022
Ran in:
All your y values are NaNs, which is why nothing is showing up. Check your equations.
I'd also suggest not using globals, as that can have some unintended consequences. I'd try something like nested functions. I did look into the why a little bit, but it was going to take more time that I had.
main
function main
%%%%%%%%%%%%%%%%%%%%%%%%%%% Cycle Duration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_normal=500; %Time for normal operation [s]
t_sw=30; %Time for switching operation [s]
t_cycle=t_normal+t_sw; %Total cycle time [s]
step=5; %Evaluate data every 5 sec
N_cycle=20; %No of cycles to be evaluated
%%%%%%%%%%%%%%%%%%%%%%%%%%% Operating conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tcw_in=30; %Temperature of cold water in [C]
Thw_in=85; %Temperature of hot water in [C]
T_chw_in=14; %Temperature of chiled water in [C]
m_dot_cw_bed=1.7; %Temperature of bed cooling water [C]
m_dot_cw_cond=1.3; %Temperature of condenser cooling water [C]
m_dot_hw_bed=1.7; %Temperature of bed heating water [C]
m_dot_chw=0.7; %Temperature of chilled water flow in evaporator [C]
%%%%%%%%%%%%%% Physical and thermal properties of adsorbate-water %%%%%%%%%%%%%
Cpw=4180; %Specific heat of water [J/kgK]
Cpv=1890; %Specific heat of water vapor [J/kgK]
Lv=2500*10^3; %Latent heat of vaporization
%%%%%%%%%%%%%% Physical and thermal properties of adsorbent %%%%%%%%%%%%%%%%
Cpads=924; %Specific heat capacity of silica gel, [J/kg.K]
Mads=47; %Mass of silica, [kg]
Rp=1.6*10^-4; %Radious of silica gel particle, [m]
Dso=2.54*10^-4; %Pre-exponential constant, (m^2 S^-1)
Ea=4.2*10^4; %Activation energy of surface diffusion, (J/mol)
R=8.314; %Ideal gas constant, (j/mol.K)
Q=2.8*10^6; %Heat of adsorption, (J/kg)
%%%%%%%%%%%%%%%%%%%%%%% SBK mode constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ao=-6.5314;
A1=0.072452;
A2=-0.23951*10^-3;
A3=0.25493*10^-6;
Bo=-15.587;
B1=0.15915;
B2=-0.50612*10^-3;
B3=0.5329*10^-6;
%%%%%%%%%%Physical and thermal properties of adsorber/desorber%%%%%%%
MhexCphex_bed=77719.4; %Mass capacitance of heat exchanger (Adsorber/desorber), [J/K]
A_bedA=2.46; %Surface area of adsorber, [m^2]
U_bedA=1602.56; %Heat transfer coefficient of desorber, [W/m^2.K]
UA_bedA=U_bedA*A_bedA;
A_bedD=2.46;
U_bedD=1724.14;
UA_bedD=A_bedD*U_bedD;
%%%%%%%%%Physical and thermal Properties of Condenser%%%%%%%%%%%%%%%%
A_cond=3.73; %Heat transfer area of condenser, [m^2]
Mhex_cond=24.28; %Mass of heat exchanger (Condenser), [kg]
U_cond=4115.23; %Heat transfer coefficient of condenser, [W/m^2.K]
Mw_cond=20; %Mass of water in Condenser [kg]
Cphex_cond=386; %Specific heat of condenser heat exchanger tubes [J/kgK]
%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%%
A_evap=1.91; %Surface area of evaporator, (m^2)
Mhex_evap=12.45; %Mass of heat exchanger (Evaporator), (kg)
U_evap=2557.54; %Heat transfer coefficient of evaporator, (W/m^2.K)
Cphex_evap=386; %Specific heat of heat exchanger tube of condenser and evaporator [J/kgK]
Mw_evap_max=50; %Initial amount of water in evaporator [kg]
%%%%%%%%%%%Initial values%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Time_i=0;
t=Time_i;
T_bedA_i=27;
T_bedD_i=27;
T_evap_i=27;
T_cond_i=27;
w_bedA_i=0.05;
w_bedD_i=0.05;
Mw_evap_i=Mw_evap_max;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=0;
for M=1:N_cycle %Loop for cycle iteration
t_round=0; %flag determinator
Q_D=0;
Q_E=0;
for L=1:step:t_cycle %Loop for inside cycle iteration
S=S+1;
t=t+step;
Time(S)=t;
t_round=t_round+step;
if (t_round<=t_normal)
FLAG1=1; %Normal cycle operatoon
else
FLAG1=0; %Switching operation (preheating/precooling)
end
Timerange=[Time_i t];
Initial=[w_bedA_i T_bedA_i w_bedD_i T_bedD_i T_cond_i T_evap_i Mw_evap_i];
options=odeset('RelTol',1e-4,'AbsTol',1e-4);
Y=ode45(@dYdt,Timerange,Initial,options); %solving all the differential equation
Y_t=deval(Y,t); %Y_t = deval(Y,t) evaluates the solutions of the differential equations at the points contained in t.
w_bedA(S)=Y_t(1);
T_bedA(S)=Y_t(2);
w_bedD(S)=Y_t(3);
T_bedD(S)=Y_t(4);
T_cond(S)=Y_t(5);
T_evap(S)=Y_t(6);
Mw_evap(S)=Y_t(7);
%%% I had to modify this next part to recreate 4 variables
if(FLAG1==1)
T_w_bedA_out=BedHex(T_bedA(S),Tcw_in,m_dot_cw_bed);
T_w_bedD_out=BedHex(T_bedD(S),Thw_in,m_dot_hw_bed);
Thw_bed_in=Thw_in;
Thw_bed_out=T_w_bedD_out;
Tcw_bed_in=Tcw_in;
Tcw_bed_out=T_w_bedA_out;
elseif(FLAG1==0)
T_w_bedA_out=BedHex(T_bedA(S),Thw_in,m_dot_hw_bed);
T_w_bedD_out=BedHex(T_bedD(S),Tcw_in,m_dot_cw_bed);
Thw_bed_in=Thw_in;
Thw_bed_out=T_w_bedA_out;
Tcw_bed_in=Tcw_in;
Tcw_bed_out=T_w_bedD_out;
end
T_w_cond_out=CondHex(T_cond(S),Tcw_in,m_dot_cw_cond);
T_chw_out=EvapHex(T_evap(S),T_chw_in,m_dot_chw);
T_cw_cond_in=Tcw_in;
T_cw_cond_out=T_w_cond_out;
%%%%%%%%%%%Adsorption/Desorption switching%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if((mod(M,2)==1)||(M==1))
T_bed2(S)=T_bedD(S);
W_bed2(S)=w_bedD(S);
T_bed1(S)=T_bedA(S);
W_bed1(S)=w_bedA(S);
else
T_bed1(S)=T_bedD(S);
W_bed1(S)=w_bedD(S);
T_bed2(S)=T_bedA(S);
W_bed2(S)=w_bedA(S);
end
Tcond_in(S)=T_cw_cond_in;
Tcond_out(S)=T_cw_cond_out;
T_cw_in(S)=Tcw_bed_in;
T_cw_out(S)=Tcw_bed_out;
T_hw_in(S)=Thw_bed_in;
T_hw_out(S)=Thw_bed_out;
T_chill_in(S)=T_chw_in;
T_chill_out(S)=T_chw_out;
dq_D=(m_dot_hw_bed*step*Cpw*(Thw_bed_in-Thw_bed_out))/t_cycle;
Q_D=Q_D+dq_D;
dq_E=(m_dot_chw*Cpw*step*(T_chw_in-T_chw_out))/t_cycle;
Q_E=Q_E+dq_E;
%Storing the newly calculated values as the new initial
%values for the next round of caluculation (ncycle value)
w_bedA_i=Y_t(1);
T_bedA_i=Y_t(2);
w_bedD_i=Y_t(3);
T_bedD_i=Y_t(4);
T_cond_i=Y_t(5);
T_evap_i=Y_t(6);
Mw_evap_i=Y_t(7);
end
Q_evap(M)=Q_E;
Q_des(M)=Q_D;
COP(M)=Q_evap(M)/Q_des(M);
SCP(M)=Q_evap(M)/(2*Mads);
cycle(M)=M;
w_bedA_i=Y_t(3);
T_bedA_i=Y_t(4);
w_bedD_i=Y_t(1);
T_bedD_i=Y_t(2);
end
figure(1)
plot(Time,T_bed1,Time,T_bed2,Time,T_cond,Time,T_evap);
figure(2)
plot(Time,W_bed1,Time,W_bed2);
figure(3)
plot(Time,T_cw_out,Time,T_hw_out,Time,T_chill_out,Time,Tcond_out);
figure(4)
plot(cycle,SCP);
figure(5)
plot(cycle,COP);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Solver %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=dYdt(t,Y)
dy=zeros(7,1); %dy is a 7*1 matrix
%Storing initial values from main function to local variables
w_bedA=Y(1);
T_bedA=Y(2);
w_bedD=Y(3);
T_bedD=Y(4);
T_cond=Y(5);
T_evap=Y(6);
Mw_evap=Y(7);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bed 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(FLAG1==1)
T_w_bedA_in=Tcw_in;
m_dot_w_bedA=m_dot_cw_bed;
elseif(FLAG1==0)
T_w_bedA_in=Thw_in;
m_dot_w_bedA=m_dot_hw_bed;
end
T_w_bedA_out=BedHex(T_bedA,T_w_bedA_in,m_dot_w_bedA);
if(FLAG1==1)
P_bedA=10^(8.07131-(1730.63/(T_evap+233.426)));
Ps_bedA=10^(8.07131-(1730.63/(T_bedA+233.426)));
w_star_bedA=(Ao+(A1*(T_bedA+273.16))+(A2*((T_bedA+273.16)^2))+(A3*((T_bedA+273.16)^3)))*(P_bedA/Ps_bedA)^(Bo+(B1*(T_bedA+273.16))+(B2*((T_bedA+273.16)^2))+(B3*((T_bedA+273.16)^3)));
w_bed_const1=(15*Dso)/(Rp^2);
elseif(FLAG1==0)
P_bedA=10^(8.07131-(1730.63/(T_bedA+233.426)));
Ps_bedA=10^(8.07131-(1730.63/(T_bedA+233.426)));
w_star_bedA=(Ao+(A1*(T_bedA+273.16))+(A2*((T_bedA+273.16)^2))+(A3*((T_bedA+273.16)^3)))*(P_bedA/Ps_bedA)^(Bo+(B1*(T_bedA+273.16))+(B2*((T_bedA+273.16)^2))+(B3*((T_bedA+273.16)^3)));
w_bed_const1=0;
end
T_bed_const1=FLAG1*Mads*Q*(15*Dso/(Rp^2));
T_bed_const2=FLAG1*Mads*Cpv*(T_evap-T_bedA)*(15*Dso/(Rp^2));
T_bed_const3=m_dot_w_bedA*Cpw*(T_w_bedA_in-T_w_bedA_out);
T_bed_const4=(Mads*Cpads)+(MhexCphex_bed);
T_bed_const5=Mads*Cpw;
if(FLAG1==1)
T_w_bedD_in=Thw_in;
m_dot_w_bedD=m_dot_hw_bed;
elseif(FLAG1==0)
T_w_bedD_in=Tcw_in;
m_dot_w_bedD=m_dot_cw_bed;
end
T_w_bedD_out=BedHex(T_bedD,T_w_bedD_in,m_dot_w_bedD);
if(FLAG1==1)
P_bedD=10^(8.07131-(1730.63/(T_cond+233.426)));
Ps_bedD=10^(8.07131-(1730.63/(T_bedD+233.426)));
w_star_bedD=(Ao+(A1*(T_bedD+273.16))+(A2*((T_bedD+273.16)^2))+(A3*((T_bedD+273.16)^3)))*(P_bedD/Ps_bedD)^(Bo+(B1*(T_bedD+273.16))+(B2*((T_bedD+273.16)^2))+(B3*((T_bedD+273.16)^3)));
w_bed_const2=(15*Dso)/(Rp^2);
elseif(FLAG1==0)
P_bedD=10^(8.07131-(1730.63/(T_bedD+233.426)));
Ps_bedD=10^(8.07131-(1730.63/(T_bedD+233.426)));
w_star_bedD=(Ao+(A1*(T_bedD+273.16))+(A2*((T_bedD+273.16)^2))+(A3*((T_bedD+273.16)^3)))*(P_bedD/Ps_bedD)^(Bo+(B1*(T_bedD+273.16))+(B2*((T_bedD+273.16)^2))+(B3*((T_bedD+273.16)^3)));
w_bed_const2=0;
end
T_bed_const6=FLAG1*Mads*Q*(15*Dso/(Rp^2));
T_bed_const7=m_dot_w_bedD*Cpw*(T_w_bedD_in-T_w_bedD_out);
T_bed_const8=(Mads*Cpads)+(MhexCphex_bed);
T_bed_const9=Mads*Cpw;
T_w_cond_in=Tcw_in;
m_dot_w_cond=m_dot_cw_cond;
T_w_cond_out=CondHex(T_cond,T_w_cond_in,m_dot_w_cond);
T_cond_const1=-FLAG1*Mads*(15*Dso/(Rp^2));
T_cond_const2=-FLAG1*Mads*Cpv*(T_bedD-T_cond)*(15*Dso/(Rp^2));
T_cond_const3=m_dot_w_cond*Cpw*(T_w_cond_in-T_w_cond_out);
T_cond_const4=Mw_cond*Cpw+Mhex_cond*Cphex_cond;
T_chw_out=EvapHex(T_evap,T_chw_in,m_dot_chw);
T_evap_const1=-FLAG1*Mads*Lv*(15*Dso/(Rp^2));
T_evap_const2=-FLAG1*Mads*Cpw*(T_cond-T_evap)*(15*Dso/(Rp^2));
T_evap_const3=m_dot_chw*Cpw*(T_chw_in-T_chw_out);
T_evap_const4=Mw_evap*Cpw+Mhex_evap*Cphex_evap;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dy(1)=w_bed_const1*(-Ea/(R*(Y(2)+273.16)))*(w_star_bedA-Y(1));
dy(2)=((T_bed_const1+T_bed_const2)*exp(-Ea/(R*(Y(2)+273.16)))*(w_star_bedA-Y(1))+T_bed_const3)/(T_bed_const4+(T_bed_const5*Y(1)));
dy(3)=w_bed_const2*(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3));
dy(4)=(T_bed_const6*exp(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3))+T_bed_const7)/(T_bed_const8+(T_bed_const9*Y(3)));
dy(5)=((T_cond_const1+T_cond_const2)*exp(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3))+T_cond_const3)/T_cond_const4;
dy(6)=((T_evap_const1+T_evap_const2)*exp(-Ea/(R*(Y(2)+273.16)))*(w_star_bedA-Y(1))+T_evap_const3)/T_evap_const4;
dy(7)=-Mads*FLAG1*((w_bed_const1*(w_star_bedA-Y(1))*exp(-Ea/(R*(Y(2)+273.16))))+(w_bed_const2*exp(-Ea/(R*(Y(4)+273.16)))*(w_star_bedD-Y(3))));
end
function [T_w_bed_out]=BedHex(T_bed,T_w_bed_in,m_dot_w_bed)
if(T_w_bed_in==Tcw_in)
T_w_bed_out=T_bed+(T_w_bed_in-T_bed)*exp((-UA_bedA)/(m_dot_w_bed*Cpw));
elseif(T_w_bed_in==Thw_in)
T_w_bed_out=T_bed+(T_w_bed_in-T_bed)*exp((-UA_bedD)/(m_dot_w_bed*Cpw));
end
end
function T_w_cond_out=CondHex(T_cond,T_w_cond_in,m_dot_w_cond)
T_w_cond_out=T_cond+(T_w_cond_in-T_cond)*exp((-U_cond*A_cond)/(m_dot_w_cond*Cpw));
end
function T_chw_out=EvapHex(T_evap,T_chw_in,m_dot_chw)
T_chw_out=T_evap+(T_chw_in-T_evap)*exp((-U_evap*A_evap)/(m_dot_chw*Cpw));
end
end
  1 comentario
Star Strider
Star Strider el 1 de Feb. de 2022
I would suggest putting all the globals into a structure, and then passing the structure as an additional parameter.

Iniciar sesión para comentar.

Categorías

Más información sobre Heat and Mass Transfer en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by