3. 扭振数学模型：单自由度扭振模型、多自由度扭振模型、模态分析基础

各位工程师朋友，咱们今天聊聊扭振的数学模型。说实话，我刚入行那会儿，觉得扭振这东西特别玄乎——不就是轴在转吗？后来在项目里吃过亏，才明白这玩意儿不搞清楚，传动链设计就是纸上谈兵。

扭振模型，说白了就是描述“轴怎么扭”的数学工具。你想想看，风机传动链那么长，从叶片到发电机，中间经过齿轮箱、联轴器，每个部件都在扭。它们不是铁板一块，而是有弹性的。弹性体就会振动，这就是扭振的根源。

3.1 单自由度扭振模型

先从最简单的说起。单自由度，就是只考虑一个转动惯量、一根轴。比如你只盯着发电机转子看，忽略其他部件。

它的运动方程长这样：

J·θ'' + C·θ' + K·θ = T(t)

其中：

• J — 转动惯量，单位 kg·m²。说白了就是“转起来有多费劲”。
• C — 阻尼系数，单位 N·m·s/rad。它消耗能量，让振动慢慢停下来。
• K — 扭转刚度，单位 N·m/rad。轴越硬，K越大。
• θ — 角位移，单位 rad。就是扭了多少角度。
• T(t) — 外部扭矩，随时间变化。

我在项目里遇到过一件事。有次做某型2MW风机的扭振分析，客户说“单自由度够用了”。结果算出来的共振频率和实测差了15%。后来发现，发电机转子的转动惯量确实占主导，但联轴器的柔性不可忽略。单自由度模型把联轴器当成刚性了，自然不准。

3.2 多自由度扭振模型

实际的风机传动链，怎么可能只有一个自由度？叶片、轮毂、低速轴、齿轮箱、高速轴、发电机……每个部件都有自己的惯量和刚度。它们串联在一起，就构成了多自由度系统。

多自由度模型的方程是矩阵形式：

[J]{θ''} + [C]{θ'} + [K]{θ} = {T(t)}

这里：

• [J] — 惯量矩阵，对角阵。每个对角元素对应一个部件的转动惯量。
• [C] — 阻尼矩阵。通常不是对角阵，因为阻尼会耦合。
• [K] — 刚度矩阵。非对角元素表示部件之间的弹性耦合。
• {θ} — 角位移向量。每个元素对应一个自由度的转角。

举个例子，一个三自由度系统（比如叶片-齿轮箱-发电机），刚度矩阵长这样：

[K] = [K1,  -K1,   0;
       -K1, K1+K2, -K2;
        0,   -K2,   K2]

你看，非对角元素 -K1 和 -K2 就是耦合项。它们把三个自由度“绑”在一起，一个地方扭了，其他地方跟着动。

多自由度模型的求解，通常用数值方法。比如Newmark-β法、Wilson-θ法。我自己习惯用Newmark-β，稳定性和精度都不错。代码实现也不复杂：

% Newmark-β法求解多自由度扭振
% 输入: J, C, K, T(t), 时间步长dt
% 输出: 角位移θ(t), 角速度θ'(t), 角加速度θ''(t)

function [theta, theta_dot, theta_ddot] = newmark_beta(J, C, K, T, dt, beta, gamma)
    n = length(J);          % 自由度数量
    N = length(T);          % 时间步数
    
    % 初始条件
    theta = zeros(n, N);
    theta_dot = zeros(n, N);
    theta_ddot = J \ (T(:,1) - C*theta_dot(:,1) - K*theta(:,1));
    
    % 有效刚度矩阵
    K_eff = K + (gamma/(beta*dt))*C + (1/(beta*dt^2))*J;
    
    for i = 2:N
        % 有效载荷
        F_eff = T(:,i) + J*(theta(:,i-1)/(beta*dt^2) + theta_dot(:,i-1)/(beta*dt)) ...
                + C*(theta(:,i-1)*gamma/(beta*dt) + theta_dot(:,i-1)*(gamma/beta - 1) ...
                + theta_ddot(:,i-1)*dt*(gamma/(2*beta) - 1));
        
        % 求解位移
        theta(:,i) = K_eff \ F_eff;
        
        % 更新速度和加速度
        theta_ddot(:,i) = (theta(:,i) - theta(:,i-1))/(beta*dt^2) ...
                         - theta_dot(:,i-1)/(beta*dt) ...
                         - theta_ddot(:,i-1)*(1/(2*beta) - 1);
        theta_dot(:,i) = theta_dot(:,i-1) + (1-gamma)*dt*theta_ddot(:,i-1) ...
                        + gamma*dt*theta_ddot(:,i);
    end
end

3.3 模态分析基础

模态分析，说白了就是找系统的“固有频率”和“振型”。为什么要找？因为当外部激励频率接近固有频率时，系统会发生共振，振幅急剧增大。在风机里，这可能导致齿轮断齿、联轴器损坏、甚至轴断裂。

模态分析从无阻尼自由振动方程开始：

[J]{θ''} + [K]{θ} = {0}

假设解的形式为 {θ} = {φ}·sin(ωt)，代入后得到特征方程：

([K] - ω²[J]){φ} = {0}

这是一个广义特征值问题。ω²就是特征值，对应固有频率的平方。{φ}是特征向量，对应振型。

举个例子，一个两自由度系统（比如低速轴和高速轴），惯量矩阵和刚度矩阵为：

[J] = [J1, 0; 0, J2]
[K] = [K1, -K1; -K1, K1+K2]

解特征方程，得到两个固有频率：

ω1² = 0  （刚体模态，整个系统一起转）
ω2² = K1/J1 + K2/J2 + K1/J2  （弹性模态，两轴相对扭转）

对应的振型：

• 第一阶：{φ1} = [1, 1]ᵀ，两个自由度同向等幅运动。
• 第二阶：{φ2} = [1, -J1/J2]ᵀ，两个自由度反向运动。

我曾经处理过一个案例。某3MW风机在并网瞬间，高速轴扭矩波动异常大。模态分析发现，第二阶固有频率（约85Hz）正好和电网的某次谐波频率接近。虽然没到完全共振，但已经放大了扭矩波动。后来我们调整了联轴器的刚度，把固有频率挪到了95Hz以上，问题就解决了。

3.4 知识体系总览

下面这张图，是我自己总结的扭振数学模型知识体系。从单自由度到多自由度，再到模态分析，层层递进。你把它印在脑子里，做扭振分析就有方向了。

嗯，这张图把今天的内容串起来了。单自由度是基础，多自由度是实战，模态分析是工具。三者缺一不可。

最后说一句心里话：模型再漂亮，最终要落到工程上。我见过太多人沉迷于复杂的数学模型，却忽略了最基本的物理假设。比如，齿轮箱的扭振模型里，齿轮啮合刚度到底取多少？轴承的阻尼怎么给？这些参数如果拍脑袋，再精密的算法也是白搭。

所以，我的习惯是：先做实验模态分析，拿到真实的固有频率和阻尼比，再去调模型参数。模型和实验对上了，才敢用它做载荷预测。这个习惯，让我少走了很多弯路。