a06.for

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

C
      K=3
      SIGP(K)=SIGP(K) + C(K,K)*(EP(K) - TEP(K))
      IST=4
      IF (ITYP2D.GE.2) IST=2
      DO 57 I=1,IST
      DO 57 J=1,IST
      IF (I.EQ.3 .OR. J.EQ.3) GO TO 57
      SIGP(I)=SIGP(I) + C(I,J)*(EP(J) - TEP(J))
   57 CONTINUE
C
      EVI=-(STRAIN(1) + STRAIN(2) + STRAIN(4))
C
      IF (MODEL.EQ.4) GO TO 58
      TEP(1)=EP(1)
      TEP(2)=EP(2)
      TEP(4)=EP(4)
   58 IF (ANGLE.LT.3.61D2)
     1  CALL CRAKID (STRESS,STRAIN,PGRAV,CRKSTR,RKLD,RKUN,GLD,SP33,
     2               ANGLE,EP,NUMCRK,MODEL,3)
C
      IF (ANGLE.LT.3.61D2) ANG=ANGLE
      CALL DCRACK (C,STRESS,ANG,MODEL,ITYP2D,NUMCRK,1,2)
      IF (ANGLE.GT.3.61D2) GO TO 60
C
      GO TO 65
C
C
C
C
   50 EVI=-(STRAIN(1) + STRAIN(2) + STRAIN(4))
      EPSMM=-EVI/3.
      EPS11=STRAIN(1)-EPSMM
      EPS22=STRAIN(2)-EPSMM
      EPS33=STRAIN(4)-EPSMM
      EPS12=STRAIN(3)
C
      PEE=3.*RK*(EPSMM - EMM) + TMM
C
      S11=2.*G*(EPS11-E11)+T11
      S22=2.*G*(EPS22-E22)+T22
      S33=2.*G*(EPS33-E33)+T33
      S12=   G*(EPS12-E12)+T12
C
      STRESS(1)=S11+PEE
      STRESS(2)=S22+PEE
      STRESS(4)=S33+PEE
      STRESS(3)=S12
      IF (ITYP2D .GE. 2) STRESS(4)=0.D0
C
C
C
C
   60 CALL PRNCPL (STRESS,STRAIN,ANG,RKLD,RKUN,GLD,SP33,EP,PGRAV,
     1             MODEL,1)
      IF (MODEL.EQ.4 .AND. ICRACK.LT.1) GO TO 65
C
C
      IF (MODEL.NE.5) GO TO 59
      TEP(1)=EP(1)
      TEP(2)=EP(2)
      TEP(4)=EP(4)
      IF (PGRAV.NE.1.D2) GO TO 59
      NUMCRK=3
      CALL DCRACK (C,STRESS,ANGLE,MODEL,ITYP2D,NUMCRK,2,2)
      FALSTR=SIGCP
      IF (KPRI.EQ.0) GO TO 65
      GO TO 70
C
C
   59 IF (KPRI.EQ.0) GO TO 65
C
      IF (MODEL.GT.4) GO TO 150
      IF (P1 .GT. PGRAV) ANGLE=ANG
      IF (STRESS(4) .GT. PGRAV) ANGLE=ANG - 361.
      IF (ANGLE .GT. 3.61D2) GO TO 65
      NUMCRK=0
      IF (P1 .GT. PGRAV) NUMCRK=1
      IF (P2.GT.PGRAV) NUMCRK=2
      IF (NUMCRK.EQ.2) ANGLE=-ANGLE
      CRKSTR(1)=EP(1)
      CRKSTR(2)=EP(2)
      CRKSTR(3)=EP(4)
      IF (STRESS(4) .GT. PGRAV) NUMCRK=NUMCRK + 4
      IF (NUMCRK.EQ.5) ANGLE=ANG - 722.
      IF (NUMCRK.EQ.6) ANGLE=-ANG - 361.
C
      CALL DCRACK (C,STRESS,ANGLE,MODEL,ITYP2D,NUMCRK,1,2)
      ILFSET=1
C
C
C
C
   65 IF (KPRI.NE.0) GO TO 70
      IF (IPRI.NE.0 .OR. IPS.EQ.0) GO TO 66
      IF (IPT.NE.1) GO TO 66
      IF (IPRNT.EQ.1) write(66,2000)
      write(66,2035) NEL
   66 IF (MODEL.EQ.4 .AND. ICRACK.LT.1) GO TO 67
      IF (ANGLE.LT.3.61D2) GO TO 165
      IF (MODEL.GT.4) GO TO 150
      PTOT=-PGRAV + P1
      IF (PTOT.GT.0.D0) ANGPRI=ANG
      IF (STRESS(4) .GT. PGRAV) ANGPRI=ANG
      IF (P1.GT.PGRAV) NUMCRK=1
      IF (P2.GT.PGRAV) NUMCRK=2
      IF (STRESS(4).GT.PGRAV) NUMCRK=NUMCRK + 4
      GO TO 160
C
C
  150 IF (P1.LT.0.D0 .AND. P3.LT.0.D0) GO TO 155
      IF (P1.GT.FALSTR) NUMCRK=1
      IF (P2.GT.FALSTR) NUMCRK=2
      IF (P3.GT.FALSTR) NUMCRK=NUMCRK + 4
  155 IF (P2.LT.SIGCP .OR. P3.LT.SIGCP) NUMCRK=3
      IF (NUMCRK.GT.0) ANGPRI=ANG
      IF (NUMCRK.EQ.3) FALSTR=SIGCP
C
      IF (KPRI.EQ.0) GO TO 160
      IF (NUMCRK.EQ.0) GO TO 70
      ANGLE=ANG
      IF (NUMCRK.EQ.2) ANGLE=-ANGLE
      IF (NUMCRK.GE.4) ANGLE=ANGLE - 361.
      IF (NUMCRK.EQ.5) ANGLE=ANG - 722.
      CRKSTR(1)=EP(1)
      CRKSTR(2)=EP(2)
      CRKSTR(3)=EP(4)
C
      IF (NUMCRK.NE.3) GO TO 159
      PGRAV=1.D2
      CRKSTR(1)=SIGCP
      CRKSTR(2)=EPSCP
C
  159 CALL DCRACK (C,STRESS,ANGLE,MODEL,ITYP2D,NUMCRK,1,2)
      ILFSET=1
      GO TO 70
C
  160 IF (ANGPRI.GT.3.61D2) GO TO 67
      CALL DCRACK (C,STRESS,ANGPRI,MODEL,ITYP2D,NUMCRK,1,2)
      GO TO 166
C
  165 ANGPRI=ANGLE
      IF (ANGPRI.LT.-5.41D2) ANGPRI=ANGPRI + 722.
      IF (ANGPRI .LT. (-1.8D2)) ANGPRI=ANGPRI + 361.
      IF (ANGPRI .LT. 0.D0) ANGPRI=-ANGPRI
      IF (ANGPRI.GT.1.8D2) ANGPRI=ANGPRI - 180.
C
  166 NF=NUMCRK + 1
      IPHCRA=NUMCRK
      PHANGL=ANGPRI
      IF (INDNL .EQ. 2) CALL CAUCHY
      IF (IPRI.NE.0 .OR. IPS.EQ.0) GO TO 170
      IF (NF .GT. 5) NF=NF - 3
      IF (NUMCRK .EQ. 6) GO TO 68
      GO TO (61,62,63,64,62) ,NF
   61 write(66,2040) IPT,STRESS(4),(STRESS(I),I=1,3),P1,P2,ANGPRI,
     1               FALSTR
      GO TO 72
   62 write(66,2041) IPT,STRESS(4),(STRESS(I),I=1,3),P1,P2,ANGPRI,
     1               FALSTR
      GO TO 72
   63 write(66,2042) IPT,STRESS(4),(STRESS(I),I=1,3),P1,P2,ANGPRI,
     1               FALSTR
      GO TO 72
   64 write(66,2043) IPT,STRESS(4),(STRESS(I),I=1,3),P1,P2,ANGPRI,
     1               FALSTR
      GO TO 72
   68 write(66,2045) IPT,STRESS(4),(STRESS(I),I=1,3),P1,P2,ANGPRI,
     1               FALSTR
      GO TO 72
   67 IF (INDNL .EQ. 2) CALL CAUCHY
      IPHCRA=10
      PHANGL=ANG
      IF (IPRI.NE.0 .OR. IPS.EQ.0) GO TO 170
      write(66,2044) IPT,STRESS(4),(STRESS(I),I=1,3),P1,P2,ANG,
     1               FALSTR
   72 IF (MODEL.EQ.5 .AND. ITHERM.EQ.1) write(66,2046) EVGRAV
C
C
  170 IF (JNPORT.EQ.0 .OR. KPLOTE.NE.0) GO TO 70
      IF (ISVE.EQ.0) GO TO 70
      IF (JNPORT.EQ.1)
     1   WRITE (IBPORT     ) 'OUTPUT-2',NEL,IPT,IPHCRA,
     2                       (STRESS(I),I=1,4),(STRAIN(I),I=1,4),P1,P2,
     3                       PHANGL,FALSTR,EVGRAV
      IF (JNPORT.EQ.2)
     1   WRITE (IFPORT,9000) 'OUTPUT-2',NEL,IPT,IPHCRA,
     2                       (STRESS(I),I=1,4),(STRAIN(I),I=1,4),P1,P2,
     3                       PHANGL,FALSTR,EVGRAV
C
 9000 FORMAT ( A,/,3I10,/,(4E20.13) )
C
C
C
C
   70 SIG(1)=STRESS(1)
      SIG(2)=STRESS(2)
      SIG(3)=STRESS(3)
      SIG(4)=STRESS(4)
C
      EPS(1)=STRAIN(1)
      EPS(2)=STRAIN(2)
      EPS(3)=STRAIN(3)
      EPS(4)=STRAIN(4)
      IF (KPRI .EQ. 0) RETURN
      IF (ICOUNT.EQ.3) RETURN
      IF (MODEL.EQ.4) GO TO 71
      IF (PGRAV.NE.1.D2) GO TO 69
      E=0.D0
      GO TO 100
C
   69 SBAR=((SIG(1)-SIG(2))**2+(SIG(1)-SIG(4))**2+(SIG(2)-SIG(4))**2)/6.
      EVI=SQRT(SBAR + SIG(3)**2) + ALFA*(SIG(1) + SIG(2) + SIG(4))
      IF (IKAS.EQ.1 .AND. ILFSET.EQ.1) EVMAX=EVI
   71 IF (EVI.GT.EVMAX) EVMAX=EVI
C
C
C
C
C
      IF (IEQREF.EQ.1) GO TO 85
C
      IF (EVI.LT.EVMAX) GO TO 75
C
C
   74 IF (MODEL.EQ.4) EVTOT=EVMAX + EVGRAV
      IK=J
      IKAS=1
      ITEM=2
      GO TO 5
C
   75 IF (IKAS) 100,85,85
C
C
   85 IF (MODEL.EQ.4) GO TO 90
      E=EKK(1)
      VNU=EKK(2)
      MOD45=1
      GO TO 100
C
   90 RKUNLO=RKUN(I) + RATIO * (RKUN(J) - RKUN(I))
      G=G*(RKUNLO/RK)
      RK=RKUNLO
C
C
  100 IF (MODEL.EQ.4) GO TO 95
      IF (E.LE.0.D0) E=EKK(1)*STIFAC
      RK=E/(3.*(1. - 2.*VNU))
      G=E/(2.*(1. + VNU))
      IF (PGRAV.EQ.1.D2) GO TO 101
      IF (MOD45.EQ.2) GO TO 110
   95 IF (MODEL.EQ.4 .AND. ICRACK.LT.1) GO TO 101
      IF (ANGLE.LT.3.61D2) GO TO 110
C
C
  101 DUM=0.6666667 * G
      A1=RK + 2.*DUM
      B1=RK - DUM
C
      C(1,1)=A1
      C(1,2)=B1
      C(1,3)=0.D0
      C(2,1)=B1
      C(2,2)=A1
      C(2,3)=0.D0
      C(3,1)=0.D0
      C(3,2)=0.D0
      C(3,3)=G
C
C
      C(1,4)=B1
      C(2,4)=B1
      C(3,4)=0.D0
      C(4,1)=B1
      C(4,2)=B1
      C(4,3)=0.D0
      C(4,4)=A1
C
      IF (ITYP2D.LT.2) GO TO 200
C
      DO 105 I=1,2
      A=C(I,4)/C(4,4)
      DO 105 J=I,2
      C(I,J)=C(I,J) - C(4,J)*A
  105 C(J,I)=C(I,J)
C
      GO TO 200
C
C
  110 IF (ANGLE.LT.3.61D2) ANG=ANGLE
      CALL DCRACK (C,STRESS,ANG,MODEL,ITYP2D,NUMCRK,2,1)
C
  200 IF (ITHERM.GT.0 .AND. IEQUIT.EQ.1)
     1CALL THERM2 (NODS,NOD5,EVGRAV,TEMPV2,YZ,EKK,ITYP2D,2)
C
      IF (IUPDT.EQ.0) RETURN
      DO 210 I=1,15
  210 WA(I)=DUMWA(I)
      RETURN
C
C
 2000 FORMAT (1X,7HELEMENT,7X,9HSTRESS-XX,5X,9HSTRESS-YY,5X,
     1        9HSTRESS-ZZ,5X,9HSTRESS-YZ,4X,10H  SIGMA-P1,4X,
     2        10H  SIGMA-P2,2X,5HANGLE,4X,7HFAILSTR/,1X,7HNUM/IPT,
     3        97X,7H(PGRAV),/)
 2005 FORMAT(49H **STOP - PLANE STRESS CANNOT BE USED WITH CDMOD  )
 2015 FORMAT (///37H * * *     STOP OF SOLUTION     * * *//,
     1        4X,22HFOR ELEMENT GROUP NO =,I2,13H ELEMENT NO =,I5/
     2        4X,29HTHE CURRENT VOLUMETRIC STRAIN,E14.6/
     3        4X,25HEXCEEDS THE TABLE MAXIMUM,E14.6//,
     4           37H * * *   END OF ERROR MESSAGE   * * *///)
 2035 FORMAT (I4)
 2040 FORMAT (5X,I3,2X,6E14.6,F7.2,E11.3,5X,8H  CLOSED)
 2041 FORMAT (5X,I3,2X,6E14.6,F7.2,E11.3,5X,8H1 CRACK )
 2042 FORMAT (5X,I3,2X,6E14.6,F7.2,E11.3,5X,8H2 CRACKS)
 2043 FORMAT (5X,I3,2X,6E14.6,F7.2,E11.3,5X,8H CRUSHED)
 2044 FORMAT (5X,I3,2X,6E14.6,F7.2,E11.3,5X,8HNO CRACK)
 2045 FORMAT (5X,I3,2X,6E14.6,F7.2,E11.3,5X,8H3 CRACKS)
 2046 FORMAT (15X,14HTEMPERATURE = ,E14.6)
C
      END
      SUBROUTINE PRNCPL (STR,EPS,ANG,SP1,SP31,SP32,SP33,EPSL,PGRAV,
     1                   MODEL,KKK)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /CRACK/ GAMMA,STIFAC,SHEFAC,SIGMAT,SIGMAC,EPSC,SIGMAU,EPSU,
     1               RAM5,RBM5,RCM5,ICRACK,MOD45
      COMMON /CONCR2/ BETA,GAMA,RKAPA,ALFA,SIGP(4),TEP(4),EP(4),YP(3),
     1                E,VNU,RK,G,E12,E14,E24,EPSCP,SIGCP,FALSTR,ILFSET
      COMMON /MATMOD/ STRESS(4),STRAIN(4),D(4,4),IPT,NEL,IPS,ISVE
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,IREF,
     1             IEQUIT,IPRI,KPLOTN,KPLOTE
C
      DIMENSION STR(4),EPS(4),SP1(*),SP31(*),SP32(*),SP33(*),EPSL(4),
     1          SIG(3)
      IF (KKK.GT.1) GO TO 10
C
C
      AA=(STR(1) + STR(2))*0.5
      BB=(STR(1) - STR(2))*0.5
      CC=SQRT(BB*BB + STR(3)*STR(3))
      SIGP(1)=AA + CC
      SIGP(2)=AA - CC
      SIGP(3)=0.D0
      SIGP(4)=STR(4)
      ANG=4.5D1
      IF (STR(3).EQ.0.D0) ANG=0.1D-3
      IF (ABS(BB).LT.0.1D-6) GO TO 10
      DUM=ABS(STR(3)/BB)
      ANG=57.296*ATAN(DUM)
C
      IF (BB.GT.0.D0 .AND. STR(3).GT.0.D0) ANG=ANG
      IF (BB.LT.0.D0 .AND. STR(3).GT.0.D0) ANG=180. - ANG
      IF (BB.LT.0.D0 .AND. STR(3).LE.0.D0) ANG=180. + ANG
      IF (BB.GT.0.D0 .AND. STR(3).LE.0.D0) ANG=360. - ANG
      ANG=ANG/2.
C
C
   10 PI=4.D0*ATAN(1.D0)
      GAM=ANG*PI/180.
      SG=SIN(GAM)
      CG=COS(GAM)
C
      EPSL(1)=EPS(1)*CG*CG + EPS(2)*SG*SG + EPS(3)*SG*CG
      EPSL(2)=EPS(1)*SG*SG + EPS(2)*CG*CG - EPS(3)*SG*CG
      EPSL(3)=0.D0
      EPSL(4)=EPS(4)
      IF (KKK.EQ.2) RETURN
      IF (MODEL.EQ.4) RETURN
      IF (PGRAV.EQ.1.D2) RETURN
      SIGCP=SIGMAC
      EPSCP=EPSC
      FALSTR=SIGMAT
C
C
      SIG(1)=SIGP(1)
      SIG(2)=SIGP(2)
      SIG(3)=SIGP(4)
  114 IS=0
      DO 116 I=1,2
      IF (SIG(I + 1).LE.SIG(I)) GO TO 116
      IS=IS + 1
      TEMP=SIG(I + 1)
      SIG(I + 1)=SIG(I)
      SIG(I)=TEMP
  116 CONTINUE
      IF (IS.GT.0) GO TO 114
      P1=SIG(1)
      P2=SIG(2)
      P3=SIG(3)
C
C
      IF (P3.GE.0.D0) GO TO 180
      IF (P1.LT.0.D0) GO TO 135
      IF (P2.LT.0.D0) GO TO 130
      FALSTR=SIGMAT*(1. - P3/SIGMAC)
      GO TO 180
C
  130 SP31I=SP31(1)
      SP32I=SP32(1)
      SP33I=SP33(1)
      TEMP=0.D0
      GO TO 150
C
  135 TEMP=P1/SIGMAC
      DO 140 I=2,6
      J=I - 1
      IF (TEMP.LT.SP1(I)) GO TO 145
  140 CONTINUE
      write(66,3000) NG,NEL,IPT,P1
      write(66,3100) P1,P2,P3
      STOP
C
  145 DSP=SP1(I) - SP1(J)
      DSPI=TEMP - SP1(J)
      FRAC=DSPI/DSP
      SP31I=SP31(J) + FRAC*(SP31(I) - SP31(J))
      SP32I=SP32(J) + FRAC*(SP32(I) - SP32(J))
      SP33I=SP33(J) + FRAC*(SP33(I) - SP33(J))
C
  150 RATIO=P2/SIGMAC
      IF (RATIO.GT.BETA*SP32I) GO TO 160
      SLOPE=(SP32I - SP31I)/(BETA*SP32I - TEMP)
      SIGCP=SP31I*SIGMAC + SLOPE*(P2 - P1)
      IF (P1.GT.0.D0) SIGCP=SP31I*SIGMAC + SLOPE*P2
      GO TO 170
  160 SLOPE=(SP33I - SP32I)/(SP33I - BETA*SP32I)
      SIGCP=SP32I*SIGMAC + SLOPE*(P2 - BETA*SP32I*SIGMAC)
C
  170 IF (SIGCP.LT.SIGMAC) EPSCP=GAMA*EPSC*SIGCP/SIGMAC
      IF (P1.GE.0.D0) FALSTR=SIGMAT*(1. - P2/SIGCP)*(1. - P3/SIGCP)
      IF (FALSTR.LT.0.001*SIGMAT) FALSTR=SIGMAT
      IF (P1.LT.0.D0) FALSTR=SIGCP
C
  180 RETURN
C
 3000 FORMAT (//37H * * *     STOP OF SOLUTION     * * *//,
     1        4X,22HFOR ELEMENT GROUP NO =,I2,13H ELEMENT NO =,I5,
     2        9H AT IPT =,I5/
     3        4X,41HTHE CURRENT VALUE OF PRINCIPAL STRESS 1 =,E15.5,/
     4        4X,49HIS LARGER THAN THE MAXIMUM INPUT VALUE FOR SP1(6))
 3100 FORMAT (//4X,40HTHE CURRENT PRINCIPAL/CRACK STRESSES ARE,3E15.5/,
     1          4X,35HERROR OCCURED IN SUBROUTINE PRNCPL.  //,
     2          37H * * *   END OF ERROR MESSAGE   * * *///)
C
      END
      SUBROUTINE CRAKID (STR,EPS,PGRAV,CRKSTR,SP1,SP31,SP32,SP33,
     1                   ANGLE,EPSL,NUMCRK,MODEL,KKK)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /CRACK/ GAMMA,STIFAC,SHEFAC,SIGMAT,SIGMAC,EPSC,SIGMAU,EPSU,
     1               RAM5,RBM5,RCM5,ICRACK,MOD45
      COMMON /CONCR2/ BETA,GAMA,RKAPA,ALFA,SIGP(4),TEP(4),EP(4),YP(3),
     1                E,VNU,RK,G,E12,E14,E24,EPSCP,SIGCP,FALSTR,ILFSET
      COMMON /MATMOD/ STRESS(4),STRAIN(4),D(4,4),IPT,NEL,IPS,ISVE
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,IREF,
     1             IEQUIT,IPRI,KPLOTN,KPLOTE