Created a aquifer system with 3 reservoirs and cant get the switch to another reservoir to work

Hello,

First time using this so sorry for any mistakes. Im trying to create a aquifer system below the ground with 3 reservoirs. One reservoir is used for heat output, one is charged with heat and one doesnt do anything. After reservoir 1 is heated to 50 degrees Celcius im trying to make the system use reservoir 1 for its heat output instead of reservoir 2. Currently it does correctly start to heat the coldest reservoir, but does not start to use hottest reservoir for its heat output. This can be seen in figure 2. Could someone tell me what im doing wrong?

Kind regards,

This is my code:

% Parameters

T_inflow = 60+273.15; % °K

density_water = 1000; % kg/m^3

c_p = 4186; % J/kg/K

heat_exchanger_eff = 0.8;

rho_water = 998; % kg/m^3

% Parameters Aquifer in general

phi = 0.3; % Porosity

k = 1e-12; % Permeability (m^2)

k_soil = 0.001; % Thermal Conductivity of Soil (W/(m*K))

T_soil= 5 + 273.15; % °K

% Parameters Aquifer 1

T_initial_1 = 10+273.15; % °K

L1 = 100; % Length (m)

W1 = 100; % Width (m)

H1 = 10; % Height (m)

V_max_1 = L1 * W1 * H1; % Volume of the First Aquifer (m^3)

A_1 = 2*L1*W1+2*W1*H1+2*L1*H1; % Area of the First Aquifer (m^2)

m_1 = density_water * V_max_1; % Mass of First Aquifer (kg)

% Parameters Aquifer 2

T_initial_2 = 50+273.15; % °K

L2 = 500; % Length (m)

W2 = 50; % Width (m)

H2 = 10; % Height (m)

V_max_2 = L2 * W2 * H2; % Volume of the second Aquifer (m^3)

A_2 = 2*L2*W2+2*W2*H2+2*L2*H2; % Area of the second Aquifer (m^2)

m_2 = density_water * V_max_2; % Mass of second Aquifer (kg)

% Parameters Aquifer 3

T_initial_3 = 30+273.15; % °K

L3 = 300; % Length (m)

W3 = 100; % Width (m)

H3 = 8; % Height (m)

V_max_3 = L3 * W3 * H3; % Volume of the second Aquifer (m^3)

A_3 = 2*L3*W3+2*W3*H3+2*L3*H3; % Area of the second Aquifer (m^2)

m_3 = density_water * V_max_3; % Mass of second Aquifer (kg)

%% Daily Simulation

% Parameters for flow rate

inflow_start_time = 8; % h

inflow_end_time = 16; % h

outflow_start_time1 = 0; % h

outflow_end_time1 = 8; % h

outflow_start_time2 = 16; % h

outflow_end_time2 = 24; % h

peak_flow_rate = 200; % kg/s (this is the maximum, or peak, flow rate)

outflow_rate = 40; % kg/s (this is the maximum, or peak, flow rate)

mass_flow_in_rate = @(t) peak_flow_rate * max(0, sin(pi * ((t/3600) - inflow_start_time) / (inflow_end_time - inflow_start_time))).* ((t/3600) >= inflow_start_time & (t/3600) <= inflow_end_time);

mass_flow_out_rate = @(t) outflow_rate * (((t/3600) >= outflow_start_time1 & (t/3600) <= outflow_end_time1) | ((t/3600) >= outflow_start_time2 & (t/3600) <= outflow_end_time2));

% Convert times from hours to seconds

inflow_start_time = inflow_start_time * 60 * 60;

inflow_end_time = inflow_end_time * 60 * 60;

outflow_start_time1 = outflow_start_time1 * 60 * 60;

outflow_end_time1 = outflow_end_time1 * 60 * 60;

outflow_start_time2 = outflow_start_time2 * 60 * 60;

outflow_end_time2 = outflow_end_time2 * 60 * 60;

% Determine the initial reservoir based on the coldest temperature

[~, initial_reservoir] = min([T_initial_1, T_initial_2, T_initial_3]);

[~, initial_draining_reservoir] = max([T_initial_1, T_initial_2, T_initial_3]);

initial_conditions = [T_initial_1; T_initial_2; T_initial_3; initial_reservoir; initial_draining_reservoir];

% Time span

tspan = 0:60:24*3600; % for a 24-hour simulation, with outputs every minute

% Run ode45

[t, Y] = ode45(@(t, y) ates_ode_function(t, y, density_water, c_p, V_max_1, V_max_2, V_max_3, A_1, A_2, A_3, T_soil, T_inflow, k_soil, mass_flow_in_rate, mass_flow_out_rate, L1, W1, H1, L2, W2, H2, L3, W3, H3, m_1, m_2, m_3, heat_exchanger_eff), tspan, initial_conditions);

% Extract state variables

T1 = Y(:,1) - 273.15;

T2 = Y(:,2) - 273.15;

T3 = Y(:,3) - 273.15;

Q_in = arrayfun(mass_flow_in_rate, t);

Q_out = arrayfun(mass_flow_out_rate, t);

% Create the plot

figure(1);

subplot(2,1,1);

plot(t/3600, T1, 'r', 'LineWidth', 2); % Plot temperature in reservoir 1

hold on;

plot(t/3600, T2, 'b', 'LineWidth', 2); % Plot temperature in reservoir 2

plot(t/3600, T3, 'y', 'LineWidth', 2); % Plot temperature in reservoir 3

xlabel('Time (hours)');

ylabel('Temperature (°C)');

legend('T in Reservoir 1', 'T in Reservoir 2');

title('Temperature vs Time');

hold off;

subplot(2,1,2);

plot(t/3600, Q_in, 'r', 'LineWidth', 2);

hold on;

plot(t/3600, Q_out, 'b', 'LineWidth', 2);

xlabel('Time (hours)');

ylabel('Mass Flow Rate (kg/s)');

legend('Inflow Rate', 'Outflow Rate');

%% Weekly simulation

% Initialize storage for the entire week's data

t_week = [];

Y_week = [];

tspan_single_day = 0:60:24*3600; % for a 24-hour simulation, with outputs every minute

for day = 1:100 % Loop over each day of the week

% Run ode45 for a single day

[t, Y] = ode45(@(t, y) ates_ode_function(t, y, density_water, c_p, V_max_1, V_max_2, V_max_3, A_1, A_2, A_3, T_soil, T_inflow, k_soil, mass_flow_in_rate, mass_flow_out_rate, L1, W1, H1, L2, W2, H2, L3, W3, H3, m_1, m_2, m_3, heat_exchanger_eff), tspan, initial_conditions);

% Store this day's data

t_week = [t_week; t + (day-1)*24*3600]; % Adjust time for the week

Y_week = [Y_week; Y];

% Use the final state as the initial condition for the next day

initial_conditions = Y(end, :);

end

% Extract state variables from the weekly data

T1_week = Y_week(:,1) - 273.15;

T2_week = Y_week(:,2) - 273.15;

T3_week = Y_week(:,3) - 273.15;

% Calculate flow rates for plotting using the weekly data

Q_in_week = arrayfun(@(t) mass_flow_in_rate(mod(t, 24*3600)), t_week);

Q_out_week = arrayfun(@(t) mass_flow_out_rate(mod(t, 24*3600)), t_week);

% Create the weekly plot

figure(2); % Create a new figure for the weekly plot

subplot(2,1,1);

plot(t_week/(3600*24), T1_week, 'r', 'LineWidth', 2); % Plot temperature in reservoir 1

hold on;

plot(t_week/(3600*24), T2_week, 'b', 'LineWidth', 2); % Plot temperature in reservoir 2

plot(t_week/(3600*24), T3_week, 'y', 'LineWidth', 2); % Plot temperature in reservoir 2

xlabel('Time (days)');

ylabel('Temperature (°C)');

legend('T in Reservoir 1', 'T in Reservoir 2', 'T in Reservoir 3');

title('Weekly Temperature vs Time');

hold off;

subplot(2,1,2);

plot(t_week/(3600*24), Q_in_week, 'r', 'LineWidth', 2);

hold on;

plot(t_week/(3600*24), Q_out_week, 'b', 'LineWidth', 2);

xlabel('Time (days)');

ylabel('Mass Flow Rate (kg/s)');

legend('Weekly Inflow Rate', 'Weekly Outflow Rate');

%% Functions

function dydt = ates_ode_function(t, y, rho_water, c_p, V1, V2, V3, A1, A2, A3, T_soil, T_inflow, k_soil, mass_flow_in_rate, mass_flow_out_rate, L1, W1, H1, L2, W2, H2, L3, W3, H3, m_1, m_2, m_3, heat_exchanger_eff)

T1 = y(1);

T2 = y(2);

T3 = y(3);

current_reservoir = round(y(4)); % The current reservoir being heated

draining_reservoir = round(y(5));

% Set idx_in to current_reservoir initially

idx_in = current_reservoir;

idx_out = draining_reservoir;

if y(current_reservoir) >= (50+273.15)

[sorted_temps, sorted_idx] = sort([T1, T2, T3]);

% Modification: Update the reservoir being drained to the current reservoir if it reaches 50°C.

draining_reservoir = current_reservoir;

% Find the next coldest reservoir that hasn't reached 50°C yet.

idx_in = sorted_idx(1); % Default to the coldest

for i = 1:3

if sorted_temps(i) < (50+273.15)

idx_in = sorted_idx(i);

break;

end

end

% If the current reservoir is below 50°C, continue heating it.

if y(current_reservoir) < (50+273.15)

current_reservoir = idx_in;

end

end

idx_out = draining_reservoir; % Use the draining_reservoir to determine which one to drain

Q_in = zeros(1,3);

Q_out = zeros(1,3);

Q_in(current_reservoir) = mass_flow_in_rate(t);

Q_out(idx_out) = mass_flow_out_rate(t);

% Update flow rates based on the switching logic

switch idx_out

case 1

Q_out1 = mass_flow_out_rate(t);

Q_out2 = 0;

Q_out3 = 0;

case 2

Q_out1 = 0;

Q_out2 = mass_flow_out_rate(t);

Q_out3 = 0;

case 3

Q_out1 = 0;

Q_out2 = 0;

Q_out3 = mass_flow_out_rate(t);

end

switch idx_in

case 1

Q_in1 = mass_flow_in_rate(t);

Q_in2 = 0;

Q_in3 = 0;

case 2

Q_in1 = 0;

Q_in2 = mass_flow_in_rate(t);

Q_in3 = 0;

case 3

Q_in1 = 0;

Q_in2 = 0;

Q_in3 = mass_flow_in_rate(t);

end

% Heat loss

heat_loss1_top = k_soil * L1 * W1 * (T1 - T_soil);

heat_loss1_bottom = k_soil * L1 * W1 * (T1 - T_soil);

heat_loss1_sides = k_soil * 2 * (L1 + W1) * H1 * (T1 - T_soil);

heat_loss1 = heat_loss1_top + heat_loss1_bottom + heat_loss1_sides;

heat_loss2_top = k_soil * L2 * W2 * (T2 - T_soil);

heat_loss2_bottom = k_soil * L2 * W2 * (T2 - T_soil);

heat_loss2_sides = k_soil * 2 * (L2 + W2) * H2 * (T2 - T_soil);

heat_loss2 = heat_loss2_top + heat_loss2_bottom + heat_loss2_sides;

heat_loss3_top = k_soil * L3 * W3 * (T3 - T_soil);

heat_loss3_bottom = k_soil * L3 * W3 * (T3 - T_soil);

heat_loss3_sides = k_soil * 2 * (L3 + W3) * H3 * (T3 - T_soil);

heat_loss3 = heat_loss3_top + heat_loss3_bottom + heat_loss3_sides;

% ODEs for temperature

dT1dt = (Q_in1 * c_p * heat_exchanger_eff * (T_inflow - T1) + Q_out1 * c_p * heat_exchanger_eff * (T_soil - T1) - heat_loss1) / (m_1 * c_p);

dT2dt = (Q_in2 * c_p * heat_exchanger_eff * (T_inflow - T2) + Q_out2 * c_p * heat_exchanger_eff * (T_soil - T2) - heat_loss2) / (m_2 * c_p);

dT3dt = (Q_in3 * c_p * heat_exchanger_eff * (T_inflow - T3) + Q_out3 * c_p * heat_exchanger_eff * (T_soil - T3) - heat_loss3) / (m_3 * c_p);

dydt = [dT1dt; dT2dt; dT3dt; 0; 0];

end

%% Logic system

% Initially, pick the coldest reservoir to heat and the hottest one to extract heat.

% Continue heating that reservoir until it reaches 50°C and draining the hottest one.

% Once it reaches 50°C, switch to the next coldest reservoir that hasn't reached 50°C yet.

% Switch to the hottest reservoir for the output instead of the previous one.

% If all reservoirs are above 50°C, continue heating the least hot reservoir.

Answer 1

[SOUMNATH PAUL]

Hi,

I understand you are facing issues with the aquifer system not using the hottest reservoir for heat output, In my opinion we should try to adjust the logic within "ates_ode_function" that determines which reservoir is currently being drained. The code you provided already includes a mechanism for choosing the reservoir to be heated, but it seems like the logic for choosing the draining reservoir is not being updated correctly.

Here I have tried to revise the code:

% Parameters
T_inflow = 60+273.15; % °K
density_water = 1000; % kg/m^3
c_p = 4186; % J/kg/K
heat_exchanger_eff = 0.8;
rho_water = 998; % kg/m^3
% Parameters Aquifer in general 
phi = 0.3;               % Porosity
k = 1e-12;               % Permeability (m^2)
k_soil = 0.001;            % Thermal Conductivity of Soil (W/(m*K))
T_soil= 5 + 273.15;      % °K
% Parameters Aquifer 1 
T_initial_1 = 10+273.15;  % °K
L1 = 100;                 % Length (m)
W1 = 100;                 % Width (m)
H1 = 10;                  % Height (m)
V_max_1 = L1 * W1 * H1;          % Volume of the First Aquifer (m^3)
A_1 = 2*L1*W1+2*W1*H1+2*L1*H1;   % Area of the First Aquifer (m^2)
m_1 = density_water * V_max_1;   % Mass of First Aquifer (kg)
% Parameters Aquifer 2 
T_initial_2 = 50+273.15; % °K
L2 = 500;                 % Length (m)
W2 = 50;                  % Width (m)
H2 = 10;                  % Height (m)
V_max_2 = L2 * W2 * H2;          % Volume of the second Aquifer (m^3)
A_2 = 2*L2*W2+2*W2*H2+2*L2*H2;   % Area of the second Aquifer (m^2)
m_2 = density_water * V_max_2;   % Mass of second Aquifer (kg)
% Parameters Aquifer 3 
T_initial_3 = 30+273.15; % °K
L3 = 300;                 % Length (m)
W3 = 100;                  % Width (m)
H3 = 8;                  % Height (m)
V_max_3 = L3 * W3 * H3;          % Volume of the second Aquifer (m^3)
A_3 = 2*L3*W3+2*W3*H3+2*L3*H3;   % Area of the second Aquifer (m^2)
m_3 = density_water * V_max_3;   % Mass of second Aquifer (kg)
%% Daily Simulation
% Parameters for flow rate
inflow_start_time = 8; % h
inflow_end_time = 16; % h
outflow_start_time1 = 0; % h
outflow_end_time1 = 8; % h
outflow_start_time2 = 16; % h
outflow_end_time2 = 24; % h
peak_flow_rate = 200; % kg/s (this is the maximum, or peak, flow rate)
outflow_rate = 40; % kg/s (this is the maximum, or peak, flow rate)
mass_flow_in_rate = @(t) peak_flow_rate * max(0, sin(pi * ((t/3600) - inflow_start_time) / (inflow_end_time - inflow_start_time))).* ((t/3600) >= inflow_start_time & (t/3600) <= inflow_end_time);
mass_flow_out_rate = @(t) outflow_rate * (((t/3600) >= outflow_start_time1 & (t/3600) <= outflow_end_time1) | ((t/3600) >= outflow_start_time2 & (t/3600) <= outflow_end_time2));
% Convert times from hours to seconds
inflow_start_time = inflow_start_time * 60 * 60;
inflow_end_time = inflow_end_time * 60 * 60;
outflow_start_time1 = outflow_start_time1 * 60 * 60;
outflow_end_time1 = outflow_end_time1 * 60 * 60;
outflow_start_time2 = outflow_start_time2 * 60 * 60;
outflow_end_time2 = outflow_end_time2 * 60 * 60;
% Determine the initial reservoir based on the coldest temperature
[~, initial_reservoir] = min([T_initial_1, T_initial_2, T_initial_3]);
[~, initial_draining_reservoir] = max([T_initial_1, T_initial_2, T_initial_3]);
initial_conditions = [T_initial_1; T_initial_2; T_initial_3; initial_reservoir; initial_draining_reservoir];
% Time span
tspan = 0:60:24*3600;  % for a 24-hour simulation, with outputs every minute
% Run ode45
[t, Y] = ode45(@(t, y) ates_ode_function(t, y, density_water, c_p, V_max_1, V_max_2, V_max_3, A_1, A_2, A_3, T_soil, T_inflow, k_soil, mass_flow_in_rate, mass_flow_out_rate, L1, W1, H1, L2, W2, H2, L3, W3, H3, m_1, m_2, m_3, heat_exchanger_eff), tspan, initial_conditions);
% Extract state variables
T1 = Y(:,1) - 273.15;
T2 = Y(:,2) - 273.15;
T3 = Y(:,3) - 273.15;
Q_in = arrayfun(mass_flow_in_rate, t);
Q_out = arrayfun(mass_flow_out_rate, t);
% Create the plot
figure(1);
subplot(2,1,1);
plot(t/3600, T1, 'r', 'LineWidth', 2);  % Plot temperature in reservoir 1
hold on;
plot(t/3600, T2, 'b', 'LineWidth', 2);  % Plot temperature in reservoir 2
plot(t/3600, T3, 'y', 'LineWidth', 2);  % Plot temperature in reservoir 3
xlabel('Time (hours)');
ylabel('Temperature (°C)');
legend('T in Reservoir 1', 'T in Reservoir 2');
title('Temperature vs Time');
hold off;
subplot(2,1,2);
plot(t/3600, Q_in, 'r', 'LineWidth', 2);
hold on;
plot(t/3600, Q_out, 'b', 'LineWidth', 2);
xlabel('Time (hours)');
ylabel('Mass Flow Rate (kg/s)');
legend('Inflow Rate', 'Outflow Rate');

%% Weekly simulation
% Initialize storage for the entire week's data
t_week = [];
Y_week = [];
tspan_single_day = 0:60:24*3600;  % for a 24-hour simulation, with outputs every minute
for day = 1:100  % Loop over each day of the week
    % Run ode45 for a single day
    [t, Y] = ode45(@(t, y) ates_ode_function(t, y, density_water, c_p, V_max_1, V_max_2, V_max_3, A_1, A_2, A_3, T_soil, T_inflow, k_soil, mass_flow_in_rate, mass_flow_out_rate, L1, W1, H1, L2, W2, H2, L3, W3, H3, m_1, m_2, m_3, heat_exchanger_eff), tspan, initial_conditions);
    
    % Store this day's data
    t_week = [t_week; t + (day-1)*24*3600];  % Adjust time for the week
    Y_week = [Y_week; Y];
    
    % Use the final state as the initial condition for the next day
    initial_conditions = Y(end, :);
end
% Extract state variables from the weekly data
T1_week = Y_week(:,1) - 273.15;
T2_week = Y_week(:,2) - 273.15;
T3_week = Y_week(:,3) - 273.15;
% Calculate flow rates for plotting using the weekly data
Q_in_week = arrayfun(@(t) mass_flow_in_rate(mod(t, 24*3600)), t_week);
Q_out_week = arrayfun(@(t) mass_flow_out_rate(mod(t, 24*3600)), t_week);
% Create the weekly plot
figure(2);  % Create a new figure for the weekly plot
subplot(2,1,1);
plot(t_week/(3600*24), T1_week, 'r', 'LineWidth', 2);  % Plot temperature in reservoir 1
hold on;
plot(t_week/(3600*24), T2_week, 'b', 'LineWidth', 2);  % Plot temperature in reservoir 2
plot(t_week/(3600*24), T3_week, 'y', 'LineWidth', 2);  % Plot temperature in reservoir 2
xlabel('Time (days)');
ylabel('Temperature (°C)');
legend('T in Reservoir 1', 'T in Reservoir 2', 'T in Reservoir 3');
title('Weekly Temperature vs Time');
hold off;
subplot(2,1,2);
plot(t_week/(3600*24), Q_in_week, 'r', 'LineWidth', 2);
hold on;
plot(t_week/(3600*24), Q_out_week, 'b', 'LineWidth', 2);
xlabel('Time (days)');
ylabel('Mass Flow Rate (kg/s)');
legend('Weekly Inflow Rate', 'Weekly Outflow Rate');

%% Functions
function dydt = ates_ode_function(t, y, rho_water, c_p, V1, V2, V3, A1, A2, A3, T_soil, T_inflow, k_soil, mass_flow_in_rate, mass_flow_out_rate, L1, W1, H1, L2, W2, H2, L3, W3, H3, m_1, m_2, m_3, heat_exchanger_eff)
    T1 = y(1);
    T2 = y(2);
    T3 = y(3);
    current_reservoir = round(y(4)); % The current reservoir being heated
    draining_reservoir = round(y(5)); % The current reservoir being drained

    % Check if any reservoir has reached the target temperature of 50°C
    % If so, update the draining reservoir to be the one that has reached the target temperature
    if T1 >= (50+273.15)
        draining_reservoir = 1;
    end
    if T2 >= (50+273.15)
        draining_reservoir = 2;
    end
    if T3 >= (50+273.15)
        draining_reservoir = 3;
    end

    % Determine the next reservoir to be heated based on the coldest temperature
    [~, idx_in] = min([T1, T2, T3]);

    % Set the current_reservoir to the coldest one if it hasn't reached 50°C yet
    if y(current_reservoir) < (50+273.15)
        current_reservoir = idx_in;
    end
    idx_out = draining_reservoir;  % Use the draining_reservoir to determine which one to drain
    Q_in = zeros(1,3);
    Q_out = zeros(1,3);
    Q_in(current_reservoir) = mass_flow_in_rate(t);
    Q_out(idx_out) = mass_flow_out_rate(t);
    % Update flow rates based on the switching logic
    switch idx_out
        case 1
            Q_out1 = mass_flow_out_rate(t);
            Q_out2 = 0;
            Q_out3 = 0;
        case 2
            Q_out1 = 0;
            Q_out2 = mass_flow_out_rate(t);
            Q_out3 = 0;
        case 3
            Q_out1 = 0;
            Q_out2 = 0;
            Q_out3 = mass_flow_out_rate(t);
    end
    switch idx_in
        case 1
            Q_in1 = mass_flow_in_rate(t);
            Q_in2 = 0;
            Q_in3 = 0;
        case 2
            Q_in1 = 0;
            Q_in2 = mass_flow_in_rate(t);
            Q_in3 = 0;
        case 3
            Q_in1 = 0;
            Q_in2 = 0;
            Q_in3 = mass_flow_in_rate(t);
    end
    % Heat loss
    heat_loss1_top = k_soil * L1 * W1 * (T1 - T_soil);
    heat_loss1_bottom = k_soil * L1 * W1 * (T1 - T_soil);
    heat_loss1_sides = k_soil * 2 * (L1 + W1) * H1 * (T1 - T_soil);
    heat_loss1 = heat_loss1_top + heat_loss1_bottom + heat_loss1_sides;
    heat_loss2_top = k_soil * L2 * W2 * (T2 - T_soil);
    heat_loss2_bottom = k_soil * L2 * W2 * (T2 - T_soil);
    heat_loss2_sides = k_soil * 2 * (L2 + W2) * H2 * (T2 - T_soil);
    heat_loss2 = heat_loss2_top + heat_loss2_bottom + heat_loss2_sides;
    
    heat_loss3_top = k_soil * L3 * W3 * (T3 - T_soil);
    heat_loss3_bottom = k_soil * L3 * W3 * (T3 - T_soil);
    heat_loss3_sides = k_soil * 2 * (L3 + W3) * H3 * (T3 - T_soil);
    heat_loss3 = heat_loss3_top + heat_loss3_bottom + heat_loss3_sides;
    % ODEs for temperature
    dT1dt = (Q_in1 * c_p * heat_exchanger_eff * (T_inflow - T1) + Q_out1 * c_p * heat_exchanger_eff * (T_soil - T1) - heat_loss1) / (m_1 * c_p);
    dT2dt = (Q_in2 * c_p * heat_exchanger_eff * (T_inflow - T2) + Q_out2 * c_p * heat_exchanger_eff * (T_soil - T2) - heat_loss2) / (m_2 * c_p);
    dT3dt = (Q_in3 * c_p * heat_exchanger_eff * (T_inflow - T3) + Q_out3 * c_p * heat_exchanger_eff * (T_soil - T3) - heat_loss3) / (m_3 * c_p);
    dydt = [dT1dt; dT2dt; dT3dt; 0; 0];
end

Hope it helps!

Regards,

Soumnath