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.

• 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.

• 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.

• 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.

##### 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
```