a20.for

ADINA84有限元编程学习的好例子

      SUBROUTINE EL3D10
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1            ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /DIMEL/ N101,N102,N103,N104,N105,N106,N107,N108,N109,N110,
     1               N111,N112,N113,N114,N120,N121,N122,N123,N124,N125
      COMMON /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
      COMMON /DPR/ ITWO
      COMMON A(1)
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
C
      REAL A
C
      DIMENSION IA(1)
C
      EQUIVALENCE (NPAR(10),NINT),(A(1),IA(1)),(NPAR(17),NCON),
     1            (NPAR(7),MXNODS),(NPAR(11),NINTZ)
C
C
C
C
C
      IDW=47*ITWO
      NPT=NINT*NINT*NINTZ
C
C
      MATP=IA(N107 + NEL - 1)
C
C
      NM=N111 + (MATP - 1)*NCON*ITWO
C
C
      IF(IND.NE.0) GO TO 100
      NN=N112+(NEL-1)*(IDW*NPT+MXNODS)
      CALL IEPC3(A(NN),A(NN + IDW*NPT),A(NN),IDW,A(N6A + ITWO),A(NM))
      GO TO 200
C
C
  100 NN=N112+(NEL-1)*(IDW*NPT+MXNODS)+(IPT-1)*IDW
      NN1=NN
      NN2=NN+6*ITWO
      NN3=NN+12*ITWO
      NN4=NN+18*ITWO
      NN5=NN+24*ITWO
      NN6=NN+25*ITWO
      NN7=NN + 26*ITWO
      NN8=NN + 32*ITWO
      NN9=NN + 44*ITWO
      NN10=NN + 45*ITWO
      NN11=NN + 46*ITWO
C
C
      KK=N112+(NEL-1)*(IDW*NPT+MXNODS)+IDW*NPT
C
C
      ND9DIM=MXNODS-8
      LL=N108+(NEL-1)*ND9DIM
C
C
      CALL EPC3(A(NM),A(NN1),A(NN2),A(NN3),A(NN4),A(NN5),A(NN6),A(NN7),
     1          A(NN8),A(NN9),A(NN10),A(NN11),A(KK),A(LL),A(N6A + ITWO),
     2          A(N6B + ITWO))
C
  200 CONTINUE
C
      RETURN
      END
      SUBROUTINE IEPC3(WA,IWA,IIWA,IIDW,TEMPV1,PROP)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1            ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /EM3D/ NOD(27),NODM(27),NOD9M(19)
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),E1,E2,E3
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /DPR/ ITWO
      COMMON /THREED/ DE,IELD,IELX,IEL,NPT,IDW,NND9,ISOCOR
      COMMON /SOLPM3/ TREF,XINTP,XCON1,XCON2,XPARM1,XPARM2,TOLIL,TOLPC,
     1                TOL1,TOL2,TOL3,TOL4,TOL5,TOL6,TCHK,CRPCON(8),DTMN,
     2                SUBDD,RNGL,RNGU,DTT,TOL7,
     3                KCRP,NITE,NALG,IINTP,NPTS,ITCHK
C
      DIMENSION WA(47,*),IWA(*),IIWA(IIDW,*),TEMPV1(*),H(27),
     1          XDM1(3,27),XDM2(3,3),XDM3(3,1),PROP(16,*),PROP1(5)
C
      EQUIVALENCE (NPAR(10),NINT),(NPAR(11),NINTZ)
C
      NPT=NINT*NINT*NINTZ
      IINTP=1
      NPTS=INT(PROP(9,7))
      TOLMT=1.D-2
C
      TOLL=TOLMT*ABS(PROP(1,1))
      IF (TOLL.EQ.0.D0) TOLL=TOLMT
      TOLU=TOLMT*ABS(PROP(NPTS,1))
      IF (TOLU.EQ.0.D0) TOLU=TOLMT
C
      RNGL=PROP(1,1) - TOLL
      RNGU=PROP(NPTS,1) + TOLU
C
C           TO ZERO
C
   15 DO 20 J=1,NPT
      DO 20 I=1,45
   20 WA(I,J)=0.D0
C
C
      II=0
      DO 25 K=1,27
      IF(NODM(K).EQ.0) GO TO 25
      II=II+1
      IWA(II)=NODM(K)
   25 CONTINUE
C
C
      IPT=0
      DO 30 LX=1,NINT
      E1=XG(LX,NINT)
      DO 30 LY=1,NINT
      E2=XG(LY,NINT)
      DO 30 LZ=1,NINTZ
      E3=XG(LZ,NINTZ)
      IPT=IPT+1
      CALL FUNCT(E1,E2,E3,H,XDM1,NOD9M,XDM2,XDUM,XDM3,IINTP)
      TEMP1=0.D0
      DO 35 K=1,IEL
      KK=IWA(K)
   35 TEMP1=TEMP1+H(K)*TEMPV1(KK)
      WA(45,IPT)=TEMP1
C
C
      CALL MTITP3(PROP,TEMP1,PROP1)
      YS1=PROP1(3)
   30 WA(25,IPT)=YS1
C
C           TO ONE
C
      KJ=45*ITWO + 1
      KJJ=46*ITWO + 1
      DO 40 I=1,NPT
      IIWA(KJ,I)=1
   40 IIWA(KJJ,I)=1
C
      RETURN
C
      END
      SUBROUTINE EPC3(PROP,SIG,EPS,EPSP,EPSC,YLD,EPSTR,ALFA,ORIG,TMPOLD,
     1                IPEL,NORG,NDS,NOD9M,TEMPV1,TEMPV2)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /THREED/ DE,IELD,IELX,IEL,NPT,IDW,NND9,ISOCOR
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1            ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
      COMMON /TEMPRT/ TEMP1,TEMP2,ITEMPR,ITP96,N6A,N6B
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),Z1,Z2,Z3
      COMMON /SOLPM3/ TREF,XINTP,XCON1,XCON2,XPARM1,XPARM2,TOLIL,TOLPC,
     1               TOL1,TOL2,TOL3,TOL4,TOL5,TOL6,TCHK,CRPCON(8),DTMN,
     2               SUBDD,RNGL,RNGU,DTT,TOL7,
     3               KCRP,NITE,NALG,IINTP,NPTS,ITCHK
      COMMON /CONST/ DT,DTA,CONS(21),DTOD,IOPE
      COMMON /PORT/ INPORT,JNPORT,NPUTSV,IBPORT,IFPORT,JDC,JVC,JAC
C
      DIMENSION PROP(16,*),SIG(*),EPS(*),EPSC(*),ALFA(*),ORIG(6,*),
     1          ORIGD(6,2),NDS(*),NOD9M(*),TEMPV1(*),TEMPV2(*),
     2          H(27),XDM1(3,27),XDM2(3,3),XDM3(3,1),
     3          DELSIG(6),DELEPS(6),DEPS(6),EPSP1(6),EPSP2(6),
     4          STRSS1(6),STRSS2(6),STRSSM(6),EPSC1(6),EPSC2(6),
     5          EPSCM(6),DPSC(6),ALFA1(6),ALFA2(6),ALFAM(6),
     6          EPS1(6),EPS2(6),PROP1(5),DEPST(6),PROP2(5),PROPM(5)
      DIMENSION STRSSD(6),DPSP(6),EPST2(6),EPSP(6),CEP(6,6),EPST1(6),
     1          DSTSS(6)
C
      EQUIVALENCE (NPAR(15),MODEL),(NPAR(3),INDNL),
     1            (STRESS(1),STRSS2(1))
C
      CHARACTER*2 STATE(2)
      DATA STATE /' E' , '*P'/
C
C
C
C
      INDEX=1
      ISUB=1
      ISUBM=INT(PROP(13,7))
C
      IF(IPT.GT.1) GO TO 5
C
C
      DTT=DT
      IF (MODEX.EQ.2 .AND. KSTEP.EQ.1 .AND. IEQUIT.EQ.1 .AND.
     1   KPRI.NE.0) DTT=DTOD
C
      SUBDD=5.D0
      DTMN=DTT/PROP(13,7)
      IINTP=1
C
      DO 2 J=1,8
    2 CRPCON(J)=PROP(J,7)
C
      NPTS=INT(PROP(9,7))
      TREF=PROP(10,7)
      KCRP=INT(PROP(11,7))
      XINTP=PROP(12,7)
      NITE=INT(PROP(14,7))
      NALG=INT(PROP(15,7))
      TOLIL=PROP(16,7)
      TOLPC=PROP(1,8)
C
      XCON1=2.D0/3.D0
      XCON2=1.D0/3.D0
C
      ITCHK=1
      IF(NITE.LT.6) ITCHK=0
C
C
      XPARM1=1.0 - XINTP
      XPARM2=XINTP
C
C     SET TOLERANCES **
C
      TOL1=TOLIL*TOLIL
      TOL4=5.D-6
      TOL5=1.D-20
      TOL2=TOL5*TOL5
      TOL3=2.0*TOL4
      TOL6=0.1D0
      TOL7=2.D0
      TOLMT=1.D-2
      TCHK=DTT*(1.0 - TOL4)
C
      TOLL=TOLMT*ABS(PROP(1,1))
      IF (TOLL.EQ.0.D0) TOLL=TOLMT
      TOLU=TOLMT*ABS(PROP(NPTS,1))
      IF (TOLU.EQ.0.D0) TOLU=TOLMT
C
      RNGL=PROP(1,1) - TOLL
      RNGU=PROP(NPTS,1) + TOLU
C
C
    5 DO 10 I=1,6
      EPS1(I)=EPS(I)
      EPSP1(I)=EPSP(I)
      EPSP2(I)=EPSP(I)
      ALFA1(I)=ALFA(I)
      ALFA2(I)=ALFA(I)
      EPSC1(I)=EPSC(I)
      EPSC2(I)=EPSC(I)
      DPSC(I)=0.D0
      EPST1(I)=0.D0
      EPST2(I)=0.D0
      DEPST(I)=0.D0
      STRSS2(I)=SIG(I)
   10 STRSS1(I)=SIG(I)
C
      YLD1=YLD
      EPSTR1=EPSTR
      EPSTR2=EPSTR
      ECSTR1=0.D0
      CRSRM=0.D0
      TMP1=TMPOLD
      IPELD=IPEL
      NORGD=NORG
      TAU=0.D0
      ESTM=0.D0
C
      DO 20 I=1,6
      DO 20 J=1,2
   20 ORIGD(I,J)=ORIG(I,J)
C
C
      DO 25 J=1,6
   25 DELEPS(J)=STRAIN(J) - EPS(J)
C
C
C
      CALL FUNCT(Z1,Z2,Z3,H,XDM1,NOD9M,XDM2,XDUM,XDM3,IINTP)
      TEMP1=0.D0
      TEMP2=0.D0
C
      DO 30 K=1,IEL
      KK=NDS(K)
      TEMP2=TEMP2 + H(K)*TEMPV2(KK)
   30 TEMP1=TEMP1 + H(K)*TEMPV1(KK)
C
      CTEMP=TEMP2
C
C
      IF (IEQUIT.EQ.1 .AND. KPRI.NE.0) CTEMP=TEMP1
C
      DELTMP=CTEMP - TMPOLD
C
C
      CALL EMAT3(TMPOLD,PROP,PROP1,A1,B1,C1,D1,E1,F1,1)
C
      YM1=PROP1(1)
      ET1=PROP1(4)
      YS1=PROP1(3)
C
      EET1=YM1*ET1/(YM1 - ET1)
C
C
   40 DELT=DTMN
      IF(KCRP.EQ.0.AND.NALG.EQ.2) DELT=DTT
C
C
   60 XFAC=(TAU + DELT)/DTT
      DO 65 J=1,6
      EPS2(J)=EPS(J) + XFAC*DELEPS(J)
   65 DEPS(J)=EPS2(J) - EPS1(J)
C
C
      TMP2=TMPOLD + XFAC*DELTMP
      TMPM=XPARM1*TMP1 + XPARM2*TMP2
C
      CALL EMAT3(TMP2,PROP,PROP2,A2,B2,C2,D2,E2,F2,2)
C
C
      ALPHA2=PROP2(5)
C
      EPST2(1)=ALPHA2*(TMP2 - TREF)
      EPST2(2)=EPST2(1)
      EPST2(3)=EPST2(1)
C
C
      IF(KCRP.EQ.0) GO TO 95
C
   70 DO 75 J=1,6
   75 STRSSM(J)=XPARM1*STRSS1(J) + XPARM2*STRSS2(J)
C
C
      DO 80 J=1,6
   80 DPSC(J)=0.D0
      CRSRM=0.D0
C
      CALL EFST3(ESTM,SXM,SYM,SZM,SXYM,SXZM,SYZM,STRSSM)
      IF(ESTM.LE.TOL5.AND.INDEX.GT.1) GO TO 95
C
      DO 90 J=1,6
   90 EPSCM(J)=XPARM1*EPSC1(J) + XPARM2*EPSC2(J)
C
      CALL CREEP3(DELT,DPSC,TMPM,EPSCM,ORIGD,NORGD,STRSSM,
     1            GAMA,CRSRM,PTIME,ESTM,SXM,SYM,SZM,SXYM,SXZM,SYZM,
     2            FF,RR,GG,FP,INDEX,ECSTRM)
C
      IF(INDEX.EQ.1) ECSTR1=ECSTRM
C
C
   95 CALL SIGMA3(STRSS2,EPS2,EPSP2,EPST2,EPSC1,EPSC2,DPSC,GAMA,PTIME,
     1            CRSRM,FF,RR,GG,FP,ESTM,SXM,SYM,SZM,SXYM,SXZM,SYZM,
     2            DELT,B2,C2,D2)
C
C
C
  100 IF(KCRP.EQ.0) GO TO 215
      IF (XPARM2.EQ.0.D0) GO TO 205
      IF(ITCHK.EQ.1) GO TO 120
C
      INDEX=INDEX + 1
      IF(INDEX.LE.NITE) GO TO 70
      GO TO 205
C
C     STRESS VECTOR **
C
  120 IF(INDEX - 4) 122,135,125
C
  122 INDEX=INDEX + 1
      GO TO 70
C
  125 DNORM2=0.D0
      DO 130 J=1,6
  130 DNORM2=DNORM2 + (STRSS2(J) - STRSSD(J))*(STRSS2(J) - STRSSD(J))
C
C
  135 SNORM=0.D0
      DO 140 J=1,6
  140 SNORM=SNORM + STRSS2(J)*STRSS2(J)
C
C
      IF(INDEX.GT.5) GO TO 155
      SNORM2=SNORM
      IF(INDEX.EQ.4) SNORM1=SNORM2
      IF(INDEX.EQ.5) DNORM1=DNORM2
      INDEX=INDEX + 1
C
      DO 150 J=1,6
  150 STRSSD(J)=STRSS2(J)
      GO TO 70
C
C
C
C
  155 IF(DNORM2.LE.DNORM1) GO TO 185
C
C     TOLERANCE BAND
C
      XTOL=TOL3*SNORM1
      IF(SNORM1.LE.TOL2) XTOL=TOL2
      IF(DNORM1.LE.XTOL.AND.DNORM2.LE.XTOL) GO TO 205
C
C     (NALG .EQ. 2) *
C
      DELT=DELT*(SQRT(DNORM1/DNORM2))/SUBDD
      IF(NALG.EQ.2.AND.ISUB.LT.ISUBM) GO TO 170
C
      write(66,3004)
      write(66,3002) NEL,IPT,ISUB,TAU,DELT
      STOP
C
C
  170 INDEX=1
C
      DO 180 I=1,6
      STRSS2(I)=STRSS1(I)
  180 EPSC2(I)=EPSC1(I)
C
      GO TO 60
C
C
  185 XTOL=TOL1*SNORM1
      IF(SNORM1.LE.TOL2) XTOL=TOL2
      IF(DNORM1.LE.XTOL) GO TO 205
C
C     NO CONVERGENCE
C
  190 INDEX=INDEX + 1
      IF(INDEX.LE.NITE) GO TO 195
C
      write(66,3001)
      write(66,3011) NEL,IPT,ISUB,TAU,DELT
      STOP
C
  195 DNORM1=DNORM2
      SNORM1=SNORM2
      SNORM2=SNORM
C
      DO 200 J=1,6
  200 STRSSD(J)=STRSS2(J)
      GO TO 70
C
C
C
C
C
  205 IF (NALG.EQ.1) GO TO 215
      DECSTR=CRSRM*DELT
      IF(DECSTR.LE.TOL5 .OR. ECSTR1.LE.TOL5) GO TO 215
C
      CHECK=DECSTR/(ECSTR1*TOLPC)
      IF (CHECK.LE.1.1D0) GO TO 215
C
C
      DELT=DELT/CHECK
      IF(NALG.EQ.2.AND.ISUB.LT.ISUBM) GO TO 208
C
      write(66,3006)
      write(66,3002) NEL,IPT,ISUB,TAU,DELT
      STOP
C
C
  208 INDEX=1
C
      DO 210 I=1,6
      EPSC2(I)=EPSC1(I)
  210 STRSS2(I)=STRSS1(I)
C
      GO TO 60
C
C
C     PLASTICITY) **
C
  215 CALL EFST3(EST2,SX2,SY2,SZ2,SXY2,SXZ2,SYZ2,STRSS2)
C
      DO 218 I=1,6
  218 DELSIG(I)=STRSS2(I) - STRSS1(I)
C
      CALL EFST3(EST,DX,DY,DZ,DXY,DXZ,DYZ,DELSIG)
C
C
      IF(MODEL.EQ.10) GO TO 220
C
      SX2=SX2 - ALFA1(1)
      SY2=SY2 - ALFA1(2)
      SZ2=SZ2 - ALFA1(3)
      SXY2=SXY2 - ALFA1(4)
      SXZ2=SXZ2 - ALFA1(5)
      SYZ2=SYZ2 - ALFA1(6)
C
C
  220 YM2=PROP2(1)
      ET2=PROP2(4)
      YS2=PROP2(3)
C
      EET2=YM2*ET2/(YM2 - ET2)
C
      DYLD=YS2 - YS1
      IF(MODEL.EQ.10) DYLD=DYLD + (EET2 - EET1)*EPSTR1
      YLD2=YLD1 + DYLD
C
  225 RA=DX*DX + DY*DY + DZ*DZ + 2.0*(DXY*DXY + DXZ*DXZ + DYZ*DYZ)
      FTA=SX2*SX2 + SY2*SY2 + SZ2*SZ2 + 2.0*(SXY2*SXY2 + SXZ2*SXZ2 +
     1    SYZ2*SYZ2)
C
C
      IF (RA.EQ.0.D0.AND.TMP1.EQ.TMP2) GO TO 228
C
      FTB=XCON1*YLD2*YLD2
      IF(FTA.GT.FTB) GO TO 250
C
C
      IPELD=1
  228 TAU=TAU + DELT
C
C
      IF(NALG.EQ.2) GO TO 230
C
C     NALG .EQ. 1 *
C
      IF(ISUB.EQ.ISUBM) GO TO 245
      GO TO 235
C
C    NALG .EQ. 2 *
C
  230 IF(TAU.GE.TCHK .OR. KCRP.EQ.0) GO TO 245
      IF(DECSTR.LE.TOL5) GO TO 232
C
      DELT=DELT*TOLPC*(1.0 + (ECSTR1/DECSTR))
      IF(TAU + DELT.GE.TCHK) DELT=DTT - TAU
      GO TO 233
C