Pythonによるバネマスダンパ系のシミュレーション

はじめに

本記事ではバネマスダンパ系をPythonでシミュレーションする方法を説明します。
バネマスダンパ系は自動車のサスペンションやロボットアームなどあらゆる振動現象に適用可能な基本的な力学システムです。
実際の業務ではMATLAB・Simulinkを使うケースが多いと思いますが

• MATLABが使えない環境の方
• 自分の書いたコードでシミュレーションしたい

という方向けに記事を書いてみました。

環境

• Python 3.10.11

シミュレーション対象

バネ定数$ k $のバネ、質量$ m $の物体、減衰係数$ c $のダンパが接続された普通のバネマスダンパ系です。
位置$ x $と速度$ v $に初期値を与え自由振動させるモデルを解いていきます。

運動方程式を定義する

バネマスダンパ系の運動方程式は以下のようになります。一般的なシステムなので導出は省略です。

\frac{d^2 x}{d t^2}+c \frac{d x}{d t}+k x=0

位置$ x $と速度$ v $を状態変数とする状態方程式に変換しておきます。

\left\{\begin{array}{c}
\frac{d x}{d t}=v \\
\frac{d v}{d t}=-\frac{k}{m} x-\frac{c}{m} v
\end{array}\right.

オイラー法の適用（前進オイラー法）

上記の状態方程式に前進オイラー法を適用し、数値積分によりプログラムで解ける状態にします。
位置$ x $は以下の式で更新していけば求まります。

x_{i+1}=x_i+\frac{d x}{d t} \Delta t=x_i+v_i \Delta t

続いて速度$ v $は以下の式で更新します。

v_{i+1}=v_i+\frac{d v}{d t} \Delta t=v_i+\left(-\frac{k}{m} x-\frac{c}{m} v_i\right) \Delta t

これで実装する準備ができました。

Pythonによる実装

spring_mass_damper_system.py

import numpy as np
import matplotlib.pyplot as plt

# 定数の設定
m = 1.0  # 質量[kg]
k = 2.0  # バネ定数[N/m]
c = 0.5  # 減衰係数[Ns/m]
timestep = 0.01  # 時間ステップ[sec]
total_time = 10.0  # 総時間[sec]

# 初期条件
x_init = 1.0  # 初期変位[m]
v_init = 0.0  # 初期速度[m/s]

# 時間配列の生成
time_steps = int(total_time / timestep)
t = np.linspace(0, total_time, time_steps)

# 変位と速度の配列
x = np.zeros(time_steps)
v = np.zeros(time_steps)

# 初期条件の設定
x[0] = x_init
v[0] = v_init

# オイラー法による数値積分
for i in range(1, time_steps):
    x[i] = x[i-1] + v[i-1] * timestep
    v[i] = v[i-1] + (-(k/m) * x[i-1] - (c/m) * v[i-1]) * timestep

# グラフの描画
plt.figure(figsize=(10, 5))
plt.plot(t, x, label='Displacement (x)')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Spring-Mass-Damper System')
plt.legend()
plt.grid(True)
plt.show()

シミュレーション結果

このプログラムを実行すると変位$ x $の時間経過がグラフに表示されます。
減衰しながら振動している様子が確認できます。

パラメータスタディ

パラメータを変更して結果を比較してみます。
振動減衰、臨界減衰、過減衰となるように減衰係数$ c $を設定してシミュレーションすると以下のようになります。

このコードも一応貼っておきます。

parameter_stady.py

import numpy as np
import matplotlib.pyplot as plt

# 定数の設定
m = 1.0  # 質量[kg]
k = 2.0  # バネ定数[N/m]
timestep = 0.01  # 時間ステップ[sec]
total_time = 10.0  # 総時間[sec]

# 初期条件
x_init = 1.0  # 初期変位[m]
v_init = 0.0  # 初期速度[m/s]

# 時間配列の生成
time_steps = int(total_time / timestep)
t = np.linspace(0, total_time, time_steps)

# 変位と速度の配列を初期化
x_underdamped = np.zeros(time_steps)
v_underdamped = np.zeros(time_steps)
x_critically_damped = np.zeros(time_steps)
v_critically_damped = np.zeros(time_steps)
x_overdamped = np.zeros(time_steps)
v_overdamped = np.zeros(time_steps)

# 初期条件の設定
x_underdamped[0] = x_init
v_underdamped[0] = v_init
x_critically_damped[0] = x_init
v_critically_damped[0] = v_init
x_overdamped[0] = x_init
v_overdamped[0] = v_init

# 減衰係数の設定
c_underdamped = 0.5  # 振動減衰
c_critically_damped = 2 * np.sqrt(m * k)  # 臨界減衰
c_overdamped = 5.0  # 過減衰

# 数値積分（オイラー法）によるシミュレーション
for i in range(1, time_steps):
    # 振動減衰
    x_underdamped[i] = x_underdamped[i-1] + v_underdamped[i-1] * timestep
    v_underdamped[i] = v_underdamped[i-1] + (-(k/m) * x_underdamped[i-1] - (c_underdamped/m) * v_underdamped[i-1]) * timestep
    
    # 臨界減衰
    x_critically_damped[i] = x_critically_damped[i-1] + v_critically_damped[i-1] * timestep
    v_critically_damped[i] = v_critically_damped[i-1] + (-(k/m) * x_critically_damped[i-1] - (c_critically_damped/m) * v_critically_damped[i-1]) * timestep
    
    # 過減衰
    x_overdamped[i] = x_overdamped[i-1] + v_overdamped[i-1] * timestep
    v_overdamped[i] = v_overdamped[i-1] + (-(k/m) * x_overdamped[i-1] - (c_overdamped/m) * v_overdamped[i-1]) * timestep

# グラフの描画
plt.figure(figsize=(12, 6))
plt.plot(t, x_underdamped, label='Underdamped (c = 0.5)')
plt.plot(t, x_critically_damped, label='Critically Damped (c = 2*sqrt(m*k))')
plt.plot(t, x_overdamped, label='Overdamped (c = 5.0)')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Spring-Mass-Damper System')
plt.legend()
plt.grid(True)
plt.show()

おわりに

この記事では、バネマスダンパ系をPythonでシミュレーションする方法を紹介しました。
実際に手を動かしてみることで、力学システムの理解が深まると思いますので参考にしていただけると幸いです。