newmark法计算响应的收敛问题

laneliu

用Fortran编写的程序，不收敛？请各位大侠看看，帮助分析一下。

程序是王勖成书中的程序

质量、刚度矩阵的存储形式已经由原来的半带宽存储改为全矩阵存储的形式，特殊需要

质量、刚度矩阵应该没有问题，因为计算出来的固有频率是正确的

类似的问题在论坛中的讨论不少，但是也没有找到比较好的结论

C =============================== SUB: 0602  =============================
      SUBROUTINE  NEWMARK(GMM,GKK,DAMP,GP,U0,V0,A0)
C ========================================================================
C     SOLVE DYNAMIC RESPONSE BY NEWMARK METHOD
C ========================================================================
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/COM1/MND,NUMEL,NUMPT,MBAND,NMATI
      COMMON/COM2/NF,NFSTR,MSOLV,MPROB,MTYPE,NVA
      COMMON/COM3/MND2,NUMPT2
      COMMON/DYN/XMUV,OMEGA,CC1,CC2,TT,DT,ALPHA,DELTA
      DIMENSION  GMM(NUMPT2,NUMPT2),GKK(NUMPT2,NUMPT2),
     $           AW(NUMPT2,NUMPT2),A0(NUMPT2),B(NUMPT2),
     $           DAMP(NUMPT2,NUMPT2),GP(NUMPT2),U0(NUMPT2),V0(NUMPT2),
     $           U1(NUMPT2),V1(NUMPT2),A1(NUMPT2)

      CLOSE(31,STATUS='DELETE')
        CLOSE(32,STATUS='DELETE')
      OPEN(31,FILE='00OUT_NMK.DAT',STATUS='UNKNOWN')
        OPEN (32, FILE='00NOD_DISP_NMK.DAT',STATUS='UNKNOWN')
        WRITE(32,201)
      WRITE(31,*) 'DYNAMIC RESPONSE RESULT BY NEWMARK METHOD'
      WRITE(31,*) 'TOTAL TIME=', TT
      WRITE(*,101)
101   FORMAT(/6X,' #  SOLVE DYNAMIC RESPONSE BY NEWMARK METHOD  #')
      WRITE(*,102)
102   FORMAT(/6X, '% OUTPUT DYNAMIC DISPLACEMENT IN FILE<OUT_NMK>')
C **************************************************************************
C     INITIAL COMPUTATIONS
C **************************************************************************
      C0=1.0/(ALPHA*DT*DT)
      C1=DELTA/(ALPHA*DT)
      C2=1.0/(ALPHA*DT)
      C3=0.5/ALPHA-1.0
      C4=DELTA/ALPHA-1.0
      C5=DT/2.0*(DELTA/ALPHA-2.0)
      C6=DT*(1.0-DELTA)
      C7=DELTA*DT
C **********************************************************************
C     COMPUTE  K'=K+C0*M+C1*C
C **********************************************************************
        B=GP
      AW=GKK+C0*GMM+C1*DAMP
C **********************************************************************
C     COMPUTATIONS FOR EACH TIME STEP
C **********************************************************************
      DO 300 Y=0,TT,DT
         GP=B
           IF(OMEGA.GT.0.0)  AP=SIN(OMEGA*Y)
         IF(OMEGA.LT.0.0)  AP=COS(OMEGA*Y)
         IF(OMEGA.NE.0.0)  GP=GP*AP

         DO 50 I=1,NUMPT2
         DO 50 J=1,NUMPT2          ! 计算t+Δt时刻的有效载荷
            GP(I)=GP(I)+GMM(I,J)*(C0*U0(J)+C2*V0(J)+C3*A0(J))+
     $            DAMP(I,J)*(C1*U0(J)+C4*V0(J)+C5*A0(J))
50       CONTINUE      

         CALL AGAUS(AW,GP,NUMPT2,L) !计算t+Δt时刻的位移
           IF(L.EQ.0) WRITE(*,80)
80       FORMAT(/6X, "NEWMARK方程组系数矩阵奇异")

           U1=GP                        ! 当前时刻的位移
           A1=C0*(U1-U0)-C2*V0-C3*A0    ! 当前时刻的速度
           V1=V0+C6*A0+C7*A1            ! 当前时刻的加速度

           U0=U1                        !下一时刻的位移、速度、加速度
           V0=V1
           A0=A1  
C ******** OUTPUT NODE DISPLACEMENT ******************************************                    
           NSTEP=Y/DT+1
         WRITE(31,61) NSTEP,DT,Y+DT
         WRITE(31,71)
         WRITE(31,72) (U1(I),I=1,NUMPT2)
C ********* OUTPUT REFERED NODE DISPLACEMENT *********************************
         NDISP=5        
         WRITE(32,200) Y+DT,U1((NDISP-1)*NF+1),U1((NDISP-1)*NF+2),
     $                      U1((NDISP-1)*NF+3),U1((NDISP-1)*NF+4)      
300   CONTINUE
C ****************************************************************************
61    FORMAT(2X,'NO. OF STEP=',I5,3X,'STEP LENGTH =',E15.6,3X,
     $        'AT TIME=',E15.6)
71    FORMAT(2X,'DISPLACEMENT:'/2(13X,'X-',13X,'Y-'))
72    FORMAT(4X,4E16.8)
201   FORMAT(8X,'TIME=',12X,'X-',15X,'Y-',12X,'TX-',15X,'TY-')
200   FORMAT(1X, 5E16.8)
      CLOSE (31)
      RETURN
      END

[ 本帖最后由 laneliu 于 2008-5-15 17:50 编辑 ]

laneliu

问题已经解决，编程有问题。

xiaowenwen1107

laneliu,你好
你能跟我说一下你的问题出现在哪里么？
我现在做毕业设计，出现的也是这方面的问题呢，急死了！
谢谢！

laneliu

newmark算法的思路还是比较清晰简单的，我的问题是程序编制中除了问题。

AW=GKK+C0*GMM+C1*DAMP的AW是不变的矩阵，而在调用

CALL AGAUS(AW,GP,NUMPT2,L)时，每一次的AW都从AGAUS中传回来，发生了变化，所以程序计算发散。

稍加修改就不发散了。