Help with MATLAB ODE45 to solve Heat Transfer Model

9 次查看（过去 30 天）
Zakaria
Zakaria 2024-12-9
评论： Torsten 2024-12-11
Ran in:
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

请先登录，再回答此问题。

回答（1 个）

Cris LaPierre
Cris LaPierre 2024-12-9
编辑：Cris LaPierre 2024-12-9
Ran in:
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
  8 个评论
显示 6更早的评论
Zakaria
Zakaria 2024-12-11
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
Torsten 2024-12-11
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by