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

Steady state Analysis:

the 2D heat conduction equation for steady state is :

`(del^2T)/(del^2x)+(del^2T)/(del^2y) = 0`

now discretizing this equation by central differencing

equation wrt x becomes as

`(del^2T)/(del^2x) = (T(i-1,j)-2*T(i,j)+T(i+1,j))/(Deltax^2)`

equation wrt y becomes as

`(del^2T)/(del^2y)= (T(i,j-1)-2*T(i,j)+T(i,j+1))/(Deltay^2)`

hence the final equation is

`T(i,j) = (T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1))*(1/(K*dx^2))`, dx = dy as nx = ny

where `k = 2*(dx^2+dy^2)/(dx^2*dy^2)`

CODE :

clear all
close all
clc
% steady state 

nx = 30;
ny = nx;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
term= 1/(k*dx^2); % as nx = ny ,dx = dy, 1/(k*dx^2) = 1/(k*dy^2)

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T+T_L)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,ny) = (T_B+T_R)/2;
T(nx,1) = (T_B+T_L)/2;
Told = T;
T_pre_val = Told;
[X Y] = meshgrid(x,y);
iterative_solver = 1:3;

% jacobi method
for iterative_solver = 1
    tic;
    if iterative_solver ==1
        jacobi_iter = 1;
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    H = Told(i-1,j)+Told(i+1,j);
                    V = Told(i,j-1)+Told(i,j+1);
                    T(i,j) = (H+V)*term;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            jacobi_iter = jacobi_iter+1;
        end
        F= toc;
    end
    %plot
    figure(1)
    [m n] = contourf(x,y,T);
    clabel(m,n);
    colormap(jet);
    colorbar;
    title(sprintf(\'steady state jacobi iteration number = %d Total time = %fs\',jacobi_iter,F));
    xlabel(\'X-axis\');
    ylabel(\'Y-axis\');
    fprintf(\'no of iterations = %d\',jacobi_iter);
end
% gauss siedel method
nx = 30;
ny = nx;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
term= 1/(k*dx^2); % as nx = ny ,dx = dy, 1/(k*dx^2) = 1/(k*dy^2)

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T+T_L)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,ny) = (T_B+T_R)/2;
T(nx,1) = (T_B+T_L)/2;
Told = T;
T_pre_val = Told;
[X Y] = meshgrid(x,y);

for iterative_solver = 2
    tic;
    if iterative_solver ==2
        gs_iter = 1;
        while(error>tol)
             for i = 2:nx-1
                 for j = 2:ny-1;
                     H = T(i-1,j)+T(i+1,j);
                     V = T(i,j-1)+T(i,j+1);
                     T(i,j) = (H+V)*term;
                 end
             end
             error = max(max(abs(Told -T)));
             Told = T;
             gs_iter = gs_iter+1;
        end
        F = toc;
    end
    % plot
    figure(2)
    [m n] = contourf(x,y,T);
    clabel(m,n);
    colormap(jet);
    colorbar;
    title(sprintf(\'steady state gauss siedel method iteration number = %d Total time = %fs\',gs_iter,F));
    xlabel(\'X-axis\');
    ylabel(\'Y-axis\');
    fprintf(\'no of iterations = %d\',gs_iter);
end
% sor method
nx = 30;
ny = nx;
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = linspace(0,1,ny);
dy = y(2)-y(1);
error = 9e9;
tol = 1e-4;
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
term= 1/(k*dx^2); % as nx = ny ,dx = dy, 1/(k*dx^2) = 1/(k*dy^2)

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T+T_L)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,ny) = (T_B+T_R)/2;
T(nx,1) = (T_B+T_L)/2;
Told = T;
T_pre_val = Told;
[X Y] = meshgrid(x,y);
alpha = 1.8;

for iterative_solver = 3
    tic;
    if iterative_solver == 3
        sor_iter = 1;
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    H = T(i-1,j)+T(i+1,j);
                    V = T(i,j-1)+T(i,j+1);
                    T(i,j) = (1-alpha)*Told(i,j)+alpha*(H+V)*term;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            sor_iter = sor_iter +1;
        end
        F= toc;
    end
    % plot
    figure(3)
    [m n] = contourf(x,y,T);
    clabel(m,n);
    colormap(jet);
    colorbar;
    title(sprintf(\'steady state sor method iterations = %d Total time = %fs\',sor_iter,F));
    xlabel(\'X-axis\');
    ylabel(\'Y-axis\');
    fprintf(\' no of iterations = %d\',sor_iter);
end



    
    

graphs/plots of steady state analysis:

for jacobi method:

for gauss-siedel method

for sor method

Result:

from the above graphs it is observed that SOR method requires least time to compute with minimum number of iterations , hence it is the fastest convergence method.

Unsteady state analysis :

the 2D heat conduction equation for unsteady state is

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

now time derivative is being discretied by forward differencing and the remaining terms by central differencing.

after simplification final equation for explicit solver is :

`T^(n+1)(i,j) = T^n(i,j)+k1*(T^n(i-1,j)-2*T^n(i,j)+T^n(i+1,j))+k2*(T^n(i,j-1)-2*T^n(i,j)+T^n(i,j+1))`

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

final equation for implicit solver is :

`T^(n+1)(i,j) = (T^n(i,j)+k3*(T^n(i-1,j)+T^n(i+1,j))+k4*(T^n(i,j-1)+T^n(i,j+1)))/(1+2*k3+2*k4)`

where `k3 = alpha*(dt)/(dx^2) , k4 = alpha*(dt)/(dy^2)`

Solving using explicit method:

code:

clear all
close all
clc


nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
error = 9e9;
tol = 1e-4;
dt = 10;
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.88e-5;
T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T+T_L)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,ny) = (T_B+T_R)/2;
T(nx,1) = (T_B+T_L)/2;

k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
[X Y] = meshgrid(x,y);
Told = T;
T_pre_val = Told;

for iterative_solver = 1
    tic;
    if iterative_solver ==1
        exp_iter = 1;
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    term1 = Told(i,j);
                    term2 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
                    term3 = k2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
                    T(i,j) = term1+term2+term3;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            exp_iter = exp_iter+1;
        end
        F = toc;
    end
    % plot 
    [M N] = contourf(x,y,T);
    clabel(M,N);
    colormap(jet);
    colorbar;
    title(sprintf(\'unsteady heat conduction by solving explicitly iteration number = %d Total time = %fs\',exp_iter,F));
    xlabel(\'X-axis\');
    ylabel(\'Y-sxis\');
    fprintf(\'no of iterations = %d\', exp_iter);
end

    
          

graphs for transient flow (explicit method):

Solving for transient flow (implicit method):

Code:

clear all
close all
clc

nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
error = 9e9;
tol = 1e-4;
dt = 2000;
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.88e-5;
T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T+T_L)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,ny) = (T_B+T_R)/2;
T(nx,1) = (T_B+T_L)/2;

k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
[X Y] = meshgrid(x,y);
Told = T;
T_pre_val = Told;

iterative_solver = 1:3;
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;


for iterative_solver = 1
    tic;
    if iterative_solver == 1
        jacobi_iter = 1;
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    H = (Told(i-1,j)+Told(i+1,j));
                    V = (Told(i,j-1)+Told(i,j+1));
                    T(i,j) = term1*T_pre_val(i,j)+H*term2+V*term3;
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            jacobi_iter = jacobi_iter+1;
        end
    end
   F = toc;
   % plot
   figure(1)
   [M N] = contourf(x,y,T);
   clabel(M,N);
   colormap(jet);
   colorbar;
   title(sprintf(\'jacobi unsteady iteration number = %d, Total time = %fs\', jacobi_iter,F));
   xlabel(\'X-axis\');
   ylabel(\'Y-axis\');
   fprintf(\'no of iterations is = %d\',jacobi_iter);
end

% updated values
nx= 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
error = 9e9;
tol = 1e-4;
dt = 2000;
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.88e-5;
T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T+T_L)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,ny) = (T_B+T_R)/2;
T(nx,1) = (T_B+T_L)/2;

k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
[X Y] = meshgrid(x,y);
Told = T;
T_pre_val = Told;

iterative_solver = 1:3;
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;

% gauss siedel method
for iterative_solver = 2;
    tic;
    if iterative_solver == 2
        gs_iter = 1;
        while(error>tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    H = T(i-1,j)+T(i+1,j);
                    V = T(i,j-1)+T(i,j+1);
                    T(i,j) = T_pre_val(i,j)*term1+H*term2+V*term3;
                end
            end
            error = max(max(abs(Told - T)));
            Told = T;
            gs_iter = gs_iter+1;
        end
        F = toc;
    end
    % plot 
    figure(2)
    [M N] = contourf(x,y,T);
    clabel(M,N);
    colormap(jet);
    colorbar;
    title(sprintf(\'gauss siedel iteration number = %d, Total time = %fs\',gs_iter,F));
    xlabel(\'X-axis\');
    ylabel(\'Y-axis\');
    fprintf(\'no of iterations = %d\',gs_iter);
end


nx = 30;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
error = 9e9;
tol = 1e-4;
dt = 2000;
dx = x(2)-x(1);
dy = y(2)-y(1);
alpha = 1.88e-5;
T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T+T_L)/2;
T(1,ny) = (T_T+T_R)/2;
T(nx,ny) = (T_B+T_R)/2;
T(nx,1) = (T_B+T_L)/2;

k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
[X Y] = meshgrid(x,y);
Told = T;
T_pre_val = Told;

iterative_solver = 1:3;
term1 = (1+2*k1+2*k2)^(-1);
term2 = k1*term1;
term3 = k2*term1;
alpha1 = 1.8;

for iterative_solver = 3
    tic;
    if iterative_solver == 3
        sor_iter = 1
        while (error>tol)
             for i = 2:nx-1
                 for j = 2:ny-1
                     H = alpha1*T(i-1,j)+Told(i+1,j)+(1-alpha1)*Told(i-1,j);
                     V = alpha1*T(i,j-1)+Told(i,j+1)+(1-alpha1)*Told(i,j-1);
                     T(i,j) = T_pre_val(i,j)*term1+H*term2+V*term3;
                 end
             end
             error = max(max(abs(Told - T)));
             Told = T;
             sor_iter = sor_iter +1;
        end
        F = toc;
    end
    % plot
    figure(3)
    [M N] = contourf(x,y,T);
    clabel(M,N);
    colormap(jet);
    colorbar;
    title(sprintf(\'sor unsteady iteration number = %d, Total time = %fs\',sor_iter,F));
    xlabel(\'X-axis\');
    ylabel(\'Y-axis\');
    fprintf(\'no of iterations = %d\',sor_iter);
end

graphs/plots for transient flow implicit method :

for jacobi method:

for gauss method:

for sor method :

Result:

even in this case sor method is taking least time to compute and that to with minimum number of iterations.So SOR method is cost and time effective method when compared to jacobi and gauss seidel method.And the accuracy is much higher for implicit methods when comared to explicit.