- #1

lukas

- 2

- 0

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):

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?

Last edited: