a10.for

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

      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
      IF(ITYP2D.LT.2) GO TO 40
C
      ECBM=1.0 - PRM*PRM
      ECAM=YMM/ECBM
      ECCM=1.0 + PRM*PRM
C
      COEF1=(ECBM*YMD + 2.0*YMM*PRM*PRD)/(YMM*ECBM*ECBM)
      COEF2=(YMM*PRD*ECCM + YMD*PRM*ECBM)/(YMM*ECBM*ECBM)
C
C
C
C
   40 IF(MODEL.EQ.10) GO TO 45
C
      SXM=SXM - ALFAM(1)
      SYM=SYM - ALFAM(2)
      SXYM=SXYM - ALFAM(3)
      SZM=SZM - ALFAM(4)
C
C     PLANE STRESS **
C
   45 IF(ITYP2D.LT.2) GO TO 50
C
      WP1=0.5*ECAM*((SXM + PRM*SYM)*(DEPS(1) - DPSC(1) - DPST) + (SYM +
     1    PRM*SXM)*(DEPS(2) - DPSC(2) - DPST)) + CM*SXYM*(DEPS(3) -
     2    DPSC(3))
C
      WP2=0.5*COEF1*(SXM*STRSSM(1) + SYM*STRSSM(2) - PRM*(SXM*
     1    STRSSM(2) + SYM*STRSSM(1)))
C
      WP2=WP2 + 0.5*COEF2*(SXM*STRSSM(2) + SYM*STRSSM(1) - PRM*(SXM*
     1    STRSSM(1) + SYM*STRSSM(2))) + (CD/CM)*SXYM*STRSSM(3)
C
      DENMP=(2.0*XCON2*XCON2*YLDM*YLDM*EETM) + 0.5*ECAM*(SXM*SXM +
     1       SYM*SYM + 2.0*PRM*SXM*SYM) + 2.0*CM*SXYM*SXYM
C
      GO TO 60
C
C
   50 WP1=CM*(SXM*(DEPS(1) - DPSC(1) - DPST) + SYM*(DEPS(2) - DPSC(2) -
     1    DPST) + SXYM*(DEPS(3) - DPSC(3)) + SZM*(DEPS(4) - DPSC(4) -
     2    DPST))
C
      WP2=(0.5*CD/CM)*(SXM*STRSSM(1) + SYM*STRSSM(2) + 2.0*SXYM*
     1     STRSSM(3) + SZM*STRSSM(4))
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)=2.0*XLAMDA*SXYM
      DPSP(4)=XLAMDA*SZM
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
      SXYM=SXYT
      SZM=SXT
C
      RETURN
C
C
   90 YLD1=YLDM
C
      SX1=SXM
      SY1=SYM
      SXY1=SXYM
      SZ1=SZM
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)=   - GAMA*SXY1
      CEP(1,4)=B2 - GAMA*SZ1
C
      GAMA=GAMA1*SY1
      CEP(2,1)=CEP(1,2)
      CEP(2,2)=A2 - GAMA*SY1
      CEP(2,3)=   - GAMA*SXY1
      CEP(2,4)=B2 - GAMA*SZ1
C
      GAMA=GAMA1*SXY1
      CEP(3,1)=CEP(1,3)
      CEP(3,2)=CEP(2,3)
      CEP(3,3)=C2 - GAMA*SXY1
      CEP(3,4)=   - GAMA*SZ1
C
      GAMA=GAMA1*SZ1
      CEP(4,1)=CEP(1,4)
      CEP(4,2)=CEP(2,4)
      CEP(4,3)=CEP(3,4)
      CEP(4,4)=A2 - GAMA*SZ1
C
      IF(ITYP2D.LT.2) RETURN
C
C
      DO 100 I=1,3
      FAC=CEP(I,4)/CEP(4,4)
      DO 100 J=I,3
      CEP(I,J)=CEP(I,J) - CEP(4,J)*FAC
  100 CEP(J,I)=CEP(I,J)
C
      RETURN
C
      END
      SUBROUTINE EMAT2(TMP,PROP,PROPI,XA1,XB1,C1,D1,E1,F1,KKK)
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
C
      DIMENSION PROP(16,*),PROPI(*)
C
      EQUIVALENCE (NPAR(5),ITYP2D)
C
C
C
C
C
      CALL MTITP2(PROP,TMP,PROPI)
      YM=PROPI(1)
      PR=PROPI(2)
C
C
      A2=YM/(1.0 + PR)
      C1=A2*0.5
      A2=A2/(1.0 - 2.0*PR)
      B2=A2*PR
      A2=A2 - B2
      A1=A2
      B1=B2
      XA1=A1
      XB1=B1
      D1=PR/(PR - 1.0)
      E1=1.0/YM
      F1=-PR*E1
C
      IF(ITYP2D.LT.2) GO TO 30
C
C        PLANE STRESS **
C
   20 A1=YM/(1.0 - PR*PR)
      B1=PR*A1
C
C
   30 IF(KKK.EQ.1) RETURN
C
      DO 40 I=1,4
      DO 40 J=1,4
   40 C(I,J)=0.D0
C
      C(1,1)=A1
      C(1,2)=B1
      C(2,1)=B1
      C(2,2)=A1
      C(3,3)=C1
C
      IF(ITYP2D.GE.2) RETURN
C
      C(1,4)=B1
      C(2,4)=B1
      C(4,1)=B1
      C(4,2)=B1
      C(4,4)=A1
C
      RETURN
C
      END
      SUBROUTINE MTITP2(PROP,TMP,PROPI)
C
C
C
C
C        YOUNGS MODULUS
      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
      COMMON /MATMOD/ STRESS(4),STRAIN(4),C(4,4),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 MTITP2))
C
      END
      SUBROUTINE CREEP2(DDT,DEPSC,TEMPD,EPSC,ORIG,NORG,STRESS,GAMA,
     1                  STRNR,PTIME2,EST,SX,SY,SXY,SZ,F,R,G,FP,INDEX,
     2                  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 /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 DEPSC(4),ORIG(4,*),EPSC(*),STRESS(*)
C
C
C
      IMAX=50
      ETOL1=5.D-3
      ETOL4=5.D-6
      ETOL5=1.D-20
C
C
   10 CALL CYCRP2(ECSTR,EPSC,ORIG,NORG,STRESS,INDEX)
      IF(EST.LE.TOL5) RETURN
C
      IF (KCRP.GE.2) GO TO 20
C
C
      CALL CRPLW2 (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 CRPLW2 (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*2.0*SXY
      DEPSC(4)=C1*SZ
C
      RETURN
C
 2000 FORMAT(//,69H    ERROR   NEWTON ITERATION FAILED TO CONVERGE   (SU
     1BROUTINE CREEP2))
C
      END
      SUBROUTINE CYCRP2(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(4,*)
C
C
C
      IF(INDEX.GT.1) GO TO 50
C
C
C
      CALL EFCSTR(EPSD,ORIG,ORIG,2)
C
C
      DUM=0.D0
      DO 15 I=1,4
   15 DUM=DUM+(EPSC(I)-ORIG(I,NORG))*STRESS(I)
      IF(DUM.GE.0.D0) GO TO 50
C
C
      CALL EFCSTR(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,4
   20 DUM=DUM+(EPSC(I)-ORIG(I,NN))*STRESS(I)
      IF(DUM.GE.0.D0) GO TO 25
C
C
      CALL EFCSTR(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,4
   48 ORIG(I,NORG)=EPSC(I)
C
C
   50 CALL EFCSTR(ECSTR,EPSC,ORIG,NORG)
C
      RETURN
C
      END
      SUBROUTINE CRPLW2 (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 /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
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))
C
      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
      END
      SUBROUTINE EFCSTR(ECSTR,EPSC,ORIG,NORG)
C
C
C
C     CURRENT ORIGIN
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,
     1               SUBDD,RNGL,RNGU,DTT,TOL7,
     3               KCRP,NITE,NALG,IINTP,NPTS,ITCHK,IST,ISR
C
      DIMENSION DEPSC(4),EPSC(*),ORIG(4,*)
C
C
C
      DO 10 I=1,4
   10 DEPSC(I)=EPSC(I)-ORIG(I,NORG)
C
      ECSTR=SQRT(XCON1*(DEPSC(1)*DEPSC(1) + DEPSC(2)*DEPSC(2) +
     1           DEPSC(4)*DEPSC(4)) + XCON2*(DEPSC(3)*DEPSC(3)))
C
      RETURN
C
C*FILE END
      END