Help programing 2D conduction heat transfer in time, using finite diference method (forward euler for time, centered euler for space).

25 Points Grid. r = x, z = y.

The formulas are correct, but i don't kwon what i am doing wrong with the matlab programming. Please help >.<

%Body properties

H = 0.3; % Height [m]

D = 0.2; % Diameter [m]

Alpha = 1.5; %Difusion heat transfer coefficient [m^2/s]

% Factores para Diferencia Finita

dr = D/8; 

dz = H/8; 

dt = 7; 

tp = 3600; %Total Process Time. [s]

tc = 2000; %Proces tranfer from heating to cooling. [s]

Ti = 18; %Initial body temperature. [°C]

Tc = 120; %Heating Temperature. [°C]

Tf = 8; %Cooling Temperature. [°C]

%Grid Initialization

T(:,:,1)= Ti*ones(5,5);

time(1) = 0;

time(2) = dt;

t=1;

while time(t+1)<=tp 

  for i=[5,4,3,2,1]

    for j=[5,4,3,2,1]

      if i==5&&time(t+1)<=tc

                              T(i,j,t+1) = Tc;

                          elseif j==5&&time(t+1)<=tc

                              T(i,j,t+1) = Tc;

                          elseif i==5&&time(t+1)>tc

                              T(i,j,t+1) = Tf;

                          elseif j==5&&time(t+1)>tc

                              T(i,j,t+1) = Tf;

                          elseif i>1&&i<5&&j>1&&j<5

                              T(i,j,t+1) = Alpha*dt*((T(i+1,j,t)-2*T(i,j,t)+T(i-1,j,t))/(dr)^2 + (1/((i-1)*dr))*((T(i+1,j,t)-T(i,j,t))/(2*dr)) + (T(i,j+1,t)-2*T(i,j,t)+T(i,j-1,t))/(dz)^2) + T(i,j,t); 

                          elseif i==1&&j>1&&j<5

                              T(i,j,t+1) = Alpha*dt*(4*(T(i+1,j,t)-T(i,j,t))/(dr)^2 + (T(i,j+1,t)-2*T(i,j,t)+T(i,j-1,t))/(dz)^2) + T(i,j,t); 

                          elseif i<1&&i>5&&j==1

                              T(i,j,t+1) = Alpha*dt*((T(i+1,j,t)-2*T(i,j,t)+T(i-1,j,t))/(dr)^2 + (1/((i-1)*dr))*((T(i+1,j,t)-T(i,j,t))/(2*dr)) + 2*(T(i,j+1,t)-T(i,j,t))/(dz)^2) + T(i,j,t);

                          elseif i==1&&j==1

                              T(i,j,t+1) = 2*Alpha*dt*(2*(T(i+1,j,t)-T(i,j,t))/(dr)^2 + (T(i,j+1,t)-T(i,j,t))/(dz)^2) + T(i,j,t);

      end

    end

  end

  t = t+1;

  time(t+1) = time(t)+dt;

end

[Nataly Challapa]

Hello! My question is about the graphic, how can I do a graphic in 3D of this solution? it is possible? because this method would be a solution for the heat equation. Is not it right?

Answer 1

[Dr. Seis]

Do you get any errors... or the result is just incorrect?

One problem: Your second to last "elseif" statement has "i<1&&i>5&&j==1" when I think you meant "i>1&&i<5&&j==1" since "i" can never be both less than 1 and greater than 5.

Answer 2

[Walter Manns]

Thank for your helpfor finding that mistake.

I'm still getting errors though. Fot t=1 my matrix goes ok. For t=2 it should start calculating, but it only change the frontier values.

From 3 and more it just goes to hell :(

Ac

val(:,:,1) =

    18    18    18    18    18
    18    18    18    18    18
    18    18    18    18    18
    18    18    18    18    18
    18    18    18    18    18

val(:,:,2) =

    18    18    18    18   120
    18    18    18    18   120
    18    18    18    18   120
    18    18    18    18   120
   120   120   120   120   120

val(:,:,3) =

1.0e+006 *

    0.0000    0.0000    0.0000    0.7616    0.0001
    0.0000    0.0000    0.0000    0.7616    0.0001
    0.0000    0.0000    0.0000    0.7616    0.0001
    1.9992    1.9992    1.9992    2.7608    0.0001
    0.0001    0.0001    0.0001    0.0001    0.0001

[Dr. Seis]

Let's start with your value at val(4,4,3)... I assume this value should actually be somewhere between 18 and 120 after 1 time step, right?

[Dr. Seis]

Can you provide the non-code version of the equation that you based your "i>1&&i<5&&j>1&&j<5" case?

Answer 3

[Walter Manns]

Grant, thank you so much for helping me.

I finnaly found what i did wrong. The alpha i was using was huge in comparison to a real one (it should be in the order of 10^(-7)), so my ecuations went unstable. After reducing it to the real order, i started geting good results :D

Thank you again :D