用模態(tài)疊加法求固有頻率資料

1、 用模態(tài)疊加法求固有頻率模態(tài)分析法（振型疊加法）原理對于n個(gè)自由度系統(tǒng)，其在廣義坐標(biāo)系下的運(yùn)動微分方程為M & k x F （t）（ 1-1）設(shè)在t=0時(shí)，有初始條件：x（0） x 0 和 X0）g0通過求解特征值問題，可得系統(tǒng)的固有頻率和振型向量niu i (i 1,2,L ,n)i 1 u i (i 1,2,L ,n)以正則振型矩陣作為變換矩陣，令(a)代入方程（1-1），并前乘以正則振型矩陣的轉(zhuǎn)置T，得Z&T F(t)(b)2n12n2O2 nnP(t)F(t)是正則坐標(biāo)系下的激勵(lì)。(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），為對應(yīng)第i個(gè)正則坐標(biāo)的激勵(lì)。對于方程（1-2 ）是一組n個(gè)獨(dú)立的方程，每個(gè)方程和單自由度系統(tǒng)的強(qiáng)迫振動相同，因此可按單自由度系統(tǒng)中的方法獨(dú)立地求解每個(gè)方程。則由杜哈美積分得方程（1-2 ）的通zi (t)zi0 cos nit也 sin nitni 0( )sinni(t )d i 1,2,L ,nni式中z0和&0是第i個(gè)正則坐標(biāo)的初始位移和初始速度。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),將正則坐標(biāo)的解T M X0z變換到原廣義坐標(biāo)x，就得到方程（1-1 ）的解。（2 ）在許多工程問題中系統(tǒng)的自由度很多，要想求出系統(tǒng)的所有固有頻率和振型向量，計(jì)算成本很大，有時(shí)甚至是不可能的。由于激勵(lì)的高頻成分很微弱，或者由于系統(tǒng)的高頻振動 沒有激發(fā)出來，總之系統(tǒng)的響應(yīng)中只有較低的幾階振型分量。因此，使用振型疊加法可以使計(jì)算大大地簡化。例如，若系統(tǒng)為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.