AnyMath 文档
Notebook

多质量系统扭转振动研究

本例将展示多质量机械系统扭转振动的建模和分析。 该模型演示了脉冲作用如何激发系统自身的模态,以及共振效应如何导致振动的显着放大。

扭转振动发生在具有多个惯性质量通过弹性轴连接的系统中。 每个这样的系统具有一组固有频率(模式),在该频率下发生振动的共振放大。

初始激励:1000Nm的脉冲时刻,持续时间为0.01s。

模型图:

使用具有弹性键的多质量体系来实现扭转振动模型。

torsional_oscillations--1759313337141.png

定义加载和运行模型的函数:

In [ ]:
m = "torsional_oscillations"
function start_model_engee(m)
    try
        engee.close(m, force=true) # 如果模型处于打开状态,则关闭该模型
    catch err # 以防没有模型关闭。
        # 模型没有打开
    end
    
    try
        m = engee.load("$((@__DIR__))/$m.engee") # 加载模型
        engee.run(m) # 启动模型
    catch err # 如果模型未加载
        m = engee.load("$((@__DIR__))/$m.engee") # 加载模型
        engee.run(m) # 启动模型
    end
end
Out[0]:
start_model_engee (generic function with 1 method)

定义用于分析的函数:

固有频率分析FFT函数:

In [ ]:
import Pkg; 
Pkg.add(["DSP", "Peaks"])
In [ ]:
using FFTW, DSP

function analyze_natural_frequencies(time_data, speed_data)
    # 为FFT准备数据
    dt = time_data[2] - time_data[1]
    fs = 1/dt  # 采样率
    
    # 应用窗口函数以获得更好的分辨率
    windowed_data = speed_data .* hanning(length(speed_data))
    
    # FFT计算
    fft_result = fft(windowed_data)
    
    # 频率轴
    frequencies = fftfreq(length(fft_result), fs)
    
    # 振幅谱(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)

固有频率峰值搜索功能:

In [ ]:
using Peaks

function find_natural_frequencies(frequencies, magnitude_db, min_height=20)
    # 搜索光谱正部分的峰
    pos_freq_idx = frequencies .> 0
    pos_freqs = frequencies[pos_freq_idx]
    pos_magnitudes = magnitude_db[pos_freq_idx]
    
    # 搜寻本地千里马
    peaks, _ = findmaxima(pos_magnitudes)
    
    # 高度过滤
    high_peaks = peaks[pos_magnitudes[peaks] .> min_height]
    
    natural_freqs = pos_freqs[high_peaks]
    amplitudes = pos_magnitudes[high_peaks]
    
    return natural_freqs, amplitudes
end
[ Info: Precompiling MakieExt [c3e534c0-8be5-5387-af52-275fbf2a048d] 
GKS: cannot open display - headless operation mode active
GKS: could not find font bold.ttf
WARNING: could not import Makie.documented_attributes into MakieExt
ERROR: LoadError: UndefVarError: `documented_attributes` not defined in `MakieExt`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
 [1] top-level scope
   @ /usr/local/ijulia-core/packages/MakieCore/hXATT/src/recipes.jl:286
 [2] include
   @ ./Base.jl:557 [inlined]
 [3] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::Nothing)
   @ Base ./loading.jl:2881
 [4] top-level scope
   @ stdin:6
in expression starting at /usr/local/ijulia-demos/packages/Peaks/frdrm/ext/MakieExt.jl:1
in expression starting at stdin:6
Error: Error during loading of extension MakieExt of Peaks, use `Base.retry_load_extensions()` to retry.
  exception =
   1-element ExceptionStack:
   Failed to precompile MakieExt [c3e534c0-8be5-5387-af52-275fbf2a048d] to "/user/.packages/compiled/v1.11/MakieExt/jl_L1Hxn7".
   Stacktrace:
     [1] error(s::String)
       @ Base ./error.jl:35
     [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool; flags::Cmd, cacheflags::Base.CacheFlags, reasons::Dict{String, Int64}, loadable_exts::Vector{Base.PkgId})
       @ Base ./loading.jl:3174
     [3] (::Base.var"#1110#1111"{Base.PkgId})()
       @ Base ./loading.jl:2579
     [4] mkpidlock(f::Base.var"#1110#1111"{Base.PkgId}, at::String, pid::Int32; kwopts::@Kwargs{stale_age::Int64, wait::Bool})
       @ FileWatching.Pidfile /usr/local/julia/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:95
     [5] #mkpidlock#6
       @ /usr/local/julia/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:90 [inlined]
     [6] trymkpidlock(::Function, ::Vararg{Any}; kwargs::@Kwargs{stale_age::Int64})
       @ FileWatching.Pidfile /usr/local/julia/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:116
     [7] #invokelatest#2
       @ ./essentials.jl:1057 [inlined]
     [8] invokelatest
       @ ./essentials.jl:1052 [inlined]
     [9] maybe_cachefile_lock(f::Base.var"#1110#1111"{Base.PkgId}, pkg::Base.PkgId, srcpath::String; stale_age::Int64)
       @ Base ./loading.jl:3698
    [10] maybe_cachefile_lock
       @ ./loading.jl:3695 [inlined]
    [11] _require(pkg::Base.PkgId, env::Nothing)
       @ Base ./loading.jl:2565
    [12] __require_prelocked(uuidkey::Base.PkgId, env::Nothing)
       @ Base ./loading.jl:2388
    [13] #invoke_in_world#3
       @ ./essentials.jl:1089 [inlined]
    [14] invoke_in_world
       @ ./essentials.jl:1086 [inlined]
    [15] _require_prelocked
       @ ./loading.jl:2375 [inlined]
    [16] _require_prelocked
       @ ./loading.jl:2374 [inlined]
    [17] run_extension_callbacks(extid::Base.ExtensionId)
       @ Base ./loading.jl:1544
    [18] run_extension_callbacks(pkgid::Base.PkgId)
       @ Base ./loading.jl:1576
    [19] run_package_callbacks(modkey::Base.PkgId)
       @ Base ./loading.jl:1396
    [20] __require_prelocked(uuidkey::Base.PkgId, env::String)
       @ Base ./loading.jl:2399
    [21] #invoke_in_world#3
       @ ./essentials.jl:1089 [inlined]
    [22] invoke_in_world
       @ ./essentials.jl:1086 [inlined]
    [23] _require_prelocked
       @ ./loading.jl:2375 [inlined]
    [24] macro expansion
       @ ./loading.jl:2314 [inlined]
    [25] macro expansion
       @ ./lock.jl:273 [inlined]
    [26] __require(into::Module, mod::Symbol)
       @ Base ./loading.jl:2271
    [27] #invoke_in_world#3
       @ ./essentials.jl:1089 [inlined]
    [28] invoke_in_world
       @ ./essentials.jl:1086 [inlined]
    [29] require(into::Module, mod::Symbol)
       @ Base ./loading.jl:2260
    [30] eval
       @ ./boot.jl:430 [inlined]
    [31] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
       @ Base ./loading.jl:2734
    [32] softscope_include_string(m::Module, code::String, filename::String)
       @ SoftGlobalScope /usr/local/ijulia-core/packages/SoftGlobalScope/u4UzH/src/SoftGlobalScope.jl:65
    [33] execute_request(socket::ZMQ.Socket, msg::IJulia.Msg)
       @ IJulia /usr/local/ijulia-core/packages/IJulia/czQmB/src/execute_request.jl:86
    [34] #invokelatest#2
       @ ./essentials.jl:1055 [inlined]
    [35] invokelatest
       @ ./essentials.jl:1052 [inlined]
    [36] eventloop(socket::ZMQ.Socket)
       @ IJulia /usr/local/ijulia-core/packages/IJulia/czQmB/src/eventloop.jl:8
    [37] (::IJulia.var"#15#17")()
       @ IJulia /usr/local/ijulia-core/packages/IJulia/czQmB/src/eventloop.jl:38
@ Base loading.jl:1550
Out[0]:
find_natural_frequencies (generic function with 2 methods)

使用脉冲效应开始仿真

In [ ]:
start_model_engee(m);

通过在仿真环境中运行模型,在信号可视化窗口中,您可以在频域中看到FFT(快速傅立叶变换)后的速度信号的图形。:

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

通过分析该图,可以检测系统的固有频率。

让我们在这个交互式脚本中构建一个类似的图表。 为此,请将仿真数据输出到变量。

输出数据进行分析:

In [ ]:
# 获得所有惯性的角速度
omega1 = collect(simout["torsional_oscillations/FFT按速度。1"]);

# 分配用于分析的时间和数据
time_data = omega1[:,1];
speed_data = omega1[:,2];  # 分析第一惯性

执行FFT分析:

In [ ]:
# 用于确定固有频率的FFT分析
frequencies, magnitude_db = analyze_natural_frequencies(time_data, speed_data);
natural_freqs, amplitudes = find_natural_frequencies(frequencies, magnitude_db);

# 找到的固有频率的输出
println("检测到的固有频率:")
for (i, (freq, amp)) in enumerate(zip(natural_freqs, amplitudes))
    println("$(i)-I模式:$(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)

所执行分析的可视化:

In [ ]:
using Plots
using Printf

function create_torsional_analysis_plots(time_data, speed_data, frequencies, magnitude_db)
    gr()
    
    # 将所有数据转换为实数
    real_frequencies = real.(frequencies)
    real_magnitude_db = real.(magnitude_db)
    real_time_data = real.(time_data)
    real_speed_data = real.(speed_data)
    
    # 时间回应时间表
    p1 = plot(real_time_data, real_speed_data, 
              title="脉冲后的扭转振动",
              xlabel="时间,从", 
              ylabel="角速度,rad/s",
              linewidth=2,
              linecolor=:blue)
    
    # FFT频谱图
    pos_idx = (real_frequencies .> 0) .& (real_frequencies .< 300)
    p2 = plot(real_frequencies[pos_idx], real_magnitude_db[pos_idx] .- 30,
              title="FFT固有频率分析", 
              xlabel="频率,赫兹",
              ylabel="振幅,dBm",
              xlim=(0, 200),
              linewidth=2,
              linecolor=:red)
    
    # 发现的固有频率的分配
    experimental_freqs = natural_freqs
    experimental_amps = amplitudes .- 30
    
    scatter!(p2, experimental_freqs, experimental_amps, 
            markersize=8, 
            markercolor=:green,
            markerstroke=:black,
            label="固有频率")
    
    # 频率签名
    for (i, (freq, amp)) in enumerate(zip(experimental_freqs, experimental_amps))
        annotate!(p2, freq+5, amp+3, text(@sprintf("%.1f赫兹", freq), 8, :black))
    end
    
    # 组合图形
    combined_plot = plot(p1, p2, layout=(2,1), size=(600,800))
    
    return combined_plot
end

# 创建图形
analysis_plots = create_torsional_analysis_plots(time_data, speed_data, frequencies, magnitude_db)
Out[0]:
No description has been provided for this image

固有频率下的共振现象

现在让我们检查作用扭矩源的不同频率如何影响第一惯性块的角速度。

扭矩的幅值为1Nm。

为此,让我们运行一个模型,其中正弦信号用作影响源,而不是脉冲发生器。

模型图:

torsional_oscillations_sin--1759313194565.png

使用变量**omega**,我们在每次仿真运行之前确定正弦信号的频率。

In [ ]:
n = "torsional_oscillations_sin"
Out[0]:
"torsional_oscillations_sin"
In [ ]:
omega = 10
start_model_engee(n)
w = collect(simout["△n/速度。1"]);
torque = collect(simout["$n/时刻"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="时间,从", ylabel="转速,rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="时间,从", ylabel="扭矩矩,N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
No description has been provided for this image
In [ ]:
omega = 20
start_model_engee(n)
w = collect(simout["△n/速度。1"]);
torque = collect(simout["$n/时刻"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="时间,从", ylabel="转速,rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="时间,从", ylabel="扭矩矩,N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
No description has been provided for this image
In [ ]:
omega = 30
start_model_engee(n)
w = collect(simout["△n/速度。1"]);
torque = collect(simout["$n/时刻"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="时间,从", ylabel="转速,rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="时间,从", ylabel="扭矩矩,N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
No description has been provided for this image
In [ ]:
omega = natural_freqs[1] # 第一固有频率
start_model_engee(n)
w = collect(simout["△n/速度。1"]);
torque = collect(simout["$n/时刻"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="时间,从", ylabel="转速,rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="时间,从", ylabel="扭矩矩,N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
No description has been provided for this image
In [ ]:
omega = 60
start_model_engee(n)
w = collect(simout["△n/速度。1"]);
torque = collect(simout["$n/时刻"]);
time_data1 = w[:,1];
speed_data1 = w[:,2];
time_torque = torque[:,1]
torque_data = torque[:,2]
p1 = plot(time_data1, speed_data1, xlabel="时间,从", ylabel="转速,rad/s", color="red")
p2 = plot(time_torque, torque_data, xlabel="时间,从", ylabel="扭矩矩,N*m")
plot(p1,p2, layout=(2,1))
Out[0]:
No description has been provided for this image

可以看出,与其他频率不同,暴露于固有频率会导致扭矩和旋转惯量的倍数增加。 观察到"跳动"现象,这往往伴随着共振。

结论:

在这个例子中,演示了一个用于分析多质量系统扭转振动的综合模型。

该模型可适用于各种扭转系统:涡轮发电机组,轧机驱动器,船用轴线和其他工业机构。