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

Question

0 Stimmen

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;