1 equation from a system of coupled ODE provides output which is used as input in the other equation of the coupled ODE.

2 views (last 30 days)
TM Abir Ahsan
TM Abir Ahsan on 28 Dec 2021
Edited: Torsten on 28 Dec 2021
Hotwater enters in one chamber and the temperature of the water at outlet is found by solving an ODE and a linear equation. I am trying to solve the equation by system of ode and system of linear equation how ever I think my code is wrong. Can anybody suggest me what I can do? Here is my code.
I have added a image for detailed visualization.
clear all
close all
clc
global Cps Ms Cpw Mhex Cphex Qs Ea R Dso Qs Uads Udes Aad m Rp tcycle Thwdin Tcond Teva
Cps=924; %Specific heat capacity of silica gel, (J/kg.K)
Ms=9; %Mass of silica, (kg)
Cpw=4200; %Specific heat capcity of water, (J/kg.k)
Mhex=65; %Mass of heat exchanger (Adsorber/desorber), (kg)
Cphex=386; %Heat capacity of heat excahgner, (J/kg.K)
Qs=2.8*10^6; %Heat of adsorption (+)/heat of desorption (-), (J/kg)
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)
Uads=189.60; %Heat transfer coefficient of adsorber, (W/m^2.K)
Udes=203.99; %Heat transfer coefficient of desorber, (W/m^2.K)
Aad=2.11; %Surface area of adsorber/desorber, (m^2)
m=7/60; %Flow rate of circulating water (kg/s)
Rp=5.6*10^-4; %Radious of silica gel particle, (m)
tcycle=30;
Thwdin=45;
Tcond=26.857;
Teva=4.4621;
global Cps Ms Cpw Mhex Cphex Qs Ea R Dso Qs Uads Udes Aad m Rp tcycle Thwdin Tcond Teva
%Initial conditions
Tads=29.598;
Tdes=83.936;
Wads=0.086543;
Wdes=0.062419;
Initial=[Wdes Tdes Wads Tads];
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);
Tcwaout = 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);
Tcwaout(i) = dsoln_actual(2);
Tcwain(i) = dsoln_actual(3);
end
Tcwainout=transpose(Tcwain);
Soln = [Soln,Thwdout,Tcwaout,Tcwainout];
colNames = {'Time','Wdes','Tdes','Wads','Tads','Thwdout', 'Tcwaout','Tcwain'};
Result = array2table([t,Soln],'VariableNames',colNames)
plot(t,Soln(:,2),'m',t,Soln(:,4),'g',t,Thwdout,'-p',t,Tcwaout,'--r',t,Tcwainout,'r')
hold on
for i=1:tcycle
Thwdin__out(i)=Thwdin;
end
plot(t,Thwdin__out,'b')
legend('TAdsorber1','TAdsorber2','TWater-Ad1-out','TWater-Ad2-out','TWater-Ad2-in','TWater-Ad1-in')
xlabel('Time (s)');
ylabel('Temperature (^oC)')
xticks(0:10:tcycle)
function [dYdt] =ODEdw(t,Soln,iflag)
global Cps Ms Cpw Mhex Cphex Qs Ea R Dso Qs Uads Udes Aad m Rp tcycle Thwdin Tcond Teva
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%%%%%%%%%%%%%%%%%%%%%%%%%%%
UAads=Aad*Uads;
UAdes=Aad*Udes;
%%%%%%%%%%%%%%%Estimated Parameters%%%%%%%%%%%
Wdes=Soln(1); Tdes=Soln(2); Wads=Soln(3); Tads=Soln(4);
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)/(m*Cpw));
Tcwain=Thwdout;
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))); %Vapour pressure inside of the adsorber, (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*(Cps+Cpw*Wads)+Mhex*Cphex); %Left hand side of adsorber equation
Tcwaout=Tads+(Thwdout-Tads)*exp((-UAads)/(m*Cpw));
if iflag==1
dYdt=[Kdes*(Weqd-Wdes); ((-Qs*Ms*Ades)+(m*Cpw*(Thwdin-Thwdout)))/TLdes;
Kads*(Weqa-Wads); (((Qs*Ms*Aads)+(m*Cpw*(Thwdout-Tcwaout)))/TLads)]
else
dYdt=[Thwdout Tcwaout Tcwain]
end
end
  1 Comment
Torsten
Torsten on 28 Dec 2021
So what exactly is your problem ?
From what you write, the code works, but you get unsatisfactory results. Then it will be difficult for us to help since the problem is related to your model, not related to MATLAB.
If you try to include an ODE to use THWDOUT for calculating TCWAOUT, you can do this in your function ODEdw.

Sign in to comment.

Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by