a20.for

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

  232 DELT=DTT - TAU
C
C     LARGE ENOUGH
C
  233 IF(ISUB + 1.LT.ISUBM.OR.TAU + DELT.GE.TCHK) GO TO 235
C
      write(66,3007)
      write(66,3011) NEL,IPT,ISUBM,TAU,DELT
      STOP
C
C
  235 ISUB=ISUB + 1
      INDEX=1
      YS1=YS2
      EET1=EET2
      YLD1=YLD2
      TMP1=TMP2
C
      DO 240 J=1,6
      EPS1(J)=EPS2(J)
      STRSS1(J)=STRSS2(J)
  240 EPSC1(J)=EPSC2(J)
C
      GO TO 60
C
C
  245 YLDC=YS2
      IF(MODEL.EQ.10) YLDC=EET2*EPSTR2 + YLDC
C
      GO TO 440
C
C
C            FOLLOW
C
C
  250 CALL EFST3(EST1,SX1,SY1,SZ1,SXY1,SXZ1,SYZ1,STRSS1)
C
C
      IF(MODEL.EQ.10) GO TO 255
C
      SX1=SX1 - ALFA1(1)
      SY1=SY1 - ALFA1(2)
      SZ1=SZ1 - ALFA1(3)
      SXY1=SXY1 - ALFA1(4)
      SXZ1=SXZ1 - ALFA1(5)
      SYZ1=SYZ1 - ALFA1(6)
C
  255 RB=SX1*DX + SY1*DY + SZ1*DZ + 2.0*(SXY1*DXY + SXZ1*DXZ + SYZ1*DYZ)
      RD=SX1*SX1 + SY1*SY1 + SZ1*SZ1 + 2.0*(SXY1*SXY1 + SXZ1*SXZ1 +
     1   SYZ1*SYZ1)
      RE=RB - (XCON1*YLD1*DYLD)
      RF=RA - (XCON1*DYLD*DYLD)
      RG=RD - (XCON1*YLD1*YLD1)
C
      IF (IPELD.EQ.2 .OR. RG.GE.0.D0) GO TO 270
C
C
  260 IF(ABS(RF).GT.TOL5) GO TO 265
      RATIO=-RG/(2.0*RE)
      GO TO 275
C
  265 RATIO=(-RE + SQRT(RE*RE - RF*RG))/RF
      GO TO 275
C
C     OF SUBDIVISION *
C
  270 RATIO=0.D0
      IF (RF.GT.TOL5.AND.RE.LT.0.D0) RATIO=-2.0*RE/RF
C
C
  275 IF(RATIO.GE.(-TOL6) .AND. RATIO.LE.(1.0 + TOL6)) GO TO 280
C
      write(66,3012)
      write(66,3013) NEL,IPT,ISUB,TAU,IPELD,RA,RB,RD,RE,RF,RG,RATIO
      STOP
C
  280 IF (RATIO.GT.1.D0) RATIO=1.D0
      IF (RATIO.LT.0.D0) RATIO=0.D0
      IPELD=2
C
C
      TAU=TAU + RATIO*DELT
      TMP1=TMPOLD + DELTMP*TAU/DTT
      YLD1=YLD1 + RATIO*DYLD
      IF(RATIO.GT.TOL5) ISUB=ISUB + 1
      IF(ISUB.GT.ISUBM) ISUBM=ISUBM + 1
C
C
      DO 285 J=1,6
      EPS1(J)=EPS1(J) + RATIO*DEPS(J)
      EPSC1(J)=EPSC1(J) + RATIO*DPSC(J)
      EPSC2(J)=EPSC1(J)
      STRSS1(J)=STRSS1(J) + RATIO*DELSIG(J)
  285 STRSS2(J)=STRSS1(J)
C
C
  288 CALL EMAT3(TMP1,PROP,PROP1,A1,B1,C1,D1,E1,F1,1)
C
C
      ALPHA1=PROP1(5)
      EPST1(1)=ALPHA1*(TMP1 - TREF)
C
C
  290 XNWDT=DTT - TAU
      DELT=XNWDT/FLOAT(ISUBM - ISUB + 1)
      INDEX=1
      DEPSTR=0.D0
      DESTR=1.D0
      DENOM=0.D0
C
C
  292 TMP2=TMPOLD + DELTMP*(TAU + DELT)/DTT
      TMPM=XPARM1*TMP1 + XPARM2*TMP2
C
C
      CALL EMAT3(TMP2,PROP,PROP2,A2,B2,C2,D2,E2,F2,2)
C
      DO 295 J=1,5
  295 PROPM(J)=XPARM1*PROP1(J) + XPARM2*PROP2(J)
C
      YMM=PROPM(1)
      ETM=PROPM(4)
C
      EETM=YMM*ETM/(YMM - ETM)
      CEM=XCON1*EETM
C
C
      ALPHA2=PROP2(5)
C
      EPST2(1)=ALPHA2*(TMP2 - TREF)
      EPST2(2)=EPST2(1)
      EPST2(3)=EPST2(1)
C
      DPST=EPST2(1) - EPST1(1)
C
C
  300 XFAC=(TAU + DELT)/DTT
      DO 305 J=1,6
      EPS2(J)=EPS(J) + XFAC*DELEPS(J)
  305 DEPS(J)=EPS2(J) - EPS1(J)
C
C
  310 DO 315 J=1,6
  315 STRSSM(J)=XPARM1*STRSS1(J) + XPARM2*STRSS2(J)
C
      CALL EFST3(ESTM,SXM,SYM,SZM,SXYM,SXZM,SYZM,STRSSM)
C
C
      EPSTRM=XPARM1*EPSTR1 + XPARM2*EPSTR2
C
C
      IF(MODEL.EQ.10) GO TO 320
C
      DO 318 J=1,6
  318 ALFAM(J)=XPARM1*ALFA1(J) + XPARM2*ALFA2(J)
      GO TO 324
C
C
C
  320 YLDM=ESTM
      GO TO 328
C
C
  324 DO 325 J=1,6
  325 DSTSS(J)=STRSSM(J) - ALFAM(J)
C
      CALL EFST3(YLDM,SXT,SYT,SZT,SXYT,SXZT,SYZT,DSTSS)
C
C
  328 IF(KCRP.EQ.0) GO TO 335
C
      DO 330 J=1,6
  330 DPSC(J)=0.D0
      CRSRM=0.D0
C
      DO 332 J=1,6
  332 EPSCM(J)=XPARM1*EPSC1(J) + XPARM2*EPSC2(J)
C
      CALL CREEP3(DELT,DPSC,TMPM,EPSCM,ORIGD,NORGD,STRSSM,
     1            GAMA,CRSRM,PTIME,ESTM,SXM,SYM,SZM,SXYM,SXZM,SYZM,
     2            FF,RR,GG,FP,INDEX,ECSTRM)
C
      IF(INDEX.EQ.1) ECSTR1=ECSTRM
C
C
  335 CALL EPMAT3(STRSSM,ALFAM,EPSTRM,DEPS,DPSC,DPST,CEP,XLAMDA,
     1            PROP1,PROP2,PROPM,YLDM,1,A2,B2,C1,C2,DPSP,SXM,SYM,SZM,
     2            SXYM,SXZM,SYZM,INDEX,EETM)
C
  342 DO 344 J=1,6
  344 EPSP2(J)=EPSP1(J) + DPSP(J)
C
      IF (XLAMDA.LE.0.D0) GO TO 345
C
C
      DEPSTR=SQRT(XCON1*(DPSP(1)*DPSP(1) + DPSP(2)*DPSP(2) + DPSP(3)*
     1            DPSP(3)) + XCON2*(DPSP(4)*DPSP(4) + DPSP(5)*DPSP(5) +
     2            DPSP(6)*DPSP(6)))
      EPSTR2=EPSTR1 + DEPSTR
C
C
      IF(MODEL.EQ.10) GO TO 345
C
      ALFA2(1)=ALFA1(1) + CEM*DPSP(1)
      ALFA2(2)=ALFA1(2) + CEM*DPSP(2)
      ALFA2(3)=ALFA1(3) + CEM*DPSP(3)
      ALFA2(4)=ALFA1(4) + 0.5*CEM*DPSP(4)
      ALFA2(5)=ALFA1(5) + 0.5*CEM*DPSP(5)
      ALFA2(6)=ALFA1(6) + 0.5*CEM*DPSP(6)
C
C
  345 CALL SIGMA3(STRSS2,EPS2,EPSP2,EPST2,EPSC1,EPSC2,DPSC,GAMA,PTIME,
     1            CRSRM,FF,RR,GG,FP,ESTM,SXM,SYM,SZM,SXYM,SXZM,SYZM,
     2            DELT,B2,C2,D2)
C
C
C
C
  350 IF (XPARM2.EQ.0.D0) GO TO 400
      IF(ITCHK.EQ.1) GO TO 358
C
      INDEX=INDEX + 1
      IF(INDEX.LE.NITE) GO TO 310
      GO TO 400
C
C     STRESS VECTOR **
C
  358 IF(INDEX - 4) 360,366,362
C
  360 INDEX=INDEX + 1
      GO TO 310
C
  362 DNORM2=0.D0
      DO 365 J=1,6
  365 DNORM2=DNORM2 + (STRSS2(J) - STRSSD(J))*(STRSS2(J) - STRSSD(J))
C
C
  366 SNORM=0.D0
      DO 368 J=1,6
  368 SNORM=SNORM + STRSS2(J)*STRSS2(J)
C
C
      IF(INDEX.GT.5) GO TO 375
      SNORM2=SNORM
      IF(INDEX.EQ.4) SNORM1=SNORM2
      IF(INDEX.EQ.5) DNORM1=DNORM2
      INDEX=INDEX + 1
C
      DO 370 J=1,6
  370 STRSSD(J)=STRSS2(J)
      GO TO 310
C
C
C
C
  375 IF(DNORM2.LE.DNORM1) GO TO 390
C
C     TOLERANCE BAND
C
      XTOL=TOL3*SNORM1
      IF(SNORM1.LE.TOL2) XTOL=TOL2
      IF(DNORM1.LE.XTOL.AND.DNORM2.LE.XTOL) GO TO 400
C
C     (NALG .EQ. 2) *
C
      DELT=DELT*(SQRT(DNORM1/DNORM2))/SUBDD
      IF(NALG.EQ.2.AND.ISUB.LT.ISUBM) GO TO 380
C
  378 write(66,3008)
      write(66,3002) NEL,IPT,ISUB,TAU,DELT
      STOP
C
C
  380 INDEX=1
      DEPSTR=0.D0
      DESTR=1.D0
      DENOM=0.D0
      EPSTR2=EPSTR1
C
      DO 385 I=1,6
      STRSS2(I)=STRSS1(I)
      EPSP2(I)=EPSP1(I)
      ALFA2(I)=ALFA1(I)
  385 EPSC2(I)=EPSC1(I)
C
      GO TO 292
C
C
  390 XTOL=TOL1*SNORM1
      IF(SNORM1.LE.TOL2) XTOL=TOL2
      IF(DNORM1.LE.XTOL) GO TO 400
C
C     NO CONVERGENCE
C
      INDEX=INDEX + 1
      IF(INDEX.LE.NITE) GO TO 395
C
      write(66,3003)
      write(66,3011) NEL,IPT,ISUB,TAU,DELT
      STOP
C
  395 DNORM1=DNORM2
      SNORM1=SNORM2
      SNORM2=SNORM
C
      DO 398 J=1,6
  398 STRSSD(J)=STRSS2(J)
      GO TO 310
C
C
C
C
C
  400 IF (NALG.EQ.1) GO TO 410
      DECSTR=CRSRM*DELT
      DESTR=DECSTR + DEPSTR
      DENOM=ECSTR1 + EPSTR1
      IF(DESTR.LE.TOL5.OR.DENOM.LE.TOL5) GO TO 410
C
      CHECK=DESTR/(DENOM*TOLPC)
      IF (CHECK.LE.1.1D0) GO TO 410
C
C
      DELT=DELT/CHECK
      IF(NALG.EQ.2.AND.ISUB.LT.ISUBM) GO TO 405
C
      write(66,3009)
      write(66,3002) NEL,IPT,ISUB,TAU,DELT
      STOP
C
C
  405 INDEX=1
      DEPSTR=0.D0
      DESTR=1.D0
      DENOM=0.D0
      EPSTR2=EPSTR1
C
      DO 408 I=1,6
      STRSS2(I)=STRSS1(I)
      EPSP2(I)=EPSP1(I)
      ALFA2(I)=ALFA1(I)
  408 EPSC2(I)=EPSC1(I)
C
      GO TO 292
C
C
  410 YM2=PROP2(1)
      ET2=PROP2(4)
      YS2=PROP2(3)
C
      EET2=YM2*ET2/(YM2 - ET2)
C
C     PLASTIC STRAIN **
C
C
      YLDC=YS2
C
C
      IF(MODEL.EQ.10) YLDC=EET2*EPSTR2 + YLDC
C
C
C
      IF(MODEL.EQ.11) GO TO 412
      CALL EFST3(YLD2,SX2,SY2,SZ2,SXY2,SXZ2,SYZ2,STRSS2)
      GO TO 414
C
C
  412 DO 413 J=1,6
  413 DSTSS(J)=STRSS2(J) - ALFA2(J)
C
      CALL EFST3(YLD2,SXT,SYT,SZT,SXYT,SXZT,SYZT,DSTSS)
C
C
C
  414 CHECK = ABS(YLD2 - YLDC) / MIN(YLD2,YLDC)
      IF(CHECK.LE.TOL7) GO TO 415
C
      write(66,3014)
      write(66,3015) NEL,IPT,ISUB,TAU,YLD2,YLDC
      STOP
C
C         HARDENING ***
C
C
C
  415 TAU=TAU + DELT
C
C
      IF(NALG.EQ.2) GO TO 416
C
C     NALG .EQ. 1 *
C
      IF(ISUB.EQ.ISUBM) GO TO 440
      GO TO 425
C
C     NALG .EQ. 2 *
C
  416 IF(TAU.GE.TCHK) GO TO 440
      IF(DESTR.LE.TOL5) GO TO 420
C
      DELT=DELT*TOLPC*(1.0 + (DENOM/DESTR))
      IF(TAU + DELT.GE.TCHK) DELT=DTT - TAU
      GO TO 422
C
  420 DELT=DTT - TAU
C
C     LARGE ENOUGH
C
  422 IF(ISUB + 1.LT.ISUBM.OR.TAU + DELT.GE.TCHK) GO TO 425
C
      write(66,3010)
      write(66,3011) NEL,IPT,ISUBM,TAU,DELT
      STOP
C
C
  425 ISUB=ISUB + 1
      INDEX=1
      DEPSTR=0.D0
      DESTR=1.D0
      DENOM=0.D0
      TMP1=TMP2
      EPSTR1=EPSTR2
      EPST1(1)=EPST2(1)
C
      A1=A2
      B1=B2
      C1=C2
      D1=D2
      E1=E2
      F1=F2
C
      DO 430 J=1,6
      EPS1(J)=EPS2(J)
      STRSS1(J)=STRSS2(J)
      EPSC1(J)=EPSC2(J)
      ALFA1(J)=ALFA2(J)
  430 EPSP1(J)=EPSP2(J)
C
      DO 435 J=1,5
  435 PROP1(J)=PROP2(J)
C
      GO TO 292
C
C
  440 IF(IUPDT.NE.0) GO TO 455
C
      DO 445 J=1,6
      SIG(J)=STRSS2(J)
      ALFA(J)=ALFA2(J)
      EPSC(J)=EPSC2(J)
      EPSP(J)=EPSP2(J)
  445 EPS(J)=STRAIN(J)
C
      YLD=YLD2
      EPSTR=EPSTR2
      TMPOLD=TMP2
      IPEL=IPELD
      NORG=NORGD
C
      DO 450 I=1,6
      DO 450 J=1,2
  450 ORIG(I,J)=ORIGD(I,J)
C
C
C
  455 IF(ICOUNT.EQ.3) RETURN
C
C
      IF(KPRI.EQ.0) GO TO 600
C
C
      IF (IEQUIT.EQ.1) GO TO 474
      GO TO (471,472) IPELD
C
  471 CALL EMAT3 (TEMP2,PROP,PROP2,A2,B2,C2,D2,E2,F2,2)
      RETURN
C
  472 DO 473 I=1,6
      DEPS(I) = EPS2(I)  - EPS1(I)
      DPSC(I) = EPSC2(I) - EPSC1(I)
      DPSP(I) = EPSP2(I) - EPSP1(I)
      STRSSM(I) = XPARM1*STRSS1(I) + XPARM2*STRSS2(I)
      IF (MODEL.NE.11) GO TO 473
      ALFAM(I)  = XPARM1*ALFA1(I)  + XPARM2*ALFA2(I)
      DSTSS(I)  = STRSSM(I) - ALFAM(I)
  473 CONTINUE
C
      DO 476 I=1,5
  476 PROPM(I) = XPARM1*PROP1(I) + XPARM2*PROP2(I)
C
      YMM = PROPM(1)
      ETM = PROPM(4)
      EETM = YMM*ETM/(YMM-ETM)
C
      DPST = EPST2(1) - EPST1(1)
      EPSTRM = XPARM1*EPSTR1 + XPARM2*EPSTR2
C
      CALL EFST3 (YLDM,SXM,SYM,SZM,SXYM,SXZM,SYZM,STRSSM)
      IF (MODEL.EQ.11)
     1CALL EFST3 (YLDM,SXT,SYT,SZT,SXYT,SXZT,SYZT,DSTSS)
C
      CALL EPMAT3 (STRSSM,ALFAM,EPSTRM,DEPS,DPSC,DPST,CEP,XLAMDA,
     1             PROP1,PROP2,PROPM,YLDM,2,A2,B2,C1,C2,DPSP,SXM,SYM,
     1             SZM,SXYM,SXZM,SYZM,1,EETM)
C
      IF (XLAMDA.GE.0.0D+00) GO TO 477
      CALL EMAT3 (TEMP2,PROP,PROP2,A2,B2,C2,D2,E2,F2,2)
      RETURN
C
  477 DO 478 I=1,6
      DO 478 J=1,6
  478 C(I,J) = CEP(I,J)
      RETURN
C
  474 CONTINUE
C
C
C
C
      DO 458 J=1,6
      EPSC1(J)=EPSC2(J)
      EPSP1(J)=EPSP2(J)
      EPST1(J)=EPST2(J)
      STRSS1(J)=STRSS2(J)
  458 ALFA1(J)=ALFA2(J)
C
      YLD1=YLD2
      EPSTR1=EPSTR2
C
C
      A1=A2
      B1=B2
      C1=C2
      D1=D2
      E1=E2
      F1=F2
C
      DO 460 I=1,5
  460 PROP1(I)=PROP2(I)
C
      CALL EMAT3(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,3