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

Objective

To solve the 2D conduction equations using Transient solver and Steady-state solver using iterative techniques Jacobi, Gauss seidal and SOR.

Assumptions made

• The 2D model is assumed to be a perfect square
• nx=ny
• The initial temperature is assumed to be 298K(@ normal atmospheric temperature)
• Error=9e9

Boundary conditions:

Left: 400K

Right: 800K

Top: 600K

Bottom: 900K

Tolerance= 1e-4

Unsteady without thermal energy generation

Steady-state analysis

In the steady-state analysis, the rate of change of temperature with respect to time is negligible

In the above problem, both k1 and k2 are equal since nx=ny

here nx=20

Matlab code

Steady State Second-Order 2D PDE

clear all

close all

clc

%inputs

nx=20;

ny=nx;

%meshing

x=linspace(0,1,nx);

y=linspace(0,1,ny);

y1=fliplr(y);

dx=x(2)-x(1);

dy=y(2)-y(1);

%Error and Tolerance

error=9e9;

tolerance=1e-4;

%defining constants

k1=dy^2/(2*(dx^2+dy^2))

k2=dx^2/(2*(dx^2+dy^2))

% Boundary condition

T=298*ones(nx,ny); % normal atmospheric temperture

T(1,:)=600;

T(20,:)=900;

T(:,1)=400;

T(:,20)=800;

T_old=T;

%gathering input from user

iterative_solver=input(\'Choose a method to solve \"Steady state second order pde\" : 1) Jacobi Method 2)Gauss Seidal 3)SOR(gauss)\');

if iterative_solver==1

jacobi(x,y,dx,dy,T,T_old,nx,ny,k1,k2,error,tolerance,y1) % Jacobi method steady state

elseif iterative_solver==2

gauss_seidal(x,y,dx,dy,T,T_old,nx,ny,k1,k2,error,tolerance,y1) % Gauss seidal steady state

elseif iterative_solver==3

sor(x,y,dx,dy,T,T_old,nx,ny,k1,k2,error,tolerance,y1)%SOR steady state

else fprintf(\'Number is not valid\');

end

Jacobi Method

function jacobi(x,y,dx,dy,T,T_old,nx,ny,k1,k2,error,tolerance,y1) %FUNCTION 1

[X,Y]=meshgrid(x,y1)

count=1

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

T(i,j)=(k1*(T_old(i-1,j)+T_old(i+1,j)))+(k2*(T_old(i,j+1)+T_old(i,j-1))); %standard 5-point formula

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1;

end

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h)

text=sprintf(\'Jacobian Method: %d\',count)

title(text);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.03);

end

Gauss_seidal

function gauss_seidal(x,y,dx,dy,T,T_old,nx,ny,k1,k2,error,tolerance,y1) %FUNCTION 2

[X,Y]=meshgrid(x,y1)

count=1

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

T(i,j)=k1*(T(i-1,j)+T(i+1,j))+k2*(T(i,j+1)+T(i,j-1));

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1;

end

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h)

text=sprintf(\'Jacobian Method: %d\',count);

title(text);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.03);

end

SOR

function sor(x,y,dx,dy,T,T_old,nx,ny,k1,k2,error,tolerance,y1) %FUNCTION 3

[X,Y]=meshgrid(x,y1)

count=1

a=1.5;

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

Tsor(i,j)=k1*(T(i-1,j)+T(i+1,j))+k2*(T(i,j+1)+T(i,j-1));

T(i,j)=(T_old(i,j)*(1-a))+(a*Tsor(i,j));

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1;

end

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h)

text=sprintf(\'SOR: %d\',count)

title(text);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.03);

end

Output

From the above image, we can conclude that the iteration number for Jacobi method is more when compared to the Gauss seidal and SOR method.

Transient-state analysis

here nx is assumed to be 80

dt=1e-2

t=200

nt=t/dt

Matlab code

Transient state Second Order 2D PDE

Explicit

clear all

close all

clc

%inputs

nx=80;

ny=nx;

%meshing

x=linspace(0,1,nx);

y=linspace(0,1,ny);

y1=fliplr(y);

dx=x(2)-x(1);

dy=y(2)-y(1);

dt=1e-2;

% Time

t=200

nt= t/dt;

alpha=1e-4;

%Error and Tolerance

error=9e9;

tolerance=1e-4;

%defining constants

k3=(dt*alpha)/(dx^2);

k4=(dt*alpha)/(dy^2);

% Boundary condition

T=298*ones(nx,ny); % normal atmospheric temperture

T(1,:)=600;

T(80,:)=900;

T(:,1)=400;

T(:,80)=800;

T_old=T;

iterative_solver2=input(\'1)Jacobi 2) Gauss Seidal 3)SOR\')

if iterative_solver2==1

jacobi_Transient_Explicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k3,k4)

elseif iterative_solver2==2

gauss_seidal_Transient_Explicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k3,k4)

elseif iterative_solver2==3

SOR_Transient_Explicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k3,k4)

else sprintf(\'Number is not valid\')

end

Jacobi Transient Explicit method

function jacobi_Transient_Explicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k3,k4) %FUNCTION 4

[X,Y]=meshgrid(x,y1);

tic

count=1;

for k=1:nt

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

T(i,j)=T_old(i,j)+k3*(T_old(i-1,j)-(2*T_old(i,j))+T_old(i+1,j))+k4*(T_old(i,j+1)-(2*T_old(i,j))+T_old(i,j-1)); %standard 5-point formula

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1

end

end

time=toc;

output=sprintf(\'jacobi explicit transient method=%d\',time,count);

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h);

title(output);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.03)

Gauss seidal transient method

function gauss_seidal_Transient_Explicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k3,k4) %FUNCTION 5

[X,Y]=meshgrid(x,y1);

tic

count=1;

for k=1:nt

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

T(i,j)=T_old(i,j)+k3*(T(i-1,j)-(2*T(i,j))+T(i+1,j))+k4*(T(i,j+1)-(2*T(i,j))+T(i,j-1)); %standard 5-point formula

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1

end

end

time=toc;

output=sprintf(\'Gauss seidal explicit transient method=%d\',time,count);

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h);

title(output);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.03)

SOR Transient explicit

function SOR_Transient_Explicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k3,k4) %FUNCTION 6

[X,Y]=meshgrid(x,y1);

a=1.4;

tic

count=1;

for k=1:nt

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

T(i,j)=a*(T_old(i,j)+k3*(T(i-1,j)-(2*T(i,j))+T(i+1,j))+k4*(T(i,j+1)-(2*T(i,j))+T(i,j-1)))+((1-a)*T_old(i,j)); %standard 5-point formula

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1;

end

end

time=toc;

output=sprintf(\'Sor explicit transient method=%d\',time,count);

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h);

title(output);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.003)

Output

In the transient method, the temperature variation depends upon the time hence Time marching is carried on to solve this equation. In the explicit method, previous values are taken as ie nth value is taken as an input

Implicit method

Matlab code

Implicit

clear all

close all

clc

%inputs

nx=80;

ny=nx;

%meshing

x=linspace(0,1,nx);

y=linspace(0,1,ny);

y1=fliplr(y);

dx=x(2)-x(1);

dy=y(2)-y(1);

dt=1e-2;

% Time

t=200

nt= t/dt;

alpha=1e-4;

%Error and Tolerance

error=9e9;

tolerance=1e-4;

%defining constants

k1=(dt*alpha)/(dx^2);

k2=(dt*alpha)/(dy^2);

term1= (1+(2*k1)+(2*k2))^(-1);

term2= k1*term1;

term3=k2*term1

% Boundary condition

T=298*ones(nx,ny); % normal atmospheric temperture

T(1,:)=600;

T(80,:)=900;

T(:,1)=400;

T(:,80)=800;

T_old=T;

iterative_solver2=input(\'1)Jacobi 2) Gauss Seidal 3)SOR\')

if iterative_solver2==1

jacobi_Transient_Implicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k1,k2,term1,term2,term3)

elseif iterative_solver2==2

gauss_seidal_Transient_Implicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k1,k2,term1,term2,term3)

elseif iterative_solver2==3

SOR_Transient_Implicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k1,k2,term1,term2,term3)

else sprintf(\'Number is not valid\')

end

Jacobi Transient Implicit method

function jacobi_Transient_Implicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k1,k2,term1,term2,term3) %FUNCTION 4

[X,Y]=meshgrid(x,y1);

tic

count=1;

for k=1:nt

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

H= T_old(i-1,j)+T_old(i+1,j);

V= T_old(i,j-1)+T_old(i,j+1);

T(i,j)=(T_old(i,j)*term1)+(H*term2)+(V*term3); %standard 5-point formula

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1

end

end

time=toc;

output=sprintf(\'jacobi implicit transient method=%d\',time,count);

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h);

title(output);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.3)

Gauss seidal Transient Implicit method

function gauss_seidal_Transient_Implicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k1,k2,term1,term2,term3) %FUNCTION 5

[X,Y]=meshgrid(x,y1);

tic

count=1;

for k=1:nt

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

H= T(i-1,j)+T(i+1,j);

V= T(i,j-1)+T(i,j+1);

T(i,j)=(T_old(i,j)*term1)+(H*term2)+(V*term3); %standard 5-point formula

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1

end

end

time=toc;

output=sprintf(\'Gauss seidal implicit transient method=%d\',time,count);

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h);

title(output);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.03)

Sor Transient Implicit method

function SOR_Transient_Implicit(x,y,dx,dy,T,T_old,nx,ny,error,tolerance,y1,t,nt,alpha,k1,k2,term1,term2,term3) %FUNCTION 6

[X,Y]=meshgrid(x,y1);

a=1.4;

tic

count=1;

for k=1:nt

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

H= T(i-1,j)+T(i+1,j);

V= T(i,j-1)+T(i,j+1);

Tgs(i,j)=(T_old(i,j)*term1)+(H*term2)+(V*term3);

T(i,j)=((1-a)*T_old(i,j))+(a*Tgs(i,j));%standard 5-point formula

end

end

error=max(abs(T_old(i,j)-T(i,j)));

T_old=T;

count=count+1

end

end

time=toc;

output=sprintf(\'Gauss seidal implicit transient method=%d\',time,count);

[c,h]=contourf(X,Y,T);

colorbar;

colormap(jet);

clabel(c,h);

title(output);

xlabel(\'xaxis\');

ylabel(\'yaxis\');

pause(0.003)

Output

In the implicit method, values are taken at n+1th value as an input

Conclusion

From the above methods either Steady or Transient, we can conclude that the number of iterations is less for SOR then followed by Gauss seidal then Jacobi method.