a20.for

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

      EPST2(J)=EPST
  465 DEPST(J)=DPST
C
C
      INDEX=1
      CALL EFST3(EST1,SX1,SY1,SZ1,SXY1,SXZ1,SYZ1,STRSS1)
C
      IF(KCRP.EQ.0) GO TO 480
C
      DO 470 J=1,6
  470 DPSC(J)=0.D0
C
      IF(EST1.LE.TOL5) GO TO 480
      TMPM=XPARM1*TEMP1 + XPARM2*TEMP2
C
      CALL CREEP3(DT,DPSC,TMPM,EPSC1,ORIGD,NORGD,STRSS1,GAMA,CRSR1,
     1            PTIME,EST1,SX1,SY1,SZ1,SXY1,SXZ1,SYZ1,FF,RR,GG,FP,
     2            INDEX,ECSTR1)
C
      DO 475 I=1,6
  475 EPSC2(I)=EPSC1(I) + DPSC(I)
C
C
  480 IF(IPELD.EQ.2) GO TO 490
C
C
      DO 482 J=1,6
  482 STRESS(J)=0.D0
C
      DO 485 I=1,6
      DO 485 J=1,6
  485 STRESS(I)=STRESS(I) + C(I,J)*(STRAIN(J) - EPSP2(J) - EPSC2(J) -
     1          EPST2(J))
C
      RETURN
C
C
C
C
  490 DO 494 J=1,5
  494 PROPM(J)=XPARM1*PROP1(J) + XPARM2*PROP2(J)
C
C
  500 YMM=PROPM(1)
      ETM=PROPM(4)
      EETM=YMM*ETM/(YMM - ETM)
C
      CALL EPMAT3(STRSS1,ALFA1,EPSTR1,DELEPS,DPSC,DPST,CEP,XLAMDA,
     1            PROP1,PROP2,PROPM,YLD1,2,A2,B2,C1,C2,DPSP,SX1,SY1,SZ1,
     2            SXY1,SXZ1,SYZ1,INDEX,EETM)
C
      DO 505 J=1,6
      STRESS(J)=0.D0
  505 EPSP2(J)=EPSP1(J) + DPSP(J)
C
      IF (IEQREF.EQ.1 .OR. XLAMDA.LT.0.D0) GO TO 520
C
C
      DO 510 I=1,6
      DO 510 J=1,6
  510 STRESS(I)=STRESS(I) - CEP(I,J)*(DPSC(J) + DEPST(J)) + C(I,J)*
     1          (STRAIN(J) - EPSP1(J) - EPSC1(J) - EPST1(J) - DPSP(J))
C
      DO 515 I=1,6
      DO 515 J=1,6
  515 C(I,J)=CEP(I,J)
C
      RETURN
C
C
  520 DO 535 I=1,6
      DO 535 J=1,6
  535 STRESS(I)=STRESS(I) + C(I,J)*(STRAIN(J) - EPSP2(J) - EPSC2(J)
     1          - EPST2(J))
C
      RETURN
C
C
C
  600 FT=YLD2
      IF(IPELD.EQ.1) FT=SQRT(1.5*FTA)
C
C
  605 IF(INDNL.EQ.2) CALL CAUCH3
C
  610 IF (IPRI.NE.0 .OR. IPS.EQ.0) GO TO 730
      IF (IPS.LT.0) GO TO 700
C
C
      IF(IPT.GT.1) GO TO 620
C
C     PRINT HEADING *
C
      write(66,2000)
C
C
      write(66,2005) NEL
C
C
  620 write(66,2200) IPT,STATE(IPELD),(STRSS2(J),J=1,6)
      write(66,2700) TEMP2,EPSTR2,ISUB,FT,YLD2,YLDC
C
      GO TO 730
C
C
  700 IF(IPT.GT.1) GO TO 720
C
C     PRINT HEADING *
C
      write(66,2000)
C
C
      write(66,2005) NEL
C
C
  720 write(66,2200) IPT,STATE(IPELD),(STRSS2(J),J=1,6)
      write(66,2300) (STRAIN(J),J=1,6)
      write(66,2400) (EPSP2(J),J=1,6)
      write(66,2500) (EPSC2(J),J=1,6)
      write(66,2600) (EPST2(J),J=1,6)
      write(66,2700) TEMP2,EPSTR2,ISUB,FT,YLD2,YLDC
C
C
  730 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 740
      IF (ISVE.EQ.0) GO TO 740
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-3',NEL,IPT,IPELD,ISUB,
     2                       (STRESS(I),I=1,6),(STRAIN(I),I=1,6),
     3                       (EPSP2(I),I=1,6),(EPSC2(I),I=1,6),
     4                       (EPST2(I),I=1,6),TEMP2,EPSTR2,FT,YLD2,YLDC
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9000) 'OUTPUT-3',NEL,IPT,IPELD,ISUB,
     2                       (STRESS(I),I=1,6),(STRAIN(I),I=1,6),
     3                       (EPSP2(I),I=1,6),(EPSC2(I),I=1,6),
     4                       (EPST2(I),I=1,6),TEMP2,EPSTR2,FT,YLD2,YLDC
C
 9000 FORMAT ( A,/,4I10,/,(4E20.13) )
C
C
  740 CONTINUE
      RETURN
C
 2000 FORMAT (1X,7HELEMENT,4X,6HSTRESS,2X,13HSTRESS/STRAIN,14X,2HXX,
     1        13X,2HYY,13X,2HZZ,13X,2HXY,13X,2HXZ,13X,2HYZ,/,1X,
     2        7HNUM/IPT,5X,5HSTATE,2X,10HCOMPONENTS)
 2005 FORMAT (/,1X,I3)
 2200 FORMAT (6X,I2,2X,A2,6HLASTIC,2X,6HSTRESS,9X,6(E14.6,1X))
 2300 FORMAT (20X,12HSTRAIN-TOTAL,3X,6(E14.6,1X))
 2400 FORMAT (25X,7HPLASTIC,3X,6(E14.6,1X))
 2500 FORMAT (27X,5HCREEP,3X,6(E14.6,1X))
 2600 FORMAT (25X,7HTHERMAL,3X,6(E14.6,1X))
 2700 FORMAT (20X,14HTEMPERATURE = ,E14.6,1X,
     1        29HACCUM. EFF. PLASTIC STRAIN = ,E14.6,1X,
     2        25HNUMBER OF SUBDIVISIONS = ,I5,/,20X,
     3        12HEFF STRESS = ,E14.6,1X,
     4        13HYLD STRESS = ,E14.6,1X,
     5        41HYLD STRESS(TMP.,ACC. EFF. PLAS. STRN.) = ,E14.6,/)
 3001 FORMAT(//,68H    ERROR   STRESS LOOP NO. 1 FAILED TO CONVERGE   (S
     1UBROUTINE EPC3))
 3002 FORMAT(//,5X,10HELEMENT = ,I5,2X,20HINTEGRATION POINT = ,I5,2X,
     1       13HSUBDIVISION = ,I5,/,5X,
     2       38HTIME RELATIVE TO START OF TIME STEP = ,E14.6,2X,
     3       40HAPPROXIMATE REQUIRED SUBDIVISION SIZE = ,E14.6)
 3003 FORMAT(//,68H    ERROR   STRESS LOOP NO. 2 FAILED TO CONVERGE   (S
     1UBROUTINE EPC3))
 3004 FORMAT(//,115H    ERROR   SUBDIVISION SIZE REQUIRED TO ELIMINATE D
     1IVERGENCE IS TOO SMALL IN STRESS LOOP NO. 1   (SUBROUTINE EPC3))
 3006 FORMAT(//,129H    ERROR   SUBDIVISION SIZE REQUIRED TO SATISFY INE
     1LASTIC STRAIN TOLERANCE IS TOO SMALL IN STRESS LOOP NO. 1   (SUBRO
     2UTINE EPC3))
 3007 FORMAT(//,114H    ERROR   MAXIMUM NUMBER OF SUBDIVISIONS WILL NOT
     1REACH END OF TIME STEP IN STRESS LOOP NO. 1 (SUBROUTINE EPC3)    )
 3008 FORMAT(//,115H    ERROR   SUBDIVISION SIZE REQUIRED TO ELIMINATE D
     1IVERGENCE IS TOO SMALL IN STRESS LOOP NO. 2   (SUBROUTINE EPC3))
 3009 FORMAT(//,129H    ERROR   SUBDIVISION SIZE REQUIRED TO SATISFY INE
     1LASTIC STRAIN TOLERANCE IS TOO SMALL IN STRESS LOOP NO. 2   (SUBRO
     2UTINE EPC3))
 3010 FORMAT(//,114H    ERROR   MAXIMUM NUMBER OF SUBDIVISIONS WILL NOT
     1REACH END OF TIME STEP IN STRESS LOOP NO. 2 (SUBROUTINE EPC3)    )
 3011 FORMAT(//,5X,10HELEMENT = ,I5,2X,20HINTEGRATION POINT = ,I5,2X,
     1       14HSUBDIVISION = ,I5,/,5X,
     2       38HTIME RELATIVE TO START OF TIME STEP = ,E14.6,2X,
     3       24HLAST SUBDIVISION SIZE = ,E14.6)
 3012 FORMAT(//,70H    ERROR   INCORRECT VALUE CALCULATED FOR *RATIO*
     1(SUBROUTINE EPC3)        )
 3013 FORMAT(//,5X,10HELEMENT = ,I5,1X,20HINTEGRATION POINT = ,I5,1X,
     1       14HSUBDIVISION = ,I5,1X,
     2       38HTIME RELATIVE TO START OF TIME STEP = ,E14.6,1X,
     3        7HIPEL = ,I2,/,5X,
     4        5HRA = ,E14.6,1X,5HRB = ,E14.6,1X,5HRD = ,E14.6,
     5        5HRE = ,E14.6,1X,5HRF = ,E14.6,1X,5HRG = ,E14.6,/,5X,
     6        8HRATIO = ,E14.6,//)
 3014 FORMAT(//,101H    ERROR   DIFFERENCE BETWEEN THE TWO MEASURES OF Y
     1IELD STRESS EXCEEDS TOLERANCE   (SUBROUTINE EPC3))
 3015 FORMAT(//,5X,10HELEMENT = ,I5,2X,20HINTEGRATION POINT = ,I5,2X,
     1       14HSUBDIVISION = ,I5,2X,
     2       38HTIME RELATIVE TO START OF TIME STEP = ,E14.6,/,5X,
     3       35HYIELD STRESS (FROM STRESS STATE) = ,E14.6,2X,
     4       48HYIELD STRESS (FROM TEMP. AND ACCUM. EFF. PLASTIC   ,
     5       10HSTRAIN) = ,E14.6,//)
C
      END
      SUBROUTINE SIGMA3(STRSS2,EPS2,EPSP2,EPST2,EPSC1,EPSC2,DPSC,GAMA,
     1                  PTIME,STRNR,F,R,G,FP,EST,SX,SY,SZ,SXY,SXZ,SYZ,
     2                  DELT,XB2,XC2,XD2)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      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 /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
C
      DIMENSION STRSS2(*),EPS2(*),EPSP2(*),EPST2(*),EPSC1(*),EPSC2(*),
     1          DPSC(*),TSTRSS(6),RH(6),TC(6,6),TTC(6,6)
C
C
C
C
      DO 10 J=1,6
      TSTRSS(J)=STRSS2(J)
   10 RH(J)=0.D0
C
      DO 20 I=1,6
      DO 20 J=1,6
   20 TC(I,J)=0.D0
C
C
      DO 25 I=1,6
      DO 25 J=1,6
   25 RH(I)=RH(I) + C(I,J)*(EPS2(J) - EPSP2(J) - EPSC1(J) - DPSC(J)
     1      - EPST2(J))
C
      IF (XPARM2.GT.0.D0.AND.KCRP.GT.0.AND.EST.GT.TOL5) GO TO 35
C
      DO 30 I=1,6
      EPSC2(I)=EPSC1(I) + DPSC(I)
   30 STRSS2(I)=RH(I)
C
      RETURN
C
C
   35 A0=CRPCON(1)
      A1=CRPCON(2)
      A2=CRPCON(3)
      A3=CRPCON(4)
      A4=CRPCON(5)
      A5=CRPCON(6)
      A6=CRPCON(7)
      A7=CRPCON(8)
C
      IF (KCRP - 2) 60,40,70
C
C     CREEP LAW NO. 1 **
C
   60 GG=1.5*(A1 - A2)*GAMA / (A2*EST*EST)
      GO TO 50
C
C     CREEP LAW NO. 2 **
C
   40 C2=A4 - 1.0
      C1=1.5*A2*A4*(A3**(-A4))
      C3=1.5*A1
      C4=1.5*A6
C
      D1=EXP(-R*PTIME)*F
      D2=EST**C2
C
      DTPSP=(D1*(C3 - PTIME*C1*D2) - C3*F - PTIME*C4*G)/STRNR
C
      GG=(1.5*D1*(C1*D2 - R*(PTIME*C1*D2 + R*DTPSP - C3)) + 1.5*C4*G -
     1   1.5*GAMA)/(EST*EST)
      GO TO 50
C
C     CREEP LAW NO. 3 **
C
   70 C1=A2*(A2 - 1.)
      C2=A2 - 1.
      C3=A3*A4*(A4 - 1.)
      C4=A4 - 2.
      C5=A5*A6*(A6 - 1.)
      C6=A6 - 2.
C
      HP=C1*PTIME**C2 + C3*PTIME**C4 + C5*PTIME**C6
C
      C7=A0*A1
      C8=A1 - 1.
C
      GG=1.5*((1.5*C7*EST**C8*FP - R*HP*F/(FP*EST))*G - GAMA)/(EST*EST)
C
   50 COEF=2.0*XC2*XPARM2*DELT
C
      TC(1,1)=COEF*(GG*SX*SX + XCON1*GAMA)
      TC(1,2)=COEF*(GG*SX*SY - XCON2*GAMA)
      TC(1,3)=COEF*(GG*SX*SZ - XCON2*GAMA)
      TC(1,4)=COEF*(2.0*GG*SX*SXY)
      TC(1,5)=COEF*(2.0*GG*SX*SXZ)
      TC(1,6)=COEF*(2.0*GG*SX*SYZ)
C
      TC(2,1)=TC(1,2)
      TC(2,2)=COEF*(GG*SY*SY + XCON1*GAMA)
      TC(2,3)=COEF*(GG*SY*SZ - XCON2*GAMA)
      TC(2,4)=COEF*(2.0*GG*SY*SXY)
      TC(2,5)=COEF*(2.0*GG*SY*SXZ)
      TC(2,6)=COEF*(2.0*GG*SY*SYZ)
C
      TC(3,1)=TC(1,3)
      TC(3,2)=TC(2,3)
      TC(3,3)=COEF*(GG*SZ*SZ + XCON1*GAMA)
      TC(3,4)=COEF*(2.0*GG*SZ*SXY)
      TC(3,5)=COEF*(2.0*GG*SZ*SXZ)
      TC(3,6)=COEF*(2.0*GG*SZ*SYZ)
C
      TC(4,1)=TC(1,4)
      TC(4,2)=TC(2,4)
      TC(4,3)=TC(3,4)
      TC(4,4)=COEF*(2.0*GG*SXY*SXY + GAMA)
      TC(4,5)=COEF*(2.0*GG*SXY*SXZ)
      TC(4,6)=COEF*(2.0*GG*SXY*SYZ)
C
      TC(5,1)=TC(1,5)
      TC(5,2)=TC(2,5)
      TC(5,3)=TC(3,5)
      TC(5,4)=TC(4,5)
      TC(5,5)=COEF*(2.0*GG*SXZ*SXZ + GAMA)
      TC(5,6)=COEF*(2.0*GG*SXZ*SYZ)
C
      TC(6,1)=TC(1,6)
      TC(6,2)=TC(2,6)
      TC(6,3)=TC(3,6)
      TC(6,4)=TC(4,6)
      TC(6,5)=TC(5,6)
      TC(6,6)=COEF*(2.0*GG*SYZ*SYZ + GAMA)
C
      DO 55 I=1,6
      DO 55 J=1,6
   55 TTC(I,J)=TC(I,J)/(2.0*XC2)
C
C
   80 DO 95 I=1,6
      DO 95 J=1,6
   95 RH(I)=RH(I) + TC(I,J)*STRSS2(J)
C
      DO 110 J=1,6
  110 TC(J,J)=TC(J,J) + 1.0
C
C
      CALL EQSOL3(TC,RH,1,6)
C
      DO 130 J=1,6
  130 STRSS2(J)=RH(J)
C
C
      DO 160 I=1,6
      DO 150 J=1,6
  150 DPSC(I)=DPSC(I) + TTC(I,J)*(STRSS2(J) - TSTRSS(J))
  160 EPSC2(I)=EPSC1(I) + DPSC(I)
C
      RETURN
C
      END
      SUBROUTINE EQSOL3(A,R,KEY,N)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      DIMENSION A(6,6),R(6),ICOL(6),TR(6),IROW(6)
C
C
C
      EPS=1.D-35
      NN=N-1
      IF(KEY.EQ.2) GO TO 110
C
      DO 10 I=1,N
      IROW(I)=I
   10 ICOL(I)=I
C
      DO 100 J=1,NN
C
C
      PIVI=0.D0
      DO 20 I=J,N
      DO 20 K=J,N
      IF(ABS(A(I,K)).LE.ABS(PIVI)) GO TO 20
      IMAX=I
      KMAX=K
      PIVI=A(I,K)
   20 CONTINUE
C
      IF(ABS(PIVI).GT.EPS) GO TO 25
      write(66,7000)
      write(66,7010) J
      STOP
C
C
   25 IF(KMAX.EQ.J) GO TO 40
C
      ISAVE=ICOL(KMAX)
      ICOL(KMAX)=ICOL(J)
      ICOL(J)=ISAVE
C
      DO 30 JJ=1,N
      SAVE=A(JJ,KMAX)
      A(JJ,KMAX)=A(JJ,J)
   30 A(JJ,J)=SAVE
C
C
   40 IF(IMAX.EQ.J) GO TO 85
C
      ISAVE=IROW(IMAX)
      IROW(IMAX)=IROW(J)
      IROW(J)=ISAVE
C
      DO 50 K=1,N
      SAVE=A(J,K)
      A(J,K)=A(IMAX,K)
   50 A(IMAX,K)=SAVE
C
      SAVE=R(J)
      R(J)=R(IMAX)
      R(IMAX)=SAVE
C
C
   85 I2=J + 1
      DO 90 K2=I2,N
      XMULT=A(K2,J)/PIVI
      A(K2,J)=XMULT
      DO 90 J2=I2,N
   90 A(K2,J2)=A(K2,J2) - XMULT*A(J,J2)
C
  100 CONTINUE
      GO TO 150
C
C
C
  110 DO 120 J=1,N
      K=IROW(J)
  120 TR(J)=R(K)
C
      DO 130 J=1,N
  130 R(J)=TR(J)
C
  150 DO 200 J=1,NN
      I2=J + 1
      DO 200 K2=I2,N
      XMULT=A(K2,J)
  200 R(K2)=R(K2) - XMULT*R(J)
C
C
      DO 320 I=1,N
      TSUM=0.D0
      II=N + 1 - I
      IF(II.EQ.N) GO TO 315
      JJ=II + 1
C
      DO 310 K=JJ,N
  310 TSUM=TSUM + A(II,K)*R(K)
C
  315 R(II)=(R(II) - TSUM)/A(II,II)
  320 CONTINUE
C
C
      DO 330 J=1,N
      K=ICOL(J)
  330 TR(K)=R(J)
C
      DO 340 J=1,N
  340 R(J)=TR(J)
C
      RETURN
C
 7000 FORMAT(///,64H   ERROR   UNABLE TO OBTAIN NON-ZERO PIVOT   (SUBROU
     1TINE EQSOL3))
 7010 FORMAT(/,10X,15HPIVOT NUMBER = ,I5)
C
      END
      SUBROUTINE EFST3(EST,SX,SY,SZ,SXY,SXZ,SYZ,STRESS)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      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 STRESS(6)
C
C
C
      SM=(STRESS(1) + STRESS(2) + STRESS(3))*XCON2
C
      SX=STRESS(1)-SM
      SY=STRESS(2)-SM
      SZ=STRESS(3)-SM
      SXY=STRESS(4)
      SXZ=STRESS(5)
      SYZ=STRESS(6)
C
      EST=SQRT(1.5*(SX*SX + SY*SY + SZ*SZ + 2.0*(SXY*SXY + SXZ*SXZ +
     1         SYZ*SYZ)))
C
      RETURN
C
      END
      SUBROUTINE EPMAT3(STRSSM,ALFAM,EPSTRM,DEPS,DPSC,DPST,CEP,XLAMDA,
     1                  PROP1,PROP2,PROPM,YLDM,KEY,A2,B2,C1,C2,DPSP,
     2                  SXM,SYM,SZM,SXYM,SXZM,SYZM,INDEX,EETM)
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 /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      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 /TPLAS3/ EETD,PRM,CM,CD,YSD
C
      DIMENSION STRSSM(*),ALFAM(*),DEPS(*),DPSC(*),PROP1(*),
     1          PROP2(*),PROPM(*),CEP(6,6),DPSP(6)
C
      EQUIVALENCE (NPAR(15),MODEL)
C
C
C
C
C