In Projects

MATLAB FEA 2D Steady State Heat Transfer

This program is a thermal Finite Element Analysis (FEA) solver for steady state heat transfer across 2D plates. The program numerically solves the steady state conduction problem using the Finite Difference Method. After the results are calculated, the program displays a color contour plot of the temperature of the plate after a given time interval.

Results are shown for a 2D fin with specified boundry conditions.

The model assumes 2D conduction with constant properties throughout the simulation. The model takes advantage of the vertical line of thermal symmetry at x = 0.020 m by adding an adiabatic boundary condition along the line of symmetry. The model assumes a specified convection boundary condition on x = 0 and y = Ly. The model assumes constant temperature boundary condition along the fin base.

MATLAB Code
        % Program:
        % Steady_Conduction_With_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