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

8 views (last 30 days)
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

Accepted Answer

Cris LaPierre
Cris LaPierre on 2 Jan 2022
Edited: Cris LaPierre on 3 Jan 2022
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.
  3 Comments

Sign in to comment.

More Answers (1)

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

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by