How to solve a thermal model? NEED HELP!!

Hi everyone,

I'm working on a transient thermal analysis using a thermal network model in MATLAB, and I'm encountering an issue I can't quite resolve. I want to solve this energy balance:

System Description:

The model consists of 16 nodes, where node 16 represents ambient air. Internal heat is generated at various nodes (1–15), and heat can only exit the system through node 15, which is thermally connected to the ambient (heat sink). The goal is to simulate temperature evolution over time in this network.

Model Setup:

S – a [16×16] thermal conductance matrix, including node 16.

S_eff – the [15×15] reduced matrix, excluding the ambient node.

Q_gen – a [15×1] vector of internal heat generation values.

MC – a [15×1] vector of thermal capacities (mass × specific heat) for nodes 1–15.

G_air – a [15×1] vector, where G_air(i) = 1 / R_i-amb, modeling conductance to ambient.

T_air – a fixed ambient temperature (293.15 K).

Transient Solver:

I'm using a forward Euler integration scheme. My understanding is that the term G_air .* (T_air) captures the heat leaving each node to the environment (with the only non-zero value corresponding to node 15), while -S_eff * T_nodes represents heat transfer between internal nodes. The energy balance at each time step is:

MC .* (dT/dt) = Q_gen + G_air .* (T_air) - S_eff * T

MATLAB Code:

T0 = 298.15;   % Initial temperature in K (25 °C)
T_air = 293.15;
dt = 0.01; tol = 0.1; max_steps = 1e6;
T = T0 * ones(15,1);  % Nodes 1–15
T_hist = T;
Q_eff = Q_gen + G_air * T_air;
time_vector = 0;
for k = 1:max_steps
    dTdt = (Q_eff - S_eff * T) ./ MC;
    dE = abs(dTdt .* MC);
    if all(dE < tol)
        fprintf('Converged at step %d (t = %.2f s)\n', k, k * dt);
        break;
    end
    T = T + dt * dTdt;
    T_hist(:, end+1) = T;
    time_vector(end+1) = k * dt;
end

Problem:

The steady-state solution, calculated as:

T_nodes = S_eff \ (Q_gen + G_air * T_air);

returns a physically reasonable result (e.g., ~174°C). However, the transient solution diverges downward, reaching unrealistic values (e.g., −157°C), as if the system were cooling despite internal heating.

I’ve attached my simulation results and input matrices in case they’re helpful. If anyone can spot what’s going wrong or has suggestions on better practices for transient simulation in MATLAB, I’d really appreciate it!

Thanks so much in advance!

[Torsten]

dTdt = (Q_eff - S_eff * T - G_air .* T) ./ MC;

is not in accordance with your energy balance

MC .* (dT/dt) = Q_gen + G_air .* (T_air) - S_eff * T

[Giselle]

That was a typo.. fixed it.

[Torsten]

I guess your ODEs for T_i are stiff, and you use an explicit integration method (Euler forward). Did you try a smaller time step ?

[Giselle]

I have tried reducing the time step, as well using ode23t. both results in reaching steady state at a very low temperature, like the curve is mirrored over the horizon.

This is result of the current explicit method..

✅ Converged at step 324444 (t = 3244.44 s)

--- Final Transient Node Temperatures ---

Node 1: -194.00 °C

Node 2: -184.41 °C

Node 3: -177.69 °C

Node 4: -170.84 °C

Node 5: -161.72 °C

Node 6: -165.92 °C

Node 7: -157.23 °C

Node 8: -151.76 °C

Node 9: -151.76 °C

Node 10: -154.44 °C

Node 11: -154.44 °C

Node 12: -151.65 °C

Node 13: -151.65 °C

Node 14: -157.22 °C

Node 15: -147.38 °C

Answer 1

[Torsten]

Params = load('Params.mat');
S_matrix = load('S_matrix.mat');
MC = Params.MC;
Q = Params.Q;
S = S_matrix.S;
T0 = 298.15;
T_air = 293.15;  
y0 = [T0*ones(17,1);T_air];
y0(1) = (Q(1) - S(1,3:18)*y0(3:18))/S(1,1);
y0(2) = (Q(2) - S(2,3:18)*y0(3:18))/S(2,2);
M = eye(18);
M(1,1) = 0;
M(2,2) = 0;
options = odeset('Mass',M);
tspan = [0 20000];
[T,Y] = ode15s(@(t,y)fun(t,y,MC,Q,S),tspan,y0,options);
plot(T,Y-273.15)
grid on

Y(end,1:10)-273.15

ans = 1×10

348.6030 354.9855 334.1085 332.4272 331.4959 330.3602 326.7618 327.1214 324.9385 320.9913

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Y(end,11:18)-273.15

ans = 1×8

320.9863 319.3709 319.3712 317.5164 317.5152 325.9185 311.2224 20.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function dTdt = fun(t,T,MC,Q,S)
  n = numel(T);
  dTdt = zeros(n,1);
  dTdt(1) = Q(1) - S(1,:)*T;
  dTdt(2) = Q(2) - S(2,:)*T;
  dTdt(3:n-1) = (Q(3:n-1) - S(3:n-1,:)*T)./MC(3:n-1);
  dTdt(n) = 0; 
end

[Giselle]

Amazing.. Thank you very much.. So, my approach was wrong after all..

[Torsten]

The last row in the S-matrix seems to indicate that the air is heated up by node 17, but I ignored this and set the air temperature to be constant over time. If this is not wanted, you will have to set an MC value for air.

[Giselle]

No, that is correct.. the room temperature is kept constant by the ventilation.

I saw that too, and was confused why it is + rather than - (air as heat sink). Should I somehow force the last line in Q to be - ?

[Torsten]

The last line of the S-matrix would give

MC(air)*dT_air/dt = -(-11.3153*T17 + 11.3153*T_air) = 11.3153*(T17-T_air)

and this is in principle correct because the air cools node 17 and acts as a heat sink (but would at the same time heat up by this setting).

What do you mean by "force the last line in Q to be -" ?

[Giselle]

I meant to put + in my defined function dTdt = (Q - S(1:17, 1:17) * T - S(1:17, 18) * T_air) ./ MC, so it becomes dTdt = (Q - S(1:17, 1:17) * T + S(1:17, 18) * T_air) ./ MC. But since you correctly mentioned it is correct in principal, I let the function to be as it should dTdt = (Q - S(1:17, 1:17) * T - S(1:17, 18) * T_air) ./ MC.

Also, if you don't mind, I need to ask another question..

The housing also transfers heat through radiation, and since radiation depends on surface temperature, I need to update the last thermal resistance at each time step. I implemented it as follows:

load('Params.mat', 'Q', 'MC');
load('S_matrix.mat', 'S');
Properties
SA_H=1.1412;
l_ch_h=[0.1140 0.2240];
U_air=2;
use_forced_convection = true;  % <— Toggle this to false for natural convection
T0 = 298.15 * ones(17, 1);
T_air = 293.15;
tspan = [0, 10000];  
S_base = S;
M = diag(MC);            

opts = odeset('Mass', M, 'RelTol', 1e-6, 'AbsTol', 1e-6);

[t_sol, T_sol] = ode15s(@(t, T) thermalODE_fullDAE(t, T, Q, MC, S_base, T_air, Stef_bolt,  ...
       eps, SA_H, k_liq, l_ch_h, Pr, nu, beta, use_forced_convection, U_air), tspan, T0, opts);

figure;
plot(t_sol, T_sol, 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Temperature (K)');
title('Transient Temperature of Nodes 1–17');
grid on;
hold on;

% Annotate steady-state (final) temperatures
final_time = t_sol(end);
final_temps = T_sol(end, :);

for n = 1:17
    text(final_time, final_temps(n), ...
        sprintf('%.1f K', final_temps(n)), ...
        'FontSize', 8, 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');
end

legend(arrayfun(@(n) sprintf('Node %d', n), 1:17, 'UniformOutput', false), ...
       'Location', 'eastoutside');

function dTdt = thermalODE_fullDAE(t, T, Q, MC, S_base, ...
    T_air, Stef_bolt, eps, SA_H, k_liq, l_ch_h, Pr, nu, beta, ...
    use_forced, U_air)

    T17 = T(17);  % Housing node

    % Get air properties (index 2)
    k_air   = k_liq(2);
    l_char  = l_ch_h(2);
    Pr_air  = Pr(2);
    nu_air  = nu(2);
    beta_air = beta(2);
    % --- Radiation ---
    R_rad = 1 / ((Stef_bolt * eps * SA_H) * (T17^2 + T_air^2) * (T17 + T_air));
    G_rad = 1 / R_rad;

    % --- Convection ---
    if use_forced
        Re = U_air * l_char / nu_air;
        Nu = 0.27 * Re^0.63;  % Tangential forced convection
    else
        Gr = 9.81 * beta_air * (T17 - T_air) * l_char^3 / (nu_air^2);
        Nu = 0.28 * (Gr * Pr_air)^0.3;  % Natural convection, vertical plate
    end

    R_conv = l_char / (SA_H * k_air * Nu);
    G_conv = 1 / R_conv;

    % --- Update S matrix ---
    G_total = G_rad + G_conv;
    S = S_base;  % Copy so original is not modified
    G_old = -S(17,18);
    S(17,18) = -G_total;
    S(18,17) = -G_total;
    S(17,17) = S(17,17) - G_old + G_total;
    S(18,18) = S(18,18) - G_old + G_total;

    % Compute dT/dt
    dTdt = (Q - S(1:17, 1:17) * T - S(1:17, 18) * T_air) ./ MC;
end

however, I get this error now "

Error using daeic12 (line 167)

Need a better guess y0 for consistent initial conditions.

Error in ode15s (line 304)

[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in ODE15S_DAE_V4 (line 21)

[t_sol, T_sol] = ode15s(@(t, T) thermalODE_fullDAE(t, T, Q, MC, S_base, T_air, Stef_bolt,..."

Can you by any chance decipher the error?!

[Torsten]

SA_H, l_ch_h and U_air cannot be resolved.

[Giselle]

those are system constants :

SA_H=1.1412;
l_ch_h=[0.1140 0.2240];
U_air=2;

[Torsten]

You already set M = diag(MC).

So

dTdt = (Q - S(1:17, 1:17) * T - S(1:17, 18) * T_air) ./ MC;

has to be changed to

dTdt = Q - S(1:17, 1:17) * T - S(1:17, 18) * T_air;

Otherwise you would divide by MC.^2 (and additionally by 0 for nodes 1 and 2).

[Giselle]

Thank you again for your time and big help :)

Answer 2

[Torsten]

I get these curves with your equations and your data.

T0 = 298.15;   
T_air = 293.15;
G_air = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 22.6305575112];
MC = [137.8208259751 591.4641353992 182.006441341 2436.9295445096 454.7099992049 640.8641472606 1679.8407170664 245.2936372241 245.2936372241 237.4649788763 237.4649788763 237.4649788763 237.4649788763 2235.552 13540.63392];
Q_gen = [359.1893251004 346.2558099057 410.9467641385 373.4019314744 48.447 0 22.0506266673 66.2445241332 66.0841062751 32.8397915689 32.8534513026 4.7328994143 4.6852295455 1527.6932754538 0];
S_eff = [66.7227944533383	0	0	0	-21.125528529005	0	0	0	0	0	0	0	0	-22.1510837715006	0
0	85.6332899505597	0	0	0	-6.26902680203049	0	0	0	0	0	0	0	-56.7623197632202	0
0	0	103.735901694769	0	0	-63.486611502084	0	0	0	0	0	0	0	-23.148051397717	0
0	0	0	104.875162688071	0	0	-14.1390959198035	0	0	0	0	0	0	-75.197227789859	0
-21.125528529005	0	0	0	136.868718804119	0	0	-10.3490606606084	-10.3490606606084	0	0	0	0	-95.0450689538968	0
0	-6.26902680203049	-63.486611502084	0	0	159.296129885689	0	0	0	-15.2814860879717	-15.2814860879717	0	0	-58.9775194056312	0
0	0	0	-14.1390959198035	0	0	175.473307903045	0	0	0	0	-15.2814860879717	-15.2814860879717	-130.771239807298	0
0	0	0	0	-10.3490606606084	0	0	32.2823856766301	0	0	0	0	0	-6.008559285848	-15.9247657301737
0	0	0	0	-10.3490606606084	0	0	0	32.2823856766301	0	0	0	0	-6.008559285848	-15.9247657301737
0	0	0	0	0	-15.2814860879717	0	0	0	41.0438730036644	0	0	0	-3.99063940270331	-21.7717475129894
0	0	0	0	0	-15.2814860879717	0	0	0	0	41.0438730036644	0	0	-3.99063940270331	-21.7717475129894
0	0	0	0	0	0	-15.2814860879717	0	0	0	0	39.3004065856245	0	-2.24717298466343	-21.7717475129894
0	0	0	0	0	0	-15.2814860879717	0	0	0	0	0	39.3004065856245	-2.24717298466343	-21.7717475129894
-22.1510837715006	-56.7623197632202	-23.148051397717	-75.197227789859	-95.0450689538968	-58.9775194056312	-130.771239807298	-6.008559285848	-6.008559285848	-3.99063940270331	-3.99063940270331	-2.24717298466343	-2.24717298466343	807.096354173603	-320.551099938051
0	0	0	0	0	0	0	-15.9247657301737	-15.9247657301737	-21.7717475129894	-21.7717475129894	-21.7717475129894	-21.7717475129894	-320.551099938051	462.118178961545];
y0 = T0*ones(15,1);
tspan = [0 3600];
[T,Y] = ode15s(@(t,y)fun(t,y,T_air,G_air,MC,Q_gen,S_eff),tspan,y0);
plot(T,Y)
grid on

Y(end,1:10)

ans = 1×10

79.1482 88.7339 95.4541 102.3082 111.4276 107.2298 115.9228 121.3899 121.3849 118.7081

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Y(end,11:15)

ans = 1×5

118.7085 121.4966 121.4954 115.9237 125.7667

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function dTdt = fun(t,T,T_air,G_air,MC,Q_gen,S_eff)
  n = size(T,1);
  dTdt = zeros(n,1);
  for i = 1:n
   dTdt(i) = (Q_gen(i) + G_air(i)*T_air - S_eff(i,:)*T)/MC(i);
  end
end

[Giselle]

Yes, you're solving for [K], and I'm converting to [C], but the trend is similar. My issue is that I'm heating up the system, but I’m seeing a downward trend—I should be seeing an increase in temperature. And I don't understand why..

[Torsten]

And I don't understand why..

But you don't heat it up sufficiently because obviously, Q_gen(i) < S_eff(i,:)*T right from the beginning. If you multiply Q_gen by a factor of 10, temperatures will rise.

[Giselle]

Aah, you mean my solution and solver are correct, but my S matrix is arranged in a way that causes the system to dump more heat than it receives..

So, I should double check my S matrix, rather than checking my methodology..

[Torsten]

Yes. But it surprises me that you didn't get the same final temperatures from both the steady-state and the time-dependent computation.

With the data given, I get:

T0 = 298.15;   
T_air = 293.15;
G_air = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 22.6305575112];
MC = [137.8208259751 591.4641353992 182.006441341 2436.9295445096 454.7099992049 640.8641472606 1679.8407170664 245.2936372241 245.2936372241 237.4649788763 237.4649788763 237.4649788763 237.4649788763 2235.552 13540.63392];
Q_gen = [359.1893251004 346.2558099057 410.9467641385 373.4019314744 48.447 0 22.0506266673 66.2445241332 66.0841062751 32.8397915689 32.8534513026 4.7328994143 4.6852295455 1527.6932754538 0];
S_eff = [66.7227944533383	0	0	0	-21.125528529005	0	0	0	0	0	0	0	0	-22.1510837715006	0
0	85.6332899505597	0	0	0	-6.26902680203049	0	0	0	0	0	0	0	-56.7623197632202	0
0	0	103.735901694769	0	0	-63.486611502084	0	0	0	0	0	0	0	-23.148051397717	0
0	0	0	104.875162688071	0	0	-14.1390959198035	0	0	0	0	0	0	-75.197227789859	0
-21.125528529005	0	0	0	136.868718804119	0	0	-10.3490606606084	-10.3490606606084	0	0	0	0	-95.0450689538968	0
0	-6.26902680203049	-63.486611502084	0	0	159.296129885689	0	0	0	-15.2814860879717	-15.2814860879717	0	0	-58.9775194056312	0
0	0	0	-14.1390959198035	0	0	175.473307903045	0	0	0	0	-15.2814860879717	-15.2814860879717	-130.771239807298	0
0	0	0	0	-10.3490606606084	0	0	32.2823856766301	0	0	0	0	0	-6.008559285848	-15.9247657301737
0	0	0	0	-10.3490606606084	0	0	0	32.2823856766301	0	0	0	0	-6.008559285848	-15.9247657301737
0	0	0	0	0	-15.2814860879717	0	0	0	41.0438730036644	0	0	0	-3.99063940270331	-21.7717475129894
0	0	0	0	0	-15.2814860879717	0	0	0	0	41.0438730036644	0	0	-3.99063940270331	-21.7717475129894
0	0	0	0	0	0	-15.2814860879717	0	0	0	0	39.3004065856245	0	-2.24717298466343	-21.7717475129894
0	0	0	0	0	0	-15.2814860879717	0	0	0	0	0	39.3004065856245	-2.24717298466343	-21.7717475129894
-22.1510837715006	-56.7623197632202	-23.148051397717	-75.197227789859	-95.0450689538968	-58.9775194056312	-130.771239807298	-6.008559285848	-6.008559285848	-3.99063940270331	-3.99063940270331	-2.24717298466343	-2.24717298466343	807.096354173603	-320.551099938051
0	0	0	0	0	0	0	-15.9247657301737	-15.9247657301737	-21.7717475129894	-21.7717475129894	-21.7717475129894	-21.7717475129894	-320.551099938051	462.118178961545];
y0 = T0*ones(15,1);
T = fsolve(@(y)fun(y,T_air,G_air,MC,Q_gen,S_eff),y0);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

T(1:10).'

ans = 1×10

79.1481 88.7338 95.4539 102.3080 111.4275 107.2297 115.9226 121.3897 121.3848 118.7080

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

T(11:15).'

ans = 1×5

118.7083 121.4965 121.4952 115.9235 125.7665

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function res = fun(T,T_air,G_air,MC,Q_gen,S_eff)
  n = size(T,1);
  res = zeros(n,1);
  for i = 1:n
   res(i) = Q_gen(i) + G_air(i)*T_air - S_eff(i,:)*T;
  end
end

[Giselle]

Thank you for your interest and for bringing this up.

Let me explain the issue I’m facing in more detail.

I’m working on a thermal model of a mechanical system, where the losses act as internal heat sources. Two of these sources come from friction between solid bodies. It’s common practice to model this type of interaction using a massless node to represent the friction interface, which is thermally connected to two bulk nodes (representing the bodies in contact) through finite thermal resistances.

The issue arises when solving the transient system: these friction nodes are massless, so their thermal capacity (MC) is zero, which leads to division by zero in the differential equation formulation. For steady-state analysis, this isn’t a problem—since mass doesn’t affect steady-state temperature, I can solve the full 18-node system without modifying the MC vector. However, for the transient case, I needed to avoid the divide-by-zero issue.

To work around this, I redistributed the heat generated at the friction nodes to the two adjacent bulk nodes. Specifically, if Q was the heat at the massless node, I split it into Q1 and Q2 such that Q1 + Q2 = Q, and assigned these directly to the connected nodes. This allowed me to remove the massless nodes and reduce the system, ensuring that all remaining nodes had non-zero thermal capacity.

Now here’s where it gets confusing. When I run the steady-state solver on both the original (18-node) and the reduced system (15-node), I get different temperature results—even though the total heat input (Q) is the same in both cases and connectivity is preserved. I don’t understand why these results don’t match.

On top of that, I’m observing another issue that seems physically impossible. Heat transfer is driven by temperature difference—energy flows from hot to cold. My environment is at 293 K, and I’m adding heat to the system. But instead of heating up, the system cools down to around 100 K. That implies heat is somehow flowing from the cooler system to the warmer environment, which violates the second law of thermodynamics.

Even if the amount of heat added isn’t enough to significantly raise the temperature, the system should still not drop below the environment temperature. I’m defining positive Q as heat input, and using the basic energy balance Q = m·C·(T − T_env), which implies T should be greater than T_env. But for some unknown reason, I’m getting T much lower than T_env.

Does this make sense to you? If you are interested, I can share the original S, Q, and MC matrix as well..

[Torsten]

MC = 0 means that you have an algebraic equation for temperature in the friction nodes. Thus you have a differential-algebraic system of equations: 15 (16 ?) differential equations and 2 algebraic equations for the friction nodes. MATLAB ODE integrators solve differential-algebraic systems by defining a mass matrix and using the option "Mass" in the ode-options.

Thus you could try solving the original 18-node time-dependent system this way.

[Giselle]

Thank you for the clarification. I’m not very familiar with this approach. In my setup, I have two massless nodes (nodes 1 and 2), which are connected to nodes 3–4 and 5–6, respectively.

What I’ve been doing so far is redistributing the heat input (Q) from the friction nodes to the connected bulk nodes and then solving the reduced system using an ODE solver. Once I have the temperatures for nodes 3 to 6, I use the relation Q=ΔT/R to algebraically compute the temperatures of nodes 1 and 2.

However, I’m still facing the issue of a downward temperature trend—where the system temperature drops significantly below ambient despite positive heat input. That’s the part I can’t physically explain, and I’m not sure if it’s due to the redistribution approach or an issue with how I’m solving the system.

Could you please elaborate more on how to properly implement the mass matrix and set up the differential-algebraic system using MATLAB’s ODE solver with the "Mass" option?

[Giselle]

Here is the original Q, MC, and S, if you are interested.

node 18 is the air, which is not included in the Q and MC matrix.