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

In this project, we will solve steady and unsteady 2D heat conduction problem using iterative solvers such as Jacobi, gauss sidle and successive over relaxation (SOR) methods.

For solving the steady state equation of 2D heat conduction:

`∂²t/∂x²+∂²t/∂y²=0`

By using CDS method

`∂²t/∂x²=(T(i-1,j)-2T(i,j)+2(i+1,j))/∆x²`

`∂²t/∂y²=(T(i,j-1)-2T(i,j)+2(i,j+1))/∆y²`

for better understanding

`T(i-1,j)=T_L`

`T(i+1,j)=T_R`

`T(i,j-1)=T_T`

`T(i,j+1)=T_B`

`T(i,j)=T_P`

Therefore the equation becomes

`(T(i-1,j)+2T(i,j)+2(i+1,j))/∆x²+(T(i,j-1)+2T(i,j)+2(i,j+1))/∆y²=0`

i.e `T_P=(1/K)+((T_L+T_R)/(∆x^2 )+(T_T+T_B)/(∆y^2 ))`

Where K is,

`K=(dx^2+dy^2)/(dx^2.dy^2 )`

Code For Steasdy State_Jacobi Method:

clear all
close all
clc
% code to solve a second order 2D ode
% d^2T/dx^2 + d^2T/dy^2;
% length of domain
L=1;
% define no.of nodes
n_x=10;
n_y=n_x;
% gridspacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=x(2)-x(1);
dy=y(2)-y(1);
% dx=L/(n_x-1);
% dy=L/(n_y-1);
% Assign Boundary Conditions
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
% copy of teperature
Told=T;
% value of K
k=2*((dx^2+dy^2)/((dx^2)*(dy^2)));
tol=1e-4;
error=7e7;
% using Jacobi method
jacobi_iteration=1;
% Caluculate simulation time
tic
while(error>tol)
    for i=2:n_x-1;
        for j=2:n_y-1;
            term_1=(1/k)*(((Told(i-1,j)+Told(i+1,j))/(dx^2)));
            term_2=(1/k)*(((Told(i,j-1)+Told(i,j+1))/(dy^2)));
            T(i,j)=term_1+term_2;
          
        end
    end
    % update the error
%     c=abs(Told-T);
%     error=max(c,[],'all');
error=max(max(abs(Told-T)));
    % update the Told
    Told=T;
    jacobi_iteration=jacobi_iteration+1;
     
end

iteration_time=toc;
figure(1);
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
title_name=sprintf('no of iterations=%d iteration time=%d',jacobi_iteration,iteration_time);
title(title_name)

Output:

`Code For Steady State_GS Method`:

clear all
close all
clc
% code to solve a second order 2D ode
% d^2T/dx^2 + d^2T/dy^2;
% length of domain
L=1;
% define no.of nodes
n_x=10;
n_y=n_x;
% gridspacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
% dx=x(2)-x(1);
% dy=y(2)-y(1);
dx=L/(n_x-1);
dy=L/(n_y-1);
% Assign Boundary Conditions
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
% copy of teperature
Told=T;
% value of K
k=2*((dx^2+dy^2)/((dx^2)*(dy^2)));
tol=1e-4;
error=7e7;
% using Gauss shedel method
gs_iteration=1;
% Caluculate simulation time
tic
while(error>tol)
    for i=2:n_x-1;
        for j=2:n_y-1;
            term_1=(1/k)*(((T(i-1,j)+T(i+1,j))/(dx^2)));
            term_2=(1/k)*(((T(i,j-1)+T(i,j+1))/(dy^2)));
            T(i,j)=term_1+term_2;
          
        end
    end
    % update the error
    c=abs(Told-T);
    error=max(c,[],'all');
    % update the Told
    Told=T;
   gs_iteration=gs_iteration+1;
     
end
iteration_time=toc;
figure(1);
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
title_name=sprintf('no of iterations=%d iteration time=%d',gs_iteration,iteration_time);
title(title_name)

Output:

`Code For Steady State_SOR Method`:

clear all
close all
clc
% code to solve a second order 2D ode
% d^2T/dx^2 + d^2T/dy^2;
% length of domain
L=1;
% define no.of nodes
n_x=10;
n_y=n_x;
% gridspacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
% dx=x(2)-x(1);
% dy=y(2)-y(1);
dx=L/(n_x-1);
dy=L/(n_y-1);
% Assign Boundary Conditions
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
% copy of teperature
Told=T;
% value of K
k=2*((dx^2+dy^2)/((dx^2)*(dy^2)));
tol=1e-4;
error=7e7;
% relaxationn factor omega is(>1-over relaxaction,<1-under relaxaction)
omega=1.02;
% using Gauss shedel method
sor_iteration=1;
% Caluculate simulation time
tic
while(error>tol)
    for i=2:n_x-1;
        for j=2:n_y-1;
            term_1=(1/k)*(((T(i-1,j)+T(i+1,j))/(dx^2)));
            term_2=(1/k)*(((T(i,j-1)+T(i,j+1))/(dy^2)));
            T(i,j)=(1-omega)*Told(i,j)+omega*(term_1+term_2);
          
        end
    end
    % update the error
    c=abs(Told-T);
    error=max(c,[],'all');
    % update the Told
    Told=T;
   sor_iteration=sor_iteration+1;
     
end
iteration_time=toc;
figure(1);
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
title_name=sprintf('no of iterations=%d iteration time=%d',sor_iteration,iteration_time);
title(title_name)

Output:

The time required for convergence is0.00100248sec and iteration carried out are 112 which is the lowest number of iterations in steady state solution.Also,the SOR method is faster to reach convergence among other iterative solvers. Also in the SOR method, by increasing the relaxation factor omega 1.01 to 1.5.it is seen that the number of iterations reduces from 114 to 30 iterations. But relaxation factor increases above the 1.5 the number of iterations starts increasing.

For solving the Transient case of 2D heat conduction problem:

`∂T/∂t=α((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 ))`

By forward differencing for the temporal terms LHS

`∂T/∂t=(T_((i,j))^(n+1)-T_((i,j))^n)/∆t`

Using central differencing for the Spacial terms

`∂²T/∂x²=(T(i-1,j)-2T(i,j)+2(i+1,j))/∆x²`

`∂²T/∂y²=(T(i,j-1)-2T(i,j)+2(i,j+1))/∆y²`

The Implicit Form of the equation becomes

`T_P^(n+1)=(T_P^n+K_1 〖(T_L+T_R)〗^(n+1)+K_2 〖(T_T+T_B)〗^(n+1))/(1+〖2K〗_1+2K_2 )`

Where K1 & K2 are,

`K_1=α(∆t/(∆x^2 ))`

`K_2=α(∆t/(∆y^2 ))`

The Explicit Form of the equation becomes

`T_P^(n+1)=(T_P^n+K_1 〖(T_L+T_R)〗^n+K_2 〖(T_T+T_B)〗^n)/(1+〖2K〗_1+2K_2 )`

Code for Explicit form of transient case:

clear all
close all
clc
% solve the unstedy/transient equation d^2T/dt^2=alpha((d^2T/dx^2)+(d^2T/dy^2))
%  domain length
L=1;
% define no.of node points
n_x=10;
n_y=n_x;
% grid spacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[X,Y]=meshgrid(x,y);
dx=L/(n_x-1);
dy=L/(n_y-1);
% boundary conditions
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
% define the temperature terms
alpha=1.1;
nt=1000;
dt=0.0001;
% assign k values
k1=alpha*(dt/(dx^2));
k2=alpha*(dt/(dy^2));
% Read CFL Value
CFL=k1+k2;
% copy of teperature
Told=T;
% explicit calculations
explicit_iteration=1;
% time loop
tic
for t=1:nt
%     spacial loop
for i=2:n_x-1
    for j=2:n_y-1
        term_1=Told(i,j);
        term_2=k1*((Told(i-1,j))-(2*Told(i,j))+(Told(i+1,j)));
        term_3=k2*((Told(i,j-1))-(2*Told(i,j))+(Told(i,j+1)));
        T(i,j)=term_1+term_2+term_3;
    end
end
% update old value of T
Told=T;
explicit_iteration=explicit_iteration+1;
end
iteration_time=toc;
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
colorbar
xlabel('X axix');
ylabel('Y axis');
set(gca,'YDIR','reverse');
title_text=sprintf('no of iterations=%d Iteration time=%d',explicit_iteration,iteration_time);
title(title_text);

Output:

The iterations required to converge the simulation is 1001 and time is 0.083sec. The explicit scheme the results become unstable when dt is crosses 0.002sec and we don’t obtain any solution for it. Also, the explicit scheme takes less time than any other implicit scheme methods.

Implicit form of transient case:

(A)Code For Transient State_Jacobi Method:

clear all
close all
clc
% length of domain
L=1;
% define no.of nodes
n_x=10;
n_y=n_x;
% gridspacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[x,y]=meshgrid(x,y);
dx=L/(n_x-1);
dy=L/(n_y-1);
% define the temperature terms
alpha=1.1;
nt=1000;
dt=0.0001;
% Assign Boundary Conditions
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
% assign k values
k1=alpha*(dt/(dx^2));
k2=alpha*(dt/(dy^2));
% Read CFL Value
CFL=k1+k2;
error=7e7;
tol=1e-4;
% copy of teperature
Told=T;
T_prv_dt=Told;
term_1=1/(1+2*k1+2*k2);
term_2=k1*term_1;
term_3=k2*term_1;
% jacobi method
jacobi_iter=1;
% time loop
tic
for t=1:nt
    error=7e7;
    tol=1e-4;
%     converge loop
while(error>tol)
%     spacial loop
for i=2:n_x-1
    for j=2:n_y-1
        H=(Told(i-1,j))+(Told(i+1,j));
        V=(Told(i,j-1))+(Told(i,j+1));
        T(i,j)=term_1*((T_prv_dt(i,j))+(term_2*H)+(term_3*V));
    end
end
c=abs(Told-T);
    error=max(c,[],'all');
Told=T;
jacobi_iter=jacobi_iter+1;
end
T_prv_dt=T;
end
iteration_time=toc;
figure(1);
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
set(gca,'YDIR','reverse');
title_text=sprintf('no.of iterations=%d iteration time=%d',jacobi_iter,iteration_time);
title(title_text);

Output:

The time required for convergence is 0.0278 sec and iteration carried out are3917 which is the highest number of iterations in transient state solution.

(B)Code For Transient State_GS Method:

clear all
close all
clc
% length of domain
L=1;
% define no.of nodes
n_x=10;
n_y=n_x;
% gridspacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[x,y]=meshgrid(x,y);
dx=L/(n_x-1);
dy=L/(n_y-1);
% alpha=1;
alpha=1.1;
nt=1000;
dt=0.0001;
% dt=0.001;
% Assign Boundary Conditions
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
% assign k values
k1=alpha*(dt/(dx^2));
k2=alpha*(dt/(dy^2));
% Read CFL Value
CFL=k1+k2;
error=7e7;
tol=1e-4;
% copy of teperature
Told=T;
T_prv_dt=Told;
term_1=1/(1+2*k1+2*k2);
term_2=k1*term_1;
term_3=k2*term_1;
% jacobi method
gs_iter=1;
% time loop
tic
for t=1:nt
    error=7e7;
    tol=1e-4;
%     converge loop
while(error>tol)
%     spacial loop
for i=2:n_x-1
    for j=2:n_y-1
        H=(T(i-1,j))+(T(i+1,j));
        V=(T(i,j-1))+(T(i,j+1));
        T(i,j)=term_1*((T_prv_dt(i,j))+(term_2*H)+(term_3*V));
    end
end
c=abs(Told-T);
    error=max(c,[],'all');
Told=T;
gs_iter=gs_iter+1;
end
T_prv_dt=T;
end
iteration_time=toc;
figure(1);
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
set(gca,'YDIR','reverse');
title_text=sprintf('no.of iterations=%d iteration time=%d',gs_iter,iteration_time);
title(title_text);

Output:

The time required for convergence is 0.02013 sec and iteration carried out are 3527 which is the moderate number of iterations in transient state solution.

(C)Code For Transient State_SOR Method:

clear all
close all
clc
% length of domain
L=1;
% define no.of nodes
n_x=10;
n_y=n_x;
% gridspacing
x=linspace(0,L,n_x);
y=linspace(0,L,n_y);
[x,y]=meshgrid(x,y);
dx=L/(n_x-1);
dy=L/(n_y-1);
% Define the temperature terms
alpha=1.1;
nt=1000;
dt=0.0001;
% Assign Boundary Conditions
T=ones(n_x,n_y);
T(:,1)=400;
T(:,10)=800;
T(1,:)=600;
T(10,:)=900;
T(1,1)=(600+400)/2;
T(1,10)=(600+800)/2;
T(10,10)=(800+900)/2;
T(10,1)=(900+400)/2;
k1=alpha*(dt/(dx^2));
k2=alpha*(dt/(dy^2));
% Read CFL Value
CFL=k1+k2;
%  Relaxation factor(optimum value 1 to 1.02)
omega=1.02;
error=7e7;
tol=1e-4;
% copy of teperature
Told=T;
T_prv_dt=Told;
term_1=1/(1+2*k1+2*k2);
term_2=k1*term_1;
term_3=k2*term_1;
% sor method
sor_iter=1;
% time loop
tic
for t=1:nt
    error=7e7;
    tol=1e-4;
%     converge loop
while(error>tol)
%     spacial loop
for i=2:n_x-1
    for j=2:n_y-1
        H=(T(i-1,j))+(T(i+1,j));
        V=(T(i,j-1))+(T(i,j+1));
%         gs=(term_1*(T_prv_dt(i,j)))+(term_2*H)+(term_3*V);
        gs=term_1*((T_prv_dt(i,j))+(term_2*H)+(term_3*V));
        T(i,j)=((1-omega)*Told(i,j))+(omega*(gs));
        
    end
end
error=max(max(abs(Told-T)));
% c=abs(Told-T);
%     error=max(c,[],'all');
Told=T;
sor_iter=sor_iter+1;
end
T_prv_dt=T;
end
iteration_time=toc;
contourf(x,y,T);
[a,b]=contourf(x,y,T);
colormap(jet);
clabel(a,b);
xlabel('X axis');
ylabel('Y axis');
set(gca,'YDIR','reverse');
colorbar
set(gca,'YDIR','reverse');
title_text=sprintf('no.of iterations=%d iteration time=%d',sor_iter,iteration_time);
title(title_text);

Output:

The time required for convergence is 0.0239 sec and iteration carried out are 3298 which is the lowest number of iterations in transient state solution.

Conclusions:

>The steady state solution is fastest schemes to complete the solution.in the view of above Successive Over-Relaxation iterative methods turn out to be the fastest among jacobian, and GS by running the above programme for each case.

>Explicit scheme forms to be the fastest among the transient method. It requires very less number of iterations and computation time to converge compared to implicit scheme.

>Relaxation factor (omega) must be maintained within a range of 1.1-1.3. More or less than the range it may become over-relaxed or under-relaxed.