回転数変化を模擬した周波数応答解析で次数分析するPythonコード

モータやエンジンといった回転機は構造の共振に当たらないように加振力を設計する必要があります。ここではPythonで回転数変化を模擬した周波数応答解析を行い、実験の次数トラッキング解析と比較できるシミュレーション結果を得る事を目標とします。

こんにちは。wat(@watlablog)です。ここでは覚えたての周波数応答解析で回転数変化を模擬した加振力の計算してみます！

※本記事は過去に学んだ内容の応用、小Tipsです！

目次（項目クリックでジャンプできます）

本記事の対象者

加振力の変化を考慮した周波数応答解析をしたい人

前回の記事「振動モデルを直接法の周波数応答解析で解くPythonコード」ではPythonを使って振動モデルの周波数応答解析を行い、各周波数における振動応答のレベルを計算できるようになりました。

加振力の変化を考慮した周波数応答解析のPythonコード

動作環境

このページに記載のコードは以下のPC環境で検証しました。参考までに。

PC

Windows	OS	Windows10 64bit
	CPU	2.4[GHz]
	メモリ	4[GB]

Mac	OS	macOS Catalina 10.15.7
	CPU	1.4[GHz]
	メモリ	8[GB]

Python環境

Python環境は以下です。外部ライブラリはNumpyと結果可視化用にmatplotlibを使っています。

Python	Python 3.7.7
PyCharm (IDE)	PyCharm CE 2020.1
Numpy	1.19.0
matplotlib	3.2.2

モデル定義

今回も理解する事を優先させて、2自由度振動系を定義します。

# シミュレーションモデル
def model():
    # 質量・剛性の集中定数を設定する
    m1 = 0.05
    m2 = 0.01
    k1 = 500000
    k2 = 200000

    # 質量マトリクスと剛性マトリクスを作成
    M = np.array([[m1, 0],
                  [0, m2]])
    K = np.array([[k1 + k2, -k2],
                  [-k2, k2]])
    return M, K

# シミュレーションモデル

def model():

# 質量・剛性の集中定数を設定する

m1 = 0.05

m2 = 0.01

k1 = 500000

k2 = 200000

# 質量マトリクスと剛性マトリクスを作成

M = np.array([[m1, 0],

[0, m2]])

K = np.array([[k1 + k2, -k2],

[-k2, k2]])

return M, K

変化する加振力を作成

変化する加振力セットは以下の関数で作成します。dataとして時間波形を作成した後にフーリエ変換して周波数波形にするという作り方をしています。周波数応答解析時にそのまま使えるように周波数分解能やサンプリングレートをサンプリングの定理に則って設定しています。

# 擬似加振力スイープ周波数波形の作成
def freq_wave(freq_l, freq_u, step, f_base, f_max):
    samplerate = freq_u * 2.56                              # サンプリングレート[Hz]は上限周波数の2.56倍
    dt = 1 / samplerate                                     # 時間刻み[s]
    end = 1 / step                                          # 周波数分解能を周波数応答解析のステップにあわせる[s]
    t = np.arange(0, end + dt, dt)                          # 時間軸[s]

    n = int((freq_u - freq_l) / step)                       # 計算ステップ数
    f_increment = (f_max - f_base) / n                      # ステップ当たりの回転数[rps]変化

    amp_list = []
    rpms = []
    for i in range(n):
        freq_i = f_base + (i + 1) * f_increment
        rpms.append(freq_i * 60)                            # 基本周波数をRPMに変換

        #data = signal.sawtooth(2 * np.pi * freq_i * t)     # のこぎり波
        data = np.sin(2 * np.pi * freq_i * t)               # 正弦波

        amp, freq = calc_fft(data, samplerate)              # フーリエ変換
        amp_sub = amp[int(freq_l / step):int(freq_u)]       # 周波数応答解析で使用する範囲のみ振幅を抽出
        amp_list.append(amp_sub)

    freq_sub = freq[int(freq_l / step):int((freq_u + step) / step)]    # 周波数応答解析で使用する範囲のみ周波数軸を抽出

    return amp_list, freq_sub, rpms

# 擬似加振力スイープ周波数波形の作成

def freq_wave(freq_l, freq_u, step, f_base, f_max):

samplerate = freq_u * 2.56 # サンプリングレート[Hz]は上限周波数の2.56倍

dt = 1 / samplerate # 時間刻み[s]

end = 1 / step # 周波数分解能を周波数応答解析のステップにあわせる[s]

t = np.arange(0, end + dt, dt) # 時間軸[s]

n = int((freq_u - freq_l) / step) # 計算ステップ数

f_increment = (f_max - f_base) / n # ステップ当たりの回転数[rps]変化

amp_list = []

rpms = []

for i in range(n):

freq_i = f_base + (i + 1) * f_increment

rpms.append(freq_i * 60) # 基本周波数をRPMに変換

#data = signal.sawtooth(2 * np.pi * freq_i * t) # のこぎり波

data = np.sin(2 * np.pi * freq_i * t) # 正弦波

amp, freq = calc_fft(data, samplerate) # フーリエ変換

amp_sub = amp[int(freq_l / step):int(freq_u)] # 周波数応答解析で使用する範囲のみ振幅を抽出

amp_list.append(amp_sub)

freq_sub = freq[int(freq_l / step):int((freq_u + step) / step)] # 周波数応答解析で使用する範囲のみ周波数軸を抽出

return amp_list, freq_sub, rpms

補足：実用的な加振力データを作るには

今回はサンプルとしてPython上で理想的な波形を生成していますが、実際の設計では製品にどのような外力が加わるのかをよく分析して設定する事が重要です。

例えばモータであれば、モータ内部のステータやロータに働く電磁力を別途磁界解析を行なって得るといった事が実用的には考えられます。

他にも機構部品がぶつかったり摩擦によって力が発生したり…解析の目的に応じて外力の詳細度や種類を検討する必要があります。

周波数応答解析の関数を修正

「振動モデルを直接法の周波数応答解析で解くPythonコード」は関数の内部で加振波形を作っていましたが、今回は関数の引数excitationとして波形を渡す形で書いてみました。

# 周波数応答解析を計算する関数
def frequency_response_analysis(M, K, C, steps, exc_node, res_node, excitation):
    steps = steps * 2 * np.pi                                           # 周波数を角振動数に変換
    n_dofs = len(M)                                                     # 自由度の総数
    F = np.zeros((len(steps), n_dofs))                                  # 周波数毎の外力ベクトルを初期化

    # 角振動数分解能毎に計算するループ
    X = []
    for i in range(len(steps)-1):
        F[i][exc_node] = excitation[i]                                  # 外力ベクトルの更新
        matrix = (-steps[i] ** 2) * M + (1j * steps[i] * C) + K         # 逆行列にする前の行列を計算しておく
        inv_matrix = np.linalg.inv(matrix)                              # 逆行列計算
        xn = np.dot(inv_matrix, F[i])[res_node]                         # 方程式を解き、複素変位ベクトルを求める
        xn = np.sqrt(xn.real ** 2 + xn.imag ** 2)                       # 振幅成分を計算する
        X.append(xn)
    steps = steps / (2 * np.pi)
    return steps, X

# 周波数応答解析を計算する関数

def frequency_response_analysis(M, K, C, steps, exc_node, res_node, excitation):

steps = steps * 2 * np.pi # 周波数を角振動数に変換

n_dofs = len(M) # 自由度の総数

F = np.zeros((len(steps), n_dofs)) # 周波数毎の外力ベクトルを初期化

# 角振動数分解能毎に計算するループ

X = []

for i in range(len(steps)-1):

F[i][exc_node] = excitation[i] # 外力ベクトルの更新

matrix = (-steps[i] ** 2) * M + (1j * steps[i] * C) + K # 逆行列にする前の行列を計算しておく

inv_matrix = np.linalg.inv(matrix) # 逆行列計算

xn = np.dot(inv_matrix, F[i])[res_node] # 方程式を解き、複素変位ベクトルを求める

xn = np.sqrt(xn.real ** 2 + xn.imag ** 2) # 振幅成分を計算する

X.append(xn)

steps = steps / (2 * np.pi)

return steps, X

メイン関数の部分で先ほど作成した加振力セットを1行ずつ抽出し、周波数応答解析を実行する事で加振力が変化する分析を行います。

# 周波数応答解析を実行
X_list = []
for i in range(len(amp)):
    freq_axis, X = frequency_response_analysis(M=M,
                                               K=K,
                                               C=C,
                                               steps=freq,
                                               exc_node=1,
                                               res_node=1,
                                               excitation=amp[i])
    X_list.append(X)

# 周波数応答解析を実行

X_list = []

for i in range(len(amp)):

freq_axis, X = frequency_response_analysis(M=M,

K=K,

C=C,

steps=freq,

exc_node=1,

res_node=1,

excitation=amp[i])

X_list.append(X)

全コード(コピペ用)

以下にコピペ用全コードを示します。減衰の考え方等は一緒ですが、フーリエ変換関数を追加したり、前回記事とはちょこちょこ変更していますので、コード内コメントもご確認下さい。

import numpy as np
from scipy import signal
from scipy import fftpack
from matplotlib import pyplot as plt

# シミュレーションモデル
def model():
    # 質量・剛性の集中定数を設定する
    m1 = 0.05
    m2 = 0.01
    k1 = 500000
    k2 = 200000

    # 質量マトリクスと剛性マトリクスを作成
    M = np.array([[m1, 0],
                  [0, m2]])
    K = np.array([[k1 + k2, -k2],
                  [-k2, k2]])
    return M, K

# Rayleigh dampingのパラメータを決定して減衰マトリクスを生成する関数
def generate_rayleigh_damping(f1, f2, zeta1, zeta2, M, K):
    # 指定周波数を角振動数に変換
    w1 = 2 * np.pi * f1
    w2 = 2 * np.pi * f2

    # alphaとbetaを使って減衰行列を計算
    alpha = (2 * w1 * w2 * (zeta2 * w1 - zeta1 * w2)) / (w1 ** 2 - w2 ** 2)
    beta = (2 * (zeta1 * w1 - zeta2 * w2)) / (w1 ** 2 - w2 ** 2)
    C = np.array(alpha * M + beta * K)

    return C

# 周波数応答解析を計算する関数
def frequency_response_analysis(M, K, C, steps, exc_node, res_node, excitation):
    steps = steps * 2 * np.pi                                           # 周波数を角振動数に変換
    n_dofs = len(M)                                                     # 自由度の総数
    F = np.zeros((len(steps), n_dofs))                                  # 周波数毎の外力ベクトルを初期化

    # 角振動数分解能毎に計算するループ
    X = []
    for i in range(len(steps)-1):
        F[i][exc_node] = excitation[i]                                  # 外力ベクトルの更新
        matrix = (-steps[i] ** 2) * M + (1j * steps[i] * C) + K         # 逆行列にする前の行列を計算しておく
        inv_matrix = np.linalg.inv(matrix)                              # 逆行列計算
        xn = np.dot(inv_matrix, F[i])[res_node]                         # 方程式を解き、複素変位ベクトルを求める
        xn = np.sqrt(xn.real ** 2 + xn.imag ** 2)                       # 振幅成分を計算する
        X.append(xn)
    steps = steps / (2 * np.pi)
    return steps, X

# 固有値解析を計算する関数
def eigen(M, K):
    # 質量マトリクスの逆行列を計算
    M_inv = np.linalg.inv(M)

    # 固有値と固有ベクトルを計算
    omega, v = np.linalg.eig(np.dot(M_inv, K))

    # 固有値の順番を昇順にソート
    omega_sort = (1 / (2 * np.pi)) * np.sqrt(np.sort(omega))

    # 固有値のソート時のインデックスを取得
    sort_index = np.argsort(omega)

    # 固有値に対応する固有ベクトルをソート
    v_sort = []
    for i in range(len(sort_index)):
        v_sort.append(v.T[sort_index[i]])
    v_sort = np.array(v_sort)
    return omega_sort, v_sort, sort_index

# フーリエ変換をする関数
def calc_fft(data, samplerate):
    spectrum = fftpack.fft(data)                                     # 信号のフーリエ変換
    amp = np.sqrt((spectrum.real ** 2) + (spectrum.imag ** 2))       # 振幅成分
    amp = amp / (len(data) / 2)                                      # 振幅成分の正規化（辻褄合わせ）
    phase = np.arctan2(spectrum.imag, spectrum.real)                 # 位相を計算
    phase = np.degrees(phase)                                        # 位相をラジアンから度に変換
    freq = np.linspace(0, samplerate, len(data))                     # 周波数軸を作成
    return amp, freq

# 擬似加振力スイープ周波数波形の作成
def freq_wave(freq_l, freq_u, step, f_base, f_max):
    samplerate = freq_u * 2.56                              # サンプリングレート[Hz]は上限周波数の2.56倍
    dt = 1 / samplerate                                     # 時間刻み[s]
    end = 1 / step                                          # 周波数分解能を周波数応答解析のステップにあわせる[s]
    t = np.arange(0, end + dt, dt)                          # 時間軸[s]

    n = int((freq_u - freq_l) / step)                       # 計算ステップ数
    f_increment = (f_max - f_base) / n                      # ステップ当たりの回転数[rps]変化

    amp_list = []
    rpms = []
    for i in range(n):
        freq_i = f_base + (i + 1) * f_increment
        rpms.append(freq_i * 60)                            # 基本周波数をRPMに変換

        #data = signal.sawtooth(2 * np.pi * freq_i * t)     # のこぎり波
        data = np.sin(2 * np.pi * freq_i * t)               # 正弦波

        amp, freq = calc_fft(data, samplerate)              # フーリエ変換
        amp_sub = amp[int(freq_l / step):int(freq_u)]       # 周波数応答解析で使用する範囲のみ振幅を抽出
        amp_list.append(amp_sub)

    freq_sub = freq[int(freq_l / step):int((freq_u + step) / step)]    # 周波数応答解析で使用する範囲のみ周波数軸を抽出

    return amp_list, freq_sub, rpms

# モデルの質量マトリクスと剛性マトリクスを求める関数を実行
M, K = model()

# 減衰マトリクスを計算する関数を実行
C = generate_rayleigh_damping(f1=10, f2=1000, zeta1=0.2, zeta2=0.1, M=M, K=K)

# サンプルの加振力周波数波形を作成
freq_l = 10
freq_u = 1000
step = 5
f_base = 10
f_max = 500
amp, freq, rpms = freq_wave(freq_l=freq_l, freq_u=freq_u, step=step, f_base=f_base, f_max=f_max)

# 周波数応答解析を実行
X_list = []
for i in range(len(amp)):
    freq_axis, X = frequency_response_analysis(M=M,
                                               K=K,
                                               C=C,
                                               steps=freq,
                                               exc_node=1,
                                               res_node=1,
                                               excitation=amp[i])
    X_list.append(X)

# 固有値解析を実行(確認用)
omega_sort, v_sort, index= eigen(M, K)
print('Natural frequency[Hz]=', omega_sort)

# ここからグラフ描画
# グラフをオブジェクト指向で作成する。
fig = plt.figure()
ax1 = fig.add_subplot(111)

# データをプロットする。
im = ax1.pcolormesh(rpms,
                    freq_axis,
                    np.array(X_list).T,
                    cmap='jet',
                    vmin=0, vmax=np.max(X_list))

# カラーバーを設定する。
cbar = fig.colorbar(im)
cbar.set_label('Displacement[m]')

# 軸設定する。
ax1.set_xlabel('Rev.[rpm]')
ax1.set_ylabel('Frequency[Hz]')

# グラフを表示する。
plt.show()
plt.close()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

import numpy as np

from scipy import signal

from scipy import fftpack

from matplotlib import pyplot as plt

# シミュレーションモデル

def model():

# 質量・剛性の集中定数を設定する

m1 = 0.05

m2 = 0.01

k1 = 500000

k2 = 200000

# 質量マトリクスと剛性マトリクスを作成

M = np.array([[m1, 0],

[0, m2]])

K = np.array([[k1 + k2, -k2],

[-k2, k2]])

return M, K

# Rayleigh dampingのパラメータを決定して減衰マトリクスを生成する関数

def generate_rayleigh_damping(f1, f2, zeta1, zeta2, M, K):

# 指定周波数を角振動数に変換

w1 = 2 * np.pi * f1

w2 = 2 * np.pi * f2

# alphaとbetaを使って減衰行列を計算

alpha = (2 * w1 * w2 * (zeta2 * w1 - zeta1 * w2)) / (w1 ** 2 - w2 ** 2)

beta = (2 * (zeta1 * w1 - zeta2 * w2)) / (w1 ** 2 - w2 ** 2)

C = np.array(alpha * M + beta * K)

return C

# 周波数応答解析を計算する関数

def frequency_response_analysis(M, K, C, steps, exc_node, res_node, excitation):

steps = steps * 2 * np.pi # 周波数を角振動数に変換

n_dofs = len(M) # 自由度の総数

F = np.zeros((len(steps), n_dofs)) # 周波数毎の外力ベクトルを初期化

# 角振動数分解能毎に計算するループ

X = []

for i in range(len(steps)-1):

F[i][exc_node] = excitation[i] # 外力ベクトルの更新

matrix = (-steps[i] ** 2) * M + (1j * steps[i] * C) + K # 逆行列にする前の行列を計算しておく

inv_matrix = np.linalg.inv(matrix) # 逆行列計算

xn = np.dot(inv_matrix, F[i])[res_node] # 方程式を解き、複素変位ベクトルを求める

xn = np.sqrt(xn.real ** 2 + xn.imag ** 2) # 振幅成分を計算する

X.append(xn)

steps = steps / (2 * np.pi)

return steps, X

# 固有値解析を計算する関数

def eigen(M, K):

# 質量マトリクスの逆行列を計算

M_inv = np.linalg.inv(M)

# 固有値と固有ベクトルを計算

omega, v = np.linalg.eig(np.dot(M_inv, K))

# 固有値の順番を昇順にソート

omega_sort = (1 / (2 * np.pi)) * np.sqrt(np.sort(omega))

# 固有値のソート時のインデックスを取得

sort_index = np.argsort(omega)

# 固有値に対応する固有ベクトルをソート

v_sort = []

for i in range(len(sort_index)):

v_sort.append(v.T[sort_index[i]])

v_sort = np.array(v_sort)

return omega_sort, v_sort, sort_index

# フーリエ変換をする関数

def calc_fft(data, samplerate):

spectrum = fftpack.fft(data) # 信号のフーリエ変換

amp = np.sqrt((spectrum.real ** 2) + (spectrum.imag ** 2)) # 振幅成分

amp = amp / (len(data) / 2) # 振幅成分の正規化（辻褄合わせ）

phase = np.arctan2(spectrum.imag, spectrum.real) # 位相を計算

phase = np.degrees(phase) # 位相をラジアンから度に変換

freq = np.linspace(0, samplerate, len(data)) # 周波数軸を作成

return amp, freq

# 擬似加振力スイープ周波数波形の作成

def freq_wave(freq_l, freq_u, step, f_base, f_max):

samplerate = freq_u * 2.56 # サンプリングレート[Hz]は上限周波数の2.56倍

dt = 1 / samplerate # 時間刻み[s]

end = 1 / step # 周波数分解能を周波数応答解析のステップにあわせる[s]

t = np.arange(0, end + dt, dt) # 時間軸[s]

n = int((freq_u - freq_l) / step) # 計算ステップ数

f_increment = (f_max - f_base) / n # ステップ当たりの回転数[rps]変化

amp_list = []

rpms = []

for i in range(n):

freq_i = f_base + (i + 1) * f_increment

rpms.append(freq_i * 60) # 基本周波数をRPMに変換

#data = signal.sawtooth(2 * np.pi * freq_i * t) # のこぎり波

data = np.sin(2 * np.pi * freq_i * t) # 正弦波

amp, freq = calc_fft(data, samplerate) # フーリエ変換

amp_sub = amp[int(freq_l / step):int(freq_u)] # 周波数応答解析で使用する範囲のみ振幅を抽出

amp_list.append(amp_sub)

freq_sub = freq[int(freq_l / step):int((freq_u + step) / step)] # 周波数応答解析で使用する範囲のみ周波数軸を抽出

return amp_list, freq_sub, rpms

# モデルの質量マトリクスと剛性マトリクスを求める関数を実行

M, K = model()

# 減衰マトリクスを計算する関数を実行

C = generate_rayleigh_damping(f1=10, f2=1000, zeta1=0.2, zeta2=0.1, M=M, K=K)

# サンプルの加振力周波数波形を作成

freq_l = 10

freq_u = 1000

step = 5

f_base = 10

f_max = 500

amp, freq, rpms = freq_wave(freq_l=freq_l, freq_u=freq_u, step=step, f_base=f_base, f_max=f_max)

# 周波数応答解析を実行

X_list = []

for i in range(len(amp)):

freq_axis, X = frequency_response_analysis(M=M,

K=K,

C=C,

steps=freq,

exc_node=1,

res_node=1,

excitation=amp[i])

X_list.append(X)

# 固有値解析を実行(確認用)

omega_sort, v_sort, index= eigen(M, K)

print('Natural frequency[Hz]=', omega_sort)

# ここからグラフ描画

# グラフをオブジェクト指向で作成する。

fig = plt.figure()

ax1 = fig.add_subplot(111)

# データをプロットする。

im = ax1.pcolormesh(rpms,

freq_axis,

np.array(X_list).T,

cmap='jet',

vmin=0, vmax=np.max(X_list))

# カラーバーを設定する。

cbar = fig.colorbar(im)

cbar.set_label('Displacement[m]')

# 軸設定する。

ax1.set_xlabel('Rev.[rpm]')

ax1.set_ylabel('Frequency[Hz]')

# グラフを表示する。

plt.show()

plt.close()

単一次数の結果

上記コードの実行結果は下図です。共振である400[Hz]付近に次数が当たった部分で振幅大となっています。

高調波次数も追加してみる

続いて、先ほどの加振力セットの関数内のdata部分コメント蘭を切り替えて、のこぎり波バージョンにしてみます。のこぎり波は基本周波数の高調波が出るので、今回のように高調波を出したい時によく使います。

# 擬似加振力スイープ周波数波形の作成
def freq_wave(freq_l, freq_u, step, f_base, f_max):
    samplerate = freq_u * 2.56                              # サンプリングレート[Hz]は上限周波数の2.56倍
    dt = 1 / samplerate                                     # 時間刻み[s]
    end = 1 / step                                          # 周波数分解能を周波数応答解析のステップにあわせる[s]
    t = np.arange(0, end + dt, dt)                          # 時間軸[s]

    n = int((freq_u - freq_l) / step)                       # 計算ステップ数
    f_increment = (f_max - f_base) / n                      # ステップ当たりの回転数[rps]変化

    amp_list = []
    rpms = []
    for i in range(n):
        freq_i = f_base + (i + 1) * f_increment
        rpms.append(freq_i * 60)                            # 基本周波数をRPMに変換

        data = signal.sawtooth(2 * np.pi * freq_i * t)     # のこぎり波
        #data = np.sin(2 * np.pi * freq_i * t)               # 正弦波

        amp, freq = calc_fft(data, samplerate)              # フーリエ変換
        amp_sub = amp[int(freq_l / step):int(freq_u)]       # 周波数応答解析で使用する範囲のみ振幅を抽出
        amp_list.append(amp_sub)

    freq_sub = freq[int(freq_l / step):int((freq_u + step) / step)]    # 周波数応答解析で使用する範囲のみ周波数軸を抽出

    return amp_list, freq_sub, rpms

# 擬似加振力スイープ周波数波形の作成

def freq_wave(freq_l, freq_u, step, f_base, f_max):

samplerate = freq_u * 2.56 # サンプリングレート[Hz]は上限周波数の2.56倍

dt = 1 / samplerate # 時間刻み[s]

end = 1 / step # 周波数分解能を周波数応答解析のステップにあわせる[s]

t = np.arange(0, end + dt, dt) # 時間軸[s]

n = int((freq_u - freq_l) / step) # 計算ステップ数

f_increment = (f_max - f_base) / n # ステップ当たりの回転数[rps]変化

amp_list = []

rpms = []

for i in range(n):

freq_i = f_base + (i + 1) * f_increment

rpms.append(freq_i * 60) # 基本周波数をRPMに変換

data = signal.sawtooth(2 * np.pi * freq_i * t) # のこぎり波

#data = np.sin(2 * np.pi * freq_i * t) # 正弦波

amp, freq = calc_fft(data, samplerate) # フーリエ変換

amp_sub = amp[int(freq_l / step):int(freq_u)] # 周波数応答解析で使用する範囲のみ振幅を抽出

amp_list.append(amp_sub)

freq_sub = freq[int(freq_l / step):int((freq_u + step) / step)] # 周波数応答解析で使用する範囲のみ周波数軸を抽出

return amp_list, freq_sub, rpms

結果は以下です。複数次数が発生していますが、共振に当たる400[Hz]付近で振幅レベルが増大している結果を得る事ができました。

ちなみに800[Hz]ちょいにも共振があり、薄ら明るくなっている事に気づきましたでしょうか（2自由度なので）？

全コード内には固有値解析コードも付いていますので、是非パラメータを色々変えて遊んでみて下さい。

まとめ

本記事では「振動モデルを直接法の周波数応答解析で解くPythonコード」で学んだ周波数応答解析の応用として、小修正で加振力を変化させてみました。

特に新しい式やアルゴリズムは登場していませんが、学んだ技術でこんな事ができるのは非常に面白いと思いました。

ちょっと考えた事がすぐ実現できるのはPythonのメリットだと思いますので、是非皆さんも始めてみて下さい！

周波数応答解析は便利なので、加振力を色々変えてみてシミュレーションしてみるのも面白いですね！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

本記事の対象者

加振力の変化を考慮した周波数応答解析をしたい人

加振力の変化を考慮した周波数応答解析のPythonコード

動作環境

PC

Python環境

モデル定義

変化する加振力を作成

補足：実用的な加振力データを作るには

周波数応答解析の関数を修正

全コード(コピペ用)

単一次数の結果

高調波次数も追加してみる

まとめ

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

本記事の対象者

加振力の変化を考慮した周波数応答解析をしたい人

加振力の変化を考慮した周波数応答解析のPythonコード

動作環境

PC

Python環境

モデル定義

変化する加振力を作成

補足：実用的な加振力データを作るには

周波数応答解析の関数を修正

全コード(コピペ用)

単一次数の結果

高調波次数も追加してみる

まとめ

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル