I can't believe that the infinite series really is an analytical solution for your PDE that has such a complicated boundary condition.
Can you provide a link for the underlying PDE of your infinite series solution ?
Comparing Finite difference method answer with analytical for 2d plate heat transfer
5 views (last 30 days)
Show older comments
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.
The matrices created using FDM then needs to be compared with a matrices created using an analytical solution (formula added)
I have the code for the FDM matrix and graph but struggling to correctly enter the code for the analytical solution to be able to compare the 2.
Any help is much appreciated
clc %% CLear Command Window
clear %% Clear Variables from code
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
h=0.1;
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); %% T 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 or equal to 2/3 of x complete following formula
T(1,c)=12.5*d;
else %% If variable 'd' is greater than or equal to 2/3 of x complete following formula
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
dt=0.1*h^2/alpha;
Fo= (alpha*dt)/(h^2); %% Calculation for Fouriers Number
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 Further Matrices
max_change_threshold = 0.01; %% Maximum change allowed for equilibrium, Change this value as needed
time = 0; %% Initialize time
while true %% Iterate until thermal equilibrium is reached
T_new = T;
for i = 2:(rows-1) %% Calculate new temperatures for interior points
for j = 2:(columns-1)
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
fprintf('Thermal equilibrium reached at time = %f seconds\n', time); % Display Thermal equilibrium message has been reached to user
fprintf('Temperature from Numerical Solution = \n%T')
disp(T)
contourf(0:h:x,0:h:y,T.')
colorbar
%--------------------------------------------------------------------------------------------------------------------------%
% Analytical Solution
AT=zeros(rows,columns); %% AT 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 or equal to 2/3 of x complete following formula
AT(1,c)=12.5*d;
else %% If variable 'd' is greater than or equal to 2/3 of x complete following formula
AT(1,c)= -25*d+150;
end
end
while true %% Iterate until thermal equilibrium is reached
AT_new = AT;
for i = 2:(rows-1) %% Calculate new temperatures for interior points
for j = 2:(columns-1)
for n = 1:200
AT_new(i,j) = AT_new+(sin(2 * n * pi / 3) / (n^2 * sinh(n * pi))) * sin(n * pi * x / 6) * sin(n * pi * y / 6);
end
end
end
end
fprintf('Temperature from Analytical Solution = \n%AT')
disp(AT)
%--------------------------------------------------------------------------------------------------------------------------%
% Plot the final temperature distribution
figure(1);
contourf(flip(T));
colormap("jet");
colorbar;
title('Temperature Distribution');
xlabel('X Axis');
ylabel('Y Axis');
axis equal tight
5 Comments
Torsten
on 16 Feb 2024
Edited: Torsten
on 16 Feb 2024
X = 0:0.1:6;
nx = numel(X);
Y = 0:0.1:6;
ny = numel(Y);
mmax = 100;
tol = 1e-6;
T = zeros(nx,ny);
for ix = 1:nx
x = X(ix);
for iy = 1:ny
y = Y(iy);
error = Inf;
im = 1;
while error > tol && im <= mmax
summand = sin(2*pi*im/3)/(im^2*sinh(im*pi))*sin(im*pi*x/6)*sinh(im*pi*y/6);
error = abs(summand);
T(ix,iy) = T(ix,iy) + summand;
im = im+1;
end
T(ix,iy) = 450/pi^2*T(ix,iy);
end
end
contourf(X,Y,T.')
colorbar
Answers (1)
William Rose
on 16 Feb 2024
Edited: William Rose
on 16 Feb 2024
[Edit: add a line of code to display the maximum temperature for each value of N.]
The boundary conditions are T=0 on three sides (x=0, x=6, y=6). On the y=0 edge, T=0 to 50 deg (linear ramp) from x=0 to x=4, then T ramps down from 50 to 0 deg, from x=4 to x=6. Therefore we expect to see T=0 along 3 edges, and at all four corners. We expect to see t=50 at (x,y)=(4,0).
Your plot has "x" running vertically, which is fine of course - just so you know.
Analytical:
The script below computes and plots the analytical temperature distribution using summations of N=5, 10, and 200 terms. You can see small differences between the plots, near the hottest point.
x=0:.06:6; y=0:.08:6;
% I like having different number of elements along x and y for grids,
% to protect against transposing without realizing.
N=[5,10,200]; % elements in summation
T=zeros(length(y),length(x),length(N));
for i=1:length(N)
T(:,:,i)=temp(N(i),x,y); % analytical solution for temperature
figure; contourf(x,y,T(:,:,i)) % plot temperature
axis equal; xlabel('X'); ylabel('Y')
title(['Analytical Temperature Solution, N=',num2str(N(i))])
clim([0,50]); colorbar % use same color range for all plots
end
fprintf('Maximum temperature for N=%d,%d,%d: %.1f,%.1f,%.1f.\n',...
N,max(T,[],[1 2]))
function T=temp(N,x,y)
T=zeros(length(y),length(x));
for n=1:N
T=T+(sin(2*pi*n/3)*sinh(n*pi*y'/6)*sin(n*pi*x/6))/(n^2*sinh(pi*n));
end
T=450*T/pi^2;
end
@Torsten's contour plot shows the non-zero temperature boundary along the edge y=6, which disagrees with the boundary conditions in the document you posted. I get the same result as Torsten. Therefore I think the analytical solution has an error.
Good luck.
6 Comments
Torsten
on 25 Feb 2024
Why do you use this while loop to compute AT ? Your analytical solution does not depend on t.
And what are x and y in the expression
AT_new(i,j) = AT(i,j) + (sin(2*pi*n/3)/(n^2*sinh(n*pi))*sin(n*pi*x/6)*sinh(n*pi*y/6));
They are both set to 6 and don't change in the for-loops although they had to.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)