How do I make/solve a diffrential equation with two variables

Hello, I'm trying to solve this model (see sketch). This leads to a diffrential equation with three unkonws being T1,T2 and Ts. Ive gotten it to a point where I can solve by stating a fixed value for T2 and computing Ts by taking the average of T1 and T2. This is wrong because Ts is time and place dependent and T2 will evventualy heat up as well. So the next step is implementing place dependency for T in my model and.

This is where I'm stuck. I have very little experience with matlab and have no idea on how to make this work.

This is the current code I have:

clear all
close all
clc

% Solve thermal first order difrential equation rho*V*cp DT/dt = sum of
% heat fluxes and plot temprature in function of time

rho = 0.0022 ;       %density in g/mm³
k=1.38        ;      %thermal conductivty in mW/(mm K)
cp=0.772  ;
h=0.01;             % heat transfer coeficient in mW/(mm²K)
Tinf=22;            % nominal temprature of surouniding
L=16.8;          %lenght of block in mm height is the direction conduction moves in
H=1.31;         %height of the block in mm
W=7;             %width of the block in mm
Aw=H*W;        % surface trough which the conduction flows in mm²
As=W*L;         % surface where convecgtion takes place in mm²
V=Aw*L;         % Volume of the block
C=rho*cp*V;

Qin=0.0216;             %heat input in W

nP=3;                   %number of periods you want to plot
P=14.4;                 %Period
scan=0.05;              %input time       
t=[0 P*nP];      

sw=@(t) sawtooth(t*2*pi/P)/2+0.5<(scan/P);  %make a square wave
OPT=odeset('MaxStep',1e-3);

[t,T]=ode45(@(t,T) (Qin*sw(t)-(1e-3*k*Aw*(T-Tinf)/L)-(1e-3*h*As*((T+Tinf)/2-Tinf)))/C,t,Tinf,OPT);

figure;
hold on
plot(t,sw(t)*Qin)
plot(t,T-Tinf)          % in order to put everyting on the zero axis

I would love to hear feedback and suggestions.

Thanks in advance,

Robin

[James Tursa]

Please post the differential equations you are trying to solve, along with any constants and initial conditions.

[Robin Muylle]

@James Tursa. This is the equation that matlab is solving: dT/dt=Qin*sw(t)-(1e-3*k*Aw*(T-Tinf)/L)-(1e-3*h*As*((T+Tinf)/2-Tinf)))/C. With initial conditions of T0 being 22°C.

Ive come to understand that in order to solve this problem with place dependency I'll need to solve for a partial diffrential equation. Which would be C dT/dt=a d²T/dl² +sum of heat in-outputs

(C and a are constants)

The same initial conditions apply and for boundry conditions I need to make the left side equal to the heat input and the right side no heat flux allowed (isolate).

I'm trying to make this work at the moment.

[Torsten]

Use "pdepe" for this task.

Best wishes

Torsten.

[Robin Muylle]

Hey Torsten,

I'm trying to do that. But I get very strange results. The temprature doesn't seem to be affected by time and it is highest at the end of the block. Ive tried flipping the boundry conditions aroud (changing pl to pr etc..) But that has no effect. Do you see what's wrong?

clear all
close all
clc

global rho k cp h Tinf L H W Aw As V C a Qin sw 

rho = 0.0022 ;       %density in g/mm³
k=1.38        ;      %thermal conductivty in mW/(mm K)
cp=0.772  ;
h=0.01;             % heat transfer coeficient in mW/(mm²K)
Tinf=22;            % nominal temprature of surouniding
L=16.8;          %lenght of block in mm height is the direction conduction moves in
H=1.31;         %height of the block in mm
W=7;             %width of the block in mm
Aw=H*W;        % surface trough which the conduction flows in mm²
As=W*L;         % surface where convecgtion takes place in mm²
V=Aw*L;         % Volume of the block
C=rho*cp*V;
a=k/(rho*cp);

Qin=0.0216;             %heat input in W

nP=3;                   %number of periods you want to plot
P=14.4;                 %Period
scan=0.05;              %input time

sw=@(t) sawtooth(t*2*pi/P)/2+0.5<(scan/P);  %make a square wave
OPT=odeset('MaxStep',1e-3);

t=linspace(0,P*nP,10);
l=linspace(0,L,10);

sol=pdepe(0,@pdefun,@pdefunic,@pdefunbc,l,t);
T=sol(:,:,1);
surf(t,l,T)
xlabel('time');
ylabel('lenght');

%
function[C,f,s]=pdefun(l,t,T,DTDl);
global rho cp V k sw Qin 
C=rho*cp*V;
a=k/(rho*cp);
f=-a*DTDl;
s=Qin*sw(t);
end
%
function T0 = pdefunic(l)
global Tinf
T0=Tinf;   
end
%
function [pl,ql,pr,qr]=pdefunbc(ll,Tl,lr,Tr,t)
global Qin Aw k sw
pr=Qin*sw(t)/(Aw*k);
qr=1;
pl=0;
ql=1;
end