Sample Code of Projectile Motion. Unknown constant that's used.

Disclamer: I am using the following code to help guide me through a projectile motion problem.

function projectile(theta,v0,dt,tf)
t = 0:dt:tf;
vx = v0*cos(theta);
vy = v0*sin(theta);
vx(1) = vx;
vy(1) = vy;
v = sqrt((vx^2)+(vy^2));
v(1) = v;
x(1) = 0;
y(1) = 0;
B2divm = 0.00004;
g = 9.81;
for i = 1:length(t)-1
    v = sqrt((vx(i)^2)+(vy(i)^2));
    x(i+1) = x(i) + vx(i)*dt;
    vx(i+1) = vx(i) - (B2divm*v*vx(i))*dt;
    y(i+1) = y(i) + vy(i)*dt;
    vy(i+1) = vy(i) - g*dt - (B2divm*v*vy(i))*dt;
end
plot(x,y,'r')
end

I follow and reasonably understand the code, except for the use of the constant 'B2divm=0.0004'.

I strongly suspect that 'B2divm' above has something to do with Drag and air density, but I just don't know.

Anybody know what it means? How it's derived?

Thanks!!

Answer 1

[Jim Riggs]

You are probably right about that.

As you can see, the position states (X & Y) are updated by adding a delta position which is computed as V*dt, i.e.

x(i+1) = x(i) - vx(i)*dt;
y(i+1) = y(i) - vy(i)*dt;

You would expect a similar expression for the velocity update.

vx(i+1) = vx(i) + ax(i)*dt
vy(i+1) = vy(i) + ay(i)*dt

The delta velocity is the acceleration times the time step,

deltaV = a*dt. So, a reaonable expression for the acceleratio of a projectile is the sum of (weight + drag force) divided by the mass. Note that the constant name is B2divm - maybe "B2 divided by mass". This suggests that "B2" is the sum of the forces acting on the projectile. I see that there is a g*dt in the Y axis, therefore "B2" represents only the aerodynamic drag force.

Drag Force = 0.5*(air density)*Velocity^2*(Reference Aera)*(Drag Coefficient)

Note that the Drag Force is computed using the total velocity (a scalar value). We want to determine the amount of the force in the X and Y directions, so the X component is (Drag Force) * cos(gamma) and

the Y component is (Drag Force)*sin(gamma) where gamma is the flight path angle (i.e. direction of the velocity vector).

We substitude cos(gamma) = Vx/V and sin(gamma) = Vy/V, now we have:

Fx = (Drag Force) * Vx/V

Fy = (Drag Force) * Vy/V

or

Fx = 0.5*(air density) * V^2 * (Reference Area) * (Drag Coefficient) * Vx / V

= 0.5*(air density) * (Reference area) * (Drag Coefficient) * V * Vx

Now the acceleration in the X direction due to air resistance = Fx / m, so

Ax(from drag) = {0.5 * (air density) * (Reference Area) * (drag Coefficient) / m } * V * Vx

For short trajectories, everything inside the curly brackets can be considered as constant, and this would be the "B2 divided by m" term, B2divm.

If the trajectory has any significant variation in altitude, the air density term must be computed inside the for loop for each pass.

[threex5x7]

Wow!

Thank you very much for taking the time to explain this to me!!

Very clear and meaningful!!

[Jim Riggs]

My Pleasure!