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

SSheetanshu Deepak25 Dec 2021
MATLAB
% Steady State 2D Heat conduction equation
clc 
clear all
close all

% Input
nx = 10;   % number of grids points in x direction
ny = nx;   % number of grids points in y direction
Told = magic(10);

% Temperature in Kelvin
T_l = 400;
T_r = 800;
T_b = 900;
T_t = 600;
T(1,1) = (T_t+T_l)/2;
Told(nx,1) = (T_l+T_b)/2;
Told(1,ny) = (T_r+T_t)/2;
Told(nx,ny) = (T_r+T_b)/2;
Told(1,2:ny-1) = 600;
Told(nx,2:ny-1) = 900;
Told(2:nx-1,ny) = 800;
Told(2:nx-1,1) = 400;
Told(2:nx-1,2:ny-1) = 300;

Tol = 0.0001;
% Length of domain
L = 1;
x = linspace(0,L,nx);
y = x;
dx = x(2)-x(1);
dy = dx;
T = Told; % HERE WE ARE ASSUMING SOME OF VALUES AT SOME NODES BECAUSE WHILE CALCULATING TEMPERATUE THERE WILL SOME % NODES WHICH ARE YET TO COMPUTE SO MATLAB AUTOMATICALLY VALUE FROM HERE
w = 1.1
k = 2*((((dx)^2) + (dy)^2)/(((dx)^2)*((dy)^2)));
iterative_solver = 1 %Jacobi
%iterative_solver = 2 % Gauss Seidal
%iterative_solver = 3  % SOR

%Jacobi Method
% Temp profile by the iterative solver
error = 10;
 itt=1;
 if iterative_solver==1;
     while(error>Tol)
         iter = 1;
         for i = 2:nx-1;
             it=19;
             for j = 2:ny-1;
                 A = ((Told(i-1,j) + Told(i+1,j))/((dx)^2));
                 B = ((Told(i,j+1) + Told(i,j-1))/((dy)^2));
                 T(i,j)= (1/k)*(A+B)
                 it=it+1;
             end
             iter=iter+1;
         end
         T(1,2:ny-1)=600;
         T(1,1)=500;
         T(nx,2:ny-1)=900;
         T(2:ny-1,1)=400;
         T(nx,1)=650;
         T(2:ny-1,10)=800;  
         T(1,ny)=700;
         T(nx,ny)=850;
         error = max(max(abs(Told-T)));
         Told=T;
         itt= itt+1
     end
     contourf(x,y,Told) % it plot the one dependent variable
     colorbar % it dark the color
     title_text = sprintf('Jacobi Method')
    title(title_text)
     set(gca,'YDIR','reverse') % it reverse the axis
     xlabel('length of domain(x)')
     ylabel('length of domain(y)')
 end
% Temp profile by the iterative solver

% Gauss-Siedal
if iterative_solver ==2
    while(error>Tol)
        iter = 1;
        for i = 2:nx-1;
            it=1;
            for j = 2:ny-1;
                A = ((T(i-1,j) + Told(i+1,j))/((dx)^2));
                B = ((Told(i,j+1) + T(i,j-1))/((dy)^2));
                T(i,j) = (1/k)*(A+B);
                it=it+1;
            end
            iter=iter+1;
        end
        T(1,2:ny-1)=600;
         T(1,1)=500;
         T(nx,2:ny-1)=900;
         T(2:ny-1,1)=400;
         T(nx,1)=650;
         T(2:ny-1,10)=800;  
         T(1,ny)=700;
         T(nx,ny)=850;
        error = max(max(abs(Told-T)));
        Told=T;
        itt= itt+1
    end
    contourf(x,y,Told)
    colorbar
    title_text = sprintf('Gauss Method')
    title(title_text)
    set(gca,'YDIR','reverse')
    xlabel('length of domain(x)')
    ylabel('length of domain(y)')
end
% SOR
if iterative_solver == 3;
    % Temp profile by the iterative solver
    error = 10;
    itt=1;
    while(error>Tol)
        iter = 1;
        for i = 2:nx-1;
            it=19;
            for j = 2:ny-1;
                A = ((T(i-1,j) + Told(i+1,j))/((dx)^2));
                B = ((Told(i,j+1) + T(i,j-1))/((dy)^2));
                T(i,j) = (1-w)*Told(i,j)+w*((1/k)*(A+B));
                it=it+1;
            end
            iter=iter+1;
        end
        T(1,2:ny-1)=600;
         T(1,1)=500;
         T(nx,2:ny-1)=900;
         T(2:ny-1,1)=400;
         T(nx,1)=650;
         T(2:ny-1,10)=800;  
         T(1,ny)=700;
         T(nx,ny)=850;
        error = max(max(abs(Told-T)));
        Told=T;
        itt= itt+1
    end
    contourf(x,y,Told)
    colorbar;
    title_text = sprintf('SOR Method')
    title(title_text)
    set(gca,'YDIR','reverse')
    colormap(jet);
    xlabel('length of domain(x)')
    ylabel('length of domain(y)')
end

Result

It is evident form the above plot is that SOR take very less time to get converge

% Unsteady state 2D Heat transfer equation(Explicit approach)
clc 
clear all
close all

% Input
nx = 10;   % number of grids points in x direction
ny = nx;   % number of grids points in y direction
Told = magic(10);

T_l = 400;
T_r = 800;
T_b = 900;
T_t = 600;
Told(1,1) = (T_t+T_l)/2;
Told(nx,1) = (T_l+T_b)/2;
Told(1,ny) = (T_r+T_t)/2;
Told(nx,ny) = (T_r+T_b)/2;
Told(1,2:ny-1) = 600;
Told(nx,2:ny-1) = 900;
Told(2:nx-1,ny) = 800;
Told(2:nx-1,1) = 400;
Told(2:nx-1,2:ny-1) = 300;

Tol = 0.0001;
% Length of domain
L = 1;
x = linspace(0,L,nx);
y = x;
dx = x(2)-x(1);
dy = dx;
%T = Told; % THIS IS DEFIND THAT FIRST TIME STEPS HAS VALUE WHICH IS AUTOMATICALLY TAKE FROM HERE.
w = 1.1
%k = 2*((((dx)^2) + (dy)^2)/(((dx)^2)*((dy)^2)));
K1 = 0.2;% K1= (thermal diffusitivity*time steps)/(dx)2)
K2 = K1;
Total_time = 200  % integration would happen uo to 200 sec, dt(time steps), nt(number of time steps through which integartion would happen)
dt = 0.1 % assumed
nt =100/dt;

% Here no Iterative method required becasue all parameters are known RHS
% Because when value at  nodes is find out for next time steps then we have
% all values at all nodes for previous time steps then it means all
% parameter are known to us so no iterative method is required here

p=1
for t = 1:100/dt;
    
    % Space loop
    for i = 2:nx-1;
        for j = 2:ny-1;
            H=(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
            V=(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
            T(i,j)=K1*H + K2*V+Told(i,j) % here Told(i,j) reprsents the value at i,j node at previous time steps
        end 
    end
   T(1,2:ny-1)=600;
   T(1,1)=500;
   T(nx,2:ny-1)=900;
   T(2:ny-1,1)=400;
   T(nx,1)=650;
   T(2:ny-1,10)=800;  
   T(1,ny)=700;
   T(nx,ny)=850;;
    %error = max(max(abs(Told-T)))
    Told=T;
    % FOR EACH TIME STEPS, THE SPACE LOOP WILL RUN AND NEXT TIME STEPS
    % WOULD USE THE VALUE OF PREVIOUS NODE CALCULATED AT PREVIOUS TIME STEP
    % TO COMPUTE THE NEXT NODE AT NEXT TIME STEPS.
    p=p+1
end
contourf(Told)
txt = sprintf('Iteration Number=%d',itt);
colormap(jet);
colorbar
set(gca,'YDIR','reverse')
xlabel('length of domain(x)')
ylabel('length of domain(y)')

Result

Since all are known at RHS so no iterative method required to solove.

% Unsteady state 2D Heat transfer equation(IMplicit approach)
clc 
clear all
close all

% Input
nx = 10;   % number of grids points in x direction
ny = nx;   % number of grids points in y direction
Told = magic(10);

T_l = 400;
T_r = 800;
T_b = 900;
T_t = 600;
Told(1,1) = (T_t+T_l)/2;
Told(nx,1) = (T_l+T_b)/2;
Told(1,ny) = (T_r+T_t)/2;
Told(nx,ny) = (T_r+T_b)/2;
Told(1,2:ny-1) = 600;
Told(nx,2:ny-1) = 900;
Told(2:nx-1,ny) = 800;
Told(2:nx-1,1) = 400;
Told(2:nx-1,2:ny-1) = 300;

Tol = 0.0001;
T = Told;
Tpre = Told;
% Length of domain
L = 1;
x = linspace(0,L,nx);
y = x;
dx = x(2)-x(1);
dy = dx;
%T = Told; % THIS IS DEFIND THAT FIRST TIME STEPS HAS VALUE WHICH IS AUTOMATICALLY TAKE FROM HERE.
w = 1.1
%k = 2*((((dx)^2) + (dy)^2)/(((dx)^2)*((dy)^2)));
K1 = 0.2;% K1= (thermal diffusitivity*time steps)/(dx)2)
K2 = K1;
Total_time = 100  % integration would happen uo to 200 sec, dt(time steps), nt(number of time steps through which integartion would happen)
dt = 0.1 % assumed
nt =100/dt;
% Here no Iterative method required becasue all parameters are known RHS
% Because when value at  nodes is find out for next time steps then we have
% all values at all nodes for previous time steps then it means all
% parameter are known to us so no iterative method is required here
error = 10;
  it = 1
   p=1
   f=1
   q=1
 %iterative_solver = 1 %Jacobi
%iterative_solver = 2 % Gauss Seidal
iterative_solver = 3  % SOR
if   iterative_solver == 1

for t = 1:100/dt
    % Space loop
     error = 0.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));
            term1=(1+(2*K1)+(2*K2))^-1;
            term2=K1*term1;
            term3=K2*term1;
            T(i,j)=(Tpre(i,j)*term1)+(H*term2)+(V*term3);
        end
        
        
    end
    
    T(1,2:ny-1)=600;
    T(1,1)=500;
    T(nx,2:ny-1)=900;
    T(2:ny-1,1)=400;
    T(nx,1)=650;
    T(2:ny-1,10)=800;  
    T(1,ny)=700;
    T(nx,ny)=850;
    error = max(max(abs(Told-T)));
    Told=T
    p=p+1;
    end
    it= it+1;
    Tpre = T;
end
contourf(x,y,Told)
colorbar;
%title_text = sprintf('SOR Method')
%txt = sprintf('Iteration Number(SOR method)=%d',itt);
% title(title_text)
set(gca,'YDIR','reverse')
colormap(jet);
xlabel('length of domain(x)')
ylabel('length of domain(y)')
end
if iterative_solver==2
    for t = 1:100/dt
    % Space loop
     error = 0.1
   
    while(error>Tol)
     
    for i = 2:nx-1
      
        for j = 2:ny-1
            H=(T(i-1,j)+Told(i+1,j));
            V=(T(i,j-1)+Told(i,j+1));
            term1=(1+(2*K1)+(2*K2))^-1;
            term2=K1*term1;
            term3=K2*term1;
            T(i,j)=(Tpre(i,j)*term1)+(H*term2)+(V*term3);
        end
        
        
    end
    
    T(1,2:ny-1)=600;
    T(1,1)=500;
    T(nx,2:ny-1)=900;
    T(2:ny-1,1)=400;
    T(nx,1)=650;
    T(2:ny-1,10)=800;  
    T(1,ny)=700;
    T(nx,ny)=850;
    error = max(max(abs(Told-T)));
    Told=T
    p=p+1;
    end
    it= it+1;
    Tpre = T;
end
contourf(x,y,Told)
colorbar;
%title_text = sprintf('SOR Method')
%txt = sprintf('Iteration Number(SOR method)=%d',itt);
% title(title_text)
set(gca,'YDIR','reverse')
colormap(jet);
xlabel('length of domain(x)')
ylabel('length of domain(y)')
end
if iterative_solver==3
    for t = 1:100/dt
    % Space loop
     error = 0.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));
            term1=(1+(2*K1)+(2*K2))^-1;
            term2=K1*term1;
            term3=K2*term1;
            T(i,j)=w*((Tpre(i,j)*term1)+(H*term2)+(V*term3))+(1-w)*Told(i,j)
        end
        
        
    end
    
    T(1,2:ny-1)=600;
    T(1,1)=500;
    T(nx,2:ny-1)=900;
    T(2:ny-1,1)=400;
    T(nx,1)=650;
    T(2:ny-1,10)=800;  
    T(1,ny)=700;
    T(nx,ny)=850;
    error = max(max(abs(Told-T)));
    Told=T
    p=p+1
    end
    it= it+1;
    Tpre = T;
end
contourf(x,y,Told)
colorbar;
%title_text = sprintf('SOR Method')
%txt = sprintf('Iteration Number(SOR method)=%d',itt);
% title(title_text)
set(gca,'YDIR','reverse')
colormap(jet);
xlabel('length of domain(x)')
ylabel('length of domain(y)')
end

Result- It has seen that SOR takes less time as compare to the other to get converge.

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.