Index in position 2 is invalid. Array indices must be positive integers or logical values.
2 次查看(过去 30 天)
显示 更早的评论
%% Material parameters
data_rho = [7840,0.44];
data_Cp = [490,0.0733333];
data_k = [13.1947,0.0126919];
%% Geometry
L = 0.03;
%% Boundary conditions
T0 = 1000;
Tinfinity = 20 + 273;
htc_top = 50;
htc_sides = 100;
htc_btm = 20;
Emissivity = 0.0;
time = 0;
%% thermocouple positions
tpx = [0.015 0.027 0.027];
tpy = [0.015 0.015 0.027];
%% Simulation control
n = 3;
CFL = 0.5;
itPlot = 10;
%% Discretize space
dx = L/(n-1);
x = 0:dx:L;
y = L:-dx:0; % T(1,1) is the top left
%% Initialize variables
T = ones(n,n)*T0;
T_new = zeros(n,n);
dtMat = zeros(n,n);
it = 0;
%% Physical constants
stefan = 5.670367e-8;
%% Record predictions
Temp_thermocouples = interp2(x,y,T,tpx,tpy);
SaveNumber = 1;
TimeVector(SaveNumber) = time;
TempSavedMatrix(SaveNumber,:) = Temp_thermocouples;
%% Model calculation
while max(max(T)) > 300 + 273.15
%% Thermophysical properties
rho = data_rho(1) + data_rho(2)*T;
k0 = data_k(1) + data_k(2)*T;
Cp = data_Cp(1) + data_Cp(2)*T;
%% Calculate temperature dependent parameters
alpha = k0./rho./Cp;
Br = stefan*Emissivity*dx./k0;
Biot_top = htc_top*dx./k0;
Biot_bottom = htc_btm*dx./k0;
Biot_sides = htc_sides*dx./k0;
%% time step
% interior
for i = 2: n-1
for j = 2: n-1
dtMat(i,j) = CFL*dx*dx/4./alpha(i,j);
end
end
% top left
i = 1;
j = 1;
dtRadTerm = 0;
dtTerm = 2+ Biot_top(i,j) + Biot_sides(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
for j = 2: n-1
% top edge
i = 1;
dtRadTerm = 0;
dtTerm = 2+ Biot_top(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
end
% top right
i = 1;
j = n;
dtRadTerm = 0;
dtTerm = 2+ Biot_top(i,j) + Biot_sides(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
for i = 2: n-1
% left edge
j = 1;
dtRadTerm = 0;
dtTerm = 2+ Biot_sides(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
end
for i = 2: n-1
% right edge
j = n;
dtRadTerm = 0;
dtTerm = 2+ Biot_sides(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
end
% btm left
i = n;
j = 1;
dtRadTerm = 0;
dtTerm = 2+ Biot_bottom(i,j) + Biot_sides(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
% btm edge
for j = 2: n-1
i = n;
dtRadTerm = 0;
dtTerm = 2+ Biot_bottom(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
end
% btm right
i = n;
j = n;
dtRadTerm = 0;
dtTerm = 2+ Biot_bottom(i,j) + Biot_sides(i,j)+ dtRadTerm;
dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
dt = min(dtMat(:));
%% Fourier number
Fo = alpha*dt/(dx*dx);
%% Calculate the new T vector T^(p+1)
% top left
i = 1;
j = 1;
part1 = (1-4*Fo(i,j)*T(i,j));
part2 = 2*Fo(i,j)*T(i+1,j) + T(i,j+1);
part3 = 2*Fo(i,j)*(Biot_top(i,j)+Biot_sides(i,j))*(Tinfinity-T(i,j));
T_new(i,j) = part1+part2+part3;
for j = 2: n-1
% top surface
i = 1;
part1 = (1-4*Fo(i,j)*T(i,j));
part2a = 2*Fo(i,j)*T(i+1,j);
part2b = Fo(i,j)*T(i,j+1);
part2c = Fo(i,j)*T(i,j-1);
part3 = 2*Fo(i,j)*(Biot_top(i,j))*(Tinfinity-T(i,j));
T_new(i,j) = part1+part2a+part2b+part2c+part3;
end
% top right
i = 1;
j = n;
part1 = (1-4*Fo(i,j)*T(i,j));
part2 = 2*Fo(i,j)*T(i+1,j) + T(i,j-1);
part3 = 2*Fo(i,j)*(Biot_top(i,j)+Biot_sides(i,j))*(Tinfinity-T(i,j));
T_new(i,j) = part1+part2+part3;
for i = 2: n-1
% left surface
j = 1;
part1 = (1-4*Fo(i,j)*T(i,j));
part2a = 2*Fo(i,j)*T(i,j+1);
part2b = Fo(i,j)*T(i+1,j);
part2c = Fo(i,j)*T(i-1,j);
part3 = 2*Fo(i,j)*(Biot_sides(i,j))*(Tinfinity-T(i,j));
T_new(i,j) = part1+part2a+part2b+part2c+part3;
end
% right surface
j = n;
part1 = (1-4*Fo(i,j)*T(i,j));
part2a = 2*Fo(i,j)*T(i,j-1);
part2b = Fo(i,j)*T(i+1,j);
part2c = Fo(i,j)*T(i-1,j);
part3 = 2*Fo(i,j)*(Biot_sides(i,j))*(Tinfinity-T(i,j));
T_new(i,j) = part1+part2a+part2b+part2c+part3;
% bottom left
i = n;
j = 1;
part1 = (1-4*Fo(i,j)*T(i,j));
part2 = 2*Fo(i,j)*T(i-1,j) + T(i,j+1);
part3 = 2*Fo(i,j)*(Biot_bottom(i,j)+Biot_sides(i,j))*(Tinfinity-T(i,j));
T_new(i,j) = part1+part2+part3;
% bottom edge
%
i = n;
part1 = (1-4*Fo(i,j)*T(i,j));
part2a = 2*Fo(i,j)*T(i-1,j);
part2b = Fo(i,j)*T(i,j+1);
part2c = Fo(i,j)*T(i,j-1);
part3 = 2*Fo(i,j)*(Biot_bottom(i,j))*(Tinfinity-T(i,j));
T_new(i,j) = part1+part2a+part2b+part2c+part3;
issue is with
part2c =
Fo(i,j)*T(i,j-1);
2 个评论
Jan
2021-2-4
@John McCann: Please post the complete error message and only the relevcant part of the code. Then the readers do not have to guess, where the problem occurs.
回答(1 个)
Jan
2021-2-4
The debugger helps you to find the problem. Type in the command window:
dbstop if error
Then run your code again. Matlab will stop when the error occurs and you can examine the values of the locally used variables.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Clocks and Timers 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!