下三角行列をQR分解したい

モチベーション

三角行列は通常の正方行列と比べ、要素数が\(\frac{n(n+1)}{2}\)(\(n\):行列のサイズ)と少ないため扱いやすい。

それなら通常は数値的に解くことの多いQR分解を、解析的に解いた場合にきれいに解くことができるかも

\(2\times 2\)下三角行列の場合

まずは、最も簡単な場合から解いていく。\(2\times 2\)の下三角行列の各要素を次式のように定義する。

\begin{align} \left[\begin{array}{cc}l_{11}&0\\l_{21}&l_{22}\end{array}\right]\end{align}

QR分解の方法としては、「グラム・シュミット法」「ハウスホルダー変換」「ギブンス変換」があるが、解析的に解くにはギブンス変換が良さそう。ギブンス変換はQR分解したい行列を左から回転行列をかけることで行う方法である。

まずは次の角度を求める。

\begin{align} \theta_{21}&=\tan^{-1}\left(\frac{-l_{21}}{l_{11}}\right)\\\sin\theta_{21}&=-\frac{l_{21}}{\alpha_{21}}\\\cos\theta_{21}&=\frac{l_{11}}{\alpha_{21}}\\\alpha_{21}&=\sqrt{l_{11}^2+l_{21}^2}\end{align}

次のように行列を演算する。なお\(\sin\theta\)を\(s\theta\)、\(\cos\theta\)を\(c\theta\)と表記する。

\begin{align} \left[\begin{array}{cc}c\theta_{21}&-s\theta_{21}\\s\theta_{21}&c\theta_{21}\end{array}\right] \left[\begin{array}{cc}l_{11}&0\\l_{21}&l_{22}\end{array}\right]=\left[\begin{array}{cc}\alpha_{21}&\frac{l_{21}l_{22}}{\alpha_{21}}\\0&\frac{l_{11}l_{22}}{\alpha_{21}}\end{array}\right] \end{align}

以上のように\(2\times 2\)下三角行列のQR分解が得られた。

\(3\times 3\)下三角行列の場合

\(3\times 3\)の下三角行列の各要素を次式のように定義する。

\begin{align} \left[\begin{array}{ccc}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\end{array}\right]\end{align}

まずは先ほどと同様の手順で、左上\(2\times 2\)の下三角行列のQR分解を求める。

\begin{align} \theta_{21}&=\tan^{-1}\left(\frac{-l_{21}}{l_{11}}\right)\\s\theta_{21}&=-\frac{l_{21}}{\alpha_{21}}\\c\theta_{21}&=\frac{l_{11}}{\alpha_{21}}\\\alpha_{21}&=\sqrt{l_{11}^2+l_{21}^2}\end{align}

次のように行列を演算する。

\begin{align} \left[\begin{array}{ccc}c\theta_{21}&-s\theta_{21}&0\\s\theta_{21}&c\theta_{21}&0\\0&0&1\end{array}\right] \left[\begin{array}{ccc}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\end{array}\right]=\left[\begin{array}{ccc}\alpha_{21}&\frac{l_{21}l_{22}}{\alpha_{21}}&0\\0&\frac{l_{11}l_{22}}{\alpha_{21}}&0\\l_{31}&l_{32}&l_{33}\end{array}\right] \end{align}

次に(3,1)要素をゼロにする。ここで(3,2)要素をゼロにしようとすると、(2,1)要素が非ゼロになってしまう。

まずは次の角度を求める。

\begin{align} \theta_{31}&=\tan^{-1}\left(\frac{-l_{31}}{\alpha_{21}}\right)\\s\theta_{31}&=-\frac{l_{31}}{\alpha_{31}}\\c\theta_{31}&=\frac{\alpha_{21}}{\alpha_{31}}\\\alpha_{31}&=\sqrt{\alpha_{21}^2+l_{31}^2}=\sqrt{l_{11}^2+l_{21}^2+l_{31}^2}\end{align}

次のように行列を演算する。

\begin{align} \left[\begin{array}{ccc}c\theta_{31}&0&-s\theta_{31}\\0&1&0\\s\theta_{31}&0&c\theta_{31}\end{array}\right] \left[\begin{array}{ccc}\alpha_{21}&\frac{l_{21}l_{22}}{\alpha_{21}}&0\\0&\frac{l_{11}l_{22}}{\alpha_{21}}&0\\l_{31}&l_{32}&l_{33}\end{array}\right]=\left[\begin{array}{ccc}\alpha_{31}&\frac{l_{21}l_{22}+l_{31}l_{32}}{\alpha_{31}}&\frac{l_{31}l_{33}}{\alpha_{31}}\\0&\frac{l_{11}l_{22}}{\alpha_{21}}&0\\0&\frac{\alpha_{21}^2l_{32}-l_{21}l_{31}l_{22}}{\alpha_{21}\alpha_{31}}&\frac{\alpha_{21}l_{33}}{\alpha_{31}}\end{array}\right] \end{align}

一方、回転行列の積は次式のようになる。

\begin{align} \left[\begin{array}{ccc}c\theta_{31}&0&-s\theta_{31}\\0&1&0\\s\theta_{31}&0&c\theta_{31}\end{array}\right] \left[\begin{array}{ccc}c\theta_{21}&-s\theta_{21}&0\\s\theta_{21}&c\theta_{21}&0\\0&0&1\end{array}\right]=\left[\begin{array}{ccc}c\theta_{21}c\theta_{31}&-s\theta_{21}c\theta_{31}&-s\theta_{31}\\s\theta_{21}&c\theta_{21}&0\\c\theta_{21}s\theta_{31}&-s\theta_{21}s\theta_{31}&c\theta_{31}\end{array}\right]\end{align}

この時点で大分ややこしい…

最後に、(3,2)要素をゼロにするために次の角度を求める。

\begin{align} \theta_{32}&=\tan^{-1}\left(\frac{l_{21}l_{31}l_{22}-\alpha_{21}^2l_{32}}{\alpha_{31}l_{11}l_{22}}\right)\\s\theta_{31}&=\frac{l_{21}l_{31}l_{22}-\alpha_{21}^2l_{32}}{\alpha_{32}}\\c\theta_{31}&=\frac{\alpha_{31}l_{11}l_{22}}{\alpha_{32}}\\\alpha_{32}&=\sqrt{\alpha_{31}^2l_{11}^2l_{22}^2+\left(l_{21}l_{31}l_{22}-\alpha_{21}^2l_{32}\right)^2}\end{align}

次のように行列を演算する。

\begin{align} \left[\begin{array}{ccc}1&0&0\\0&c\theta_{32}&-s\theta_{32}\\0&s\theta_{32}&c\theta_{32}\end{array}\right] \left[\begin{array}{ccc}\alpha_{31}&\frac{l_{21}l_{22}+l_{31}l_{32}}{\alpha_{31}}&\frac{l_{31}l_{33}}{\alpha_{31}}\\0&\frac{l_{11}l_{22}}{\alpha_{21}}&0\\0&\frac{\alpha_{21}^2l_{32}-l_{21}l_{31}l_{22}}{\alpha_{21}\alpha_{31}}&\frac{\alpha_{21}l_{33}}{\alpha_{31}}\end{array}\right]=\left[\begin{array}{ccc}\alpha_{31}&\frac{l_{21}l_{22}+l_{31}l_{32}}{\alpha_{31}}&\frac{l_{31}l_{33}}{\alpha_{31}}\\0&\frac{\alpha_{32}}{\alpha_{21}\alpha_{31}}&\frac{\alpha_{21}l_{33}\left(\alpha_{21}^2l_{32}-l_{21}l_{31}l_{22}\right)}{\alpha_{31}\alpha_{32}}\\0&0&\frac{\alpha_{21}l_{11}l_{22}l_{33}}{\alpha_{32}}\end{array}\right] \end{align}

一方、回転行列の積は次式のようになる。

\begin{align} \left[\begin{array}{ccc}1&0&0\\0&c\theta_{32}&-s\theta_{32}\\0&s\theta_{32}&c\theta_{32}\end{array}\right]\left[\begin{array}{ccc}c\theta_{21}c\theta_{31}&-s\theta_{21}c\theta_{31}&-s\theta_{31}\\s\theta_{21}&c\theta_{21}&0\\c\theta_{21}s\theta_{31}&-s\theta_{21}s\theta_{31}&c\theta_{31}\end{array}\right]\\=\left[\begin{array}{ccc}c\theta_{21}c\theta_{31}&-s\theta_{21}c\theta_{31}&-s\theta_{31}\\s\theta_{21}c\theta_{32}-c\theta_{21}s\theta_{31}s\theta_{32}&c\theta_{21}c\theta_{32}+s\theta_{21}s\theta_{31}s\theta_{32}&-c\theta_{31}s\theta_{32}\\s\theta_{21}s\theta_{32}+c\theta_{21}s\theta_{31}c\theta_{32}&c\theta_{21}s\theta_{32}-s\theta_{21}s\theta_{31}c\theta_{32}&c\theta_{31}c\theta_{32}\end{array}\right]\end{align}

ちょっともうきついかも…

\(4\times 4\)下三角行列の場合

\(4\times 4\)の下三角行列の各要素を次式のように定義する。

\begin{align} \left[\begin{array}{cccc}l_{11}&0&0&0\\l_{21}&l_{22}&0&0\\l_{31}&l_{32}&l_{33}&0\\l_{41}&l_{42}&l_{43}&l_{44}\end{array}\right]\end{align}

途中までは、\(3\times 3\)の場合と同様の手順で下記のように得られる。

\begin{align}\left[\begin{array}{ccc}c\theta_{21}c\theta_{31}&-s\theta_{21}c\theta_{31}&-s\theta_{31}&0\\s\theta_{21}&c\theta_{21}&0&0\\c\theta_{21}s\theta_{31}&-s\theta_{21}s\theta_{31}&c\theta_{31}&0\\0&0&0&1\end{array}\right]\left[\begin{array}{cccc}l_{11}&0&0&0\\l_{21}&l_{22}&0&0\\l_{31}&l_{32}&l_{33}&0\\l_{41}&l_{42}&l_{43}&l_{44}\end{array}\right]\\=\left[\begin{array}{ccc}\alpha_{31}&\frac{l_{21}l_{22}+l_{31}l_{32}}{\alpha_{31}}&\frac{l_{31}l_{33}}{\alpha_{31}}&0\\0&\frac{l_{11}l_{22}}{\alpha_{21}}&0&0\\0&\frac{\alpha_{21}^2l_{32}-l_{21}l_{31}l_{22}}{\alpha_{21}\alpha_{31}}&\frac{\alpha_{21}l_{33}}{\alpha_{31}}&0\\l_{41}&l_{42}&l_{43}&l_{44}\end{array}\right] \end{align}

次に(4,1)要素をゼロにする。

\begin{align} \theta_{41}&=\tan^{-1}\left(\frac{-l_{41}}{\alpha_{31}}\right)\\s\theta_{41}&=-\frac{l_{41}}{\alpha_{41}}\\c\theta_{41}&=\frac{\alpha_{31}}{\alpha_{41}}\\\alpha_{41}&=\sqrt{\alpha_{31}^2+l_{41}^2}=\sqrt{l_{11}^2+l_{21}^2+l_{31}^2+l_{41}^2}\end{align}

次のように行列を演算する。

\begin{align} \left[\begin{array}{cccc}c\theta_{41}&0&0&-s\theta_{41}\\0&1&0&0\\0&0&1&0\\s\theta_{41}&0&0&c\theta_{31}\end{array}\right] \left[\begin{array}{ccc}\alpha_{31}&\frac{l_{21}l_{22}+l_{31}l_{32}}{\alpha_{31}}&\frac{l_{31}l_{33}}{\alpha_{31}}&0\\0&\frac{l_{11}l_{22}}{\alpha_{21}}&0&0\\0&\frac{\alpha_{21}^2l_{32}-l_{21}l_{31}l_{22}}{\alpha_{21}\alpha_{31}}&\frac{\alpha_{21}l_{33}}{\alpha_{31}}&0\\l_{41}&l_{42}&l_{43}&l_{44}\end{array}\right]\\=\left[\begin{array}{ccc}\alpha_{41}&\frac{l_{21}l_{22}+l_{31}l_{32}+l_{41}l_{42}}{\alpha_{41}}&\frac{l_{31}l_{33}+l_{41}l_{43}}{\alpha_{41}}&\frac{l_{41}l_{44}}{\alpha_{41}}\\0&\frac{l_{11}l_{22}}{\alpha_{21}}&0&0\\0&\frac{\alpha_{21}^2l_{32}-l_{21}l_{31}l_{22}}{\alpha_{21}\alpha_{31}}&\frac{\alpha_{21}l_{33}}{\alpha_{31}}&0\\0&\frac{\alpha_{31}^2l_{42}-l_{41}\left(l_{21}l_{22}+l_{31}l_{32}\right)}{\alpha_{31}\alpha_{41}}&\frac{\alpha_{31}^2l_{43}-l_{31}l_{41}l_{33}}{\alpha_{31}\alpha_{41}}&\frac{\alpha_{31}l_{44}}{\alpha_{41}}\end{array}\right]\end{align}

ちょっとこれ以上解析的に解くのはきつい。ただ、第1行はきれいになるので、その部分の計算だけは簡単に計算できそう。

やはりQR分解は数値的に解くのが良さそうです。

ナイキストの安定判別法の解釈

背景

制御工学における安定性を検証するために、ナイキストの安定判別法が良く用いられる。ナイキストの安定判別法のメリットとして、

• モデル化が必要なく、周波数応答さえ取得できれば安定性の判別が可能
• 連続系と離散系、無駄時間が複合するようなシステムでも判定可能
• 安定余裕という指標を用いることで、安定性の強さの定量化が可能

といったように多くある。しかしこの安定判別の方法が、

• 一巡伝達関数の極に不安定極がない場合、一巡伝達関数の0[Hz]から∞[Hz]（離散系の場合はナイキスト周波数）における周波数応答を複素平面状に描画した線図（ナイキスト線図）が、-1+0jの点を左側に見て進めば安定
• 不安定極が\(P\)個ある場合、一巡伝達関数の-∞[Hz]から∞[Hz]のナイキスト線図が、-1+0jの点を反時計回りにP回回れば安定

といったような、パッと見わかりづらい方法である。この方法の最大の問題点は描画した線図を人の目が見て判断する必要があることである。つまり、もしプログラムで計測した周波数特性の安定判別を行う際は、このナイキスト安定判別法を数式等に直す必要があります。

不安定な極を持たない場合

まずは次のような一巡伝達関数を考えてみます。

\begin{align}G(s)=K\cfrac{s-1}{s^2+6s+5}\end{align}

\(K=1\)の時のナイキスト線図は下図のようになります。Octaveでは実線が正の周波数領域を表していることに注意しよう。

K=1のナイキスト線図

上記のようにナイキスト線図が-1+0jを回らないため、この系の閉ループ系は安定である。

次に\(K=10\)の時のナイキスト線図は下図のようになります。

K=10のナイキスト線図

0Hzの時は\(G(j0)=-2\)であるため、-1+0jを右に見ながら原点に収束しており、この系の閉ループ系は不安定である。

この問題をプログラムで判定させるために、一巡伝達関数に1+0jを加えた伝達関数を考える。

\begin{align}F(s)=1+G(s)=1+K\cfrac{s-1}{s^2+6s+5}\end{align}

つまりナイキスト線図を1だけ左にずらした線図、-1を原点に持つ線図に変更する。こうすることで原点から見た位相が時計回りに回るかどうかで判定できるようになるため、ボード線図を書いたときの位相から判定できるようになる。

\(K=1\)の時の\(F(s)\)のボード線図は下図のようになります。

K=1のボード線図

このように位相が0degを出発して、0degに収束しており回転していないため安定であるといえる。

次に\(K=10\)の時の\(F(s)\)のボード線図は下図のようになります。

K=10のボード線図

上図は位相が180degを出発して、0degに収束しており半回転している。周波数領域を負の領域まで考えた場合、位相が360degを出発して0degを収束していることになるため、反時計回りに1周しており、不安定である。

このように一巡伝達関数に1をプラスした位相線図から安定性を判別できる。以降では別のパターンで考える。

不安定な極を持つ場合

次のような一巡伝達関数を考えてみます。不安定極が1つあるため、ナイキスト線図は反時計回りに1周する必要があります。

\begin{align}G(s)=K\cfrac{s+1}{s^2+4s-5}\end{align}

\(K=1\)の時のナイキスト線図と\(F(s)\)のボード線図は下図のようになります。

K=1のナイキスト線図

K=1のボード線図

位相線図の位相差が0degであり、不安定極を持つのに回転していないためこの系の閉ループ系は不安定である。

\(K=10\)の時のナイキスト線図と\(F(s)\)のボード線図は下図のようになります。

K=10のナイキスト線図

K=10のボード線図

位相線図の位相差が180degであり、反時計回りとなるためには周波数が増えることで位相が増えていく方向である。そのため、閉ループ系は安定である。

不安定な極と原点に極を持つ場合

次のような一巡伝達関数を考えてみます。不安定極（原点極を含む）が2つあるため、ナイキスト線図は反時計回りに2周する必要があります。

\begin{align}G(s)=K\cfrac{s+1}{s(s^2+4s-5)}\end{align}

\(K=1\)の時のナイキスト線図と\(F(s)\)のボード線図は下図のようになります。

K=1のナイキスト線図

K=1のボード線図

虚軸上に極を持つ場合、その周波数におけるナイキスト線図の不連続性は右側に回転して繋がっていると考えることとなっている。

位相線図の位相差が90degしかなく、不安定極を持つのに回転していないためこの系の閉ループ系は不安定である。

\(K=10\)の時のナイキスト線図と\(F(s)\)のボード線図は下図のようになります。

K=10のナイキスト線図

K=10のボード線図

位相線図の位相差が270degではあるが、虚軸上の極を持つ場合はむりやり右側に回転して繋がっていると考えるため、周波数0Hzの時の位相は90degだけ引いて0degとなる。そうなると位相差が360degであるため、0Hz~∞Hzで1周回っており、-∞Hz~0Hzで1周回っていることになるため、全体として2周回っており安定である。

まとめ

このように計測した周波数特性に対して、1+0jを加えた特性の位相差を見ることで安定判別が可能となる。しかしこの場合でも不安定極の個数が事前にわかっている必要があるため、完全にモデルなしで安定判別ができるわけではないことに注意したい。

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

背景

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

モータの基礎方程式

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

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

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

アナロジー視点

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

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

\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のようになる。

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

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