メカトロニックなカメ

メカトロニクス技術者になりたいカメです

モータの仕様書を眺めてみる

背景

実際にモータを購入した使用する場合に仕様書を必ず眺めるでしょう。しかしモータメーカー各社によって仕様書の記載方法が違うため、見方がわからないことがまあまああるでしょう。特に制御を行う上で重要なのがトルク定数と逆起電力定数なのだが、それについての記載も不十分なことが多いため着目してみよう。 

モータの基礎方程式

まずはモータの基礎的な回路方程式と運動方程式を次式に示す。

\begin{align}v&=Ri+L\cfrac{di}{dt}+K_e\omega\\K_ti&=J\cfrac{d\omega}{dt}+C\omega \end{align}

\(v\)は電圧[V]、\(i\)は電流[A]、\(\omega\)はモータの回転角速度[rad/s]、\(R,L\)はモータの電気抵抗[\(\Omega\)]と自己インダクタンス[H]、\(J,C\)はモータの慣性モーメント[kgm\(^2\)]と粘性抵抗[Nms/rad]、\(K_e,K_t\)はモータの逆起電力定数[Vs/rad]とトルク定数[Nm/A]を表す。

注意してもらいたいのは、上記で示した単位でしかこの方程式が満たされないことである。また、この単位であればモータの逆起電力定数とトルク定数は一致するため、その点も確認していこう。

実際のモータ仕様書

maxon motor

実際のモータ仕様書を見ながらモータの基礎方程式と見比べてみましょう。とりあえずここではmaxon motorのDCモータを見てみます。

https://www.maxongroup.co.jp/medias/sys_master/root/8846534246430/20-JP-99.pdf

maxonのモータの仕様書にはトルク定数が示されていますが、単位が[mNm/A]となっていることに注意しましょう。しかし逆起電力定数はありません。それっぽいものとして回転数定数[rpm/V]というものがあります。この定数を逆起電力定数と同じ次元に置き換えるために、次式のような計算をします。

\begin{align}\cfrac{60}{2\pi\times回転数定数} \end{align}

 実際に回転数定数1330[rpm/V]に対し上式を適用してみますと、0.0072[Vs/rad]となります。トルク定数が7.19[mNm/A]つまり0.00719[Nm/A]であるため、上式で示した計算結果は逆起電力定数を表してることになります。

トルク定数

次に、他のデータからトルク定数を計算してみましょう。トルク定数を計算するのに用いるものとして、定格もしくは最大連続トルクと電流を用います。maxonのデータでは最大連続電流0.577[A]を与えたときの最大連続トルクが4.06[mNm]であるため、この比4.06*1e-3/0.577=0.00704[Nm/A]となり、トルク定数とほぼ一致します。また停動トルクと起動電流の関係からも同様に計算できます(10.5*1e-3/1.46=0.00719[Nm/A])。この若干のずれは計測条件の違いであると予想されます。

回転数定数

次に回転数定数1330[rpm/V]をほかのデータから計算してみましょう。無負荷回転時というのは逆起電力と印加電圧が平衡状態になることであるため、ここから回転数定数が計算できます。つまり無負荷回転数7890[rpm]/公称電圧6[V]=1315[rpm/V]と計算され、ほぼ一致します。しかし実際には摩擦などで逆起電力と印加電圧が完全に平衡状態にならず、少し電流が流れます。そのためモータに加えられる電圧は公称電圧から抵抗による電圧降下を引いた電圧であるため、無負荷電流14.7[mA]と端子間抵抗4.10[\(\Omega\)]を用いると、無負荷回転数7890[rpm]/(公称電圧6[V]-端子間抵抗4.10[\(\Omega\)]×無負荷電流14.7[mA]*1e-3)=1328[rpm/V]となり、回転数定数1330[rpm/V]とより近い値になります。

効率

今度は最大効率81[%]について見てみましょう。効率を計算するためには出力[W}を計算する必要があります。入力する電力[W]は公称電圧[V]×最大連続電流[A]を用いて6[V]×0.577[A]=3.462[W]と計算されます。次に出力パワーは最大連続トルク[Nm]×最大連続トルク時の回転数[rad/s]を用いて計算できます。しかしここでは単位が異なるため、そこに注意して計算すると、4.06*1e-3[Nm]×4830/60*2π[rad/s]=2.054[W]となります。よって最大トルク時の効率は2.054/3.462=59%となります。また、抵抗による熱損失も概算することができ、抵抗4.1[\(\Omega\)]×(0.577[A])\(^2\)=1.37[W]であり、このモータの定格回転時の損失のほとんどが抵抗による熱損失であることも確認できる。

最大効率となる回転数を求めるために、まずは最大効率となる電流を計算してみましょう。このモータの損失のほとんどが抵抗による熱損失であることから、入力電力に対し19%が抵抗による損失であることが推定されます。つまり4.1[\(\Omega\)]×(I[A])\(^2\)=0.19×6[V]×I[A]を満たす電流は0.278[A]となる。この時の発生トルクは、0.278[A]×7.19e-3[Nm/A]=0.0020[Nm]となる。効率が81%であるため、0.0020[Nm]×\(\omega\)/60*2π[rad/s]=0.81×6[V]×0.278[A]の等式から最大効率時の回転数は6451[rpm]となる。

粘性係数

機械的時定数も記載されていますので、ここにも着目してみます。運動方程式から機械的時定数\(t_m\)[s]は慣性モーメント\(J\)[kgm\(^2\)]と粘性係数\(C\)[Nms/rad]の比\(t_m=J/C\)から求められる。慣性モーメントと機械的時定数から粘性係数が概算できます。単位に注意して計算すると1.12*1e-3*1e-4[kgm\(^2\)]/8.87*1e-3[s]=1.263e-05[Nms/rad]のようになります。

澤村電気工業

次に澤村電気工業のDCモータを見てみます。

http://www.sawamura.co.jp/pdf/ss32g.pdf

澤村電気工業では逆起電力定数3.6[V/krpm]とトルク定数0.034[Nm/A]が設定されており、制御設計には親切な仕様書となっています。しかし逆起電力定数の単位が[V/krpm]となっているため、この単位[V/krpm]を基礎方程式の単位[Vs/rad]に合わせるために、60/(2π)×1000をかけると3.6[V/krpm]×60/(2π)/1000=0.0344[Vs/rad]となり、トルク定数[Nm/A]と一致していることが確認できました。

トルク定数

定格電流が2.3[A]で定格トルクが0.054[Nm]であるため、トルク定数を計算すると0.054[Nm]/2.3[A]=0.0235[Nm/A]となります。このように仕様書記載のトルク定数0.034[Nm/A]と乖離があります。

逆起電力定数

定格電圧12[V]で無負荷回転速度3300[rpm]であり、この時の電流が0.4[A]で電機子抵抗1.5[\(\Omega\)]であるため、逆起電力定数[V/krpm]は(12[V]-1.5[\(\Omega\)]×0.4[A])/3.3[krpm]=3.455[V/krpm]であるため、仕様書記載の逆起電力定数3.6[V/krpm]と若干の乖離があります。

無負荷回転速度以外にも定格回転速度を用いても同様に計算できます。計算すると(12[V]-1.5[\(\Omega\)]×2.3[A])/2.5[krpm]=3.42[V/krpm]である。

 

効率

効率の記載はないですが、定格時の効率を計算してみましょう。入力電力は12[V]×2.3[A]=27.6[W]であり、出力パワーは0.054[Nm]*2500/60*2*π[rad/s]=19.79[W]であるため、効率は19.79/27.6=71.7%である。この時の抵抗による熱損失も計算することができ、1.5[\(\Omega\)]×(2.3[A])\(^2\)=7.94[W]であり、このモータの定格回転時の損失のほとんどが抵抗による熱損失であることも確認できる。

粘性係数

慣性モーメントと機械的時定数から粘性係数が概算できます。単位に注意して計算すると0.11*1e-4[kgm\(^2\)]/12*1e-3[s]=9.17e-04[Nms/rad]のようになります。

 

マブチモータ

最後にマブチモータのDCモータを見てみます。

https://www.mabuchi-motor.co.jp/motorize/branch/motor/pdf/fa_130ra.pdf

この仕様書にはトルク定数も逆起電力定数も記載されていない。無負荷時(NO LOAD)や最大効率時(AT MAXIMUM EFFICIENCY)、停動時(STALL)のデータから計算してみましょう。

トルク定数

最大効率時電流が0.66[A]でトルクが0.59[mNm]であるため、トルク定数を計算すると0.59*1e-3[Nm]/0.66[A]=0.000894[Nm/A]となります。また停動時の電流が2.2[A]でトルクが2.55[mNm]であるため、トルク定数は2.55*1e-3[Nm]/2.2[A]=0.0012[Nm/A]となります。このように計測条件によって変動することもあるため注意しましょう。

逆起電力定数

定格電圧1.5[V]で無負荷回転速度9100[rpm]であり、この時の電流が0.2[A]であるが、抵抗値の記載はない。ひとまず抵抗値が0とした場合の、逆起電力定数[V/krpm]は1.5[V]/(9100/60*2π)[rad/s]=0.0016[Vs/rad]であるため、トルク定数と近しい値にはなる。一方最大効率時のデータもあるため、無負荷時と最大効率時で回路方程式が2つ得られ、未知定数が抵抗と逆起電力定数と2つであるため、連立方程式を解くことができる。つまり次式のように表せる。

\begin{align}1.5-R*0.2&=K_e(9100/60*2\pi)\\1.5-R*0.66&=K_e(6990/60*2\pi) \end{align}

上記連立方程式を解くと抵抗値は0.6868[\(\Omega\)]、逆起電力定数は0.0014[Vs/rad]と求まる。 

効率

最大効率時の出力パワーは0.43[W]と記載されているが、これはトルク0.59*1e-3[Nm]×回転数6990/60*2π[rad/s]から計算されたものである。この時の入力電力は1.5[V]×0.66[A]=0.99[W]であるため、この時の効率は43.4%である。抵抗による熱損失は0.6868[\(\Omega\)]×(0.66[A])\(^2\)=0.30[W]であることもわかる。

 

まとめ

このように、各社によってモータの仕様書の記載方法が違うため、ややこしく感じるが、冷静にそれぞれのデータからトルク定数や逆起電力定数を割り出すことができるため、購入する際はこの点に注意しよう。それにしても表記のバラバラっぷりにはいつも困らされます…

なぜ伝達関数にs=jωを代入するのか

背景

1自由度振動系の運動方程式は次式のように表される。

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

上式をラプラス変換すると以下のような伝達関数が得られる。

 \begin{align}G(s)=\cfrac{X(s)}{F(s)}=\cfrac{1}{ms^2+cs+k} \end{align}

この運動方程式の周波数特性を求める際に、伝達関数の\(s\)に\(j\omega\)を代入すればよいと知られている。

 \begin{align}G(j\omega)&=\cfrac{1}{m(j\omega)^2+cj\omega+k}\\&=\cfrac{1}{k-m\omega^2+cj\omega} \end{align}

この時のゲイン(振幅)と位相差は単に絶対値と偏角を求めればよい。

 \begin{align}\left|G(j\omega)\right|&=\cfrac{1}{\sqrt{(k-m\omega^2)^2+(c\omega)^2}}\\\angle G(j\omega)&=\tan^{-1}\cfrac{-c\omega}{k-m\omega^2} \end{align}

間違ってはいないのだが、なぜ\(s=j\omega\)とすればよいのか不透明である。この点について考えてみる。 

フーリエ変換との関係性

フーリエ変換ラプラス変換には密接な関係がある。フーリエ変換の公式は次式である。

\begin{align}\mathcal{F}[f(t)]=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt \end{align}

一方ラプラス変換は次式である。

\begin{align}\mathcal{L}[f(t)]=\int_{0}^{\infty}f(t)e^{-s t}dt \end{align}

フーリエ変換ラプラス変換の違いは積分範囲が\([-\infty, \infty]\)、\([0,\infty]\)であることと、被積分関数が\(e^{-j\omega t}\)、\(e^{-st}\)であることである。雑に言ってしまえば、ラプラス変換の\(s\)を\(j\omega\)にしてしまえばフーリエ変換っぽくなるわけである。

そしてフーリエ変換は周波数特性を求める変換であるため、ラプラス変換によって得られた伝達関数の\(s\)を\(j\omega\)に変えてしまえば、周波数特性が得られるというわけである。

しかし、この方法ではフーリエ変換が周波数特性を表すものであることを理解している必要がある。次に別の説明を試みる。

微分方程式ベースで考える

周波数特性を求める方法として微分方程式に対して、\(f=Fe^{j\omega t}\)を代入する手法がある。例えば1自由度振動系に対して代入し、その解が\(x=Xe^{j\omega t}\)とすると、次式が得られる。

 \begin{align} X\left(-m\omega^2+jc\omega+k\right)e^{j\omega t}=Fe^{j\omega t}\end{align}

上式を満たす\(X\)は次式のようになる。

 \begin{align} X=\cfrac{F}{-m\omega^2+jc\omega+k}\end{align}

つまり、外力振幅\(F\)と変位振幅\(X\)の比は次式のように表され、これはまさに伝達関数に\(s=j\omega\)を代入したものである。

上記では一自由度振動系を例にしたが、一般的な伝達関数微分方程式)でも同様のことがいえる。

次のような伝達関数を考えてみる。

 \begin{align} G(s)=\cfrac{Y(s)}{U(s)}=\cfrac{b_{n-1}s^{n-1}+b_{n-2}s^{n-2}+\cdots+b_0}{s^n+a_{n-1}s^{n-1}+\cdots+a_0}\end{align}

上式の伝達関数は次の微分方程式で表される。

 \begin{align} \cfrac{d^n y}{dt^n}+a_{n-1}\cfrac{d^{n-1} y}{dt^{n-1}}+\cdots+a_0y=b_{n-1}\cfrac{d^{n-1} u}{dt^{n-1}}+b_{n-2}\cfrac{d^{n-2} u}{dt^{n-2}}+\cdots+b_0u\end{align}

上式に、\(u=Ue^{j\omega t}\)を代入し、その解が\(y=Ye^{j\omega t}\)とすると、次式が得られる。

 \begin{align} \left((j\omega)^n+a_{n-1}(j\omega)^{n-1}+\cdots+a_0\right)Y=\left(b_{n-1}(j\omega)^{n-1}+b_{n-2}(j\omega)^{n-2}+\cdots+b_0\right)U\end{align}

 \begin{align} \cfrac{Y}{U}=\cfrac{b_{n-1}(j\omega)^{n-1}+b_{n-2}(j\omega)^{n-2}+\cdots+b_0}{(j\omega)^n+a_{n-1}(j\omega)^{n-1}+\cdots+a_0}\end{align}

これで、ラプラス変換して得られた伝達関数に\(s=j\omega\)を代入することで周波数特性が得られることを示された。

ラプラス逆変換で考える

周波数特性を求めるために\(f=\sin(\omega t)\)を代入することを考える。この\(f\)をラプラス変換すると\(F(s)\)は次式のようになることが知られている。

\begin{align}F(s)=\cfrac{\omega}{s^2+\omega^2} \end{align}

次のような伝達関数を考える。

 \begin{align}\cfrac{X(s)}{F(s)}=\cfrac{a}{s+b} \end{align}

この時、\(X(s)\)は次式のように表せる。

 \begin{align}X(s)=\cfrac{a}{s+b}\cfrac{\omega}{s^2+\omega^2} \end{align}

上式を部分分数分解すると次式のようになる。

 \begin{align}X(s)&=\cfrac{A}{s+b}+\cfrac{Bs+C}{s^2+\omega^2}\\&=\cfrac{A(s^2+\omega^2)+(Bs+C)(s+b)}{(s+b)(s^2+\omega^2)}\\&=\cfrac{(A+B)s^2+(Bb+C)s+A\omega^2+Cb}{(s+b)(s^2+\omega^2)}\end{align}

係数を比較すると\(A,B,C\)は以下のように得られる。

 \begin{align}A&=\cfrac{a\omega}{\omega^2+b^2}\\B&=\cfrac{-a\omega}{\omega^2+b^2}\\C&=\cfrac{ab\omega}{\omega^2+b^2}\end{align}

 つまり部分分数分解した式は次式のようになる。

 \begin{align}X(s)&=\cfrac{a\omega}{\omega^2+b^2}\cfrac{1}{s+b}+\cfrac{a\omega}{\omega^2+b^2}\cfrac{-s+b}{s^2+\omega^2}\\&=\cfrac{a\omega}{\omega^2+b^2}\cfrac{1}{s+b}-\cfrac{a\omega}{\omega^2+b^2}\cfrac{s}{s^2+\omega^2}+\cfrac{ab}{\omega^2+b^2}\cfrac{\omega}{s^2+\omega^2}\end{align}

上式を逆ラプラス変換すると、次式のようになる。

 \begin{align}x(t)&=\cfrac{a\omega}{\omega^2+b^2}e^{-bt}-\cfrac{a\omega}{\omega^2+b^2}\cos(\omega t)+\cfrac{ab}{\omega^2+b^2}\sin(\omega t)\\&=\cfrac{a\omega}{\omega^2+b^2}e^{-bt}+\cfrac{a}{\sqrt{\omega^2+b^2}}\sin\left(\omega t+\tan^{-1}\left(-\frac{\omega}{b}\right)\right)\end{align}

上式の右辺第一項は過渡応答であるため、周波数特性は右辺第二項になります。

一方、伝達関数に\(s=j\omega\)を代入すると次式のようになる。

 \begin{align}G(j\omega)&=\cfrac{a}{j\omega+b}\\\left|G(j\omega)\right|&=\cfrac{a}{\sqrt{\omega^2+b^2}}\\\angle G(j\omega)&=\tan^{-1}\left(-\cfrac{\omega}{b}\right) \end{align}

上記のように一致することが確認できた。

 

上記では1次の場合で証明したが、2次以上でも複雑になるが同様に証明できる(はず)。

 

 

モータブレーキの原理

前準備

1自由度回転振動系の運動方程式は次式のように表される。

\begin{align}J\cfrac{d^2\theta}{dt^2}+C\cfrac{d\theta}{dt}+K\theta=N \end{align}

 

モータブレーキとは

モータブレーキは永久磁石式モータが装着された運動系に対し、モータの端子間を短絡することでブレーキ力を発生させることである。この原理の数式は次式のようになる。

\begin{align}L\cfrac{di}{dt}+Ri=e_v \end{align}

\(e_v\)[V]は運動系の速度に起因される起電力で、速度起電力(motional electromotive force)や逆起電力(back electromotive force)と呼ばれる。速度起電力は基本的に運動の速度に比例しており、次式のように表記する。\(\phi_e\)は速度起電力定数や逆起電力定数と呼ばれる。

\begin{align}e_v = \phi_e \cfrac{d\theta}{dt}\end{align}

一方、モータに電流が流れるとローレンツ力\(f_e\)が発生し、運動系に力を及ぼす。このことを数式で表すと次式のようになる。 \(\phi_t\)はトルク定数と呼ばれる。

\begin{align}N_e = -\phi_t i\end{align}

なお、速度起電力定数\(\phi_e\)[V/(rad/s)]とトルク定数\(\phi_t\)[Nm/A]は原理的に同じ値になります。ただし単位が異なる場合があるため注意が必要である。具体的には速度起電力定数の単位の種類として[V/(deg/s)]や[V/Hz]、[Vrms/rpm]、トルク定数として[Nm/Arms]など様々なパターンがあるため、使用する際は注意しよう。

外力モーメント以外にローレンツ力が加わった場合の運動方程式は次式のようになる。

\begin{align}J\cfrac{d^2\theta}{dt^2}+C\cfrac{d\theta}{dt}+K\theta=N+N_e \end{align}

最後に、以上の4つの方程式をまとめると次式のようになる。

\begin{align}J\cfrac{d^2\theta}{dt^2}+C\cfrac{d\theta}{dt}+K\theta+\phi_ti=N\\L\cfrac{di}{dt}+Ri=\phi_e \cfrac{d\theta}{dt} \end{align}

上式を1つにまとめるために、状態方程式を用いると次式のようになる。

\begin{align}\left[\begin{array}{ccc}1&0&0\\0&J&0\\0&0&L\end{array}\right]\cfrac{d}{dt}\left(\begin{array}{c}\theta\\\cfrac{d\theta}{dt}\\i\end{array}\right)=\left[\begin{array}{ccc}0&1&0\\-K&-C&-\phi_t\\0&\phi_e&-R\end{array}\right]\left(\begin{array}{c}\theta\\\cfrac{d\theta}{dt}\\i\end{array}\right)+\left(\begin{array}{c}0\\N\\0\end{array}\right) \end{align}

別の方法としてラプラス変換を用いて、1つにまとめると次式のようになる。\(\tilde{\theta}\)や\(\tilde{N}\)はラプラス変換された変数である。

\begin{align}\left(Js^2+Cs+K+\cfrac{\phi_e\phi_ts}{Ls+R}\right)\tilde{\theta}=\tilde{N} \end{align}

性能の良いモータはインダクタンス\(L\)が抵抗\(R\)に比べて十分に小さい、つまり電気系時定数\(L/R\)が小さい。そのためラプラス変換された式は次式のように近似することができる。

\begin{align}\left(Js^2+\left(C+\cfrac{\phi_e\phi_t}{R}\right)s+K\right)\tilde{\theta}=\tilde{N} \end{align}

上式より、モータによる効果は減衰力を増やしていることを意味しており、抵抗値が小さいほど、速度起電力定数やトルク定数が大きいほどその効果は大きい。これがモータブレーキと呼ばれる所以である。

アナロジー視点

以前に(1自由度振動系のお話 その4 - メカトロニックなカメ)、電気系は機械系に変換(アナロジー)することができると述べた。ここでもアナロジーを行ってみよう。

上述したようにモータブレーキでは抵抗の逆数は減衰に相当する。以前の記事では逆数ではない点に注意されたい。

\begin{align}L\cfrac{di}{dt}+\cfrac{1}{c'}i=\phi_e \cfrac{d\theta}{dt} \end{align}

上式を整理すると次のようになる。

\begin{align}c'\cfrac{di}{dt}+\cfrac{1}{L}i=\cfrac{c'\phi_e}{L} \cfrac{d\theta}{dt} \end{align}

つまり電流の1階微分\(\cfrac{di}{dt}\)は角速度に相当し、電流は角度に相当する。このことからインダクタンスの逆数は剛性係数に相当することもわかる。

\begin{align}c'\cfrac{d\theta'}{dt}+k'\theta'=k'c'\phi_e \cfrac{d\theta}{dt} \end{align}

 上式から主系の角速度\( \cfrac{d\theta}{dt}\)が0の時は、アナロジーされた減衰力や剛性力は発生しない。つまり減衰とばねが直列に繋がれたモデルであると推察できる。実際に減衰とばねが直列に繋がれたモデルは力のつり合いの関係から次式のように表せる。

\begin{align}f=k\theta_0=c\cfrac{d(\theta-\theta_0)}{dt} \\c\cfrac{d\theta_0}{dt}+k\theta_0=c\cfrac{d\theta}{dt} \\\cfrac{d\theta_0}{dt}+\cfrac{k}{c}\theta_0=\cfrac{d\theta}{dt}\end{align}

再度アナロジーされた式に戻ると、以下のように変形される。 

\begin{align}c'\cfrac{d\theta'}{dt}+k'\theta'=k'c'\phi_e \cfrac{d\theta}{dt}\\\cfrac{1}{k'\phi_e}\cfrac{d\theta'}{dt}+\cfrac{1}{c'\phi_e}\theta'= \cfrac{d\theta}{dt} \end{align}

つまり、\(\theta_0=\cfrac{\theta'}{k'\phi_e}, \cfrac{k}{c}=\cfrac{k'}{c'}\)の関係があり、推察通り減衰とばねが直列に繋がれたモデルが主系に繋がれていることがわかった。

 

ヒステリシス減衰におけるモード分解と周波数特性

前準備

\(\boldsymbol{M}\)を質量行列、\(\boldsymbol{C}\)を減衰行列、\(\boldsymbol{K}\)を剛性行列、\(\boldsymbol{x}\)を位置ベクトル、\(\boldsymbol{f}\)を外力ベクトルとしたとき多自由度振動系の運動方程式は次式のように表される。

\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}

このとき、剛性行列に相当する部分がゴムのような高分子化合物である場合、ヒステリシス減衰(複素剛性、構造減衰とも)が発生する。ヒステリシス減衰を考慮した運動方程式は次式のようになる。

\begin{align} \boldsymbol{M}\cfrac{d^2\boldsymbol{x}}{dt^2}+\boldsymbol{C}\cfrac{d\boldsymbol{x}}{dt}+\boldsymbol{K}\left(\boldsymbol{I}+i\boldsymbol{G}\right)\boldsymbol{x}=\boldsymbol{f} \end{align}

\(\boldsymbol{G}\)は損失係数で構成された行列、\(i\)は虚数単位である。

モード分解

ヒステリシス減衰がない場合と同様の手順でモード分解を行う。まずは次式の一般固有値問題の解である変換行列\(\boldsymbol{T}\)を求める。

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

 変換行列\(\boldsymbol{T}\)と変位ベクトルを\(\boldsymbol{x}=\boldsymbol{T}\boldsymbol{\phi}\)のように変換すると、以下のようにモード分解された運動方程式が得られる。

\begin{align} \boldsymbol{T}^\textrm{T}\boldsymbol{M}\boldsymbol{T}\cfrac{d^2\boldsymbol{\phi}}{dt^2}+\boldsymbol{T}^\textrm{T}\boldsymbol{C}\boldsymbol{T}\cfrac{d\boldsymbol{\phi}}{dt}+\boldsymbol{T}^\textrm{T}\boldsymbol{K}\left(\boldsymbol{I}+i\boldsymbol{G}\right)\boldsymbol{T}\boldsymbol{\phi}=\boldsymbol{T}^\textrm{T}\boldsymbol{f}\\\boldsymbol{\bar{M}}\cfrac{d^2\boldsymbol{\phi}}{dt^2}+\boldsymbol{\bar{C}}\cfrac{d\boldsymbol{\phi}}{dt}+\left(\boldsymbol{\bar{K}}+i\boldsymbol{T}^\textrm{T}\boldsymbol{G}\boldsymbol{T}\right)\boldsymbol{\phi}=\boldsymbol{T}^\textrm{T}\boldsymbol{f}\end{align}

ここで損失係数行列である\(\boldsymbol{G}\)は対角化できるか、2自由度振動系のを例に考えてみる。2自由度振動系の剛性行列は次のように表せる。

\begin{align}\boldsymbol{K}=\left[\begin{array}{cc}k_1+k_{12}&-k_{12}\\-k_{12}&k_{12}+k_2\end{array}\right]\end{align}

それぞれの剛性係数に損失係数が含まれていると考えると、ヒステリシス減衰を含んだ剛性行列は次のようになる。

\begin{align}\boldsymbol{K}\left(\boldsymbol{I}+i\boldsymbol{G}\right)&=\left[\begin{array}{cc}k_1(1+ig_1)+k_{12}(1+ig_{12})&-k_{12}(1+ig_{12})\\-k_{12}(1+ig_{12})&k_{12}(1+ig_{12})+k_2(1+ig_{2})\end{array}\right]\\&=\left[\begin{array}{cc}k_1+k_{12}&-k_{12}\\-k_{12}&k_{12}+k_2\end{array}\right]+i\left[\begin{array}{cc}g_1+g_{12}&-g_{12}\\-g_{12}&g_{12}+g_2\end{array}\right]\end{align}

上式より損失係数行列と剛性行列は必ずしも比例関係を満たさないことがわかる。よって減衰行列と同様、損失係数行列が質量行列や剛性行列と比例関係にあるという仮定を満たせば、モード分解が可能である。

周波数特性

上記で表したヒステリシス減衰は周波数特性上で近似された式である。実際に周波数特性がどうなるか確認してみよう。簡単のためここでは1自由度振動系の場合を考える。

 \begin{align} m\cfrac{d^2x}{dt^2}+c\cfrac{dx}{dt}+k\left(1+ig\right)x=f \end{align}

周波数特性を考える際に上記の微分方程式を解く方法もあるが、別の方法を用いる。まず外力に\(f=Fe^{i\omega t}\)を代入する。これはオイラーの等式を用いると次式のようになる。 

\begin{align} f=Fe^{i\omega t}=F\cos \omega t+Fi\sin \omega t\end{align}

上式の虚部は実際には発生しえない力である。しかし周波数特性を考える際に、位相差に相当する虚部を抽出するために必要な操作である。上式を代入し、その時の解が\(x=Xe^{i\omega t}\)とすると、次式が得られる。

 \begin{align} X\left(-m\omega^2+ic\omega+k\left(1+ig\right)\right)e^{i\omega t}=Fe^{i\omega t}\end{align}

上式を満たす\(X\)は次式のようになる。

 \begin{align} X=\cfrac{F}{-m\omega^2+ic\omega+k\left(1+ig\right)}\end{align}

上式は複素数を含んでおり、やや扱いが難しい。そこで次のような手順で変形を加えていく。ここで式の見通しをよくするため、一旦分母は\(a+ib\)とする。

 \begin{align} X&=\cfrac{F}{a+ib}\\&=\cfrac{F}{a^2+b^2}(a-ib)\\&=\cfrac{F}{a^2+b^2}\sqrt{a^2+b^2}\left(\cfrac{a}{\sqrt{a^2+b^2}}-i\cfrac{b}{\sqrt{a^2+b^2}}\right)\\&=\cfrac{F}{\sqrt{a^2+b^2}}\left(\cos \Theta-i\sin \Theta\right)\\&=\cfrac{F}{\sqrt{a^2+b^2}}e^{-i\Theta}\end{align}

元々の解が\(x=Xe^{i\omega t}\)であるため、解\(x\)は次式のようになる。

 \begin{align} x&=\cfrac{F}{\sqrt{a^2+b^2}}e^{-i\Theta}e^{i\omega t}\\&=\cfrac{F}{\sqrt{a^2+b^2}}e^{i\left(\omega t-\Theta\right)}\end{align}

このように外力\(f=Fe^{i\omega t}\)を加えたときの解は、上式のように大きさが\(\cfrac{1}{\sqrt{a^2+b^2}}\)倍され、位相が\(\Theta=\cos^{-1}\cfrac{a}{\sqrt{a^2+b^2}}\)だけ遅れる応答が得られる。つまり無次元化すると振幅倍率は次のようになる。

 \begin{align} \cfrac{kX}{F}&=\cfrac{k}{\sqrt{\left(k-m\omega^2\right)^2+\left(c\omega+kg\right)^2}}\\&=\cfrac{1}{\sqrt{\left(1-\alpha^2\right)^2+\left(2\zeta\alpha+g\right)^2}}\end{align}

\(\alpha=\omega/\omega_0, \omega_0=\sqrt{k/m}, \zeta=c/2\sqrt{mk}\)である。ヒステリシス減衰がない場合と違い、減衰項に\(g\)が無次元化周波数\(\alpha\)に依存せずに加わっていることに大きな違いがある。これを図に表すと図1のようになる。

f:id:meckam:20210208234844j:plain

損失係数を変化させたときの応答

図1からわかるように損失係数が大きくなると、振幅倍率はピークが小さくなるだけでなく低周波も小さくなっていくことが特徴である。一方減衰係数が大きくなる場合ではピークが小さくなるのみで低周波は変化しない。ここがヒステリシス減衰を含む振動系の特徴である。

いろいろな減衰

前準備

\(\boldsymbol{M}\)を質量行列、\(\boldsymbol{C}\)を減衰行列、\(\boldsymbol{K}\)を剛性行列、\(\boldsymbol{x}\)を位置ベクトル、\(\boldsymbol{f}\)を外力ベクトルとしたとき多自由度振動系の運動方程式は次式のように表される。

\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}

上式のように表される減衰は速度に比例した減衰であり、最も一般的な減衰モデルである。しかし、減衰は運動する物体の形状や減衰力を生み出す流体の粘性などで大きく変化する。それらについてまとめてみよう。

速度比例減衰

前準備で示した減衰は速度に比例した減衰であるため、速度比例減衰と呼べる。この減衰の場合、モデルは線形微分方程式であるため非常に扱いが容易である。

\begin{align} \boldsymbol{f}_c=-\boldsymbol{C}\cfrac{d\boldsymbol{x}}{dt} \end{align}

この減衰モデルで表せる条件は、運動の周りの流体が層流であることである。つまりレイノルズ数が小さいことが条件である。レイノルズ数\(\textrm{Re}\)は代表長さ\(L\)、代表速度\(U\)、動粘性係数\(\nu\)、流体密度\(\rho\)、粘性係数\(\mu\)を用いると次式のように表されます。

\begin{align} \textrm{Re}=\cfrac{LU}{\nu}=\cfrac{LU\rho}{\mu}\end{align}

上式からさらに言い換えると速度が小さくて、動粘性係数が大きい場合にこの減衰が適用可能である。

速度二乗減衰

レイノルズ数\(\textrm{Re}\)が大きくて、物体の周りの流体が乱流となる場合、減衰力は速度の二乗に比例するとされている。 

\begin{align} \boldsymbol{f}_c=-\boldsymbol{C}\left|\cfrac{d\boldsymbol{x}}{dt}\right|\cfrac{d\boldsymbol{x}}{dt} \end{align}

速度が大きい範囲、動粘性係数が小さい場合にはこの減衰力モデルを用いることとなる。振動を扱う際は微小振動に着目することが多いため、この減衰力モデルを用いることは少ない。しかし水中など動粘性係数が小さい場合にはその限りではない。

 

クーロン摩擦減衰

上記で示した速度比例減衰や速度二乗減衰は流体と物体との摩擦を意味しており、摩擦と減衰は言葉は違うが、厳密には同じである。クーロン摩擦は一般的に動摩擦と静摩擦と呼ばれるもので、垂直抗力を一定とすると次式のように表される。

\begin{align} \boldsymbol{f}_c=\left\{\begin{array}{ll}-\boldsymbol{F}_d\textrm{sgn} \left(\cfrac{d\boldsymbol{x}}{dt}\right)&\textrm{if }   \cfrac{d\boldsymbol{x}}{dt}\neq 0\\-\boldsymbol{F}_s&\textrm{if }\cfrac{d\boldsymbol{x}}{dt}=0\end{array}\right. \end{align}

厳密には上記よりも複雑な式となるが、上記で充分なことが多い。

 

レイリー減衰

比例粘性減衰とも呼ばれ、質量行列と剛性行列に比例した減衰である。

\begin{align} \boldsymbol{f}_c=-\left(\alpha\boldsymbol{M}+\beta\boldsymbol{K}\right)\cfrac{d\boldsymbol{x}}{dt} \end{align}

この減衰はモード分解する際に非常に都合の良い減衰力モデルである。解析用の減衰力モデルといえる。

 

ヒステリシス減衰

履歴減衰、構造減衰、複素剛性とも呼ばれ、ゴムなど高分子化合物でよくみられる特性である。このモデルの特徴として、変位に比例した減衰力が働くことであり、ゴムを伸縮させたときに発生する熱を表すモデルである。上記のように一般的な数式で表すことができず、周波数特性の時に次式のような近似式を用いられる。

\begin{align} \boldsymbol{f}_c=-\boldsymbol{K}\left(\boldsymbol{I}+i\boldsymbol{G}\right)\boldsymbol{x} \end{align}

 

剛体モードを持つ振動系へのモード分解

前準備

\(\boldsymbol{M}\)を質量行列、\(\boldsymbol{C}\)を減衰行列、\(\boldsymbol{K}\)を剛性行列、\(\boldsymbol{x}\)を位置ベクトル、\(\boldsymbol{f}\)を外力ベクトルとしたとき多自由度振動系の運動方程式は次式のように表される。

\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{T}\)を求める。

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

 変換行列\(\boldsymbol{T}\)と変位ベクトルを\(\boldsymbol{x}=\boldsymbol{T}\boldsymbol{\phi}\)のように変換すると、以下のようにモード分解された運動方程式が得られる。

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

上記のモード分解では\(\boldsymbol{\bar{M}}^{-1}\boldsymbol{\bar{C}}\)が対角行列になるとは限らない。しかし実用上は減衰が小さいため、\(\boldsymbol{\bar{M}}^{-1}\boldsymbol{\bar{C}}\)の非対角項を無視しても問題ないことが多い。

剛体モード(Rigid body mode)とは

上記で扱ってきたモードというのは、厳密には振動モードと呼ばれるものであり、モードには他にも剛体モードというものがある。剛体モードとは、数学的にはモード分解やモード形式によって得られた固有値(もしくは固有角周波数)が0となるモードを差し、実用上では復元力(ばね力などの元の位置に戻す力)が働かないモードであり、力を加えると動き続ける状態のことを意味している。

例えば、人が歩行する際の移動は剛体モードであり、人がジャンプしたときは地面に引き寄せられるので振動モードである。メカトロニクスの分野で剛体モードが良く用いられるのが、電磁モータの回転である。モータの回転はトルクを与え続けると常に回転し続け、回転角度は無限に大きくなっていく。このように回転する物体を制御する際には剛体モードの考えは不可欠である。

また、電磁モータにトルクを与える際に、一定のトルクを与えることができればモータは滑らかな回転をするが、実際には一定のトルクを加えることは困難である。そのためトルクに脈動が生じ、モータの回転も振動的になる。この振動が故障の問題となることが多い。

自動車を例にするとさらに顕著である。ガソリン等の内燃機関の多くはロータ(回転する物体)が2回転する間に一度だけガソリンを燃焼させてトルクを生み出している。そのためトルク脈動が顕著となるため、自動車特有の揺れを発生させている。

以上のことを理由に、今回は剛体モードが含まれた多自由度振動系を考える。

剛体モードを含んだ運動方程式

剛体モードを含んだ2自由度振動系の運動方程式は次のようになる。

\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_{12}&-k_{12}\\-k_{12}&k_{12}\end{array}\right]\left(\begin{array}{c}x_1\\x_2\end{array}\right)=\left(\begin{array}{c}f\\0\end{array}\right) \end{align}

上式をモード分解してみよう。モード分解の際は減衰行列を無視して、下記の一般化固有値問題を考える。

\begin{align} \left[\begin{array}{cc}m_1&0\\0&m_2\end{array}\right]\boldsymbol{T}\boldsymbol{\Omega}^2=\left[\begin{array}{cc}k_{12}&-k_{12}\\-k_{12}&k_{12}\end{array}\right]\boldsymbol{T}\end{align}

固有値行列と固有ベクトル行列は次のようになる。

\begin{align} \boldsymbol{\Omega}^2&=\left[\begin{array}{cc}0&0\\0&k_{12}\left(\cfrac{1}{m_1}+\cfrac{1}{m_2}\right)\end{array}\right]\\\boldsymbol{T}&=\left[\begin{array}{cc}1&1\\1&-\cfrac{m_1}{m_2}\end{array}\right]\end{align}

よってモード分解を行うと次式のようになる。

\begin{align} \cfrac{d^2\phi_1}{dt^2}+\bar{c}_1\cfrac{d\phi_1}{dt}=\bar{f}_1\\\cfrac{d^2\phi_2}{dt^2}+2\zeta_2\omega_2\cfrac{d\phi_2}{dt}+\omega_2^2=\bar{f}_2 \end{align}

一つのモードでは固有角周波数は持たず、外乱\(\bar{f}_1\)が一定値ならば、モード変位\(\phi_1\)は常に増加(または減少)し続け、モード速度\(\frac{d\phi_1}{dt}\)はモード減衰\(\bar{c}_1\)に応じた値に収束する。

また、モード変位\(\phi_2\)は振動モードであり、ステップ上の外乱を印加したときに振動しながら、ある値に収束する。

まとめ

振動モードのみの場合でも、剛体モードを含む場合でもモード分解を行うことで、モードの分離ができる。複雑な振動系ほどモード分解が有効である。今回はモード分解を行ったが、減衰が大きい場合はモード形式を用いても同様にモード分解が行える。

モード分解の数値例

前準備

\(\boldsymbol{M}\)を質量行列、\(\boldsymbol{C}\)を減衰行列、\(\boldsymbol{K}\)を剛性行列、\(\boldsymbol{x}\)を位置ベクトル、\(\boldsymbol{f}\)を外力ベクトルとしたとき多自由度振動系の運動方程式は次式のように表される。

\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{T}\)を求める。

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

 変換行列\(\boldsymbol{T}\)と変位ベクトルを\(\boldsymbol{x}=\boldsymbol{T}\boldsymbol{\phi}\)のように変換すると、以下のようにモード分解された運動方程式が得られる。

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

上記のモード分解では\(\boldsymbol{\bar{M}}^{-1}\boldsymbol{\bar{C}}\)が対角行列になるとは限らない。しかし実用上は減衰が小さいため、\(\boldsymbol{\bar{M}}^{-1}\boldsymbol{\bar{C}}\)の非対角項を無視しても問題ないことが多い。

減衰が無視できないときはモード形式(Eigenstructure decomposition)を用いるとよい。モード形式を求めるために、まず次式で示すように運動方程式状態方程式に変換する。

\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}

次に、以下に示す一般化固有値問題固有値\(\boldsymbol{\Lambda}\)と固有ベクトル行列\(\boldsymbol{V}\)を求める。

\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}

実数な固有値に対応した固有ベクトルを\(\boldsymbol{v}_i^r\)、複素数固有値に対応した固有ベクトルを\(\boldsymbol{v}_j^c\)としたとき、次の変換行列を考える。

\begin{align}&\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}

上記のように、元の運動方程式をモード減衰比\(\zeta_i\)とモード固有角周波数\(\omega_i\)に分解が可能である。

減衰が小さい場合の数値例

次のような運動方程式を考える。

\begin{align} \left[\begin{array}{cc}3&0\\0&2\end{array}\right]\cfrac{d^2}{dt^2}\left(\begin{array}{c}x_1\\x_2\end{array}\right)&+ \left[\begin{array}{cc}0.1&-0.2\\-0.2&0.5\end{array}\right]\cfrac{d}{dt}\left(\begin{array}{c}x_1\\x_2\end{array}\right)\\&+\left[\begin{array}{cc}30&-10\\-10&10\end{array}\right]\left(\begin{array}{c}x_1\\x_2\end{array}\right)=\left(\begin{array}{c}f\\0\end{array}\right) \end{align}

状態方程式で表記すると次式のようになる。

\begin{align} \left[\begin{array}{cccc}1&0&0&0\\0&1&0&0\\0&0&3&0\\0&0&0&2\end{array}\right]\cfrac{d}{dt}\left(\begin{array}{c}x_1\\x_2\\v_1\\v_2\end{array}\right)= \left[\begin{array}{cccc}0&0&1&0\\0&0&0&1\\-30&10&-0.1&0.2\\10&-10&0.2&-5\end{array}\right]\left(\begin{array}{c}x_1\\x_2\\v_1\\v_2\end{array}\right)+\left(\begin{array}{c}0\\0\\f\\0\end{array}\right) \end{align}

モード分解

次式の一般化固有値問題の解\(\boldsymbol{T}\)を求める。固有値問題を解く際はOctaveを使用した。

\begin{align} \left[\begin{array}{cc}3&0\\0&2\end{array}\right]\boldsymbol{T}\boldsymbol{\Omega}^2=\left[\begin{array}{cc}30&-10\\-10&10\end{array}\right]\boldsymbol{T}\\\boldsymbol{T}=\left[\begin{array}{cc}-0.5037&-0.2822\\0.3456&-0.6169\end{array}\right]\\\boldsymbol{\Omega}^2=\left[\begin{array}{cc}12.2871&0\\0&2.7129\end{array}\right]\end{align}

上記変換行列を用いると、運動方程式は次式のようになる。なお、質量行列が単位行列となっており、この変換行列は正規化の特性も含まれている。

\begin{align} \cfrac{d^2\boldsymbol{\phi}}{dt^2}&+ \left[\begin{array}{cc}0.1547&-0.135\\-0.135&0.1286\end{array}\right]\cfrac{d\boldsymbol{\phi}}{dt}\\&+\left[\begin{array}{cc}12.29&0\\0&2.713\end{array}\right]\boldsymbol{\phi}=\left(\begin{array}{c}-0.5037\\-0.2822\end{array}\right)f\\\left(\begin{array}{c}x_1\\x_2\end{array}\right)&=\left[\begin{array}{cc}-0.5037&-0.2822\\0.3456&-0.6169\end{array}\right] \boldsymbol{\phi}\end{align}

 完全にモードごとに分解するために、減衰項の非対角項(連成項)を0とする必要がある。連成項を0とすることで次式のように完全に分離が可能となる。

\begin{align} \cfrac{d^2\phi_1}{dt^2}+0.1547\cfrac{d\phi_1}{dt}+12.29\phi_1=-0.5037f\\\cfrac{d^2\phi_2}{dt^2}+0.1286\cfrac{d\phi_2}{dt}+2.713\phi_1=-0.2822f\\\left(\begin{array}{c}x_1\\x_2\end{array}\right)=\left[\begin{array}{cc}-0.5037&-0.2822\\0.3456&-0.6169\end{array}\right] \left(\begin{array}{c}\phi_1\\\phi_2\end{array}\right)\end{align}

この系のモード減衰比とモード固有角周波数は以下のようになる。

\begin{align} \zeta_1&=0.022069\\\zeta_2&=0.039042\\\omega_1&=3.5053\\\omega_2&=1.6471\end{align}

モード形式

次式の一般化固有値問題の解\(\boldsymbol{V}\)を求める。なお上付きバーは複素共役を表す。

\begin{align} \left[\begin{array}{cccc}1&0&0&0\\0&1&0&0\\0&0&3&0\\0&0&0&2\end{array}\right]\boldsymbol{V}\boldsymbol{\Lambda}&= \left[\begin{array}{cccc}0&0&1&0\\0&0&0&1\\-30&10&-0.1&0.2\\10&-10&0.2&-5\end{array}\right]\boldsymbol{V}\\\boldsymbol{V}&= \left[\begin{array}{cccc}\boldsymbol{v}_1&\boldsymbol{\bar{v}}_1&\boldsymbol{v}_2&\boldsymbol{\bar{v}}_2\end{array}\right]\\\boldsymbol{v}_1=\left(\begin{array}{c}-0.010955+0.28062i\\-0.014844-0.19252i\\0.98332+0.016677i\\-0.67287+0.066842i\end{array}\right)&, \boldsymbol{v}_2=\left(\begin{array}{c}0.0064523+0.27488i\\-0.018709-0.60148i\\0.45242-0.028338i\\-0.99207-0.0079267i\end{array}\right)\end{align}

次の変換行列\(\boldsymbol{T}\)を求める。

\begin{align} \boldsymbol{T}&= \left[\begin{array}{cccc}\textrm{Re}(\boldsymbol{v}_1)&\textrm{Im}(\boldsymbol{v}_1)&\textrm{Re}(\boldsymbol{v}_2)&\textrm{Im}(\boldsymbol{v}_2)\end{array}\right]\\&=\left[\begin{array}{cccc}-0.010955&0.28062&0.0064523&0.27488\\-0.014844&-0.19252&-0.018709&-0.60148\\0.98332&0.016677&0.45242&-0.028338\\-0.67287&0.066842&-0.99207&-0.0079267\end{array}\right]\end{align}

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

\begin{align} \cfrac{d}{dt}\left(\begin{array}{c}\phi_1\\\phi_2\\\psi_1\\\psi_2\end{array}\right)=&\left[\begin{array}{cccc}-0.0772&-3.5011&0&0\\3.5011&-0.0772&0&0\\0&0&-0.0644&-1.6474\\0&0&1.6474&-0.0644\end{array}\right]\left(\begin{array}{c}\phi_1\\\phi_2\\\psi_1\\\psi_2\end{array}\right)\\&+\left(\begin{array}{c}0.2585\\-0.004226\\0.1757\\0.01049\end{array}\right)f\\\left(\begin{array}{c}x_1\\x_2\\v_1\\v_2\end{array}\right)=&\left[\begin{array}{cccc}-0.010955&0.28062&0.0064523&0.27488\\-0.014844&-0.19252&-0.018709&-0.60148\\0.98332&0.016677&0.45242&-0.028338\\-0.67287&0.066842&-0.99207&-0.0079267\end{array}\right]\left(\begin{array}{c}\phi_1\\\phi_2\\\psi_1\\\psi_2\end{array}\right)\end{align}

 モード別に分けることで次式が得られる。

\begin{align} &\cfrac{d^2\phi_1}{dt^2}+0.1545\cfrac{d\phi_1}{dt}+12.26\phi_1=0.2585\cfrac{df}{dt}+0.03477f\\&\cfrac{d^2\phi_2}{dt^2}+0.1545\cfrac{d\phi_2}{dt}+12.26\phi_1=-0.004226\cfrac{df}{dt}+0.9048f\\&\cfrac{d^2\psi_1}{dt^2}+0.1288\cfrac{d\psi_1}{dt}+2.718\psi_1=0.1757\cfrac{df}{dt}-0.005967f\\&\cfrac{d^2\psi_2}{dt^2}+0.1288\cfrac{d\psi_2}{dt}+2.718\psi_1=0.01049\cfrac{df}{dt}+0.2901f\end{align}

この系のモード減衰比とモード固有角周波数は以下のようになる。

\begin{align} \zeta_1&=0.022059\\\zeta_2&=0.039074\\\omega_1&=3.5020\\\omega_2&=1.6486\end{align}

モード分解による式と比べ若干ずれていることが確認できる。

周波数特性

 モード分解した場合とモード形式に変換した場合の周波数特性を図1に示します。赤で示したモード分解による特性が反共振の部分で若干ずれているがおおむね一致していることがわかる。また、モード形式は元のモデルに対して状態変数を変えただけなので、完全に一致している。

f:id:meckam:20210201225043p:plain

外力から変位\(x_1\)までの周波数特性(青:元のモデル、赤:モード分解、黄:モード形式)

減衰が大きい場合の数値例

次に減衰が大きい場合の運動方程式を考える。先ほどの運動方程式に対して、減衰項が大きくなっていることがわかる。

\begin{align} \left[\begin{array}{cc}3&0\\0&2\end{array}\right]\cfrac{d^2}{dt^2}\left(\begin{array}{c}x_1\\x_2\end{array}\right)&+ \left[\begin{array}{cc}0.5&-1\\-1&25\end{array}\right]\cfrac{d}{dt}\left(\begin{array}{c}x_1\\x_2\end{array}\right)\\&+\left[\begin{array}{cc}30&-10\\-10&10\end{array}\right]\left(\begin{array}{c}x_1\\x_2\end{array}\right)=\left(\begin{array}{c}f\\0\end{array}\right) \end{align}

モード分解

先ほどと同様の手順で同時対角化を行うと次式が得られる。質量行列と剛性行列が先ほどと同じであるため、この時変換行列は先ほどと同じものになる。

\begin{align} \cfrac{d^2\boldsymbol{\phi}}{dt^2}&+ \left[\begin{array}{cc}3.461&-5.472\\-5.472&9.206\end{array}\right]\cfrac{d\boldsymbol{\phi}}{dt}\\&+\left[\begin{array}{cc}12.29&0\\0&2.713\end{array}\right]\boldsymbol{\phi}=\left(\begin{array}{c}-0.5037\\-0.2822\end{array}\right)f\\\left(\begin{array}{c}x_1\\x_2\end{array}\right)&=\left[\begin{array}{cc}-0.5037&-0.2822\\0.3456&-0.6169\end{array}\right] \boldsymbol{\phi}\end{align}

モード形式

こちらも同様の手順で同時対角化を行うと次式が得られる。最終的な状態方程式のみ示す。

\begin{align} \cfrac{d}{dt}\left(\begin{array}{c}\phi_1\\\phi_2\\\psi_1\\\psi_2\end{array}\right)=&\left[\begin{array}{cccc}-12.09&0&0&0\\0&-0.2808&0&0\\0&0&-0.1496&-3.131\\0&0&3.131&-0.1496\end{array}\right]\left(\begin{array}{c}\phi_1\\\phi_2\\\psi_1\\\psi_2\end{array}\right)\\&+\left(\begin{array}{c}0.002338\\0.01398\\-0.3348\\-0.001632\end{array}\right)f\end{align}

上式より、減衰が大きいため実数の固有値が得られたことから、この系の振動モードは1つのみであることがわかる。

周波数特性

先ほどと同様に周波数特性(図2)を観察すると、赤で示したモード分解による特性には全く振動的な応答が得られなくて、うまくモード分解ができていないことが確認できる。

f:id:meckam:20210201230824p:plain

減衰が大きい場合のボード線図(青:元のモデル、赤:モード分解、黄:モード形式)

 

まとめ

減衰が小さい場合にはモード分解で充分にシステムの挙動を再現することが可能であるが、減衰が大きい場合ではモード分解は全く適切ではない。モード分解した際に計算されるモード減衰比が1に近い、1より大きい場合には再考する必要があることに注意せよ。