a20.for

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

      SXT=SXM
      SYT=SYM
      SZT=SZM
      SXYT=SXYM
      SXZT=SXZM
      SYZT=SYZM
C
      IF(INDEX.GT.1) GO TO 40
C
      YMM=PROPM(1)
      PRM=PROPM(2)
      ETM=PROPM(4)
C
      YMD=PROP2(1) - PROP1(1)
      PRD=PROP2(2) - PROP1(2)
      ETD=PROP2(4) - PROP1(4)
      YSD=PROP2(3) - PROP1(3)
C
      CD=C2 - C1
      CM=0.5*YMM/(1.0 + PRM)
      EETD=((YMM*YMM*ETD) - (ETM*ETM*YMD))/((YMM - ETM)*(YMM - ETM))
C
C
C
   40 IF(MODEL.EQ.10) GO TO 50
C
      SXM=SXM - ALFAM(1)
      SYM=SYM - ALFAM(2)
      SZM=SZM - ALFAM(3)
      SXYM=SXYM - ALFAM(4)
      SXZM=SXZM - ALFAM(5)
      SYZM=SYZM - ALFAM(6)
C
   50 WP1=CM*(SXM*(DEPS(1) - DPSC(1) - DPST) + SYM*(DEPS(2) - DPSC(2) -
     1   DPST) + SZM*(DEPS(3) - DPSC(3) - DPST) + SXYM*(DEPS(4) -
     2   DPSC(4)) + SXZM*(DEPS(5) - DPSC(5)) + SYZM*(DEPS(6) -
     3   DPSC(6)))

C
      WP2=(0.5*CD/CM)*(SXM*STRSSM(1) + SYM*STRSSM(2) + SZM*STRSSM(3)
     1     + 2.0*(SXYM*STRSSM(4) + SXZM*STRSSM(5) + SYZM*STRSSM(6)))
C
      DENMP=(XCON1*YLDM*YLDM)*(CM + EETM*XCON2)
C
   60 WP=WP1 + WP2
C
      IF(MODEL.EQ.11) GO TO 65
C
C
      XLAMDA=(WP - (YLDM*XCON2)*(EPSTRM*EETD + YSD))/DENMP
      WPP=XLAMDA
      GO TO 70
C
C
   65 XLAMDA=(WP - (YLDM*YSD*XCON2))/DENMP
      WPP=XLAMDA
C
C
   70 IF (XLAMDA.GT.0.D0) GO TO 75
      XLAMDA=0.D0
      GO TO 80
C
   75 IF(KEY.EQ.2 .AND. IEQREF.NE.1) XLAMDA=XLAMDA - WP1/DENMP
C
   80 DPSP(1)=XLAMDA*SXM
      DPSP(2)=XLAMDA*SYM
      DPSP(3)=XLAMDA*SZM
      DPSP(4)=2.0*XLAMDA*SXYM
      DPSP(5)=2.0*XLAMDA*SXZM
      DPSP(6)=2.0*XLAMDA*SYZM
C
      XLAMDA=WPP
C
      IF (KEY.EQ.2 .AND. IEQREF.NE.1 .AND. XLAMDA.GE.0.D0) GO TO 90
C
      SXM=SXT
      SYM=SYT
      SZM=SZT
      SXYM=SXYT
      SXZM=SXZT
      SYZM=SYZT
C
      RETURN
C
C
   90 YLD1=YLDM
C
      SX1=SXM
      SY1=SYM
      SZ1=SZM
      SXY1=SXYM
      SXZ1=SXZM
      SYZ1=SYZM
C
      YIELD=YLD1*YLD1/(3.0*C2)
      GAMA1=1.0/(YIELD*(1.0 + EETM/(3.0*CM)))
C
      GAMA=GAMA1*SX1
      CEP(1,1)=A2 - GAMA*SX1
      CEP(1,2)=B2 - GAMA*SY1
      CEP(1,3)=B2 - GAMA*SZ1
      CEP(1,4)=   - GAMA*SXY1
      CEP(1,5)=   - GAMA*SXZ1
      CEP(1,6)=   - GAMA*SYZ1
C
      GAMA=GAMA1*SY1
      CEP(2,1)=CEP(1,2)
      CEP(2,2)=A2 - GAMA*SY1
      CEP(2,3)=B2 - GAMA*SZ1
      CEP(2,4)=   - GAMA*SXY1
      CEP(2,5)=   - GAMA*SXZ1
      CEP(2,6)=   - GAMA*SYZ1
C
      GAMA=GAMA1*SZ1
      CEP(3,1)=CEP(1,3)
      CEP(3,2)=CEP(2,3)
      CEP(3,3)=A2 - GAMA*SZ1
      CEP(3,4)=   - GAMA*SXY1
      CEP(3,5)=   - GAMA*SXZ1
      CEP(3,6)=   - GAMA*SYZ1
C
      GAMA=GAMA1*SXY1
      CEP(4,1)=CEP(1,4)
      CEP(4,2)=CEP(2,4)
      CEP(4,3)=CEP(3,4)
      CEP(4,4)=C2 - GAMA*SXY1
      CEP(4,5)=   - GAMA*SXZ1
      CEP(4,6)=   - GAMA*SYZ1
C
      GAMA=GAMA1*SXZ1
      CEP(5,1)=CEP(1,5)
      CEP(5,2)=CEP(2,5)
      CEP(5,3)=CEP(3,5)
      CEP(5,4)=CEP(4,5)
      CEP(5,5)=C2 - GAMA*SXZ1
      CEP(5,6)=   - GAMA*SYZ1
C
      GAMA=GAMA1*SYZ1
      CEP(6,1)=CEP(1,6)
      CEP(6,2)=CEP(2,6)
      CEP(6,3)=CEP(3,6)
      CEP(6,4)=CEP(4,6)
      CEP(6,5)=CEP(5,6)
      CEP(6,6)=C2 - GAMA*SYZ1
C
      RETURN
C
      END
      SUBROUTINE EMAT3(TMP,PROP,PROPI,A1,B1,C1,D1,E1,F1,KKK)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /MTMD3D/ C(6,6),STRESS(6),STRAIN(6),IPT,NEL,IPS,ISVE
C
      DIMENSION PROP(16,*),PROPI(*)
C
C
C
C
C
      CALL MTITP3(PROP,TMP,PROPI)
      YM=PROPI(1)
      PR=PROPI(2)
C
C
      A1=YM/(1.0 + PR)
      C1=A1*0.5
      A1=A1/(1.0 - 2.0*PR)
      B1=A1*PR
      A1=A1 - B1
      D1=PR/(PR - 1.0)
      E1=1.0/YM
      F1=-PR*E1
C
C
   30 IF(KKK.EQ.1) RETURN
C
      DO 40 I=1,6
      DO 40 J=1,6
   40 C(I,J)=0.D0
C
      C(1,1)=A1
      C(1,2)=B1
      C(1,3)=B1
      C(2,1)=B1
      C(2,2)=A1
      C(2,3)=B1
      C(3,1)=B1
      C(3,2)=B1
      C(3,3)=A1
      C(4,4)=C1
      C(5,5)=C1
      C(6,6)=C1
C
      RETURN
C
      END
      SUBROUTINE MTITP3(PROP,TMP,PROPI)
C
C
C
C
C        YOUNGS MODULUS
      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 PROP(16,*),PROPI(*)
C
C
C
C
C
    5 IF(TMP.GE.RNGL) GO TO 10
      write(66,3001)
      STOP
C
   10 L=0
      DO 20 K=2,NPTS
      L=L + 1
      DUM=PROP(K,1)
      IF(K.EQ.NPTS) DUM=RNGU
      IF(TMP.GT.DUM) GO TO 20
      GO TO 25
   20 CONTINUE
      write(66,3001)
      STOP
C
   25 XRATIO=(TMP - PROP(L,1))/(PROP(L + 1,1) - PROP(L,1))
C
C
      IF (XRATIO.GT.1.D0) XRATIO=1.D0
      IF (XRATIO.LT.0.D0) XRATIO=0.D0
C
C
C
C
      DO 30 M=2,6
   30 PROPI(M - 1)=PROP(L,M) + XRATIO*(PROP(L + 1,M) - PROP(L,M))
C
      RETURN
C
 3001 FORMAT(//,92H    ERROR   TEMPERATURE OUTSIDE RANGE OF MATERIAL PRO
     1PERTY TEMPERATURES  (SUBROUTINE MTITP3))
C
      END
      SUBROUTINE CREEP3(DDT,DEPSC,TEMPD,EPSC,ORIG,NORG,STRESS,GAMA,
     1                  STRNR,PTIME2,EST,SX,SY,SZ,SXY,SXZ,SYZ,F,R,G,FP,
     2                  INDEX,ECSTR)
C
C
C
C     CYCLIC CREEP
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 /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      COMMON /ADINAI/ OPVAR(7),TSTART,IRINT,ISTOTE
      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 DEPSC(6),ORIG(6,*),EPSC(*),STRESS(*)
C
C
C
      IMAX=50
      ETOL1=5.D-3
      ETOL4=5.D-6
      ETOL5=1.D-20
C
C
   10 CALL CYCRP3(ECSTR,EPSC,ORIG,NORG,STRESS,INDEX)
      IF(EST.LE.TOL5) RETURN
C
      IF (KCRP.GE.2) GO TO 20
C
C
      CALL CRPLW3 (EST,ECSTR,STRN,STRNR,DDT,TEMPD,F,R,G,FP)
      GO TO 60
C
C
C
C
C
   20 PTIME1=FLOAT(KSTEP)*DTT + TSTART
      IF (KCRP.EQ.3) PTIME1=PTIME1 + 0.5*DDT
      IF (ECSTR.LE.ETOL5 .AND. KCRP.EQ.3) PTIME1=1.D-10
      PTIME2=PTIME1
C
      KOUNT=1
   25 CALL CRPLW3 (EST,ECSTR,STRN,STRNR,PTIME2,TEMPD,F,R,G,FP)
C
      IMOD=0
      FUNCT=STRN-ECSTR
      DELTA=FUNCT/STRNR
      IF (ECSTR.LE.ETOL5 .AND. KCRP.EQ.3) GO TO 60
      IF (ECSTR.EQ.0.D0) DELTA=PTIME2
C
C
      IF ((PTIME2 - DELTA).GE.0.D0) GO TO 30
      DELTA=0.5*PTIME2
      IMOD=1
C
   30 IF(KOUNT.GT.1) GO TO 40
      PTIME2=PTIME1 - DELTA
      DNORM1=ABS(DELTA)
      KOUNT=KOUNT + 1
      GO TO 25
C
C
   40 DNORM2=ABS(DELTA)
      IF(IMOD.EQ.1) GO TO 50
      IF(DNORM2.LE.DNORM1) GO TO 45
C
C     TOLERANCE BAND **
C
      XTOL=ETOL4*PTIME1
      IF(PTIME1.LE.ETOL5) XTOL=ETOL5
      IF(DNORM2.LE.XTOL.AND.DNORM1.LE.XTOL) GO TO 60
      GO TO 50
C
   45 XTOL=ETOL1*PTIME1
      IF(PTIME1.LE.ETOL5) XTOL=ETOL5
      IF(DNORM1.LE.XTOL) GO TO 60
C
C     NO CONVERGENCE *
C
   50 KOUNT=KOUNT + 1
      IF(KOUNT.LE.IMAX) GO TO 55
      write(66,2000)
      STOP
C
   55 PTIME1=PTIME2
      PTIME2=PTIME2 - DELTA
      DNORM1=DNORM2
      GO TO 25
C
C
   60 GAMA=1.5*STRNR/EST
      C1=GAMA*DDT
C
  100 DEPSC(1)=C1*SX
      DEPSC(2)=C1*SY
      DEPSC(3)=C1*SZ
      DEPSC(4)=C1*2.0*SXY
      DEPSC(5)=C1*2.0*SXZ
      DEPSC(6)=C1*2.0*SYZ
C
      RETURN
C
 2000 FORMAT(//,69H    ERROR   NEWTON ITERATION FAILED TO CONVERGE   (SU
     1BROUTINE CREEP3))
C
      END
      SUBROUTINE CYCRP3(ECSTR,EPSC,ORIG,NORG,STRESS,INDEX)
C
C
C
C
C
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
C
      DIMENSION STRESS(*),EPSC(*),ORIG(6,*)
C
C
C
      IF(INDEX.GT.1) GO TO 50
C
C
C
      CALL EFCST3(EPSD,ORIG,ORIG,2)
C
C
      DUM=0.D0
      DO 15 I=1,6
   15 DUM=DUM+(EPSC(I)-ORIG(I,NORG))*STRESS(I)
      IF (DUM.GE.0.D0) GO TO 50
C
C
      CALL EFCST3(ECSTR,EPSC,ORIG,NORG)
C
C     NEW VALUES **
C
      IF(ECSTR.GT.EPSD) GO TO 40
C
C
      IF(NORG.EQ.2) GO TO 18
   17 NN=2
      GO TO 19
   18 NN=1
   19 DUM=0.D0
C
      DO 20 I=1,6
   20 DUM=DUM+(EPSC(I)-ORIG(I,NN))*STRESS(I)
      IF (DUM.GE.0.D0) GO TO 25
C
C
      CALL EFCST3(TECSTR,EPSC,ORIG,NN)
      IF(ECSTR.GE.TECSTR) RETURN
C
C
   25 IF(NORG.EQ.2) GO TO 35
   30 NORG=2
      GO TO 50
   35 NORG=1
      GO TO 50
C
C
   40 IF(NORG.EQ.2) GO TO 42
   41 NORG=2
      GO TO 45
   42 NORG=1
C
   45 DO 48 I=1,6
   48 ORIG(I,NORG)=EPSC(I)
C
C
   50 CALL EFCST3(ECSTR,EPSC,ORIG,NORG)
C
      RETURN
C
      END
      SUBROUTINE EFCST3(ECSTR,EPSC,ORIG,NORG)
C
C
C
C     CURRENT ORIGIN
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 DEPSC(6),EPSC(*),ORIG(6,*)
C
C
C
      DO 10 I=1,6
   10 DEPSC(I)=EPSC(I)-ORIG(I,NORG)
C
      ECSTR=SQRT(XCON1*(DEPSC(1)*DEPSC(1) + DEPSC(2)*DEPSC(2) +
     1           DEPSC(3)*DEPSC(3)) + XCON2*(DEPSC(4)*DEPSC(4)
     2           +DEPSC(5)*DEPSC(5)+DEPSC(6)*DEPSC(6)))
C
      RETURN
C
      END
      SUBROUTINE CRPLW3 (STRESS,ECSTR,STRAIN,STRNR,TIME,TEMPD,F,R,G,FP)
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
C
C
      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) 10,50,60
C
C
   10 IF (A2.GE.1.D0) GO TO 20
C
      RTTOL=2.D1
      EX1=1./(1.-A2)
      EX2=A1*EX1
      EX3=A2*EX1
      EX4=-RTTOL*EX3
C
C
      ECMIN=(A0**EX1)*(STRESS**EX2)*(A2**EX3)*(10.0**EX4)
      IF(ECSTR.LE.ECMIN) GO TO 40
C
   20 IF (ECSTR.EQ.0.D0) GO TO 40
      EX5=1.0/A2
      EX6=A1*EX5
      EX7=(A2-1.)/A2
      EX8=ECSTR**EX7
C
   30 STRNR=(A0**EX5)*(STRESS**EX6)*A2*EX8
C
      RETURN
C
   40 STRAIN=A0*(STRESS**A1)*(TIME**A2)
      STRNR=STRAIN/TIME
C
      RETURN
C
C
   50 F=A0*EXP(A1*STRESS)
      R=A2*((STRESS/A3)**A4)
      G=A5*EXP(A6*STRESS)
      STRAIN=F*(1.-EXP(-R*TIME)) + (G*TIME)
      STRNR=F*R*EXP(-R*TIME) + G
C
      RETURN
C
C
   60 F=A0*STRESS**A1
      R=TIME**A2 + A3*TIME**A4 + A5*TIME**A6
      G=EXP(-A7/(TEMPD +273.16))
      STRAIN=F*R*G
C
      A2M1=A2 - 1.
      A4M1=A4 - 1.
      A6M1=A6 - 1.
      A34 =A3*A4
      A56 =A5*A6
C
      FP=A2*TIME**A2M1 + A34*TIME**A4M1 + A56*TIME**A6M1
C
      STRNR=F*FP*G
C
      RETURN
C
C*FILE END
      END