仿真干货丨模态叠加法的基本原理及其实现方法（附Matlab代码）
将物理坐标x通过坐标变换x=Φq转换到模态空间对系统原振动微分方程进行解耦，转换为模态坐标q下的(用模态坐标表示的)系统振动微分方程。
对于系统的振动微分方程：
这里质量矩阵和刚度矩阵均为实对称矩阵，阻尼矩阵也为实对称矩阵。
由于此常微分方程组中各个微分方程相互耦合，无法单独求解。因此，要想办法通过坐标变换，将此微分方程组解耦——即各个微分方程相互独立，可以单独求解。
令x=Φq代入上式有：
当方程可以解耦时，应有：
式中Mp、Cp和Kp均为如下形式的对角阵：
即：
其中：
公式(8)为n阶微分方程组，各个坐标qi(i=1，2，...，n)之间相互独立，可以单独求解，即实现了原方程组的解耦。
问题是，如何求得此Φ?
将(4)式变为：
代入(6)式得：
其中Mp和Kp均为对角阵，从而有：
其中对角线上元素为：
于是式(11)变为：
式(14)为矩阵理论中的广义特征值问题。
如果将矩阵Φ按列分开：
其中：
由分块矩阵乘法：
则式(14)就变为如下熟悉的特征值形式：
之所以叫“广义”是因为M的存在。如果矩阵M为单位阵，那么就是常规的特征值问题。
于是，原方程(1)的解耦问题最终归结为求解特征值问题(14)，相应的特征值方程为：
一旦求出这个多项式的根λ，代回式(14)，便可以确定变换矩阵Φ 。
习惯上我们让方程的第一项系数为1，即：
问题是，如何找到这样的Φ，使得Mp=E，我们令此时的Φ为ΦN ，则有：
问题是ΦN如何求?
首先引入振型的正交性
方法一：
由前面的论述
可见，若λi≠λj，则有：
因为：
可得：
因此有：
由这两个式子可见：任意2个振型之间，既有对M的正交性，又有对K的正交性，它们统称为振型的正交性。
方法二：
设振系的第i个与第j个振型向量分别为Φi和Φj，按照振型方程(19)有：
可得：
因为M和K都是对称矩阵，故将上式转置后得：
所以有：
若λi≠λj，则有正交关系：
同样有：
由这两个式子可见：任意2个振型之间，既有对M的正交性，又有对K的正交性，它们统称为振型的正交性。因此，由振型的正交性有：
即
故对于ΦN有：
因此，必须有：
令：
代入上式有：
因此：
所以：
此时，由式(14)式得：
因此，对原系统振动微分方程做如下变换：
(1) 令x=ΦNqN代入上式有：
(2) 可得：
即：
上式中qN称为正则坐标，fNp为正则激励力。
对于动力学方程
当阻尼矩阵Cp可对角化时，其分量形式为：
两边同除以mi有：
其中：
则有：
(之所以这样做是仿效一元二次方程的形式)上式化为：
当采用正则振型解耦时，
可化为：
(注意：由于已经关于质量正则化，故这里的mi=1)
上述方程的解为：
其中：
利用x=ΦNqN可求得物理坐标下的响应。
注意，我们已知的初始条件是物理坐标下的，而解耦方程是模态坐标下的方程，因此，需要将物理坐标下的初始条件转换到模态坐标下去。
如何将初始条件变换到模态坐标(正则坐标)下？
我们知道x=ΦNqN，因此似乎可以直接进行变换。但实际上，由于进行矩阵的求逆运算计算量非常大，而且数值计算中会产生误差，因此，我们不采用此方法进行变换。
将其两端同时左乘则有：
由于：
因此有：
1、求解广义特征值问题，求得特征值λr和特征向量Φr。
广义特征值问题为KΦ=MΦΛ，其特征方程为：
在MATLAB中利用如下语句求得特征值和特征向量：
[V，D]=eig(K，M)；
%特征向量矩阵和特征值矩阵
其中D为特征值矩阵(对角阵)，V为对应的特征向量矩阵。
2、对特征值和特征向量进行排序
求得固有频率ωi和振型Φi后，进而得到振型矩阵Φ(习惯上令振型Φi中第一个元素为1)。
我们习惯上对特征值按从小到大的顺序进行排序，这就要求对得到的特征值进行排序，其对应的特征向量也要进行相应的排序。关于排序的算法，有很多，这里采用冒泡法进行排序。
由于：
可求得圆频率，进而可得到固有频率：
V1=zeros(n，n)；
VN=zeros(n，n)；
for i=1:n
V1(:,i)=V(:,i)/V(1,i)；
%求振型矩阵，各列除以各列的第一个元素
end
3、振型矩阵Φ关于主质量归一化，得到正则振型矩阵ΦN。
由前面论述，可利用：
求得正则振型矩阵Φ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；
%结果为特征值对角阵，与理论相符合
CNp=VN'*C*VN；
%若CNp为对角阵，则原方程组可以解耦
4、将初始条件变换到正则坐标上。
MATLAB中如下实现：
x=zeros(n，1)；
%物理坐标系下位移向量
xd=zeros(n，1)；
%物理坐标系下速度向量
x0=[0 0]'；
%物理坐标系下初位移向量
xd0=[0 0]'；
%物理坐标系下初速度向量
qN=zeros(n，1)；
%模态坐标系下位移向量
qNd=zeros(n，1)；
%模态坐标系下速度向量
%模态坐标系下初位移
qN0=VN'*M*x0；
%物理坐标系下初位移转换到模态坐标系下初位移
%模态坐标系下初速度
qNd0=VN'*M*xd0；
%物理坐标系下初速度转换到模态坐标系下初速度
5、求振系在正则坐标系下的响应。
MATLAB中代码如下：
deltat=0.01；
%正则坐标系下的响应公式中卷积积分t的分段
t0=0；tf=1；
%动力学响应的起止时间
t=t0:deltat:tf；
qNt=zeros(n，length(t))；
%用于存储各时刻模态坐标系下的位移
m=length(t)；
%%%%%%%%%%%%%%%
f=[-50000*1.6*sin(57.5*t)-10*1.6*57.5*cos(57.5*t)；
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)；
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))；
end
6、将正则坐标系下的响应再变回到物理坐标下。
MATLAB中如下实现：
xt=VN*qNt；
%求解物理坐标x
前面讨论的模态叠加法是基于系统的特征值没有重根。
利用振型解耦的前提是已经有了全部的振型向量，如果所有特征值互异，那么肯定存在对应的特征向量，它们对刚度矩阵和质量矩阵正交。
重频特征向量的非唯一性：如果特征方程有两个根λ1=λ2，那么数学上可以证明确实存在Φ1和Φ2两个特征向量(至少就对称的刚度矩阵和质量矩阵，这是成立的)。但是正交性却无法自动保证。(陈奎孚机械振动基础P214)
可以证明：对应于重根λ有无穷多振型向量。
我记得线性代数中指出什么时候即使有重根也会有n个相互独立的特征向量的。方法是施密特正交化。
如果已经知道了重频所对应的两个独立振型向量，则总可以通过线性代数中的格拉姆——施密特正交化过程找到两个正交向量。例如有两个独立振型向量Φ1和Φ2，那么先保留Φ1，而第二个向量假定为：
它肯定也是对应于同一频率的特征向量，现在目标是要调整μ使得上诉变量与已经选定的Φ1也加权正交，也就是：
即：
因为质量矩阵总是正定的，所以上式第二项非系数部分肯定非0，故可以解出：
这样就保证了加权的正交性。
三个重根的证明：
则应有：
即：
由于
故有：
