1. 弾道計算の理論 #

1. 要約 #

ニュートンの運動方程式の範囲内で、エアガンから飛び出した BB弾は下記の２つの方程式で記述できます。

• BB弾の運動方程式 \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}

• BB弾の回転の減衰を表す方程式 \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}

ここで、 $t$は時間, $\mathbf{r}$は BB弾の位置ベクトル,
$\mathbf{v}, v\equiv|\mathbf{v}|$は BB弾の速度ベクトルとその大きさ,
$\mathbf{V}\equiv \mathbf{v}-\mathbf{u}, V\equiv|\mathbf{V}|$は風と共に動く座標系から見た BB弾の速度ベクトルとその大きさ,
$\mathbf{u}$は風の速度ベクトル,
$\mathbf{e}_z$は鉛直上方向の単位ベクトル,
$\boldsymbol{\omega}, \omega\equiv |\boldsymbol{\omega}|$は BB弾の角速度ベクトルとその大きさ,
$\mathbf{e}_r$はBB弾の中心を原点とする動径方向の単位ベクトル,
$\mathbf{U}\equiv \mathbf{u}-\mathbf{v}-\boldsymbol{\omega}\times R\mathbf{e}_r, U\equiv|\mathbf{U}|$はBB弾の表面におけるBB弾から見た空気の速度ベクトルとその大きさ を表します。

回転の減衰を簡単に扱いたい場合、下記のように近似できます。

\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}

ここで、$c$は$0\sim 1$の適当な定数です(エアガンの計算を行う場合は$0.5$で良いと思います)。

その他の物理量は下記の表のとおりです。下記の物理量を正確に計算したい場合は、空気の特徴量をご覧ください。

2. 定式化 #

2.1. 座標系 #

真っすぐ正面に立ってエアガンを構えたとき、x, y, z 軸を下記のように考えるデカルト座標系で考えます。

• 正面方向を$x$軸正の方向
  （背中の方向は$x$軸負の方向）
• 左方向を$y$軸正の方向
  （右方向を$y$軸負の方向）
• 鉛直上向き方向を$z$軸正の方向
  （地球の中心方向は$z$軸負の方向）

座標の原点は、自分の足と地面に接する点です。 また、$\mathbf{r}=(x, y, z)$を BB弾の位置ベクトルとし、それぞれの方向の単位ベクトルを$\mathbf{e}_x = (1, 0, 0), \mathbf{e}_y = (0, 1, 0), \mathbf{e}_z = (0, 0, 1)$と表します。

まず無風状態の運動方程式を考え、最後に風がある状態の運動方程式を立てていきます。

2.2. BB弾に働く力 #

BB弾の運動方程式は、下記の通り３つの力$\mathbf{F}_g, \mathbf{F}_d, \mathbf{F}_L$から構成されると考えます。

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

ここで、$\mathbf{F}_g$は重力、$\mathbf{F}_d$は空気抵抗力、$\mathbf{F}_L$は回転によって生じる揚力を表します。 BB弾の問題を考えるうえで、コリオリ力は非常に小さいので無視します。

念のため、コリオリ力の大きさ$F_\text{Coriolis}$を計算しておきますと、

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

です。目に見えて運動に影響する重力の大きさは$F_g\approx 2\times 10^{-3}$ですので、その影響の 1/1000 程度な訳です。十分無視できるでしょう。

2.2.1. 重力 #

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

重力について説明は不要でしょう。

2.2.2. 空気抵抗 #

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

完全な球体の抗力係数$C_d$をある程度の範囲でフィッティングさせようとする場合、下記の式が良い近似となります¹。

\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}

ただし、$R_e<2\times 10^{5}$の層流領域のみ有効です。BB弾の場合、$R_e\approx 2\times 10^4$ですので、十分使用可能です。

また$C_M,G_M,H_M$は、マッハ数$M_a=|V|/c$, ($V$は流体中の粒子の相対速度, $c$は空気の音速) を用いて以下の通り与えられます。

\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. 回転による揚力 #

回転による揚力は、回転する円柱の揚力に関する理論であるクッタ・ジュコーフスキーの定理から導きます²。

クッタ・ジュコーフスキーの定理によると、揚力$F_L$の大きさは、循環$\Gamma$を用いて

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

で与えられます。ここで、$\Gamma$は対象の外側での閉経路に対する線積分であり、

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

です。ここで、$\mathbf{v}_f$は対象の周りに生じている流体の速度ベクトルです。

この$\Gamma$は球体の周りに生じる流体上で行わなければならないため、なかなか難しいです。 そのため、対象の周りの循環は球体の循環に比例するだろう、と仮定します。 特にこの仮定に根拠があるわけではないですが、的外れではないと信じて進めていきます。

この仮定の下では、循環$\Gamma$は球体の循環$\Gamma_\text{sp}$の定数倍で書けるとします。つまり、

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

とします。ここで$C_l$は定数で、$\Gamma_\text{sp}$は球(BB弾)の循環です。

半径$R$で角速度$\omega$を持つ球の循環は \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}

と計算できます。ここで$\Gamma(r)$は半径$r$の円盤に対する循環で$\Gamma(r)=2\pi r^2\omega$となります。$l, \phi$は下記の右図の通りです。

以上より、揚力$F_L$の大きさは

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

と書けます。方向も含めれば、下記のようになります。

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

実測結果と比較したところ、BB弾の場合、$C_l\approx 0.12$となります。

2.2.4. 回転の減衰 #

流体中に板を置けば、板は流体から力を受け引きずられるはずです。 BB弾が回転をしているのならば、回転方向とは反対の向きに力のモーメントを受け回転が減衰するはずです。

物体が角速度ベクトル$\mathbf{\omega}$、慣性モーメント$I$で回転する場合、その運動方程式は力のモーメント$\mathbf{N}$を用いて

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

と書けます。 BB弾のような球体の場合、$I=2mR^2/5$で書けます。

$\mathbf{N}$を求める取っ掛かりは、“境界層流れ"という話の中で出てくる摩擦抗力$D_f$です。 流れている流体に置かれた平板に生じる、流体からの摩擦抗力の大きさ$D_f$は

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

で与えられます。ここで、$\rho$は流体の密度、$U$は流体との相対速度、$S_{\text{c}}$は流体と物体の接触面積を表します。また $C_f=C_f(R_e)$は摩擦抗力係数で、下記の値を持ちます³⁴。

\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}

$D_f$が球体に対して使用可能か否かは不明ですが、使用できると信じて計算します。 また、風もある瞬間においては空間で一様であると仮定します。

方針としては、微小区間に働く力のモーメント$\Delta\mathbf{N}=\mathbf{r}\times \Delta\mathbf{f}$を考えて、それを積分することで式\eqref{eqdecay}の右辺を求めます。 まず、微小面積$\Delta S$に働く微小摩擦抗力$\Delta \mathbf{D}_f$は、平板に働く摩擦抗力から

\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}

と書けると仮定します ($U\equiv |\mathbf{U}|$)。ここで、$\mathbf{r}’$は球の中心を原点とする位置ベクトルとします。この$\mathbf{D}_f$ が$\Delta \mathbf{f}$となります。

つまり、球体と流体が接する微小面積$\Delta S$が位置$\mathbf{r’}$で指定されるとき、その位置における球体から見た流体の相対速度$\mathbf{U}(\mathbf{r’})$が分かれば良い、ということになります。

座標は下図のように取ります。

すると位置$\mathbf{r}’$における流体との相対速度$\mathbf{U}(\mathbf{r’})$は、

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

と書けます。ここで$\mathbf{v}’(\mathbf{r}’)$は$\mathbf{r}’$における球体の速度、$\mathbf{v}=\mathbf{v}(t)$はBB弾が平行移動する速度⁵を表します。

また流体からの力は、球体の表面でのみ発生します。つまり、球の中心を原点とする球面座標系において$\mathbf{r’}=R\mathbf{e}_{r’}$であるという条件が入ります。この下では、

\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}

と書くことができます。ここで、$\mathbf{e}_{r’} $は動径方向の単位ベクトルです。

以上から、球体表面の微小面積に働く力のモーメント$\Delta \mathbf{N}$は

\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}

と書けます。半径$R$の球面$\Omega$に渡る表面積分を実施して、球体全体に働く力のモーメント$\mathbf{N}$を求めれば、

\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}

となります。 以上で、球の角速度の変化を示す式\eqref{eqdecay}の右辺を得ました。

2.3. BB弾の運動方程式 #

これまでの議論をまとめますと、BB弾の運動方程式は下記の二本の方程式で記述できることが分かりました。

• BB弾の運動方程式 \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}

• BB弾の回転の減衰を表す方程式 \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}

ここで、BB弾の運動方程式には風の影響が入っていません。流れる空気中の粒子の運動を見るためには、相対速度$\mathbf{V}=\mathbf{v}-\mathbf{u}$を導入すれば良いです。BB弾が飛ぶ時間はおおよそ2秒ほどであり、$\mathbf{u}$は時間に依存しない定ベクトルとして扱います。

以上より、風を含む運動方程式は下記のように書き変えられます。

• BB弾の運動方程式 \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}

• BB弾の回転の減衰を表す方程式 \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}

これが今までの仮定を考えたBB弾の運動方程式です。

この方程式をそのまま解けば良いですが、\eqref{eom1x}は別として\eqref{eom2x}の右辺の積分が面倒です。そのため\eqref{eom2x}の積分を近似することをこれから考えましょう。

2.4. 式の近似 #

\eqref{eom2x}は考えてきた仮定を満たす厳密な式ですが、実際に計算するのは非常に面倒です。
特に右辺の積分は厄介です。その上、回転の減衰もそんなに正確な描像ではないのに、精度よく積分を求める必要はないとも感じます。 そのため荒っぽいけれども実用に足り得る近似を実施しましょう。

ここからはかなりグレーな近似が含まれますので、気に入らない方は\eqref{eom2x}から出発して貴方の近似式を作成していただければと思います。

2.4.1. 積分の除去 #

今、具体的に近似したい部分は下記の積分で、先に結論を言いますと下記のように近似します。

\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}

すると、この近似の下で方程式は

\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}

となります。ここで$c$は$[0,1]$のどこかの値です。例えば$c=1/2$でも良いです⁶。また、途中で$U’\equiv (|\mathbf{u} - \mathbf{v}(t)|^2 + (c R\omega)^2)^{1/2}$と置いています⁷。

更に、式\eqref{w_norm}の右辺を見ると回転の角速度の大きさ$\omega$だけに依存する関数だと分かるので、球面座標系で考えれば微分方程式

\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}

を得ます。詳細はAppendixを参照してください。 以上から、この近似の下では$\omega$以外の球面座標の$\theta, \varphi$方向の成分の変化量はゼロという結論が得られるため、$\boldsymbol{\omega}$の大きさ以外は変化しないという結論が得られます。

上で使用した近似を説明します。使用した近似は、下記の$U(\mathbf{r}’)$の部分で、このような式変形が可能であると仮定しました。

\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}

\eqref{approxU}において、まず第三項は積分を実施する場合に平均的にゼロになるだろうとして無視しました。 そして、第二項の大きさは、正確には$|\boldsymbol{\omega}\times \mathbf{r}’|=\omega R \cdot c(\mathbf{r}’),~~c(\mathbf{r}’)=\sin\theta_{\boldsymbol{\omega}, \mathbf{r}’}$ですが、適当な角度で 良いだろうと安易に考えて$c(\mathbf{r}’)=c$と定数と考えています。

2.4.2. 運動方程式（近似1） #

項2.4.1の近似を仮定した場合、BB弾の運動方程式は下記の二本の方程式で記述できます。

• 仮定

  • BB弾から見た空気の速度は、BB弾表面の至る所でおおよそ平均的に扱える

• BB弾の運動方程式 \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}

• BB弾の回転の減衰 \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}

この近似下でBB弾の回転の減衰は、単に回転の角速度の大きさが徐々に小さくなるだけで、角速度ベクトルの方向は初期状態から変化しません。

2.4.3. 微分方程式の積分 #

さて、ここまでくるともっと近似したくなります。もう少しで解けてしまいそうだからです。

まず、近似としてBB弾の速度の変化は、回転の角速度の変化に対して十分ゆっくり変化すると仮定します。 つまり、

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

が成立すると仮定します。以降の式変形の途中では\eqref{approx1}下で$\mathbf{v}$を定数として扱いますが、最後に$\mathbf{v}$の時間依存性を復活させます。

仮定\eqref{approx1}の下、微分方程式は

\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}

と書けます。 式\eqref{approx2}の方程式は、適当な変数変換によって下記の微分方程式に帰着できます。

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

ここで$\gamma$は定数です。解く方法は割愛しますが、式\eqref{approx4}の解は、

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

です。ここで$c’$は積分定数であり、$\sinh(x)$は双曲線正弦関数 (hyperbolic sine, ハイパボリックサイン)です。 結局、\eqref{approx2}の解は、

\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}

となります。最後に$\mathbf{v}$の時間依存性を復活させて

\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}

が近似的な解となります。再掲しますが、$c$は$0\sim 1$の範囲のどこかの値で、例えば$c=1/2$です。

2.4.4. 運動方程式（近似2） #

下記の近似を仮定した場合、BB弾の運動方程式は下記の二本の方程式で記述できます。

• 仮定

  • BB弾から見た空気の速度は、BB弾表面の至る所でおおよそ平均的に扱える
  • BB弾自体の速度変化$\frac{d\mathbf{v}}{dt}$が、回転の角速度の変化$\frac{d\boldsymbol{\omega}}{dt}$に比べて小さい

• BB弾の運動方程式 \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}

• BB弾の回転の減衰 \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}

ここで、$\omega(t=t_0)=\omega_0$は初期条件を表します。また、$\omega=|\boldsymbol{\omega}|$であり、$\sinh^{-1}(x)$は逆双曲線関数で、$\sinh^{-1}(x)=\ln(x+\sqrt{x^2+1})$と書けます。この近似下で、回転の減衰は単に回転の角速度の大きさが徐々に小さくなるだけで、角速度ベクトルの方向は変化しません。

2.4.5. 近似の評価 #

回転の減衰を近似した場合の評価を行います。 つまり、近似の度合いに応じでどのくらい正確に近似できているかを評価します。正解を式\eqref{eom2x}に従った減衰として、近似式\eqref{App1.B}, \eqref{App2.B}がどの程度有効なのか比較してみましょう。それぞれの近似式は\eqref{eom1x}と連立させて実際に解いていきます。

本稿では、エアガンでBB弾を射出する場合を主に考えるので、パラメータも適したものに固定します。 ここで評価する条件は以下の通りです。

• 条件
  • BB弾: 直径 6mm で 0.25g の弾を使用し, 0.90J (84.9m/s), 300回転/sで高さ1mから地面と平行に射出
  • 環境: 重力加速度 9.8m/s², 温度 20℃, 湿度60%, 気圧 1013.25hPa

2.4.5.1. 通常の計算 #

速度の大きさ $v(t)$	回転角速度の大きさ $\omega(t)$

空気抵抗がある状態で計算した結果です。計算範囲は地面に到達するまでです。

この場合、\eqref{approx1}の近似は成立しないため、近似式\eqref{App2.B}は正確に記述できません。この近似は無理があるということです。よって、回転を示す右の図において、正解の黒線\eqref{eom2x}との差が大きく出ています。

一方で、近似式\eqref{App1.B}は悪くない結果を与えています。多少の差はあれど、近似\eqref{approxU}は有効に働いているのだろうと予測できます。そもそも回転の減衰を計算する仮定が正確ではない（と思う）ため、二重積分を実施する労力よりも、近似式\eqref{App1.B}を採用して回転の減衰を計算した方が良いと思います。

2.4.5.2. 空気抵抗が無い場合で計算 #

速度の大きさ $v(t)$	回転角速度の大きさ $\omega(t)$

続いて、空気抵抗を無くした状態で計算した結果です。この結果は、近似式\eqref{App2.B}が有効な場合を示すために実施しています。 計算範囲は4秒間です。この条件では、マグヌス力が減少してBB弾の軌道が頂点を迎えるより前の間が計算されています。

つまり、マグヌス力と重力が大体釣り合い、速度がほとんど一定になるので、近似の条件\eqref{approx1}が満たされます。 そのため、近似式\eqref{App2.B}は近似式\eqref{App1.B}と良く一致するはずです。 実際に右の図を見ると、近似式\eqref{App1.B}とよく一致していることが分かります。 また、正解の黒線\eqref{eom2x}との差も悪くないと思われます。

以上のことから、近似式\eqref{App1.B}を採用するのが、計算コスト・精度の面から良いと思われます。

3. 結果 #

3.1. ゼロインを40mに合わせた場合 #

ここでは、ゼロインを40mに合わせた場合で、重さの違いによる飛び方の違いを比較します。 計算結果は全て式\eqref{eom1x}と\eqref{eom2x}を連立させて実施した計算です。

図の説明	グラフ
距離 vs 高さ
距離 vs 速さ
距離 vs 回転

• 0.25gの方が、上下方向の振れ幅が小さい
  • 重いBB弾の方が、軌道が安定する。

3.2. 同じ回転量で射出した場合 #

ここでは、射出条件を合わせた場合で、重さの違いによる飛び方の違いを比較します。 計算結果は全て式\eqref{eom1x}と\eqref{eom2x}を連立させて実施した計算です。

図の説明	グラフ
距離 vs 高さ
距離 vs 速さ
距離 vs 回転

• 0.25gの方が、上下方向の振れ幅が小さい

3.3. 同条件で近似具合を比較 #

ここでは、射出条件を合わせた場合で、近似によってどのくらい軌道が歪められるかを比較します。

図の説明	グラフ
距離 vs 高さ
距離 vs 速さ
距離 vs 回転

• 回転量を一定にすると、軌道は合わない。特に、\eqref{App2.B}のずれは大きい。
  • 回転量を変えれば、近似に関わらず一致する軌道が見つかるかもしれません。

3.4. その他の比較 #

その他の計算結果や詳細は、ハイパー道楽様のBB弾の弾道学の結果が見やすいので、そちらをご覧ください。 当時の方程式と変わった点が下記の2点有りますが、結果を大きく変えるものではありません。

• 空気抵抗の抗力係数を求める式を変更
  • 当時の抗力係数$C_d$は、F. A. Morrison⁸による式を使用しておりました。この式によるとエアガンから射出するBB弾の速度域では$C_d \approx 0.40$となります。しかしながら、@rakugaki_camさんの調査⁹により、BB弾の速度域では小さく見積もりすぎていることが分かりました(ありがとうございます！)。 今回使用した新しい式では、乱流域こそ表現できませんがBB弾の速度域では精度よく求められ、$C_d \approx 0.45$となります。この差が若干見えています。
• 回転の減衰の近似方法を変更
  • 運動方程式は当時のもので、本稿の物と見た目上異なっております。これは、球の回転の減衰\eqref{eom2x}から近似式を作る過程で異なる近似を採用した結果であり、そこから導かれる結果は無視できる差(1~5%位?)しか生じません。

4. 実測結果との比較 #

大体以前の結果と一致するので、端折ります。 詳細は、
(https://slpr.sakura.ne.jp/qp/bullet-course)
や
(https://slpr.sakura.ne.jp/qp/bullet-course2)
をご覧ください。

5. Appendix #

5.1. 球面座標系における時間微分 #

\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}

球面座標系における$r, \theta, \varphi$方向の単位ベクトルは

\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}

となります。球面座標系において、任意の位置ベクトル$\mathbf{r}$は$\mathbf{r}=r \mathbf{e}_r$と書くことができます。そのため、球面座標系における$\mathbf{r}$の時間微分は

\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}

と書けます。$\frac{d \mathbf{e}_r}{dt}$は$(x, y, z)$座標で微分して適当に$\mathbf{e}_\theta, \mathbf{e}_\varphi $まとめることで得られます。

6. 修正履歴 #

• 移行：2023 年 7 月 6 日
• 初版：2015 年 1 月 29 日 (https://slpr.sakura.ne.jp/qp/bullet-course)

7. 参考文献 #