1D Heat Conduction using explicit Finite Difference Method
169 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello I am trying to write a program to plot the temperature distribution in a insulated rod using the explicit Finite Central Difference Method and 1D Heat equation. The rod is heated on one end at 400k and exposed to ambient temperature on the right end at 300k. I am using a time of 1s, 11 grid points and a .002s time step. When I plot it gives me a crazy curve which isn't right. I think I am messing up my initial and boundary conditions. Here is my code.
L=1;
t=1;
k=.001;
n=11;
nt=500;
dx=L/n;
dt=.002;
alpha=k*dt/dx^2;
T0(1)=400;
for j=1:nt
for i=2:n
T1(i)=T0(i)+alpha*(T0(i+1)-2*T0(i)+T0(i-1));
end
T0=T1;
end
plot(x,T1)
2 commentaires
KSSV
le 15 Déc 2016
First place, it is not giving any curve..there is a error in your code. Please recheck your code once.
Réponse acceptée
michio
le 15 Déc 2016
Modifié(e) : michio
le 15 Déc 2016
It seems your initial condition and boundary conditon (x = L) are missing in the code. Try
L=1;
t=1;
k=.001;
n=11;
nt=500;
dx=L/n;
dt=.002;
alpha=k*dt/dx^2;
T0=400*ones(1,n);
T1=300*ones(1,n);
T0(1) = 300;
T0(end) = 300;
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
plot(T1)
http://jp.mathworks.com/matlabcentral/fileexchange/59916-simple-heat-equation-solver could make a good example for you.
9 commentaires
Torsten
le 4 Déc 2022
Yes, code should be
L=1;
t=1;
k=.001;
n=11;
nt=10000;
dx=L/n;
dt=.002;
alpha=k*dt/dx^2;
T0=300*ones(1,n+1);
T0(1) = 400;
T1 = T0;
for j=1:nt
for i=2:n
T1(i)=T0(i)+alpha*(T0(i+1)-2*T0(i)+T0(i-1));
end
T0=T1;
end
plot((0:n)*L/n,T1)
Plus de réponses (3)
youssef aider
le 12 Fév 2019
here is one, you can just change the boundaries
clear
clc
clf
% domain descritization
alpha = 0.05;
xmin = 0;
xmax = 0.2;
N = 100;
dx = (xmax-xmin)/(N-1);
x = xmin:dx:xmax;
dt = 4.0812E-5;
tmax = 1;
t = 0:dt:tmax;
% problem initialization
phi0 = ones(1,N)*300;
phiL = 230;
phiR = phiL;
% solving the problem
r = alpha*dt/(dx^2) % for stability, must be 0.5 or less
for j = 2:length(t) % for time steps
phi = phi0;
for i = 1:N % for space steps
if i == 1 || i == N
phi(i) = phiL;
else
phi(i) = phi(i)+r*(phi(i+1)-2*phi(i)+phi(i-1));
end
end
phi0 = phi;
plot(x,phi0)
shg
pause(0.05)
end
0 commentaires
Voir également
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!