モード分解とモード形式とは

前準備

次式のような2自由度の振動系を考える。

\begin{align} m_1\cfrac{d^2x_1}{dt^2}+c_1\cfrac{dx_1}{dt}+k_1x_1(t)+c_{12}\cfrac{d(x_1-x_2)}{dt}+k_{12}(x_1-x_2)=f_1(t)\\m_2\cfrac{d^2x_2}{dt^2}+c_2\cfrac{dx_2}{dt}+k_2x_2(t)+c_{12}\cfrac{d(x_2-x_1)}{dt}+k_{12}(x_2-x_1)=f_2(t) \end{align}

 上式を行列形式で書くと次式のようになる。

\begin{align} \left[\begin{array}{cc}m_1&0\\0&m_2\end{array}\right]\cfrac{d^2}{dt^2}\left(\begin{array}{c}x_1\\x_2\end{array}\right)&+\left[\begin{array}{cc}c_1+c_{12}&-c_{12}\\-c_{12}&c_{12}+c_2\end{array}\right]\cfrac{d}{dt}\left(\begin{array}{c}x_1\\x_2\end{array}\right)\\&+\left[\begin{array}{cc}k_1+k_{12}&-k_{12}\\-k_{12}&k_{12}+k_2\end{array}\right]\left(\begin{array}{c}x_1\\x_2\end{array}\right)=\left(\begin{array}{c}f_1\\f_2\end{array}\right) \end{align}

上式をまとめると次式のようになる。

\begin{align} \boldsymbol{M}\cfrac{d^2\boldsymbol{x}}{dt^2}+\boldsymbol{C}\cfrac{d\boldsymbol{x}}{dt}+\boldsymbol{K}\boldsymbol{x}=\boldsymbol{f} \end{align}

\(\boldsymbol{M}\)を質量行列、\(\boldsymbol{C}\)を減衰行列、\(\boldsymbol{K}\)を剛性行列、\(\boldsymbol{x}\)を位置ベクトル、\(\boldsymbol{f}\)を外力ベクトルを表す。この時、作用反作用の関係から質量行列、減衰行列、剛性行列は対称行列となることに留意する。

上式は3自由度の振動系だったり、もっと複雑な構造物の運動方程式に対しても、同様の形式で表される。

モード分解(Mode decomposition)

モード分解を考える際に、初めに無減衰系を考える。

無減衰系\((\boldsymbol{C}=\boldsymbol{0})\)

無減衰系の多自由度振動系の運動方程式は次のようになる。

\begin{align} \boldsymbol{M}\cfrac{d^2\boldsymbol{x}}{dt^2}+\boldsymbol{K}\boldsymbol{x}=\boldsymbol{f} \end{align}

質量行列\(\boldsymbol{M}\)と剛性行列\(\boldsymbol{K}\)の同時対角化を考える。通常の（同時）対角化問題では行列の相似変換\(\boldsymbol{B}=\boldsymbol{P}^{-1}\boldsymbol{A}\boldsymbol{P}\)を用いるが、ここでは合同変換\(\boldsymbol{B}=\boldsymbol{P}^\textrm{T}\boldsymbol{A}\boldsymbol{P}\)を考える。

下記の一般化固有値問題の解である固有ベクトルを並べた行列\(\boldsymbol{T}\)を考える。

\begin{align} \boldsymbol{M}\boldsymbol{T}\boldsymbol{\Omega}^2=\boldsymbol{K}\boldsymbol{T}\end{align}

\(\boldsymbol{\Omega}^2\)は固有値を対角成分に持つ対角行列（固有値行列）である。上式から左から\(\boldsymbol{T}^\textrm{T}\)をかけることで、対角行列\(\boldsymbol{\bar{M}}(=\boldsymbol{T}^\textrm{T}\boldsymbol{M}\boldsymbol{T})\)、\(\boldsymbol{\bar{K}}(=\boldsymbol{T}^\textrm{T}\boldsymbol{K}\boldsymbol{T})\)を得ることができる。その証明を簡単に示す。

• \(\boldsymbol{M}\)はコレスキー分解可能であるため、三角行列\(\boldsymbol{U}\)を用いて\(\boldsymbol{M}=\boldsymbol{U}^\textrm{T}\boldsymbol{U}\)のように分解する。\begin{align} \boldsymbol{U}^\textrm{T}\boldsymbol{U}\boldsymbol{T}\boldsymbol{\Omega}^2=\boldsymbol{K}\boldsymbol{T}\end{align}
• 固有ベクトル行列\(\boldsymbol{T}\)を三角行列\(\boldsymbol{U}\)を用いて変換する。\begin{align} \boldsymbol{W}=\boldsymbol{U}\boldsymbol{T}\end{align}
• 新しい固有ベクトル\(\boldsymbol{W}\)を用いると固有方程式は次のようになる。\begin{align} \boldsymbol{U}^\textrm{T}\boldsymbol{W}\boldsymbol{\Omega}^2=\boldsymbol{K}\boldsymbol{U}^{-1}\boldsymbol{W}\end{align}
• 上式の左から\(\boldsymbol{U}^\textrm{-T}\)をかけることで左辺は固有ベクトルと固有値の積で表されるようになり、ただの固有値問題となる。\begin{align} \boldsymbol{W}\boldsymbol{\Omega}^2=\boldsymbol{U}^\textrm{-T}\boldsymbol{K}\boldsymbol{U}^{-1}\boldsymbol{W}\end{align}\(\boldsymbol{K}\)が対称行列であるため、\(\boldsymbol{U}^\textrm{-T}\boldsymbol{K}\boldsymbol{U}^{-1}\)は対称行列となるので、この固有値問題の固有ベクトル行列である\(\boldsymbol{W}\)は直行行列である。
• \(\boldsymbol{T}^\textrm{T}\boldsymbol{M}\boldsymbol{T}\)と\(\boldsymbol{T}^\textrm{T}\boldsymbol{K}\boldsymbol{T}\)を\(\boldsymbol{W}\)を用いて表すと次式のように対角行列が得られる。\begin{align} \boldsymbol{T}^\textrm{T}\boldsymbol{M}\boldsymbol{T}&=\boldsymbol{W}^\textrm{T}\boldsymbol{U}^{-\textrm{T}}\boldsymbol{M}\boldsymbol{U}^{-1}\boldsymbol{W}\\&=\boldsymbol{W}^\textrm{T}\boldsymbol{U}^{-\textrm{T}}\boldsymbol{U}^\textrm{T}\boldsymbol{U}\boldsymbol{U}^{-1}\boldsymbol{W}\\&=\boldsymbol{W}^\textrm{T}\boldsymbol{W}=\boldsymbol{\bar{M}}\\\boldsymbol{T}^\textrm{T}\boldsymbol{K}\boldsymbol{T}&=\boldsymbol{W}^\textrm{T}\boldsymbol{U}^{-\textrm{T}}\boldsymbol{K}\boldsymbol{U}^{-1}\boldsymbol{W}\\&=\boldsymbol{W}^\textrm{T}\boldsymbol{W}\boldsymbol{\Omega}^2=\boldsymbol{\bar{K}}\end{align}
• 上式から対角化した質量行列と剛性行列には次の関係式がある。\begin{align} \boldsymbol{\bar{K}}=\boldsymbol{\bar{M}}\boldsymbol{\Omega}^2\end{align}このことから、\(\boldsymbol{\Omega}^2\)は一般化固有値問題の固有値を対角成分に持つものであると同時に、その固有値の平方根は（無減衰）固有角周波数であることがわかる。

以上のことから、変換行列\(\boldsymbol{T}\)によって変換された変位ベクトル\(\boldsymbol{x}=\boldsymbol{T}\boldsymbol{\bar{x}}\)を用いると運動方程式は次式のように変換される。なお、この変換された変位ベクトル\(\boldsymbol{\bar{x}}\)をモード変位と呼ぶ。

\begin{align} \boldsymbol{\bar{M}}\cfrac{d^2\boldsymbol{\bar{x}}}{dt^2}+\boldsymbol{\bar{K}}\boldsymbol{\bar{x}}=\boldsymbol{T}^\textrm{T}\boldsymbol{f} \end{align}

対角化された質量行列と剛性行列の対角成分をモード質量、モード剛性と呼ばれる。上式を変換すると次式も得られる。

\begin{align} \cfrac{d^2\boldsymbol{\bar{x}}}{dt^2}+\boldsymbol{\Omega^2}\boldsymbol{\bar{x}}=\boldsymbol{\bar{M}}^{-1}\boldsymbol{T}^\textrm{T}\boldsymbol{f} \end{align}

通常、固有値行列\(\boldsymbol{\Omega^2}\)は左上から右下にかけて小さい順に並べ、最左上の固有値の平方根（固有角周波数）を1次モード固有角周波数と呼び、あとは順に2次3次と呼ばれる。それに対応した位置にあるモード質量、モード剛性もn次モード質量、n次モード剛性と呼ばれる。

例えば2自由度振動系の場合、左辺を分解すると次のように表すことができ、\(\boldsymbol{\bar{x}}\)の各要素は錬成していない形となる。これがモード分解である。

\begin{align} \cfrac{d^2}{dt^2}\left(\begin{array}{c}\bar{x}_1\\\bar{x}_2\end{array}\right)+\left[\begin{array}{cc}\omega_1^2&0\\0&\omega_2^2\end{array}\right]\left(\begin{array}{c}\bar{x}_1\\\bar{x}_2\end{array}\right)=\boldsymbol{\bar{M}}^{-1}\boldsymbol{T}^\textrm{T}\boldsymbol{f} \end{align}

比例減衰系\((\boldsymbol{C}=\alpha\boldsymbol{M}+\beta\boldsymbol{K})\)

減衰行列が質量行列や剛性行列に比例する場合、同じ手順でモード分解が可能である。同じ変換行列を用いると次式のようになる。

\begin{align} \boldsymbol{T}^\textrm{T}\boldsymbol{M}\boldsymbol{T}\cfrac{d^2\boldsymbol{\bar{x}}}{dt^2}+\boldsymbol{T}^\textrm{T}\boldsymbol{C}\boldsymbol{T}\cfrac{d\boldsymbol{\bar{x}}}{dt}+\boldsymbol{T}^\textrm{T}\boldsymbol{K}\boldsymbol{T}\boldsymbol{\bar{x}}=\boldsymbol{T}^\textrm{T}\boldsymbol{f}\\\boldsymbol{\bar{M}}\cfrac{d^2\boldsymbol{\bar{x}}}{dt^2}+\boldsymbol{T}^\textrm{T}\left(\alpha\boldsymbol{M}+\beta\boldsymbol{K}\right)\boldsymbol{T}\cfrac{d\boldsymbol{\bar{x}}}{dt}+\boldsymbol{\bar{K}}\boldsymbol{\bar{x}}=\boldsymbol{T}^\textrm{T}\boldsymbol{f}\\\boldsymbol{\bar{M}}\cfrac{d^2\boldsymbol{\bar{x}}}{dt^2}+\left(\alpha\boldsymbol{\bar{M}}+\beta\boldsymbol{\bar{K}}\right)\cfrac{d\boldsymbol{\bar{x}}}{dt}+\boldsymbol{\bar{K}}\boldsymbol{\bar{x}}=\boldsymbol{T}^\textrm{T}\boldsymbol{f}\\\cfrac{d^2\boldsymbol{\bar{x}}}{dt^2}+\left(\alpha\boldsymbol{I}+\beta\boldsymbol{\Omega}^2\right)\cfrac{d\boldsymbol{\bar{x}}}{dt}+\boldsymbol{\Omega}^2\boldsymbol{\bar{x}}=\boldsymbol{\bar{M}}^{-1}\boldsymbol{T}^\textrm{T}\boldsymbol{f} \end{align}

先ほどと同様、2自由度振動系の場合、左辺を分解すると次のように表すことができ、\(\boldsymbol{\bar{x}}\)の各要素は錬成していない形となる。

\begin{align} \cfrac{d^2}{dt^2}\left(\begin{array}{c}\bar{x}_1\\\bar{x}_2\end{array}\right)+\left[\begin{array}{cc}\alpha+\beta\omega_1^2&0\\0&\alpha+\beta\omega_2^2\end{array}\right]\cfrac{d}{dt}\left(\begin{array}{c}\bar{x}_1\\\bar{x}_2\end{array}\right)+\left[\begin{array}{cc}\omega_1^2&0\\0&\omega_2^2\end{array}\right]\left(\begin{array}{c}\bar{x}_1\\\bar{x}_2\end{array}\right)=\boldsymbol{\bar{M}}^{-1}\boldsymbol{T}^\textrm{T}\boldsymbol{f} \end{align}

しかし、この減衰行列が質量行列や剛性行列に比例するという仮定を満たすことは実際にはない。そのため通常は\(\boldsymbol{T}^\textrm{T}\boldsymbol{C}\boldsymbol{T}\)に非対角項（連成項）が0とならない。しかし、多くの場合でこの非対角項は対角項と比べて小さくなるため、モード分解を行った後、非対角項を0に近似するのが実用的である。

モード形式(Eigenstructure decomposition)

上記のように減衰行列が質量行列や剛性行列に比例しない場合、具体的には意図的に減衰を付加した場合など減衰が大きい系に対しては、完全にモード分解ができない。そこで異なるモード分解の手法としてモード形式について説明する。モード形式という名前はメカトロニクス屋さん御用達のソフトウェア、MATLABで扱っている名称であるが、あまり馴染みがない。英語ではEigenstructure decompositionという名称のようだが、こちらもそこまでメジャーではない。

正準状態空間実現 - MATLAB canon - MathWorks 日本

モード形式を用いるには、まず運動方程式を1階の微分方程式（状態方程式）に変換する必要がある。位置ベクトルの一階微分である速度ベクトル\(\boldsymbol{v}=\frac{d}{dt}\boldsymbol{x}\)を用いると次式のようになる。 

\begin{align} \left[\begin{array}{cc}\boldsymbol{I}&\boldsymbol{0}\\\boldsymbol{0}&\boldsymbol{M}\end{array}\right]\cfrac{d}{dt}\left(\begin{array}{c}\boldsymbol{x}\\\boldsymbol{v}\end{array}\right)=\left[\begin{array}{cc}\boldsymbol{0}&\boldsymbol{I}\\-\boldsymbol{K}&-\boldsymbol{C}\end{array}\right]\left(\begin{array}{c}\boldsymbol{x}\\\boldsymbol{v}\end{array}\right)+\left(\begin{array}{c}\boldsymbol{0}\\\boldsymbol{f}\end{array}\right) \end{align}

次の一般化固有値問題の固有値と固有ベクトルを求める。

\begin{align} \left[\begin{array}{cc}\boldsymbol{I}&\boldsymbol{0}\\\boldsymbol{0}&\boldsymbol{M}\end{array}\right]\boldsymbol{V}\boldsymbol{\Lambda}=\left[\begin{array}{cc}\boldsymbol{0}&\boldsymbol{I}\\-\boldsymbol{K}&-\boldsymbol{C}\end{array}\right]\boldsymbol{V} \end{align}

各固有値のうち、実数固有値\(\lambda_{1,2,...}\)に対応した固有ベクトルを\(\boldsymbol{v}^r_{1,2,..}\)、複素固有値\(-\zeta_{1,2,...}\omega_{1,2,...}\pm\omega_{1,2,...}\sqrt{1-\zeta_{1,2,...}^2}\)に対応した固有ベクトルを\(\boldsymbol{v}^c_{1,2,..}\)と表したとき、次の変換行列\(\boldsymbol{T}\)と新しい状態ベクトル\(\boldsymbol{\phi}, \boldsymbol{\psi}\)を考える。

\begin{align}\boldsymbol{\Lambda}&=\textrm{diag}\left(\lambda_1,\lambda_2,...,-\zeta_1\omega_1\pm\omega_1\sqrt{1-\zeta_1^2},-\zeta_2\omega_2\pm\omega_2\sqrt{1-\zeta_2^2},...\right)\\\boldsymbol{T}&=\left[\boldsymbol{v}^r_1,\boldsymbol{v}^r_2,...,\textrm{Re}\left(\boldsymbol{v}^c_1\right),\textrm{Im}\left(\boldsymbol{v}^c_1\right),\textrm{Re}\left(\boldsymbol{v}^c_2\right),\textrm{Im}\left(\boldsymbol{v}^c_2\right),...\right]\\\left(\begin{array}{c}\boldsymbol{x}\\\boldsymbol{v}\end{array}\right)&=\boldsymbol{T}\left(\begin{array}{c}\boldsymbol{\phi}\\\boldsymbol{\psi}\end{array}\right)\end{align}

変換行列\(\boldsymbol{T}\)と状態ベクトル\(\boldsymbol{\phi}, \boldsymbol{\psi}\)を用いると状態方程式は次式のようになる。

\begin{align} \cfrac{d}{dt}\left(\begin{array}{c}\boldsymbol{\phi}\\\boldsymbol{\psi}\end{array}\right)&=\boldsymbol{T}^{-1}\left[\begin{array}{cc}\boldsymbol{0}&\boldsymbol{I}\\-\boldsymbol{M}^{-1}\boldsymbol{K}&-\boldsymbol{M}^{-1}\boldsymbol{C}\end{array}\right]\boldsymbol{T}\left(\begin{array}{c}\boldsymbol{\phi}\\\boldsymbol{\psi}\end{array}\right)+\boldsymbol{T}^{-1}\left(\begin{array}{c}\boldsymbol{0}\\\boldsymbol{M}^{-1}\boldsymbol{f}\end{array}\right)\\&=\left[\begin{array}{cc}\boldsymbol{\Lambda}&\boldsymbol{0}\\\boldsymbol{0}&\boldsymbol{D}\end{array}\right]\left(\begin{array}{c}\boldsymbol{\phi}\\\boldsymbol{\psi}\end{array}\right)+\left(\begin{array}{c}\boldsymbol{f}\\\boldsymbol{g}\end{array}\right) \end{align}

\begin{align} \boldsymbol{\Lambda}&=\textrm{diag}\left(\lambda_1,\lambda_2,...\right)\\\boldsymbol{D}&=\textrm{blkdiag}\left(\boldsymbol{D}_1,\boldsymbol{D}_2,...\right)\\\boldsymbol{D}_i&=\left[\begin{array}{cc}-\zeta_i\omega_i&\omega_i^d\\-\omega_i^d&-\zeta_i\omega_i\end{array}\right]\\\omega_i^d&=\omega_i\sqrt{1-\zeta_i^2} \end{align}

上式のように、実数固有値の成分\(\boldsymbol{\Lambda}\)と複素固有値の成分\(\boldsymbol{D}\)とに分けることができ、複素固有値の成分も実部と虚部で分けて表すことができる。\(\zeta_i\)と\(\omega_i\)はそれぞれモード減衰比、モード固有角周波数である。そしてそれぞれのモードは他のモードと分離されている。上式の行列形式を分解すると下記のようになる。

\begin{align} \cfrac{d\phi_i}{dt}&=\lambda_i\phi_i+f_i\\\cfrac{d\psi_{2j-1}}{dt}&=-\zeta_j\omega_j\psi_{2j-1}+\omega_j^d\psi_{2j}+g_{2j-1}\\\cfrac{d\psi_{2j}}{dt}&=-\zeta_j\omega_j\psi_{2j}-\omega_j^d\psi_{2j-1}+g_{2j}\end{align}

\(\psi_{2j-1}\)と\(\psi_{2j}\)が混在しているが、1つにまとめると次式のように表せる。

\begin{align} \cfrac{d^2\psi_{2j-1}}{dt^2}+2\zeta_j\omega_j\cfrac{d\psi_{2j-1}}{dt}+\omega_j^2\psi_{2j-1}&=\cfrac{dg_{2j-1}}{dt}+\zeta_j\omega_jg_{2j-1}+\omega_j^dg_{2j}\\\cfrac{d^2\psi_{2j}}{dt^2}+2\zeta_j\omega_j\cfrac{d\psi_{2j}}{dt}+\omega_j^2\psi_{2j}&=\cfrac{dg_{2j}}{dt}+\zeta_j\omega_jg_{2j}-\omega_j^dg_{2j-1}\end{align}

まとめ

基本的にモード分解を必要とするモード解析の分野では、減衰が付加されていない無減衰系を対象とする場合が多いため、前半で示した比例減衰系を仮定したモード分解で充分である。

しかし、減衰が付加されているような系では単純にモード分解が行えないため、後半で示したモード形式のように分解することで、厳密なモード減衰比を求めることができるため有用である。

１自由度振動系のお話 その４

１自由度振動系のお話 その１ - メカトロニックなカメ（無次元化について）

１自由度振動系のお話 その２ - メカトロニックなカメ（周波数応答・過渡応答について）

１自由度振動系のお話 その３ - メカトロニックなカメ（状態方程式について）

の続き

復習

質量を\(m\)[kg]、ばね定数を\(k\)[N/m]、減衰係数を\(c\)[Ns/m]、位置を\(x(t)\)[m]、外力を\(f(t)\)[N]とすると、一自由度マスばねダンパ系の運動方程式は下記のようになります。

\begin{align} m\cfrac{d^2x}{dt^2}+c\cfrac{dx}{dt}+kx(t)=f(t) \end{align}

 上式を無次元化した運動方程式は次式のようになります。

\begin{align} \cfrac{d^2\tilde{x}}{d\tilde{t}^2}+2\zeta\cfrac{d\tilde{x}}{d\tilde{t}}+\tilde{x}(\tilde{t})=\tilde{f}(\tilde{t}) \end{align} 

回転運動に関する運動方程式(Equation of rotational motion)

ここまでは直線運動に関する運動方程式に着目してきました。次に回転運動、つまり回転したり、ねじれたりする運動に着目します。回転運動に関する一自由度振動系は次のように表されます。

\begin{align} J\cfrac{d^2\theta}{dt^2}+D\cfrac{d\theta}{dt}+K\theta(t)=N(t) \end{align}

慣性モーメントを\(J\)[kgm\(^2\)]、ねじり剛性係数を\(K\)[Nm/rad]、ねじり減衰係数を\(D\)[Nms/rad]、角度を\(\theta(t)\)[rad]、外力モーメントを\(N(t)\)[Nm]を表します。

さて、

１自由度振動系のお話 その１ - メカトロニックなカメと上記の運動方程式を同じ手順で無次元化を行います。

まず初めに角度\(\theta(t)\)を無次元化します。といっても角度はすでに無次元化されています。なぜなら角度の単位[rad]は、円の半径とその円弧の長さの比であるからである。ここが直線運動の運動方程式と大きく違う点である。

次に外力モーメントの無次元化を行います。ここで角度単位の[rad]は無次元単位であるため、ねじり剛性係数と外力モーメントは同じ次元である。そのため、無次元化された外力モーメント\(\tilde{N}(t)=N(t)/K\)[rad]を用いると、運動方程式は次式のようになる。

\begin{align} \cfrac{J}{K}\cfrac{d^2\theta}{dt^2}+\cfrac{D}{K}\cfrac{d\theta}{dt}+\theta(t)=\tilde{N}(t) \end{align}

最後に直線運動と同様に時間の無次元化を行います。固有角周波数\(\omega\)[rad/sec]は\(\sqrt{K/J}\)で表されるので、無次元化された時間\(\tilde{t}=\omega t\)を用いると回転運動に関する運動方程式は次式のようになります。

\begin{align}\cfrac{d^2\theta}{d\tilde{t}^2}+\cfrac{D}{\sqrt{JK}}\cfrac{d\theta}{d\tilde{t}}+\theta(\tilde{t})=\tilde{N}(\tilde{t}) \end{align}

減衰比を\(\zeta=\cfrac{D}{2\sqrt{JK}}\)とすると、下記のように直線運動と同じ運動方程式が得られます。

\begin{align}\cfrac{d^2\theta}{d\tilde{t}^2}+2\zeta\cfrac{d\theta}{d\tilde{t}}+\theta(\tilde{t})=\tilde{N}(\tilde{t}) \end{align}

さて、ここで気にしたい点として両辺の単位が[rad]となっていることである。このままでは外力モーメントの単位が[rad]であるため、直感的にわかりづらい。そこで正規化を行う。正規化とは規格化とも呼ばれ、運動方程式などの数式（数理モデル）においては、ある変数（定数）を基準となる量で割ることで、新しい変数（定数）が１を基準とした運動方程式に変形される。１を基準とするため数値を用いて計算する際も非常にわかりやすくなるというメリットがある。実は直線運動の運動方程式でも変数\(x\)を基準量\(x_{st}\)で割っているため、正規化を行っていたということである。

正規化にはまず基準角度を決める必要がある。通常は基準角度として、この場合では使用する角度範囲の最大値を基準角度とすることが多い。ここで基準の角度を\(\theta_0\)としたとき、正規化された回転運動の運動方程式は次のようになる。

\begin{align}\cfrac{d^2\tilde{\theta}}{d\tilde{t}^2}+2\zeta\cfrac{d\tilde{\theta}}{d\tilde{t}}+\tilde{\theta}(\tilde{t})=\tilde{N_0}(\tilde{t}) \end{align}

正規化角度が\(\tilde{\theta}(\tilde{t})=\theta(\tilde{t})/\theta_0\)、正規化外力モーメントが\(\tilde{N_0}(\tilde{t})=N(\tilde{t})/\theta_0\)と表される。これで直線運動の運動方程式と同じ方程式が得られた。

電気系の回路方程式(Circuit equation)

メカトロニクスを考える上で、機械系の運動方程式だけを考えるだけでは不十分である。ここでは電気系の回路方程式に着目する。今回はもっとも簡単なRLS直列回路、つまり抵抗とインダクタとコンデンサが直列につながれた回路に対し、電圧を印加するモデルを考えます。この時の回路方程式は次のようになります。

\begin{align} L\cfrac{di}{dt}+Ri(t)+\cfrac{1}{C}\int i dt=v(t) \end{align}

インダクタンスを\(L\)[H]、抵抗値を\(R\)[\(\Omega\)]、キャパシタンスを\(C\)[F]、電流を\(i(t)\)[A]、電圧を\(v(t)\)[V]を表します。この方程式は積分微分方程式であり、若干扱いが難しい。そこで、コンデンサに貯まる電荷\(q=\int i dt\)を用いることで、次式のように微分方程式が得られる。

\begin{align} L\cfrac{d^2 q}{dt^2}+R\cfrac{d q}{dt}+\cfrac{1}{C}q(t)=v(t) \end{align}

こうなると運動方程式と同じ形であることがわかる。基準電荷\(q_0\)、無次元化電荷\(\tilde{q}(t)=q(t)/q_0\)、無次元化電圧\(\tilde{v}(t)=Cv(t)/q_0\)、固有角周波数\(1/\sqrt{LC}\)、無次元化時間\(\tilde{t}=\omega t\)、減衰比\(\zeta = R\sqrt{\frac{C}{L}}\)を用いると、次のような回路方程式が得られる。

\begin{align} \cfrac{d^2 \tilde{q}}{d\tilde{t}^2}+2\zeta\cfrac{d \tilde{q}}{d\tilde{t}}+\tilde{q}(\tilde{t})=\tilde{v}(\tilde{t}) \end{align}

もちろん無次元化された回路方程式も運動方程式と同じ形状となる。

まとめ

このように、機械系の運動方程式と電気系の回路方程式は同じ方程式で表される。言い換えると機械系の運動を回路で表現することができ、逆に電気系の応答を機械で表現することができる。このことを日本語ではアナロジーと呼ばれ、英語ではMechanical–electrical analogiesと呼ばれている。

具体的には上記の数式を用いると、インダクタは質量、抵抗は減衰係数、コンデンサはばね定数の逆数、位置は電荷、速度は電流、外力は電圧に相当する。

この現象は非常に興味深く、電気を専門としている人は機械系を電気回路に変換してしまえば、慣れている領域で考えることができるため、扱いやすくなる。逆もまた然り。メカトロニクスを考えるうえで非常に重要な性質である。

しかしここで挙げた例では線形な要素のみであり、すべての要素に対して機械電気アナロジーを適用できるわけではない。例えばトランジスタのような半導体を機械で表すことは困難である。そのため、あくまでこの機械電気アナロジーは、電気や機械を専門としている者の橋渡しをする役目であり、メカトロニクス技術者にとっては万能な方法論ではないことを注意しておこう。

１自由度振動系のお話 その３

１自由度振動系のお話 その１ - メカトロニックなカメ（無次元化について）

１自由度振動系のお話 その２ - メカトロニックなカメ（周波数応答・過渡応答について）

の続き

復習

質量を\(m\)[kg]、ばね定数を\(k\)[N/m]、減衰係数を\(c\)[Ns/m]、位置を\(x(t)\)[m]、外力を\(f(t)\)[N]とすると、一自由度マスばねダンパ系の運動方程式は下記のようになります。

\begin{align} m\cfrac{d^2x}{dt^2}+c\cfrac{dx}{dt}+kx(t)=f(t) \end{align}

 上式を無次元化した運動方程式は次式のようになります。

\begin{align} \cfrac{d^2\tilde{x}}{d\tilde{t}^2}+2\zeta\cfrac{d\tilde{x}}{d\tilde{t}}+\tilde{x}(\tilde{t})=\tilde{f}(\tilde{t}) \end{align} 

この系の過渡応答には次式のようになります。

 \begin{align}\tilde{x}(\tilde{t})=C_1e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_2e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\end{align}

\(\zeta<1\)の時、この系の過渡応答には次式のように振動的な応答が発生します。

 \begin{align}\tilde{x}(\tilde{t})=X'e^{-\zeta\tilde{t}}\sin\left(\sqrt{1-\zeta^2}\tilde{t}+\Theta'\right)\end{align}

この過渡応答における振動の周波数比\(\alpha_{damp}\)は無次元化減衰固有周波数と呼ばれ次式で表される。

  \begin{align}\alpha_{damp}=\sqrt{1-\zeta^2}\end{align}

状態方程式(State equation)

制御工学において、高次の微分方程式を1階の微分方程式の形で表した式を状態方程式と呼ばれ、非常に基礎的な考え方である。また、非線形な挙動を含む複雑なシステムは一般に厳密な解を求めることができないため、ソフトウェアを用いたシミュレーションの分野では、複雑なシステムを1階の微分方程式に変換する。そしてルンゲ=クッタ法等の数値積分手法を用いて、シミュレーションを行う。というように1階の微分方程式に変換することは非常に重要であるので、その手順を簡単に示す。

まずは、一自由度マスばねダンパ系の2階の運動方程式を1階の微分方程式に変換するために、次のような新しい変数\(\tilde{v}(\tilde{t})\)を導入します。

 \begin{align}\tilde{v}(\tilde{t})&=\cfrac{d\tilde{x}}{d\tilde{t}}\end{align}

上式を見てわかる通り、\(\tilde{v}(\tilde{t})\)は単なる無次元化速度である。この無次元化速度を用いると、運動方程式は次式のようになる。

\begin{align} \cfrac{d\tilde{v}}{d\tilde{t}}+2\zeta\tilde{v}(\tilde{t})+\tilde{x}(\tilde{t})=\tilde{f}(\tilde{t}) \end{align} 

これで1階の微分方程式になりました。しかしこの方程式には2つの変数\(\tilde{v}, \tilde{x}\)があり、多変数な微分方程式となってしまい扱いが困難になります。そこで、また新しい変数のベクトル\(\tilde{\boldsymbol{w}}\)を次式のように用意します。

 \begin{align}\tilde{\boldsymbol{w}}(\tilde{t})&=\left(\begin{array}{c} \tilde{x}(\tilde{t}) \\\tilde{v}(\tilde{t})\end{array}\right)\end{align}

この新しい変数ベクトル（状態変数ともいう）を用いると、次式のような状態方程式となります。

\begin{align} \cfrac{d\tilde{\boldsymbol{w}}}{d\tilde{t}}=\left[\begin{array}{cc} 0&1\\-1&-2\zeta\end{array}\right]\tilde{\boldsymbol{w}}(\tilde{t})+\left(\begin{array}{c} 0 \\\tilde{f}(\tilde{t})\end{array}\right)\end{align}

 さて、この状態方程式の行列部の固有値を求めてみましょう。固有値の求め方は、次式を満たすような\(\lambda\)を求めることです。

\begin{align} \left[\begin{array}{cc} 0&1\\-1&-2\zeta\end{array}\right]\boldsymbol{z}=\lambda\boldsymbol{z}\end{align}

求められた固有値\(\lambda_1,\lambda_2\)は次式のようになる。

\begin{align} \lambda_{1,2}=-\zeta\pm\sqrt{\zeta^2-1}\end{align}

この値は、まさに過渡応答における指数部に相当します。

対角化(Diagonalization)

固有値には非常に重要な性質を持つことが分かった。そこで対角化して固有値を抜き出す作業を行ってみる。まずは固有値に対応した固有ベクトルを求めてよう。固有ベクトルは次式を満たす\(\boldsymbol{z}_0\)である。

\begin{align} \left[\begin{array}{cc} 0&1\\-1&-2\zeta\end{array}\right]\boldsymbol{z}_0=\lambda_{1,2}\boldsymbol{z}_0\end{align}

\begin{align} \boldsymbol{z}_{1,2}=\left(\begin{array}{c} 1\\\lambda_{1,2}\end{array}\right)=\left(\begin{array}{c} 1\\-\zeta\pm\sqrt{\zeta^2-1}\end{array}\right)\end{align}

 状態方程式の対角化のために、固有ベクトルを横に並べた変換行列\(\boldsymbol{A}\)を導入します。

\begin{align} \boldsymbol{A}&=\left[\begin{array}{cc} 1&1\\\lambda_1&\lambda_2\end{array}\right]=\left[\begin{array}{cc} 1&1\\-\zeta+\sqrt{\zeta^2-1}&-\zeta-\sqrt{\zeta^2-1}\end{array}\right]\\\boldsymbol{A}^{-1}&=-\cfrac{1}{\sqrt{\zeta^2-1}}\left[\begin{array}{cc} \lambda_2&-1\\-\lambda_1&1\end{array}\right]\end{align}

対角化のために状態変数を次のように変換する。

\begin{align} \tilde{\boldsymbol{w}}(\tilde{t})=\boldsymbol{A}\tilde{\boldsymbol{w}}_d(\tilde{t})\end{align}

この状態変数を用いると状態方程式は次のようになります。

\begin{align} \cfrac{d\tilde{\boldsymbol{w}}_d}{d\tilde{t}}&=\boldsymbol{A}^{-1}\left[\begin{array}{cc} 0&1\\-1&-2\zeta\end{array}\right]\boldsymbol{A}\tilde{\boldsymbol{w}}_d(\tilde{t})+\boldsymbol{A}^{-1}\left(\begin{array}{c} 0 \\\tilde{f}(\tilde{t})\end{array}\right)\\&=\left[\begin{array}{cc} \lambda_1&0\\0&\lambda_2\end{array}\right]\tilde{\boldsymbol{w}}_d(\tilde{t})-\cfrac{1}{\sqrt{\zeta^2-1}}\left(\begin{array}{c} -1 \\1\end{array}\right)\tilde{f}(\tilde{t})\end{align}

上式は固有値でまとめることができて、非常にすっきりした形となったといえる。しかし、減衰比\(\zeta\)が1より小さい場合、各係数に虚数や複素数が出てきてしまい、扱いが非常に複雑になる。そこで対角化ではなくモード形式への変換を行ってみよう。

モード形式（Mode form）

減衰比が1未満\(\zeta<1\)の時、固有値は次のように表すこともできる。

\begin{align} \lambda_{1,2}=-\zeta\pm i \alpha_{damp}\end{align}

次のような変換行列\(\boldsymbol{B}\)を導入します。 

\begin{align} \boldsymbol{B}&=\left[\begin{array}{cc} \textrm{Re}(\boldsymbol{z}_{1,2})&\textrm{Im}(\boldsymbol{z}_{1,2})\end{array}\right]\\&=\left[\begin{array}{cc} 1&0\\-\zeta&\alpha_{damp}\end{array}\right]\\\boldsymbol{B}^{-1}&=\cfrac{1}{\alpha_{damp}}\left[\begin{array}{cc} \alpha_{damp}&0\\\zeta&1\end{array}\right]\end{align}

先ほどと同様の手順で状態変数を次のように変換する。

\begin{align} \tilde{\boldsymbol{w}}(\tilde{t})=\boldsymbol{B}\tilde{\boldsymbol{w}}_m(\tilde{t})\end{align}

この状態変数を用いると状態方程式は次のようになります。

\begin{align} \cfrac{d\tilde{\boldsymbol{w}}_m}{d\tilde{t}}&=\boldsymbol{B}^{-1}\left[\begin{array}{cc} 0&1\\-1&-2\zeta\end{array}\right]\boldsymbol{B}\tilde{\boldsymbol{w}}_m(\tilde{t})+\boldsymbol{B}^{-1}\left(\begin{array}{c} 0 \\\tilde{f}(\tilde{t})\end{array}\right)\\&=\left[\begin{array}{cc} -\zeta&\alpha_{damp}\\-\alpha_{damp}&-\zeta\end{array}\right]\tilde{\boldsymbol{w}}_m(\tilde{t})+\cfrac{1}{\alpha_{damp}}\left(\begin{array}{c} 0 \\1\end{array}\right)\tilde{f}(\tilde{t})\end{align}

上式のように1自由度振動系で重要な性質を表す減衰比\(\zeta\)と無次元化減衰固有周波数\(\alpha_{damp}\)を抜き出すことに成功しました。

まとめ

以上のように、1自由度振動系を表す2階の微分方程式を1階の微分方程式に置き換える方法にはいくつかある。最後にこれらの特徴をまとめると、

• 変換なし：元の運動方程式から求めやすい。\(\tilde{\boldsymbol{w}}(\tilde{t})\)に物理的意味（位置、速度）があるため、わかりやすい。1階の微分方程式の形が直感にそぐわない
• 対角化：固有値を用いた対角行列が得られてすっきりする。\(\tilde{\boldsymbol{w}}_d(\tilde{t})\)に物理的意味がない。減衰比が1未満の時は虚数が現れる。
• モード形式：減衰比が1未満の時に、減衰比と無次元化減衰固有周波数に分けて表記できて、見通しが良い。\(\tilde{\boldsymbol{w}}_m(\tilde{t})\)に物理的意味がない。

１自由度振動系のお話 その２

１自由度振動系のお話 その１ - メカトロニックなカメの続き

復習

質量を\(m\)[kg]、ばね定数を\(k\)[N/m]、減衰係数を\(c\)[Ns/m]、位置を\(x(t)\)[m]、外力を\(f(t)\)[N]とすると、一自由度マスばねダンパ系の運動方程式は下記のようになります。

\begin{align} m\cfrac{d^2x}{dt^2}+c\frac{dx}{dt}+kx(t)=f(t) \end{align}

 上式を無次元化した運動方程式は次式のようになります。

\begin{align} \cfrac{d^2\tilde{x}}{d\tilde{t}^2}+2\zeta\frac{d\tilde{x}}{d\tilde{t}}+\tilde{x}(\tilde{t})=\tilde{f}(\tilde{t}) \end{align}

 この運動方程式に正弦波の外力\(\tilde{f}_0\sin(\alpha\tilde{t})\)を与えたときの応答は次式のようになります。

 \begin{align}\tilde{x}(\tilde{t})=&\cfrac{\tilde{f}_0}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\sin\left(\alpha\tilde{t}+\cos^{-1}\left(\cfrac{1-\alpha^2}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\right)\right)\\&+C_1e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_2e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\end{align}

周波数応答(Frequency response)

今回は右辺第一項の周期的な成分にのみ着目します。というのも第二項、第三項に含まれる\(-(\zeta\pm\sqrt{\zeta^2-1})\tilde{t}\)の実部は必ず負であるため、\(e^{-A\tilde{t}}(A>0)\)は時間とともに0に収束します。よって、十分に長い時間周期的な外乱を与えたときは第一項しか残りません。

 \begin{align}\tilde{x}(\tilde{t})&=\cfrac{\tilde{f}_0}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\sin\left(\alpha\tilde{t}+\cos^{-1}\left(\cfrac{1-\alpha^2}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\right)\right)\\&=\tilde{f}_0X(\alpha)\sin\left(\alpha\tilde{t}+\Theta(\alpha)\right)\end{align}

この式は元の外乱\(\tilde{f}_0\sin(\alpha\tilde{t})\)に対して、振幅が\(X(\alpha)\)倍、位相が\(\Theta(\alpha)\)だけずれていることを意味します。\(X(\alpha)\)の部分は振幅倍率とも呼ばれます。何に対しての倍率かというと、周波数がゼロの外乱、つまり\(\alpha=0\)で定数な外乱を与えたときの位置\(\tilde{x}\)を1としたとき倍率です。

この振幅倍率と位相差を無次元化周波数\(\alpha\)の周波数応答と呼び、図１のように振幅倍率と位相差を縦に並べたものを、制御工学ではボード(Bode)線図と呼びます。

ボード線図\((\zeta = 0.1)\)

 図１より、振幅倍率が無次元化周波数が１付近で極大値を取っています。次にこの振幅倍率が最大となる無次元化周波数を求めてみましょう。

振幅倍率\(X(\alpha)\)が最大となる条件は\((1-\alpha^2)^2+4\zeta^2\alpha^2\)が最小となるときです。つまり振幅倍率が最大となる無次元化周波数\(\alpha_{xpeak}\)は次式のように表されます。

 \begin{align}\left.\cfrac{d}{d\alpha}\left((1-\alpha^2)^2+4\zeta^2\alpha^2\right)\right|_{\alpha=\alpha_{xpeak}}&=0\\\left.\cfrac{d\alpha^2}{d\alpha}\cfrac{d}{d\alpha^2}\left((1-\alpha^2)^2+4\zeta^2\alpha^2\right)\right|_{\alpha=\alpha_{xpeak}}&=0\\2\alpha_{xpeak}\left(-2(1-\alpha_{xpeak}^2)+4\zeta^2\right)&=0\\\alpha_{xpeak}^2&=1-2\zeta^2\end{align}

つまり\(\zeta<1/\sqrt{2}\approx0.707\)の時、振幅倍率\(X(\alpha)\)は次式の無次元化周波数で極大値を持ち、その時の振幅倍率\(X(\alpha_{xpeak})\)と位相差\(\Theta(\alpha_{xpeak})\)は次式のように得られる。

 \begin{align}\alpha_{xpeak}&=\sqrt{1-2\zeta^2}\\X(\alpha_{xeak})&=\cfrac{1}{2\zeta\sqrt{(1-\zeta^2)}}\\\Theta(\alpha_{xpeak})&=\cos^{-1}\left(\cfrac{\zeta}{\sqrt{1-\zeta^2}}\right)\end{align}

また、位相差\(\Theta(\alpha)\)が90degとなるような、無次元化周波数は\(\alpha=1\)の時であることは、簡単に確認することができる。つまり、位相差が90degずれるのは、振動系の固有周波数のときのみである。

このように、振幅倍率が極大となる周波数は一般に固有周波数と異なり、位相差が90degずれる周波数は固有周波数と一致するといった特徴が見受けられた。

過渡応答(Transient response)

先ほどまでは、一定周波数の外乱を与えたときの定常応答を観察しました。今度は定常応答以外の成分、過渡応答に着目します。過渡応答成分のみ抜き出すと次式のようになります。

 \begin{align}\tilde{x}(\tilde{t})=C_1e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_2e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\end{align}

\(\zeta>1\)の場合、この過渡応答は単なる指数関数の重ね合わせになるため、振動的にならない。逆に\(\zeta<1\)の場合、指数部に虚数が現れるためオイラーの等式を用いて整理すると下記のように表されます。

 \begin{align}\tilde{x}(\tilde{t})&=e^{-\zeta\tilde{t}}\left(C_1'\sin\left(\sqrt{1-\zeta^2}\tilde{t}\right)+C_2'\cos\left(\sqrt{1-\zeta^2}\tilde{t}\right)\right)\\&=X'e^{-\zeta\tilde{t}}\sin\left(\sqrt{1-\zeta^2}\tilde{t}+\Theta'\right)\end{align}

このように、無次元化周波数\(\alpha_{damp}=\sqrt{1-\zeta^2}\)で振動しながら、\(e^{-\zeta\tilde{t}}\)の減衰率で過渡応答が0に収束していくことがわかる。この減衰しながら振動する周波数を（無次元化された）減衰固有周波数と呼ばれます。この時の振幅倍率\(X(\alpha_{damp})\)と位相差\(\Theta(\alpha_{damp})\)は次式のように得られる。

  \begin{align}X(\alpha_{damp})&=\cfrac{1}{\zeta\sqrt{(4-3\zeta^2)}}\\\Theta(\alpha_{damp})&=\cos^{-1}\left(\cfrac{\zeta}{\sqrt{4-3\zeta^2}}\right)\end{align}

少し脱線するが、 この振動的な過渡応答が1周期ごとにどれだけ振幅が減衰するかを表す指標に対数減衰率というものがある。過渡応答の1周期は\(T=1/\sqrt{1-\zeta^2}\)であるため、\(\tilde{t}=\tilde{t}_0\)の時の振幅と\(\tilde{t}=\tilde{t}_0+T\)の時の振幅の比の自然対数を対数減衰率\(\delta\)と呼び、次式のように表せる。

 \begin{align}\delta&=\ln\left(\cfrac{\tilde{x}(\tilde{t}_0)}{\tilde{x}(\tilde{t}_0+T)}\right)\\&=\ln\left(\cfrac{X'e^{-\zeta\tilde{t}_0}\sin\left(\sqrt{1-\zeta^2}\tilde{t}_0+\Theta'\right)}{X'e^{-\zeta(\tilde{t}_0+T)}\sin\left(\sqrt{1-\zeta^2}(\tilde{t}_0+T)+\Theta'\right)}\right)\\&=\cfrac{\zeta}{\sqrt{1-\zeta^2}}\end{align}

この対数減衰率は実験的に求めやすく、減衰が小さい振動系の場合\(\delta\approx\zeta\)とみなすことができるため、減衰比の推定によく用いられる。

速度・加速度の場合

最後に位置の一階微分である速度、二階微分である加速度についても着目してみる。周波数応答成分は次式のようになる。

 \begin{align}\cfrac{d\tilde{x}}{d\tilde{t}}&=\tilde{f}_0\alpha X(\alpha)\cos\left(\alpha\tilde{t}+\Theta(\alpha)\right)\\\cfrac{d^2\tilde{x}}{d\tilde{t}^2}&=-\tilde{f}_0\alpha^2 X(\alpha)\sin\left(\alpha\tilde{t}+\Theta(\alpha)\right)\end{align}

微分するごとに振幅倍率\(X(\alpha)\)に無次元化周波数\(\alpha\)が掛けられていき、位相差は変化しない。また過渡応答成分は次式のようになる。

 \begin{align}\cfrac{d\tilde{x}}{d\tilde{t}}&=C_3e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_4e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\\\cfrac{d^2\tilde{x}}{d\tilde{t}^2}&=C_5e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_6e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\end{align}

過渡応答成分は初期値に関する定数が変化するのみで、時間に関する項は変化しない、つまり振動の周期は変化しないことがわかる。

微分することで振幅倍率が変わるため、振幅倍率が最大となる周波数も変化する。計算の過程は省略するが、速度の振幅倍率が最大となる無次元周波数\(\alpha_{vpeak}\)と加速度の振幅倍率が最大となる無次元化周波数\(\alpha_{apeak}\)は次式のようになる。

 \begin{align}\alpha_{vpeak}&=1\\\alpha_{apeak}&=\cfrac{1}{\sqrt{1-2\zeta^2}}\end{align}

まとめ

以上のことより、一自由度振動系の無次元化された（無減衰）固有周波数\(\alpha_0\)と、周波数応答における位置、速度、加速度に関する振幅倍率が最大となる無次元化周波数\(\alpha_{xpeak}\),\(\alpha_{vpeak}\),\(\alpha_{apeak}\)、無次元化された減衰固有周波数\(\alpha_{damp}\)には、次の関係式が存在する。

 \begin{align}\alpha_{xpeak}\leq\alpha_{damp}\leq\alpha_0=\alpha_{vpeak}\leq\alpha_{apeak}\\\sqrt{1-2\zeta^2}\leq\sqrt{1-\zeta^2}\leq1=1\leq\cfrac{1}{\sqrt{1-2\zeta^2}}\end{align}

同じ一自由度振動系の中にも4つの特徴的な周波数が見られるといった不思議な現象が発見できました。

１自由度振動系のお話 その１

振動工学や制御工学を学ぶ際に必ず出てくる制御対象として、マスばねダンパで構成された1自由度振動系。今回はここを深堀してみます。

質量を\(m\)[kg]、ばね定数を\(k\)[N/m]、減衰係数を\(c\)[Ns/m]、位置を\(x(t)\)[m]、外力を\(f(t)\)[N]とすると、下記のような運動方程式が得られます。

\begin{align} m\cfrac{d^2x}{dt^2}+c\frac{dx}{dt}+kx(t)=f(t) \end{align}

今回はこの式を使っていろいろやっていきます。

無次元化(Nondimensionalization)

 無次元化とは、例えば位置には[m]の次元を持っていますが、この次元をなくす（無次元にする）ことです。無次元化のメリットは

• 式に含まれるパラメータの数が減って見通しが良くなる
• 数値計算の際に桁落ちが生じづらい
• 例えば船や飛行機の試験では、大きすぎて実機を使うわけにはいかないので、無次元化した方程式が同じになるような条件にすることで、実機を模擬した試験が可能になる
• 具体的な質量や位置など有次元の情報は、企業秘密も含まれているため、それを隠して学会や特許等で公開するために使う

といったところです。非常に実用的なのですが、学生はあまり意識することが少ないかと思います。覚えておいて損はないでしょう。デメリットとして考えられるのは、直感にそぐわない式になるため慣れが必要なことでしょうか。

無次元化の方法にはいくつかありますが、今回は私の方法を紹介します。まず、変数、ここでは位置\(x(t)\)[m]と外力\(f(t)\)[N]、そして時間\(t\)[sec]を無次元化します。初めに同じ次元を持つ新しい定数を用意します。この同じ次元を持つ新しい定数の用意の仕方が、一番重要で慣れと経験が必要になります。

マスばねダンパ系の場合、位置に関する新しい定数（基準量）として\(x_{st}\)[m]を用意します。そして無次元化された位置\(\tilde{x}(t)=\frac{x(t)}{x_{st}}\)を導入すると、運動方程式は下記のように変わります。

\begin{align} m\cfrac{d^2\tilde{x}}{dt^2}+c\frac{d\tilde{x}}{dt}+k\tilde{x}(t)=\cfrac{f(t)}{x_{st}} \end{align}

次に外力の無次元化を行います。先ほど新しい定数\(x_{st}\)を用意したので、これを用いて外力と同じ次元の定数を作りたいと思います。そこで運動方程式を眺めていると、\(kx_{st}\)が同じ次元になることが見えてきます。その情報を使って、無次元化された外力\(\tilde{f}(t)=\frac{f(t)}{kx_{st}}\)を導入すると、運動方程式は下記のように変わります。

\begin{align} \cfrac{m}{k}\cfrac{d^2\tilde{x}}{dt^2}+\cfrac{c}{k}\frac{d\tilde{x}}{dt}+\tilde{x}(t)=\tilde{f}(t) \end{align}

これで両辺が無次元化されました。最後に時間\(t\)を無次元化します。ここも外力と同様に、時間と同じ次元を持つ定数を用意します。そこで振動工学の知見を持つ人には馴染みの固有角周波数を用います。固有角周波数\(\omega\)[rad/sec]は\(\sqrt{\cfrac{k}{m}}\)で表されるので、無次元化された時間\(\tilde{t}=\omega t\)を導入すると、運動方程式は下記のように変わります。

\begin{align} \cfrac{d^2\tilde{x}}{d\tilde{t}^2}+\cfrac{c}{\sqrt{mk}}\frac{d\tilde{x}}{d\tilde{t}}+\tilde{x}(\tilde{t})=\tilde{f}(\tilde{t}) \end{align}

最後にこれも振動工学で良く用いられる減衰比\(\zeta=\cfrac{c}{2\sqrt{mk}}\)という無次元量を用いることで、運動方程式は下記のように無次元化されます。

\begin{align} \cfrac{d^2\tilde{x}}{d\tilde{t}^2}+2\zeta\frac{d\tilde{x}}{d\tilde{t}}+\tilde{x}(\tilde{t})=\tilde{f}(\tilde{t}) \end{align}

\(m\)、\(k\)、\(c\)と３つあった定数が\(\zeta\)に集約されており、非常に見通しの良い形になりました。また、物体の振動は減衰比\(\zeta\)のみによって特徴づけられるということがわかりました。

このような無次元化は他の様々なところでも用いられており、例えば流体力学の代表的なナビエストークス方程式を無次元化すると、レイノルズ数などが現れます。特に模型実験などを行う流体力学の分野では無次元化は頻繁に使うので、必須知識といえるでしょう。

解析解(Analytical solution)

次に、この式の解析解（厳密解ともいう）を示します。解法は省略しますが、この運動方程式は2階線形非同次（非斉次）微分方程式ですので、同次微分方程式の一般解（余関数）は\(\zeta\neq1\)の時に

\begin{align} \tilde{x}(\tilde{t})=C_1e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_2e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}} \end{align}

となります。とりあえず今回は重根となる状況は無視しておきます。

特殊解も含めた一般解は定数変化法を用いることで下記のように得られます。

 \begin{align} \tilde{x}(\tilde{t})=\cfrac{1}{2\sqrt{\zeta^2-1}}\left(\int e^{(\zeta+\sqrt{\zeta^2-1})\tilde{t}}\tilde{f}(\tilde{t})d\tilde{t}e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}-\int e^{(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\tilde{f}(\tilde{t})d\tilde{t}e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\right) \end{align}

なお、\(\tilde{f}(\tilde{t})=\tilde{f}_0\)のように定数の場合、一般解は下記のようになる。\(C_1, C_2\)は初期位置と初期速度に関する定数である。

 \begin{align}\tilde{x}(\tilde{t})=\tilde{f}_0+C_1e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_2e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\end{align}

無次元化された初期位置を\(\tilde{x}_0=x_0/x_{st}\)、初期速度を\(\tilde{v}_0=v_0/(\omega x_{st})\)とすると、\(C_1, C_2\)は次のように表すことができます。

 \begin{align}C_1&=-\cfrac{(\zeta-\sqrt{\zeta^2-1})(\tilde{x}_0-\tilde{f}_0)+\tilde{v}_0}{2\sqrt{\zeta^2-1}}\\C_2&=\cfrac{(\zeta+\sqrt{\zeta^2-1})(\tilde{x}_0-\tilde{f}_0)+\tilde{v}_0}{2\sqrt{\zeta^2-1}}\end{align}

 また、\(\tilde{f}(\tilde{t})=\tilde{f}_0\sin(\alpha \tilde{t})\)のような正弦波の場合、一般解は下記のようになる。

 \begin{align}\tilde{x}(\tilde{t})&=\cfrac{(1-\alpha^2)\tilde{f}_0}{(1-\alpha^2)^2+4\zeta^2\alpha^2}\sin(\alpha\tilde{t})-\cfrac{2\zeta\alpha\tilde{f}_0}{(1-\alpha^2)^2+4\zeta^2\alpha^2}\cos(\alpha\tilde{t})+C_1e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_2e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\\&=\cfrac{\tilde{f}_0}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\sin\left(\alpha\tilde{t}+\cos^{-1}\left(\cfrac{1-\alpha^2}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\right)\right)+C_1e^{-(\zeta+\sqrt{\zeta^2-1})\tilde{t}}+C_2e^{-(\zeta-\sqrt{\zeta^2-1})\tilde{t}}\end{align}

なお、\(\alpha\)は固有角周波数で無次元化された周波数比である。

このように、無次元化された運動方程式の解析解が得られました。無次元化したことによって、パラメータの数が減っているため計算ミスが少なくなりそうな気がしませんか？もし有次元の解析解を得たい場合は、無次元化された解析解に元の変数を代入してあげればよいです。

\begin{align}x(t)=\cfrac{1}{k}\cfrac{f_0}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\sin\left(\alpha\omega t+\cos^{-1}\left(\cfrac{1-\alpha^2}{\sqrt{(1-\alpha^2)^2+4\zeta^2\alpha^2}}\right)\right)+C_1x_0e^{-(\zeta+\sqrt{\zeta^2-1})\omega t}+C_2x_0e^{-(\zeta-\sqrt{\zeta^2-1})\omega t}\end{align}

長くなったので、今回はここまで

メカトロニクスって何だろう

メカトロニクス(Mechatronics)とは

メカトロニクスと聞いて、どんなことを思い浮かべるでしょうか？

• メカニクス(Mechanics)とエレクトロニクス(Electronics)の合成語
• 機械工学、電気工学、電子工学、情報工学の技術を融合させた分野
• 機械を電気で動かすこと
• 機械と電気とソフトを組み合わせた製品
• 工場自動化
• 機電一体

 様々な解釈や認識がありますが、メカトロニクスは主に下記の4要素で構成されています。

1. 制御したい対象の情報（位置や速度、画像情報、温度、圧力等々）を
2. センサで計測し
3. 計測した情報をある目的に応じて処理して
4. アクチュエータで制御対象を思い通りに操作する

以下に簡単な例を挙げてみます。

エアコン（Air Conditioner）

 今やどの家庭にも1台はあるような必需品ですが、メカトロニクスの技術が使われています。

上記の4要素に照らし合わせてみると、

1. 室内の温度を
2. 温度センサで計測し
3. 計測温度を設定温度と比較して、その偏差を減らすように
4. 室外機にある圧縮機を動かすことで冷気を作り、ファンを回して風速や風向を制御する

他にも多機能なエアコンでは、赤外線センサやカメラによって家具の配置等の部屋情報を取得して、最適な風速・風向を作り出すなどの機能もメカトロニクスの一種ですね。

自動車（Vehicle）

 日本が誇る一大産業の自動車もメカトロニクス技術が盛り込まれています。特にモータで動かす電気自動車と、モータとエンジンの両方を持つハイブリッド車は、メカトロニクス技術の結晶といってよいでしょう。逆にエンジンのみをもつ自動車は機械技術の結晶とも言えます。

エンジン車にも用いられている自動車のメカトロニクスの一つに、電動パワーステアリングというものがあります。

電動パワーステアリングは、進路を変更するためにハンドルを回す際に、腕の力だけではタイヤを回せないため、モータの力で回してあげよう、といった装置です。これを4要素にまとめてみますと、

1. ハンドルを回す力を
2. トルクセンサで計測し
3. その力を何倍にも増幅するように
4. モータで回す力をアシストする

といった感じでしょう。

もちろん、近年話題の自動運転機能もメカトロニクスの一つです。あまりにも盛り込み過ぎて紹介しきれません…

HDD（Hard Disk Drive）

 最期に、ここまで比較的大きなものを見てきたが、小さいものとしてPCの記憶装置であるハードディスクドライブ、HDDを紹介します。

HDDは、磁性体が塗布された回転する円盤（プラッタ）に対して、磁気ヘッドを近づけることで情報を読み取っている。

磁気ヘッドをプラッタの所定の位置に高速に動かす必要があるが、その際に振動が発生してしまい、ヘッドが所定の位置に留まるまでに時間がかかる場合がある。

ここでメカトロニクス的な技術が使われている。これを4要素にまとめてみると、

1. 磁気ヘッドを現在位置と速度を
2. 位置センサで計測（もしくは推定）し
3. 振動を効率的に抑えるように
4. 磁気ヘッドを動かすモータへの電流を制御する

このように、様々な製品でメカトロニクスな技術が用いられています。電気で動く製品すべてにメカトロニクス的な技術が使われているわけではないですが、探してみると面白いでしょう。