%Unsteady-Implicit-Gauss_seidel
Tuioj = To;%Initialising Tuio as To for iterations
Tuij = To;%Initialising Tui as To for iterations
Tuipj=To;%Initialising Tuip as To for iterations
k1=CFLx;
k2 = CFLy;
term=1/(1+(2*k1)+(2*k2));
ecnj = 9e9;% Reinitializing error for unsteady implicit
for nt =1:500
while ecnj>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
li1=(Tuioj(p,m-1)+Tuioj(p,m+1));%Temperture gradient along X direction
li2 = (Tuioj(p-1,m)+Tuioj(p+1,m));%Temperture gradient along Y direction
Tuij(p,m)=(term*Tuipj(p,m))+(k2*term*li1)+(k1*term*li2);%Temperature at the node of interest
end
end
ecnj = max(abs(Tuioj-Tuij));%Calculating the error between the initial value and calculated value
Tuioj=Tuij;% Updating old iteration value
end
Tuipj=Tuij;% Updating previous time step value
ecnj=9e9; % Reinitializing error for next time loop
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tuij);
title('Unsteady-Implicit-Jacobi');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Unsteady-Implicit-Gauss_seidel
Tuios = To;%Initialising Tuio as To for iterations
Tuis = To;%Initialising Tui as To for iterations
Tuips=To;%Initialising Tuip as To for iterations
k1=CFLx;
k2 = CFLy;
term=1/(1+(2*k1)+(2*k2));
ecns = 9e9;% Reinitializing error for unsteady implicit
for nt =1:500
while ecns>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
li1=(Tuis(p,m-1)+Tuis(p,m+1));%Temperture gradient along X direction
li2 =(Tuis(p-1,m)+Tuis(p+1,m));%Temperture gradient along Y direction
Tuis(p,m)=(term*Tuips(p,m))+term*(1-w)*Tuios(p,m)+w*((k2*term*li1)+(k1*term*li2));%Temperature at the node of interest
end
end
ecns = max(abs(Tuios-Tuis));%Calculating the error between the initial value and calculated value
Tuios=Tuis;% Updating old iteration value
end
Tuips=Tuis;% Updating previous time step value
ecns=9e9; % Reinitializing error for next time loop
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tuis);
title('Unsteady-Implicit-SOR');
xlabel('X-direction,m');
ylabel('Y-direction,m');
clear all
close all
clc
%Initiating physical domain and initial input value
lx = 10; %Length of the physical domain m
ly = 10;%breadth of the physical domain m
dx = 1 ;%element size in x-direction m
dy = 1;%element size in y-direction m
nx = lx/dx; % number of elements in X-direction
ny = ly/dy; % number of elements in X-direction
a =1.4;%
dt = 0.1;%time step
CFLx = (a*dt)/(dx^2); %Courant number in x
CFLy = (a*dt)/(dy^2);%Courant number in y
%Initialising temperature field Ti = 27+273 = 300 K (Ambient temperature)
Ti = 300*ones(nx,ny);
%Initialising boundary condition in temperature field Tf
Tf =Ti;
%Initialising boundary condition in temperature field Tf
for i = 1:nx
if i == 1
for j = 1:ny
Ti(i,j) = 900;
end
else if i==nx
for j = 1:ny
Ti(i,j) = 600;
end
end
end
end
for j = 1:ny
if j == 1
for i = 1:nx
Ti(i,j) = 400;
end
else if j==ny
for i = 1:nx
Ti(i,j) = 800;
end
end
end
end
Ti(1,1) = (900+400)/2;
Ti(1,ny) = (800+900)/2;
Ti(nx,1) = (400+600)/2;
Ti(nx,ny) = (800+600)/2;
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Ti);
title('Initial condition');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Defining the error criteria ec
To= Ti;
Tn =Ti;
ecn = 9e9*ones(nx,ny);%Define ecn for initialising while loop
ec = 1e-4*ones(nx,ny);% Tolerance below which the the result converges
k=2*((dx^2)+(dy^2)/((dx^2)*(dy^2)));%Calculte term k for constants
%Jacobi
Toj = To;%Initialising Toj as To for iterations
Tj =To;%Initialising Tj as To for iterations
while max(ecn)>=ec(1,1) % Convergence criteria
for l=2:nx-1
for m=2:ny-1
lj1=(Toj(l,m-1)+Toj(l,m+1))/(dx^2);%Temperture gradient along X direction
lj2 = (Toj(l-1,m)+Toj(l+1,m))/(dy^2);%Temperture gradient along Y direction
Tj(l,m)=(1/k)*(lj1+lj2);%Temperature at the node of interest
end
end
ecn = max(abs(Toj-Tj)); %Calculating the error between the initial value and calculated value
Toj =Tj;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tj);
title('Jacobi');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Guass Seidel
Tog =To;%Initialising Tog as To for iterations
Tg = To;%Initialising Tg as To for iterations
ecn = 9e9; % Reinitializing error for gs
while max(ecn)>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
lg1=(Tg(p,m-1)+Tg(p,m+1))/(dx^2);%Temperture gradient along X direction
lg2 = (Tg(p-1,m)+Tg(p+1,m))/(dy^2);%Temperture gradient along Y direction
Tg(p,m)=(1/k)*(lg1+lg2);%Temperature at the node of interest
end
end
ecn = max(abs(Tog-Tg));%Calculating the error between the initial value and calculated value
Tog=Tg;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tg);
title('Gauss Seidel');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%SOR
Tos =To;%Initialising Tos as To for iterations
Ts = To;%Initialising Ts as To for iterations
w = 0.6;%Initialising a weighted error for SOR
ecn = 9e9; % Reinitializing error for sor
while max(ecn)>=ec(1,1)% Convergence criteria
for l=2:nx-1
for m=2:ny-1
ls1=(Tos(l,m-1)+Tos(l,m+1))/(dx^2);%Temperture gradient along X direction
ls2 = (Tos(l-1,m)+Tos(l+1,m))/(dy^2);%Temperture gradient along Y direction
Ts(l,m)=(1-w)*Tos(l,m)+w*((1/k)*(ls1+ls2));%Temperature at the node of interest
end
end
ecn = max(abs(Tos-Ts));%Calculating the error between the initial value and calculated value
Tos =Ts;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Ts);
title('SOR');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Unsteady
Tot = To;%Initialising Tot as To for iterations
Tt = To;%Initialising Tt as To for iterations
Tp=To;%Initialising Tp as To for retaining the initial field values
k1=CFLx;
k2 = CFLy;
for nt =1:1500
for i=2:nx-1
for j = 2:ny-1
t1 = (Tp(i+1,j)-2*Tp(i,j)+Tp(i-1,j));%Temperture gradient along X direction
t2 = (Tp(i,j+1)-2*Tp(i,j)+Tp(i,j-1));%Temperture gradient along Y direction
Tn(i,j)= Tp(i,j)+k1*t1+k2*t2;%Temperature at the node of interest
end
end
Tp=Tn;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tn);
title('Unsteady-Explicit');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Unsteady-Implicit-Gauss_seidel
Tuiog = To;%Initialising Tuio as To for iterations
Tuig = To;%Initialising Tui as To for iterations
Tuipg=To;%Initialising Tuip as To for iterations
k1=CFLx;
k2 = CFLy;
term=1/(1+(2*k1)+(2*k2));
ecng = 9e9;% Reinitializing error for unsteady implicit
for nt =1:500
while ecng>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
li1=(Tuig(p,m-1)+Tuig(p,m+1));%Temperture gradient along X direction
li2 = (Tuig(p-1,m)+Tuig(p+1,m));%Temperture gradient along Y direction
Tuig(p,m)=(term*Tuipg(p,m))+(k2*term*li1)+(k1*term*li2);%Temperature at the node of interest
end
end
ecng = max(abs(Tuiog-Tuig));%Calculating the error between the initial value and calculated value
Tuiog=Tuig;% Updating old iteration value
end
Tuipg=Tuig;% Updating previous time step value
ecng=9e9; % Reinitializing error for next time loop
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tuig);
title('Unsteady-Implicit-Gauss-Siedel');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Unsteady-Implicit-Gauss_seidel
Tuioj = To;%Initialising Tuio as To for iterations
Tuij = To;%Initialising Tui as To for iterations
Tuipj=To;%Initialising Tuip as To for iterations
k1=CFLx;
k2 = CFLy;
term=1/(1+(2*k1)+(2*k2));
ecnj = 9e9;% Reinitializing error for unsteady implicit
for nt =1:500
while ecnj>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
li1=(Tuioj(p,m-1)+Tuioj(p,m+1));%Temperture gradient along X direction
li2 = (Tuioj(p-1,m)+Tuioj(p+1,m));%Temperture gradient along Y direction
Tuij(p,m)=(term*Tuipj(p,m))+(k2*term*li1)+(k1*term*li2);%Temperature at the node of interest
end
end
ecnj = max(abs(Tuioj-Tuij));%Calculating the error between the initial value and calculated value
Tuioj=Tuij;% Updating old iteration value
end
Tuipj=Tuij;% Updating previous time step value
ecnj=9e9; % Reinitializing error for next time loop
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tuij);
title('Unsteady-Implicit-Jacobi');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Unsteady-Implicit-Gauss_seidel
Tuios = To;%Initialising Tuio as To for iterations
Tuis = To;%Initialising Tui as To for iterations
Tuips=To;%Initialising Tuip as To for iterations
k1=CFLx;
k2 = CFLy;
term=1/(1+(2*k1)+(2*k2));
ecns = 9e9;% Reinitializing error for unsteady implicit
for nt =1:500
while ecns>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
li1=(Tuis(p,m-1)+Tuis(p,m+1));%Temperture gradient along X direction
li2 =(Tuis(p-1,m)+Tuis(p+1,m));%Temperture gradient along Y direction
Tuis(p,m)=(term*Tuips(p,m))+term*(1-w)*Tuios(p,m)+w*((k2*term*li1)+(k1*term*li2));%Temperature at the node of interest
end
end
ecns = max(abs(Tuios-Tuis));%Calculating the error between the initial value and calculated value
Tuios=Tuis;% Updating old iteration value
end
Tuips=Tuis;% Updating previous time step value
ecns=9e9; % Reinitializing error for next time loop
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tuis);
title('Unsteady-Implicit-SOR');
xlabel('X-direction,m');
ylabel('Y-direction,m');
clear all
close all
clc
%Initiating physical domain and initial input value
lx = 10; %Length of the physical domain m
ly = 10;%breadth of the physical domain m
dx = 1 ;%element size in x-direction m
dy = 1;%element size in y-direction m
nx = lx/dx; % number of elements in X-direction
ny = ly/dy; % number of elements in X-direction
a =1.4;%
dt = 0.1;%time step
CFLx = (a*dt)/(dx^2); %Courant number in x
CFLy = (a*dt)/(dy^2);%Courant number in y
%Initialising temperature field Ti = 27+273 = 300 K (Ambient temperature)
Ti = 300*ones(nx,ny);
%Initialising boundary condition in temperature field Tf
Tf =Ti;
%Initialising boundary condition in temperature field Tf
for i = 1:nx
if i == 1
for j = 1:ny
Ti(i,j) = 900;
end
else if i==nx
for j = 1:ny
Ti(i,j) = 600;
end
end
end
end
for j = 1:ny
if j == 1
for i = 1:nx
Ti(i,j) = 400;
end
else if j==ny
for i = 1:nx
Ti(i,j) = 800;
end
end
end
end
Ti(1,1) = (900+400)/2;
Ti(1,ny) = (800+900)/2;
Ti(nx,1) = (400+600)/2;
Ti(nx,ny) = (800+600)/2;
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Ti);
title('Initial condition');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Defining the error criteria ec
To= Ti;
Tn =Ti;
ecn = 9e9*ones(nx,ny);%Define ecn for initialising while loop
ec = 1e-4*ones(nx,ny);% Tolerance below which the the result converges
k=2*((dx^2)+(dy^2)/((dx^2)*(dy^2)));%Calculte term k for constants
%Jacobi
Toj = To;%Initialising Toj as To for iterations
Tj =To;%Initialising Tj as To for iterations
while max(ecn)>=ec(1,1) % Convergence criteria
for l=2:nx-1
for m=2:ny-1
lj1=(Toj(l,m-1)+Toj(l,m+1))/(dx^2);%Temperture gradient along X direction
lj2 = (Toj(l-1,m)+Toj(l+1,m))/(dy^2);%Temperture gradient along Y direction
Tj(l,m)=(1/k)*(lj1+lj2);%Temperature at the node of interest
end
end
ecn = max(abs(Toj-Tj)); %Calculating the error between the initial value and calculated value
Toj =Tj;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tj);
title('Jacobi');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Guass Seidel
Tog =To;%Initialising Tog as To for iterations
Tg = To;%Initialising Tg as To for iterations
ecn = 9e9; % Reinitializing error for gs
while max(ecn)>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
lg1=(Tg(p,m-1)+Tg(p,m+1))/(dx^2);%Temperture gradient along X direction
lg2 = (Tg(p-1,m)+Tg(p+1,m))/(dy^2);%Temperture gradient along Y direction
Tg(p,m)=(1/k)*(lg1+lg2);%Temperature at the node of interest
end
end
ecn = max(abs(Tog-Tg));%Calculating the error between the initial value and calculated value
Tog=Tg;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tg);
title('Gauss Seidel');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%SOR
Tos =To;%Initialising Tos as To for iterations
Ts = To;%Initialising Ts as To for iterations
w = 0.6;%Initialising a weighted error for SOR
ecn = 9e9; % Reinitializing error for sor
while max(ecn)>=ec(1,1)% Convergence criteria
for l=2:nx-1
for m=2:ny-1
ls1=(Tos(l,m-1)+Tos(l,m+1))/(dx^2);%Temperture gradient along X direction
ls2 = (Tos(l-1,m)+Tos(l+1,m))/(dy^2);%Temperture gradient along Y direction
Ts(l,m)=(1-w)*Tos(l,m)+w*((1/k)*(ls1+ls2));%Temperature at the node of interest
end
end
ecn = max(abs(Tos-Ts));%Calculating the error between the initial value and calculated value
Tos =Ts;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Ts);
title('SOR');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Unsteady
Tot = To;%Initialising Tot as To for iterations
Tt = To;%Initialising Tt as To for iterations
Tp=To;%Initialising Tp as To for retaining the initial field values
k1=CFLx;
k2 = CFLy;
for nt =1:1500
for i=2:nx-1
for j = 2:ny-1
t1 = (Tp(i+1,j)-2*Tp(i,j)+Tp(i-1,j));%Temperture gradient along X direction
t2 = (Tp(i,j+1)-2*Tp(i,j)+Tp(i,j-1));%Temperture gradient along Y direction
Tn(i,j)= Tp(i,j)+k1*t1+k2*t2;%Temperature at the node of interest
end
end
Tp=Tn;% Updating old iteration value
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tn);
title('Unsteady-Explicit');
xlabel('X-direction,m');
ylabel('Y-direction,m');
%Unsteady-Implicit-Gauss_seidel
Tuio = To;%Initialising Tuio as To for iterations
Tui = To;%Initialising Tui as To for iterations
Tuip=To;%Initialising Tuip as To for iterations
k1=CFLx;
k2 = CFLy;
term=1/(1+(2*k1)+(2*k2));
ecn = 9e9;% Reinitializing error for unsteady implicit
for nt =1:500
while ecn>=ec(1,1)% Convergence criteria
for p=2:nx-1
for m=2:ny-1
li1=(Tui(p,m-1)+Tui(p,m+1));%Temperture gradient along X direction
li2 = (Tui(p-1,m)+Tui(p+1,m));%Temperture gradient along Y direction
Tui(p,m)=(term*Tuip(p,m))+(k2*term*li1)+(k1*term*li2);%Temperature at the node of interest
end
end
ecn = max(abs(Tuio-Tui));%Calculating the error between the initial value and calculated value
Tuio=Tui;% Updating old iteration value
end
Tuip=Tui;% Updating previous time step value
ecn=9e9; % Reinitializing error for next time loop
end
%plotting a graph for Temperature gradient in X-Y field
figure
x=linspace(1,lx,nx);
y=linspace(1,ly,ny);
[X Y] = meshgrid(x,y);
contourf(X,Y,Tui);
title('Unsteady-Implicit');
xlabel('X-direction,m');
ylabel('Y-direction,m');