Heat Transfer: Matlab 2D Conduction Question

Charles (view profile)

on 27 Mar 2012
Latest activity Answered by ashwath suresh

ashwath suresh (view profile)

on 16 Feb 2015

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

Charles

Charles (view profile)

on 27 Mar 2012

hold on i think im onto something. be right back

Charles

Charles (view profile)

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

Image Analyst (view profile)

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.

Products

No products are associated with this question.

Charles (view profile)

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

Charles

Charles (view profile)

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

Ahmed Hakim (view profile)

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 (view profile)

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