ENGINEERING

Home Engineering
  • Engineering
  • Senior Project
  • Projects
  • 3 Doors Game Simulation
  • Woodworking: Bedstand
  • Woodworking: Desk
  • Numbers Game
  • Path Planning
  • Snake Video Game AI
  • Minimum Spanning Tree
  • Heat Map
  • MATLAB: Rocket Simulation
  • MATLAB: Bike Suspension
  • MATLAB: Thermal FEA SS
  • MATLAB: Thermal FEA Tr
  • Notes
  • Statistics
  • Steoscopic Estimates
  • Cubed Roots
  • Other
  • Free Fall Distance
  • Speed of Sound
  • Travel
  • Solar System Scale
  • Videos
  • Yellowstone Geysers
  • Geysers and Waterfalls
  • Yellowstone and Glacier
  • Grinnell Glacier
  • Yellowstone Geysers
  • Photography

    Subjects

    MATLAB FEA 2D Steady-State Heat Transfer



    This program is a thermal Finite Element Analysis (FEA) solver for steady-state heat transfer involving 2D plates. The program numerically solves the steady-state two dimensional conduction problem using the finite difference method and creates the color contour plot.


    The model assumes steady-state 2D conduction with constant properties. The model takes advantage of vertical line of thermal symmetry at x = 0.020 m by adding adiabatic BC on line of symmetry. Assumes specified convection BC on x = 0 and y = Ly. Assumes constant temperature BC along fin base.



    MATLAB Source 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