Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

The 2d heat conduction equation is `(delT)/(delt)+alpha*((del^2T)/(delx)+(del^2T)/(dely))=0`

T=temperature along x and y direction, t=time, aplha=thermal diffusivity.

A) Steady State :

let us first consider the 2D state state part of equation : (del^2T)/(delx)+(del^2T)/(dely)=0'

The steady state 2D heat conduction equation is solved using Iterative techiques:

1) Jacobi Method

2)Gauss Siedel Method

3)SOR Method

The Steady state 2D heat conduction equation is discritized using central difference scheme

in X and Y direction.

we get,

T(i,j)=1/k*(T(L)+T(R)/DX^2)+1/K*(T(T)+T(B)/DY^2)

Where K=2*(dx^2+dy^2/dx^2*dy^2)

B) Unsteady State:

The 2D unsteady state heat conduction equation :

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

The above equation is discritized using forward diffrence scheme on the left hand side for space derivative and forward difference for time derivative. Depending on selection of nodes in time direction the equation can be solved by Implicit method and explicit method:

1)explicit method:

we condider that the value of T is known at all grid point at time level n and hence we can find values of t at n+1 using our time marching solution for given equation.

T(i,j)^(n+1)=T(i,j)+k1*(T(L)-2*T(P)+T(R)/D(x)^2)+k2*(T(T)-2*T(P)+T(B)/D(y)^2)

here k1 and k2 is CFL number

the equation solved directly withou using iterarivre methods usng space loop and convergence loop.

2)Implicit Method:

The unknowns are obtained by means of a simultenous solution of the difference equation applied at all grid points arrayed at given time level.

T(i,j)^(n+1)=(1/1+2k1+2k2)*[T(i,j)^n+k1*(T(L)+T(R))^(n+1)+k2*(T(T)+T(L))^(n+1)]

The equation is solve using iterartive methods:

1) Jacobi Method

2)Gauss Siedel Method

3)SOR Method

1)The steady state 2D heat conduction equation is solved using Iterative techiques:

clear all
close all
clc
alpha=1.1;%relaxation factor
nx=10;
ny=10;
lx=1;
ly=1;
x=linspace(0,lx,nx);
y=linspace(0,ly,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
k=(2*dx^2+2*dy^2)/(dx^2*dy^2);
u=300*ones(nx,ny);
u(1,:)=600;
u(nx,:)=900;
u(:,1)=400;
u(:,ny)=800;
uold=u;
error=1;
tol=1e-4;
iterative_solver=input('iterative_solver=');
if iterative_solver==1;
  jacobi_iter=1;
while (error>tol);
  for i=2:nx-1;
    for j=2:ny-1;
      u(i,j)=(1/k)*(((uold(i+1,j)+uold(i-1,j))/(dx^2))+((uold(i,j+1)+uold(i,j-1))/(dy^2)));
    endfor
   endfor
    error=max(max(abs(uold-u)));
    jacobi_iter=jacobi_iter+1
    uold=u;
    
  end  
  contourf(x,y,u, 'showtext','on')
  colorbar;
  title_text=sprintf('iteration number= %d',jacobi_iter);
     title(title_text)
end

if iterative_solver==2;
  gs_iter=1;
while (error>tol);
  for i=2:nx-1;
    for j=2:ny-1;
      u(i,j)=(1/k)*(((u(i+1,j)+u(i-1,j))/(dx^2))+((u(i,j+1)+u(i,j-1))/(dy^2)));
    endfor
   endfor
    error=max(max(abs(uold-u)));
    gs_iter=gs_iter+1
    uold=u;
    
  end  
  contourf(x,y,u,'showtext','on')
  colorbar;
  title_text=sprintf('iteration number= %d',gs_iter);
     title(title_text)
  
end

if iterative_solver==3;
  sor_iter=1;
while (error>tol);
  for i=2:nx-1;
    for j=2:ny-1;
      u(i,j)=(1-alpha)*uold(i,j)+alpha*((1/k)*(((u(i+1,j)+u(i-1,j))/(dx^2))+((u(i,j+1)+u(i,j-1))/(dy^2))));
    endfor
   endfor
    error=max(max(abs(uold-u)));
    sor_iter=sor_iter+1
    uold=u;
    
  end  
  contourf(x,y,u,'showtext','on')
  colorbar;
  title_text=sprintf('iteration number= %d',sor_iter);
     title(title_text)
end

OUTPUT:

1) Jacobi:

2)Gauss Siedel:

3)SOR:

2) Transient state (Explicit):

clear all
close all
clc
nx=10;
ny=10;
xl=1;
yl=1;
x=linspace(0,xl,nx);
y=linspace(0,yl,ny);
dx=x(2)-x(1);
dy=x(2)-x(1);
u=300*ones(nx,ny);
u(1,:)=600;
u(nx,:)=900;
u(:,1)=400;
u(:,ny)=800;
uold=u;
t=1.4;
dt=0.001;
alpha=1.4;  %thermal diffusivity;
ts=t/dt;
k1= alpha*dt/dx;
k2= alpha*dt/dy;
for n=1:ts
  for i=2:nx-1
    for j=2:ny-1
      a=k1*(uold(i-1,j)-2*uold(i,j)+uold(i+1,j));
      b=k2*(uold(i,j+1)-2*uold(i,j)+uold(i,j-1));
      u(i,j)=uold(i,j)+a+b;     
    endfor
  endfor
  error=max(max(abs(uold-u)))
  uold=u;
endfor
contourf(x,y,u,'showtext','on')
colorbar
xlabel('length')
ylabel('width')

OUTPUT:

3)Transient State (Implicit):

clear all
close all
clc
lx=1;
ly=1;
nx=10;
ny=10;
x=linspace(0,lx,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
u=300*ones(nx,ny);
%boundary conditions
u(1,:)=600;
u(nx,:)=900;
u(:,1)=400;
u(:,ny)=800;
uold=u;
u_prev=u; % u at previsous time step
alpha=1.4;  %thermal diffusivity
dt=0.1;   %time size
t=40;      %time step
ts=t/dt;
omega=1.1;  %relaxation factor 
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=1/(1+2*k1+2*k2);
term2=k1*term1;
term3=k2*term1;
tol=1e-4;
error=9e9;
iterative_solver=input('iterative_solver=');

%jacobi method;
if iterative_solver==1;
  jacobi_iter=1;
  for k=1:ts
    error=1
    while (error>tol)
      for i=2:nx-1
        for j=2:ny-1
          h=(uold(i-1,j)+uold(i+1,j));
          v=(uold(i,j-1)+uold(i,j+1));
          u(i,j)=(u_prev(i,j)*term1)+(h*term2)+(v*term3);
        endfor
      endfor
      error=max(max(abs(uold-u)));
      uold=u;
      jacobi_iter=jacobi_iter+1;
    end
    u_prev=u;
    contourf(x,y,u,'showtext','on')
  colorbar
  title_text=sprintf('iteration number=%d',jacobi_iter)
  title(title_text)
  pause(0.0003)
  end
 end 
 
 %gauss siedel
  if iterative_solver==2;
  gs_iter=1;
  for k=1:ts
    error=1;
    while (error>tol)
      for i=2:nx-1
        for j=2:ny-1
          h=(u(i-1,j)+u(i+1,j));
          v=(u(i,j-1)+u(i,j+1));
          u(i,j)=(u_prev(i,j)*term1)+(h*term2)+(v*term3);
        endfor
      endfor
      error=max(max(abs(uold-u)));
      uold=u;
      gs_iter=gs_iter+1;
    end
    u_prev=u;
    contourf(x,y,u,'showtext','on')
  colorbar
  title_text=sprintf('iteration number=%d',gs_iter)
  title(title_text)
  pause(0.0003)
  end
 end     
 
%sor methods
  if iterative_solver==3;
  sor_iter=1;
  for k=1:ts
    error=1;
    while (error>tol)
      for i=2:nx-1
        for j=2:ny-1
          h=(u(i-1,j)+u(i+1,j));
          v=(u(i,j-1)+u(i,j+1));
          u(i,j)= (1-omega)*uold(i,j)+omega*((u_prev(i,j)*term1)+(h*term2)+(v*term3));
        endfor
      endfor
      error=max(max(abs(uold-u)));
      uold=u;
      sor_iter=sor_iter+1;
    end
    u_prev=u;
    contourf(x,y,u,'showtext','on')
  colorbar;
  title_text=sprintf('iteration number=%d',sor_iter)
  title(title_text)
  pause(0.0003)
  end
  
 end  

OUTPUT:

1)Jacobi:

2)Gauss Siedel:

3)SOR:

Conclusion:

1) For steady state heat conduction we see that a steady state is reached but time taken by Jacobi is more than gauss siedel and Sor as number of iteration is more for jacobi>Gauss siedel>SOR.This is because in jacobi we finf unknowns at each grid point hence simulation time is more whereas in gauss siedel and SOR we update the value from previous time steps. hence solution is reached in lesser no of iteration hence simulation time is also less.

2) Transient state (explicit): we can use iterative methods but we can even directly solve the equation as we have done using time loop space loop and convergence loop. The simulation time is less than implicit methods as we simply solve the equation at known grid points.

3) Transient state (Implicit): in implicit approach the jacobian method takes more no of iteration and simulation time as it takes information both current time ste and later time step hence it increases computational timing as compared to explicit method. but in implicit we see that gauss sidel method takes less number of iterartion than jacobi and SOR takes less number of iterartion as compared to both jacobi and gauss sidel as it takes relaxation factor alpha.