拟合衰减振动模型，估算阻尼比和阻尼系数

flyfish

衰减振动模型

在自由振动系统中，阻尼振动可以用以下公式描述：
$x(t)=x_0{e}^{-\zeta {\omega }_{n}t}\cos ({\omega }_{d}t+\varphi )$
其中：

$x(t)$ 是时间 $t$ 时的位移（Displacement at time $t$）

$x_0$ 是初始位移（Initial displacement）

$\zeta$ 是阻尼比（Damping ratio）
希腊字母 $\zeta$，英文表示为 “zeta”

${\omega }_n$ 是无阻尼固有频率（Undamped natural frequency）
希腊字母 $\omega$，英文表示为 “omega”

$t$ 是时间（Time）

${\omega }_d$ 是阻尼振动频率（Damped natural frequency），其公式为：
${\omega }_d={\omega }_n\sqrt{1-{\zeta }^2}$

$\varphi$ 是相位角（Phase angle）
希腊字母 $\varphi$，英文表示为 “phi”

公式说明

1. 初始位移 $x_0$ (Initial displacement) :
   系统开始自由振动时的位移。

2. 阻尼比 $\zeta$ (Damping ratio) :
   衡量系统阻尼程度的无量纲参数。
   $\zeta$ 越大，阻尼越强。

3. 无阻尼固有频率 ${\omega }_n$ (Undamped natural frequency) :
   系统在没有阻尼情况下的固有振动频率。
   单位：弧度每秒（radians per second, rad/s）。

4. 阻尼振动频率 ${\omega }_d$ (Damped natural frequency) :
   系统在有阻尼情况下的实际振动频率。

计算公式： ${\omega }_d={\omega }_n\sqrt{1-{\zeta }^2}$

当 $\zeta <1$ 时，系统会振荡，${\omega }_d$ 表示振荡频率。

5. 相位角 $\varphi$ (Phase angle) :
   描述系统初始位移和初始速度之间的关系。

衰减振动的描述

在阻尼系统中，振动会逐渐衰减，幅度随着时间指数下降。振动系统的运动可以分解为以下几个部分：

1. 指数衰减部分 $e^{-\zeta {\omega }_{n}t}$ :
   描述振幅随时间的衰减。

2. 余弦振荡部分 $\cos ({\omega }_{d}t+\varphi )$ :
   描述系统的振荡行为，频率为 ${\omega }_d$，初始相位为 $\varphi$。

代码实现

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

# 定义衰减振动函数
def damped_vibration(t, x0, zeta, omega_n, phi):
    omega_d = omega_n * np.sqrt(1 - zeta**2)
    return x0 * np.exp(-zeta * omega_n * t) * np.cos(omega_d * t + phi)

# 合成实验数据
t_data = np.linspace(0, 10, 1000)  # 时间数据
x0_true = 1.0  # 初始位移
zeta_true = 0.1  # 真正的阻尼比
omega_n_true = 2 * np.pi  # 真正的无阻尼固有频率
phi_true = 0  # 真正的相位角
x_data = damped_vibration(t_data, x0_true, zeta_true, omega_n_true, phi_true)  # 生成理想数据

# 加入噪声
x_data_noisy = x_data + 0.05 * np.random.normal(size=t_data.shape)  # 加入随机噪声

# 拟合数据
initial_guess = [x0_true, zeta_true, omega_n_true, phi_true]  # 初始猜测
params, params_covariance = opt.curve_fit(damped_vibration, t_data, x_data_noisy, p0=initial_guess)  # 拟合模型

# 提取拟合参数
x0_est, zeta_est, omega_n_est, phi_est = params

# 假设质量为1 kg
m = 1.0
gamma_est = 2 * zeta_est * omega_n_est * m

# 输出结果
print(f"Estimated Initial Displacement (x0): {x0_est}")
print(f"Estimated Damping Ratio (zeta): {zeta_est}")
print(f"Estimated Natural Frequency (omega_n): {omega_n_est}")
print(f"Estimated Damping Coefficient (gamma): {gamma_est} N·s/m")

# 绘图比较
plt.figure(figsize=(10, 6))
plt.plot(t_data, x_data_noisy, label='Noisy Data', alpha=0.5)
plt.plot(t_data, damped_vibration(t_data, *params), label='Fitted Curve', color='red')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.legend()
plt.show()

t_data: 时间数据，生成一个从0到10秒的1000个时间点的数组。
x0_true: 初始位移（True initial displacement），设定为1.0米。
zeta_true: 真正的阻尼比（True damping ratio），设定为0.1。
omega_n_true: 真正的无阻尼固有频率（True undamped natural frequency），设定为 2π（即每秒一个完整的振动周期）。
phi_true: 真正的相位角（True phase angle），设定为0。

Estimated Initial Displacement (x0): 0.9969144391835894
Estimated Damping Ratio (zeta): 0.09975574475118448
Estimated Natural Frequency (omega_n): 6.2712971489170615
Estimated Damping Coefficient (gamma): 1.2511958352924026 N·s/m

每个结果的单位及其读法：

1. Estimated Initial Displacement (x0) :
   单位 : 米 (meters, m)

读法 : “Estimated Initial Displacement is value meters”

2. Estimated Damping Ratio (zeta) :
   单位 : 无单位（damping ratio 是一个无量纲参数）

读法 : “Estimated Damping Ratio is value”

3. Estimated Natural Frequency (omega_n) :
   单位 : 弧度每秒 (radians per second, rad/s)

读法 : “Estimated Natural Frequency is value radians per second”

4. Estimated Damping Coefficient (gamma) :
   单位 : 牛顿·秒每米 (Newton-seconds per meter, N·s/m)

读法 : “Estimated Damping Coefficient is value Newton-seconds per meter”