Graphing the movement of a projectile under drag

Question

I am trying to write a program to plot the movement of a projectile under the influence of air resistance.

The inputs into the program are the initial velocities of the projectile [$i_x, i_z$] in $x$ and $z$ respectively and the mass in kg ($m$) of the object.

After that, my approach was to calculate the time of flight of the projectile using the formula:

$t_f = (2v_0sin\theta)/g$

I know $v_0 = \sqrt{i_x{}^2, i_z{}^2}$

From there I calculated $\theta = arctan(i_z/i_x)$

Therefore I can calculate $t_f$.

Then the terminal velocity $v_t = m g/c$ where $c$ is momentum, $c = mv_0$. As such, $v_t = 9.8/v0$.

Then I continuously calculated the $x$ and $z$ coordinates by incrementing the time in intervals until the time surpassed $t_f$. The extra result can be fixed by establishing a range of $z \geqslant 0$

We can do this since we know $x(t) = ((v_0v_tcos\theta)/g)(1-e^{-gt/v_t})$ and $z(t) = (v_t/g)(v_0sin\theta + v_t)(1 - e^{-gt/v_t}) - v_tt$

I have tried this and it doesn't work, it produces a graph like this http://prntscr.com/poh0jo for inputs $[35, 40]$ and $1.0$. The shape of the curve is wrong and $z$ goes way beyond $0$ which implies that $t_f$ may be wrong.

I'm also concerned about $v_t$. It doesn't explicitly state that $c$ is momentum. But $v_t$ is in $m/s$, and so I want the right side of the equation to be in $m/s$. Since $v_t=mg/c$, we have $kgm/s^2$ on the top and as such $c$ must be in $kgm/s$ which is momentum. So is $c=mv_0$? This is also confusing since it means we are not actually using $m$ as it implies $v_t = g/v_0$.

Here is my source: http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node29.html#e5.23u

EDIT: Is $c$ the drag constant? If so, does that mean I need to supply more information such as the area?

EDIT2: Should I just use $v_t = \sqrt{(2mg)/(pAC_d)}$ where $C_d$ is the drag coefficient, $p$ is the mass density of air ($1.225 kg/m^3$) and $A$ is the projected area of the projectile? I can figure all of these out and $C_d = (2F_d)/(pu^2A)$ where $F_d$ is the drag force and $u$ is the flow speed of the object relative to the fluid. However, wikipedia then says I need $C_d$ to calculate $F_d$ which means I can't calculate either of them... Not totally sure how to calculate $u$ either.

EDIT3:

$F_d = \dfrac{1}{2}C_dpu^2A$. This will produce a vector since $u$ is a vector of the velocities in $x$ and $z$. As such, we will have $F_{d_x}$ and $F_{d_z}$. We can calculate this since $u = \begin{bmatrix} v_x(0) \\ v_z(0)\end{bmatrix}$, $C_d = 0.5$, $p = 1.225$ and $A = \pi r^2$ where $r$ is known.

Therefore, $\sum f_x = -F_{d_x}$. This means $-F_{d_x} = ma_x(t)$. So, $a_x(t) = \dfrac{-F_{d_x}}{m}$.

This also means $\sum f_z = -F_{d_z} - mg$. And so $-F_{d_z} - mg = ma_z(t)$. Meaning $a_z(t) = \dfrac{-F_{d_z}}{m} - g$.

We know $v_x(0)$ and $v_z(0)$ and therefore we can calculate $a_x(0)$ and $a_z(0)$. This therefore means we can calculate $v_x(0 + \Delta t)$ and $v_z(0 + \Delta t)$ and as such we can calculate $x(0 + \Delta t)$, $z(0 + \Delta t)$.

Then we calculate the new $F_d$ given the new velocity vector. So we can calculate the new accelerations, velocities and therefore coordinates. The formula for $a_z(t)$ will change to $a_z(t) = \dfrac{F_{d_z}}{m} - g$ once $v_z(t) \leq 0$. Then we just keep repeating this until $z(t) < 0$.

Is this correct?

EDIT4:

I have plugged this into my code, and using the values $u = \begin{bmatrix} 60 \\ 20\end{bmatrix}$, $m = 1.0$, $r = 0.5$, $\Delta t = 0.0005$ I have gotten this graph:

Clearly this is not quite correct...

EDIT5:

Should I be using $F_d = \dfrac{1}{2}C_dpu|u|A$? Where $u$ is the velocity vector and $|u| = \sqrt{u_x^2 + u_z^2}$? That seems to be what they're doing here https://pdfs.semanticscholar.org/3fb8/577794f3eb802de98aadc06b0a1120a00c02.pdf. And it results in this graph instead:

Answer 1

In order to graph the projectile motion of a sphere, we must be able to calculate the forces acting on said sphere at a given time. One of these forces is drag, $F_d(t)$. $$F_d(t) = \dfrac{1}{2} \cdot C_d \cdot \rho \cdot v(t) \cdot |v(t)| \cdot A$$

Where $C_d$ is the drag coefficient, $\rho$ is the density of air, $v(t)$ is the velocity vector at time $t$ and $A$ is the cross sectional area of the sphere. Drag is a two dimensional vector such that $$F_d(t) = \begin{bmatrix} F_d^x (t) \\ F_d^z (t)\end{bmatrix}$$

We can therefore work out the sum of forces acting in $x$ and $z$ by using $\sum f = ma$. In this case, $\sum f^x (t) = -F_d^x (t) = ma^x (t)$. As such, $$a^x (t) = \dfrac{-F_d^x (t)}{m}$$

When considering the $z$ component, $\sum f^z (t) = -F_d^z (t) - mg = ma^z (t)$. Therefore, $$a^z(t) = \dfrac{-F_d^z (t)}{m} - g$$

According to Euler's method, we can add a small amount of time $\delta$ (say $0.0005$) and calculate the coordinates at $t + \delta$ where $t$ is the current time. Therefore,

$$v^x (t + \delta) = v^x (t) + a^x (t) \cdot \delta$$ $$v^z (t + \delta) = v^z (t) + a^z (t) \cdot \delta$$

We can then use these formulae to work out the coordinate (displacement) vector, $s(t)$ at time $t$.

$$s^x (t + \delta) = s^x (t) + v^x (t) \cdot \delta$$ $$s^z (t + \delta) = s^z (t) + v^z (t) \cdot \delta$$

We now almost have all of the information we need to begin calculating the coordinates of the motion of the sphere. Although, we must first define some values. $C_d = 0.5$ as this is the standard drag coefficient for a sphere, $\rho = 1.225$ as this is the density of air, let the mass $m = 1.0$ and $A = \pi r^2$ where $r = 0.5$. We know that $c(0) = \begin{bmatrix} 0 \\ 0\end{bmatrix}$, and we can define the initial velocities as $u(0) = \begin{bmatrix} 50 \\ 30\end{bmatrix}$.

We can therefore calculate $F_d(0)$ which will allow us to calculate $a^x (0)$, $a^z(0)$ meaning we can calculate $v(0 + \delta)$, and therefore $c(0 + \delta)$. This process should be repeated until $s^z(t) \leq 0$, as this means we have reached or gone below the $z$ axis. Once this has occurred, we can plot our coordinates and set a range of $z \geq 0$, and our graph will be complete.

Note that we have ignored the $y$ axis as the sphere is being projected directly along the $x$ axis, and so changes only occur in $x$ and $z$.

Thanks to @David K for his help.

---

Below is a proof for the Euler's method equations.

Consider the acceleration vector $a(t)$. Since $a(t) = \dfrac{dv}{dt}$, $a(t) \approx \dfrac{\Delta v}{\Delta t}$ = $\dfrac{v(t + \delta) - v(t)}{\delta}$. Therefore, $$v(t + \delta) = v(t) + a(t) \cdot \delta$$

We can then use this formula to work out the displacement vector, $s(t)$ at time $t + \delta$ which will give us our new coordinates. Consider the velocity vector $v(t)$. As $v(t) = \dfrac{ds}{dt}$, $v(t) \approx \dfrac{\Delta s}{\Delta t} = \dfrac{s(t + \delta) - s(t)}{\delta}$. Therefore, $$s(t + \delta) = s(t) + v(t) \cdot \delta$$

Answer 2

The formula $t_f=(2v_0\sin\theta)/g$ is a formula for projectile motion without air resistance, from one point to another point on a level, non-rotating flat surface. You might want to write $\theta_0$ instead of $\theta,$ because this is the angle of the velocity at the instant the projectile is launched or fired. It is not a function of time.

In the problem this formula is supposed to solve, $v_f = v_0$ by conservation of energy.

In the page you linked to (which is taking some kind of air resistance into account), $c$ is simply an arbitrary constant that can be found by experimentation. It has nothing to do with the initial conditions of the projectile, least of all the initial momentum or $v_0.$ Moreover, let's not gloss over the disclaimer in the first paragraph on that page:

This is not a particularly accurate model of the drag force due to air resistance

Your implied formula for $F_d,$ $$ F_d = \tfrac12 C_d \rho u^2 A, \tag1$$ where $\rho$ is the air density, $u$ is the relative flow speed, and $A$ is the frontal area of the projectile, is much more realistic than the formula $F_d = c v$ implied by the linked page. The constant $C_d,$ however, is not something you can easily derive by simple formulas; it can be approximated by some fancy fluid dynamics, but in reality I think it is generally established experimentally by putting objects of various shapes in an airflow and measuring the forces on them. (Yes, $C_d$ is dependent on the shape of your projectile.)

If you are plotting the trajectory of an artillery shell, which typically travels at supersonic speeds (at least initially), you have an extra complication because Equation $(1)$ is no longer valid if you regard $C_d$ as a constant. But if you're using this to plot the trajectory of a ping-pong ball launched from a homemade catapult, the equation is probably accurate enough. And if you're using a ping-pong ball you can just look up $C_d$ for a sphere. (Actually the $C_d$ of a sphere depends on a number of factors; it's only approximately constant over some intervals of speed. The value $C_d = 0.47$ is probably good for the example of a ping-pong ball.)

For terminal velocity, set $u = v_t$ and $F_d = mg$ in Equation $(1)$ and solve for $v_t$; I think you solved this correctly in Edit 2.

But if you use Equation $(1)$ for the drag on the projectile, however, the equations $x(t) = ((v_0v_t\cos\theta)/g)(1-e^{-gt/v_t})$ and $z(t) = (v_t/g)(v_0\sin\theta + v_t)(1 - e^{-gt/v_t}) - v_tt$ are no longer valid, because they assume the force is directly proportional to speed and not to the square of speed. This brings us to the secret of why that writer chose to use the (unrealistic) assumption of directly proportional drag rather than $v^2$ drag: the differential equation for motion in the $v^2$ case does not have a closed-form solution. That is, you can't just write a formula with a fixed number of terms, plug $t$ into the formula, and come out with an accurate result. You have to use what are called numerical methods for the more realistic problem, which means you would be computing the position and velocity of the projectile at many points along the trajectory, using the position and velocity at the previous point (or multiple previous points) to get the position and velocity at the next point.

If all you want to do is to reproduce what you saw on that web page, just make an arbitrary guess about $v_t$ and then plug this and the initial data $v_0$ and $\theta$ into the formulas for $x(t)$ and $z(t)$.

For the more realistic problem involving Equation $(1),$ the simplest method might be the Euler method. Given the position ($x(t),z(t)$) and velocity ($v_x(t),v_z(t)$) of the projectile at time $t$, you compute the forces on the projectile (gravity and drag) to get the components of acceleration ($a_x(t),a_z(t)$), and then you estimate the position and velocity at time $t+\Delta t$ as follows: \begin{align} x(t+\Delta t) &= x(t) + v_x(t)\Delta t, \\ z(t+\Delta t) &= z(t) + v_z(t)\Delta t, \\ v_x(t+\Delta t) &= v_x(t) + a_x(t)\Delta t, \\ v_z(t+\Delta t) &= v_z(t) + a_z(t)\Delta t. \end{align}

You can start at the time the projectile is launched and proceed forward step by step in this way until you reach a time when $z(t)$ becomes negative, then interpolate between that step and the previous step to estimate where and when the projectile lands. With a small enough $\Delta t,$ this might be reasonably accurate. You can even try re-running the calculation with a smaller value of $\Delta t$ in each step and see how that affects the final result.

---

Note that drag is a force and is therefore a vector, much like velocity, with $x$ and $z$ components. In fact it is a vector in the exact opposite direction from velocity. Equation $(1)$ gives only the magnitude of the vector.

One way to get the components of force is to figure out the angle at which the vector points; let's call that angle $\phi_d.$ If the angle is measured counterclockwise from the positive direction of the $x$ axis, then the components follow the usual formulas: the $x$ component is $(F_d)_x = F_d \cos(\phi_d)$ and the $z$ component is $(F_d)_z = F_d \sin(\phi_d).$ Knowing that the force is exactly opposite the velocity, you could first compute the angle of the velocity, which you already know how to do: it's $\arctan\frac{v_z}{v_x}.$ Then $\phi_d = \arctan\frac{v_z}{v_x} + 180^\circ$; or if you're OK measuring your angles in radians, which is usually a good idea if you're using trig functions in software, $\phi_d = \arctan\frac{v_z}{v_x} + \pi.$ Plugging this into the component formulas you would get \begin{align} (F_d)_x &= F_d \cos\left(\arctan\frac{v_z}{v_x} + \pi\right),\\ (F_d)_z &= F_d \sin\left(\arctan\frac{v_z}{v_x} + \pi\right). \end{align}

If I were doing this, however, I would exploit the fact that I already have components of the velocity, and knowing that the force is exactly opposite, I know that $(F_d)_z/(F_d)_x = v_z/v_x,$ although the signs of $(F_d)_z$ and $(F_d)_x$ are opposite the signs of $v_z$ and $v_x.$ I also know that $v_x^2+v_z^2 = v^2$ (the square of the magnitude of the velocity vector) and $(F_d)_z^2 + (F_d)_x^2 = F_d^2.$ Putting all this together with a little algebra, I would find that

\begin{align} (F_d)_x &= - v_x\frac{F_d}{\sqrt{v_x^2 + v_z^2}},\\ (F_d)_z &= - v_z\frac{F_d}{\sqrt{v_x^2 + v_z^2}}, \end{align}

which I think is a lot simpler than all those trig functions. (You might try proving this to yourself; you can do it with similar triangles, or you can use trig identities on the other equations for the components of drag.)

You should find that the $z$ component of drag is downward (the same direction as the force of gravity) when the projectile is rising and upward (opposing gravity) when the projectile is falling. But the $x$ component is always negative (assuming you launched in a positive $x$ direction); the only horizontal force on the projectile is drag and it always points in the direction we came from. So the horizontal velocity $v_x$ just keeps decreasing; if you launch the projectile high enough, with an initial velocity greater than the terminal velocity, $v_x$ will eventually decrease almost to zero, and then you'll be falling almost straight down at the terminal velocity. But if you launch at less than terminal velocity the projectile will never reach this state before it hits the ground; likewise if you launch at a low angle, it will hit the ground before $v_x$ gets near zero.

---

There are other methods such as Heun's method and Runge-Kutta methods that involve more calculations per step but tend to give more accurate results for the same number of calculations. There are also multistep methods and symplectic methods that might be useful for projectile motion. If you want to explore those possibilities you could look these things up and post new questions here if you encounter things that don't make sense to you.

Any kind of numerical method like these would be an impressive undertaking for someone in high school, I think.