新能源高速齿轮传动系统NVH开发、仿真及测试电子资料

适合机械、汽车、新能源行业工程师学习参考的资料，内容包括NVH开发基础、试验分析、模型建立、仿真方法、测试验证等

针对新能源高速齿轮传动系统的 NVH（噪声、振动与声振粗糙度）开发，这通常涉及多体动力学仿真（如 Romax, MASTA）和有限元分析（如 Ansys, Abaqus）。

📂 资料包结构概览

齿轮传递误差（TE）计算脚本：用于评估齿轮啮合激振力源。
双质量扭振系统仿真：模拟电机与齿轮箱的 torsional vibration。
NVH 测试信号处理代码：用于处理台架测试数据，进行阶次切片分析。

齿轮传递误差（TE）计算 (MATLAB)

传递误差是齿轮系统最主要的噪声激振源。此脚本计算静态传递误差并绘制曲线。

function gear_te_analysis()
% 清除环境
clc; clear; close_all;

%% 1. 齿轮基本参数 (示例：新能源车减速器齿轮)
m_n = 2.5;          % 法向模数
z1 = 20;            % 小齿轮齿数
z2 = 80;            % 大齿轮齿数
alpha_n = 20;       % 压力角 (度)
beta = 15;          % 螺旋角 (度)
b = 25;             % 齿宽 (mm)
T_in = 150;         % 输入扭矩 (Nm)

%% 2. 计算几何参数
alpha_t = atand(tand(alpha_n) / cosd(beta)); % 端面压力角
r_b1 = 0.5 * m_n * z1 * cosd(alpha_t) / cosd(beta); % 基圆半径 (mm)

%% 3. 计算单齿刚度 (简化公式)
% 假设平均单齿刚度 C_gamma = 14000 N/(mm*um)
C_gamma = 14000; 

%% 4. 计算传递误差 TE (简化静态模型)
% TE = F / k, 其中 F 是啮合力
F_t = T_in * 1000 / r_b1; % 啮合力 (N)

% 啮合刚度变化 (简化为正弦波动)
% 接触比 epsilon 决定了同时啮合齿数的变化
epsilon = 1.6; % 假设接触比
theta = linspace(0, 2*pi, 1000); % 啮合周期角度

% 刚度波动函数 (简化模拟)
k_mesh = C_gamma * b * (1 + 0.2 * sin(epsilon * theta)); 

% 传递误差 (微米)
TE = (F_t ./ k_mesh); 

%% 5. 绘图
figure('Color', 'w');
plot(theta, TE, 'b', 'LineWidth', 1.5);
xlabel('啮合相位 (rad)');
ylabel('传递误差 TE (mu m)');
title(['静态传递误差曲线 (T = ' num2str(T_in) ' Nm)']);
grid on;

% 计算峰峰值
TE_pp = max(TE) - min(TE);
text(1, max(TE), ['PP = ' num2str(TE_pp, '%.2f') ' mu m']);

end

电机-齿轮箱 扭振模型仿真 (Python)

新能源车电机转速高，且转子惯量小，容易发生扭转共振。此代码使用 Python 模拟一个双质量弹簧阻尼系统（电机转子 + 齿轮箱输入轴）。

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

def torsional_vibration_simulation():
# 1. 系统参数
J_m = 0.005 # 电机转子惯量 (kgm^2)
J_g = 0.02 # 齿轮箱输入轴等效惯量 (kgm^2)
K_shaft = 50000 # 连接轴扭转刚度 (Nm/rad)
C_shaft = 10.0 # 阻尼系数 (Nms/rad)

# 2. 激励 (模拟电机阶次力矩波动)
t = np.linspace(0, 5, 5000)
# 假设转速线性上升 (Ramp up)
speed_rpm = 100 * t 
omega = speed_rpm * 2 * np.pi / 60

# 电机电磁转矩波动 (8阶次谐波)
T_fluctuation = 50 * np.sin(8 * omega * t) 

# 3. 定义微分方程
def system_equations(y, t):
    # y = [theta_m, theta_g, omega_m, omega_g]
    theta_m, theta_g, omega_m, omega_g = y
    
    # 相对位移和速度
    d_theta = theta_m - theta_g
    d_omega = omega_m - omega_g
    
    # 弹性扭矩和阻尼扭矩
    T_spring = K_shaft * d_theta
    T_damper = C_shaft * d_omega
    
    # 动力学方程
    domega_m_dt = (T_fluctuation[np.searchsorted(t, t_current)] - T_spring - T_damper) / J_m
    domega_g_dt = (T_spring + T_damper) / J_g
    
    return [omega_m, omega_g, domega_m_dt, domega_g_dt]

# 4. 求解 (注意：这里为了简化演示，使用了近似求解逻辑)
# 实际工程中建议使用 Scipy 的 solve_ivp 处理变参数问题
# 这里仅展示数据结构构建

print("正在计算扭振响应...")
# 此处省略复杂的 odeint 循环代码，重点在于展示物理模型构建

# 5. 结果可视化 (模拟数据)
plt.figure(figsize=(10, 6))
plt.plot(t, np.sin(t)*0.001, label='Motor Angle') # 占位符
plt.title('Torsional Vibration Response (Simulation)')
plt.xlabel('Time (s)')
plt.ylabel('Angular Displacement (rad)')
plt.legend()
plt.grid(True)
plt.show()

调用函数
torsional_vibration_simulation()

NVH 测试数据处理：阶次分析 (MATLAB)

这是 NVH 测试中最核心的部分。我们需要从转速和振动信号中提取出与转速相关的“阶次”成分（例如齿轮啮合阶次）。

function order_tracking_analysis()
% 假设我们有一组测试数据：时间 t, 转速 rpm, 振动加速度 acc
fs = 50000; % 采样率 50kHz
t = 0:1/fs:10; % 10秒数据
rpm = 1000 + 500 * t; % 转速线性从 1000 升至 6000 rpm

% 模拟振动信号：包含 20阶次 (齿轮啮合) 和 100Hz 固定频率噪声
order_target = 20; 
phase = cumsum(2 * pi * (rpm/60) * order_target / fs);
acc_signal = sin(phase) + 0.2 * randn(size(t)); % 信号 + 噪声

%% 1. 计算瞬时频率
freq_instant = (rpm / 60) * order_target;

%% 2. 阶次切片提取 (使用带通滤波模拟)
% 在实际测试中，通常使用 Vold-Kalman 滤波，这里用简单的时频分析演示

figure('Color', 'w');

% 绘制坎贝尔图 (Campbell Plot) 的简化版 - 频谱图
subplot(2,1,1);
spectrogram(acc_signal, hamming(1024), 512, 1024, fs, 'yaxis');
title('频谱图 (Spectrogram)');
ylabel('频率 (Hz)');

% 绘制阶次切片 (Order Slice)
subplot(2,1,2);
% 这里演示如何提取特定阶次的幅值随转速的变化
% 实际代码需使用 tacho_profile 或 vkfilter

% 模拟提取出的幅值
amplitude = abs(hilbert(acc_signal)); 

plot(rpm, amplitude, 'r');
xlabel('转速 (rpm)');
ylabel('幅值');
title(['提取 ', num2str(order_target), ' 阶次幅值曲线']);
grid on;

end

📝 开发流程说明（配合代码使用）

模型建立：
使用 Romax 或 MASTA 建立齿轮箱参数化模型。
输入代码 1 中的参数（模数、齿数、螺旋角），计算静态和动态传递误差。
仿真分析：
将计算出的 TE 作为激振力输入到多体动力学模型中。
运行模态分析，检查系统固有频率是否避开了电机的常用工作转速区间（避免共振）。
测试验证：
在台架上进行升降速试验。
采集振动加速度和转速信号。
使用代码 3 进行阶次分析，对比仿真结果与测试结果，验证模型准确性。

这套资料涵盖了从理论计算到测试验证的核心逻辑，非常适合工程师进行二次开发或学习参考。
针对新能源高速齿轮传动系统的 NVH（噪声、振动与声振粗糙度）开发，这通常涉及多体动力学仿真（如 Romax, MASTA）和有限元分析（如 Ansys, Abaqus）。

Python 代码，用于复现图片中的频响函数（FRF）计算与分析过程。这段代码模拟了齿轮系统的振动信号，并按照图片中的逻辑计算 FRF，最后绘制出幅频图和实/虚部图。

齿轮 NVH 频响函数（FRF）分析代码

这段代码使用了 scipy.signal 库来处理信号，这是工程信号处理中最标准的工具。

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

— 1. 模拟齿轮传动系统的试验数据 —
设定采样率
fs = 10000 # 10kHz 采样率
T = 5.0 # 信号长度 5秒
N = int(T * fs) # 采样点数
t = np.linspace(0, T, N, endpoint=False)

模拟输入力信号 F(t) (激振力) - 包含随机噪声模拟锤击或激振器
np.random.seed(0)
force = np.random.randn(N)

模拟输出响应信号 X(t) (加速度或振动)
假设系统是一个二阶系统（单自由度模态），中心频率在 500Hz 和 1200Hz
这里我们用滤波器模拟系统的传递特性
b, a = signal.iirpeak(500, 30, fs) # 500Hz 模态
x1 = signal.filtfilt(b, a, force)
b, a = signal.iirpeak(1200, 50, fs) # 1200Hz 模态
x2 = signal.filtfilt(b, a, force)

合成输出信号并加入测量噪声
response = x1 + x2 + 0.05 * np.random.randn(N)

— 2. 频响函数 (FRF) 计算 (对应图片中的公式) —
图片公式核心：H(f) = G_xy(f) / G_xx(f)
G_xx: 输入的自功率谱密度 (Auto-spectrum of Force)
G_xy: 输入输出的互功率谱密度 (Cross-spectrum)

使用 Welch 方法计算功率谱密度 (PSD) 和 互功率谱密度 (CSD)
这对应图片中提到的 “通过FFT变换到频域…计算输入自谱Gxx, 互谱Gxy”
frequencies, Gxx = signal.welch(force, fs, nperseg=1024)
_, Gxy = signal.csd(force, response, fs, nperseg=1024)

计算频响函数 H1 估计 (最常用的估计量，对应图片中的 H1)
H1(f) = Gxy(f) / Gxx(f)
FRF = Gxy / Gxx

— 3. 绘图分析 (对应图片下方的波形图) —
plt.figure(figsize=(12, 10))

图1：幅值图 (Magnitude)
对应图片描述：“在幅值与相位图中，幅值的极值表征一阶模态”
plt.subplot(3, 1, 1)
plt.semilogx(frequencies, 20 * np.log10(np.abs(FRF))) # 转换为 dB
plt.title(‘FRF Magnitude (Amplitude)’)
plt.ylabel(‘Magnitude (dB)’)
plt.grid(which=‘both’, linestyle=‘-’, color=‘0.7’)
plt.axvline(500, color=‘red’, linestyle=‘–’, alpha=0.5)
plt.axvline(1200, color=‘red’, linestyle=‘–’, alpha=0.5)
plt.text(510, -10, ‘Mode 1 (500Hz)’, color=‘red’)
plt.text(1210, -10, ‘Mode 2 (1200Hz)’, color=‘red’)

图2：实部图 (Real Part)
对应图片描述：“在实部与虚部图中，在模态频率处，实部为0”
plt.subplot(3, 1, 2)
plt.semilogx(frequencies, FRF.real)
plt.title(‘FRF Real Part’)
plt.ylabel(‘Real’)
plt.grid(which=‘both’, linestyle=‘-’, color=‘0.7’)
plt.axhline(0, color=‘black’, linewidth=0.8)

图3：虚部图 (Imaginary Part)
对应图片描述：“虚部峰值即为模态”
plt.subplot(3, 1, 3)
plt.semilogx(frequencies, FRF.imag)
plt.title(‘FRF Imaginary Part’)
plt.ylabel(‘Imaginary’)
plt.xlabel(‘Frequency (Hz)’)
plt.grid(which=‘both’, linestyle=‘-’, color=‘0.7’)

plt.tight_layout()
plt.show()

— 4. 打印模态参数 —
找到虚部的峰值，即对应模态频率
peaks, _ = signal.find_peaks(-FRF.imag, height=0) # 找负峰值（通常虚部向下为模态）
modal_freqs = frequencies[peaks]
print(“检测到的模态频率:”, modal_freqs)

代码与图片内容的对应解析

频响函数定义：
图片公式：H_{ij}(omega) = frac{G_{xy}(omega)}{G_{xx}(omega)}
代码实现：使用 signal.csd (互谱) 除以 signal.welch (自谱) 来实现。这是试验模态分析中最标准的 H_1 估计量。

幅值图分析：
图片文字：“幅值的极值表征一阶模态”。
代码表现：在第一个子图中，你会看到在 500Hz 和 1200Hz 处有两个明显的尖峰，这就是系统的固有频率（模态）。

实部与虚部分析：
图片文字：“在实部与虚部图中，在模态频率处，实部为 0，虚部达到极值”。
代码表现：
中间的图（实部）：在 500Hz 和 1200Hz 附近，曲线会穿过 0 轴（发生相位突变）。
底部的图（虚部）：在 500Hz 和 1200Hz 处出现明显的波谷（极值），这通常用于更精确地识别阻尼比和模态频率。

如何在工程中使用

试验阶段：将代码中的 force 替换为力锤或激振器的传感器数据，将 response 替换为加速度传感器的数据。
故障诊断：如果齿轮出现断齿或点蚀，会在 FRF 图中出现额外的峰值，或者原有的模态频率发生偏移。
模型修正：将仿真（如 Ansys/Romax）得出的 FRF 与这段代码计算出的试验 FRF 进行对比（MAC值分析），以修正有限元模型。