Add Harmonic Balance Method

[jguarato]

Summary
This PR introduces Harmonic Balance Analysis to ROSS, allowing users to compute steady-state periodic responses of rotors under harmonic excitation using the Harmonic Balance Method (HBM).

Features

• New method: run_harmonic_balance_response()
• Computes steady-state harmonic responses directly in the frequency domain
• Supports multiple excitation nodes, magnitudes, phases, and harmonics
• Provides tools for post-processing: time-domain reconstruction, and frequency response plots

Notes

• The execution time using run_harmonic_balance_response() was significantly reduced compared to the traditional run_time_response() method — decreasing ~ 90 % for the same rotor and excitation conditions.
• This highlights the efficiency of the Harmonic Balance approach for steady-state analysis.

Example
Consider the following forces:

$$ F_{x,29}(t) = A_1 \cos(\omega t + p_1) + A_2 \cos(2\omega t + p_2) + A_3 \cos(3\omega t + p_3), $$

$$ F_{y,29}(t) = A_1 \sin(\omega t + p_1) + A_2 \sin(2\omega t + p_2) + A_3 \sin(3\omega t + p_3), $$

$$ F_{x,33}(t) = m \ e \ \omega^2 \cos(\omega t), $$

$$ F_{y,33}(t) = m \ e \ \omega^2 \sin(\omega t). $$

# Rotor operating speed (rad/s)
speed = 200.0

# Time vector for reconstruction
t = np.arange(0, 10, 1e-4)

# Harmonic excitation parameters
A1, A2, A3 = 1.0, 10.0, 5.0  # Force amplitudes (N)
p1, p2, p3 = 0.0, 0.0, 0.0  # Phase angles (rad)
m, e = 0.2, 0.01  # Unbalance mass (kg) and eccentricity (m)

# Define probe for observation
probe = rs.Probe(15, rs.Q_(45, "deg"))

# Run Harmonic Balance Analysis
results = rotor.run_harmonic_balance_response(
    speed=speed,
    t=t,
    harmonic_forces=[
        {
            "node": 29,
            "magnitudes": [A1, A2, A3],
            "phases": [p1, p2, p3],
            "harmonics": [1, 2, 3],
        },
        {"node": 33, "magnitudes": [m * e * speed**2], "phases": [0], "harmonics": [1]},
    ],
    n_harmonics=3,  # Number of harmonics considered in the HBM solution
)

probes = [rs.Probe(15, 0), rs.Probe(15, rs.Q_(90, "deg"))]

# Frequency-domain
fig1 = results.plot(probes)
fig2 = results.plot_deflected_shape()

# Time-domain reconstruction of the steady-state response
time_response = results.get_time_response()
fig3 = time_response.plot_1d(probes)

Results considering only unbalance response