Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Objective is to solve the steady and unsteady 2D heat conduction equation by using iterative solvers {Jacobi , Gauss Seidel and Successive Over relaxation (SOR) method
Domain of interest is a square surface
Boundary conditions are given as follows:
Temperature at left surface = 400K
Temperature at right surface = 600K
Temperature at top surface = 800K
Temperature at bottom surface = 900K
First I have discretised the domain in (20 by 20) matrix
Assuming the temperature at all the other nodes to be 300K
Also the temperature at corner points are the average of connecting surfaces.
A generalised 3D transient heat conduction equation equation looks like as follows:
Since there is no internal heat generation term(assumption) and we are dealing with only 2D case ,the 3rd and fourth term of left hand side of the generalised equation becomes zero.
So, for steady state equation we are left with
Now we can solve the above equation by discretization
here we will use central differencing technique
So,solving the equation using Jacobi method
Using Gauss Seidel
Using SOR method {which is the enhanced version of Gauss Seidel method}
The MATLAB code is written here to show the results;
close all
clear all
clc
%%Solving the steady state 2D heat conduction
% Input parameters
% No of grid points
nx = 20;
ny = nx;
omega = 1.2;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = dx;
%Absolute error criteria
error = 9e9;
tolerance = 1e-4;
%% Defining boundary conditions
T_L = 400; %temperature of left surface
T_T = 600; %temperature of top surface
T_R = 800; %temperature of right surface
T_B = 900; %temperature of bottom surface
T = 200*ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:nx-1) = T_B;
%calculating average temperature at corners
T(1,1) = (T_T+T_L)/2;
T(nx,ny) = (T_R+T_B)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,1) = (T_L+T_B)/2;
%%assigning values to T
T_initial = T;
T_old = T;
%calculating by Jacobi method
iterative_solver = 1;
if iterative_solver == 1
jacobi_iteration = 1;
while (error>tolerance)
for i = 2:nx-1
for j=2:ny-1
T(i,j)=0.25*(T_old(i-1,j)+T_old(i+1,j)+T_old(i,j-1)+T_old(i,j+1));
end
end
error = max(max(abs(T_old-T)));
T_old = T;
jacobi_iteration = jacobi_iteration + 1;
end
end
%plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('No of Jacobi iterations(Implicit)=%d',jacobi_iteration));
pause(0.03)
%calculating by Gauss seidel method
iterative_solver = 2;
if iterative_solver == 2
Gauss_seidel_iteration = 1;
while (error>tolerance)
for i = 2:nx-1
for j=2:ny-1
T(i,j)=0.25*(T(i-1,j)+T_old(i+1,j)+T(i,j-1)+T_old(i,j+1));
end
end
error = max(max(abs(T_old-T)));
T_old = T;
Gauss_seidel_iteration = Gauss_seidel_iteration + 1;
end
end
%plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('No of Gauss seidel iterations(Implicit)=%d',Gauss_seidel_iteration));
pause(0.03)
%calculating by SOR method
iterative_solver = 3;
if iterative_solver == 3
SOR = 1;
while (error>tolerance)
for i = 2:nx-1
for j=2:ny-1
T(i,j)=omega*0.25*(T(i-1,j)+T_old(i+1,j)+T(i,j-1)+T_old(i,j+1)) + (1-omega)*T_old(i,j);
end
end
error = max(max(abs(T_old-T)));
T_old = T;
SOR = SOR + 1;
end
end
%plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('No of SOR iterations(Implicit)=%d',SOR));
pause(0.03)
%% Solving the unsteady state 2D heat conduction equation
2D unsteady state equation looks like as follows:
Upon discretizing the above equation we get:
Now the three schemes has been applied to solve this equation
And the MATLAB code and results are given below:
close all
clear all
clc
%%Solving the unsteady state 2D heat conduction
% Input parameters
% No of grid points
nx = 20;
ny = nx;
nt = 1200;
omega = 1.2;
dt = 1e-3;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = dx;
%Absolute error criteria
error = 9e9;
tolerance = 1e-4;
%% Defining boundary conditions
T_L = 400; %temperature of left surface
T_T = 600; %temperature of top surface
T_R = 800; %temperature of right surface
T_B = 900; %temperature of bottom surface
T = 200*ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:nx-1) = T_B;
%calculating average temperature at corners
T(1,1) = (T_T+T_L)/2;
T(nx,ny) = (T_R+T_B)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,1) = (T_L+T_B)/2;
%%assigning values to T
T_initial = T;
T_old = T;
%%calcualtion of 2D steady heat equation by Jacobi method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
jacobi_iteration = 1;
for k = 1:nt
error = 9e9;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
term_1 = 1/(1+(2*k1) + (2*k2));
term_2 = k1*term_1;
term_3 = k2*term_1;
H = (T_old(i-1,j))+(T_old(i+1,j));
V = (T_old(i,j+1))+(T_old(i,j-1));
T(i,j) = (T_initial(i,j)*term_1)+(H*term_2)+(V*term_3);
end
end
error = max(max(abs(T_old-T)));
T_old = T;
jacobi_iteration = jacobi_iteration + 1;
end
T_initial = T;
%plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('No of unsteady Jacobi iterations(Explicit)=%d',jacobi_iteration));
pause(0.03)
end
%%calcualtion of 2D steady heat equation by Gauss seidel method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
Gauss_seidel_iteration = 1;
for k = 1:nt
error = 9e9;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
term_1 = 1/(1+(2*k1) + (2*k2));
term_2 = k1*term_1;
term_3 = k2*term_1;
H = (T(i-1,j))+(T_old(i+1,j));
V = (T(i,j+1))+(T_old(i,j-1));
T(i,j) = (T_initial(i,j)*term_1)+(H*term_2)+(V*term_3);
end
end
error = max(max(abs(T_old-T)));
T_old = T;
Gauss_seidel_iteration = Gauss_seidel_iteration + 1;
end
T_initial = T;
%plotting
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('No of unsteady Gauss seidel iterations(Explicit)=%d',Gauss_seidel_iteration));
pause(0.03)
end
close all
clear all
clc
%%Solving the unsteady state 2D heat conduction
% Input parameters
% No of grid points
nx = 20;
ny = nx;
nt = 1200;
omega = 1.2;
dt = 1e-3;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2)-x(1);
dy = dx;
%Absolute error criteria
error = 9e9;
tolerance = 1e-4;
%% Defining boundary conditions
T_L = 400; %temperature of left surface
T_T = 600; %temperature of top surface
T_R = 800; %temperature of right surface
T_B = 900; %temperature of bottom surface
T = 200*ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:nx-1) = T_B;
%calculating average temperature at corners
T(1,1) = (T_T+T_L)/2;
T(nx,ny) = (T_R+T_B)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,1) = (T_L+T_B)/2;
%%assigning values to T
T_initial = T;
T_old = T;
T_end = T;
%%calcualtion of 2D steady heat equation by SOR method
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
SOR_iteration = 1;
for k = 1:nt
error = 9e9;
while(error>tolerance)
for i=2:nx-1
for j=2:ny-1
term_1 = 1/(1+(2*k1) + (2*k2));
term_2 = k1*term_1;
term_3 = k2*term_1;
H = (T(i-1,j))+(T_old(i+1,j));
V = (T(i,j+1))+(T_old(i,j-1));
T(i,j) = (T_initial(i,j)*term_1)+(H*term_2)+(V*term_3);
T_end(i,j) = ((1-omega)*T_old(i,j)) + (T(i,j)*omega);
end
end
error = max(max(abs(T_old-T)));
T_old = T_end;
SOR_iteration = SOR_iteration + 1;
end
%plotting
T_initial = T_end;
figure(1)
contourf(x,y,T)
clabel(contourf(x,y,T))
colorbar
colormap(jet)
set(gca,'ydir','reverse')
xlabel('X-axis')
ylabel('Y-axis')
title(sprintf('No of unsteady SOR iterations(Explicit)=%d',SOR_iteration));
pause(0.03)
end
