基于Python的非线性振动控制模型与MR阻尼器集成技术咨询

Got it, let's break down how to integrate your MR Damper model into the existing state-space system. Here's a practical, step-by-step approach with modified code to make it work:

First, Let's Align the System Logic #

Your original linear model uses a constant damping coefficient c, but the MR damper produces a nonlinear force dependent on relative velocity, relative displacement, and input current. We need to replace the linear damping term with this nonlinear force and solve the resulting nonlinear ODE (since signal.lsim only works for linear systems).

Step 1: Redefine the Nonlinear System Equations #

Your single-DOF system (with linear damping) was:
$$m\ddot{x}1 + c(\dot{x}1 - \dot{u}) + k(x_1 - u) = 0$$
Replace the linear damping force with the MR damper force $F{mr}$:
$$m\ddot{x}1 + F{mr} + k(x_1 - u) = 0$$
Where $F{mr}$ uses:

• Relative velocity: $v = \dot{x}_1 - \dot{u}$ (mass velocity minus base excitation velocity)
• Relative displacement: $x_{rel} = x_1 - u$ (mass displacement minus base excitation displacement)
• Input current $i$ (you can fix this or make it time-varying)

Step 2: Rewrite as First-Order ODEs #

Define state variables:

• $x_0 = x_1$ (mass displacement)
• $x_1 = \dot{x}_1$ (mass velocity)

The first-order equations become:
$$
\begin{cases}
\dot{x}_0 = x_1 \
\dot{x}1 = \frac{ -k(x_0 - u(t)) - F{mr} }{ m }
\end{cases}
$$

Step 3: Modified Code with MR Damper Integration #

Here's the full working code, with comments explaining key changes:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# --- Reuse your existing helper functions ---
def rms(v):
    return np.sqrt(np.mean(np.square(v)))

def transmissibility(_in, _out):
    return rms(_out) / rms(_in)

# --- MR Damper Function (unchanged, but we'll use relative displacement directly) ---
def MR(v, x_rel, i):
    """
    @v: relative velocity (mass velocity - base velocity)
    @x_rel: relative displacement (mass displacement - base displacement)
    @i: input current (A)
    return: MR damper force (N)
    """
    b1 = 3415.7
    b2 = 93.324
    b3 = 74.487
    F0 = b1 * (i**2) + b2 * i + b3
    
    b4 = 2534.1
    b5 = 19.55
    b6 = 643.1
    C1 = b4 * (i**2) + b5*i + b6
    
    beta = 50
    p1 = 4
    p2 = 0.2
    
    part1 = F0 * np.tanh(beta * (v + p1 * x_rel))
    part2 = C1 * (v + p2 * x_rel)
    return part1 + part2

# --- Nonlinear System ODE for solve_ivp ---
def nonlinear_system(t, x, k, m, u_func, u_dot_func, i):
    x0, x1 = x  # x0 = mass displacement, x1 = mass velocity
    u_t = u_func(t)
    u_dot_t = u_dot_func(t)
    
    # Calculate relative motion
    x_rel = x0 - u_t
    v_rel = x1 - u_dot_t
    
    # Get MR damper force
    F_mr = MR(v_rel, x_rel, i)
    
    # ODEs
    dx0_dt = x1
    dx1_dt = (-k * x_rel - F_mr) / m
    return [dx0_dt, dx1_dt]

# --- Function to compute transmissibility for a given frequency ---
def compute_transmissibility(f, t, amplitude, k, m, i):
    # Base excitation: u(t) = amplitude * sin(2πft)
    u_func = lambda t: amplitude * np.sin(2 * np.pi * f * t)
    # Derivative of excitation (velocity)
    u_dot_func = lambda t: amplitude * 2 * np.pi * f * np.cos(2 * np.pi * f * t)
    
    # Solve nonlinear ODE
    sol = solve_ivp(
        fun=nonlinear_system,
        t_span=[t[0], t[-1]],
        y0=[0, 0],  # initial conditions: displacement=0, velocity=0
        t_eval=t,
        args=(k, m, u_func, u_dot_func, i)
    )
    
    # Mass displacement is sol.y[0], excitation is u_func(t)
    mass_displacement = sol.y[0]
    excitation = u_func(t)
    return transmissibility(excitation, mass_displacement)

# --- Main Simulation Setup ---
# Parameters (adjust these to match your original model)
k = 10000  # stiffness (N/m), example value
m = 100    # mass (kg), example value
i_mr = 0.2 # MR damper current (A)
amplitude = 0.05  # excitation amplitude (m)
freqs = np.linspace(1, 20, 50)  # frequency range (Hz)
t = np.linspace(0, 10, 1000)    # time vector (long enough to reach steady state)

# Compute transmissibility for each frequency
trans_mr = [compute_transmissibility(f, t, amplitude, k, m, i_mr) for f in freqs]

# Plot results
plt.figure(figsize=(12, 8))
plt.plot(freqs, trans_mr, label=f'MR Damper (I={i_mr}A)', linewidth=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Transmissibility')
plt.title('Transmissibility with MR Damper')
plt.legend()
plt.grid(True)
plt.show()

Key Notes & Next Steps #

1. Parameter Alignment: Make sure k and m match the values from your original linear model to compare results.
2. Current Flexibility: You can modify i_mr to be a time-varying function (e.g., i_func(t)) if your MR damper has changing current input—just update the nonlinear_system function to accept it.
3. Validation: To check if integration works, set the MR parameters to approximate a linear damper (e.g., tweak F0 to 0 and C1 to your original c value) and compare with your original linear plot.
4. Steady-State Check: The time vector t should be long enough for the system to reach steady state (we used 10s, adjust based on your system's natural frequency).

内容的提问来源于stack exchange，提问作者S.R