Finite difference temperature distribution with TDMA

A stainless steel (thermal conductivity, k = 20 W/mK) fin has a circular cross-sectional area (diameter, D=4 cm) and length L=16 cm. The wing is attached to a wall with a temperature of 200C. The ambient temperature of the fluid surrounding the fin is 20C and the heat transfer coefficient is h=10 W/m2K. Fin tip is insulated. Determine the following by applying finite differences with TDMA.

1. Temperatures inside the fin and the amount of heat radiated from the fin. Plot the temperature distribution.

2. Discuss the effect of the number of element dimensions for N = 5, 10, 50, and 100. Plot the temperature distribution.

The temperature distribution from the analitic solution is like this: T1=150.54, T2=111.107, T3=78.671, T4=50.741, T5=25.172, my code does not give these values.

clc
clear all
L=0.16;             %lenght 
N=6;                %number of nodes
dx=L/(N-1);         %lenght between nodes
T=zeros(N+1,1);
Tb=200;             %the wall temperature (boundary condition)
k=6;                %number of iterations

for j=1:1:k
    T(1,1)=Tb;
    T(5,1)=T(7,1);  %the second boundary condition due to the isulation on the tip of the fin
  for i=2:1:N
    T(i,1)=(T(i+1,1)+T(i-1,1)+1.536)/2.0768;    %FDE narrowed down
  end
  T(N,1)=T(N-1,1);
  plot(T);
  hold on
end
hold off

If I can get some help with my code.

[Torsten]

What is TDMA ?

[Cubii4]

It is a tridiagonal matrix algorithm and to solve it we use the Thomas algorithm.

[Torsten]

And what equation (including the boundary conditions) are you trying to solve with your code ?

[Cubii4]

The first boundary condition is that the fin tip is insulated (Neumann condition) becaus of that (if N is the number of nodes) : (when there are five nodes).

The second boundary condition is that the wall temperature (Tb) is equal to the temperature at the the zeroth node (T0), T0=Tb.

The ODE (from the general fin equation) is :

The FDE is :

I hope I explained with enough details.

[Torsten]

What is P/A ? 4/D ?

[Cubii4]

I do not really know what P/A and 4/D is. If you could explain a little bit I would be happy to answer it then.

[Torsten]

You defined D as the diameter of the fin.

P/A appears in your general fin equation.

[Cubii4]

Yes of course I'm sorry, P/A would be 6/D (A=2*4 and P=(2+4)*2).

[Torsten]

Isn't 6/D for a sphere and not for a cylinder ?

I tried to compute the analytical solution for your problem, but could not reproduce the temperatures you included. Maybe you can tell me where I'm wrong.

syms x T(x)
eqn = diff(T,x,2)*20 == 2/0.02*10*(T-20);
dT = diff(T,x);
conds = [T(0)==200,dT(0.16)==0];
sol(x) = dsolve(eqn,conds)

sol(x) = 

double(sol(0.16))  % temperature at insulated side 

ans = 125.1865

double(sol(0)) % temperature at wall side

ans = 200

[Cubii4]

You are right, I just noticed that I did not used the equation for a cylinder (P=π*D).

How did you get the 2/0.02*10*(T-20) part, I couldn't quite understand it.

[Torsten]

How did you get the 2/0.02*10*(T-20) part, I couldn't quite understand it.

4/D * h * (T-T_inf)

[Cubii4]

I did the analytical solution on paper with P/A as 4/D and I got T1=155.2548, T2=117.4346, T3=84.6031, T4=55.0792, T5=27.3514. Something must be off because T5 to be 125.1865 as Your analytical solution shows is too much I think.

[Torsten]

If the equation is

k * d^2T/dx^2 = 4/D * h * (T-T_inf)
T(0) = 200
dT/dx (0.16) = 0

with the parameters for D, T_inf, k and h you posted, then I believe in MATLAB's "dsolve" solution.

[Cubii4]

I did it again and I got the same analytical solution as you did:

T1=171.3794, T2=150.5095, T3=136.3216, T4=128.0894, T5=125.3914

Now my FDE is

I tried this equation in my code but it is still not giving my the same as the analytical solution.

[Cubii4]

clc
clear all
L=0.16;
N=6;
dx=L/(N-1);
T=zeros(N+1,1);
Tb=200;
k=1000;

for i=1:1:k
    T(1,1)=Tb;
    for j=2:1:N
    T(j,1)=(T(j+1,1)+T(j-1,1)+1.024)/2.0512;
    end
    T(N+1)=T(N-1);
    plot(T);
    hold on
end
hold off
for j=1:N
    disp(T(j,1));
end

I just corrected the code like this and it gave me the temperature results as the analytic solutions. I can't thank you enough, is there anything else I have to correct in my code.

[Cubii4]

Thank you very much I really appreciated Your help. Can you write it as an answer so I can accept it.

Answer 1

[Torsten]

You solve the system iteratively. But you were told to solve it with the Thomas-Algorithm. So apply this algorithm to the matrix A below.

k = 20;
T_inf = 20;
h = 10;
D = 0.04;
L = 0.16;
T_at_wall = 200;
N = 6;
dL = L/(N-1);
A = zeros(N+1);
b = zeros(N+1,1);
A(1,1) = 1.0;
b(1) = T_at_wall;
for i = 2:N
    A(i,i-1) = k/dL^2;
    A(i,i) = -(2*k/dL^2 + 4/D*h);
    A(i,i+1) = k/dL^2;
    b(i) = -4/D*h*T_inf;
end
A(N+1,N-1) = -1.0;
A(N+1,N+1) = 1.0;
b(N+1) = 0;
T = A\b;
T = T(1:N)

T = 6×1

200.0000 171.3794 150.5095 136.3216 128.0894 125.3914

[EGE]

can you please explain the line after "end"? I could not figure out what it is. Thanks in advance.

[Torsten]

It's the discretized insulated Neumann condition at x=L: 0 = dT/dx ~ T_(n+1)-T_(n-1).

And the point (n+1) is a ghost node in the discretization.