Engee documentation
Notebook

Investigation of torsional oscillations of a multi-mass system

This example will show the modeling and analysis of torsional vibrations of a multi-mass mechanical system. The model demonstrates how an impulse action excites the system's own modes, and how a resonant effect leads to a significant amplification of vibrations.

Torsional vibrations occur in systems with several inertial masses connected by elastic shafts. Each such system has a set of natural frequencies (modes) at which resonant amplification of vibrations occurs.

Initial excitation: a pulse moment of 1000 Nm with a duration of 0.01 s.

Model diagram:

To implement the torsional vibration model, a multi-mass system with elastic bonds is used.

torsional_oscillations--1759313337141.png

Defining the function to load and run the model:

In [ ]:
m = "torsional_oscillations"
function start_model_engee(m)
    try
        engee.close(m, force=true) # closing the model if it is open
    catch err # in case there is no model to close.
        # The model was not opened
    end
    
    try
        m = engee.load("$((@__DIR__))/$m.engee") # loading the model
        engee.run(m) # launching the model
    catch err # in case the model is not loaded
        m = engee.load("$((@__DIR__))/$m.engee") # loading the model
        engee.run(m) # launching the model
    end
end
Out[0]:
start_model_engee (generic function with 2 methods)

Defining functions for analysis:

Natural frequency analysis FFT function:

In [ ]:
using FFTW, DSP

function analyze_natural_frequencies(time_data, speed_data)
    # Preparing data for FFT
    dt = time_data[2] - time_data[1]
    fs = 1/dt  # Sampling rate
    
    # Using the window function for better resolution
    windowed_data = speed_data .* hanning(length(speed_data))
    
    # FFT Calculation
    fft_result = fft(windowed_data)
    
    # Frequency axis
    frequencies = fftfreq(length(fft_result), fs)
    
    # Amplitude spectrum in dB
    magnitude_db = 20 * log10.(abs.(fft_result) .+ 1e-10)
    
    return frequencies, magnitude_db
end
Out[0]:
analyze_natural_frequencies (generic function with 1 method)

Natural frequency peak search function:

In [ ]:
using Peaks

function find_natural_frequencies(frequencies, magnitude_db, min_height=20)
    # Search for peaks in the positive part of the spectrum
    pos_freq_idx = frequencies .> 0
    pos_freqs = frequencies[pos_freq_idx]
    pos_magnitudes = magnitude_db[pos_freq_idx]
    
    # Search for local maxima
    peaks, _ = findmaxima(pos_magnitudes)
    
    # Height filtering
    high_peaks = peaks[pos_magnitudes[peaks] .> min_height]
    
    natural_freqs = pos_freqs[high_peaks]
    amplitudes = pos_magnitudes[high_peaks]
    
    return natural_freqs, amplitudes
end
Out[0]:
find_natural_frequencies (generic function with 2 methods)

Starting a simulation with a pulse effect

In [ ]:
start_model_engee(m);

By running the model in the simulation environment, in the Signal Visualization window, you can see a graph of the velocity signal after the FFT (fast Fourier transform) in the frequency domain.:

newplot - 2025-10-01T115035.780.png

By analyzing this graph, it is possible to detect the natural frequencies of the system.

Let's build a similar graph in this interactive script. To do this, we need to output simulation data to variables.

Output of data for analysis:

In [ ]:
# Obtaining the angular velocities of all inertia
omega1 = collect(simout["torsional_oscillations/FFT by speed.1"]);

# Allocating time and data for analysis
time_data = omega1[:,1];
speed_data = omega1[:,2];  # Analyzing the first inertia

Performing FFT analysis:

In [ ]:
# FFT analysis for determining natural frequencies
frequencies, magnitude_db = analyze_natural_frequencies(time_data, speed_data);
natural_freqs, amplitudes = find_natural_frequencies(frequencies, magnitude_db);

# Output of the found natural frequencies
println("Detected natural frequencies:")
for (i, (freq, amp)) in enumerate(zip(natural_freqs, amplitudes))
    println("$(i)-I mode: $(round(freq, digits=2)) Hz ($(round(amp.-30, digits=2)) dBm)")
end
Обнаруженные собственные частоты:
1-я мода: 54.16 Гц (63.55 dBm)
2-я мода: 102.51 Гц (25.96 dBm)
3-я мода: 141.2 Гц (30.16 dBm)
4-я мода: 165.7 Гц (11.16 dBm)

Visualization of the analysis performed:

In [ ]:
using Plots
using Printf

function create_torsional_analysis_plots(time_data, speed_data, frequencies, magnitude_db)
    gr()
    
    # Convert all data to real numbers
    real_frequencies = real.(frequencies)
    real_magnitude_db = real.(magnitude_db)
    real_time_data = real.(time_data)
    real_speed_data = real.(speed_data)
    
    # Time response schedule
    p1 = plot(real_time_data, real_speed_data, 
              title="Torsional vibrations after an impulse",
              xlabel="Time, from", 
              ylabel="Angular velocity, rad/s",
              linewidth=2,
              linecolor=:blue)
    
    # FFT Spectrum graph
    pos_idx = (real_frequencies .> 0) .& (real_frequencies .< 300)
    p2 = plot(real_frequencies[pos_idx], real_magnitude_db[pos_idx] .- 30,
              title="FFT natural frequency analysis", 
              xlabel="Frequency, Hz",
              ylabel="Amplitude, dBm",
              xlim=(0, 200),
              linewidth=2,
              linecolor=:red)
    
    # Allocation of found natural frequencies
    experimental_freqs = natural_freqs
    experimental_amps = amplitudes .- 30
    
    scatter!(p2, experimental_freqs, experimental_amps, 
            markersize=8, 
            markercolor=:green,
            markerstroke=:black,
            label="Natural frequencies")
    
    # Frequency signatures
    for (i, (freq, amp)) in enumerate(zip(experimental_freqs, experimental_amps))
        annotate!(p2, freq+5, amp+3, text(@sprintf("%.1f Hz", freq), 8, :black))
    end
    
    # Combining graphs
    combined_plot = plot(p1, p2, layout=(2,1), size=(600,800))
    
    return combined_plot
end

# Creating graphs
analysis_plots = create_torsional_analysis_plots(time_data, speed_data, frequencies, magnitude_db)
Out[0]:

The phenomenon of resonance at natural frequency

Now let's check how the different frequencies of the acting torque source affect the angular velocity of the first inertia block.

The amplitude of the torque is 1 Nm.

To do this, let's run a model where a sinusoidal signal is used as the source of the effect, not a pulse generator.

Model diagram:

torsional_oscillations_sin--1759313194565.png

Using the variable omega, we determine the frequency of the sinusoidal signal before each simulation run.

In [ ]:
n = "torsional_oscillations_sin"
Out[0]:
"torsional_oscillations_sin"
In [ ]:
omega = 10
start_model_engee(n)
w = collect(simout["$n/Speed.1"]);
torque = collect(simout["$n/Moment"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Time, from", ylabel="Rotation speed, rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="Time, from", ylabel="Torque moment, N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
In [ ]:
omega = 20
start_model_engee(n)
w = collect(simout["$n/Speed.1"]);
torque = collect(simout["$n/Moment"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Time, from", ylabel="Rotation speed, rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="Time, from", ylabel="Torque moment, N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
In [ ]:
omega = 30
start_model_engee(n)
w = collect(simout["$n/Speed.1"]);
torque = collect(simout["$n/Moment"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Time, from", ylabel="Rotation speed, rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="Time, from", ylabel="Torque moment, N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
In [ ]:
omega = natural_freqs[1] # The first natural frequency
start_model_engee(n)
w = collect(simout["$n/Speed.1"]);
torque = collect(simout["$n/Moment"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Time, from", ylabel="Rotation speed, rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="Time, from", ylabel="Torque moment, N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
In [ ]:
omega = 60
start_model_engee(n)
w = collect(simout["$n/Speed.1"]);
torque = collect(simout["$n/Moment"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="Time, from", ylabel="Rotation speed, rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="Time, from", ylabel="Torque moment, N*m")
plot(p1,p2, layout=(2,1))
Out[0]:

It can be seen that, unlike other frequencies, exposure to natural frequency causes a multiple increase in torque and rotational inertia. The phenomenon of "beating" is observed, which is often accompanied by resonance.

Conclusions:

In this example, a comprehensive model for analyzing torsional vibrations of a multi-mass system was demonstrated.

The model can be adapted for various torsion systems: turbo generator sets, rolling mill drives, marine shaft lines and other industrial mechanisms.