a10.for

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

      CALL EMAT2(TEMP2,PROP,PROP2,A2,B2,C2,D2,E2,F2,2)
C
C
      ALPHA2=PROP2(5)
      EPST=ALPHA2*(TEMP2 - TREF)
      DPST=EPST - EPST1(1)
C
      DO 465 J=1,4
      EPST2(J)=EPST
  465 DEPST(J)=DPST
C
      EPST2(3)=0.D0
      DEPST(3)=0.D0
C
C
      INDEX=1
      CALL EFST(EST1,SX1,SY1,SXY1,SZ1,STRSS1)
C
      IF(KCRP.EQ.0) GO TO 480
C
      DO 470 J=1,4
  470 DPSC(J)=0.D0
C
      IF(EST1.LE.TOL5) GO TO 480
      TMPM=XPARM1*TEMP1 + XPARM2*TEMP2
C
      CALL CREEP2(DT,DPSC,TMPM,EPSC1,ORIGD,NORGD,STRSS1,GAMA,CRSR1,
     1            PTIME,EST1,SX1,SY1,SXY1,SZ1,FF,RR,GG,FP,INDEX,ECSTR1)
C
      DO 475 I=1,4
  475 EPSC2(I)=EPSC1(I) + DPSC(I)
C
C
  480 IF(IPELD.EQ.2) GO TO 490
C
C
      DO 482 J=1,4
  482 STRESS(J)=0.D0
C
      DO 485 I=1,IST
      DO 485 J=1,IST
  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 EPMAT2(STRSS1,ALFA1,EPSTR1,DELEPS,DPSC,DPST,CEP,XLAMDA,
     1            PROP1,PROP2,PROPM,YLD1,2,A2,B2,C1,C2,DPSP,SX1,SY1,
     2            SXY1,SZ1,INDEX,EETM)
C
      DO 505 J=1,4
      STRESS(J)=0.D0
  505 EPSP2(J)=EPSP1(J) + DPSP(J)
C
      IF(IEQREF.EQ.1 .OR. XLAMDA.LT. 0.0) GO TO 520
C
C
      DO 510 I=1,IST
      DO 510 J=1,IST
  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,IST
      DO 515 J=1,IST
  515 C(I,J)=CEP(I,J)
C
      RETURN
C
C
  520 DO 535 I=1,IST
      DO 535 J=1,IST
  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 CAUCHY
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(4),(STRSS2(J),J=1,3)
      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(4),(STRSS2(J),J=1,3)
      IF(ITYP2D.GE.2) STRAIN(4)=EPS2(4)
      write(66,2300) STRAIN(4),(STRAIN(J),J=1,3)
      write(66,2400) EPSP2(4),(EPSP2(J),J=1,3)
      write(66,2500) EPSC2(4),(EPSC2(J),J=1,3)
      write(66,2600) EPST2(4),(EPST2(J),J=1,3)
      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 (ITYP2D.GE.2) STRAIN(4)=EPS2(4)
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-2',NEL,IPT,IPELD,ISUB,
     2                       (STRESS(I),I=1,4),(STRAIN(I),I=1,4),
     3                       (EPSP2(I),I=1,4),(EPSC2(I),I=1,4),
     4                       (EPST2(I),I=1,4),TEMP2,EPSTR2,FT,YLD2,YLDC
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9000) 'OUTPUT-2',NEL,IPT,IPELD,ISUB,
     2                       (STRESS(I),I=1,4),(STRAIN(I),I=1,4),
     3                       (EPSP2(I),I=1,4),(EPSC2(I),I=1,4),
     4                       (EPST2(I),I=1,4),TEMP2,EPSTR2,FT,YLD2,YLDC
C
 9000 FORMAT ( A,/,4I10,/,(4E20.13) )
C
C
  740 CONTINUE
C
      RETURN
C
 2000 FORMAT (1X,7HELEMENT,4X,6HSTRESS,2X,13HSTRESS/STRAIN,14X,2HXX,
     1        13X,2HYY,13X,2HZZ,13X,2HYZ,/,1X,7HNUM/IPT,5X,5HSTATE,
     2        2X,10HCOMPONENTS)
 2005 FORMAT (/,1X,I3)
 2200 FORMAT (6X,I2,2X,A2,6HLASTIC,2X,6HSTRESS,9X,4(E14.6,1X),
     1        2(1X,E14.6))
 2300 FORMAT (20X,12HSTRAIN-TOTAL,3X,4(E14.6,1X))
 2400 FORMAT (25X,7HPLASTIC,3X,4(E14.6,1X))
 2500 FORMAT (27X,5HCREEP,3X,4(E14.6,1X))
 2600 FORMAT (25X,7HTHERMAL,3X,4(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 EPC2))
 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 EPC2))
 3004 FORMAT(//,115H    ERROR   SUBDIVISION SIZE REQUIRED TO ELIMINATE D
     1IVERGENCE IS TOO SMALL IN STRESS LOOP NO. 1   (SUBROUTINE EPC2))
 3006 FORMAT(//,129H    ERROR   SUBDIVISION SIZE REQUIRED TO SATISFY INE
     1LASTIC STRAIN TOLERANCE IS TOO SMALL IN STRESS LOOP NO. 1   (SUBRO
     2UTINE EPC2))
 3007 FORMAT(//,114H    ERROR   MAXIMUM NUMBER OF SUBDIVISIONS WILL NOT
     1REACH END OF TIME STEP IN STRESS LOOP NO. 1 (SUBROUTINE EPC2)    )
 3008 FORMAT(//,115H    ERROR   SUBDIVISION SIZE REQUIRED TO ELIMINATE D
     1IVERGENCE IS TOO SMALL IN STRESS LOOP NO. 2   (SUBROUTINE EPC2))
 3009 FORMAT(//,129H    ERROR   SUBDIVISION SIZE REQUIRED TO SATISFY INE
     1LASTIC STRAIN TOLERANCE IS TOO SMALL IN STRESS LOOP NO. 2   (SUBRO
     2UTINE EPC2))
 3010 FORMAT(//,114H    ERROR   MAXIMUM NUMBER OF SUBDIVISIONS WILL NOT
     1REACH END OF TIME STEP IN STRESS LOOP NO. 2 (SUBROUTINE EPC2)   )
 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 EPC2)     )
 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 EPC2))
 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 SIGMA2(STRSS2,EPS2,EPSP2,EPST2,EPSC1,EPSC2,DPSC,GAMA,
     1                  PTIME,STRNR,F,R,G,FP,EST,SX,SY,SXY,SZ,DELT,XB2,
     2                  XC2,XD2)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /MATMOD/ STRESS(4),STRAIN(4),C(4,4),IPT,NEL,IPS,ISVE
      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 STRSS2(*),EPS2(*),EPSP2(*),EPST2(*),EPSC1(*),EPSC2(*),
     1          DPSC(*),TSTRSS(4),RH(4),TC(4,4),TTC(4,4)
C
C
C
      XFAC1=-XD2
C
      DO 10 J=1,4
      TSTRSS(J)=STRSS2(J)
   10 RH(J)=0.D0
C
      DO 20 I=1,4
      DO 20 J=1,4
   20 TC(I,J)=0.D0
C
C
      DO 25 I=1,IST
      DO 25 J=1,IST
   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,4
      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)
C
      GO TO 50
C
C     CREEP LAW NO. 3 **
C
   70 C1=A2*(A2 - 1.)
      C2=A2 - 2.
      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*(2.0*GG*SX*SXY)
      TC(1,4)=COEF*(GG*SX*SZ - XCON2*GAMA)
C
      TC(2,1)=TC(1,2)
      TC(2,2)=COEF*(GG*SY*SY + XCON1*GAMA)
      TC(2,3)=COEF*(2.0*GG*SY*SXY)
      TC(2,4)=COEF*(GG*SY*SZ - XCON2*GAMA)
C
      TC(3,1)=TC(1,3)
      TC(3,2)=TC(2,3)
      TC(3,3)=COEF*(2.0*GG*SXY*SXY + GAMA)
      TC(3,4)=COEF*(2.0*GG*SXY*SZ)
C
      TC(4,1)=TC(1,4)
      TC(4,2)=TC(2,4)
      TC(4,3)=TC(3,4)
      TC(4,4)=COEF*(GG*SZ*SZ + XCON1*GAMA)
C
      DO 55 I=1,4
      DO 55 J=1,4
   55 TTC(I,J)=TC(I,J)/(2.0*XC2)
C
C
   80 IF(IST.EQ.4) GO TO 90
C
C
      DO 85 I=1,2
      DO 85 J=1,3
   85 TC(I,J)=TC(I,J) - XFAC1*TC(4,J)
C
   90 DO 95 I=1,IST
      DO 95 J=1,IST
   95 RH(I)=RH(I) + TC(I,J)*STRSS2(J)
C
      DO 110 J=1,IST
  110 TC(J,J)=TC(J,J) + 1.0
C
C
      CALL EQSOL2(TC,RH,1,IST)
C
      DO 130 J=1,4
  130 STRSS2(J)=RH(J)
C
C
      DO 160 I=1,4
      DO 150 J=1,IST
  150 DPSC(I)=DPSC(I) + TTC(I,J)*(STRSS2(J) - TSTRSS(J))
  160 EPSC2(I)=EPSC1(I) + DPSC(I)
C
      RETURN
C
      END
      SUBROUTINE EQSOL2(A,R,KEY,N)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      DIMENSION A(4,4),R(4),ICOL(4),TR(4),IROW(4)
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 EQSOL2))
 7010 FORMAT(/,10X,15HPIVOT NUMBER = ,I5)
C
      END
      SUBROUTINE EFST(EST,SX,SY,SS,SZ,STRESS)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      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 STRESS(4)
C
      SM=(STRESS(1) + STRESS(2) + STRESS(4))*XCON2
      SX=STRESS(1)-SM
      SY=STRESS(2)-SM
      SS=STRESS(3)
      SZ=STRESS(4)-SM
      EST=SQRT(1.5*(SX*SX + SY*SY + 2.0*SS*SS + SZ*SZ))
C
      RETURN
C
      END
      SUBROUTINE EPMAT2(STRSSM,ALFAM,EPSTRM,DEPS,DPSC,DPST,CEP,XLAMDA,
     1                  PROP1,PROP2,PROPM,YLDM,KEY,A2,B2,C1,C2,DPSP,SXM,
     2                  SYM,SXYM,SZM,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 /MATMOD/ STRESS(4),STRAIN(4),C(4,4),IPT,NEL,IPS,ISVE
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      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 /TPLAS2/ ECAM,ECBM,ECCM,EETD,COEF1,COEF2,PRM,CM,CD,
     1                YSD
C
      DIMENSION STRSSM(*),ALFAM(*),DEPS(*),DPSC(*),PROP1(*),PROP2(*),
     1          PROPM(*),CEP(4,4),DPSP(4)
C
      EQUIVALENCE (NPAR(5),ITYP2D),(NPAR(15),MODEL)
C
C
C
C
C
      SXT=SXM
      SYT=SYM
      SXYT=SXYM
      SZT=SZM
C
      IF(INDEX.GT.1) GO TO 40
C
      YMM=PROPM(1)
      PRM=PROPM(2)
      ETM=PROPM(4)
C
      YMD=PROP2(1) - PROP1(1)