Simple heat transfer problem

Hi everyone!

I'm having trouble solving a very simple problem. Basically two blocks at 2 different initial temperatures T1 and T2 in contact. We assume that they exchange heat only via their contact surface, the rest being insulated. We assume a very high thermal conductivity inside the material, so homogeneous temperature in the blocks. The heat equation therefore reads:

T1,t= -U*A*(T2-T1)

T2,t= U*A*(T2-T1), where T1,t and T2,t are the temporal derivatives of the temperatures T1 and T2, and T1 > T2.

I tried solving this with ode23, ode45 and also solving it with matrices, but all give very odd results (clearly wrong). We obviously some sort of exponential decay for T1 and an increase for T2, both going to an equilibrium temperature in between T1 and T2. In my case, both diverge.

Here is my code. I don't know what I'm doing wrong. Thank you !

syms T_c(t) T_hs(t)

%propreties of component

V_c=1;

rho_c=9000;

cp_c=450;

%properties of heat sink

V_hs=1;

rho_hs=1000;

cp_hs=4200;

%heat transfer coefficient

h=1000; % [W/m2K]

A=1; % [m2]

a=h*A/(rho_c*V_c*cp_c);

b=h*A/(rho_hs*V_hs*cp_hs);

A=[a -a; -b b];

T=[T_c; T_hs];

odes= diff(T) == A*T

C=T(0) == [60; 30];

[T_cSol(t), T_hsSol(t)] = dsolve(odes,C);

T_cSol(t) = simplify(T_cSol(t))

T_hsSol(t) = simplify(T_hsSol(t))

clf

fplot(T_cSol)

hold on

fplot(T_hsSol)

grid on

legend('T_cSol','T_hsSol','Location','best')

Answer 1

[Sulaymon Eshkabilov]

Hi,

There are a couple of points that you've overlooked. (1) is the sign of U; (2) is the simulation time. This is a slow process and thus, longer simulation is needed.

clearvars
syms T1(t) T2(t)
dT1 = diff(T1, t);
dT2 = diff(T2, t);
%propreties of component
V_c=1;
rho_c=9000;
cp_c=450;
%properties of heat sink
V_hs=1;
rho_hs=1000;
cp_hs=4200;
%heat transfer coefficient
h=1000;  % [W/m2K]
A=1; % [m2]
a=h*A/(rho_c*V_c*cp_c);
b=h*A/(rho_hs*V_hs*cp_hs);
A=[a -a; -b b];
U=1;               % Specify: U
T = [T1; T2];
dT = diff(T, t);
C=T(0) ==[60; 30];
EQN = dT==- U*A*T;  % (-) SIGN problem was here
C1=T1(0) ==60; C2=T2(0)==30; % [60; 30];
TS = dsolve(EQN, C1, C2);
SOLUTION = [TS.T1; TS.T2];   % To see analytical solution formulation
figure
fplot(TS.T1, [0, 10000])     % Longer simulation time
hold on
fplot(TS.T2, [0, 10000])     % Longer simulation time
grid on
legend('T_{cSol}','T_{hsSol}','Location','best')
xlabel('$t$', 'interpreter', 'latex')
ylabel('$T_{cSol}(t), T_{hsSol}(t)$', 'interpreter', 'latex')

Good luck.

Answer 2

[suresh]

clearvars
syms T1(t) T2(t)
dT1 = diff(T1, t);
dT2 = diff(T2, t);

% Properties of component
V_c = 1;
rho_c = 9000;
cp_c = 450;

% Properties of heat sink
V_hs = 1;
rho_hs = 1000;
cp_hs = 4200;

% Heat transfer coefficient
h = 1000; % [W/m2K]
A_surf = 1; % surface area [m2]

a = h * A_surf / (rho_c * V_c * cp_c);
b = h * A_surf / (rho_hs * V_hs * cp_hs);

% Define the system matrix
M = [a, -a; -b, b];

% Define variables vector
T = [T1; T2];
dT = diff(T, t);

% System of ODEs
eqns = dT == M * T;

% Initial conditions
conds = [T1(0) == 60, T2(0) == 30];

% Solve ODE system
S = dsolve(eqns, conds);

% Extract solutions
T1_sol = S.T1;
T2_sol = S.T2;

% Plotting solutions
figure
fplot(T1_sol, [0, 10000])
hold on
fplot(T2_sol, [0, 10000])
grid on
legend('T_1(t)', 'T_2(t)', 'Location', 'best')
xlabel('Time t')
ylabel('Temperature (°C)')
title('Heat Transfer Between Two Blocks')

Answer 3

[suresh]

Your MATLAB symbolic code attempts to solve the coupled ODE system representing heat transfer between two blocks exchanging heat through their contact surface only. The model and equations are:dT1dt=−UA(T2−T1)\frac{dT_1}{dt} = -U A (T_2 - T_1)dtdT1​​=−UA(T2​−T1​)dT2dt=UA(T2−T1)\frac{dT_2}{dt} = U A (T_2 - T_1)dtdT2​​=UA(T2​−T1​)

where UA=hA(1ρcVccpc+1ρhsVhscphs)U A = h A \left( \frac{1}{\rho_c V_c c_{p_c}} + \frac{1}{\rho_{hs} V_{hs} c_{p_{hs}}} \right)UA=hA(ρc​Vc​cpc​​1​+ρhs​Vhs​cphs​​1​) but you've defined coefficients separately and built matrix AAA accordingly.

Verification of your ODE system and solution approach:

1. Physical equations as given:

dT1dt=−UA(T2−T1)dT2dt=+UA(T2−T1)\frac{dT_1}{dt} = - U A (T_2 - T_1) \\ \frac{dT_2}{dt} = + U A (T_2 - T_1)dtdT1​​=−UA(T2​−T1​)dtdT2​​=+UA(T2​−T1​)

1. Your matrix A is:

A=[a−a−bb]A = \begin{bmatrix} a & -a \\ -b & b \end{bmatrix}A=[a−b​−ab​]

where

a=hAρcVccpc,b=hAρhsVhscphsa = \frac{hA}{\rho_c V_c c_{p_c}}, \quad b = \frac{hA}{\rho_{hs} V_{hs} c_{p_{hs}}}a=ρc​Vc​cpc​​hA​,b=ρhs​Vhs​cphs​​hA​

and

ddt[T1T2]=−U⋅A[T1T2]\frac{d}{dt} \begin{bmatrix} T_1 \\ T_2 \end{bmatrix} = - U \cdot A \begin{bmatrix} T_1 \\ T_2 \end{bmatrix}dtd​[T1​T2​​]=−U⋅A[T1​T2​​]

But this does not match your original scalar equations, because those involve the difference (T2−T1)(T_2 - T_1)(T2​−T1​) and the matrix AAA you constructed doesn’t reflect that form directly.

Re-deriving the matrix system consistent with scalar ODEs:

From the scalar equations:

dT1dt=−UA(T2−T1)=UA(T1−T2)\frac{dT_1}{dt} = - U A (T_2 - T_1) = U A (T_1 - T_2)dtdT1​​=−UA(T2​−T1​)=UA(T1​−T2​)dT2dt=UA(T2−T1)\frac{dT_2}{dt} = U A (T_2 - T_1)dtdT2​​=UA(T2​−T1​)

If you write in vector form T=[T1;T2]\mathbf{T} = [T_1; T_2]T=[T1​;T2​]:

dTdt=UA[−111−1]T\frac{d\mathbf{T}}{dt} = U A \begin{bmatrix} -1 & 1 \\ 1 & -1 \end{bmatrix} \mathbf{T}dtdT​=UA[−11​1−1​]T

where

A=hA(1ρcVccpc+1ρhsVhscphs)A = h A \left( \frac{1}{\rho_c V_c c_{p_c}} + \frac{1}{\rho_{hs} V_{hs} c_{p_{hs}}} \right)A=hA(ρc​Vc​cpc​​1​+ρhs​Vhs​cphs​​1​)

Alternatively, as two separate constants:

a=hAρcVccpc,b=hAρhsVhscphsa = \frac{h A}{\rho_c V_c c_{p_c}}, \quad b = \frac{h A}{\rho_{hs} V_{hs} c_{p_{hs}}}a=ρc​Vc​cpc​​hA​,b=ρhs​Vhs​cphs​​hA​

then

dT1dt=−a(T2−T1)\frac{dT_1}{dt} = - a (T_2 - T_1)dtdT1​​=−a(T2​−T1​)dT2dt=b(T2−T1)\frac{dT_2}{dt} = b (T_2 - T_1)dtdT2​​=b(T2​−T1​)

Note the different denominators imply a≠ba \neq ba=b, so the system is:

ddt[T1T2]=[a−a−bb][T1T2]\frac{d}{dt} \begin{bmatrix} T_1 \\ T_2 \end{bmatrix} = \begin{bmatrix} a & -a \\ -b & b \end{bmatrix} \begin{bmatrix} T_1 \\ T_2 \end{bmatrix}dtd​[T1​T2​​]=[a−b​−ab​][T1​T2​​]

No negative sign out front here. So the ODE system is:

dTdt=AT\frac{d\mathbf{T}}{dt} = A \mathbf{T}dtdT​=AT