Sources Esope | Cast3M

projvm
C PROJVM    SOURCE    GG250959  17/09/20    21:16:15     9554           
      SUBROUTINE PROJVM(SIG,EPS,EPST,EPSTAR,DEPS,PENTE,SN,
     . PANTIN,SELAS,ITER,SO,SSTAR,SI,NK,IBOU,COEF,JKI,
     . KERRE,ecou,necou)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
C         CALCULE L ECOULEMENT EN CAS DE CRITERE DE VON MISES
C
      DIMENSION SIG(*),EPS(*),COEF(*)
C
      SEGMENT ECOU
C**     COMMON/ECOU/TEST,ALFAH,
       REAL*8 TEST, ALFAH,
     1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
     2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),W1(6),W2(6),
     1 DALPHA(6),DSIGO(6),E(12),XINV(3),
     2 SIPLAD(6),DSIGP0(6),TETM,TETI
      ENDSEGMENT
C      COMMON/ECOU/TEST,ALFAH,
C     1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
C     2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),W1(6),W2(6),
C     1 DALPHA(6),DSIGO(6),E(12),XINV(3),
C     2 SIPLAD(6),DSIGP0(6),TETM,TETI
C
       SEGMENT NECOU
C      COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
       INTEGER NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
     . ITYP,IFOUR,IFLUAG,
     . ICINE,ITHER,IFLUPL,ICYCL,IBI,
     . JFLUAG,KFLUAG,LFLUAG,
     . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
      ENDSEGMENT
C      COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
C     . ITYP,IFOUR,IFLUAG,
C     . ICINE,ITHER,IFLUPL,ICYCL,IBI,
C     . JFLUAG,KFLUAG,LFLUAG,
C     . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
C
      DATA ITMAX/200/
C
      JEBOUC=0
      PENTEJ=PENTE
      DEPSJ=DEPS
2020  JEBOUC=JEBOUC+1
      PENTE=PENTEJ
      DEPS=DEPSJ
      IF(JEBOUC.GT.1) WRITE(6,77881) PENTE,DEPS
77881 FORMAT('0 VALEURS INITIALES  PENTE=',1PE16.5,2X,
     . 'DEPS=',1PE16.5/)
C
C
C          CAS PLASTIQUE
      PENTE1=PENTE
      EPST=EPSTAR+DEPS
C CONTRAINTE SUR LA COURBE DE TRACTION CORRESPONDANTE
      IF(ICINE.EQ.1) GO TO 81
C         CAS ISOTROPE
      CALL TRACTI(SN,EPST,SIG,EPS,NCOURB,2,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GO TO 57
      ENDIF
C
C         CAS MIXTE CALCUL DU NOUVEAU RAYON DE LA SPHERE
      IF(ICINE.EQ.2) SN=SN-PANTIN*EPST
      GO TO 82
C           CAS CINEMATIQUE
   81 SN=SELAS
   82 CONTINUE
C CONTRAINTE MOYENNE
      SM=0.5*(SN+SELAS)
C 2222 CONTINUE
C
C ITERATIONS INTERNES CALCUL DE DELTA EPSILON
C
      ITER=0
      DEPS0=DEPS
      DEPSI=DEPS0
      SO=SN
      DEPS00=DEPS0
      JYCOMP=0
      UNSSM=1./SM
      SIGT(1)=SIGEL(1)+DSIGP(1)
      SI=SSTAR
      IF(JEBOUC.GT.1) WRITE(6,77882) SN,SM,SI
77882 FORMAT('0 SN =',1PE16.5,2X,'SM=',1PE16.5,2X,
     . 'SI=',1PE12.5/)
 555  ITER=ITER+1
      IF(JEBOUC.GT.1) WRITE(6,77883) ITER
77883 FORMAT('0 ----->  ITER=',I3/)
      DEPSI1=DEPSI
      SI1=SI
      STEST=TEST*SN
      Z=0.
      DO 52 IB=NK,IBOU
 
      Y=UNSSM*E(IB)
      X=Y*DEPS
 
      IF(X.EQ.0.) GO TO 51
      EXPMX=EXP(-X)
      UNSX=1./X
      UNSX2=UNSX*UNSX
      X2=X*X
      IF (ABS(X).LT.1.D-17) GO TO 50
      SIGT(IB)=(SIGEL(IB)-DSIGP(IB)*UNSX)*EXPMX+DSIGP(IB)*UNSX
      Z=Z+COEF(JKI+IB)*SIGT(IB)*(EXPMX*(DSIGP(IB)*UNSX2-SIGEL(IB)+
     1DSIGP(IB)*UNSX)-DSIGP(IB)*UNSX2)*E(IB)
      GO TO 52
   50 CONTINUE
      SIGT(IB)=SIGEL(IB)*(1.-X+0.5*X2)+DSIGP(IB)*(1-0.5*X)
      Z=Z+COEF(JKI+IB)*SIGT(IB)*E(IB)*(DSIGP(IB)*((X/3.)-0.5)-
     1 SIGEL(IB)*(1-X+0.5*X2))
      GO TO 52
C         UNE VALEUR PROPRE EST NULLE
   51 SIGT(IB)=SIGEL(IB)+DSIGP(IB)
 52   CONTINUE
C NOUVEL DEPS
      SI=VNMISD(SIGT,ITYP,ALFAH,CVNMSD)
      X=PENTE-Z*UNSSM/SI
      ZOB1 = UNSSM/SI
      ZOB2 = ZOB1 * Z
      IF(JEBOUC.GT.1) WRITE(6,77888) PENTE,Z,UNSSM
77888 FORMAT('0 PENTE =',1PE16.5,2X,'Z=',1PE16.5,2X,'UNSSM=',1PE16.5/)
      DEPSI=(SI-SN)/X
      DEPS=DEPS+DEPSI
      EPST=DEPS+EPSTAR
      IF(JEBOUC.GT.1) WRITE(6,77884) SI,SN,X,DEPSI,DEPS
77884 FORMAT('0 SI =',1PE16.5,2X,'SN=',1PE16.5,2X, 
     .'X=',1PE16.5,2X, 'DEPSI=',1PE16.5,2X,'DEPS=',1PE16.5/)
      IF(ICINE.EQ.1) GO TO 580
C            ON CALCULE LA CONTRAINTE SUR LA COURBE DE TRACTION A
C            LA FIN DU PAS
      EPSTP=EPST
C           ON CONTROLE SI LE DEPS EST POSITIF
C           SI DEPS EST NEGATIF ON PREND EPSTP =EPST A L ITERATION PRECE
      IF(DEPS.LT.0.) EPSTP=EPSTP-DEPSI
      CALL TRACTI(SN,EPSTP,SIG,EPS,NCOURB,2,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GO TO 57
      ENDIF
C        ON CALCULE LA PENTE A LA COURBE DE TRACTION A LA FIN DU PAS
      CALL TRACTI(PENTE1,EPSTP,SIG,EPS,NCOURB,1,IBI)
      IF(ICINE.EQ.2) SN=SN-EPST*PANTIN
  580 CONTINUE
      IF(ITER.GE.ITMAX) GO TO 56
C            ON N A PAS DEPASSE LE NOMBRE MAX D ITERATIONS INTERNES
      IF(DEPS) 54,500,500
C             DEPS EST POSITIF
  500 CONTINUE
      PENTE=PENTE1
C           ON PREND COMME NOUVELLE PENTE LA PENTE A LA FIN DU PAS
      IF(ABS(SI-SN)-STEST) 57,57,555
C               A T ON CONVERGE
   56 SSTST=ABS(SI-SN)/SN
      KERRE=2
        write(6,fmt='('' kerre projvm'',i6)')kerre
      GO TO 57
   54 CONTINUE
C                 DEPS EST NEGATIF
C         LA PENTE A LA FIN DU PAS EST-ELLE SUPERIEURE A LA PENTE
C          AU DEBUT DANS LES CAS ISOTROPE ?
      IF(ICINE.EQ.1) GO TO 543
C
C   MILL 26/3/93
      IF(PENTE.EQ.0.D0.AND.DEPS.LE.0.D0) GOTO 540
C
      X=ABS(PENTE)
      IF(X.LT.ABS(PENTE1)) X=ABS(PENTE1)
      X=ABS((PENTE-PENTE1)/(X+1.D-10))
      IF(X.LT..2D00) GO TO 543
C                LA PENTE VARIE DE PLUS DE 20%
C               ON REGARDE SI LA COURBE EST CONCAVE
      IF(PENTE1.LT.PENTE) GO TO 543
      EPST=EPST-DEPS
      DEPS=DEPS-DEPSI
C               ON PREND  PENTE =.5*(PENTE AU DEBUT + PENTE A LA FIN )
      PENTE=0.5D00*(PENTE+PENTE1)
C              ON CALCULE  LE NOUVEAU DEPSI
      X=PENTE-Z*UNSSM/SI
      DEPSI=(SI-SN)/X
 
      DEPS=DEPS+DEPSI
      IF(DEPS.GT.0.D00) GO TO 542
C           DEPS EST ENCORE NEGATIF
      GO TO 540
  543 CONTINUE
C        LE POINT SI EST AU DESSUS DE LA SURFACE
C        ON RECOMMENCE AVEC LE PREMIER DEPS DIVISE PAR 2
C        LE POINT SI EST EN DESSOUS ON ASSIMILE LA TGT A LA CORDE
      EPST=EPST-DEPS
      DEPS=DEPS-DEPSI
      IF(SI.GT.SN) GO TO 540
      X=PENTE-(SI-SI1)/DEPSI1
      GO TO 541
540   CONTINUE
      IF(JYCOMP.GT.0) DEPS00=0.5*DEPS00
      JYCOMP=JYCOMP+1
      IF(JEBOUC.GT.1) WRITE(6,77887) JYCOMP
77887 FORMAT('0 JYCOMP = ',I4/)
      DEPS=0.5*DEPS00
      EPSOM=EPSTAR+DEPS
      DEPSI=DEPS
      CALL TRACTI(SN,EPSOM,SIG,EPS,NCOURB,2,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GO TO 57
      ENDIF
      SM=0.5D00*(SN+SELAS)
      UNSSM=1.D00/SM
      GO TO 542
  541 CONTINUE
      DEPSI=(SI-SN)/X
      DEPS=DEPS+DEPSI
      IF(DEPS.LE.0.D00) GO TO 540
  542 CONTINUE
      EPST=EPSTAR+DEPS
      GO TO 555
   57 CONTINUE
      IF( DEPS.LT.0.D0) THEN
        IF(JEBOUC.EQ.1) GO TO 2020
      ENDIF
      RETURN
      END