1. Theory of Airsoft ballistics #

1. Summary #

Within the scope of Newton’s equations of motion, a BB pellet launched from an airsoft gun can be described by the following two equations:

• Equation of motion for the BB pellet \begin{align} m\frac{d^2 \mathbf{r}}{dt^2}= -mg\mathbf{e}_z-\frac{1}{2}C_d \rho \pi R^2 V\mathbf{V} -C_l \frac{4}{3}\pi R^3 2\rho \mathbf{V}\times\boldsymbol{\omega} \label{eom} \end{align}

• Damping of the rotation of the BB pellet: \begin{align}\label{eom2} I\frac{d\boldsymbol{\omega}}{dt}= \frac{1}{2}C_f \rho R^3 \int_\Omega d\Omega ~U(\Omega) [\mathbf{e}_{r}\times\mathbf{U}(\Omega)] \end{align}

Here, $t$ denotes time, and $\mathbf{r}$ denotes the position vector of the BB pellet,
$\mathbf{v}$ and $v\equiv|\mathbf{v}|$ denote the velocity vector of the BB pellet and its magnitude respectively,
$\mathbf{V}\equiv \mathbf{v}-\mathbf{u}$, and $V\equiv|\mathbf{V}|$ denote the velocity vector of the BB pellet and its magnitude as seen from a coordinate system moving with the wind,
$\mathbf{u}$ is the wind velocity vector,
$\mathbf{e}_z$ represents the unit vector in the vertical upward direction,
$\boldsymbol{\omega}$, and $\omega\equiv |\boldsymbol{\omega}|$ denote the angular velocity vector of the BB pellet and its magnitude respectively,
$\mathbf{e}_r$ is the unit vector in the radial direction with the center of the BB pellet as the origin,
$\mathbf{U}\equiv \mathbf{u}-\mathbf{v}-\boldsymbol{\omega}\times R\mathbf{e}_r$, and $U\equiv|\mathbf{U}|$ denote the velocity vector of the air as seen from the surface of the BB pellet and its magnitude respectively,
$\Omega$ is the area at the surface of the sphere, $d\Omega = \sin\theta d\theta d\varphi$.

If you want to handle the decay of rotation simply, it can be approximated as follows:

\begin{align} I\frac{d\boldsymbol{\omega}}{dt}\approx -\frac{1}{2}C_f \rho R^2 \frac{8\pi R^2}{3}\left(|\mathbf{u}-\mathbf{v}(t)|^2+(cR\omega )^2 \right)^{1/2}\boldsymbol{\omega} \end{align}

Here, $c$ is an appropriate constant ranging from $0$ to $1$ (for airsoft calculations, $0.5$ is sufficient).

The other physical quantities are as follows. If you want to accurately calculate the following physical quantities, please refer to the Characteristics of Air.

2. Formulation #

2.1. Coordinate #

We consider a Cartesian coordinate system, assuming you are standing straight and holding an airgun, with the x, y and z axes defined as follows:

• The direction straight ahead as the positive x-axis direction
  (The direction behind is the negative x-axis direction)
• The direction to the left as the positive y-axis direction
  (The direction to the right is the negative y-axis direction)
• The direction upward vertically as the positive z-axis direction
  (The direction towards the center of the earth is the negative z-axis direction)

The origin of the coordinates is the point where your feet meet the ground. Furthermore, we define $\mathbf{r}=(x, y, z)$ as the position vector of the BB pellet and the unit vectors in each direction as $\mathbf{e}_x = (1, 0, 0), \mathbf{e}_y = (0, 1, 0), \mathbf{e}_z = (0, 0, 1)$.

We first consider the equation of motion in a windless condition and then formulate the equation of motion when there is wind.

2.2. Forces Acting on the BB Pellet #

The equation of motion of the BB pellet is considered to be composed of three forces, $\mathbf{F}_g, \mathbf{F}_d, \mathbf{F}_L$, as follows:

\begin{align} m\frac{d^2\mathbf{r}}{dt^2} = \mathbf{F}_g + \mathbf{F}_d + \mathbf{F}_L \end{align}

Here, $\mathbf{F}_g$ represents gravity, $\mathbf{F}_d$ represents air resistance, and $\mathbf{F}_L$ represents the lift generated by rotation. In considering the problem of the BB pellet, the Coriolis force is very small and can be ignored.

For the record, if we calculate the magnitude of the Coriolis force, $F_\text{Coriolis}$, we get

$$ F_\text{Coriolis} = 2m|\omega_\text{earth}|v \approx 3\times 10^{-6} $$

The magnitude of gravity, which visibly affects motion, is $F_g\approx 2\times 10^{-3}$, so the Coriolis force is about 1/1000 of that. It can be safely ignored.

2.2.1. Gravity #

\begin{align} \mathbf{F}_g = -m g \mathbf{e}_z \end{align}

There’s no need to explain about gravity.

2.2.2. Air Resistance #

\begin{align} \mathbf{F}_d = -\frac{1}{2}C_d(R_e) \rho \pi R^2 |\mathbf{v}|\mathbf{v} \end{align}

When attempting to fit the drag coefficient $C_d$ of a perfect sphere within a certain range, the following equation provides a good approximation¹.

\begin{align} C_d(R_e) =\frac{24}{R_e}(1+0.15 R_e^{0.687})H_M +\frac{0.42 C_M}{1+(\frac{42500}{R_e^{1.16C_M}})+\frac{G_M}{R_e^{0.5}}} \label{dragCoeff} \end{align}

However, it is only valid in the laminar flow region where $R_e<2\times 10^{5}$. In the case of a BB pellet, $R_e\approx 2\times 10^4$, so it can be used without problems.

Also, $C_M,G_M,H_M$ are given as follows using the Mach number $M_a=|V|/c$, ($V$ is the relative speed of particles in the fluid, $c$ is the speed of sound in air).

\begin{align} C_M &= \Biggl\{ \begin{aligned} &1.65+0.65\tanh(4M_a-3.4) &&\mbox{for } M_a\lt 1.5\\ &2.18-0.13\tanh(0.9M_a-2.7)&&\mbox{for } M_a\gt 1.5 \end{aligned} \\ G_M &= \Biggl\{ \begin{aligned} &166M_a^3+3.29M_a^2-10.9M_a+20 &&\mbox{for } M_a\lt 0.8\\ &5+40M_a^{-3} &&\mbox{for } M_a\gt 0.8 \end{aligned} \\ H_M &= \Biggl\{ \begin{aligned} &0.0239M_a^3+0.212M_a^2-0.074M_a+1 &&\mbox{for } M_a\lt 1\\ &0.93+\frac{1}{3.5+M_a^5} &&\mbox{for } M_a\gt 1 \end{aligned} \end{align}

2.2.3. Lift Due to Rotation #

The lift due to rotation is derived from the Kutta-Joukowski theorem, a theory about the lift of a rotating cylinder².

According to the Kutta-Joukowski theorem, the magnitude of lift $F_L$ is given by using the circulation $\Gamma$,

\begin{align} F_L=\rho \Gamma v. \end{align}

Here, $\Gamma$ is a line integral with respect to a closed path on the outside of the object, and is

\begin{align} \Gamma \equiv \oint_s \mathbf{v}_f\cdot d\mathbf{s}. \end{align}

Here, $\mathbf{v}_f$ is the velocity vector of the fluid generated around the object.

This $\Gamma$ must be performed on the fluid generated around the sphere, which is quite difficult. Therefore, we assume that the circulation around the object will be proportional to the circulation of the sphere. There is no particular basis for this assumption, but we believe it is not far off and proceed.

Under this assumption, we will write the circulation $\Gamma$ as a constant multiple of the circulation $\Gamma_\text{sp}$ of the sphere. In other words,

\begin{align} \Gamma=C_l \Gamma_\text{sp} \end{align}

Here, $C_l$ is a constant and $\Gamma_\text{sp}$ is the circulation of the sphere (BB pellet).

The circulation of a sphere with radius $R$ and angular velocity $\omega$ can be calculated as

\begin{align} \Gamma_{\text{sp}} &= \int_{-R}^{R} \Gamma(r) dl \nonumber \\ &= 2\pi \omega R^3 \int_0^\pi \sin^3\phi d\phi \nonumber \\ &= \frac{4}{3}\pi R^3 \cdot 2\omega \end{align}

Here, $\Gamma(r)$ is the circulation for a disk of radius $r$, and $\Gamma(r)=2\pi r^2\omega$. $l, \phi$ are as shown in the figure on the right below.

Therefore, the magnitude of the lift $F_L$ can be written as

\begin{align} F_L = \rho \Gamma v = C_l \frac{4}{3}\pi R^3 2\omega \rho v \end{align}

Including the direction, it becomes as follows.

\begin{align} \mathbf{F}_L = -C_l \frac{4}{3}\pi R^3 2\rho \mathbf{v}\times\boldsymbol{\omega} \end{align}

Comparing with the actual measurement results, in the case of a BB pellet, $C_l\approx 0.12$.

2.2.4. Damping of Rotation #

If a plate is placed in a fluid, it should be dragged by the force from the fluid. If a BB pellet is rotating, it should receive a moment of force in the direction opposite to the rotation, causing the rotation to dampen.

If an object is rotating with an angular velocity vector $\mathbf{\omega}$ and a moment of inertia $I$, its motion equation can be written as

\begin{align} I\frac{d\boldsymbol{\omega}}{dt}&=\mathbf{N} \label{eqdecay} \end{align}

using the moment of force $\mathbf{N}$. In the case of a sphere like a BB pellet, it can be written as $I=2mR^2/5$.

The starting point to find $\mathbf{N}$ is the frictional drag $D_f$ that comes up in the context of “boundary layer flow”. The magnitude of $D_f$ from the fluid, which arises on a flat plate placed in a flowing fluid, is given by

\begin{align} D_f=C_f\frac{\rho U^2}{2}S_{\text{c}} \end{align}

Here, $\rho$ represents the density of the fluid, $U$ the relative speed with the fluid, and $S_{\text{c}}$ the contact area between the fluid and the object. In addition, $C_f=C_f(R_e)$ is the friction drag coefficient, which has the following values³⁴:

\begin{align} C_f &= \left\{ \begin{aligned} &\frac{1.328}{\sqrt{R_e}} &&(R_e < 5.3\times 10^{5})\\ &\frac{0.455}{(\log_{10}R_e)^{2.58}}-\frac{1700}{R_e} &&(5.3 \times 10^{5} < R_e < 2\times 10^{8})\\ &\frac{0.455}{(\log_{10}R_e)^{2.58}} &&(2\times 10^8 < R_e)\\ \end{aligned} \right. \end{align}

It is unclear whether $D_f$ is usable for a sphere, but we believe it can be used and calculate accordingly. Also, we assume that the wind is uniform in space at a certain moment.

As a strategy, we consider the moment of force $\Delta\mathbf{N}=\mathbf{r}\times \Delta\mathbf{f}$ acting on a small section and calculate the right-hand side of Eq. \eqref{eqdecay} by integrating it. First, we assume that the infinitesimal frictional drag $\Delta \mathbf{D}_f$ acting on the infinitesimal area $\Delta S$ can be written as

\begin{align} \Delta \mathbf{D}_{f}=\frac{\rho C_f}{2} U^2(\mathbf{r’})\frac{\mathbf{U}(\mathbf{r’})}{U(\mathbf{r’})}\Delta S \end{align}

from the frictional drag acting on a flat plate ($U\equiv |\mathbf{U}|$). Here, $\mathbf{r}’$ is a position vector with the center of the sphere as the origin. This $\mathbf{D}_f$ becomes $\Delta \mathbf{f}$.

In other words, when a small area $\Delta S$ where the sphere and fluid make contact is specified at position $\mathbf{r’}$, all we need to know is the relative speed $\mathbf{U}(\mathbf{r’})$ of the fluid from the sphere at that position.

We take coordinates as shown in the figure below.

Then, the relative speed $\mathbf{U}(\mathbf{r’})$ of the fluid at position $\mathbf{r’}$ can be written as

\begin{align} \mathbf{U}(\mathbf{r}’) &= \mathbf{u} - \mathbf{v}’(\mathbf{r}’) \\ &=\mathbf{u} - \boldsymbol{\omega}\times \mathbf{r}’ - \mathbf{v}. \end{align}

Here, $\mathbf{v}’(\mathbf{r}’)$ represents the speed of the sphere at $\mathbf{r}’$, and $\mathbf{v}=\mathbf{v}(t)$ represents the speed at which the BB pellet is moving in parallel⁵.

Also, the force from the fluid only occurs on the surface of the sphere. This means that there is a condition that $\mathbf{r’}=R\mathbf{e}_{r’}$ in the spherical coordinate system with the center of the sphere as the origin. Under this condition, we can write

\begin{align} \mathbf{r}’ &= R\mathbf{e}_{r’} \\ &= R \begin{pmatrix} \sin\theta \cos\varphi \\ \sin\theta \sin\varphi \\ \cos\theta \end{pmatrix} \end{align}

where $\mathbf{e}_{r’}$ is the unit vector in the radial direction.

From the above, the moment of force $\Delta \mathbf{N}$ acting on the infinitesimal area of the sphere surface can be written as

\begin{align} \Delta\mathbf{N}&=\mathbf{r’}\times \Delta\mathbf{f} \nonumber \\ &=\frac{1}{2}C_f \rho U(\mathbf{r’}) [\mathbf{r’} \times\mathbf{U}(\mathbf{r’})] \Delta S \end{align}

By performing the surface integral over the sphere of radius $R$, $\Omega$, and determining the moment of force $\mathbf{N}$ acting on the whole sphere, we get

\begin{align} \mathbf{N} = \frac{1}{2}C_f \rho R^3 \int_\Omega ~U(\Omega) [\mathbf{e}_{r’}\times\mathbf{U}(\Omega)] \end{align}

With this, we have obtained the right-hand side of the Eq. \eqref{eqdecay} that represents the change in angular velocity of the sphere.

2.3. Equation of Motion of the BB Pellet #

Summing up our discussion so far, we found that the equation of motion of the BB pellet can be described by the following two equations:

• Equation of motion of the BB pellet \begin{align} m\frac{d^2 \mathbf{r}}{dt^2}= -mg\mathbf{e}_z-\frac{1}{2}C_d \rho \pi R^2 v\mathbf{v} -C_l \frac{4}{3}\pi R^3 2\rho \mathbf{v}\times\boldsymbol{\omega} \end{align}

• Equation representing the damping of the BB pellet’s rotation \begin{align} I\frac{d\boldsymbol{\omega}}{dt}= \frac{1}{2}C_f \rho R^3 \int_\Omega d\Omega ~U(\Omega) [\mathbf{e}_{r}\times\mathbf{U}(\Omega)] \end{align}

Note that the effect of the wind is not included in the equation of motion of the BB pellet. To see the motion of particles in flowing air, we can introduce the relative velocity $\mathbf{V}=\mathbf{v}-\mathbf{u}$. The flight time of the BB pellet is approximately 2 seconds, and $\mathbf{u}$ is treated as a constant vector not dependent on time.

Thus, the equation of motion including the wind can be rewritten as follows:

• Equation of motion of the BB pellet \begin{align} \tag{Ex. A}\label{eom1x} m\frac{d^2 \mathbf{r}}{dt^2}= -mg\mathbf{e}_z-\frac{1}{2}C_d \rho \pi R^2 V\mathbf{V} -C_l \frac{4}{3}\pi R^3 2\rho \mathbf{V}\times\boldsymbol{\omega} \end{align}

• Damping of the rotation of the BB pellet: \begin{align} \tag{Ex. B}\label{eom2x} I\frac{d\boldsymbol{\omega}}{dt}= \frac{1}{2}C_f \rho R^3 \int_\Omega d\Omega ~U(\Omega) [\mathbf{e}_{r}\times\mathbf{U}(\Omega)] \end{align}

This is the equation of motion of the BB pellet considering our assumptions so far.

Although this equation should be solved as it is, the integration on the right-hand side of \eqref{eom2x} is troublesome, separate from Eq. \eqref{eom1x}. Therefore, let’s consider approximating the integration of Eq. \eqref{eom2x} from here.

2.4. Approximation of the Equations #

Though Eq. \eqref{eom2x} is a exact equation satisfying our previous assumptions, it is quite cumbersome to calculate in practice. The integration on the right hand side is particularly tricky. Given the imprecision of our depiction of rotational decay, it seems unnecessary to compute the integral precisely. Thus, we shall proceed with a rough yet practical approximation.

Be aware that the following discussion includes some rather loose approximations. If these are unsatisfactory, you’re welcome to devise your own approximation starting from Eq. \eqref{eom2x}.

2.4.1. Removal of the Integral #

The specific part we want to approximate is the integral below. Jumping to the conclusion, we’ll approximate it as follows:

\begin{align} \int_\Omega ~U(\Omega) [\mathbf{e}_r\times\mathbf{U}(\Omega)] d\Omega &= \int_0^{\pi} d\theta\int_0^{2\pi} d\varphi ~U(\mathbf{r}’) [\mathbf{e}’\times\mathbf{U}(\mathbf{r}’)]\sin\theta\nonumber \\ &\approx U’ \int_0^{\pi} d\theta\int_0^{2\pi} d\varphi [\mathbf{e}’\times\mathbf{U}(\mathbf{r}’)]\sin\theta \nonumber \\ &= U’ \left(-\frac{8\pi R}{3}\boldsymbol{\omega} \right)\nonumber \\ &=-\frac{8\pi R}{3}\left[|\mathbf{u}-\mathbf{v}(t)|^2 + (cR\omega )^2\right]^{1/2}\boldsymbol{\omega} \end{align}

With this approximation, the equation becomes

\begin{align}\label{w_norm} I\frac{d\boldsymbol{\omega}}{dt}\approx -\frac{1}{2}C_f \rho R^2 \frac{8\pi R^2}{3}\left(|\mathbf{u}-\mathbf{v}(t)|^2+(cR\omega )^2 \right)^{1/2}\boldsymbol{\omega}. \end{align}

where $c$ is a value somewhere in $[0,1]$. For instance, $c=1/2$ would be acceptable⁶. In the middle of the process, we have set $U’\equiv (|\mathbf{u} - \mathbf{v}(t)|^2 + (c R\omega)^2)^{1/2}$⁷.

Further, from the right hand side of Eq. \eqref{w_norm}, we realize it is a function depending solely on the magnitude of angular velocity $\omega$. Considering this in spherical coordinates, we derive the differential equation

\begin{align} \label{approx3} I\frac{d\omega}{dt}=-\frac{1}{2}C_f \rho R^2 \frac{8\pi R^2}{3}\left(|\mathbf{u}-\mathbf{v}|^2+(cR\omega )^2 \right)^{1/2}\omega \end{align}

Please refer to the Appendix for details.
From the above, under this approximation, we conclude that the changes in $\theta, \varphi$ components of spherical coordinates are zero, other than the magnitude of $\boldsymbol{\omega}$. Thus, we deduce that the magnitude of $\boldsymbol{\omega}$ doesn’t change except for its magnitude.

The approximation used in the above argument is explained here. The approximation was used on the part of $U(\mathbf{r}’)$ below, and it was assumed that such a transformation of the equation would be possible.

\begin{align} {U}(\mathbf{r}’) &= |\mathbf{u} - \mathbf{v}(t) - \boldsymbol{\omega}\times \mathbf{r}’| \nonumber\\ &= [|\mathbf{u} - \mathbf{v}(t)|^2 + |\boldsymbol{\omega}\times \mathbf{r}’|^2 - 2 (\mathbf{u} - \mathbf{v}(t))\cdot(\boldsymbol{\omega}\times \mathbf{r}’)]^{1/2} \nonumber\\ &\approx [|\mathbf{u} - \mathbf{v}(t)|^2 + (c R \omega )^2]^{1/2}\label{approxU} \end{align}

In Eq. \eqref{approxU}, first, the third term is expected to be zero on average when performing the integration, so it’s ignored. The magnitude of the second term is exactly $|\boldsymbol{\omega}\times \mathbf{r}’|=\omega R \cdot c(\mathbf{r}’),~~c(\mathbf{r}’)=\sin\theta_{\boldsymbol{\omega}, \mathbf{r}’}$, but a suitable constant angle was assumed and considered as $c(\mathbf{r}’)=c$.

2.4.2. Equation of motion (Approximation 1) #

Assuming the approximation from section 2.4.1, the equation of motion of the BB pellet can be expressed as the following pair of equations:

• Assumption

  • The velocity of the air as seen from the BB pellet can be treated roughly as an average over the entire surface of the BB pellet.

• Equation of motion of the BB pellet \begin{align} \tag{App1. A}\label{App1.A} m\frac{d^2 \mathbf{r}}{dt^2}= -mg\mathbf{e}_z-\frac{1}{2}C_d \rho \pi R^2 V\mathbf{V} -C_l \frac{4}{3}\pi R^3 2\rho \mathbf{V}\times\boldsymbol{\omega} \end{align}

• Damping of the rotation of the BB pellet: \begin{align} \tag{App1. B}\label{App1.B} I\frac{d\omega}{dt}=-\frac{1}{2}C_f \rho R^2 \frac{8\pi R^2}{3}\left(|\mathbf{u}-\mathbf{v}(t)|^2+(cR\omega )^2 \right)^{1/2}\omega \end{align}

Under this approximation, the rotational decay of the BB pellet is simply that the magnitude of the rotational speed gradually decreases, while the direction of the angular velocity vector does not change from the initial state.

2.4.3. Integration of Differential Equations #

At this point, we’re eager to develop a more accurate approximation. It seems we’re just on the brink of a solution.

Initially, we make an approximation that the change in velocity of the BB pellet is sufficiently slow compared to the change in angular velocity. In other words, we assume that the following equation holds:

\begin{align} \label{approx1} \mathbf{v}(t)\approx \mathbf{v} \end{align}

During the intermediate steps of our calculations, we treat $\mathbf{v}$ as a constant under the assumption \eqref{approx1}, but ultimately reintroduce the time-dependence of $\mathbf{v}$.

Under the assumption \eqref{approx1}, the differential equations can be written as:

\begin{align} \label{approx2} I\frac{d\omega}{dt}= -\frac{1}{2}C_f \rho R^2 \frac{8\pi R^2}{3}\left(|\mathbf{u}-\mathbf{v}|^2+(cR\omega )^2 \right)^{1/2}\omega \end{align}

With an appropriate variable substitution, Eq. \eqref{approx2} can be reduced to the following differential equation:

\begin{align} \label{approx4} \frac{dx}{dt} = -\gamma (1+x^2)^{1/2} x,~~(\gamma, x \geq 0) \end{align}

Here, $\gamma$ is a constant. The method to solve this is not discussed here, but the solution to Eq. \eqref{approx4} is:

\begin{align} \label{approx5} x(t) = \frac{1}{\sinh(\gamma t + c’)} \end{align}

Here, $c’$ is the constant of integration, and $\sinh(x)$ is the hyperbolic sine function. Consequently, the solution to Eq. \eqref{approx2} is:

\begin{align} \label{approx6} \omega(t) \approx \frac{|\mathbf{u}-\mathbf{v}|}{cR} \frac{1}{\sinh\left(\frac{1}{2I}C_f \rho R^2 \frac{8\pi R^2}{3}|\mathbf{u}-\mathbf{v}|\cdot t + c’\right)} \end{align}

Finally, reintroducing the time-dependence of $\mathbf{v}$ results in an approximate solution:

\begin{align} \omega(t) \approx \frac{|\mathbf{u}-\mathbf{v}(t)|}{cR} \frac{1}{\sinh\left(\gamma(t)\cdot t + c’\right)} \\ \gamma(t) \equiv \frac{1}{2I}C_f \rho R^2 \frac{8\pi R^2}{3}|\mathbf{u}-\mathbf{v}(t)| \end{align}

2.4.4. Equations of Motion (Approximation 2) #

Assuming the approximations listed below, the equations of motion for the BB pellet can be represented by the following two equations:

• Assumptions:

  • The speed of the air as viewed from the BB pellet can be treated as approximately uniform across the surface of the BB pellet.
  • The change in the velocity of the BB pellet, $\frac{d\mathbf{v}}{dt}$, is small compared to the change in angular velocity, $\frac{d\boldsymbol{\omega}}{dt}$.

• Equations of motion for the BB pellet: \begin{align} \tag{App2. A}\label{App2.A} m\frac{d^2 \mathbf{r}}{dt^2}= -mg\mathbf{e}_z-\frac{1}{2}C_d \rho \pi R^2 V\mathbf{V} -C_l \frac{4}{3}\pi R^3 2\rho \mathbf{V}\times\boldsymbol{\omega} \end{align}

• Damping of the rotation of the BB pellet: \begin{align} \tag{App2. B}\label{App2.B} \omega(t) &= \frac{|\mathbf{u}-\mathbf{v}(t)|}{cR} \frac{1}{\sinh\left(\gamma(t)\cdot t + c’\right)} \\ &\gamma(t) = \frac{1}{2I}C_f \rho R^2 \frac{8\pi R^2}{3}|\mathbf{u}-\mathbf{v}(t)| \tag{App2. Ba}\\ &c’ = \sinh^{-1}\left(\frac{|\mathbf{u}-\mathbf{v}(t_0)|}{cR\omega_0}\right)-\gamma(t_0)\cdot t_0\tag{App2. Bb} \end{align}

Here, $\omega(t=t_0)=\omega_0$ represents the initial condition. Also, $\omega=|\boldsymbol{\omega}|$, and $\sinh^{-1}(x)$ is the inverse hyperbolic sine function, which can be written as $\sinh^{-1}(x)=\ln(x+\sqrt{x^2+1})$. Under this approximation, the decay of the rotation simply means the magnitude of the angular velocity decreases gradually, with no change in the direction of the angular velocity vector.

2.4.5. Evaluation of Approximations #

We will evaluate the case of approximating the rotational decay. Specifically, we will determine how accurately the approximation can be made according to the degree of approximation. Assuming the correct answer is the decay according to Eq. \eqref{eom2x}, let’s compare how effective the approximation formulas \eqref{App1.B} and \eqref{App2.B} are. Each of these approximation formulas will be actually solved in conjunction with Eq. \eqref{eom1x}.

In this paper, as we mainly consider the case of firing BB pellets with an air gun, the parameters will also be fixed to suitable ones. The conditions for evaluation here are as follows:

• Conditions
  • BB pellets: Using a pellet of 0.25g with a diameter of 6mm, launched horizontally from a height of 1m at 0.90J (84.9m/s), with 300 rotations/s
  • Environment: Gravitational acceleration 9.8m/s², Temperature 20℃, Humidity 60%, Atmospheric pressure 1013.25hPa

2.4.5.1. Regular Calculations #

Velocity magnitude $v(t)$	Angular velocity magnitude $\omega(t)$

This is the result of calculations in a state with air resistance. The calculation range is until reaching the ground.

In this case, the approximation \eqref{approx1} is not valid, so the approximation \eqref{App2.B} cannot accurately describe it. It means that this approximation is unreasonable. Therefore, there is a large difference from the correct black line \eqref{eom2x} in the right figure indicating rotation.

On the other hand, the approximation \eqref{App1.B} provides a not bad result. Although there is some difference, it can be predicted that the approximation \eqref{approxU} is effectively working. Since the assumption to calculate the rotational decay is not accurate (I think), it would be better to adopt the approximation formula \eqref{App1.B} and calculate the rotational decay rather than doing the effort of double integration.

2.4.5.2. Calculations in Case of No Air Resistance #

Velocity magnitude $v(t)$	Angular velocity magnitude $\omega(t)$

Next, this is the result of calculations in a state without air resistance. This result is done to show the case where the approximation formula \eqref{App2.B} is effective. The calculation range is for 4 seconds. Under this condition, the Magnus force decreases and the trajectory of the BB pellet is calculated before it reaches its peak.

In other words, the Magnus force and gravity balance out roughly, the speed is almost constant, so the approximation condition \eqref{approx1} is met. Therefore, the approximation formula \eqref{App2.B} should match well with the approximation formula \eqref{App1.B}. In fact, looking at the right figure, you can see that it matches well with the approximation formula \eqref{App1.B}. The difference from the correct black line \eqref{eom2x} also seems not bad.

From the above, it is considered that it is good to adopt the approximation formula \eqref{App1.B} in terms of calculation cost and accuracy.

3. Results #

3.1. When Zeroing is set at 40m #

Here, we compare the difference in the way they fly due to differences in weight when zeroing is set to 40m. All calculated results are calculations performed by combining equations \eqref{eom1x} and \eqref{eom2x}.

Figure description	Graph
Distance vs Height
Distance vs Velocity
Distance vs Rotation

• The 0.25g has smaller swings in the vertical direction
  • The heavier BB pellet has a more stable trajectory.

3.2. Comparisons with Same Amount of Spin #

In this section, we compare the trajectories of BBs with different weights, assuming the same launch conditions. The calculations were carried out by solving equations \eqref{eom1x} and \eqref{eom2x} simultaneously.

Graph Description	Graph
Distance vs Height
Distance vs Velocity
Distance vs Rotation

• The BB with a weight of 0.25g shows less variation in vertical oscillation.
  

3.3. Comparison of Approximations Under the Same Conditions #

Here we compare how the trajectory is distorted by the approximation when the launch conditions are the same.

Graph Description	Graph
Distance vs Height
Distance vs Velocity
Distance vs Rotation

• When the amount of spin is kept constant, the trajectories do not match. Particularly, the deviation of \eqref{App2.B} is significant.
  • If the spin is changed, there might be a trajectory that matches regardless of the approximation.

3.4. Other Comparisons #

For other calculation results and details, please see the results of BB pellet ballistics by Hyper Doraku in Japanese, which are easy to understand. There are two differences from the equations used at that time, but they do not significantly affect the results.

• The formula to determine the drag coefficient of air resistance has been changed.

  • The drag coefficient $C_d$ used at that time was based on the formula by F. A. Morrison⁸. According to this formula, the speed range of BB pellets fired from airsoft is $C_d \approx 0.40$. However, thanks to the investigation by @rakugaki_cam, it was found that the speed range of BB pellets is underestimated. The new formula used this time cannot express the turbulent flow region, but it can accurately estimate the speed range of BB pellets, resulting in $C_d \approx 0.45$. This difference is slightly visible.

• The method of approximating the attenuation of spin has been changed.

  • The equation of motion is different from the one used at that time and the one in this paper. This is due to the different approximations used in creating the approximate equation from the attenuation of the sphere’s rotation \eqref{eom2x}. The results derived from this cause negligible differences (about 1~5%).

4. Comparison with Actual Measurement Results #

As the results largely match those from previous resultss, I’ll omit detailed comparisons here. For more details, please refer to
(https://slpr.sakura.ne.jp/qp/bullet-course)
and
(https://slpr.sakura.ne.jp/qp/bullet-course2)
in Japanese.

5. Appendix #

5.1. Time Derivative in Spherical Coordinate System #

\begin{align} \left\{ \begin{aligned} x&=r\sin\theta\cos\varphi \\ y&=r\sin\theta\sin\varphi \\ z&=r\cos\theta \end{aligned} \right. \end{align}

The unit vectors in the $r, \theta, \varphi$ directions in the spherical coordinate system are:

\begin{align} \left\{ \begin{aligned} \mathbf{e}_r &= \mathbf{e}_x\sin\theta\cos\varphi + \mathbf{e}_y \sin\theta\sin\varphi + \mathbf{e}_z \cos\theta \\ \mathbf{e}_\theta &= \mathbf{e}_x\cos\theta\cos\varphi + \mathbf{e}_y \cos\theta\sin\varphi - \mathbf{e}_z \sin\theta\\ \mathbf{e}_\varphi &= -\mathbf{e}_x\sin\varphi + \mathbf{e}_y\cos\varphi \end{aligned} \right. \end{align}

In spherical coordinates, any position vector $\mathbf{r}$ can be written as $\mathbf{r}=r \mathbf{e}_r$. Therefore, the time derivative of $\mathbf{r}$ in spherical coordinates can be written as: \begin{align} \frac{d\mathbf{r}}{dt} &= \frac{d r\mathbf{e}_r}{dt} \\ &=\frac{dr}{dt}\mathbf{e}_r + r\frac{d \mathbf{e}_r}{dt} \\ &=\frac{dr}{dt}\mathbf{e}_r + \frac{d\theta}{dt}\mathbf{e}_\theta + \frac{d\varphi}{dt}\sin\theta \mathbf{e}_\varphi \end{align}

The value $\frac{d \mathbf{e}_r}{dt}$ can be obtained by differentiating in $(x, y, z)$ coordinates and conveniently grouping into $\mathbf{e}_\theta, \mathbf{e}_\varphi$.

6. Revision History #

• Migration: July 6, 2023
• First Edition: January 29, 2015 (https://slpr.sakura.ne.jp/qp/bullet-course)

7. References #