二輪車の運動特性(1) ベンチマーク自転車の固有値解析

これまでの記事では、二輪車の運動特性を理解するための準備として、倒立振子を題材に固有値・固有ベクトルと挙動の関係を確認してきました。本記事から、いよいよ対象を二輪車そのものに移し、安定性解析モデルに対して固有値・固有ベクトルを求めていきます。

倒立振子では自由度が1つしかないため、固有ベクトルは単に状態変数（ロール角と角速度）の比を表すだけでしたが、ベンチマーク自転車モデルでは舵角の自由度が加わり、各状態変数がどのような比率で組み合わさるかという「運動のパターン（モード形状）」として固有ベクトルを解釈できるようになります。

すなわち、固有ベクトルは単なる状態の比ではなく、各モードにおける運動の現れ方そのものを表す量になります。

ベンチマーク自転車モデルの運動方程式については、以前の記事で紹介しました。今回はこれを状態方程式に変換し、実際に固有値・固有ベクトルを得る手順について解説します。

二輪車安定性解析モデルのシステム行列

ベンチマーク自転車モデルの運動方程式は、こちらの記事で紹介したように、次の形で表されます。

$$\bold{M}\ddot{\bm{q}} + v\bold{C_1}\dot{\bm{q}} + (g\bold{K_0} + v^2\bold{K_2})\bm{q} = \bm{f} \\$$

この式は二階の連立常微分方程式ですが、固有値解析を行うために倒立振子と同様の手法で一階の連立微分方程式に変換します。

まず、質量行列 $\bold{M}$ は通常正則（逆行列を持つ）であるため、運動方程式は次のように書き換えることができます。

$$\ddot{\bm{q}} + v\bold{M^{-1}C_1}\dot{\bm{q}} + \bold{M^{-1}}(g\bold{K_0} + v^2\bold{K_2})\bm{q} = \bold{M^{-1}}\bm{f}$$

さらに、$\dot{\bm{q}}=\bm{p}$ とおき、これを微分すると $\ddot{\bm{q}}=\dot{\bm{p}}$ となります。これらを元の式に代入して整理すると、次の連立微分方程式が得られます。

$$\begin{bmatrix} \dot{\bm{q}} \\ \dot{\bm{p}} \end{bmatrix} = \begin{bmatrix} \bold{O} & \bold{E} \\ -(g\bold{M^{-1}K_0}+v^2\bold{M^{-1}K_2}) & -v\bold{M^{-1}C_1} \end{bmatrix} \begin{bmatrix} \bm{q} \\ \bm{p} \end{bmatrix} + \begin{bmatrix} \bold{O} \\ \bold{M^{-1}} \end{bmatrix} \bm{f}$$

ここで、状態ベクトルを $\bm{x} = [\bm{q}\:\bm{p}]^T = [\bm{q}\:\dot{\bm{q}}]^T$ とおくと、この系は

$$\dot{\bm{x}} = \bold{A}(v)\bm{x} + \bold{B}\bm{f}$$

と表されます。システム行列 $\bold{A}(v)$ は

$$\bold{A}(v) = \begin{bmatrix} \bold{O} & \bold{E} \\ -(g\bold{M^{-1}K_0}+v^2\bold{M^{-1}K_2}) & -v\bold{M^{-1}C_1} \end{bmatrix}$$

で与えられ、速度 $v$ に依存することがわかります。

この $\bold{A}(v)$ の固有値・固有ベクトルを求めることで、各速度における安定性と運動パターン（モード）を調べることができます。

ただし、この行列はパラメータ $v$ を含む高次の行列であるため、固有値を解析的に求めることは一般に困難です。そのため実際には、具体的なパラメータを代入し、速度を変化させながら数値的に固有値問題を解くことになります。

ベンチマーク自転車の固有値

ここでは、Meijaardらの論文^[1]に記載されているベンチマーク自転車の諸元を使用し、Scilabによって車速 $v=0-10\,\mathrm{m/s}$ における固有値を計算しました。各車速における固有値の実部および虚部をプロットすると次のようになり、論文と同様の結果が得られていることが確認できます。また、固有値計算とプロットを実行するScilabスクリプトを末尾に示します。

固有値の実部を実線、虚部を破線で表し、複素共役な固有値のペアを同じ色で表しています。車速$v_d$以下の領域では複素数の固有値はありませんが、固有値の連続的な変化を追跡するため、同じ色で表しました。

グレーに塗った車速 $v_w < v < v_c$ の領域では全ての固有値の実部が負であるため、システムは漸近安定となります。つまりこの速度領域では、手放し運転などライダーが操作を行わなくても二輪車は直立して直進走行を保つことができる自己安定な領域となっており、二輪車が倒れずに走ることができることを表しています。

一方、$v < v_w$ の領域では赤線で描いたweaveの固有値実部が正となり、また$v_c < v$の領域では緑線のcapsize固有値が正となって、不安定となります。これらのモードはどのような挙動を示すのかを理解するため、対応する固有ベクトルで確認していきます。

ベンチマーク自転車の固有モード

まず、車速$v=5\,\mathrm{m/s}$での固有値・固有ベクトルの計算結果を例として、二自由度系における固有ベクトルの構造を確認します。末尾のScilabスクリプトを実行後、Scilabコンソール上で

sl=BicycleSys(5,9.81,M,K0,K2,C1);
[V D]=spec(sl.A);

を実行すると、車速 $v=5\,\mathrm{m/s}$ における固有値・固有ベクトルを得ることができます。

車速 $v=5\,\mathrm{m/s}$ での固有値は、$\lambda = [-14.078,\:-0.32287,\:-0.77534+4.4649i,\:-0.77534-4.4649i]$ で、それぞれに対応する固有ベクトルは以下のように計算されました。

$$\begin{bmatrix} \phi \\ \delta \\ \dot{\phi} \\ \dot{\delta} \end{bmatrix} = \begin{bmatrix} -0.00016187 \\ -0.070852 \\ 0.0022788 \\ 0.99748 \end{bmatrix},\; \begin{bmatrix} -0.87481 \\ -0.37457 \\ 0.28245 \\ 0.12094 \end{bmatrix},\; \begin{bmatrix} 0.00034845 - 0.13152i \\ -0.029204 - 0.16818i \\ 0.58694 + 0.10353i \\ 0.77353 \end{bmatrix},\; \begin{bmatrix} 0.00034845 + 0.13152i \\ -0.029204 + 0.16818i \\ 0.58694 - 0.10353i \\ 0.77353 \end{bmatrix}$$

Scilabで計算された固有ベクトルは、ベクトルの大きさが1になるように正規化されていますが、固有ベクトルは方向のみを表し、大きさは不定なので、見やすくするためにここでは $\delta=1$ として $\phi$ との比を表すことにします。

また、固有モードに対応する解は $e^{\lambda t}$ に比例するため、速度成分 $\dot{\phi},\:\dot{\delta}$ は変位成分 $\phi,\:\delta$ の $\lambda$ 倍として表されます。この関係を用いると、速度成分については変位と固有値の積で表すことができ、以下のようになります。

$$\begin{bmatrix} \phi \\ \delta \\ \dot{\phi} \\ \dot{\delta} \end{bmatrix} = \begin{bmatrix} 0.0022846 \\ 1 \\ 0.0022846*(-14.078) \\ 1*(-14.078) \end{bmatrix},\; \begin{bmatrix} 2.3355 \\ 1 \\ 2.3355*(-0.32287) \\ 1*(-0.32287) \end{bmatrix},\; \begin{bmatrix} 0.75879 + 0.13384i \\ 1 \\ (0.75879 + 0.13384i)*(-0.77534 + 4.4649i) \\ 1*(-0.77534 + 4.4649i) \end{bmatrix},\; \begin{bmatrix} 0.75879 - 0.13384i \\ 1 \\ (0.75879 - 0.13384i)*(-0.77534 - 4.4649i) \\ 1*(-0.77534 - 4.4649i) \end{bmatrix}$$

したがって、固有ベクトルの変位成分を見れば、そのモードにおける各自由度の比率、すなわち運動パターンを読み取ることができます。ここからは、各モードごとにその特徴を見ていきます。

キャスタリングモード

一つ目のモードの固有値は $\lambda=-14.078$ で非常に安定であり、固有値のグラフ上では青の実線で表されいます。このモードは常に安定で、速度が増加するにつれて実部がより負の方向へ移動し、安定性が増加します。

また、固有ベクトルは $\delta$ の成分が支配的であり、操舵系がほぼ単独で応答するモードであることがわかります。これはステアリングが直進方向に収束するモードであることを表しており、ショッピングカートなどの車輪（caster）が自動的に進行方向を向くことになぞらえて、キャスタリングモードと呼ばれます。

速度ごとに固有ベクトルの振幅比を算出すると以下のようなグラフとなり、ほぼ速度によらず常に舵角 $\delta$ が卓越していることがわかります。

キャプサイズモード

二つ目のモードの固有値は $\lambda=-0.32287$ であり、固有値のグラフ上では緑の実線に対応します。このモードは低速域では安定ですが、速度が装荷して $v_c$ を超えると固有値の実部が正に転じ、不安定になることがわかります。

固有ベクトルはロール角 $\phi$ の成分が支配的であり、速度ごとの振幅比をキャスタリングモードのときとは逆に $\phi$ を基準として正規化すると以下のようなグラフとなります。速度が増加して不安定になる領域でロール角 $\phi$ が卓越し、その寄与がさらに大きくなることがわかります。このモードによる挙動は、低速では車体が直立状態に維持される一方、高速ではわずかに不安定になってゆっくり倒れていくことを表しています。この不安定なときの挙動は船が転覆する現象（capsize）と似ていることから、キャプサイズモードと呼ばれています。

ウィーブモード

三つ目のモードの固有値は $\lambda=-0.77534 \pm 4.4649i$ であり、共役複素数となっていて振動モードであることがわかります。固有値のグラフ上では実部が赤の実線、虚部が赤の破線で表されています。このモードはキャプサイズとは逆に高速では安定ですが、$v < v_w$の低速域では不安定になることがわかります。

複素固有値の場合、固有ベクトルも複素数となり、複素平面上での極形式で表すことで振幅（絶対値）と位相（偏角）の情報が得られます。そこで舵角 $\delta$ を基準としたロール角 $\phi$ の振幅比と位相差を速度ごとにプロットすると、次のようなグラフが得られます。

振幅比はおおむね0.6から0.8であり、ロール角の振幅は舵角よりやや小さいものの、舵角もロール角も左右に振動するモードであることがわかります。タイヤの接地拘束により、舵角に応じてヨーレート（旋回運動）が発生するため、車体全体が左右に傾きながら蛇行するような挙動を示します。このようにうねうねと蛇行する挙動（weave）からウィーブモードと呼ばれています。

位相をみると、ウィーブモードが安定な$v_w < v$の領域では舵角とロール角がおおむね同位相であるのに対し、$v < v_w$になるとロール角の位相が進む（舵角の位相が遅れる）ようになることが確認できます。このことから、操舵（旋回）によって生じる復元力（遠心力）の発生がロール運動に対して遅れることで、復元作用が十分に働かなくなり、不安定になるといえます。

車速が$v < v_d$になると、固有値が2つの実数根に分かれ、操舵とロールの連成が失われます。
舵角が固定されている場合には車体全体が倒立振子として振る舞いますが、操舵系が自由に動く場合には操舵系も含めた二重倒立振子と同等の挙動を示します。

まとめ

二輪車の挙動には、キャスタリング、キャプサイズ、ウィーブの3つの特徴的なモードがあり、特定の速度域（$v_w < v < v_c$）では自己安定性を示して倒れずに走行できることが、固有値解析によって説明できることがわかりました。また、高速ではキャプサイズモードが、低速ではウィーブモードが不安定になることが、固有値解析によって理論的に示されました。

また、ウィーブモードの安定・不安定の境界となる速度$v_w$はウィーブ速度、キャプサイズモードの安定・不安定の境界となる速度$v_c$はキャプサイズ速度と呼ばれます。これらの速度は設計諸元から直接算出することが可能であり、安定性設計における重要な指標の一つです。

次の記事では、これらの速度について解説します。

benchmark_bicycle.sce

//================================================
// Eigenvalue Analysis for the benchmark bicycle
//================================================
clear;
// Define Bicycle Specifications
function bicycle=DefBicycleSpec()
    bicycle.w = 1.02; // Wheel base
    bicycle.c = 0.08; // Trail
    bicycle.lambda = %pi/10; // Steer axis tilt
    bicycle.g = 9.81; // Gravity
    //// Rear WHeel R
    bicycle.rR = 0.3; // Radius
    bicycle.mR = 2; // Mass
    bicycle.IRxx = 0.0603; // Mass moment of inertia xx
    bicycle.IRyy = 0.12; // Mass moment of inertia yy
    //// Rear Body and frame assembly B
    bicycle.xB = 0.3; // Position centre of mass x
    bicycle.zB = -0.9; //Position centre of mass z
    bicycle.mB = 85; // Mass
    bicycle.IBxx = 9.2; // Mass moment of inertia xx
    bicycle.IByy = 11; // Mass moemnt of inertia yy
    bicycle.IBzz = 2.8; // Mass moement of inertia zz
    bicycle.IBxz = 2.4; // Mass product of inertia xz
    //// Front Handlebar and fork assembly H
    bicycle.xH = 0.9; // Position centre of mass x
    bicycle.zH = -0.7; // Position centre of mass z
    bicycle.mH = 4; // Mass
    bicycle.IHxx = 0.05892; // Mass moment of inertia xx
    bicycle.IHyy = 0.06; // Mass moment of inertia yy
    bicycle.IHzz = 0.00708; // Mass moment of inertia zz
    bicycle.IHxz = -0.00756; // Mass product of inertia xz
    //// Front wheel F
    bicycle.rF = 0.35; //Radius
    bicycle.mF = 3; // Mass
    bicycle.IFxx = 0.1405; // Mass moment of inertia xx
    bicycle.IFyy = 0.28; // Mass moemnt of inertia yy
endfunction
//
// Create Bicycle Matrices
function[M, K0, K2, C1]=BicycleMatrices(b)
    // Total system T properties
    mT = b.mR + b.mB + b.mH + b.mF; // Mass
    xT = (b.xB*b.mB + b.xH*b.mH + b.w*b.mF)/mT; // Position centre of mass x
    zT = (-b.rR*b.mR + b.zB*b.mB + b.zH*b.mH - b.rF*b.mF)/mT; // Position centre of mass z
    ITxx = b.IRxx + b.IBxx + b.IHxx + b.IFxx + b.mR*b.rR^2 + b.mB*b.zB^2 + b.mH*b.zH^2 + b.mF*b.rF^2;
    ITxz = b.IBxz + b.IHxz - b.mB*b.xB*b.zB - b.mH*b.xH*b.zH + b.mF*b.w*b.rF;
    ITzz = b.IRxx + b.IBzz + b.IHzz + b.IFxx + b.mB*b.xB^2 + b.mH*b.xH^2 + b.mF*b.w^2;
    // front assembly A properties
    mA = b.mH + b.mF; // Mass
    xA = (b.xH*b.mH + b.w*b.mF)/mA; // Position centre of mass x
    zA = (b.zH*b.mH - b.rF*b.mF)/mA; // Position centre of mass z
    IAxx = b.IHxx + b.IFxx + b.mH*(b.zH - zA)^2 + b.mF*(b.rF + zA)^2;
    IAxz = b.IHxz - b.mH*(b.xH - xA)*(b.zH - zA) + b.mF*(b.w - xA)*(b.rF + zA);
    IAzz = b.IHzz + b.IFxx + b.mH*(b.xH - xA)^2 + b.mF*(b.w - xA)^2;
    uA = (xA - b.w - b.c)*cos(b.lambda) - zA*sin(b.lambda);
    IAll = mA*uA^2 + IAxx*sin(b.lambda)^2 + 2*IAxz*sin(b.lambda)*cos(b.lambda) + IAzz*cos(b.lambda)^2;
     IAlx = -mA*uA*zA + IAxx*sin(b.lambda) + IAxz*cos(b.lambda);
     IAlz = mA*uA*xA + IAxz*sin(b.lambda) + IAzz*cos(b.lambda);
     // Other parameters
     mu = (b.c/b.w)*cos(b.lambda);
     SR = b.IRyy/b.rR;
     SF = b.IFyy/b.rF;
     ST = SR + SF;
     SA = mA*uA + mu*mT*xT;
     // Coefficient Matrices
     M = [ITxx IAlx+mu*ITxz; IAlx+mu*ITxz IAll+2*mu*IAlz+mu^2*ITzz];
     K0 = [mT*zT -SA; -SA -SA*sin(b.lambda)];
     K2 = [0 (ST-mT*zT)*cos(b.lambda)/b.w; 0 (SA+SF*sin(b.lambda))*cos(b.lambda)/b.w];
     C1 = [0 mu*ST+SF*cos(b.lambda)+ITxz*cos(b.lambda)/b.w-mu*mT*zT; -(mu*ST+SF*cos(b.lambda)) IAlz*cos(b.lambda)/b.w+mu*(SA+ITzz*cos(b.lambda)/b.w)];
endfunction
//
// Create Bicycle System
function sl=BicycleSys(v,g,M,K0,K2,C1)
    MI = inv(M);
    gMIK0 = g*MI*K0;
    MIK2 = MI*K2;
    MIC1 = MI*C1;
    A = [zeros(2,2) eye(2,2); -(gMIK0 + v^2*MIK2) -v*MIC1];
    B = [zeros(2,2); MI];
    C = [eye(2,2) zeros(2,2)];
    D = zeros(2,2);
    sl = syslin('c', A,B,C,D);
endfunction
//
// Main routine
// Eigenvalue Calculation
bicycle = DefBicycleSpec();
[M K0 K2 C1]=BicycleMatrices(bicycle);
v=[0:0.1:0.68428307889246,0.7:0.1:10];
i=0;
for vi=v
    i=i+1;
    sl=BicycleSys(vi,bicycle.g,M,K0,K2,C1);
    [V D]=spec(sl.A);
    D=diag(D).';
    sort_key=[abs(imag(D)); real(D)];
    [dmy idx]=gsort(sort_key,"lc","i");
    E(i,:)=D(idx);
end
//
// Draw Eigenvalue Plots
f=scf();
f.axes_size=[600 600];
// add origin axis
xpoly([v(1),v($)],[0,0],'lines');
e = gce();
e.thickness = 2;
e.line_style = 1;
e.foreground = color("black");
// plot eigenvalues
plot(v,[E(:,1),E(:,2),real(E(:,3)),imag(E(:,3))]);
xlabel("v [m/s]"),xgrid(),ylabel("Eigenvalue");
// adjust the imaginary line style and the color to the real one
h1 = gce().children;
h1(1).foreground=h1(2).foreground;
h1(1).line_style=3;
// add the plot for the conjugate complex data
plot(v,[real(E(:,4)),imag(E(:,4))]);
h2 = gce().children;
h2(1).line_style=3;
h2.foreground=h1(2).foreground;
// ajust the y scale
replot([0 -10 10 10]);
// add legend
l=legend(h1(4:-1:1),['castering','capsize','weave(Re)','weave(Im)']);

脚注

1. J. P. Meijaard, et al., "Linearized dynamics equations for the balance and steer of a bicycle: a benchmark and review", Proc. R. Soc. A 463:1955-1982, 2007