Heat Transfer: Matlab 2D Conduction Question

1. The problem statement, all variables and given/known data

Having trouble with code as seen by the gaps left where it asks me to add things like the coefficient matrices. Any Numerical Conduction matlab wizzes out there?

A long square bar with cross-sectional dimensions of 30 mm x 30 mm has a specied temperature on each side, The temperatures are:

Tbottom = 100 C
Ttop = 150 C
Tleft = 250 C
Tright = 300 C

Assuming isothermal surfaces, write a software program to solve the heat equation to determine the two-dimensional steady-state spatial temperature distribution within the bar. Your analysis should use a finite difference discretization of the heat equation in the bar to establish a system of equations:

2. Relevant equations

AT = C

Must use Gauss-Seidel Method to solve the system of equations

3. The attempt at a solution

    clear all
  close all
  %Specify grid size
  Nx = 10;
  Ny = 10;
  %Specify boundary conditions
  Tbottom = 100
  Ttop = 150
  Tleft = 250
  Tright = 300
  % initialize coefficient matrix and constant vector with zeros
  A = zeros(Nx*Ny);
  C = zeros(Nx*Ny,1);
  % initial 'guess' for temperature distribution
  T(1:Nx*Ny,1) = 100;
  % Build coefficient matrix and constant vector
  % inner nodes
  for n = 2:(Ny-1)
  for m = 2:(Nx-1)
  i = (n-1)*Nx + m;
  A(i,i+Nx) = 1;
  A(i,i-Nx) = 1;
  A(i,i+1) = 1;
  A(i,i-1) = 1;
  A(i,i) = -4;
  end
  end
  % Edge nodes
  % bottom
  for m = 2:(Nx-1)
  %n = 1
  i = m;

A(i,i+Nx) = 1;
A(i,i+1) = 1;
A(i,i-1) = 1;
A(i,i) = -4;
C(i) = -Tbottom;

end
%top:
for m = 2:(Nx-1)
% n = Ny
i = (Ny-1)*Nx + m;

A(i,i-Nx) = 1;
A(i,i+1) = .5;
A(i,i-1) = .5;
A(i,i) = -5;
C(i) = -Ttop;

end
%left:
for n=2:(Ny-1)
%m = 1
i = (n-1)*Nx + 1;
A(i,i+Nx) = .5;
A(i,i+1) = 1;
A(i,i-Nx) = .5;
A(i,i) = -2;
end
%right:
for n=2:(Ny-1)
%m = Nx
i = (n-1)*Nx + Nx;
A(i,i+Nx) = 1;
A(i,i-1) = 1;
A(i,i-Nx) = 1;
A(i,i) = -4;
C(i) = -Tright;

DEFINE COEFFICIENT MATRIX AND CONSTANT VECTOR ELEMENTS HERE
end
% Corners
%bottom left (i=1):
i=1
A(i,Nx+i) = 1;
A(i,2) = 1;
A(i,1) = -4;
C(i) = -(Tbottom + Tleft);
%bottom right:
i = Nx
A(i,Nx+i) = 1;
A(i,2) = 1;
A(i,1) = -4;
C(Nx) = -(Tbottom + Tright);
%top left:
i = (Ny-1)*Nx + 1;
A(i,i+1) = .5
A(i,i) = 
%top right:
i = Nx*Ny;
DEFINE COEFFICIENT MATRIX AND CONSTANT VECTOR ELEMENTS HERE
%Solve using Gauss-Seidel
residual = 100;
iterations = 0;
while (residual > 0.0001) % The residual criterion is 0.0001 in this example
7% You can test different values
iterations = iterations+1
%Transfer the previously computed temperatures to an array Told
Told = T;
%Update estimate of the temperature distribution
INSERT GAUSS-SEIDEL ITERATION HERE
%compute residual
deltaT = abs(T - Told);
residual = max(deltaT);
end
iterations % report the number of iterations that were executed
%Now transform T into 2-D network so it can be plotted.
delta_x = 0.03/(Nx+1)
delta_y = 0.03/(Ny+1)
for n=1:Ny
for m=1:Nx
i = (n-1)*Nx + m;
T2d(m,n) = T(i);
x(m) = m*delta_x;
y(n) = n*delta_y;
end
end
T2d
surf(x,y,T2d)
figure
contour(x,y,T2d)

[Walter Roberson]

Please be more specific about the difference between what you observe and what you expect.

[Charles]

I am not sure what to expect, as this is my first time working with matlab. I just edited the code adding in what I think is the correct coefficient matrices, do the values look right?

[Walter Roberson]

I am not familiar with heat transfer equations, so I have no idea if the coefficients "look right".

I _am_ familiar with many MATLAB error messages, and with tracing back _specific_ problems to causes (given inputs), but there isn't much I can help with unless you can clarify

Having trouble with code as seen by the gaps left where it asks me to add things like the coefficient matrices"

[Geoff]

Most people here are probably not involved in your course, so we don't necessarily know what the correct values are. We are here to help with MatLab questions or issues, not to verify your results. I expect your results can be used to reproduce the initial conditions. That is verification. May I also suggest (if you doubt the correctness of your code) that you set breakpoints at relevant places (F12), run it, and check your intermediate values for sanity. Use F10 to step through code one line at a time, F5 to resume, Shift-F5 to cancel....

[Charles]

hold on i think im onto something. be right back

[Charles]

okay I cranked it out below, its working, but the distributions of the heat looks kind of odd. It is so overwhelming hotter at the top. Idk if anyone knows how to check for wrong coefficient matrice values or something, please give me a tip :)

[Image Analyst]

You mean like how to use the debugger? Like how to stop at breakpoints and how you can "check for wrong coefficient matrice values or something"? You can click in the margin to set a breakpoint and, after you've stopped at it, click in a variable and type control-D to check it's value in the variable editor, or look in the workspace panel.

[bimal]

What if we have values for coefficient of conductivity (k) and convection coefficient (h). How can we modify codes

[Jonathan Ayala]

Hello, any luck with modifying the code with values for K and h? I am working on a similar problem, Thanks.

[ARJUN MODIA]

If we work on steady-state problem then the values of h & k will not affect your results. @ Jonathan and @ Bimal.

Answer 1

[Charles]

clear all
close all
%Specify grid size
Nx = 10;
Ny = 10;
%Specify boundary conditions
Tbottom = 100
Ttop = 150
Tleft = 250
Tright = 300
% initialize coefficient matrix and constant vector with zeros
A = zeros(Nx*Ny);
C = zeros(Nx*Ny,1);
% initial 'guess' for temperature distribution
T(1:Nx*Ny,1) = 100;
% Build coefficient matrix and constant vector
% inner nodes
for n = 2:(Ny-1)
for m = 2:(Nx-1)
i = (n-1)*Nx + m;
A(i,i+Nx) = 1;
A(i,i-Nx) = 1;
A(i,i+1) = 1;
A(i,i-1) = 1;
A(i,i) = -4;
end
end
% Edge nodes
% bottom
for m = 2:(Nx-1)
%n = 1
i = m;

A(i,i+Nx) = 1;
A(i,i+1) = 1;
A(i,i-1) = 1;
A(i,i) = -4;
C(i) = -Tbottom;

end
%top:
for m = 2:(Nx-1)
% n = Ny
i = (Ny-1)*Nx + m;

A(i,i-Nx) = 1;
A(i,i+1) = 1;
A(i,i-1) = 1;
A(i,i) = -4;
C(i) = -Ttop;

end
%left:
for n=2:(Ny-1)
%m = 1
i = (n-1)*Nx + 1;
A(i,i+Nx) = 1;
A(i,i+1) = 1;
A(i,i-Nx) = 1;
A(i,i) = -4;
end
%right:
for n=2:(Ny-1)
%m = Nx
i = (n-1)*Nx + Nx;
A(i,i+Nx) = 1;
A(i,i-1) = 1;
A(i,i-Nx) = 1;
A(i,i) = -4;
C(i) = -Tright;

end
% Corners
%bottom left (i=1):
i=1;
A(i,Nx+i) = 1;
A(i,2) = 1;
A(i,1) = -4;
C(i) = -(Tbottom + Tleft);
%bottom right:
i = Nx;
A(i,i+Nx) = 1;
A(i,i-1) = 1;
A(i,i) = -4;
C(i) = -(Tbottom + Tright);
%top left:
i = (Ny-1)*Nx + 1;
A(i,i+1) = 1;
A(i,i) = -4;   
A(i,i-Nx) = 1;
C(i) = -(Ttop + Tleft);
%top right:
i = Nx*Ny;
A(i,i-1) = 1;
A(i,i) = -4;
A(i,i-Nx) = 1;
C(i) = -(Tright + Ttop);
%Solve using Gauss-Seidel
residual = 100;
iterations = 0;
while (residual > 0.0001) % The residual criterion is 0.0001 in this example
% You can test different values
iterations = iterations+1;
%Transfer the previously computed temperatures to an array Told
Told = T;
%Update estimate of the temperature distribution

for n=1:Ny
for m=1:Nx
i = (n-1)*Nx + m;
Told(i) = T(i);
end
end
% iterate through all of the equations
for n=1:Ny
for m=1:Nx
i = (n-1)*Nx + m;
%sum the terms based on updated temperatures
sum1 = 0;
for j=1:i-1
sum1 = sum1 + A(i,j)*T(j);
end
%sum the terms based on temperatures not yet updated
sum2 = 0;
for j=i+1:Nx*Ny
sum2 = sum2 + A(i,j)*Told(j);
end
% update the temperature for the current node
T(i) = (1/A(i,i)) * (C(i) - sum1 - sum2);
end
end
residual = max(T(i) - Told(i));
end

%compute residual
deltaT = abs(T - Told);
residual = max(deltaT);

iterations; % report the number of iterations that were executed
%Now transform T into 2-D network so it can be plotted.
delta_x = 0.03/(Nx+1);
delta_y = 0.03/(Ny+1);
for n=1:Ny
for m=1:Nx
i = (n-1)*Nx + m;
T2d(m,n) = T(i);
x(m) = m*delta_x;
y(n) = n*delta_y;
end
end
T2d;
surf(x,y,T2d)
figure
contour(x,y,T2d)

[Charles]

okay I cranked it out, its working but the distributions of the heat looks kind of odd. It is so overwhelming hotter at the top. Idk if anyone knows how to check for wrong coefficient matrice values or something, please give me a tip :)

[Frank Mota]

the code wont run for me , is there something im doing wrong?

[Max]

Charles: It's been a while since you posted this, and perhaps you've found the error - but you forgot to include the C vector update in your left edge, which is why the distribution looked odd!

Answer 2

[Ahmed Hakim]

Very nice Code

I would like to use SOR method for finding the optimum omega...can u help me?

thanks

Answer 3

[ashwath suresh]

Hi could you please explain how codes for inner nodes and edge nodes are given? why is the value for A(i,i)=-4