I am working on simulating an AC just like what's been shown in this link. I have managed to convert it into an AC by flipping some equations and re-initializing some values. Now, I need to make it in Matlab to able to let my system work smoothly and faster without having to run simulink for every iteration. I did some reading on ODE45 and decided to use it, firstly, I converted the model in the webpage to a function as the following:
function dTroom_dt = odes_thermal(Troom,Tout)
TinIC = 26;
r2d = 180/pi;
% -------------------------------
% Define the house geometry
% -------------------------------
% House length = 30 m
lenHouse = 30;
% House width = 10 m
widHouse = 10;
% House height = 4 m
htHouse = 4;
% Roof pitch = 40 deg
pitRoof = 40/r2d;
% Number of windows = 6
numWindows = 6;
% Height of windows = 1 m
htWindows = 1;
% Width of windows = 1 m
widWindows = 1;
windowArea = numWindows*htWindows*widWindows;
wallArea = 2*lenHouse*htHouse + 2*widHouse*htHouse + ...
2*(1/cos(pitRoof/2))*widHouse*lenHouse + ...
tan(pitRoof)*widHouse - windowArea;
% -------------------------------
% Define the type of insulation used
% -------------------------------
% Glass wool in the walls, 0.2 m thick
% k is in units of J/sec/m/C - convert to J/hr/m/C multiplying by 3600
kWall = 0.038*3600; % hour is the time unit
LWall = .2;
RWall = LWall/(kWall*wallArea);
% Glass windows, 0.01 m thick
kWindow = 0.78*3600; % hour is the time unit
LWindow = .01;
RWindow = LWindow/(kWindow*windowArea);
% -------------------------------
% Determine the equivalent thermal resistance for the whole building
% -------------------------------
Req = RWall*RWindow/(RWall + RWindow);
% c = cp of air (273 K) = 1005.4 J/kg-K
c = 1005.4;
% -------------------------------
% Enter the temperature of the heated air
% -------------------------------
% The air exiting the heater has a constant temperature which is a heater
% property. THeater =20 deg C
THeater = 20;
% Air flow rate Mdot = 1 kg/sec = 3600 kg/hr
Mdot = 3600; % hour is the time unit
% -------------------------------
% Determine total internal air mass = M
% -------------------------------
% Density of air at sea level = 1.2250 kg/m^3
densAir = 1.2250;
M = (lenHouse*widHouse*htHouse+tan(pitRoof)*widHouse*lenHouse)*densAir;
%%
dQlosses_dt= (THeater-Troom)*Mdot*c;
dQlosses_dt=(Troom-Tout)/Req;
dTroom_dt=(1/(M*c))*(dQlosses_dt);
end
Then I just called it using the following lines:
clear;
t=1:1:1440;
Troom_0=25;
[t,Troom] = ode45('odes_thermal',[t(1) t(3)],Troom_0);
figure(1);
plot(t,Troom);
when I run it this is what I get:
My question is, how do I incorporate the status or on off features in it. In simulink they were using a relay that will produce a status signal that would be 1 or 0 depending on the temperature, this is how they did it:
I know that this can be represented by "dQlosses_dt= (THeater-Troom)Mdotc*status;" but I don't know how to update status for every time sample as there is no way to check Troom at every time slot.Can you guide me in this, or at least from where to start?
Edit: following Pacta's comment, I tried to make it into a multivariant system which has two outputs Troom and Qlosses. Here are the new functions, but I am getting an error (too many input arguments in odes_thermal2)
clear;
t=1:1:1440;
Troom_0=25;Qlosses_0=10;
X0=[Troom_0;Qlosses_0];
[t,X] = ode45('odes_thermal2',[t(1) t(end)],X0);
Troom=X(1);
Qlosses=X(2);
figure(1);
plot(t,Troom);
ylabel('Troom');
xlabel('t');
And the ODE function has the following changes:
function dX_dt = odes_thermal2(X)
TinIC = 26;
r2d = 180/pi;
and
X(2)= (THeater-Troom)*Mdot*c;
dQlosses_dt=(Troom-Tout)/Req;
X(1)=(1/(M*c))*(dQlosses_dt);
Thanks
