はじめに
強制振動では、外力により変位がどのように変化するのかを調査することができる。
特に、系固有の周波数と外力に含まれる振動数が一致する場合（共振）は、他の周波数成分の場合と比べて応答が大きく変化すると考えられる。
そこで、今回はPythonを用いて1質点,1バネ系に、固有周波数と一致する周波数を持つ方形波の外力を加えた場合、どのように応答値（変位）が変化するのかを調査する。
問題設定
1質点($m$),1バネ($k$)系に外力F(t)を加えて応答値である変位$x$（バネの伸び）をみる。力学モデルは以下の通り。
m\ddot{x}=-kx +F(t)
ここで、系の固有周波数と一致する周波数の方形波の外力$F(t)$を加える。
$F(t)$は方形波である。したがって、フーリエ級数展開を用いると、以下のように表現することができる。
F(t)=\frac{4}{\pi}\sum_{n=1}^{\infty} \frac{1}{2n-1} \sin (2n-1)\theta
ただし、$\theta=\omega t=2\pi f t$とする。
したがって、固有周波数と一致する振動数の外力を受けた場合、系は共振し、変位は大きく変化すると考えられる。
ただし、今回の系は減衰項を持たないため、共振時には変位振幅が時間とともに発散し続ける点に注意が必要である。
数値計算
以下のようなプログラムを作成して、数値計算を実施した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
b=10
#フーリエ変換のための設定
sampling_rate = 1000 # サンプリング周波数（Hz）
n=sampling_rate
dt=b/(n-1) # 実際の時間刻み（tと整合させる）
#質量
m=1
#方形波周波数
f=2
#バネ定数（方形波の周波数fに固有周波数を一致させ、共振させる）
k=m*(2*math.pi*f)**2
omega=2*math.pi*f
t=np.linspace(0,b,n)
#方形波の生成
wave = np.sign(np.sin(2 * np.pi * f * t))
x=np.zeros(n)
v=np.zeros(n)
a=np.zeros(n)
F_ary=np.zeros(n)
for i in range(n):
if i==0:
x[i]=0
v[i]=0
a[i]=0
else:
F=wave[i]
F_ary[i]=F
a[i]=F/m-k*x[i-1]/m
v[i]=v[i-1]+a[i]*dt
x[i]=x[i-1]+v[i]*dt
plt.plot(t,x,label="変位")
plt.plot(t,v,label="速度")
plt.plot(t,a,label="加速度")
plt.plot(t,F_ary,label="外力")
plt.xlabel("時間[s]")
plt.ylabel("変位、速度、加速度、外力")
plt.legend()
plt.savefig("強制振動_方形波.png")
plt.show()
# フーリエ変換のための信号を用意（変位）
signal=x
# FFT実行
fft_result = np.fft.fft(signal)
# 周波数軸の計算
freqs = np.fft.fftfreq(len(signal), dt)
# FFT結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title('Original Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.subplot(2, 1, 2)
plt.stem(freqs[:len(freqs)//2], np.abs(fft_result)[:len(freqs)//2], 'b', markerfmt=" ", basefmt="-b")
plt.title('FFT Result')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.tight_layout()
plt.savefig("強制振動_方形波_FFT.png")
plt.show()
結果
結果は以下のようになった。
応答値の時間応答
外力を加えた際の変位、速度、加速度の応答を調査した。
結果、速度、加速度は外力の波形の影響を受けやすい。
一方、変位は方形波の周波数成分と一致する固有周波数によって共振し、振幅が時間とともにほぼ線形に増大し続ける挙動を示した（例えば t=2sで振幅約0.10、t=10sで振幅約0.50）。
これは、今回の系が減衰項を持たない共振系であるため、外力の投入したエネルギーが散逸せず蓄積し続けることによるものである。
変位のフーリエ変換
そこで、変位について、フーリエ変換を施すことで、どの周波数帯の影響が支配的になるか調査した。
結果、固有周波数（=方形波の周波数）に一致する成分が支配的になることが分かった。
まとめ
今回はPythonを用いて、強制振動の変位を調査した。
具体的には、入力として外力に、系の固有周波数と一致する周波数の方形波を加え、変位の応答値を調査した。
結果、変位については固有周波数と一致する周波数の振動が支配的になり、かつ減衰のない系であるため振幅が時間とともに発散し続けることが分かった。
このようなことから、質量とバネの系には、固有周波数と外力の周波数が一致した場合にのみ応答が大きく成長するという周波数選択性(フィルタのようなもの)があると考えられる。
ただし、今回のように減衰がない場合は定常的な応答には収束せず、共振周波数では振幅が発散し続ける点には注意が必要である。
参考文献
- フーリエ級数展開と方形波
- Numpyでのフーリエ変換