Solving the steady and unsteady 2D heat conduction problem
OBJECTIVE
Solving the steady and unsteady 2D heat conduction problem
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
We will be using iterative methods for solving implicit equations.
1. Jacobi
2. Gauss-seidel
3. Successive over-relaxation
absolute error criteria are 1e-4.
Explicit method
When a computation is solved in terms of a known quantity for an unknown quantity. The stability criterion must be followed while using the explicit method. Steps must be chosen properly for stable solutions.
Implicit method
when a computation is made for a dependent variable are defined by a set of equations and a matrix involved to find an unknown variable. This computation is called the implicit method. This method is unconditionally stable as a result no stability criterion is required to obtain stable results. This method is low in accuracy and high in computational time as it uses simultaneous equations and new values to find approximate values.
2D Heat conduction equation
`(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))`
the time derivative term `(delT)/(delt)` is hyperbolic in nature which is approximated by forward difference scheme and the spatial derivative sare elliptical therefore descritized by central differencing scheme.
the `(delT)/(delt)` term is solve dby time marching method.
Discretization of spatial derivatives is done by central differencing that takes information from neighbouring points to obtain a solution.
The heat conduction equation will be given as
where n is the current time step and n+1 is the next timestep.
i,i+1 are neighbouring points in x-direction and j,j+1 in y-direction.
`Deltax` is grid spacing in the x-direction and `Deltay`in the y-direction.
Iterative solver
1) Jacobi solver
In this method, the temperature at n+1th iteration is found using old values in the nth iteration.
2) Gauss Siedel
n this method, the solver uses n+1 already computed iterations in solving n+1 remaining computations. In this way, old values are overwritten by newer values for iterations. Unlike Jacobi G-S requires one storage vector that needs to be overwritten with each computation. G-S converges faster than Jacobi.
3) Selective Over Relaxation
Solver relaxation factor is used to increase speed and reduce iterations.
Formula:`T^(n+1)=T^n(1-omega)+omega*T(GS) or T(JACOBI)`
`omega` is the relaxation factor. Its values vary from `0<omega<2`
When omega is between 1 to 2 it is known as an over relaxed state.
`T^n` is the value of previous iterations and T(GS) or T(JACOBI) are values for n+1 iteration computed by gauss-seidel or Jacobi iterations.
STEADY-STATE HEAT CONDUCTION EQUATION.
and diffusivity is constant.
Hence the equation becomes
Discretization for steady-state
Solving for T(i,j)
Solving 2-D heat conduction equation implicitly.
Formula for Jacobi-iterative solver
Formula for Gauss-Siedel iterative solver
`(Tp)^n+1=K_1(TL+TR)^n+K_2(T T+TB)^(n+1)`
Formula for Successive-Overrelaxation iterative solver
TRANSIENT STATE EQUATION
`(delT)/(delt)=alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))`
Discritization
Solving transient equation explicitly we get,
`Tp^(n+1)=1/K((1-2.K1-2.K2).Tp^(n)+K1(TL+TR)^(n+1)+K2(T T+TB)^(n+1))`
`K1=alpha.((Deltat)/(Deltax^2))`
`K2=alpha.((Deltat)/(Deltay^2))`
Stability Criteria:
`K1+K2<=0.5` Stability criteria for 2-D heat conduction equation when solving explicitly.
Solving transient equation implicitly we get,
`Tp^(n+1)=1/K(Tp^(n)+K1(TL+TR)^(n+1)+K2(T T+TB)^(n+1))`
`K=1+2.K1+2.K2`
`K1=alpha.((Deltat)/(Deltax^2))`
`K2=alpha.((Deltat)/(Deltay^2))`
Formula for Jacobi-iterative solver
Formula for Gauss-Siedel iterative solver
The formula for Successive-Overrelaxation iterative solver
Main Code for steady-state
clear all
close all
clc
%2D steady state problem
% input 1 for Jacobi 2 for GS 3 for SOR
n=1
nx=10;
ny=10;
l=1 %length of domain
x=linspace(0,l,nx);
y=linspace(l,0,ny)
dx=l/(nx-1);
dy=l/(ny-1);
error=9e9;
tol=1e-4;
%BC
TL=400;
TR=800;
TT=600;
TB=900;
% T matrix
T=300*ones(nx,ny);
T(1,:)=TT;
T(end,:)=TB;
T(:,1)=TL;
T(:,end)=TR;
%edge refinement
T(1,1)=0.5*(900+400);
T(nx,1)=0.5*(900+800);
T(nx,ny)=0.5*(800+600);
T(1,ny)=0.5*(400+600);
% ITERATIVE SOLVERS
Told=T;%(storing values for furthur calcuations)
Tini=T;
% defining K
K=(dx^2*dy^2)/(2*(dx^2+dy^2));
[X,Y]=meshgrid(x,y);
%error, tolerance
error=1;
tol=1e-4;
% JACOBI
jac_iter=0;
if n==1
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=((Told(i+1,j)+Told(i-1,j))/dx^2)
term2=((Told(i,j-1)+Told(i,j+1))/dy^2)
T(i,j)=K*(term1+term2);
end
end
error=max(max(abs(T-Told)));
Told=T;
jac_iter=jac_iter+1;
figure(1)
contourf(flipud(T),'ShowText','ON')
pause(0.005)
title(sprintf('Steady State Jacobi iterations=%d',jac_iter));
colorbar
colormap(jet)
end
end
%GS
if n==2
Told=T;%(storing values for furthur calcuations)
Tini=T;
%error, tolerance
error=1;
tol=1e-4;
gs_iter=0;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=((Told(i+1,j)+T(i-1,j))/dx^2)
term2=((T(i,j-1)+Told(i,j+1))/dy^2)
T(i,j)=K*(term1+term2);
end
end
error=max(max(abs(T-Told)));
Told=T;
gs_iter=gs_iter+1;
figure(2)
contourf(flipud(T),'ShowText','ON')
pause(0.005)
title(sprintf('Steady State GS iterations=%d',gs_iter));
colorbar
colormap(jet)
end
end
%SOR
omega=1.2
if n==3
Told=T;%(storing values for furthur calcuations)
Tini=T;
%error, tolerance
error=1;
tol=1e-4;
sor_iter=0;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=((Told(i+1,j)+T(i-1,j))/dx^2)
term2=((T(i,j-1)+Told(i,j+1))/dy^2)
T(i,j)=K*(term1+term2);
T(i,j)=Told(i,j)*(1-omega)+(omega*T(i,j));
end
end
error=max(max(abs(T-Told)));
Told=T;
sor_iter=sor_iter+1;
figure(3)
contourf(flipud(T),'ShowText','ON')
pause(0.005)
title(sprintf('Steady State GS iterations=%d',sor_iter));
colorbar
colormap(jet)
end
end
EXPLANATION:
1) Specifying 10 grid points in the x and y-direction.
2) Specifying length in x and y domain in our case it is 1 for x and y directions.
3) Defining spatial domain for grid and length.
4) Defining temperature with guess values as 300K for the number of grid points defined.
5) Applying BC in problem.
6) The accuracy is measured with computed solution and compared with tolerance for the desired accuracy.
7) The copy of initial temperatures is made by Told and Tini which is used for evaluating n+1 time step implicitly and explicitly.
9) Calculating K,K1(term1),K2(term2).
10) Told stores the value of the preceding time step.
11) T stores the value for the n+1 step.
12) error is computed and updated with each iteration we input only the single(max) and positive value of error from each computation because of the while loop condition.
Main code for transient implicit analysis
clear all
close all
clc
%2D transient state problem
% input 1 for Jacobi 2 for GS 3 for SOR
n=2
nx=10;
ny=10;
nt=500;
dt=0.001
l=1 %length of domain
x=linspace(0,l,nx);
y=linspace(l,0,ny)
dx=abs(x(1)-x(2));
dy=abs(y(1)-y(2));
%BC
TL=400;
TR=800;
TT=600;
TB=900;
% T matrix
T=300*ones(nx,ny);
T(1,:)=TT;
T(end,:)=TB;
T(:,1)=TL;
T(:,end)=TR;
% thermal diffusivity alpha=1.2
alpha=1.4;
%edge refinement
T(1,1)=0.5*(900+400);
T(nx,1)=0.5*(900+800);
T(nx,ny)=0.5*(800+600);
T(1,ny)=0.5*(400+600);
% ITERATIVE SOLVERS
Told=T;%(storing values for furthur calcuations)
Tini=T;
%k1 and k2
K1=alpha*(dt/(dx^2));
K2=alpha*(dt/(dy^2));
%error,tolerance
error=1;
tol=1e-5;
conv=1;
% JACOBI
jac_iter=0;
jac_t_iter=0;
if n==1
while(conv>tol)
while(error>tol)
for j=2:ny-1
for i=2:nx-1
term1=1/(1+(2*K1)+(2*K2));
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)=(Tini(i,j)*term1)+(h*term2)+(v*term3);
end
end
error=max(max(abs(T-Told))); %calculating error
Told=T;
jac_iter=jac_iter+1;
end
conv=max(max(abs(T-Tini)));
Tini=T; % Updating Tinitial
jac_t_iter=jac_t_iter+1
error=1;
end
figure(1)
contourf(flipud(T),'ShowText','ON')
title(sprintf(' Jacobi iterations=%d n time iterations=%d',jac_iter,jac_t_iter));
colorbar
colormap(jet)
end
% GS
gs_iter=0;
gs_t_iter=0;
if n==2
conv=1;
error=1;
while(conv>tol)
while(error>tol)
for j=2:ny-1
for i=2:nx-1
term1=1/(1+(2*K1)+(2*K2));
term2=K1*term1;
term3=K2*term1;
h=(Told(i+1,j)+T(i-1,j));
v=(T(i,j-1)+Told(i,j+1));
T(i,j)=(Tini(i,j)*term1)+(h*term2)+(v*term3);
end
end
error=max(max(abs(T-Told))); %calculating error
Told=T;
gs_iter=gs_iter+1;
end
conv=max(max(abs(T-Tini)))
gs_t_iter=gs_t_iter+1
error=1;
Tini=T; % Updating Tinitial
end
figure(2)
contourf(flipud(T),'ShowText','ON')
title(sprintf('GS iterations=%d n time iterations=%d',gs_iter,gs_t_iter));
colorbar
colormap(jet)
end
% SOR
sor_iter=0;
sor_t_iter=0;
omega=1.2; % over relaxation factor
if n==3
conv=1;
error=1;
while(conv>tol)
while(error>tol)
for j=2:ny-1
for i=2:nx-1
term1=1/(1+(2*K1)+(2*K2));
term2=K1*term1;
term3=K2*term1;
h=(Told(i+1,j)+T(i-1,j));
v=(T(i,j-1)+Told(i,j+1));
Tgs(i,j)=(Tini(i,j)*term1)+(h*term2)+(v*term3);
T(i,j)=Told(i,j)*(1-omega)+(Tgs(i,j)*omega);
end
end
error=max(max(abs(T-Told))); %calculating error
Told=T;
sor_iter=sor_iter+1;
end
conv=max(max(abs(T-Tini)))
error=1;
sor_t_iter=sor_t_iter+1;
Tini=T; % Updating Tinitial
end
figure(3)
contourf(flipud(T),'ShowText','ON')
title(sprintf('sor iterations=%d n time iterations=%d',sor_iter,sor_t_iter));
colorbar
colormap(jet)
end
The second while loop functions as similar to the first while loop. In the second while loop, it finds the converged value of temperature in a particular time step.
Main code transient explicit analysis
clear all
close all
clc
%2D transient explicit problem
n=1
nx=10;
ny=10;
%time steps
nt=500;
dt=0.001;
l=1 %length of domain
x=linspace(0,l,nx);
y=linspace(l,0,ny)
dx=abs(x(1)-x(2));
dy=abs(y(1)-y(2));
%BC
TL=400;
TR=800;
TT=600;
TB=900;
% T matrix
T=300*ones(nx,ny);
T(1,:)=TT;
T(end,:)=TB;
T(:,1)=TL;
T(:,end)=TR;
% thermal diffusivity alpha=1.2
alpha=1.4;
%edge refinement
T(1,1)=0.5*(900+400);
T(nx,1)=0.5*(900+800);
T(nx,ny)=0.5*(800+600);
T(1,ny)=0.5*(400+600);
% ITERATIVE SOLVERS
Told=T;%(storing values for furthur calcuations)
Tini=T;
%k1 and k2
K1=alpha*(dt/(dx^2));
K2=alpha*(dt/(dy^2));
%error,tolerance
error=1;
tol=1e-4;
iter=0;
% EXplicit solution
for k=1:nt
while(error>tol)
for j=2:ny-1
for i=2:nx-1
term1=(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
term2=(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
T(i,j)=Told(i,j)+(K1*term1)+(K2*term2);
end
end
error=max(max(abs(T-Told))); %calculating error
Told=T;
iter=iter+1;
end
end
figure(1)
contourf(flipud(T),'ShowText','ON')
title(sprintf('Transient explicit iterations=%d',iter));
colorbar
colormap(jet)
CONCLUSION
1) Successive overrelaxation method convergence rate is higher than Jacobi and gauss-seidel in both steady and unsteady implicit methods. The number of steps is less compared to Jacobi and gauss-seidel iterations.
Below table clear an idea of the number of iterations.
| Analysis type | Jacobi | Gauss-Siedel | SOR |
| steady | 207 | 111 | 75 |
| transient | 3716 | 2922 | 2809 |
2)Implicit method can be used when convergence is needed in less computation time as they are unconditionally stable. With an increase in time step the convergence rate is increased which leads to equilibrium in less number of iterations.
3)Explicit method results in time-accurate solutions as it is governed by stability criterion. Hence explicit time rate leads to convergence faster. But if the stability criterion is not satisfied then it may lead to an unstable solution. As timestep size increases solution rate convergence rate decreases and iterations increase.