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

clear all
close all
clc

% initial conditions
L = 1;
nx = 10;   % no. of grid pts. along x direction
ny = 10;  % no. of grid pts. along y direction
x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);

% BC
T = ones(nx,ny)*298;   %temp
T(1,:)= 600;
T(:,1)= 400;
T(end,:) = 900;
T(:,end) = 800;
T(1,1) = 500;
T(end,1) =650;
T(1,end) = 700;
T(end,end)= 850;
error = 1;
tol= 1e-4;
Told = T;
aplha = 1.4;


%  JACOBI

jacobi_iter = 0;
c = (dx^2)*(dy^2)/(2*((dx^2)+(dy^2)));
while(error>tol)
    jacobi_iter = jacobi_iter+1  ;
    for i = 2:nx-1
        for j = 2:ny-1
            T(i,j) = ((Told(i-1,j)+Told(i+1,j))/(dx^2))+((Told(i,j-1)+Told(i,j+1))/(dy^2));
            T(i,j) = T(i,j)*c;
        end
    end
    error = max(max(abs(Told-T)));
    Told = T;

end

figure(1)
contourf(T,'ShowText', 'ON')
set(gca, 'ydir', 'reverse')
colormap(jet)
title_text = sprintf('no. of iterations = %d', jacobi_iter);
title(title_text)
xlabel('X-Axis')
ylabel('Y-Axis')
colorbar
legend('2D temeprature profile- Jacobi method')

% GAUSS SEIDEL METHOD

T = ones(nx,ny)*298;
T(1,:)= 600;
T(:,1)= 400;
T(end,:) = 900;
T(:,end) = 800;
T(1,1) = 500;
T(end,1) =650;
T(1,end) = 700;
T(end,end)= 850;
tol= 1e-4;
Told = T;
gs_iter = 0;
error = 1;
while(error>tol)
    gs_iter = gs_iter+1  ;       % over relaxation factor
    for i = 2:nx-1
        for j = 2:ny-1
            T(i,j) = ((T(i-1,j)+T(i+1,j))/(dx^2))+(T(i,j-1)+T(i,j+1))/(dy^2);
            T(i,j) = T(i,j)*c;
        end
    end
    error = max(max(abs(Told-T)));
    Told = T;

end
figure(2)
contourf(T,'ShowText', 'ON')
set(gca, 'ydir', 'reverse')
colormap(jet)
title_text = sprintf('no. of iterations = %d', gs_iter);
title(title_text)
xlabel('X-Axis')
ylabel('Y-Axis')
colorbar
legend('2D temeprature profile- gauss seidel method')


% Successive over-relaxation

c = (dx^2)*(dy^2)/(2*((dx^2)+(dy^2)));
T = ones(nx,ny)*298;  %temp
T(1,:)= 600;
T(:,1)= 400;
T(end,:) = 900;
T(:,end) = 800;
T(1,1) = 500;
T(end,1) = 650;
T(1,end) = 700;
T(end,end)= 850;
tol= 1e-4;
Told = T;
sor_iter = 0;
error = 1;
omega = 1.2;
while(error>tol)
    for i = 2:nx-1
        for j = 2:ny-1
            T(i,j) = (1-omega)*Told(i,j)+omega*c*(((T(i-1,j)+T(i+1,j))/(dx^2))+((T(i,j-1)+T(i,j+1))/(dy^2)));

        end
    end
    error = max(max(abs(Told-T)));
    Told = T;
    sor_iter = sor_iter+1;

end
figure(3)
contourf(T,'ShowText', 'ON')
set(gca, 'ydir', 'reverse')
colormap(jet)
title_text = sprintf('no. of iterations = %d', sor_iter);
title(title_text)
xlabel('X-Axis')
ylabel('Y-Axis')
colorbar
legend('2D temeprature profile - Successive over-relaxation method')


% Unsteady State Solution

dt = 0.01;
alpha = 1.4;


% EXPLICIT METHOD
% JACOBI
L = 1;
nx = 10;
ny = 10;
x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);
T = ones(nx,ny)*298;
T(:,1)= 400;
T(end,:) = 900;
T(:,end) = 800;
T(1,1) = 500;
T(end,1) =650;
T(1,end) = 700;
T(end,end)= 850;
error = 1;
tol= 1e-4;
Told = T;
T_prev_dt = Told;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

jacobi_iter = 1;
for nt = 1:1400
    error = 1;
    while (error>tol)
        for  i = 2:(nx-1)
            for j = 2:(ny-1)
                term1 = 1/(1+(2*k1)+(2*k2));
                term2 = k1*term1;
                term3 = k2*term1;
                H = (Told(i-1,j)+Told(i+1,j));
                V = (Told(i,j-1)+Told(i,j+1));
                T(i,j) = ((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        jacobi_iter = jacobi_iter+1;

    end
    T_prev_dt = T;

end
figure(4)
contourf(T,'ShowText', 'ON')
set(gca, 'ydir', 'reverse')
colormap(jet)
title_text = sprintf('no. of iterations = %d', jacobi_iter);
title(title_text)
xlabel('X-Axis')
ylabel('Y-Axis')
colorbar
legend('2D temeprature profile -EXPLICIT Jacobi method')


% Gauss seidel Method
L = 1;
nx = 10;
ny = 10;
x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);
T = ones(nx,ny)*298;
T(1,:)= 600;
T(:,1)= 400;
T(end,:) = 900;
T(:,end) = 800;
T(1,1) = 500;
T(end,1) =650;
T(1,end) = 700;
T(end,end)= 850;
error = 1;
tol= 1e-4;
Told = T;
gs_iter = 1;
for nt = 1:1400
    error = 1;
    while (error>tol)
        for  i = 2:(nx-1)
            for j = 2:(ny-1)
                term1 = 1/(1+(2*k1)+(2*k2));
                term2 = k1*term1;
                term3 = k2*term1;
                H = (Told(i-1,j)+Told(i+1,j));
                V = (Told(i,j-1)+Told(i,j+1));
                T(i,j) = ((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
            end
        end
        error = max(max(abs(Told-T)));
        gs_iter = gs_iter+1;
        Told = T;

    end
    T_prev_dt = T;
end
figure(5)
contourf(T,'ShowText', 'ON')
set(gca, 'ydir', 'reverse')
colormap(jet)
title_text = sprintf('no. of iterations = %d',gs_iter);
title(title_text)
xlabel('X-Axis')
ylabel('Y-Axis')
colorbar
legend('2D temeprature profile - EXPLICIT gauss seidel method')

% Successive over-relaxation
L = 1;
nx = 10;
ny = 10;
x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);
T = ones(nx,ny)*298;
T(1,:)= 600;
T(:,1)= 400;
T(end,:) = 900;
T(:,end) = 800;
T(1,1) = 500;
T(end,1) =650;
T(1,end) = 700;
T(end,end)= 850;
error = 1;
tol= 1e-4;
Told = T;
alpha = 10;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
omega = 1.05;
sor_iter = 1;
for nt = 1:1400
    error = 1;
    while (error>tol)
        for  i = 2:(nx-1)
            for j = 2:(ny-1)
                term1 = 1/(1+(2*k1)+(2*k2));
                term2 = k1*term1;
                term3 = k2*term1;
                H = (T(i-1,j)+Told(i+1,j));
                V = (T(i,j-1)+Told(i,j+1));
                T(i,j) =Told(i,j)*(1-omega)+(omega*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3)));
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        sor_iter = sor_iter+1;

    end
end
figure(6)
contourf(T,'ShowText', 'ON')
set(gca, 'ydir', 'reverse')
colormap(jet)
title_text = sprintf('no. of iterations = %d', sor_iter);
title(title_text)
xlabel('X-Axis')
ylabel('Y-Axis')
colorbar
legend('2D temeprature profile- Successive over-relaxation method')




%     Transient state solution with explicit

L = 1;
nx = 10;   % no. of grid pts. along x direction
ny = 10;  % no. of grid pts. along y direction
x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
alpha = 1.4;
dy = y(2)-y(1);
dt = 0.001;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

T = ones(nx,ny)*298;   %temp
T(1,:)= 600;
T(:,1)= 400;
T(end,:) = 900;
T(:,end) = 800;
T(1,1) = 500;
T(end,1) =650;
T(1,end) = 700;
T(end,end)= 850;
error = 1;
tol= 1e-4;
Told = T;
explicit_iter=1;
for nt = 1:1400
    for  i = 2:(nx-1)
        for j = 2:(ny-1)
            term1 =(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
            term2 =(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
            T(i,j) =((Told(i,j)+(k1*term1)+(k2*term2)));
        end
    end

    Told = T;
    explicit_iter = explicit_iter+1;
    T_prev_dt = T;
end
figure(7)
contourf(T,'ShowText', 'ON');
set(gca, 'ydir', 'reverse')
colormap(jet)
title_text = sprintf('simulation time = %d',dt*nt);
title(title_text);
xlabel('X-Axis')
ylabel('Y-Axis')
colorbar
legend('Transient state heat conduction');

Governing Equation -

2D conduction equation in cartesion coordinates is given below -

d^2T/dx^2 + d^2/dy^2 = 0 steady state

d^2T/dx^2 + d^2/dy^2 = ( 1/aplha)dT/dt unsteady state

Descretisation for steady state --

d^2T/dx^2 =( T L - 2T P+T R)/dx^2 + (T T -2T P + T B)/dy^2 = 0

T P = T i,j

T L = T i-1,j

T R = T i+1,j

T T = T i,j+1

T B = T i,j-1

Descretisation for unsteady state --

( T L - 2T P+T R)/dx^2 + (T T -2T P + T B)/dy^2 = ( 1/aplha)dT/dt

Obtain temperature distribution for a thin square plate by solving 2D heat conduction equation. To solve for 2D heat conduction equation for steady-state and unsteady-state by explicit and implicit method with given boundry conditions.

By implicit method

Solve using different iterative solvers i.e Jacobi,Gauss-seidal and successive over relaxation(SOR) for both steady and unsteady(transient) state heat conduction equation.

The heat equation is a partial differential equation that describes how the distribution of heat envolves over time in a solid medium , as it spontaneously flows from places where it is higher towards places where it is lower.

In STEADY_STATE condition the amount of heat entering any region of an object is equal to the amount of heat coming out, Steady-state is where the fluid properties does not change with respect to time.

Jacobi iterative solver --

In this method the tenperature at n+1th iteration is found using old values of the temperatures defined at nth iteration. Using this approximation the iteration procedure goinon continue untill the desire accuracy of the solution is obtained.

Gauss seidel iteration solver --

In this type of iteration ,temperature at n+1 iteration used values of n+1 iteration that have already been computed. by this way old approximated values are overwritten with new approximation values for the iterations performed. unlike jacobi, gauss seidal require one storage vector which can be overwritten with each computation. This helps for the solution to coverage faster than the jacobi method.

Successive over relaxation iterative solver ---

In this type of iterative solver relaxation factor is used which helps in increasing the convergence rate

Steady state

in steay state the amount of heat entering any region of an object is equal to the amount of heat coming out

Explicit and impicit methods are approaches used in numerical analysis for obtaining numerical approximations to the solution of time-dependent ordinary and partial differential equations , as it required in computer simulations of physical processes. Explicit method calculate the state of a system at a lter time from the state of the system at the current time, while implicit methods find a solution by solving an equation involving both the current state of the system and the later one.

Coding procedure ----

1. Compute Area is to be a square

2. Assume nx = ny (no. of points in x-direction is equal to no. of points in y-direction)

3. Assume nx = ny = 10

4. Boundary conditions are

Top = 600k , Bottom = 900k , Left = 400k , Right = 800k

5. Various method to be employed in solving the equation are given below-

(a). VARIOUS STATES : steady & unsteady-state conditions

(b). VARIOUS SCHEMES : jacobi, Gauss-seidal & Successive over-relaxation(SOR)

(c). VARIOUS METHODS TO BE USED IN UNSTEADY-STATE : Implicit & Explicit

6. After initializing the input parameters calculate average temperature at corners of the plate.

7. Assigning originals Values to T because we are going to calculate heat conduction equation for jacobi,gauss seidal and successive over relaxation.

8. Putting the values of k1 , k2 and error & tolerance value and start exicuting while loop and within the while loop apply for loop for the x and y direction and inside the while loop obtained a value of maximum error in the temperature and update the value of Told as a T and closed the loop .

9. Now, plotting a graph for jacobi giving a ' contourf ' command used for display the contour lines within the graph. and 'set' command for reverse the parameters of temperature in the graph. 'colormap(jet)' command for showing the best specific color in the figure for better analysis the color for different temperature range. 'title_text = sprintf' command for showing the no. of iteration used in the specific iterative solver.'colorbar' command used for showing the bar for different temperature range and enclosed with the legend of that specific iterative solver .

in addition it it, the same fashion goes to SOR and gauss seidel method for steady state. in the case of unsteady(transient) state in Explicit we have updating the value of T_prev = T and using a omega = 1.4(relaxation factor) in the SOR for both steady and unsteady case.

Solving Transient state Equation-

Transient state 2D steady state Heat conduction equation is solved in both ways i.e explicitly and implicitly.

Solvng this equation explicitly - expilicit methods calculate the state of a system at a later time from the state of the sytem at the current time.

In this method, each differential equation approach contains only one unknown and can be explicitly for this one unknown using this method. Here in this equation time derivative terms are calculated using the first order differential equation and whereas the term 2 has second order derivative terms and elliptical in nature so information needs to be taken from all the points equally, So, therefore a central differencing scheme is applied for the discretising the term.

Conclusion --

1. when comapring steady state and unsteady state ,it is proven that steady state is the most successful one as it takes least no. of iterations and least time to compute.

2. SOR iteration method is the most effective iteration method as it took least no. of iteration of covering solution in both the schemes(EXPLICT & IMPLICIT) comapred to the jacobi and Gauss-seidel.