newmark法计算响应的收敛问题
用Fortran编写的程序，不收敛？请各位大侠看看，帮助分析一下。程序是王勖成书中的程序
质量、刚度矩阵的存储形式已经由原来的半带宽存储改为全矩阵存储的形式，特殊需要
质量、刚度矩阵应该没有问题，因为计算出来的固有频率是正确的
类似的问题在论坛中的讨论不少，但是也没有找到比较好的结论
C =============================== SUB: 0602=============================
SUBROUTINENEWMARK(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
DIMENSIONGMM(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 COMPUTEK'=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,你好
你能跟我说一下你的问题出现在哪里么？
我现在做毕业设计，出现的也是这方面的问题呢，急死了！
谢谢！ newmark算法的思路还是比较清晰简单的，我的问题是程序编制中除了问题。
AW=GKK+C0*GMM+C1*DAMP的AW是不变的矩阵，而在调用
CALL AGAUS(AW,GP,NUMPT2,L)时，每一次的AW都从AGAUS中传回来，发生了变化，所以程序计算发散。
稍加修改就不发散了。
页:
[1]