Pythonで多自由度振動系の過渡応答計算をする時はSciPyを使おう

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

このブログではしばらく自前のルンゲ・クッタ法を使って振動解析を行ってきました。しかしPythonで常微分方程式の数値計算をする時はSciPyのodeintを使った方が良いです。ここでは巨人の肩に乗るつもりでSciPy/odeintで振動問題を解くための例題を紹介します。

こんにちは。wat(@watlablog)です。ここではようやく使い始めたSciPyのodeintで振動問題の過渡応答解析をするコード例を紹介します

scipy.integrate.odeintを使った方が良い理由

圧倒的に高速

Pythonで書いた4次のルンゲ・クッタ法の問題点

当WATLABブログでは2020年2月に以下の記事で多自由度振動系を4次のルンゲ・クッタ法で数値計算するコードを紹介しました。この記事はTwitterへのDMやブログのコメントで「多自由度の振動を解析したい」と要望があって書いたものです。

Advertisements

それより前は1自由度の記事でSciPyodeintを使っていましたが、当時はスキル不足からどうしても多自由度化のコードが書けず、odeintで多自由度の問題を解くのは保留していました。

0からスクラッチで書くルンゲ・クッタ法のコードであればブラックボックスにはならないだろうという理由だけです。

ところが自前のルンゲ・クッタ法コードはとにかく計算が遅いです。Pythonでforループを使っていることもあり、数十秒の計算を時間刻み1e-4[s]程度で解こうものなら結構待ちます。NumbaをインストールしてJITを適用しても気休め程度でした。

さらに自前の手法は方程式がいわゆるスティッフな系硬い方程式)だとすぐに発散してしまいました。これは数値解析を専門にしている人であれば、さらにアルゴリズムを改良して解決していくと思いますが、筆者にはハードルが高すぎました。
ちなみにNumbaもPython高速化に関するライブラリです。こちらを参考にしてください。

odeintはFORTRANライブラリを使用している

SciPyのodeintは公式ドキュメントによるとFORTRANodepackというライブラリを使用しているそうです。数値解析の分野では現役でFORTRANが使われていますが、やはり速度はFORTRANかC言語系(最近ではJuliaも?)が強いです。

Solve a system of ordinary differential equations using lsoda from the FORTRAN library odepack.

scipy.integrate.odeint:https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

車輪の再発明で新しいバグを仕込む心配が少ない

先ほどのルンゲ・クッタ法の記事内コメント欄を確認していただければわかりますが、記事公開後、色々な人にご指摘をいただきコードを修正してきました。この辺はブログをやっていてよかった点ですが、やはり自前で用意したコードはバグが潜むリスクを多く持ちます。

0からアルゴリズムを書くのは勉強としてやってよかったと思いますが、ルンゲ・クッタ法等は既に先人達がパッケージングして誰でも使える状態にしてくれています。そちらの方がより多くのレビューを受けているのでバグの可能性は低くなるでしょう(ないとは言えませんが)。

今回は久しぶりにSciPyのodeintを使ってみたら、多自由度でも特に悩まず使えたのでその方法を紹介します。

scipy.integrate.odeintで振動解析をするPythonコード例

動作環境

この記事は以下のPython環境でコードの動作を確認しました。

Python Python 3.9.6
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.21.1
Scipy 1.4.1
Pandas 1.0.5
matplotlib 3.4.3
pillow 7.1.2

1自由度非減衰自由振動

odeintを使った振動解析コードですが、せっかくなので1自由度非減衰系の自由振動から順を追って紹介します。

1自由度非減衰系の自由振動は下図のモデルで表現されます。

1自由度非減衰系振動モデル

続いて運動方程式を式(1)に示します。

\[ m\ddot{x}+kx=0 \tag{1} \]

式(1)をodeintで計算するコードを以下に示します。
メイン部分(コードの下部分)で以下の流れとなっています。

  • モデルを定義する。
  • 解析時間を設定する。
  • 初期条件を設定する。
    初期条件は1D配列にする必要があるためvarとまとめています。
  • odeintを呼び出して計算する。
    関数fを呼び出し、varとt、その後ろに変数をまとめてargs=()で設定します。
    呼び出された関数f内部で初期条件で指定した1D配列変数を分離させ、運動方程式に組み込みます。
    結果をreturnに設定しますが、returnは1D配列になるようにします。
    この1D配列にするというのがodeintを使う時のキーポイントになります。
  • 結果のフォーマットを揃える。
    計算されたvarは列方向が変数、行方向が時間であるため、転置してから指標で抽出することで時間波形を取得できます。
  • グラフに描画する。

運動方程式は「Pythonで1自由度非減衰系の自由振動シミュレーション」で説明したように変数変換を使って表現しており、今後もこのスタイルで書いていきます。
このコードを実行すると以下の結果を得ます。非減衰の自由振動なのでずっと振動する結果となります。

1自由度非減衰系の自由振動結果

1自由度減衰自由振動

下図のようにダッシュポットを追加し、減衰をモデル化すると運動方程式は式(2)となります。

1自由度減衰系振動モデル
\[ m\ddot{x}+c\dot{x}+kx=0 \tag{2} \]

コードを以下に示します。パラメータの設定と運動方程式部分を修正しただけです。

減衰振動が確認できました。減衰係数cの値次第で過減衰、臨界減衰、減衰振動の3パターンの振る舞いをしますが、詳細は「Pythonで1自由度減衰系の自由振動シミュレーション」の記事をご確認ください。

1自由度減衰振動の結果例

1自由度減衰強制振動(サインスイープ加振)

次は強制振動問題を扱ってみます。SciPyのodeintで強制振動を扱う時は、odeintで呼びだす関数fの中で外力を定義する必要があります。
強制振動の運動方程式は式(3)です。外力が追加されました。

\[ m\ddot{x}+c\dot{x}+kx=f \tag{3} \]

ここではサインスイープ信号チャープ信号)を扱ってみましょう。サインスイープは式(4)で与えます。\(f_{0}\)と\(f_{1}\)の順番を変更することでアップスイープやダウンスイープを切り替えることが可能です。

\[ \alpha = \frac{f_{1} - f_{0}}{t_{max}}\\ f = A \sin \left \{ 2 \pi (f_{0} t + \frac{\alpha}{2}t^{2}) \right \} \tag{4} \]

コードを以下に示します。せっかくサインスイープで外力を与えるので、周波数応答関数を計算してみます。周波数応答関数については「Pythonで伝達関数!2つの信号の周波数応答関数を計算する方法」の記事を参照してください。

結果はこちら。変位の波形で途中から振幅が大きくなっているのはサインスイープの周波数がちょうど共振にあたっているからです。周波数応答関数は約5Hzにピークを保ちますが、この系の固有振動数が \(f_{n}=\frac{1}{2\pi}\sqrt{\frac{k}{m}}=5.03\)[Hz]から説明できます。

1自由度強制加振(サインスイープ)

1自由度減衰強制振動(インパルス加振)

外力をインパルス加振に変更した場合のコードを以下に示します。odeintで外力を指定する時は時間の関数にする必要がありますが、プログラミングであれば時間\(t\)をif文の判定に使ってインパルス形状を作成することも可能です(一般的なやり方は他にあるかも?)。

こちらが結果です。外力はインパルス加振になりましたが、モデルのパラメータは変更していないので固有振動数は一緒です。

1自由度強制加振(インパルス加振)

2自由度減衰強制振動(サインスイープ加振)

ここからが本題。odeintで多自由度を表現する方法を紹介します。まずは2自由度から。
今回扱うモデルは以下の2自由度振動系です。

2自由度減衰振動系

運動方程式の形は式(5)で示す通り1自由度系と同じですが、質量・剛性・減衰は行列形式、変位・速度・加速度・外力はベクトル形式になっています。

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

行列とベクトルの中身も書いた運動方程式は式(6)となります。この運動方程式の導出方法は「Pythonで計算するために多自由度振動系を行列形式にする方法」や「csvの値から直列ばねマスモデルを自動生成するPythonコード」をご確認ください。

\[ \begin{bmatrix} m_{1} &0 \\ 0 &m_{2} \end{bmatrix} \begin{bmatrix} \ddot{x}_{1}\\ \ddot{x}_{2} \end{bmatrix} + \begin{bmatrix} c_{1}+c_{2} &-c_{2} \\ -c_{2} &c_{2} \end{bmatrix} \begin{bmatrix} \dot{x}_{1}\\ \dot{x}_{2} \end{bmatrix} + \begin{bmatrix} k_{1}+k_{2} &-k_{2} \\ -k_{2} &k_{2} \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{2} \end{bmatrix} = \begin{bmatrix} f_{1}\\ f_{2} \end{bmatrix} \tag{6} \]

以下に例題のコードを示します。ここではxやFの値を[-1]で抽出しています(最後の質点を加振し、最後の質点の結果を表示)。

多自由度を扱う時の注意点ですが、初期値や結果を1D配列で扱うことが挙げられます(ブログ開始間もない当時の筆者はこれに気づかなかった…)。初期値varや結果outputをnp.concatnate()で結合した形でodeintに渡しています。

計算結果としてのvarはmatplotlibで表示させるために、前処理としてそれぞれの変数(変位と速度)に分離しています。
変位と速度は同じ数だけあるので、半分で分離すれば必ず分けることが可能です。

こちらが結果です。周波数応答関数の裾野がぐにゃぐにゃしているのは応答信号が減衰しきっていないためのリーケージエラーです。この応答曲線はサインスイープ加振だったとしても共振を過ぎると振幅が減衰しているので、本来は指数ウィンドウを使った上で周波数応答関数を処理するのですが、今回は面倒なので割愛(←)。

指数ウィンドウは過去にPythonで作ってみたことがあります。ご参考までに。
Pythonで窓関数が無い場合は?指数窓を自作してみる!

2自由度振動系の結果(サインスイープ)

N自由度減衰強制振動(片端固定)

次はN自由度振動系を表現してみましょう。まずは質点群片方の端部を固定した場合です。
N自由度になっても運動方程式自体は式(5)で表現されています。

N自由度片端固定振動モデル

行列やベクトルの中身を可視化しながらN自由度に拡張すると、片端固定の場合の運動方程式は式(7)となります。

\[ \begin{bmatrix} m_{1} &0 &0&\cdots&\cdots&0 \\ 0 &m_{2}&0&\cdots &\cdots&0\\ 0&0&\ddots&&&\vdots \\ \vdots&\vdots&&\ddots&&\vdots \\ \vdots &\vdots &&&m_{n-1}&0\\ 0&0&\cdots &\cdots &0&m_{n} \end{bmatrix} \begin{bmatrix} \ddot{x}_{1}\\ \ddot{x}_{2} \\ \vdots \\ \vdots \\ \ddot{x}_{n-1} \\ \ddot{x}_{n} \end{bmatrix} + \begin{bmatrix} c_{1}+c_{2} &-c_{2} &0&\cdots&\cdots&0 \\ -c_{2} &c_{2}+c_{3}&-c_{3}&0&\cdots &0\\ 0&-c_{3}&\ddots&&&\vdots \\ 0&0&&\ddots&&\vdots\\ \vdots &\vdots &&&c_{n-1}+c_n&-c_{n}\\ 0&0&\cdots&\cdots &-c_{n} &c_{n} \end{bmatrix} \begin{bmatrix} \dot{x}_{1}\\ \dot{x}_{2} \\ \vdots \\ \vdots \\ \dot{x}_{n-1} \\ \dot{x}_{n} \end{bmatrix} + \begin{bmatrix} k_{1}+k_{2} &-k_{2} &0&\cdots&\cdots&0 \\ -k_{2} &k_{2}+k_{3}&-k_{3}&0&\cdots &0\\ 0&-k_{3}&\ddots&&&\vdots \\ 0&0&&\ddots&&\vdots\\ \vdots &\vdots &&&k_{n-1}+k_n&-k_{n}\\ 0&0&\cdots&\cdots &-k_{n} &k_{n} \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{2} \\ \vdots \\ \vdots \\ x_{n-1} \\ x_{n} \end{bmatrix} = \begin{bmatrix} f_{1}\\ f_{2}\\ \vdots \\ \vdots \\ f_{n-1}\\ f_{n}\\ \end{bmatrix} \tag{7} \]

N自由度の場合、コードの中に長い行列を書くのが非常に煩わしいと感じます。
ここでは表計算ソフトで質量・剛性・減衰をcsvファイルとして作成するパターンで書いた以下のコードを使ってシミュレーションを行います(もちろん2自由度の時のようにコードに直接書いても問題ありません)。

csvファイルをダウンロードできるようにしておきますので、検証でお使いください。

適当に自由度を設定し、例題のファイルは19自由度となっていました。コードを実行すると以下の結果が得られます。ピークがたくさんありますが、最初のピークが最も大きい結果でした。

N自由度強制振動の例

N自由度の質点挙動をアニメーションで確認する方法

これまでは任意質点1つだけの時間波形や周波数応答関数をプロットしていましたが、それとは別に全体の挙動を確認するコードも書いてみました。

plot_animation()という関数を作成し、interval(間引き数)とend_time(終了時間:シミュレーション上でアニメーションを打ち切る時間)を設定しました。

シミュレーション結果xやvは始めから全ての質点に関する全時間データが入っていますので、それをforループで呼び出してプロットしているだけです。

このコードを実行すると、まず始めに先ほどの時間波形と周波数応答関数が表示され、それを閉じるとアニメーション作成が開始されます。指定した時間でループが終了するとcantilever.gifというGIF動画が作成されます。

片端固定のN自由度振動系解析結果(挙動)

びよんびよんしています。

N自由度減衰強制振動(両端固定)

中央加振

次は先ほどのN自由度振動系において、両端固定の場合を検討します。

N自由度振動系(両端固定)

片端固定の運動方程式が書けていれば両端固定の運動方程式を導出するのは簡単です。
下図のように固定したい自由度の行と列を削除することで固定(拘束)を表現できます。

そのため両端固定N自由度振動系の運動方程式は式(7)を少し変形させて式(8)となります。

\[ \begin{bmatrix} m_{1} &0 &0&\cdots&\cdots&0 \\ 0 &m_{2}&0&\cdots &\cdots&0\\ 0&0&\ddots&&&\vdots \\ \vdots&\vdots&&\ddots&&\vdots \\ \vdots &\vdots &&&\ddots&0\\ 0&0&\cdots &\cdots &0&m_{n-1} \end{bmatrix} \begin{bmatrix} \ddot{x}_{1}\\ \ddot{x}_{2} \\ \vdots \\ \vdots \\ \vdots \\ \ddot{x}_{n-1} \end{bmatrix} + \begin{bmatrix} c_{1}+c_{2} &-c_{2} &0&\cdots&\cdots&0 \\ -c_{2} &c_{2}+c_{3}&-c_{3}&0&\cdots &0\\ 0&-c_{3}&\ddots&&&\vdots \\ 0&0&&\ddots&&\vdots\\ \vdots &\vdots &&&\ddots&-c_{n-1}\\ 0&0&\cdots&\cdots &-c_{n-1} &c_{n-1}+c_{n} \end{bmatrix} \begin{bmatrix} \dot{x}_{1}\\ \dot{x}_{2} \\ \vdots \\ \vdots \\ \vdots \\ \dot{x}_{n-1} \end{bmatrix} + \begin{bmatrix} k_{1}+k_{2} &-k_{2} &0&\cdots&\cdots&0 \\ -k_{2} &k_{2}+k_{3}&-k_{3}&0&\cdots &0\\ 0&-k_{3}&\ddots&&&\vdots \\ 0&0&&\ddots&&\vdots\\ \vdots &\vdots &&&\ddots&-k_{n-1}\\ 0&0&\cdots&\cdots &-k_{n-1} &k_{n-1}+k_{n} \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{2} \\ \vdots \\ \vdots \\ \vdots \\ x_{n-1} \end{bmatrix} = \begin{bmatrix} f_{1} \\ f_{2} \\ \vdots \\ \vdots \\ \vdots \\ f_{n-1}\\ \end{bmatrix} \tag{8} \]

コード上ではgenerate_matrix()関数の中で、np.delete()K_all, C_all, M_allに対して追加して式変形を表現しました。また両端支持の振動系なのでインパルス加振を与える位置はf()関数の中でforce_posを設定し、ちょうど中間の質点に印加するようにしています。

時間波形と周波数応答関数を示します。先ほどとちょっとピークの形が変わりましたね。

両端支持N自由度振動系の時間波形と周波数応答関数(中央加振)

コードを実行すると動画はdouble-supported-beam.gifが作成されます。

動画で可視化してみると一目瞭然です。運動方程式を少し変えただけで両端の変位が小さくなりました。当たり前の結果ですが面白いですね。

両端支持N自由度振動系中央加振の挙動

びよんびよんしています!

端部加振

f()関数内のforce_posを以下の内容に修正して、加振点を左端にしてみましょう。

動画のみ紹介しますが、左端から波が伝わっていく様子がわかりました。

両端支持N自由度振動系端部加振の挙動

端部サインスイープ加振

せっかくなのでサインスイープ加振も試してみましょう。既に紹介したコードを組み合わせるだけですが、コピペで動作できるようにした方が便利かもしれないので全コードを以下に示します。

縄跳びのシミュレーションのようになりました。

両端支持N自由度振動系端部サインスイープ加振の挙動

N自由度Duffing方程式の強制振動アニメーション

最後はN自由度Duffing方程式を可視化してみましょう。
Duffing方程式(式(9))は「PythonでDuffing振動子を解析してアトラクターを見る」の記事でアトラクターを確認しましたが、今回はN自由度ということで以下の式(10)と表現しました。

\[ \ddot{x} + \delta \dot{x} + \alpha x + \beta x^{3} = \gamma \cos(\omega t + \phi) \tag{9} \]
\[ \mathbf{M}\mathbf{\ddot{x}} + \mathbf{C}\mathbf{\dot{x}} + \mathbf{K_{\alpha}}\mathbf{x} + \mathbf{K_{\beta}}\mathbf{x}^{3} =\mathbf{F} \tag{10} \]

コード例を以下に示します。式(10)の3乗項にかかる係数は簡単に既存剛性行列を1000倍するだけにしました。