Simulation of steady and unsteady conduction
AIM :- To solve 2D conduction equation
conduction is the mode of heat transfer in solids. As the molecules of the solid gains thermal energy they start to vibrate along a mean position and since thy are closely packed they collide with adjacent molecules transfering their kinetic energy. and heat propogates through the solid.
The conduction process is governed by the \'Fourier\'s law of conduction\'. It is another statement of the first law of thermodynamics that is law of conservation of energy. Conduction occurs whenever there is an existance of thermal gradient or temprature difference. As the heat proogates through the solid it sets up a temprature field. that is distribution of temprature along space coordinates.
The conduction process can be mathematically modelled as a partial differential equation which gives the temprature field as the time goes on. the conduction equation is parabolic in nature the conduction could be steady state or transient in nature
steady state conduction : it is a form of conduction in which the temprature differences which are causing the conduction are constant over a period of time and thermal equilibrium has been reached. this condition makes the partial derivative of temprature with resoect to time equal to zero. the partial derivatives with respect to space may or may not have a non zero value
`(del^2T)/(delx^2) +(del^2T)/(dely^2) = 0`
on discreetization we get the following numerical equation
`T_p = [(T_R +T_L)/(deltax^2) + (T_T +T_B)/(deltay^2)] *1/k`
where k = `2(deltax^2 + deltay^2)/(deltax^2*deltay^2)`
T_p = point over which calcultion has to be done
T_R = point in forward direction on grid
T_L = point in backward direction on grid
T_T = point in top direction on grid
T_b = point in bottom direction on grid
unsteady state conduction : Also called transient condition conduction. basically the temprature gradients due to which conduction is occuring are time dependent. as every time the temprature gradient changes the thermal equilibriumis disturbed and the system tends to move to a new equilibrium position. this type of system is more common in nature. The conduction depends on the constant of diffusion of the solid
`(deltaT)/(deltat) + alpha((del^2T)/(delx^2) +(del^2T)/(dely^2)) = 0`
on discreetization the following numerical equation is obtained
`T_P^(n+1) = [T_P^(n) + (k_1(T_L + T_R) +k_2(T_l + T_B))^(n+1)]/(1 + 2k_1 +2k_2)`
where `k_1 = alpha*(deltat)/(deltax^2)`
and `k_2 = alpha*(deltat)/(deltay^2)`
the equations obtained above can be solved by various methods to obtain the temprature field set up in the plate due to the initial conditions
the methods which can be used are:-
- Time marching method : this method can be applied when there is a time dependency of the system as such it can be applied only for the transient condition conduction. this is an explicit type method. in an explicit method each difference equation contains only one unknown and can be solved for this unknown(`T_P^(n+1)`)
- implicit methods :in these methods the difference equations contains multiple unknowns which can be solved simultaneously to get the desired result.(`T_R^(n+1) , T_L^(n+1) , T_T^(n+1) ,T_B^(n+1)`) are the unknowns in this method.
these unknowns can be solved simultaneously by using matrix manipulation techniques like jacobi, gauss siedel, successive over relaxation method
different iterative methods used are:
- jacobi method : an initial guess solution is selected and its value is plugged in. the process is repeated until convergence is reached
- gauss siedel method : a variation of jacobi. The difference is the calculated value in current iteration of previous elements is used for better and faster convergence
- SOR method : Basically weighted gauss method. a factor alpha (>1) is assigned to the gauss method. it results in even faster convergence to the true solution
CODE :
%% for solving the steady state equation input 1
%% for solving the unsteady state equation input 2
clear all
close all
clc
% assuming square domain
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = 1/(nx-1);
dy = 1/(ny-1);
%error and tolerance criteria
error = 9e9;
tol = 1e-4;
[x,y] = meshgrid(x,y);
scheme = input(\'enter 1 for steady or 2 for unsteady\');
if scheme == 1
solver = input(\'jacobi = 1 gausss siedel = 2 SOR = 3\');
T = steady_state(nx,ny,dx,dy,error,tol,x,y,solver)
endif
if scheme == 2
condition = input(\'explicitly = 1 implicitly = 2\');
if condition == 1
T = explicit(nx,ny,dx,dy,x,y);
endif
if condition == 2
solver = input(\'jacobi = 1 gausss siedel = 2 SOR = 3\')
T = implicit(nx,ny,dx,dy,error,tol,x,y,solver);
endif
endif
main code defines the domain and specifies the error and tolerance criteria for the solution. the code asks for user input for the conditions weather transient or steady. depending upon the selection code asks user for implicit or explicit methods of solution and then asks for the type of iterative method to solve.
\'scheme\' is responsible to specify steady or unsteady conditions. for steady state scheme = 1 and for unsteady state scheme = 2
\'condition\' is responsible to specify implicit or explicit methods of calculation. for explicit condition = 1 and for implicit condition = 2
\'solver\' specifies the type of iterative method employed. 1 for jacobi, 2 for gauss seidel, 3 for SOR
Temprature function: defines the boundary and initial conditions. the initial conditions is to be assumed at NTP
function out = temprature_matrix(nx,ny)
T = 300*ones(nx,ny);
T_l = 400;
T_t = 600;
T_r = 800;
T_b = 900;
T(:,end) = T_r;
T(1,:) = T_t;
T(:,1) = T_l;
T(end,:) = T_b;
T(1,1) = 0.5*(T_l+T_t);
T(end,1) = 0.5*(T_l+T_b);
T(1,end) = 0.5*(T_t+T_r);
T(end,end) = 0.5*(T_r+T_b);
out = T;
endfunction
k : constatant used in steady states calculation
%for steady state heat conduction equation
function out = k(dx,dy)
out = 2*(dx^2 +dy^2)/((dx^2)*(dy^2));
endfunction
steady state : based on the user input the function takes in values and calls the respective iterative solver for the steady state condition
%function for solving the steady state 2D conduction equation
function out = steady_state(nx,ny,dx,dy,error,tol,x,y,solver)
scheme = 1;
if solver == 1
out = jacobi(nx,ny,dx,dy,error,tol,x,y,scheme);
endif
if solver == 2;
out = gauss_siedel(nx,ny,dx,dy,error,tol,x,y,scheme);
endif
if solver == 3;
out = sor(nx,ny,dx,dy,error,tol,x,y,scheme);
endif
endfunction
implicit : based on the user input the function takes in values and calls the respective iterative solver for the unsteady state condition
% function for solving the 2D conduction equation implicitly
function out = implicit(nx,ny,dx,dy,error,tol,x,y,solver)
scheme = 2;
if solver == 1 %%for using jacobi method
out = jacobi(nx,ny,dx,dy,error,tol,x,y,scheme);
endif
if solver == 2 %%for using gauss siedel method
out = gauss_siedel(nx,ny,dx,dy,error,tol,x,y,scheme);
endif
if solver == 3 %%for using SOR method
out = sor(nx,ny,dx,dy,error,tol,x,y,scheme)
endif
endfunction
explicit : takes in the required inputs and solves the unsteady equation by a time marching method and plots a temprature distribution curve
%solution by explicit method
function out = explicit(nx,ny,dx,dy,x,y)
% defining the tempature matrix
T = temprature_matrix(nx,ny);
t = T;
t_p = T;
%defining trhe variables
a = 1.48e-4; %thermal diffusivity
dt = 0.5; %time steps
k1 = a*dt/dx^2;
k2 = a*dt/dy^2;
k3 = (1+2*k1+2*k2)^-1 ;
% Since the equation is solved explicitly time marching approach is used
tic
counter = 1;
for p = 1:2000 %time loop
for i = 2:nx-1 %
for j = 2:ny-1
H = (t(i+1,j) + t(i-1,j) - 2*t(i,j));
V = (t(i,j+1) + t(i,j-1) - 2*t(i,j));
T(i,j) = (t(i,j) + k1*H + k2*V);
endfor
endfor
t = T;
out = T;
counter = counter + 1
endfor
iteration_time = toc;
figure(1);
[m n] = contourf(x,y,T);
clabel(m ,n);
colormap(jet);
xlabel (\'xdomain\');
ylabel(\'ydomain\');
title1 = sprintf(\'EXPLICIT SOLUTION OF TRANSIENT EQUATION\');
title2 = sprintf(\'no of iterations = %d\',counter);
title3 = sprintf(\'time elapsed = %f sec\',iteration_time);
colorbar;
title({title1, title2, title3});
endfunction
jacobi : based on the user selection solves either steady or unsteady state equation by jaobi method and then plots a temprature distribution plot
function out = jacobi(nx,ny,dx,dy,error,tol,x,y,scheme)
% defining the temprature matrix
T = temprature_matrix(nx,ny);
t = T;
t_p = T;
% defining the variables
a = 1.45e-4;
dt = 0.5;
k1 = a*dt/dx^2;
k2 = a*(dt/dy^2);
k3 = (1+2*k1+2*k2)^-1 ;
% iteration counter
jacobi_counter = 1;
if (scheme == 1) %for solving the steady state equation
tic
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
H = (t(i+1,j) + t(i-1,j))/dx^2;
V = (t(i,j+1) + t(i,j-1))/dy^2;
T(i,j) = (H + V)/k(dx,dy);
endfor
endfor
error = max(max(abs(T-t)));
t = T;
out = T;
jacobi_counter = jacobi_counter +1;
title1 = sprintf(\'SOLUTION OF STEADY STATE BY JACOBI METHOD\');
end
iteration_time = toc;
end
if (scheme == 2) %for solving the unsteady state equation implicitly
tic
for p = 1:200
while(error > tol)
for i = 2:nx-1
for j = 2:ny-1
H = t(i-1,j) + t(i+1,j);
V = t(i,j-1) + t(i,j+1);
T(i,j) = k3*t(i,j) + k1*k3*H +k2*k3*V;
end
end
error = max(max(abs(T - t)));
t = T;
jacobi_counter = jacobi_counter + 1;
end
t_p = T;
end
iteration_time = toc;
title1 = sprintf(\'SOLUTION OF TRANSIENT STATE IMPLICITLY BY JACOBI METHOD\');
out = T;
end
figure(1);
[m n] = contourf(x,y,T);
clabel(m ,n);
colormap(jet);
xlabel (\'xdomain\');
ylabel(\'ydomain\');
title2 = sprintf(\'no of iterations = %d\',jacobi_counter);
title3 = sprintf(\'time elapsed = %f sec\',iteration_time);
colorbar;
title({title1, title2, title3});
endfunctiongauss siedel : based on the user selection solves either steady or unsteady state equation by gauss siedel method and then plots a temprature distribution plot
function out = gauss_siedel(nx,ny,dx,dy,error,tol,x,y,scheme)
%defining the temprature matrix
T = temprature_matrix(nx,ny);
t = T;
t_p = T;
% defining the variables
a = 1.45e-4;
dt = 0.5;
k1 = a*dt/dx^2;
k2 = a*dt/dy^2;
k3 = (1+2*k1+2*k2)^-1 ;
%iteration counter
counter = 1;
if (scheme == 1) %for steady state conduction
tic
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
H = (T(i-1,j) + t(i+1,j))/dx^2;
V = (T(i,j-1) + t(i,j+1))/dy^2;
T(i,j) = (H+V)/k(dx,dy);
endfor
endfor
error = max(max(abs(T-t)));
t = T;
out = T;
counter = counter + 1
end
iteration_time = toc;
title1 = sprintf(\'STEADY EQUATION IMPLICITLY SOLVED by GAUSS SEIDEL method\');
end
if (scheme == 2 ) %for unsteady state conduction solved implicitly
tic
for p = 1:20000
while(error > tol)
for i = 2:nx - 1
for j = 2:ny - 1
H = T(i-1,j) + t(i+1,j);
V = T(i,j-1) + t(i,j+1);
T(i,j) = k3*t(i,j) + k1*k3*H +k2*k3*V;
end
end
error = max(max(abs(T - t)));
t = T;
out = T;
counter = counter + 1
end
t_p = t;
end
iteration_time = toc;
title1 = sprintf(\'TRANSIENT EQUATION IMPLICITLY SOLVED by GAUSS SEIDEL method\');
end
figure(1);
[m n] = contourf(x,y,T);
clabel(m ,n);
colormap(jet);
title2 = sprintf(\'no of iterations = %d\',counter);
title3 = sprintf(\'time elapsed = %f sec\',iteration_time);
colorbar;
title({title1, title2, title3});
endfunction
SOR : based on the user selection solves either steady or unsteady state equation by SOR method and then plots a temprature distribution plot
function out = sor(nx,ny,dx,dy,error,tol,x,y,scheme)
% defining the temprature matrix
T = temprature_matrix(nx,ny);
t = T;
t_p = T;
% defining the variables
a = 1.45e-4;
dt = 0.5;
k1 = a*dt/dx^2;
k2 = a*(dt/dy^2);
k3 = (1+2*k1+2*k2)^-1 ;
alpha = 1.5;
% number of iterations
sor_counter = 1;
if (scheme == 1 ) % for solving the steady state equation
tic
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
H = (t(i+1,j) + T(i-1,j))/dx^2;
V = (t(i,j+1) + T(i,j-1))/dy^2;
T(i,j) = (1-alpha)*t(i,j)+ alpha*(H + V)/k(dx,dy);
endfor
endfor
error = max(max(abs(T-t)));
t = T;
out = T;
sor_counter = sor_counter +1
end
iteration_time = toc
title1 = sprintf(\'STEADY STATE EQUATION IMPLICITLY SOLVED by SOR method\');
endif
if (scheme == 2) % for solving the unsteady equation implicitly
tic
for p = 1:20000
while(error > tol)
for i = 2:nx - 1
for j = 2:ny - 1
H = T(i-1,j) + t(i+1,j);
V = T(i,j-1) + t(i,j+1);
T(i,j) = (1-alpha)*t(i,j) + alpha*(k3*t(i,j) + k1*k3*H +k2*k3*V);
end
end
error = max(max(abs(T - t)));
t = T;
out = T;
sor_counter = sor_counter +1
end
t_p = T;
end
iteration_time = toc;
title1 = sprintf(\'TRANSIENT EQUATION IMPLICITLY SOLVED by SOR method\');
endif
figure(1);
[m n] = contourf(x,y,T);
clabel(m ,n);
colormap(jet);
title2 = sprintf(\'no of iterations = %d\',sor_counter);
title3 = sprintf(\'time elapsed = %f sec\',iteration_time);
colorbar;
title({title1, title2, title3});
endfunction
plots:
steady state:
- jacobi no of iterations : 208
- gaussno of iterations : 112
- SORno of iterations : 29
unsteady state :
explicit :
no of iterations :2001
implicit:
jacobi :
no of iterations : 6491
gauss :
no of iterations : 6430
SOR :
no of iterations : 4457
CONCLUSION :
- for explicit approach dt is restricted to be equal or less than a certain value given by the stability criterion. above which the solution will go unstable. this makes dt very small and calculations has to be performed over large time interval
- for most implicit methods dt stability can be maintaned over a large value of dt without compromising with the stability of the solution
- but as dt increases in an implicit method so does the truncation error
- generally for calculation time wise implicit>explicit
- for iterative methods sor>gauss>jacobi