a06.for

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

C
      DIMENSION STR(4),EPS(4),CRKSTR(3),SP1(*),SP31(*),SP32(*),SP33(*),
     1          EPSL(4),SIG(3)
C
C
      PI=4.D0*ATAN(1.D0)
      TANG=ANGLE
      IF (TANG.LT.-5.41D2) TANG=TANG + 722.
      IF (TANG .LT. (-1.8D2)) TANG=TANG + 361.
      IF (TANG.GE.1.8D2) TANG=TANG - 180.
      GAM=2.*ABS(TANG)*PI/180.
      SG=SIN(GAM)
      CG=COS(GAM)
      IF (KKK.EQ.2) GO TO 10
      IF (KKK.EQ.3) GO TO 107
C
      R11=(STR(1) + STR(2))*0.5
      R12=(STR(1) - STR(2))*0.5
      SIGP(1)=R11 + R12*CG + STR(3)*SG
      SIGP(2)=R11 - R12*CG - STR(3)*SG
      SIGP(3)=-R12*SG + STR(3)*CG
      SIGP(4)=STR(4)
C
   10 D11=(EPS(1) + EPS(2))*0.5
      D12=(EPS(1) - EPS(2))*0.5
      EPSL(1)=D11 + D12*CG + EPS(3)*SG*0.5
      EPSL(2)=D11 - D12*CG - EPS(3)*SG*0.5
      EPSL(3)=-D12*SG + EPS(3)*CG
      EPSL(4)=EPS(4)
      IF (KKK.EQ.2) RETURN
C
C
  107 NUMCRK=1
      IF (TANG.LT.0.D0) NUMCRK=2
      IF (ANGLE.GE.1.8D2) NUMCRK=0
      IF (ANGLE.LT.-1.81D2) NUMCRK=4
      IF (ANGLE.LT.-3.61D2) NUMCRK=6
      IF (ANGLE.LT.-5.41D2) NUMCRK=5
      NCROLD=NUMCRK
      FALSTR=PGRAV
      IF (MODEL.EQ.4) GO TO 109
      FALSTR=SIGMAT
      SIGCP=SIGMAC
      EPSCP=EPSC
      SIG(1)=SIGP(1)
      SIG(2)=SIGP(2)
      SIG(3)=SIGP(4)
C
  109 NF=NUMCRK + 1
      GO TO (110, 111, 180, 185, 195, 190, 190), NF
  110 IF (SIGP(1).LT.FALSTR .AND. SIGP(2).LT.FALSTR .AND. SIGP(4).LT.
     1    FALSTR) GO TO 112
      IF (SIGP(1).GT.FALSTR .OR. SIGP(2).GT.FALSTR) NUMCRK=1
      ANGLE=TANG
      CRKSTR(1)=EPSL(1)
      CRKSTR(2)=EPSL(2)
      CRKSTR(3)=EPSL(4)
      IF (SIGP(4).LT.FALSTR) GO TO 112
      NUMCRK=NUMCRK + 4
      ANGLE=ANGLE - 361.
      IF (NUMCRK.EQ.5) ANGLE=ANGLE - 361.
      GO TO 112
C
  111 IF (EPSL(1).LT.0.D0 .AND. EPSL(1).LT.CRKSTR(1)) NUMCRK=0
      IF (NUMCRK.EQ.0) ANGLE=ANGLE + 180.
      IF (SIGP(2).LT.FALSTR) GO TO 113
      NUMCRK=2
      CRKSTR(2)=EPSL(2)
      ANGLE=-ANGLE
  113 IF (SIGP(4).LT.FALSTR) GO TO 112
      NUMCRK=NUMCRK + 4
      CRKSTR(3)=EPSL(4)
      ANGLE=ANGLE - 361.
      IF (NUMCRK.EQ.5) ANGLE=ANGLE - 361.
C
C
  112 IF (MODEL.EQ.4) RETURN
C
C
  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
      IF (P3.GE.0.D0) GO TO 200
      IF (P2.LT.0.D0) GO TO 115
      IF (P3.GT.SIGMAC) GO TO 200
      GO TO 185
  115 IF (P1.LT.0.D0) GO TO 135
      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)
  170 IF (SIGCP.LT.SIGMAC) EPSCP=GAMA*EPSC*SIGCP/SIGMAC
      IF (P3.LE.SIGCP) GO TO 185
      GO TO 200
C
  180 IF (EPSL(1).LT.0.D0 .AND. EPSL(1).LT.CRKSTR(1)) NUMCRK=1
      IF (EPSL(2).LT.0.D0 .AND. EPSL(2).LT.CRKSTR(2)) NUMCRK=NUMCRK - 1
      IF (NUMCRK.LT.2) ANGLE=ABS(ANGLE)
      IF (NUMCRK.EQ.0) ANGLE=ANGLE + 180.
      IF (SIGP(4).LT.FALSTR) GO TO 112
      NUMCRK=NUMCRK + 4
      CRKSTR(3)=EPSL(4)
      ANGLE=ANGLE - 361.
      IF (NUMCRK.EQ.5) ANGLE=ANGLE - 361.
      GO TO 112
C
  185 NUMCRK=3
      PGRAV=1.D2
      FALSTR=SIGCP
      CRKSTR(1)=SIGCP
      CRKSTR(2)=EPSCP
      GO TO 200
C
  190 IF (EPSL(2).LT.0.D0 .AND. EPSL(2).LT.CRKSTR(2) .AND. NUMCRK.EQ.6)
     1           NUMCRK=5
      IF (EPSL(1).LT.0.D0 .AND. EPSL(1).LT.CRKSTR(1) .AND. NUMCRK.EQ.5)
     1           NUMCRK=4
      IF (EPSL(4).LT.0.D0 .AND. EPSL(4).LT.CRKSTR(3)) NUMCRK=NUMCRK - 4
      IF (NCROLD.NE.5) GO TO 112
      IF (SIGP(2).GE.FALSTR) NUMCRK=6
      IF (NUMCRK.EQ.6) ANGLE=-(ANGLE + 722.) - 361.
      GO TO 112
C
  195 CONTINUE
      R11=(SIGP(1) + SIGP(2))*0.5
      R12=(SIGP(1) - SIGP(2))*0.5
      RAD=SQRT(R12*R12 + SIGP(3)*SIGP(3))
      SIG(1)=R11 + RAD
      SIG(2)=R11 - RAD
      IF (MODEL.EQ.5 .AND. SIG(2).LT.0.D0) FALSTR=SIGMAT*(1. - SIG(2)
     1                                                          /SIGMAC)
      IF (SIG(1).GT.FALSTR) NUMCRK=NUMCRK + 1
      IF (SIG(2).GT.FALSTR) NUMCRK=NUMCRK + 1
      IF (NUMCRK.EQ.4) GO TO 198
C
      S11=(EPSL(1) + EPSL(2))*0.5
      S12=(EPSL(1) - EPSL(2))*0.5
      RAD=SQRT(S12*S12 + EPSL(3)*EPSL(3)/4.)
      EPS2=S11 + RAD
      EPS3=S11 - RAD
      IF (NUMCRK.GT.4) CRKSTR(1)=EPS2
      IF (NUMCRK.GT.4) CRKSTR(2)=EPS3
      SIGP(1)=SIG(1)
      SIGP(2)=SIG(2)
      EPSL(1)=EPS2
      EPSL(2)=EPS3
      EPSL(3)=0.D0
      ANG=4.5D1
      IF (ABS(R12).GT.0.1D-5) ANG=ATAN2(SIGP(3),R12)*28.648
      ANG=ANG + TANG
      IF (ANG.LT.0.D0) ANG=ANG + 180.
      IF (NUMCRK.EQ.6) ANG=-ANG
      ANGLE=ANG
      SIGP(3)=0.D0
  198 IF (EPSL(4).LT.0.D0 .AND. EPSL(4).LT.CRKSTR(3)) NUMCRK=NUMCRK - 4
      IF (NUMCRK.EQ.5) ANGLE=ANG - 722.
      IF (NUMCRK.EQ.6) ANGLE=ANG - 361.
      GO TO 112
C
  200 IF (KKK.EQ.3 .AND. NCROLD.NE.NUMCRK) ILFSET=1
      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 CRAKID.  //,
     2          37H * * *   END OF ERROR MESSAGE   * * *///)
      END
      SUBROUTINE DCRACK (C,SIG,ANGLE,MODEL,ITYP2D,NUMCRK,ILOCAL,KKK)
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
      DIMENSION C(4,4),SIG(4),D(4,4),T(4,4),DSIG(4)
C
      IF (KKK.EQ.2) GO TO 42
C
C
      IF (MOD45.EQ.2) GO TO 8
      DUM=2.0*G/3.0
      A1=RK + 2.*DUM
      B1=RK - DUM
C
      C(1,1)=A1
      C(1,2)=B1
      C(1,3)=0.D0
      C(1,4)=B1
      C(2,2)=A1
      C(2,3)=0.D0
      C(2,4)=B1
      C(3,3)=G
      C(3,4)=0.D0
      C(4,4)=A1
      GO TO 10
C
C
    8 A1=(1. + VNU)*(1. - 2.*VNU)
      B1=(1. - VNU)/A1
      D1=VNU/A1
      C1=1./(2.*(1. + VNU))
      C(1,1)=B1*YP(1)
      C(1,2)=E12*D1
      C(1,3)=0.D0
      C(1,4)=E14*D1
      C(2,2)=B1*YP(2)
      C(2,3)=0.D0
      C(2,4)=E24*D1
      C(3,3)=C1*E12
      C(3,4)=0.D0
      C(4,4)=B1*YP(3)
C
   10 IF (NUMCRK.EQ.0 .OR. NUMCRK.EQ.3) GO TO 35
C
C
      C(3,3)=C(3,3)*SHEFAC
      IF (MOD45.EQ.2) GO TO 13
      A2=4.*G*(RK+G/3.)/(RK+4.0*G/3.)
      B2=2.*G*(RK-2.*G/3.)/(RK+4.0*G/3.)
      C2=A2
      D2=(9.*RK*G)/(3.*RK + G)
      R2=A2
      S2=B2
      T2=A2
      W2=G
      Z2=D2
      GO TO 14
C
   13 A1=1./(1. - VNU*VNU)
      B1=VNU*A1
      A2=A1*YP(2)
      B2=B1*E24
      C2=A1*YP(3)
      D2=YP(3)
      R2=A1*YP(1)
      S2=B1*E12
      T2=A2
      W2=C1*E12
      Z2=E12
C
   14 GO TO (15, 20, 35, 25, 30, 33), NUMCRK
C
   15 DO 16 I=1,4
   16 C(1,I)=C(1,I)*STIFAC
      C(2,2)=A2
      C(2,4)=B2
      C(4,4)=C2
      GO TO 35
C
   20 DO 21 I=1,2
      DO 21 J=I,4
   21 C(I,J)=C(I,J)*STIFAC
      C(4,4)=D2
      GO TO 35
C
   25 DO 26 I=1,4
   26 C(I,4)=C(I,4)*STIFAC
      C(1,1)=R2
      C(1,2)=S2
      C(2,2)=T2
      C(3,3)=W2
      GO TO 35
C
   30 DO 31 I=1,4
      DO 31 J=I,4
   31 C(I,J)=C(I,J)*STIFAC
      C(2,2)=Z2
      C(3,3)=W2*SHEFAC
      GO TO 35
C
   33 DO 34 I=1,4
      DO 34 J=I,4
   34 C(I,J)=C(I,J)*STIFAC
      C(3,3)=W2*SHEFAC
C
   35 DO 40 I=1,3
      DO 40 J=I,4
   40 C(J,I)=C(I,J)
      IF (ILOCAL.EQ.1) GO TO 90
C
   42 IF (NUMCRK.NE.3) GO TO 43
      IF (KKK.EQ.1) GO TO 90
      IF (ILOCAL.EQ.2) GO TO 97
C
C
   43 PI=4.D0*ATAN(1.D0)
      TANGLE=ANGLE
      IF (TANGLE.LT.-5.41D2) TANGLE=TANGLE + 722.
      IF (TANGLE .LT. (-1.8D2)) TANGLE=TANGLE + 361.
      IF (TANGLE.GT.1.8D2) TANGLE=TANGLE - 180.
      GAM=ABS(TANGLE)*PI/180.
      SG = SIN(GAM)
      CG = COS(GAM)
      T(1,1) =  CG**2
      T(1,2) =  SG**2
      T(1,3) =  CG* SG
      T(1,4) =  0.D0
      T(2,1) =  T(1,2)
      T(2,2) =  T(1,1)
      T(2,3) = -T(1,3)
      T(2,4) =  0.D0
      T(3,1) =  T(2,3)* 2.0
      T(3,2) = -T(3,1)
      T(3,3) =  T(1,1)- T(1,2)
      T(3,4) =  0.D0
      T(4,1) =  0.D0
      T(4,2) =  0.D0
      T(4,3) =  0.D0
      T(4,4) =  1.D0
      IF (KKK.EQ.2) GO TO 97
C
C
C
      DO 60 IR=1,4
      DO 60 IC=1,4
      D(IR,IC) = 0.D0
      DO 50 IN=1,4
   50 D(IR,IC) = D(IR,IC) + T(IN,IR)* C(IN,IC)
   60 CONTINUE
C
C
      DO 80 IR=1,4
      DO 80 IC=IR,4
      C(IR,IC) = 0.D0
      DO 70 IN=1,4
   70 C(IR,IC) = C(IR,IC) + D(IR,IN)* T(IN,IC)
   80 C(IC,IR)=C(IR,IC)
C
C
C
   90 IF (ITYP2D.LT.2) GO TO 97
      DO 95 I=1,3
      A=C(I,4)/C(4,4)
      DO 95 J=I,3
      C(I,J)=C(I,J) - C(4,J)*A
   95 C(J,I)=C(I,J)
C
   97 IF (KKK.LT.2) RETURN
      IF (MODEL.EQ.4 .AND. ICRACK.LT.2) GO TO 140
      ETATAU=0.D0
      IF (SHEFAC .GT. 0.1D-2) ETATAU=1.D0
C
C
      DO 91 I=1,4
   91 DSIG(I)=0.D0
      IF (NUMCRK.NE.3) GO TO 98
      EPSUP=EPSU*EPSCP/EPSC
      IF (TEP(1).GT.EPSUP .AND. TEP(2).GT.EPSUP .AND. TEP(4).GT.EPSUP)
     1    GO TO 93
      DO 92 I=1,4
   92 SIGP(I)=0.D0
      GO TO 98
C
   93 IF (ILOCAL.EQ.1) GO TO 140
      RETURN
C
C
   98 NF=NUMCRK + 1
      GO TO (140,120,110,155,100,100,100), NF
  100 SIGP(4)=0.D0
      IF (NUMCRK - 5) 140,120,110
  110 SIGP(2)=0.D0
  120 SIGP(3)=SIGP(3)*ETATAU
      SIGP(1)=0.D0
C
C
  140 DO 150 IR=1,4
      DO 150 IC=1,4
  150 DSIG(IR)=DSIG(IR) + T(IC,IR)*SIGP(IC)
C
  155 DO 160 I=1,4
  160 SIG(I)=DSIG(I)
C
      RETURN
C
      END
      SUBROUTINE THERM2 (NODS,NOD5,EVGRAV,TEMPV2,YZ,
     1                   PROP,ITYP2D,KKK)
C
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /GAUSS/ XG(6,6),WGT(6,6),EVAL2(9,2),EVAL3(27,3),R,S,T
      COMMON /MATMOD/ STRESS(4),STRAIN(4),C(4,4),IPT,NEL,IPS,ISVE
      COMMON /TODIM/ BET,THIC,DE,IEL,NND5
C
      DIMENSION NODS(*),NOD5(*),TEMPV2(*),YZ(*),PROP(*)
      DIMENSION H(9),P(2,9),XJ(2,2)
C
      IF (KKK.EQ.2) GO TO 200
C
      DSTR=PROP(3)*(EVGRAV-PROP(39))
      DO 100 I=1,2
  100 STRAIN(I)=STRAIN(I) - DSTR
      IF (ITYP2D.LT.2) STRAIN(4)=STRAIN(4) - DSTR
      RETURN
C
C
C
C
  200 CALL FUNCT2 (R,S,H,P,NOD5,XJ,DET,YZ,NEL,1)
      TEMP2=0.D0
      DO 220 K=1,IEL
      KK=NODS(K)
  220 TEMP2=TEMP2 + H(K)*TEMPV2(KK)
      DELT=TEMP2 - EVGRAV
      DSTR=PROP(3)*DELT
      DO 250 I=1,2
      DO 250 J=1,2
  250 STRESS(I)=STRESS(I) - C(I,J)*DSTR
      EVGRAV=TEMP2
      IF (ITYP2D.GE.2) RETURN
      STRESS(1)=STRESS(1) - C(1,4)*DSTR
      STRESS(2)=STRESS(2) - C(2,4)*DSTR
      STRESS(4)=STRESS(4) - (C(4,1)+C(4,2)+C(4,4))*DSTR
      RETURN
C*FILE END
      END