Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Steady state analysis & Transient State Analysis
we will implement the following methods for solving implicit equations.
Jacobi
Gauss-seidel
Successive over-relaxation
Steady and unsteady states
Steady state : In this fluid state in which there is no change in with respective fluid property and no time is known as steady state condition
Solving steady state heat conduction using three different iteravitive
- Jacobi
- Gaussc seidel
- successive over relaxation
The equation is derived from heat equation and by centrel differencing meathod
2D conduction equation is :
The 1st equation is jacobi method
The 2nd equation is Gauss seidel
The 3rd equation is successive over relaxation method
Main code
% script for solving 2d steady state heat equation
close all
clear all
clc
%Jacobi method
%inputs
L=1 %Length
nx=30; %Grids along x axis
ny=nx; %Grids along y axis
x=linspace(0,L,nx); %spacing along x axis
y=linspace(0,L,ny); %spacing along y axis
dx=x(2)-x(1);
dy=y(2)-y(1);
error=9e9;
tol=1e-4; %tolerance limit
nt=1:400; %change in time
alpha=1.8;
%BC
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=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(nx,2:ny-1)=T_B;
%calculating the average of a matrix
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;
[X,Y]=meshgrid(x,y);
Told=T;
tic;
%calculating of 2D steady state heat equation
iterative_solver=1:3 %(enter the iteration solver );
%%jacobi iteration
%time loop
for iterative_solver=1
tic;
if iterative_solver==1
%space loop
jacobi_iter=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_iter=jacobi_iter+1;
end
T_prev_dt=T;
end
F=toc;
%plot
figure(1)
[M,N]=contourf(X,Y,T);
clabel(M,N);
colormap(jet)
colorbar
title_text=sprintf('jacobi steady iter number = %d, Total time F=%fs',jacobi_iter,F);
title(title_text);
xlabel('X-axis')
ylabel('Y-asix')
fprintf('No of iter in jacobi=%d\n',jacobi_iter);
pause(0.03)
end
%Gauss seidel
%gs iter
%inputs
L=1;
nx=30;
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;
alpha=1.8;
%BC
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=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(nx,2:ny-1)=T_B;
%calculating the average of a matrix
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;
[X,Y]=meshgrid(x,y);
Told=T;
tic;
%calculating of 2D steady state heat equation
%%GS iteration
for iterative_solver=2
tic;
if iterative_solver==2
gs_iter=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_iter=gs_iter+1;
end
end
F=toc;
%plot
figure(2)
[M,N]=contourf(X,Y,T);
clabel(M,N);
colormap(jet)
colorbar
title_text=sprintf('GS iter number = %d, Total time F=%fs',gs_iter,F);
title(title_text);
xlabel('X-axis')
ylabel('Y-asix')
fprintf('No of iter in GS=%d\n',gs_iter);
pause(0.03)
end
%SOR iter
%inputs
%solving iteration 3
L=1;
nx=30;
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;
alpha=1.8;
%BC
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=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(nx,2:ny-1)=T_B;
%calculating the average of a matrix
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;
[X,Y]=meshgrid(x,y);
Told=T;
tic;
%calculating of 2D steady state heat equation
%%SOR iteration
for iterative_solver=3
tic;
if iterative_solver==3
sor_iter=1
while(error>tol)
for i=2:nx-1
for j=2:ny-1
T(i,j)=((1-alpha)*Told(i,j)+(alpha*0.25*(T(i,j-1)+Told(i,j+1)+T(i-1,j)+Told(i+1,j))));
end
end
error=max(max(abs(Told-T)));
Told=T;
sor_iter=sor_iter+1;
end
end
F=toc;
%plot
figure(3)
[M,N]=contourf(X,Y,T);
clabel(M,N);
colormap(jet)
colorbar
title_text=sprintf('SOR iter number = %d, Total time F=%fs',sor_iter,F);
title(title_text);
xlabel('X-axis')
ylabel('Y-asix')
fprintf('No of iter in SOR=%d\n',sor_iter);
pause(0.03)
end
Unsteady state (Transient state )
It is a term in which a fluid property changes with respective time
Main code
% script for solving 2d unsteady(transient) state heat equation
close all
clear all
clc
%jacobi
%inputs
L=1 %domain length
nx=30; %Grid point along x axis
ny=nx; %Grid point along y axis
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
tol=1e-4;
alpha=1.8;
dt = 1e-3;
nt=1:400;
%BC
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=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(nx,2:ny-1)=T_B;
%calculating the average of a matrix
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;
%we know the value of constant k
k1=(alpha*dt)/(dx^2);
k2= (alpha*dt)/(dy^2);
%[X,Y]=meshgrid(x,y);
Told=T;
T_prev_dt = Told;
% calculation 2D transient heat conduction equation iteration solver
iterative_solver=input('enter the iterative_solver=')
%jacobic iteration
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
tic;
jacobi_iter=1;
if iterative_solver==1
%time loop
for k=1:nt
error=9e9;
%space loop
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)= (term1*T_prev_dt(i,j))+term2*H+term3*V;
end
end
error = max(max(abs(Told-T)));
Told = T;
jacobi_iter= jacobi_iter+1;
end
T_prev_dt=T;
end
F=toc;
%plot
figure(1)
[M,N]=contourf(x,y,T);
clabel(M,N);
colormap(jet)
colorbar
title_text=sprintf('jacobi unsteady iter number = %d, Total time F=%fs',jacobi_iter,F);
title(title_text);
xlabel('X-axis')
ylabel('Y-asix')
fprintf('No of iter in jacobi=%d\n',jacobi_iter);
end
%GS Iteration solver
%inputs
L=1
nx=30;
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;
alpha=1.8;
dt = 1e-3;
%BC
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=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(nx,2:ny-1)=T_B;
%calculating the average of a matrix
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;
%we know the value of constant k
k1=(alpha*dt)/(dx^2);
k2= (alpha*dt)/(dy^2);
[X,Y]=meshgrid(x,y);
Told=T;
T_prev_dt = Told;
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
gs_iter=1;
if iterative_solver==2
%time loop
for k=1:nt
error=9e9;
while(error>tol);
%space loop
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,1+j));
T(i,j)= term1*T_prev_dt(i,j)+term2*H+term3*V;
end
end
error= max(max(abs(Told-T)));
Told=T;
gs_iter=gs_iter+1;
end
T_prev_dt+T;
end
F=toc;
%plot
figure(2)
[M,N]=contourf(X,Y,T);
clabel(M,N);
colormap(jet)
colorbar
title_text=sprintf('GS unsteady iter number = %d, Total time F=%fs',gs_iter,F);
title(title_text);
xlabel('X-axis')
ylabel('Y-asix')
fprintf('No of iter in jacobi=%d\n',gs_iter);
end
%sor iteration solver
%inputs
L=1
nx=30;
ny=nx;
x=linspace(0,L,nx);
y=linspace(0,L,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
tol=1e-4;
alpha=1.8;
dt = 1e-3;
nt=400;
%BC
T_L=400;
T_R=800;
T_T=600;
T_B=900;
T=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(nx,2:ny-1)=T_B;
%calculating the average of a matrix
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;
%we know the value of constant k
k1=(alpha*dt)/(dx^2);
k2= (alpha*dt)/(dy^2);
%[X,Y]=meshgrid(x,y);
Told=T;
T_prev_dt = Told;
%SOR iteration
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
tic;
sor_iter=1
if iterative_solver==3
%time loop
for k=1:nt
error=9e9;
while(error>tol);
%space time
for i = 2:nx-1
for j = 2:ny-1;
H = (alpha*T(i-1,j)+(1-alpha)*Told(i-1,j)+Told(i+1,j));
V = (alpha*T(i,j-1)+(1-alpha)*Told(i,j-1)+Told(i,j+1));
T(i,j)= T_prev_dt(i,j)*term1+(term2*H+term3*V);
end
end
error= max(max(abs(Told-T)));
Told=T;
sor_iter=sor_iter+1;
end
T_prev_dt=T;
end
F=toc;
%plot
figure(3)
[M,N]=contourf(X,Y,T);
clabel(M,N);
colormap(jet)
colorbar
title_text=sprintf('SOR unsteady iter number = %d, Total time F=%fs',sor_iter,F);
title(title_text);
xlabel('X-axis')
ylabel('Y-asix')
fprintf('No of iter in SOR=%d\n',sor_iter);
end