用模态叠加法求固有频率.

用模态叠加法求固有频率一、 模态分析法（振型叠加法）原理对于n 个自由度系统，其在广义坐标系下的运动微分方程为[]{}[]{}{}()M x k x F t += （1-1）设在t=0时，有初始条件：{}{}(0)0x x = 和 {}{}(0)0x x =通过求解特征值问题，可得系统的固有频率和振型向量{},(1,2,,)u i n ni iω= {}}(1,2,,)u i n i i ϕ=以正则振型矩阵[]ϕ作为变换矩阵，令{}[]{}x z ϕ= （a ）代入方程（1-1），并前乘以正则振型矩阵的转置T ϕ⎡⎤⎣⎦，得[][][]{}[][][]{}[]{}()T T T M z k z F t ϕϕϕϕϕ+= （b ）∵ [][][][]T M I ϕϕ=[][][][]21222n T k n nn ωωϕϕω⎡⎤⎢⎥⎢⎥=Λ=⎢⎥⎢⎥⎢⎥⎣⎦令 {}[]{}()()T P t F t ϕ= ---- 是正则坐标系下的激励。

则方程（b ）为{}[]{}{}()z z P t +Λ= （c ）展开后，得2()11112()22222()z z P t n z z P t n z z P t n nn n nωωω⎧+=⎪⎪+=⎨⎪⎪+=⎩ （1-2）式中 {}{}()()(1,2,,)T P t F t i n i i ϕ== ，为对应第i 个正则坐标的激励。

对于方程（1-2）是一组n 个独立的方程，每个方程和单自由度系统的强迫振动相同，因此可按单自由度系统中的方法独立地求解每个方程。

则由杜哈美积分得方程（1-2）的通解()0()cos sin 01()sin ()1,2,,0z i z t z t t i i ni ni nit P t d i n ni i ni ωωωτωττω=+⎰+-=式中0z i 和 0z i 是第i 个正则坐标的初始位移和初始速度。

∵ {}[]{}x z ϕ=∴ {}[]{}00x z ϕ= （d ）和 {}[]{}00x z ϕ= （e ）用 [][]T M ϕ 前乘以式（d ）两端，得[][]{}[][][]{}00T T M x M z ϕϕϕ=∴ {}[][]{}00T z M x ϕ=同理，有 {}[][]{}00T z M x ϕ=写成分量形式{}[]{},(1,2,,)00T z M x i n i i ϕ=={}[]{},(1,2,,)00T z M x i n i i ϕ==最后，由方程（a ），将正则坐标的解{}z 变换到原广义坐标{}x ，就得到方程（1-1）的解。

（2）在许多工程问题中系统的自由度很多，要想求出系统的所有固有频率和振型向量，计算成本很大，有时甚至是不可能的。

由于激励的高频成分很微弱，或者由于系统的高频振动没有激发出来，总之系统的响应中只有较低的几阶振型分量。

因此，使用振型叠加法可以使计算大大地简化。

例如，若系统为n 自由度，且只需考虑前p （p<<n ）阶固有特性，即只需计算出系统前p 阶固有频率和振型向量：{},(1,2,,)i p ni i ωϕ=则振型矩阵 []{}{}{}12p ϕϕϕϕ⎡⎤=⎢⎥⎣⎦是一个n ×p 的矩阵。

作如下的线性变换 {}[]{}11x z n p n p ϕ=⨯⨯⨯代入方程（1-1），并前乘以正则振型矩阵的转置[]T ϕ，得[][][]{}[][][]{}[]{}()T T T M z k z F t ϕϕϕϕϕ+= （b ）[][][][]T M I p n n n n p p p ϕϕ=⨯⨯⨯⨯[][][][]21222n T n k p n n n n p p p pp p pωωϕϕω⎡⎤⎢⎥⎢⎥==Λ⎢⎥⨯⨯⨯⨯⎢⎥⎢⎥⎣⎦⨯ {}[]{}()()11T P t F t p n p n ϕ=⨯⨯⨯∴ {}[]{}{}()111z z P t p p p p p +Λ=⨯⨯⨯⨯展开后，得2()11112()22222()z z P t n z z P t n z z P t ppp p p ωωω⎧+=⎪⎪+=⎨⎪⎪+=⎩ 则原广义坐标下的响应为{}[]{}{}{}{}121112x z z z z p n p p n p ϕϕϕϕ==+++⨯⨯⨯这种处理方法称为振型截断法。

二、振型叠加法计算步骤1. 建立广义坐标系下的运动微分方程[]{}[]{}{}()M x k x F t +=给出广义坐标系下初始条件：{}0x 和 {}0x2. 计算固有频率和振型向量{},(1,2,,)u i n ni i ω=计算主质量： {}[]{}(1,2,,)T M u M u i n a i i i== 得正则振型向量： {}}(1,2,,)u i n i iϕ= 3. 初始条件变换{}[]{}(1,2,,)00T z M x i n i i ϕ== {}[]{}(1,2,,)00T z M x i n i i ϕ== 4. 激励变换{}{}()()(1,2,,)T P t F t i n i i ϕ==5. 计算正则坐标系下的响应正则坐标系下运动微分方程2()(1,2,,)z z P t i n i ni i i ω+==利用杜哈美积分写出正则坐标系下的响应0()cos sin 01()sin ()(1,2,,)0z i z t z t t i i ni ni nit P t d i n ni i ni ωωωτωττω=+⎰+-=6. 写出原广义坐标系下的响应{}[]{}{}{}{}1212x z z z z n n ϕϕϕϕ==+++三、源程序******************************* Program to generate the response of an N-DOF structure * to any arbitrary applied loading using the* MODE SUPERPOSITION METHOD******************************* Program Author : Sanjoy Chakraborty* Auburn University******************************program mode_suse msimslimplicit real *8 (a-h,o-z)real *8 mr,krcharacter *12 flin,foutdimension mr(10,1),cr(10,1),kr(10,1),phi(10,10),tt(20,10),1 ft(20,10),fr(6000,10),q(6000,10),q1(6000,10),q2(6000,10),2 qt(10,1),q1t(10,1),q2t(10,1),nfor(10),xmin(10,1),vmin(10,1)data ndim/10/******************************** Assign I/O Devices and read in system parameters*******************************write(*,50)50format(///,' Program for Eigensolution and Mode Superposition'1 ,//,' Program Author : Sanjoy Chakraborty, Auburn University'2 ,//////,' Enter Input File Name : ')read(*,'(a)')flinopen(1,file=flin,status='old')write(*,'(//a)')' Enter Output File Name : 'read(*,'(a)')foutopen(5,file=fout,status='unknown')call input1(mr,cr,kr,phi,tt,ft,nd,nfor,ndim,delt,ttot,xmin,vmin)open(2,file='x.out',status='unknown')open(3,file='x1.out',status='unknown')open(4,file='x2.out',status='unknown')******************************** Generate Modal Force Vectors at discrete time intervals* fr(i,j) contains the modal force vector Fr at time step 'i'* j=1,nd corresponds to the DOF at which fr is applied * {fr} = [Phi]'.{ft}* where, {ft} = physical force vector at time step 'i'* [Phi] = modal transformation matrix for given system* [Phi]' = transpose of [Phi]*******************************call force1(tt,ft,delt,phi,nd,fr,nfor,ndim,ttot)******************************** Generate MODAL time history for all DOF's* q(i,j), q1(i,j), q2(i,j) contain the modal* displacements, velocities and accelerations resp.* i --> corresponds to time step number* j --> corresponds to a particular DOF******************************** Direct Integration of the uncoupled equations of motion* are carried out using the Newmark-B method*******************************do 10 i=1,ndcall newb(i,mr(i,1),cr(i,1),kr(i,1),fr,q,q1,q2,delt,ndim,ttot,1 xmin(i,1),vmin(i,1))10continue******************************** Transform MODAL time histories back to PHYSICAL Coords.* {x(t)} = [Phi].{q(t)} ..etc.* Print Output*******************************ic=1do 20 t=0.,ttot,deltdo 30 i=1,ndqt(i,1)=q(ic,i)q1t(i,1)=q1(ic,i)q2t(i,1)=q2(ic,i)30continuecall matmul(phi,qt,nd,mr,ndim)write(2,43)t,(mr(nn,1),nn=1,nd)43format(f8.4,10e14.6)do 40 i=1,ndq(ic,i)=mr(i,1)40continuecall matmul(phi,q1t,nd,mr,ndim)write(3,43)t,(mr(nn,1),nn=1,nd)do 41 i=1,ndq1(ic,i)=mr(i,1)41continuecall matmul(phi,q2t,nd,mr,ndim)write(4,43)t,(mr(nn,1),nn=1,nd)do 42 i=1,ndq2(ic,i)=mr(i,1)42continueic=ic+120continue*************************call res_spec(ic-1,q,q1,q2,delt,ndim,nd)************************** Generate values of MAXIMUM RESPONSE* @ each DOF*************************stopend*************************subroutine input1(m,c,k,phi,t,f,nd,nfor,ndim,delt,ttot,xmin,1 vmin)************************** [m1] is the system mass matrix (input)* [k1] is the system stiffness matrix (input)* [m] contains the modal mass parameters (generated)* [c] contains the modal damping parameters (generated)* [k] contains the modal stiffness parameters (generated)* alpha, beta = parameters to generate proportional damping matrix * [c] = alpha[m] + beta[k]* [phi] contains the modal transformation matrix* {t} contains the times at which the load vector* {f} is specified* nd = the number of degrees of freedom* nsol = 1 for EIGENSOLUTION ONLY* nsol = 2 for EIGENSOLUTION + MODE SUPERPOSITION ANALYSIS* delt = step size* ttot = total time for which response will be evaluated************************** This program uses the IMSL Math Subroutine DGVCSP to evaluate* the eigenvalues and eigenvectors of the dynamic system*************************implicit real *8 (a-h,o-z)real *8 m,k,m1,k1external dgvcspdimension m(ndim,1),c(ndim,1),k(ndim,1),phi(ndim,ndim),t(20,10)dimension nfor(10),f(20,10),m1(10,10),c1(10,10)dimension k1(10,10),omega(10)dimension phit(10,10),eva(10),evec(10,10)dimension temp1(10,10),temp2(10,10),temp3(10,10),xin(10,1),1 vin(10,1),xmin(10,1),vmin(10,1)read(1,*)nsol,nddo 90005 i=1,ndread(1,*)(m1(i,j),j=1,nd)90005continuedo 90006 i=1,ndread(1,*)(k1(i,j),j=1,nd)90006continue********************do 825 i=1,nddo 825 j=1,ndtemp1(i,j)=m1(i,j)temp2(i,j)=k1(i,j)825continuec Call the IMSL routine DGVCSP to evaluate the eigenvalues and c eigenvectors of the system.call dgvcsp(nd,temp2,ndim,temp1,ndim,eva,evec,ndim)do 903 i=1,ndomega(i)=sqrt(eva(nd-i+1))903continuec write(*,*)(omega(i),i=1,nd)c pausec stopdo 904 i=1,nddo 904 j=1,ndphi(i,j)=evec(i,nd-j+1)904continue******************************* Normalize Eigenvectors wrt 1st row values******************************do 826 j=1,ndtemp=phi(1,j)do 826 i=1,ndphi(i,j)=phi(i,j)/temp826continue******************************* Print EIGEN OUTPUT* Terminate if NSOL=1******************************write(5,905)905format(/,' E I G E N S O L U T I O N',//, 1' Mode # Omega(rad/sec)',/,2' ----------------------------')do 906 i=1,ndwrite(5,907)i,omega(i)907format(i6,e22.8)906continuewrite(5,912)912format(//,' Mode Shapes',/)do 908 i=1,ndwrite(5,'(10E14.6)')(phi(i,j),j=1,nd)908continueif(nsol.eq.1)stop****************************************read(1,*)delt,ttotread(1,*)alpha,betacall trans(phi,phit,nd,ndim)call matmul1(phit,m1,temp1,nd,ndim)call matmul1(temp1,phi,temp2,nd,ndim)do 909 i=1,ndm(i,1)=temp2(i,i)909continuecall matmul1(phit,k1,temp1,nd,ndim)call matmul1(temp1,phi,temp2,nd,ndim)do 910 i=1,ndk(i,1)=temp2(i,i)910continuedo 911 i=1,ndzeta=0.5*(alpha/omega(i)+beta*omega(i))c(i,1)=2.0*m(i,1)*omega(i)*zeta911continuedo 921 i=1,nddo 921 j=1,ndc1(i,j)=alpha*m1(i,j)+beta*k1(i,j)921continuewrite(5,922)nd922format(//,' S Y S T E M P R O P E R T I E S',1 //,' Degrees of Freedom =',i3,//,' Mass Matrix',/)do 913 i=1,ndwrite(5,914)(m1(i,j),j=1,nd)914format(10e15.6)913continuewrite(5,915)915format(//,' Stiffness Matrix',/)do 916 i=1,ndwrite(5,914)(k1(i,j),j=1,nd)916continuewrite(5,917)917format(//,' Damping Matrix',/)do 918 i=1,ndwrite(5,914)(c1(i,j),j=1,nd)918continuewrite(5,923)(m(i,1),i=1,nd)923format(//,' Modal Masses',//,10e15.6)write(5,924)(k(i,1),i=1,nd)924format(//,' Modal Stiffnesses',//,10e15.6)write(5,925)(c(i,1),i=1,nd)925format(//,' Modal Damping',//,10e15.6)read(1,*)(xin(i,1),i=1,nd)read(1,*)(vin(i,1),i=1,nd)do 851 i=1,nddo 851 j=1,ndtemp1(i,j)=0.0temp2(i,j)=0.0temp3(i,j)=0.0851continuedo 852 i=1,ndtemp1(i,i)=1.0/m(i,1)852continuecall matmul1(temp1,phit,temp2,nd,ndim)call matmul1(temp2,m1,temp3,nd,ndim)call matmul(temp3,xin,nd,xmin,ndim)call matmul(temp3,vin,nd,vmin,ndim)********************************** Input of forcing Function*********************************do 90011 i=1,ndnfor(i)=2t(1,i)=0.0t(2,i)=ttot+100.f(1,i)=0.0f(2,i)=0.090011continue90014read(1,*,err=90012)i,nfor(i)do 90013 j=1,nfor(i)read(1,*)t(j,i),f(j,i)90013continuegoto 9001490012continuereturnend**************************subroutine force1(tt,f,dt,p,nd,fr,nfor,ndim,ttot)**************************implicit real *8 (a-h,o-z)dimension tt(20,10),f(20,10),fr(6000,ndim),p(ndim,ndim),1 f1(10,1),pt(10,10),nfor(ndim),ft(10,1)do 90030 i=1,nddo 90030 j=1,ndpt(i,j)=p(j,i)90030continueic=1do 90020 t=0.,ttot,dt******************************** Interpolate to obtain value of forcing function Ft* at time = t*******************************do 90022 nf=1,nddo 90024 n=1,nfor(nf)-1if(t.ge.tt(n,nf).and.t.le.tt(n+1,nf))thenft(nf,1)=f(n,nf)+(f(n+1,nf)-f(n,nf))*1 (t-tt(n,nf))/(tt(n+1,nf)-tt(n,nf))goto 90022endif90024continue90022continue******************************** Generate MODAL force vector at time = t*******************************call matmul(pt,ft,nd,f1,ndim)do 90026 ii=1,ndfr(ic,ii)=f1(ii,1)90026continueic=ic+190020continuereturnend*******************************subroutine newb(ii,m,c,k,f,x,x1,x2,del,ndim,ttot,xin,vin) ******************************** Direct Integration of the uncoupled equations of motion * using the NEWMARK-b method******************************implicit real *8 (a-h,o-z)real *8 m,k,ksdimension f(6000,ndim),x(6000,ndim),x1(6000,ndim),x2(6000,ndim)******************************* Assign Initial Conditions (@ t=0.0)******************************x(1,ii)=xinx1(1,ii)=vinx2(1,ii)=(f(1,ii)-c*x1(1,ii)-k*x(1,ii))/mks=k+2.0*c/del+4.0*m/del**2******************************* Obtain time history at discrete time intervals of del******************************ic=1do 101 t=0.,ttot,deldps=f(ic+1,ii)+m*(4.0*x(ic,ii)/del**2+4.0*x1(ic,ii)1 /del+x2(ic,ii))+c*(2.0*x(ic,ii)/del+x1(ic,ii))x(ic+1,ii)=dps/ksx2(ic+1,ii)=(4.0/del**2)*(x(ic+1,ii)-x(ic,ii)-1 del*x1(ic,ii))-x2(ic,ii)x1(ic+1,ii)=(2.0/del)*(x(ic+1,ii)-x(ic,ii))-x1(ic,ii)c write(5,105)ic,t,f(ic+1,ii),f(ic,ii),dps,dxc105 format('**',i6,e17.7,/,5x,2f15.9,/,5x,2e20.10)ic=ic+1101continuereturnend*****************************subroutine matmul(a,b,nd,c,ndim)******************************* Used to multiply two matrices* [C] = [A].{B}*****************************implicit real *8 (a-h,o-z)dimension a(ndim,ndim),b(ndim,1),c(ndim,1)do 9030 i=1,nddo 9030 j=1,1c(i,j)=0.0do 9030 ii=1,ndc(i,j)=c(i,j)+a(i,ii)*b(ii,j)9030continuereturnend*****************************subroutine res_spec(n,x,x1,x2,del,ndim,nd)****************************** Obtain the maximum value of the response* (displacement, velocity or acceleration)* at each DOF for the time span considered****************************implicit real *8 (a-h,o-z)dimension x(6000,ndim),x1(6000,ndim),x2(6000,ndim)dimension xl(10),x1l(10),x2l(10),t(10),t1(10),t2(10)do 9070 i=1,ndxl(i)=x(1,i)x1l(i)=x1(1,i)x2l(i)=x2(1,i)t(i)=0.t1(i)=0.t2(i)=0.9070continuedo 9060 i=2,ndo 9060 j=1,ndif(abs(x(i,j)).gt.abs(xl(j)))thenxl(j)=x(i,j)t(j)=float(i-1)*delendifif(abs(x1(i,j)).gt.abs(x1l(j)))thenx1l(j)=x1(i,j)t1(j)=float(i-1)*delendifif(abs(x2(i,j)).gt.abs(x2l(j)))thenx2l(j)=x2(i,j)t2(j)=float(i-1)*delendif9060continuewrite(5,9068)9068format(///,' R E S P O N S E S P E C T R A',///, 1' DOF Maxm. @ t Maxm. @ t ', 2' Maxm. @ t',/,3' X Vel. ',4' Accln.',/,5' ----------------------------------------------------', 6'----------------')do 9065 i=1,ndwrite(5,9066)i,xl(i),t(i),x1l(i),t1(i),x2l(i),t2(i)9066format(i5,e15.7,f7.3,e15.7,f7.3,e15.7,f7.3)9065continuereturnend*************************************subroutine matmul1(a,b,c,nd,ndim)*************************************implicit real*8 (a-h,o-z)dimension a(10,10),b(10,10),c(10,10)do 80010 i=1,nddo 80010 j=1,ndc(i,j)=0.0do 80010 k=1,ndc(i,j)=c(i,j)+a(i,k)*b(k,j)80010continuereturnend***************************************subroutine trans(p,pt,nd,ndim)***************************************implicit real*8 (a-h,o-z)dimension p(10,10),pt(10,10)do 80001 i=1,nddo 80001 j=1,ndpt(i,j)=p(j,i)80001continuereturnend***********************************************。