第8回 軌道シミュレータver.3.2 -重力・揚力・抗力全方向-

エクセルで作った軌道シミュレータver.3.1を全方向速度に対する空気力へ拡張しver.3.2にアップデートします。

数式は、ver3.0の抗力、ver3.1の揚力と同じ形の式をy,z方向速度についても拡張したものになります。

では式を立てます。少し長くなります。


軌道計算の運動方程式

ボールに作用する力を重量で割ると加速度が求まります。その加速度を時間積分すると速度になり、速度を時間積分するとボールの位置が求まります。

ボールに作用する力は「重力」、「抗力」、「揚力」の３つです。

[計算式]

x軸をホーム方向、y軸を一塁方向、z軸を上空方向と定義。

ボール回転軸は、z軸からx-y平面に向かう角度（真上から水平に向かう角度）をθs、x軸からy軸に向かう角度（ホーム方向から一塁側へ向かう角度）をΦsと定義。
今回は4シームとしてθs=110度、Φs=-80度とします。

回転軸により揚力を各方向成分へ振り分けます。

抗力係数は暫定的に回転軸の向きによらず一定とします。

　空気力による各方向への加速度

　各方向の加速度を行列式でまとめて表すと以下のようになります。

　各速度成分により発生する各方向への抗力、揚力をボールの重量で割ったものです。

　

　ここで、sgn(dx/dt)はx方向速度(dx/dt)の向きを表す符号で、dx/dt>0（投手から
　ホーム方向へ飛翔）ならsgn(dx/dt)=1, dx/dt<0(ホームから投手方向飛翔)なら
　sgn=-1です。
　二乗すると速度のプラスマイナスが分からなくなってしまうため、これを付けて区別
　します。
　例えば投手が投げるトップスピンのカーブと捕手が投げ返すバックスピンの球は、
　三塁側から見ればどちらも同じ右回りの回転です。
　しかしx方向速度が反対のため前者は下に変化、後者は上にホップと揚力の方向が反
　対になります。
　これを正しく再現するために速度のプラスマイナスの区別が必要となります。

　Sgn(dy/dt), sgn(dz/dt)も同様です。

x方向

上記行列式のx方向加速度を通常の式の形に書き下すと以下のようになります。

y方向

上記行列式のy方向加速度を通常の式の形に書き下すと以下のようになります。

z方向

上記行列式のx方向加速度を通常の式の形に書き下し、重力加速度を追加すると以下のようになります。

ここで、

v0:リリース時の球速、x0,y0,z0:リリース位置、
θ:上向きリリース角度、Φ:横向きリリース角度、t:リリース後経過時間

[計算式おわり]

ボールの空間位置x,y,zはtを代入すれば求まる形の式にできないため、エクセルで数値積分の方法で求めます。

具体的には、ti+1秒後の速度=ti秒後の速度 + ti秒後の加速度×(ti+1 - ti)、
ti+1秒後の位置=ti秒後の位置 + ti秒後の速度×(ti+1 - ti)、といった具合です。



エクセルへの数式入力

では、エクセルへ入力していきます。

単位系はSI単位系です。

 [エクセル入力]

まず初期条件を入力。

ver3.1と同様にしました。

ボールの回転軸角度は完全なバックスピンではなく、ジャイロ成分やサイドスピン(シュート)成分が混じっていることを考慮してθs=110度,φs=-80度とします。

揚力係数CD、抗力係数CDの値は個人で計算や実測をするのが難しいのでネットで拾った値を使います。




次に時間tを少しずつ増加させて複数入力します。
ver.3.1と同じく、0.02秒おきにしました。

数値積分をするので時間間隔が小さいほど計算結果は正確になります。
そして時間t列のとなりに、d2x/dt、dx/dt、x、d2y/dt、dy/dt、y、d2z/dt、dz/dt、z用の列を9列挿入します。





X方向

dx/dt列の一番上のセルに、時間t=0における速度を入力します。(④式)



x列の一番上のセルに、時間t=0におけるx位置を入力します。(⑤式)



そしてd2x/dt2列に、dx/dtとdy/dtとdz/dtからd2x/dt2を計算する式を入力します。(①式)

再びdx/dt列に戻り、t>0のセルに数式を入力します。
②式の計算を数値積分で行うために下図のように数式を入力します。



再びx列に戻り、xのt>0について③式の計算を数値積分で行うために下図のように数式を入力します。


x位置の計算式はこれで完成です。


Ｙ方向

dy/dt列の一番上のセルに、時間t=0における速度を入力します。(⑨式)



y列の一番上のセルに、時間t=0におけるy位置を入力します。(⑩式)



そしてd2y/dt2列に、dx/dtとdy/dtとdz/dtからd2y/dt2を計算する式を入力します。(⑥式)



再びdy/dt列に戻り、t>0のセルに数式を入力します。
⑦式の計算を数値積分で行うために下図のように数式を入力します。



再びy列に戻り、yのt>0について⑧式の計算を数値積分で行うために下図のように数式を入力します。



y位置の計算式はこれで完成です。


Z方向

dz/dt列の一番上のセルに、時間t=0における速度を入力します。(⑮式)



z列の一番上のセルに、時間t=0におけるz位置を入力します。(⑮式)



そしてd2z/dt2列に、dx/dtとdy/dtとdz/dtからd2yzdt2を計算する式を入力します。
重力加速度gも忘れないよう注意です。(⑪式)



再びdz/dt列に戻り、t>0のセルに数式を入力します。
⑫式の計算を数値積分で行うために下図のように数式を入力します。



再びz列に戻り、zのt>0について⑬式の計算を数値積分で行うために下図のように数式を入力します。z位置の計算式はこれで完成です。





完成です!




これで各時刻におけるボール位置x,y,zの値が計算できました。

[エクセル入力おわり]

軌道計算結果のプロット

では、計算されたx,y,z値を散布図でグラフ上にプロットします。

[エクセルグラフ化]

空気力(揚力、抗力)を全方向速度に対して完全に考慮した場合の軌道計算結果は、以下のようになりました。
グラフ中の点は0.02秒ごとのボール位置を表します。

空気力をx方向速度成分に対してのみ考慮し、y,ｚ方向速度に対する空気力はなしとして計算しているver.3.1の軌道計算結果と重ねてみました。



しかし差が小さすぎてわかりません。完全に線が重なっており見分けがつきません。

ホームベース上における両者の軌道の差はわずか4mmです。

[エクセルグラフ化おわり]

ピッチングではx方向の速度成分が大部分であり、y,z方向速度成分は小さいためそれらによる空気力はあってもなくてもほぼ変わりません。

ver.3.1のようにy,z方向の速度成分による空気力は省略しても、通常のピッチングにおける軌道計算をするなら十分な精度が得られるということです。

ここで行っている計算はエクセルで一瞬で終わるようなものなのでどちらでもよいのですが、サーバーでの計算時間が数百時間に及ぶような大規模シミュレーションではあってもなくても大差がない部分を上手く省略し効率化を図ることは重要です。


遠投軌道計算結果による影響比較

y,z方向速度成分による空気力あり、なしの場合の差が明確になるよう、遠投の軌道を計算してみます。
遠投では、ピッチングに比べボール軌道が山なりとなりz方向速度成分が大きくなりそれによる空気量の影響も無視できない程度に大きくなること、また飛翔時間が長いため空気力を受ける時間が長くなることにより、両者の軌道の差が顕著になります。

[計算結果2]

リリース角度を上向き35度に変更します。
それ以外の条件はピッチングの時と同様で球速140km/h、抗力係数CD=0.33、
揚力係数CL=0.20とします。

計算された両者の軌道はこのようになりました。

ver.3.2ではボールが上昇していくときz方向速度に対しても抗力が働きより減速が大きいため、ver.3.1に比べ頂点の高さが6m程減少しています。
100m過ぎあたりのボールが落下していく段階では、ver.3.2方はz方向速度に対して揚力が働きまたz方向への落下速度に対して抗力のブレーキが働くため、ver3.1に比べるとやや水平に近い軌道になっています。

全体としては軌道前半の抗力による影響の方が大きく、ver3.2はver.3.1に比べ飛距離が10m減少しています。

(*) ver.3.2において140km/hでの飛距離115mは、実際と比べるとやや飛び過ぎているようです。今回使用した抗力係数CDの値0.33が、実際はもう少し大きいことが原因として考えられます。

[計算結果2おわり]

*****

今回でボールに作用する「重力」、「抗力」、「揚力(マグヌス力)」の3つの力を、全方向の速度に対して計算できる軌道シミュレータが完成しました。

加速度、速度からの時間積分を数値積分の方法で実行していますが、Microsoft Officeのエクセルや類似の表計算ソフトで簡単にできますので、よかったら作ってみてくださいね。

では、また。

（よかったら押してください。）

　 　
　にほんブログ村