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