Golf Ball launcher code

% Apps Lab

% Lab 1

% Sam Georgal

clc

clear all

close all

Ball_D = .04; % diameter (m)

Ball_Du = 0.0001; % ball diamter uncertainty

Ball_M = 0.00275; % Mass (kg)

Ball_Mu = 0.00005; % Ball mass uncertainty

Pipe_L = 1.524; % Pipe Length in (m)

Pipe_Lu = 0.0254; % Pipe length uncertainty

Pipe_ID = 0.040386; % Pipe Inner Diameter (m)

Patm = 101325; % Pressure outisde of the pipe (Pa)

Patmu = 4; % atmospheric Pressure uncertainty

Air_D = 1.225; % Air density (kg/m^3)

dt = 0.00001; % Time step in seconds

Tape_s = 130000; % Tape breaking threshold (Pa)

Tape_su = 25000; % Tape threshold uncertainty

i = 1; % initilize variable

Air_k = 1.4; % air k value

Pmax = 10; % max pressure calculator

K = 0.8; % Loss coefficient

f = 0.02; % Friction factor

for k = 1:Pmax

  

   for j = 1:9

       Pd(1) = k*1000; % Downstream pressure (Pa)

       Pu(1) = Patm; % set initial pressure upstream (Pa)

       Ball_v(1) = 0; % Set initial ball velocity to 0 (m/s)

       Ball_x(1) = 0.02; % Set initial ball position to the ball radius (m)

       t(1) = 0; % Create time vector

       i = 1;

      

       % sets uncertainty values

       switch j

           case 2 % Mass upper bound

               Ball_M = Ball_M + Ball_Mu;

           case 3 % Mass lower bound

               Ball_M = Ball_M - 2*Ball_Mu;

           case 4 % Pressure lower bound

               Ball_M = Ball_M + Ball_Mu;

               Patm = Patm - Patmu;

           case 5 % Pressure upper bound

               Patm = Patm + 2*Patmu;

           case 6 % Tape upper bound

               Tape_s = Tape_s + Tape_su;

               Patm = Patm - Patmu;

           case 7 % Tape lower bound

               Tape_s = Tape_s - 2 * Tape_su;

           case 8 % Pipe upper bound

               Pipe_L = Pipe_L + Pipe_Lu;

               Tape_s = Tape_s + Tape_su;

           case 9 % Pipe Lower bound

               Pipe_L = Pipe_L - 2*Pipe_Lu;

       end

  

       Ball_A = pi*((Ball_D/2)^2); % Ball cross sectional area

       Ball_V = (4/3)*pi*((Ball_D/2)^3); % Ball volume

     

       while Ball_x(i) < Pipe_L

          

           Air_Vu = Ball_x(i)*pi*((Pipe_ID/2)^2) - (Ball_V/2); % Air upstream volume

           Air_Mu = Air_Vu*Air_D; % Air mass upstream

      

  

           % Calculate pressure downstream 

           Pd(i) = Pd(1)*((Pipe_L-Ball_x(1))/( ...

               Pipe_L-Ball_x(i)))^(Air_k); 

          

           % Calculate pressure upstream 

           Pu(i)= Patm - (0.5*Air_D*(Ball_v(i)^2))-(K*0.5*Air_D*(Ball_v(i)^2))-(f*0.5*(Ball_x(i)/Pipe_ID)*Air_D*(Ball_v(i)^2));

          

           % Test if tape has ruptured and downstream pressure is released

           if ((Pd(i)-Patm) >= Tape_s) || (Ball_x(i) >= Pipe_L-0.02)

               Pd(i) = Patm;

           end

  

  

           Fp(i) = Pu(i) * Ball_A; % Force from air pressure (N)

           Fv(i) = Pd(i) * Ball_A; % Force from air in pipe (N)

           Ft(i) = Fp(i)-Fv(i); % Sum of forces in x direction (N)

           Accel(i) = Ft(i)/(Ball_M+Air_Mu); % Ball acceleration in x direction (m/s^2) 

           dv = Accel(i)*dt; % Change in velocity (m/s)

           dY = Ball_v(i)*dt; % Change in position (m) 

          

           Ball_v(i+1) = Ball_v(i) + dv; % Find velocity of next instance

           Ball_x(i+1) = Ball_x(i) + dY; % Find positon at next instance

         

           t(i+1) = t(i) + dt;

           i = i + 1;

      

       end

       switch j

           case 1 % control

               Ball_vn(k) = Ball_v(end);

           case 2 % mass upper

               Ball_vmu = Ball_v(end) - Ball_vn(k);

           case 3 % Mass lower

               Ball_vml = Ball_v(end) - Ball_vn(k);

           case 4 % Pressure lower

               Ball_vPl = Ball_v(end) - Ball_vn(k);

           case 5 % Pressure upper

               Ball_vPu = Ball_v(end) - Ball_vn(k);

           case 6 % Tape upper

               Ball_vtu = Ball_v(end) - Ball_vn(k);

           case 7 % Tape lower

               Ball_vtl = Ball_v(end) - Ball_vn(k);

           case 8 % Pipe upper

               Ball_vpu = Ball_v(end) - Ball_vn(k);

           case 9 % Pipe lower

               Ball_vpl = Ball_v(end) - Ball_vn(k);

       end

  

   Pd = [];

   Pu = [];

   Ball_v = [];

   Ball_x = [];

   t = [];

   i = [];

   Accel = [];

   Fp = [];

   Ft = [];

   Fv = [];

   Vaxis(k) = k;

   end

   Ball_vu(k) = Ball_vn(k) + sqrt(Ball_vmu^2+ ...

       Ball_vPl^2+Ball_vtu^2 + Ball_vpl^2); % Ball velocity upper bound

   Ball_vl(k) = Ball_vn(k) - sqrt(Ball_vml^2+ ...

       Ball_vPu^2+Ball_vtl^2 + Ball_vpu^2); % Ball velocity lower bound

end

plot(Vaxis,Ball_vn,'b')

hold on

plot(Vaxis,Ball_vu,'r')

plot(Vaxis,Ball_vl,'r')

title('Predicted Exit velocities')

xlabel('Pressure (kPa)')

ylabel('Exit Velocity (m/s)')

axis tight

grid on

legend('Velocity','Error Bounds')