airdragnewtonian-mechanics
I'm simulating snowfall in a discrete-time physics engine, using this formula for the drag force on a snowflake:
$F_d=CAv^2$ (where $A$ is the flake's cross-sectional area and $C$ is a catch-all constant). I need to calculate the drag impulse over a given $dt$, and that's where I'm stuck.
I've read this similar question, but it's been so long since I've dealt with differential equations that I can't follow the accepted answer well enough to adapt it for my purposes.
$C_d$ is a function of speed via the Reynolds number. See here and here.
For some numeric values for the $C_d$ vs. $\rm Re$ use the following from here of which you take the log values to interpolate.
The kinematic viscosity of Air is $\nu = 14.8\; \rm{cSt} = 14.8 \cdot 10^{−6}\; \rm{m^2/s}$. At a speed of $v=20\;\rm{m/s}$ a ball with diameter $d=5\,\rm{cm} = 0.05\;\rm{m}$ has Reynolds number of $\rm{Re} = \frac{v d}{\nu} = 67400 $.
When you look at the $C_d$ vs. $\rm Re$ graph you get $C_d = 0.47$. The area of the ball is $A=\pi \frac{d^2}{4} = 0.001964\;\rm{m^2}$ so the drag force is
$F_d = \frac{1}{2} \rho A C_d v^2 = \frac{1}{2} (1.2\;\rm{kg/m^3}) (0.001964\;\rm{m^2}) (0.47) (20\;\rm{m/s})^2 = 0.2215\;\rm{N}$
The force F is taken to be positive if it is upward and negative if it is downward. So -mg means that the gravitational force is downward. The term -kv means that, if the body is moving upward (positive v) the drag force is downward, and if the body is moving downward (negative v), the drag force is upward. This is completely consistent with our expectations.
Best Answer
If you model the acceleration of a particle for small changes in speed as
$$ a(v) = A_0 + A_1 v + A_2 v^2 $$
with appropriate constants (like $A_0$ with gravity and $A_2$ with drag) then you can calculate the position $\Delta x$ and velocity change $\Delta v$ from the current speed $v_1$ ( $v \rightarrow v_1 + \Delta v$) as
$$ \Delta t = \int_{v_1}^{v_1+\Delta v} \frac{1}{a(v)}\,{\rm d} v \\ \Delta x = \int_{v_1}^{v_1+\Delta v} \frac{v}{a(v)}\,{\rm d} v $$
The result is well defined, but quite complex. You can simplify things a bit for simulations by only considering 2nd order terms of $\Delta v$. Simplifications can be made then using the constant $a_1 = a(v_1) = A_0 + A_1 v_1 + A_2 v_2^2$
$$ \Delta t = \frac{v-v_1}{a_1} - \frac{(v-v_1)^2}{2 a_1^2} (A_1 + 2 A_2 v_1) \\ \Delta x = \frac{v^2-v_1^2}{2 a_1} - \frac{v_1 (v-v_1)^2}{2 a_1^2} (A_1 + 2 A_2 v_1) $$
Example
Initial Conditions: $$a(v) = 10 - 1 v - \frac{1}{5} v^2 \\ v_1 = 3 \\ a_1 = 10 - 1 3 - \frac{9}{5} = \frac{26}{5} = 5.2 $$ Step $\Delta v = 0.5$, $v = v_1 + \Delta v = 3.5$ with exact solution $\Delta t = 0.10632$, $\Delta x = 0.3430$, and approximate solution
$$ \Delta t = \frac{5 (v-3)}{26} + \frac{55 (v-3)^3}{1352} = 0.10632 \\ \Delta x = \frac{5 (v^2-9)}{52} + \frac{165 (v-3)^3}{1352} = 0.3430 $$
so the simulation (for example) would move from $$ \begin{Bmatrix} t_1=10 \\ x_1=1 \\ v_1 = 3 \\ a_1 = 5.2 \end{Bmatrix} \rightarrow \begin{Bmatrix} t_2=10.1053 \\ x_2=1.3430 \\ v_2 = 3.5 \\ a_2 = 4.05 \end{Bmatrix}$$
Impulse Calculation
Now you can turn the solution around and state that
$$ \Delta v = \frac{A_0 + A_1 v_1 + A_2 v_1^2}{A_1 + 2 A_2 v_1} \left( 1 - \sqrt{ 1 - 2 (A_1 + 2 A_2 v_1) \Delta t} \right) $$
which is a change is velocity over some time step, with impulse $J = m \Delta v$.