Sources Esope | Cast3M

vparf
C VPARF     SOURCE    CHAT      05/01/13    04:09:39     5004
      SUBROUTINE VPARF(YOUNG,POISSON,Y,CK,CN,RSIG0,RDSIG,RDEFP,
     &                RSIGF,DT)
C
C ==================================================================
C CE SOUS-PROGRAMME REALISE L'INTEGRATION DE LA LOI DE VISCOPLASTICITE
C
C ENTREES:
C -------
C YOUNG
C POISSON
C Y  limite elastique
C CK coefficient de viscosite
C CN puissance de la loi d'ecoulement
C RSIG0(6)  = CONTR. AU DEBUT DU PAS D'INTEGRATION
C RDSIG(6) = INCREMENT DES CONTR. CALCULE ELASTIQUEMENT
C                A PARTIR DE L'INCREMENT DES DEFORM. TOTALES
C
C SORTIES:
C -------
C RSIGF(NSTRS)= CONTR. A LA FIN DU PAS D'INTEGRATION
C RDEFP(NSTRS)= INCREMENT DES DEFORM. PLASTIQUES A LA FIN DU PAS
C              D'INTEGRATION
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION RSIG0(6),RDSIG(6),RDEFP(6),RSIGF(6),RSIGTR(6)
      DIMENSION COEF(6)
      F(X,A,B,C,D,E,G) = X + 2.D0*G*B*((X-A)/C)**D - E
c      a c'est y
c      b c'est mu
c      c c'est k
c      d c'est cn
c      e c'est seqtr
c      g c'est dt
c
c  executable
*
      DATA COEF/1.D0,1.D0,1.D0,2.D0,2.D0,2.D0/
*
      DO 5 I=1,6
          RDEFP(I) = 0.D0
          RSIGF(I) = 0.D0
  5   CONTINUE
      CMU = YOUNG / 2.D0 / (1.d0 +POISSON)
      DO 10 I=1,6
        RSIGTR(I) = RSIG0(I) + RDSIG(I)
  10  CONTINUE
c
c   partie hydrostatique du trial stress
      PTR = (RSIGTR(1) +  RSIGTR(2) + RSIGTR(3))/3.D0
c
c   on calcul la partie deviatorique du trial stress
      DO 20 I=1,3
         RSIGTR(I) = RSIGTR(I) - PTR
  20  CONTINUE
c
c   contrainte equivalente du trial stress
      SEQTR = 0.D0
      DO 30 I=1,6
          SEQTR = SEQTR + RSIGTR(I)*RSIGTR(I)*COEF(I)
  30  CONTINUE
      SEQTR = SQRT( 1.5D0 * SEQTR)
c
      IF ( SEQTR .GT. Y ) THEN
c          le materiau subit une evolution nonlineaire
c
c       ordre de grandeur de l'equation pour la convergence
        ODG = ABS (Y-SEQTR)
c      resolution de l'equation f(x) = 0
        X1 = Y
        X2 = SEQTR
c        Y1 = F(X1,Y,CMU,CK,CN,SEQTR,DT)
c        Y2 = F(X2,Y,CMU,CK,CN,SEQTR,DT)
c
        DO 40 I=1,200
c          resolution par dicotomie
c         X3 = X1 + abs((X2 - X1)*Y1/ (Y2-Y1))
         X3 = (X1 +X2) /2.d0
         Y3 =  F(X3,Y,CMU,CK,CN,SEQTR,DT)
         IF ( Y3 .GT. (1.D-4*ODG) ) THEN
           X2 = X3
c           Y2 = Y3
         ELSE IF ( Y3 .LT. (-1.D-4*ODG) ) THEN
           X1 = X3
c           Y1 = Y3
         ELSE
c       on a converge
           GOTO 50
         ENDIF
  40    CONTINUE
        CALL ERREUR(268)
c        write (*,1000) young,poisson,y,ck,cn
c 1000   format('E=',1X,G12.7,1X,'Nu=',G12.7,1X,'Y=',G12.7/
c     & 'K=',G12.7,1X,'N=',G12.7)
c        WRITE(*,1001) X1,X2
c        WRITE(*,1001) Y1,Y2
c 1001   format('X1=',G12.7,' X2=',G12.7)
c        Print *, 'AIE CA NE CONVERGE PAS'
  50    CONTINUE
*
        AMPLI = X3 / SEQTR
        DO 60 I=1,6
          RSIGF(I) = RSIGTR(I) * AMPLI
  60    CONTINUE
        AMPLI = (1.D0/AMPLI -1.D0) / 2.D0 / CMU
        DO 70 I=1,6
          RDEFP(I) = RSIGF(I) * AMPLI * COEF(I)
  70    CONTINUE
c       de deviatorique a complet
        DO 80 I=1,3
          RSIGF(I) = RSIGF(I) + PTR
  80    CONTINUE
c        verification
c       zero = rsigtr(1) - rsigf(1)+ptr - 2.d0*cmu*rdefp(1)
c       if ( (abs(zero)) .GT. 1.d-3 ) then
c         print *, 'ouille',zero
c       endif
      ELSE
c    l'evolution est elastique
        DO 90 I=1,6
          RSIGF(I) = RSIGTR(I)
  90    CONTINUE
        DO 100 I=1,3
          RSIGF(I) = RSIGF(I) + PTR
  100    CONTINUE
      ENDIF
c
      RETURN
      END