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

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

Assumptions:

1. The domain is unit square.

2. nx=ny, no of grid points along the x-direction is equal to the no of grid points along the y- direction.

3. Temperature specify on the boundaries are

(i) left 400K (ii) top 600K (iii) right 800K (iv) bottom 900K

1.Transient Case:

When a direct computation of the dependent variables can be made in terms of known quantities, the computation is said to be explicit. When the dependent variables are defined by coupled sets of equation, and either a matrix or iterative technique is needed to obtain the solution, the numerical method is said to be implicit.

1.1 Solving the transient 2D heat conduction equation explicitly:

The governing equation for 2D heat conduction equation is

`frac{partialT}{partialt}=alpha(frac{partial^2T}{partialx^2}+frac{partial^2T}{partialy^2})`

Discretize the right hand side by using central differencing scheme and left hand side by time marching scheme.

`frac{T_(i,j)^(n+1)-T_(i,j)^n}{Deltat}= alpha(frac{T_(i-1,j)-2T_(i,j)+T_(i+1,j)}{Deltax^2}+frac{T_(i,j-1)-2T_(i,j)+T_(i,j+1)}{Deltay^2})^n`

rearrange the above equation

`T_(i,j)^(n+1)= T_(i,j)^n+alphaDeltat(frac{T_(i-1,j)-2T_(i,j)+T_(i+1,j)}{Deltax^2}+frac{T_(i,j-1)-2T_(i,j)+T_(i,j+1)}{Deltay^2})^n`

again rearrange

`T_(i,j)^(n+1)= T_(i,j)^n+k_1(T_(i-1,j)-2T_(i,j)+T_(i+1,j))^n+k_2(T_(i,j-1)-2T_(i,j)+T_(i,j+1))^n`

where `k_1=frac{alpha deltat}{deltax^2}` `k2= frac{alpha delta t}{delta y^2}`

1.2 Solving the transient 2D heat conduction equation implicitly:

`frac{partialT}{partialt}=alpha(frac{partial^2T}{partialx^2}+frac{partial^2T}{partialy^2})`

After discretization

`frac{T_(i,j)^(n+1)-T_(i,j)^n}{Deltat}= alpha(frac{T_(i-1,j)-2T_(i,j)+T_(i+1,j)}{Deltax^2}+frac{T_(i,j-1)-2T_(i,j)+T_(i,j+1)}{Deltay^2})^(n+1)`

rearrange the above equation, k1 and k2 are defined above.

`T_(i,j)^(n+1)= T_(i,j)^n+k_1(T_(i-1,j)-2T_(i,j)+T_(i+1,j))^(n+1)+k_2(T_(i,j-1)-2T_(i,j)+T_(i,j+1))^(n+1)`

Rearrange the above terms

`T_(i,j)^(n+1)+2k_1T_(i,j)^(n+1)+2k_2T_(i,j)^(n+1)= T_(i,j)^(n)+k_1(T_(i-1,j)+T_(i+1,j))^(n+1)+k_2(T_(i.j-1)+T_(i,j+1))^(n+1)`

`T_(i,j)^(n+1)(1+2k_1+2k_2)= T_(i.j)^(n)+k_1(T_(i-1,j)+T_(i+1,j))^(n+1)+k_2(T_(i.j-1)+T_(i,j+1))^(n+1)`

`T_(i,j)^(n+1)=frac{ T_(i,j)^(n)+k_1(T_(i-1,j)+T_(i+1,j))^(n+1)+k_2(T_(i.j-1)+T_(i,j+1))^(n+1)}{1+2k_1+2k_2}`

The above equation is solved using jacobi, Gauss seidel, and successive over relaxation iterative methods.

2. Steady state case:

`frac{partial^2T}{partialx^2}+frac{partial^2T}{partialy^2}=0`

`frac{T_(i-1,j)-2T_(i,j)+T_(i+1,j)}{Deltax^2}+frac{T_(i,j-1)-2T_(i,j)+T_(i,j+1)}{Deltay^2}=0`

since `Deltax=Deltay`in this case then the above equation becomes

`T_(i,j)=frac{T_(i-1,j)+T_(i+1,j)+T_(i,j-1)+T_(i,j+1)}{4}`

This is solved by using Jacobi, Gauss seidel and successive over relaxation iterative technique.

Main Matlab Code:

% Mid term project 2D heat conduction Solver- Transient & steady state solver
clc
clear
clc
global nx ny tol nt k1 k2 omega %declaring global variables to be used in all functions
L=1;
nx = 10;     % no of grid points along x-axis
ny=nx;       %no of grid points along y-axis
tol = 1e-4;  %Tolerance
T = ones(nx,ny);  %initializing Temperature matrix
% Assigning boundary conditions
T_T = 600;  %top edge
T_B = 900;  %bottom edge
T_L = 400;   %left edge
T_R = 800;   %right edge
      % corner adjustment
T(1,1)=(900+400)/2;  %corner left and bottom
T(nx,1)= (800+900)/2;  %corner right and bottom
T(nx,ny)= (800+600)/2;  %corner right and top
T(1,ny) = (600+400)/2;  %corner top and left
T(2:nx-1,ny)= T_T;
T(2:nx-1,:)= T_B;
T(1,2:ny-1)= T_L;
T(nx,2:ny-1)= T_R;
%alpha = 1.9e-5;  %diffusivity constant
alpha = 1.4;
dt = 1/1000;      %time step
nt=800;            
omega = 1.5;
x = linspace(0,L,nx);
y= linspace(0,L,ny);
dx = x(2)-x(1);
dy = y(2)-y(1);
T_old =T;  %Define a matrix to contain old values of T
T_prev_dt = T_old; %Define a matrix to contain a values of T at a previous time step
% Defining constants to be used in the equation
k1=  (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
                     %Transient Analysis
 for i=1:4  % loop for calling different functions and plotting the solutions
     if(i==1)
         figure(i)
         [time_taken,T,iter] = Jacobi(T,T_old,T_prev_dt); % solve by jacobi implicit
     end
     if(i==2)
         figure(i)
         [time_taken,T, iter] = Gaussseidel(T,T_old,T_prev_dt); % solve by gauss seidel implicit
     end
     if(i==3)
         figure(i)
         [time_taken,T, iter] = SOR(T,T_old,T_prev_dt); % solve by SOR implicit
     end
         if(i==4)
         figure(i)
         [time_taken,T, iter] = Explicit(T,T_old); % solve by Explicit method
         end
     [M,N]=contourf(x,y,T);   %create a filled contour plot
     clabel(M,N);           %label the contour
     colormap(jet);         %set the colormap to jet
     colorbar;              %Display the colorbar
    xlabel(\'X-axis\')
    ylabel(\'Y-axis\')
    title_1 = sprintf(\'Number of Iterations = %d\',iter);
    title_2 = sprintf(\'Simulation Time = %d\',time_taken);
    title({\'Transient scheme\',;title_1;title_2});
        end
                    %steady state Analysis
    for i=5:7  % loop for calling different functions and plotting the solutions
     if(i==5)
         figure(i)
         [time_taken,T,iter] = JacobiSS(T,T_old); % solve by jacobi explicit
     end
     if(i==6)
         figure(i)
         [time_taken,T, iter] = SeidelSS(T,T_old); % solve by gauss seidel explicit
     end
     if(i==7)
         figure(i)
         [time_taken,T, iter] = SORSS(T,T_old); % solve by SOR implicit
     end
     [M,N]=contourf(x,y,T); %create a filled contour plot
     clabel(M,N);           %label the contour
     colormap(jet);         %set the colormap to jet
     colorbar;              %Display the colorbar
     xlabel(\'X-axis\')
    ylabel(\'Y-axis\')
    title_1 = sprintf(\'Number of Iterations = %d\',iter);
    title_2 = sprintf(\'Simulation Time = %d\',time_taken);
    %title({title_1;title_2});
    title({\'Steady state scheme\',;title_1;title_2});
    end

Functions codes are displayed below which is called in above program

Jacobi function code for transient implicit scheme:

function[time_taken1,T, jacobi_iter] = Jacobi(T,T_old,T_prev_dt)
global nx ny tol nt k1 k2 
tic;
jacobi_iter = 1;  %initializing the counter variable
for time = 1:nt  %time loop
    error=9e9;   %initializing the error
    while(error>tol)  %convergence loop
        for i=2:nx-1     %nodal loop
            for j=2:ny -1 %nodal loop
                term1 = (1+2*k1+2*k2)^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                H = T_old(i-1,j)+T_old(i+1,j);
                V = T_old(i,j-1)+T_old(i,j+1);
                T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
            end
        end
         error = max(max(abs(T_old-T)));  %calculate error
    T_old = T;  %update old values
    jacobi_iter = jacobi_iter+1;   %increment counter
    end
    T_prev_dt=T;   %update previous time step value
end
toc;
time_taken1= toc;
end
 

Gauss Seidel function code for transient implicit scheme:

function[time_taken2,T,GS_iter] = Gaussseidel(T,T_old,T_prev_dt)
global nx ny tol nt k1 k2 
tic;
GS_iter = 1;  %initializing the counter variable
for time = 1:nt  %time loop
    error=9e9;   %initializing the error
    while(error>tol)  %convergence loop
        for i=2:nx-1     %nodal loop
            for j=2:ny -1 %nodal loop
                term1 = (1+2*k1+2*k2)^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                H = T(i-1,j)+T_old(i+1,j);
                V = T(i,j-1)+T_old(i,j+1);
                T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
            end
        end
         error = max(max(abs(T_old-T)));  %calculate error
    T_old = T;  %update old values
    GS_iter = GS_iter+1;   %increment counter
    end
    T_prev_dt=T;   %update previous time step value
end
toc;
time_taken2= toc;
end

SOR function code for transient implicit scheme:

function[time_taken3,T,SOR_iter] = SOR(T,T_old,T_prev_dt)
global nx ny tol nt k1 k2 omega
tic;
SOR_iter = 1;  %initializing the counter variable
for time = 1:nt  %time loop
    error=9e9;   %initializing the error
    while(error>tol)  %convergence loop
        for i=2:nx-1     %nodal loop
            for j=2:ny -1 %nodal loop
                term1 = (1+2*k1+2*k2)^-1;
                term2 = k1*term1;
                term3 = k2*term1;
                H = T(i-1,j)+T_old(i+1,j);
                V = T(i,j-1)+T_old(i,j+1);
                T(i,j) = (T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
                T(i,j) = ((T_old(i,j))*(1-omega))+omega*T(i,j);
            end
        end
         error = max(max(abs(T_old-T)));  %calculate error
    T_old = T;  %update old values
    SOR_iter = SOR_iter+1;   %increment counter
    end
    T_prev_dt=T;   %update previous time step value
end
toc;
time_taken3= toc;
end

Results discussion for implicit schemes:

Jacobi contours display for transient implicit scheme:

Gauss Seidel Contours display for transient implicit scheme:

SOR Contours display for transient implicit scheme:

Above results discussion of implicit schemes:

As we can see in above results that we can obtained less number of iterations in the case of SOR implicit scheme as compared to the Jacobi and Gauss seidel.

no of iterations obtained in Jacobi method:2974

no of iterations obtained in Gauss seidel method:2402

no of iterations obtained in SOR method:2059

Bascially the range of omega lies between 0 and 2. At the value of omega =1, Gauss seidel and SOR obtained the same number of iterations.

For different values of omega the sor number of iterations results are given below.

for omega =1, no of iterations in SOR implicit=2402

for omega=1.1, no of iterations in SOR implicit=2059

for omega=1.2, no of iterations in SOR implicit=2073

for omega=1.3, no of iterations in SOR implicit=2526

for omega=1.4, no of iterations in SOR implicit=3154

for omega=1.5, no of iterations in SOR implicit=4045

for omega=1.6, no of iterations in SOR implicit=5433

for omega=1.7, no of iterations in SOR implicit=7949

for omega =1.8, no of iterations in SOR implicit=13516

for omega=1.9, no of iterations in SOR implicit=26836

in the above contours we see that the red colour shows that the heat is diffused slowly to the place where the temperature is not too much and blue colour shows that the temperature is not too much high at this place. in this rectangular plate corner temperaure is adjust by taking the average of temperatures. in the above graphs behaviour of solutions are correct.

Function code for transient explicit schemes:

function  [time_taken4,T,counter] = Explicit(T,T_old)
global nx ny tol nt k1 k2 
tic;
counter = 0;  %initializing the counter variable
for time = 1:nt  %time loop
    error=9e9;   %initializing the error
    while(error>tol)  %convergence loop
        for i=2:nx-1     %nodal loop
            for j=2:ny -1 %nodal loop
                H = T_old(i+1,j)+T_old(i-1,j)-2*T_old(i,j);
                V = T_old(i,j-1)+T_old(i,j+1)-2*T_old(i,j);
                T(i,j) = T_old(i,j)+(k1*H)+(k2*V);
            end
        end
         error = max(max(abs(T_old-T)));  %calculate error
    T_old = T;  %update old values
    counter = counter+1;   %increment counter
    end
   toc;
time_taken4= toc;
end

Contours display for transient explicit scheme:

In explicit case there is no need to apply iterative scheme because we have only one unknown and all other values are known.so we need to solve only one equation.

Jacobi function for steady state scheme:

function  [time_taken5,T, jacobi_iter] = JacobiSS(T,T_old)
global nx ny tol 
tic;
jacobi_iter = 1;  %initializing the counter variable
    error=9e9;   %initializing the error
    while(error>tol)  %convergence loop
        for i=2:nx-1     %nodal loop
            for j=2:ny -1 %nodal loop
                H = T_old(i-1,j)+T_old(i+1,j);
                V = T_old(i,j-1)+T_old(i,j+1);
                T(i,j) = 0.25*(H+V);
            end
        end
         error = max(max(abs(T_old-T)));  %calculate error
    T_old = T;  %update the old values
    jacobi_iter = jacobi_iter+1;   %increment counter
    end
toc;
time_taken5= toc;
end

Gauss Seidel function code for steady state scheme:

function  [time_taken6,T, seidel_iter] = SeidelSS(T,T_old)
global nx ny tol 
tic;
seidel_iter = 1;  %initializing the counter variable
    error=9e9;   %initializing the error
    while(error>tol)  %convergence loop
        for i=2:nx-1     %nodal loop
            for j=2:ny -1 %nodal loop
                H = T(i-1,j)+T_old(i+1,j);
                V = T(i,j-1)+T_old(i,j+1);
                T(i,j) = 0.25*(H+V);
            end
        end
         error = max(max(abs(T_old-T)));  %calculate error
    T_old = T;  %update the old values
    seidel_iter = seidel_iter+1;   %increment counter
    end
toc;
time_taken6= toc;
end

SOR function code for steady state scheme:

function [time_taken7,T, sor_iter] = SORSS(T,T_old)
global nx ny tol omega
tic;
sor_iter = 1;  %initializing the counter variable
    error=9e9;   %initializing the error
    while(error>tol)  %convergence loop
        for i=2:nx-1     %nodal loop
            for j=2:ny -1 %nodal loop
                H = T(i-1,j)+T_old(i+1,j);
                V = T(i,j-1)+T_old(i,j+1);
                T(i,j) = 0.25*(H+V);
                T(i,j) = (T_old(i,j)*(1-omega))+(omega*T(i,j));
            end
        end
         error = max(max(abs(T_old-T)));  %calculate error
    T_old = T;  %update the old values
    sor_iter = sor_iter+1;   %increment counter
    end
toc;
time_taken7= toc;
end

Jacobi Contours display for steady state scheme:

Gauss seidel Contours display for steady state scheme:

SOR Contours display for steady state scheme:

Steady State results analysis:

In steady state scheme time derivative term is zero, so it requires less computational time. In SOR no of iterations is also too much less as compared to Jacobi and gauss seidel.

no of iterations in jacobi=235

no of iterations in Gauss seidel=135

no of iterations in SOR=103

To get the optimal value of omega i do experiment for differnent values of omega, whose details are given below:

at omega=1, no of iterations in seidel is equal to the no of iterations in SOR

omega=1.1, no of iterations in SOR=103

omega=1.2, no of iterations in SOR=84

omega=1.3, no of oterations in SOR=67

omega=1.4, no of iterations in SOR=52

omwga=1.5, no of iterations in SOR=36

omega=1.6, no of iterations in SOR=33

omega=1.7, no of iterations in SOR=44

omega=1.8, no od iterations in SOR=70

omega =1.9, no of iterations in SOR=145