Conduction equation using a Transient solver and a Steady state solver using Iterative techniques Jacobi Gauss Seidal SOR

REPORT

AIM

The main objective is to solve conduction equation using a Transient solver and a Steady state solver using Iterative techniques (Jacobi,Gauss Seidal,SOR)

2 Dimensional Steady State Heat Conduction Equation is as follows

`(partial^2T)/(partialx^2)+(partial^2T)/(partialy^2)=0` eq1

discritize eq1 as

`frac{(Tl-2Tp+Tr)}(trianglex^2)+frac{(Tt-2Tp+Tb)}(triangley^2)=0`

where

l stands for left node

p stands for present node

t stands for top node

b stands for left node

r stands for right node

simplify the above equation as

`frac{(Tl+Tr)}(trianglex^2)+frac{(Tt+Tb)}(triangley^2)=frac{(2Tp)}(trianglex^2)+frac{(2Tp)}(triangley^2`

Taking 2Tp common on R.H.S of above equation

`frac{(Tl+Tr)}(trianglex^2)+frac{(Tt+Tb)}(triangley^2)=2Tp(frac{1}(trianglex^2)+frac{1}(triangley^2))`

The above equation can be writteen as

`TpK=frac{(Tl+Tr)}(trianglex^2)+frac{(Tt+Tb)}(triangley^2)` eq2

where

`K=(2trianglex^2+2triangley^2)/(trianglex^2triangley^2)`

eq 2 can be written as

`Tp=1/k((Tl+Tr)/(trianglex^2)+(Tt+Tb)/(triangley^2))`

As nx = ny [Number of points along the x direction is equal to the number of points along the y direction].

so can take same step size along x and y direction as nx=ny

△x=△y=h

putting values △x and △y, k reduces to

`K=(2h^2+2h^2)/h^4`

`K=4/h^2`

Putting the value of k and replacing

p by i,j

l by i-1,j

t by i,j+1

b by i,j-1

The finnal form of eq 2 will be

`T(i,j)=[T(i-1,j)+T(i+1,j)+T(i,j+1)+T(i,j-1)]/4` eq3

JACOBI METHOD

`T^(k+1)(i,j)=[T^k(i-1,j)+T^k(i+1,j)+T^k(i,j+1)+T^k(i,j-1)]/4` eq4

GAUS SIEDLE METHOD

`T^(k+1)(i,j)=[T^(k+1)(i-1,j)+T^k(i+1,j)+T^k(i,j+1)+T^(K+1)(i,j-1)]/4` eq5

SOR METHOD

`T^(k+1)(i,j)=w.[T^(k+1)(i-1,j)+T^k(i+1,j)+T^k(i,j+1)+T^((k+1)(i,j-1))]/4+(1-w)T^k(i,j)` eq6

where k represents itterations

Selection of w for SOR Method

we know that 0<w<2

but for the best choice i selected different values of w and noted the number of itterations

when w=1 no. of itterations = 841

Note that w=1 gives gauss siedle method

when w=1.1 no. of itterations = 689

when w=1.2 no. of itterations = 553

when w=1.3 no. of itterations = 441

when w=1.4 no. of itterations = 329

when w=1.5 no. of itterations = 209

when w=1.6 no. of itterations = 249

when w=1.7 no. of itterations = 369

when w=1.8 no. of itterations = 585

when w=1.9 no. of itterations = 1193

So, we can see that the best choice for w=1.5, as it is giving minimum nuumber of itterations

Selection Of Temperature At Corners

At the corner points i took average of both temperatures

At the left bottom corner i took (400+900)/2

At the left top corner i took (400+600)/2

At the right bottom corner i took (800+900)/2

At the right top corner i took (600+800)/2

CODE

%solve  2nd  order 2 diimensional pde
%(d^2t/dx^2) + (d^2t/dy^2)
%x=0 t = 400 k
%x=1 t = 800 k
%y=0 t = 900 k
%y=1 t = 600 k
clear all
close all
clc
% inputs
w=1.5;
alpha=1.4;
dt=0.01;
%k1=(dt*alpha)/(dx^2);

L=1;
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
t=zeros(length(x),length(y));
dx=L/(nx-1);
k1=(dt*alpha)/(dx^2);
% boundry conditions
t(1,1)= (900+400)/2; % corner joining left  and botom
t(2:nx-1,1) = 900;% bottom
t(1,2:ny-1)= 400; %left
t(nx,2:ny-1)=800; %right
t(2:nx-1,ny)=600; %top
t(nx,1)=(800+900)/2;  %corner right n bottom
t(nx,ny)=(800+600)/2;%corner  right and  top
t(1,ny)=(600+400)/2; %corner top and left
told=t;
error = 9e9;
tol=1e-4;
iterative_solver=2;
% jacobi method
if iterative_solver==1;
     jacobi_iter=1;
while (error>tol)
for j=2:ny-1
    for i=2:nx-1
        t(i,j)=0.25*(told(i-1,j)+told(i+1,j)+told(i,j+1)+told(i,j-1));
    end

 error=max(max(abs(told-t)))
  
told=t;
 jacobi_iter= jacobi_iter+1;

 pause(0.03);
end
end
surf(x,y,t);
title_text= sprintf(\'iteration number = %d\',jacobi_iter);
title(title_text)
xlabel(\'x axis\');
ylabel(\'y axis\');
zlabel(\'temperature\');
end
% gauss seidle
if iterative_solver==2;
     gaus_siedle_iter=1;
while (error>tol)
for j=2:ny-1
    for i=2:nx-1
        t(i,j)=0.25*(t(i-1,j)+told(i+1,j)+told(i,j+1)+t(i,j-1));
    end


 error=max(max(abs(told-t)))
told=t;

 gaus_siedle_iter= gaus_siedle_iter+1;
   
end

surf(x,y,t);
title_text=sprintf(\'iteration number = %d\',gaus_siedle_iter);
    title(title_text)
xlabel(\'x axis\');
ylabel(\'y axis\');
zlabel(\'temperature\');
end
end
% SOR METHOD

if iterative_solver==3;
     sor_iter=1;
while (error>tol)
for j=2:ny-1
    for i=2:nx-1
        t(i,j)=w*0.25*(t(i-1,j)+told(i+1,j)+told(i,j+1)+t(i,j-1))+((1-w)*told(i,j));
    end

 error=max(max(abs(told-t)))
  
told=t;
 sor_iter= sor_iter+1;
 %sprintf(\'iteration number = %d\',sor_iter);
 pause(0.03);
end
end
surf(x,y,t);
title_text=sprintf(\'iteration number = %d\',sor_iter);
title(title_text)
xlabel(\'x axis\');
ylabel(\'y axis\');
zlabel(\'temperature\');
end

Solved equation 4 for jaccobi method. equation 5 for Gauss siedle method and equation 6 for SOR method. First given inputs. initialize the arrays and given bounry conditions. Applied loops and soved the desired equations. In the end graph is obtained using \"surf\" comand. i used this instead of contorf comand to observe the boundry conditions. as graph is lifted upwards because average is appplied on the corners.

The code for jaccobi, guass siedle and SOR is given bellow

RESULTS

Jaccobi method give results at 1217 number of itterations

Gauss Siedle method gives results at 841 numbber of itterations

SOR methood gives results at 209 number of itterations.

so, we can say that SOR is best method amoung all three. And Gauss Siedle is better then Jaccobi, as we use latest values for T(i-1,j) and T(i,j-1).

GRAPH OF JACCOBI METHOD

GRAPH OF GAUSS SIEDLE METHOD

GRAPH OF SOR METHOD

NOTE:

Graphs are lifted upward from the corners as i have taken average value at the corners.

2 Dimensional Un Steady State Heat Conduction Equation is as follows

`(delT)/(delt)=alpha(nabla^2T)`

`(Tp^((k+1))-Tp^k)/(trianglet)=alpha[(Tl-2Tp+Tr)/(trianglex^2)+(Tt-2Tp+Tb)/(triangley^2)]` eqA

I) Implicit Case

II) Explicit Case

I) For Implicit Case

For implicit case the above equation A will get the form

`(Tp^((k+1))-Tp^k)/(trianglet)=alpha[(Tl-2Tp+Tr)/(trianglex^2)+(Tt-2Tp+Tb)/(triangley^2)]^(k+1)`

`Tp^((k+1))-Tp^k=alpha*trianglet[(Tl-2Tp+Tr)/(trianglex^2)+(Tt-2Tp+Tb)/(triangley^2)]^(k+1)`

so can take same step size along x and y direction as nx=ny

△x=△y=h

putting values △x and △y, k reduces to

`Tp^((k+1))-Tp^k=(alpha*trianglet)/h^2[Tl-2Tp+Tr+Tt-2Tp+Tb)]^(k+1)``

put

`k1=(trianglet*alpha)/h^2`

`Tp^((k+1))-Tp^k=k1[Tl-2Tp+Tr+Tt-2Tp+Tb)]^(k+1`

`Tp^((k+1))-Tp^k=k1[Tl+Tr+Tt+Tb)]^(k+1)-4k1*Tp^(k+1)`

`Tp^((k+1))(1+4*K1)=k1[Tl+Tr+Tt+Tb)]^(k+1)+Tp^k`

`Tp^((k+1)=(k1)/(1+4*k1)*[(Tl+Tr+Tt+Tb)^(k+1)]+(Tp^k)/(1+4*k1)`

the final form of equation will be as folllows

`Tp^((k+1)=(k1)/(1+4*k1)*[Tl^(k+1)+Tr^(k+1)+Tt^(k+1)+Tb^(k+1)]+(Tp^k)/(1+4*k1)`

NOTE: n represents itterations

JACOBI METHOD

`(n+1)_(Tp^(k+1)) =(k1)/(1+4*k1)[(n_(Tl^(k+1)))+(n_(Tr^(k+1)))+( n_(Tt^(k+1)))+(n_(Tb^(k+1))]+(Tp^k)/(1+4*k1)` eq7

GAUSS SIEDLE METHOD

`(n+1)_(Tp^(k+1)) =(k1)/(1+4*k1)[(n+1)_(Tl^(k+1)))+(n_(Tr^(k+1)))+( n_(Tt^(k+1)))+(n+1)_(Tb^(k+1))]+(Tp^k)/(1+4*k1` eq8

SOR METHOD

`(n+1)_(Tp^(k+1)) =w*{(k1)/(1+4*k1)[(n+1)_(Tl^(k+1)))+(n_(Tr^(k+1)))+( n_(Tt^(k+1)))+(n+1)_(Tb^(k+1))]+(Tp^k)/(1+4*k1]}+(1-w)(n_Tp^(k+1))` eq9

where

l stands for left node

`Tl=T_(i-1,j)`

p stands for present node

`Tp=T_(i,j)`

t stands for top node

`Tt=T_(i,j+1)`

b stands for left node

`Tb=T_(i,j-1)`

r stands for right node

`Tr=T_(i+1,j)`

Selection of w for SOR Method

we know that 0<w<2

but for the best choice i selected different values of w and noted the number of itterations

when w=1 no. of itterations = 10459

Note that w=1 gives gauss siedle method

when w=1.1 no. of itterations = 6364

when w=1.2 no. of itterations = 5104

when w=1.28 no. of itterations = 4141

when w=1.3 no. of itterations = 4204

when w=1.4 no. of itterations = 4897

when w=1.5 no. of itterations = 6004

when w=1.6 no. of itterations = 7849

when w=1.7 no. of itterations = 10981

when w=1.8 no. of itterations = 17272

when w=1.9 no. of itterations = 36451

So, we can see that the best choice for w=1.5, as it is giving minimum nuumber of itterations

CODE

Solved equation 7, 8 and 9 for jaccobi, gauss siedle and SOR methods. Given inputs. initialized the arrays and given boundry conditions. averaging is used at the corners. applied loops and solved the desired equations. ploted a graph using conourf command and observed temperature distribution using colorbar command

%solve  2nd  order 2 diimensional pde
%(d^2t/dx^2) + (d^2t/dy^2)
%x=0 t = 400 k
%x=1 t = 800 k
%y=0 t = 900 k
%y=1 t = 600 k
clear all
close all
clc
% inputs
w=1.28;
L=1;
nx=11;
dx=L/(nx-1);
alpha=1.4;
dt=0.01;
k1=(alpha*dt)/dx^2;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
t=zeros(length(x),length(y));
 % boundry conditions
t(1,1)= (900+400)/2; % corner joining left  and botom
t(2:nx-1,1) = 900;% bottom
t(1,2:ny-1)= 400; %left
t(nx,2:ny-1)=800; %right
t(2:nx-1,ny)=600; %top
t(nx,1)=(800+900)/2;  %corner right n bottom
t(nx,ny)=(800+600)/2;%corner  right and  top
t(1,ny)=(600+400)/2; %corner top and left
told=t;
tpre=t;
iterative_solver=3;
% jacobi method
jacobi_iter=1;
if iterative_solver==1;
    %time loop
for k=1:100
       error = 1;
tol=1e-4;
while (error>tol)
for i=2:nx-1
    for j=2:ny-1
        t(i,j)=tpre(i,j)/(1+4*k1)+(k1/(1+4*k1))*(told(i-1,j)+told(i+1,j)+told(i,j+1)+told(i,j-1));
    end

 error=max(max(abs(told-t)));
  
    told=t;

 jacobi_iter= jacobi_iter+1;
end
end
fprintf(\'TIME: %i ERROR:  %f5.6 Iteratiion:  %i TOTAL ITTER %i\\n\',k,error,jacobi_iter)

contourf(x,y,t);
title_text= sprintf(\'iteration number = %d\',jacobi_iter);
title(title_text)
xlabel(\'x axis\');
ylabel(\'y axis\');
zlabel(\'temperature\');

tpre=t;
colorbar
end
end

% gauss seidle
if iterative_solver==2;
    gaus_siedle_iter=1;
    
     %time loop
    for k=1:100
        error = 1;
tol=1e-4;
    
while (error>tol)
for j=2:ny-1
    for i=2:nx-1
        t(i,j)=tpre(i,j)/(1+4*k1)+(k1/(1+4*k1))*(t(i-1,j)+told(i+1,j)+told(i,j+1)+t(i,j-1));
    end


 error=max(max(abs(told-t)));
told=t;

 gaus_siedle_iter= gaus_siedle_iter+1;
end
end
fprintf(\'TIME: %i ERROR:  %f5.6 Iteratiion:  %i\\n\',k,error,gaus_siedle_iter)
contourf(x,y,t);
title_text=sprintf(\'iteration number = %d\',gaus_siedle_iter);
    title(title_text)
xlabel(\'x axis\');
ylabel(\'y axis\');
zlabel(\'temperature\');
tpre=t;
colorbar
end
end
% SOR METHOD

if iterative_solver==3;
      sor_iter=1;
     %time loop
    for k=1:100
        
        error = 1;
tol=1e-4;
while (error>tol)
for j=2:ny-1
    for i=2:nx-1
        %M=tpre(i,j)/(1+4*k1)+(k1/(1+4*k1))*(t(i-1,j)+told(i+1,j)+told(i,j+1)+t(i,j-1));
        t(i,j)=tpre(i,j)/(1+4*k1)+(k1/(1+4*k1))*(t(i-1,j)+told(i+1,j)+told(i,j+1)+t(i,j-1));
        t(i,j)=w*t(i,j)+(1-w)*told(i,j);
       % t(i,j)=w*(tpre(i,j)/(1+4*k1)+(k1/(1+4*k1))*(t(i-1,j)+told(i+1,j)+told(i,j+1)+t(i,j-1)))+(1-w)*told(i,j);
    end

 error=max(max(abs(told-t)));
  
told=t;
 sor_iter= sor_iter+1;
 %sprintf(\'iteration number = %d\',sor_iter);
end
end
fprintf(\'TIME: %i ERROR:  %f5.6 Iteratiion:  %i\\n\',k,error,sor_iter)
contourf(x,y,t);
title_text=sprintf(\'iteration number = %d\',sor_iter);
title(title_text)
xlabel(\'x axis\');
ylabel(\'y axis\');
zlabel(\'temperature\');
tpre=t;
colorbar
end
end

GRAPHS

OBSERVATIONS

It is observed that jacobi method takes more itterations to converge and Gauss Siedle is better then jacobi mehtod. SoR method gives good results for perfect choice of w.

II) Explicit Case

For explicit case the equation A will get the form

`(Tp^((k+1))-Tp^k)/(trianglet)=alpha[(Tl-2Tp+Tr)/(trianglex^2)+(Tt-2Tp+Tb)/(triangley^2)]^k`

`(Tp^((k+1))-Tp^k)=(alpha*trianglet)[(Tl-2Tp+Tr)/(trianglex^2)+(Tt-2Tp+Tb)/(triangley^2)]^k`

so can take same step size along x and y direction as nx=ny

△x=△y=h

putting values △x and △y, k reduces to

`(Tp^((k+1))-Tp^k)=(alpha*trianglet)/h^2[(Tl-4Tp+Tr+Tt+Tb]^k`

`Tp^((k+1))=(alpha*trianglet)/h^2[(Tl-4Tp+Tr+Tt+Tb]^k+Tp^k`

`Tp^((k+1))=k1[(Tl-4Tp+Tr+Tt+Tb]^k+Tp^k`

where

`k1=(alpha*trianglet)/h^2`

`Tp^(k+1)=k1*Tl^k-4*k1*Tp^k+k1*Tr^k+k1*Tt^k+k1*Tb^k+Tp^k`

The final form of equation becomes