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

STEADY AND UNSTEADY 2D HEAT CONDUCTION PROBLEM

SShreyas G11 Aug 2022

STEADY AND UNSTEADY 2D HEAT CONDUCTION PROBLEM:

Given data:

  • The domain under consideration is a unit square
  • Assume nx = ny
  • Left = 400K
  • Right = 800K
  • Top = 600K
  • Bottom = 900K

Assumptions used in explanation:

  • `T_(i+1,j) = T_R`
  • `T_(i-1,j) = T_L`
  • `T_(i,j+1) = T_T`
  • `T_(i,j-1) = T_B`
  • `T_(i,j) = T_P` (Current point)

STEADY-STATE ANALYSIS:

The steady-state equation is given by:

`(del^2T)/(del x^2) + (del^2 T)/(del x^2) = 0`

Discretzing the above equation using CDS we get:

`(T_L - 2T_P + T_R) / (Delta x^2) + (T_T - 2T_P + T_B) / (Delta y^2) = 0`

Rearranging the terms,

`(2T_P) / (Delta x^2) + (2T_P)/(Delta y^2) = [(T_L + T_R) / (Delta x^2)] + [(T_T + T_B) / (Delta y^2)]`

But we have nx = ny we can consider `Delta x = Deltay`, hence the equation becomes:

`T_P = 0.25 (T_L + T_R + T_B + T_T)`

Code for Steady state 2D heat conduction:

clear all
close all
clc

%No of grid points
nx = 50;
ny = nx;

%Initializing temperate at all node points
u = 400*ones(nx,ny);

%Boundary Conditions
u(:,ny)=800;
u(1,:)=600;
u(:,1)=400;
u(nx,:)=900;

%Making a copy of u
uold = u;

%Grid generation
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = x;
dy = y(2)-y(1);

iterative_solver = 1;
tol = 1e-5;
error = 9e9;
[X,Y] = meshgrid(x,y);
tic;

%Jacobi Method
if(iterative_solver == 1)
    jacobi_iter = 1;
    while(error>tol)
        for i = 2:nx-1
            for j = 2:ny-1
            u(i,j) = 0.25*(uold(i,j-1)+uold(i,j+1)+uold(i-1,j)+uold(i+1,j));
            end
        end
        error = max(abs(uold(:)-u(:)));
        uold = u;
        jacobi_iter = jacobi_iter+1;
    end
end

%Gauss Seidel Method
if(iterative_solver == 2)
    gs_iter = 1;
    while(error>tol)
        for i = 2:nx-1
            for j = 2:ny-1
            u(i,j) = 0.25*(u(i,j-1)+uold(i,j+1)+u(i-1,j)+uold(i+1,j));
            end
        end
        error = max(abs(uold(:)-u(:)));
        uold = u;
        gs_iter = gs_iter+1;
    end
end

%SOR(Succive over relaxation) Method
f = 2/(1+sin(pi*dx));
if(iterative_solver == 3)
    SOR_iter = 1;
    while(error>tol)
        for i = 2:nx-1
            for j = 2:ny-1
            u(i,j) = ((1-f)*uold(i,j)+(f*0.25*(u(i,j-1)+uold(i,j+1)+u(i-1,j)+uold(i+1,j))));
            end
        end
        error = max(abs(uold(:)-u(:)));
        uold = u;
        SOR_iter = SOR_iter+1;
    end
end

time = toc;
figure(1)
[C,h] = contourf(X,Y,u);
colorbar;
colormap(jet);
clabel(C,h);

if(iterative_solver == 1)
    title(sprintf('Temperature Distribution of Steady state heat conduction using Jacobi Method \n No. of Iterations: %d \n Computation Time: %0.3f s \n \n',jacobi_iter,time));
end

if(iterative_solver == 2)
    title(sprintf('Temperature Distribution of Steady state heat conduction using Gauss Seidel Method \n No. of Iterations: %d \n Computation Time: %0.3f s \n \n',gs_iter,time));
end

if(iterative_solver == 3)
    title(sprintf('Temperature Distribution of Steady state heat conduction using SOR Method \n No. of Iterations: %d \n Computation Time: %0.3f s \n Relaxation factor: %0.3f \n \n',SOR_iter,time,real(f)));
end

OUTPUT for Steady state 2D heat conduction:

Fig 1. Temperature Distribution of Steady state heat conduction using Jacobi Method

Fig 2. Temperature Distribution of Steady state heat conduction using Gauss Seidel Method

Fig 3. Temperature Distribution of Steady state heat conduction using SOR method.

From the plots obtained we can observe that the number of iterations and the computational time required is varied in the following pattern:

SOR>Gauss Seidel>Jacobi Method.

TRANSIENT STATE ANALYSIS:

2D Transient Heat Equation is given by:

`(del T) / (del t) + alpha [(del^2 T) / (del x^2) + (del^2 T) / (del y^2)] = 0`

Appling CDS to the above equation we get:

`[T_P^(n+1) - T_P^n] + alpha[(T_R -2T_P + T_L)^n / (Delta x^2) + (T_T - 2T_P + T_B)^n / (Delta y^2) ]`

But `Delta x = Delta y`

Solving and Rearranging the terms we get:

`T_p^(n+1) = k (T_L^n + T_R^n + T_T^n + T_B^n) + [1 - 4k] T_P^n`

where `k = alpha (Delta T) / (Delta x^2) = alpha (Delta T) / (Delta y^2)`

Code for 2D Transient heat conduction - Explicit Method:

clear all 
close all
clc

%No of grid points
nx = 50;
ny = nx;

%No of time steps
nt = 2000;

%Assuming value of c
c = 1.656e-4;

%Initializing temperate at all node points
u = 400*ones(nx,ny);

%Boundary Conditions
u(:,ny)=800;
u(1,:)=600;
u(:,1)=400;
u(nx,:)=900;

%Making a copy of u
uold = u;

%Grid generation
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = x;
dy = y(2)-y(1);
dt=0.1;

% k = (c*dt)/(dx^2);
k=0.26;
iterative_solver = 3;
tol = 1e-5;
error = 9e9;
[X,Y] = meshgrid(x,y);
tic;
for p = 1:nt
    for i = 2:nx-1
        for j = 2:ny-1
            u(i,j) = k*(uold(i-1,j) + uold(i,j-1) + uold(i+1,j) + uold(i,j+1)) + (1-4*k)*uold(i,j);
        end
    end
    error=max(abs(uold(:)-u(:)));
    uold=u;
end
time = toc;
figure(1)
[C,h] = contourf(X,Y,u);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('Unsteady state heat conduction (Explicit Method) \n Time: %0.3f \n Computation Time: %0.3f \n \n',dt*nt,time));

OUTPUT for 2D Transient heat conduction - Explicit Method:

IMPLICIT METHOD:

Implicit Equation is given by:

`[T_P^(n+1) - T_P^n] + alpha[(T_R -2T_P + T_L)^(n+1) / (Delta x^2) + (T_T - 2T_P + T_B)^(n+1) / (Delta y^2) ]`

Simplifying and Rearranging the terms we get:

`T_p^(n) = [1 + 4k] T_P^(n+1) - k (T_L^(n+1) + T_R^(n+1) + T_T^(n+1) + T_B^(n+1))`

where `k = alpha (Delta T) / (Delta x^2) = alpha (Delta T) / (Delta y^2)`

Code for 2D Transient heat conduction - Implicit Method:

clear all
close all
clc

%No of grid points
nx = 50;
ny = nx;
dt=0.1;

%No of time steps
nt = 2000; 

%Grid generation
x = linspace(0,1,nx);
dx = x(2)-x(1);
y = x;
dy = y(2)-y(1);
[X,Y] = meshgrid(x,y);


%Assuming value of c
c = 1.656e-4;

k=(c*dt)/(dx^2);

%Initializing temperate at all node points
u = 400*ones(nx,ny);
c=0;

%Boundary Conditions
u(2:ny,nx)=800;
u(1,2:nx)=600;
u(1:ny-1,1)=400;
u(ny,1:nx-1)=900;

%Making a copy of u
uold = u;

t1=u;
r1=u;
k1=k/2;
k2=k1;

a=zeros(nx*ny, nx*ny);
b=zeros(nx*ny, 1);

for i = 2:nx-1
    for j = 2:ny-1
        a(j+ny*(i-1),j+ny*(i-1))=1+2*k1+2*k2;
        a(j+ny*(i-1),j+ny*(i-1)-1)=-k1;
        a(j+ny*(i-1),j+ny*(i-1)+1)=-k1;
        a(j+ny*(i-1),j+ny*(i-2))=-k2;
        a(j+ny*(i-1),j+ny*(i))=-k2;
        b(j+ny*(i-1))=u(i,j);
    end
end

%L/R Boundary
for j = 2:ny
    a(j*ny,j*ny)=1;
    b(j*ny)=800;
    a(nx*(j-2)+1,nx*(j-2)+1)=1;
    b(nx*(j-2)+1)=400;
end

%T/B Boundary
for r=2:ny
    a(r,r)=1;
    a((ny*ny)-r+1,(ny*ny)-r+1)=1;
    b(r)=600;
    b((ny*ny)-r+1)=900;
end
tic;
for q=1:nt
    x1=inv(a)*b;
    b=x1;
    for i = 2:nx-1
        for j =2:ny-1
            u(i,j)=b(j+ny*(i-1));
        end
    end
end

time=toc;
figure(1)
[C,h] = contourf(X,Y,u);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('Unsteady state heat conduction (Implicit Method) \n Time: %0.3f \n Computation Time: %0.3f \n \n',dt*nt,time));

OUTPUT for 2D Transient heat conduction - Implicit Method:

From the plots obtained we observe that the time required for computing the Implicit Method is more than the Explicit method but the Implicit method generates more stable results than Explicit method.

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.