Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Heat Transfer: Matlab 2D Conduction Question

Asked by Charles on 27 Mar 2012

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)

7 Comments

Charles on 27 Mar 2012

hold on i think im onto something. be right back

Charles on 28 Mar 2012

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 on 28 Mar 2012

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.

Charles

Products

No products are associated with this question.

2 Answers

Answer by Charles on 28 Mar 2012
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)

1 Comment

Charles on 28 Mar 2012

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 :)

Charles
Answer by Ahmed Hakim on 17 Nov 2012

Very nice Code

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

thanks

0 Comments

Ahmed Hakim

Contact us