How can get the value for each time step for a variable inside of a function?

I need to get the value of some variables (Defined as golobal variables in the attached file) for each time steps. These variables are calculated inside of the fucntion for each time steps. Based on my code I am geeting the value only for the last time step. How can i get the value for each time steps

I have attached the files (main file. function file)

Answer 1

[Stephan]

Hi,

you try this code (all in one file as nested function). The values you want to have in every iteration are now added to the solution and also saved in the results table. This should be what you want.

Result = main
function Result = main
%Initial conditions
Wads=0.045;
Tads=80;
Wdes=0.035;
Tdes=27;
Teva=7;
Tcond=27;
Mweva=0.65;         
Initial=[Wads Tads Wdes Tdes Mweva Teva Tcond 0 0 0 0];           
tspan=[1:1:800];
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,Soln]=ode45(@(t,Soln)ODEdw(t,Soln),tspan,Initial,opts);

colNames = {'Time','Wads','Tads','Wdes','Tdes','Mweva', 'Teva','Tcond',...
    'Tcwaout', 'Thwdout', 'Tchwout', 'Tcwcout'};
Result = array2table([tspan',Soln],'VariableNames',colNames);

function [dYdt] =ODEdw(~,Soln)
Ea=4.2*10^4;         %Activation energy of surface diffusion, (J/mol)
R=8.314;             %Ideal gas constant, (j/mol.K)
Dso=2.54*10^-4;       %Pre-exponential constant, (m^2 S^-1)
Qs=2.8*10^6;         %Heat of adsorption, (J/kg)
%%%%%%%%%%%Operating conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Thwdin=85;            %Inlet temperature of hot water in desorber, (oC)
Tcwain=27;            %Inlet temperature of cooling water in adsorber, (oC)
Tcwcin=27;            %Inlet temperature of cooling water in condenser, (oC)
Tchwin=14;            %Inlet temperature of chilling water in evaporator, (oC)
mcwa=8/60;            %Mass flow rate of cooling water in adsorber, (kg/s)
mcwc=20/60;           %Mass flow rate of cooling water in condenser, (kg/s)
mchw=1.6/60;          %Mass flow rate of chilling water in evaporator, (kg/s)
mhwd=7/60;            %Mass flow rate of hot water in desorber, (kg/s)
%%%%%%%%%Physical and thermal properties of adsorbent%%%%%%%%%%
Cps=924;              %Specific heat capacity of silica gel, (J/kg.K)
Ms=9;                 %Mass of silica, (kg)
Rp=5.6*10^-4;         %Radious of silica gel particle, (m)
ps=700;              %density of silica gel, (kg/m^3)
Por=0.3;              %Porosity of silica gel,

%%%%%%Physical and thermal properties of adsorbate%%%%%%%%%%%%%
Cpw=4200;             %Specific heat capcity of water, (J/kg.k)
Cpv=1900;             %Specific heat capacity of water vapor (J/kg.K)
M=18.02;              %Molar mass of water, (g/mol)
Lv=2500*10^3;         %Lantent heat of vaporization (J/kg) 
%%%%%%%%%Physical and thermal Properties of Condenser%%%%%%%%%%%%%
Acond=5.76;            %Surface area of condenser, (m^2)
Mhexcond=19.87;        %Mass of heat exchanger (Condenser), (kg)
Ucond=486.87;          %Heat transfer coefficient of condenser, (W/m^2.K)
UAcond=Acond*Ucond;
Mwcond=1.0;         %Mass of water in Condenser (kg)
%%%%%%%%%%Physical and thermal properties of adsorber/desorber   %%%%%%%
Mhex=65;             %Mass of heat exchanger (Adsorber/desorber), (kg)
Cphex=386;           %Heat capacity of heat excahgner, (J/kg.K)
Aad=2.11;            %Surface area of adsorber/desorber, (m^2)
Uads=189.60;         %Heat transfer coefficient of adsorber, (W/m^2.K)
Udes=203.99;         %Heat transfer coefficient of desorber, (W/m^2.K)
UAads=Aad*Uads;
UAdes=Aad*Udes;
Vt=0.5;              %volume of adsorber/desorber tank, (m^3)
%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%
Aeva=0.34;           %Surface area of evaporator, (m^2)
Mhexeva=16.1;        %Mass of heat exchanger (Evaporator), (kg)
Ueva=302.58;         %Heat transfer coefficient of evaporator, (W/m^2.K)
UAeva=Aeva*Ueva;

Wads=Soln(1);    
Tads=Soln(2);  
Wdes=Soln(3);    
Tdes=Soln(4);
Mweva=Soln(5);
Teva=Soln(6);
Tcond=Soln(7);
Kads=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tads+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwa=133.32*exp(18.3-(3820/((Teva+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)
Psga=133.32*exp(18.3-(3820/((Tads+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqa=0.346*(Pwa/Psga)^0.625;
Aads=(Kads*(Weqa-Wads));
TLads=(Ms*(Cpw+Cpw*Wads)+Mhex*Cphex); %Left hand side of adsorber equation
Tcwaout=Tads+(Tcwain-Tads)*exp((-UAads)/(mcwa*Cpw));
   

Kdes=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tdes+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwd=133.32*exp(18.3-(3820/((Tcond+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)
Psgd=133.32*exp(18.3-(3820/((Tdes+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqd=0.346*(Pwd/Psgd)^0.625;
Ades=(Kdes*(Weqd-Wdes));
TLdes=(Ms*(Cps+Cpw*Wdes)+Mhex*Cphex); %Left hand side of desorber equation
Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(mhwd*Cpw));
  

TLeva=(Mweva*Cpw+Mhexeva*Cphex); %Left hand side of evaporator equation
Tchwout=Teva+(Tchwin-Teva)*exp((-UAeva)/(mchw*Cpw));
  


TLcond=(Mwcond*Cpw+Mhexcond*Cphex); %Left hand side of evaporator equation
Tcwcout=Tcond+(Tcwcin-Tcond)*exp((-UAcond)/(mcwc*Cpw));

dYdt=[Kads*(Weqa-Wads);((Qs*Ms*Aads)+(Ms*Cpv*Aads*(Teva-Tads))+mcwa*Cpw*(Tcwain-Tcwaout))/TLads;
    Kdes*(Weqd-Wdes); ((Qs*Ms*Ades)+(mcwa*Cpw*(Thwdin-Thwdout)))/TLdes; 
    -Ms*((Kads*(Weqa-Wads))+(Kdes*(Weqd-Wdes)));((-Lv*Ms*Aads)+Ades*Ms*Cpw*(Tcond-Teva)+mchw*Cpw*(Tchwin-Tchwout))/TLeva;
    (Ades*Ms*Cpv*(Tdes-Tcond)+(Lv*Ms*Ades)+mcwc*Cpw*(Tcwcin-Tcwcout))/TLcond; Tcwaout; Thwdout; Tchwout; Tcwcout];
end
end

Best regards

Stephan

[Md Shamim Ahamed]

Thanks for your answer. The code is working, but there is problem. The vaiables (Tcwaout Thwdout Tchwout Tcwcout) should not assign initial value becuase these value would be calculated based on the initial value (t=0) and for next time steps these variaboles would be calcualted based on the final solution ODE paramters (Tads, Tdes, Teva, Tcond) from ODE45, herefore the results are not correct. Do you have any suggestions.

[Md Shamim Ahamed]

Hi

Based on code, its displaying the cumulated (summation of last interation and cureent interation value) value for these variable (Tcwaout Thwdout Tchwout Tcwcout) . Any help is appreciated.

Kind regrads

[Stephan]

In every step of the integration the values of Tcwaout, Thwdout, Tchwout, Tcwcout are calculated by using the actual values of Tads, Tdes, Teva, Tcond. This seems to be correct for me.

The only issue seems to be the initial value - is it possible to calculate the initial value for Tcwaout, Thwdout, Tchwout, Tcwcout by using the initial values of the other Tads, Tdes, Teva, Tcond? Would not this solve the problem?

If not, what is the expected behavior? How should it be? How would you do it with pencil and paper?

[Md Shamim Ahamed]

Hi Stephan

You are right that its doing integration the values of Tcwaout, Thwdout, Tchwout, Tcwcout, but this is now what I want. I want these variable will be calculated for each time steps for the final value of Tads, Tdes, Teva, Tcond.

In your code the value of Tcwaout, Thwdout, Tchwout, Tcwcout are increaing until the sumilation stop. However, it could increase or decrease denepding on the value of Tads, Tdes, Teva, Tcond

Like I have Tads=20, Tdes=30, Teva=14. Tcond=27

Now these variables (Tcwaout, Thwdout, Tchwout, Tcwcout) would be calcualted based on the following equations:

Tcwaout=Tads+(Tcwain-Tads)*exp((-UAads)/(mcwa*Cpw));

Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(mhwd*Cpw));

Tchwout=Teva+(Tchwin-Teva)*exp((-UAeva)/(mchw*Cpw));

Tcwcout=Tcond+(Tcwcin-Tcond)*exp((-UAcond)/(mcwc*Cpw));

[Stephan]

I see - Torstens solution is doing exactly what you want. I allowed myself to solve the last problem, that Result is not given to workspace. I posted this in Torstens answer, because it were his thoughts that solved the problem finally. I think now you have a working solution doing what you want it to do.

Answer 2

[Torsten]

Result = main
function Result = main
%Initial conditions
Wads=0.045;
Tads=80;
Wdes=0.035;
Tdes=27;
Teva=7;
Tcond=27;
Mweva=0.65;         
Initial=[Wads Tads Wdes Tdes Mweva Teva Tcond];           
tspan=[1:1:800];
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
iflag = 1;
[t,Soln]=ode45(@(t,Soln)ODEdw(t,Soln,iflag),tspan,Initial,opts);

iflag = 2;
Tcwaout = zeros(numel(t),1); 
Thwdout = zeros(numel(t),1); 
Tchwout = zeros(numel(t),1); 
Tcwcout = zeros(numel(t),1);
for i=1:numel(t)
  t_actual =t(i);
  sol_actual = Soln(i,:);
  [~,dsol_actual]=ODEdw(t,sol_actual,iflag);
  Tcwaout(i) = dsol_actual(1);
  Thwdout(i) = dsol_actual(2);
  Tchwout(i) = dsol_actual(3);
  Tcwcout(i) = dsol_actual(4); 
end
Soln = horzcat(Soln,Tcwaout,Thwdout,Tchwout,Tcwcout);
colNames = {'Time','Wads','Tads','Wdes','Tdes','Mweva', 'Teva','Tcond',...
    'Tcwaout', 'Thwdout', 'Tchwout', 'Tcwcout'};
Result = array2table([tspan',Soln],'VariableNames',colNames);

function [dYdt] =ODEdw(~,Soln,iflag)
Ea=4.2*10^4;         %Activation energy of surface diffusion, (J/mol)
R=8.314;             %Ideal gas constant, (j/mol.K)
Dso=2.54*10^-4;       %Pre-exponential constant, (m^2 S^-1)
Qs=2.8*10^6;         %Heat of adsorption, (J/kg)

%%%%%%%%%%%Operating conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Thwdin=85;            %Inlet temperature of hot water in desorber, (oC)
Tcwain=27;            %Inlet temperature of cooling water in adsorber, (oC)
Tcwcin=27;            %Inlet temperature of cooling water in condenser, (oC)
Tchwin=14;            %Inlet temperature of chilling water in evaporator, (oC)
mcwa=8/60;            %Mass flow rate of cooling water in adsorber, (kg/s)
mcwc=20/60;           %Mass flow rate of cooling water in condenser, (kg/s)
mchw=1.6/60;          %Mass flow rate of chilling water in evaporator, (kg/s)
mhwd=7/60;            %Mass flow rate of hot water in desorber, (kg/s)

%%%%%%%%%Physical and thermal properties of adsorbent%%%%%%%%%%

Cps=924;              %Specific heat capacity of silica gel, (J/kg.K)
Ms=9;                 %Mass of silica, (kg)
Rp=5.6*10^-4;         %Radious of silica gel particle, (m)
ps=700;              %density of silica gel, (kg/m^3)
Por=0.3;              %Porosity of silica gel,



%%%%%%Physical and thermal properties of adsorbate%%%%%%%%%%%%%

Cpw=4200;             %Specific heat capcity of water, (J/kg.k)
Cpv=1900;             %Specific heat capacity of water vapor (J/kg.K)
M=18.02;              %Molar mass of water, (g/mol)
Lv=2500*10^3;         %Lantent heat of vaporization (J/kg) 

%%%%%%%%%Physical and thermal Properties of Condenser%%%%%%%%%%%%%

Acond=5.76;            %Surface area of condenser, (m^2)
Mhexcond=19.87;        %Mass of heat exchanger (Condenser), (kg)
Ucond=486.87;          %Heat transfer coefficient of condenser, (W/m^2.K)
UAcond=Acond*Ucond;
Mwcond=1.0;         %Mass of water in Condenser (kg)

%%%%%%%%%%Physical and thermal properties of adsorber/desorber   %%%%%%%

Mhex=65;             %Mass of heat exchanger (Adsorber/desorber), (kg)
Cphex=386;           %Heat capacity of heat excahgner, (J/kg.K)
Aad=2.11;            %Surface area of adsorber/desorber, (m^2)
Uads=189.60;         %Heat transfer coefficient of adsorber, (W/m^2.K)
Udes=203.99;         %Heat transfer coefficient of desorber, (W/m^2.K)
UAads=Aad*Uads;
UAdes=Aad*Udes;
Vt=0.5;              %volume of adsorber/desorber tank, (m^3)

%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%

Aeva=0.34;           %Surface area of evaporator, (m^2)
Mhexeva=16.1;        %Mass of heat exchanger (Evaporator), (kg)
Ueva=302.58;         %Heat transfer coefficient of evaporator, (W/m^2.K)
UAeva=Aeva*Ueva;

Wads=Soln(1);    
Tads=Soln(2);  
Wdes=Soln(3);    
Tdes=Soln(4);
Mweva=Soln(5);
Teva=Soln(6);
Tcond=Soln(7);
Kads=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tads+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwa=133.32*exp(18.3-(3820/((Teva+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)
Psga=133.32*exp(18.3-(3820/((Tads+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqa=0.346*(Pwa/Psga)^0.625;
Aads=(Kads*(Weqa-Wads));
TLads=(Ms*(Cpw+Cpw*Wads)+Mhex*Cphex); %Left hand side of adsorber equation
Tcwaout=Tads+(Tcwain-Tads)*exp((-UAads)/(mcwa*Cpw));  
Kdes=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tdes+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwd=133.32*exp(18.3-(3820/((Tcond+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)
Psgd=133.32*exp(18.3-(3820/((Tdes+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqd=0.346*(Pwd/Psgd)^0.625;
Ades=(Kdes*(Weqd-Wdes));
TLdes=(Ms*(Cps+Cpw*Wdes)+Mhex*Cphex); %Left hand side of desorber equation
Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(mhwd*Cpw));
TLeva=(Mweva*Cpw+Mhexeva*Cphex); %Left hand side of evaporator equation
Tchwout=Teva+(Tchwin-Teva)*exp((-UAeva)/(mchw*Cpw));  
TLcond=(Mwcond*Cpw+Mhexcond*Cphex); %Left hand side of evaporator equation
Tcwcout=Tcond+(Tcwcin-Tcond)*exp((-UAcond)/(mcwc*Cpw));
if iflag==1
  dYdt=[Kads*(Weqa-Wads);((Qs*Ms*Aads)+(Ms*Cpv*Aads*(Teva-Tads))+mcwa*Cpw*(Tcwain-Tcwaout))/TLads;
    Kdes*(Weqd-Wdes); ((Qs*Ms*Ades)+(mcwa*Cpw*(Thwdin-Thwdout)))/TLdes; 
    -Ms*((Kads*(Weqa-Wads))+(Kdes*(Weqd-Wdes)));((-Lv*Ms*Aads)+Ades*Ms*Cpw*(Tcond-Teva)+mchw*Cpw*(Tchwin-Tchwout))/TLeva;
    (Ades*Ms*Cpv*(Tdes-Tcond)+(Lv*Ms*Ades)+mcwc*Cpw*(Tcwcin-Tcwcout))/TLcond];
else
  dYdt = [Tcwaout; Thwdout; Tchwout; Tcwcout]
end
end
end

[Stephan]

Line 30 - horzcat error - fixed this way:

Soln = horzcat(Soln,Tcwaout',Thwdout',Tchwout',Tcwcout');

Line 24 - Too many output arguments error - fixed this way:

[dsol_actual]=ODEdw(t,sol_actual,iflag);

But the result seems not to be correct. All zeros except the first and the last line.

[Torsten]

This code works for me:

function test

%Initial conditions

Wads=0.045;

Tads=80;

Wdes=0.035;

Tdes=27;

Teva=7;

Tcond=27;

Mweva=0.65;         

Initial=[Wads Tads Wdes Tdes Mweva Teva Tcond];           

tspan=[1:1:800];

opts = odeset('RelTol',1e-6,'AbsTol',1e-6);

iflag = 1;

[t,Soln]=ode45(@(t,Soln)ODEdw(t,Soln,iflag),tspan,Initial,opts);



iflag = 2;

Tcwaout = zeros(numel(t),1); 

Thwdout = zeros(numel(t),1); 

Tchwout = zeros(numel(t),1); 

Tcwcout = zeros(numel(t),1);

for i=1:numel(t)

  t_actual =t(i);

  sol_actual = Soln(i,:);

  [dsol_actual]=ODEdw(t_actual,sol_actual,iflag);

  Tcwaout(i) = dsol_actual(1);

  Thwdout(i) = dsol_actual(2);

  Tchwout(i) = dsol_actual(3);

  Tcwcout(i) = dsol_actual(4); 

end

Soln = [Soln,Tcwaout,Thwdout,Tchwout,Tcwcout];

colNames = {'Time','Wads','Tads','Wdes','Tdes','Mweva', 'Teva','Tcond',...
    'Tcwaout', 'Thwdout', 'Tchwout', 'Tcwcout'};

Result = array2table([t,Soln],'VariableNames',colNames);

end

function [dYdt] =ODEdw(~,Soln,iflag)

Ea=4.2*10^4;         %Activation energy of surface diffusion, (J/mol)

R=8.314;             %Ideal gas constant, (j/mol.K)

Dso=2.54*10^-4;       %Pre-exponential constant, (m^2 S^-1)

Qs=2.8*10^6;         %Heat of adsorption, (J/kg)



%%%%%%%%%%%Operating conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%



Thwdin=85;            %Inlet temperature of hot water in desorber, (oC)

Tcwain=27;            %Inlet temperature of cooling water in adsorber, (oC)

Tcwcin=27;            %Inlet temperature of cooling water in condenser, (oC)

Tchwin=14;            %Inlet temperature of chilling water in evaporator, (oC)

mcwa=8/60;            %Mass flow rate of cooling water in adsorber, (kg/s)

mcwc=20/60;           %Mass flow rate of cooling water in condenser, (kg/s)

mchw=1.6/60;          %Mass flow rate of chilling water in evaporator, (kg/s)

mhwd=7/60;            %Mass flow rate of hot water in desorber, (kg/s)



%%%%%%%%%Physical and thermal properties of adsorbent%%%%%%%%%%



Cps=924;              %Specific heat capacity of silica gel, (J/kg.K)

Ms=9;                 %Mass of silica, (kg)

Rp=5.6*10^-4;         %Radious of silica gel particle, (m)

ps=700;              %density of silica gel, (kg/m^3)

Por=0.3;              %Porosity of silica gel,







%%%%%%Physical and thermal properties of adsorbate%%%%%%%%%%%%%



Cpw=4200;             %Specific heat capcity of water, (J/kg.k)

Cpv=1900;             %Specific heat capacity of water vapor (J/kg.K)

M=18.02;              %Molar mass of water, (g/mol)

Lv=2500*10^3;         %Lantent heat of vaporization (J/kg) 



%%%%%%%%%Physical and thermal Properties of Condenser%%%%%%%%%%%%%



Acond=5.76;            %Surface area of condenser, (m^2)

Mhexcond=19.87;        %Mass of heat exchanger (Condenser), (kg)

Ucond=486.87;          %Heat transfer coefficient of condenser, (W/m^2.K)

UAcond=Acond*Ucond;

Mwcond=1.0;         %Mass of water in Condenser (kg)



%%%%%%%%%%Physical and thermal properties of adsorber/desorber   %%%%%%%



Mhex=65;             %Mass of heat exchanger (Adsorber/desorber), (kg)

Cphex=386;           %Heat capacity of heat excahgner, (J/kg.K)

Aad=2.11;            %Surface area of adsorber/desorber, (m^2)

Uads=189.60;         %Heat transfer coefficient of adsorber, (W/m^2.K)

Udes=203.99;         %Heat transfer coefficient of desorber, (W/m^2.K)

UAads=Aad*Uads;

UAdes=Aad*Udes;

Vt=0.5;              %volume of adsorber/desorber tank, (m^3)



%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%



Aeva=0.34;           %Surface area of evaporator, (m^2)

Mhexeva=16.1;        %Mass of heat exchanger (Evaporator), (kg)

Ueva=302.58;         %Heat transfer coefficient of evaporator, (W/m^2.K)

UAeva=Aeva*Ueva;



Wads=Soln(1);    

Tads=Soln(2);  

Wdes=Soln(3);    

Tdes=Soln(4);

Mweva=Soln(5);

Teva=Soln(6);

Tcond=Soln(7);

Kads=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tads+273.16)));   %Overall mass tarsnfer coeffcient of adsorber

Pwa=133.32*exp(18.3-(3820/((Teva+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)

Psga=133.32*exp(18.3-(3820/((Tads+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)

Weqa=0.346*(Pwa/Psga)^0.625;

Aads=(Kads*(Weqa-Wads));

TLads=(Ms*(Cpw+Cpw*Wads)+Mhex*Cphex); %Left hand side of adsorber equation

Tcwaout=Tads+(Tcwain-Tads)*exp((-UAads)/(mcwa*Cpw));  

Kdes=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tdes+273.16)));   %Overall mass tarsnfer coeffcient of adsorber

Pwd=133.32*exp(18.3-(3820/((Tcond+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)

Psgd=133.32*exp(18.3-(3820/((Tdes+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)

Weqd=0.346*(Pwd/Psgd)^0.625;

Ades=(Kdes*(Weqd-Wdes));

TLdes=(Ms*(Cps+Cpw*Wdes)+Mhex*Cphex); %Left hand side of desorber equation

Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(mhwd*Cpw));

TLeva=(Mweva*Cpw+Mhexeva*Cphex); %Left hand side of evaporator equation

Tchwout=Teva+(Tchwin-Teva)*exp((-UAeva)/(mchw*Cpw));  

TLcond=(Mwcond*Cpw+Mhexcond*Cphex); %Left hand side of evaporator equation

Tcwcout=Tcond+(Tcwcin-Tcond)*exp((-UAcond)/(mcwc*Cpw));

if iflag==1

  dYdt=[Kads*(Weqa-Wads);((Qs*Ms*Aads)+(Ms*Cpv*Aads*(Teva-Tads))+mcwa*Cpw*(Tcwain-Tcwaout))/TLads;

    Kdes*(Weqd-Wdes); ((Qs*Ms*Ades)+(mcwa*Cpw*(Thwdin-Thwdout)))/TLdes; 

    -Ms*((Kads*(Weqa-Wads))+(Kdes*(Weqd-Wdes)));((-Lv*Ms*Aads)+Ades*Ms*Cpw*(Tcond-Teva)+mchw*Cpw*(Tchwin-Tchwout))/TLeva;

    (Ades*Ms*Cpv*(Tdes-Tcond)+(Lv*Ms*Ades)+mcwc*Cpw*(Tcwcin-Tcwcout))/TLcond];

else

  dYdt = [Tcwaout; Thwdout; Tchwout; Tcwcout]

end

end

[Md Shamim Ahamed]

Hi

Thanks for your time to help me. You code is working correctly. But, can you please inform me that why RESULT is not showing in Workspace?

kind reagrds

[Stephan]

Use:

Result = Main
function Result = Main
%Initial conditions
Wads=0.045;
Tads=80;
Wdes=0.035;
Tdes=27;
Teva=7;
Tcond=27;
Mweva=0.65;         
Initial=[Wads Tads Wdes Tdes Mweva Teva Tcond];           
tspan=[1:1:800];
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
iflag = 1;
[t,Soln]=ode45(@(t,Soln)ODEdw(t,Soln,iflag),tspan,Initial,opts);
iflag = 2;
Tcwaout = zeros(numel(t),1); 
Thwdout = zeros(numel(t),1); 
Tchwout = zeros(numel(t),1); 
Tcwcout = zeros(numel(t),1);
for i=1:numel(t)
  t_actual =t(i);
  sol_actual = Soln(i,:);
  [dsol_actual]=ODEdw(t_actual,sol_actual,iflag);
  Tcwaout(i) = dsol_actual(1);
  Thwdout(i) = dsol_actual(2);
  Tchwout(i) = dsol_actual(3);
  Tcwcout(i) = dsol_actual(4); 
end
Soln = [Soln,Tcwaout,Thwdout,Tchwout,Tcwcout];
colNames = {'Time','Wads','Tads','Wdes','Tdes','Mweva', 'Teva','Tcond',...
    'Tcwaout', 'Thwdout', 'Tchwout', 'Tcwcout'};
Result = array2table([t,Soln],'VariableNames',colNames);
end
function [dYdt] =ODEdw(~,Soln,iflag)
Ea=4.2*10^4;         %Activation energy of surface diffusion, (J/mol)
R=8.314;             %Ideal gas constant, (j/mol.K)
Dso=2.54*10^-4;       %Pre-exponential constant, (m^2 S^-1)
Qs=2.8*10^6;         %Heat of adsorption, (J/kg)
%%%%%%%%%%%Operating conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Thwdin=85;            %Inlet temperature of hot water in desorber, (oC)
Tcwain=27;            %Inlet temperature of cooling water in adsorber, (oC)
Tcwcin=27;            %Inlet temperature of cooling water in condenser, (oC)
Tchwin=14;            %Inlet temperature of chilling water in evaporator, (oC)
mcwa=8/60;            %Mass flow rate of cooling water in adsorber, (kg/s)
mcwc=20/60;           %Mass flow rate of cooling water in condenser, (kg/s)
mchw=1.6/60;          %Mass flow rate of chilling water in evaporator, (kg/s)
mhwd=7/60;            %Mass flow rate of hot water in desorber, (kg/s)
%%%%%%%%%Physical and thermal properties of adsorbent%%%%%%%%%%
Cps=924;              %Specific heat capacity of silica gel, (J/kg.K)
Ms=9;                 %Mass of silica, (kg)
Rp=5.6*10^-4;         %Radious of silica gel particle, (m)
ps=700;              %density of silica gel, (kg/m^3)
Por=0.3;              %Porosity of silica gel,
%%%%%%Physical and thermal properties of adsorbate%%%%%%%%%%%%%
Cpw=4200;             %Specific heat capcity of water, (J/kg.k)
Cpv=1900;             %Specific heat capacity of water vapor (J/kg.K)
M=18.02;              %Molar mass of water, (g/mol)
Lv=2500*10^3;         %Lantent heat of vaporization (J/kg) 
%%%%%%%%%Physical and thermal Properties of Condenser%%%%%%%%%%%%%
Acond=5.76;            %Surface area of condenser, (m^2)
Mhexcond=19.87;        %Mass of heat exchanger (Condenser), (kg)
Ucond=486.87;          %Heat transfer coefficient of condenser, (W/m^2.K)
UAcond=Acond*Ucond;
Mwcond=1.0;         %Mass of water in Condenser (kg)
%%%%%%%%%%Physical and thermal properties of adsorber/desorber   %%%%%%%
Mhex=65;             %Mass of heat exchanger (Adsorber/desorber), (kg)
Cphex=386;           %Heat capacity of heat excahgner, (J/kg.K)
Aad=2.11;            %Surface area of adsorber/desorber, (m^2)
Uads=189.60;         %Heat transfer coefficient of adsorber, (W/m^2.K)
Udes=203.99;         %Heat transfer coefficient of desorber, (W/m^2.K)
UAads=Aad*Uads;
UAdes=Aad*Udes;
Vt=0.5;              %volume of adsorber/desorber tank, (m^3)
%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%
Aeva=0.34;           %Surface area of evaporator, (m^2)
Mhexeva=16.1;        %Mass of heat exchanger (Evaporator), (kg)
Ueva=302.58;         %Heat transfer coefficient of evaporator, (W/m^2.K)
UAeva=Aeva*Ueva;
Wads=Soln(1);    
Tads=Soln(2);  
Wdes=Soln(3);    
Tdes=Soln(4);
Mweva=Soln(5);
Teva=Soln(6);
Tcond=Soln(7);
Kads=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tads+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwa=133.32*exp(18.3-(3820/((Teva+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)
Psga=133.32*exp(18.3-(3820/((Tads+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqa=0.346*(Pwa/Psga)^0.625;
Aads=(Kads*(Weqa-Wads));
TLads=(Ms*(Cpw+Cpw*Wads)+Mhex*Cphex); %Left hand side of adsorber equation
Tcwaout=Tads+(Tcwain-Tads)*exp((-UAads)/(mcwa*Cpw));  
Kdes=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tdes+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwd=133.32*exp(18.3-(3820/((Tcond+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)
Psgd=133.32*exp(18.3-(3820/((Tdes+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqd=0.346*(Pwd/Psgd)^0.625;
Ades=(Kdes*(Weqd-Wdes));
TLdes=(Ms*(Cps+Cpw*Wdes)+Mhex*Cphex); %Left hand side of desorber equation
Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(mhwd*Cpw));
TLeva=(Mweva*Cpw+Mhexeva*Cphex); %Left hand side of evaporator equation
Tchwout=Teva+(Tchwin-Teva)*exp((-UAeva)/(mchw*Cpw));  
TLcond=(Mwcond*Cpw+Mhexcond*Cphex); %Left hand side of evaporator equation
Tcwcout=Tcond+(Tcwcin-Tcond)*exp((-UAcond)/(mcwc*Cpw));
if iflag==1
  dYdt=[Kads*(Weqa-Wads);((Qs*Ms*Aads)+(Ms*Cpv*Aads*(Teva-Tads))+mcwa*Cpw*(Tcwain-Tcwaout))/TLads;
    Kdes*(Weqd-Wdes); ((Qs*Ms*Ades)+(mcwa*Cpw*(Thwdin-Thwdout)))/TLdes; 
    -Ms*((Kads*(Weqa-Wads))+(Kdes*(Weqd-Wdes)));((-Lv*Ms*Aads)+Ades*Ms*Cpw*(Tcond-Teva)+mchw*Cpw*(Tchwin-Tchwout))/TLeva;
    (Ades*Ms*Cpv*(Tdes-Tcond)+(Lv*Ms*Ades)+mcwc*Cpw*(Tcwcin-Tcwcout))/TLcond];
else
  dYdt = [Tcwaout; Thwdout; Tchwout; Tcwcout];
end
end

[Md Shamim Ahamed]

Thanks Toresten and Stephan

Answer 3

[Md Shamim Ahamed]

Hi Toresten and Stepphan

Thanks for your help. The code was working perfectly. Now, I am facing a problem after changing the order of ODE euqations. Can you please help where is the problem?

I have attached both codes (Changed and Original) with this comment.

Kind Reagrds

Answer 4

[Md Shamim Ahamed]

Hi

now, its works. i made a sily mistake

thanks