How do I fix my P_ups function?

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

Unable to perform assignment because the left and right sides have a different number of elements.

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

Answer 1

[Walter Roberson]

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.

[Walter Roberson]

We see here one of the reasons that working with globals is so trouble-prone.

[Anthony]

I just declared 'f' global but it still gives me the same error.

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
global f
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

Unable to perform assignment because the left and right sides have a different number of elements.

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

[Sam Chak]

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