MATLAB: Numerical approximation of projectile motion with air resistance

air resistancenumerical approximationprojectile motion

For a school project, I need to estimate the maximum distance of a projectile without neglecting air resistance. I'm following the procedure outlined here:

http://wps.aw.com/wps/media/objects/877/898586/topics/topic01.pdf

On the second page it shows a nice, step by step process to find a numerical approximation. I am trying to reproduce the trajectory of the baseball that is shown on the last page in order to verify my model. However, the plot shows that the baseball will travel over 100 meters, while my model shows that it will travel about 80. Here is my code:

clear 
clc
close all
%Constants for Baseball
  m=.145 %kg
  A=pi*(.0366*2)^2/4 %m^2
  C=.5 %Drag Coefficient of a sphere
  rho= 1.2 %kg/m^3 (density of air)
  D=rho*C*A/2
  g=9.81 %m/s^2 (acceleration due to gravity)
    %Initial Conditions
    delta_t= .001 %s
    x(1)=0
    y(1)=0
    v=50 %m/s
    theta=35 %deg
    vx=v*cosd(theta)
    vy=v*sind(theta)
    t(1)=0 
    %Start Loop
    i=1
    while min(y)> -.001
        ax=-(D/m)*v*vx;
        ay=-g-(D/m)*v*vy;
        vx=vx+ax*delta_t;
        vy=vy+ay*delta_t;
        x(i+1)=x(i)+vx*delta_t+.5*ax*delta_t^2;
        y(i+1)=y(i)+vy*delta_t+.5*ay*delta_t^2;
        t(i+1)=t(i)+delta_t;
        i=i+1;
    end
    plot(x,y)
    xlabel('x distance (m)')
    ylabel('y distance (m)')
    title('Projectile Path')

If anyone has a few minutes to take a look at this, I'd really appreciate it. I can't seem to find the problem. Thank you!

Best Answer

You are forgetting to update v inside your loop.

v = sqrt(vx^2 + vy^2);

Right now it is sitting at v = 50 the whole time.

Related Solutions

MATLAB: Sample Code of Projectile Motion. Unknown constant that’s used.

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.

MATLAB: Problem in plotting in a nested while loop in a for loop

Replace the code after count = count+1; with the following to get separate figures (though the curves ae all the same!):

    figure
    plot(x,y)
    xlabel('x distance (m)')
    ylabel('y distance (m)')
    title(['Projectile Path    Vx = ' num2str(Vx(aa))] )
    
    
    %hold on

I think your while loop logic needs modifying to the following in order to get different curves:

        flag = true;
    while flag
        v = sqrt(vx^2 + vy^2);
        ax=-(D/m)*vx^2;
        ay=-g-(D/m)*vy^2;
        vx=vx+ax*delta_t;
        vy=vy+ay*delta_t;
        x(nn+1)=x(nn)+vx*delta_t+.5*ax*delta_t^2;
        y(nn+1)=y(nn)+vy*delta_t+.5*ay*delta_t^2;
        t(nn+1)=t(nn)+delta_t;
        nn=nn+1;
        if y(nn)<=0
            x(nn) = NaN; y(nn) = NaN; t(nn) = NaN;
            flag = false;
        end
        
    end

If you want all the curves to appear on the same figure then keep your original plot commands.

Also, I suspect your ay term should be:

     ay=-g-(D/m)*vy*abs(vy);

as the drag will oppose gravity when the projectile is coming down (i.e. when vy is negative).

Related Question