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

弾道計算, 物理学

エアガンの集弾限界

コメントする

理想的なエアガンとBB弾、それに風などが無い理想的な環境があったとします。

バレルとBB弾に直径差があることによって生じるばらつきは、

　30m先 | 直径3cm
　50m先 | 直径10cm

となります。この直径以下にばらつきを抑えることは構造的に不可能です（60cmのバレルで0.9Jで射出する場合）。
つまり、このばらつきよりも大きくばらつく場合、バレルとBB弾の直径差以外にばらつく原因が存在します。

バレルとBB弾に直径差があることによって生じるこのばらつきを軽減するには、
　・バレルの直径を小さくする
　・バレルを長くする
　・バレル内部の素材とBB弾の反発が起こりにくいものにする
にすることで軽減されますが、他の要因が場合はこの限りではありません。

ここでの理想的とは、

・球の直径のばらつきは無い
・いつも同じ初速で撃ち出せる
・無風

とします。この場合でも同じ点に集弾することはありません。
なぜなら、バレルの直径とBB弾の直径が異なる為です。この差によってどの位集弾性が悪くなるのか、見積もりましょう。

これは構造的な問題であり、ばらつきの原因の中で取り除く事ができない１つの原因です。

ばらつく原因として、以下の三つが考えられますが、まず本稿では理想状態のばらつきを考えます。

ばらつきの種類 理想状態のばらつき（本稿） 製品誤差によるばらつき セッティングによるばらつき
BB弾重さ – あり なし
BB弾の大きさ なし あり なし
回転のばらつき なし なし あり
手振れ なし なし あり

ばらつきが生じる原因

何故ばらつきが生じるかを考えましょう。
全くの理想であれば、銃口から同じ初速、角度、回転量で射出されたBB弾にばらつきが生じることはありません。
しかし、現実にはバレルの大きさとBB弾の大きさに差があります。これによって銃口から射出する時に進行方向と垂直な平面に速度を持つことになります。
銃口から出てきたBB弾の速度\(\mathbf{v}\)を以下のように表現します。

ここで、\(v_{z}, \mathbf{e}_{z}\)はバレルの方向に沿う速度、単位ベクトルで、
\(v_{x,y}, \mathbf{e}_{x,y}\)はバレルの方向に垂直な面への速度、単位ベクトルです。
すなわち、壁面に垂直な方向に対する速度ベクトルの大きさ\(v_\perp\)は、\(v_\perp=\sqrt{v_x^2+v_y^2}\)と表されます。

ばらつきが無くまっすぐ飛ぶのであれば\(v_{\perp}=0\)であり、そうでなければ\(v_{\perp}\ne 0\)です。

BB弾のばらつきが生じる原因は、射出方向に垂直な方向（横方向）に有限の速度が生じている、と仮定します。
この横方向の速度が生じる原因の一つは、ピストンでBB弾を空気で押す際に圧力が一定ではないとか、回転を掛けるためのゴムで生じる、などいろいろ考えられます。
ここでは、パッキンを通過して、もうこれ以上横揺れを増やすような原因が生じないのだ、と仮定します。壁の反射によって変化はあるものの、速度は減衰する一方であるとします。

それでは、定式化をしていきましょう。

定式化

横方向の速度はそれほど大きくないだろう、と予想するので空気抵抗は考えません。この条件のもと考えていきます。

バレルの直径\(d_s\)とBB弾の直径\(d_B\)の差があることによって、バレルに垂直な方向にBB弾が自由に進める距離\(l\)は

と書けます。\(n\)回目の衝突から\(n+1\)回目の衝突までにかかる時間\(t_n\)は、その間のBB弾の垂直方向の速度\(v_n\)に依存して、

と書けます。よって、\(N\)回衝突するのにかかる時間\(t\)はそれぞれの衝突までに掛かる時間を足し合わせればよいので、

です。バレルの内壁との衝突が起こることによってBB弾の速度の変化するとします。するとBB弾とバレルとの反発係数\(r\)を用いれば、

と表現されます。初速度を\(v_0\)と書いてこれを代入すれば、

と書くことが出来ます。ただし、\(r\ne 1\)です。\(r=1\)の場合は\(\displaystyle
t=\frac{l}{v_0}(N+1)
\)
です。\(r\ne 1\)の場合に\(N\)について変形すれば、

と書くことが出来ます。
初速\(v_0\)、反発係数\(r\), 銃口に到達するまでに掛かる時間\(t\)が分かれば、銃口から出た時のバレルに垂直な方向の速度\(v_\perp\)は式(5b)より、

と書けます。具体的に妥当な値を入れて射出時の、進行方向に垂直な方向の速度を計算してみましょう。

具体的な衝突回数の見積もり

60cmのバレル長を持ち、0.20gで0.9J程度の球が射出される場合、セットされた位置から射出までにかかる時間は約\(t=0.015\mathrm{[s]}\)と分かっていますのでこれを利用します（バレル内部の計算より）。
バレルBB弾の直径差は

とします。\(r\)はどうにかして求めることにして、\(v_0\)はとりあえず変数としましょう。

反発係数の見積もり

\(r\)を見積もります。
反発係数の見積もりは、実験して大雑把な値を見積もります。
\(r\)を見積もるために、家にある固い材質のものとBB弾とを衝突させて、高さを計測しました。以下の実測結果を得ました。

  • BB弾-フローリング 30cmから落として15cm上がる
  • BB弾-アクリル板 30cmから落として10cm上がる
  • BB弾-アルミ 10cmから落として6cmまで上がる

高さ\(h_0\)から静かに落として、反発した後に極大値をとるときの高さ\(h_1\)が分かっているとき、反発係数\(r\)は

と求められることが分かっているので、反発係数は

  • BB弾-フローリング　　\(r\approx 0.71\)
  • BB弾-アクリル板　　\(r\approx 0.57\)
  • BB弾-アルミ　　\(r\approx 0.77\)

となります。バレル内部はどちらかといえば金属に近いと思いますので、\(r\approx 0.75\)と仮定して計算を進めていきます。

横軸に進行方向に垂直な方向の初速度…つまり、ホップを掛けるためのゴムなどで、進行方向ではない方向の速度のことを意味します…、縦軸に射出時の垂直な方向の速度をプロットしました。横方向の速度が仮に50m/sになっていても、10m/sになっていようともあまり変わらないことが分かります。なので、10m/sと仮定しましょう。
ばらつきの上限を与えるにはよい指標です（※1）。

仮に\(r=0.75, v_0=10\mathrm{[m/s]}\)とすると、\(N\approx 20.6\)と計算されます。つまり、20.6回バレルに衝突してから飛び出していくわけです。実際には衝突回数は整数しかありえないので、小数点以下を切り捨てて20回衝突が起こって飛び出していきます。
ただし計算を行う上では衝突回数を切り捨てると不連続性が発生しますので、小数点の衝突回数を認めることにします。

この場合、銃口から飛び出す際の、進行方向に垂直な方向の速度\(v_\perp\)は

を得ます。この垂直方向の速度は上下左右に振れる可能性がありますが、最も触れる場合は横方向でしょう。

30, 50m先の広がり具合を考えます。
着弾までの時間\(t_{30m}, t_{50m}\)はそれぞれ\(0.5, 2 [s]\)分かっているので（弾道計算（BB弾）の結果）、着弾時の直径\(d_{30m}, d_{50m}\)は

と求められます。ここで垂直方向の速度が非常に遅いので、空気抵抗は無視して考えています。BB弾の重さは関係ありません。軽いものほど空気抵抗を受けて減衰するので、ぶれは小さくなります。

空気抵抗を考慮すると、これよりは小さくなるということです。

実測定との比較

さて、30mチャレンジというものがあります30mチャレンジ公式ランキング　2016年度。
この記録によると、30m先で直径15cm程度の範囲に収まるようです。
すなわち、バレル-BB弾の直径差による集弾性の悪化よりも大きな影響を及ぼす要因がある、ということです。

恐らくは回転量が一定ではないとか、BB弾の重さにばらつきがあったり、銃身の振動があったりという要因のほうが大きいということでしょう。

逆に言えば、スタンダードな大きさのバレルとBB弾を用いる限り、30m先で3cm以内に収めることは不可能、ということです。

注釈

※1
\(v_0\to \infty\)で射出時の速度\(v_\perp\)は発散しますが対数の発散です。非常にゆっくり発散するため、\(v_0=10\)であろうが\(v_0=50\)であろうがほとんど変化しません。どこかでカットオフ（これ以上はとらない上限）となるような値を取れればよいかもしれません。
豆知識ですが、対数の発散をほとんど無視する目的としてカットオフを設けるやり方は、量子力学の繰り込み論が有名ですね。だからと言って、エアガンの計算で量子力学が現れる！などは言わないでください。全く関係ありません。

弾道計算

流速チューン？設定におけるホップ量と初速の変化

コメントする

流速チューンと呼ばれるエアガンのチューニングがあるそうですね。

流速チューンとは
　・重いピストンヘッドを使う
　・強いばねを使う
　・短いバレルにする
というチューニングを行う事のようです[1]。この結果、
　・ホップを掛けるほど初速が上がる
　・破裂音が鳴る
という効果が得られるらしいです。これらの効果を物理的な言葉に置き換えます。

「ホップを掛ける」ということは、BB弾がホップを掛ける為のゴムの所で止まりそこを抜けるために力が必要、ということです。
すなわち、内圧がある程度高まらないと動き出さないということを意味します。
「破裂音が鳴る」ということは、射出時に外気圧と内圧の差が大きいことを意味します。
また、短いバレルにするのは初速が上がりすぎるためだそうです。

以上の振る舞いがシミュレーションで再現できるでしょうか。
グラフの横軸を射出に必要な内圧にとり、バレル長を10cmに設定した時の射出速度と射出時の内圧をプロットしたのが以下の図です。

確かによく言われている通り、BB弾が動き出すまでに力が必要なほど（ホップを掛けるほど）、初速が上がる振る舞いが再現できています。
ある程度のホップの強さからは頭打ちになるようです。このある程度の強さが実際のエアガンで掛けられる強さなのかはわかりません。
射出時の内圧も3気圧と非常に高くなります。外気圧は1気圧なので、これだけ差があれば破裂音が鳴るのも納得でしょう。

流速チューン？のパラメータ設定は以下の通りです。通常の設定と比べ、ピストンヘッド+ばねの重さを22→30g, ばね定数を431→800kg s-2に変更しました。
————————-
BB弾直径 5.95mm
BB弾重さ 0.25g
バレルの直径 6.05mm
バレル長 10cm
シリンダー直径 24mm
ピストンヘッド+ばねの重さ 30g
ばね自然長 150mm
ばね定数 800kg s-2
————————-
流速チューンだ!と言っている方々のブログ等を拝見しましたが、なかなかばねの強さやピストンヘッドの重さ等を詳細に詳しく述べているページはありませんでしたので、上記パラメータはなんとなく設定しています。
また、流速チューン？とはてなマークを付けているのは、具体的なパラメータ例が無いことや、定義が個人によって違うからです。

追記）
通常のばねや、ピストンヘッドの重さでもホップを掛けるほど初速が上がるという傾向は変わらないことが分かりました。
すなわち、上記の傾向はバレル長が短いことによる普遍的な傾向だ、ということです。バレルが短ければホップを掛けるほど初速が上がる、という意味です。

上の計算は以下のパラメータで計算を行いました。
————————-
BB弾直径 5.95mm
BB弾重さ 0.25g
バレルの直径 6.05mm
バレル長 10cm
シリンダー直径 24mm
ピストンヘッド+ばねの重さ 22g
ばね自然長 150mm
ばね定数 431kg s-2
————————-

以上の計算は、下記ページで配布しているエクセルシート(ver1.3以降)にて計算することが出来ます。
弾道計算のコード(Excel)

参考資料

[1]流速チューンの原理と製作例、注意点

プログラミングと数値計算, 弾道計算

弾道計算のコード(Excel)

コメントする

Excelのマクロ(VBA)で弾道計算を行うコードを作成しました。

https://slpr.sakura.ne.jp/qp/supplement_data/BBullet_course_excel/BBbullet_course_ver1.41.xlsm
上記プログラムの利用の際は連絡先と著作権についてをご覧ください。

空気抵抗を修正したβ版を作成しました。参考設定は機能しません。
https://slpr.sakura.ne.jp/qp/supplement_data/BBullet_course_excel/BBbullet_course_ver1.5_beta.xlsm

以下の微分方程式を計算します。

エクセルシート”bullet”
\(
\begin{align}
m\frac{d^2 \vec{r}}{dt^2}&=
-mg\vec{k}-\frac{1}{2}C_d \rho \pi R^2 |\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}|} \\
\frac{d\mathbf{\omega}}{dt}&=N/I \\
N&=\frac{\rho C_f}{2}R\frac{4\pi R^2}{2}\left(|v_{up}|v_{up}+|v_{down}|v_{down}\right)
\end{align}
\)
回転の減衰は本来、弾道計算(BB弾)の理論で導いた通り積分を行わなければなりませんが、簡単にするためにBB弾を真四角の箱として近似して導いています。
この近似は本来の積分値と1パーセント程度の違いしか出てこないので、BB弾の計算には非常に良い近似です。
回転の減衰を取り入れたくない場合、入力パラメータの”回転の減衰”項目をNOにすることでそのように設定が出来ます。
ゼロインを固定した時の弾道の探索機能は実装していないので出来ません。
Excelの設定項目では、Fortranで実装したコードから求めたパラメータをあらかじめ用意しておきました。
設定されたものと書いてあっても、高精度に積分計算をして求めた結果と上の近似をもちいた本エクセルの結果は違くなります。数10cm程度のずれは出てきますので、その点は注意してください。

エクセルシート”barrel”
\(
\begin{align}
m\frac{d^2}{dt^2}x(t)&=\left[P(t)-P_0\right]S_{BB}-\frac{1}{2}C_d \rho \pi R^2 |v(t)|v(t) \\
m_s\frac{d^2}{dt^2}x_0(t)&=-k\left[x_0(t)-x_B-l\right]-[P(t)-P_0]S_s -F_f\label{bbin2} \\
\frac{d}{dt}V(t)&=v(t)S_{\text{barrel}}-v_0(t)S_s+c'(S_{barrel}-\pi R^2)\mathrm{sgn}(P(t)-P_0)\sqrt{\frac{2}{\rho}|P(t)-P_0|}
\end{align}
\)

内部は以下の構造をしていると想定しています。詳細はバレル内部でのBB弾の方程式をご覧ください。

計算手法はどちらも陽的4次ルンゲクッタ法を用いています。

使い方

  1. ファイルのダウンロード

    2. 最新
    https://slpr.sakura.ne.jp/qp/supplement_data/BBullet_course_excel/BBbullet_course_ver1.5_beta.xlsm

    をクリックでダウンロード、もしくは
    右クリックして、”名前を付けてリンク先を保存”からダウンロードしてください。

  2. 計算パラメータの入力

    3. 入力する場所は、色付きのセルです。

    • 弾道計算について


      • 出力結果は大きく分けて２つの項目に分かれています。
      １つ目は定数、射出高さに着弾、地面に着弾の項目は重要と思われるパラメータの抜粋、２つ目は弾道の時間ごとのデータです。
      グラフは、２つ目のデータを元に表示しているだけですので、必要に応じて追加、変更を行ってください。

    • バレル内部の計算について


      • 出力結果は大きく分けて２つの項目に分かれています。
      １つ目は定数、内圧が大気圧に等しい時に射出される時と指定したバレル長におけるパラメータの抜粋、２つ目はBB弾、ピストンの位置、速度の時間ごとのデータです。
      グラフは、２つ目のデータを元に表示しているだけですので、必要に応じて追加、変更を行ってください。

    手入力で入力することも可能ですが、参考設定と書かれたボタンを押すことで、予め設定されたパラメータを読み込ませることもできます。

  3. 計算開始

    4. 　　入力パラメータを入力した後、計算開始ボタンを押すことで計算が開始されます。

注意点

常識的なパラメータの範囲内でしか動作確認をしておりません。
この結果はあくまで一例です。本シミュレーションを信用して設計・組立を行い、本シミュレーションと結果が違う事があります。その場合私は責任を負いませんので、その点をご了承して使用してください。
BB弾で出来るからと言って、銃弾等の計算はそれなりの信用しかできないことを考慮してください。

追記

バレル内部の計算において、抜弾抵抗の考慮はしていません。すなわちホップが全く掛かっていない状況の計算です。
もし、ホップを掛けた状況を再現したいのであれば、以下のような考えで出来るでしょう。

まず、ホップを掛けた状態でBB弾をバレル内部から引き出す際に必要な力をどうにかして実験的に求めます。
このホップを掛けた際に必要な力\(F\)は、内部圧力\(P\)とBB弾の面積\(S_{BB}\)を用いて
\(
F=PS_{BB}
\)
と書けるはずです。すなわち、内圧がこの力に達した時にBB弾が射出されるとして考え、その内部圧力に達した際にBB弾の運動方程式を解き始める、という方法で可能なはずです。

参考文献

[1] 弾道計算の理論 -シキノート https://slpr.sakura.ne.jp/qp/bullet-course/
[2] BB弾の弾道学 -ハイパー道楽 http://www.hyperdouraku.com/colum/ballistics/index.html

弾道計算

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

コメントする

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

まとめ

目的

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

方法

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

結論

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

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

[adsense1]

きっかけ

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

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

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

方針

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

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

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

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

理論上の問題

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

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

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

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

実験上の問題

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

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

[adsense2]

実験と理論の比較

理論上の式を数値計算で解きます。数値計算で用いる数式は弾道計算(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弾の回転量について(実験との比較)

2件のコメント

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

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

[adsense1]

対象とする問題

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\)が良いと思います。

[adsense2]

結論

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

・BB弾の運動方程式
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\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)\)として考えます。
１．どんな回転量でも\(C_l\)は正である必要があります(負のマグヌス力は今回の場合発生しないと考える)。なので、
\(
C_l(\omega) > 0
\)
でなければなりません。
２．回転が非常に0に近い時、球体の循環と周りで発生する循環は一致するはずです。
なので、
\(
C_l(\omega\to 0)=1
\)
であることが予想されます。
３．回転が非常に強い場合、周辺の空気は限られるため発生する揚力は弱くなっていくはずです。
ただし、回転が無限大の時、揚力自体は発散すると予想できるので、係数\(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)
弾道計算のコード(Excel)
バレル内部でのBB弾の方程式←今ここ

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

空気の漏れも考えます。

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

内部圧力変化と力

ピストン-BB弾間の圧力変化による力\(F_V(t)\)は、
\(
F_V(t)=S\cdot P(t)
\)
と書けます。ここで\(P(t)\)はピストン-BB弾間の内部圧力です。
\(P(t)\)は、断熱変化を仮定すると、ピストン-BB弾間体積\(V(t)\)を用いて
\(
P(t)V^{7/5}(t)=\text{const}
\)
の関係があります。体積\(V(t)\)は
\(
V(t)=x(t)S_b-x_0(t)S_s
\)
と書けます。ここで、ピストンの位置による時間変化を\(x_0(t)\)、BB弾の位置を\(x(t)\)、バレルの断面積を\(S_b\)、シリンダーの断面積を\(S_s\)としました。
この式は時刻に依存せずに決まるので、時刻\(t=0\)の時に大気圧\(P_0\)で、体積\(V_0\)であれば、任意の時刻での圧力\(P(t)\)は
\(
\displaystyle P(t)=P_0 \left(\frac{V_0}{V(t)}\right)^{7/5}
\)
と書けます。

空気の漏れについて

十分短い時間の間、時刻\(t\)において内部圧力\(P(t)\), 空気の密度\(\rho\)とすると、漏れ出る空気の速度\(v(t)\)はベルヌーイの定理から
\(
\begin{align}
P(t)&=\frac{1}{2}\rho v(t)^2+P_0 \nonumber \\
&\to v(t)=\pm \sqrt{\frac{2}{\rho}(P(t)-P_0)}
\end{align}
\)
が成立するとします。厳密には、ピストンの速度は漏れ出る空気の速度に対して無視できるほど小さい、という仮定の下で成立します。漏れ出る流量\(Q\)は、バレル-BB弾の隙間の断面積\(A\)、実験とのズレを調節する無次元の流量係数\(c’\)を用いて\(Q=c’ A v(t)\)と書けます。
また、断熱過程を圧力-体積の関係を用いたいので、漏れ出る空気はBB弾-ピストン間の体積の増加として扱います。

フルシリンダーの場合の運動方程式

その他、空気抵抗による力を入れると、BB弾の位置\(x(t)\), ピストンの位置\(x_0(t)\), 体積\(V(t)\)の運動方程式は
\(
\begin{align}
m\frac{d^2}{dt^2}x(t)&=\left[P(t)-P_0\right]S_b-\frac{1}{2}C_d \rho \pi R^2 |v(t)|^2\cdot\frac{v(t)}{|v(t)|} \label{bbin1}\\
m_s\frac{d^2}{dt^2}x_0(t)&=-k\left[x_0(t)-x_B-l\right]-[P(t)-P_0]S_s -F_f\label{bbin2} \\
\frac{d}{dt}V(t)&=v(t)S_b-v_0(t)S_s+c'(S_b-\pi R^2)\mathrm{sgn}(P(t)-P_0)\sqrt{\frac{2}{\rho}|P(t)-P_0|} \label{bbin3}
\end{align}
\)
と導くことが出来ます。

ここで、\(P_0\)は大気圧、\(\eta\)は粘性率、\(R\)はBB弾の半径、\(C_d\)は抗力係数、\(\rho\)は空気の密度, \(S_b\)はバレルの断面積、\(S_b\)はシリンダーの断面積、\(k\)はピストンのばね定数、\(l\)はばねの自然長、\(F_f\)はピストン-シリンダーの摩擦、\(v(t),v_0(t)\)はそれぞれBB弾の速度、ピストンの速度を意味します。

空気抵抗に関する詳細は球体の空気抵抗と係数をご覧ください。

ここで１つ言えることは、外部と内部の圧力が一定になる最適なバレル長というものが存在する、ということです。

2019/01/20 追記)
本の頒布日より1年以上たちました。
計算結果は載せないつもりでしたが、載せたくなりましたので少しだけ載せておきます。
まず、実測データと上で立てたモデルの運動方程式の比較をします。
比較対象はインターネット上で公開されていた二つのデータです。

実測データとの比較


1つ目(図の←)はさばなび様で公開されていた記事
https://www.saba-navi.com/2015/10/29/laboratory_work_barrel_cut/で実測されていたデータで、バレル長の長さと初速の関係です。ただし、既にリンク切れのようです。
ピストンの重さの記述が無かったので、典型的とされている重さ25gを仮定して計算しています。

2つ目(図の→)は石岡様のホームページ
○電動ガンバレル、シリンダの組み合わせによる初速実験
の結果との比較です。

どちらの計算結果もまぁまぁ合っていることから、私が仮定したモデルは見当はずれなものではないということが分かるかと思います。

バレル内部の様子

さて、上のモデルが正しいとした時、運動方程式を解いてバレルの動きや内圧を考えてみましょう。
解いてみますと、こういった図が得られます。

図は、東京マルイM16を想定した時のBB弾の位置、速度、ピストンの位置、内圧の時間変化です。
計算のパラメータは、
BB弾の重さ0.20g
ピストンの重さ24g
ばね定数431N/m
ばねの自然長150mm
押し切られた時のばねの長さ120mm
BB弾の半径5.95mm
バレルの直径6.05mm
とした時の結果です。
この場合、生み出せる最高初速は96m/s（バレル長63cmの時）。
例えばバレル長が50cmの時、ピストンが動いてからBB弾が射出されるまでにかかる時間は0.013秒ということが分かります。
すなわち、上のピストンの重さやBB弾の重さ、ばねの強さ、シリンダー容量の場合は最適なバレル長は63cmであるということです。
また、内圧の上昇によるピストンのブレーキなども計算できていることが分かるかと思います。

また、漏れ出る空気の量は\(3000[\mathrm{mm^3}]\)だということが見積もられました。この量はM16の場合は本来のシリンダー容量の15%程度です。

典型的なバレル長である50cmの場合、BB弾は射出されるまでに約0.013秒かかります。内圧の上昇によるブレーキがかからず、同じ過程でばねが引かれると仮定すると、ピストンは約0.010秒で戻ると考えられます。
理論上の効率的な射出サイクルの限界は

（1発）÷（0.013秒＋0.010秒）≒ 43発/秒

という事です。現在販売されている東京マルイのハイサイクルは25発/秒なので、理論上あと1.7倍早くすることが可能です。

BB弾の重さを0.25gにした結果や、ばねの強さを変えた時の結果は
弾道計算本として計算結果をまとめた本に載せているので、そちらをご参照ください。本の詳細については弾道計算本の自家通販をご覧ください。
また時間が経ったらここに載せる…かもしれません。

補遺

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

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

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

参考文献

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

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

弾道計算, 物理学

球体の抗力係数

1件のコメント

説明

本稿では、完全な球体に働く抗力係数について述べます。
ここでは概要だけを説明します。詳細は以下のPDFをご覧ください。

https://slpr.sakura.ne.jp/qp/supplement_data/drag_coefficient_air/drag_coefficient.pdf

概要

速度\(v(=|\mathbf{v}|)\)で動く半径\(R\)の完全な球体に働く空気抵抗力の大きさ\(F_d(=|\mathbf{F}_d|)\)は,
\begin{align}
F_d =\frac{1}{2} C_d(R_e) \rho \pi R^2 v^2
\end{align}
と書けます. 空気抵抗力\(\mathbf{F}_d\)は, 球体の速度方向と反対に働くので, 方向を含めれば
\begin{align}
\mathbf{F}_d =-\frac{1}{2} C_d(R_e) \rho \pi R^2 v^2\frac{\mathbf{v}}{v}
\end{align}
と書けます.

つまり、重力加速度\(g\)の重力下における質量\(m\)の完全な球体に対する運動方程式は
\begin{align}
m\frac{d^2}{dt^2}\mathbf{r} =-mg\mathbf{k}-\frac{1}{2} C_d(R_e) \rho \pi R^2 v^2\frac{\mathbf{v}}{v}
\end{align}
となります。ここで、\(\mathbf{k}\)は重力の方向の単位ベクトル（多くの場合では鉛直上向きを\(z\)軸の正の方向、単位ベクトル\(\mathbf{k}\)とするので、負号が付いています）です。
速度が非常に遅い極限\((R_e\to 0)\)においては、\( C_d(R_e)|_{R_e\to 0}\to \frac{24}{R_e}\)となります。

完全な球体の抗力係数\(C_d\)における近似式は、レイノルズ数を\(R_e\)として以下の通り与えられます。

  • \(R_e\lt 2\times 10^5\)層流領域のみ
    \begin{align}
    C_d(R_e)&=\left[\frac{24}{R_e}(1+0.15 R_e^{0.687})\right]+\frac{0.42}{1+(\frac{42500}{R_e^{1.16}})},~~~\mbox{式(1)}
    \end{align}
    ただし、\(R_e \lt 2\times 10^5~~\mbox{and}~~K_n \lt 0.01 ~~\mbox{and}~~ M_a\to 0 \)
  • \(R_e\lt 2\times 10^5\)の層流領域のみ、マッハ数\(M_a\)依存性込み
    \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}}} ,~~~\mbox{式(2)}
    \end{align}
    ここで、\(C_M, G_M, H_M\)は以下の通り与えられます。
    \(\displaystyle
    \begin{eqnarray}
    C_M=\left\{~~
    \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}
    \right.
    \end{eqnarray}
    \)

  • \(\displaystyle
    \begin{eqnarray}
    G_M=\left\{~~
    \begin{aligned}
    &166M_a^3+3.29M_a^2-10.9M_a+20~~~\mbox{for} ~~~M_a\lt 0.8\\
    &5+40M_a^{-3}~~~\hspace{9.4em}\mbox{for} ~~~ M_a\gt 0.8
    \end{aligned}
    \right.
    \end{eqnarray}
    \)

    \(\displaystyle
    \begin{eqnarray}
    H_M=\left\{~~
    \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}~~~\hspace{8.8em}\mbox{for} ~~~ M_a\gt 1
    \end{aligned}
    \right.
    \end{eqnarray}
    \)

  • 定性的、広範囲
    \begin{align}
    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{0.25\frac{R_e}{10^6}}{1+\frac{R_e}{10^6}}\right),~~~\mbox{式(3)}
    \end{align}

対象物が流体中を高速で動かなければ、層流領域で実験値を良く合う式(1)（図中の赤線）, 式(2)で十分でしょう。
温度なども含めたいならば、PDFで説明しているレイノルズ数、マッハ数を置き換えて利用することを勧めます。
また精度はあまりいらない状況で、定性的な理解だけを対象とするならば、乱流領域の振る舞いを含む式(3)（図中の水色線）を用いると良いです。

以降、詳細は
https://slpr.sakura.ne.jp/qp/supplement_data/drag_coefficient_air/drag_coefficient.pdf
をご覧ください。

サバイバルゲーム, 弾道計算, 物理学

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

9件のコメント

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

この結果に付随して、
0.25gで0.8Jの軌道は
0.20gで1.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.20gの1.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乗に比例して強くなるため、それに見合うようなの飛距離の伸びは得られません。

G&G バイオBB弾 0.25g 1KG 4000発入り 【G07126】

東京マルイ 0.25g　ベアリングバイオBB弾 生分解バイオ 1300発入り【３袋セット】

補足

「流速チューン」というものがあるそうですね。
流速チューンの問題点と改善案について -週休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弾)のコード

17件のコメント

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)←今ここ
弾道計算のコード(Excel)
バレル内部での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
を参照しました。

改良版

これから最適な回転量、角度を探すプログラムに更新したいと思いますが、とりあえず
改良したプログラムを載せておきます。
input
main.f90