Calculate temperature of 2d plate over time using finite difference method and matrix

I have an assignment where I need to calculate the temperature distrubtion over time of a 2D aluminium plate until thermal equalibrium is reached.

temperature calculations are based on the explicit finite difference method, calculations will be done considering the region of interest covered by a uniform grid. The plate will be split up into a grid, user inputs their grid spacing requirment.

matrices will be used to show the temperature changes

initial matrix will be plate at start of calculations all filled with 0's except for boundary conditions.

equation for obtaining further matrices is attached in photo

i have included code so far which initialises first matrix and then sets boundary conditions.

user is asked to enter their grid spacing and time step and checks they are ok.

I am struggling with how to enter code to create second matrix and so on until thermal equalibrium is reached

any help please

clc
clear
close all                                                   %% Close all figure windows
%--------------------------------------------------------------------------------------------------------------------------%
% Define values for code
x=6;                                                        %% X Variable equal to
y=6;                                                        %% Y Variable equal to
k= 143;                                                     %% Conductivity of Aluminium
p= 2.8*10^3;                                                %% Density value
c= 795;                                                     %% Specific heat value
alpha= k/(p*c);                                             %% Calculate Alpha using above values
%--------------------------------------------------------------------------------------------------------------------------%
% User enters their grid spacing 
h=input('Enter your material grid spaceing here - ');       %% User enters their grid spaceing
while mod(6, h) ~= 0 ||  h<0 || h>=6                        %% Check is users grid spacing is valid by checking divisible into 6
    fprintf('Your grid spacing is not valid \n');           %% If grid spaceing entered is not valid tell user
    h=input('Enter your material grid spaceing here - ');
end
%--------------------------------------------------------------------------------------------------------------------------%
% Calculate Matrix dimensions
rows= (y/h)+1;                                                   %% Matrix rows is x(sheet width) / user input for h       
columns= (x/h)+1;                                                %% Matrix columns is y(Sheet Length) / user input h 
%--------------------------------------------------------------------------------------------------------------------------%
%% Initialise Matrix
T=zeros(rows,columns);                                      %% M is Matrix, filled with zeros until conditions are determined 
%--------------------------------------------------------------------------------------------------------------------------%
%% Ammend Matrix with boundary conditions
for c=1:columns
    d=(c-1)*h;                                              %% Variable 'd' is equal to c-1*user grid spacing
    if d<=2/3*x                                             %% If variable 'd' is less than 2/3 of x
        T(1,c)=12.5*d
    else
        T(1,c)= -25*d+150
    end
end
%--------------------------------------------------------------------------------------------------------------------------%
% User enters their time step value, this gets checked to ensure its valid
dt=input('Enter your time step here - ');                   %% User enters their time step here
Fo= (alpha*dt)/(h^2)                                        %% Calculation for **********
while Fo<0 || Fo>0.25
    fprintf('Your time step is not valid \n');              %% If time step entered is not valid tell user
    dt=input('Enter your time step here - ');
end
%--------------------------------------------------------------------------------------------------------------------------%
%Creating second Matrix

Answer 1

[Shubh]

Hi Christopher,

To create the subsequent matrices for temperature distribution, you'll need to implement the finite difference method as per the equation you've provided. The equation describes how to update the temperature at each grid point Tp​ based on the temperatures of the neighboring points and the current temperature, using a factor Fo​.

Here's a complete MATLAB code that performs the temperature distribution calculations using the explicit finite difference method:

clc
clear
close all % Close all figure windows

% Define values for code
x = 6; % X Variable equal to
y = 6; % Y Variable equal to
k = 143; % Conductivity of Aluminium
p = 2.8e3; % Density value
c = 795; % Specific heat value
alpha = k / (p * c); % Calculate Alpha using above values

% User enters their grid spacing
h = input('Enter your material grid spacing here - ');
while mod(6, h) ~= 0 || h < 0 || h >= 6
    fprintf('Your grid spacing is not valid \n');
    h = input('Enter your material grid spacing here - ');
end

% Calculate Matrix dimensions
rows = (y / h) + 1;
columns = (x / h) + 1;

% Initialise Matrix
T = zeros(rows, columns); % T is the temperature matrix, initially filled with zeros

% Amend Matrix with boundary conditions
for c = 1:columns
    d = (c - 1) * h;
    if d <= 2 / 3 * x
        T(1, c) = 12.5 * d;
    else
        T(1, c) = -25 * d + 150;
    end
end

% User enters their time step value, this gets checked to ensure its valid
dt = input('Enter your time step here - ');
Fo = (alpha * dt) / (h ^ 2);
while Fo < 0 || Fo > 0.25
    fprintf('Your time step is not valid \n');
    dt = input('Enter your time step here - ');
end

% Calculate further matrices until thermal equilibrium is reached
equilibrium_reached = false;
while ~equilibrium_reached
    T_new = T; % Create a new matrix for the updated temperatures
    for i = 2:rows-1
        for j = 2:columns-1
            T_new(i, j) = Fo * (T(i+1, j) + T(i-1, j) + T(i, j+1) + T(i, j-1) ...
                               - (4 - 1/Fo) * T(i, j));
        end
    end
    % Update the temperature matrix with the new values
    T = T_new;
    
    % Check for thermal equilibrium (no significant change in temperatures)
    if max(max(abs(T - T_new))) < 1e-4
        equilibrium_reached = true;
    end
    
    % You can add a condition to stop the simulation after a certain number of iterations
    % or a maximum time to prevent infinite loops in case equilibrium is never reached.
end

% Plot the final temperature distribution
imagesc(T);
colorbar;
title('Temperature Distribution');
xlabel('X Axis');
ylabel('Y Axis');

In this script, the "while" loop calculates the temperature at each grid point based on the current and neighboring grid points until thermal equilibrium is reached, which is when the maximum change in temperature between two successive iterations is less than a small threshold, here chosen as 1e-4. Adjust the threshold according to the precision required for your assignment. The "imagesc" function is used to visualize the final temperature distribution. Please ensure to validate the stopping criterion to avoid infinite loops.

Answer 2

[Saarthak Gupta]

Hi Christopher,

I understand you wish to calculate the dynamic temperature distribution of a body until equilibrium.

Using the explicit finite difference method, you will need to iteratively update the temperature values in the matrix based on the finite difference approximation of the heat equation, as given in the formula. The matrix represents the temperature at each grid point, and with each time step, you will create a new matrix representing the updated temperatures.

Refer to the following code:

clc
clear
close all                                                   %% Close all figure windows
%--------------------------------------------------------------------------------------------------------------------------%
% Define values for code
x=6;                                                        %% X Variable equal to
y=6;                                                        %% Y Variable equal to
k= 143;                                                     %% Conductivity of Aluminium
p= 2.8*10^3;                                                %% Density value
c= 795;                                                     %% Specific heat value
alpha= k/(p*c);                                             %% Calculate Alpha using above values
%--------------------------------------------------------------------------------------------------------------------------%
% User enters their grid spacing 
h=input('Enter your material grid spaceing here - ');       %% User enters their grid spaceing
while mod(6, h) ~= 0 ||  h<0 || h>=6                        %% Check is users grid spacing is valid by checking divisible into 6
    fprintf('Your grid spacing is not valid \n');           %% If grid spaceing entered is not valid tell user
    h=input('Enter your material grid spaceing here - ');
end
%--------------------------------------------------------------------------------------------------------------------------%
% Calculate Matrix dimensions
rows= (y/h)+1;                                                   %% Matrix rows is x(sheet width) / user input for h       
columns= (x/h)+1;                                                %% Matrix columns is y(Sheet Length) / user input h 
%--------------------------------------------------------------------------------------------------------------------------%
%% Initialise Matrix
T=zeros(rows,columns);                                      %% M is Matrix, filled with zeros until conditions are determined 
%--------------------------------------------------------------------------------------------------------------------------%
%% Ammend Matrix with boundary conditions
for c=1:columns
    d=(c-1)*h;                                              %% Variable 'd' is equal to c-1*user grid spacing
    if d<=2/3*x                                             %% If variable 'd' is less than 2/3 of x
        T(1,c)=12.5*d
    else
        T(1,c)= -25*d+150
    end
end
%--------------------------------------------------------------------------------------------------------------------------%
% User enters their time step value, this gets checked to ensure its valid
dt=input('Enter your time step here - ');                   %% User enters their time step here
Fo= (alpha*dt)/(h^2)                                        %% Calculation for **********
while Fo<0 || Fo>0.25
    fprintf('Your time step is not valid \n');              %% If time step entered is not valid tell user
    dt=input('Enter your time step here - ');
end
%--------------------------------------------------------------------------------------------------------------------------%

% Maximum change allowed for equilibrium
max_change_threshold = 1e-5; % Change this value as needed

% Initialize time
time = 0;

% Iterate until thermal equilibrium is reached
while true
    T_new = T;
    
    % Calculate new temperatures for interior points
    for i = 2:(rows-1)
        for j = 2:(columns-1)
            % A simplified version of the attached formula
            T_new(i,j) = T(i,j) + Fo*(T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1) - 4*T(i,j));
        end
    end
    
    % Check for thermal equilibrium
    max_change = max(abs(T_new - T),[],"all");
    if max_change < max_change_threshold
        % Exit loop if thermal equilibrium is reached
        break; 
    end
    
    % Update the temperature matrix and time for the next iteration
    T = T_new;
    time = time + dt;
end

disp(T);
fprintf('Thermal equilibrium reached at time = %f seconds\n', time);

Please note that this code assumes that the boundary conditions do not change over time. If they do, you will need to reapply them in each iteration. Also, the ‘max_change_threshold’ is a small value you choose to determine when the system is close enough to equilibrium (i.e., when changes between iterations are negligible).

Refer to the following resources for further reference:

https://www.visualslope.com/Library/FDM-for-heat-transfering.pdf