Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
AIM: Solving the steady and unsteady 2D heat conduction problem
OBJECTIVE:
- Steady-state analysis & Transient State Analysis
Solve the 2D heat conduction equation by using the point iterative techniques that were taught in the class. The Boundary conditions for the problem are as follows;
Top Boundary = 600 K
Bottom Boundary = 900 K
Left Boundary = 400 K
Right Boundary = 800 K
You will implement the following methods for solving implicit equations.
- Jacobi
- Gauss-seidel
- Successive over-relaxation
Your absolute error criteria are 1e-4.
Part 1 -
Steady-state heat conduction
Input conditions to solve the 2D heat conduction at a steady-state using an implicit technique
Length = 1m
Grid points along with x-axis = 100
Grid points along with y-axis = 100
Temperature at left = 400K
Temperature at Right = 500K
Temperature at Top = 600K
Temperature at Bottom = 900K
ωω = 1.25
2D steady state heat conduction equation and its solver
implicit iterative solver
∂2T∂x2+∂2T∂y2=0∂2�∂�2+∂2�∂�2=0
∂2T∂x2=T(i−1,j)+2T(i,j)+T(i+1,j)δx2∂2�∂�2=�(�-1,�)+2�(�,�)+�(�+1,�)��2
∂2T∂y2=T(i,j−1)+2T(i,j)+T(i,j+1)δy2∂2�∂�2=�(�,�-1)+2�(�,�)+�(�,�+1)��2
1 jacobi
T(i,j)=(T(i−1,j)+T(i+1,j)k⋅dx2+T(i,j−1)+T(i,j−1)k⋅dy2)old�(�,�)=(�(�-1,�)+�(�+1,�)�⋅��2+�(�,�-1)+�(�,�-1)�⋅��2)���
Where k=2(∂x2+∂y2)∂x2⋅∂y2�=2(∂�2+∂�2)∂�2⋅∂�2
- Gauss Seidel
T(i,j)=(T(i−1,j)old+T(i+1,j)k⋅dx2+T(i,j−1)old+T(i,j−1)k⋅dy2)�(�,�)=(�(�-1,�)���+�(�+1,�)�⋅��2+�(�,�-1)���+�(�,�-1)�⋅��2)
k=2(∂x2+∂y2)∂x2⋅∂y2�=2(∂�2+∂�2)∂�2⋅∂�2
- Sucessive Over Relexation
T(i,j)=(1−ω)T(i,j)old+ω(T(i−1,j)old+T(i+1,j)k⋅dx2+T(i,j−1)old+T(i,j−1)k⋅dy2)�(�,�)=(1-�)�(�,�)���+�(�(�-1,�)���+�(�+1,�)�⋅��2+�(�,�-1)���+�(�,�-1)�⋅��2)
k=2(∂x2+∂y2)∂x2⋅∂y2�=2(∂�2+∂�2)∂�2⋅∂�2
MATLAB CODE
% 2D heat conduction at steady state using implicit technique
clear all
close all
clc
L=1;%Length of the domain
omega=1.25;% SOR factor
nx=100;%grid size
ny=nx;
x=linspace(0,L,nx);
dx=x(2)-x(1);
y=linspace(0,L,ny);
dy=y(2)-y(1);
T=300*ones(nx,ny);
T(:,1)=400;
T(1,:)=900;
T(:,end)=800;
T(end,:)=600;
T(1,1)=(900+400)/2;
T(1,end)=(900+800)/2;
T(end,1)=(400+600)/2;
T(end,end)=(800+600)/2;
Told=T;
k1=2*(dx^2+dy^2)/(dx^2*dy^2);
tol=1e-4;
error=9e9;
solver=input('Enter the number for the solver method: n Jacobi=1 n Gauss Seidel=2 n SOR=3 n Solver Type: ');
if solver == 1
tic
jacobi_iter=1;
while (error>tol)
for j=2:ny-1
for i=2:nx-1
H=((Told(i-1,j)+Told(i+1,j))/(k1*dx^2));
V=((Told(i,j-1)+Told(i,j+1))/(k1*dy^2));
T(i,j)=H+V;
end
end
error=max(max(abs(Told-T)))
Told=T;
jacobi_iter=jacobi_iter+1
timecounter=toc;
figure(1)
[C,h]=contourf(x,y,T,'ShowText','on');
colorbar
colormap(jet)
xlabel('xx');
ylabel('yy');
clabel(C,h)
title(sprintf('2d Steady State Heat Conduction n Method:Jacobi n No of iterations=%d n Time:%f s',jacobi_iter,timecounter))
end
end
if solver == 2
tic
gs_iter=1;
while (error>tol)
for j=2:ny-1
for i=2:nx-1
H=((T(i-1,j)+Told(i+1,j))/(k1*dx^2));
V=((T(i,j-1)+Told(i,j+1))/(k1*dy^2));
T(i,j)=H+V;
end
end
error=max(max(abs(Told-T)))
Told=T;
gs_iter=gs_iter+1
timecounter_1=toc;
figure(1)
[C,h]=contourf(x,y,T,'ShowText','on');
colorbar
colormap(jet)
xlabel('xx');
ylabel('yy');
clabel(C,h)
title(sprintf('2d Steady State Heat Conduction n Method:Gauss Seidel n No of iterations=%d n Time:%f s',gs_iter,timecounter_1))
end
end
if solver == 3
tic
sor_iter=1;
while (error>tol)
for j=2:ny-1
for i=2:nx-1
H=((T(i-1,j)+Told(i+1,j))/(k1*dx^2));
V=((T(i,j-1)+Told(i,j+1))/(k1*dy^2));
T(i,j)=(Told(i,j)*(1-omega))+(omega*(H+V));
end
end
error=max(max(abs(Told-T)))
Told=T;
sor_iter=sor_iter+1
timecounter_2=toc;
figure(1)
[C,h]=contourf(x,y,T,'ShowText','on');
colorbar
colormap(jet)
xlabel('xx');
ylabel('yy');
clabel(C,h)
title(sprintf('2d Steady State Heat Conduction Method:SOR No of iterations=%d Time:%f s',sor_iter,timecounter_2))
end
end
RESULTS
Part 2 - Transient state heat conduction
Input condition to solve 2D heat conduction under transient state using implicit & explicit technique
Length = 1m
Grid points along with x-axis = 51
Grid points along with y-axis = 51
Temperature at left = 400K
Temperature at Right = 500K
Temperature at Top = 600K
Temperature at Bottom = 900K
ωω = 1.3
α� = 4
2D heat conduction transient - state equation and its solver
∂T∂t=α(∂2T∂x2+∂2T∂y2)∂�∂�=�(∂2�∂�2+∂2�∂�2)
∂2T∂x2=T(i−1,j)+2T(i,j)+T(i+1,j)δx2∂2�∂�2=�(�-1,�)+2�(�,�)+�(�+1,�)��2
∂2T∂y2=T(i,j−1)+2T(i,j)+T(i,j+1)δy2∂2�∂�2=�(�,�-1)+2�(�,�)+�(�,�+1)��2
- Explicit
T(i,j)n=(1−2k1−2k2)⋅T(i,j)n−1+k1⋅(T(i−1,j)+T(i.j−1))n−1+k2⋅(T(i,j−1)+T(i.j+1))n−1�(�,�)�=(1-2�1-2�2)⋅�(�,�)�-1+�1⋅(�(�-1,�)+�(�.�-1))�-1+�2⋅(�(�,�-1)+�(�.�+1))�-1
Where k1=α⋅d2tdx2�1=�⋅�2���2
k2=α⋅d2tdy2�2=�⋅�2���2
- Implicit - jacobi
term1=11+2k1+2k2����1=11+2�1+2�2
term2=k1⋅term1����2=�1⋅����1
term3=k2⋅term1����3=�2⋅����1
T(i,j)n=term1⋅T(i,j)n−1+term2⋅(T(i−1,j)+T(i+1,j))old+term3⋅(T(i,j−1)+T(i,j+1))old�(�,�)�=����1⋅�(�,�)�-1+����2⋅(�(�-1,�)+�(�+1,�))���+����3⋅(�(�,�-1)+�(�,�+1))���
- Implicit - gauss seidel
term1=11+2k1+2k2����1=11+2�1+2�2
term2=k1⋅term1����2=�1⋅����1
term3=k2⋅term1����3=�2⋅����1
T(i,j)n=term1⋅T(i,j)n−1+term2⋅(T(i−1,j)+T(i+1,j)old)+term3⋅(T(i,j−1)+T(i,j+1)old)�(�,�)�=����1⋅�(�,�)�-1+����2⋅(�(�-1,�)+�(�+1,�)���)+����3⋅(�(�,�-1)+�(�,�+1)���)
- Implicit - SOR
term1=11+2k1+2k2����1=11+2�1+2�2
term2=k1⋅term1����2=�1⋅����1
term3=k2⋅term1����3=�2⋅����1
T(i,j)n=(1−ω)T(i,j)old+ω(term1⋅T(i,j)n−1+term2⋅(T(i−1,j)+T(i+1,j)old)+term3⋅(T(i,j−1)+T(i,j+1)old))�(�,�)�=(1-�)�(�,�)���+�(����1⋅�(�,�)�-1+����2⋅(�(�-1,�)+�(�+1,�)���)+����3⋅(�(�,�-1)+�(�,�+1)���))
%% Transient state 2D heat conduction using Explicit technique
clear all
close all
clc
%% intializing variables
Nx = 51;
Ny = 51;
x = linspace (0,1,Nx);
y = linspace (0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 0.0001;
dt = 1;
tol = 1e-4;
%% Temperature gride
T = 298*ones(Nx,Ny);
T(1,:) = 900; %bottom temp
T(end,:) = 600; %Top temp
T(:,1) = 400; %Left temp
T(:,end) = 800; %right temp
%edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
%% CFL number calculation and adaptive Time-step control
CFL_number = (alpha*dt)*((1/dx^2)+(1/dy^2));
%% Calculation of temperature distribution explicitly
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = T;
tic;
counter = 1;
for z = 1:2500
% Explicit Iteration Scheme
for i = 2:Nx-1
for j = 2:Ny-1
term1 = (1-2*k1-2*k2)*Told(i,j);
term2 = k1*(Told(i-1,j)+Told(i+1,j));
term3 = k2*(Told(i,j-1)+Told(i,j+1));
T(i,j) = term1+term2+term3;
end
end
Told = T;
counter = counter+1;
time_counter = toc;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Unstedy state heat conduction Method:Explicit Method Number of iterations: %d(Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f', counter-1,z*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
Result
Transient state heat conduction using implicit method
%% 2D transient state heat conduction equation solved using implicit method
clear all
close all
clc
%% Initializing Variable
Nx = 51;
Ny = 51;
x = linspace (0,1,Nx);
y = linspace (0,1,Ny);
dx = 1/(Nx-1);
dy = 1/(Ny-1);
alpha = 0.0001;
dt = 1;
tol = 1e-4;
%% Temperature gride
T = 298*ones(Nx,Ny);
T(1,:) = 900; %bottom temp
T(end,:) = 600; %Top temp
T(:,1) = 400; %Left temp
T(:,end) = 800; %right temp
%edge refining
T(1,1) = (900+400)/2;
T(1,end) = (900+800)/2;
T(end,1) = (400+600)/2;
T(end,end) = (800+600)/2;
%% CFL number calculation and adaptive Time-step control
CFL_number = (alpha*dt)*((1/dx^2)+(1/dy^2));
%% Calculation of temperature distribution using implicit method
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
[xx,yy] = meshgrid(x,y);
Told = T;
T_prev_dt = T;
term1 = 1/(1+2*k1+2*k2);
term2 = k1*term1;
term3 = k2*term1;
tic;
counter = 1;
solver = input('Enter the number for the solver method: n Jacobi=1 n Gauss Seidal=2 n SOR=3 n Solver type: ');
% Jacobin Method
if solver == 1
for k = 1:2500
error = 9e9;
while(tol<error)
for j = 2:Ny-1
for i = 2:Nx-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)+(H*term2)+(V*term3);
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter+1;
end
T_prev_dt = T;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
time_counter = toc;
title(sprintf('2D Unsteady state heat conduction Method:Implicit Jacobi Number of iterations: %d (Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
end
% Gauss Seidel Method
if solver == 2
for k = 1:2500
error = 9e9;
while (tol<error)
for j = 2:Ny-1
for i = 2:Nx-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)+(H*term2)+(V*term3);
end
end
error = max(max(abs(Told-T)));
Told = T;
counter = counter+1;
end
T_prev_dt = T;
figure (1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
time_counter = toc;
title(sprintf('2D Unsteady state heat conduction Method:Implicit Gauss Seidel Number of iterations: %d (Time: %0.2fs) Computation Time: %0.4f CFL Number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
end
%SOR method
if solver == 3
for k=1:2500
error = 9e9;
while(tol<error)
for j = 2:Ny-1
for i = 2:Nx-1
H= (T(i-1,j)+Told(i+1,j));
V= (T(i,j-1)+Told(i,j+1));
T(i,j) =(1-1.3)*(Told(i,j))+(1.3*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3)));
end
end
error=max(max(abs(Told-T)));
Told = T;
counter = counter+1;
end
T_prev_dt = T;
figure(1)
[C,h] = contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
time_counter = toc;
title(sprintf('2D Unsteady state heat conduction n Method:Implicit SOR n Number of iterations: %d (Time: %0.2fs) n Computation Time: %0.4f n CFL number: %0.4f',counter-1,k*dt,time_counter,CFL_number));
xlabel('xx');
ylabel('yy');
end
end
RESULT
Conclusion
PART 1: Steady-state heat conduction
From the above result, the iteration for convergence of plot is very high for the Jacobi method when compare with gauss- Seidel and SOR. By increasing the grid size we can drastically reduce the iteration which as well gives a convergence plot with stability.
PART 2 - transient- state heat conduction
The above result shows the convergence of transient heat equation using an explicit and implicit technique
- Unsteady state problems can be solved by both implicit and Explicit method
- For stiff equations required lesser timesteps and larger computational time to solve using the explicit method.
Among implicit solvers, SOR gives the result in less time compared to other techniques. It gives stable results irrespective of the CFL number