Verifying my code for this problem, are these the right solution and graphs?

11 views (last 30 days)
clc; clear; close all;
% Parameters
L = 10; % mm
W = 4; % mm
dx = 1; % mm
dy = 1; % mm
T_inf = 22; % Ambient temperature (°C)
Tb = 50; % Base plate temperature (°C)
h = 630; % Convective heat transfer coefficient (W/m^2.K)
k = 20; % Thermal conductivity (W/m.K)
% Grid size
Nx = L/dx + 1;
Ny = W/dy + 1;
% Initialize temperature field
T = ones(Nx, Ny) * T_inf; % Start with ambient temperature
% Set bottom boundary condition
T(:,1) = Tb;
% Convection boundary condition coefficient
conv_coeff = h * dy / k;
% Iteration parameters
tol = 1e-4;
max_iter = 5000;
error = 1;
iter = 0;
% Jacobi Iteration
while error > tol && iter < max_iter
T_old = T;
for i = 2:Nx-1
for j = 2:Ny-1
T(i,j) = (T_old(i+1,j) + T_old(i-1,j) + T_old(i,j+1) + T_old(i,j-1)) / 4;
end
end
% Apply convective boundary conditions
T(1, :) = (conv_coeff * T_inf + T_old(2, :)) / (1 + conv_coeff); % Left boundary
T(Nx, :) = (conv_coeff * T_inf + T_old(Nx-1, :)) / (1 + conv_coeff); % Right boundary
T(:, Ny) = (conv_coeff * T_inf + T_old(:, Ny-1)) / (1 + conv_coeff); % Top boundary
% Compute error
error = max(max(abs(T - T_old)));
iter = iter + 1;
end
fprintf('Converged in %d iterations\n', iter);
Converged in 57 iterations
% (b) Plot temperature profile along vertical lines at x=1mm, x=4mm, x=6mm
x_positions = [2, 5, 7]; % Corresponding to x=1mm, x=4mm, x=6mm
figure;
hold on;
for x_index = x_positions
plot(T(x_index, :), 'LineWidth', 2);
end
xlabel('Vertical position (y-axis)');
ylabel('Temperature (°C)');
legend('x=1mm', 'x=4mm', 'x=6mm');
title('Temperature Profile along Vertical Lines');
grid on;
hold off;
% (c) Plot temperature contour
figure;
contourf(T', 20);
colorbar;
xlabel('x (mm)');
ylabel('y (mm)');
title('Temperature Contour of the Rib');
% (d) Refine the grid with dx=dy=0.5mm
dx = 0.5; dy = 0.5;
Nx = L/dx + 1;
Ny = W/dy + 1;
T_fine = ones(Nx, Ny) * T_inf;
T_fine(:,1) = Tb;
conv_coeff = h * dy / k;
error = 1;
iter = 0;
while error > tol && iter < max_iter
T_old = T_fine;
for i = 2:Nx-1
for j = 2:Ny-1
T_fine(i,j) = (T_old(i+1,j) + T_old(i-1,j) + T_old(i,j+1) + T_old(i,j-1)) / 4;
end
end
T_fine(1, :) = (conv_coeff * T_inf + T_old(2, :)) / (1 + conv_coeff);
T_fine(Nx, :) = (conv_coeff * T_inf + T_old(Nx-1, :)) / (1 + conv_coeff);
T_fine(:, Ny) = (conv_coeff * T_inf + T_old(:, Ny-1)) / (1 + conv_coeff);
error = max(max(abs(T_fine - T_old)));
iter = iter + 1;
end
fprintf('Refined grid converged in %d iterations\n', iter);
Refined grid converged in 204 iterations
% Compare temperature profile at x=4mm for both grids
x_index_coarse = 5; % Corresponding to x=4mm in the coarse grid
x_index_fine = 9; % Corresponding to x=4mm in the fine grid
figure;
plot(T(x_index_coarse, :), 'r--', 'LineWidth', 2);
hold on;
plot(T_fine(x_index_fine, :), 'b-', 'LineWidth', 2);
xlabel('Vertical position (y-axis)');
ylabel('Temperature (°C)');
legend('Coarse grid (1mm)', 'Fine grid (0.5mm)');
title('Comparison of Temperature Profiles at x=4mm');
grid on;
hold off;
  1 Comment
Saif
Saif on 18 Mar 2025
% I will go and check, thanks for noticing , someone else in my team wrote another code, any chance someone can take a look at that pls, thanks
clc;
clear;
close all;
% Define Parameters
k = 20; % Thermal conductivity (W/m.K)
h = 630; % Convective heat transfer coefficient (W/m^2.K)
T_inf = 295.15; % Ambient temperature (Kelvin)
T_base = 323.15; % Base temperature (Kelvin)
L = 10e-3; % Length of the rib (m)
W = 4e-3; % Width of the rib (m)
% Define grid spacing (Δx = Δy = 1mm)
dx = 1e-3;
dy = 1e-3;
% Compute number of grid points
nx = L/dx + 1;
ny = W/dy + 1;
% Initialize matrices
T = ones(nx, ny) * T_inf; % Initial guess
T(1, :) = T_base; % Bottom boundary at base temperature
% Jacobi Method Parameters
max_iter = 10000; % Maximum iterations
MinError = 1e-6; % Convergence criteria
error = inf;
iter = 0;
% Coefficient for convective boundary
beta = h * dx / k;
% Solve using Jacobi method
while error > MinError && iter < max_iter
T_old = T;
% Update interior points using Jacobi iteration
for i = 2:nx-1
for j = 2:ny-1
T(i, j) = (T_old(i+1, j) + T_old(i-1, j) + T_old(i, j+1) + T_old(i, j-1)) / 4;
end
end
% Apply convective boundary conditions (top, left, and right)
for j = 2:ny-1
T(nx, j) = (T_old(nx-1, j) + 2 * beta * T_inf) / (1 + 2 * beta); % Right
end
for i = 2:nx-1
T(i, ny) = (T_old(i, ny-1) + 2 * beta * T_inf) / (1 + 2 * beta); % Top
T(i, 1) = (T_old(i, 2) + 2 * beta * T_inf) / (1 + 2 * beta); % Left
end
% Apply convective boundary to top-right and top-left corners
T(nx, ny) = (T_old(nx-1, ny) + T_old(nx, ny-1) + 2 * beta * T_inf) / (2 + 2 * beta);
T(nx, 1) = (T_old(nx-1, 1) + T_old(nx, 2) + 2 * beta * T_inf) / (2 + 2 * beta);
% Compute error for convergence
error = max(max(abs(T - T_old)));
iter = iter + 1;
end
% Part B: Temperature Profile at Selected x Locations
x_values = [1, 4, 6] * 1e-3; % In meters
figure;
hold on;
for i = 1:length(x_values)
x_idx = round(x_values(i) / dx) + 1;
plot(linspace(0, W, ny), T(x_idx, :), 'LineWidth', 2);
end
legend('x = 1 mm', 'x = 4 mm', 'x = 6 mm');
xlabel('Distance along y (m)');
ylabel('Temperature (K)');
title('Temperature Profile at Selected x Locations');
grid on;
% Part C: Temperature Contour Plot
figure;
contourf(linspace(0, L, nx), linspace(0, W, ny), T', 20, 'LineWidth', 0.5);
colorbar;
xlabel('X (m)');
ylabel('Y (m)');
title('Temperature Contour Plot');
axis equal;
% Part D: Solve for Δx = Δy = 0.5 mm
dx2 = 0.5e-3;
dy2 = 0.5e-3;
% Compute number of grid points
nx2 = L/dx2 + 1;
ny2 = W/dy2 + 1;
% Initialize matrices
T2 = ones(nx2, ny2) * T_inf; % Initial guess
T2(1, :) = T_base; % Bottom boundary at base temperature
% Reset error and iteration count
error = inf;
iter = 0;
beta2 = h * dx2 / k;
% Solve using Jacobi method
while error > MinError && iter < max_iter
T_old = T2;
for i = 2:nx2-1
for j = 2:ny2-1
T2(i, j) = (T_old(i+1, j) + T_old(i-1, j) + T_old(i, j+1) + T_old(i, j-1)) / 4;
end
end
for j = 2:ny2-1
T2(nx2, j) = (T_old(nx2-1, j) + 2 * beta2 * T_inf) / (1 + 2 * beta2); % Right
end
for i = 2:nx2-1
T2(i, ny2) = (T_old(i, ny2-1) + 2 * beta2 * T_inf) / (1 + 2 * beta2); % Top
T2(i, 1) = (T_old(i, 2) + 2 * beta2 * T_inf) / (1 + 2 * beta2); % Left
end
T2(nx2, ny2) = (T_old(nx2-1, ny2) + T_old(nx2, ny2-1) + 2 * beta2 * T_inf) / (2 + 2 * beta2);
T2(nx2, 1) = (T_old(nx2-1, 1) + T_old(nx2, 2) + 2 * beta2 * T_inf) / (2 + 2 * beta2);
error = max(max(abs(T2 - T_old)));
iter = iter + 1;
end
% Compare Temperature Profile at x = 4mm for Δx=1mm and Δx=0.5mm
x_position = 4e-3;
x_idx1 = round(x_position / dx) + 1;
x_idx2 = round(x_position / dx2) + 1;
figure;
plot(linspace(0, W, ny), T(x_idx1, :), 'r-', 'LineWidth', 2);
hold on;
plot(linspace(0, W, ny2), T2(x_idx2, :), 'b--', 'LineWidth', 2);
legend('Δx=1mm', 'Δx=0.5mm');
xlabel('Distance along y (m)');
ylabel('Temperature (K)');
title('Comparison of Temperature Profile at x=4mm');
grid on;

Sign in to comment.

Answers (1)

Torsten
Torsten on 18 Mar 2025
Edited: Torsten on 18 Mar 2025
Look at the contour plot. The temperature at x = 0 should be at T_b = 50°C. This is not the case in your graphics. Maybe you interchanged the axes ?
And why do the temperature profiles in your last plot start at 1 and end at 9 ? And why does the red curve end at 5 ? Both curves should start at 0 and end at 4.
I suggest creating
x = linspace(0,L,Nx)
y = linspace(0,W,Ny)
in both cases (coarse and fine) so that you can make plots over the correct x- and y-grids.

Categories

Find more on Introduction to Installation and Licensing in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!