Blade Structural Dynamics Flap/Lag/Torsion (12-state)

def _aero_forcing(t):
    """Per-revolution aerodynamic excitation at 1/rev and 2/rev."""
    psi = _OMEGA * t  # azimuth angle
    f = np.zeros(6)
    q_dyn = 0.5 * _RHO_AIR * (_OMEGA * 0.75 * _R_BLADE)**2 * _CHORD * _R_BLADE
    f[0] = _F_AERO_AMP[0] * q_dyn * np.sin(psi)           # 1/rev flap
    f[1] = _F_AERO_AMP[1] * q_dyn * np.sin(2.0 * psi)     # 2/rev flap
    f[2] = _F_AERO_AMP[2] * q_dyn * (1.0 + 0.3 * np.cos(psi))  # lag: steady + 1/rev
    f[3] = _F_AERO_AMP[3] * q_dyn * np.cos(2.0 * psi)     # 2/rev lag
    f[4] = _F_AERO_AMP[4] * q_dyn * np.cos(psi)            # 1/rev torsion
    f[5] = _F_AERO_AMP[5] * q_dyn * np.cos(3.0 * psi)      # 3/rev torsion
    return f

def rhs_blade_structural(t, y):
    """Rotating blade flap-lag-torsion dynamics with centrifugal stiffening."""
    q = y[:6]
    qd = y[6:12]

    # Effective stiffness with centrifugal stiffening
    omega_eff_sq = _OMEGA_NR**2 + _K_CF * _OMEGA**2

    # Stiffness force: -omega_eff^2 * q
    f_stiffness = -omega_eff_sq * q

    # Damping force: -2 * zeta * omega_eff * qdot
    omega_eff = np.sqrt(omega_eff_sq)
    f_damping = -2.0 * _ZETA * omega_eff * qd

    # Cross-coupling forces (Coriolis-like from rotation)
    f_coupling_stiffness = -_COUPLING_FL @ q * _OMEGA**2
    f_coupling_velocity = -_COUPLING_FT @ qd * _OMEGA

    # Aerodynamic forcing
    f_aero = _aero_forcing(t)

    # EOM: M * qdd = f_stiffness + f_damping + f_coupling + f_aero
    qdd = (f_stiffness + f_damping + f_coupling_stiffness
           + f_coupling_velocity + f_aero) / _M_GEN

    return np.concatenate([qd, qdd])