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

模态叠加法算法理论及其编程实现作者：谢胜龙以下是本人硕士期间学习总结，供大家学习参考，如有错误，请多包涵。

参考资料： 陈奎孚《机械振动基础》、《少齿差星轮型减速器弹性动力学分析》一、系统动力学方程的解耦——坐标变换，将物理坐标x 通过坐标变换=x Φq 转换到模态空间对系统原振动微分方程进行解耦，转换为模态坐标q 下的（用模态坐标表示的）系统振动微分方程。

对于系统的振动微分方程：++=MxCx Kx f （1） 这里质量矩阵和刚度矩阵均为实对称矩阵，阻尼矩阵也为实对称矩阵。

由于此常微分方程组中各个微分方程相互耦合，无法单独求解。

因此，要想办法通过坐标变换，将此微分方程组解耦——即各个微分方程相互独立，可以单独求解。

令=x Φq 带入上式有：++=M ΦqC Φq K Φq f （2） 两边同乘以TΦ得：++=T T T T ΦM ΦqΦC Φq ΦK Φq Φf （3） 当方程可以解耦时，应有=T p ΦM ΦM （4） =T p ΦC ΦC （5） =T p ΦK ΦK （6）式中p M 、p C 和p K 均为如下形式的对角阵：120000=0n m m m ⎛⎫ ⎪⎪⎪⎪⎝⎭p M ，120000=00n c c c ⎛⎫⎪⎪⎪⎪⎝⎭p C，120000=0n k k k ⎛⎫⎪⎪⎪⎪⎝⎭p K（7）即：++=p p p p M qC q K q f （8） 其中=T p f Φf （9）公式（8）为n 阶微分方程组，各个坐标i q （1,2,,i n = ）之间相互独立，可以单独求解，即实现了原方程组的解耦。

问题是，如何求得此Φ？二、坐标变换矩阵Φ的求解——求解广义特征值问题将（4）式变为11--=T p ΦM ΦM，带入（6）式得11--=p p M ΦM K ΦK （10）1-=p p K ΦM ΦM K （11）其中p M 和p K 均为对角阵，1-p p M K 可合并为一个对角阵：121000000n λλλ-⎛⎫⎪ ⎪Λ==⎪⎪⎝⎭p p M K（12）其中对角线上元素为ii ik m λ=（1,2,,i n = ）（13） 于是式（11）变为=ΛK ΦM Φ（14）式（14）为矩阵理论中的广义特征值问题。

如果将矩阵Φ按列分开12[,,,]n φφφ=Φ （15）其中12r r r rn φφφφ⎧⎫⎪⎪⎪⎪=⎨⎬⎪⎪⎪⎪⎩⎭（1,2,,r n = ）（16） 由分块矩阵乘法1212120000[,,,][,,,]00n n n λλφφφφφφλ⎛⎫⎪ ⎪= ⎪⎪⎝⎭K K K M M M（17） 121122[,,,][,,,]n n n φφφλφλφλφ=K K K M M M （18）则式（14）就变为如下熟悉的特征值形式：111222n n nφλφφλφφλφ===K M K M K M（19）之所以叫“广义”是因为M 的存在。

如果矩阵M 为单位阵，那么就是常规的特征值问题。

于是，原方程（1）的解耦问题最终归结为求解特征值问题（14），相应的特征值方程为0λ-=K M （20）一旦求出这个多项式的根λ，带回式（14），便可以确定变换矩阵Φ。

三、归一化——关于主质量归一化习惯上我们让方程的第一项系数为1，即==Tp M ΦM ΦE问题是，如何找到这样的Φ，使得=p M E ，我们令此时的Φ为N Φ，则有==Tp N N M ΦM ΦE 问题是N Φ如何求？ 首先引入振型的正交性 方法一： 由前面的论述11112122122212121122[,,,]000000T T T T n T T TTn n T T TTn nn n n T TTn n φφφφφφφφφφφφφφφφφφφφφφφφφφφφφφ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦T p M M M M M M M ΦM ΦM M M M M M M可见，若i j λλ≠，则有0Ti j φφ=M ，因为j j j φλφ=K M ，两边左乘T i φ可得T T i j j i j φφλφφ=K M ，因此有0Ti j φφ=K 。

由这两个式子可见：任意2个振型之间，既有对M 的正交性，又有对K 的正交性，它们统称为振型的正交性。

方法二：设振系的第i 个与第j 个振型向量分别为i φ和j φ，按照振型方程（19）有i i i φλφ=K M （3-1） j j j φλφ=K M （3-2）用Tj φ左乘式（3-1）有：T Tj i i j i φφλφφ=K M （3-3）用Ti φ左乘式（3-2）有：T T i j j i j φφλφφ=K M （3-4）因为M 和K 都是对称矩阵，故将式（3-4）转置后得：T T j i j j i φφλφφ=K M （3-5）式（3-3）减式（3-4）得()0T i j j i λλφφ-=M若i j λλ≠，则有正交关系0Ti j φφ=M 将其带入式（3-3）得正交关系0T i j φφ=K由这两个式子可见：任意2个振型之间，既有对M 的正交性，又有对K 的正交性，它们统称为振型的正交性。

因此，由振型的正交性有：1111212212221212111222[,,,]000000000000T T T T n T T T Tn n T T TTn n nn n T TTn n n m m m φφφφφφφφφφφφφφφφφφφφφφφφφφφφφφ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦T p M M M M M M M ΦM ΦM M M M M M M即Tr r r m φφ=M 故对于N Φ，有11112122122212121122[,,,]000000T T T TN N N N N N Nn T T T TN N N N N N Nn N N Nn T T T TNn Nn N Nn N Nn Nn T N N TNN TNn Nn φφφφφφφφφφφφφφφφφφφφφφφφφφφφφφ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦T p N N M M M M M M M ΦM ΦM M M M M M ME因此，必须有1TNr Nr φφ=M （1,2,,r n = ）令Nr r r φμφ=带入上式有1Tr r r r μφμφ=M ，即21r r m μ=因此，r μ==因此，121122[,,,][,,,]N N Nn n n φφφμφμφμφ==N Φ此时，由式（14）=ΛK ΦM Φ得===TTp N N N N K ΦK ΦΦM ΦΛΛ。

因此，对原系统振动微分方程做如下变换： （1）令=N N x Φq 带入上式有：++=N N N N N N M ΦqC Φq K Φq f （2）两边同乘以TN Φ得：++=T T T TN N N N N N N N N N ΦM Φq ΦC Φq ΦK Φq Φf即：++=N p N N Np qC q Λq f 上式中N q 称为正则坐标，Np f 为正则激励力。

四、解耦方程的求解对于动力学方程++=p p p p M qC q K q f ，当阻尼矩阵p C 可对角化时，其分量形式为： i i i i i i pi m qc q k q f ++= （1,2,,i n =……） 两边同除以i m 有：pi i i i i i i i if c kqq q m m m ++= 其中2i ni i k m ω=，令2i i ni ic m ξω=，则有2i i i ni c m ξω==（之所以这样做是仿效一元二次方程的形式）上式化为22pi i i ni i ni i if q q q m ξωω++=当采用正则振型解耦时，++=N p N N Np q C q Λq f ，22pi i i ni i ni i if q q q m ξωω++=可化为：22Ni i ni Ni ni Ni Ni q q q f ξωω++= （注意：由于已经关于质量正则化，故这里的i m =1）上述方程的解为：()0(0)(0)1()[(0)cos sin ()sin ()]i ni i ni tt t Ni i ni Ni Ni Ni di di Ni di didiqq q t e q t t f e t d ξωξωτξωωωτωττωω---+=++-⎰其中，di ni ω=利用=N N x Φq 可求得物理坐标下的响应。

注意，我们已知的初始条件是物理坐标下的，而解耦方程是模态坐标下的方程，因此，需要将物理坐标下的初始条件转换到模态坐标下去。

如何将初始条件变换到模态坐标（正则坐标）下？我们知道=N N x Φq ，因此，利用=-1N N q Φx 似乎可以直接进行变换。

但实际上，由于进行矩阵的求逆运算计算量非常大，而且数值计算中会产生误差，因此，我们不采用此方法进行变换。

将其两端同时左乘-1TN N M ΦM ，则有=-1T-1TN N N N N N M ΦMx M ΦM Φq由于=-1N M E ，=TN N N ΦM ΦM ，因此有=TN N q ΦMx 。

五、模态叠加法步骤1、求解广义特征值问题，求得特征值r λ和特征向量r φ。

广义特征值问题为=ΛK ΦM Φ，其特征方程为：0λ-=K M在MATLAB 中利用如下语句求得特征值和特征向量： [V ,D]=eig(K,M); %特征向量矩阵和特征值矩阵 其中D 为特征值矩阵（对角阵），V 为对应的特征向量矩阵。

2、对特征值和特征向量进行排序，求固有频率i ω和振型i φ，进而得到振型矩阵Φ（习惯上令振型i φ中第一个元素为1）。

我们习惯上对特征值按从小到大的顺序进行排序，这就要求对得到的特征值进行排序，其对应的特征向量也要进行相应的排序。

关于排序的算法，有很多，这里采用冒泡法进行排序。

利用i ω=2ii f ωπ=求得固有频率。

V1=zeros(n,n); VN=zeros(n,n); for i=1:nV1(:,i)=V(:,i)/V(1,i); %求振型矩阵，各列除以各列的第一个元素 end3、振型矩阵Φ关于主质量归一化，得到正则振型矩阵N Φ。

由前面论述，可利用Nr r r φμφ=，r μ==N Φ。

Mp=V1'*M*V1; %主质量矩阵Kp=V1'*K*V1; %主刚度矩阵 Cp=V1'*C*V1; for i=1:nVN(:,i)=V1(:,i)/sqrt(Mp(i,i)); %求正则振型矩阵 end验证程序：MNp=VN'*M*VN; %结果为1对角阵，与理论相符合 KNp=VN'*K*VN; %结果为特征值对角阵，与理论相符合 CNp=VN'*C*VN; %若CNp 为对角阵，则原方程组可以解耦 4、将初始条件变换到正则坐标上。