steady-state-conduction-finite-differences.m

% Program:
% steady-state-conduction-finite-differences.m
% Steady-state 2D conduction solver using finite difference method.
%
% Description:
% Numerically solves the steady-state two dimensional conduction problem
% using the finite difference method and plots color contour plot. Assumes 
% steady-state 2D conduction with constant properties. Takes advantage of 
% vertical line of thermal symmetry at x = 0.020 m by adding adiabatic BC 
% on line of symmetry. Implements specified convection BC on x = 0 and
% y = Ly. Implements constant temperature BC along fin base.
% 
% Variable List:
% T = Temperature (deg. Celsius)
% T1 = Boundary condition temperature 1 (deg. Celsius)
% T2 = Boundary condition temperature 2 (deg. Celsius)
% Tinf = Ambient fluid temperature (deg. Celsius)
% theta = Non-dimensionalized temperature difference = (T-T1)/(T2-T1)
% Lx = Plate length in x-direction (m)
% Ly = Plate length in y-direction (m)
% AR = Aspect ratio of Ly / Lx to ensure dx = dy
% h = Convection coefficient (W/m^2K)
% k = Thermal conductivity (W/mK)
% Bi = Finite-difference Biot number
% Nx = Number of increments in x-direction
% Ny = Number of increments in y-direction
% dx = Increment size in x-direction (m)
% dy = Increment size in y-direction (m)
% dT = Temperature step between contours
% tol = Maximum temperature difference for convergence (deg. Celsius)
% pmax = Maximum number of iterations
% Told = Stores temperature matrix for previous time step
% diff = Maximum difference in T between iterations (deg. Celsius)
% i = Current column
% j = Current row
% p = Current iteration
% v = Sets temperature levels for contours
% x =  Create x-distance node locations
% y =  Create y-distance node locations
% Nc = Number of contours for plot

clear, clc      % Clear command window and workspace

Lx = .020;      % Plate half-length in x-direction (m)
Ly = .200;      % Plate length in y-direction (m)
Nx = 14;        % Number of increments in x-direction

AR = Ly/Lx;     % Aspect ratio of Ly /Lx to ensure dx = dy
Ny = AR*Nx;     % Number of increments in y-direction
dx = Lx/Nx;     % Increment size in x-direction (m)
dy = Ly/Ny;     % Increment size in y-direction (m)

T1 = 200;       % BC temperature at end of fin (deg. Celsius)
T2 = 100;       % BC temperature at base of fin (deg. Celsius)
Tinf = T2;      % Ambient fluid temperature (deg. Celsius)
h = 500;        % Convection coefficient (W/m^2K)
k = 50;         % Thermal conductivity (W/m^2K)
Bi = h*dx/k;    % Finite-difference Biot number

T = T1*ones(Nx+1,Ny+1);     % Initialize T matrix to T1 everywhere
T(1:Nx+1,1) = T1;           % Initialize base of fin to T1 BC

tol = 10^-6;                % Max temp delta to converge (deg. Celsius)
pmax = 10*10^6;             % Maximum number of iterations

x = 0:dx:Lx;                % Create x-distance node locations
y = 0:dy:Ly;                % Create y-distance node locations

for p = 1:pmax              % Loop through iterations
    Told = T;               % Store previous T array as Told for later
    
    for j = 2:Ny            % Loop through rows
        for i = 2:Nx        % Loop through interior columns
            
            % Calculates convection BC along left side
            if i == 2
                T(1,j) = (2*T(2,j)+T(1,j-1)+T(1,j+1)+2*Bi*Tinf)/(4+2*Bi); 
            end
            
            % Calculates interior node temperatures
            T(i,j) = .25*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)); 

            
            % Calculates adiabatic BC along right side
            if i == Nx
                T(Nx+1,j) = .25*(2*T(Nx,j)+T(Nx+1,j-1)+T(Nx+1,j+1)); 
            end
            
        end
    end

    for i = 2:Nx	% Loop through interior columns at top row
        
        % Calculates top left corner, conv/conv BC's            
        if i == 2
            T(1,Ny+1) = (T(1,Ny)+T(2,Ny+1)+2*Bi*Tinf)/(2+2*Bi); 
        end
        
        % Calculates convection BC along top
        T(i,Ny+1) = (2*T(i,Ny)+T(i-1,Ny+1)+T(i+1,Ny+1)+2*Bi*Tinf)...
            /(4+2*Bi);
        
        % Calculates top right, conv/adiabatic BC's
        if i == Nx
            T(Nx+1,Ny+1) = (T(Nx+1,Ny)+T(Nx,Ny+1)+Bi*Tinf)/(2+Bi); 
        end
       
    end
    
    diff = max(max(abs(T - Told)));     % Max difference between iterations
    fprintf('Iter = %8.0f - Dif. = %10.6f deg. C\n', p, diff);
    if (diff < tol)     % Exit iteration loop because of convergence
        break
    end
end

fprintf('Number of iterations = \t %8.0f \n\n', p)  % Print time steps

if (p == pmax)          % Warn if number of iterations exceeds maximum
    disp('Warning: code did not converge')
    fprintf('\n')                               
end
 
disp('Temperatures in brick in deg. C = ') 
for j = Ny+1:-1:1             	% Loop through each row in T
    fprintf('%7.1f', T(:,j))	% Print T for current row
    fprintf('\n')
end

Nc = 50;                        % Number of contours for plot
dT = (T2 - T1)/Nc;              % Temperature step between contours
v = T1:dT:T2;                   % Sets temperature levels for contours
colormap(jet)                   % Sets colors used for contour plot
contourf(x*100, y*100, T',v, 'LineStyle', 'none') 
colorbar                        % Adds a scale to the plot
axis equal tight                % Makes the axes have equal length
title('Contour Plot of Temperature in deg. C')  
xlabel('x (cm)') 
ylabel('y (cm)')
pause(0.001)                    % Pause between time steps to display graph