Having trouble modeling temperature flux through a building

I am trying to model the temperature flux through a simple two story cabin but I keep running into issues either with preallocating my temperature arrays or in my loop to model the temperature with Euler's method. The model I'm trying to make is of a cabin that has a fireplace that outputs a certain amount of energy on its first floor. There are three flux equations: the flux from the 1st floor to the outside, from the 1st floor to the 2nd floor, and from the second floor to the outside. These flux equations use the Newton Law of cooling to create two 1st order ODEs for the temperatures for the 1st and 2nd floors.

The flux equations are:

The temperature ODEs are: and

In my code I have included all my constants shown here in the equations.

% Identify Constants

k1 = 0.2; % W/(m^2*K)

k2 = 3; % range from 1 to 10 W/(m^2*K)

k3 = 0.5; % W/(m^2*K)

Qin = 2500; % Max of 2500, W/(m^2*k)

Ain = 0.5; % m^2

A1 = 30; % Square footage of the 1st floor, m^2

A2 = 30; % Square footage of 2nd floor, m^2

rho_air = 1.2; % kg/m^3

Cair = 1005; % 1005 J/(kg*K)

V1 = 90; % volume of lower floor

V2 = 38.97; % volume of upper floor

h = 1; % stepsize

t = linspace(0,2*86400);

N = t/h;

% Initialize Arrays

T1 = zeros(1,N);

T2 = zeros(1,N);

% IC

T1(1) = 5; % C

T2(1) = 7; % C

T_out = -10*sin((2*pi.*t)/86400); % C

% Functions

Q1 = k1*(T1 - T_out);

Q2 = k2*(T1 - T2 + 5);

Q3 = k3*(T2 - T_out);

dT1 = (Ain*Qin)/(Cair*rho_air*V1) - (A1*Q1)/(Cair*rho_air*V1) + (A1*Q2)/(Cair*rho_air*V1);

dT2 = (Ain*Qin)/(Cair*rho_air*V1) + (A1*Q2)/(Cair*rho_air*V1) + (A2*Q3)/(Cair*rho_air*V2);

for i = 1:N-1

t(i+1) = t(i) + h;

T1(i+1) = T1(i) + h*dT1(T1(i),t(i));

T2(i+1) = T2(i) + h*dT2(T2(i),t(i));

end

plot(t,T1,'r-')

plot(t,T2,'g-')

Any help with getting this thing going would be greatly appreciated.

Answer 1

[Davide Masiello]

The following code works

% Identify Constants
k1 = 0.2; % W/(m^2*K)
k2 = 3; % range from 1 to 10 W/(m^2*K)
k3 = 0.5; % W/(m^2*K)
Qin = 2500; % Max of 2500, W/(m^2*k)
Ain = 0.5; % m^2
A1 = 30; % Square footage of the 1st floor, m^2
A2 = 30; % Square footage of 2nd floor, m^2
rho_air = 1.2; % kg/m^3
Cair = 1005; % 1005 J/(kg*K)
V1 = 90; % volume of lower floor
V2 = 38.97; % volume of upper floor
h = 1; % stepsize
t = 0:h:2*86400;
N = length(t);
% Initialize Arrays
T1 = zeros(1,N);
T2 = zeros(1,N);

% IC
T1(1)  = 5; % C
T2(1) = 7; % C
T_out = -10*sin((2*pi.*t)/86400); % C

for i = 1:N-1
    Q1 = k1*(T1(i) - T_out(i));
    Q2 = k2*(T1(i) - T2(i) + 5);
    Q3 = k3*(T2(i) - T_out(i));
    dT1 = (Ain*Qin)/(Cair*rho_air*V1) - (A1*Q1)/(Cair*rho_air*V1) - (A1*Q2)/(Cair*rho_air*V1);
    dT2 = (Ain*Qin)/(Cair*rho_air*V1) + (A1*Q2)/(Cair*rho_air*V1) + (A2*Q3)/(Cair*rho_air*V2);
    T1(i+1) = T1(i) + h*dT1;
    T2(i+1) = T2(i) + h*dT2;
end

figure(1)
plot(t,T1,'r-',t,T2,'g-')
legend('T1','T2','Location','best')