1D Heat Equation Explicit Scheme (fixed)

28 visualizaciones (últimos 30 días)
Gabriel Quattromani
Gabriel Quattromani el 19 de Abr. de 2022
Editada: Gabriel Quattromani el 30 de Abr. de 2022
clear all
p=2050; % density(kg/m^3)
k=0.52; % conductivity(W/m-K)
c=1840; % specific heat capacity(J/kg-K)
T_ini=20; % initial temperature(K)
T_inf=-15; % external temperature(K)
alpha=k/(p*c); % Let thermal diffusivity be alpha
t=60*60*60*24; %time in seconds
DELTA_t=100; % change in time per timestep
M=t/DELTA_t; % number of timesteps
N=300; % number of steps of depth
L=3; % maximum depth being analyzed
DELTA_x=L/(N); % change in depth per step
t_diff=t/M;
time=0:t_diff:t;
%explicit scheme
% Boundary Conditions
T(1:N,1:M+1)=T_ini;
T(1,1:M+1)=T_inf;
CFL=(alpha*DELTA_t)/(DELTA_x^2); % Courant number, which should be less than 0.5
% Finding temperature for next timestep for each depth
for j=1:M
for i=2:N-1
T(i,j+1)=T(i,j)+((alpha*DELTA_t)/((DELTA_x)^2))*(T(i+1,j)-2*T(i,j)+T(i-1,j));
end
end
Z=N*(0.1/L); % t=0.1 seconds in terms of the time step
xm=(69/300)*L;
plot(time,T(Z,:), 'b',time,T(69,:), 'r')
xlabel('Time')
ylabel('Temperature')
grid on
set(gca,'fontsize',16)
legend('0.1','xm')

Respuesta aceptada

Torsten
Torsten el 19 de Abr. de 2022
You use
for i=0:0.1:50
...
T(i,j) = ...
But only positive integer values are possible as array indices. Thus something like T(0,j), T(0.1,j) etc is not allowed.
Additionally, you need a boundary condition at x=xstart and x=xend. What are these conditions ?
Delta_x and Delta_t are scalar values. Thus addressing them in the if-statement
if (alpha*DELTA_t(j))/(DELTA_x(i)^2)<=0.5
as vector is wrong.

Más respuestas (1)

Alan Stevens
Alan Stevens el 19 de Abr. de 2022
With
for i=0:0.1:50
you are trying to index variables with zero. However, Matlab's indexing starts at 1.

Productos


Versión

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by