Projectile motion of a cannonball with varying drag

Hey everyone,

https://imgur.com/a/yUTmYXh

For a project, I need to simulate the projectile motion of a "Paris Cannon". I need to consider the air resistance changing according to the altitude with the following equation:

where y0 = 1000m.

And the drag is:

So far, my code looks like this:

clc
clear
m=94;
A=pi*(.21)^2;
C=.12;
p0= 1.225;
g=9.81;
%Initial Conditions
dt= .001;
x(1)=0;
y(1)=0;
v=10000; 
angle=20 ;
vx=v*cosd(angle);
vy=v*sind(angle);
t(1)=0 ;
i=1;
y0 = 1000;
    while min(y)> -.001
        p = p0*exp(-y/y0);
        Drag = p*C*A/2*v^2;
        ax=-(Drag/m)*(vx*vx+vy*vy)^0.5*vx;
        ay=-g-(Drag/m)*(vx*vx+vy*vy)^0.5*vy;
        vx=vx+ax*dt;
        vy=vy+ay*dt;
        x(i+1)=x(i)+vx*dt+.5*ax*dt^2;
        y(i+1)=y(i)+vy*dt+.5*ay*dt^2;
        t(i+1)=t(i)+dt;
        i=i+1;
    
    end
      plot(x,y,'-')
      xlabel('x distance (m)')
      ylabel('y distance (m)')
      pause(0.001);

I cannot get the correct simulation here. If I take drag as constant and place it in initial conditions, it works fine. But when I put it in while loop, the whole simulation gives me absolute nonsense. I would appriciate your help.

Thanks.

Answer 1

[Michael Madelaire]

Hi,

I am not sure about the equations you use. But here is my attempt.

Try running it, and if it is what you are looking for, maybe I can explain how I did it.

m=94;
A=pi*(.21)^2;
C=.12;
p0= 1.225;
g=9.81;
%Initial Conditions
dt= .001;
x=0;
y=0;
v=10000; 
angle=20 ;
vx=v*cosd(angle);
vy=v*sind(angle);
t=0 ;
y0 = 1000;
i=1;
% Half step
p = p0*exp(-y/y0);
D = C*(p*v^2)/2*A;
ax = -D/m;
ay = -g;
while y(i) >= 0
    i = i+1;
    vx(i) = vx(i-1)+ax(i-1)*dt;
    vy(i) = vy(i-1)+ay(i-1)*dt;
    
    x(i) = x(i-1)+(vx(i-1)+vx(i))/2*dt;
    y(i) = y(i-1)+(vy(i-1)+vy(i))/2*dt;
    
    p = p0*exp(-y(i)/y0);
    v = sqrt(vx(i)^2+vy(i)^2);
    D = C*(p*v^2)/2*A;
    ax(i) = -D/m;
    ay(i) = -g;
end

figure;
plot(x,y,'-')
xlabel('x distance (m)')
ylabel('y distance (m)')

[Ozan Terzioglu]

Hello and thank you,

In your case, you are completely eliminating the air resistance since you took ax -D/m and ay =-g. When drag is changing

ax and ay should change with the following;

        ax=-(Drag/m)*(vx*vx+vy*vy)^0.5*vx;
        ay=-g-(Drag/m)*(vx*vx+vy*vy)^0.5*vy;
       

But i cannot implement it because drag changes as y changes and vice versa.

It find it hard to explain myself right now. I have uploaded the question here so that you can check it out.

https://imgur.com/a/yUTmYXh

[Michael Madelaire]

That cannot be the correct equations the units do not add up.

(m/s^2)*sqrt(m^2/s^2)*(m/s) = (m/s^2)*(m/s)*(m/s) = (m/s^2)*(m^2/s^2)

= m^3/s^4

And the drag D is a force parallel to x counteracting the movement. Given newtons 2'nd law F=m*a => a=F/m=D/m

Then a is the deacceleration due to air drag, assuming the eqaution for D is correct.

%% Edit

I have included the difference, here ax=0

[Michael Madelaire]

I realized that you have to consider the drag in both components, duh.

m=94;
A=pi*(.21)^2;
C=.12;
p0= 1.225;
g=9.81;
%Initial Conditions
dt= .001;
x=0;
y=0;
v=10000; 
angle=20 ;
vx=v*cosd(angle);
vy=v*sind(angle);
t=0 ;
y0 = 1000;
i=1;
% Half step
p = p0*exp(-y/y0);
Dx = C*(p*vx^2)/2*A;
Dy = C*(p*vy^2)/2*A;
ax = -Dx/m;
ay = -g-Dy/m;
while y(i) >= 0
    i = i+1;
    vx(i) = vx(i-1)+ax(i-1)*dt;
    vy(i) = vy(i-1)+ay(i-1)*dt;
    
    x(i) = x(i-1)+(vx(i-1)+vx(i))/2*dt;
    y(i) = y(i-1)+(vy(i-1)+vy(i))/2*dt;
    
    p = p0*exp(-y(i)/y0);
    Dx = C*(p*vx(i)^2)/2*A;
    Dy = C*(p*vy(i)^2)/2*A;
    ax(i) = -Dx/m;
    if vy(i)>0
        ay(i) = -g-Dy/m;
    else
        ay(i) = -g+Dy/m;
    end
end
figure;
plot(x,y,'-')
xlabel('x distance (m)')
ylabel('y distance (m)')

[Ozan Terzioglu]

Thanks a lot, I will examine it deeply.

Another question, can I change this code into a Runge-Kutta method?

[Michael Madelaire]

Yes, but you need to find what differential equations you are solving and to what order of Runge-Kutta you want to use.