Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
% Program:
% steady-state-conduction-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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment