1D Heat Conduction using explicit Finite Difference Method

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)

[KSSV]

First place, it is not giving any curve..there is a error in your code. Please recheck your code once.

[Derek Shaw]

I am not sure I understand correctly. In the above I wrote this equation to be iterated

With Boundary conditions

and Initial Conditions

with T0=400k and TL=Ti=300k

I am not sure how to set these boundary conditions in the code. Or if there is a curve I need to derive before doing the iterations.

Answer 1

[michio]

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.

[Derek Shaw]

Although your code didn't plot the graph properly for me it helped me figure out how to define the boundaries. I was able to plot it properly by switching around a few things. Thanks for your help!

[michio]

Glad to know that you figure things out. One additional tip is vectorization instead of for-loop, ie.

 for i=2:n-1
  T1(i) = T0(i) + alpha*(T0(i+1)-2*T0(i)+T0(i-1));
 end

is equivalent to

 T1(2:n-1) = T0(2:n-1) + alpha*(T0(3:n)-2*T0(2:n-1)+T0(1:n-2));

The later runs much faster.

[Kunpeng Liao]

Are the time step and grid spacing missing? I guess it would be:

if true
  T1(2:n-1) = T0(2:n-1) + alpha*dt/dx^2*(T0(3:n)-2*T0(2:n-1)+T0(1:n-2));
end

[Torsten]

Look at the definition of "alpha" in the code ...

Best wishes

Torsten.

[Siddarth Banerjee]

@Derek Shaw Can you please share the code. I am also struck in the same and cannot figure out the correct code for plotting..

Thanks

[Torsten]

 L=1;
 t=1;
 k=.001;
 n=11;
 nt=10000;
 dx=L/n;
 dt=.002;
 alpha=k*dt/dx^2;
 T0=400*ones(1,n+1);
 T1=300*ones(1,n+1);
 T0(1) = 300;
 T0(end) = 300;

 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)

[SYML2nd]

Hi @Torsten I have a question that regards this code.

If one end should be at 400K and the other end is at 300K.

Why, before the for cycle, is T0 = [300,400,..,400,300]?

For me it should have been T0=[400,300,...,300,300]. If T(0,t)=T0=400 K, then this boundary condition should be fixed during the whole calculation. Can you explain?

[Torsten]

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)

[SYML2nd]

Thank you @Torsten

Answer 2

[youssef aider]

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

Answer 3

[okncyln]

implicit?

Answer 4

[Alugunuri]

how to write code with neumqn BCs ?