Heat Transfer: Matlab 2D Conduction Question

173 views (last 30 days)
Charles
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)
  10 Comments
ARJUN MODIA
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.

Sign in to comment.

Answers (3)

Charles
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)
  3 Comments
Max
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!

Sign in to comment.


Ahmed Hakim
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
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

Community Treasure Hunt

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

Start Hunting!

Translated by