バンクを走行する自転車の簡単な数理モデル
現在、私は簡単な競輪のレースシミュレーションプログラムを書いている。そこで使っている、バンクを単独で走行する自転車の運動の簡単な数理モデルを紹介する。
力学モデル
まず、単純に自転車に乗った選手を質量Mの質点とみなす。
選手がペダルをこいでその力がチェーンやギアを介して車輪を回し、車輪が摩擦によって走路に力を及ぼしてその反作用で自転車が動く。これを単純化して何らかの外力Faが質点に働いていると考える。Faは一般に時間の関数である。
走行する自転車には空気抵抗によって走るのと反対方向に力が加わる。この大きさは自転車の速度vの2乗に比例している。また、空気抵抗に比べたらかなり小さいが、走路の転がり摩擦とかもある。これは速度に関係なく一定である。両者の和をFrとすると、
Fr = Ca*v^2 + Cg ... (1)
と書ける。Ca、Cgは定数である。
カーブを曲がるには進行方向と垂直な方向の力が必要である。平地を自転車で走る場合は、その力はタイヤと地面の摩擦で得られるが、摩擦力には限りがあるので全速で急なカーブを曲がることができない。速すぎるとすべってコースアウトしてしまう。だから、速度を調節したり、なるべく曲率半径の大きいコース取りをしたりと考える必要がある。
しかし、バンクにはカントがあるので十分な向心力が簡単に得られる。よって選手はカーブでも速度を落としたりコースを変更することはない。実際には直線よりもカーブは若干走りにくいのだろうが本モデルではそれは無視して、単純にカーブも直線であるとみなしてしまう。直線上を運動する1次元の問題と考える。
よって、Newtonの運動方程式は速度vを変数として
Fa - Fr = M*dv/dt ... (2)
となる。
生理学モデル
選手がペダルをこぐ力には限界があり、またこぎ続けると疲れて出せる力が小さくなってくる。人体にはエネルギが化学的に蓄えられていて、それを消費して自転車を動かす仕事をしている。蓄えられたエネルギ、出せる最大の力、疲れたときの力の落ち方を説明するモデルを考える。私はよく知らないが、運動生理学なんかでこういうモデルが詳しく研究されていると思うのだが、ここでは非常に簡単に以下のように考える。
人体には使えるエネルギEがあるとする。ある瞬間に出せる仕事率すなわちパワーの最大値は、Ceを定数としてE/Ceであるとする。エネルギが減ってくると最大パワーも小さくなることを最も簡単に表現したモデルである。
パワー=力x速度だからdE/dt=Fa*vとなり、
Fa = (1/v)*(dE/dt) ... (3)
と選手が与える力Faが決まる。
また、選手が常に最大パワーを出すとすると、
dE/dt = -E/Ce
となり、これは良く知られた微分方程式で解は、
E = E0 * e^(-t/Ce)
である。E0は時刻0でのエネルギ、eは自然対数の底である。定数Ceは時間の次元を持ち、Ceだけ時間が経過するとエネルギがあらかた無くなる。
以上でモデルができた。選手が持つエネルギEが時間の関数として与えられると、運動方程式(2)と(1)、(3)式より速度vが決まる。速度が分かれば位置も速度を積分して求まる。
実測値での検証
自転車競技者向けのCS-netというwebサイトがあって、そこのスポーツバイオメカニクスの解説に速度などの実測データがある。これを元にモデルのパラメタを大雑把に決めてみる。
まず、抵抗力(1)式に関する定数は自転車運動の抵抗とパワーにある図15から、Ca=0.195(N*秒^2/m^2)、Cg=3.0(N)と読める。このグラフは選手と自転車の質量の合計が79(kg)の場合のものなので、以後も質量Mはこの値に大雑把に決めてしまう。
あと選手固有のパラメタとして、初期エネルギE0、定数Ceがある。これらには実走行中の速度、トルク、パワー(SRMパワーメータによる測定)にあるグラフを使う。そこの図8で静止状態から全力で250m走ったときの速度とパワーの変化が分かる。
パワーは60mくらいまで上昇し以後下降している。しかし、このモデルだとパワーはE/Ceなのでその最高値はEが最大のときすなわち時刻0での値E0/Ceとなり、時間がたつと単調に減少する。いきなり合わない。
しょうがないので、図8でパワーが上昇する60m付近までのデータを無視する。60m以降のパワーが単調に減少する部分のみを考える。しかしこれはそれほど無茶なことではない。ペダル踏力とクランク軸回りのパワー変化(30秒全力漕ぎの場合)によると、ペダルをこぐ運動は速度が遅いときには効率が悪いらしい。大きな力を出しても回転に役にたたない無駄な力が多いのだ。だから、消費されたエネルギの多くの部分が運動エネルギにならないのでこのモデルでは合わないと判断できる。
60m走った時点での残っているエネルギをE0だと考える。図8からこの付近のパワーは1560(W)くらいだと読めるので、E0/Ce=1560(W)なる関係があることになる。これを満たすE0、Ceの組を適当に決めて、運動方程式(2)の解を計算機を使った数値計算で求めた。すると、パワーの変化の実測値とよく合うE0、Ceの組が分かった。しかし、それで速度の変化を見てみると合わない。計算解は最高速度が実測より1.0(m/秒)ほど小さくなってしまう。パワーが合うなら速度も合いそうなうなものである。そこで空気抵抗をちょっと小さくしてみるとうまくいった。(1)式のCaを小さくしたのである。
Ca=0.16(N*秒^2/m^2)、E0=28080(J)、Ce=18.0(秒)としたときの、位置とパワーの関係を数値解と実測値でプロットする。
同様に、位置と速度は下になる。
まあまあよく合っている。速度がピークを過ぎた後のデータも欲しいところだし、Caを変えたのが妥当かという問題もあるが、元々大雑把な計算なのであまり気にしてもしょうがない。少なくとも本モデルは大はずれではないだろう。
TU = 0.05 MASS = 79.0 C_GROUND = 3.0 def sim(x, v, energy, c_e, c_air) t = 0.0 500.times do fr = c_air * v * v + C_GROUND pw = energy / c_e fa = 1/v * pw dv = TU * (fa - fr) / MASS v += dv x += TU * v energy -= TU * pw printf("%5.2f %5.2f %5.2f %8.2f %5.2f\n", t, x, v, energy, pw) t += TU end end ce = 18.0 e0 = 1560 * ce sim(60.0, 13.8, e0, ce, 0.16)
