Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM: Solving the steady state and unsteady state 2D heat conduction by using the point iterative technique.
- Jacobian method
- Gauss siedol method
- Successive-over Relaxation method
With absolute error criteria 1e-4.
- STEADY STATE ANALYSIS
The 2D steady state heat conduction equation with no heat generation is
Discretizing the above equation
Where
TP=Ti,j">TP=Ti,j
TL=Ti-1,j">TL=Ti−1,j
TR=Ti+1,j">TR=Ti+1,j
TT=Ti,j+1">TT=Ti,j+1
TB=Ti,j-1">TB=Ti,j−1
On further simplication we get
Where,
To solve these we use the ieterative solvers.
1. Jacobian iterative solver :
The 2D steady state heat conduction equation for this method is
2. Gauss-siedel iterative solver :
The 2D steady state heat conduction equation for this method is
3. Successive-over Relaxation(SOR) iterative solver :
The 2D steady state heat conduction equation for this method is
for scaling jacobian method,where Tj is the result of jacobian method,
and for scaling Gauss-Seidel method,where Tj is the result of Gauss-Seidel method,
when omega = 0 we get new value is equal to old value which gives wrong result.
when omega goes from 0 to 1 we are switching from no update to an update using gauss siedal or jacobi method.
when omega >1 then we are taking the update from gauss siedal method and scaling up the convergence rate. The value we give is called over relaxition factor.
when omega <1 ,The value we give is called under relaxition factor.This is used when the solution given by gauss siedel is over relaxed we give omega value less than 1 to get the converged solution.
MAIN PROGRAM
% solving the steady state 2D heat conduction problem
close all
clear all
clc
%inputs
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=30;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=ones(nx,ny);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%creating a copy of T
[X, Y]=meshgrid(x,y);
Told=T;
% jacobi method
for iterative_solver=1;
tic
if iterative_solver==1;
jacobi_iterations=1;
while (error >tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=0.25*(Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1));
end
end
error=max(max(abs(Told-T)));
Told=T;
jacobi_iterations = jacobi_iterations+1;
end
end
F=toc;
%plotting
figure(1);
[P,Q]=contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('Jacobi method,Iteration number=%d,computation time F=%fs',jacobi_iterations,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf('jacobi=%d/n',jacobi_iterations);
pause(0.03);
end
% for Gauss siedel method
for iterative_solver = 2;
% update values
%inputs
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=30;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=300*ones(nx,ny);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
tic
if iterative_solver ==2;
GS_iteration = 1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j)=0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1));
end
end
error = max(max(abs(Told-T)));
Told = T
GS_iteration = GS_iteration + 1;
end
end
F = toc;
%plotting
figure(2);
[P,Q]=contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('Gauss seidol method,Iteration number=%d,computation time F=%fs',GS_iteration,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf('jacobi=%d/n',GS_iteration);
pause(0.03);
end
% SOR method
for iterative_solver = 3;
% update value
% inputs
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=30;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=ones(nx,ny);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%creating a copy of T
[X, Y]=meshgrid(x,y);
Told=T;
tic;
if iterative_solver ==3;
SOR_iteration = 1;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=((1-omega)*(Told(i,j))+(omega*0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1))));
end
end
error = max(max(abs(Told-T)));
Told = T;
SOR_iteration = SOR_iteration+1;
end
end
F = toc;
%plotting
figure(3);
[P,Q]=contourf(X,Y,T);
clabel(P,Q);
colormap(jet);
colorbar;
title_text=sprintf('SOR method,Iteration number=%d,computation time F=%fs',SOR_iteration,F);
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
fprintf('jacobi=%d/n',SOR_iteration);
pause(0.03);
end
NO. OF ITERATION NUMBER IN JACOBI METHOD = 1884
NO. OF ITERATION NUMBER IN GAUSS SIEDEL METHOD = 954
NO. OF ITERATION NUMBER IN SOR METHOD = 93
By comparing above plots i.e jacobi,gs & sor method, sor method is solving faster with less iterations and compution time when compared to two other methds.
SOR is taking 93 iterations and computation time =0.009769s
Gauss seidel is taking 954 iterations and computation time=0.04521s
Jacobi method is taking 1884 iterations and computation time=0.11530s
Solving the unsteady/transient 2D heat conduction
1. Explicit Method
MAIN PROGRAM
% solving the steady 2D transient/unsteady state heat conduction problem
% explicit
clear all
close all
clc
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=51;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1e-3; %convective coefficient
t=200;
%CFL number
M = alpha*dt/(dx^2);
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=300*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
nt=t/dt;
%2d transient unsteady state heat conduction equation in explicit method
tic;
for k=1:nt
for i=2:nx-1
for j=2:ny-1
term1=Told(i,j);
term2=k1*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j));
term3=k2*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1));
T(i,j)=term1+term2+term3;
end
end
Told=T;
end
time=toc;
%plottings
figure(1);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('Explicit method,computation time =%fs,M=%f',time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.3);
2. Implicit Method
MAIN PROGRM
% solving the steady 2d transient unsteady state heat conduction problem
% implicit
clear all
close all
clc
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=20;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1; %convective coefficient
t=100;
%CFL number
M = alpha*dt/(dx^2);
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=300*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
nt=t/dt;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
T_prev_dt=Told;
%2d transient unsteady state heat conduction equation in explicit method
nt=t/dt;
iterative_solver=1;
tic;
jacobi_iterations=1;
for v=1:nt
if iterative_solver==1
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=100;
while(error >tol)
for i=2:nx-1
for j=2:ny-1
H=(Told(i-1,j)+Told(i+1,j));
V=(Told(i,j-1)+Told(i,j+1));
T(i,j)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
jacobi_iterations=jacobi_iterations+1;
end
T_prev_dt=T;
end
time=toc;
%plottings
figure(1);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('Jacobi method,Iteration number=%d,computation time =%fs,M=%f',jacobi_iterations,time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.03);
end
%for gs method
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=20;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1; %convective coefficient
t=200;
% %CFL number
M = alpha*dt/(dx^2);
%define boundary conditions
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=300*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
% %calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
nt=t/dt;
% %creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
T_prev_dt=Told;
%2d transient unsteady state heat conduction equation in explicit method
iterative_solver=2;
tic;
gs_iterations=1;
for w=1:nt
if iterative_solver==2
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=100;
while(error >tol)
for i=2:nx-1
for j=2:ny-1
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
gs_iterations=gs_iterations+1;
end
T_prev_dt=T;
end
F=toc;
%plottings
figure(2);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('Gauss siedal,Iteration number=%d,computation time =%fs,M=%f',gs_iterations,time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.03);
end
%for SOR method
L=1; %length of domain and it is a square
%nx,ny no of grid points
nx=20;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4;
dt=0.1;
alpha = 1; %convective coefficient
t=100;
% %CFL number
M = alpha*dt/(dx^2);
%calculation of relaxation factor
omega=2/(1+sin(pi*dx));
%define boundary conditions
T_L=400;% T_R=800;
T_T=600;
T_B=900;
T=300*ones(nx,ny);
% temperature at boundary conditions
T(:,1)=T_L;
T(:,end)=T_R;
T(end,:)=T_T;
T(1,:)=T_B;
%calculation of temperature distribution explicity using iterative
%solvers
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
nt=t/dt;
%creating a copy of T
[X Y]=meshgrid(x,y);
Told=T;
T_prev_dt=Told;
%2d transient unsteady state heat conduction equation in explicit method
iterative_solver=3;
tic;
SOR_iterations=1;
for Z=1:nt
if iterative_solver==3
term1=(1+2*k1+2*k2)^(-1);
term2=k1*term1;
term3=k2*term1;
error=100;
while(error >tol)
for i=2:nx-1
for j=2:ny-1
H=(omega*T(i-1,j)+(1-omega)*Told(i-1,j)+Told(i+1,j));
V=(omega*T(i,j-1)+(1-omega)*Told(i,j-1)+Told(i,j+1));
T(i,j)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
SOR_iterations=SOR_iterations+1;
end
T_prev_dt=T;
end
F=toc;
%plottings
figure(3);
[P,Q] = contourf(X,Y,T);
clabel(P,Q,"fontsize",12);
colormap(jet);
colorbar;
title_text=sprintf('SOR,Iteration number=%d,computation time =%fs,M=%f',SOR_iterations,time,M');
title(title_text);
xlabel('Xaxis');
ylabel('Yaxis');
pause(0.3);
end
Jacobo method iteration number = 4200
Gauss siedel method iteration number = 3798
SOR method iteration number = 1588
From the above plots you can conclude that explicit method takes less time to execute but when it comes to stability it mainly depend on CFL condition if this condition is not satisfied the solution gets blow up.Where as in the implicit methods jacobian method gives the solution in 4200 iterations where as gauss method takes 3798 iterations but using SOR method you get the solution in 1588 iterations but the computational tine high for SOR and then jacobian but for the gauss method is is faster than other two.Here in SOR method if you use omega value more or less than 1.06 you can see the iteration count increases when it crosses 1.71 then you see that the solution got blew up.so its value should be given less than 2.you can also observe that for any values of alpha dx and dt the solution will be stable, like explicit it is not dependent of CFL condition.