Assignment 2 (Date – 20/01/2025)
Submitted in partial fullfilment of the course
ME-864
COMPUTATIONAL FLUID DYNAMICS
in
THERMAL ENGINEERING
DEPARTMENT OF MECHANICAL ENGINEERING (M.TECH)
NATIONAL INSTITUTE OF TECHNOLOGY, KARNATAKA
SURATHKAL, MANGALORE-575025
Submitted By – Submitted To –
Ratnesh Kumar Sharma Dr.Ranjith M.
Roll No.-242TH013
Resistration No. - 2420373
Matlab code
1:solve fin equation by TDMA method and validate the results by analytical
solution.
% Parameters
tic;
L =1.0;
n =9;
dx =0.1;
% Boundary conditions
T_left =300;
T_right =400;
% Coefficients for the discretized system
A = -1;
B = 2.0005;
C = -1;
D = 0.1;
% Initialize arrays for the tridiagonal system
a = A * ones(1, n);
b = B * ones(1, n);
c = C * ones(1, n);
d = D * ones(1, n);
% Apply boundary conditions to the RHS
d(1) = d(1) - A * T_left;
d(end) = d(end) - C * T_right;
% TDMA (Thomas Algorithm)
for i = 2:n
r = a(i) / b(i-1);
b(i) = b(i) - r * c(i-1);
d(i) = d(i) - r * d(i-1);
end
T_interior = zeros(1, n);
T_interior(end) = d(end) / b(end);
for i = n-1:-1:1
T_interior(i) = (d(i) - c(i) * T_interior(i+1)) / b(i);
end
% Include boundary temperatures
T = [T_left, T_interior, T_right];
% Print the temperature values at each grid point
fprintf('Temperature Distribution:\n');
fprintf('Node\tPosition (m)\tTemperature (k)\n');
x = linspace(0, L, n + 2);
for i = 1:(n + 2)
fprintf('%d\t%.2f\t\t%.2f\n', i, x(i), T(i));
end
% Define the equation
temperature_profile = @(x) 200+(266.223 * exp(0.2236 * x)) + (-166.223 * exp(-
0.2236 * x));
% Create an array of x values from 0 to 1
x_values = 0:dx:1;
% Calculate temperatures for each x value
temperatures = arrayfun(temperature_profile, x_values);
% Create a table
result_table = table(x_values', temperatures', 'VariableNames', {'x', 'T'});
% Display the table
disp('Table for x, T, and:');
disp(result_table);
% Save the table to a CSV file
writetable(result_table, 'result_table.csv');
disp('Table saved to result_table.csv');
% Plot the results
plot(x, T, '-o', 'LineWidth', 2);
hold on
plot(x_values, temperatures, 'b', 'LineWidth', 2); % Plot T in blue
xlabel('Distance (x)');
ylabel('Temperature (T)');
legend('TDMA','Analytical temperatures');
title('Temperature Profile');
grid on;
xlabel('Position (m)');
ylabel('Temperature (°C)');
title('fin equation (TDMA)');
toc;
grid on;
Table -1:
Node Position Temperatures Temperatures
(m) (TDMA) (Analytical)
1 0.00 300.00 300.00
2 0.10 309.69 309.70
3 0.20 319.44 319.45
4 0.30 329.25 329.26
5 0.40 339.12 339.13
6 0.50 349.07 349.07
7 0.60 359.08 359.09
8 0.70 369.18 369.19
9 0.80 379.36 379.37
10 0.90 389.63 389.63
11 1.00 400.00 400.00
Elapsed time is 0.099084 seconds.
Jacobi Method
% Parameters
tic
L = 1.0;
n = 9;
dx = 0.1;
% Boundary conditions
T_left = 300;
T_right = 400;
% Coefficients
A = -1;
B = 2.0005;
C = -1;
D = 0.1;
% Initialize coefficients and arrays
a = A * ones(1, n);
b = B * ones(1, n);
c = C * ones(1, n);
d = D * ones(1, n);
% Apply boundary conditions to d
d(1) = d(1) - A * T_left;
d(end) = d(end) - C * T_right;
% Initialize temperature array
T_interior = zeros(1, n);
T_new = T_interior;
% Convergence criteria
Error_limit = 1e-5;
error = inf;
iteration = 0;
% iteration table
fprintf('Iteration\tMax Error (Absolute)\tMax Error (Percentage)\n');
% Jacobi iteration
while error > Error_limit
for i = 1:n
if i == 1
T_new(i)=(d(i)-c(i)*T_interior(i+1))/b(i);
elseif i == n
T_new(i)=(d(i)-a(i)*T_interior(i-1))/b(i);
else
T_new(i)=(d(i)-a(i)*T_interior(i-1)-c(i)*T_interior(i+1))/b(i);
end
end
abs_error = max(abs(T_new - T_interior));
percent_error = max(abs((T_new - T_interior) ./ T_new)) * 100;
% Display the iteration results
fprintf('%d\t\t%.5f\t\t\t%.5f%%\n',iteration,abs_error,percent_error);
error = abs_error;
T_interior = T_new;
iteration = iteration + 1;
end
% Include boundary temperatures
T = [T_left, T_interior, T_right];
% Display final results
fprintf('\nConverged after %d iterations.\n', iteration);
fprintf('Temperature Distribution:\n');
fprintf('Node\tPosition (m)\tTemperature (K)\n');
x = linspace(0, L, n + 2);
for i = 1:(n + 2)
fprintf('%d\t%.2f\t\t%.2f\n', i,x(i),T(i));
end
toc
gauss-seidel method:
% Parameters
tic:
L = 1.0;
n = 9;
dx = 0.1;
% Boundary conditions
T_left = 300;
T_right = 400;
% Coefficients
A = -1;
B = 2.0005;
C = -1;
D = 0.1;
% Initialize coefficients and arrays
a = A * ones(1, n);
b = B * ones(1, n);
c = C * ones(1, n);
d = D * ones(1, n);
d(1) = d(1) - A * T_left;
d(end) = d(end) - C * T_right;
% Initialize temperature array
T_interior = linspace(T_left, T_right, n);
% Convergence criteria
Error_limit = 1e-5;
error = inf;
iteration = 0;
% Display header for iteration table
fprintf('Iteration\tMax Error (Absolute)\tMax Error (Percentage)\n');
% Gauss-Seidel iteration
while error > Error_limit
T_old = T_interior;
for i = 1:n
if i == 1
T_interior(i) = (d(i) - c(i) * T_interior(i + 1)) / b(i);
elseif i == n
T_interior(i) = (d(i) - a(i) * T_interior(i - 1)) / b(i);
else
T_interior(i)=(d(i)-a(i)*T_interior(i-1)-c(i)*T_interior(i+1))/b(i);
end
end
abs_error = max(abs(T_interior - T_old));
error = max(abs((T_interior - T_old) ./ T_interior)) * 100;
% Display the iteration results
fprintf('%d\t\t%.6f\t\t\t%.6f%%\n', iteration,abs_error,percent_error);
error = abs_error;
iteration = iteration + 1;
end
% Include boundary temperatures
T = [T_left, T_interior, T_right];
% Display final results
fprintf('\nConverged after %d iterations.\n', iteration);
fprintf('Temperature Distribution:\n');
fprintf('Node\tPosition (m)\tTemperature (K)\n');
x = linspace(0, L, n + 2); % Grid positions
for i = 1:(n + 2)
fprintf('%d\t%.2f\t\t%.4f\n', i, x(i), T(i));
end
toc;
Successive-over relaxation (SOR) method.
% Parameters
tic;
L = 1.0;
n = 9;
dx = 0.1;
% Boundary conditions
T_left = 300;
T_right = 400;
% Coefficients
A = -1;
B = 2.0005;
C = -1;
D = 0.1;
% Initialize coefficients and arrays
a = A * ones(1, n);
b = B * ones(1, n);
c = C * ones(1, n);
d = D * ones(1, n);
d(1) = d(1) - A * T_left;
d(end) = d(end) - C * T_right;
% Initialize temperature array
T_interior = linspace(T_left, T_right, n);
% Convergence criteria
Error_limit = 1e-5;
error = inf;
iteration = 0;
% Relaxation factor
w_opt= 1.53;
% iteration table
fprintf('Iteration\tMax Error (Absolute)\tMax Error (Percentage)\n');
% SOR iteration
while error > Error_limit
T_old = T_interior;
for i = 1:n
if i == 1
T_new=(d(i)-c(i)*T_interior(i + 1))/b(i);
elseif i == n
T_new=(d(i)-a(i) * T_interior(i-1))/b(i);
else
T_new=(d(i)-a(i)*T_interior(i-1)-c(i)*T_interior(i+1))/b(i);
end
T_interior(i)=T_interior(i)+w_opt*(T_new - T_interior(i));
end
abs_error = max(abs(T_interior - T_old));
percent_error = max(abs((T_interior - T_old) ./ T_interior)) * 100;
% Display the iteration results
fprintf('%d\t\t%.6f\t\t\t%.6f%%\n',iteration,abs_error,percent_error);
error = abs_error;
iteration = iteration + 1;
end
% Include boundary temperatures
T = [T_left, T_interior, T_right];
% Display final results
fprintf('\nConverged after %d iterations.\n', iteration);
fprintf('Temperature Distribution:\n');
fprintf('Node\tPosition (m)\tTemperature (K)\n');
x = linspace(0, L, n + 2);
for i = 1:(n + 2)
fprintf('%d\t%.2f\t\t%.5f\n', i, x(i), T(i));
end
toc;
Table-2 : Temperatures calculated by all three Methods:
Node Position Temperatures Temperatures Temperatures
(K) (K) (K)
(m)
(Jacobi) (Gauss-seidel) (SOR)
1 0.00 300.00 300.0000 300.0000
2 0.10 309.69 309.6940 309.6940
3 0.20 319.44 319.4429 319.4428
4 0.30 329.25 329.2515 329.2514
5 0.40 339.12 339.1247 339.1246
6 0.50 349.07 349.0674 349.0674
7 0.60 359.08 359.0847 359.0847
8 0.70 369.18 369.1816 369.1815
9 0.80 379.36 379.3630 379.3629
10 0.90 389.63 389.6340 389.6340
11 1.00 400.00 400.0000 400.0000
Table-3: Relaxation factor Vs Number of iterations:
Relaxation Number of Relaxation Number of Iterations
factor (w) Iterations factor(w)
0.1 1231 1.2 75
0.2 666 1.3 61
0.3 454 1.4 48
0.4 369 1.5 34
0.5 269 1.53 (w_opt) 28
0.6 218 1.55 29
0.7 180 1.6 30
0.8 151 1.7 40
0.9 127 1.8 67
1.0 107 1.9 132
1.1 90 1.95 279
Relaxation Factor Vs Number of Iterations
1400
1200
1000
Number of Iterations
800
600
400
200
0
0 0.5 1 1.5 2 2.5
Relaxation Factor
Table-4:
Iterative Number of Elapsed
Methods Iterations time(seconds)
Jacobi 304 0.007145
Gauss-seidel 107 0.005849
SOR 28 0.005307
Inferences:
The study involved solving the fin equation using four numerical methods:
TDMA, Jacobi, Gauss-Seidel, and Successive Over-Relaxation (SOR). Clearly,
SOR is a very efficient method if the optimum value of the relaxation factor is
known.The results indicate that all methods produced accurate temperature
profiles, as confirmed by their alignment with the analytical solution. Among the
iterative methods, SOR demonstrated the highest efficiency, converging in just
28 iterations with an optimal relaxation factor (w_opt) of 1.53, compared to 107
iterations for Gauss-Seidel and 304 for Jacobi. The computational time also
reflected this efficiency, with SOR requiring the least time (0.0053 seconds)
compared to Jacobi (0.0071 seconds) and Gauss-Seidel (0.0058 seconds).
Additionally, the relaxation factor directly influenced the convergence speed,
with higher values of w generally leading to faster convergence up to an optimal
point. These findings highlight the practical advantages of using optimized
iterative methods like SOR for solving heat transfer problems efficiently.
