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:
Best Answer
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$$