「弾道計算」カテゴリーアーカイブ

水中下でのBB弾の弾道計算

水中にBB弾が入った時の軌道は計算できるのでしょうか?
もしも水中での軌道が分かったのならば、そこから回転量、初速が分かるはずです。

まとめ


目的

  • ハイスピードカメラを用いず、空気中で50mに及ぶ軌道を追わずに回転量を知ることが出来るか?
  • 弾道計算(BB弾)の理論の検証

方法

  • 水槽に向けてエアガンを撃った軌道とBB弾の弾道計算理論による軌道と比較

結論

  • 水中の軌道は定性的に一致
  • 定量的に合うのは到達距離、合わないのは水中での揚力
  • 水中の軌道から回転量を知るには理論的に不十分
  1. まとめ
  2. きっかけ
  3. 方針
    1. 理論上の問題
    2. 実験上の問題
  4. 実験と理論の比較
  5. 理論で回転の減衰を入れないと?
  6. 参考文献

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算←今ここ


スポンサーリンク


きっかけ


つい先日、こんな動画を見つけました。
夏休み 自由研究? BB弾水に撃つと?検証してみた。 マック堺 エアガンレビュー

水中の軌道じゃないですか。面白そう!
計算できない理由は無い。
式もある、プログラムもある。ならやってみよう。

がきっかけです。何故計算しようと思ったか?素晴らしい理由はありません。そこに物理があるからです。

方針


厳密さはそこまで求めません。

これには理由が2つあり、

  1. 水中を想定して作った理論ではないため、水の揺らぎが与える影響が不明。
  2. 実験がされた動画の縮尺、距離、回転量を見積もれない。

だからです。詳しい理由は以下の通り。

理論上の問題

空気も水も同じ流体ですが、その特性は大きく違います。

大気中の密度は低く、BB弾の密度と比べても0.1%未満です。そのため、大気の揺らぎがBB弾の影響することは無く、無視できるだろうと予想できます。
しかし、水ではBB弾の半分程度の密度になるために無視はできません。これは、0.1gのBB弾が詰まった容器に0.2gのBB弾を打ち込むようなものです。なので、容器が揺らぐとBB弾もそれにつられて動きます。

なので、現状のBB弾の弾道計算理論では考慮していない効果が現れるはずです。

荒っぽいですが、水は全く揺らいでいないし、揺ぐこともないと考えましょう。
この近似の下で、計算と実際の軌道が程度一致するのかを確かめていきます。

実験上の問題

実験の動画は軌道がはっきりわかるくらいのハイスピードカメラが使われています。
しかし、動画の時間、位置スケールが分からない点と回転量が分からないため、大雑把な値しか見積もれません。
せめて、同じ回転量で大気中の軌道があれば何とかなったかもしれませんが、無い以上は仕方ありません。

よって、方針としては
位置はBB弾のサイズから見積もり、
時間の情報は軌道だけに注目することで無くし、
回転量は数値計算する時の回転量を変えていって、重なる軌道があるかどうか?

で考えで行きます。

スポンサーリンク


実験と理論の比較


理論上の式を数値計算で解きます。数値計算で用いる数式は弾道計算(BB弾)の理論で導いた、空気抵抗、揚力、回転の減衰を考慮した運動方程式です。

まず、マック堺様による実験動画[1]からサイズ、入射角などを見積もります。
スクリーンショットをとって、BB弾の位置を等間隔に撮っていったのが下の画像になります。

入射角の推定方法などは画像を見ていただければわかると思います。

画像から
入射角:51[度](→53[度]と推定)
入射位置からの到達距離:18~22[cm]
と分かりました。
重さ、初速は動画上で0.20g, 90m/sと紹介されていたのでそれを用います。

水中ではカーブするため、水中の2点間から導いた角度よりも大きくなることは確実です。
53度という値は数値計算で得られた軌道から53度位が良さそうだ、という推定をしました。

上記パラメータを利用し、
水の密度を\(1000[\mbox{kg/m}^3]\),
粘度を\(10^-3[\mbox{Pa}\cdot \mbox{s}]\)
として回転量を変えていった時の軌道は以下のようになります。

実験の状況から、横軸は拡大されて見えているはずなので、BB弾を打ち込んだ位置の到達距離は、推定した一番短い距離を利用しています。
図中の赤い大きな点が実験のBB弾の位置で、線と点で表した4つの線が数値計算です。4つの黒、青、緑、紫は回転量の違いを表し、順に100、150、200、250回転/秒で計算した結果です。
この回転数は妥当なものです。200回転の場合で大気中を想定し同じパラメータで計算しますと

のようになり、ホップが効かな過ぎず効き過ぎずの妥当な軌道になっています。

水の中を全く想定していない理論だとはいえ、良い結果だと思います。
特に、到達距離が数値計算でも20cm前後であることが定量的に再現できたのは大きな成果だと思います。また、水中でも揚力によって上昇することが定性的に再現できました。
揚力について、定量的に異なっています。実験では水面-最下点の半分以上上昇しているように見えますが、計算ではそこまで上昇しません。

空気が入っていることが大きな原因だと思いますが、はっきりとはわかりません。
先ほども述べた通り、水が揺らぐ効果は全く入っていないためです。
もしかするとエアガンを水中に沈め、空気が入らない状態で発射すればよいと思いますが、エアガンが壊れるので現実的ではありません。

銃口に膜でも作って、その膜を通してBB弾の衝撃をうまく伝わらせれば何とかなるかもしれませんが、今度は回転が無くなってしまいます。
ホップのゴムパッキンのすぐ後ろに膜でも作れば何とか…?

理論で回転の減衰を入れないと


回転の減衰を入れないとどうなるでしょうか?
大気中の場合はさほど重大な問題ではありませんでしたが、水中では以下のようになりました。200回転/秒で計算しています。

回転の減衰を入れない計算は図中、水色の線であり、明らかに物理が破綻しています。
回転の減衰を入れた計算は図中、緑色の線であり、直感と合い、良い振る舞いをしていることが分かります。

参考文献


[1]夏休み 自由研究? BB弾水に撃つと?検証してみた。 マック堺 エアガンレビュー
このページでのみ、動画のスクリーンショットに変更を加えた画像の使用許諾済み。

[2]水・空気の物性 密度 粘度 動粘度

[3]弾道計算(BB弾)の理論 -シキノート

BB弾の回転量について(実験との比較)

BB弾の回転量について,係数\(C_l\)は
\(
\begin{align}
C_l&\approx 0.12\\
\end{align}
\)

で近似できる、と実験データから求まりました。
実験データ(ハイスピードカメラによる回転量計測及び弾道データ)を提供してくださった
ガンジニア、石岡様(http://gungineer.matrix.jp/)に感謝いたします。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)←今ここ
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式


スポンサーリンク


対象とする問題


BB弾の運動方程式を求める事が出来ましたが、BB弾の回転量に関して計算と実験との間に4~5倍近い差が生じていました。この原因は、BB弾の回転量と発生する循環を一緒にしていたためです。式で表せば運動方程式は

・BB弾の運動方程式
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-C_l\left(S_p\right)\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
\)

・BB弾の回転角速度の減衰を記述する方程式
\(
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
\)

と書けます。今まで係数\(C_l\left(S_p\right)\)を1にしていました。
この係数を求める事が今回の目標です。

良く使われるパラメータにスピンパラメータ\(S_p\)なるものが存在するようで[1]、今回はこれに倣います。
これは回転球の表面での速度\(R\omega\)と飛翔速度\(v\)に関する無次元量で,
\(
S_p=R\omega/v
\)

と書かれます。
BB弾の場合、実験[2]より、おおよそ150[回転/秒, rps]ということが分かりました。
飛翔速度はおおよそ\(80[{\rm m/s}]\)と考えると、スピンパラメータの値の範囲はおおよそ
\(
\displaystyle S_p=\frac{3\times 10^{-3}{\rm [m]}\cdot 2\pi\cdot 150{\rm [/s]}}{80{\rm [m/s]}}\sim 0.0353\cdots
\)

という量になります。よって、\(C_l(S_p)\)を\(S_p=0\)回りでテーラー展開して、
\(
C_l(S_p)\sim c_0+c_1 S_p+c_2 S_p^2 +O(S_p^3)
\)

と考えられます(補足1)。

\(C_l(S_p)\)の導出


係数\(C_l(S_p)\)を以下のように求めます。

ある既知の回転量で射出されたBB弾が描く軌道

その既知の回転量を初期条件として運動方程式を解いた場合の軌道

の分散を最小にするパラメータをモンテカルロ法によって求めます。モンテカルロ法などと言っていますが、これは\(c_0,c_1\)に乱数を適当に振って良さげなものを求める、ということです。

ここでは実験から得られたある3つの回転量
case1:117回転/秒
case2:147回転/秒
case3:191回転/秒
の軌道との比較を行います。

実験のデータ点\(x_i\)と数値計算でのデータ点\(X_i\)との差の2乗の平方根をデータ点数\(N\)で割ったものを\(\Delta\)と定義します。

この\(\Delta\)が最も小さくなるパラメータの組み合わせが今求めたいものです。
縦軸に\(\Delta\)、横軸に各々のパラメータ\(c_0,c_1,c_2\)を取ったものをグラフにすると以下のようになります。
Delta_param2

各々のcase1,case2,case3で、\(\Delta\)を最小にするパラメータで軌道と終端速度の妥当性を見てみますとこうなります。

velocity_expthe5

速度が合わない原因としては恐らく、終端速度の計測機器による不具合でしょう。

なぜなら、初期速度を実験から得られたデータ値にした際、case3に対して同様にフィッティングするとばらつき具合を示す\(\Delta\)が1.9という値になってしまい、明らかに悪くなるからです。

実験のデータ点に+3[m/s]にするとおおよそ妥当なことが分かります。
velocity_expthe3
きっと計測機器の初期値がずれていただけでしょう。

\(\Delta\)の図から、最小地点3点での平均を取ればいいです。すると
\(
c_0\approx 0.18,\;\;\; c_1\approx -1.4 \;\; c_2=?
\)

と求める事が出来ました。\(c_2\)の項が求まりませんでしたので、\(c_1\)で打ち切ります。

考察


\(c_1\)で打ち切ることを考えると少し問題が生じます。今、\(c_1\)はどうやら負のようです。
と言うことは、これをこのまま高速回転する状況に使うと、回転が増えると逆に落ちやすくなる事を意味しています。\(S_p\)が上のように限定されている場合では精度よく求められますが、広い範囲ではちょっと上の級数展開式は使えません。

現実的な問題で考えましょう。
物理的に正しくなるように、係数\(C_l(S_p)\)は常に正の値でなければならないと仮定します。なので、上の級数展開はよろしくありません。
単純に定数として求めてしまいましょう。すなわち、
\(
\begin{align}
C_l(S_p)&\approx c_0
\end{align}
\)

とします。球自体の循環だけに依存するとして広い範囲で(とりあえず)使える定数を求めます。
Delta_param_const2
そうして実験データから\(\Delta\)を最小にする\(c_0\)を求めてやると、
\(
C_0=0.12
\)

と求まりました。
実際に\(C_l(S_p)=0.12\)で実験との比較をしてみますと、こうなります。
constcl012_c
定数として考えてもそんなに悪くないと思います。
case1,case3 を用いて、\(C_l(S_p)=0.12\)の妥当性を考えてみます。
\(C_l=0.12\)として、大体どの程度の回転量の誤差に収まっているかを示したのが以下の図です。
click1020pm20
数値計算上で求められた回転量は現実の回転量と10%程度の誤差に収まっていることが分かりました。

どれが良いかは一概には言えません。
回転量を知らなくていいなら\(C_l(S_p)=1\)で良いと思います。
回転数を大雑把にでも知りたければ\(C_l(S_p)=0.12\)を、
\(S_p\)が小さい、すなわちホップをあまりかけないで0.9J近くで飛ばす限定された状況では\(C_l(S_p)= 0.18-1.4S_p\)が良いと思います。

スポンサーリンク

結論


エアガンから射出された、BB弾の従う方程式は

・BB弾の運動方程式
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-C_l\left(S_p\right)\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
\)

・BB弾の回転角速度の減衰を記述する方程式
\(
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
\)

ここで、
\(
\begin{align}
&C_l(S_p)\sim 0.18-1.4S_p+O(S_p^2) \\
&S_p=R\omega/v
\end{align}
\)

もしくは
\(
C_l(S_p)\sim 0.12
\)

で近似できる。
という結論を得ました。

最後に、実験データ(ハイスピードカメラによる回転量計測[2]及び弾道データと速度データ)を提供してくださった
ガンジニア、石岡大樹様(http://gungineer.matrix.jp/)に心より感謝申し上げます。
以下の実験時の条件は石岡様が行っていただいた条件です。

実験[2]とデータ採取時の条件


エアガン:東京マルイ VSR-10 G-SPEC
BB弾:ギャロップ 0.3g
ハイスピードカメラ:PHOTRON MH4-10K、10000fps
データ点の採取方法:Stealth-Target ST17

室温、湿度、気圧の計測
1つのデータ点につき、30回計測

補足1


私の予測ですが、\(C_l\)の展開方法について述べておきます。
論文[1]によると、回転による揚力を私のようには書いておらず、
揚力を\(L\)と記述し、
\(
\displaystyle L=\frac{1}{2}C_l\rho U^2 A
\)

のように記述しています。ここで、\(\rho\)は空気の密度、\(U\)は球の速度、\(A\)は球の直径断面積です。
[1]では\(L\)を実験から直接求めて\(C_l\)を求めています。
この時に\(C_l\)をスピンパラメータ\(S_p\)で展開して実験データをフィッティングしているので私の計算と事情が違うのかもしれません。

[1]を頼らないことにしましょう。
物理的な予測から求めていきます。回転量依存だと考え、\(C_l(\omega)\)として考えます。
1.どんな回転量でも\(C_l\)は正である必要があります(負のマグヌス力は今回の場合発生しないと考える)。なので、
\(
C_l(\omega) > 0
\)

でなければなりません。
2.回転が非常に0に近い時、球体の循環と周りで発生する循環は一致するはずです。
なので、
\(
C_l(\omega\to 0)=1
\)

であることが予想されます。
3.回転が非常に強い場合、周辺の空気は限られるため発生する揚力は弱くなっていくはずです。
ただし、回転が無限大の時、揚力自体は発散すると予想できるので、係数\(C_L\)は揚力の発散を抑えてはなりません。なので、
\(
C_l(\omega\to\infty) = A\omega^{k},\;\;(A=\mbox{const},\;\; -1 \lt k )
\)

であるはずです。

物理的な考察から、係数\(C_l\)が満たすべき条件はこんなものでしょうか。
BB弾の回転量\(\omega\)が程々の場合は予想が出来ません。どこかにピークがあるのかもしれませんし、緩やかに減少していくだけかもしれません。
適当な関数を持ってくるのも手ですが、私は嫌います。
そんなことをやるのなら定数でいいだろうとして、
\(C_l=0.12\)と言う結論となりました。

参考文献

[1]鳴尾 丈司, 溝田 武人, 下園 仁志, “一様気流中で高速回転するゴルフボールの空気力測定と飛しょう実験” 日本機械学会論文集 B編 Vol. 70 (2004) No. 697 P 2371-2377
[2]石岡大樹、ホップアップの回転数測定 TM VSR-10 G-SPEC

上下幅が一定の時の軌道と、撃ち上げ撃ち下し時の軌道のずれについて

今まではゼロインする距離を固定した時の比較を行ってきましたが、
通常軌道を見て適正ホップだ、とか決定するかと思います。
この条件ならば上下幅を一定にしなければ適切な比較とは言えないと思ったのでその場合の結果です。

また、山のフィールドで
斜面の下から上に向かって打つとき、
斜面の上から下へ向かって打つとき
の着弾点は上にずれることを示します。

上下幅が一定の時の軌道


弾道の上下の幅を、0.25gのBB弾で0.9Jで射出した時の弾道の上下幅を基準にして重さを変えた時、弾道がどうなるかを計算すると以下の通りになります。
変更するBB弾の重さは販売されているものを基準とし、
0.12g, 0.20g, 0.21g, 0.23g, 0.25g, 0.28g, 0.30g
にしています。
グラフ中の点は0.05秒ごとの位置をプロットしており、図中の上の数本の線は水平距離のみを見た時のものです。
orbit_updown

静止画では
重さと軌道の関係_c
のようになります。

飛距離に対して到達速度は0.12g以外ではほとんど変わりません。
30m付近でおおよそ同じ時間に到達していることが分かります。
それ以降は重いほど早く到達し、かつ遠くまで飛んでいることが分かるかと思います。

やはり銃が許す限り、BB弾は重ければ重いほど良い、ということです。

撃ち上げ、撃ち下し時の軌道のずれについて


山のフィールドにおいて、斜面下から上、または上から下に向かって撃つ時、軌道はどう変わるのでしょうか。
答えは、照準よりも上に飛んでいきます。

以下の画像の状況を考えます。
角度_c

エネルギーは0.90J, 重さは0.25gを考えます。
下から上に向かって打つときの角度を正にとって\(+\theta\)
上から下に向かって打つときの角度を負にとって\(-\theta\)と表すことにします。
また、平坦な場所でゼロイン50mで最小の振れ幅になる軌道で比較します。

それぞれの角度\(\theta\)の向きで射出した後、撃った角度\(\theta\)で座標を回転させて主観的な視点で見た時にどうなるかを考えます。

すると、以下の結果を得ます。
撃ち上げ撃ち下し時_軌道のずれ_c
横軸は水平の距離ではなく、その角度で見た時の相対的な距離です。
\(\theta=0\)の時とそれぞれの角度の時を比較すると、角度をつける場合はほとんどの場合で上向きにずれることが分かります。
すなわち、角度がある時はシューティングレンジで合わせた照準よりも上に着弾するのです。

何故でしょうか。ちゃんと理由があります。
上向きにでも下向きにでも、角度があって撃つとき、BB弾に掛かる見かけ上の重力が減少します。
gravity_updown_c
なので角度がある時は重力が小さく感じるので、角度が無いところで調節した時よりも上の方に飛んでいくことになります。
なので、山のフィールドで斜面があるときは、若干下方向を、すなわち相手の膝や腰辺りを狙えばいいのです。

バレル内部でのBB弾の方程式

バレル内部でのBB弾の運動方程式です。

目的は、

BB弾は何秒間バレル内部に存在しているのか?
バレルが長いと減速になりうるのか?

を知ることです。
ただし、今は時間が無いため運動方程式を提示する、にとどめます。

弾道計算(BB弾)の理論
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式←今ここ


電動ガンの場合です。
バレル内部1
空気抵抗は初速0[m/s]から90[m/s]前後にまで加速されるわけですから、空気抵抗の粘性抵抗、慣性抵抗どちらの項も無視することはできないでしょう。

いくつか仮定があります。

  1. 気体の状態方程式\(PV=nRT\)が成立している(空気を理想気体として扱う, 補遺1)
  2. BB弾とバレルの間に空気の漏れは無い(補遺2)

ということです。
この条件下では、バレル内部の運動方程式は以下のように導くことができます。

ピストン-BB弾間の体積\(V\)はピストンの位置とBB弾の位置の2つに依存します。
なので、ピストンの位置による時間変化を\(x_0(t)\)、BB弾の位置を\(x(t)\)、バレルの断面積を\(S\)とすると
\(
V(t)=S\cdot (x(t)-x_0(t))
\)

と記述されます。
気体の状態方程式より、
\(
PV=nRT
\)

であり、右辺を定数\(c_0\)と置き、圧力\(P(t)\)をBB弾に掛かる力F(t)、すなわち\(P(t)=\frac{F(t)}{S}\)とすれば、
\(
\begin{align}
\frac{F(t)}{S}\cdot S\cdot (x(t)-x_0(t))=c_0 \\
F(t)=\frac{c_0}{x(t)-x_0(t)}
\end{align}
\)
と書けます。よって、空気抵抗と共に書けば、運動方程式は
\(
\displaystyle m\frac{d^2}{dt^2}x(t)=\left(\frac{c_0}{x(t)-x_0(t)}-P_0\cdot S\right)-6\pi \eta R \cdot v(t)-\frac{1}{2}C_d \rho \pi R^2 \cdot v^2(t)
\)

と導けます。ここで、\(P_0\)は大気圧、\(\eta\)は粘性率、\(R\)はBB弾の半径、\(C_d\)は抗力係数、\(\rho\)は空気の密度を意味します。
空気抵抗に関する詳細は球体の空気抵抗と係数をご覧ください。

あとはこの運動方程式を解いて、BB弾が出てくるバレル位置において適当な速度になるような物理的な意味を十分に持つ\(c_0, x_0(t)\)を入れて計算すればいいのです。

ここで1つ言えることは、外部と内部の圧力が一定になる最適なバレル長というものが存在する、ということです。
このバレル長を使えば、BB弾が射出されたときに、銃口で気圧差は無く、それによって生じる不確定な要素を最小にすることができるでしょう。

そのバレル長とは、式(4)の右辺が0になる時であり、計算によって求めることができます。
いくつになるかは実際求めてみないと想像できません。数cmであるのか数十mであるのか、はたまた30cmとかその位なのか・・・。

ちなみに実銃ではほぼ延々に火薬によって非常に多くのガスが生成され、加速されるためこのようなことは起こらないでしょう(あるかもしれませんが、とてつもなく長いバレルでなければ実現しないと思われます)。

補遺


補遺1
理想気体として扱える条件は、
低い気圧(分子の数が少なく、衝突等が無視できる)かつ高い温度(分子間力が分子の運動エネルギーに比べて無視できる)
です。
どうやら実在気体では10気圧以下、室温以上でこの条件は良く満たされ、理想気体とのずれは1%以内のようです[1]。
[1]を引用すると、

一般に,沸点の低い酸素・窒素・水素・ヘリウム等は,室温またはそれ以上の温度で10atm以下の圧力の場合,理想気体の値の1%以内で理想気体に近い性質を示す。

とあります。ピストンで空気が圧縮されたとき、BB弾とピストン間の圧力が10気圧以上にならなければ良い近似だと言えるでしょう。

補遺2
BB弾とバレル内側との間における空気の漏れはどの程度あるのかはわかりません。

バレルの内径:6.08mm (東京マルイのものはこれらしいですが、信頼できるページが見当たりませんでした。)
BB弾の直径:5.95mm [2]

完全に理想気体として、時間を考慮しないで扱うのであれば、少しでも空間が空いていればどれだけ早く押し込もうとも気圧は一定になろうとするため、BB弾に掛かる力はゼロです。
なので時間を考慮した方程式が必要でしょう。
連続の式
\(
\displaystyle \frac{\partial \rho}{\partial t}+\nabla \cdot \vec{j}=0
\)

を用いれば時間と共に漏れ出る空気量を知ることが出来そうです。
ただ、まずは簡単のため、空気が漏れない場合をまず考えましょう。

簡単に空気の漏れを表現するには、あたかもピストン-BB弾間の体積\(V\)が時間と共に増加していく、と考えることである程度の効果を入れることができると考えられます。

実際に式を解いてから、空気の漏れが起こる場合と起こらない場合でどれほどの差が出てくるのか?を考えることでこの効果を入れるべきか入れないべきかを考えるのがいいでしょう。

良く空気の漏れがーと言われているので本当かどうか理論的に検証するべきでしょう。

参考文献


[1] 実在気体 -第2節 気体の状態方程式
[2]PERFECT HIT -TOKYO MARUI

『集弾性アップへの道』 BB弾とバレル内部に隙間があることが写真で確認できます。

球体の空気抵抗と係数

まとめ


半径\(R\)の完全な球体に働く空気抵抗力\(F_d\)は
\(
\displaystyle F_d = 6\pi R \eta v + \frac{1}{2} C_d \rho \pi R^2 v^2
\)

です。方向も含めるのであれば
\(
\displaystyle \vec{F}_d = -6\pi R \eta |\vec{v}| \frac{\vec{v}}{|\vec{v}|} – \frac{1}{2} C_d \rho \pi R^2 |\vec{v}|^2 \frac{\vec{v}}{|\vec{v}|}
\)

となります。

ここで、

  • \(\eta\ \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) : 空気の粘性率(Viscosity) ※動粘度ではないです。
  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 空気の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 球体の速度の大きさ
  • \(R\ \mathrm{[m]}\) : 球の半径

を表します。
\(C_d\)は抗力係数(Drag Coefficient)と呼ばれ、厳密にはレイノルズ数\(R_e\)の関数\(C_d=C_d(R_e)\)となります。
レイノルズ数\(R_e\)とは無次元量で、
\(\displaystyle R_e=\frac{v \rho L}{\eta}\)
という量です。ここで\(L\)は、物体の大きさを表し、半径\(R\)の球の場合は直径\(L=2R\)に相当する量です。

球の場合の抗力係数\(C_d(R_e)\)をフィッティングしたものは、論文[2]より、
\(
\displaystyle \small C_d=\frac{24}{R_e}+\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
+ \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}} + \left(\frac{R_e^{0.80}}{461000}\right)
\)

と表されます。
Cdグラフ

20℃で1気圧(101325[Pa])、湿度0%の空気中では、
\(
\begin{align}
\eta &= 18.2\times 10^{-6}\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]} \\
\rho &= 1.205\mathrm{[kg \cdot m^{-3}]}
\end{align}
\)

という値になります。

温度、圧力、湿度を考慮する場合

温度\(T[\mathrm{K}]\)、圧力\(P[\mathrm{Pa}]\)の時,
空気の粘性率\(\eta\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\)は、
\(\displaystyle \eta\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}=1.487 \times 10^{-6}\cdot \left(\frac{T\mathrm{[K]}^{1.5}}{T\mathrm{[K]}+117}\right)\)
であり、
空気の密度\(\rho\ \mathrm{[kg \cdot m^{-3}]}\)は、
\(\displaystyle \rho\ \mathrm{[kg \cdot m^{-3}]}=0.0034856447\cdot \frac{P\mathrm{[Pa]}}{T\mathrm{[K]}-0.670}\)
です。

この2つの式を[1]が正しいと思って確かめた温度の範囲-10℃~40℃で大体3桁くらい一致します。
この式は下に説明してある空気の場合を代入し、求めたものです。

2016/03/13 追)
湿度がある場合、空気の密度\(\rho_{w}\mathrm{[kg \cdot m^{-3}]}\)は湿度0%の時と比べて減少します。
温度\(t\mathrm{[{}^{\circ}C]}\)、湿度M(湿度60%の時、\(M=0.60\))、気圧\(P\mathrm{[Pa]}\)の場合、空気の密度\(\rho_{w}\mathrm{[kg \cdot m^{-3}]}\)は
\(
\displaystyle \rho_{w}\mathrm{[kg \cdot m^{-3}]}=
\left(0.0034856447\cdot\frac{P\mathrm{[Pa]}}{t\mathrm{[{}^{\circ}C]}+272.48}\right)\cdot
\left(1-M\cdot 0.378\cdot\frac{100\cdot 6.1078\times 10^{\frac{7.5\cdot t\mathrm{[{}^{\circ}C]}}{t\mathrm{[{}^{\circ}C]}+237.3}}}{P\mathrm{[Pa]}}\right)
\)

と記述されます。
ここで用いた関係式は
[5]に圧力pの水蒸気を含んだ空気の密度の関係式

[6]、[7]の飽和水蒸気圧を求めるTetensの式
を用いています。また、atmによる単位系では[8]に記述があります。


スポンサーリンク



スポンサーリンク

風の抵抗力の詳細


これと同じことは弾道計算(BB弾)の理論に記述してあります。
その中の風の抵抗力を部分を抜き出したものになります。

風(流体)の抵抗力の大きさ\(F_d\)は一般には次元解析により
\(
\displaystyle F_d=\frac{\eta^2}{\rho}\sum_{n=1}^{\infty}K_n \left(\frac{v \rho L}{\eta}\right)^n
\)

として与えられます[3]。ここで

  • \(\eta\ \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) : 流体の粘性率(Viscosity) ※動粘度ではないです。
  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 流体の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 物体の速度の大きさ
  • \(L\ \mathrm{[m]}\) : 物体の大きさ(半径\(R\)の球の場合、\(L\)は直径\(L=2R\)に相当する量)
  • \(\frac{v \rho L}{\eta}\) : レイノルズ数(\(R_e\),Reynolds number)と呼ばれる無次元量

です。次元解析から具体的な定数の値\(K_n\)を求めることはできません。

速度の3乗以上に比例する項というのは見かけません。
おそらく、物体の動きを記述するために必要なパラメータがこの2つを考慮すれば含まれるため、\(n=3\)以上の項は冗長になるからでしょう。つまり、第一項には空気の密度\(\rho\)が含まれておらず、第二項まで入れることで\(\eta,\rho\)があらわに含まれるようになるからだと思います。

\(n=3\)以上の高次の効果は、抗力係数\(C_d\)に押し込めていると思います。

上の式の\(n=1,n=2\)をそのまま書き下すと
\(F_d=K_1\eta L v +K_2\rho L^2 v^2\)
となりますが、慣例として同じ次元を持つように式を変形させて
\(F_d=C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right)\)
の形で記述されます。ここでSは速度に垂直な、物体の断面積です。

\(C_1\)は物体が完全な球で、小さく、その速度がゆっくりである場合(レイノルズ数が2以下くらい)、流体の運動を記述する方程式、
ストークス方程式(Stokes’ law、ナビエ・ストークス方程式のレイノルズ数が小さい極限の方程式)
から導くことができて\(C_1=3\pi\)であることが理論的に導けます。

\(C_d\)は抗力係数(Drag Coefficient)と呼ばれ、
これはレイノルズ数の関数となり、
\(C_d=C_d(R_e)\)
となります。
この抗力係数は非常に厄介らしいです。物体の形状に依存したりするらしく、球の場合でも、どうやら実験的に求められているようです。
Wikipediaの抗力係数のページに若干の記述があります。
球の場合の抗力係数の関数\(C_d(R_e)\)をフィッティングした関数は論文[2]より、
\(
\displaystyle \small C_d=\frac{24}{R_e}+\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
+ \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}} + \left(\frac{R_e^{0.80}}{461000}\right)
\)

として書けるそうです。実際にプロットしてみるとこんな感じです。
Cdグラフ
確かに抗力係数のグラフと同じようになります。

よって半径Rの球の場合、風の抵抗力\(F_d\)の大きさは結局、
\(
\begin{align}
F_d &= C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right) \\
&= 3\pi \cdot 2R\cdot \eta v +\frac{1}{2} C_d \rho \pi R^2 v^2 \\
&= 6\pi R \eta v + \frac{1}{2} C_d \rho \pi R^2 v^2
\end{align}
\)
として与えられます。

各定数の値
名称 捕捉
空気の粘性率(Viscosity)\(\eta\) \(\small 18.2\times 10^{-6} \\ \small \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) 空気中20度において。また、粘性率は数十気圧から数気圧の間圧力には依存しません。
が、温度に依存し、ある温度\(T_1[K]\)での粘性率\(\eta_1\)が分かっているならば、
温度\(T_2[K]\)での粘性率\(\eta_2\)はサザランドの定数Cを用いて、
\(
\eta_2=\eta_1 \left(\frac{T_1[K]+C}{T_2[K]+C}\right)\left(\frac{T_2[K]}{T_1[K]}\right)^{3/2}
\)
と記述されます。空気の場合、C=117です[4]。
乾燥した空気の密度\(\rho\) \(\small 1.205\mathrm{[kg \cdot m^{-3}]}\) 乾燥した空気の密度は温度と圧力に依存します。温度\(t[℃]\)、圧力をH[torr]とすると
\(
\rho=\frac{1.293}{1+0.00367t[℃]}\cdot \frac{H[torr]}{760}
\)
です[5]。

参考文献


[1]水・空気の物性 密度 粘度 動粘度
[2]Faith A. Morrison, “Data Correlation for Drag Coefficient for Sphere,” Department of Chemical Engineering, Michigan Technological University, Houghton, MI
[3]伊東敏雄著『な~るほど!の力学』学術図書出版社 (1994) p.220
[4]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.378
[5]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.376
[6]相対湿度の月別平年値 -理科年表オフィシャルサイト
[7]飽和水蒸気量 -wikipedia
[8] 空気 -wikipedia

ちなみに、
\(1 [\mathrm{torr}] = 133.322368 [\mathrm{Pa}]\\ 1[\mathrm{Pa}]=1[\mathrm{N\cdot m^{-2}}]=1[\mathrm{kg \cdot m^{-1}\cdot s^{-2}}]\)
という関係です。

弾道計算(BB弾)の結果2、違う重さでの比較

本稿の主要な結果は、様々なパラメータでの、BB弾の軌道の詳細なデータです。

この結果に付随して、
0.25g0.8Jの軌道は
0.20g1.2Jの軌道と同じであることが分かりました。

上下振れ幅が最小になる軌道における、重さとエネルギー、ゼロイン位置を様々にとった時のそれぞれのデータです。
弾道計算を数値計算によって行い、結果を考察することが本稿の目的となります。
※本稿では規定のエネルギーを超える場合のデータがありますが、このデータは全て数値シミュレーションであるため、こういったエアガンを持っている訳ではありません。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ←今ここ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算


様々なパラメータの最適な軌道


ここで載せるデータは、以下の組み合わせです。

BB弾の重さ: 0.20g , 0.25g
エネルギー:0.60J 0.65J 0.70J 0.75J 0.80J 0.85J 0.90J
ゼロイン位置:25m 30m 35m 40m 45m 50m 55m 60m(0.25gのみ)
各々でどのような軌道を描くかの計算結果を載せます。
破線はホップ無しの軌道に対応しています。

↓これは低画質です。高画質版はこのリンク(4.3MB)を踏んでください。
弾道詳細データ_高画質2 - コピー_c


最適な軌道に合わせるためには下向きに初速を与えます。
その時、下向きに何度傾けて撃てばいいか、のデータはこちらです。

ゼロインと下向き角度


前提


まず、規定速度うんぬんは別にしまして、事実を述べます。

  1. BB弾は重いほど弾はまっすぐ飛ぶ
  2. 射出速度が早いほど弾はまっすぐ飛ぶ
  3. ホップを掛けるほど弾は上下に揺れる

これらは紛れもない事実です。
よって、BB弾を飛ばす最高の条件とは、
「射出速度を早くし、ホップは余りかけないで重いBB弾を使うこと」
となります。
だからこそ、射出速度を速めようとして規定速度の話になります。

最適な軌道とは?


BB弾を飛ばすうえで最適な軌道とはどんな軌道でしょう?
それはホップによる上下方向の揺れを最小限に抑える軌道です。
この最適な軌道とは、BB弾の重さとエネルギーとゼロイン位置を決めた時の理想的な軌道、ということです。

ゼロイン位置と上下の振れ幅とは何かは、下の画像をご覧ください。
0.25,0.8_compressed

重さとエネルギーを決めても、どの角度で射出すればいいか、回転数はいくつか、など他のパラメータが残ります。
それを決めるため、ホップによる上下方向の揺れを最小限に抑える軌道を定めます。

この軌道を ”最適な軌道” と呼ぶことにします。

重さの違いによる軌道の具体的な影響


今度は、ゼロイン位置を固定します。
そして、BB弾の重さとエネルギーを変化させたとき、どのような軌道をたどるか見てみましょう。
zeroin50_2_cc

この結果から分かることは2点あります。

1つ目は、エネルギーを上げていくと軌道の変化が小さくなる
2つ目は、重いBB弾の軌道ほどまっすぐ飛ぶ

ということです。

エネルギーを上げると軌道の変化が小さくなる

0.20gの場合でゼロインを50mに合わせ、エネルギーを0.1Jずつ増やしていったとき、上下振れ幅がどのくらい減少していくのか見てみましょう。

0.60J→0.70J … 14.2cm
0.70J→0.80J … 10.4cm
0.80J→0.90J … 8.0cm
0.90J→1.00J … 6.2cm
1.00J→1.10J … 5.0cm
1.10J→1.20J … 4.2cm
1.20J→1.30J … 3.6cm
1.30J→1.40J … 3.0cm

となります。エネルギーを変えるよりも重さを変えるべきです。
それだけで軌道は大きく変わります。

例えば、0.80Jで射出できるエアガンで、使うBB弾を0.25gで射出すると、この時、軌道は0.20g1.2Jの軌道にほぼ一致します。
m020J12m025J08_2_c

まぁ、軌道自体は一致するのですが、残念ながら到達時間は一致しません。
約0.1秒(まばたき程度)、0.20g1.2Jのほうが早くなります。
zeroin50_2_time_c

追記)
0.25gで0.90Jは、0.20gの1.35Jの軌道におおよそ相当します。
0.30gで0.80Jは、0.20gの1.60J, 0.25gの1.05Jの軌道におおよそ相当します。
0.30gで0.90Jは、0.20gの1.80J, 0.25gの1.20Jの軌道におおよそ相当します。

0.25gで0.8Jでうまく合わせられたら50m飛ばしても上下の振れ幅40cmです。
上に20cm、下に20cmしかずれません。
また、エネルギーを上げても空気抵抗の強さは速度の2乗に比例して強くなるため、それに見合うようなの飛距離の伸びは得られません。


補足


「流速チューン」というものがあるそうですね。
流速チューンの問題点と改善案について -週休5日

流速チューン比較テスト その1「初速変化」 -Gunsmithバトン
にあるように。

本稿ではこれらのことについては全く触れていません。
というのも、銃口から射出された直後の初速と回転数の2つだけが軌道に影響するためであり、発射するまでの過程に何があったか?などどうでもいいのです。
数値計算の観点から言いますと、通常の軌道よりも「まっすぐ遠くまで飛ぶ」場合、初速が上がっていること以外には考えられません
BB弾を飛んでいる最中で加速させる要因など無いのです。


※1
なぜBB弾を使う限り初速を稼ぐことが無駄なのか?
これは空気抵抗が関係します。
空気抵抗の大きさは弾道計算(BB弾)の理論で書いたように、BB弾の半径\(R\)と速度\(v\)にのみ依存します。
空気抵抗力を半径\(R\)と速度\(v\)の関数として\(F_d(R,v)\)、質量\(m\)とすると、その運動方程式は
\(
m\frac{d^2\vec{r}}{dt^2}=F_d(R,v)
\)
であり、両辺をmで割れば
\(
\frac{d^2\vec{r}}{dt^2}=F_d(R,v)/m
\)
となります。BB弾の半径や、温度、空気は変えられないので、可変なパラメータは質量と速度のみです。
そして右辺だけに注目すれば、質量が大きいほど空気抵抗力があたかも小さくなるのです。
すなわち、空気抵抗力を減らそうとするにはBB弾の重さを変えるほかない。
また、初速を変えても空気抵抗は速度の2乗で効いてくるため、皮肉なことに速度が早いほど空気抵抗力が強くなっていくのです。
よって、高いエネルギーでは初速を上げてもそれに見合っただけの結果は得ることはできないのです。

弾道計算(BB弾)のコード

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
弾道計算(BB弾)のfortran90での計算コードです。


2016/03/18 以降のコードに関して
コードは改変、複製、配布、営利、非営利目的に使用していただいて構いません。
ただし著作権を放棄している訳ではありません。
また、このコードを利用した場合、クレジットの表記をお願いします。


Linuxで、OSがUbuntuやLinuxMintの場合はfortranコンパイラは

sudo apt-get install gfortran

でインストールできます。
intel®fortranコンパイラとgfortranコンパイラで正常に動くことは確かめています。

もろもろのファイルをmain.f90に記述した後、
gfortranコンパイラなら

gfortran main.f90
./a.out

intel®fortranコンパイラなら

ifort main.f90
./a.out

で実行、ファイル”orbit.txt”が生成されます。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)←今ここ
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算

できること


・BB弾の弾道計算(重力、慣性抵抗、粘性抵抗、ホップ有、回転の減衰を考慮)ができます。
式で表すと
BB弾の運動方程式(式1)
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
\)

とBB弾の回転角速度の減衰を記述する方程式(式2)
\(
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
\)
の2つを考慮したものになります。

・ゼロインを決めた時の最適な軌道
・上下方向が最小の振れ幅になる軌道
も出力することができます。

※BB弾の回転の減衰を考えるか考えないかは、inputパラメータ”omgdecay”により制御できます。
回転の減衰に疑問がある場合はこれを利用してください。

注意事項


・とんでもないパラメータでは計算できない場合がありますので、そこら辺はご了承ください。
・僕のプログラミング知識ではこれが限界です。
僕はこのプログラムを使用したことによって生じた、いかなる場合でも一切の責任を負いません。

主なルーチンと説明

計算理論


式1の時間発展はルンゲ-クッタ-フェールベルグ法による適応刻み幅制御を行っています。デフォルトでは8桁の精度で(※1)一致するようになっています。モジュールGBL内のtol=1d-8を変更することで精度を制御することができます。
式2の右辺の積分は、デフォルトでは15次のガウス=クロンロッド求積法を2回用いています。今回の問題に対してはあまり良い方法とは言えません。ですが、5-6桁はロンバーグ積分法による結果と一致するので、実用上問題ないと判断し、計算速度を優先しました。
式2の回転減衰の効果はルンゲクッタ法を行う際の連立方程式に組み込んでいません。これだけはオイラー法で行っています。なので、1ステップ目の刻み幅をどうするかで、本来はずれないはずですが、変更するとそれが伝搬してゼロイン位置がずれます。デフォルトでは全ての1ステップ目の刻み幅は統一してあるので上記問題は起こりません。本当に厳密にやるならば、コードを変更して減衰を連立方程式内に組み込んでください。
2016/03/13更新
コードを変更して減衰を連立方程式内に組み込みました。計算は若干遅くなりますが、上記問題はもう起こり得ません。

※1 相対誤差もしくは絶対誤差で、です。この理由は、相対誤差に統一してしまうと値が非常に小さいとき、機械誤差によっていつまでたっても収束しないためです。なので、元になる値が1より小さいときは絶対誤差を、1より大きいときは相対誤差を取るようにしています。

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
2017/10/07更新

inputファイルとコード


inputというファイル名のファイルを作り、中に

&input
m=0.25d-3,          ! Weight of bullet [kg]
a=3d-3,             ! Radious of bullet [m]
energy=0.90d0,      ! Energy of bullet [J]
g=9.80665d0,        ! Gravity [m s^-2]
temperature=20d0,   ! Temperature [degree of Celsius]
pressure=101325d0,  ! Atmospheric pressure [Pa]
moisture=0.6d0,     ! Value of moisture 0~1
omgdecay="yes",     ! Introduce omega decay. "yes" or "no"
search_zeroin="yes", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0,       !   +-- if yes, set zeroin length [m]
search_updown="yes",! Search good vz and omgy. "yes" or "no"
updownz=0.20d0,     !   +-- if yes, set up-down henght [m]
theta=-1.7d0,
rx=0.d0,
ry=0.d0,  
rz=1.0d0,
vy=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi
omgy=-340d0, !*2pi
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0,       ! Output interval [s]
ik=0, ! ik=0 exact, ik=1 approximation
&end

と記述してください。このinputファイルの例では、

重さ(m)0.25g
エネルギー(energy)0.90J
重力加速度(g)9.80665m/(s^2)
温度(temperature)20度
大気圧(pressure)101325Pa
湿度(moisture)60%
回転角速度の減衰(omgdecay)有り
ゼロインによる最適な弾道の探索(search_zeroin)しない
↑がyesなら、ゼロインの位置(zeroinx)50m
上下幅による最適な弾道の探索(search)する
↑がyesなら、上下幅(updownz)0.1m
銃口の位置(rz)1m
y方向の速度(vy)0m/s
射出角度(theta)-0.5度
風速(ux=uy=uz=0)0m/s
ホップ(omgy)-210回転/s
出力ファイル名の接頭語(outputfile)orbit
出力の時間間隔(stept)0.01秒ごとに出力
回転の減衰を表す積分の近似方法(ik=0 : 積分の直接計算、ik=1 : 積分を近似)

となっています。
回転の減衰を考えたくない時は、omgdecay=”no”としてください。

そして下のコードを展開し、コピー&ペーストして一つのファイル(例えばmain.f90という名前)に入れてください。
(カーソルをmodule GBL…のmに合わせた後、一番下までスクロールしてshift+クリックで楽に選択できます。)
個人的にファイルをダウンロードして…は嫌いなのでこういう形で公開しています。

2016/03/13 更新
2016/03/18 更新
2016/03/20 更新
2016/07/08 更新
2016/07/23 更新
2016/07/25 更新
2016/08/06 更新
ソースコードは1523行あります。


使い方


・fortranコンパイラ(ifortやgfortran)
・gnuplot
を想定しています。

2017/02/04追)
Bernadotte66さんの報告より、windows上でgfortranコンパイラを用いる際、gfortranコンパイラがver5.1以前のものである場合エラーが出てしまうようです。ver6.2以上のコンパイラを推奨します。

Linuxでは、gfortranコンパイラver4.8.4で問題ないことは確かめています。

それでコンパイルして実行してください。
一連の流れとしては、以下の通りです。

$ gfortran main.f90
$ ./a.out
&INPUT
 M=  2.5000000000000001E-004,
 A=  3.0000000000000001E-003,
 ENERGY= 0.90000000000000002     ,
 G=  9.8066499999999994     ,
 TEMPERATURE=  20.000000000000000     ,
 PRESSURE=  101325.00000000000     ,
 MOISTURE= 0.59999999999999998     ,
 OMGDECAY="yes                                             ",
 SEARCH_ZEROIN="yes                                             ",
 ZEROINX=  50.000000000000000     ,
 SEARCH_UPDOWN="yes                                             ",
 UPDOWNZ= 0.20000000000000001     ,
 THETA=-0.50000000000000000     ,
 RX=  0.0000000000000000     ,
 RY=  0.0000000000000000     ,
 RZ=  1.0000000000000000     ,
 VY=  0.0000000000000000     ,
 UX=  0.0000000000000000     ,
 UY=  0.0000000000000000     ,
 UZ=  0.0000000000000000     ,
 OMGX=  0.0000000000000000     ,
 OMGY= -210.00000000000000     ,
 OMGZ=  0.0000000000000000     ,
 STEPT=  1.0000000000000000E-002,
 OUTPUTFILE="orbit                                           ",
 IK=          0,
 /
Reach ground at  0.980E+00
===============================
 search vz and wy by zeroin -->   50.00000
 ---- false_position
 Sequence 1 of 25 : omgy/2pi  -767.96000000 : zeroin   100.00000000
 Sequence 2 of 25 : omgy/2pi  -735.92000000 : zeroin   100.00000000
 Sequence 3 of 25 : omgy/2pi  -703.88000000 : zeroin   100.00000000
 Sequence 4 of 25 : omgy/2pi  -671.84000000 : zeroin   100.00000000
 Sequence 5 of 25 : omgy/2pi  -639.80000000 : zeroin   100.00000000
 Sequence 6 of 25 : omgy/2pi  -607.76000000 : zeroin   100.00000000
 Sequence 7 of 25 : omgy/2pi  -575.72000000 : zeroin    17.22935680
 Sequence 8 of 25 : omgy/2pi  -543.68000000 : zeroin    15.16813299
 Sequence 9 of 25 : omgy/2pi  -511.64000000 : zeroin    12.90696782
 Sequence 10 of 25 : omgy/2pi  -479.60000000 : zeroin    10.42356698
 Sequence 11 of 25 : omgy/2pi  -447.56000000 : zeroin     7.68923851
 Sequence 12 of 25 : omgy/2pi  -415.52000000 : zeroin     4.66717441
 Sequence 13 of 25 : omgy/2pi  -383.48000000 : zeroin     1.30910154
 Sequence 14 of 25 : omgy/2pi  -351.44000000 : zeroin    -2.44916818
  +--- Root find, go to conversion phase --+  
Conversion =>       0.00033361%
Conversion =>       0.00902931%
Conversion =>       0.24389215%
Conversion =>       6.77148580%
Conversion =>     104.78145034%
 ---- false_position
 Sequence 1 of 1 : omgy/2pi  -366.89192543 : zeroin    -0.57095843
  +--- Root find, go to conversion phase --+  
Conversion =>      0.00075051%
Conversion =>      0.08150935%
Conversion =>      8.77790138%
Conversion =>    579.19923386%
  -371.80132940800553       -2.5094934272766087    
Reach ground at  0.170E+01
===============================
 search vz and wy by up-down height -->    0.20000
 ---- false_position
 Sequence 1 of 25 : omgy/2pi  -959.96000000 : zeroin    19.21215678
 Sequence 2 of 25 : omgy/2pi  -919.92000000 : zeroin    17.82618788
 Sequence 3 of 25 : omgy/2pi  -879.88000000 : zeroin    16.40387340
 Sequence 4 of 25 : omgy/2pi  -839.84000000 : zeroin    14.95914779
 Sequence 5 of 25 : omgy/2pi  -799.80000000 : zeroin     9.80000000
 Sequence 6 of 25 : omgy/2pi  -759.76000000 : zeroin     9.80000000
 Sequence 7 of 25 : omgy/2pi  -719.72000000 : zeroin     9.80000000
 Sequence 8 of 25 : omgy/2pi  -679.68000000 : zeroin     9.80000000
 Sequence 9 of 25 : omgy/2pi  -639.64000000 : zeroin     9.80000000
 Sequence 10 of 25 : omgy/2pi  -599.60000000 : zeroin     9.80000000
 Sequence 11 of 25 : omgy/2pi  -559.56000000 : zeroin     1.45287636
 Sequence 12 of 25 : omgy/2pi  -519.52000000 : zeroin     1.10739833
 Sequence 13 of 25 : omgy/2pi  -479.48000000 : zeroin     0.80320434
 Sequence 14 of 25 : omgy/2pi  -439.44000000 : zeroin     0.54089117
 Sequence 15 of 25 : omgy/2pi  -399.40000000 : zeroin     0.32062896
 Sequence 16 of 25 : omgy/2pi  -359.36000000 : zeroin     0.14208536
 Sequence 17 of 25 : omgy/2pi  -319.32000000 : zeroin     0.00427472
 Sequence 18 of 25 : omgy/2pi  -279.28000000 : zeroin    -0.09464709
  +--- Root find, go to conversion phase --+  
Conversion =>       0.00183551%
Conversion =>       0.00989716%
Conversion =>       0.05297720%
Conversion =>       0.28399231%
Conversion =>       1.52644398%
Conversion =>       8.26417579%
Conversion =>      39.67268127%
Conversion =>     179.58494188%
 ---- false_position
 Sequence 1 of 1 : omgy/2pi  -312.86039405 : zeroin    -0.01405310
  +--- Root find, go to conversion phase --+  
Conversion =>      0.00065868%
Conversion =>      0.03115036%
Conversion =>      1.49063057%
Conversion =>    102.25144566%
  -317.78761751909860       -1.5473682498931858    
Reach ground at  0.149E+01
     2.681[CPU sec]
$ gnuplot

    G N U P L O T
    Version 4.6 patchlevel 4    last modified 2013-10-02
    Build System: Linux x86_64

    Copyright (C) 1986-1993, 1998, 2004, 2007-2013
    Thomas Williams, Colin Kelley and many others

    gnuplot home:     http://www.gnuplot.info
    faq, bugs, etc:   type "help FAQ"
    immediate help:   type "help"  (plot window: hit 'h')

Terminal type set to 'wxt'
gnuplot> plot "orbit.txt" u 2:4 w l
gnuplot> replot "orbit_h.txt" u 2:4 w l
gnuplot> replot "orbit_opt.txt" u 2:4 w l

gnuplot上での出力が
orbitGK
のようになったら成功です。

出力ファイルは最大で3つ生成(inputパラメータ “search_??????” )に依存します。
orbit.txtは、inputによって与えられたパラメータでそのまま計算したもの
orbit_opt.txtはそのエネルギー、重さで上下振れ幅が最小でゼロイン位置がinputのzeroinxのもの
orbit_h.txtはそのエネルギー、重さで上下振れ幅が最小で上下振れ幅がinputのupdownzのもの
となります。

出力ファイルはそれぞれ14列あり、
時刻t[秒],位置x,y,z[m]:x,y,z方向の速度[m/s]:BB弾の回転数x,y,z[回転/s]:x,y,z方向の風速:エネルギー[J]
になっています。
出力ファイルの中身は下のように書いてあるので、見てください。

# m:    0.25000E-3[kg]
# g:       9.80665[m s^{-2}]
# a:    0.30000E-2[m]
# temperature:      20.00000[degree Celsius]
# moisture:       0.60000[no-dimension]
# pressure:  101325.00000[Pa]
# eta:  0.18197E-4[kg m^{-1} s^{-1}]
# rho:     1.20123[kg m^{-3}]
# tol:  0.10000E-7
# eps:  0.10000E-8
#         t[s]         x[m]         y[m]         z[m]      vx[m/s]      vy[m/s]      vz[m/s]    wx[rot/s]    wy[rot/s]    wz[rot/s]   windx[m/s]   windy[m/s]   windz[m/s]    Energy[J]
     0.000000     0.000000     0.000000     1.000000    84.849583     0.000000    -0.740471     0.000000  -210.000000     0.000000     0.000000     0.000000     0.000000     0.900000
     0.010000     0.838863     0.000000     0.992911    82.937695     0.000000    -0.678265     0.000000  -209.261019     0.000000     0.000000     0.000000     0.000000     0.859890
    ...

実銃のパラメータでの計算


実銃は特に想定していませんが、計算ができると言えばできます。
実銃と比較して考慮されていない点は3つあり、

  1. 滑らかな球として考えているので抗力係数\(C_d\)が異なる
  2. 銃弾の歳差運動が考慮されていない(球ではなく細長いので)
  3. ライフリングによる回転が異なる
  4. コリオリ力を考慮に入れていない

1番目はかなり重大です。
私のプログラムを使うと1000メートル地点でのエネルギーの値が2倍ほど大きく出ます。
ライフリングによる横回転だけなら問題なく考慮できます。バックスピンで入っていたものを変えればいいだけです。
ただし、風が無い、またはコリオリ力を考慮しない状況では計算上は横回転があっても無くても変わりません。
風とコリオリ力がある時にライフリングによる弾道のずれが影響します。
2番目はどんな力が働いたりするか私には見当がつきません。専門の方じゃないと分かるまでに時間がかかると思います。
コリオリ力はプログラムで入れればいいです。これはそこまで辛くないでしょう。
ライフリングの効果は入れられるには入れられますが、ちょっと変な結果になるので使わないほうが良いでしょう。
なのでこの効果は入れていません。
パッと思いつく限り、違う点は以上です。

よって、本プログラムによる実銃を想定した計算とは、
風が無く、赤道上で実銃と同じ口径の滑らかな球が実銃のエネルギーで射出されたときの弾道
となります。

ちなみに、.338ラプアマグナム弾の場合を想定したinputパラメータはこちら
角度を使って計算しているのでそこのとこのプログラムは各自直してください。

&input
! For .338 lapua magnum( but treated as sphire)
m=16.2d-3,          ! Weight of bullet [kg]
a=4.29d-3,      ! Radious of bullet [m]
energy=6562.16d0,      ! Energy of bullet [J]
g=9.80665d0,        ! Acceleration of Gravity [m s^-2]
temperature=20d0,   ! Temperature [degree of Celsius]
pressure=101325d0,  ! Atmospheric pressure [Pa]
moisture=0.6d0,     ! Value of moisture 0~1
omgdecay="no",     ! Introduce omega decay. "yes" or "no"
search_zeroin="no", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0,       !   +-- if yes, set zeroin length [m]
search_updown="no",! Search good vz and omgy. "yes" or "no"
updownz=0.10d0,     !   +-- if yes, set up-down henght [m]
theta=0.6d0,        ! default, theta is neglected. change program main.
rx=0.d0,
ry=0.d0,  
rz=1.0d0,
vy=0.d0,
vz=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi [omgx*2*pi rotate/s]
omgy=0.d0, !*2pi [omgy*2*pi rotate/s]
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0,       ! Output interval [s]
&end

上記パラメータで計算するとこんな軌道が得られます(点は0.05秒ごとの軌跡)
lapua338_c

.338ラプアマグナム弾のデータは↓
.338 ラプアマグナム弾 / .338 Lapua Magnum -MEDIAGUN DATABASE
を参照しました。

弾道計算(BB弾)の結果

弾道計算(BB弾)の数値計算結果を中心にまとめます。


2016/07/09 追)
このページに置いて、「○○回転」の「回転」と書いてあるところは実際の回転量の約0.12倍です。
実際の回転量はこのページに書いてある回転量を8.3倍して下さい。

(例えば”20回転/秒”とあったら、これは”166回転/秒”と言うように読み変えてください)

目次

  1. 目的
  2. 結論
  3. 設定した条件
  4. 実測データとの比較
  5. BB弾の様々な軌道
  6. 「良い軌道」の詳細なデータ
  7. 動画と静止画
  8. 「良い軌道」へ調節するには
  9. 温度の違いについて
  10. 僕が実際に試した時の報告
  11. お礼
  12. アイディア募集
  13. 参考文献

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果←今ここ
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算

目的


屋外で行われるサバイバルゲームで優位に立つために、BB弾は重い球が良いのか、軽い球が良いのかを数値計算によって,良い軌道を描く弾道の軌道のパラメータを探索して、現実で調節するための方法を模索します。

結論


サバゲーを行うにあたり、風があったり、遠くを狙う必要のある屋外では

BB弾は重い方が良い

という結論が得られました。
BB弾は重ければ重いほどまっすぐ、遠くまで飛ぶ。また、当たった時に判定がされやすい事になります。
屋内などの風が無く、接近戦(30m以内)が主になる場合ではほぼ差は見られないため、お値段的に軽い球を使うのがいいと思います。

重さの違いにありがちな話は本当?


「0.20gのほうが着弾が早い」
   →ほとんど嘘
    初速は確かに0.20gのほうが早いですが、最大でも15mの地点で0.01秒程度の違いしかありません。30m以上では0.25gのほうが早く着弾します。

「0.25gはまっすぐ飛ぶ」
   →本当
    →風の影響を受けにくいため、横にぶれにくいだけでなく、縦にもぶれが少くなります。

「0.20gのほうが遠くに飛ぶ」
   →
    →0.25gのほうが遠くまで飛びます。

「0.25gのほうがヒットを取りやすい」
   →本当
    →弾の持つエネルギーは0.25gの方がどの地点でも大きくなるため、当たった時に判定されやすいことになります。

まとめると、
0.20gが優れているのはコストだけ
ということになります。


スポンサーリンク



以下、詳細なデータを載せます。

設定した条件


ニュートン力学の範囲内で計算します。
考慮した力と条件は、

  1. 重力
  2. 空気の慣性抵抗力、粘性抵抗力
  3. 回転による揚力(Magnus力)
  4. BB弾の回転の減衰

です。コリオリ力は重力に比べ3桁程度小さい力なので考慮に入れていません。
空気の粘性率は、空気の密度は気温20度、乾燥した空気中のものを使用しています。
僕が調べた範囲で他のページではBB弾の回転の減衰が考慮されていませんでした。本稿ではそれを考慮したことが一つの特徴です。
運動方程式の導出や、理論に関してはこちらの弾道計算(BB弾)の理論をご覧ください。

実測データとの比較


計算が実測のデータと合っているのかをまず初めに検証します。
東京マルイ G36CにおけるノーマルホップとG-HOPの弾道特性実測データと、本稿の数値計算結果との比較を行います。
実験と比較(結果用)_c
ここで数値計算結果(緑線)との比較対象は紫色の点です。
※実測データにおいて、BB弾の回転数は測定されていなかったため、最高到達高さを元に回転数を決定しています。

全体的に非常に良い一致を見せていることが分かるかと思います。50m,60m地点においては10cm程度の差しかないのではないでしょうか。この実測データの再現、という点においては十分に信頼できる数値計算結果になっていると思います。
ただし、打ち出し直後の浮き上がり具合に若干の差があります。
このずれの原因として、実測データの高低さの加味が十分ではない、ということが考えられます。
実測時ではなく、改めて高低差を加味した、と実測データのページにあるので、場所が完全に同じではないかもしれないと考えられます。
風で揺らいでいる、ということは無いでしょう。もしも風で揺らいでいるとしたら、遠くで差が出るはずです。射出直後でこれだけ揺らいで遠くで揺らがないというのは納得ができません。

本稿の数値計算は実測データと十分に一致している、と確認したうえで、いろいろとパラメータを変えて計算してみました。

BB弾の様々な軌道


まず、一番気になるBB弾の軌道を見ます。
エネルギーを0.90Jで固定して、パラメータを様々に変え0.20gと0.25gのBB弾で比較してみます。
その結果がこちらです。黒文字が0.20g, 赤文字が0.25gであることを意味しています。
弾道軌道09J_c
BB弾が描く良い軌道とは

  • まっすぐ飛ぶ
  • 遠くまで飛ぶ

軌道であることだ、と考えます(左右のブレはここでは除きます)。

すなわち「良い軌道」とは上下振れ幅が最小になる軌道のことだと言えるでしょう。
上の画像の例でいうならば
0.20g → BB弾が35回転/sで回転し、下向きに1.2°方向に射出されたとき
0.25g → BB弾が42回転/sで回転し、下向きに1.2°方向に射出されたとき

がそれに近い軌道でしょう。

0.25gでこの「良い軌道」に合わせることができたなら、50m先の敵までは、喉~胸あたりに照準を合わせることで確実にヒットが得られるでしょう。
また、60mの距離を狙いたければ敵の頭の1.26m上を狙うことで一応届きます。

もう少し、この「良い軌道」について調べてみましょう。

「良い軌道」の詳細なデータ


ここでは、重さの違う0.20gと0.25gとを比較するために、ゼロインを50[m]に合わせ、上下方向の振れ幅が最小になる時で比較します。

弾道詳細_ゼロイン50m2_c
ゼロインを50mに合わせた時、0.20gでは60cmの振れ幅、0.25gでは40cm程度の振れ幅になります。
0.20gよりも0.25gの方が小さい振れ幅です。予想通りといえば予想通りの結果ですね。

次に着弾時刻を見てみましょう。30m地点までは0.20gの方が若干早く着弾し、それ以降は0.25gのほうが早く着弾します。50mの遠距離で比較しますと、0.20gと0.25gとの着弾時刻の差は約0.13秒。まばたき程度の差があります。
ちなみに15mでの着弾の時間差は0.01秒です。とても微々たるもので、この差は認識できないと思われます。よって着弾時刻も上下の振れ幅も少ない、遠距離にも対応できる0.25gを使うのがいいでしょう。

また、エネルギーは発射した瞬間から着弾するまで、どこでも0.25gの方が高いため、遠くでも判定がされやすいことになります。

遠距離を狙いたい場合は着弾速度、上下振れ幅に影響するだけでなく、判定もされやすい重い球を使いましょう。


2016/03/13
計算コードを変更し、動画化しました。点の間隔は0.04秒で、gifアニメの速度の時間間隔と実際の時間間隔は同じにしてあります(すなわち、gifアニメの1秒と現実の1秒が同じ)。
bullet_orbit
こちらは静止画。
orbit_st_c

ゼロイン距離と上下方向振れ幅との関係


さて、上下方向の振れ幅を最小限に抑えながらゼロインを適当な位置に合わせたい場合について考えます。
遠くにゼロインを合わせたい時、飛ばせばと飛ばすほどホップを強くかけなければなりません。その結果、上にも下にもブレることになります。
この関係性がどうなっているかを調べますと、以下のようになります。
zeroin_width2_c

エネルギー0.9Jの場合で、上下の振れ幅を最小に抑えるパラメータで考えますと、段々と上下方向の振れ幅が大きくなっていることがわかります。
これは空気抵抗によって弾速が減少するためにその分、ホップによって時間を稼がなければならないことに由来します。
例えば、
0.20gのBB弾で50mにゼロインを合わせた時は上下振れ幅が60cm(射出位置から30cm下がり、ホップによって射出位置から30cm上昇し、ゼロインを迎える)ですし、
0.25gのBB弾で50mにゼロインを合わせた時は上下振れ幅が40cm(射出位置から20cm下がり、ホップによって射出位置から20cm上昇し、ゼロインを迎える)です。
この重さの違いによる上下方向の振れ幅は顕著に表れます
60mなんかにゼロインを合わせようとすると、0.20gでは上下方向に1.6mずれることになり、0.25gでは上下方向に1mずれることになります。
もしかしたら、0.20gで60mに合わせた場合、フィールドで立って撃たない限り、先に地面についてしまってダメになるでしょう。

これらから、遠距離ではますます重いBB弾を使ったほうが良いという結果が得られるわけです。

※ただし、中距離に当てるつもりがない場合はホップをある程度かけて上向きに射出したほうがいいでしょう。

「良い軌道」へ調節するには


さて、サバゲーに行き、シューティングレンジで調節をすることを考えます。
BB弾の回転数など知っているわけがありません。また、角度1.2°などそんなこと分かるわけがありません。
さらに、BB弾がどこまで飛んだかは大抵の場合は10m単位でしか分かるはずがないのです。
せめてジュール数(弾速)は知っているとします。この状況で良い軌道へ調節することを考えましょう。

0.25g, 0.9J(初速84.85m/s)の場合

ここでは、0.25g, 0.9J(初速84.85[m/s])の場合を考えます。

50mにゼロインを合わせることを考えます。
調節方法_c

この場合の手順は、
手順1,    射出する銃口の高さ50m地点にある同じ高さにある的に照準を合わせる。

手順2,    ホップ無しで、50m先の的に照準に合わせたままBB弾が地面に落ちた時の位置が、おおよそ
    \(10\times(\mbox{射出高さ[m]}) +10[\mathrm{m}]\)
    になるように調節する。

手順3,    その照準に合わせたまま、徐々にホップを掛けていき、50[m]の的に当たるまでホップを掛けて調節する。

です。

詳しく説明すると、
手順1では角度を決めるための基準点を探しています。
手順2では角度をホップ無しの場合で、到達距離を見ることによって決定します。何点か計算したら、射出する高さと良い角度の時の落下地点は、\(10\times(\mbox{射出高さ[m]}) +10[\mathrm{m}]\)という関係でした。
なので、地面からの銃口の位置が分かれば、良い角度での射出は、到達距離を見れば分かるのです。
※空気抵抗のない場合、\((\mbox{落下距離})\propto \sqrt{射出高さ}\)になります。ただし、短い区間であれば線形補間で良いでしょう。ちゃんと線形補間すると、\(10.40\times(\mbox{射出高さ}[\mathrm{m}])+9.87[\mathrm{m}]\)となるのですが、実用上は簡単な\(10\times(\mbox{射出高さ[m]}) +10[\mathrm{m}]\)で十分でしょう。

手順3では、その照準のまま、だんだんホップを掛けていき、50mに合わせればいいのです。

0.20g, 0.533J(初速73m/s)の場合

もう一つ具体例、0.20g, 73m/sの場合を考えます。
これは、僕が持っている電動のmp7a1のなにもカスタムしていない時の値です。
調節方法m0200533J_c
50mに合わせると上に40cm、下に40cmも上下するので、無理に狙わず、40mで合わせています。
(※上のグラフとはエネルギーが異なるため、”ゼロイン距離と上下方向振れ幅との関係”で紹介したグラフとは上下振れ幅が異なっています。)
この時のパラメータは、エネルギー0.533J, 下向きに角度1.81度、回転数42回転/sです。
mp7a1はこちら。

2丁持っている場合


2丁銃を持っている場合、精密射撃用と遠距離用と分けることができます。
30m内では上下方向のブレが6cm程度になるので、近接用と割り切りってしまい、
遠距離では上のパラメータで調節しておけば、どんな距離でも気持ち良く狙えます。
例えば、基本的に上のパラメータで調節した銃を使っておき、相手が隠れた時など、壁の小さい隙間を狙って近距離用で精密射撃を行う、なんてことができるかと思います。

温度の違いについて


夏場と冬場ではどのくらい弾道が変わるのでしょうか?
0.28g, 0.9Jで、温度10℃、20℃、30℃で撃った時の軌道の違いを考えます。
紫が10℃、緑が20℃、青が30℃です。
基準は20℃です。20℃でゼロイン50mの時の最小の振れ幅に合わせた時に他の温度ではどうなるかを考えます。
弾道の温度依存性_c
結論として、10℃変わるごとに約1m飛距離が変わります。
なので夏場の方が冬場よりも2m程度よく飛ぶことになります。

僕が実際に試した時の感想


2015年12月、内部をカスタムしていないmp7a1でopsに行き、実際に撃ってきて理論通りの軌道になるだろうか、と試してきました。
レーザーサイト、その他もろもろは持っていないので考察ではなく、感想としてみてください。
僕個人の感想をまとめると、

  • 0.25gと0.20gの差は明らかに感じられる。特に集弾性が素晴らしい。
  • 最小の振れ幅になる軌道は確かにその通りになるが、少なくとも起伏の激しいフィールドでは使いにくい。

でした。
0.20gと0.25gの差は明らかだと感じました。軌道に対しては体感では少しだけ分かる程度でしたが、集弾性がとにかくいいのです。着弾点でばらけることがほぼ無いのです。実際にどの程度影響していたか分かりませんが、0.20gから0.25gに変えた時、キル数が上がりました。
とても気持ちよく戦えるので、初心者の方を誘ったりする際は是非0.25gを使わせて、サバゲーマーの仲間入りをさせてください。

一度下がってから上がる、これは理論上では一番良い軌道であることは確かなのですが、起伏の激しいフィールドでは、この”一度下がる”が非常に厄介です。一度下がるため、そこに障害物(例えば斜面等越しに当てるとき)があったら当たらないのです。少し上向きに撃とうものなら強ホップのため、敵の頭上を通り越します。
結果、敵から反撃を受け、こちらからは当てられない状況が生まれました。
起伏の激しいフィールドでは、30~35m程度まで数cmの振れ幅でまっすぐに進む軌道が一番良いと思います。
それよりも遠い場合は山なりに何発か撃てばいいので、非常に適していると感じました。

今度は平坦なフィールドで試してみたいと思います。

スポンサーリンク

お礼


@hthsi アイディア提供。実際にサバゲーをやっている上で気になること、○○の違いが知りたい!といったアイディアを弾道計算を行う初期の頃からいっぱいいただきました。重さの違いにありがちな話はほんと?の部分があるといいよね、等の分かりやすい説明は@hthsiの提案です。ありがとうございます!
また、サバゲーるにて告知をしていただきました。
0.25g弾の優位性について

募集


2015/07/06
現在も計算を行い、まとめている最中です。気になること、面白いアイディア、計算に関する質問等がありましたらコメントでもtwitter(@sig_colon)でも、メールでも構いませんのでお知らせください。
計算するまでもない質問の場合はコメントで返信、特に面白そうな場合は、綺麗にまとめて本稿に追加いたします。

2016/02/07
Tech Report  ハイスピードカメラによるエアソフトガンの考察
より、ガスガンと電動ガンの銃口付近での空気の流れが見て取れます。
非常に素晴らしく、実際の状況なので信頼できます。15000fpsの高速度カメラだそうです。

ひとつ、気になるのは報告されているBB弾の回転数が私の導き出した値と2倍近く違うことです。
回転数が違うことによって、数値計算結果の回転数以外に影響はありません。ただのホップの強さの目安と見ているので、もしも間違いがあったとしてもホップの強さと回転数との換算が違うだけです。結果には影響しません。
管理人様にお尋ねしたところ、回転数は管理人様自身が導き出したものではなく、考察程度としてみるのがよい、と返答いただきました。残念ながら確証が取れないので回転数の議論はどなたかやっていただけるまで待つことにします。
動画があり、非常に説得力のあるサイトなので、是非ともご覧いただくことをお勧めします。

参考文献


東京マルイ G36CにおけるノーマルホップとG-HOPの弾道特性

2016/03/22 追)
ガンジニア様のハイスピードカメラによる映像


より、本サイト内で表記している回転数にミスがある疑いが濃厚になってきました。
サイト内に「○○回転している」や「○○回転/s」と言う表記があるのはたぶんです。
これから検証していきますのでお待ちください。

2016/07/09追)この問題は解決しました。詳細はBB弾の回転量について(実験との比較)をご覧ください。

コメント返信用、画像

2015/07/27
>ウエダさん宛て
縦軸:高さ[m]
横軸:到達距離[m]
ウエダさんパラメータ_c

弾道計算(BB弾)の理論

エアガンから発射したBB弾がどのように飛んでいくか?
BB弾の弾道計算理論パートです。

完全な計算のためには圧縮ナビエ・ストークス方程式を解いて…ということをしなければなりませんが、そこまでは考えません。
射出後のBB弾がどのような軌道を描き進むのか、ニュートンの運動方程式の範囲で考えていき,定量的な一致を目標にしていきます。

結果として得られる2つの運動方程式は以下の2つです。
・BB弾の運動方程式
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-C_l \frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
\)

・BB弾の回転角速度の減衰を記述する方程式
\(
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
\)

の二つにより構成されます。

実験データから、\(C_l=0.12\)で良く記述できることが分かりました。
詳細はBB弾の回転量について(実験との比較)をご覧ください。

結論を言うと、サバゲーで勝ちたければ、重いBB弾を使うべきです。
飛距離を伸ばすためには
BB弾の重さ、銃口から出てくるときの初速回転数以外に飛距離を伸ばすことに関係する変数はありません。
0.20,0.25gの違いと回転減衰
※初速と共に載せました。この初速の値は実験データがあるサイトと同じ値にしています。

目次

  1. 先行研究の紹介
  2. BB弾に働く力
  3. BB弾の運動方程式
  4. 結果
  5. 流体中に置かれた球体の回転の減衰
  6. 参考文献

弾道計算に関するその他ページ
弾道計算(BB弾)の理論←今ここ
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算

先行研究の紹介


誰かやっているかなーっと調べてみると、やっぱりやってる人がいました。

解析解(慣性抵抗とマグヌス力を考慮?)を用いての、定性的な一致に終わっています。
BB弾の弾道計算まとめ -チン太さんのノート(2013年)
(※2016/01/16 正しさの確認がとれないので一応の解析解、ということにしておきます。なぜなら、2次元の空気抵抗を考慮した微分方程式でさえ解析解が存在しないはずだからです[12]。)

本稿では可能な限りの影響を考慮し、定量的に一致させることを目標にしていきます。

BB弾に働く力


質量mのBB弾の場合、運動方程式は言葉で書けば、
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=(\mbox{重力})+(\mbox{風の抵抗力})+(\mbox{回転による揚力})
\)

で十分だと思います。
コリオリ力は今の場合、重力に比べて3桁程度小さい力であるため考慮に入れる必要はないでしょう。
また、角速度ベクトルの歳差運動はあるかもしれませんが、見積もる方法が分からなかったため考慮していません。

追記)
各々の力の大きさを簡単に計算してみました。BB弾の速さを100m/s, 重さ0.2g, 166回転/sと仮定

重力: \(1.96\times 10^{-3}\)
粘性抵抗力: \(1.03\times 10^{-4}\)
慣性抵抗力: \(8.52\times 10^{-2}\)
Magnus力: \(1.71\times 10^{-3}\)
コリオリ力(最大): \(2m|\vec{\omega_{earth}}||\vec{v}|=2.8\times 10^{-6}\)

であるため、コリオリ力はほとんど無視できるんじゃないでしょうか。

重力\(F_g\)


場所による重力加速度の違いはあるかもしれませんが、無視できるでしょう。
なので重力\(F_g\)の大きさは
\(F_g=mg\)
でしょう。ここで

  • \(m\ \mathrm{[kg]}\) : 物体の質量
  • \(g\ \mathrm{[m\cdot s^{-2}]}\) : 重力加速度

です。

風の抵抗力\(F_d\)


空気抵抗と呼ばれているものです。
空気抵抗を計算するには物体の形状が分からなければ求めることはできません。
速度に比例する項は理論的に求めることができますが、速度の2乗に比例する空気抵抗の係数(抗力係数)はレイノルズ数の関数となります。

ここでは、半径\(R\)の完全な球体に働く空気抵抗を求め、それを導出、計算します。

風(流体)の抵抗力の大きさ\(F_d\)は一般には次元解析により
\(
\displaystyle F_d=\frac{\eta^2}{\rho}\sum_{n=1}^{\infty}K_n \left(\frac{v \rho L}{\eta}\right)^n
\)

として与えられます[1]。ここで

  • \(\eta\ \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) : 流体の粘性率(Viscosity) ※動粘度ではないです。
  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 流体の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 物体の速度の大きさ
  • \(L\ \mathrm{[m]}\) : 物体の大きさ(半径\(R\)の球の場合、\(L\)は直径\(L=2R\)に相当する量)
  • \(\frac{v \rho L}{\eta}\) : レイノルズ数(\(R_e\),Reynolds number)と呼ばれる無次元量

です。次元解析から具体的な定数の値\(K_n\)を求めることはできません。
通常の場合は\(n=1\)(粘性抵抗)と\(n=2\)(慣性抵抗)のみ考慮されます。
速度の3乗以上に比例する項というのは見かけません。
おそらく、これから説明する抗力係数\(C_d\)にその効果を全部押し込めているためでしょう。

上の式の\(n=1,n=2\)をそのまま書き下すと
\(F_d=K_1\eta L v +K_2\rho L^2 v^2\)
となりますが、慣例として同じ次元を持つように式を変形させて
\(F_d=C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right)\)
の形で記述されます。ここでSは速度に垂直な、物体の断面積です。

\(C_1\)は物体が完全な球で、小さく、その速度がゆっくりである場合(レイノルズ数が2以下くらい)、流体の運動を記述する方程式、
ストークス方程式(Stokes’ law、ナビエ・ストークス方程式のレイノルズ数が小さい極限の方程式)
から導くことができて\(C_1=3\pi\)であることが理論的に導けます。

\(C_d\)は抗力係数(Drag Coefficient)と呼ばれ、
これはレイノルズ数の関数となり、
\(C_d=C_d(R_e)\)
となります。
この抗力係数は非常に厄介らしいです。物体の形状に依存したりするらしく、球の場合でも、どうやら実験的に求められているようです。
Wikipediaの抗力係数のページに若干の記述があります。
球の場合の抗力係数の関数\(C_d(R_e)\)をフィッティングした関数は論文[2]より、
\(
\displaystyle C_d=\frac{24}{R_e}+\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
+ \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}} + \left(\frac{R_e^{0.80}}{461000}\right)
\)

として書けるそうです。実際にプロットしてみるとこんな感じです。
Cdグラフ
確かに抗力係数のグラフと同じようになります。

よって半径Rの球の場合、風の抵抗力\(F_d\)の大きさは結局、
\(
\begin{align}
F_d &= C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right) \\
&= 3\pi \cdot 2R\cdot \eta v +\frac{1}{2} C_d \rho \pi R^2 v^2 \\
&= 6\pi R \eta v + \frac{1}{2} C_d \rho \pi R^2 v^2
\end{align}
\)
として与えられます。

※BB弾をエアガンから、無風状態で発射する場合、抵抗は

粘性抵抗力: \(1.03\times 10^{-6}|\vec{v}|\)
慣性抵抗力: \(8.52\times 10^{-6}|\vec{v}|^2\)

となります。\(|\vec{v}|\cong 10\sim 100 \mathrm[m\cdot s^{-1}]\)くらいであるため、粘性抵抗力は慣性抵抗力のおおよそ0.1%位ということになります。

回転による揚力\(F_L\)


回転による揚力の発生はマグヌス効果(Magnus effect)と呼ばれています。
空気力学で、回転する円柱の揚力に関する理論クッタ・ジュコーフスキーの定理(Kutta–Joukowski theorem)によって回転による揚力は記述されます。
発生する揚力\(F_L\)の大きさは
\(
F_L=\rho \Gamma v
\)

と書けます[3]。ここで

  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 流体の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 物体の速度の大きさ
  • \(\Gamma\ \mathrm{[m^2\cdot s^{-1}]}\) : 循環

を意味します。循環\(\Gamma\)は
\(
\displaystyle \Gamma=\oint_s \vec{v}\cdot d\vec{s}
\)

という量です。
ここで、\(\oint\)は球体の外側での閉経路に対する線積分であり、球体そのものの循環ではありません[9]。

私は流体力学は良く知らないで余り突っ込まないでおきますが、簡単に仮定します。
実際に発生する循環\(\Gamma\)を無次元の係数\(C_l\)と球体の循環\(\Gamma^{sp}\)を用いて、
\(
\Gamma=C_l\Gamma^{sp}
\)

で書ける、とします。
(実際には迎え角なども関係するようなのでこれは厳密には正しくはありませんが、定性的に説明するので良いでしょう、という仮定しているだけにすぎません。)

球の循環を考えます。
半径\(R\mathrm{[m]}\)の角速度\(\omega\mathrm{[s^{-1}]}\)を持つ円盤に対して循環\(\Gamma^{ci}\)を計算すると、
\(
\begin{align}
\Gamma^{ci} &= \oint_s \vec{v}\cdot d\vec{s} \\
&= \int_0^{2\pi}R\omega\cdot Rd\theta \\
&= 2\pi R^2 \omega
\end{align}
\)
と計算できます。しかし、これから考えるのは球です。
半径\(r\)の円盤に対して循環は、\(\Gamma^{ci}=\Gamma^{ci}(r)=2\pi r^2 \omega\)と書くことができるので、積分して球の循環\(\Gamma^{sp}\)が求められます。
球に対する循環の計算
図のように考えれば、計算は、
\(
\begin{align}
\Gamma^{sp} &= \int_{-R}^{R} \Gamma(r) dl \\
\end{align}
\)

です。\(l=R\cos\phi, \ \ r=R\sin{\phi}\)の変数変換を行って計算すると、
\(
\begin{align}
\Gamma^{sp} &= \int_{-R}^{R} \Gamma(r) dl \\
&= 2\pi \omega R^3 \int_0^\pi \sin^3\phi d\phi \\
&= \frac{4}{3}\pi R^3 \cdot 2\omega
\end{align}
\)
と求められるため、揚力\(F_L\)は無次元の係数\(C_l\)を用いて
\(
\begin{align}
F_L &= C_l \rho \Gamma^{sp} v \\
&= C_l \frac{4}{3}\pi R^3 2\omega \rho v
\end{align}
\)
と書き表せるとします。

2016/07/09 追)
係数\(C_l\)は実験[11]とその弾道データより、
\(
C_l=0.12
\)

で良く近似される、という事が分かりました。詳しくは
BB弾の回転量について(実験との比較)
をご覧ください。

2016/07/09 追)
このページのここより以降の結果に対しては書いてある回転量に8.3倍してください。そうすると現実での、実際の回転量になります。
(例えば”20回転/秒”とあったら、これは”166回転/秒”と言うように読み変えてください)



BB弾の運動方程式


重力\(F_g\), 空気抵抗\(F_d\)と回転による揚力\(F_L\)を取り入れたBB弾に対する運動方程式は方向も考慮に入れて、
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}
+\left\{-6\pi \eta R |\vec{v}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{v}|^2\right\}\frac{\vec{v}}{|\vec{v}|}
-C_l\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{v}|\frac{\vec{v}\times\vec{\omega}}{|\vec{v}\times\vec{\omega}|}
\)

と書けます。上の式の場合、風速の影響は入っていません。
風速の影響を取り入れるのは風と共に動く座標系で考えてから元の座標系に戻ればいいです。
すなわち、地表に固定した座標に対して風の速度ベクトル(定ベクトル)を\(\vec{u}\), 物体の速度ベクトルを\(\vec{v}(=d\vec{r}/dt)\), とし、
相対的な速度ベクトルを\(\vec{V}=\vec{v}-\vec{u}\)と定義すると、風がある中での運動方程式は、
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-C_l\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
\)

と書けます。

それぞれの定数の具体的な値を考えます。

各定数の値
名称 記号 捕捉
重力加速度\(g\) \(9.80665\mathrm{[m\cdot s^{-2}]}\)
空気の粘性率(Viscosity)\(\eta\) \(18.2\times 10^{-6}\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) 空気中20度において。また、粘性率は数十気圧から数気圧の間圧力には依存しません。
が、温度に依存し、ある温度\(T_1[K]\)での粘性率\(\eta_1\)が分かっているならば、
温度\(T_2[K]\)での粘性率\(\eta_2\)はサザランドの定数Cを用いて、
\(
\eta_2=\eta_1 \left(\frac{T_1[K]+C}{T_2[K]+C}\right)\left(\frac{T_2[K]}{T_1[K]}\right)^{3/2}
\)
と記述されます。空気の場合、C=117です[4]。
乾燥した空気の密度\(\rho\) \(1.205\mathrm{[kg \cdot m^{-3}]}\) 乾燥した空気の密度は温度と圧力に依存します。温度\(t[℃]\)、圧力をH[torr]とすると
\(
\rho=\frac{1.293}{1+0.00367t[℃]}\cdot \frac{H[torr]}{760}
\)
です[5]。
BB弾の半径\(R\) \(3.0\times 10^{-3}\mathrm{[m]}\)

数値計算は刻み幅制御ルンゲ・クッタ法を用います。
コードはこちら→弾道計算(BB弾)のコード


スポンサーリンク


結果


結果をまとめます。
※ここ以降に表示されている回転数に対しては8.3倍してください。

回転数と飛距離依存性

回転数依存性

最終的な到達距離はBB弾の質量にほとんど依存しないことが分かります。
また、BB弾の重さにより同じ回転数ではホップが効きすぎたりということがあるので、重さによる調整は欠かせないですね。

理想的な弾道の軌跡

0.25gの場合、BB弾の回転数25回転くらいが一番よさそうな軌道として考えます。
この時打つ際に若干、下方向に向けて打つことでまっすぐ飛ぶ領域を増やせそうです。どのくらい下向きに打てばいいんでしょうか。その結果がこれです。
下向き射出角度依存性
※追)上で述べた通り、回転数はあてにしないでください。実際の回転数は上の3~4倍したものです。
この結果から、一つの最適な軌道を与えるのは回転数25回転/s, 角度0.2°で射出した時ですね。その時の上下方向のブレは25m位飛んで3cmくらい。
また、一つ言えるのは一番いい軌道は角度\(\theta\)を変化させたときに上下方向のブレが少ないとき、と言えると思います。
なのでこの値が最適ではないです。もう少し調べる必要があるわけですね。ホップアップの回転数に依存して飛距離と角度が決定されています。

BB弾の重さ0.20gと0.25gの違い

速さ,到達距離の時間依存性
回転なし、風なしの場合に到達するまでの時間を見ると30[m]まで0.20gのBB弾のほうが早く到達するようですね。

重さと風の影響

よく飛ぶであろう0.25gは25回転/s, 0.20gは20回転/sとして計算をしています。速度は0.9Jのエネルギーを持つようにしています。
風の影響
0.2gと0.25gの違いと風の影響は余りなさそうな気がします。例えば30[m]地点ではどちらも5[m/s]の風が吹いていたとしても1.0~1.3[m]のずれに収まっています。
この30cmが問題になるのかは…状況次第ですね。

横向きに打ち出した時の曲がり具合

横向き射出2
意外と曲がらないものですね。大体2[m]です。
理由は、この計算では打ち出す時は銃を完全にまっすぐにして固定した状態なのでMagnus力の影響が0から段々と強くなります。
Magnus力の強さは回転軸に垂直の速度成分が速いほど加速度的に増えていくものなので、はじめ銃を曲がる方向に向けて撃てば
曲がる方向の初速度+更なるMagnus力でもっと曲がります。
なので物陰にいる敵を倒す(狙っちゃだめですが)場合は若干敵の方向へ向けて回転による効果を得られるようにしましょう。

銃を逆さまにして上から狙う

上から落として狙う
じゃあ銃を逆さまにして上向きに撃てば障害物の上から狙えるんじゃ・・・?
適正なホップ(m=0.20[g],20回転)の場合、2[m]位の障害物があれば上の黄色Θ=0.64°の軌道でよさそうです。
しかし、障害物を超えて1m落下するのに5[m]位離れてしまいます。なので、敵が障害物から5[m]程度離れている時に撃てば当てられますね。

2015/4/28 追)

実際の弾道のデータとの比較

BB弾の実際の飛距離をデータとして掲載しているページを見つけました。ここです↓
東京マルイ G36CにおけるノーマルホップとG-HOPの弾道特性
測定は大変なものなので緻密なデータとはいきませんが、大体の傾向と、この計算がどのくらいよいものなのか、はわかりそうです。
BB弾の回転数は上のページを参考にし、最大到達点を目安に決定しました。
0.2gのデータは良いデータではなさそうだったのであきらめ、0.25gの結果のみで弾道の軌道を比較しましょう。
数値計算では無風状態でBB弾の重さ0.25g, 初速77.51[m],高さ1.1[m]から射出します。
回転数はおおよそ41回転の時、最高到達点が2.4[m]程度と数値計算によって分かったためそれを使います。その時のデータは以下のようになりました。
実際のデータと比較
比較対象はピンク色の線です。割と良い一致をしています。
原因を考察します。どの程度の統計性があるかは分かりませんが、写真を見る限り、数十発程度の平均値であろうと予想できます。
なので、毎回50m地点で向かい風が吹いて失速し、それが数値計算とのずれだ、などという状況は考えにくそうです。
一番考えられるのはBB弾の回転角速度の減衰でしょう。次の節で取り扱います。

そのほか、この計算に一致しない理由としては、温度による粘性抵抗のずれ、
若干の向かい風がある(計算よりも早くホップにより上昇しているため)、横方向の考慮、測定ミスetc…上げればきりがありません。
しかし、おおよその傾向はグッドです。ホップが効いて下がり始めるまではほとんど同じ軌跡を辿っています。

ちなみに、地面に着弾するまでの時間は計算上では2.27秒でした。

追)2015/05/31

流体中に置かれた球体の回転の減衰

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


最終的に、最初に角運動量がz軸方向の成分しかない場合、角速度の時間変化\(\frac{d\omega_z}{dt}\)は、
\(
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
\)
と求められます。(ここで\(I\)は慣性モーメント、\(C_f\)は摩擦抗力係数、\(u\)は物体の、風との相対速度を表します。)
これを計算し、BB弾の角速度に反映させます。

ここから自信が無いので、あっているか、物理の背景が正しいかは保証できません。
まず、回転を記述する運動方程式は、
\(
\displaystyle \frac{d\vec{L}}{dt}=\vec{N}
\)

慣性モーメント\(I\)を用いると、\(\vec{L}=I\vec{\omega}\)であるため、
\(
\begin{align}
I\frac{d\vec{\omega}}{dt}&=\vec{N} \\
I\frac{d^2\vec{\theta}}{dt^2}&=\vec{N}
\end{align}
\)
ここで\(\frac{d\vec{\theta}}{dt}=\vec{\omega}\)と置いています。

求めるべきは右辺の力のモーメント\(\vec{N}\)です。

摩擦抗力\(D_f\)、摩擦抗力係数\(C_f\)というものがありました[6,7,10]。
境界層流れ、という分野に分類されるそうです。
流れている流体に対して摩擦によって平板に生まれる、摩擦抗力\(D_f\)は
\(
\displaystyle D_f=-C_f\frac{\rho U^2}{2}S_{contact}
\)

と記述でき、\(C_f\)は摩擦抗力係数、\(\rho\mathrm{[kg\cdot m^{-3}]}\)は流体の密度、
\(U\mathrm{[m\cdot s^{-1}]}\)は流体との相対速度、\(S_{contact}\)は流体と物体の接触面積を記述しています。
\(C_f\)はレイノルズ数\(R_e\)の関数であり、
層流領域(\(R_e\lt 5\times 10^5 \sim 10^6\))では[10]
\(
\displaystyle C_f=\frac{1.328}{\sqrt{R_e}}
\)
乱流領域(\(5\times 10^{5}\lt R_e \lt 5\times 10^6\))では[6]
\(
\displaystyle C_f=\frac{0.074}{{R_e}^{1/5}}
\)
乱流領域(\(10^{7}\lt R_e\))では[6]
\(
\displaystyle C_f=\frac{0.455}{(\log_{10}{R_e})^{2.58}}
\)

を用いるのが推奨されるようです。
専門家ではないのでこれらの係数の導出等々は知りません。
エアガンの場合、レイノルズ数はおおよそ\(3\times 10^{4}\)程度なので層流領域だと見積もられます。

この摩擦抗力を元に、空気抵抗によって回転が減衰することを調べてみましょう。
ただし、ここでは2つの仮定をします。

  1. 風は角度\(\varphi, (0\le\varphi\lt 2\pi)\)のみに依存する
  2. 位置\(\vec{r}\)で速度\(\vec{v}(\vec{r})\)を持つとき、その微小面積\(\Delta S\)に働く摩擦抗力\(\Delta D_\vec{f}\)を
    \(\displaystyle \Delta D_\vec{f}=-\frac{\rho C_f}{2}|\vec{v}(\vec{r})|^2\frac{\vec{v}(\vec{r})}{|\vec{v}(\vec{r})|}\Delta S\)
    と仮定

これらの仮定の下、話を進めます。
微小区間に働く力のモーメント\(\Delta\vec{N}=\vec{r}\times \Delta\vec{f}\)を考えてから、全空間で積分を行います。
そのために位置\(\vec{r}\)の点における速度\(\vec{v}\)を導出します。
減衰_回転と風_c
摩擦抗力は相対速度\(\vec{v}\)によって記述されると考えられるので、回転によって生じる速度\(\vec{v}_{rot}\)とその地点での風の速度\(\vec{v}_{win}\)を使えば、
\(
\vec{v}=\vec{v}_{win}-\vec{v}_{rot}
\)

と書くことができます。
まず、\(\vec{v}_{rot}\)は、
\(
\vec{v}_{rot}=\vec{r}\times \vec{\omega}
\)

と表現されます。また、\(\vec{v}_{win}\)は、
\(
\vec{v}_{win}=(\vec{u}\cdot \vec{e}_{\varphi})\vec{e}_{\varphi}
\)

であると考えられます(仮定①)。ここで\(\vec{e}_{\varphi}\)は回転方向の単位ベクトルを表します。
よって、位置\(\vec{r}\)の相対速度\(\vec{v}\)は、
\(
\vec{v}=(\vec{u}\cdot \vec{e}_{\varphi})\vec{e}_{\varphi}-(\vec{r}\times \vec{\omega})
\)

となります。
微小区間に働く力のモーメント\(\Delta\vec{N}\)は、
\(
\begin{align}
\Delta\vec{N}&=\vec{r}\times \Delta\vec{f} \\
&=-\vec{r}\times\left(\frac{\rho C_f}{2}|\vec{v}|^2\frac{\vec{v}}{|\vec{v}|}\Delta S\right) \\
&=-\frac{\rho C_f}{2}|\vec{v}|(\vec{r}\times\vec{v})\Delta S
\end{align}
\)
と書かれます(仮定②を使用)。
今、3次元極座標を考えて積分すると位置\(R\)での微小面積はヤコビアンより、\(R^2\sin(\theta)d\varphi d\theta\)なので、
(動径方向は\(r=R\)のみ)、
\(
\displaystyle \vec{N}=-\frac{\rho C_f}{2}R^2 \int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |\vec{v}|\cdot(\vec{r}\times\vec{v}) \sin{\theta}
\)

と表せます。
あとはこれを解いて剛体の運動方程式に当てはめれば、回転の減衰が表せます。


具体的な形を当てはめていきましょう。
半径\(R\)の球体に対して、3次元極座標で考えます。すると、位置\(\vec{r}\)は、
\(
\vec{r}=
\left( \begin{array}{c}
R\sin\theta \cos\varphi \\
R\sin\theta \sin\varphi \\
R\cos\theta
\end{array} \right)
\)
であり、回転方向の単位ベクトル\(\vec{e}_{\varphi}\)は、
\(
\vec{e}_{\varphi}=
\left( \begin{array}{c}
-\sin\varphi \\
\cos\varphi \\
0
\end{array} \right)
\)
です。
回転角速度\(\vec{\omega}=(0,0,\omega)\)、風の方向\(\vec{u}=(u,0,0)\)
であると考えれば、
\(
\vec{r}\times \vec{\omega}=
\left( \begin{array}{c}
R\omega \sin\theta\sin\varphi \\
-R\omega \sin\theta\cos\varphi \\
0
\end{array} \right)
\)
であり、
\(
(\vec{u}\cdot \vec{e}_{\varphi})\vec{e}_{\varphi}=
\left( \begin{array}{c}
u \sin\varphi\sin\varphi \\
-u \sin\varphi\cos\varphi \\
0
\end{array} \right)
\)
となります。よって、相対速度\(\vec{v}\)は
\(
\vec{v}=(\vec{u}\cdot \vec{e}_{\varphi})\vec{e}_{\varphi}-(\vec{r}\times \vec{\omega})=
\left( \begin{array}{c}
(u\sin\varphi-R\omega\sin\theta)\sin\varphi \\
-(u\sin\varphi-R\omega\sin\theta)\cos\varphi \\
0
\end{array} \right)
\)
と計算できます。
故に、\(\vec{r}\times\vec{v}\)は、
\(
\vec{r}\times\vec{v}=
\left( \begin{array}{c}
R\cos\theta(u\sin\varphi-R\omega\sin\theta)\cos\varphi \\
R\cos\theta(u\sin\varphi-R\omega\sin\theta)\sin\varphi \\
-R\sin\theta(u\sin\varphi-R\omega\sin\theta)
\end{array} \right)
\)
です。
具体的に、球体に働く力のモーメント\(\vec{N}\)を求めるためにやらなければならない積分は、
\(
\displaystyle \vec{N}=-\frac{\rho C_f}{2}R^2\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta|\vec{v}|(\vec{r}\times\vec{v})\sin\theta \\
\displaystyle =-\frac{\rho C_f}{2}R^2
\left( \begin{array}{c}
R\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin\theta\cos\theta\cos\varphi \\
R\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin\theta\cos\theta\sin\varphi \\
-R\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{array} \right)
\)
です。まず絶対値の符号を取るために場合分けを行います。
\(
|u\sin\varphi-R\omega\sin\theta| \\
=|R\omega|\cdot|\frac{u}{R\omega}\sin{\varphi}-\sin\theta|
\)
は、\(\frac{u}{R\omega}\sin{\varphi}\)に関して3つの領域に分けられます。

  1. \(\displaystyle \frac{u}{R\omega}\sin{\varphi}\le 0\)の領域
  2. \(\displaystyle 0\lt \frac{u}{R\omega}\sin{\varphi}\lt 1\)の領域
  3. \(\displaystyle 1\le \frac{u}{R\omega}\sin{\varphi}\)の領域

です。

\(0\lt \frac{u}{R\omega}\sin{\varphi}\lt 1\)が厄介なので特に考えましょう。
\(|\vec{v}|=|u\sin\varphi-R\omega\sin\theta|\)に対して、図のように考えます。
場合分け

例えば力のモーメントのx成分の積分は
場合分け_積分
となります。

こう表したのですが、x,y成分の被積分関数は\(\pi/2\)を中心に奇関数です。z成分だけは\(\pi/2\)を中心に偶関数となります。
なのでx,y成分は計算するまでもなくゼロになり、z成分だけ残ります。
これはもともとの回転角速度\(\omega\)がz成分しか持っていないから、という概要に合っています。

z成分の\(\theta\)に関する積分までは求めました。
シータ方向z成分_c
となります。
この後更に\(\varphi\)でこれを積分するのですが、僕はできませんでした。多分楕円積分に近いものが必要になるかと式の形から思われます。
\(\varphi\)方向に関しては数値計算でやってしまおうと思います。

2016/08/07追)
上記積分を近似して簡略化します。
回転している球体を上半分と下半分に分けて、球体を2つの平面に近似してしまいます。
上半分の速度\(v_{up}\)(風のベクトルと球体の回転方向が一致)と下半分の速度\(v_{down}\)(風のベクトルと球体の回転方向が反対方向)はそれぞれ
\(
\begin{align}
v_{up}&=u\sin_{\varphi_c}-R\omega\sin_{\theta_c} \\
v_{down}&=-u\sin_{\varphi_c}-R\omega\sin_{\theta_c}
\end{align}
\)
と書けます。ここで、\(\varphi_c,\theta_c\)は球を近似するに最適なある特定の角度を表します。
すると、近似した力のモーメントは
\(
\displaystyle N=\frac{\rho C_f}{2}R\frac{4\pi R^2}{2}\left(|v_{up}|v_{up}+|v_{down}|v_{down}\right)
\)

と導くことが出来ます。

実際の弾道のデータとの比較(回転の減衰考慮)


2015/06/07追記)
数値計算によって\(\varphi\)方向の数値計算を行い、比較してみます。
減衰ありの場合で実験と比較_c
回転数は、最大到達高さからおおよその値を導いています。
回転の減衰を考慮に入れた計算結果は緑の線です。
素晴らしいですね。最高です!
回転の減衰によってMagnus力が減少し、到達距離が短くなっています。
以前の計算結果である赤線の時に問題になった、最後のあたりのBB弾の伸びが違う、といった辺りがちゃんと改善されました。


0.20,0.25gの違いと回転減衰
なぜ質量が軽いほうが回転の減衰が早いのか?
今回の回転の減衰の導出では、表面積によってのみ、回転の減衰を引き起こす力のモーメントが変化します。なので、質量が重ければ思いほど、その分慣性モーメント\(I\)が大きいことになり、回転量は変化しにくいことになります。

スポンサーリンク


ちなみに、Magnus力のする仕事はゼロです。
なぜなら、Magnus力の働く方向は粒子の速度に垂直なはずで、方向は\(\vec{v}\times \vec{\omega}\)です。
その仕事は、係数を無視して、
\(
\begin{align}
&\int_{\vec{r}_{(t)}}^{\vec{r}_{(t+\Delta t)}}(\vec{v}\times\vec{\omega})\cdot d\vec{r} \\
&=\int_t^{t+\Delta t} (\vec{v}\times\vec{\omega})\cdot \vec{v} dt \\
&=\int_t^{t+\Delta t} 0{dt}
&=0
\end{align}
\)
となるからです。
・・・ということは、Magnus力が仕事をすることはあり得ません。

fortran90によるコードはこちらへ↓
弾道計算(BB弾)のコード

参考文献


[1]伊東敏雄著『な~るほど!の力学』学術図書出版社 (1994) p.220
[2]Faith A. Morrison, “Data Correlation for Drag Coefficient for Sphere,” Department of Chemical Engineering, Michigan Technological University, Houghton, MI (2015/02/08 アクセス確認)
[3]NASA, Ideal Lift of a Spinning Ball(2015/02/08 アクセス確認)
[4]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.378
[5]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.376

[6]第9章 境界層と物体まわりの流れ
[7]Skin friction drag

[8]流体力学講話・つまみ食い(その5)

[9]鳴尾 丈司, 溝田 武人, 下園 仁志, “一様気流中で高速回転するゴルフボールの空気力測定と飛しょう実験” 日本機械学会論文集 B編 Vol. 70 (2004) No. 697 P 2371-2377

[10]第3章 境界層理論, 3.5.4摩擦抵抗係数 \(C_{D_f}\)

[11]石岡、ホップアップの回転数測定 TM VSR-10 G-SPEC

[12]放物運動 -初等力学の未解決問題:速度の2乗に比例する空気抵抗を持つ放物運動