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

Solving the 2D heat conduction equation using iterative techniques:

The iterative techniques used in this are as followed:

1. Jacobi

2. Gauss-seidel

3. Successive over-relaxation

First, we are solving for steady-state equation and then Unsteady/Transient state for Explicit and Implicit.

The steady-state equation for 2D:

The equation for steady-state is

`(d^2 T)/dx^2 + (d^2T)/dy^2=0`

For a unit square of length 1m, the equation is being solved by given temperatures into follows,

`Tp = 1/k((TL+TR)/dx^2) +1/k((T T+TB)/dy^2)`

`TP = temperature at point, T(i,j)`,

`TL =temperature at Left,400 degrees, T(i-1,j)`

`TR=temperature at Right,800 degrees,T(i+1,j)`

`T T=temperature at Top,600 degrees,T(i,j-1)`

`TB= temperature at Bottom,900 degrees, T(i,j+1)`

By coding this equation in matlab in convergence loop, we get as follows

%Steady state using iterative methods

clear all
close all
clc


nx=10;
ny=nx;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = x;
dy=dx;
error =9e9;
tol = 1e-4;
 w = 1.5; %relaxation factor for SOR method
T = 298*ones(nx,ny); 


T(:,1)= 400; % Left Temperature
T(1,:)= 600; % Top Temperature
T(:,end)= 800; %Right Temperature
T(end,:)=900; %bottom Temperature


Told= T;
k= 2*(dx^2+dy^2)/(dx^2*dy^2);



jacobi_SS (nx,ny,error,tol,T,Told,k,x,y,dx,dy)
GS_SS(nx,ny,error,tol,T,Told,k,x,y,dx,dy)
SOR_SS(nx,ny,error,tol,T,Told,k,x,y,dx,dy,w)

Function for Jacobi:

%jacobi
function jacobi_SS(nx,ny,error,tol,T,Told,k,x,y,dx,dy)
[X,Y]=meshgrid(x,y);
    jacobi_iter = 1;
    tic;
while  error>tol
    
    for i = 2:nx-1
        for j = 2:ny-1
            T(i,j) = 1/k*((Told(i-1,j)+Told(i+1,j))/dx^2 + (Told(i,j-1)+Told(i,j+1))/dy^2);
        end
    end
    error = max(abs(Told(:) - T(:)));
    Told=T;
    jacobi_iter = jacobi_iter+1;
    toc;
    time = toc;
    
    %plot
    figure(1)
    contourf(X,Y,T)
    set(gca,'YDIR','reverse')
    colormap(jet)
    %clabel(C,h)
    colorbar
    xlabel('X');
    ylabel('Y');
    title(sprintf('Jacobi using steady state \n number of iterations = %d \n computaional time = %f s',jacobi_iter,time))
   
end
end

Function for Gauss-seidel:

%GS
function GS_SS(nx,ny,error,tol,T,Told,k,x,y,dx,dy)
[X,Y]=meshgrid(x,y);
    GS_iter = 1;
    tic;
while  error>tol
    
    for i = 2:nx-1
        for j = 2:ny-1
            T(i,j) = 1/k*((T(i-1,j)+Told(i+1,j))/dx^2 + (T(i,j-1)+Told(i,j+1))/dy^2);
        end
    end
    error = max(abs(Told(:) - T(:)));
    Told=T;
    GS_iter = GS_iter+1;
    toc;
    time=toc;
    
    %plot
    figure(2)
    contourf(X,Y,T)
    set(gca,'YDIR','reverse')
    colormap(jet)
    %clabel(C,h)
    colorbar
    xlabel('X');
    ylabel('Y');
    title(sprintf('GS using steady state \n number of iterations = %d\n computaional time = %f s ',GS_iter,time))
   
end

Function for Successive over relaxation:

%SOR
function SOR_SS(nx,ny,error,tol,T,Told,k,x,y,dx,dy,w)
[X,Y]=meshgrid(x,y);
tic;
    SOR_iter = 1;
while  error>tol
    
    for i = 2:nx-1
        for j = 2:ny-1
            T(i,j) = (1-w)*Told(i,j) +w/k*((T(i-1,j)+Told(i+1,j))/dx^2 + (T(i,j-1)+Told(i,j+1))/dy^2);
        end
    end
    error = max(abs(Told(:) - T(:)));
    Told=T;
    SOR_iter = SOR_iter+1;
    toc;
    time=toc;
    
    %plot
    figure(3)
    contourf(X,Y,T)
    set(gca,'YDIR','reverse')
    colormap(jet)
    %clabel(C,h)
    colorbar
    xlabel('X');
    ylabel('Y');
    title(sprintf('SOR using steady state \n number of iterations = %d \n computaional time = %f s',SOR_iter,time))
   
end

Output:

The following output for these 3 iterative techniques used to solve 2D steady state is shown below

as we could see, the computation time dropped drastically using Successive over-relaxation iterative compared to other two iterative methods.

The number of iterations is also low in number for SOR iterative since omega helps to shoot the value closer to the solution fast compare to other iterative technique.

relaxation factor(i.e., omega) is 1.5 for the least amount of iterations and solution blows off if >2. Thus, SOR iterative is preferred for solving 2D Steady State.

Solving 2D Transient equation using an explicit method:

In this, we do not need to use iterative techniques.

The equation for Transient equation is:

`(dT)/dt = alpha*((d^2 T)/dx^2 + (d^2T)/dy^2)`

The equation when simplified in the following will be

`T^(n+1)P = T^nP + k1*(TL-2*TP+TR)^n + k2*(T T-2*TP+TB)^n`

where, `k1= alpha*dt/dx^2`

`k2= alpha*dt/dy^2`

and CFL number will be <0.5, `CFL= CFLx+CFLy<=0.5`

and alpha will be 1.4 or 1.5

In this solution alpha is 1.5,

The code using Explicit method is

%Transient state using explicit methods

clear all
close all
clc

nx=11;   %no. of grid points
ny=nx;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = x;
dy=dx;
error =9e9;
tol = 1e-4;
T =298*ones(nx,ny);
nt=1500;     %Time step
dt=0.001;  
alpha = 1.5; %Thermal diffusivity


T(:,1)= 400; % Left Temperature
T(1,:)= 600; % Top Temperature
T(:,end)= 800; %Right Temperature
T(end,:)=900; %bottom Temperature


Told= T;
k1= alpha*dt/dx^2;
k2= alpha*dt/dy^2;
CFL=k1+k2;

[X,Y]=meshgrid(x,y);
iter = 1;
tic;
for k = 1:nt 
    for i = 2:nx-1
        for j = 2:ny-1
            h=Told(i,j);
            l= k1*(Told(i-1,j)- 2*Told(i,j)+Told(i+1,j));
            m= k2*(Told(i,j-1)- 2*Told(i,j)+Told(i,j+1));
            T(i,j)=h+l+m;
        end
    end

error = max(abs(Told(i,j) - T(i,j)));
Told=T; 
iter=iter+1;
      
end

toc
time=toc;

%plot
figure(1)
contourf(X,Y,T)
set(gca,'YDIR','reverse')
colormap(jet)
colorbar
xlabel('X')
ylabel('Y');
title(sprintf('Transient State using explicit method\n number of iterations = %d\n computaional time =%f s \n CFL number = %f',iter,time,CFL))

The output from this code is:

we could see the CFL value is < 0.5 and is stable.if it is >0.5, the solution will be unstable and blown up.

The iterations are 1500 as given and computational time taken is 0.016 sec.

Solving 2D Transient equation using an Implicit method:

We have to solve the Implicit method using three iterative techniques.

The equation for will be same but we will be solving the equation for n+1th power whereas n in the explicit method.

`T^(n+1)P = 1/(1+2*k1+2*k2) *[T^nP + k1*(TL-2*TP+TR)^(n+1) + k2*(T T-2*TP+TB)^(n+1)]`

The following code to solve this equation is

close all
clear all
clc



nx=100;%nodes
ny=nx;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = x;
dy=dx;
nt = 1500;
dt=0.01;
error =9e9;  %error value
tol = 1e-4; %tolerance
w = 1.2; %omega for SOR
alpha = 1.4; %Thermal Diffusivity
T = 298*ones(nx,ny); 

T(:,1)= 400; % Left Temperature
T(1,:)= 600; % Top Temperature
T(:,end)= 800; %Right Temperature
T(end,:)=900; %bottom Temperature


Told= T;
T_prev= T;
k1= alpha*dt/dx^2;
k2= alpha*dt/dy^2;


Jacobi_Transient_implicit(nx,ny,error,tol,T,Told,T_prev,k1,k2,x,y,nt)
GS_Transient_implicit(nx,ny,error,tol,T,Told,T_prev,k1,k2,x,y,nt)
SOR_Transient_implicit(nx,ny,error,tol,T,Told,T_prev,k1,k2,x,y,w,nt)

Function for Jacobi iterative:

%Jacobi for implicit method
function Jacobi_Transient_implicit(nx,ny,error,tol,T,Told,T_prev,k1,k2,x,y,nt)
[X,Y]=meshgrid(x,y);
jacobi_iter = 1;
tic;
for k =1:nt
    while (error>tol)
    
    for i = 2:nx-1
        for j = 2:ny-1
           term1= (1+2*k1+2*k2)^-1;
           term2= k1*term1;
           term3= k2*term1;
           T(i,j)= term1*T_prev(i,j) + term2*(Told(i-1,j)+Told(i+1,j))+term3*(Told(i,j-1)+Told(i,j+1));
        end
    end
    error = max(abs(Told(i,j) - T(i,j)));
    Told=T;
    jacobi_iter = jacobi_iter+1;
    
    end
    T_prev = T;
       
end
    toc;
    time=toc;

    %plot
    figure(1)
    contourf(X,Y,T)
    set(gca,'YDIR','reverse')
    colormap(jet)
    colorbar
    grid on
    xlabel('X');
    ylabel('Y');
    title(sprintf('jacobi using unsteady state by implicit method\n number of iterations = %d\n computaional time = %f s ',jacobi_iter,time))
end

Function for Gauss-sidel

%GS for implicit method
function GS_Transient_implicit(nx,ny,error,tol,T,Told,T_prev,k1,k2,x,y,nt)
[X,Y]=meshgrid(x,y);
GS_iter = 1;
tic;
for k=1:nt
    while (error>tol)
    
    for i = 2:nx-1
        for j = 2:ny-1
           term1= (1+2*k1+2*k2)^-1;
           term2= k1*term1;
           term3= k2*term1;
           T(i,j)= term1*T_prev(i,j) + term2*(T(i-1,j)+Told(i+1,j))+term3*(T(i,j-1)+Told(i,j+1));
        end
    end
    error = max(abs(Told(i,j) - T(i,j)));
    Told=T;
    GS_iter = GS_iter+1;
    
    end
    T_prev = T;
       
end
    toc;
    time=toc;

    %plot
    figure(2)
    contourf(X,Y,T)
    set(gca,'YDIR','reverse')
    colormap(jet)
    colorbar
    grid on
    xlabel('X');
    ylabel('Y');
    title(sprintf('GS iterative using unsteady state by implicit method\n number of iterations = %d\n computaional time = %f s ',GS_iter,time))
end

Function for Successive over relaxation:

%SOR for implicit method
function SOR_Transient_implicit(nx,ny,error,tol,T,Told,T_prev,k1,k2,x,y,w,nt)
[X,Y]=meshgrid(x,y);
SOR_iter = 1;
tic;
for k=1:nt
    while (error>tol)
    
    for i = 2:nx-1
        for j = 2:ny-1
           term1= (1+2*k1+2*k2)^-1;
           term2= k1*term1;
           term3= k2*term1;
           T(i,j)= (1-w)*Told(i,j)+w*(term1*T_prev(i,j) + term2*(T(i-1,j)+Told(i+1,j))+term3*(T(i,j-1)+Told(i,j+1)));
        end
    end
    error = max(abs(Told(i,j) - T(i,j)));
    Told=T;
    SOR_iter = SOR_iter+1;
    
    end
    T_prev = T;
       
end
    toc;
    time=toc;

    %plot
    figure(3)
    contourf(X,Y,T)
    set(gca,'YDIR','reverse')
    colormap(jet)
    colorbar
    grid on
    xlabel('X');
    ylabel('Y');
    title(sprintf('SOR iterative using unsteady state by implicit method\n number of iterations = %d\n computaional time = %f s ',SOR_iter,time))
end

The output for the iteratives used to solve 2D transient equation are:

As we could see from the plots, the SOR iterative has the least number of iterations and also the computational time compared to the other two iterative.