0
$\begingroup$

I'm trying to simulate the motion of a projectile taking into account aerodynamic drag. The code works perfectly if the aerodynamic drag is zero. However, if the drag coefficient or velocity is too high, the trajectory starts to bend towards the negative x axis. I used this source to help me create the simulation.

Here's the code:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

cd = 0.1
v = 150
ang = 40 # in degrees
t = [0]  # list to keep track of time
vx = [v * np.cos(ang / 180 * np.pi)] 
vy = [v * np.sin(ang / 180 * np.pi)]
x = [0]  # list for x and y position
y = [0]
g = 9.8
M = 1
density = 1.2
area = 1
# Drag force
drag = area * density * cd * 0.5 * v ** 2
# Acceleration x and y
ax = [-(drag * np.cos(ang / 180 * np.pi)) / M]
ay = [-g - (drag * np.sin(ang / 180 * np.pi) / M)]
# time-step
dt = 0.01
f1 = plt.figure()
ax1 = f1.add_subplot(111)
line, = ax1.plot([], [], '-')
def animate(i):
    t.append(t[i] + dt) 
    vx.append(vx[i] + dt * ax[i])  # Update the velocity
    vy.append(vy[i] + dt * ay[i])
    # Update position
    x.append(x[i] + dt * vx[i])
    y.append(y[i] + dt * vy[i])
    # Calculate updated velocity
    vel = np.sqrt(vx[-1] ** 2 + vy[-1] ** 2)  # magnitude of velocity
    drag = area * density * 0.5 * cd * vel ** 2
    ax.append(-(drag * np.cos(ang / 180 * np.pi)) / M)
    ay.append(-g - (drag * np.sin(ang / 180 * np.pi) / M))
    line.set_xdata(x[:-1])
    line.set_ydata(y[:-1])
    return line,

ax1.set_ylabel("y (meters)")
ax1.set_xlabel("x (meters)")
ax1.set_xlim(0, 90)
ax1.set_ylim(0, 90)
ax1.autoscale(False)
ani = animation.FuncAnimation(f1, animate, interval=0.1, blit=True)
plt.show()

This is the output:

enter image description here

It looks pretty accurate up until the peak of the trajectory, where eventually it starts bending inward which is obviously incorrect. Why is this happening? Is there an error in the script, or is it simply the result of Euler's method being inaccurate?

$\endgroup$
1
  • $\begingroup$ Something looks wrong with the direction of the drag force - when the projectile starts falling vertically (near the point x=40,y=10) the net force (gravity+drag) can be only vertical; but according to the picture there is some component of the net force directed in the negative X direction. $\endgroup$ – Maxim Umansky Dec 16 '20 at 22:48
3
$\begingroup$

You are not updating the angle in your simulation as the velocity vector is changing. Here is a modified version of your 'animate' function and the result.

def animate(i):
    t.append(t[i] + dt) 
    vx.append(vx[i] + dt * ax[i])  # Update the velocity
    vy.append(vy[i] + dt * ay[i])
    vel = np.sqrt(vx[-1] ** 2 + vy[-1] ** 2)  # magnitude of velocity
    cos = vx[-1]/vel # note: should deal with vel=0 properly
    sin = vy[-1]/vel
    # Update position
    x.append(x[i] + dt * vx[i])
    y.append(y[i] + dt * vy[i])
    # Calculate updated velocity
    drag = area * density * 0.5 * cd * vel ** 2
    ax.append(-(drag * cos) / M)
    ay.append(-g - (drag * sin / M))
    line.set_xdata(x[:-1])
    line.set_ydata(y[:-1])
    return line,

and the result:

enter image description here

$\endgroup$
3
  • $\begingroup$ Thanks. It works perfectly if the values of the drag coefficient and density are low. If the cd or density value is too high, the graph glitches out, just a straight line appears and a RuntimeWarning: error. I believe this is because the vel value becomes 0 if it goes too slow, and as a result it is divided by zero. However, the issue gets solved if I input a low value for the initial velocity. $\endgroup$ – Star Man Dec 17 '20 at 0:54
  • $\begingroup$ Yes, I slipped a comment in the code that suggested dealing with vel=0 properly. A simple ‘if’ statement would be a crude but effective way to handle it. $\endgroup$ – rpm2718 Dec 17 '20 at 1:04
  • $\begingroup$ Actually, that straight line sounds like a different problem - that might be numerical instability if your deceleration is too high, causing the velocity change between time steps to be too large. Try decreasing the time increment. $\endgroup$ – rpm2718 Dec 17 '20 at 13:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.