transient-conduction-finite-differences.m

% 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

Hi Kevin How are you
I appreciate your sharing the code...I just want to let you know I modified your code for the implicit solution, which as no stability condition for time step, notice I tested the solution for dT = 1500 ...just replace the following code line in the calculation

% Calculates temperatures for new time step
T(i,j) = (Fo*(T(i-1,j) + T(i+1,j) + T(i,j-1)...
                + T(i,j+1)) +Told(i,j))/((1+4*Fo));