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

YYash Deshmukh30 Nov 2021
% function code for jacobi method

function T = jacobi_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told) 

jacobi_iter= 1;
while (start_error>error_tol)
   
    for i=2:nx-1
        for j=2:ny-1
            term1 = (Told(i-1,j)+Told(i+1,j))/(k*(dx^2));
            term2 = (Told(i,j-1)+Told(i,j+1))/(k*(dy^2));
            T(i,j) = term1+term2;
        end
    end
    start_error=max(max(abs(Told-T))); 
    Told=T;
    jacobi_iter=jacobi_iter+1;
   figure(1)
    [c,d]=contourf(xx,yy,T);
    colormap(jet)
    colorbar
    title_text=sprintf('jacobi iteration number=%d @steady state',jacobi_iter);
    title(title_text)
    xlabel('x axis')
    ylabel('y axis')
    clabel(c,d)
    pause(0.003)
fprintf('No. of iteration in jacobi =%d',jacobi_iter) 
end
 figure(1)
    [c,d]=contourf(xx,yy,T);
    colormap(jet)
    colorbar
    title_text=sprintf('jacobi iteration number=%d @steady state',jacobi_iter);
    title(title_text)
    xlabel('x axis')
    ylabel('y axis')
    clabel(c,d)
    pause(0.003)
fprintf('No. of iteration in jacobi =%d',jacobi_iter)
end 
% function code for gauss seidal
function gauss_siedal= gauss_iter(start_error,error_tol,x,y,nx,ny,dx,dy,T,k,Told) 
GS_iter = 1;
while(start_error>error_tol)
    for i=2:nx-1
        for j=2:ny-1
            T(i,j)=((((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k);
        end
    end
    start_error=max(max(abs(Told-T)));  
    Told=T;
    GS_iter=GS_iter+1;
    gauss_siedal=figure(2);
    [c,d]=contourf(x,y,T); 
    colormap(jet)
    colorbar
    title_text=sprintf('GS iteration number=%d @steady state',GS_iter);
    title(title_text)
    xlabel('x-axis')
    ylabel('y-axis')
    clabel(c,d)
    pause(0.003)
end
fprintf('No. of iteration in GS=%d',GS_iter)
end
% function code for  Successive over relaxation
function T = sor_iter(start_error,error_tol,x,y,nx,ny,dx,dy,T,k,Told) 
sor_iter=1;
omega=1.5; 
while(start_error>error_tol)
    for i=2:nx-1
        for j=2:ny-1
            T_gs=((((T(i-1,j)+Told(i+1,j))/dx^2)+((T(i,j-1)+Told(i,j+1))/dy^2))/k);
            T(i,j)=(omega*(T_gs)+(1-omega)*(Told(i,j)));
        end
    end
    start_error=max(max(abs(Told-T))); 
    Told=T;
    sor_iter=sor_iter+1;
    figure(3)
    [c,d]=contourf(x,y,T); 
    colormap(jet)
    colorbar
    title_text=sprintf('sor iteration number=%d @steady state',sor_iter);
    title(title_text)
    xlabel('x-axis')
    ylabel('y-axis')
    clabel(c,d)
    pause(0.003)
end
fprintf('No. of iteration in sor=%d',sor_iter)
end
% main code
%solving steady state equation using different solver
clear all
close all
clc

% Input
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1); 
[xx,yy]=meshgrid(x,y);

%Boundary condition
T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
T(1,1) = 650;
T(1,end) = 850;
T(end,1) = 500;
T(end,end) = 700;
Told=T;

start_error=9e9;
error_tol=1e-4;
k=2*((dx^2+dy^2)/(dx^2*dy^2));

iterative_solver=input('enter number with respect to solver:');
if iterative_solver==1
    tic;
    T_jacobi = jacobi_iter(start_error,error_tol,xx,yy,nx,ny,dx,dy,T,k,Told); 
    time_jacobi_in_sec=toc;
elseif iterative_solver==2
    tic;
    T_gauss=gauss_iter(start_error,error_tol,x,y,nx,ny,dx,dy,T,k,Told);
    time_gauss_in_sec=toc;
elseif iterative_solver==3
    tic;
    T_sor=sor_iter(start_error,error_tol,x,y,nx,ny,dx,dy,T,k,Told); 
    time_sor_in_sec=toc;
end

Code for Transient explicit 

%function code for explicit

function T = jacobi_iter_tran_expli(xx,yy,nx,ny,T,k1,k2,Told)

jacobi_iter=1;

for k=1:1200
    for i=2:nx-1
        for j=2:ny-1
            term1 = Told(i,j);
            term2 = k1*(Told(i+1,j)-2*Told(i,j)+Told(i-1,j));
            term3 = k2*(Told(i,j+1)-2*Told(i,j)+Told(i,j-1));
            T(i,j) = term1+term2+term3;
        end
    end
    start_error=max(max(abs(Told-T)));
    Told=T;
    jacobi_iter = jacobi_iter+1;
end
figure(1)
[c,d] = contourf(xx,yy,T);
colormap(jet)
colorbar
title_text=sprintf('jacobi iteration number = %d @transient state explicitly',jacobi_iter);
title(title_text)
xlabel('x axis')
ylabel('y axis')
clabel(c,d)
pause(0.003)
fprintf('No. of iteration in jacobi = %d\n',jacobi_iter);
end
% main code for transient heat conduction equation explicitly 

clear all
close all
clc

% inputs
nx=10;
ny=nx;
x=linspace(0,1,nx);
y=linspace(0,1,nx);
dx=x(2)-x(1);
dy=y(2)-y(1); 
dt=1e-3;
alpha=1.4;
[xx,yy]=meshgrid(x,y);

%Boundary condition
T=ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;
Told=T;

% rounding off corner value 
T(1,1)=(400+900)/2;
T(1,ny)=(900+800)/2;
T(nx,1)=(400+600)/2;
T(nx,ny)=(600+800)/2;

error_tol=1e-4;
start_error=9e9;
k1=alpha*(dt/dx^2);
k2=alpha*(dt/dy^2);

% defining solver
iterative_solver=1;
tic;
if iterative_solver==1
    T_jacobi=jacobi_iter_tran_expli(xx,yy,nx,ny,T,k1,k2,Told);
end
time_jacobi_in_sec=toc;

Code for Transient implicit 

% function code for jacobi (implicit)

function T=jacobi_iter_tran_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)

jacobi_iter=1;
for k = 1:1200
    start_error=9e9;
    while(start_error > error_tol)
        for i=2:nx-1
            for j=2:ny-1
                term1 = (1+(2*k1)+(2*k2))^(-1);
                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)=((Tprev(i,j)*term1)+term2*H+term3*V);
            end
        end
        start_error=max(max(abs(Told-T)));
        Told=T;
        jacobi_iter=jacobi_iter+1;
    end
    Tprev=T;
end
time_jacobi=toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('x axis')
ylabel('y axis')
colorbar
title_text=sprintf('Solution for 2D heat conduction equation solved implicitly using jacobi with no of iteration =%d and time to solve =%fsec',jacobi_iter,time_jacobi);
title(title_text)
fprintf('No of iteration in jacobi=%d and time in sec=%f',jacobi_iter,time_jacobi);
end

% function code for gauss siedal (implicit)

function T=gs_iter_tran_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev)
gs_iter=1;
for k=1:1200
    start_error=9e9;
    while(start_error > error_tol)
    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)+Told(i+1,j));
            V=(T(i,j-1)+Told(i,j+1));
            T(i,j)=(Tprev(i,j)*term1+term2*H+term3*V);
        end
    end
    start_error=max(max(abs(Told-T)));
    Told=T;
    gs_iter=gs_iter+1;
end
Tprev=T;
end
time_gs=toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('x axis')
ylabel('y axis')
colorbar
title_text=sprintf('Solution for 2D heat conduction equation solved implicitly using gauss siedal witn no of iteration =%d and time to solve =%fsec',gs_iter,time_gs);
title(title_text)
fprintf('No of iteration in gauss siedal=%d and time in sec =%f',gs_iter,time_gs);
end

% fucntion code for Successive over relaxation (implicit)

function T=sor_iter_tran_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega)
sor_iter=1;
for k=1:1200
    start_error=9e9;
    while(start_error > error_tol)
    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)+Told(i+1,j));
            V=(T(i,j-1)+Told(i,j+1));
            T(i,j)=(omega*(Tprev(i,j)*term1+term2*H+term3*V)+(1-omega)*Tprev(i,j));
        end
    end
    start_error=max(max(abs(Told-T)));
    Told=T;
    sor_iter=sor_iter+1;
    end
    Tprev=T;
end
time_sor=toc;
contourf(xx,yy,T,'ShowText','on')
colormap(jet)
colorbar
xlabel('x axis')
ylabel('y axis')
colorbar
title_text=sprintf('Solution for 2D heat conduction equation solved implicitly using SOR witn no of iteration =%d and time to solve =%fsec',sor_iter,time_sor);
title(title_text)
fprintf('No of iteration in SOR=%d and time in sec =%f',sor_iter,time_sor);
end
% main code implicit 
clear all
close all
clc

% inputs
nx=10;
ny=nx;
L=1;
x=linspace(0,1,nx);
y=linspace(0,1,nx);

dx=x(2)-x(1);
dy=y(2)-y(1);
dt=1e-3;
alpha=1.4;
[xx,yy] = meshgrid(x,y);

% boundary condition
T = ones(nx,ny);
T(1,:)=900;
T(end,:)=600;
T(:,1)=400;
T(:,end)=800;

% rounding off corner value 
T(1,1)=(400+900)/2;
T(1,ny)=(900+800)/2;
T(nx,1)=(400+600)/2;
T(nx,ny)=(600+800)/2;
Told=T;
Tprev=Told;
error_tol=1e-4;
start_error=9e9;
k1 = alpha*(dt/dx^2);
k2 = alpha*(dt/dy^2);
omega = 1.2;
CFL_1 = alpha*(dt/dx^2);
CFL_2 = alpha*(dt/dy^2);
CFL = CFL_1+CFL_2;

% defining solver

iterative_solver=input('enter number with respect to solver:');

if iterative_solver==1
    tic;
    T=jacobi_iter_tran_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
elseif iterative_solver==2
    tic;
    T_gs=gs_iter_tran_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev);
elseif iterative_solver==3
    tic;
    T_sor=sor_iter_tran_impli(error_tol,xx,yy,nx,ny,T,k1,k2,Told,Tprev,omega);
end

General 3D heat conduction equation has both time derivative and spatial derivative terms. 3D heat conduction equation is given below.

`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2)+(del^2T)/(delz^2))=0`

The equation above has two terms where one term represent time derivative and other term is spatial derivative and alpha is the thermal diffusivity.

2D heat conduction is same as 3D but z is neglected which is given below.

`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`

Steady state equation is under condition which is the rate of change of tempreature remain zero and the time derivative term become zero.

Steady state 2D heat equation is given below:-

`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`

Transient state equation in this conditon both the term are considered and equation is below:-

`(delT)/(delt)+alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`

Equation can be solved two ways by using three iterative method such as jacobi,Gauss-seidal and successive over relaxation:-

1) Explicit Method: -

Explicit method is used in numerical analysis for obtaining numerical approximations to the solution of the time dependent ordinary or partial differential equations as required in computer simulation of physical process. It calculate the state of a system at later time from the state of system at current time. In this method each diffrential equation approach contain only one unknown using this method.

2) Implicit Method: -

Implicit approach is used in numercial analysis for obtaining numerical approximations to the solution of the time dependent ordinary or partial differential equations as required in computer simulation of physical process. Implicit method find the solution by solving an equation involving both the current state of the system and later one. Mathematically if Y(t) is the current system state and Y(T+dT) is the state at later time dt then for an explicit method Y(T+dt)=F(Y(t)).

Conclusion: -

In steady state equation SOR method yeild the result very quickly comapre to other two metghod whearas jacobi take more number of iteration and more number of iteration and more to solve unknown. Therefore SOR is the most efficient when comapred to other. For transient the time required to solve is more when to steadyn state equation. So it is recomended to solve the 2d heat conduction equation with steady state equation. When transient is solved implicitly jacobi take less time when compare to the other two method whearas SOR is more stable compared to gauss seidal.

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.