BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

VVenkatavamsikrishna Mahavratayajula1 Nov 2020
%% Steady state Jacobian iterative technique

clear all
close all
clc
% inputs
L=1 % Length of domain (square);
%nx,ny= no of grid points;
nx=10;
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;
T=300*ones(nx,ny);
% Define boundary conditions
T_t =600;
T_b =900;
T_l =400;
T_r =800;
%Edge refinement
T(1,:)=T_t
T(end,:)=T_b
T(:,1)=T_l
T(:,end)=T_r
Told=T;
iter = 0;
% Jacobi method
##Tp=(T1+Tr/dx^2*k)+(Tt+Tb/dx^2*k);
%Convergence loop
while error>tol
% space loop
for i=2:nx-1
for j=2:ny-1
T1 = Told(i-1,j)+Told(i+1,j); 
T2 = Told(i,j-1)+Told(i,j+1);
T(i,j)= 0.25*(T1+T2);
end
end
error=max(max(abs(Told-T)));
Told=T;
iter = iter+1;
end
figure(1)
contourf(x,fliplr(y),T,'ShowText','ON');
colorbar;
colormap(jet)
disp(iter);


%% Steady state Gauss siedel iterative technique
clear all
close all
clc
% inputs
L=1  % Length of domain (square);
%nx,ny= no of grid points;
nx=10;
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;
T=300*ones(nx,ny);
% Define boundary conditions
T_t =600;
T_b =900;
T_l =400;
T_r =800;
%Edge refinement
T(1,:)=T_t
T(end,:)=T_b
T(:,1)=T_l
T(:,end)=T_r
Told=T;
Tintial=T
iter = 0;
% Gauss seisel
%Convergence loop
while error>tol
% space loop
for i=2:nx-1
 for j=2:ny-1
T1 = T(i-1,j)+Told(i+1,j);   
T2 = T(i,j-1)+Told(i,j+1);
T(i,j)= 0.25*(T1+T2);
end
end
error=max(max(abs(Told-T)));
Told=T;
iter = iter+1;
end
figure(1)
contourf(x,fliplr(y),T,'ShowText','ON');
colorbar;
colormap(jet)
disp(iter);


## Steady state s.o.r method ##
clear all
close all
clc
% inputs
L=1  % Length of domain (square);
%nx,ny= no of grid points;
nx=10;
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;
T=300*ones(nx,ny);
% Define boundary conditions
T_t =600;
T_b =900;
T_l =400;
T_r =800
%Edge refinement
T(1,:)= T_t
T(end,:)= T_b
T(:,1)= T_l
T(:,end)= T_r
% Relaxation factor
rf=0.1;
% S.O.r
Told = T;
Tintial=Told;
iter=0;
%Convergence loop
while error>tol
% space loop
for i=2:nx-1
for j=2:ny-1
T1 = Told(i-1,j)+Told(i+1,j);   
T2 = Told(i,j-1)+Told(i,j+1);
T(i,j)= ((1-rf)*Told(i,j)+(rf*0.25*(T1+T2)));
end
end
error=max(max(abs(Told-T)));
Told=T;
iter=iter+1;
end
figure(1)
contourf(x,fliplr(y),T,'ShowText','ON');
colorbar;
colormap(jet)
disp(iter);

%%Unsteady Explicit 

clear all
close all
clc
% inputs
L=1  % Length of domain (square);
%nx,ny= no of grid points;
nx=10;
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.4;
dt=1e-3;
T=300*ones(nx,ny);
% Define boundary conditions
T_t =600;
T_b =900;
T_l =400;
T_r =800;
%Edge refinement
T(1,:)=T_t
T(end,:)=T_b
T(:,1)=T_l
T(:,end)=T_r
Tprev=T;
iter=0;
% constants
K1=(alpha*dt)/(dx^2);
K2=(alpha*dt)/(dy^2);
cfl = K1+K2;
%time loop
while error>tol
% space loop
for i=2:nx-1
 for j=2:ny-1
term2 = K1*((Tprev(i-1,j))+Tprev(i+1,j));   
term3 = K2*(Tprev(i,j-1)+Tprev(i,j+1));
T(i,j)= Tprev(i,j)*(1-2*K1-2*K2)+term2+term3;
end
end
iter=iter+1;
error=max(max(abs(Tprev-T)));
Tprev=T;
end
figure(1)
contourf(x,fliplr(y),T,'ShowText','ON');
colorbar;
colormap(jet)
disp(iter);

%% unsteady implicit jacobian

clear all
close all
clc
% inputs
% unsteady 2D heat conduction
% by implicit method
L=1 % Length of domain (square);
nx=10;
ny=nx;
rf=1.5;
dt=0.1;
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.4;
nt=1000;
% Define boundary conditions
T=300*ones(nx,ny);
T_t =600;
T_b =900;
T_l =400;
T_r =800;
%Edge refinement
T(1,:)=T_t
T(end,:)=T_b
T(:,1)=T_l
T(:,end)=T_r
Told=T;
Tprev=Told;
titer=0;
siter = titer;
% constants
K1=(alpha*dt)/(dx^2);
K2=(alpha*dt)/(dy^2);
cfl = K1+K2;
% Jacobian
%Time loop
for k=1:nt
%Convergence loop
error = 9e9;
while error>tol
% space loop
for i=2:nx-1
for j=2:ny -1
term2 = K1*(Told(i-1,j)+Told(i+1,j));
term3 = K2*(Told(i,j-1)+Told(i,j+1));
T(i,j)= (Tprev(i,j)+term2+term3)/(1+2*(K1+K2));
end
end
error=max(max(abs(Told-T)));
Told=T;
siter=siter+1;
end
conv = max(max(abs(Tprev-T)));
if conv<tol
break;
end
Tprev=T;
titer=titer+1;
end

figure(1)
contourf(x,fliplr(y),T,'ShowText','ON');
colorbar;
colormap(jet)
pause(0.3);
disp(siter);
disp(titer);

%% Unsteady implicit Gauss seidel

clear all
close all
clc
% inputs
% unsteady 2D heat conduction
% by implicit method
L=1 % Length of domain (square);
nx=10;
ny=nx;
rf=0.1;
dt=0.1;
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.4;
nt=1000;
% Define boundary conditions
T=300*ones(nx,ny);
T_t =600;
T_b =900;
T_l =400;
T_r =800;
%Edge refinement
T(1,:)=T_t
T(end,:)=T_b
T(:,1)=T_l
T(:,end)=T_r
Told=T
Tprev=T
titer=0;
siter = titer;
% constants
K1=(alpha*dt)/(dx^2);
K2=(alpha*dt)/(dy^2);
cfl = K1+K2;
% Gauss seidel
%Convergence loop
for k=1:nt
error = 9e9;
while error>tol
% space loop
for i=2:nx-1
for j=2:ny-1
term2 = K1*((Told(i-1,j))+(T(i+1,j)));
term3 = K2*((Told(i,j-1))+(T(i,j+1)));
T(i,j)= (Tprev(i,j)+term2+term3)/(1+2*(K1+K2));
end
end
error=max(max(abs(Told-T)));
Told=T;
siter=siter+1;
end
conv = max(max(abs(Tprev-T)));
if conv<tol
break;
end
Tprev=T;
titer=titer+1;
end
figure(1)
contourf(x,fliplr(y),T,'ShowText','ON');
colorbar;
colormap(jet)
pause(0.3);
disp(titer);
disp(siter);

%% Unsteady implicit S.O.R

clear all
close all
clc
% inputs
% unsteady 2D heat conduction
% by implicit method
L=1 % Length of domain (square);
nx=10;
ny=nx;
rf=0.1;
dt=0.1;
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.4;
nt=1000;
% Define boundary conditions
T=300*ones(nx,ny);
T_t =600;
T_b =900;
T_l =400;
T_r =800;
%Edge refinement
T(1,:)=T_t
T(end,:)=T_b
T(:,1)=T_l
T(:,end)=T_r
Tprev=T;
Told=T;
titer=0;
siter = titer;
% constants
K1=(alpha*dt)/(dx^2);
K2=(alpha*dt)/(dy^2);
cfl = K1+K2;
% S.O.I
%time loop
for k=1:nt
%Convergence loop
error = 9e9;
while error>tol
% space loop
for i=2:nx-1
for j=2:ny-1
term2 = K1*((Told(i-1,j))+(Told(i+1,j)));
term3 = K2*((Told(i,j-1))+(Told(i,j+1)));
T(i,j)= ((1-rf)* Told(i,j))+rf*((Tprev(i,j)+term2+term3)/(1+2*(K1+K2)));
end
end
error=max(max(abs(Told-T)));
Told=T;
siter=siter+1;
end
conv = max(max(abs(Tprev-T)));
if conv<tol
break;
end
Tprev=T;
titer=titer+1;
end
figure(1)
contourf(x,fliplr(y),T,'ShowText','ON');
colorbar;
colormap(jet)
pause(0.3);
disp(siter);
disp(titer);

##Explination

Steady state analysis

1)Change in properties with respect to time remain constant known as Steady state analysis

2) 2D steady state equation `(d^2t)/(dx^2) + (d^2t)/(dx^2) =0`

3)Assumed it as square domain (nx=ny)

4) Temperature at point given by (Tp)=`1/k((Tl+Tr) / (dx^2)) +1/k((Tt+Tb) / (dy^2))`

5) convergeloop-spaceloop-Update Told

6)converge loop (while error>tolerence=run) else break or end loop

Iterative techniques

1)Jacobian

2)Gauss seidel

3) suceesive over relaxation

1)Jacobian:

1)Temperature at a point(T(i,j))=`1/k((Tl(old)+Tr(old)) / (dx^2)) +1/k((Tt (old)+Tb (old)) / (dy^2))`

2)since Tl= T(i-1,j) :Tr = T(I+1,j) : Tt = T(j+1,i) : Tb = T(j-1,i)

3)Update Told values

4) error in jacobian= 0.000094868

5) No of iterations = 207

2)Gauss seidel:

1)Temperature at a point(T(i,j))=`1/k((Tl+Tr(old)) / (dx^2)) +1/k((Tt +Tb (old)) / (dy^2))`

2)since Tl= T(i-1,j) :Tr = T(I+1,j) : Tt = T(j+1,i) : Tb = T(j-1,i)

3)Update Told values

4) error in Gauss seidel= 0.000095865

5) No of iterations = 111

3)Steady S.O.R:

1)Temperature at a point(T(i,j))=`(1-omega)*Told+omega*(Jacobian or gauss seidel)

2)since Tl= T(i-1,j) :Tr = T(I+1,j) : Tt = T(j+1,i) : Tb = T(j-1,i)

3) where omega is relaxation factor omega<1 (under relaxation) omega>1 (over relaxation)

4)Update Told values

5) error in Gauss seidel= 0.000099972

6) No of iterations depends on omega values

No of iterations omega
1730

0.1

920

 0.2
634 0.3
487 0.4
396 0.5
334 0.6
289 0.7
255 0.8
228 0.9
207 1.0

Un-steady/Transiant state analysis :

1)Change in properties with respect to time are not constant/varying known as Unsteady state analysis

2) 2D Transiant state equation `α*((d^2t)/(dx^2) + (d^2t)/(dx^2)) =((dt)/(dτ))

4) Temperature at point given by `((Tp^(n+1)-Tp^n)/(dt)) = α*((Tl-2*(Tp)+Tr)/(dx^2) + (Tt-2*(Tp)+Tb)/(dy^2))^((n)or(n+1))`

5)Time loop- convergeloop-spaceloop-Update Told, Tprev

6)converge loop (while error>tolerence=run) else break or end loop

7)

Steady state analysis:

##Explination

Steady state analysis

7) EXplicit solution(n) is independent of other variables but implicit solution(n+1) is dependant on variable`

1)Unsteady explicit:

1)unsteady explicitT(i,j) given by

(Tp^(n+1)-Tp^n)/(dt)) = α*((Tl-2*(Tp)+Tr)/(dx^2) + (Tt-2*(Tp)+Tb)/(dy^2))^((n)

`Tp^(n+1)=Tp^n+K1*(Tl-2*Tp+Tr)^n+K1*(Tt-2*Tp+Tb)^n`

`Where K1=α*dt/dx^2` `K2=α*dt/dx^2`

Explict refers to diffusion problem given by CFL number

`CFl<=;=0.5 Stable CFl>=0.5 unstable & CFlx+CFly<=0.5`

2) since Tl= T(i-1,j) :Tr = T(I+1,j) : Tt = T(j+1,i) : Tb = T(j-1,i)

3)Update Told values

4) error in unsteady explicit= 0.000098836

5) No of iterations = 433

Implicit:

Expression for implicit

1)((Tp^(n+1)-Tp^n)/(dt)) = α*((Tl-2*(Tp)+Tr)/(dx^2) + (Tt-2*(Tp)+Tb)/(dy^2))^(n+1)

`Tp^(n+1) = 1/(1+2k1+2k2) [Tp^n+K1(Tl+Tr)^(n+1) +K2 (Tt+Tb)^(n+1)`

`K1=(α*dt)/(dx^2) K2=(α*dt)/(dy^2)`

`term1=1/(1+2k1+2k2`

Term2=K1*term1

Term3=K2*Term1

`Tp^(n+1)=Term1*Tp^n+Term2*(Tl+Tr)^(n+1)+Term3*(Tt+Tb)^n+1`

1)Jacobian:

1)Terem2=K1*( Tl(old)+Tr(old)))

2) Term3 = K2*(Tt (old)+Tb (old))

2)since Tl= T(i-1,j) :Tr = T(I+1,j) : Tt = T(j+1,i) : Tb = T(j-1,i)

T(i,j)=(Tprev(i,j)+term2+term3)/(1+2*(K1+K2))

3)Update Told,Tprev values

4) error in jacobian= 0.000094406

5) No of iterations = 840

2)Gauss seidel:

`Terem2=K1*( Tl(old)+Tr`

`Term3 = K2*(Tt +Tb (old)`

2)since Tl= T(i-1,j) :Tr = T(I+1,j) : Tt = T(j+1,i) : Tb = T(j-1,i)

T(i,j)=(Tprev(i,j)+term2+term3)/(1+2*(K1+K2))

3)Update Told ,Tprev values

4) error in Gauss seidel= 0.000094406

5) No of iterations = 840

3)Steady S.O.R:

1)term2 = K1*((Told(i-1,j))+(Told(i+1,j)))
2)term3 = K2*((Told(i,j-1))+(Told(i,j+1)))

3)Temperature at a point(T(i,j))=`(1-omega)*Told+omega*(Jacobian or gauss seidel)

4)since Tl= T(i-1,j) :Tr = T(I+1,j) : Tt = T(j+1,i) : Tb = T(j-1,i)

5)T(i,j)= ((1-rf)* Told(i,j))+rf*((Tprev(i,j)+term2+term3)/(1+2*(K1+K2)));

3) where omega is relaxation factor omega<1 (under relaxation) omega>1 (over relaxation)

4)Update Told,Tprev values

5) error in Gauss seidel= 0.000099287 (at rf=0.1)

6) No of iterations depends on omega values

No of iterations omega

iterations Rf
5935

0.1
3332 0.2
2369 0.3
1855 0.4
1531 0.5
1310 0.6
1147 0.7
1021 0.8
921 0.9
840 1.0
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.