Phil Dias Midterm Project 2D Conduction
Problem:
Solve the conduction equation using a Transient solver and a Steady state solver using Iterative techniques (Jacobi,Gauss Seidal,SOR)
1. Assume that the domain is a unit square.
2. Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
3. Boundary conditions for steady and transient case
1. Left:400K
2. Top:600K
3. Right:800K
4. Bottom:900K
Main Octave Code:
% 2D Heat Conduction Steady and Unsteady State
clear all
close all
clc
% Initial Conditions
L = 1; % Length [m]
nx = 10; % Nodes in x direction
ny = nx; % Nodes in y direction
% Mesh setup
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = L/(nx-1);
dy = L/(ny-1);
% Tolerance criteria
error = 9e9;
tol = 1e-4;
t = 5 ;
dt = 1e-2;
alpha = 0.3;
% Steady State Constant
k = 2*(1/(dx^2)+1/(dy^2));
%Unsteady State Constants
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);
% Defining BC in [K]
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
% Assigning Boundary Temperatures
T(:,1) = T_L;
T(:,end) = T_R;
T(1,:) = T_T;
T(end,:) = T_B;
% Boundary Temperatures
T(1,1) = (T_L + T_T)/2;
T(1,end) = (T_B + T_R)/2;
T(end,1) = (T_T + T_R)/2;
T(end,end) = (T_B + T_L)/2;
% Steady State (SS) Loop Selection
for iteration_selection = [1 2 3]
if iteration_selection == 1
time = jacobi_SS(nx,ny,x,y,dx,dy,T,tol,k);
else if iteration_selection == 2
time = GS_SS(nx,ny,x,y,dx,dy,T,tol,k);
else if iteration_selection == 3
time = SOR_SS(nx,ny,x,y,dx,dy,T,tol,k);
end
end
end
end
% Transient State Loop Selection
% Implicit Method
for iteration_selection = [1 2 3]
if iteration_selection == 1
time = jacobi_implicit(nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2);
else if iteration_selection == 2
time = GS_implicit(nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2);
else if iteration_selection == 3
time = SOR_implicit(nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2);
end
end
end
end
% Explicit Method
for iteration_selection = [1 2 3]
if iteration_selection == 1
time = jacobi_explicit(nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2);
else if iteration_selection == 2
time = GS_explicit(nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2);
else if iteration_selection == 3
time = SOR_explicit(nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2);
end
end
end
end
Steady State Functions:
% Jacobi Method Steady State
function [out] = jacobi_SS (nx,ny,x,y,dx,dy,T,tol,k)
meshgrid(x,y);
T_old = T;
jacobi_iteration = 1;
error = 9e9;
tic
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T_old(i-1,j) + T_old(i+1,j)) / (k*(dx^2));
T2 = (T_old(i,j-1) + T_old(i,j+1)) / (k*(dy^2));
T(i,j) = (T1 + T2);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
jacobi_iteration = jacobi_iteration + 1;
end
iteration_time = toc;
% Plot for Jacobi Method Steady State
figure(1)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'Jacobi Method Iterations = %d\', jacobi_iteration);
title2 = sprintf(\'Jacobi Method Time = %f sec\', iteration_time);
title3 = sprintf(\'Jacobi Method Error = %d\', error);
colorbar
title({\'Jacobi Method Steady State\'; title1;title2;title3});
out = iteration_time
endfunction
% Gauss-Seidel Method Steady State
function [out] = GS_SS (nx,ny,x,y,dx,dy,T,tol,k)
meshgrid(x,y);
T_old = T;
GS_iteration = 1;
error = 9e9;
tic
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T(i-1,j) + T_old(i+1,j)) / (k*(dx^2));
T2 = (T(i,j-1) + T_old(i,j+1)) / (k*(dy^2));
T(i,j) = (T1 + T2);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
GS_iteration = GS_iteration + 1;
end
iteration_time = toc;
% Plot for Gauss-Seidel Method Steady State
figure(2)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'GS Method Iterations = %d\', GS_iteration);
title2 = sprintf(\'GS Method Time = %f sec\', iteration_time);
title3 = sprintf(\'GS Method Error = %d\', error);
colorbar
title({\'GS Method Steady State\'; title1;title2;title3});
out = iteration_time
endfunction
% Successive Over Relxation (SOR) Method Steady State
function [out] = SOR_SS (nx,ny,x,y,dx,dy,T,tol,k)
% Over Relxation Factor
ORF = 1.5;
meshgrid(x,y);
T_old = T;
SOR_iteration = 1;
error = 9e9;
tic
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T(i-1,j) + T_old(i+1,j)) / (k*(dx^2));
T2 = (T(i,j-1) + T_old(i,j+1)) / (k*(dy^2));
T(i,j) = (1-ORF)*T_old(i,j) + ORF*(T1 + T2);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
SOR_iteration = SOR_iteration + 1;
end
iteration_time = toc;
% Plot for SOR Method Steady State
figure(3)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'SOR Method Iterations = %d\', SOR_iteration);
title2 = sprintf(\'SOR Method Time = %f sec\', iteration_time);
title3 = sprintf(\'SOR Method Error = %d\', error);
colorbar
title({\'SOR Method Steady State\'; title1;title2;title3});
out = iteration_time
endfunctionImplicit Method for Transient State Functions:
% Jacobi Method Transient State Implict Form
function [out] = jacobi_implicit (nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2)
meshgrid(x,y);
T_old = T;
Tp = T;
nt = t/dt;
jacobi_iteration = 1;
tic
for m = 1:nt
error = 9e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T_old(i-1,j) + T_old(i+1,j));
T2 = (T_old(i,j-1) + T_old(i,j+1));
T(i,j) = (Tp(i,j) + (k1*T1) + (k2*T2))/ (1 + 2*k1 +2*k2);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
jacobi_iteration = jacobi_iteration + 1;
end
Tp = T;
end
iteration_time = toc;
% Plot for Implicit Jacobi Method Transient State
figure(4)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'Implict Jacobi Method Iterations = %d\', jacobi_iteration);
title2 = sprintf(\'Implict Jacobi Method Time = %f sec\', iteration_time);
title3 = sprintf(\'Implict Jacobi Method Error = %d\', error);
colorbar
title({\'Implicit Jacobi Method Transient State\'; title1;title2;title3});
out = iteration_time
endfunction
% Implict Gauss-Seidel Method Transient State
function [out] = GS_implicit (nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2)
meshgrid(x,y);
T_old = T;
Tp = T;
nt = t/dt;
GS_iteration = 1;
tic
for m = 1:nt
error = 9e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T(i-1,j) + T_old(i+1,j));
T2 = (T(i,j-1) + T_old(i,j+1));
T(i,j) = (Tp(i,j) + (k1*T1) + (k2*T2))/ (1 + 2*k1 +2*k2);
end
end
error = max(max(abs(T_old - T)));
T_old = T;
GS_iteration = GS_iteration + 1;
end
Tp = T;
end
iteration_time = toc;
% Plot for Implict GS Method Transient State
figure(5)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'Implict GS Method Iterations = %d\', GS_iteration);
title2 = sprintf(\'Implict GS Method Time = %f sec\', iteration_time);
title3 = sprintf(\'Implict GS Method Error = %d\', error);
colorbar
title({\'Implict GS Method Transient State\'; title1;title2;title3});
out = iteration_time
endfunction
% Implict SOR Method Transient State
function [out] = SOR_implicit (nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2)
ORF = 1.09;
meshgrid(x,y);
T_old = T;
Tp = T;
nt = t/dt;
SOR_iteration = 1;
tic
for m = 1:nt
error = 9e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T(i-1,j) + T_old(i+1,j));
T2 = (T(i,j-1) + T_old(i,j+1));
T(i,j) = (1-ORF)*T_old(i,j) + (ORF*((Tp(i,j) + (k1*T1) + (k2*T2))/ (1 + 2*k1 +2*k2)));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
SOR_iteration = SOR_iteration + 1;
end
Tp = T;
end
iteration_time = toc;
% Plot for Implict GS Method Transient State
figure(6)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'Implict SOR Method Iterations = %d\', SOR_iteration);
title2 = sprintf(\'Implict SOR Method Time = %f sec\', iteration_time);
title3 = sprintf(\'Implict SOR Method Error = %d\', error);
colorbar
title({\'Implict SOR Method Transient State\'; title1;title2;title3});
out = iteration_time
endfunction
Explicit Method for Transient State Functions:
% Explicit Jacobi Method Transient State
function [out] = jacobi_explicit (nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2)
meshgrid(x,y);
T_old = T;
Tp = T;
nt = t/dt;
jacobi_iteration = 1;
tic
for m = 1:nt
error = 9e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T_old(i-1,j) + T_old(i+1,j));
T2 = (T_old(i,j-1) + T_old(i,j+1));
T(i,j) = (Tp(i,j)*(1 - 2*k1 -2*k2) + (k1*T1) + (k2*T2));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
jacobi_iteration = jacobi_iteration + 1;
end
Tp = T;
end
iteration_time = toc;
% Plot for Explicit Jacobi Method Transient State
figure(7)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'Explict Jacobi Method Iterations = %d\', jacobi_iteration);
title2 = sprintf(\'Explict Jacobi Method Time = %f sec\', iteration_time);
title3 = sprintf(\'Explict Jacobi Method Error = %d\', error);
colorbar
title({\'Explicit Jacobi Method Transient State\'; title1;title2;title3});
out = iteration_time
endfunction
% Explict Gauss-Seidel Method Transient State
function [out] = GS_explicit (nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2)
meshgrid(x,y);
T_old = T;
Tp = T;
nt = t/dt;
GS_iteration = 1;
tic
for m = 1:nt
error = 9e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T(i-1,j) + T_old(i+1,j));
T2 = (T(i,j-1) + T_old(i,j+1));
T(i,j) = ((Tp(i,j)*(1 - 2*k1 - 2*k2)) + (k1*T1) + (k2*T2));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
GS_iteration = GS_iteration + 1;
end
Tp = T;
end
iteration_time = toc;
% Plot for Explict GS Method Transient State
figure(8)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'explict GS Method Iterations = %d\', GS_iteration);
title2 = sprintf(\'explict GS Method Time = %f sec\', iteration_time);
title3 = sprintf(\'explict GS Method Error = %d\', error);
colorbar
title({\'Explict GS Method Transient State\'; title1;title2;title3});
out = iteration_time
endfunction
% Explict SOR Method Transient State
function [out] = SOR_explicit (nx,ny,x,y,dx,dy,T,t,dt,tol,k,k1,k2)
ORF = 1.3;
meshgrid(x,y);
T_old = T;
Tp = T;
nt = t/dt;
SOR_iteration = 1;
tic
for m = 1:nt
error = 9e9;
while (error > tol)
for i = 2:nx-1
for j = 2:ny-1
T1 = (T(i-1,j) + T_old(i+1,j));
T2 = (T(i,j-1) + T_old(i,j+1));
T(i,j) = (1-ORF)*T_old(i,j) + (ORF*((Tp(i,j)*(1 - 2*k1 - 2*k2)) + (k1*T1) + (k2*T2)));
end
end
error = max(max(abs(T_old - T)));
T_old = T;
SOR_iteration = SOR_iteration + 1;
end
Tp = T;
end
iteration_time = toc;
% Plot for Explict GS Method Transient State
figure(9)
[a,b] = contourf(x,y,T);
colormap(jet)
clabel(a,b);
xlabel(\'X Domain\');
ylabel(\'Y Domain\');
title1 = sprintf(\'Explict SOR Method Iterations = %d\', SOR_iteration);
title2 = sprintf(\'Explict SOR Method Time = %f sec\', iteration_time);
title3 = sprintf(\'Explict SOR Method Error = %d\', error);
colorbar
title({\'Explict SOR Method Transient State\'; title1;title2;title3});
out = iteration_time
endfunction
Results:
Steady State:
In the State State from, the SOR method uses the least amount of iterations, 30, and is therefore most efficient compared to Jacobi, 217, and GS, 117, methods. Variations of the Over Relation Factors (ORF) were used to determine its most efficient value.
Implicit Transient State:
Since the Transient state equation is more complex, it require many more iterations in all three methods. A trend can be noticed in comparison to the State State results in that the SOR method took the least iterations.
Explicit Transient State:
When comparing Implicit and Explicit solutions, it is noted that Explicit required a significantly lower amount of iterations. The trend contunues for this solutions and they show SOR to be the most efficient method with the least number of iterations. Since the ORF can be easily changed, the SOR methods gains an advantage over the preceding methods and is the fastest and most accurate method.