Governing Physics and Discretization

The problem focuses on steady-state energy transport combining conduction mechanisms with convective boundary interactions. The differential form describing thermal equilibrium within the control volume is:

$$\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right) + \dot{q} = 0$$

Convective boundaries follow a Robin condition where surface flux balances environmental exchange:

$$-k \frac{\partial T}{\partial n} = h(T_{surface} - T_{ambient})$$

The computational mesh employs a structured Cartesian grid with $N_x$ and $N_y$ divisions, defining cell widths $\Delta x$ and $\Delta y$.

Algorithmic Implemantation

Domain Configuration and Initialization

% Physical constants
span_x = 0.1; span_y = 0.1;
res_x = 50; res_y = 50;
conductivity = 200;     % W/m*K
film_coeff = 25;        % W/m^2*K
temp_surroundings = 25; % Celsius

% Spatial step sizes
dx = span_x / (res_x - 1);
dy = span_y / (res_y - 1);

[x_grid, y_grid] = meshgrid(linspace(0, span_x, res_x), linspace(0, span_y, res_y));

Global Stiffness Matrix Construction

Linearizing the discretized equations results in a sparse system $A\mathbf{x} = \mathbf{b}$.

node_count = res_x * res_y;
SystemMatrix = spalloc(node_count, node_count, 9 * node_count);
LoadVector = zeros(node_count, 1);

% Traverse internal control volumes
for i = 2:(res_x-1)
    for j = 2:(res_y-1)
        center_idx = sub2ind([res_x, res_y], i, j);
        
        % Diffusion coefficients
        c_w = conductivity / dy^2;
        c_e = conductivity / dy^2;
        c_s = conductivity / dx^2;
        c_n = conductivity / dx^2;
        
        A_diag = -(c_w + c_e + c_s + c_n);
        
        SystemMatrix(center_idx, center_idx) = A_diag;
        SystemMatrix(center_idx, sub2ind([res_x, res_y], i+1, j)) = c_e;
        SystemMatrix(center_idx, sub2ind([res_x, res_y], i-1, j)) = c_w;
        SystemMatrix(center_idx, sub2ind([res_x, res_y], i, j+1)) = c_n;
        SystemMatrix(center_idx, sub2ind([res_x, res_y], i, j-1)) = c_s;
        
        LoadVector(center_idx) = 0; % Source term placeholder
    end
end

Boundary Condition Enforcement

Specific rows are overwritten to enforce Dirichlet or Neumann-type constraints dierctly into the matrix.

% West Wall: Fixed Temperature
for j = 1:res_y
    idx = sub2ind([res_x, res_y], 1, j);
    SystemMatrix(idx, :) = 0;
    SystemMatrix(idx, idx) = 1;
    LoadVector(idx) = 100; % Prescribed temperature
end

% Top Wall: Convective Exchange
% Adjusting diagonal for heat transfer coefficient
for i = 1:res_x
    idx = sub2ind([res_x, res_y], i, res_y);
    SystemMatrix(idx, idx) = 1 + (film_coeff * dx / conductivity);
    LoadVector(idx) = (film_coeff * dx / conductivity) * temp_surroundings;
    
    % Zero out adjacent neighbors to decouple boundary from interior logic here
    SystemMatrix(idx, [sub2ind([res_x, res_y], i-1, res_y), ... 
                        sub2ind([res_x, res_y], i+1, res_y)]) = 0;
end

Solver and Convergence Loop

An iterative relaxation scheme resolves the linear system until residuals fall below tolerance.

% Initial state
TempField = ones(res_x, res_y) * 25;

max_steps = 10000;
convergence_tol = 1e-6;

for iter = 1:max_steps
    OldState = TempField(:);
    
    % Update core nodes using Gauss-Seidel updates
    for i = 2:(res_x-1)
        for j = 2:(res_y-1)
            lin_id = sub2ind([res_x, res_y], i, j);
            neighbor_sum = (TempField(i+1,j) + TempField(i-1,j))*conductivity/dy^2 + ...
                           (TempField(i,j+1) + TempField(i,j-1))*conductivity/dx^2;
            divisors = 2*conductivity*(1/dy^2 + 1/dx^2);
            TempField(i,j) = neighbor_sum / divisors;
        end
    end
    
    % Check residual magnitude
    error_norm = max(abs(TempField(:) - OldState));
    if error_norm < convergence_tol
        break;
    end
end

Output Generation

Visualizing the scalar field validates the numerical convergence and physical behavior.

figure('Color','w');
contourf(x_grid, y_grid, TempField, 20, 'LineStyle','none');
colormap(parula);
xlabel('X Coordinate (m)'); ylabel('Y Coordinate (m)');
title('Steady State Thermal Distribution');
colorbar; axis equal tight;

Post-processing computes temperature gradients to derive local heat flux vectors along both axes for energy balance verification.