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

 AIM: 

To solve 2D heat conduction equation for steady and transient state using various iterative techniques such as Jacobi, Gauss Seidel and Successive over relaxation methods using MATLAB. 

Problem Define; 

Heat conduction equation is define by the pde 

`(delT)/(delt)=alpha.(del^2T)/(delx)+(del^2T)/(dely)+(del^2T)/(delz)`

which is the heat distribution in 3D space with the period of time. 

where alpha is known to be \'thermal diffusivity\' (m^2/s) 

Given that, the number of grid points in x direction is equal to y direction and the domain is assumed to be an unit squre. 

Boundary conditions given 

1. Left (TL) = 400K  2. Right (TR) = 800K  3. Top(TT) = 600K 4. Bottom (TB) = 900K

First let\'s define what are the iterative techniques we are using here: 

1. Jacobi

In this iteration technique we assume the initial values at all nodes to calculate the new values, which will updated once we find the all the nodes values and this process will goes on untill we arrive with an converge solution. For different simulateneous equations with many unknowns first we have to make the equations into Diagonally dominant matrix form for converging results. 

2. Gauss-seidel

In this iteration we assume the values obtained in Jacobi iteration for the next equation immediately without waiting for completion of first cycle. 

3. Successive Over-relaxation (SOR)

This method is an updated version of previous two. In addition we use the relaxation factor to ease the process. For this analysis we use the relaxation facor to be greater than 1 and less than 2.

`(delT)/(delt)=alpha.(del^2T)/(delx)+(del^2T)/(dely)`

For steady state there are no change in temperature hence the Steady state becomes

 `(del^2T)/(delx)+(del^2T)/(dely)=0`

For solving steady state 2D using different iterative techniques 

The steady state equation is discretised using 2nd order CDS: 

`(T_L-2T+T_R)/(Deltax^2)+(T_B-2T+T_T)/(Deltay^2)=0`

On simplifying,  

assume,`k=2(Deltax^2+Deltay^2)/(Deltax^2.Deltay^2)`

we can write 

`T_P=(1)/(k)((T_L+T_R)/(Deltax^2)+(T_B+T_T)/(Deltay^2))`

MATLAB CODE FOR STEADY STATE: 

clear all
close 
clc

%% inputs
nx = 20;   % number of grids in x direction
ny = nx;   % number of grids in y direction

tolerance = 1e-5; % tolerance
error = 9e9; % initial guess for the error
T=5;  % time limit
dt = 1e-2;   % time step
x = linspace(0,1,nx);    % x axis mesh
y = linspace(0,1,ny);    % y axis mesh

% grid spacing
dx= x(2)-x(1);
dy= y(2)-y(1); 

[xx,yy] = meshgrid(x,y);

% matrix difination
t= 300*ones(nx,ny);   % Initail temperatures guess
t(end,:) = 900;  % top face
t(1,:) = 600; % Bottom face
t(:,1) = 400;  % left face
t(:,end) = 800;   % right face
k = (2*(dx^2+dy^2)/(dx^2*dy^2));

t_old = t;
t_p = t;

%% Solving method
method = 'which method ?'
solver = input(method);  % \'1\' for jacobi
                         % \'2\' for Gauss-sidel
                         % \'3\' for SOR
%% JACOBI SOLVER

if solver ==1
    iter_1= 1;
    tic;
    % Convergence loop
    while(error>tolerance)
        for i = 2:(nx-1)   % loop of equations
            for j = 2:(ny-1)  % loop of summation
                TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                TB = t_old(i,j-1);
                TT = t_old(i,j+1);
                
                H= (TL +TR)/dx^2;
                V = (TB +TT)/dy^2;
                
                t(i,j) = (H +V)/k;
            end
        end
        time = toc;
        error=max(max(abs(t_old-t)));
        t_old = t;
        iter_1=iter_1 +1;
        
        [x,y]= contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('Steady state using javobi = %d', iter_1);
        title_text2 = sprintf('Simulation time  +%f s', time);
        title_text3 = sprintf('Final error = %f' , error);
        title({title_text1; title_text2;title_text3})
    end
end

%% Gauss-siedel method
if solver ==2
    iter_2= 1;
    tic;
    % Convergence loop
    while(error>tolerance)
        for i = 2:(nx-1)   % loop of equations
            for j = 2:(ny-1)  % loop of summation
                %TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                %TB = t_old(i,j-1);
                TT = t_old(i,j+1);
                
                H= (t(i-1,j) + TR)/dx^2;
                V= (t(i,j-1) + TT)/dy^2;
                
                t(i,j) = (H +V)/k;
            end
        end
        time2 = toc;
        error=max(max(abs(t_old-t)));
        t_old = t;
        iter_2=iter_2 +1;
        time= toc;
        
        [x,y]= contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('Steady state using GS = %d', iter_2);
        title_text2 = sprintf('Simulation time  +%f s', time);
        title_text3 = sprintf('Final error = %f' , error);
        title({title_text1; title_text2;title_text3})
    end

end


%% Successive over relaxation (SOR) method
alpha = 0.3;   % thermal diffusivity
omega = 1.2;   % SOR value, should be greater than 1 and less than 2
if solver ==3
    iter_3= 1;
    tic;
    % Convergence loop
    while(error>tolerance)
        for i = 2:(nx-1)   % loop of equations
            for j = 2:(ny-1)  % loop of summation
                %TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                %TB = t_old(i,j-1);
                TT = t_old(i,j+1);
                
                H= (t(i-1,j) + TR)/dx^2;
                V= (t(i,j-1) + TT)/dy^2;
                tgs(i,j) = (H +V)/k;
                t(i,j) = (1-omega)*(t_old(i,j))+omega*tgs(i,j);
            end
        end
        time3 = toc;
        error=max(max(abs(t_old-t)));
        t_old = t;
        iter_3=iter_3 +1;
        time= toc;
        
        [x,y]= contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('Steady state using SOR = %d', iter_3);
        title_text2 = sprintf('Simulation time  +%f s', time3);
        title_text3 = sprintf('Final error = %f' , error);
        title({title_text1; title_text2;title_text3})
    end

end

output:

Jacobi

Gauss-seidel

Successive over-relaxation

Solving Transient state (Implicit) 

Implicit method uses the parameters which are inter-dependent with each other in the same time step. Most of the mordern day discretization uses this method because of its usefulness in hiher time steps and stable characteristics. 

`T^(n+1)=T^n+k1(T_L-2T+T_R)^(n+1)+k1(T_B-2T+T_T)^(n+1)`

MATHS CODE:

clear all
close all
clc

%% inputs
nx = 20;   % number of grids in x direction
ny = nx;   % number of grids in y direction
alpha = 1.4;   % thermal diffusivity
omega = 1.2;   % SOR value, should be greater than 1 and less than 2

tolerance = 1e-5; % tolerance
error = 9e9; % initial guess for the error
T=5;  % time limit
dt = 1e-2;   % time step
x = linspace(0,1,nx);    % x axis mesh
y = linspace(0,1,ny);    % y axis mesh

% grid spacing
dx= x(2)-x(1);
dy= y(2)-y(1); 

[xx,yy] = meshgrid(x,y);

% matrix difination
t= 300*ones(nx,ny);   % Initail temperatures guess
t(end,:) = 900;  % top face
t(1,:) = 600; % Bottom face
t(:,1) = 400;  % left face
t(:,end) = 800;   % right face
k1= alpha*dt/dx^2;
k2= alpha*dt/dy^2;
k = 1+(2*k1)+(2*k2);

t_old = t;
t_p = t;

%% Solving method
method = 'which method ? '
solver = input(method);  % \'1\' for jacobi
                         % \'2\' for Gauss-sidel
                         % \'3\' for SOR
%% JACOBI SOLVER

if solver ==1
    iter_1= 1;
    tic;
    % Convergence loop
    while(error>tolerance)
        for i = 2:(nx-1)   % loop of equations
            for j = 2:(ny-1)  % loop of summation
                TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                TB = t_old(i,j-1);
                TT = t_old(i,j+1);
                
                H= TL +TR;
                V = TB +TT;
                t(i,j) = ((t_p(i,j)+H*k1 +V*k2)/k);
            end
        end
        time = toc;

        error=max(max(abs(t_old-t)));
        t_old = t;
        iter_1=iter_1 +1;
        
        [x,y]= contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('transient state Implicit using javobi = %d', iter_1);
        title_text2 = sprintf('Simulation time  +%f s', time);
        title_text3 = sprintf('Final error = %f' , error);
        title({title_text1; title_text2;title_text3})
    end
end

%% Gauss-siedel method
if solver ==2
    iter_2= 1;
    tic;
    % Convergence loop
    while(error>tolerance)
        for i = 2:(nx-1)   % loop of equations
            for j = 2:(ny-1)  % loop of summation
                %TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                %TB = t_old(i,j-1);
                TT = t_old(i,j+1);
                H = t(i-1,j) +TR;
                V = t(i,j-1) +TT;
                t(i,j) = ((t_p(i,j)+H*k1 +V*k2)/k);

            end
        end
        time2 = toc;

        error=max(max(abs(t_old-t)));
        t_old = t;
        iter_2=iter_2 +1;
        time= toc;
        
        [x,y]= contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('transient state Implicit using GS = %d', iter_2);
        title_text2 = sprintf('Simulation time  +%f s', time);
        title_text3 = sprintf('Final error = %f' , error);
        title({title_text1; title_text2;title_text3})
    end

end


%% Successive over relaxation (SOR) method
if solver ==3
    iter_3= 1;
    tic;
    % Convergence loop
    while(error>tolerance)
        for i = 2:(nx-1)   % loop of equations
            for j = 2:(ny-1)  % loop of summation
                TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                TB = t_old(i,j-1);
                TT = t_old(i,j+1);
                H = t(i-1,j) +TR;
                V = t(i,j-1) +TT;
                tgs(i,j) = ((t_p(i,j)+H*k1 +V*k2)/k);
                t(i,j) = (1-omega)*(t_old(i,j))+omega*tgs(i,j);
            end
        end
        
        time3 = toc;

        error=max(max(abs(t_old-t)));
        t_old = t;
        iter_3=iter_3 +1;
        time= toc;
        
        [x,y]= contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('transient state Implicit using SOR = %d', iter_3);
        title_text2 = sprintf('Simulation time  +%f s', time3);
        title_text3 = sprintf('Final error = %f' , error);
        title({title_text1; title_text2;title_text3})
    end

end

OUTPUT:

Jacobi

Gauss-seidel

Successive over-relaxation

Transient State 2D (Explicit): 

In the explicit method the next values is estimated based on the already exsisted parameters, which means the n+1 nodes is calculated by nth step parameters.  

Here, the time derivative is discretised by FDS first order while space derivative is discretised by 2nd order CDS; 

`(delT)/(delt)=alphagrad^2T`

`(delT)/(delt)=(T^(n+1)-T^m)/(Deltat)`

`T^(n+1)=T^n+alphaDeltat((T_L-2T+T_R)/(Deltax^2)+(T_B-2T+T_T)/(Deltay^2))`

Assuming `k1=alpha((Deltat)/(Deltax^2)),k2=alpha((Deltat)/(Deltay^2))`

`T^(n+1)+T^n+k1(T_L-2T^n+T_R)+k2(T_B-2T^n+T_T)`

MATHS CODE:

clear all
close all
clc

%% inputs
nx = 10;   % number of grids in x direction
ny = nx;   % number of grids in y direction
alpha = 0.3;   % thermal diffusivity
omega = 1.2;   % SOR value, should be greater than 1 and less than 2

tolerance = 1e-4; % tolerance
error = 9e9; % initial guess for the error
T=5;  % time limit in sec
dt = 0.01;   % time step
nt = T/dt;  
x = linspace(0,1,nx);    % x axis mesh
y = linspace(0,1,ny);    % y axis mesh

% grid spacing
dx= x(2)-x(1);
dy= y(2)-y(1); 

[xx,yy] = meshgrid(x,y);

% matrix difination
t= 300*ones(nx,ny);   % Initail temperatures guess
t(end,:) = 900;  % top face
t(1,:) = 600; % Bottom face
t(:,1) = 400;  % left face
t(:,end) = 800;   % right face
k1= alpha*dt/dx^2;
k2= alpha*dt/dy^2;
k = 1-(2*k1)-(2*k2);   % Explicit constant term

t_old = t;
t_p = t;

%% Solving method
method = 'which method ?   '
solver = input(method);  % \'1\' for jacobi
                         % \'2\' for Gauss-sidel
                         % \'3\' for SOR
%% JACOBI SOLVER
if solver ==1
iter_1 = 1;
nt = T/dt;
tic
for m = 1:nt
    while (error>tolerance)
        for i = 2:nx-1
            for j = 2:ny-1
                TL = t_old(i-1,j);
                TR = t_old(i+1,j);
                TB = t_old(i,j-1);
                TT = t_old(i,j+1);
                H = TL +TR;
                V = TB +TT;
                t(i,j) = t_p(i,j)*k + k1*H +k2*V;
            end
        end
        error = max(max(abs(t-t_old)));
        t_old = t;
        iter_1 = iter_1 +1;
    end
    t_p = t;
end
time1 = toc;
[x,y]= contourf(xx,yy,t);
        colormap(jet)
        colorbar
        clabel(x,y)
        pause(0.03)
        title_text1 = sprintf('transient state Explicit using Jacobi = %d', iter_1);
        title_text2 = sprintf('Simulation time  +%f s', time1);
        title_text3 = sprintf('Final error = %f' , error);
        title({title_text1; title_text2;title_text3})
end
%% Gauss-siedel method
if solver ==2
    iter_2 = 1;
    nt = T/dt;
    tic
    for m = 1:nt
        while (error>tolerance)
            for i = 2:nx-1
                for j = 2:ny-1
                    %TL = t_old(i-1,j);
                    TR = t_old(i+1,j);
                    %TB = t_old(i,1j-1);
                    TT = t_old(i,j+1);
                    H = t(i-1,j) +TR;
                    V = t(i,j-1) +TT;
                    t(i,j) = t_p(i,j)*k +H*k1 +V*k2;
                end
            end
            error = max(max(abs(t-t_old)));
            t_old = t;
            iter_2 = iter_2 +1;
        end
        t_p = t;
    end
    time2 = toc;
    [x,y]= contourf(xx,yy,t);
    colormap(jet)
    colorbar
    clabel(x,y)
    pause(0.03)
    title_text1 = sprintf('transient state Explicit using GS = %d', iter_2);
    title_text2 = sprintf('Simulation time  +%f s', time2);
    title_text3 = sprintf('Final error = %f' , error);
    title({title_text1; title_text2;title_text3})
end

%% Successive over relaxation (SOR) method
if solver ==3
    iter_3 = 1;
    tic
    for m = 1:nt
        while (error>tolerance)
            for i = 2:nx-1
                for j = 2:ny-1
                    %TL = t_old(i-1,j);
                    TR = t_old(i+1,j);
                    %TB = t_old(i,j-1);
                    TT = t_old(i,j+1);
                    H = t(i-1,j) +TR;
                    V = t(i,j-1) +TT;
                    tgs(i,j) = t_p(i,j)*k +H*k1 +V*k2;
                    t(i,j) = (1-omega)*(t_old(i,j))+omega*tgs(i,j);
                end
            end
            error = max(max(abs(t-t_old)));
            t_old = t;
            iter_3 = iter_3 +1;
        end
        t_p = t;
    end
    time3 = toc;
    [x,y]= contourf(xx,yy,t);
    colormap(jet)
    colorbar
    clabel(x,y)
    pause(0.03)
    title_text1 = sprintf('transient state Explicit using SOR = %d', iter_3);
    title_text2 = sprintf('Simulation time  +%f s', time3);
    title_text3 = sprintf('Final error = %f' , error);
    title({title_text1; title_text2;title_text3})
end

OUTPUT:

Jacobi

 

Gauss-seidel

 

Successive over-relaxation

For Explicit Scheme, we cannot use the same time steps and length of mesh domain, as it yields to higher CFL number, which leads to unstable solution unconverging. For the convenient I have taken the time steps and length domain for both x and y axis mentioned in the code section.  

Results and Observations: 

Comparing with the Different iterations approaches and states, we can conclude the following things 

The computational time and iteration required are SOR for both steady and unsteady states.  Steady states are less complicated to analyze because time marching properties change has been avoided, on the other hand Transient approach deals with more sophisticated and complicated unknown and consideration of time marching phenomena changes such as diffusion rates.  Between Explicit and Implicit results, the latter one does not depend greatly on Courant number for its stability and convergence, on the other hand in the Explicit approach information from a given cell or mesh element must propagate to its immediate neighbors, hence the low CFL number needs to be considered.