在上篇文章《理论分享 | 模态分析（2）：理解各阶模态对结构振动的真实贡献》中，介绍了动力学分析中的基本概念，包括模态解耦原理、模态参与因子以及模态有效质量等关键参数，并详细说明了相关计算公式。特别地，针对不含旋转自由度的实体单元，提到旋转方向的模态参与因子是基于虚拟旋转位移场（单位旋转角度）进行计算的，并提到在完备振型条件下，旋转方向各阶模态有效质量之和等于结构绕该轴的转动惯量这一结论。本文在此基础上通过理论推导来说明，并将编程实现的各方向模态参与因子与模态有效质量计算结果与商业软件进行对比验证。

1、旋转方向的模态参与因子和模态有效质量

1.1 虚拟旋转位移场的构造原理

模态参与因子是“给定激励向量”在模态坐标系下的投影，故要计算旋转方向的模态参与因子，首先要构造出旋转方向的“单位激励向量”，也就是让结构旋转单位角度（弧度制）。以绕X轴旋转为例，如上图所示，绕着X轴旋转1弧度，各个节点的位移可以表示为$[0,-z_i,y_i]$，即与节点本身的坐标相关，故最终整个结构绕X轴旋转方法激励向量可以表示为：

同理，绕Y轴和Z轴旋转的单位激励向量可以类似构造。

1.2 理论推导：旋转方向模态有效质量之和与转动惯量的关系：

以下给出了ROTX旋转方向模态有效质量之和等于绕X轴转动惯量的推导过程，其余两个方向是类似的：

与平动方向的模态有效质量之和等于结构的质量类似。从物理意义角度解读如下：结构发生单位平动位移或单位旋转角度的总能量守恒，平动或转动的总能量（${\gamma }^T\mathbf{M}\gamma$或${\theta }_x^T\mathbf{M}{\theta }_x$，用单位位移替代速度）会等于投影坐标系下的能量总和（${\sum }_i({\Phi }_i\mathbf{M}\gamma {)}^2$或${\sum }_i({\Phi }_i\mathbf{M}{\theta }_x{)}^2$）,类似于平面曲线运动中，结构的总动能等于X向、Y向动能之和。

2、数值计算与软件验证

作为理论分享的验证，针对上篇文章《理论分享 | 模态分析（2）：理解各阶模态对结构振动的真实贡献》中提到的悬臂梁案例进行了编程实现，组装总体刚度矩阵和质量矩阵进行模态分析，为确保分析过程的一致性，相应单元刚度矩阵和单元质量矩阵是从ANSYS中导出的，并自行组装成总体矩阵，对比结果如下：

简洁起见，上述仅给出了前十阶固有频率，X向、ROTX向模态参与因子和模态有效质量的对比。可以看出，除了数值误差导致的部分接近0的数（小于1E-10）存在差异，其它完全一致。

关键部分代码如下：

%% 模态求解
nEig = 10;  % 求解10阶
matFun = @(x) dK\(M(free, free)*x);
warning off;[eivecs(:, :), D(:, :)] = eigs(matFun, length(free), nEig+4, 'la', ...
    'Tolerance', 1e-12, 'MaxIterations', 200); % 10阶特征值计算
freqModal = zeros(ndof, nEig);
gama_num = zeros(nEig, 1);
[mu, ii] = sort(diag(D(:, :)), 'descend'); % 圆频率平方倒数
eivSort = eivecs(:, ii(1:nEig));
gama_num(:) = mu(1:nEig);
freqModal(free, :) = eivSort./sqrt(diag(eivSort'*M(free, free)*eivSort)'); % 振型向量归一化
lambda = 1./mu(1:nEig); % 圆频率
frequency = sqrt(lambda)/2/pi; % 固有频率

%% 模态参与因子与模态有效质量计算
% X向
lambda1 = zeros(ndof,1);lambda1(1:3:end) = 1; % X向单位激励向量
xParticipationFactor = freqModal'*M*lambda1; % X向模态参与因子
xEffectiveMass = xParticipationFactor.^2; % X向模态有效质量
% Y向
lambda2 = zeros(ndof,1);lambda2(2:3:end) = 1; % Y向单位激励向量
yParticipationFactor = freqModal'*M*lambda2; % Y向模态参与因子
yEffectiveMass = yParticipationFactor.^2;% Y向模态有效质量
% Z向
lambda3 = zeros(ndof,1);lambda3(3:3:end) = 1; % Z向单位激励向量
zParticipationFactor = freqModal'*M*lambda3; % Z向模态参与因子
zEffectiveMass = zParticipationFactor.^2;% Z向模态有效质量
% ROTX向
lambda4 = zeros(ndof,1); % ROTX向虚拟旋转位移场
lambda4(2:3:end) = -node_xyz_all(:,3);
lambda4(3:3:end) = node_xyz_all(:,2);
RotxParticipationFactor = freqModal'*M*lambda4; % ROTX向模态参与因子
RotxEffectiveMass = RotxParticipationFactor.^2;% ROTX向模态有效质量
% ROTY向
lambda5 = zeros(ndof,1); % ROTY向虚拟旋转位移场
lambda5(1:3:end) = node_xyz_all(:,3);
lambda5(3:3:end) = -node_xyz_all(:,1);
RotyParticipationFactor = freqModal'*M*lambda5; % ROTY向模态参与因子
RotyEffectiveMass = RotyParticipationFactor.^2;% ROTY向模态有效质量
% ROTZ向
lambda6 = zeros(ndof,1); % ROTZ向虚拟旋转位移场
lambda6(1:3:end) = -node_xyz_all(:,2);
lambda6(2:3:end) = node_xyz_all(:,1);
RotzParticipationFactor = freqModal'*M*lambda6; % ROTZ向模态参与因子
RotzEffectiveMass = RotzParticipationFactor.^2; % ROTZ向模态有效质量

3、总结

• 本文进一步对旋转方向各阶模态有效质量之和等于结构绕该轴的转动惯量进行了推导证明，并从能量守恒的角度尝试进行了深层含义的解读；
• 通过Matlab具体编程实现了最近两篇文章中分享的模态参与因子、模态有效质量的相关计算，并与软件ANSYS进行了对比验证。

最后

有需要上述完整代码的读者请联系作者！

ps：验证过程中，发现对于相同的一个单元（均为SOLID185单元，节点编号完全一致），ANSYS经典界面（APDL）导出的单元刚度矩阵与Workbench导出的结果始终存在一定差异。找了很久发现两者对于SOLID185的默认计算方法不同，ANSYS经典界面（APDL）中采用的是全积分，而workbench采用的是简化增强应变计算，可以通过SOLID185单元的KEYOPT(2)进行设置，具体介绍如下图所示：

-完-
公众号已同步更新：仿真优化工匠
欢迎联系！