模态叠加法算法理论及其编程实现

M N ，因此有 qN
ΦNT Mx 。
五、模态叠加法步骤
1、求解广义特征值问题，求得特征值 r 和特征向量r 。
5
广义特征值问题为 KΦ MΦ ，其特征方程为： K M 0
在 MATLAB 中利用如下语句求得特征值和特征向量： [V,D]=eig(K,M); %特征向量矩阵和特征值矩阵 其中 D 为特征值矩阵（对角阵），V 为对应的特征向量矩阵。
Q(i,:)=deltat*conv(fN(i,:),h);
qNt(i,:)=exp(-ksi(i)*Wn(i)*t).*(qN0(i)*cos(Wd(i)*t)+(qNd0(i)+ksi(i)*Wn(i)*qN0(i))/Wd(i)*sin(Wd(i)*t)+Q(i,1:1/deltat+1));
fi
i 2
求得固有频率。
V1=zeros(n,n);
VN=zeros(n,n);
for i=1:n
V1(:,i)=V(:,i)/V(1,i); %求振型矩阵，各列除以各列的第一个元素
end
3、振型矩阵 Φ 关于主质量归一化，得到正则振型矩阵 ΦN 。
由前面论述，可利用Nr rr ， r
1 mr
1 rT Mr
0 ，因为 K j
j Mj
，两边左乘 iT
可得iT K j
T
ji
M
j
，因此有
iT
K
j
0。
由这两个式子可见：任意 2 个振型之间，既有对 M 的正交性，又有对 K 的正交性，它们统称为振型的正交性。
方法二：
设振系的第 i 个与第 j 个振型向量分别为i 和 j ，按照振型方程（19）有
Ki i Mi （3-1）
将其带入式（3-3）得正交关系
iT K j 0
3
由这两个式子可见：任意 2 个振型之间，既有对 M 的正交性，又有对 K 的正交性，它们统称为振型的正交性。
因此，由振型的正交性有：
Mp
ΦT
MΦ
12TT
M[1,2 ,,n ]
12TT
M1 M1
1T M2 2T M2
1T 2T
Mn Mn
对于系统的振动微分方程：
Mx Cx Kx f （1）
这里质量矩阵和刚度矩阵均为实对称矩阵，阻尼矩阵也为实对称矩阵。 由于此常微分方程组中各个微分方程相互耦合，无法单独求解。因此，要想办法通过坐标变换，将此微分方程组解 耦——即各个微分方程相互独立，可以单独求解。
令 x Φq 带入上式有：
MΦq CΦq KΦq f （2）
因此， r
1 mr
1 rT Mr
因此， ΦN [N1,N 2 ,,Nn ] [11, 22 ,, nn ]
此时，由式（14） KΦ MΦ 得 K p ΦNT KΦN ΦNT MΦN Λ Λ。
因此，对原系统振动微分方程做如下变换：
（1）令 x ΦN qN 带入上式有：
MΦN qN CΦN qN KΦN qN f
我们知道 x ΦN qN ，因此，利用 qN ΦN-1 x 似乎可以直接进行变换。但实际上，由于进行矩阵的求逆运算计算量
非常大，而且数值计算中会产生误差，因此，我们不采用此方法进行变换。
将其两端同时左乘
M
Φ -1 T
NN
M
，则有
M
Φ -1 T
NN
Mx
M
Φ -1 T
NN
MΦN qN
由于
M
-1 N
E ， ΦNT MΦN
（2）两边同乘以 ΦNT 得：
即：
ΦNT MΦN qN ΦNT CΦN qN ΦNT KΦN qN ΦNT f
4
qN C pqN ΛqN fNp
上式中 qN 称为正则坐标， fNp 为正则激励力。
四、解耦方程的求解
对于动力学方程 M pq C pq K pq f p ，当阻尼矩阵 C p 可对角化时，其分量形式为：
K j j M j （3-2）
用
T j
左乘式（3-1）有：
T j
Ki
i
T j
Mi
（3-3）
用iT 左乘式（3-2）有：
iT K j
T
ji
M
j
（3-4）
因为 M 和 K 都是对称矩阵，故将式（3-4）转置后得：
T j
Ki
j
T j
M
i
（3-5）
式（3-3）减式（3-4）得
(i
j
)
T j
M i
0
若 i j ，则有正交关系iT M j 0
50000*1.6*sin(57.5*t)+10*1.6*57.5*cos(57.5*t)];
fN=VN'*f; %fN 为正则激励力
Q=zeros(n,2*m-1); %动力学响应中的卷积积分部分
%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n
h=1/Wd(i)*sin(Wd(i)*t).*exp(-ksi(i)*Wn(i)*t);
qNi
(t)
e [q init Ni
(0)
cos di t
q Ni
(0)
iniqNi di
(0)
sin
dit
1 di
t 0
f ( )eini (t ) Ni
sin di (t
)d ]
MATLAB 中代码如下：
deltat=0.01;
%正则坐标系下的响应公式中卷积积分 t 的分段
t0=0;tf=1;
%结果为特征值对角阵，与理论相符合
CNp=VN'*C*VN;
%若 CNp 为对角阵，则原方程组可以解耦
4、将初始条件变换到正则坐标上。
qN ΦNT Mx
qN (0) ΦNT Mx (0)
MATLAB 中如下实现：
x=zeros(n,1);
%物理坐标系下位移向量
xd=zeros(n,1);
%物理坐标系下速度向量
q Ni
(0)
iniqNi di
(0)
sin
dit
1 di
t 0
f ( )eini (t ) Ni
sin di (t
)d ]
其中， di 1 i2ni
利用 x ΦN qN 可求得物理坐标下的响应。
注意，我们已知的初始条件是物理坐标下的，而解耦方程是模态坐标下的方程，因此，需要将物理坐标下的初始条件转 换到模态坐标下去。如何将初始条件变换到模态坐标（正则坐标）下？
求得正则振型矩阵 ΦN 。
Mp=V1'*M*V1;
%主质量矩阵
Kp=V1'*K*V1;
%主刚度矩阵
Cp=V1'*C*V1;
for i=1:n VN(:,i)=V1(:,i)/sqrt(Mp(i,i)); %求正则振型矩阵
end 验证程序：
MNp=VN'*M*VN;
%结果为 1 对角阵，与理论相符合
KNp=VN'*K*VN;
于是，原方程（1）的解耦问题最终归结为求解特征值问题（14），相应的特征值方程为
K M 0 （20）
一旦求出这个多项式的根 ，带回式（14），便可以确定变换矩阵 Φ 。
三、归一化——关于主质量归一化
2
习惯上我们让方程的第一项系数为 1，即 M p ΦT MΦ E
问题是，如何找到这样的 Φ ，使得 M p E ，我们令此时的 Φ 为 ΦN ，则有 M p ΦNT MΦN E
耦。
问题是，如何求得此 Φ ？ 二、坐标变换矩阵 Φ 的求解——求解广义特征值问题
1
将（4）式变为 ΦT M pΦ1M 1 ，带入（6）式得
M pΦ1M 1KΦ K p （10）
KΦ
MΦM
1 p
K
p
（11）
其中
M
p
和
K
p
均为对角阵，
M
1 p
K
p
可合并为一个对角阵：
其中对角线上元素为
M
1 p
K
p
两边同乘以 ΦT 得：
当方程可以解耦时，应有
ΦT MΦq ΦTCΦq ΦT KΦq ΦT f （3） ΦT MΦ M p （4）
ΦTCΦ C p （5）
ΦT KΦ K p （6）
式中 M p 、 C p 和 K p 均为如下形式的对角阵：
M
p
=
m1 0
0
m2
0 0
，
C
miqi ciqi kiqi f pi （ i 1, 2,……, n ）
两边同除以 mi 有：
qi
ci mi
q i
ki mi
qi
f pi mi
其中 ki mi
2 ni
，令
ci mi
2ini ，则有 i
ci 2mini
2
ci mi ki
（之所以这样做是仿效一元二次方程的形式） 上式化为
qi
2iniqi
模态叠加法算法理论及其编程实现
作者：谢胜龙
以下是本人硕士期间学习总结，供大家学习参考，如有错误，请多包涵。参考资料： 陈奎孚《机械振动基础》、《少齿差星轮型减速器弹性动力学分析》
一、系统动力学方程的解耦——坐标变换，将物理坐标 x 通过坐标变换 x Φq 转换到模态空间对系统原振动微分方程 进行解耦，转换为模态坐标 q 下的（用模态坐标表示的）系统振动微分方程。
1 0Leabharlann 02 00
（12）
0
0
n
i
ki mi
（i