Solving the Steady and Transient state 2-D Equation
OBJECTIVE: To solve the steady and unsteady 2D heat conduction equation problem using point iterative techniques.
Conduction: Conduction is a mode of heat transfer between high energy particle and a low energy particle. It takes place in a stationary medium due to the presence of temperature gradient.
Conduction is of two types:
1. Steady state conduction: The temperature within the system does not change with time.
2. Unsteady (or) Transient state conduction: The temperature within the system changes with time.
Given Data:
i) The domain is a unit square.
ii) nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
iii) Boundary conditions for steady and transient case
1. Left:400K
2. Top:600K
3. Right:800K
4. Bottom:900K
iv) Initial condition i.e. T=300K.
Assumptions:
1. No convection taking place
2. No internal heat generation taking place
Problem Solving:
1. Steady state conduction equation will be solved using:
a) Jacobi method
b) Gauss seidel method
c) Successive over-relaxation method
2. Unsteady state conduction equation will be solved using:
a) Implicit method- In this method is solved using:
i) Jacobi method
ii) Gauss seidel method
iii) Successive over-relaxation method
b) Explicit method
The 2D transient heat conduction equation is as follows:
`(delT)/(delt)-alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`---------- equation (1)
1. SOLVING THE STEADY STATE EQUATION:
To obtain the steady state equation, we substitute `(delT)/(delt)=0` in equation (1) then we have
`(del^2T)/(delx^2)+(del^2T)/(dely^2)=0`
Since this is an elliptic profile, we descritize each term using central differencing scheme then we have
`(T_(i+1,j)-2T_(i,j)+T_(i-1,j))^n/(/_\\x^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))^n/(/_\\y^2)=0`
`=>T_(i,j)*(2(/_\\x^2+/_\\y^2)/(/_\\x^2*/_\\y^2))=(T_(i-1,j)+T_(i+1,j))/(/_\\x^2)+(T_(i,j+1)+T_(i,j-1))/(/_\\y^2`
let `k=2*(/_\\x^2+/_\\y^2)/(/_\\x^2*/_\\y^2`
`T_(i,j)=(1/k)*((T_(i-1,j)+T_(i+1,j))/(/_\\x^2)+(T_(i,j+1)+T_(i,j-1))/(/_\\y^2))`
The above equation is our solution equation.
Now representing the equation in jacobi method form, we have:
`T_(i,j)^(n+1)=(1/k)*((T_(i-1,j)+T_(i+1,j))/(/_\\x^2)+(T_(i,j+1)+T_(i,j-1))/(/_\\y^2))^n`
In jacobi method, we solve the equation using the assumed values which are the old values.
Representing the equation in Gauss siedel form, we have:
`T_(i,j)^(n+1)=(1/k)*((T_(i-1,j)+T_(i+1,j))/(/_\\x^2)+(T_(i,j+1)+T_(i,j-1))/(/_\\y^2))^(n+1)`
In Gauss siedel method we solve the equation using the updated values which are the new values.
when we are computing at (i,j) utilising (i+1) and (j+1) value, then their values at nth time step and (n+1)th time step will be the same.
Representing the equation in Successive over relaxation form, we have:
`T_(i,j)^(n+1)=T_(i,j)^n(1-omega)+omega*(1/k)*((T_(i-1,j)+T_(i+1,j))/(/_\\x^2)+(T_(i,j+1)+T_(i,j-1))/(/_\\y^2))^(n+1)`
SOR method is an enhancer for Gauss siedel or jacobi method which utilises the new values for solving the equation.
Here `omega`is the over relaxation factor and is taken as `omega=1.4`
1.1 PROGRAM CODE FOR STEADY STATE HEAT CONDUCTION:
clear all
close all
clc
nx = 10; % No.of nodes in x direction
ny = nx; % No.of nodes in y direction
x= linspace(0,1,nx); % length of x-domain
y = linspace(0,1,ny); % length of y-domain
dx = x(2)-x(1);
dy = y(2)-y(1);
T = 300*ones(nx,ny); % Giving the initial condition i.e. T=300K
% Initialising the boundary conditions
T(:,1) = 400;
T(1,:) = 600;
T(:,ny) = 800;
T(nx,:) = 900;
% Temperature at Corner points
T(1,1) = 500;
T(1,nx) = 700;
T(nx,1) = 650;
T(nx,ny) = 850;
Told = T; % storing T values in Told
dt = 1e-2; % time step
error = 100;
tol = 1e-4;
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
% taking iterative method as input from user
iterative_method=input(\'enter the iterative method = \');
% SOLVING USING JACOBI METHOD
if(iterative_method==1)
tic;
jacobi_iter = 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) = (term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told = T; % Updating values
jacobi_iter = jacobi_iter+1;
time_jacobi = toc;
end
% PLOTTING
figure(1)
contourf(x,y,T,\'ShowText\',\'on\')
colormap(jet);
colorbar;
text1 = sprintf(\'iteration number = %d\',jacobi_iter);
text2 = sprintf(\'simulation time = %f sec\',time_jacobi);
title({\'Temperature Profile of Steady State Heat Conduction Using jacobi method\';text1;text2});
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
end
% SOLVING USING GAUSS SIEDEL METHOD
if(iterative_method==2)
tic;
GS_iter = 1;
while(error>tol)
for i = 2:(nx-1)
for j = 2:(ny-1)
term1= (T(i-1,j)+Told(i+1,j))/(dx^2);
term2= (Told(i,j+1)+T(i,j-1))/(dy^2);
T(i,j)= (term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told = T; % Updating values
GS_iter = GS_iter+1;
time_GS = toc;
end
% PLOTTING
figure(2)
contourf(x,y,T,\'ShowText\',\'on\')
colormap(jet);
colorbar;
text1 = sprintf(\'iteration number = %d\',GS_iter);
text2 = sprintf(\'simulation time = %f sec\',time_GS);
title({\'Temperature Profile of Steady State Heat Conduction Using GS method\';text1;text2});
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
end
% SOLVING USING SOR METHOD
if(iterative_method==3)
tic;
SOR_iter = 1;
omega = 1.4; % Over relaxation factor
while(error>tol)
for i = 2:(nx-1)
for j = 2:(ny-1)
term1= (T(i-1,j)+Told(i+1,j))/(dx^2);
term2= (Told(i,j+1)+T(i,j-1))/(dy^2);
T(i,j) =((Told(i,j)*(1-omega)))+(omega*(term1+term2))/k;
end
end
error=max(max(abs(Told-T)));
Told = T; % Updating values
SOR_iter =SOR_iter+1;
end
time_SOR = toc;
% PLOTTING
figure(3)
contourf(x,y,T,\'ShowText\',\'on\')
colormap(jet);
colorbar;
text1 = sprintf(\'iteration number = %d\',SOR_iter);
text2 = sprintf(\'simulation time = %f sec\',time_SOR);
text3 = sprintf(\'over relaxation factor omega = %f\',omega);
title({\'Temperature Profile of Steady State Heat Conduction Using SOR method\';text1;text2;text3});
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
end
1.2 RESULTS:
1.3 CONCLUSION:
From the above table we see that SOR method gives us excellent results in terms of No.of iterations carried out and the simulation time, as it takes least no.of iterations and less simulation time to converge than jacobi and Gauss seidel methods.This is because of the use of over relaxation factor which accelerates the convergence of the iteration methods.
WHAT ARE IMPLICIT AND EXPLICIT METHODS?
Numerical solution schemes are often referred to as being explicit or implicit. When a direct computation of the dependent variables can be made in terms of known quantities, the computation is said to be explicit. When the dependent variables are defined by coupled sets of equations, and either a matrix or iterative technique is needed to obtain the solution, the numerical method is said to be implicit.
Explicit Scheme: Is one in which the differential equation is discretized in such a way that there is only one unknown (at new time level n+1) on the left hand side (LHS) of the difference equation and it is computed in terms of all other terms on the RHS which are known (at previous time level n) as shown below.
`=>T_(i,j)^(n+1)=T_(i,j)^n+k1*(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))^n+k2*(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))^n`
Implicit Scheme: Is one in which the differential equation is discretized in such a way that there are multiple unknowns at n+1 time level on the LHS of the equation and the terms on RHS are known ones at n time level as shown below.
`=>T_(i,j)^(n+1)-T_(i,j)^n=(alpha*/_\\t)*((T_(i+1,j)-2T_(i,j)+T_(i-1,j))^(n+1)/(/_\\x^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))^(n+1)/(/_\\y^2))`
2. SOLVING THE UNSTEADY STATE EQUATION USING IMPLICIT METHOD:
Let us consider the 2D conduction equation
`(delT)/(delt)-alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
when we descritize the above equation we have:
`(T_(i,j)^(n+1)-T_(i,j)^n)/(/_\\t)-alpha*((T_(i+1,j)-2T_(i,j)+T_(i-1,j))^(n+1)/(/_\\x^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))^(n+1)/(/_\\y^2))=0`
`=>T_(i,j)^(n+1)=T_(i,j)^n+(alpha*/_\\t)*((T_(i+1,j)-2T_(i,j)+T_(i-1,j))^(n+1)/(/_\\x^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))^(n+1)/(/_\\y^2))`
let `k1=(alpha*/_\\t)/(/_\\x^2)`
`k2=(alpha*/_\\t)/(/_\\y^2)`
Now on further simplification we have our implicit equation as:
`=>T_(i,j)^(n+1)=(T_(i,j)^n+k1*(T_(i+1,j)+T_(i-1,j))^(n+1)+k2*(T_(i,j+1)+T_(i,j-1))^(n+1))/(1+2k1+2k2)`
As discussed above we see that in the above equation both LHS and RHS consists of (n+1) terms which are unknown.
Now representing the equation in jacobi method form, we have:
`=>T_(i,j)^(n+1)=(T_(i,j)^n+k1*(T_(i+1,j)+T_(i-1,j))^(n)+k2*(T_(i,j+1)+T_(i,j-1))^(n))/(1+2k1+2k2)`
In jacobi method, we solve the equation using the assumed values which are the old values.
Representing the equation in Gauss siedel form, we have:
`=>T_(i,j)^(n+1)=(T_(i,j)^n+k1*(T_(i+1,j)+T_(i-1,j))^(n+1)+k2*(T_(i,j+1)+T_(i,j-1))^(n+1))/(1+2k1+2k2)`
In Gauss siedel method we solve the equation using the updated values which are the new values.
when we are computing at (i,j) utilising (i+1) and (j+1) value, then their values at nth time step and (n+1)th time step will be the same.
Representing the equation in Successive over relaxation form, we have:
`=>T_(i,j)^(n+1)=T_(i,j)^n(1-omega)+omega*((T_(i,j)^n+k1*(T_(i+1,j)+T_(i-1,j))^(n+1)+k2*(T_(i,j+1)+T_(i,j-1))^(n+1))/(1+2k1+2k2))`
SOR method is an enhancer for Gauss siedel or jacobi method which utilises the new values for solving the equation.
Here `omega`is the over relaxation factor and is taken as `omega=1.2`
2.1 PROGRAM FOR TRANSIENT IMPLICIT:
clear all
close all
clc
nx = 10; % No.of nodes in x direction
ny = nx; % No.of nodes in y direction
nt = 300; % No.of time steps
x= linspace(0,1,nx); % length of x-domain
y = linspace(0,1,ny); % length of y-domain
dx = x(2)-x(1);
dy = y(2)-y(1);
T = 300*ones(nx,ny); % Giving the initial condition i.e. T=300K
% Initialising the boundary conditions
T(:,1) = 400;
T(1,:) = 600;
T(:,ny) = 800;
T(nx,:) = 900;
% Temperature at Corner points
T(1,1) = 500;
T(1,nx) = 700;
T(nx,1) = 650;
T(nx,ny) = 850;
Told = T; % storing T values in Told
T_prev_dt = Told; % storing TOLD values in T_prev_dt
alpha = 0.3; % Initialising Thermal Diffusivity Value
dt = 0.01; % time step
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);
tol = 1e-4;
% taking iterative method as input from user
iterative_method=input(\'enter the iterative method = \');
% SOLVING USING JACOBI METHOD
if(iterative_method==1)
tic;
jacobi_iter = 1;
for p = 1:nt
error=100;
while(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) = ((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
end
end
error = max(max(abs(Told-T)));
Told = T; % Updating old values
jacobi_iter = jacobi_iter + 1;
end
T_prev_dt = T; % Updating previous time step value
end
time_jacobi = toc;
% PLOTTING
figure(1)
contourf(x,y,T,\'ShowText\',\'on\');
colormap(jet);
colorbar;
text2 = sprintf(\' simulation time=%f sec\',time_jacobi);
text3 = sprintf(\'iteration number = %d\',jacobi_iter);
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
title({\'Temperature Profile of Transient State Heat Conduction(Implicit)- jacobi\';text2;text3});
end
% SOLVING USING GAUSS SIEDEL METHOD
if(iterative_method==2)
tic;
GS_iter = 1;
for p =1:nt % Time loop
error = 100;
while(error>tol) % Convergence loop
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) = ((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3));
end
end
error = max(max(abs(Told-T)));
Told = T; % Updating old values
GS_iter = GS_iter + 1;
end
T_prev_dt = T; % Updating previous time step value
end
time_gs = toc;
% PLOTTING
figure(1)
contourf(x,y,T,\'ShowText\',\'on\');
colormap(jet);
colorbar;
text2 = sprintf(\' simulation time=%f sec\',time_gs);
text3 = sprintf(\'iteration number = %d\',GS_iter);
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
title({\'Temperature Profile of Transient State Heat Conduction(Implicit)- GS\';text2;text3});
end
% SOLVING USING SOR METHOD
if(iterative_method==3)
tic;
SOR_iter = 1;
omega = 1.2; % Over relaxation factor
for p = 1:nt % Time loop
error = 100;
while(error>tol) % Convergence loop
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) =(Told(i,j)*(1-omega)) + (omega*((T_prev_dt(i,j)*term1)+(H*term2)+(V*term3)));
end
end
error = max(max(abs(Told-T)));
Told = T; % Updating old values
SOR_iter = SOR_iter+1;
end
T_prev_dt = T; % Updating previous time step value
end
time_sor = toc;
% PLOTTING
figure(1)
text2 = sprintf(\' simulation time=%f sec\',time_sor);
text3 = sprintf(\'iteration number = %d\',SOR_iter);
contourf(x,y,T,\'ShowText\',\'on\');
colormap(jet);
colorbar;
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
title({\'Temperature Profile of Transient State Heat Conduction(Implicit)- SOR\';text2;text3});
end
2.2 RESULTS:
2.3 CONCLUSION:
Similar to the steady state discussion, in transient state also we see that (from the above table) SOR method gives us excellent results in terms of No.of iterations carried out and the simulation time as it takes least no.of iterations and less simulation time to converge than jacobi and Gauss seidel methods.This is because of the use of over relaxation factor which accelerates the convergence of the iteration methods.
3. SOLVING THE UNSTEADY STATE CONDUCTION EQUATION USING EXPLICIT METHOD:
Let us consider the 2D conduction equation
`(delT)/(delt)-alpha((del^2T)/(delx^2)+(del^2T)/(dely^2))=0`
when we descritize the above equation we have:
`(T_(i,j)^(n+1)-T_(i,j)^n)/(/_\\t)-alpha*((T_(i+1,j)-2T_(i,j)+T_(i-1,j))^n/(/_\\x^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))^n/(/_\\y^2))=0`
`=>T_(i,j)^(n+1)=T_(i,j)^n+(alpha*/_\\t)*((T_(i+1,j)-2T_(i,j)+T_(i-1,j))^n/(/_\\x^2)+(T_(i,j+1)-2T_(i,j)+T_(i,j-1))^n/(/_\\y^2))`
let `k1=(alpha*/_\\t)/(/_\\x^2)`
`k2=(alpha*/_\\t)/(/_\\y^2)`
Now on further simplification we have our equation as:
`=>T_(i,j)^(n+1)=T_(i,j)^n+k1*(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))^n+k2*(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))^n`
As discussed above we see that in the above equation both LHS consists of (n+1) term which is unknown and RHS has all known terms.
3.1 PROGRAM FOR TRANSIENT EXPLICIT METHOD:
clear all
close all
clc
nx = 10; % No.of nodes in x direction
ny = nx; % No.of nodes in y direction
nt = 300; % No.of time steps
x = linspace(0,1,nx); % length of x-domain
y = linspace(0,1,ny); % length of y-domain
dx = x(2)-x(1);
dy = y(2)-y(1);
T = 300*ones(nx,ny); % Giving the initial condition i.e. T=300K
% Initialising the boundary conditions
T(:,1) = 400;
T(1,:) = 600;
T(:,ny) = 800;
T(nx,:) = 900;
% Temperature at Corner points
T(1,1) = 500;
T(1,nx) = 700;
T(nx,1) = 650;
T(nx,ny) = 850;
Told = T; % storing T values in Told
alpha = 0.3; % Initialising Thermal Diffusivity Value
dt = 0.01; % time step
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
CFL = alpha*dt/(dx^2);
tic;
exp_iter =1;
for p = 1:nt
for i = 2:(nx-1)
for j = 2:(ny-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
Told = T; % Updating values
exp_iter = exp_iter+1;
end
exp_time = toc;
% PLOTTING
figure(1)
contourf(x,y,T,\'ShowText\',\'on\');
colormap(jet);
colorbar;
text1 = sprintf(\'simulation time = %f sec\',exp_time);
text2 = sprintf(\'Time = %d sec\',dt*nt);
title({\'Temperature Profile of Transient State Heat Conduction(Explicit)\';text1;text2});
xlabel(\'X-axis\')
ylabel(\'Y-axis\')
3.2 RESULTS:
3.3 CONCLUSION:
From the above table we see that for the same time step of 300, explicit method is converging faster than implicit method as it is taking less time for simulation.
If only steady-state results i.e accurate results are wanted, then an implicit solution scheme with over-relaxation should be used so that steady conditions will be reached as quickly as possible.
If, instead time accuracy is important, explicit methods produce greater accuracy with less computational effort than implicit methods.