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

1模态叠加法算法理论及其编程实现作者：谢胜龙以下是本人硕士期间学习总结，供大家学习参考，如有错误，请多包涵。参考资料：陈奎孚《机械振动基础》、《少齿差星轮型减速器弹性动力学分析》一、系统动力学方程的解耦——坐标变换，将物理坐标x通过坐标变换xΦq转换到模态空间对系统原振动微分方程进行解耦，转换为模态坐标q下的（用模态坐标表示的）系统振动微分方程。对于系统的振动微分方程：MxCxKxf（1）这里质量矩阵和刚度矩阵均为实对称矩阵，阻尼矩阵也为实对称矩阵。由于此常微分方程组中各个微分方程相互耦合，无法单独求解。因此，要想办法通过坐标变换，将此微分方程组解耦——即各个微分方程相互独立，可以单独求解。令xΦq带入上式有：MΦqCΦqKΦqf（2）两边同乘以TΦ得：TTTTΦMΦqΦCΦqΦKΦqΦf（3）当方程可以解耦时，应有TpΦMΦM（4）TpΦCΦC（5）TpΦKΦK（6）式中pM、pC和pK均为如下形式的对角阵：120000=00nmmmpM，120000=00ncccpC，120000=00nkkkpK（7）即：ppppMqCqKqf（8）其中TpfΦf（9）公式（8）为n阶微分方程组，各个坐标iq（1,2,,in）之间相互独立，可以单独求解，即实现了原方程组的解耦。问题是，如何求得此Φ？二、坐标变换矩阵Φ的求解——求解广义特征值问题2将（4）式变为11TpΦMΦM，带入（6）式得11ppMΦMKΦK（10）1ppKΦMΦMK（11）其中pM和pK均为对角阵，1ppMK可合并为一个对角阵：121000000nppMK（12）其中对角线上元素为iiikm（1,2,,in）（13）于是式（11）变为KΦMΦ（14）式（14）为矩阵理论中的广义特征值问题。如果将矩阵Φ按列分开12[,,,]nΦ（15）其中12rrrrn（1,2,,rn）（16）由分块矩阵乘法1212120000[,,,][,,,]00nnnKKKMMM（17）121122[,,,][,,,]nnnKKKMMM（18）则式（14）就变为如下熟悉的特征值形式：111222nnnKMKMKM（19）之所以叫“广义”是因为M的存在。如果矩阵M为单位阵，那么就是常规的特征值问题。于是，原方程（1）的解耦问题最终归结为求解特征值问题（14），相应的特征值方程为0KM（20）一旦求出这个多项式的根，带回式（14），便可以确定变换矩阵Φ。三、归一化——关于主质量归一化3习惯上我们让方程的第一项系数为1，即TpMΦMΦE问题是，如何找到这样的Φ，使得pME，我们令此时的Φ为NΦ，则有TpNNMΦMΦE问题是NΦ如何求？首先引入振型的正交性方法一：由前面的论述11112122122212121122[,,,]000000TTTTnTTTTnnTTTTnnnnnTTTnnTpMMMMMMMΦMΦMMMMMMM可见，若ij，则有0TijM，因为jjjKM，两边左乘Ti可得TTijjijKM，因此有0TijK。由这两个式子可见：任意2个振型之间，既有对M的正交性，又有对K的正交性，它们统称为振型的正交性。方法二：设振系的第i个与第j个振型向量分别为i和j，按照振型方程（19）有iiiKM（3-1）jjjKM（3-2）用Tj左乘式（3-1）有：TTjiijiKM（3-3）用Ti左乘式（3-2）有：TTijjijKM（3-4）因为M和K都是对称矩阵，故将式（3-4）转置后得：TTjijjiKM（3-5）式（3-3）减式（3-4）得()0TijjiM若ij，则有正交关系0TijM将其带入式（3-3）得正交关系0TijK4由这两个式子可见：任意2个振型之间，既有对M的正交性，又有对K的正交性，它们统称为振型的正交性。因此，由振型的正交性有：1111212212221212111222[,,,]000000000000TTTTnTTTTnnTTTTnnnnnTTTnnnmmmTpMMMMMMMΦMΦMMMMMMM即TrrrmM故对于NΦ，有11112122122212121122[,,,]000000TTTTNNNNNNNnTTTTNNNNNNNnNNNnTTTTNnNnNNnNNnNnTNNTNNTNnNnTpNNMMMMMMMΦMΦMMMMMMME因此，必须有1TNrNrM（1,2,,rn）令Nrrr带入上式有1TrrrrM，即21rrm因此，11rTrrrmM因此，121122[,,,][,,,]NNNnnnNΦ此时，由式（14）KΦMΦ得TTpNNNNKΦKΦΦMΦΛΛ。因此，对原系统振动微分方程做如下变换：（1）令NNxΦq带入上式有：NNNNNNMΦqCΦqKΦqf（2）两边同乘以TNΦ得：TTTTNNNNNNNNNNΦMΦqΦCΦqΦKΦqΦf即：5NpNNNpqCqΛqf上式中Nq称为正则坐标，Npf为正则激励力。四、解耦方程的求解对于动力学方程ppppMqCqKqf，当阻尼矩阵pC可对角化时，其分量形式为：iiiiiipimqcqkqf（1,2,,in……）两边同除以im有：piiiiiiiiifckqqqmmm其中2iniikm，令2iiniicm，则有22iiiiniiiccmmk（之所以这样做是仿效一元二次方程的形式）上式化为22piiiniiniiifqqqm当采用正则振型解耦时，NpNNNpqCqΛqf，22piiiniiniiifqqqm可化为：22NiiniNiniNiNiqqqf（注意：由于已经关于质量正则化，故这里的im=1）上述方程的解为：()0(0)(0)1()[(0)cossin()sin()]iniinitttNiiniNiNiNididiNidididiqqqteqttfetd其中，21diini利用NNxΦq可求得物理坐标下的响应。注意，我们已知的初始条件是物理坐标下的，而解耦方程是模态坐标下的方程，因此，需要将物理坐标下的初始条件转换到模态坐标下去。如何将初始条件变换到模态坐标（正则坐标）下？我们知道NNxΦq，因此，利用-1NNqΦx似乎可以直接进行变换。但实际上，由于进行矩阵的求逆运算计算量非常大，而且数值计算中会产生误差，因此，我们不采用此方法进行变换。将其两端同时左乘-1TNNMΦM，则有-1T-1TNNNNNNMΦMxMΦMΦq由于-1NME，TNNNΦMΦM，因此有TNNqΦMx。五、模态叠加法步骤1、求解广义特征值问题，求得特征值r和特征向量r。6广义特征值问题为KΦMΦ，其特征方程为：0KM在MATLAB中利用如下语句求得特征值和特征向量：[V,D]=eig(K,M);%特征向量矩阵和特征值矩阵其中D为特征值矩阵（对角阵），V为对应的特征向量矩阵。2、对特征值和特征向量进行排序，求固有频率i和振型i，进而得到振型矩阵Φ（习惯上令振型i中第一个元素为1）。我们习惯上对特征值按从小到大的顺序进行排序，这就要求对得到的特征值进行排序，其对应的特征向量也要进行相应的排序。关于排序的算法，有很多，这里采用冒泡法进行排序。利用ii求得圆频率。2iif求得固有频率。V1=zeros(n,n);VN=zeros(n,n);fori=1:nV1(:,i)=V(:,i)/V(1,i);%求振型矩阵，各列除以各列的第一个元素end3、振型矩阵Φ关于主质量归一化，得到正则振型矩阵NΦ。由前面论述，可利用Nrrr，11rTrrrmM求得正则振型矩阵NΦ。Mp=V1'*M*V1;%主质量矩阵Kp=V1'*K*V1;%主刚度矩阵Cp=V1'*C*V1;fori=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、将初始条件变换到正则坐标上。TNNqΦMx(0)(0)TNNqΦMxMATLAB中如下实现：x=zeros(n,1);%物理坐标系下位移向量xd=zeros(n,1);%物理坐标系下速度向量x0=[00]';%物理坐标系下初位移向量xd0=[00]';%物理坐标系下初速度向量qN=zeros(n,1);%模态坐标系下位移向量qNd=zeros(n,1);%模态坐标系下速度向量%模态坐标系下初位移qN0=VN'*M*x0;%物理坐标系下初位移转换到模态坐标系下初位移7%模态坐标系下初速度qNd0=VN'*M*xd0;%物理坐标系下初速度转换到模态坐标系下初速度5、求振系在正则坐标系下的响应。()0(0)(0)1()[(0)cossin()sin()]iniinitttNiiniNiNiNididiNidididiqqqteqttfetd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);%动力学响应中的卷积积分部分%%%%%%%%%%%%%%%%%%%%%%%%%%fori=1:nh=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));end6、将正则坐标系下的响应再变回到物理坐标下。NNxΦqM