野球解析で用いる物理モデル(空気抵抗編)

自己紹介+注意事項

初めまして、Gozaと申します。
神戸大学に通う学部新4年生で、専攻は物理学です。大学の講義で学んだ
"流体力学"と呼ばれるものに興味を持ちまして、自由研究として
野球の投球シミュレーターの開発やボールの研究を行っています。
note初心者であり、数式をうまく入れられないので一度overleafでコンパイルしたものをスクショして載せます。
野球経験はほぼなく、ど素人なので技術・感覚的なことはわかりません。
また、あくまでも独学で行っていますので下記で紹介する物理モデルや考察には粗さがあるということを前提とした上で、ご一読ください。
以下にこれまでのポートフォリオ(研究記録)を載せておきます。
noteで書いている物理モデルのイメージ図やそれらを実装した投球シミュレーターを導入しています。

野球ボールにかかる力について

野球ボールにかかる力について、現在4つの力で考えるのが主流になっています。

1：重力

おそらく最も有名な力の1つでしょう。ボールから地面に向かって働く力で、その値は不変です。義務教育でも習うような内容ですし、特にこれと言って語ることはありません。

2：空気抵抗力

最も有名な力その2。ボールの進行方向と逆向きに働く力であり、一般には速度の2乗に比例すると言われています。
正直2乗則を使うだけで精度は高いのですが、少し面白い法則があるのでそれを再現するために新たな物理モデルを作りました。

3：揚力(マグヌス力)

ボールの回転によって進行方向と垂直に働く力です。野球が好きな方は「マグヌス力」として馴染みがあるかもしれません。
藤川球児さんの「火の玉ストレート」に挙げられるようなストレートの浮き上がる力も揚力ですし、カーブやスライダーなど、回転している変化球であれば何でも揚力で基本説明がつきます。
詳しくは別noteで

4：縫い目の力(SSW)

野球ボールが完璧な球体でないことから及ぼされる力です。
通称Seam-Shifted-Wake(縫い目で後流が変化する現象)。
ナックルボールなど、回転数の少ない変化球で見られますが、導入が非常に難しいため無視しています。

空気抵抗について考えよう

空気抵抗の特性について

前節でも述べた通り、空気抵抗は基本的には速度の2乗則に従えば基本何とかなります。ただ調べてみるとこんなグラフが見つかりました。

レイノルズ数は速度に依存する関数です。$${C_d}$$というのは空気抵抗係数です。上のグラフはその二つの関係を表しているわけですが、速度の二乗則に従うのであれば、右肩上がりのグラフを描くはずです。
ただし、実際は
ある速度を境にガタンと抵抗が急減少します。
これを再現しようというのが今回のnoteの説明です。

物理式を整理しよう

先ほど紹介したレイノルズ数や空気抵抗力を立式します。
ここでは結果だけ書きますが、詳しい導出などはポートフォリオに書いてありますので、気になった方はご確認ください。

レイノルズ数：
レイノルズ数は慣性力と粘性力の比で表され、整理すると以下の通りになります。
速度をkm/h単位に変換して値を代入すると2つ目の式になります。

空気抵抗力：
空気抵抗力は空気の運動エネルギーとも言える動圧qを用いて表します。

Cdのモデルを再現しよう(ver1.0)

ここからは先ほどのNASAのグラフであった$${C_d}$$のモデルを再現する作業に入ります。いくつか特徴があるので分けて考えていきます。

1：急激な下降

図1のグラフでは
$${R_e \times 10^{-5}\simeq 1.0}$$
で急激に抵抗が下降していきます。
通称Drag Crisisと呼ばれるものであり、私の物理モデルでは
シグモイド関数と呼ばれるものを使います。

$${R_{e t r a n s}}$$はこのグラフの転換ポイントであり、変数$${R_{e }}$$がこの値を超えると、分母が大きくなりこの関数は無視できるほど小さくなります。
係数kがその下降の激しさを表し、今回はk=10として扱います。

2：緩やかな上昇

Drag Crisisの後ですが、緩やかに上昇していることがわかります。
これに関しては対数関数を使います。

転換ポイントである$${R_e \times 10^{-5}\simeq 1.0}$$を超えるまではこの項は0であり、Drag Crisis後から作用する関数です。

3：合体！！！

最小値である$${C_{d _{m i n}}}$$と最大値の$${C_{d _{\text {laninar }}}}$$を定義し、2つの式を合わせるとこんな感じになります。

青線の部分はDrag Crisis、赤線の部分はlogの上昇、紫がその両方が作用している形になります。
NASAのグラフを参考に

と見積もってプロットした結果が以下の通りになります。大分NASAのグラフは再現できたと思います。

Cdのモデルを再現しよう(ver2.0)

前節のものは私がnpb+を用いたデータ解析を行う前に作成した物理モデルです。ポートフォリオであげたのもここまでです。
ただデータ解析を行うに当たり、もう少し生のデータに近づけたり他の物理現象との関係を考えるべきだと思ったため改良します(2026-03)。

回転と空気抵抗の相関

データ野球を語る上で欠かせない偉大なるAlan Nathan教授の2020年の記事で、
$${C_d}$$には回転数と中程度の正の相関があるとわかったみたいです。

また、baseball-savantや、MLBのレポートから、抵抗係数が0.30を下回ることがあまりないことがわかりました。

https://baseball.physics.illinois.edu/HRReport2018.pdf

このことから、私のモデルでも回転数に依存するような空気抵抗モデルに改良するべきだと感じました。

以下ではMLBの開閉式ドーム球場で特に閉じている割合が高い、Globe Life Fieldの2025-04-01~2025-06-30までに取れたピッチングデータの中から、回転軸のわかりやすいフォーシームに絞って相関関係を調べます。

$${C_d=C_{d0}+k\omega}$$ モデル

まずは、記事にあった通り、回転数と抵抗係数の関係を調べます。
ピッチングデータでは回転数がしっかり書いてますので、それをそのまま単回帰モデルに当てはめます。

MLBモデルの回帰結果(Cd=Cd0+kω)
平均球速:144.77670579811038
Cd_0:0.3215
k:0.0000262281
相関係数r:0.2001

平均球速がやたら遅いですが、
これは初速と終速から平均を算出したものをさらに集めて平均をとったものだからです。

$${C_d=C_{d0}+kC_l}$$ モデル

揚力係数$${C_l}$$を使ったモデルです。
運動方程式から空気抵抗と揚力を逆算し、それを回帰分析します。

MLBモデルの回帰結果(Cd=Cd0+kCl)
平均球速:144.77670579811038
Cd_0:0.2995
k:0.3746804018
相関係数r:0.4783

$${C_d=C_{d0}+kC_l^2}$$ モデル

航空力学では、揚力係数の2乗に比例したモデルが使われているので、こっちでも試してみます。

MLBモデルの回帰結果(Cd=Cd0+k*Cl^2)
平均球速:144.77670579811038
Cd_0:0.3392
k:0.8703
相関係数r:0.4750

合体！！！！！

以上の相関係数の結果から揚力係数で見るのが良さそうです。
航空力学で使われている形でもあることから、$${C_d=C_{d0}+kC_l^2}$$ モデルを採用します。

ver1.0のモデルの特徴を組み合わせた上で、以下の通り立式しました。

先ほど$${C_d=C_{d0}+kC_l^2}$$ モデルで導いた平均球速を抵抗係数が最小値を取るポイントとして設定し、
調整項のオフセットを導入しました。

ここからは打球データを用いてオフセットの最適化を始めます。
取るデータは先ほどと同じ期間・球場で、回転がバックスピンになりやすいセンターフライに絞ります。

実際の飛距離と、バックスピン前提で考えた運動方程式で導いた飛距離の差の絶対値を集計し、最も差の小さいものをオフセットとします。

以上より、$${C_{d_{\text{offset}}}=0.01}$$とし、ver2.0完成とします。

Cdのモデルを再現しよう(ver3.0)

ver2.0のモデルでは、基準値を意識し$${ΔC_{dg}}$$という項を追加していましたが、
別のモデルなどで検証を行っていた際、$${ΔC_{rise}}$$や速度の塩梅によって$${ΔC_{dg}}$$が極端にマイナスを取ってしまい
飛ばないはずなのにめっちゃ飛ぶという問題が発生しました。
$${ΔC_{dg}}$$を削除し、新たなモデルを構築します。

2.0よりだいぶシンプルな形になってます。オフセットは以下の通りで
$${C_{d_{\text{offset}}}=-0.03}$$とします。

実験で使用した基本コード(Python)

from google.colab import drive
drive.mount('/content/drive')
from pybaseball import statcast
import os
import pandas as pd
mlb_3months_path="/content/drive/MyDrive/mlb_statcast_3months.csv"

print("MLBデータを収集します。")
if os.path.exists(mlb_3months_path):
  df_mlb=pd.read_csv(mlb_3months_path)
else:
  df_mlb=statcast(start_dt="2025-04-01",end_dt="2025-06-30")
  df_mlb.to_csv(mlb_3months_path,index=False)

print("収集完了しました。")

#物理モデルや計算式まとめるところ(ver3.0)
from matplotlib.lines import lineStyles
import pandas as pd
import numpy as np
from scipy.stats import linregress
from google.colab import drive
drive.mount('/content/drive')
import matplotlib.pyplot as plt
import japanize_matplotlib
from pybaseball import statcast
from scipy.integrate import solve_ivp #微分方程式ツール
from scipy import stats
csv_file = "/content/drive/MyDrive/npb_fly_inplay.csv" #集めたデータを入れてね
df_npb = pd.read_csv(csv_file)
#ドーム球場に絞る(屋外も欲しかったらコメントアウトしてね)
#df_npb=df_npb[df_npb["Stadium_Name"].str.contains("ドーム",na=False)]
#定数や使う物理量
rho = 1.225
radius = 0.0368
eta = 1.789e-5
constFactor = (rho * 2 * radius) / (eta * 3.6) #レイノルズ数変換の時にこれと速度かければ出せる
m = 0.145
g = 9.8
area = np.pi * radius**2
ft_to_m = 0.3048 #フィートからメートルへの変換
mp_to_km=1.60934 #マイルからキロメートルへの変換
cd_offset=-0.03 #最適化して出した(ヒント参照)

#statcastのピッチングデータをもとにCdの基本モデルを導出
print("【1】statcastのピッチングデータからCdの基本モデルを作ります。")
print("上のコードを実行していることを確認してください。")


#ほぼ屋根閉じてるレンジャースで尚且つバックスピンでSSWの小さい4seamを選択
df_dome=df_mlb[(df_mlb["home_team"]=="TEX")&(df_mlb["pitch_type"]=="FF")].copy()

cols=["release_speed","release_spin_rate","vx0","vy0","vz0","ax","ay","az","release_extension"]
df_physics=df_dome[cols].dropna().copy()
#dropna()で一つでも欠損(NA)がある行を削除

# 1. リリース位置と50ft地点の時間のズレを補正し、リリース時の速度(vxR, vyR, vzR)を算出
df_physics['yR'] = 60 - df_physics['release_extension']
df_physics['tR'] = (-df_physics['vy0'] - (df_physics['vy0']**2 - 2 * df_physics['ay'] * (50 - df_physics['yR']))**0.5) / df_physics['ay']

df_physics['vxR'] = df_physics['vx0'] + df_physics['ax'] * df_physics['tR']
df_physics['vyR'] = df_physics['vy0'] + df_physics['ay'] * df_physics['tR']
df_physics['vzR'] = df_physics['vz0'] + df_physics['az'] * df_physics['tR']

# 2. ホームベース（17/12 ft）までの飛行時間(tf)を算出
df_physics['tf'] = (-df_physics['vyR'] - (df_physics['vyR']**2 - 2 * df_physics['ay'] * (df_physics['yR'] - 17/12))**0.5) / df_physics['ay']

# 3. 軌道全体の「平均速度ベクトル(vxbar, vybar, vzbar)」と「平均スピード(vbar)」を算出
df_physics['vxbar'] = (2 * df_physics['vxR'] + df_physics['ax'] * df_physics['tf']) / 2
df_physics['vybar'] = (2 * df_physics['vyR'] + df_physics['ay'] * df_physics['tf']) / 2
df_physics['vzbar'] = (2 * df_physics['vzR'] + df_physics['az'] * df_physics['tf']) / 2
df_physics['vbar'] = (df_physics['vxbar']**2 + df_physics['vybar']**2 + df_physics['vzbar']**2)**0.5

# 4. 3Dベクトル計算により、マグナス力と重力(+32.174)を排除した「純粋な空気抵抗(adrag)」を抽出
df_physics['adrag'] = -(df_physics['ax']*df_physics['vxbar'] + df_physics['ay']*df_physics['vybar'] + (df_physics['az']+32.174)*df_physics['vzbar']) / df_physics['vbar']

#statcastはy軸がキャッチャーかピッチャーの方向を考慮して計算
df_physics["v_m_s"]=df_physics['vbar']*ft_to_m
df_physics["ax_m"]=df_physics["ax"]*ft_to_m
df_physics["ay_m"]=df_physics["ay"]*ft_to_m
df_physics["az_m"]=df_physics["az"]*ft_to_m

#ma_y=-1/2*rho*v^2*area*C_dよりC_dを導出
df_physics["force_common"]=2*m/(rho*area*df_physics["v_m_s"]**2)

#拾うデータがa_y<0より、絶対値化する
# 5. 完全なCdを算出！
df_physics['Cd'] = df_physics['adrag'] / (5.153E-03 * df_physics['vbar']**2)

#揚力に関しても同じようにやりたいが、上下(z軸)左右(x軸)考慮する必要があるので
#ma_x=1/2*rho*v^2*area*Cl_xとma_z=-mg+1/2*rho*v^2*area*Cl_zから三平方で導出
df_physics["Cl"]=np.sqrt(df_physics["ax_m"]**2+(df_physics["az_m"]+g)**2)*df_physics["force_common"]

#もしくはスピンパラメータから導出(相関が弱くなるので使わない(失敗例はヒント))
df_physics["omega_rad_s"]=df_physics["release_spin_rate"]*(2*np.pi/60)
df_physics["S"]=(radius*df_physics["omega_rad_s"])/df_physics["v_m_s"]
df_physics["Cl_theory"]=0.6564*(df_physics["S"]**0.7)

#あり得る値だけ考える(0.2<cd<0.6, 0<cl<0.7)
df_valid=df_physics[(df_physics["Cd"]>0.2)&(df_physics["Cd"]<0.6)
                    &(df_physics["Cl"]>0)&(df_physics["Cl"]<0.7)].copy()

print(f"{len(df_valid)}球のフォーシームを解析します。")

#平均球速
mean_speed_ms=df_valid["v_m_s"].mean()
V_REF=mean_speed_ms*3.6
RE_REF=V_REF*constFactor
#回帰分析(Cd=Cd0+k*Cl^2)を行う
df_valid["Cl_squared"]=df_valid["Cl"]**2

#最小二乗法
slope_1, intercept_1, r_value_1, p_value_1, std_error_1=linregress(df_valid["Cl_squared"], df_valid["Cd"])
K_MLB_1=slope_1#k値
CD_REF_MLB_1=intercept_1

print(f"MLBモデルの回帰結果(Cd=Cd0+k*Cl^2)")
print(f"平均球速:{V_REF}")
print(f"Cd_0:{CD_REF_MLB_1:.4f}")
print(f"k:{K_MLB_1:.4f}")
print(f"相関係数r:{r_value_1:.4f}")
#おそらく相関係数は0.5以下と若干弱い値を取るが、
#mlbのボールの作りの荒さやその他環境を考慮した上で強い相関と解釈する

#回帰分析(Cd=Cd0+kω)を行う
slope_2, intercept_2, r_value_2, p_value_2, std_error_2=linregress(df_valid["release_spin_rate"], df_valid["Cd"])
K_MLB_2=slope_2#k値
CD_REF_MLB_2=intercept_2
print(f"MLBモデルの回帰結果(Cd=Cd0+kω)")
print(f"平均球速:{V_REF}")
print(f"Cd_0:{CD_REF_MLB_2:.4f}")
print(f"k:{K_MLB_2:.10f}")
print(f"相関係数r:{r_value_2:.4f}")

#回帰分析(Cd=Cd0+kcl)を行う
slope_3, intercept_3, r_value_3, p_value_3, std_error_3=linregress(df_valid["Cl"], df_valid["Cd"])
K_MLB_3=slope_3#k値
CD_REF_MLB_3=intercept_3
print(f"MLBモデルの回帰結果(Cd=Cd0+kCl)")
print(f"平均球速:{V_REF}")
print(f"Cd_0:{CD_REF_MLB_3:.4f}")
print(f"k:{K_MLB_3:.10f}")
print(f"相関係数r:{r_value_3:.4f}")

print("【2】解析で用いる計算式を準備します。")
#C#で以前作ったプロット
def drag_crisis(Re, Cd_laminar=0.5, Cd_min=0.25, Re_trans=1.0e5, k=10, rise_factor=0.15):
    decay = (Cd_laminar - Cd_min) / (1 + np.exp(k * (Re - Re_trans) / Re_trans))
    rise = rise_factor * np.log10(np.maximum(1, Re / Re_trans))
    return Cd_min + decay + rise

#シミュレーションから得た改良版
def drag_crisis_remaster(re_val):
    Cd_laminar = 0.5
    Re_trans = 1.0e5
    k_trans = 10
    rise_factor = 0.15
    decay = (Cd_laminar-CD_REF_MLB_1) / (1 + np.exp(k_trans * (re_val - Re_trans) / Re_trans))
    rise = rise_factor * np.log10(np.maximum(1, re_val / 1.0e5))
    return decay + rise

#調整項(オフセット)を追加
def baseball_calculate_perfect(t, state, rpm, cd_offset,tilt_angle=0):
  #変数をstatcastに合わせる
  #y:キャッチャーからピッチャーへの飛距離
  #z:高さ
    y, z, vy, vz = state
    v_ms = np.sqrt(vy**2 + vz**2)
    v_kmh = v_ms * 3.6

    #Cl
    omega_rad_s = rpm * (2 * np.pi / 60)
    S = (radius * omega_rad_s) / v_ms
    cl_current = 0.6564 * (S ** 0.7)
    #回転軸を設定してシミュレーションできるように
    #0°でバックスピン、45°でフック
    cl_vertical=cl_current*np.cos(np.radians(tilt_angle))

    #Cd (Drag Crisis + オフセット + スピンペナルティ)
    re=v_kmh*constFactor

    cd_base = CD_REF_MLB_1 + cd_offset+drag_crisis_remaster(re)
    cd_total = cd_base + K_MLB_1 * (cl_current**2)

    #運動方程式
    Fd = 0.5 * cd_total * rho * area * v_ms**2
    Fl = 0.5 * cl_vertical * rho * area * v_ms**2
    ay = (-Fd * (vy/v_ms) - Fl * (vz/v_ms)) / m
    az = (-m*g - Fd * (vz/v_ms) + Fl * (vy/v_ms)) / m
    return [vy, vz, ay, az]

def hit_ground(t, state, rpm, cd_offset,tilt_angle=0): return state[1]
hit_ground.direction = -1
hit_ground.terminal = True

#引数に cd_offset を追加し、argsで微分方程式に渡す
def simulate_perfect(v0_kmh, angle_deg, rpm, cd_offset,tilt_angle=0):
    v0_ms = v0_kmh / 3.6
    rad = np.radians(angle_deg)
    y0 = [0.0, 1.0, v0_ms * np.cos(rad), v0_ms * np.sin(rad)]
    res = solve_ivp(baseball_calculate_perfect, [0, 15], y0, args=(rpm, cd_offset,tilt_angle), events=hit_ground, max_step=0.03)
    return res.y[0][-1]
  #solve_initial_value_problemは初期値を与えたらbaseball_calculateを実行して
  #返した値の積分値を行列でまとめてくれるやつ(今回ならy,z,vy,vzを持ってくる)
  #(使う関数,関数の初めの変数の範囲(t),tその次の変数(state),
  #追加で必要な変数(args(arguments)の中に入れる),条件(events),初めの変数の増加幅(max_step))
  #res.yはbaseball_calculateで帰ってくる値のグループ
  #(minimizeなど、res自体が答えを求める系のやつだったらres.xにして、シミュレート系はres.y)
  #res.y[0]はy座標、[-1]は一番最後のデータなのでz=0になった瞬間のy

print("【3】フライの条件を設定し、mlbとnpbのデータを引き出します。")
#MLB
df_hits_mlb = df_mlb[(df_mlb['home_team'] == 'TEX') &
                     (df_mlb['type'] == 'X') & #Sストライク,Bボール,Xインプレー
                     (df_mlb['hit_distance_sc'].notnull()) &
                     (df_mlb['events'] != 'home_run')].copy()

df_hits_mlb['Speed_kmh'] = df_hits_mlb['launch_speed']*mp_to_km
df_hits_mlb['Angle_deg'] = df_hits_mlb['launch_angle']
df_hits_mlb['Distance_m'] = df_hits_mlb['hit_distance_sc']*ft_to_m

df_fly_mlb = df_hits_mlb[(df_hits_mlb['Distance_m'] > 50) &
                           (df_hits_mlb['Angle_deg'] > 15) ].copy()
#1ヶ月分だけ取り出す(一旦日付を文字から数字へ)
df_fly_mlb["game_date"]=pd.to_datetime(df_fly_mlb["game_date"])

start_date = '2025-04-01'
end_date   = '2025-04-30'

df_valid_mlb=df_fly_mlb[
    (df_fly_mlb["game_date"]>=start_date)&
    (df_fly_mlb["game_date"]<=end_date)].copy()
print(f"NPBとの比較で{start_date}~{end_date}(計{len(df_valid_mlb)}件)取り出しました。")
#NPB
df_valid_npb = df_npb

print("【4】センターフライに絞り、推測の近似式の整合性を考えます。")
df_cf_fly=df_fly_mlb[(df_fly_mlb["hit_location"]==8.0)& #センターが捕球
                     (df_fly_mlb["events"]=="field_out")&
                     (df_fly_mlb["bb_type"]=="fly_ball")].copy() #フライアウト
print(f"センターフライ{len(df_cf_fly)}件を拾いました。")

print("設定完了しました、必要な関数は適宜このコードに更新してください。")