csvの値から直列ばねマスモデルを自動生成するPythonコード

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

振動シミュレーションをプログラミングで行う場合、質量行列や剛性行列の記述が面倒です。しかし、ばね要素の結合と行列の重畳を理解する事で規則性を利用できます。ここではcsvに記録した質量や剛性の値を使ってPython/Numpy形式の行列生成を行う方法を紹介します。

こんにちは。wat(@watlablog)です。ここではPythonを使って振動モデルを自動生成する方法を紹介します

なぜ振動モデルを自動生成したいのか?

本記事はほぼ筆者のメモですが、少しでも有益なページになると良いと思い、ここでコードのモチベーションを説明しておこうと思います。

スクリプトに直打ちするのはかなり面倒!

直列振動モデルの概要

このページで扱う直列振動モデルというのは、下図のようなばねマスが直列に繋がった系のことを指します。この図は2自由度モデルですが、これは振動を扱う人であれば見慣れた図だと思います。

2自由度モデル

Pythonスクリプト表現の例

上のモデルをPythonスクリプトで表現すると下のコードのようになります。

Numpyndarrayで書いていますが、括弧を複数書いたりパラメータ同士を加算したり負値にしたりカンマで区切ったり…とスクリプト上で体裁を整えるのが面倒です(面倒でない?)。

自由度を変更した時に修正するのも面倒!

先ほどは2自由度モデルでしたが、3自由度になるとPythonでは以下のように書きます。

自由度が1つ増えただけですが、行を追加するだけではなく各行全てを編集する必要があります。

自由度が増えたり減ったりするケースを検討する場合、いちいちスクリプト上でカンマを打ったり0を追加したりといった作業は意外と面倒です。

多自由度モデルを書くのも面倒!

10自由度以下はまだなんとか頑張ろうという気にもなりますが、例えば100自由度のモデルを検討する必要が出てきた場合、これは最初の記述をするだけで面倒を通り越して無理だと思います(下図はn自由度直列振動モデル)。

n自由度モデル

絶対どこかで書き間違えをしてバグの原因にもなります。ここまでの自由度を検討する場合は自動化がマストだと思います。

また、このように多自由度の振動モデルを自動生成できるようになれば、振動シミュレーションに関しても最適化計算を適用したりしやすくなります。この記事はその方面に進むための足がかりとします。

自動生成の要件

Advertisements

それでは自動生成コードを書くために、どのようにやりたいことを実現するかの要件を書いておきます。

端部固定の直列振動系に限定する

今回は前提として直列のばねマス振動系を自動生成する方法を検討します。直列や並列、節点を飛び越したばねの接続は簡単化のため考えません。

複雑なモデルの場合はModelicaを使った1DCAEツールによる物理シミュレーションの方が良いと考えます。

下図のように一方の端部を固定したモデルを生成します。

端部固定のn自由度系

csvファイルからパラメータを読み取る

csvファイルにパラメータを任意の自由度分記録しておいて、それを読み取る形式とします。以下に2自由度の場合のcsvファイル例を示します。

列方向に質量、剛性列を作成し、行方向に自由度分数値を並べます。
ここで単位は質量が[kg], 剛性が[N/m]とします。単位系は任意に変えても構いませんが、後々シミュレーション結果で固有振動数が[Hz]になるようにと考えてこの単位系にしています。

csvファイルの例

この時、端部固定を前提としているので、ファイルは以下の図に示すパラメータ位置関係となります(mとkの個数が同一になる)。

計算済の剛性行列を作る

最後は読み込んだ数値を基に、計算済の質量行列と剛性行列を作ります。
上記モデルの質量行列は対角行列となるので特に工夫は要りませんので、剛性行列の作り方に着目します。

剛性行列はばね要素の重畳と固定節点の性質から自動的に自由度分作成します。

ばね要素の結合と剛性行列の重畳について

プログラムで規則性を記述するために、ばねを直列に接続していった場合の剛性行列の変化を確認します。

まずは下図に示す単一のばね要素ですが、節点質量を0として節点周りの力関係から運動方程式を導くと式(1)となります。

\[ \begin{cases} f_{1} = k_{1}x_{1} - k_{1}x_{2}\\ f_{2} = -k_{1}x_{1} + k_{1}x_{2} \end{cases} \tag{1} \]

上式は行列形式で書くと式(2)となります。変位ベクトルの前にかかっている行列が剛性行列です。

\[ \begin{Bmatrix} f_{1} \\ f_{2} \end{Bmatrix} = \begin{bmatrix} k_{1} & -k_{1} \\ -k_{1} & k_{1} \end{bmatrix} \begin{Bmatrix} x_{1} \\ x_{2} \end{Bmatrix} \tag{2} \]

運動方程式の立て方や行列形式の詳細は「Pythonで計算するために多自由度振動系を行列形式にする方法」にまとめていますのでこちらの記事を参照してください。

次に、下図のようにばね要素を1つ追加してみます。

本モデルの運動方程式は式(3)です。外力ベクトルや変位ベクトルは単純に要素が増えただけですが、剛性行列は一部加算されている項目があります。

\[ \begin{Bmatrix} f_{1} \\ f_{2} \\ f_{3} \end{Bmatrix} = \begin{bmatrix} k_{1} & -k_{1} & 0 \\ -k_{1} & k_{1} + k_{2} & -k_{2} \\ 0 & -k_{2} & k_{2} \end{bmatrix} \begin{Bmatrix} x_{1} \\ x_{2} \\ x_{3} \end{Bmatrix} \tag{3} \]

この加算のイメージを以下に図解します。単一のばね要素は同じ形の剛性行列を持ちますが、結合点の位置(ここでは\(x_{2}\)の節点と\(x_{3}\)の節点)で剛性行列を重畳させたことを意味しています。

ばね要素の重畳

隣接節点を飛び越えた位置同士でばね要素を結合する場合は、結合させた節点位置に重畳させますが、直列モデルの場合はこの図を覚えておけば運動方程式を簡単に求めることが可能です。

固定した場合の剛性行列の変化について

先ほどまでのモデルは拘束をかけない(つまりフリーフリー)状態でしたが、今回は片方の端部を固定します。

拘束をかけると剛性行列は下図のように変化します。

固定した場合の剛性行列の変化

ばねが1つしかない場合はわかりやすく、どちらか一方を固定することで単一の剛性値となります。

これは固定することでその節点の動きがなくなるため、行列の該当要素を削除することを意味します。2自由度端部固定の場合の剛性行列は覚えている方も多いのではないでしょうか。

ここまででばね要素を結合していった場合と固定させた場合の剛性行列の変化を見てきました。

規則性がわかったので、n自由度の場合に拡張してコーディングしてみましょう。

csvファイルから質量・剛性行列を生成するコード

全コード

コードはとても簡単です。以下に全コードを示します。
まず、generate_matrix()関数の引数にファイル名を与え、pandasで読み込みます。

使ったファイルはこちら(ダウンロード可)。

読み込んだ値をそれぞれm, kとSeriesで抽出し、forループで剛性行列を作成します。forループではnp.zerosにより0で埋めた行列K_allに、単一のばね要素として作成したKを重畳させます。

重畳させる時にNumpyのスライスで位置をずらしながら、直列のばね要素の結合部で加算しています。

このK_allに対し、最後に固定を反映させるためnp.delete(操作する行列, 削除する指標, 方向(0:行, 1:列))を使います。

質量行列は1d配列を使って対角行列を作成すれば良いので、np.diag()で作成可能です。

実行結果

こちらが実行結果です。2自由度は既に上の図解中の式で検算可能ですが、正しく作れているようです。

100自由度でも確認してみる

次のファイルは100自由度の場合です。

以下が結果。どうやらうまくいっている様子。

振動シミュレーションで使えるか確認する

参考記事と使用するcsvファイル

本来目的である振動のシミュレーションで使う事ができるかどうか確認します。
今回は周波数応答解析や固有値解析による検証を行います。

周波数応答解析については「振動モデルを直接法の周波数応答解析で解くPythonコード」、固有値解析については「Pythonで多自由度系の固有値解析!固有振動数とモードを計算」の記事に詳細を紹介していますので、是非読んでみてください。

先ほどの100自由度のファイルだと固有値がかなり隣接していて効果が確認し難かったため、今回は以下の10自由度のファイルを使います。

全コード

以下が全コードです。詳細は先ほどの記事に記載した通りですが、今回は自動生成した質量と剛性行列がこれらのシミュレーションで使えるかどうかを確認します。

固有値解析結果

コンソールに出力される固有振動数の結果はこちらです。1.0[Hz]以上50[Hz]未満の範囲で分布しているようです。

周波数応答解析結果

周波数応答解析の結果はmatplotlibによるプロットで表示されます。先ほどの固有振動数の位置でピークが立っていることが確認できました。

周波数応答解析結果

以上から、今回自動生成した質量行列と剛性行列はPythonによる振動シミュレーションで十分使用できることが確認されました。

まとめ

今回はこれまでPythonコード内に直接記述していた多自由度振動系のパラメータをcsvから読み込んで行列生成してみました。

これでn自由度のパラメータをいちいちコード上に記載しなくてもよくなりそうです。

Numpyを使えば高速に行列生成ができるようです!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

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

SNSでもご購読できます。

コメントを残す

*