a10.for

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

      SUBROUTINE EL2D10
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 /MATMOD/ STRESS(4),STRAIN(4),C(4,4),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)
C
C
C
C
C
      IDW=33*ITWO
      NPT=NINT*NINT
C
C
      MATP=IA(N107 + NEL - 1)
C
C
      NM=N109 + (MATP - 1)*NCON*ITWO
C
C
      IF(IND.NE.0) GO TO 100
      NN=N110+(NEL-1)*(IDW*NPT+MXNODS)
      CALL IEPC2(A(NN),A(NN + IDW*NPT),A(NN),IDW,A(N6A + ITWO),A(NM))
      GO TO 200
C
C
  100 NN=N110 + (NEL - 1)*(IDW*NPT + MXNODS) + (IPT - 1)*IDW
      NN1=NN
      NN2=NN + 4*ITWO
      NN3=NN + 8*ITWO
      NN4=NN + 12*ITWO
      NN5=NN + 16*ITWO
      NN6=NN + 17*ITWO
      NN7=NN + 18*ITWO
      NN8=NN + 22*ITWO
      NN9=NN + 30*ITWO
      NN10=NN + 31*ITWO
      NN11=NN + 32*ITWO
C
C
      KK=N110+(NEL-1)*(IDW*NPT+MXNODS)+IDW*NPT
C
C
      ND5DIM=MXNODS-4
      LL=N111+(NEL-1)*ND5DIM
C
C
      CALL EPC2(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 IEPC2(WA,IWA,IIWA,IDW,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 /EM2D/ S(378),XM(27),B(4,18),RE(27),EDIS(27),EDISI(27),
     1              XX(27),NOD(9),NODM(9),NOD5M(5)
      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 /TODIM/ BET,THIC,DE,IEL,NND5,ISOCOR
      COMMON /SOLPM2/ 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,IST,ISR
C
      DIMENSION WA(33,*),IWA(*),IIWA(IDW,*),TEMPV1(*),H(9),
     1          XDM1(2,9),XDM2(2,2),XDM3(2,1),PROP(16,*),PROP1(5)
C
      EQUIVALENCE (NPAR(10),NINT)
C
      NPT=NINT*NINT
      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
C           TO ZERO
C
   15 DO 20 J=1,NPT
      DO 20 I=1,31
   20 WA(I,J)=0.D0
C
C
      II=0
      DO 25 K=1,9
      IF(NODM(K).EQ.0) GO TO 25
      II=II+1
      IWA(II)=NODM(K)
   25 CONTINUE
C
C
      DO 30 LX=1,NINT
      E1=XG(LX,NINT)
      DO 30 LY=1,NINT
      E2=XG(LY,NINT)
      IPT=(LX-1)*NINT+LY
      CALL FUNCT2(E1,E2,H,XDM1,NOD5M,XDM2,XDUM,XDM3,IDUM,IINTP)
      TEMP1=0.D0
      DO 35 K=1,IEL
      KK=IWA(K)
   35 TEMP1=TEMP1+H(K)*TEMPV1(KK)
      WA(31,IPT)=TEMP1
C
C
      CALL MTITP2(PROP,TEMP1,PROP1)
      YS1=PROP1(3)
   30 WA(17,IPT)=YS1
C
C           TO ONE
C
      KJ=31*ITWO + 1
      KJJ=32*ITWO + 1
      DO 40 I=1,NPT
      IIWA(KJ,I)=1
   40 IIWA(KJJ,I)=1
C
      RETURN
C
      END
      SUBROUTINE EPC2(PROP,SIG,EPS,EPSP,EPSC,YLD,EPSTR,ALFA,ORIG,TMPOLD,
     1                IPEL,NORG,NDS,NOD5M,TEMPV1,TEMPV2)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /TODIM/ BET,THIC,DE,IEL,NND5,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 /MATMOD/ STRESS(4),STRAIN(4),C(4,4),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 /SOLPM2/ 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,IST,ISR
      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(4,*),
     1          ORIGD(4,2),NDS(*),NOD5M(*),TEMPV1(*),TEMPV2(*),
     2          H(9),XDM1(2,9),XDM2(2,2),XDM3(2,1),
     3          DELSIG(4),DELEPS(4),DEPS(4),EPSP1(4),EPSP2(4),
     4          STRSS1(4),STRSS2(4),STRSSM(4),EPSC1(4),EPSC2(4),
     5          EPSCM(4),DPSC(4),ALFA1(4),ALFA2(4),ALFAM(4),
     6          EPS1(4),EPS2(4),PROP1(5),DEPST(4),PROP2(5),PROPM(5)
      DIMENSION STRSSD(4),DPSP(4),EPST2(4),EPSP(*),CEP(4,4),EPST1(4),
     1          DSTSS(4)
C
      EQUIVALENCE (NPAR(5),ITYP2D),(NPAR(15),MODEL),(NPAR(3),INDNL)
      EQUIVALENCE (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
      IST=4
      IF(ITYP2D.GE.2) IST=3
      ISR=3
      IF(ITYP2D.EQ.0) ISR=4
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 DELEPS(4)=0.D0
      DEPS(4)=0.D0
      EPS2(4)=0.D0
      DELSIG(4)=0.D0
      STRSSM(4)=0.D0
      STRSSD(4)=0.D0
C
      DO 10 I=1,4
      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,4
      DO 20 J=1,2
   20 ORIGD(I,J)=ORIG(I,J)
C
C
      DO 25 J=1,ISR
   25 DELEPS(J)=STRAIN(J) - EPS(J)
C
C
C
      CALL FUNCT2(Z1,Z2,H,XDM1,NOD5M,XDM2,XDUM,XDM3,NEL,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 EMAT2(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,ISR
      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 EMAT2(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(4)=EPST2(1)
C
C
      IF(KCRP.EQ.0) GO TO 95
C
   70 DO 75 J=1,IST
   75 STRSSM(J)=XPARM1*STRSS1(J) + XPARM2*STRSS2(J)
C
C
      DO 80 J=1,4
   80 DPSC(J)=0.D0
      CRSRM=0.D0
C
      CALL EFST(ESTM,SXM,SYM,SXYM,SZM,STRSSM)
      IF(ESTM.LE.TOL5.AND.INDEX.GT.1) GO TO 95
C
      DO 90 J=1,4
   90 EPSCM(J)=XPARM1*EPSC1(J) + XPARM2*EPSC2(J)
C
      CALL CREEP2(DELT,DPSC,TMPM,EPSCM,ORIGD,NORGD,STRSSM,
     1            GAMA,CRSRM,PTIME,ESTM,SXM,SYM,SXYM,SZM,FF,RR,GG,FP,
     2            INDEX,ECSTRM)
C
      IF(INDEX.EQ.1) ECSTR1=ECSTRM
C
C
   95 CALL SIGMA2(STRSS2,EPS2,EPSP2,EPST2,EPSC1,EPSC2,DPSC,GAMA,PTIME,
     1            CRSRM,FF,RR,GG,FP,ESTM,SXM,SYM,SXYM,SZM,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,IST
  130 DNORM2=DNORM2 + (STRSS2(J) - STRSSD(J))*(STRSS2(J) - STRSSD(J))
C
C
  135 SNORM=0.D0
      DO 140 J=1,IST
  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,IST
  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,4
      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,IST
  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,4
      EPSC2(I)=EPSC1(I)
  210 STRSS2(I)=STRSS1(I)
C
      GO TO 60
C
C
C     PLASTICITY) **
C
  215 CALL EFST(EST2,SX2,SY2,SXY2,SZ2,STRSS2)
C
      DO 218 I=1,IST
  218 DELSIG(I)=STRSS2(I) - STRSS1(I)
C
      CALL EFST(EST,DX,DY,DXY,DZ,DELSIG)
C
C
      IF(MODEL.EQ.10) GO TO 220
C
      SX2=SX2 - ALFA1(1)
      SY2=SY2 - ALFA1(2)
      SXY2=SXY2 - ALFA1(3)
      SZ2=SZ2 - ALFA1(4)
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
      FTA=SX2*SX2 + SY2*SY2 + SZ2*SZ2 + 2.0*SXY2*SXY2
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 248
      GO TO 235
C
C    NALG .EQ. 2 *
C
  230 IF(TAU.GE.TCHK.OR.KCRP.EQ.0) GO TO 248
      IF(DECSTR.LE.TOL5) GO TO 232
C
      DELT=DELT*TOLPC*(1.0 + (ECSTR1/DECSTR))
      IF(TAU + DELT.GE.TCHK) DELT=DTT - TAU