Solving 2D Heat Conduction Equation Numerically for Transient and Steady Conditions
Solving the Steady and Transient 2D Heat Conduction
The 2D heat conduction equation for Transient case is `(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))`.
For steady-state the equation becomes `((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0`
Both the equations will be solved by the finite differencing method by employing the Jacobi, Gauss-Siedel, and Successive over relaxation scheme. Lets first solve it for steady-state 2D heat conduction.
1. Steady-state 2D Heat Conduction
The equation, `((del^2T)/(delx^2) + (del^2T)/(dely^2)) = 0` can be discretized as `(T_L-2T_P+T_R)/(trianglex^2) + (T_T-2T_P+T_B)/(triangley^2) = 0` which on further solving for `T_P` is `T_P = (1/k){(T_L+T_R)/(trianglex^2) + (T_T+T_B)/(triangley^2)}` where `k = (2(trianglex^2 + triangley^2))/(trianglex^2triangley^2)`. Now converting in to i-j notation, the `T_P` that needs to be updated at every node becomes `T_(i,j) = (1/k){(T_(i-1,j)+T_(i+1,j))/(trianglex^2) + (T_(i,j+1)+T_(i,j-1))/(triangley^2)}` The `T_(i,j)` values can be calculated by any of the above mentioned scheme i.e. Jacobi, Gauss-Siedel and Successive over relaxation.
The script for solving steady-state 2D Heat Conduction is as follow
1. Discretize square domain over n grid-points, initialize the initial temperature over these grid-points incorporating the boundary conditions. Define all other input parameters.
2. Update the new `T_(i,j)` values using different schemes inside the convergence loop until the error reduced below the defined tolerance value. Plot the graph at each iteration to see how the solution converges to the steady-state solution.
Results and Conclusion of steady-state heat conduction by different schemes
The simulations are executed for different schemes for nx = 10 (nx = ny) grid-points. The results will be judged based on the number of iterations taken by different iteration scheme.
1. Using Jacobi iteration scheme
The new `T_(i,j)` value inside the convergence loop is updated as `T_(i,j) = (1/k){(Told_(i-1,j)+Told_(i+1,j))/(trianglex^2) + (Told_(i,j+1)+Told_(i,j-1))/(triangley^2)}`
2. Using Gauss-Siedel iteration scheme
The new `T_(i,j)` value inside convergence loop is updated as `T_(i,j) = (1/k){(T_(i-1,j)+Told_(i+1,j))/(trianglex^2) + (Told_(i,j+1)+T_(i,j-1))/(triangley^2)}`. In this case while updating the value at any grid-point, the value of the previous updated points are used.
3. Using Successive over relaxation (SOR) iteration scheme
The new `T_(i,j)` value inside the convergence loop is updated as `T_(i,j) = (1-w)Told_i,j + w(Tgauss)` or `T_(i,j) = (1-w)Told_i,j + w(Tjacobi)`.
The plot using the different schemes are attached below
Results
1. Iteration taken by different iteration schemes for error below specified tolerance of 1e-4 are as follows
SOR (with Jacobi for w = 0.9) > Jacobi > SOR (with Gauss-Siedel for w = 0.9) > Gauss-Siedel > SOR (with Gauss-Siedel for w = 0.9)
2. SOR with Jacabi iteration for w value above 1 or for over relaxed values gives unstable results.
Conclusion
| Steady-state | ||||||
| Jacobi | Gauss-Siedel | SOR | ||||
| w = 0.9 | w = 1.1 | |||||
| with Jacobi | with Gauss-Siedel | with Jacobi | with Gauss-Siedel | |||
| Total iterations | 207 | 111 | 228 | 135 | Unstable | 92 |
| Simulation time | 65.81 | 34.84 | 68.2 | 36.56 | 29.05 | |
| Error | 9.54E-05 | 9.63E-05 | 9.93E-05 | 9.06E-05 | 8.64E-05 |
2. Transient 2D Heat Conduction with Explicit Method
The equation, `(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))` can be discretized as `(T_P^(n+1)-T_P^n)/(trianglet) = alpha((T_L-2T_P+T_R)/(trianglex^2) + (T_T-2T_P+T_B)/(triangley^2))^n` which on further solving for n+1 time-step is `T_P^(n+1) = T_P^n(1-2k1-2k2)+k1(T_L+T_R) + k2(T_T+T_B)` where `k1 = (alphatrianglet)/(trianglex^2)` and `k1 = (alphatrianglet)/(triangley^2)`. Now converting in to i-j notation, the `T_P` that needs to be updated at every node becomes `T_(i,j)^(n+1) = T_(i,j)^n(1-2k1-2k2)+k1(T_(i-1,j)+T_(i+1,j)) + k2(T_(i,j-1)+T_(i,j+1))`. The `T_(i,j)` values than can be calculated by different iteration scheme such as Jacobi, Gauss-Siedel and Successive over relaxation.
The script for solving Transient 2D Heat Conduction using Explicit Method is as follow
1. Discretize square domain over n grid-points, initialize the initial temperature over these grid-points incorporating the boundary conditions. Define all other input parameters such as time as total time, time-step, alpha values, etc.
2. Update the new `T_(i,j)` values using different iteration schemes inside the space loop. Plot the graph inside each iteration of time-loop.
Results and Conclusion of Transient 2D Conduction with Explicit Method
The simulations are executed for different schemes for nx = 10 (nx = ny) grid-points. Other values are
alpha = 1e-2 % Thermal diffusivity
t_total= 120 sec
dt = 0.25 % time-step
N = t_total/dt % total number of iterations
tolerance = 1e-4, for comparing the error
The results will be judged based on the number of iterations taken by different iteration schemes for reaching an error below a specified limit instead of going to the total iterations. in this way, the result will be compared to see which requires less iteration. For comparing the error below tolerance value, the absolute difference between the two iteration mid-point temperature is calculated and compared.
1. Using Jacobi iteration scheme
The new `T_(i,j)` value inside the space loop is updated as `T_(i,j)^(n+1) = Told_(i,j)^n(1-2k1-2k2)+k1(Told_(i-1,j)^n+Told_(i+1,j)^n) + k2(Told_(i,j-1)^n+Told_(i,j+1)^n)`
2. Using Gauss-Siedel iteration scheme
The new `T_(i,j)` value inside space loop is updated as `T_(i,j)^(n+1) = Told_(i,j)^n(1-2k1-2k2)+k1(T_(i-1,j)^n+Told_(i+1,j)^n) + k2(T_(i,j-1)^n+Told_(i,j+1)^n)`. In this case while updating the value at any grid-point, the value of the previous updated points are used.
3. Using Successive over relaxation (SOR) iteration scheme
The new `T_(i,j)` value inside the space loop is updated as `T_(i,j) ^(n+1)= (1-w)Told_i,j^n + w(Tgauss)^n` or `T_(i,j)^(n+1) = (1-w)Told_i,j^n + w(Tjacobi)^n`.
The plot using the different schemes are attached below
Results
Iteration taken by different iteration schemes for error (difference in midpoint temperature for two consecutive iterations) below specified tolerance of 1e-4 are as follows:
1. Based on iteration taken by time loop, the iterations taken by different techniques are as:
SOR (with Jacobi for w = 0.9) > Jacobi > SOR (with Jacobi for w = 1.1) > SOR (with Gauss-Siedel for w = 0.9) > Jacobi > SOR (with Gauss-Siedel for w = 0.9)
2. SOR method with under relaxed values results in more number of iterations than over relaxed values.
3. Mid-point temperature plot, shows us how the simulation approaches the steady state and stop the time loop once the simulation approaches the steady state.
Conclusion
| Transient Explicit | ||||||
| Jacobi | Gauss-Siedel | SOR | ||||
| w = 0.9 | w = 1.1 | |||||
| with Jacobi | with Gauss-Siedel | with Jacobi | with Gauss-Siedel | |||
| Total iterations | 252 | 159 | 278 | 187 | 231 | 136 |
| Simulation time | 100.84 | 62.37 | 112.32 | 75.87 | 91.87 | 53.38 |
| Error | 9.85E-05 | 9.87E-05 | 9.96E-05 | 9.59E-05 | 9.49E-05 | 9.86E-05 |
3. Transient 2D Heat Conduction with Explicit Method
The equation, `(delT)/(delt) = alpha((del^2T)/(delx^2) + (del^2T)/(dely^2))` can be discretized as `(T_P^(n+1)-T_P^n)/(trianglet) = alpha((T_L-2T_P+T_R)/(trianglex^2) + (T_T-2T_P+T_B)/(triangley^2))^(n+1)` which on further solving for n+1 time-step is `T_P^(n+1) = A*T_P^n+B*(T_L+T_R)^(n+1) + C*(T_T+T_B)^(n+1)` where `A = 1/(1+2k1+2k2)`, `B = A*k1` and `C = A*k2`. Now converting in to i-j notation, the `T_P` that needs to be updated at every node becomes `T_(i,j)^(n+1) = A*Tprev_(i,j)^n+B(Told_(i-1,j)+Told_(i+1,j)) + C(Told_(i,j-1)+Told_(i,j+1))`. The `T_(i,j)` values than can be calculated by different iteration scheme such as Jacobi, Gauss-Siedel and Successive over relaxation.
The script for solving Transient 2D Heat Conduction using Explicit Method is as follow
1. Discretize square domain over n grid-points, initialize the initial temperature over these grid-points incorporating the boundary conditions. Define all other input parameters such as time as total time, time-step, alpha values, etc.
2. Update the new `T_(i,j)` values using different iteration schemes inside the space convergence loop.
3. Update the new `T_(i,j)` value for each new iteration in the time loop. Plot the graph inside each iteration of time-loop until the error ( difference in midpoint temperature in two successive iterations) below specified tolerance value.
Results and Conclusion of Transient 2D Conduction with Explicit Method
The simulations are executed for different schemes for nx = 10 (nx = ny) grid-points. Other values are
alpha = 1e-2 % Thermal diffusivity
t_total= 120 sec
dt = 0.25 % time-step
N = t_total/dt % total number of iterations
tolerance = 1e-4, for comparing the error
The results will be judged based on the number of iterations taken by different iteration schemes for reaching an error below a specified limit instead of going to the total iterations. in this way, the result will be compared to see which requires less iteration. For comparing the error below tolerance value, the absolute difference between the two iteration mid-point temperature is calculated and compared.
1. Using Jacobi iteration scheme
The new `T_(i,j)` value inside the convergence loop is updated as `T_(i,j)^(n+1) = A*Tprev_(i,j)^n+B(Told_(i-1,j)+Told_(i+1,j)) + C(Told_(i,j-1)+Told_(i,j+1))`
2. Using Gauss-Siedel iteration scheme
The new `T_(i,j)` value inside convergence loop is updated as `T_(i,j)^(n+1) = A*Tprev_(i,j)^n+B(T_(i-1,j)+Told_(i+1,j)) + C(T_(i,j-1)+Told_(i,j+1))`. In this case while updating the value at any grid-point, the value of the previous updated points are used.
3. Using Successive over relaxation (SOR) iteration scheme
The new `T_(i,j)` value inside the convergence loop is updated as `T_(i,j) ^(n+1)= (1-w)Told_i,j^n + w(Tgauss)^n` or `T_(i,j)^(n+1) = (1-w)Told_i,j^n + w(Tjacobi)^n`.
The plot using the different schemes are attached below
Results
Iteration taken by different iteration schemes for error (difference in midpoint temperature for two consecutive iterations) below specified tolerance of 1e-4 are as follows:
1. Based on iteration taken by convergence loop, the iterations taken by different techniques are as:
SOR (with Jacobi for w = 0.9) > Jacobi > SOR (with Jacobi for w = 1.1) > SOR (with Gauss-Siedel for w = 0.9) > Jacobi > SOR (with Gauss-Siedel for w = 0.9)
2. SOR method with under relaxed values results in more number of iterations than over relaxed values.
3. The mid-point temperature plot shows us how the simulation approaches the steady-state and stop the time loop once the simulation approaches the steady-state.
Conclusion
| Transient Implicit | ||||||
| Jacobi | Gauss-Siedel | SOR | ||||
| w = 0.9 | w = 1.1 | |||||
| with Jacobi | with Gauss-Siedel | with Jacobi | with Gauss-Siedel | |||
| Total iterations time loop | 259 | 260 | 258 | 259 | 259 | 261 |
| Total iterations space loop | 2138 | 1605 | 2406 | 1905 | 1909 | 1350 |
| Time to reach steady state | 64.75 | 65 | 64.5 | 64.75 | 64.75 | 65.25 |
| Simulation time | 135.73 | 97.54 | 133.61 | 96.52 | 139.32 | 145.93 |
| Error | 9.69E-05 | 9.74E-05 | 9.93E-05 | 9.82E-05 | 9.88E-05 | 9.72E-05 |
Below is the attached code for the above problem
%%Week_5_Mid_Term_Project
%Solving the steady and unsteady 2D heat conduction problem
clear all
close all
clc
% Inputs
Lx = 1; % domain length in x-direction
Ly = 1; % domain length in y-direction
nx = input(\'Number of grid-points along x and y = \');
ny = nx;
% Grid-point calculation
x = linspace(0,1,nx);
y = linspace(0,1,ny);
[X,Y] = meshgrid(x,y);
dx = x(2) - x(1);
dy = dx;
% Assued temperature along grid-points
T = 298*ones(nx,ny);
% Boundary conditions
T(2:nx-1,1) = 400;
T(1,2:ny-1) = 600;
T(2:nx-1,end) = 800 ;
T(end,2:ny-1) = 900;
T(1,1) = (T(1,2) + T(2,1))/2;
T(1,end) = (T(1,end-1) + T(2,end))/2;
T(end,1) = (T(end,2) + T(end-1,1))/2;
T(end,end) = (T(end,end-1) + T(end-1,end))/2;
Told = T;
Tprev = T;
% tolerance and error values
tol = 1e-4;
error = 9e9;
relerr = 9e9
solution = input(\'Enter, 1 for Steady state,2 for Transient = \');
if solution == 2;
solver = ...
input(\'Enter 1 for Transient Explicit or 2 for Transient Implicit = \');
end
iterative_solver = input(\'Enter, 1 for Jacobi,2 for Gauss-siedel, 3 for SOR = \');
if iterative_solver == 3;
% SOR omega value
w = input(\'alpha value for SOR method = \');
solution_solver = input(\'Enter 1 with Jacobi or 2 with Gauss-Siedel = \');
end
%% steady state
if solution ==1;
k = 2*[dx^2 + dy^2]/[dx^2*dy^2];
% Jacobi Method
if iterative_solver == 1;
Jacobi_iter = 1;
tic
while error > tol;
for i = 2:nx-1;
for j = 2:ny-1;
term_1 = (Told(i-1,j) + Told(i+1,j))/(k*dx^2);
term_2 = (Told(i,j-1) + Told(i,j+1))/(k*dy^2);
T(i,j) = term_1 + term_2;
end
end
error = max(max(abs((T-Told))));
Told = T;
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
set(gca,\'Ydir\',\'reverse\');
xlabel(\'x\');
ylabel(\'y\');
l = sprintf(\'2D Steady-state Conduction Temperature contour, Jacobi solver\');
title_text = ({l,sprintf(\'Iteration number = %d, error = %.2d\', Jacobi_iter,error)});
title(title_text);
M(Jacobi_iter) = getframe(gcf); % Capturing frame to make movie
pause(.05);
Jacobi_iter = Jacobi_iter + 1;
end
total_Jacobi_iter = Jacobi_iter-1
Simulation_time = toc
error
% Movie Creation
saveas(gca,\'2D steady state conduction Temperature contour, Jacobi solver.jpg\');
movie(M);
save 2D_steady_state_conduction_Temperature_contour_Jacobi_solver.avi;
videofile = VideoWriter(\'2D_steady_state_conduction_Temperature_contour_Jacobi_solver.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
% Gause-Siedel
GS_iter = 1;
if iterative_solver == 2;
tic
while error > tol
for i = 2:nx-1;
for j = 2:ny-1
term_1 = (T(i-1,j) + Told(i+1,j))/(k*dx^2);
term_2 = (T(i,j-1) + Told(i,j+1))/(k*dy^2);
T(i,j) = term_1 + term_2;
end
end
error = max(max(abs((T-Told))));
Told = T;
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
set(gca,\'Ydir\',\'reverse\');
xlabel(\'x\');
ylabel(\'y\')
l=sprintf(\'2D Steady-state conduction Temperature contour, Gause-Siedel solver\');
title_text = ({l,sprintf(\'Iteration number = %d, error = %.2d\', GS_iter,error)});
title(title_text);
M(GS_iter) = getframe(gcf); % Capturing frame to make movie
pause(0.05);
GS_iter = GS_iter + 1;
end
total_GS_iter = GS_iter-1
Simulation_time = toc
error
% Movie Creation
saveas(gca,\'2D steady state conduction Temperature contour, Gause-Siedel solver.jpg\')
movie(M);
save 2D_steady_state_conduction_Temperature_contour_Gause-Siedel_solver.avi;
videofile = VideoWriter(\'2D_steady_state_conduction_Temperature_contour_Gause-Siedel_solver.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
% SOR
SOR_iter = 1;
if iterative_solver == 3;
tic
while error > tol;
for i = 2:nx-1;
for j = 2:ny-1;
if solution_solver == 1;
term_1 = (Told(i-1,j) + Told(i+1,j))/(k*dx^2);
term_2 = (Told(i,j-1) + Told(i,j+1))/(k*dy^2);
end
if solution_solver == 2;
term_1 = (T(i-1,j) + Told(i+1,j))/(k*dx^2);
term_2 = (T(i,j-1) + Told(i,j+1))/(k*dy^2);
end
T(i,j) = (1-w)*Told(i,j) + w*(term_1 + term_2);
end
end
error = max(max(abs((T-Told))));
Told = T;
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
set(gca,\'Ydir\',\'reverse\');
xlabel(\'x\');
ylabel(\'y\');
if solution_solver == 1;
l = sprintf(\'2D Steady-state conduction Temperature contour,SOR with Jacobi solver\');
else
l = sprintf(\'2D Steady-state conduction Temperature contour, SOR with Gauss-Siedel solver\');
end
title_text = ({l,sprintf(\'Iteration number = %d, w = %.1f, error = %.2d\', SOR_iter, w, error)});
title(title_text);
set(gcf,\'position\',[200,100,700,500]);
M(SOR_iter) = getframe(gcf); % Capturing frame to make movie
pause(0.005);
SOR_iter = SOR_iter + 1;
end
total_SOR_iter = SOR_iter-1
Simulation_time = toc
error
% Movie Creation
if solution_solver == 1;
saveas(gca,\'2D steady state conduction Temperature contour, SOR solver using Jacobi.jpg\');
movie(M);
save 2D_steady_state_conduction_Temperature_contour_SOR_solver_using_Jacobi.avi;
videofile = VideoWriter(\'2D_steady_state_conduction_Temperature_contour_SOR_solver_using_Jacobi.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
else
saveas(gca,\'2D steady state conduction Temperature contour, SOR solver using Gause-Siedel.jpg\')
movie(M);
save 2D_steady_state_conduction_Temperature_contour_SOR_solver_using_Gause-Siedel.avi;
videofile = VideoWriter(\'2D_steady_state_conduction_Temperature_contour_SOR_solver_using_Gause-Siedel.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
end
end
%% Transient
alpha = 1e-2;
dt = .25;
t_total = 2*60;
N = t_total/dt;
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
%% Explicit Solver
if solver == 1;
if k1 > .25;
sprintf(\'Unstable Solution\');
else
iter = 1;
error = 9e9;
tic
while iter<= N & error > tol;
% Jacobi Method
if iterative_solver == 1;
for i = 2:nx-1;
for j = 2:ny-1
TL = Told(i-1,j);
TR = Told(i+1,j);
TB = Told(i,j-1);
TT = Told(i,j+1);
T(i,j) = (1 -2*k1 -2*k2)*Told(i,j)...
+ k1*(TL + TR) + k2*(TB + TT);
end
end
Told = T;
end
% Gause-Siedel
if iterative_solver == 2;
for i = 2:nx-1;
for j = 2:ny-1;
TL = T(i-1,j);
TR = Told(i+1,j);
TB = T(i,j-1);
TT = Told(i,j+1);
T(i,j) = (1 -2*k1 -2*k2)*Told(i,j)...
+ k1*(TL + TR) + k2*(TB + TT);
end
end
Told = T;
end
% SOR
if iterative_solver == 3;
for i = 2:nx-1;
for j = 2:ny-1;
if solution_solver == 1;
TL = Told(i-1,j);
TB = Told(i,j-1);
end
if solution_solver == 2;
TL = T(i-1,j);
TB = T(i,j-1);
end
TR = Told(i+1,j);
TT = Told(i,j+1);
T(i,j) = (1-w)*Told(i,j) + w*((1 -2*k1 -2*k2)*Told(i,j)...
+ k1*(TL + TR) + k2*(TB + TT));
end
end
Told = T;
end
aa(iter) = iter;
%time_to_reach_steady_state
tt = dt*iter;
if rem(nx,2) == 0;
Temp(iter) = (T(nx/2,ny/2) + T(nx/2+1,ny/2+1))/2;
T_prev = (Tprev(nx/2,ny/2) + Tprev(nx/2+1,ny/2+1))/2;
else
Temp(iter) = T(nx/2,ny/2);
T_prev = (Tprev(nx/2,ny/2) + Tprev(nx/2+1,ny/2+1))/2;
end
if iter > 50; % inorder to give an error time to grow
error(iter) = max(abs(Temp(iter)-T_prev));
else
error(iter) = 9e9;
end
l = sprintf(\'Transient 2D Heat Conduction Temperature contour\');
p = sprintf(\'nx=ny = %.0f,alpha = %0.3f, dt = %0.3f, Total time = %1.0f, \\n CFLx = %0.3f, CFLy = %0.3f\',nx, alpha,dt,t_total,k1,k2);
if iterative_solver == 1;
q = sprintf(\'Explicit Jacobi solver, Iteration number = %d\',iter);
q1 = sprintf(\'Explicit Jacobi solver \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\', aa(iter),tt,error(iter));
end
if iterative_solver == 2;
q = sprintf(\'Explicit Gauss-Siedel solver, Iteration number = %d\',iter);
q1 = sprintf(\'Explicit Gauss-Siedel solver \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\', aa(iter),tt, error(iter));
end
if iterative_solver == 3;
if solution_solver == 1;
q = sprintf(\'Explicit SOR with Jacobi,w = %0.2f, Iteration number = %d\',w,iter);
q1 = sprintf(\'Explicit SOR with Jacobi,w = %0.2f \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\',w , aa(iter),tt, error(iter));
end
if solution_solver == 2;
q = sprintf(\'Explicit SOR with Gauss-Siedel,w = %0.2f, Iteration number = %d\',w,iter);
q1 = sprintf(\'Explicit SOR with Gauss-Siedel,w = %0.2f \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\',w , aa(iter),tt, error(iter));
end
end
figure(1);
plot(aa,Temp,\'.\');
xlabel(\'Iteration No.\');
ylabel(\'Mid-point Tempearature\');
title_text = {l,q1};
title(title_text);
Tprev = T;
figure(2);
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
set(gca,\'Ydir\',\'reverse\');
xlabel(\'x\');
ylabel(\'y\');
title_text = {l,p,q};
title(title_text);
M(iter) = getframe(gcf); % Capturing frame to make movie
pause(0.05);
iter = iter + 1;
end
total_iter = iter-1
error(total_iter)
Simulation_time = toc
% Movie Creation
if iterative_solver == 1;
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Explicit Jacobi solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour explicit, Jacobi solver.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_explicit_Jacobi_solver.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_explicit_Jacobi_solver.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
if iterative_solver == 2;
% Movie Creation
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Emplicit Gauss-Siedel solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour explicit, Gauss-Siedel solver.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_explicit_Gauss-Siedel_solver.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_explicit_Gauss-Siedel_solver.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
if iterative_solver == 3;
if solution_solver == 1;
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Explicit SOR with Jacobi solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour explicit, SOR with Jacobi.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_explicit_SOR_with_Jacobi.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_explicit_SOR_with_Jacobi.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
else
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Implicit SOR with Gauss-Siedel solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour explicit, SOR with Gauss-Siedel.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_explicit_SOR_with_Gauss-Siedel.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_explicit_SOR_with_Gauss-Siedel.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
end
end
end
%% Implicit Solver
if solver == 2;
term_1 = 1/(1 + 2*k1 + 2*k2);
term_2 = k1*term_1;
term_3 = k2*term_1 ;
iter = 1;
err = 9e9;
counter_1 = 0;
counter = 0;
tic
while iter <= N & error > tol;
% Jacobi
if iterative_solver == 1;
relerr = 9e9;
while(relerr>tol);
for i=2:nx-1;
for j=2:ny-1;
H = (Told(i-1,j)+Told(i+1,j));
V = (Told(i,j-1)+Told(i,j+1));
T(i,j) = Tprev(i,j)*term_1 + H*term_2 + V*term_3;
end
end
counter = counter+1;
relerr = max(max(abs(Told-T)));
Told = T;
end
end
% Gauss-Siedel
if iterative_solver == 2;
relerr = 9e9;
while(relerr>tol)
for i=2:nx-1;
for j=2:ny-1;
H = (T(i-1,j)+Told(i+1,j));
V = (T(i,j-1)+Told(i,j+1));
T(i,j) = Tprev(i,j)*term_1 + H*term_2 + V*term_3;
end
end
counter = counter+1;
relerr = max(max(abs(Told-T)));
Told = T;
end
end
% SOR
if iterative_solver == 3;
relerr = 9e9;
while(relerr>tol);
for i=2:nx-1;
for j=2:ny-1;
if solution_solver == 1;
H1 = Told(i-1,j);
V1 = Told(i,j-1);
end
if solution_solver == 2;
H1 = T(i-1,j);
V1 = T(i,j-1);
end
H = (H1 + Told(i+1,j));
V = (V1 + Told(i,j+1));
T(i,j) = (1-w)*Told(i,j) + w*(Tprev(i,j)*term_1 + H*term_2 + V*term_3);
end
end
counter = counter+1;
relerr = max(max(abs(Told-T)));
Told = T;
end
end
aa(iter) = iter;
%time_to_reach_steady_state
tt = dt*iter;
if rem(nx,2) == 0;
Temp(iter) = (T(nx/2,ny/2) + T(nx/2+1,ny/2+1))/2;
T_prev = (Tprev(nx/2,ny/2) + Tprev(nx/2+1,ny/2+1))/2;
error(iter) = max(abs(Temp(iter)-T_prev));
else
Temp(iter) = T(nx/2,ny/2);
error(iter) = max(abs(Temp(iter)-Tprev(nx/2,ny/2)));
end
l = sprintf(\'Transient 2D Heat Conduction Temperature contour\');
p = sprintf(\'nx=ny = %.0f,alpha = %0.2d, dt = %0.3f, Total time = %1.0f,\\n CFLx = %0.3d, CFLy = %0.3d\',nx, alpha,dt,t_total,k1,k2);
if iterative_solver == 1;
q = sprintf(\'Implicit Jacobi solver, Iteration number = %d\',iter);
q1 = sprintf(\'Implicit Jacobi solver \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\', aa(iter),tt, error(iter));
end
if iterative_solver == 2;
q = sprintf(\'Implicit Gauss-Siedel solver, Iteration number = %d\',iter);
q1 = sprintf(\'Implicit Gauss-Siedel solver \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\', aa(iter),tt, error(iter));
end
if iterative_solver == 3;
if solution_solver == 1;
q = sprintf(\'Implicit SOR with Jacobi,w = %0.2f, Iteration number = %d\',w,iter);
q1 = sprintf(\'Implicit SOR with Jacobi,w = %0.2f \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\',w , aa(iter),tt, error(iter));
end
if solution_solver == 2;
q = sprintf(\'Implicit SOR with Gauss-Siedel,w = %0.2f, Iteration number = %d\',w,iter);
q1 = sprintf(\'Implicit SOR with Gauss-Siedel,w = %0.2f \\n Iteration to reach steady state = %d, time = %0.2f sec, error = %.2d\',w , aa(iter),tt, error(iter));
end
end
figure(1);
plot(aa,Temp,\'.\');
xlabel(\'Iteration No.\');
ylabel(\'Mid-point Tempearature\');
title_text = {l,q1};
title(title_text);
Tprev = T;
figure(2);
contourf(x,y,T);
[a,b] = contourf(x,y,T);
clabel(a,b);
colorbar;
colormap(jet);
set(gca,\'Ydir\',\'reverse\');
xlabel(\'x\');
ylabel(\'y\');
title_text = {l,p,q};
title(title_text);
M(iter) = getframe(gcf); % Capturing frame to make movie
pause(0.05);
iter = iter + 1;
end
iteration_time_loop = iter-1% iteration under time loop
iter_space_loop = counter% iteration under space loop
error = error(iteration_time_loop)
%time_to_reach_steady_state
time_to_reach_steady_state = tt
Simulation_time = toc
% Movie Creation
if iterative_solver == 1;
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Implicit Jacobi solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour, Implicit Jacobi solver.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_Implicit_Jacobi_solver.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_Implicit_Jacobi_solver.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
if iterative_solver == 2;
% Movie Creation
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Implicit Gauss-Siedel solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour, Implicit Gauss-Siedel solver.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_Implicit_Gauss-Siedel_solver.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_Implicit_Gauss-Siedel_solver.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
if iterative_solver == 3;
if solution_solver == 1;
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Implicit SOR with Jacobi solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour, Implicit SOR with Jacobi solver.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_Implicit_SOR_with_Jacobi.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_Implicit_SOR_with_Jacobi.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
else
saveas(figure(1),\'Transient 2D Heat Conduction Mid-point Temperature, Implicit SOR with Gauss-Siedel solver.jpg\');
saveas(figure(2),\'Transient 2D Heat Conduction Temperature contour Implicit, SOR with Gauss-Siedel solver.jpg\');
movie(M);
save Transient_2D_Heat_Conduction_Temperature_contour_Implicit_SOR_with_Gauss-Siedel.avi;
videofile = VideoWriter(\'Transient_2D_Heat_Conduction_Temperature_contour_Implicit_SOR_with_Gauss-Siedel.avi\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
end
end
end