% Program: | |
% transient-conduction-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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment