In Projects

MATLAB FEA 2D Transient Heat Transfer

This program is a thermal Finite Element Analysis (FEA) solver for transient heat transfer across 2D plates. The program numerically solves the transient conduction problem using the Finite Difference Method. The program displays a color contour plot of the temperature of the plate for each time step.

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

The model assumes transient 2D conduction with constant properties throughout the simulation. A stable time step must be set by taking into account the physical time constant. Selecting unstable time steps may lead to inconsistent results.

  • #
    Figure 1
    Transient Conduction

    Contour plot of temperature (deg C) for 2D, transient, constant property conduction with Lx = Ly = 1m, Nx = Ny = 20, and stable time step of 1000 seconds after a time of 184.44 hours.

  • #
    Figure 2
    Transient Conduction, Unstable Start

    Contour plot of temperature (deg C) for 2D, transient, constant property conduction with Lx = Ly = 1m, Nx = Ny = 20, and unstable time step of 1500 seconds after a time of 27.08 hours.

  • #
    Figure 3
    Transient Conduction, Unstable Progression

    Contour plot of temperature (deg C) for 2D, transient, constant property conduction with Lx = Ly = 1m, Nx = Ny = 20, and unstable time step of 1500 seconds after a time of 33.33 hours.

  • #
    Download
    MATLAB Code

    Download the MATLAB program.

MATLAB Code
        % Program:
        % Transient_Conduction_With_Finite_Differences.m
        % Transient 2D Conduction Solver using Finite Difference Method.
        %
        % Description:
        % Numerically solves the transient two dimensional conduction problem
        % using the finite difference method and plots color contour plot. Assumes 
        % transient 2D conduction with constant properties.
        % 
        % Variable List:
        % T = Temperature (deg. Celsius)
        % T1 = Boundary condition temperature 1 (deg. Celsius)
        % T2 = Boundary condition temperature 2 (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)
        % x = Create x-distance node locations
        % y = Create y-distance node locations
        % 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
        % Lmax = Maximum number of time steps before stopping
        % Told = Stores temperature array for previous time step
        % diff = Percent difference at x = y = 0.4 m compared to theoretical
        % k = Thermal conductivity (W/mK)
        % rho = Density (kg/m^3)
        % Cp = Specific heat (J/kgK)
        % alpha = Thermal diffusivity (m^2/s)
        % deltaT = Time step (s)
        % deltaT_stable_max = Maximum stable time step (s)
        % Fo = Fourier number 
        % tol = Stop when T at x = y = 0.4 m reaches within this percent of tol_ans
        % tol_ans = Steady-state temperature at x = y = 0.4m (deg. Celsius)
        % Tplot = Stores T (deg. Celsius) at x = y = 0.4 m at each time step
        % L = Loop counter
        % p = Current iteration
        % i = Current column
        % j = Current row
        % v = Sets temperature levels for contours
        % Nc = Number of contours for plot

        clear, clc              % Clear command window and workspace

        Lx = 1;                 % Plate length in x-direction (m)
        Ly = 1;                 % Plate length in y-direction (m)
        Nx = 20;                % Number of increments in x-direction

        Ny = 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 = 0;                 % BC temperature 1 (deg. Celsius)
        T2 = 100;               % BC temperature 2 (deg. Celsius)

        k = 0.72;               % Thermal conductivity (W/mK)
        rho = 1920;             % Density (kg/m^3)
        Cp = 835;               % Specific heat (J/kgK)
        alpha = k/(rho*Cp);     % Thermal diffusivity (m^2/s)

        deltaT = 1300;                          % Time step (seconds)
        deltaT_stable_max = .25*dx^2/alpha;     % Maximum stable time step
        Fo = alpha*deltaT/dx^2;                 % Fourier number 

        tol = 0.1;              % Stop when T at x = y = 0.4m reaches percentage of tol_ans
        tol_ans = 16.8568;      % Steady-state temperature at x = y = 0.4m
        Lmax = 10^4;            % Maximum number of time steps before stopping

        T = T1*ones(Nx+1,Ny+1); % Initialize T array to T1 everywhere
        T(2:Nx,Ny+1) = T2;      % Initialize top row to T2 boundary condition
        T(1,Ny+1)    = T1;      % Initialize top left
        T(Nx+1,Ny+1) = T1;      % Initialize top right
        Tplot = ones(Lmax,1);   % Initialize Tplot to allocate memory

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

        for L = 1:Lmax          % Loop through time steps
            Told = T;           % Store previous T array as Told for next time step
            
            for j = 2:Ny        % Loop through rows
                for i = 2:Nx    % Loop through columns
                    
                    % Calculates temperatures for new time step
                    T(i,j) = Fo*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1)...
                        + Told(i,j+1)) + (1-4*Fo)*Told(i,j);
                end 
            end    
            
            Tplot(L) = T(Nx/2-1,Ny/2-1); % Store T at x = y = 0.4m
            
            % Percent difference at x = y = 0.4 m compared to theoretical
            diff = abs((T(Nx/2-1,Ny/2-1) - tol_ans)/tol_ans*100); ...
                
            fprintf('Time step = %8.0f - Diff. = %10.6f percent\n', L, diff);
            if (diff < tol)     % Exit iteration loop because reached steady-state
                break 
            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, y, 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 at time = ',...
                num2str(deltaT*L/3600),' h']) 
            xlabel('x (m)') 
            ylabel('y (m)')        
            set(gca,'XTick',0:.1:Lx)    % Sets the x-axis tick mark locations
            set(gca,'YTick',0:.1:Ly)    % Sets the y-axis tick mark locations
            pause(0.001)                % Pause between time steps to display graph
                
            %if L == 55 || L == 65 || L == 80 % Chosen time steps to save plot
            %    saveas(gcf, ['Transient_Plot_Unstable_',num2str(L)], 'jpg'); % save plot
            %end    
        end

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

        if (L == Lmax)      % Warn if number of iterations exceeds maximum
            disp('Warning: Maximum time steps exceeded')   
            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