用模态叠加法求固有频率资料

1、 用模态叠加法求固有频率模态分析法（振型叠加法）原理对于n个自由度系统，其在广义坐标系下的运动微分方程为M & k x F （t）（ 1-1）设在t=0时，有初始条件：x（0） x 0 和 X0）g0通过求解特征值问题，可得系统的固有频率和振型向量niu i (i 1,2,L ,n)i 1 u i (i 1,2,L ,n)以正则振型矩阵作为变换矩阵，令(a)代入方程（1-1），并前乘以正则振型矩阵的转置T，得Z&T F(t)(b)2n12n2O2 nnP(t)F(t)是正则坐标系下的激励。(c)则方程（b）P(t)展开后，得&n舘片(t)2 z n2 2L L2 z nn nP2(t)Pn（t

2、）(1-2)式中Pj（t） T F（t） （i 1,2,L ,n），为对应第i个正则坐标的激励。对于方程（1-2 ）是一组n个独立的方程，每个方程和单自由度系统的强迫振动相同，因此可按单自由度系统中的方法独立地求解每个方程。则由杜哈美积分得方程（1-2 ）的通zi (t)zi0 cos nit也 sin nitni 0( )sinni(t )d i 1,2,L ,nni式中z0和&0是第i个正则坐标的初始位移和初始速度。z0(e)(d)&0前乘以式（两端，得T MZ0z 0T M x0同理，有&0T M &0与成分量形式%T M X0,(i 1,2,L ,n)&0T M &0,(i 1,2,L

3、 ,n)最后，由方程（a),将正则坐标的解T M X0z变换到原广义坐标x，就得到方程（1-1 ）的解。（2 ）在许多工程问题中系统的自由度很多，要想求出系统的所有固有频率和振型向量，计算成本很大，有时甚至是不可能的。由于激励的高频成分很微弱，或者由于系统的高频振动 没有激发出来，总之系统的响应中只有较低的几阶振型分量。因此，使用振型叠加法可以使计算大大地简化。例如，若系统为n自由度，且只需考虑前 p（ p corresp onds to time step nu mberj - corresp onds to a particular DOF*Direct In tegrati on of

4、the un coupled equati ons of moti onare 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 xmi n(i,1),vmi n(i,1)10con ti nue*Tran sform MODAL time histories back to PHYSICAL Coords.x(t) = Phi.q(t) .etc.Print Output *ic=1do 20 t=0.,tto

5、t,deltdo 30 i=1, ndqt(i,1)=q(ic,i) q1t(i,1)=q1(ic,i) q2t(i,1)=q2(ic,i)30con ti nuecall matmul(phi,qt ,n d,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)con ti nuecallmatmul(phi,q1t ,n d,mr, ndim)write (3,43)t,(mr(nn,1),nn=1,nd) do 41 i=1, ndq1(ic,i)=mr(

6、i,1)con ti nuecallmatmul(phi,q2t ,n d,mr, ndim)write (4,43)t,(mr(nn,1),nn=1,nd) do 42 i=1, ndq2(ic,i)=mr(i,1)con ti nueic=ic+120con ti nue*call res_spec(ic-1,q,q1,q2,delt ,ndim,nd)*Gen erate values of MAXIMUM RESPONSE each DOF*stopend*subrouti ne in put1(m,c,k,phi,t,f, nd,nfor,n dim,delt,ttot,xm in,

7、1 vmi n)*m1 is the system mass matrix (in put)k1 is the system stiffness matrix (in put)m contains the modal mass parameters (ge nerated)c contains the modal damp ing parameters (ge nerated)k contains the modal stiffness parameters (ge nerated)alpha, beta = parameters to gen erate proporti onal damp

8、 ing matrixc = alpham + betakphi contains the modal tran sformatio n matrixt contains the times at which the load vectorf is specifiednd = the nu mber of degrees of freedomnsol = 1 for EIGENSOLUTION ONLYn sol = 2 for EIGENSOLUTION + MODE SUPERPOSITION ANALYSISdelt = step sizettot = total time for wh

9、ich resp onse will be evaluated*This program uses the IMSL Math Subrouti ne DGVCSP to evaluatethe eige nvalues and eige nvectors of the dyn amic system*implicit real *8 (a-h,o-z) real *8 m,k,m1,k1 exter nal dgvcspdimension m(ndim,1),c(ndim,1),k(ndim,1),phi(ndim,ndim),t(20,10)dimension nfor(10),f(20,

10、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),xmi n(10,1),vmi n(10,1)read (1,*)nsol,nddo 90005 i=1,ndread(1,*)(m1(i,j),j=1, nd)continuedo 90006 i=1, ndread(1,*)(k1(i,j),j=1, nd)cont

11、inue*do 825 i=1, nddo 825 j=1, ndtemp1(i,j)=m1(i,j)temp2(i,j)=k1(i,j)825 con ti nuec Call the IMSL routi ne DGVCSP to evaluate the eige nvalues and c eige nvectors of the system.call dgvcsp( nd,temp2, ndim,temp1, ndim,eva,evec, ndim)do 903 i=1, ndomega(i)= sqrt (eva(nd-i+1)con ti nuecwrite(*,*)(omeg

12、a(i),i=1, nd)c pausec stopdo 904 i=1, nddo 904 j=1, ndphi(i,j)=evec(i, nd-j+1)contin ue* Normalize Eige nvectors wrt 1st row values*do 826 j=1, ndtemp=phi(1,j)do 826 i=1, ndphi(i,j)=phi(i,j)/temp826con ti nue*Print EIGEN OUTPUTTermi nate if NSOL=1*，/,write (5,905)format (/, E I G E N S O L U T I O N

13、 Mode # Omega(rad/sec)，/, )do 906 i=1, ndwrite (5,907)i,omega(i)format (i6,e22.8)con ti nuewrite (5,912)912format (/, Mode Shapes ,/)do 908 i=1, ndwrite (5, (10E14.6)(phi(i,j),j=1,nd)con ti nueif (nsol.eq.1)stop*read (1,*)delt,ttotread (1,*)alpha,betacall tran s(phi,phit ,nd,n dim)call matmul1(phit,

14、m1,temp1, nd,n dim) call matmul1(temp1,phi,temp2, nd,n dim) do 909 i=1, ndm(i,1)=temp2(i,i)con ti nuecall matmul1(phit,k1,temp1, nd,n dim) call matmul1(temp1,phi,temp2, nd,n dim) do 910 i=1, ndk(i,1)=temp2(i,i)con ti nue911921922914913915916917918923924925do 911 i=1, ndzeta=0.5*(alpha/omega(i)+beta*

15、omega(i) c(i,1)=2.0*m(i,1)*omega(i)*zetacon ti nuedo 921 i=1, nddo 921 j=1, ndc1(i,j)=alpha*m1(i,j)+beta*k1(i,j)con ti nuewrite (5,922)ndformat (/, 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)format (10e15.6)con ti nuewrit

16、e (5,915)format (/, Stiffness Matrix ,/)do 916 i=1, ndwrite (5,914)(k1(i,j),j=1,nd)con ti nuewrite (5,917)format (/, Damping Matrix,/)do 918 i=1, ndwrite (5,914)(c1(i,j),j=1,nd)con ti nuewrite (5,923)(m(i,1),i=1,nd)format (/, Modal Masses ,/,10e15.6)write (5,924)(k(i,1),i=1,nd)format (/, Modal Stiff

17、nesses ,/,10e15.6)write (5,925)(c(i,1),i=1,nd)format (/, Modal Damping ,/,10e15.6)read (1,*)(xin(i,1),i=1,nd)read (1,*)(vin(i,1),i=1,nd)do 851 i=1,nd do 851 j=1,nd temp1(i,j)=0.0 temp2(i,j)=0.0 temp3(i,j)=0.0 851 continuedo 852 i=1,nd temp1(i,i)=1.0/m(i,1) 852 continuecall matmul1(temp1,phit,temp2,n

18、d,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,nd nfor(i)=2 t(1,i)=0.0 t(2,i)=ttot+100.f(1,i)=0.0f(2,i)=0.090011 continue90014 read (1,*,err=90012)i,nfor(i)do 90013 j=1,nfor(i)read (1,*)t(j,i

19、),f(j,i)90013 continuegoto 9001490012 continuereturnend*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),f1(10,1),pt(10,10), nfor(n dim),ft(10,1) do 90030 i=1,nddo 90030 j=1, ndpt(i,j)=p(j,i)90030 continueic=1do 90020

20、 t=0.,ttot,dt*In terpolate to obta in value of forcing function Ftat time = t*do 90022 n f=1, nddo 90024 n=1, nfor(n f)-1if (t.ge.tt(n,nf).and.t.le.tt(n+1,nf)thenft( nf,1)=f( n,n f)+(f( n+1, nf)-f(n,n f)*1(t-tt (n,n f)/(tt (n+1, nf)-tt( n,nf)goto 90022en dif90024continue90022 continue*Gen erate MODA

21、L force vector at time = t*call matmul(pt,ft, nd,f1, ndim)do 90026 ii=1, ndfr(ic,ii)=f1(ii,1)90026 continueic=ic+190020 continuereturnend*subrouti ne n ewb(ii,m,c,k,f,x,x1,x2,del, ndim,ttot,xi n,vin)*Direct In tegrati on of the un coupled equati ons of moti onusi ng the NEWMARK-b method*implicit rea

22、l *8 (a-h,o-z)real *8 m,k,ksdimension f(6000,ndim),x(6000,ndim),x1(6000,ndim),x2(6000,ndim)*Assig n In itial Con diti ons ( t=0.0)*x(1,ii)=x inx1(1,ii)=vi n x2(1,ii)=(f(1,ii)-c*x1(1,ii)-k*x(1,ii)/m ks=k+2.0*c/del+4.0*m/del*2*Obta in time history at discrete time in tervals of del*ic=1do 101 t=0.,tto

23、t,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(*,i

24、6,e17.7,/,5x,2f15.9,/,5x,2e20.10)ic=ic+1101con ti nuereturnend*subrouti nematmul(a,b ,n d,c ,n dim)*Used to multiply two matricesC = 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,nd c(i,j)=c(i,j)+a(i,ii)*b(ii,j)9030 continu

25、ereturnend*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,nd xl(i)=x(1,i) x1l(i)=x1(1,i) x2l(i)=x2(1,i) t(i)=0.t1(i)=0. t2(i)=0.9070 continuedo 9060 i=2,ndo 9060 j=1,ndif ( abs(x(i,j).gt.