One of the constant input variable need to be changed based on previous step's output in ODE, Can anyone kindly help?

The following code (I have attached below) has been solved by a constant hot water input (Thwdin=85) to see a chambers temperature (Tdes) change at different time steps. A linear equation provides the outlet water temperature Thwdout. This is an open system now.

I want to make this system a recirculating system as in each previous Thwdout will be the new input water temperature (Thwdin) in the chamber (Tdes).

Can It be solved using ODE45? If yes then how? IF no, then can anybody please suggest what approach I should follow?

Sample picture of the current situation and my expected situation is given below,

clear all, close all, clc
global Thwdin Tcwain Tcwcin Tchwin mcwa mcwc mchw mhwd Cps Ms Rp Cpw Cpv Lv Acond Mhexcond Ucond Mwcond Mhex Cphex Aad Uads Udes Aeva Mhexeva Ueva tcycle tphase thr

%%%%%%%%%%%Operating conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Thwdin=85;            %Inlet temperature of hot water in desorber, (oC)
Tcwain=27;            %Inlet temperature of cooling water in adsorber, (oC)
Tcwcin=Tcwain;        %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)

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

%%%%%%%%%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)
tcycle=900;

Result = main;

Result = 900×6 table
    Time      Wdes       Tdes     Tcond     Thwdout    Tcwcout
    ____    ________    ______    ______    _______    _______

      1        0.065        30        28    52.849     27.865 
      2     0.065182    30.565    28.232     53.18     28.066 
      3     0.065367    31.127    28.446    53.508     28.251 
      4     0.065555    31.687    28.645    53.835     28.423 
      5     0.065745    32.243    28.828    54.161     28.582 
      6     0.065937    32.797    28.999    54.485     28.729 
      7     0.066132    33.348    29.157    54.807     28.866 
      8     0.066328    33.896    29.303    55.127     28.992 
      9     0.066527    34.441    29.439    55.445      29.11 
     10     0.066728    34.983    29.566    55.762      29.22 
     11      0.06693    35.521    29.684    56.077     29.322 
     12     0.067135    36.057    29.793     56.39     29.416 
     13     0.067341    36.589    29.895    56.701     29.505 
     14     0.067548    37.118     29.99     57.01     29.587 
     15     0.067757    37.644    30.079    57.318     29.663 
     16     0.067967    38.166    30.161    57.623     29.735 

Time=Result{:,1}; Tdes=Result{:,3}; Tcond=Result{:,4};Thwdout=Result{:,5};Tcwcout=Result{:,6};

%%%%%%%%%%%%%%%%%Plot%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot(Time,Tdes,'--b',Time,Tcond,'m',Time,Thwdout,'b',Time,Tcwcout,'c');  
legend('Tdes','Tcond','Thwdout','Tcwcout')
xlabel('Time (s)');ylabel('Temperature (^oC)')

function Result = main
global tcycle

%Initial conditions
Wdes=0.065; Tdes=30; Tcond=28;
        
Initial=[Wdes Tdes Tcond];        
tspan=1:1:tcycle;
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
iflag = 1;
[t,Soln]=ode45(@(t,Soln)ODEdw(t,Soln,iflag),tspan,Initial,opts);
iflag = 2;
Thwdout = zeros(numel(t),1);  
Tcwcout = zeros(numel(t),1);

    for i=1:numel(t)
      t_actual =t(i);
      Soln_actual = Soln(i,:);
      [dsoln_actual]=ODEdw(t_actual,Soln_actual,iflag);
      Thwdout(i) = dsoln_actual(1);
      Tcwcout(i) = dsoln_actual(2); 
    end

Soln = [Soln,Thwdout,Tcwcout];
colNames = {'Time','Wdes','Tdes','Tcond','Thwdout','Tcwcout'};
Result = array2table([t,Soln],'VariableNames',colNames)

end

function [dYdt] =ODEdw(t,Soln,iflag)
global Thwdin Tcwain Tcwcin Tchwin mcwa mcwc mchw mhwd Cps Ms Rp Cpw Cpv Lv Acond Mhexcond Ucond Mwcond Mhex Cphex Aad Uads Udes Aeva Mhexeva Ueva tphase thr tcycle
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)


%%%%%%%Estimated Paramters%%%%%%%%%%%%%%%%%%%%%%%%%%%
UAcond=Acond*Ucond;
UAdes=Aad*Udes;

%%%%%%%%%%%%%%%Estimated Parameters%%%%%%%%%%%
Wdes=Soln(1);Tdes=Soln(2);Tcond=Soln(3);

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)));    % Vapour pressure Inside of the adsorber, (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));
        
TLcond=(Mwcond*Cpw+Mhexcond*Cphex);                   %Left hand side of evaporator equation
Tcwcout=Tcond+(Tcwcin-Tcond)*exp((-UAcond)/(mcwc*Cpw));

         if iflag==1
           dYdt=[(Kdes*(Weqd-Wdes)); ((-Qs*Ms*Ades)+(mhwd*Cpw*(Thwdin-Thwdout)))/TLdes;
           (Ades*Ms*Cpv*(Tcond-Tdes)-(Lv*Ms*Ades)+mcwc*Cpw*(Tcwcin-Tcwcout))/TLcond];
         else
           dYdt = [Thwdout;Tcwcout];
         end
end

Answer 1

[Cris LaPierre]

I don't think your approach to solving for Thwdout and Tcwcout are doing what you want. You run ode45 a second time as a way of computing Thwdout and Tcwcout, not for solving an ode. It is unnecessary. You could simply pull the calculation out and solve for it after.

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

%%%%%%%Estimated Paramters%%%%%%%%%%%%%%%%%%%%%%%%%%%

UAcond=param.Acond*param.Ucond;

UAdes=param.Aad*param.Udes;

Thwdout=Soln(:,2)+(param.Thwdin-Soln(:,2))*exp((-UAdes)/(param.mhwd*param.Cpw));

Tcwcout=Soln(:,3)+(param.Tcwcin-Soln(:,3))*exp((-UAcond)/(param.mcwc*param.Cpw));

However, to do what you are asking, the calculation does need to be part of your odefun. To work, then, you need to express your calculations as differential equations. Fortunately, you already have that in your odefun (differential equations for Tdes and Tcond).

Use the product rule to take the derivative of Thwdout and Tcwcout and adding them to dYdt leads to the following expressoin for dYdt:

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

dy2 = ((-Qs*param.Ms*Ades)+(param.mhwd*param.Cpw*(Thwdin-Thwdout)))/TLdes;

dy3 = (Ades*param.Ms*param.Cpv*(Tcond-Tdes)-(param.Lv*param.Ms*Ades)+param.mcwc*param.Cpw*(param.Tcwcin-Tcwcout))/TLcond;

dy4 = dy2 - dy2*exp((-UAdes)/(param.mhwd*param.Cpw));

dy5 = dy3 - dy3*exp((-UAcond)/(param.mcwc*param.Cpw));

dYdt = [dy1; dy2; dy3; dy4; dy5];

Now specify Thwcin as Soln(4) and specify your initial conditions, and you should have what you want. Big note - I am not an expert, so I will defer to others who may critique this approach. Here is my updates to your code:

% https://www.mathworks.com/matlabcentral/answers/1620675-one-of-the-constant-input-variable-need-to-be-changed-based-on-previous-step-s-output-in-ode-can-an?s_tid=srchtitle
%%%%%%%%%%%Operating conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param.Thwdin=85;            %Inlet temperature of hot water in desorber, (oC)
param.Tcwain=27;            %Inlet temperature of cooling water in adsorber, (oC)
param.Tcwcin=param.Tcwain;        %Inlet temperature of cooling water in condenser, (oC)
param.Tchwin=14;            %Inlet temperature of chilling water in evaporator, (oC)
param.mcwa=8/60;            %Mass flow rate of cooling water in adsorber, (kg/s)
param.mcwc=20/60;           %Mass flow rate of cooling water in condenser, (kg/s)
param.mchw=1.6/60;          %Mass flow rate of chilling water in evaporator, (kg/s)
param.mhwd=7/60;            %Mass flow rate of hot water in desorber, (kg/s)

%%%%%%%%%Physical and thermal properties of adsorbent%%%%%%%%%%
param.Cps=924;              %Specific heat capacity of silica gel, (J/kg.K)
param.Ms=9;                 %Mass of silica, (kg)
param.Rp=5.6*10^-4;         %Radious of silica gel particle, (m)

%%%%%%Physical and thermal properties of adsorbate%%%%%%%%%%%%%
param.Cpw=4200;             %Specific heat capcity of water, (J/kg.k)
param.Cpv=1900;             %Specific heat capacity of water vapor (J/kg.K)
param.Lv=2500*10^3;         %Lantent heat of vaporization (J/kg)

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

%%%%%%%%%%Physical and thermal properties of adsorber/desorber   %%%%%%%
param.Mhex=65;              %Mass of heat exchanger (Adsorber/desorber), (kg)
param.Cphex=386;            %Heat capacity of heat excahgner, (J/kg.K)
param.Aad=2.11;             %Surface area of adsorber/desorber, (m^2)
param.Uads=189.60;          %Heat transfer coefficient of adsorber, (W/m^2.K)
param.Udes=203.99;          %Heat transfer coefficient of desorber, (W/m^2.K)

%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%
param.Aeva=0.34;            %Surface area of evaporator, (m^2)
param.Mhexeva=16.1;         %Mass of heat exchanger (Evaporator), (kg)
param.Ueva=302.58;          %Heat transfer coefficient of evaporator, (W/m^2.K)
tcycle=900;

Result2 = main(tcycle,param);

Result = 900×6 table
    Time      Wdes       Tdes     Tcond     Thwdout    Tcwcout
    ____    ________    ______    ______    _______    _______

      1        0.065        30        28    52.849     27.865 
      2     0.065182     30.31    28.231    53.031     28.065 
      3     0.065366     30.62    28.445    53.212      28.25 
      4     0.065553    30.931    28.641    53.394     28.419 
      5     0.065741    31.242    28.822    53.575     28.576 
      6     0.065931    31.553    28.989    53.757     28.721 
      7     0.066123    31.865    29.143     53.94     28.854 
      8     0.066317    32.176    29.285    54.122     28.977 
      9     0.066513    32.488    29.417    54.304     29.091 
     10      0.06671      32.8    29.538    54.486     29.196 
     11     0.066908    33.111    29.651    54.668     29.293 
     12     0.067108    33.423    29.755     54.85     29.383 
     13     0.067309    33.734    29.851    55.032     29.466 
     14     0.067512    34.045     29.94    55.214     29.544 
     15     0.067715    34.356    30.023    55.396     29.615 
     16      0.06792    34.666      30.1    55.577     29.682 

Time=Result2.Time; Tdes=Result2.Tdes; Tcond=Result2.Tcond;Thwdout=Result2.Thwdout;Tcwcout=Result2.Tcwcout;

%%%%%%%Estimated Paramters%%%%%%%%%%%%%%%%%%%%%%%%%%%
UAcond=param.Acond*param.Ucond;
UAdes=param.Aad*param.Udes;

%%%%%%%%%%%%%%%%%Plot%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(Time,Tdes,'--b',Time,Tcond,'m',Time,Thwdout,'b',Time,Tcwcout,'c');
legend('Tdes','Tcond','Thwdout','Tcwcout')
xlabel('Time (s)');ylabel('Temperature (^oC)')

function Result = main(tcycle,param)
%Initial conditions
Wdes=0.065; Tdes=30; Tcond=28; 
Thwdout = Tdes+(param.Thwdin-Tdes)*exp(-param.Aad*param.Udes/(param.mhwd*param.Cpw));
Tcwcout=Tcond+(param.Tcwcin-Tcond)*exp(-param.Acond*param.Ucond/(param.mcwc*param.Cpw));

Initial=[Wdes Tdes Tcond Thwdout Tcwcout];
tspan=1:1:tcycle;
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,Soln]=ode45(@ODEdw,tspan,Initial,opts,param);

colNames = {'Time','Wdes','Tdes','Tcond','Thwdout','Tcwcout'};
Result = array2table([t,Soln],'VariableNames',colNames)
end

function [dYdt] =ODEdw(t,Soln,param)
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)

%%%%%%%Estimated Paramters%%%%%%%%%%%%%%%%%%%%%%%%%%%
UAcond=param.Acond*param.Ucond;
UAdes=param.Aad*param.Udes;

%%%%%%%%%%%%%%%Estimated Parameters%%%%%%%%%%%
Wdes=Soln(1);Tdes=Soln(2);Tcond=Soln(3);
Thwdin = Soln(4);

Kdes=((15*Dso)/(param.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)));    % Vapour pressure Inside of the adsorber, (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=(param.Ms*(param.Cps+param.Cpw*Wdes)+param.Mhex*param.Cphex);                 %Left hand side of desorber equation
Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(param.mhwd*param.Cpw));

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

dy1 = (Kdes*(Weqd-Wdes));
dy2 = ((-Qs*param.Ms*Ades)+(param.mhwd*param.Cpw*(Thwdin-Thwdout)))/TLdes;
dy3 = (Ades*param.Ms*param.Cpv*(Tcond-Tdes)-(param.Lv*param.Ms*Ades)+param.mcwc*param.Cpw*(param.Tcwcin-Tcwcout))/TLcond;
dy4 = dy2 - dy2*exp((-UAdes)/(param.mhwd*param.Cpw));
dy5 = dy3 - dy3*exp((-UAcond)/(param.mcwc*param.Cpw));

dYdt = [dy1; dy2; dy3; dy4; dy5];
end

Just change Thwdin = Soln(4); to Thwdin = param.Thwdin; to get your original solution.

[TM Abir Ahsan]

Thanks for your time and dedication @Cris LaPierre. However, I I dont think this solution works as the Thwdin increases the temperature of the of Tdes, however if you take a look at the results at t=2, Thwdout temperature is 53 and Tdes=30.31, if Thwdout works as Thwdin in coming steps, there is no chance for Tdes or Thwdout to reach a temperature of 80C at any point down the line. The actual system had continuous input of 85C water so the Tdes is likely to increase its temperature equivalent to that of Thwdin at some point down the line. In the recirculation case Tdes will never reach any temperature near 80C. Unfortunately your solution is showing that temperature reaches to 78C.

@Torsten Hopp Can you please provide me a reference for DD23 used in similar model calcualtion? I am unfamiliar with this MATLAB algorithm.

[Cris LaPierre]

I see what you mean. It is calculating Thwdout correctly inside the odefun, but the differential equation solution for Thwdout and Tcwcout needs to be delayed one time step (at least). Said another way, the equation is using the previous value of Tdes to compute Thwdout while the differential equation will be solved using the current value of Tdes.

You might find this helpful: Tutorial on solving DDEs with DDE23

[TM Abir Ahsan]

Thank you so much mate. @Cris LaPierre

Answer 2

[Torsten]

Either you insert some delay elements (PT1 or similar) or you have to solve a delay differential equation (e.g. using dde23).

[TM Abir Ahsan]

@Torsten Hopp Can you please provide me a reference for DD23 used in similar model calcualtion? I am unfamiliar with this MATLAB algorithm.

[Torsten]

Sorry, but I'm also unfamiliar with it. You will have to start with the MATLAB examples provided.