1D Heat Conduction using Eulers Explicit discretisation
14 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Yashraj Randad
am 8 Aug. 2022
Kommentiert: William Rose
am 8 Aug. 2022
I am trying to solve a problem regarding heat conduction. I have written down a code and using the method that I know I tried to solve it but when I'm making the T0 matrix it gives me error because I have to start at 0 Kelvin and it gives me an error. Attached is an image of the question I'm trying to solve here.
%% Clearing
clc
clear
close all
%% Initialisation
L=1;
t=100;
k=0.01;
nt=500;
dx=0.05;
n=L/dx;
dt=0.5;
x=0:dx:L;
alpha=k*dt/dx^2;
T0=0*ones(1,n);
T1=20*ones(1,n);
T0(1) = 0;
T0(end) = 20;
%% Solving
for j=1:nt
for i=2:n-1
T1(i)=T0(i)+alpha*(T0(i+1)-2*T0(i)+T0(i-1));
end
T0=T1;
end
%% plotting
timeset=[0 1 5 10 100]; % time instants required
for i=1:length(timeset) % 1 to number of elements in timeset
z=timeset(i); % time instant
ntset=z/dt+1; % time instant extracted for the corresponding step
figure, plot(x,T1(:,ntset)) % plotting temperature vs distance
ylabel('Temperature (°C)') % labels y axis
xlabel('Distance (m)') % label x axis
title('Distribution of Temperature') % title of graph
grid minor % detailed view
end
0 Kommentare
Akzeptierte Antwort
William Rose
am 8 Aug. 2022
Note that length(x)=21, but your n=20. This causes problems. I recommend that you initialize a bit differently, as shown below.
I do not understand why you have vectors for T0 and T1. I recommend that T0 and T1 be scalars. I recommend that you create a single array, called T. Each row of T corresponds to one time. Each column of T corresponds to one position along the bar.
%% Initialisation
L=1.0; dx=0.05; x=0:dx:L; nx=length(x); %space
tmax=30; dt=0.05; t=0:dt:tmax; nt=length(t); %time
alpha=0.01; %thermal diffusivity, m^2/s
T0=0; %left temperature
T1=20; %right termperature
T=zeros(nt,nx); %allocate array for temperature
T(:,1)=T0; %make the left temp = T0 at all times
T(1,:)=T0; %make the initial temp=T0 at all locations
T(:,end)=T1; %make the right temp = T1 at all times
Solve the problem.
for i=2:nt
for j=2:nx-1
T(i,j)=T(i-1,j)+alpha*(dt/dx^2)*(T(i-1,j+1)-2*T(i-1,j)+T(i-1,j-1));
end
end
Plot some results.
plot(x,T(1,:),'-r',x,T(101,:),'-y',x,T(201,:),'-g',x,T(301,:),'-c',...
x,T(401,:),'-b',x,T(501,:),'-m',x,T(601,:),'-r');
grid on; title('Temperature versus position and time')
xlabel('Position (m)'); ylabel('Temperature (K)')
legend('t=0 s','t=5 s','t=10 s','t=15 s','t=20 s','t=25 s','t=30 s');
When I ran the code with dt=0.5 s, the solution oscillated, and gave results such as T=+/-10^167 after 100 seconds. This happens if the time step is too large. So I reduced the time step by a factor of 10: dt=0.05 s. Then I got reasonable results, as seen above.
Learn how to use surf() to make a nice 3D plot with axes position, time, and temperature.
2 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Geometry and Mesh finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!