Skip to content

Instantly share code, notes, and snippets.

# KevinKParsons/steady-state-conduction-finite-differences.m Secret

Last active Jul 11, 2020
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 % 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
to join this conversation on GitHub. Already have an account? Sign in to comment