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

`del T/delt = alpha*(del^2T/del x + del^2T/del y + del^2T/del z)`

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.

Transient State of equation: for 2D)

`del T/delt = alpha*(del^2T/del x + del^2T/del y)`

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

`(del^2T/del x + del^2T/del y) =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)/(Delta x^2) + (T_B- 2T +T_T)/(Delta y^2) =0`

On simplifying,

assume, `k = 2(Delta x^2 +Delta y^2)/(Delta x^2*Delta y^2`

we can write

`T_p = 1/k((T_L +T_R)/(Delta x^2) + (T_B +T_T)/(Delta y^2))`

MATLAB SCRIPT FOR STEADY STATE:

clear
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

Plots:

As from stated on the figure, the number of iteration required to converge the solution is more in Jacobi iteration scheme followed by Gauss Siedle iteraion and SOR scheme. Which obviously means the computaional time required for a computer program to run the iteration is more in Jacobi followed by GS and SOR scheme.

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)`

MATLAB SCRIPT

clear
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

Plotas:

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;

`del T/ del t= alpha nabla^2 T`

`del T/ del t= (T^(n+1) -T^n)/(Delta t)`

`T^(n+1) = T^n + alpha Delta t ((T_L -2T^n +T_R)/(Delta x^2) + (T_B -2T^n +T_T)/(Delta y^2))`

Assuming `k1 = alpha (Delta t)/(Delta x^2) , k2 = alpha (Delta t)/(Delta y^2)`

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

MATLAB SCRIPT:

clear
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

Plots:

For Explicit Scheme, we cannot use the same time steps and length of mesh domain, as it yeids 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 iteraions approaches and states, we can conclude the following things

• The computational time and iteration required are SOR<GS<JACOBI for both steady and unsteady states.
• Steady states are less complicated to analyse 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 phenomen 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 neighbours, hence the low CFL number needs to be considered.

Thank you!