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

SSaif Ansari17 Jul 2023

Aim : To perform Steady state analysis & Transient State Analysis

1. For Steady State Case

clear all
close all
clc
%% Inputs
L = 1;        %% Domain Length
nx = 50;      %% Grid Points
omega = 1.1;  %% Relaxation Factor
error = 9e9;  %% Error
tol = 1e-4;   %% Tolerance Value
ct = 0;
ny=nx;
%% Calculations 
x=linspace(0,L,nx);
dx=x(2)-x(1);
y=linspace(0,L,ny);
dy=y(2)-y(1);
%% Boundary Condition
T=300*ones(nx,ny);
T(:,1)= 400 ;    %% Left Wall
T(1,:)= 600;     %% Top Wall
T(:,end)= 800;   %% Right Wall
T(end,:)= 900;   %% Bottom Wall

T(1,1)= (400+600)/2 ; 
T(1,end)= (600+800)/2 ; 
T(end,1)= (400+900)/2 ; 
T(end,end)= (900+800)/2 ; 

Told = T;
k = 2*((dx^2+dy^2)/(dx^2*dy^2));

%% Command Prompt For Choosing Iterative Solver
prompt = "Select Iterative Solver, Jacobi = 1, Gauss-Siedel = 2, SOR = 3  ";
iterative_solver = input(prompt);

%% Iterative Solution For Jacobi Method 
if iterative_solver == 1
    tic
    Jacobi = 1;
    % Convergence Loop
    while(error>tol)
        % Space Loop
        for i = 2:nx-1
            for j = 2:ny-1
                H = ((Told(i-1,j)+ Told(i+1,j))/(k*dx^2));  
                V = ((Told(i,j-1)+ Told(i,j+1))/(k*dy^2));              
                T(i,j) = H+V;
            end
        end
        error = max(max(abs(Told - T))); 
        Told = T;
        ct = ct + 1;
        
    end
    time = toc;
    %% Plotting
    figure(1)
    contourf(x,y,T,'ShowText','on');
    colorbar
    colormap(jet)
    xlabel('x');
    ylabel('y');
    text = sprintf("Jacobi Method & Iterations For Convergence = %d & Time: %f s",ct ,time);
    title(text);
    set(gca,'YDIR','reverse');
end
%% Iterative Solution For Gauss Siedel Method
if iterative_solver == 2
    tic
    Gauss_Seidel = 2;
    % Convergence Loop
    while(error>tol)
        % Space Loop
        for i = 2:nx-1
            for j = 2:ny-1
                H = ((T(i-1,j)+ Told(i+1,j))/(k*dx^2));  
                V = ((T(i,j-1)+ Told(i,j+1))/(k*dy^2));              
                T(i,j) = H+V;
            end
        end
        error = max(max(abs(Told - T))); 
        Told = T;
        ct = ct + 1;   
    end
    time = toc;
    %% Plotting
    figure(2)
    contourf(x,y,T,'ShowText','on');
    colorbar
    colormap(jet)
    xlabel('x');
    ylabel('y');
    text = sprintf("Gauss Siedel Method & Iterations For Convergence = %d & Time: %f s",ct ,time);
    title(text);
    set(gca,'YDIR','reverse');
end
%% Iterative Solution For Successive Over-Relaxation Method
if iterative_solver == 3
    tic
    SOR = 3;
    % Convergence Loop
    while(error>tol)
        % Space Loop
        for i = 2:nx-1
            for j = 2:ny-1
                H = ((T(i-1,j)+ Told(i+1,j))/(k*dx^2));  
                V = ((T(i,j-1)+ Told(i,j+1))/(k*dy^2));                
                T(i,j) = (1-omega)*(Told(i,j)) + (omega*(H + V));
            end
        end
        error = max(max(abs(Told - T))); 
        Told = T;
        ct = ct + 1;
        
    end
    time = toc;
    %% Plotting
    figure(3)
    contourf(x,y,T,'ShowText','on');
    colorbar
    colormap(jet)
    xlabel('x');
    ylabel('y');
    text = sprintf("Successive Over-Relaxation Method & Iterations For Convergence = %d & Time: %f s",ct ,time);
    title(text);
    set(gca,'YDIR','reverse');
end

Outputs

2. For Transient State Explicit Case

clear all
close all
clc
%% Inputs
L = 1;  %% Domain Length
nx = 50; %% Grid Points
ny=nx;
%% Calculations 
x=linspace(0,L,nx);
dx=x(2)-x(1);
y=linspace(0,L,ny);
dy=y(2)-y(1);
%% Time Marching Step
t = 0.1;
nt = 1400;
dt = t/nt;
%% Coefficient Of Thermal Diffusion
alpha = 1.4; 
%% CFL Condition
CFL_x= alpha*dt/(dx)^2; 
CFL_y= alpha*dt/(dy)^2;
ct = 0;
%% Boundary Condition
T=300*ones(nx,ny);
T(:,1)= 400 ;   %% Left Wall
T(1,:)= 600;    %% Top Wall
T(:,end)= 800;  %% Right Wall
T(end,:)= 900;  %% Bottom Wall

T(1,1)= (400+600)/2 ; 
T(1,end)= (600+800)/2 ; 
T(end,1)= (400+900)/2 ; 
T(end,end)= (900+800)/2 ; 

Told = T;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
%% Loop
tic
for k = 1:nt
    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) = Told(i,j) + k1*H + k2*V;
        end
    end
    Told = T;
    ct = ct+1;
    time = toc;
end
%% Plotting
figure(1)
contourf(x,y,T,'ShowText','on');
colorbar
colormap(jet)
xlabel('x');
ylabel('y');
text = sprintf("Transient Explicit Heat Conduction Method & Iterations For Convergence = %d & Time: %f s",ct ,time);
title(text);
set(gca,'YDIR','reverse');

Output

3. For Transient State Implicit Case

clear all
close all
clc
%% Inputs
L = 1; %% Domain Length
nx = 50; %% Grid Points
omega  = 1.1; %% Relaxation Factor
error = 9e9;  %% Error
tol = 1e-4;   %% Tolerance
alpha = 1.4;  %% Coefficient Of Thermal Diffusion
ct = 0;
ny=nx;
%% Time Marching Steps
t = 0.1; 
nt = 1400; 
dt = t/nt;
%% Calculations
x=linspace(0,L,nx);
dx=x(2)-x(1);
y=linspace(0,L,ny);
dy=y(2)-y(1);
%% Boundary Condition
T=300*ones(nx,ny);
T(:,1)= 400 ;   %% Left Wall
T(1,:)= 600;    %% Top Wall
T(:,end)= 800;  %% Right Wall
T(end,:)= 900;  %% Bottom Wall

T(1,1)= (400+600)/2 ; 
T(1,end)= (600+800)/2 ; 
T(end,1)= (400+900)/2 ; 
T(end,end)= (900+800)/2 ; 

Told = T;
T_prev = T;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);

%% Command Prompt For Choosing Iterative Solver
prompt = "Select Iterative Solver, Jacobi = 1, Gauss-Siedel = 2, SOR = 3  ";
iterative_solver = input(prompt);
%% Iterative Solution For Jacobi Method  
if iterative_solver == 1
    tic
    Jacobi = 1;
    for k = 1:nt
        error = 9e9;
        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) = (term1*(T_prev(i,j)))+(term2*H)+(term3*V);
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            ct = ct+1;
        end
        T_prev = T;
    end
    time = toc;
    %% Plotting
    figure(1)
    contourf(x,y,T,'ShowText','on')
    colorbar
    colormap(jet)
    xlabel('X-axis');
    ylabel('Y-axis');
    text = sprintf("Jacobi Method & Iterations For Convergence = %d & Time: %f s",ct ,time);
    title(text);
    set(gca,'YDIR','reverse');
end
%% Iterative Solution For Gauss Siedel Method
if iterative_solver == 2
    tic
    GS = 2;
    for k = 1:nt
        error = 9e9;
        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 = (Told(i,j+1) + T(i,j-1));            
                    T(i,j) = (term1*(T_prev(i,j)))+(term2*H)+(term3*V);
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            ct = ct+1;
        end
        T_prev = T;
    end
    time = toc;
    %% Plotting
    figure(2)
    contourf(x,y,T,'ShowText','on')
    colorbar
    colormap(jet)
    xlabel('X-axis');
    ylabel('Y-axis');
    text = sprintf("Gauss Siedel Method & Iterations For Convergence = %d & Time: %f s",ct ,time);
    title(text);
    set(gca,'YDIR','reverse');
end
%% Iterative Solution For Successive Over-Relaxation Method
if iterative_solver == 3
    tic
    SOR = 3;
    for k = 1:nt
        error = 9e9;
        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 = (Told(i,j+1) + T(i,j-1));            
                    T(i,j) = (1-omega)*(Told(i,j)) + (omega*((term1*(T_prev(i,j)))+(term2*H)+(term3*V)));
                end
            end
            error = max(max(abs(Told-T)));
            Told = T;
            ct = ct+1;
        end
        T_prev = T;
    end
    time = toc;
    %% Plotting
    figure(3)
    contourf(x,y,T,'ShowText','on')
    colorbar
    colormap(jet)
    xlabel('X-axis');
    ylabel('Y-axis');
    text = sprintf("Successive Over-Relaxation Method & Iterations For Convergence = %d & Time: %f s",ct ,time);
    title(text);
    set(gca,'YDIR','reverse');
end

Output

Compare the convergence rate [Number of iterations] of the steady state and transient [explicit and implicit] simulations and justify the numerics.

Method Iterations For Convergence Time Taken To Reach Convergence(seconds)
Jacobi Method For Steady State Case 4587 0.140518
Gauss Siedel Method For Steady State Case 2469 0.087779
SOR Method For Steady State Case
 2061 0.089496
Transient Explicit Method 1400 0.082067
Jacobi Method For Transient State Case(Implicit) 16675 0.472450
Gauss Siedel Method For Transient State Case(Implicit) 11709 0.349222
SOR Method For Transient State Case(Implicit) 9621 0.326980
  • From the above analysis, we can conclude that for steady case and transient implicit case, the order of convergence is Jacobi Method > Gauss Siedel Method > SOR Method i.e Jacobi takes the most no of iterations and SOR Method takes up the least no of iterations for reaching convergence.
  • Coefficient Of Thermal Diffusion can help us in attaining convergence faster but very high Coefficient Of Thermal Diffusion values can also cause Divergence
  • SOR Method tries to extrapolate value which helps us in reaching convergence with the help of relaxation factor,but very high value of relaxation factor can also lead to divergence
  • Transient Explicit Method is converging way faster than SOR Method also, as in explicit case CFL Factor and Coefficient Of Thermal Diffusion comes into play which is helping us in reaching convergence faster.
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.