Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
Solving Steady and Unsteady State 2d Heat Conduction Equation
The two dimension heat conduction equation can be written as
`((delT)/(delt))+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
In the above equation T represents Temperature, 'alpha' represents thermal diffusivity. x and y are spatial coordinates and t is the temporal coordinate.
There are two approaches to solve the problem. These are listed below.
1) Steady State Equation:-
In this approach the change of temperature with time is zero, i.e the term `((delT)/(delt))' is equal to zero. The equation reduces to
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
In order to discretize the above equation we create a domain of unit length in x and y direction resulting in a unit square. This domain is divided 50 equally spaced grid points in the x and y direction. The Temperature boundary conditions are applied as given below
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
A second order central differencing scheme is applied to calculate the second order derivative terms in the steady state equation. The value of Temperature at point P is evaluated using the neighbourhood values, i.e TL and TR in x direction, Tt and Tb in the y direction.
The discretised equation for 2D steady state is given below
`(T_L-2T_p+T_R)/(Deltax^2) +(T_T-2T_p+T_B)/(Deltay^2)=0`
Solving for 'T_p'
`T_p=((Deltax^2*Deltay^2)/(2(Deltax^2+Deltay^2)))*(((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)))`
The term '((2(Deltax^2+Deltay^2))/(Deltax^2*Deltay^2)))' is labelled as k
The equation Reduces to
`T_p=(1/K)*(((T_L+T_R)/(Deltax^2)+(T_T+T_B)/(Deltay^2)))`
There are three iterative solvers which have been used to solve the above equation
a) Jacobi Method
In this method, given a current approximation for x,
`x^k=(x_1^k,x_2^k,x_3^k,.......x_n^k)`
the strategy of Jacobi's Method is to use the current values of 'x_2^k','x_3^k','x_n^k' to find a new value for x_1^(k+1)
The Code written for iterative solver using Jacobi Method is given below
%Jacobi Method
clear all
close all
clc
%Input Parameter
nx=50;
ny=50;
x=linspace(0,1,nx);%Grid along X
y=linspace(0,1,ny);%Grid along y
dx=x(2)-x(1);
dy=y(2)-y(1);
T=ones(nx,ny);
%Boundary Conditions
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Tinitial=T;
Told=T;
k=2*((dx^2+dy^2)/(dx^2*dy^2));
Tol=1e-4
error1=9e9
error2=9e9
error3=9e9
jacobi_iter=1;
tic
a=1;
while (error1>Tol)
for i=2:nx-1
for j=2:ny-1
m=(Told(i+1,j)+Told(i-1,j))/(dx*dx);
n=(Told(i,j+1)+Told(i,j-1))/(dy*dy);
T(i,j)=(m+n)/k ;
end
end
error1= max(max(abs(Told-T)));
ERROR1(a)=error1; %Plotting error vs iterations
Told=T;
jacobi_iter=jacobi_iter+1;
JI(a)=jacobi_iter;%Storing iterations in counter
time1=toc;
a=a+1;
end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi solver for steady state',tit1})
colorbar
pause(0.3)-Plot for Temperature Field . For a grid spacing of 50x50, dx and dy value of 0.0204
No of iterations= 4873. Simulation time 2.189e-1seconds
b) Gauss-Seidel Method
After we compute 'x_1^(k+1)' and it is a better approximation to the true value of x_1 than x_1^k we use the new value 'x_1^(k+1)' in finding 'x_2^(k+1)' . We progressively use the updated values of the variables in the equation. The result is we get faster convergence and lesser simulation time.
Code for Gauss Seidel( TO be read in continuation of Jacobi code)
%Gauss Seidel
tic;
T=Tinitial;
Told=Tinitial;
gauss_seidel=1;
b=1;
while (error2>Tol)
for i=2:nx-1
for j=2:ny-1
m=(T(i+1,j)+T(i-1,j))/(dx*dx);
n=(T(i,j+1)+T(i,j-1))/(dy*dy);
T(i,j)=(m+n)/k ;
end
end
error2= max(max(abs(Told-T)));
Told=T;
ERROR2(b)= error2;
GS(b)=gauss_seidel;
gauss_seidel=gauss_seidel+1;
b=b+1;
time2=toc;
end
figure(2)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel solver for steady state',tit2})
colorbar
pause(0.3)
-Plot for Temperature Field . For a grid spacing of 50x50, dx and dy value of 0.0204
No of iterations= 2610. Simulation time 1.789e-1seconds
-Successive Over Relaxation Method (SOR)
The SOR method of solving a linear system of equations Ax=B is derived by extrapolating the Gauss Seidel Method. The extrapolation takes a the form of weighted average between previous iterate and the computed Gauss seidel iterate successively for each component
`x_i^k=omega*x_i^k+(1-omega)*x_i^(k-1)`
'omega' is the relaxation factor
'omega>1', the equation is overelaxed
'omega=1', The equation is same as Gauss Seidel
'omega<1', The equation is under relaxed
COde for SOR (TO be read in continuation of code for Gauss Seidel)
%Successive Over Relaxation
SOR=1;
T=Tinitial
Told=Tinitial;
omega=1.5;
tic
c=1;
while (error3>Tol)
for i=2:nx-1
for j=2:ny-1
m=(T(i+1,j)+T(i-1,j))/(dx*dx);
n=(T(i,j+1)+T(i,j-1))/(dy*dy);
T(i,j)=(1-omega)*Told(i,j) +omega*((m+n)/k );
end
end
error3= max(max(abs(Told-T)));
ERROR3(c)= error3;
Told=T;
SR(c)=SOR;
SOR=SOR+1;
time3=toc;
c=c+1;
end
figure(3)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit3=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR solver for steady state',tit3})
colorbar
pause(0.3)
Plot for Temperature Field . For a grid spacing of 50x50, dx and dy value of 0.0204, omega=1.5
No of iterations= 957. Simulation time 8.3477e-02seconds
- Plot showing variation of error with iterations for the three methods
By observing the below plot we note that SOR method converges at much faster rate compared to Gauss seidel and Jacobi.
Unsteady Transient Equation
The unsteady Transient equation for heat conduction is given as
`(delT)/(delt) =alpha*((del^2T)/(delx^2) +(del^2T)/(dely^2) )`
Discretizing this equation we get
`(T_p^(n+1)-T_p^n)/(Deltat)=alpha*((T_L-2*T_p+T_R)/(Deltax^2) +(T_T-2*T_p+T_B)/(Deltay^ 2))^n`
There are two methods to solve this equation
1) Explicit:- In the explicit method, there would be only one unknown. In this case it is `(T_p^(n+1)'. All other variables are at time step n which have been calculated.
Solving for `(T_p^(n+1)', we get
`T_p^(n+1)=T_p^(n)+ ((alpha*Deltat)/(Deltax^2))*(T_L^n-2*T_p^n+T_R^n)+ ((alpha*Deltat)/(Deltay^2))*(T_T^n-2*T_p^n+T_B^n)`
In this equation we have the one unknown `(T_p^(n+1)' on the LHS. All the parameters on RHS are at time step n , a known quantity.
In the above equation we can denote as follows.
' k1=((alpha*Deltat)/(Deltax^2))'
' k2=((alpha*Deltat)/(Deltay^2))'
'k1' and 'k2' are called the CFL numbers for the equation
For the numerical explicit diffusion equation to be stable, it must satisfy the below criteria
`k1+k2<= 0.5`
In order to solve the transient explicit equation, we solve for Temperature values at each node using values from previous time step. We then run the loop for sufficiently large time steps to obtain the solution. No iterative solver is required.
Code for transient equation
clear all
close all
clc
nx=50;
ny=50;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=4e-3;
dt=0.01;
k1=(alpha*dt)/(dx)^2;
k2=(alpha*dt)/(dy)^2;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Told=T;
T(1,1)=(400+900)/2;
T(end,1)=(400+600)/2;
T(1,end)=(900+800)/2;
T(end,end)=(800+600)/2;
% k=2*((dx^2+dy^2)/(dx^2*dy^2));
Tol=1e-4;
error1=9e9;
for nt=1:15000
tic;
for i=2:nx-1
for j=2:ny-1
m=k1*(Told(i+1,j)+Told(i-1,j)-2*Told(i,j));
n=k2*(Told(i,j+1)+Told(i,j-1)-2*Told(i,j));
T(i,j)=m+n+Told(i,j) ;
end
end
time1=toc;
Told=T;
end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Simulation Time=%d s',time1)
title({'Solution for transient explicit '})
colorbar
pause(0.3)Plot for temperature profile
Parameters considered
No of points along x=50
No of points along y=50
dx=dy=0.0204
thermal diffusivity 'alpha=4e-3'
dt=0.01
CFL along x 'k1=0.0960'
CFL along y 'k2=0.0960'
'k1+k2=0.192' which is less than 0.5. We expect a solution
No of time steps= 150000
Plot
Plot when 'k1+k2 >0.5'. The solution is unstable
Implicit Method
In this method, we have more than one unknown element
`(T_p^(n+1)-T_p^n)/(Deltat)=alpha*((T_L-2*T_p+T_R)/(Deltax^2) +(T_T-2*T_p+T_B)/(Deltay^ 2))^(n+1)`
`T_p^(n+1)=T_p^(n)+ ((alpha*Deltat)/(Deltax^2))*(T_L^(n+1)-2*T_p^(n+1)+T_R^(n+1))+ ((alpha*Deltat)/(Deltay^2))*(T_T^(n+1)-2*T_p^(n+1)+T_B^(n+1))`
In the above equation we can denote as follows.
' k1=((alpha*Deltat)/(Deltax^2))'
' k2=((alpha*Deltat)/(Deltay^2))'
Collecting terms related to `T_p^(n+1)' on LHS we get
`'T_p^(n+1)*[1+2k1+2k2]= T_p^n+(k1*(T_L^(n+1)+T_R^(n+1)) +k2*(T_T^(n+1)+T_B^(n+1)) )`
Solving for `'T_p^(n+1)'
`'T_p^(n+1)*= [1+2k1+2k2]^-1*(T_p^n+k1*(T_L^(n+1)+T_R^(n+1)) +k2*(T_T^(n+1)+T_B^(n+1)) )`
For ease of coding the terms of the above equation are labelled as
'[1+2k1+2k2]^-1= Term1'
'k1*[1+2k1+2k2]^-1= k1*Term1=Term2'
'k2*[1+2k1+2k2]^-1= k2*Term1=Term3'
The above equation is soved using the 3 iterative solvers that were used previously
1. Jacobi solver
clear all
close all
clc
nx=50;
ny=50;
x=linspace(0,1,nx);%Grid Spacing along X
y=linspace(0,1,ny);%Grid Spacing along Y
dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=1.5;%Thermal Diffusivity
dt=1e-2;%Time Step
k1=(alpha*dt)/(dx)^2;%CFL along x
k2=(alpha*dt)/(dy)^2;%CFL along y
term1=(1/(1+2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
T(1,1)=(400+600)/2;
T(end,1)=(400+800)/2;
T(1,end)=(900+800)/2;
T(end,end)=(800+600)/2;
Told=T;
Tinitial=T;
Tprev=T;
Tol=1e-4;
jacobi_iter=1;
a=1
for nt=1:1400
error1=9e9;
tic;
while(error1>Tol)
for i=2:nx-1
for j=2:ny-1
H=term2*(Told(i-1,j)+Told(i+1,j));
V=term3*(Told(i,j+1)+Told(i,j-1));
T(i,j)=term1*Tprev(i,j)+ H+V;
end
end
error1=max(max(abs(Told-T)));
Error1(a)=error1;
Told=T;%Spatial Temp values updated within the space iteration
jacobi_iter=jacobi_iter+1;
JI(a)=jacobi_iter;
a=a+1;
time1=toc;
end
Tprev=T;%After convergence the T values updated
end
figure(1)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi Iter solver for implicit',tit1})
colorbar
pause(0.3)Plot for Temperature
Parameters considered
No of points along x=50
No of points along y=50
dx=dy=0.0204
thermal diffusivity 'alpha=1,5'
dt=0.01
CFL along x 'k1=0.9604'
CFL along y 'k2=0.9604'
No of Time steps = 1400
No of iterations to converge 24199, Simulation time 3.51e-05s
- Gauss seidel Method
Code ( To be read in continuation of Implicit)
%Gauss Seidel Implicit
T=Tinitial;
Told=Tinitial;
Tprev=Tinitial;
gauss_seidel=1;
b=1;
for nt=1:1400
error2=9e9;
tic;
while(error2>Tol)
for i=2:nx-1
for j=2:ny-1
H=term2*(T(i-1,j)+T(i+1,j));
V=term3*(T(i,j+1)+T(i,j-1));
T(i,j)=term1*Tprev(i,j)+ H+V;
end
end
error2=max(max(abs(Told-T)));
Error2(b)=error2;
Told=T;
GS(b)=gauss_seidel;
gauss_seidel=gauss_seidel+1;
b=b+1;
end
Tprev=T;
time2=toc;
end
figure(2)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel Iter solver for implicit',tit2})
colorbar
pause(0.3)Plot for Temperature
Parameters considered
No of points along x=50
No of points along y=50
dx=dy=0.0204
thermal diffusivity 'alpha=1,5'
dt=0.01
CFL along x 'k1=0.9604'
CFL along y 'k2=0.9604'
No of Time steps = 1400
No of Iteration to converge 14410, Simulation Time 2.7e-05 sec
-SOR
Code for SOR (To be read in continuation with Gauss Seide Implicit)
%SOR Implicit
Told=Tinitial;
Tprev=Tinitial;
T=Tinitial;
omega=1.5
SOR=1;
c=1;
for nt=1:1400
error3=9e9;
tic;
while(error3>Tol)
for i=2:nx-1
for j=2:ny-1
H=term2*(T(i-1,j)+T(i+1,j));
V=term3*(T(i,j+1)+T(i,j-1));
T(i,j)=(term1*Tprev(i,j)+ H+V)*omega+(1-omega)*Told(i,j);
end
end
error3=max(max(abs(Told-T)));
Error3(c)=error3;
Told=T;
SR(c)=SOR;
SOR=SOR+1;
c=c+1;
time3=toc;
end
Tprev=T;
end
figure(3)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR Iter solver for implicit',tit2})
colorbar
pause(0.3)
figure(4)
semilogy(JI,Error1,GS,Error2,'r',SR, Error3,'k');
legend('Jacobi','Gauss Seidel','SOR')Plot for Temperature
Parameters considered
No of points along x=50
No of points along y=50
dx=dy=0.0204
thermal diffusivity 'alpha=1,5'
dt=0.01
relaxation factor omega=1.5
CFL along x 'k1=0.9604'
CFL along y 'k2=0.9604'
No of Time steps = 1400
No of Iteration to converge 6614, Simulation Time 3.3e-05 sec
-Comparison of decrease of error with iteration for 3 Methods
From the above graph we note that in SOR method the rate of decrease of error is much higher compared to Gauss seidel and Jacobi. Jacobi takes much longer time since we continue to use the old guess values instead of the newly calculated better guess values.
%Full COde for Steady State
clear all
close all
clc
%Input Parameter
nx=50;
ny=50;
x=linspace(0,1,nx);%Grid along X
y=linspace(0,1,ny);%Grid along y
dx=x(2)-x(1);
dy=y(2)-y(1);
T=ones(nx,ny);
%Boundary Conditions
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Tinitial=T;
Told=T;
k=2*((dx^2+dy^2)/(dx^2*dy^2));
Tol=1e-4
error1=9e9
error2=9e9
error3=9e9
jacobi_iter=1;
tic
a=1;
while (error1>Tol)
for i=2:nx-1
for j=2:ny-1
m=(Told(i+1,j)+Told(i-1,j))/(dx*dx);
n=(Told(i,j+1)+Told(i,j-1))/(dy*dy);
T(i,j)=(m+n)/k ;
end
end
error1= max(max(abs(Told-T)));
ERROR1(a)=error1; %Plotting error vs iterations
Told=T;
jacobi_iter=jacobi_iter+1;
JI(a)=jacobi_iter;%Storing iterations in counter
time1=toc;
a=a+1;
end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi solver for steady state',tit1})
colorbar
pause(0.3)
%Gauss Seidel
tic;
T=Tinitial;
Told=Tinitial;
gauss_seidel=1;
b=1;
while (error2>Tol)
for i=2:nx-1
for j=2:ny-1
m=(T(i+1,j)+T(i-1,j))/(dx*dx);
n=(T(i,j+1)+T(i,j-1))/(dy*dy);
T(i,j)=(m+n)/k ;
end
end
error2= max(max(abs(Told-T)));
Told=T;
ERROR2(b)= error2;
GS(b)=gauss_seidel;
gauss_seidel=gauss_seidel+1;
b=b+1;
time2=toc;
end
figure(2)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel solver for steady state',tit2})
colorbar
pause(0.3)
%Successive Over Relaxation
SOR=1;
T=Tinitial
Told=Tinitial;
omega=1.5;
tic
c=1;
while (error3>Tol)
for i=2:nx-1
for j=2:ny-1
m=(T(i+1,j)+T(i-1,j))/(dx*dx);
n=(T(i,j+1)+T(i,j-1))/(dy*dy);
T(i,j)=(1-omega)*Told(i,j) +omega*((m+n)/k );
end
end
error3= max(max(abs(Told-T)));
ERROR3(c)= error3;
Told=T;
SR(c)=SOR;
SOR=SOR+1;
time3=toc;
c=c+1;
end
figure(3)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit3=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR solver for steady state',tit3})
colorbar
pause(0.3)
figure(4)
semilogy(JI,ERROR1,GS,ERROR2,'r',SR, ERROR3,'k');
legend('Jacobi','Gauss Seidel','SOR')
xlabel('No of iterations')
ylabel('log of absolute error')
title ('Comaparison of rate of decrease of error with iterations for three methods')
%Full Code for Transient explicit
clear all
close all
clc
nx=50;
ny=50;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=4e-3;
dt=0.01;
k1=(alpha*dt)/(dx)^2;
k2=(alpha*dt)/(dy)^2;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
Told=T;
T(1,1)=(400+900)/2;
T(end,1)=(400+600)/2;
T(1,end)=(900+800)/2;
T(end,end)=(800+600)/2;
% k=2*((dx^2+dy^2)/(dx^2*dy^2));
Tol=1e-4;
error1=9e9;
for nt=1:15000
tic;
for i=2:nx-1
for j=2:ny-1
m=k1*(Told(i+1,j)+Told(i-1,j)-2*Told(i,j));
n=k2*(Told(i,j+1)+Told(i,j-1)-2*Told(i,j));
T(i,j)=m+n+Told(i,j) ;
end
end
time1=toc;
Told=T;
end
figure(1);
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Simulation Time=%d s',time1)
title({'Solution for transient explicit '})
colorbar
pause(0.3)
%Full COde for Implicit
clear all
close all
clc
nx=50;
ny=50;
x=linspace(0,1,nx);%Grid Spacing along X
y=linspace(0,1,ny);%Grid Spacing along Y
dx=x(2)-x(1);
dy=y(2)-y(1);
alpha=1.5;%Thermal Diffusivity
dt=1e-2;%Time Step
k1=(alpha*dt)/(dx)^2;%CFL along x
k2=(alpha*dt)/(dy)^2;%CFL along y
term1=(1/(1+2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
T=ones(nx,ny);
T(1,:)=600;
T(end,:)=900;
T(:,1)=400;
T(:,end)=800;
T(1,1)=(400+600)/2;
T(end,1)=(400+800)/2;
T(1,end)=(900+800)/2;
T(end,end)=(800+600)/2;
Told=T;
Tinitial=T;
Tprev=T;
Tol=1e-4;
jacobi_iter=1;
a=1
for nt=1:1400
error1=9e9;
tic;
while(error1>Tol)
for i=2:nx-1
for j=2:ny-1
H=term2*(Told(i-1,j)+Told(i+1,j));
V=term3*(Told(i,j+1)+Told(i,j-1));
T(i,j)=term1*Tprev(i,j)+ H+V;
end
end
error1=max(max(abs(Told-T)));
Error1(a)=error1;
Told=T;%Spatial Temp values updated within the space iteration
jacobi_iter=jacobi_iter+1;
JI(a)=jacobi_iter;
a=a+1;
time1=toc;
end
Tprev=T;%After convergence the T values updated
end
figure(1)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit1=sprintf('Number of iteration=%d\n Simulation Time=%d s',jacobi_iter,time1)
title({'Jacobi Iter solver for implicit',tit1})
colorbar
pause(0.3)
%Gauss Seidel Implicit
T=Tinitial;
Told=Tinitial;
Tprev=Tinitial;
gauss_seidel=1;
b=1;
for nt=1:1400
error2=9e9;
tic;
while(error2>Tol)
for i=2:nx-1
for j=2:ny-1
H=term2*(T(i-1,j)+T(i+1,j));
V=term3*(T(i,j+1)+T(i,j-1));
T(i,j)=term1*Tprev(i,j)+ H+V;
end
end
error2=max(max(abs(Told-T)));
Error2(b)=error2;
Told=T;
GS(b)=gauss_seidel;
gauss_seidel=gauss_seidel+1;
b=b+1;
end
Tprev=T;
time2=toc;
end
figure(2)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',gauss_seidel,time2)
title({'Gauss Seidel Iter solver for implicit',tit2})
colorbar
pause(0.3)
%SOR Implicit
Told=Tinitial;
Tprev=Tinitial;
T=Tinitial;
omega=1.5
SOR=1;
c=1;
for nt=1:1400
error3=9e9;
tic;
while(error3>Tol)
for i=2:nx-1
for j=2:ny-1
H=term2*(T(i-1,j)+T(i+1,j));
V=term3*(T(i,j+1)+T(i,j-1));
T(i,j)=(term1*Tprev(i,j)+ H+V)*omega+(1-omega)*Told(i,j);
end
end
error3=max(max(abs(Told-T)));
Error3(c)=error3;
Told=T;
SR(c)=SOR;
SOR=SOR+1;
c=c+1;
time3=toc;
end
Tprev=T;
end
figure(3)
colormap(jet);
contourf(x,y,T,'showText','on')
set(gca,'YDir','reverse')
tit2=sprintf('Number of iteration=%d\n Simulation Time=%d s',SOR,time3)
title({'SOR Iter solver for implicit',tit2})
colorbar
pause(0.3)
figure(4)
semilogy(JI,Error1,GS,Error2,'r',SR, Error3,'k');
xlabel('ITeration')
ylabel('log of absolute error')
title('Comparison of error with iterations')
legend('Jacobi','Gauss Seidel','SOR')