# Heat Transfer: Matlab 2D Conduction Question

173 views (last 30 days)
Charles on 27 Mar 2012
Commented: ARJUN MODIA on 23 Jun 2020
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)
ARJUN MODIA on 23 Jun 2020
If we work on steady-state problem then the values of h & k will not affect your results. @ Jonathan and @ Bimal.

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)
Max on 7 Dec 2017
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!

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

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

### Categories

Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by