BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

SShivasharan Mendegar5 Jul 2021
BIM

main code :

%Solving the steady and unsteady 2D heat conduction problem using different
%iterative techniques 
close all
clear all
clc
%Defining global varibles used in all functions
global nx ny to1 nt k1 k2 omega
L = 1;
% no. of grid points along x- axis
nx = 15;
% no. of grid points along y- axis
ny = nx;
%Tolorence
to1 = 1e-4;
%Inatializing temperature matrix 
T = ones(nx,ny);
% Asigning Boundary Conditionns 
% Corners adjustments
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(1,1:nx) =600;
T(ny,1:nx) = 900;
T(2:ny-1,1) = 400;
T(2:ny-1,nx) = 800;

% Diffusivity constant
alpha = 1.4;
dt = 1/1000; % time steps
nt = 800;
omega = 1.5;
x = linspace(0,L,nx);
dx = x(2)-x(1);
y = linspace(0,L,ny);
dy = y(2)-y(1);

%defining matrix which contains old value of T
T_old = T;
%defining matrix to which contain old value of T at a previous time step
T_prev_dt = T_old;

%Defining Constants k1 nd k2
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

                    %TRANSIENT ANALYSIS
%for loop for calling different functions and plotting the solutions 
for i = 1:3       
    if (i==1)
        figure(i)
        [time_taken,T,iter] = Jacobi(T,T_old,T_prev_dt); %jacobi
    end
    if (i==2)
        figure(i)
        [time_taken,T,iter] = Gaussseidal(T,T_old,T_prev_dt);%gauss Seidel
    end
    if (i==3)
        figure(i)
        [time_taken,T,iter] = SOR_t(T,T_old,T_prev_dt); %SOR
    end
    
    
    [M,N]= contourf(x,y,T); %creating field contour plot 
    clabel(M,N)
    colormap(jet)           %set colormap to jet
    colorbar;
    xlabel('X- Axis')
    ylabel('Y- Axis')
    title_1 = sprintf('no. of iteration = %d',iter);
    title_2 = sprintf('Simulation Time = %d',time_taken);
    title({'Transient Scheme',;title_1;title_2});
    set (gca,'YDir','reverse')
end

                      %EXPLICIT ANALYSIS
                                  
    %CFL = alpha*(dt/dx^2);
    CFL=0.25;
    for i = 4
    if(i==4)
    figure(i)
    [time_taken,T,iter] = Explicit(T,T_old); %Explicit
    end
    
    [M,N]= contourf(x,y,T); %creating field contour plot 
    clabel(M,N)
    colormap(jet)           %set colormap to jet
    colorbar;
    xlabel('X- Axis')
    ylabel('Y- Axis')
    title_1 = sprintf('no. of iteration = %d',iter);
    title_2 = sprintf('Simulation Time = %d,CFL = %f',time_taken,CFL);
    title({'Transient Scheme',;title_1;title_2});
    set (gca,'YDir','reverse')
    end

                   %STEADY STATE ANALYSIS
   for i=5:7
       if (i==5)
        figure(i)
       [time_taken,T,iter] = JacobiSS(T,T_old);
       end
    if (i==6)
        figure(i)
        [time_taken,T,iter] = SeidelSS(T,T_old);
    end
       if (i==7)
        figure(i)
       [time_taken,T,iter] = SORSS(T,T_old);
       end
       
    [M,N]= contourf(x,y,T);
    clabel(M,N)
    colormap(jet)
    colorbar;
    xlabel('X- Axis')
    ylabel('Y- Axis')
    title_1 = sprintf('no. of iteration = %d',iter);
    title_2 = sprintf('Simulation Time = %d',time_taken);
    title({'Steady State Scheme',;title_1;title_2});
    set (gca,'YDir','reverse') 
   end

Function codes for TRANSIENT IMPLICIT METHOD

function code for jacobi iterative technique

%jaccob function code for transient implicit
function[time_taken,T,Jacobi_iter] = Jacobi(T_old,T,T_prev_dt)
global nx ny nt to1 k1 k2
tic;
Jacobi_iter = 1;          %inatializing the counter variable
for time = 1:nt            %Time loop
   error = 9e9;           %Inatializing the error
   while(error>to1)       %convergence loop
        for i = 2:nx-1     %nodal loop
            for j = 2:ny-1
                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)));
        T_old = T;
        Jacobi_iter = Jacobi_iter+1;
    end
    T_prev_dt = T;
end
toc;
time_taken = toc;
end

function code for gauss seidel 

%Gauss Seidel function code for transient implicit
function[time_taken2,T,GS_iter] = Gaussseidal(T,T_old,T_prev_dt)
global nx ny nt to1 k1 k2
tic;
GS_iter = 1;
for time = 1:nt
    error = 9e9;
    while(error>to1)
        for i = 2:nx-1
            for j = 2:ny-1
                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)));
        T_old = T;
        GS_iter = GS_iter+1;
    end
    T_prev_dt = T;
end
toc;
time_taken2 = toc;
end

function code for SOR

%SOR function code for transient implicit
function[time_taken3,T,SOR_iter] = SOR_t(T,T_old,T_prev_dt)
global nx ny nt to1 k1 k2 omega
tic;
SOR_iter = 1;
for time = 1:nt
    error = 9e9;
   while(error>to1)
        for i = 2:nx-1
            for j = 2:ny-1
                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);
                H = omega*T(i-1,j)+(1-omega)*T_old(i-1,j)+T_old(i+1,j);
                V = omega*T(i,j-1)+(1-omega)*T_old(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)));
        T_old = T;
        SOR_iter = SOR_iter+1;
    end
    T_prev_dt = T;
end
toc;
time_taken3 = toc;
end

function code for transient explicit method

%explicit function code
function[time_taken4,T,counter] = Explicit(T,T_old)
global nx ny nt k1 k2
tic;
counter = 0;
for time = 1:nt
        for i = 2:nx-1
            for j = 2:ny-1
                H = k1*(T_old(i+1,j)+T_old(i-1,j)-2*T_old(i,j));
                V = k2*(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)));
        T_old = T;
        counter = counter+1;
toc;
time_taken4 = toc;
end

function code for STEADY STATEBIMPLICIT METHOD

function code for jacobi

%jaccob SS function code for Steady implicit
function[time_taken5,T,Jacobi_iter] = JacobiSS(T,T_old)
global nx ny to1
tic;
Jacobi_iter = 1;
    error = 9e9;
    while(error>to1)
        for i = 2:nx-1
            for j = 2:ny-1
                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)));
        T_old = T;
        Jacobi_iter = Jacobi_iter+1;
   end
toc;
time_taken5 = toc;
end

function code for gauss seidal

%gauss Seidel steady state function code
function[time_taken6,T,Seidel_iter] =SeidelSS(T,T_old)
global nx ny to1
tic;
Seidel_iter = 1;
    error = 9e9;
    while(error>to1)
        for i = 2:nx-1
            for j = 2:ny-1
                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)));
        T_old = T;
        Seidel_iter = Seidel_iter+1;
   end
toc;
time_taken6 = toc;
end

function code for SOR

%SORSS function code for Steady implicit
function[time_taken7,T,sor_iter] = SORSS(T,T_old)
global nx ny to1 omega
tic;
sor_iter = 1;
    error = 9e9;
    while(error>to1)
        for i = 2:nx-1
            for j = 2:ny-1
                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)));
        T_old = T;
       sor_iter = sor_iter+1;
   end
toc;
time_taken7 = toc;
end

RESULTS :

graphs:

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.