【RESPチーム四方山話】減衰について本気出して考えてみた |

振動解析をするときの減衰

振動解析をするときに、減衰の設定は重要です。しかし、その一方で、「減衰」というものは「振動方程式のおさまりのために考慮しよう」といった値でもあります。ですので、そこまで細かく精査するものではなく、「まぁこんなものだろう」って設定される傾向が強いと思われます。そういう背景はありますが、たまには、慣用的に与えてる値が振動方程式にどのように影響するかを考えてみるのも良いと思いこの記事を作成しました。従いまして、次のような方に読んでもらいたいものです。

減衰定数を「剛性比例減衰やレーリー減衰を指定する時のパラメータであり、大きくすれば早く振動が収まる」程度に認識している。
「剛性比例減衰は高次モードの影響が大きいと過大評価になる」という表現をよく使うが、正直、よく見るグラフを見て感覚的な理解になってる気がしてる。
減衰マトリクスと友達になりたい。

減衰定数hについて

まずは減衰定数hをおさらいしておきます。一自由度系の振動系があり、質量をm、減衰係数をc、固有円振動数をωとすると、減衰定数は次のように表せます。

$$
h = \frac{c}{2 \cdot m \cdot \omega} \tag1
$$

減衰係数を基準に考えると次のようになります。

$$
c = 2 h \cdot m \cdot \omega \tag2
$$

なお、多くのソフトで立体解析などでの多自由度系の解析のパラメータとなっているので意識してない事もあるかと思いますが、減衰定数は固有周期と同じで「一自由度の振動系に関する指標」です。

比例減衰マトリクスについて

振動解析をするときは「比例減衰」で考えることが多いです。よく教科書に載っている比例減衰は以下の種類があります。

質量比例減衰マトリクス
剛性比例減衰マトリクス
レーリー減衰マトリクス
モード別減衰マトリクス

大体の本では上記の順番で具体的な式が並んでいることが多いです。よく用いられる剛性比例減衰は全体剛性マトリクスを係数倍しただけの形で表現されるのでわかりやすく、比例減衰マトリクスの説明の前半部分に記載されていることが多いです。単純に式だけが載っている場合もあり、それぞれの減衰マトリクスの関係性がわかりづらいこともありますが、「剛性比例減衰、質量比例減衰、レーリー減衰は全てモード別減衰の特別なパターンである」と考えることができます。

減衰マトリクスの対角化

比例減衰で考える範囲は、減衰マトリクスが対角化できるという仮定で構築されています。次のように固有ベクトルを前後から乗じた場合に対角成分に「モード減衰」が並びます。

$$
[U]^T [C] [U] =
\begin{bmatrix}
_sC_1 & 0 & \cdots & 0 \\
0 & _sC_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & _sC_n
\end{bmatrix} \tag3
$$

このモード減衰は減衰係数であり、式（1）に入れることで、そのモードの減衰定数となります。

モード別減衰マトリクスについて

普通の説明では簡単な「剛性比例減衰」から説明しますが、あえて「モード別減衰マトリクス」から説明します。「モード別減衰マトリクス」は各モードの減衰定数が指定した値になる減衰マトリクスです。つまり対角化した時の「モード減衰」が各次で式(2)の関係となるように作成します。

$$
\begin{eqnarray}
[U]^T [C] [U] &=&
\begin{bmatrix}
2 h_1 \cdot _s m_1 \cdot \omega_1 & 0 & \cdots & 0 \\
0 & 2 h_2 \cdot _s m_2 \cdot \omega_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & 2 h_n \cdot _s m_n \cdot \omega_n
\end{bmatrix} \\\
&=&
\color{Red}{ \underline{ \color{black}{
\begin{bmatrix}
_s m_1 & 0 & \cdots & 0 \\
0 & _s m_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & _s m_n
\end{bmatrix}}}}
\color{blue}{ \underline{ \color{black}{ \begin{bmatrix}
2 h_1 \cdot \omega_1 & 0 & \cdots & 0 \\
0 & 2 h_2 \cdot \omega_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & 2 h_n \cdot \omega_n
\end{bmatrix} }}}\\\
&=&
\color{Red}{ \underline{ \color{black}{ [U]^T[M][U] }}} \color{blue}{ \underline{ \color{black}{ [\zeta]}}}
\end{eqnarray} \tag4
$$

$[U]^T [C] [U]$の前から$[U]^{-T}$を掛けて、後ろから$[U]^{-1}$を掛けて、さらに$[U]^{-1} = ([U]^T [M] [U])^{-1} \cdot [U]^T \cdot [M] $からモード別減衰マトリクスは次のようになります。

$$
[C] = [M] \cdot [U] \cdot [\zeta] \cdot([U]^T [M] [U])^{-1} \cdot [U]^T \cdot [M] \tag5
$$

モード別減衰マトリクスを分解してみる

モード別減衰マトリクスをもう少し深堀りしてみます。先ほどの式（5）を各モードの固有ベクトル$\lbrace u_i \rbrace$で表現してみます。

$$
[C_i] = [M] \cdot \{u_i \} \cdot \frac{2 h_i \omega_i }{_sm_i} \cdot \{u_i \}^T \cdot [M] \tag6
$$

$[C_i]$はすべてのモードで存在するので、すべての次数で計算した$[C_i]$を足したものが全体減衰マトリクスになります。

$$
[C] = [C_1] + [C_2] + \cdots + [C_n] \tag7
$$

この各次の$[C_i]$は全体減衰マトリクスを「スペクトル分解」した時の各モード成分を表してます。意味合いとしては、応答速度を$i$次モードの減衰力にする作用を持っています。例えば、ある時刻で応答速度$\lbrace \dot{x} \rbrace$が生じた場合に、1次モードの減衰マトリクスが$[C_1]$にだったとします。まず、応答速度は振動系の固有モードに展開できます。

$$
\lbrace \dot{x} \rbrace =\alpha_1 \lbrace u_1 \rbrace + \alpha_2 \lbrace u_2 \rbrace + \cdots + \alpha_n \lbrace u_n \rbrace \tag8
$$

この、応答速度が$[C_1]$に作用して発生する減衰力は次のようになります。

$$
\begin{eqnarray}
[C_1] \lbrace \dot{x} \rbrace &=& [M] \cdot \frac{2h_1 \omega_1}{_sM_1} \lbrace u_1 \rbrace \color{Red}{ \underline{ \color{black}{\lbrace u_1 \rbrace^T \cdot [M] \cdot (\alpha_1 \lbrace u_1 \rbrace + \alpha_2 \lbrace u_2 \rbrace + \cdots + \alpha_n \lbrace u_n \rbrace )}}} \\\
&=& [M] \cdot \frac{2h_1 \omega_1}{_sM_1 }\lbrace u_1 \rbrace \cdot \color{Red}{ \underline{ \color{black}{\alpha_1 \cdot _sM_1}}} \\\
&=& 2h_1 \cdot \omega_1 \cdot [M] \cdot \alpha_1 \cdot \lbrace u_1 \rbrace
\end{eqnarray} \tag9
$$

赤線の部分は直交性を利用して、全体座標系の応答から1次モードの速度$\alpha_1 \cdot \lbrace u_1 \rbrace$を抽出しています。
行列$[C_i]$はすべての次数に対して存在しますが、ある次数の減衰定数を$h_i=0$のように設定したモードについては$[C_i] = [0]$となり、そのモードに対する減衰力が発生しなくなります。このようにして見ると、この記事での内容が理解しやすくなりませんか。

剛性比例減衰について

剛性比例減衰も同様に表すことができるので過程を説明します。
剛性比例減衰は全体減衰マトリクスは全体剛性マトリクスをスカラー倍したものとして定義されます。

$$
[C] = \alpha [K] \tag{10}
$$

剛性マトリクスは固有ベクトルで対角化が可能です。

$$
\begin{eqnarray}
[U]^T [K] [U] &=&
\begin{bmatrix}
\{ u_{1} \}^{T} \\ \{ u_{2} \}^{T} \\ \vdots \\ \{ u_{n} \}^{T}
\end{bmatrix}
\begin{bmatrix}
k_{11} & k_{12} & \cdots & 0 \\
k_{21} & k_{22} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & k_{nn}
\end{bmatrix}
\begin{bmatrix}
\{ u_{1} \} & \{ u_{2} \} & \cdots & \{ u_{n} \}
\end{bmatrix} \\\
&=&
\begin{bmatrix}
k_1 & 0 & \cdots & 0 \\
0 & k_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & k_n
\end{bmatrix}
\end{eqnarray}
$$

したがって、式（10）を対角化すると式（11）のようになります。

$$
\begin{eqnarray}
[U]^T [C] [U] &=&
\begin{bmatrix}
\alpha \cdot k_1 & 0 & \cdots & 0 \\
0 &\alpha \cdot k_2 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots\\
0 & 0 & \cdots & \alpha \cdot k_n
\end{bmatrix}
\end{eqnarray} \tag{11}
$$

$k_i$は$i$次モードのモード剛性です。モード剛性はモード質量などと同じで、固有ベクトルを用いて等価1自由度系に座標変換した時の剛性となります。式(11)の各モード減衰を式（1）に入れると、$i$次モードの減衰定数$h_i$は以下のように表せます。

$$
\begin{align}
h_i & = \frac{\alpha \cdot k_i}{2 \cdot m_i \cdot \omega_i} \\\
& = \frac{\alpha \cdot \omega_i^2}{2 \cdot \omega_i} \\\
& = \frac{\alpha}{2} \omega_i
\end{align} \tag{12}
$$

この式は減衰定数と固有円振動数が比例関係にあることを表してます。「質量比例減衰」も同じように計算すれば反比例関係、「レーリー減衰」も両者を足し合わせた関係になりますので是非試してみてください。

剛性比例減衰の高次モードの減衰について

式（7）の形式は比例減衰では全て成立するので、もちろん剛性比例減衰でも成り立ちます。すべての次数で式（12）の関係性になるように行列$[C_i]$を作成して、足し合わせたら剛性比例減衰となります。
高次モードでは固有円振動数が大きくなるので、大きい減衰定数で設定された$[C_i]$を含んでいます。言い方を変えれば、1次モードに対して2%とか設定したつもりでも、勝手に高次モードに大きい減衰定数で設定された$[C_i]$が設定されています。
高次モードで式（9）での計算のイメージができると、大きい減衰定数の$[C_i]$によって大きな減衰力が発生して、よく言う「高次の減衰を過大評価」する状態もわかりやすいかもしれません。
なお、モード別減衰マトリクスを作成するには減衰定数を指定するモードの固有円振動数と固有ベクトルを知る必要がありますが、立体の固有値解析は基本的にすべての次数で計算するわけではありません。その意味では、剛性比例減衰は作成のコストも軽いというメリットがあります。

〇次に対して△%の減衰を与える操作について

剛性比例減衰を定義する時は最初に「〇次の周期に対して△%の減衰を与える」という呪文を唱えます。
この呪文を言い換えると、「比例関係にある減衰定数、円振動数の傾き$\alpha$を決定する」という事になります。例えば、1次（固有円振動数が20(rad/s))に対して5%減衰を与えた場合です。
比例関係になるので、減衰定数と円振動数のグラフは次のようになります。

このグラフの傾きは0.0025です。
つまり、先ほどの式（12）でいうと次のようになります。

$$
h_i = 0.0025 \cdot \omega_i \tag{13}
$$

そして、$\frac{\alpha}{2}=0.0025$なので$\alpha=0.005$となります。
この$\alpha$を全体剛性マトリクスに乗じれば、剛性比例減衰マトリクスが完成します。
剛性比例型減衰は式（14）のように定式化されてますが、上記の計算過程を数式化していることがわかります。$\alpha$が決まればよいので、呪文を唱える時は解析モデルの1次などの固有周期である必要ではありません。このように設定するのは、想定している円振動数（例えば、1次の円振動数）での減衰定数を指定するのにわかりやすいためです。また、当然ですがグラフの傾きを指定するので、円振動数が40(rad/s)に対して10%を与えても同じ減衰を与えていることになります。

$$2 \cdot \frac{h}{ω}・[K] \tag{14} $$

閑話～剛性比例減衰であることの恩恵～

プログラムを作る側からすると、剛性比例減衰はとてもありがたいものです。
$[C_i]$は基本的にフルマトリクスとなるため、モード別減衰にはフルマトリクスとなります。
一方で、剛性マトリクスはその成分のほとんどが0であり、対角成分付近に非ゼロの値を持つため、その特性に合わせたメモリの確保や行列演算を実装しているプログラムがほとんどです。
剛性マトリクスをスカラー倍しただけという事は、剛性比例減衰マトリクスは行列の形としても同じ性質を持っているため、高速に連立方程式を解くことが可能です。行列が大規模になればなるほど、この恩恵はすさまじいものになります。
そんなに変わるものなのか？と感じるかもしれませんが、連立方程式を解くスピードは係数行列が持っている性質にかなり依存します。「この形状の係数行列に対するアルゴリズム」ということで関数が用意されるくらいなので、その性質から外れた行列を解こうとすると信じられないくらい遅くなったという経験は数値解析プログラムを開発してる人の中ではあるあるかもしれません。
RESPチームではこれに対策するため「不完全減衰マトリクス」という考え方を提案しています。過去の建築学会で投稿されているので、興味のある方はぜひ確認してみてください。RESP-D、RESP-F3Tで研究的な機能として使用することもできます。
ちなみに、剛性比例減衰で式（7）の計算をしてみると、フルマトリクス$[C_i]$を足し算した結果で最終的な$[C]$はきれいに剛性マトリクスの係数倍になるため、思わず拍手してしまいます。※1つのモードだけでも、少しずらすとフルマトリクスになってしまいます。

まとめ

今回は、減衰マトリクスについて長尺で記載してみました。総じて次のようなことを述べてました。

比例減衰マトリクスは各モードに対する減衰マトリクス$[C_i]$に分解できて、i次モードの応答に対応する減衰力を計算する行列になる。
「質量比例減衰」「剛性比例減衰」「レーリー減衰」は減衰定数と円振動数の関係性を仮定して、その関係性になるように「モード別減衰マトリクス」を作成したものであり、すべての次数に減衰定数が自動的にセットされる。
〇次に対して△%の減衰を与える操作は減衰定数と円振動数の関係性を設定している。

なんとなく「〇造なら1次に対して、△％与えよう」と慣例的に思うだけでなく、たまには、自分が設定した減衰によってどんな減衰マトリクスができているかに思いをはせてみてはいかがでしょうか。
これで減衰マトリクスと友達になれましたかね。