Pythonでレイリー減衰を実装する方法!解析して効果を見てみる

  • このエントリーをはてなブックマークに追加
Advertisements

シミュレーション分野にはレイリー減衰という減衰の考え方があります。ここではPythonでレイリー減衰の曲線を書いたり、シミュレーションで使って効果を確かめたりしながら学んでいきます。

こんにちは。wat(@watlablog)です。今回はレイリー減衰って何なの?どう使うの?というのを調べてみました

レイリー減衰の概要

色々な減衰とレイリー減衰の特徴

「減衰」と一口に言っても様々な種類があります。

振動分野で度々出てくる、速度に比例した減衰力を発生させる粘性減衰(Viscous damping)、変位に比例する構造減衰(Structural damping)、振動モードで考えるモーダル減衰(Modal damping)…と、線形の範囲でざっと考えただけでこれだけあります。

レイリー減衰(Rayleigh damping)とは、数値解析の世界から発展した減衰の考え方であり、一般に計算を安定させる目的で使用されます。

また、レイリー減衰は固有値を基準に周波数範囲で減衰を与える事ができる特徴から、運動方程式を直接解く場合にもモーダル減衰と似たような使い方が出来る事で重宝されています。

通常の粘性減衰、構造減衰(これはヒステリシス減衰とも呼ぶ)はパラメータを設定する時に周波数軸に対する減衰度合いはなかなかイメージする事が出来ませんので、レイリー減衰はその部分に関して他の減衰と比べて異色と言えます。

但し、物理的な因果関係から成り立っている減衰では無いため、使用には注意が必要です。

レイリー減衰の構成式

レイリー減衰を説明するために、はじめに運動方程式を式(1)に示します。

\[ \mathbf{M}\mathbf{\ddot{x}}+\mathbf{C}\mathbf{\dot{x}}+\mathbf{K}\mathbf{x}=\mathbf{F} \tag{1} \]

この式の減衰マトリクス\(\mathbf{C}\)に対して次式の通り質量マトリクス\(\mathbf{M}\)と剛性マトリクス\(\mathbf{K}\)にそれぞれ比例定数\(\alpha, \beta\)をかけて加算する形にしたものがレイリー減衰(式(2))です。

\[ \mathbf{C}=\alpha \mathbf{M} + \beta \mathbf{K}\tag{2} \]

さらに、レイリー減衰では比例定数\(\alpha, \beta\)と減衰比\(\zeta\)の間に、モード毎の固有角振動数\(\omega_{i}\)を使って次式(3)で表します。

\[ \zeta_{i} = \frac{1}{2}(\frac{\alpha}{\omega_{i}}+\beta \omega_{i})\tag{3} \]

この式は固有角振動数が低周波数の場合は\(\alpha\)の方(=質量マトリクスにかかる係数)が大きくなるのに対し、高周波数の場合は\(\beta\)の方(=剛性マトリクスにかかる係数)が大きくなる事を意味しています。

一般にこの\(\alpha, \beta\)は二つの固有角振動数の範囲と対応する減衰比を決めて使用するため、以下式(4)のように二つの\(\zeta\)と二つの\(\omega\)を設定します。

\[ \begin{cases} \displaystyle \zeta_{1} = \frac{1}{2}(\frac{\alpha}{\omega_{1}}+\beta \omega_{1})\\ \displaystyle \zeta_{2} = \frac{1}{2}(\frac{\alpha}{\omega_{2}}+\beta \omega_{2}) \end{cases}\tag{4} \]

この2式は連立方程式として\(\alpha, \beta\)で解く事が可能です。
一応本ブログはPythonブログなので、PythonライブラリであるSymPyに解かせてみましょう。
コードは以下。

すると以下の結果を得ます。簡単に記号ベースの連立方程式を解く事が出来ました。

LaTeX形式で書くと式(5)。SymPyは便利ですね。

\[ \begin{cases} \displaystyle \alpha = \frac{2\omega_{1}\omega_{2}(\omega_{1}\zeta_{2}-\omega_{2}\zeta_{1})}{\omega_{1}^{2}-\omega_{2}^{2}}\\ \displaystyle \beta = \frac{2(\omega_{1}\zeta_{1}-\omega_{2}\zeta_{2})}{\omega_{1}^{2}-\omega_{2}^{2}} \end{cases}\tag{5} \]

この式を使って\(\alpha, \beta\)を決めれば、固有角振動数\(\omega_{1}, \omega_{2}\)の範囲で任意の減衰比\(\zeta_{1}, \zeta_{2}\)を滑らかに設定する事ができます。

ほんと??やってみましょう!

Pythonでレイリー減衰曲線を描くコード

Advertisements

以下がPythonでレイリー減衰曲線を描画するコードです。詳細はコメントに書きましたので是非ご確認下さい。
generate_rayleigh_damping()でalphaとbetaを決定し、check_rayleigh_damping()はただグラフ表示をしているだけです。

以下が結果です。文献等でよく見る曲線を描く事ができました。

Rayleigh damping curve with alpha and beta curve

Pythonの振動解析でレイリー減衰の効果を試してみる

Advertisements

シミュレーションモデル

当ブログではいつも直列のばねマスモデルを使っていましたが、今回シミュレーションモデルはサンプルとして以下の7自由度振動系を使います。

モデルは関数として定義します。剛性行列は「Python/SymPyで連立運動方程式の剛性行列を自動生成する」で自動生成した行列を使いました。この関数の中でgenerate_rayleigh_damping()関数を実行してレイリー減衰を考慮しています。

このモデルでは以下のレイリー減衰を設定しています。

サンプルのレイリーダンピング曲線

シミュレーションコード

Runge-Kutta法

シミュレーションにはRunge-Kutta法を使います。この方法は「Pythonの4次ルンゲ・クッタ法で多自由度連成振動を解く方法」を見て頂くと詳細がわかると思います。以下のコードを「simulation.py」として別ファイルで作成しておきます。

メインコード全文

先ほどの「simulation.py」とは別に、こちらのメインコードを作成し、実行します。

シミュレーションモデルの関数もここに記載しています。ただシミュレーションするだけでなく、固有値解析、周波数応答関数(FRF)計算、グラフ表示を行っています。

固有値解析については「Pythonで多自由度系の固有値解析!固有振動数とモードを計算」を、周波数応答関数については「Pythonで伝達関数!2つの信号の周波数応答関数を計算する方法」を参照して下さい。

実行結果

以下がコードの実行結果です。インパルス応答としてFRF上に7つのピーク(7自由度分)が発生していますが、それぞれのピークの丸み、鈍り方が異なっており、減衰比がモード毎に異なるという結果を得ました。

シミュレーション結果

おまけ:減衰比を変えてみる

以下のように減衰比を変更してみます。

しっかり変更が反映されています。

減衰比を変更したシミュレーション結果

まとめ

レイリー減衰について、ちまたでよく聞くので筆者は存在は知っていましたが実際に中身を見てみたのは初めてでした。

式をみてみればレイリー減衰が質量と剛性に対する比例減衰である事、なんちゃって周波数軸で減衰を与えたい事、実は調整が難しそうだという事がよくわかりました。

本記事ではレイリー減衰の式の紹介、これまで覚えたシミュレーションと組み合わせて効果の確認をしてみました。

なんとなく、今後仕事で汎用ソフト等で設定をする時に前よりは適切に設定できそうだと思うようになりました。ただ、もともと減衰というのは奥(闇?)が深いので、できれば減衰パラメータに悩まされない人生をおくりたいものです。

レイリー減衰という特殊な減衰を調べてみました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

  • このエントリーをはてなブックマークに追加

SNSでもご購読できます。

コメント

コメントを残す

*