# 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 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)
Udes=203.99; %Heat transfer coefficient of desorber, (W/m^2.K)
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
Tdes=83.936;
Wdes=0.062419;
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];
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')
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)
%%%%%%%Estimated Paramters%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%Estimated Parameters%%%%%%%%%%%
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;
TLdes=(Ms*(Cps+Cpw*Wdes)+Mhex*Cphex); %Left hand side of desorber equation
Tcwain=Thwdout;
Pwa=133.32*exp(18.3-(3820/((Teva+273.16)-46.1))); %Vapour pressure inside of the adsorber, (Pa)
Weqa=0.346*(Pwa/Psga)^0.625;
if iflag==1
else
dYdt=[Thwdout Tcwaout Tcwain]
end
end
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.