Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Heat conduction equation is define by the pde's
`(delta T)/(delta t) = alpha *( (delta ^ 2T )/(delta x^2) + (delta ^ 2T )/(delta y^2)+ (delta ^ 2T )/(deltaz^2))`
Heat distribution in 3d space over time peroid
`alpha=` thermal diffusivity
Given that,nod grid sise is equal in x and y direction and the domain is assumed to be an unit squre.
Diffferant types iteratives techqunies
1.Jobic Method
first we compute the assume initial values at all nodes of domain and calulate new values , which will updated once we find the all nodes values and loop runs untill sol becomes converge.
2.Gauss-seidel
In this method the values obtained in Jacobi iteration are used for the next equation immediately without waiting for completion of first cycle.
3. Successive Over-relaxation (SOR)
This method is an updated version of previous two. In addition we use the relaxation factor to ease the process. For this analysis we use the relaxation facor to be greater than 1 and less than 2.
Heat conduction equation is define by the pde's
`(delta T)/(delta t) = alpha *( (delta ^ 2T )/(delta x^2) + (delta ^ 2T )/(delta y^2)+ (delta ^ 2T )/(deltaz^2))`
For 2D Steady state
`(delta ^ 2T )/(delta x^2) + (delta ^ 2T )/(delta y^2)=0`
The steady state equation is discretised using 2nd order CDS:
`(T L − 2 T + T R)/(Deltax^ 2) + (T B − 2 T + T T)/(Deltay^ 2) = 0`
On simplifying,
`T p = (1/ k)* ( (T L + T R)/(Deltax^2 ) +(T B + T T)/(Deltax^2 ) )`
`k = 2 ((Deltax^2 + Deltay^2)/( Deltax^2 Deltay^2))`
Code for Steady State
%code for the steady state (implicity)
clear all
close all
clc
%input parameters
L=1; %Domain Length
nx=20; %no.of gird size in X direction
ny=nx; %no.of gird size in X direction
x= linspace(0,L,nx); % X axis Mesh
y= linspace(0,L,ny); % Y axis Mesh
time=5; %time
omega = 1.2; % 1< omega >2
%grid range
dx = x(2)-x(1);
dy = y(2)-y(1);
% Initial Tolerance and Error guess values
error=9e9;
tol=1e-4;
[xx,yy] = meshgrid(x,y);
% Matrix defining
T=ones(nx,ny);
% Initial temp values
T_L=400;
T_R=800;
T_T=600;
T_B=900;
% defining Boundary condition
T(1,2:ny-1)=T_T; %top face
T(nx,2:ny-1)=T_B; %Bottom Face
T(2:nx-1,1)=T_L; %Left Face
T(2:nx-1,ny)=T_R; %Right Face
%definig the corner values
T(1,1)= (T_T+T_L)/2;
T(1,ny)= (T_T+T_R)/2;
T(nx,1)= (T_B+T_L)/2;
T(nx,ny)= (T_B+T_R)/2;
Told=T;
Tintial=T;
%solver Method
method= ' which method ?';
solver=input(method); % '1' for jocobi
% '1' for Gauss-sidel
% '1' for jocobi
if solver==1
iter_1=1;
tic
%convergence loop
while(error>tol)
for i=2:nx-1 % loop for x direction
for j=2:ny-1 % loop for Y dircetion
T(i,j) = (Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1))/4;
end
end
time1=toc;
error= max(max(abs(Told-T)));
Told=T;
iter_1=iter_1+1;
[x,y]= contourf(xx,yy,T);
colormap(jet)
colorbar
clabel(x,y)
title_text=sprintf('No.jocobi itertion=%d',iter_1);
title(title_text);
pause(0.000001)
end
end
%Gauss-siedel method
if solver==2
iter_2=1;
tic
%convergence loop
while(error>tol)
for i=2:nx-1 % loop for x direction
for j=2:ny-1 % loop for Y dircetion
T(i,j) = (T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1))/4;
end
end
time2=toc;
error= max(max(abs(Told-T)));
Told=T;
iter_2=iter_2+1;
[x,y]= contourf(xx,yy,T);
colormap(jet)
colorbar
clabel(x,y)
title_text=sprintf('No.Gauss Siedel itertion=%d',iter_2);
title(title_text);
pause(0.000001)
end
end
% Successive over relaxation(sor) method
if solver==3
iter_3=1;
tic
%convergence loop
while(error>tol)
for i=2:nx-1 % loop for x direction
for j=2:ny-1 % loop for Y dircetion
Tsmaple(i,j) = (T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1))/4;
T(i,j) = (1-omega)*(Told(i,j))+omega*Tsmaple(i,j);
end
end
time3=toc;
error= max(max(abs(Told-T)));
Told=T;
iter_3=iter_3+1;
[x,y]= contourf(xx,yy,T);
colormap(jet)
colorbar
clabel(x,y)
title_text=sprintf('No.of Successive over relaxation itertion=%d',iter_3);
title(title_text);
pause(0.000001)
end
endGraph
Results
form the grph we can say the no of iteration required to converge the solution more in the Jocobic iteration followed by the Gassu seidel and SOR
which means the simluation run time more in the jocobic method followed by the Gassu seidel and SOR
TRANSIENT STATE HEAT CONDUCTION:
In 2D heat conduction under unsteady state conditions
`(delta T)/(delta t) = alpha *( (delta ^ 2T )/(delta x^2) + (delta ^ 2T )/(delta y^2))`
The Transient state equation is discretised by First order FDS for the time derivative and 2nd order central differance for the space derivative.
For Explicit method,
`T^( n + 1) = T ^(n) + k 1 ( T L − 2 T + T R ) ^(n) + k 1 ( T B − 2 T + T T ) ^(n)`
For Implicit method,
`T^( n + 1) = 1/(1+2k1+2k2)[T ^(n) + k 1 ( T L + T R ) ^(n+1) + k 1 ( T B + T T ) ^(n+1)]`
n is for explicit schemes
n+1 is for implicit schemes
`k1=alpha*dt/dx^2`
`k1=alpha*dt/dy^2`
code for transient state (implicity)
%code for the unsteady state (implicity)
clear all
close all
clc
%input parameters
L=1; %Domain Length
nx=20; %no.of gird size in X direction
ny=nx; %no.of gird size in X direction
x= linspace(0,L,nx); % X axis Mesh
y= linspace(0,L,ny); % Y axis Mesh
time=5; %time
omega = 1.2; % 1< omega >2
alpha = 1.4;
dt=0.001;
nt=time/dt;
%grid range
dx = x(2)-x(1);
dy = y(2)-y(1);
% Initial Tolerance and Error guess values
error=9e9;
tol=1e-4;
[xx,yy] = meshgrid(x,y);
% Matrix defining
T=ones(nx,ny);
% Initial temp values
T_L=400;
T_R=800;
T_T=600;
T_B=900;
% defining Boundary condition
T(1,2:ny-1)=T_T; %top face
T(nx,2:ny-1)=T_B; %Bottom Face
T(2:nx-1,1)=T_L; %Left Face
T(2:nx-1,ny)=T_R; %Right Face
%definig the corner values
T(1,1)= (T_T+T_L)/2;
T(1,ny)= (T_T+T_R)/2;
T(nx,1)= (T_B+T_L)/2;
T(nx,ny)= (T_B+T_R)/2;
Told=T;
Tintial=T;
k1=alpha*dt/dx^2;
k2=alpha*dt/dy^2;
%solver Method
method= ' which method ?';
solver=input(method); % '1' for jocobi
% '2' for Gauss-sidel
% ' 3' for SOR
if solver==1
t1=(1+2*k1+2*k2)^(-1);
t2=k1*t1;
t3=k2*t1;
iter_1=1;
nt=time/dt;
tic
for m=1:nt
%convergence loop
error=1e9;
while(error>tol)
for i=2:nx-1 % loop for x direction
for j=2:ny-1 % loop for Y dircetion
h=Told(i-1,j)+Told(i+1,j);
v=Told(i,j-1)+Told(i,j+1);
T(i,j)=(t1)*(Tintial(i,j))+(h*t2)+(v*t3);
end
end
error= max(max(abs(Told-T)));
Told=T;
iter_1=iter_1+1;
end
Tintial=T;
end
time1=toc;
[x,y]= contourf(xx,yy,T);
colormap(jet)
colorbar
clabel(x,y)
title_text=sprintf('No.of jocobi itertion=%d',iter_1);
title(title_text);
%pause(0.001)
end
%Gauss-siedel method
if solver==2
t1=(1+(2*k1)+(2*k2))^(-1);
t2=k1*t1;
t3=k2*t1;
iter_2=1;
tic
for m=1:nt
error=1e9;
%convergence loop
while(error>tol)
for i=2:nx-1 % loop for x direction
for j=2:ny-1 % loop for Y dircetion
h=T(i-1,j)+Told(i+1,j);
v=T(i,j-1)+Told(i,j+1);
T(i,j)=(t1)*(Tintial(i,j))+(h*t2)+(v*t3);
end
end
error= max(max(abs(Told-T)));
Told=T;
iter_2=iter_2+1;
end
Tintial=T;
end
time2=toc;
[x,y]= contourf(xx,yy,T);
colormap(jet)
colorbar
clabel(x,y)
title_text=sprintf('No.of Gauss Siedel itertion=%d',iter_2);
title(title_text);
%pause(0.000001)
end
% Successive over relaxation(sor) method
if solver==3
t1=(1+2*k1+2*k2)^(-1);
t2=k1*t1;
t3=k2*t1;
iter_3=1;
tic
for m=1:nt
error=1e9;
%convergence loop
while(error>tol)
for i=2:nx-1 % loop for x direction
for j=2:ny-1 % loop for Y dircetion
h=T(i-1,j)+Told(i+1,j);
v=T(i,j-1)+Told(i,j+1);
Tsmaple(i,j)=(t1)*(Tintial(i,j))+(h*t2)+(v*t3);
T(i,j) = (1-omega)*(Told(i,j))+omega*Tsmaple(i,j);
end
end
error= max(max(abs(Told-T)));
Told=T;
iter_3=iter_3+1;
end
Tintial=T;
end
time3=toc;
[x,y]= contourf(xx,yy,T);
colormap(jet)
colorbar
clabel(x,y)
title_text=sprintf('No.of Successive over relaxation itertion=%d',iter_3);
title(title_text);
%pause(0.000001)
end
Graphs
Code for theTransient State 2D (Explicit):
%code for the unsteady state (Explicity)
clear all
close all
clc
%input parameters
L=1; %Domain Length
nx=20; %no.of gird size in X direction
ny=nx; %no.of gird size in X direction
x= linspace(0,L,nx); % X axis Mesh
y= linspace(0,L,ny); % Y axis Mesh
%grid range
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.4;
dt=0.0004;
%nt=time/dt;
% Initial Tolerance and Error guess values
error=9e9;
tol=1e-4;
[xx,yy] = meshgrid(x,y);
% Matrix defining
T=ones(nx,ny);
% Initial temp values
T_L=400;
T_R=800;
T_T=600;
T_B=900;
% defining Boundary condition
T(1,2:ny-1)=T_T; %top face
T(nx,2:ny-1)=T_B; %Bottom Face
T(2:nx-1,1)=T_L; %Left Face
T(2:nx-1,ny)=T_R; %Right Face
%definig the corner values
T(1,1)= (T_T+T_L)/2;
T(1,ny)= (T_T+T_R)/2;
T(nx,1)= (T_B+T_L)/2;
T(nx,ny)= (T_B+T_R)/2;
Told=T;
Tintial=T;
k1=alpha*dt/dx^2;
k2=alpha*dt/dy^2;
CFL=k1+k2; %to verify stability of the sol
iter_1=1;
tic
for m=1:1400
%convergence loop
while(error>tol)
for i=2:nx-1 % loop for x direction
for j=2:ny-1 % loop for Y dircetion
t1=Told(i-1,j)-2*Told(i,j)+Told(i+1,j);
t2=Told(i,j-1)-2*Told(i,j)+Told(i,j+1);
T(i,j)=Told(i,j)+k1*t1+k2*t2;
end
end
error= max(max(abs(Told-T)));
Told=T;
iter_1=iter_1+1;
time=toc;
end
end
[x,y]= contourf(xx,yy,T);
colormap(jet)
colorbar
clabel(x,y)
title_text=sprintf('No.of itertion=%d',iter_1);
title(title_text);
%pause(0.001)Graphs
`CFL=0.404
For the explicty problem we dont use the same time steps values which yeids to higher value of CFL number, which leads to unstable and uncovergence
Conclusion
1) The computational time and iteration required are SOR<GS<JACOBI for both steady and unsteady states.
2) the Explicit method, low CFL number needs to be considered for the covergence solutions.