Help with MATLAB ODE45 to solve Heat Transfer Model

Hi everyone,

Hi everyone,

I'm working on a project where I need to model heat transfer in a spherical product using MATLAB. My approach is based on a simple model that interprets the heat transfer in terms of heat capacities and resistances as described in this article.

The evolution of the product's average temperature is formulated as:

I've drawn inspiration from this paper where a similar problem was solved using the ode45 function, and I'm trying to replicate it with my experimental data.

However, when I run my code, I encounter the following error:

______________________________________________________________________

Error

Error using odearguments (line 95)

FUNCTION must return a column vector.

Error in ode45 (line 104)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in untitled11 (line 50)

[t, y] = ode45(@Function, time, y0);

^^^^^^^^^^^^^^^^^^^^^^^^^^

______________________________________________________________________

Could someone please help me modify the ODE function so it works correctly with ode45 ? I believe the issue might be related to indexing or how I'm handling the input data (). I would appreciate any advice or guidance.

Thanks in advance!

Here's my script:

clear; close all; clc;

global Tair_meas time Rint Rext rho cp lambda R V A hconv t i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  constants
rho = 898;          % Density (kg/m^3)
cp = 3829;          % Specific heat capacity (J/(kg*K))
lambda = 0.463;     % Thermal conductivity (W/(m*K))
R = 0.05;           % Radius (m)
A = 4 * pi * (R^2); % Surface area (m^2)
V = (4/3) * pi * R^3; % Volume (m^3)
hconv = 15;         % Convective heat transfer coefficient (W/(m^2*K))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% thermal resistances
Rint = 1 / (A * (4 * lambda / R)); % Internal resistance
Rext = 1 / (hconv * A);            % External resistance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time and measured air temperature data
time = [0.01474, 0.92313, 1.83754, 3.45581, 6.33138, 9.13477, 11.90206, ...
        14.70545, 17.50884, 20.34832, 23.15171, 26.08744, 28.89083];
Tair_meas = [5.7557, 5.78514, 4.48993, 3.05735, 4.74014, 4.82355, 4.6273, ...
             4.4654, 4.40653, 4.43596, 4.65674, 4.4654, 4.4654];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time and measured product average temperature data
time2 = [0.01474 0.85696 1.77137 3.52198 6.33138 9.0325 11.93816 14.66936...
     17.47876 20.28215 23.08553 26.08744 28.89083];
Tavg_meas = [16.85819 12.01097 9.42056 6.50143 5.81458 6.05988...
        5.95195 5.73117 5.5104 5.18169 5.12772 5.12772 5.06885];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulated product average temperature data (Article)
    time3 = [0 0.27944 0.4539 0.62234 0.79078 1.09158 1.36229 ...
        1.73527 2.13834 2.71586 3.22119 3.66035 4.23185 4.83945 5.3508...
     5.92231 6.46373 7.54659 8.4249 9.47166 10.62069 11.29446...
     12.54576 13.15336 13.9655 14.80772 16.16129 17.06968...
     18.11644 19.53618 21.12436 21.76806 22.71255 23.59086...
     24.30073 24.94443 25.55203 26.66496 27.64555 28.31932...
     28.86075 29.63679 29.97368];
    Tavg_ref = [16.77479 15.67092 14.65536 13.63489 12.45252...
     10.85313 9.91607 8.42952 7.30112 6.05988 5.12772 4.6273 3.96988...
     3.66571 3.47437 3.63627 3.94045 4.29859 4.5488 4.57333 4.5488...
     4.74014 4.76958 4.68617 4.68617 4.60277 4.51937 4.48993 4.4654...
     4.35256 4.32803 4.43596 4.26916 4.32803 4.48993 4.4654 4.51937...
     4.48993 4.51937 4.57333 4.57333 4.48993 4.48993];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial condition for the ODE
y0 = 16;
i = [1:1:length(time)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the ODE using ode45
[t, y] = ode45(@Function, time, y0);

Error using odearguments (line 95)
FUNCTION must return a column vector.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the results
figure;
hold on;
plot(time, Tair_meas, '-o', 'LineWidth', 2, 'DisplayName', 'Tair measured');
plot(t, y, '-', 'LineWidth', 2, 'DisplayName', 'Tavg model');
plot(time2, Tavg_meas, '-x', 'LineWidth', 2, 'DisplayName', 'Tair measured');
plot(time3, Tavg_ref, '-x','LineWidth', 2)
xlabel('Time (hr)');
ylabel('Temperature (°C)');
legend('Location', 'best');
set(gca, 'FontName', 'LM Roman 10', 'TickDir', 'in', 'Box', 'on', ...
    'FontWeight', 'bold', 'FontSize', 12);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function definition for the ODE
function dydt = Function(t, y)
    global  rho cp V Rint Rext Tair_meas i
    dydt = (Tair_meas(i) - y) / ((Rint + Rext) * (rho * cp * V));
end

Answer 1

[Cris LaPierre]

The odefxn is returning a row vector instead of a column vector. However, looking at your code, it should be returning a scalar. The issue is that Tair_meas(i) is not doing what you expect.

ode45 is a variable step-size solver, meaning it calls odefxn for variable times t. However, you have defined values for Tair_meas at specific time values time. Since you can't know ahead of time what values of t ode45 will use, I think the best solution here is to use interp1 to estimate the value of Tair_meas at t.

Using global variables can have unintended consequences, so I removed the global variables and instead passed extra parameters to the ode function following this example. I also renamed the function.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  constants
rho = 898;          % Density (kg/m^3)
cp = 3829;          % Specific heat capacity (J/(kg*K))
lambda = 0.463;     % Thermal conductivity (W/(m*K))
R = 0.05;           % Radius (m)
A = 4 * pi * (R^2); % Surface area (m^2)
V = (4/3) * pi * R^3; % Volume (m^3)
hconv = 15;         % Convective heat transfer coefficient (W/(m^2*K))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% thermal resistances
Rint = 1 / (A * (4 * lambda / R)); % Internal resistance
Rext = 1 / (hconv * A);            % External resistance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time and measured air temperature data
time = [0.01474, 0.92313, 1.83754, 3.45581, 6.33138, 9.13477, 11.90206, ...
        14.70545, 17.50884, 20.34832, 23.15171, 26.08744, 28.89083];
Tair_meas = [5.7557, 5.78514, 4.48993, 3.05735, 4.74014, 4.82355, 4.6273, ...
             4.4654, 4.40653, 4.43596, 4.65674, 4.4654, 4.4654];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time and measured product average temperature data
time2 = [0.01474 0.85696 1.77137 3.52198 6.33138 9.0325 11.93816 14.66936...
     17.47876 20.28215 23.08553 26.08744 28.89083];
Tavg_meas = [16.85819 12.01097 9.42056 6.50143 5.81458 6.05988...
        5.95195 5.73117 5.5104 5.18169 5.12772 5.12772 5.06885];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulated product average temperature data (Article)
    time3 = [0 0.27944 0.4539 0.62234 0.79078 1.09158 1.36229 ...
        1.73527 2.13834 2.71586 3.22119 3.66035 4.23185 4.83945 5.3508...
     5.92231 6.46373 7.54659 8.4249 9.47166 10.62069 11.29446...
     12.54576 13.15336 13.9655 14.80772 16.16129 17.06968...
     18.11644 19.53618 21.12436 21.76806 22.71255 23.59086...
     24.30073 24.94443 25.55203 26.66496 27.64555 28.31932...
     28.86075 29.63679 29.97368];
    Tavg_ref = [16.77479 15.67092 14.65536 13.63489 12.45252...
     10.85313 9.91607 8.42952 7.30112 6.05988 5.12772 4.6273 3.96988...
     3.66571 3.47437 3.63627 3.94045 4.29859 4.5488 4.57333 4.5488...
     4.74014 4.76958 4.68617 4.68617 4.60277 4.51937 4.48993 4.4654...
     4.35256 4.32803 4.43596 4.26916 4.32803 4.48993 4.4654 4.51937...
     4.48993 4.51937 4.57333 4.57333 4.48993 4.48993];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial condition for the ODE
y0 = 16;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the ODE using ode45
[t, y] = ode45(@(t,y) odefxn(t,y,time*3600,Tair_meas,Rint,Rext,rho,cp,V), time*3600, y0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the results
figure;
hold on;
plot(time, Tair_meas, '-o', 'LineWidth', 2, 'DisplayName', 'Tair measured');
plot(t/3600, y, '-', 'LineWidth', 2, 'DisplayName', 'Tavg model');
plot(time2, Tavg_meas, '-x', 'LineWidth', 2, 'DisplayName', 'Tair measured');
plot(time3, Tavg_ref, '-x','LineWidth', 2)
xlabel('Time (hr)');
ylabel('Temperature (°C)');
legend('Location', 'best');
set(gca, 'FontName', 'LM Roman 10', 'TickDir', 'in', 'Box', 'on', ...
    'FontWeight', 'bold', 'FontSize', 12);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function definition for the ODE
function dydt = odefxn(t,y,time,Tair_meas,Rint,Rext,rho,cp,V)
    Tam = interp1(time, Tair_meas,t);
    dydt = (Tam - y) / ((Rint + Rext) * (rho * cp * V));
end

[Torsten]

Multiply the time vectors by a factor of 3600. Otherwise you work with inconsistent units.

[Zakaria]

Thank you so much for your detailed answer and help! I really appreciate the time and effort you put into explaining the issues and providing suggestions.

Yes, I have tried using the interp1 function to estimate air temperature (Tair_meas) at the variable times t that ode45 uses. However, the results were the same. The product average temperature (orange line) estimated by my model is still not close to the product average temperature of the reference model (cross violet line).

[Cris LaPierre]

I missed that time was in hours. Here, ode45 needs time to be in seconds.

Per @Torsten's recommendation, I updated the code to pass time into ode45 in seconds by multiplying by 3600. That also means the output times will be in seconds, so when plotting, I divided t by 3600.

...
    [t, y] = ode45(@(t,y) odefxn(t,y,time*3600,Tair_meas,Rint,Rext,rho,cp,V), time*3600, y0)
...
plot(t/3600, y, '-', 'LineWidth', 2, 'DisplayName', 'Tavg model');

[Walter Roberson]

ode45() does not care what the time units are. Time in ode45 should be understood to be in "timeunits" rather than seconds.

[Walter Roberson]

   Tam = interp1(time, Tair_meas,t);

The resutls of interp1() with default parameters have derivatives that are discontinuous at every time entry. You need

      Tam = interp1(time, Tair_meas, t, 'spline');

to have the needed derivative continuity.

[Torsten]

ode45() does not care what the time units are. Time in ode45 should be understood to be in "timeunits" rather than seconds.

cp is in unit J/(kg*K) = W*s/(kg*K). Thus t must be in unit s, not h for the integration.

But you are right: ode45 doesn't care about time units. It's in the responsibility of the user to use consistent physical units. I recommend SI units - then such problems cannot arise.

[Zakaria]

Thank you all for your comments. I appreciate your help!

I realize now that I didn’t account for the time unit properly when plotting. The scenario I shared was just a test to understand how ode45 works. My goal is to perform a co-simulation between EnergyPlus and MATLAB.

Here’s my approach:

1. EnergyPlus calculates the internal temperature (Tamb) at a specific time step t (e.g., a single value at time t) and sends it to MATLAB.
2. MATLAB uses this Tamb to solve an ODE problem for that same time t, calculating the product temperature (Tavg).
3. Using Tavg, MATLAB then calculates the heat load (Qprod) and sends it back to EnergyPlus.
4. EnergyPlus updates to the next time step (t + Dt) and calculates the new Tamb, which is again sent to MATLAB to repeat the process.

The problem lies with ode45 in MATLAB, which requires the time to be defined as a vector (e.g., [0 1 ... t]). However, in my case, I need to solve the ODE problem for only a single time step (e.g., at t or t + Dt) rather than over an entire time vector.

My question is:

• How can I adapt ode45 (or any suitable method) to solve the ODE for just one time step at a time, aligning it with the single-value Tamb data from EnergyPlus?
• Are there alternative solvers or approaches better suited for this kind of stepwise co-simulation?

Thank you in advance for your help!

[Torsten]

ode45 is not suited to be used in such a co-simulation case.

I suggest you write your own fixed-step solver. The explicit or implicit Euler method should be a good starting point.