Untitled

function chimney_heat_transfer(T_gas, k, h_o)
    % Other Constants
    T_ambient = 15; % °C
    T_sky = 250; % K
    h_i = 75; % W/m²°C
    epsilon = 0.9;
    sigma = 5.67e-8; % Stefan-Boltzmann constant W/m²K⁴

    % Grid
    T = zeros(10, 10); % considering symmetry

    % Initial guess
    T(:,:) = T_gas;

    % Iteration settings
    max_iter = 5000;
    tol = 1e-6;

    % Jacobi iteration
    for iter = 1:max_iter
        T_old = T;
        
        % Update inner nodes
        for i = 2:9
            for j = 2:9
                T(i,j) = 0.25 * (T_old(i+1,j) + T_old(i-1,j) + T_old(i,j+1) + T_old(i,j-1));
            end
        end
        
        % Update boundary conditions
        T(1,:) = (T_old(2,:) + dx * h_i / k * T_gas) / (1 + dx * h_i / k);
        T(10,:) = (T_old(9,:) + dx * (h_o * T_ambient + epsilon * sigma * (T_sky^4 - T_old(10,:).^4) / k)) / (1 + dx * h_o / k);
        
        % Check convergence
        if max(max(abs(T - T_old))) < tol
            break
        end
    end

    % Calculate heat loss
    Q_conv = h_o * sum(T(10,:) - T_ambient) * 1; % 1 m length
    Q_rad = epsilon * sigma * sum(T(10,:).^4 - T_sky^4) * 1; % 1 m length

    % Total heat loss per unit length
    Q_total = Q_conv + Q_rad;

    % Display result
    disp(['Total heat loss per unit length: ', num2str(Q_total), ' W/m'])

    % Plot temperature distribution
    [X,Y] = meshgrid(1:10, 1:10);
    surf(X, Y, T);
    xlabel('X');
    ylabel('Y');
    zlabel('Temperature (°C)');
    title('Temperature Distribution in Chimney Wall');
end