2D Heat Conduction simulation for a plate - For steady and transient condition
2d heat conduction is to be simulated for transient and steady state.
Parameters considered for solving are as follows:
1)For transient time period of 4second is considered.
2)Alpha (diffusivity ) value of 0.05 is considered for all the conditions.
3)Equal number of nodes (10 numbers) are taken along the horizontal as well as vertical edge.
4)Absolute error criteria is taken to be 1e-4
5) The plate boundary conditions are given as
1. Left:400K
2. Top:600K
3. Right:800K
4. Bottom:900K
Equation to be solved:
`T_(i,j) ^(n+1) = T_(i,j) ^(n) + (alpha Deltat *(T_(i+1,j)-2*T_(i,j)+T_(i-1,j))/(Deltax^2)) +( alpha Deltat *(T_(i,j+1)-2*T_(i,j)+T_(i,j-1))/(Deltay^2))`
Here i,j corresponds to the respective rows and columns of the matrixes.
Both explicit and implicit methods can be used to solve the equation.
In explicit method it will be direct substitution of formula where as in implicit solver we will do iterations till the new value converge below the error criteria.
If explicit mehtod is used we can use following methos:
a) Jacobi method
b)Gauss Siedel method
c)Successive over relaxation method
Jacobi Method
Initila values are randomly chosen mostly depending on the boundary conditions.
Calculations is done over all the rows and columns with the values from inital assumption.
When one complete iteration is over for all rows and columns then the initial values are updated.
Speciality of Jacobi method is the values are updated only after one complete iteration of rows and columns.
Gauss Siedel method
Gauss Siedel method is similar to Jacobi method only difference being the values are updated within the current iteration run.
This leads to fewer number of iterations compared to Jacobi method.
Successive over relaxation method
Relaxation factor omega is taken arbitarirly. From previous studies it\'s has shown that it\'s optimum to choos a value below 2.
Here omega is chosen as
omega = 2/(1+sin(pi*dx))
For SOR we can use Jacobi of Gauss Siedel to calculate the change of values from inital time value for one iteration to next time step.
Equation to be solved for explicit equation is:
`T^(n+1) = (1-omega) T^(n) + omega * T_(gauss-side-l)^(n+1)`
Implicit Transient Solution
#2d Heat Conduction using transient state -Implicit method
clear all
close all
clc
prompt = \'Enter 1 to use Jacobi Iteration. Enter 2 to use Gauss-Siedel Iteration. Enter 3 to use SOR with Gauss_Siedel \'
method = input(prompt)
if(method~=1)||(method~=2)||(method~=3)
disp(\'Enter correct values to do iterations\');
end
#diving horizontal and verical domain to equal nodes(divisions)
nx = ny = 10;
#distance between the nodes
dx = 1/(nx-1)
dy = 1/(ny-1)
#diffusivity
alpha = 0.05;
#time stem
dt = 0.2
#relaxation value
omega = 2/(1+sin(pi*dx))
#time of calculation
initial_time =0
final_time=4
#total time of calculation
total_time = final_time-initial_time;
cycle = (final_time-initial_time) /dt
t_space = linspace(initial_time,final_time,cycle)
#constants
k1 = (alpha*dt)/(dx^2);
k2 = (alpha*dt)/(dy^2);
T = 300 * ones(ny,nx);
#boundary conditions
T(1,:) = 600
T(end,:) = 900
T(:,1) = 400
T(:,end) = 800
#Refining the edges
T(1,1) = (600+400)/2
T(1,end) = (600+800)/2
T(end,1) = (900+400)/2
T(end,end) = (900+800)/2
#tolerance allowed
tolerance = 1e-4
error = 9e9;
Told = T;
T_prev_dt = T;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,y);
#Jacobi iteration
if method == 1;
ieration_jacobi_transient =1
#time loop
for t = 1:cycle
#convergence loop
while error > tolerance
for j = 2:ny-1;
for i = 2:nx-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);
endfor
endfor
error = max(max(abs(Told - T)))
Told = T;
ieration_jacobi_transient = ieration_jacobi_transient+1
endwhile
T_prev_dt = T
#forcing time loop to enter convergence loop
error = 1e9
endfor
B = flipud(T)
[c,h] = contourf(xx,yy, B);
colormap(jet);
colorbar;
clabel(c,h);
xlabel(\'X domain\')
ylabel(\'Y domain\')
text = sprintf(\"Implicit-Jacobi Method - Transient -2D heat conduction in a plate after %d seconds using %d x %d nodes.. Iterations Taken = %d\",total_time,nx,ny,ieration_jacobi_transient)
title(text)
pause(0.3)
endif
#gauss sieldel iteration
if method == 2;
iteration_siedel_transient= 1
#time loop
for t = 1:cycle
#convergence loop
while error > tolerance
for j = 2:ny-1;
for i = 2:nx-1;
term1 = (1+(2*k1) + (2*k2))^(-1 );
term2 = k1 * term1;
term3 = k2 * term1;
Told =T;
H = (Told(i-1,j) +Told(i+1,j));
Told =T;
V = (Told(i,j-1) +Told(i,j+1));
T(i,j) = (T_prev_dt(i,j) * term1) + (H*term2) + (V*term3);
T_gauss_one = T;
endfor
endfor
error = max(max(abs(Told - T)))
Told =T;
iteration_siedel_transient = iteration_siedel_transient+1
endwhile
#forcing time loop to enter convergence loop
T_prev_dt = T
error = 1e9
endfor
B = flipud(T)
[c,h] = contourf(xx,yy,B);
colormap(jet);
colorbar;
clabel(c,h);
xlabel(\'X domain\')
ylabel(\'Y domain\')
text = sprintf(\"Implicit-Gauss SiedelMethod - Transient -2D heat conduction in a plate after %d seconds using %d x %d nodes.. Iterations Taken = %d\",total_time,nx,ny, iteration_siedel_transient)
title(text)
endif
#SOR over gauss-siedel iteration
if method == 3;
sor = 1
#time loop
for t = 1:cycle
#convergence loop
while error>tolerance
for j = 2:ny-1;
for i = 2:nx-1;
term1 = (1+(2*k1) + (2*k2))^(-1 );
term2 = k1 * term1;
term3 = k2 * term1;
Told =T;
H = (Told(i-1,j) +Told(i+1,j));
Told =T;
V = (Told(i,j-1) +Told(i,j+1));
T(i,j) = (T_prev_dt(i,j) * term1) + (H*term2) + (V*term3);
endfor
endfor
T_gauss_one = T;
T = (T_prev_dt *(1-omega) + (T_gauss_one )*omega);
error = max(max(abs(T - T_prev_dt)))
T_prev_dt = T;
sor = sor+1
endwhile
#forcing time loop to enter convergence loop
T_prev_dt = T;
error = 1e9
endfor
B = flipud(T)
[c,h] = contourf(xx,yy,B);
colormap(jet);
colorbar;
clabel(c,h);
xlabel(\'X domain\')
ylabel(\'Y domain\')
text = sprintf(\"Implicit-Successive Over Relaxation Method - Transient -2D heat conduction in a plate after %d seconds using %d x %d nodes... Iterations Taken = %d\",total_time,nx,ny,sor)
title(text)
endif
Steady heat conduction using Implicit Method
#2d Heat Conduction using Steady state -Implicit method
clc
clear all
close all
prompt = \'Enter 1 to use Jacobi Iteration. Enter 2 to use Gauss-Siedel Iteration. Enter 3 to use SOR with Gauss_Siedel \';
method = input(prompt)
if(method~=1)||(method~=2)||(method~=3)
disp(\'Enter correct values to do iterations\');
end
#diving horizontal and verical domain to equal nodes(divisions)
nx = ny = 10
#distance between the nodes
dx = dy = 1/(ny-1)
k_steady = 2*(dx^2 + dy^2)/(dx^2 * dy^2)
T = 300 * ones(ny,nx);
#boundary conditions
T(1,:) = 600;
T(end,:) = 900;
T(:,1) = 400;
T(:,end) = 800;
T(1,1) = (600+400)/2;
T(1,end) = (600+800)/2;
T(end,1) = (900+400)/2;
T(end,end) = (900+800)/2;
tolerance = 1e-4
error = 9e9;
Told = T;
T_prev_dt = T
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,y);
#jacobi method
if method == 1;
iterations_steady_jacobi=1
#convergence loop
while error > tolerance
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) = (1/k_steady) * ((H + V)/dx^2 );
#T(i,j) = (1/4) * (H + V) ;
endfor
endfor
error = max(max(abs(Told - T)))
Told = T;
iterations_steady_jacobi = iterations_steady_jacobi + 1
endwhile
B = flipud(T)
[c,h] = contourf(xx,yy, B);;
colormap(jet);
colorbar;
clabel(c,h);
text = sprintf(\"Jacobi Method - Steady y state -2D heat conduction in a plate using %d x %d nodes.. Iterations Taken = %d\",nx,ny,iterations_steady_jacobi)
title(text)
xlabel(\'X domain\')
ylabel(\'Y domain\')
endif
#gauss siedel method
if method == 2;
iterations_steady_siedel=1
#convergence loop
while error > tolerance
for j = 2:ny-1;
for i = 2:nx-1;
Told = T;
H = (Told(i-1,j) +Told(i+1,j));
Told = T;
V = (Told(i,j-1) +Told(i,j+1));
T(i,j) = (1/k_steady) * ((H + V)/dx^2 );
#T(i,j) = (1/4) * (H + V) ;
endfor
endfor
error = max(max(abs(Told - T)))
Told = T;
iterations_steady_siedel = iterations_steady_siedel + 1
endwhile
B = flipud(T)
[c,h] = contourf(xx,yy, B);
colormap(jet);
colorbar;
clabel(c,h);
text = sprintf(\"Gauss Siedel Method - Steady y state -2D heat conduction in a plate using %d x %d nodes.. Iterations Taken = %d\",nx,ny,iterations_steady_siedel)
title(text)
xlabel(\'X domain\')
ylabel(\'Y domain\')
endif
#SOR over gauss siedel
if method == 3;
iterations_steady_sor=1;
omega = 2/(1+sin(pi*dx))
#convergence loop
while error > tolerance
for j = 2:ny-1;
for i = 2:nx-1;
Told = T;
H = (Told(i-1,j) +Told(i+1,j));
Told = T;
V = (Told(i,j-1) +Told(i,j+1));
T(i,j) = (1/k_steady) * ((H + V)/dx^2 );
endfor
endfor
T_gauss_steady = T;
T = (T_prev_dt *(1-omega) + (T_gauss_steady )*omega);
error = max(max(abs(T - T_prev_dt)))
T_prev_dt = T;
iterations_steady_sor = iterations_steady_sor+1
endwhile
B = flipud(T)
[c,h] = contourf(xx,yy, B);
colormap(jet);
colorbar;
clabel(c,h)
text = sprintf(\"SOR METHOD - Seady state -2D heat conduction in a plate using %d x %d nodes.. Iterations Taken = %d\",nx,ny,iterations_steady_sor)
title(text)
xlabel(\'X domain\')
ylabel(\'Y domain\')
endif
Explicit Method
Explicit method is comparetively simple to calculate comparing to implicit methods. In this only the values of previous and current time step is requierd . Explicit mehtod is direct method and simple. But the stability of the explicit method depends on various factors.
Stability of explicit method depends on:
a)distance between nodes
b)time step
c)thermal diffusivity
For stability
`(alpha * Delta t )/(Deltax^2)` should be equal or less than 0.5`
#Explicit method - heat conduction 2d -Transient
clear all
close all
clc
method =1
nx = ny = 10
dx = 1/(nx-1)
dy = 1/(ny-1)
alpha = 0.05;
dt = 0.05;
#check for stability
stable_value = (alpha*dt)/(dx^2)
if stable_value <= 0.5
initial_time =0
final_time=4
total_time = final_time-initial_time
cycle = (final_time-initial_time) /dt
t_space = linspace(initial_time,final_time,cycle);
k1 = (alpha*dt)/(dx^2)
k2 = (alpha*dt)/(dy^2)
T = 400 * ones(ny,nx)
T(1,:) = 600
T(end,:) = 900
T(:,1) = 400
T(:,end) = 800
T(1,1) = (600+400)/2
T(1,end) = (600+800)/2
T(end,1) = (900+400)/2
T(end,end) = (900+800)/2
tolerance = 1e-4
error = 9e9;
Told = T;
T_prev_dt = T
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[xx,yy] = meshgrid(x,y);
if method == 1;
iteration_explicit_transient =1
for t = 1:cycle
for i = 2:nx-1;
for j = 2:ny-1;
term1 = T_prev_dt(i,j)
# term2 = (T_prev_dt(i+1,j)+T_prev_dt(i-1,j) - (2*T_prev_dt(i,j)))*k1
# term3 = (T_prev_dt(i,j+1)+T_prev_dt(i,j-1) - (2*T_prev_dt(i,j)))*k2
term2 = (T(i+1,j)+T(i-1,j) - (2*T(i,j)))*k1;
term3 = (T(i,j+1)+T(i,j-1) - (2*T(i,j)))*k2;
T(i,j)= term1 + (term2) + (term3);
endfor
endfor
#error = max(max(abs(Told - T)))
T_prev_dt = T
iteration_explicit_transient = iteration_explicit_transient+1
endfor
T
B = flipud(T);
[c,h] = contourf(xx,yy, B);;
colormap(jet);
colorbar;
clabel(c,h);
xlabel(\'X domain\')
ylabel(\'Y domain\')
text = sprintf(\"Explicit Method - Transient -2D heat conduction in a plate after %d seconds using %d x %d nodes.. Iterations Taken = %d\",total_time,nx,ny,iteration_explicit_transient)
title(text)
endif
else
(\'The value of alpha or time step is exceeding 0.5 so the solution would not converge. Kindly change the value of either alpha or dt\')
endifResult:
The plate\'s boundary conditions are given. We can observe the heat is conducted through the plate from the boundary regions to the central region. The bottom part near 900K shows the maximum temperature while the region close to 400K shows the lowest temperature when comparing to neighbouring regions.
Explicit methods takes the least number of iterations to solve.
But stability criteria if not fullfilled may lead to unstable/weird results.
Inferences:
a)For transient conditions explicit method takes least number of iterations.
No: iterations:
Explicit_Transient<Succesive_over_relaxation < Gauss_Siedel < Jacobi method
Explicit method has the limitation of stability criteria. So the initial time step has to be taken small.
b) For steady state calculation also
No: iterations:
Succesive_over_relaxation < Gauss_Siedel < Jacobi method