1D Heat Equation Explicit Scheme (fixed)

9 vues (au cours des 30 derniers jours)
Gabriel Quattromani
Gabriel Quattromani le 19 Avr 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')

Réponse acceptée

Torsten
Torsten le 19 Avr 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.

Plus de réponses (1)

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

Catégories

En savoir plus sur Heat and Mass Transfer dans Help Center et File Exchange

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by