弾道計算(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乗に比例する空気抵抗を持つ放物運動


コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です