質量・ばね・ダンパ系は三つの値から応答を予測できる
機械の振動を検討するとき、「固有振動数だけ計算して終わる」と設計判断を誤ることがあります。同じ固有振動数でも、減衰が小さければ振動が長く残り、減衰が大きければゆっくり原点へ戻ります。位置決め機構、除振台、支持構造、アクチュエータなどの過渡応答を比べるには、質量、ばね定数、減衰係数を同じモデルで扱う必要があります。
本記事では、並進方向の1自由度・線形・粘性減衰系を対象に、固有角振動数、減衰比、減衰固有角振動数、オーバーシュート、整定時間の目安を計算します。最後にPythonで微分方程式を数値積分し、計算式と時刻歴を同じ入力から再現します。
運動方程式と前提をそろえる
質量を $m$、粘性減衰係数を $c$、ばね定数を $k$、変位を $x$、外力を $F(t)$ とすると、運動方程式は次式です。
$$ m\ddot{x}+c\dot{x}+kx=F(t) $$
各入力のSI単位は次のとおりです。
| 量 | 記号 | SI単位 |
|---|---|---|
| 等価質量 | m | kg |
| 粘性減衰係数 | c | N·s/m |
| ばね定数 | k | N/m |
| 変位 | x | m |
| 外力 | F | N |
このモデルでは、ばね力が変位に比例し、減衰力が速度に比例すると仮定します。摩擦、がた、材料のヒステリシス、大変形、接触離脱などが支配的な装置は、この線形モデルだけでは表せません。最初にモデルの適用範囲を決め、実測との差を後から確認します。
固有振動数と減衰比を計算する
非減衰固有角振動数と固有振動数
減衰と外力を除いた系の固有角振動数は次式です。
$$ \omega_n=\sqrt{\frac{k}{m}} $$
$\omega_n$ の単位はrad/sです。Hzで表す固有振動数は次式で求めます。
$$ f_n=\frac{\omega_n}{2\pi} $$
質量を増やすと固有振動数は低下し、ばね定数を増やすと上昇します。ただし実機では、治具、ワーク、可動配線、ばね自身の有効質量などを等価質量へ含める必要があります。
臨界減衰係数と減衰比
臨界減衰係数は次式です。
$$ c_c=2\sqrt{km}=2m\omega_n $$
実際の減衰係数を臨界減衰係数で割った無次元量が減衰比です。
$$ \zeta=\frac{c}{c_c}=\frac{c}{2\sqrt{km}} $$
減衰比によって自由応答の性質を分類できます。
| 減衰比 | 分類 | 自由応答の特徴 |
|---|---|---|
| 0 | 非減衰 | 一定振幅で振動する理想モデル |
| 0より大、1未満 | 不足減衰 | 振動しながら減衰する |
| 1 | 臨界減衰 | 振動せず、境界的に原点へ戻る |
| 1より大 | 過減衰 | 振動せず、二つの指数成分で戻る |
「減衰を増やせば必ず速く整定する」とは限りません。臨界減衰を大きく超えると、振動は消えても遅い指数成分が支配し、戻りが遅くなることがあります。
不足減衰系の応答量を求める
減衰固有角振動数
$0<\zeta<1$ の不足減衰系では、実際に観測される自由振動の角振動数は次式です。
$$ \omega_d=\omega_n\sqrt{1-\zeta^2} $$
対応する周波数は $f_d=\omega_d/(2\pi)$ です。減衰が小さいと $\omega_d$ は $\omega_n$ に近くなりますが、両者を同じ値として保存すると、減衰の変更時に誤差の原因になります。
初期変位を $x_0$、初期速度を $v_0$、外力をゼロとすると、自由応答は次式です。
$$ x(t)=e^{-\zeta\omega_n t}\left[x_0\cos(\omega_d t)+\frac{v_0+\zeta\omega_n x_0}{\omega_d}\sin(\omega_d t)\right] $$
包絡線は $e^{-\zeta\omega_n t}$ に従って減少します。したがって、減衰比だけでなく積 $\zeta\omega_n$ が振幅の減衰速度を決めます。
ステップ応答のオーバーシュート
標準的な2次遅れ系で $0<\zeta<1$、かつ応答を大きく変える零点がない場合、単位ステップ応答の最大オーバーシュート率は次式で近似できます。
$$ M_p=\exp\left(-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}\right) $$
百分率で表すときは $100M_p$ とします。この式は自由振動の初期値応答ではなく、標準2次系のステップ応答に対する式です。入力形状、零点、飽和、摩擦がある系へ無条件に適用しません。
整定時間の目安
不足減衰の標準2次系では、2%整定時間の代表的な近似として次式が使われます。
$$ t_s\approx\frac{4}{\zeta\omega_n} $$
これは設計初期の比較に便利ですが、厳密な整定時間ではありません。評価帯、初期条件、入力、離散サンプリング、非線形要素によって結果は変わります。最終判定では時刻歴から「応答が許容帯へ入り、その後に出ない最初の時刻」を求めます。
数値例：12 kgの機構を4800 N/mで支持する
次の単純化した機構を考えます。
| 項目 | 値 |
|---|---|
| 等価質量 m | 12 kg |
| ばね定数 k | 4800 N/m |
| 粘性減衰係数 c | 180 N·s/m |
| 初期変位 x_0 | 0.010 m |
| 初期速度 v_0 | 0 m/s |
非減衰固有角振動数は次のとおりです。
$$ \omega_n=\sqrt{\frac{4800}{12}}=20.0\ \mathrm{rad/s} $$
固有振動数は次の値です。
$$ f_n=\frac{20.0}{2\pi}\approx3.18\ \mathrm{Hz} $$
臨界減衰係数は次のとおりです。
$$ c_c=2\sqrt{4800\times12}=480\ \mathrm{N\cdot s/m} $$
したがって減衰比は次の値です。
$$ \zeta=\frac{180}{480}=0.375 $$
$0<\zeta<1$ なので不足減衰です。減衰固有角振動数は次のようになります。
$$ \omega_d=20\sqrt{1-0.375^2}\approx18.54\ \mathrm{rad/s} $$
対応する減衰固有振動数は約2.95 Hzです。標準2次系として見た最大オーバーシュート率は約28.1%、2%整定時間の近似は次の値です。
$$ t_s\approx\frac{4}{0.375\times20}\approx0.533\ \mathrm{s} $$
この結果から、固有振動数だけを見るのではなく、許容オーバーシュートと整定時間を満たす減衰の検討が必要だと分かります。
Pythonで指標計算と時刻歴解析を行う
次のコードは、入力値から振動指標を計算し、scipy.integrate.solve_ivpで自由応答を数値積分します。微分方程式を変位と速度の二つの一階方程式へ変換しています。
from dataclasses import asdict, dataclass
from math import exp, pi, sqrt
import numpy as np
from scipy.integrate import solve_ivp
@dataclass(frozen=True)
class OscillatorInput:
mass_kg: float
damping_n_s_m: float
stiffness_n_m: float
initial_displacement_m: float
initial_velocity_m_s: float
end_time_s: float = 2.0
samples: int = 1001
def calculate_indicators(x: OscillatorInput) -> dict[str, float | str]:
if x.mass_kg <= 0 or x.stiffness_n_m <= 0:
raise ValueError("mass and stiffness must be positive")
if x.damping_n_s_m < 0:
raise ValueError("damping must be non-negative")
if x.end_time_s <= 0 or x.samples < 2:
raise ValueError("invalid time range")
omega_n = sqrt(x.stiffness_n_m / x.mass_kg)
critical_damping = 2 * sqrt(x.stiffness_n_m * x.mass_kg)
damping_ratio = x.damping_n_s_m / critical_damping
if damping_ratio < 1:
response_type = "underdamped"
omega_d = omega_n * sqrt(1 - damping_ratio**2)
overshoot_ratio = (
exp(-pi * damping_ratio / sqrt(1 - damping_ratio**2))
if damping_ratio > 0
else 1.0
)
elif damping_ratio == 1:
response_type = "critical"
omega_d = 0.0
overshoot_ratio = 0.0
else:
response_type = "overdamped"
omega_d = 0.0
overshoot_ratio = 0.0
settling_time_2pct_s = (
4 / (damping_ratio * omega_n)
if 0 < damping_ratio < 1
else float("nan")
)
return {
"natural_angular_frequency_rad_s": omega_n,
"natural_frequency_hz": omega_n / (2 * pi),
"critical_damping_n_s_m": critical_damping,
"damping_ratio": damping_ratio,
"response_type": response_type,
"damped_angular_frequency_rad_s": omega_d,
"damped_frequency_hz": omega_d / (2 * pi),
"step_overshoot_percent": 100 * overshoot_ratio,
"settling_time_2pct_s": settling_time_2pct_s,
}
def simulate_free_response(x: OscillatorInput):
def equation(_time, state):
displacement, velocity = state
acceleration = (
-x.damping_n_s_m * velocity
- x.stiffness_n_m * displacement
) / x.mass_kg
return [velocity, acceleration]
time = np.linspace(0.0, x.end_time_s, x.samples)
solution = solve_ivp(
equation,
(0.0, x.end_time_s),
[x.initial_displacement_m, x.initial_velocity_m_s],
t_eval=time,
rtol=1e-9,
atol=1e-12,
)
if not solution.success:
raise RuntimeError(solution.message)
return solution.t, solution.y[0], solution.y[1]
example = OscillatorInput(
mass_kg=12.0,
damping_n_s_m=180.0,
stiffness_n_m=4800.0,
initial_displacement_m=0.010,
initial_velocity_m_s=0.0,
)
print(asdict(example))
for name, value in calculate_indicators(example).items():
print(f"{name}: {value}")
time_s, displacement_m, velocity_m_s = simulate_free_response(example)
np.savetxt(
"mass-spring-damper-response.csv",
np.column_stack([time_s, displacement_m, velocity_m_s]),
delimiter=",",
header="time_s,displacement_m,velocity_m_s",
comments="",
)solve_ivpは初期値問題を数値積分する関数です。外力を加えたい場合は、加速度の分子へ $F(t)$ を追加します。入力波形を関数として分離すれば、ステップ、正弦波、台形加減速などを同じモデルで比較できます。
実務で入力値を決める方法
等価質量は「動くもの」を機構方向へ換算する
直線運動では可動部とワークの質量を合計します。回転体を含む場合は、伝達比や半径を使って対象座標の等価質量へ換算します。別の計算で慣性トルクとして加えた回転体を、等価質量へ再び入れると二重計上になります。
ばね定数は部品単体ではなく荷重点で評価する
コイルばねだけでなく、フレーム、軸、軸受、ボルト締結部、接触部にも弾性があります。直列要素と並列要素を整理し、応答を評価する荷重点での等価ばね定数を使います。有限要素解析から荷重と変位を得た場合、線形範囲では $k=F/x$ として等価剛性を求められます。
減衰係数は推測値を確定値として扱わない
減衰は材料、締結、潤滑、シール、摩擦、流体、温度によって変わります。設計初期は複数の減衰比で感度解析し、試作後は自由減衰波形や周波数応答から同定します。線形粘性減衰で実測波形を再現できない場合は、クーロン摩擦や速度依存性を含むモデルへ進みます。
よくある失敗
- 固有角振動数rad/sと固有振動数Hzを混同する
- 部品単体のばね定数を、機構全体の荷重点剛性として使う
- ワーク交換による質量変化を計算ケースへ入れない
- 減衰比と減衰係数を同じ量として扱う
- 不足減衰用の式を、臨界減衰や過減衰へそのまま使う
- ステップ応答のオーバーシュート式を自由応答へ流用する
- 2%整定時間の近似だけで最終合否を決める
- 摩擦、がた、飽和、接触離脱を線形粘性減衰だけで表す
- 数値計算の時間刻みや許容誤差を変えずに高周波系へ流用する
設計データと計算ロジックへ展開する
計算を再利用するなら、mass_kg、damping_n_s_m、stiffness_n_mのように物理量と単位を変数名へ含めます。結果には固有振動数だけでなく、減衰比、応答分類、評価条件、モデル版、入力値の出典を保存します。
複数案を比較する場合は、質量・剛性・減衰の組合せをCSVまたはJSONで与え、同じ関数で一括計算するとレビューしやすくなります。将来Webツール化するときも、入力検証、単位変換、計算本体、時刻歴生成、合否判定を分ければ、式の変更や試験値への置換が容易です。
まとめ
1自由度の質量・ばね・ダンパ系では、$m$、$c$、$k$から非減衰固有角振動数 $\omega_n$、臨界減衰係数 $c_c$、減衰比 $\zeta$を求められます。減衰比は不足減衰、臨界減衰、過減衰の判定に使い、不足減衰では $\omega_d$、オーバーシュート、整定時間の目安へ展開できます。
ただし、式の精度はモデルの妥当性に依存します。等価質量と荷重点剛性をそろえ、減衰は感度解析と実測で更新します。指標計算と微分方程式の時刻歴解析を同じ入力データから実行できる形にしておくと、設計変更、試験との照合、計算ツール化へつなげやすくなります。
