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

Solving the steady and unsteady 2D Heat Conduction problem

SSiddharth Redkar22 Mar 2020

\"\"

\"\"

Here ,k1 and k2 are the CFL numbers.

Code 1:

%Matlab program to solve the 2D Heat Conduction problem

clear all;
close all;
clc;

%Inputs
L=1;
nx=20;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
T=298*ones(nx,ny);

%Boundary Conditions
T_left=400;
T_top=600;
T_right=800;
T_bottom=900;

%Initializing 
T(:,1)=T_left;
T(end,:)=T_top;
T(:,end)=T_right;
T(1,:)=T_bottom;
%Calculating corner grid point values
T(1,1)=(T_left+T_top)/2;
T(1,ny)=(T_right+T_top)/2;
T(nx,1)=(T_left+T_bottom)/2;
T(nx,ny)=(T_right+T_bottom)/2;

tolerance=1e-4;
alpha=0.05;
dt=0.0125;
t=[0:dt:3];
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
%Functions to solve 2D Heat Conduction equation using various methods.
Jacobi_steady(nx,ny,tolerance,x,dx,y,dy,T);
GS_steady(nx,ny,tolerance,x,dx,y,dy,T);
SOR_steady(nx,ny,tolerance,x,dx,y,dy,T);
explicit_unsteady(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t);
Jacobi_unsteady_implicit(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t);
GS_unsteady_implicit(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t);
SOR_unsteady_implicit(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t);

Code 2:

function Jacobi_steady(nx,ny,tolerance,x,dx,y,dy,T)
    [X Y]=meshgrid(x,y);
    jacobi_iter=0;
    Told=T;
    error=9e9;
    
    tic;                                                                                                                                                                                                                            
    %Convergence loop
    while(error>tolerance)
        %Nodal loop
        for i=2:nx-1
            for j=2:ny-1
                T(i,j)=0.25*(Told(i-1,j)+Told(i+1,j)+Told(i,j-1)+Told(i,j+1));
            end
        end
        error=max(max(abs(Told-T)));
        Told=T;
        jacobi_iter=jacobi_iter+1;
    end
    
    simulation_time=toc;
    figure(1);
    [P,Q]=contourf(X,Y,T);
    clabel(P,Q);
    colorbar;
    colormap(jet);
    xlabel(\'X\');
    ylabel(\'Y\');
    title(sprintf(\'Jacobi Iteration,No of iterations: %d\',jacobi_iter));
    display(\'2D steady state solution using Jacobi iterator\');
    fprintf(\'No of iterations is %d \\n\',jacobi_iter);
    fprintf(\'Time taken for solving problem is %f \\n\',simulation_time);
    fprintf(\'The error obtained is %f \\n\',error);
    display(\'*********************************\');
    
end

Code 3:

function GS_steady(nx,ny,tolerance,x,dx,y,dy,T)
    [X Y]=meshgrid(x,y);
    GS_iter=0;
    Told=T;
    error=9e9;
    
    tic;
    %Convergence loop
    while(error>tolerance)
        %Nodal loop
        for i=2:nx-1
            for j=2:ny-1
                T(i,j)=0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1));
            end
        end
        error=max(max(abs(Told-T)));
        Told=T;
        GS_iter=GS_iter+1;
    end
    
    simulation_time=toc;
    figure(8);
    [P,Q]=contourf(X,Y,T);
    clabel(P,Q);
    colorbar;
    colormap(jet);
    xlabel(\'X\');
    ylabel(\'Y\');
    title(sprintf(\'Gauss Seidel Iteration,No of iterations: %d\',GS_iter));
    display(\'2D steady state solution using Gauss Seidel iterator\');
    fprintf(\'No of iterations is %d \\n\',GS_iter);
    fprintf(\'Time taken for solving problem is %f \\n\',simulation_time);
    fprintf(\'The error obtained is %f \\n\',error);
    display(\'*********************************\');
end

Code 4:

function SOR_steady(nx,ny,tolerance,x,dx,y,dy,T)
    [X Y]=meshgrid(x,y);
    SOR_iter=0;
    Told=T;
    error=9e9;
    omega=1.1;%relaxation factor(over-relaxed)
    
    tic;
    %Convergence loop
    while(error>tolerance)
        %Nodal loop
        for i=2:nx-1
            for j=2:ny-1
                T(i,j)=omega*(0.25*(T(i-1,j)+Told(i+1,j)+T(i,j-1)+Told(i,j+1)))+(1-omega)*Told(i,j);
            end
        end
        error=max(max(abs(Told-T)));
        Told=T;
        SOR_iter=SOR_iter+1;
    end
    
    simulation_time=toc;
    figure(9);
    [P,Q]=contourf(X,Y,T);
    clabel(P,Q);
    colorbar;
    colormap(jet);
    xlabel(\'X\');
    ylabel(\'Y\');
    title(sprintf(\'SOR Iteration,No of iterations: %d\',SOR_iter));
    display(\'2D steady state solution using SOR iterator\');
    fprintf(\'No of iterations is %d \\n\',SOR_iter);
    fprintf(\'Time taken for solving problem is %f \\n\',simulation_time);
    fprintf(\'The error obtained is %f \\n\',error);
    display(\'*********************************\');
end

Code 5:

function explicit_unsteady(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t)
    [X Y]=meshgrid(x,y);
    Told=T;
    
    
    
    tic;
    for k=1:length(t)
    %Nodal loop
    for i=2:nx-1
        for j=2:ny-1
            H=(Told(i-1,j)+Told(i+1,j)-(2*Told(i,j)));
            V=(Told(i,j-1)+Told(i,j+1)-(2*Told(i,j)));
            T(i,j)=(Told(i,j)+k1*(H)+k2*(V));
            
        end
    end
    Told=T;
    end
    simulation_time=toc;
    figure(2);
    [P,Q]=contourf(X,Y,T);
    clabel(P,Q);
    colorbar;
    colormap(jet);
    xlabel(\'X\');
    ylabel(\'Y\');
    hold on;
    title(sprintf(\'Unsteady Explicit\'));
    display(\'2D unsteady state solution,explicit\');
    fprintf(\'Time taken for solving problem is %f \\n\',simulation_time);
    display(\'*********************************\');
end

Code 6:

function Jacobi_unsteady_implicit(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t)
    [X Y]=meshgrid(x,y);
    jacobi_iter=0;
    Told=T;
    T_prev_dt=T;
    error=9e9;
    term1=(1+(2*k1)+(2*k2))^-1;
    term2=k1*term1;
    term3=k2*term1;
    
    tic;
    %Time loop
    for k=1:length(t)
    error=9e9;
    %Convergence loop
    while(error>tolerance)
        %Nodal loop
        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)=(T_prev_dt(i,j)*term1)+(term2*H)+(term3*V);
            end
        end
        error=max(max(abs(Told-T)));
        Told=T;
        jacobi_iter=jacobi_iter+1;
    end
    T_prev_dt=T;
    end
    simulation_time=toc;
    figure(5);
    [P,Q]=contourf(X,Y,T);
    clabel(P,Q);
    colorbar;
    colormap(jet);
    xlabel(\'X\');
    ylabel(\'Y\');
    hold on;
    title(sprintf(\'Jacobi Iteration(Unsteady Implicit),No of iterations: %d\',jacobi_iter));
    display(\'2D unsteady state solution using Jacobi iterator,implicit\');
    fprintf(\'No of iterations is %d \\n\',jacobi_iter);
    fprintf(\'Time taken for solving problem is %f \\n\',simulation_time);
    fprintf(\'The error obtained is %f \\n\',error);
    display(\'*********************************\');
end

Code 7:

function GS_unsteady_implicit(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t)
    [X Y]=meshgrid(x,y);
    GS_iter=0;
    T_prev_dt=T;
    Told=T;
    error=9e9;
    term1=(1+(2*k1)+(2*k2))^-1;
    term2=k1*term1;
    term3=k2*term1;
    
    tic;
    %Time loop
    for k=1:length(t)
    error=9e9;
    %Convergence loop
    while(error>tolerance)
        %Nodal loop
        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));
                T(i,j)=(T_prev_dt(i,j)*term1)+(term2*H)+(term3*V);
            end
        end
        error=max(max(abs(Told-T)));
        Told=T;
        GS_iter=GS_iter+1;
    end
    T_prev_dt=T;
    end
    simulation_time=toc;
    figure(6);
    [P,Q]=contourf(X,Y,T);
    clabel(P,Q);
    colorbar;
    colormap(jet);
    xlabel(\'X\');
    ylabel(\'Y\');
    title(sprintf(\'Gauss Scidel Iteration(Unsteady Implicit),No of iterations: %d\',GS_iter));
    display(\'2D unsteady state solution using Gauss Seidel iterator,implicit\');
    fprintf(\'No of iterations is %d \\n\',GS_iter);
    fprintf(\'Time taken for solving problem is %f \\n\',simulation_time);
    fprintf(\'The error obtained is %f \\n\',error);
    display(\'*********************************\');
end

Code 8:

function SOR_unsteady_implicit(nx,ny,tolerance,x,dx,y,dy,T,k1,k2,t)
    [X Y]=meshgrid(x,y);
    SOR_iter=0;
    Told=T;
    T_prev_dt=T;
    term1=(1+(2*k1)+(2*k2))^-1;
    term2=k1*term1;
    term3=k2*term1;
    error=9e9;
    tic;
    omega=1.1;%relaxation factor(over-relaxed)
    
    %Time loop
    for k=1:length(t)
    error=9e9;
    %Convergence loop
    while(error>tolerance)
        %Nodal loop
        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));
                T(i,j)=omega*((T_prev_dt(i,j)*term1)+(term2*H)+(term3*V))+(1-omega)*(Told(i,j));
            end
        end
        error=max(max(abs(Told-T)));
        Told=T;
        SOR_iter=SOR_iter+1;
    end
    T_prev_dt=T;
    end
    simulation_time=toc;
    figure(7);
    [P,Q]=contourf(X,Y,T);
    clabel(P,Q);
    colorbar;
    colormap(jet);
    xlabel(\'X\');
    ylabel(\'Y\');
    title(sprintf(\'Successive Over Relaxation Iteration(Unsteady Implicit),No of iterations: %d\',SOR_iter));
    display(\'2D unsteady state solution using SOR iterator,implicit\');
    fprintf(\'No of iterations is %d \\n\',SOR_iter);
    fprintf(\'Time taken for solving problem is %f \\n\',simulation_time);
    fprintf(\'The error obtained is %f \\n\',error);
    display(\'*********************************\');
end

Output:

2D steady state solution using Jacobi iterator
No of iterations is 826 
Time taken for solving problem is 0.014018 
The error obtained is 0.000099 
*********************************
2D steady state solution using Gauss Seidel iterator
No of iterations is 440 
Time taken for solving problem is 0.006982 
The error obtained is 0.000099 
*********************************
2D steady state solution using SOR iterator
No of iterations is 365 
Time taken for solving problem is 0.008552 
The error obtained is 0.000099 
*********************************
2D unsteady state solution,explicit
Time taken for solving problem is 0.006386 
*********************************
2D unsteady state solution using Jacobi iterator,implicit
No of iterations is 3265 
Time taken for solving problem is 0.056902 
The error obtained is 0.000049 
*********************************
2D unsteady state solution using Gauss Seidel iterator,implicit
No of iterations is 2283 
Time taken for solving problem is 0.046864 
The error obtained is 0.000070 
*********************************
2D unsteady state solution using SOR iterator,implicit
No of iterations is 1861 
Time taken for solving problem is 0.032384 
The error obtained is 0.000031 
*********************************

Plots:

\"\"

\"\"

\"\"

\"\"

\"\"

\"\"

\"\"

Observation:

1.From the above plots,it can be confirmed that the SOR method takes minimum time/iteration for converging solution than Gauss-Seidel method which takes less time/iteration than Jacobi method while solving the heat conduction equation implicitly using iterative solvers in transient condition.(SOR method used over relaxation factor)

2.Transient state methods takes more computation time and iterations than steady state to converge.

3.Implicit method takes more time computation time than explicit method for a transient solver.

4.Similarly,we can observe from the above plots,that the SOR method takes minimum time/iterations for converging solution than Gauss-Seidel method which takes less time/iteration than Jacobi method while solving the heat conduction equation implicitly using iterative solvers in steady state condition.(SOR method used over relaxation factor)

5. When the CFL number is less than 0.25 we can use explicit method but when CFL number>0.25 its advised to use implicit method because explicit method becomes unstable.

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.