MID TERM PROJECT - Solving the Steady and Unsteady 2D Heat Conduction Equation
REPORT:
OBJECTIVES:
- Obtain temperature distribution for a thin square plate by 2D heat conduction solution.
- Solve both steady and unsteady 2D heat conduction equation.
- The transient state is solved using both implicit and explicit methods.
- The implicit method requires iterative solvers such as jacobi method,gauss-seidal method and sucessive over relaxation method.
ASSUMPTIONS:
- The domain is a unit square.
- nx = ny [The grid points along the x and y axis are equal].
-
Boundary conditions for steady and transient case
1. Left:400K
2. Top:600K
3. Right:800K
4. Bottom:900K
- Absolute error criteria is 1e-4.
TRANSIENT CASE:
When direct computation of the dependent variables can be made in terms of known quantities, The computation is said to be explicit. When the dependent variable are define by couple of set of equations, an either a matrix or a iterative technique is used to obtain the solution this method is called implicit.
EXPLICIT METHOD:
The forward euler procedure gives the explicit form of the govening equation. By definition in an explicit approach, each difference equation contains only one unknown and can be solved explicitly.
The first step is to discretize the equation. The time derivative is discretized by using the first order forward differencing method, Whereas the spacial terms are discretized by using second order central differencing method.
The value of T at time level \'n+1\' is thus calculated at all grid points by changing the x and y coordinates in the grid. Thus once all values of T are known at time level \'n+1\'. Then we calculate T at next time step. This is carrried out till the solution converges. This method is also called as Time marching.
Since all the values are known in this method the iterative technique is not required.
SOLVING 2D HEAT CONDUCTION EQUATION EXPLICITLY
`frac{delT}{delt} = alpha(frac{del^2T}{delx^2} + frac{del^2T}{dely^2})`
Discretizing the equation.
MATLAB CODE:
clear all
close all
clc
L = 1;
nx = 30;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
tol = 1e-4;
alpha = 1.4;
dt =1/5000;
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_B + T_L)/2;
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 800;
cfl = alpha*(1/(dx^2) + 1/(dy^2))*(dt);
[X,Y] = meshgrid(x,y);
told = T;
tic;
for g = 1:nt
for i = 2:nx-1
for j = 2:ny-1
term1 = told(i,j);
term2 = k1*(told(i+1,j) - 2*told(i,j) + told(i-1,j));
term3 = k2*(told(i,j+1) - 2*told(i,j) + told(i,j-1));
T(i,j) = (term1 + term2 + term3);
end
end
told = T;
end
f = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'No.of iterations = %d \\t time taken = %fs\',nt,f);
title(title_text)
xlabel(\'x axis\')
ylabel(\'y axix\')
set(gca,\'ydir\',\'reverse\')
OUTPUT FOR TRANSIENT STATE EXPLICIT METHOD:
IMPLICIT METHOD:
The backward euler method or the crank-nicolson method gives the implicit form of the equation. In implicit approach, each difference equation contains more than one unknown, Hence it cannot be solved directly. It results in a system of algebric equations for which the unknowns can be solved simultaneously.
Crank-Nicolson form is used to solve the parabolic governing equation ion CFD.In this project , Backward euler method is used for discretization.
The different types of iterative method used are:
- Jacobi method
- Gauss-Seidal method
- SOR method
Jacobi method:
In this method , the unknown on the RHS are assumed to have an initial value, i.e T at \'n+1\' at all nodes is assumed. Using the equation , Then the LHS is calculated, On the next value of (i,j) similar method is applied to calculate LHS. Thus in the first iteration, value of T for time \'n+1\' is calculated at all nodes. The difference between T at \'n+1\' which is obtained and T at \'n+1\' which was assumed previously is taken. If it is greater than a specified tolerance value, Iterations are carried out to calculate T at \'n+1\'. But in the second iteration, T at \'n+1\' is assumed as the value of T which was obtained in the previous iteration.
Once the value of T is found out iteratively whose error is less than the tolerance value. The value of T is calculated at the next time step \'n+2\' using the iterative method.
Gauss-seidal method:
It is similar to jacobo method, But the only difference is for the first iteration whwn the value of T at a particular node is calculated, This updated value at that node is used to calculate the value at next node. But the error term is calculated similar to jacobi method, i.e T at \'n+1\' calculated minus the T at \'n+1\' previously assumed.
Once the value of T is found out iteratively whose error is less than the tolerance value.The value of T is calculated at the next time step \'n+2\' using the same iterative method.
SOR method:
This can be applied to both jacobi and gauss seidal. Here a weighing factor \'w\' is considered due to which the increment in values of T is more rigorous, thereby the number of iterations for convergence of solution is reduced.
Solving 2D heat conductiom equation implicitly:
MATLAB CODE:
clear all
close all
clc
L = 1;
nx = 30;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
tol = 1e-4;
alpha = 1.4;
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_B + T_L)/2;
dt = 1/5000;
cfl = alpha*(1/(dx^2) + 1/(dy^2))*(dt);
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 800;
[X,Y] = meshgrid(x,y);
told = T;
t_prev = told;
iterative_solver = 1:3;
for iterative_solver = 1
tic;
if iterative_solver == 1
jacobi_iter = 1;
for g = 1:nt
error = 9e9;
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(i,j)*term1) + (H*term2) + (V*term3);
end
end
error = max(max(abs(told - T)));
told = T;
jacobi_iter = jacobi_iter + 1;
end
t_prev = T;
end
end
f = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'No.of.iterations(Jacobi) = %d Simulation time = %fs\',jacobi_iter,f);
title(title_text)
xlabel(\'x axis\')
ylabel(\'y axix\')
set(gca,\'ydir\',\'reverse\')
end
L = 1;
nx = 30;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
tol = 1e-4;
alpha = 1.4;
dt = 1/5000;
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_B + T_L)/2;
cfl = alpha*(1/(dx^2) + 1/(dy^2))*(dt);
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 800;
for iterative_solver = 2
tic;
if iterative_solver == 2
gs_iter = 1;
for g = 1:nt
error = 9e9;
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 = T(i-1,j) + told(i+1,j);
V = T(i,j-1) + told(i,j+1);
T(i,j) = (t_prev(i,j)*term1) + (H*term2) + (V*term3);
end
end
error = max(max(abs(told - T)));
told = T;
gs_iter = gs_iter + 1;
end
t_prev = T;
end
end
f2 = toc;
figure(2)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'No,of,iterations(gs) = %d Simulation time = %fs\',gs_iter,f2);
title(title_text)
xlabel(\'x axis\')
ylabel(\'y axix\')
set(gca,\'ydir\',\'reverse\')
end
L = 1;
nx = 30;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
tol = 1e-4;
alpha = 1.4;
omega =1.5;
dt = 1/5000;
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_B + T_L)/2;
cfl = alpha*(1/(dx^2) + 1/(dy^2))*(dt);
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
nt = 800;
for iterative_solver = 3
omega = 1.4;
tic;
if iterative_solver == 3
sor_iter = 1;
for g = 1:nt
error = 9e9;
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 = T(i-1,j) + told(i+1,j);
V = T(i,j-1) + told(i,j+1);
T(i,j) = ((1-omega)*told(i,j)) + omega*((t_prev(i,j)*term1) + (H*term2) + (V*term3));
end
end
error = max(max(abs(told - T)));
told = T;
sor_iter = sor_iter + 1;
end
t_prev = T;
end
end
f3 = toc;
figure(3)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'No.of.iterations(sor) = %d Simulation time = %fs\',sor_iter,f3);
title(title_text)
xlabel(\'x axis\')
ylabel(\'y axix\')
set(gca,\'ydir\',\'reverse\')
end
OUTPUT FOT TRANSIENT STATE IMPLICIT METHOD:
JACOBI METHOD:
GAUSS-SEIDAL METHOD:
SOR METHOD:
STEADY STATE CASE:
In the steady state analysis, The rate of change of Temperature with respect to time is negligible. Both k1 and k2 are equal since nx=ny.
Here in the steady state there is no need to use the iterative techniques as all the values are known.
Solving 2D Heat conduction equation for steady state:
MATLAB CODE:
clear all
close all
clc
L = 1;
nx = 30;
ny = nx;
x = linspace(0,L,nx);
y = linspace(0,L,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
error = 9e9;
tol = 1e-4;
alpha = 1.8;
T_L = 400;
T_T = 600;
T_R = 800;
T_B = 900;
T = ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(nx,2:ny-1) = T_B;
T(1,1) = (T_T + T_L)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_B + T_L)/2;
[X,Y] = meshgrid(x,y)
told = T;
t_prev = told;
tic;
k = 1;
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (0.25*((told(i-1,j) + told(i+1,j) + told(i,j-1) + told(i,j+1))));
end
end
error = max(max(abs(told - T)));
told = T;
k= k+1;
end
f = toc;
figure(1)
[M,N] = contourf(X,Y,T);
clabel(M,N)
colormap(jet)
colorbar
title_text = sprintf(\'No.of iterations = %d \\t Simulation Time = %fs\',k,f);
title(title_text)
xlabel(\'x axis\')
ylabel(\'y axix\')
set(gca,\'ydir\',\'reverse\')
OUTPUT FOR STEADY STATE CASE:
CONCLUSION:
Based on the steady and transient and steady state case are given below:
Transient implicit method requires more computational time compared with transient explicit method and steady state method.
The implicit method is uncontionally stable.If the CFL number is greater than 0.50 the implicit method is preferred as we know high courant number can be used to get stable results in implicit method.
Based on the iterations:
In the transient implicit method
Jacobi method > Gauss seidal method > SOR method
It can be seen that the gauss seidal method is considered over the jacobi since it takes less iteration to converge. Meanwhile SOR method is considered over gauss seidal method since it takes less iterations to converge
When the omega value is 1 the iterations of GS and SOR are the same when the omega value is incresed the SOR iterations are decreased.
The explicit method is considered over implicit since it converges in less iteration and takes less computational time.
Implicit method > Explicit method
The steady state takes more iterations to converges than transient implicit and explicit methods.
Steady state method > Implicit method > Explicit method
Here in the steady state method the iterative techniques are not used since it has all known values.
FOR THE EXPLICIT METHOD:
The time for convergence for an explicit method is very less compared with the implicit method. But the CFL method seems to be unstable when CFL is greater than 0.5.
In explicit method there is no need iterative techniques as all the values are known.
The explicit method solution was obtained by time marching. The convergence criteria was not checked.
It was seen that the explicit method was dominated by COURANT-FRIEDRICHS-LEWY(CFL) condition which determines the stability of the solution.
The CFL condition for 2D heat condition equation is
`alpha(frac{Deltat}{Deltax}) + alpha(frac{Deltat}{Deltax}) le 0.5`
`alpha(frac{Deltat}{(Deltax)^2}) le 0.25` when `Deltax = Deltay`.