clear,pack,clc
% material properties
k=200
cp=921
rho=2500
alpha=k/(cp*rho);
% introduce discretization constants
% time
delta_t=0.1;
%space
delta_L=0.01; % this is in meters
% how many points
N=100;
% constants a1, and a2
a1=1-2*(alpha*delta_t)/(delta_L^2);
a2=(alpha*delta_t)/(delta_L^2);
% introduce initial and boundary conditions
% initial condition
T0=zeros(N,1);
T0(1:floor(N/3))=30;
% boundary condition
% this is the left point next to the bar
T0_m1=500;
% this is the right point next to the bar
T0_Np1=10;
b=zeros(N,1);
b(1)=a2*T0_m1;
b(N)=a2*T0_Np1;
% define the A matrix
A=zeros(N,N);
for i=1:N
if i==1
A(1,1)=a1; A(1,2)=a2;
end
if i==N
A(N,N)=a1; A(N,N-1)=a2;
end
if (i>=2)&&(i<=N-1)
A(i,i)=a1;
A(i,i-1)=a2;
A(i,i+1)=a2;
end
end
% number of simulation steps
sim_steps=10000;
RESULT=zeros(N,sim_steps);
for j=1:sim_steps
if j==1
RESULT(:,1)=T0; % initial temperature
else
RESULT(:,j)=A*RESULT(:,j-1)+b;
end
end
imagesc(RESULT)
% RESULT is number of points X time
% middle point plotting
plot(RESULT(50,:))
hold on
plot(RESULT(80,:),'r')
