モード分解とモード形式とは
前準備
次式のような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}
まとめ
基本的にモード分解を必要とするモード解析の分野では、減衰が付加されていない無減衰系を対象とする場合が多いため、前半で示した比例減衰系を仮定したモード分解で充分である。
しかし、減衰が付加されているような系では単純にモード分解が行えないため、後半で示したモード形式のように分解することで、厳密なモード減衰比を求めることができるため有用である。
