このページの内容は既に広く認められた理論や実験結果に基づくものではなく、趣味で個人的に考えているものです。専門家でもありません。正確性は保証できませんので参考程度に留めておいてください。

軌道シミュについて(編集中)

問題設定

テニスボールを質点として扱い,飛翔中には重力,揚力(マグヌス力?)、 抗力が作用するとする。また,地面との衝突では,鉛直方向には反発係数を, 水平方向には動摩擦と転がりへの移行を考慮した実効的な反発係数を用いる。 コート長手方向をx軸,幅方向をy軸,鉛直上向きをz軸とする。ボールの位置,速度,角速度をそれぞれ \(\boldsymbol{x}=\left(x、y、z\right)、\boldsymbol{v}=\mathrm{d}\boldsymbol{x}/\mathrm{d}t、\omega \) と表す。重力加速度は g=9.8 m・s^-2,ボールの質量と半径はm=57.7 g、 r=0.0335 mとする。流れに垂直な投影面積は\( A=\mathrm{\pi }r^{2} \)になる。 ボールの局所座標には球座標系を用いる。例えば初速度の大きさを\(V_{0}\)[m/s],仰角を\(θ\)[rad],水平面内の方位角を\(ϕ\)[rad]とすると, 初速度は \[{\boldsymbol{v}}_{0}\mathrm{[m/s]}=V_{0}\begin{pmatrix}\cos{\theta }\cos{\phi } \\ \cos{\theta }\sin{\phi } \\ \sin{\theta }\end{pmatrix}\] である。 同様に,回転軸の仰角と方位角をそれぞれ\({\theta }_{\omega },{\phi }_{\omega }\)[rad]とし,回転数をn[rpm]とすると初期角速度は \[\boldsymbol{\omega }\mathrm{[rad/s]}=\frac{2\pi n}{60}\begin{pmatrix}{\cos{\theta }}_{\omega }{\cos{\phi }}_{\omega } \\ {\cos{\theta }}_{\omega }{\sin{\phi }}_{\omega } \\ {\sin{\theta }}_{\omega }\end{pmatrix}\] \(ω\)は初期値から変わらないものとする。

飛翔中の運動方程式

抗力を\(\boldsymbol{F_{d}}\),揚力を\(\boldsymbol{F_{l}}\)、 重力加速度ベクトルを\(\boldsymbol{g}\)とすると,ボールの運動方程式は \[\frac{{\mathrm{d}}^{2}\boldsymbol{x}}{\mathrm{d}{\mathrm{t}}^{2}}=\boldsymbol{g}+\frac{1}{m}\left(\boldsymbol{F_{d}}+\boldsymbol{F_{l}}\right)\] である。空気密度を\(\rho\)、ボールの投影面積を\(A\)、抗力係数と揚力係数をそれぞれ\(C_{d}\)と\(C_{l}\)とすると,抗力と揚力はよく知られた次の式で表される。(確かこの形になることは次元解析から言えた気がする?) \[\boldsymbol{F_{d}}=\frac{1}{2}\rho AC_{d}\left|\boldsymbol{v}\right|^2\hat{\boldsymbol{D}}\] \[\boldsymbol{F_{l}}=\frac{1}{2}\rho AC_{l}\left|\boldsymbol{v}\right|^2\hat{\boldsymbol{L}}\] ここで\(\hat{\boldsymbol{D}},\hat{\boldsymbol{L}}\)はそれぞれ抗力と揚力の働く方向を示す単位ベクトルである。 以降はこれに含まれる各要素をプレイ中に感じられるパラメータを使って表現していくことにする。

流体力の向き

まず抗力の向きは非常に簡単で\(\boldsymbol{v}\)が流体でなくボール側の速度であることに気を付けて正規化するだけでいい; \[\hat{\boldsymbol{D}}=\frac{-\boldsymbol{v}}{\left|\boldsymbol{v}\right|}\] 一方で揚力の方向は、2次元平面内で回転してるのであれば鉛直上下の2択だが、3次元の場合例えばスライスを考えればわかるように方位角と天頂角の回転軸の自由度分複雑になる。(流れに対し垂直な力ではないのでもはや揚力とは言わないかもしれないが これを表す言葉を知らないので便宜上以降も揚力と呼ぶ。マグヌス力?)しかし、結局力の向きはボール表面の速度差で決まる(流線に沿ったエネルギー的に)。 その速度差は粘性流体に対するボール表面のnonslip条件によって空気と表面がその境界で同じ速度になることに由来する。これは回転軸に関して軸対称な流れを生む。この段階では速度差はない。 一方で並進は並進軸に関して軸対称な流れを加える。これらの流れを組み合わせるとボールの中心を通り回転軸と並進軸の両方を含む面に対して非対称な圧力場が生じ、高圧側から低圧側へ力が働く。 これが知りたい向きである。条件は角速度ベクトルと速度ベクトルの両方に垂直な向きといえて、それはつまりこれらの外積の向きである; \[\hat{\boldsymbol{L}}=\frac{\boldsymbol{\omega}\times \boldsymbol{v}}{\left|\boldsymbol{\omega}\times \boldsymbol{v}\right|}\]

抗力係数、揚力係数

速度の単位ベクトルを \[\hat{\boldsymbol{v}}=\frac{\boldsymbol{v}}{\left|\boldsymbol{v}\right|}\] で表すと角速度のうち速度に垂直な成分は \[{\boldsymbol{\omega }}_{\perp }=\boldsymbol{\omega }-\left(\boldsymbol{\omega }\cdot \hat{\boldsymbol{v}}\right)\hat{\boldsymbol{v}}\] である(大学1年生御用達の直交化法)。スピンパラメータはこれを用いて \[s_{p}=\frac{r\left|{\boldsymbol{\omega }}_{\perp }\right|}{\left|\boldsymbol{v}\right|}\] と定義される(らしい。モデルかもしれない)。 抗力係数 Cdと揚力係数Clはスピンパラメータの関数として \[C_{d}\left(s_{p}\right)=0.55+\frac{1}{{\left[22.5+4.2{\left(1/s_{p}\right)}^{2.5}\right]}^{0.4}}\] \[C_{l}\left(s_{p}\right)=\frac{1}{2+1/s_{p}}\] とするモデルを採用する。sp=0の場合揚力係数は0とする。

湿潤空気の密度

気温をt[℃],絶対温度を T=273.15+t、大気圧P[Pa],相対湿度H_rとする。飽和水蒸気圧はTetensの経験式 \[P_{s}=610.8\times {10}^{\left(\frac{7.5t}{t+237.3}\right)}[\mathrm{Pa}] \] により計算する。水蒸気分圧は\(H_{r}P_{s}[Pa]\)となる。このとき湿潤空気の密度は状態方程式から \[{\rho }_{\mathrm{air}}=\frac{M_{d}P-\left(M_{d}-M_{w}\right)H_{r}P_{s}}{RT}[\mathrm{g}\cdot {\mathrm{m}}^{-3}]\] で与えられる。ここで\(M_{d}=28.966\,\mathrm{g}\cdot \mathrm{mol}^{-1}、 M_{w}=18\,\mathrm{g}\cdot \mathrm{mol}^{-1}、 R=8.314\, \mathrm{J}\cdot \mathrm{mol}^{-1}\cdot \mathrm{K}^{-1}\)であり、それぞれ乾燥空気のモル質量、水のモル質量、モルベース気体定数。

最終形

以上をまとめると運動方程式は \[\frac{{\mathrm{d}}^{2}\boldsymbol{x}}{\mathrm{d}{\mathrm{t}}^{2}}=\boldsymbol{g}-\frac{{\rho }_{\mathrm{air}}A}{2m}\left(C_{d}\left|\boldsymbol{v}\right|\boldsymbol{v}-C_{l}{\left|\boldsymbol{v}\right|}^{2}\hat{\boldsymbol{L}}\right)\] シミュレーションは4次のルンゲクッタ法を用いて実行している。誤差評価はしていない。

地面との衝突

ボール中心がz=rに達したときを地面との接触とする。衝突直前と直後の速度をそれぞれ\( {\boldsymbol{v}}^{-},{\boldsymbol{v}}^{+}\) と書き,鉛直反発係数をe=:CORとする。鉛直速度は \[v_{z}^{+}=-ev_{z}^{-}\] によって反転・減衰する。 水平方向の衝突には動摩擦係数 μ=CODF を用いる。ボールの慣性モーメントは通常の球体の式に実験的なパラメータαを乗じたものを用いる。: \[I=\alpha mr^{2}、\alpha =0.55\] 入射角に対応する量を \[\tan{\beta }=\frac{\left|{\boldsymbol{v}}_{h}^{-}\right|}{\left|v_{z}^{-}\right|}\] \[{\boldsymbol{v}}_{h}=\left(v_{x}、v_{y}\right)\] と定める。滑りが継続する条件は \[\mu \lt \frac{1-s_{p}}{\left(1+1/\alpha \right)\left(1+e\right)\tan{\beta }}+\frac{D}{\left(1+\alpha \right)r}\] である。垂直抗力の作用線と重心の距離Dを用いた。これを満たす場合,水平速度に対する実効反発係数を \[e_{h}=1-\mu \left(1+e\right)\tan{\beta }\] とする。一方,条件を満たさず,転がり状態へ移行すると判定された場合は \[e_{h}=0.645\] とする。これを用いて衝突後の水平速度は次式で求められる。: \[\boldsymbol{v}_{h}^{+}=e_{h}\boldsymbol{v}_{h}^{-}\]