Velocity change under drag and other forces: small mass, large delta-t

[lukas]

I am working on simulating the movement of small spherical objects under the influence of drag and a number of other forces that, at least for now, don't depend on the object's velocity.
Normally I would sum up all these forces, including drag, to $F_\Sigma$ and update the object's speed (vector) like so: $v_1=\frac{F_\Sigma \cdot \Delta t}{m}$. This seems to work fine for very small $\Delta t$ (e.g., 1e-6 s), but with the following example parameters for only drag and gravity I start running into trouble:
$$
\begin{align*}
\rho &= 997\ \frac{\text{kg}}{\text{m}^3} \\
C_D &= 0.8 \\
r &= 30 \cdot 10^{-9}\text{ m, and } A = \pi r^2 \\
\Delta t &= 0.001\text{ s} \\
v_0 &= 0.05\ \frac{\text{m}}{\text{s}} \\
m &= 2.65 \cdot 10^{-19}\text{ kg} \\
g &= 9.81\ \frac{\text{m}}{\text{s}^2}
\end{align*}
$$
Here, in the first time step, $\Delta v$ would be $\frac{mg \cdot \Delta t}{m} - \frac{\rho A C_D v_0^2 \cdot \Delta t}{2m} = -10.6\ \frac{\text{m}}{\text{s}}$, if I'm not mistaken, which in absolute terms is much higher than the theoretical terminal velocity of 1.5e-3 m/s, to which I would expect the object to decelerate. In the next time step, the velocity is about 478 km/s.

I was thinking that integrating over the drag force or the corresponding acceleration with respect to $t$ I should be able to get a more accurate approximation of the velocity after a large-ish time step. However, since I'm not familiar with differential equations, I started looking for existing solutions. One that stood out to me in particular was this one from StackExchange user ja72. His approach was that the acceleration for a given speed change can be written as
$$
a(v) = A_0 + A_1 v + A_2 v^2,
$$
with, for example, $A_0=g$ for gravity, $A_1=0$, and $A_2=-\frac{1}{2m} \rho A C_d$ for drag. This led him to the following expression:
$$
\Delta v = \frac{A_0 + A_1 v_0 + A_2 v_0^2}{A_1 + 2 A_2 v_0} \left( 1 - \sqrt{ 1 - 2 (A_1 + 2 A_2 v_0) \Delta t} \right).
$$
This appears to work well again for very small time steps and as long as $v_0$ is not too large. If the time step is only slightly too large, the object's velocity rises above the terminal velocity in the first couple of time steps before eventually converging. If either the time step or $v_0$ is too large, the expression under the square root becomes negative.

So here is my question (in three parts):

• Is what I am trying to do possible/feasible at all? I would really prefer to not use $\Delta t$ much lower than a millisecond.
• Is there a better solution than / a better application of the solution that I found?
• Do you perhaps know of a scientific publication of such a solution?

 

[jrmichler]

Your initial time step must be such that the velocity changes by less than 10%, and preferably less than 1%. A general rule in simulations where the results diverge rapidly is to first try smaller time steps until you have at least 5 or 10 time steps per cycle or until it behaves properly. A good rule of thumb is to always make at least one run with a time step significantly shorter than the time step that seems to work.

Also, if you search drag coefficient sphere, you will note that the drag coefficient is a function of Reynolds number. Note also that this approach may not apply in your case, so also search Stoke's Law, which does apply.

 

[lukas]

Thank you, jrmichler! If there is no other way, I guess I'll have to live with longer simulation times or try some other tricks like leaving some computations for only every n-th time step.

Regarding the drag coefficient: I got the value 0.8 for a Reynolds number of 201.5 from figure 4 in
P. P. Brown and D. F. Lawler, "Sphere Drag and Settling VelocityRevisited" https://doi.org/10.1061/(ASCE)0733-9372(2003)129:3(222) (or https://pdfs.semanticscholar.org/969c/09d2677888cdc545b887334ef6b23bc4a0b6.pdf)
Are these conditions sufficient for applying Stoke's Law? Although I'm not sure about the reliability of their citation, the Wikipedia article on drag suggests Reynolds numbers less than 1, and "moving through a fluid at relatively slow speeds" seems to refer to speeds in the order of µm/s.

 

[jrmichler]

When the drag coefficient is a function of Reynolds number, it needs to be a variable in your simulation.

You could try making the time step variable. Use the acceleration to calculate the next time step for a given percent velocity change. Start with about 10% to debug the program, then reduce for more accuracy. Program termination would be based on acceleration decreasing below a set value.

"moving through a fluid at relatively slow speeds" seems to refer to speeds in the order of µm/s.

It's not velocity, it's Reynolds number. And the Reynolds number of a sub-micron (30E-9 m is 0.03 micron) size particle is really low.