Python NumPy SciPy : 不減衰系の固有振動数と固有振動モード

今回は Python で二由度不減衰系の固有振動数と固有モードを求めてみます。 二由度系を例としていますが、より多自由度の N 自由度になっても計算方法は同じです。 以下の二つのバネと二つの質点が直列に結合している不減衰系を考えます。

1506142DOF.jpg

各質点の力のつりあい式から、以下の行列形式の運動方程式が得られます。

\begin{equation} \begin{bmatrix} {m_1} & {0}\\ {0} & {m_2} \end{bmatrix} \begin{bmatrix} {\ddot{x}_1}\\ {\ddot{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} {0}\\ {0} \end{bmatrix} \nonumber \end{equation}

これを一般化した形で表せば

\begin{equation} \mathbf{M}\ddot{\boldsymbol{x}}+\mathbf{K}\boldsymbol{x}=\boldsymbol{0} \nonumber \end{equation}

となります。ここで \(\mathbf{M}\) を質量行列、\(\mathbf{K}\) を剛性行列と言います。

より一般化するために対角項に変位 \(x_i\) を並べた対角行列 \(\mathbf{x}\) を用いると、運動方程式は

\begin{equation} \begin{bmatrix} {m_1} & {0}\\ {0} & {m_2} \end{bmatrix} \begin{bmatrix} {\ddot{x}_1} & {0}\\ {0} & {\ddot{x}_2} \end{bmatrix} + \begin{bmatrix} {k_1+k_2} & {-k_2}\\ {-k_2} & {k_2} \end{bmatrix} \begin{bmatrix} {x_1} & {0}\\ {0} & {x_2} \end{bmatrix} = \begin{bmatrix} {0}\\ {0} \end{bmatrix} \nonumber \end{equation} \begin{equation} \mathbf{M}\ddot{\mathbf{x}}+\mathbf{K}\mathbf{x}=\boldsymbol{0} \nonumber \end{equation}

と表せます。一般解を複素関数の形で仮定し

\begin{equation} \mathbf{x}=\mathbf{\Phi}\mathbf{e} \nonumber \end{equation}

これを運動方程式に代入すると以下のようになります。 ここで \(\mathbf{\Phi}\) は各変位の r 次の振幅 \(\phi_i\) を列に並べた行列、\(\mathbf{e}\) は \(e^{j\omega_r t}\) を対角項とする対角行列です。

\begin{equation} \mathbf{K}\mathbf{\Phi}=\mathbf{\Omega}^2\mathbf{M}\boldsymbol{\Phi} \nonumber \end{equation}

なお \(\mathbf{\Omega}\) は r 次の振動数 \(\omega_r\) を対角項とする対角行列です。 これは

\begin{equation} \mathbf{K}\mathbf{V}=\mathbf{D}\mathbf{M}\mathbf{V} \nonumber \end{equation}

の形の固有値問題となっており、SciPy の eig 関数で計算することができます。 ここで \(\mathbf{D}\) は固有値を対角項とする対角行列、\(\mathbf{V}\) はそれに対応する右固有ベクトルを列に持つ行列 です。

SciPy の eig 関数は Matlab と同様にこの形の固有値計算ができますが、NumPy や Octave, Scilab の eig 関数は標準型の固有値問題しか解けないので、逆行列計算等が別途必要になります。

以上から、Python では二由度不減衰系の固有振動数と固有振動モードを以下のように求めることができます。

図はわかりやすくするために水平方向である振幅を縦方向にして表示しています。 1 次モードでは二つの質点が同位相で振動し、2 次モードでは逆位相で振動していることがわかります。

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

# パラメータ設定
m1 = 1
m2 = 2
k1 = 300
k2 = 500

M = np.array([[m1, 0], [0, m2]])
K = np.array([[k1+k2, -k2], [-k2, k2]])

w, vr = linalg.eig(K, M)
omega = np.sort(np.sqrt(w))
index = np.argsort(np.sqrt(w))
f = omega/(2*np.pi)
mode = np.array([[1,1], [vr[1,index[0]]/vr[0,index[0]], vr[1,index[1]]/vr[0,index[1]]]])

print("frequency[Hz] =", f)
print("mode =\n", mode)

# プロット
plt.figure()
plt.subplot(211)
plt.plot([0, 1, 2], np.insert(mode[0], 0, 0), "bo-")
plt.hlines(0, 0, 2, linestyles="-.")
plt.xlabel("First order")
plt.ylim(-1.5, 1.5)
plt.subplot(212)
plt.plot([0, 1, 2], np.insert(mode[1], 0, 0), "bo-")
plt.hlines(0, 0, 2, linestyles="-.")
plt.xlabel("Second order")
plt.ylim(-1.5, 1.5)
# plt.show()
frequency[Hz] = [ 1.39737839+0.j  4.96428689+0.j]
mode =
 [[ 1.          1.        ]
 [ 1.44582364 -0.34582364]]

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
電子書籍
線形代数(1): Pythonによる科学技術計算 実践編