How do I fix my P_ups function?
Show older comments
clear
clc
global P_atm
global P_tube
global tube_ID
global tube_L
global ball_D
global x
global v
global a_dens
global TR
global P_rupt
ball_D = 0.04; %meters
ball_m = 0.00275; %kg
tube_OD = 0.05; %meters
tube_ID = 0.04; %meters
tube_L = 1.52; %meters
P_atm = 101.4; %kPa
dt = 0.00001; %time intervals
K = 0.78; %entrance coefficient
f = 0.018; %friction coefficient
a_dens = 1.225; %kg/m^3
A_ball = pi*(ball_D/2)^2;
P_rupt = 250; %kPa
j = 1;
for P_tube = 1:1:10
TR = 0;
i = 1;
time(i) = 0;
v(i) = 0;
x(i) = 0.02;
P_ds(i) = P_tube;
P_ups(i) = 0;
while x(i) <= tube_L-ball_D/2
P_ups(i) = P_upsCalc(i);
P_ds(i) = P_dsCalc(i);
F_ups(i) = A_ball * P_ups(i) * 1000;
F_ds(i) = A_ball * P_ds(i) * 1000;
a(i) = (F_ups(i) - F_ds(i)) / ball_m;
time(i+1) = time(i) + dt;
v(i+1) = v(i) + a(i)*dt;
x(i+1) = x(i) + v(i+1)*dt;
i = i + 1;
end
Velocity(j) = v(i);
Pressure(j) = P_tube;
j = j + 1;
end
plot(Pressure, Velocity)
xlabel("Pressure (kPa)")
ylabel("Velocity (m/sec)")
function P_ds = P_dsCalc(b)
global tube_ID
global tube_L
global ball_D
global x
global P_tube
global P_atm
global P_rupt
global TR
Volume_init = pi*(tube_ID/2)^2*(tube_L - ball_D/2);
Volume_inter = pi.*(tube_ID/2)^2*(tube_L -ball_D/2 - x(b));
Volume_inter = max(Volume_inter, 1e-6);
if TR == 1
P_ds = P_atm;
return
end
P_ds = P_tube * (Volume_init / Volume_inter)^1.4;
P_ds = min(P_ds, 300);
if P_ds >= P_rupt
TR = 1;
P_ds = P_atm;
end
end
function P_ups = P_upsCalc(d)
global P_atm a_dens K f tube_L tube_ID v
v_local = v(d); % guaranteed scalar
P_b = 0.5 * a_dens * (v_local^2);
P_f = (f * tube_L * a_dens * (v_local^2)) / (2 * tube_ID);
P_k = K * (0.5 * a_dens * (v_local^2));
P_ups = P_atm - P_b - P_f - P_k;
P_ups = max(P_ups, 0);
end
Accepted Answer
Walter Roberson
about 16 hours ago
f = 0.018; %friction coefficient
f is NOT declared global at that point in the code, so f is just a normal variable.
global P_atm a_dens K f tube_L tube_ID v
but at this point in the code, inside P_upsCalc, f is declared global . As global variable f has not been initialized, global variable f in P_upsCalc defaults to [] . That [] poisons the result of the calculation, so P_upsCalc returns [] but the calling context tries to assign that [] into an indexed variable.
3 Comments
Walter Roberson
on 4 Sep 2026 at 18:14
We see here one of the reasons that working with globals is so trouble-prone.
Hi @Anthony
Did you inherit legacy code that uses global variables? Modern AI models generally avoid using global variables unless instructed to modify existing code that already relies on them. Moreover, the coding style in this code is inconsistent.
Many of the issues occur in the while loop, which continues to increment the counter variable i by one as long as the condition x(i) <= tube_L-ball_D/2 remains satisfied. If the code becomes stuck in an infinite loop, it is advisable to include a break statement within an if statement so that the loop terminates only when a specific condition is met.
Are you attempting to solve an ordinary differential equation in which the pressure varies along the tube? if the pressure varies along the tube and can be represented by a differential equation, using ode45 solver would generally be preferable to relying on global variables as a workaround.
global P_atm
global P_tube
global tube_ID
global tube_L
global ball_D
global x
global v
global a_dens
global TR
global P_rupt
global f
% -----------------
global K % previously undeclared
% -----------------
ball_D = 0.04; % meters
ball_m = 0.00275; % kg
tube_OD = 0.05; % meters
tube_ID = 0.04; % meters
tube_L = 1.52; % meters
P_atm = 101.4; % kPa
dt = 0.00001; % time intervals
K = 0.78; % entrance coefficient
f = 0.018; % friction coefficient
a_dens = 1.225; % kg/m^3
A_ball = pi*(ball_D/2)^2;
P_rupt = 250; % kPa
j = 1;
for P_tube = 1:1:10
TR = 0;
i = 1;
time(i) = 0;
v(i) = 0;
x(i) = 0.02;
P_ds(i) = P_tube;
P_ups(i)= 0;
while x(i) <= tube_L-ball_D/2
P_ups(i) = P_upsCalc(i);
P_ds(i) = P_dsCalc(i);
F_ups(i) = A_ball * P_ups(i) * 1000;
F_ds(i) = A_ball * P_ds(i) * 1000;
a(i) = (F_ups(i) - F_ds(i)) / ball_m;
time(i+1) = time(i) + dt;
v(i+1) = v(i) + a(i)*dt;
x(i+1) = x(i) + v(i+1)*dt;
i = i + 1;
% -----------------
%% Set a condition for exiting an unintentional infinite while-loop
% -----------------
if i == 10
break; % Immediately exits the while loop
end
% -----------------
end
Velocity(j) = v(i);
Pressure(j) = P_tube;
j = j + 1;
end
plot(Pressure, Velocity)
xlabel("Pressure (kPa)")
ylabel("Velocity (m/sec)")
function P_ds = P_dsCalc(b)
global tube_ID
global tube_L
global ball_D
global x
global P_tube
global P_atm
global P_rupt
global TR
Volume_init = pi*(tube_ID/2)^2*(tube_L - ball_D/2);
Volume_inter = pi.*(tube_ID/2)^2*(tube_L -ball_D/2 - x(b));
Volume_inter = max(Volume_inter, 1e-6);
if TR == 1
P_ds = P_atm;
return
end
P_ds = P_tube * (Volume_init / Volume_inter)^1.4;
P_ds = min(P_ds, 300);
if P_ds >= P_rupt
TR = 1;
P_ds = P_atm;
end
end
function P_ups = P_upsCalc(d)
global P_atm a_dens K f tube_L tube_ID v
v_local = v(d); % guaranteed scalar
P_b = 0.5 * a_dens * (v_local^2);
P_f = (f * tube_L * a_dens * (v_local^2)) / (2 * tube_ID);
P_k = K * (0.5 * a_dens * (v_local^2));
P_ups = P_atm - P_b - P_f - P_k;
P_ups = max(P_ups, 0);
end
More Answers (0)
Categories
Find more on Assembly in Help Center and File Exchange
on 5 Sep 2026 at 9:51
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
