Sources Esope | Cast3M

endom1
C ENDOM1    SOURCE    CB215821  16/04/21    21:16:37     8920
      SUBROUTINE ENDOM1(SIGG,NSTRS,VARR,XMAT,NCOMAT,XCAR,ICARA,
     1NCOURB,TRAC,NNVARI,MFR,IFOURB,DD,DDINV,BB,R,S,XMU,
     2QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,YUNGV,XNUV)
C
C
C ===================================================================
C CE SOUS-PROGRAMME EST APPELE DANS "LOIS".
C IL CALCULE LES MATR. DD ET BB, LES VECT. R,S,XMU ET QSI DS.
C LE CAS DE L'ECROUISSAGE ET DE L'ENDOMMAGEMENT ISOTROPES.
C EN REGIME ANISOTHERME.
C
C ENTREES:
C -------
C NSTRS        = NBR. DE COMPOSANTES DES CONTR. OU DES DEFORM.
C SIGG(NSTRS)  = CONTR.
C NNVARI       =  NBR. DE VARIABLES INTERNES PILOTANT LES EQUATIONS
C                DU MODELE (VAUT 2)
C VARR(2)...   = VARIABLES INTERNES PILOTANT LES EQUATIONS DU MODELE
C ...VARR(1+NNVARI) (r ET D)
C NCOMAT       = NBR. DE CARACTERISTIQUES MECANIQUES DU MATERIAU
C XMAT(NCOMAT) = CARACTERISTIQUES MECANIQUES DU MATERIAU
C MFR          = INDICE DE LA FORMULATION MECANIQUE; SEULEMENT
C                MASSIF OU COQUE POUR LES MATERIAUX ENDOMMAGEABLES
C ICARA        = NBR. DE CARACT. GEOMETRIQUES DES ELEMENTS FINIS
C XCAR(ICARA)  = CARACT. GEOMETRIQUES DES ELEMENTS FINIS
C NCOURB       = NBR. DE PTS. DE LA COURBE DE TRACTION
C TRAC(2NCOURB)= COURBE DE TRACTION
C YUNGV        = DERIVEE DU MODULE D'YOUNG / TEMPERATURE
C XNUV         = DERIVEE DU COEFFICIENT DE POISSON / TEMPERATURE
C
C SORTIES:
C -------
C DD(NSTRS,NSTRS) = MATRICE DE HOOK ENDOMMAGEE
C BB(NSTRS,NNVARI)= MATRICE DEFINIE DS. LA RELATION SUIVANTE:
C                    .          .         .
C                   SIG = DD*EPSELAS + BB*Q   AVEC:
C                   - SIG     = CONTR.
C                   - EPSELAS = DEFORM. ELASTIQUES
C                   - Q       = VAR. INT. PILOTANT LES LOIS DU MODELE
C R(NSTRS)        = LOI D'ECOULEMENT PLASTIQUE
C S(NNVARI)       = LOIS D'EVOLUTION DES VARIABLES Q
C XMU(NSTRS)      = DERIVEE DU CRITERE PAR RAPPORT AUX CONTR.
C QSI(NNVARI)     = DERIVEE DU CRITERE PAR RAPPORT AUX VARIABLES Q
C KERRE           = INDICE QUI REGIT LES ERREURS
C                 = 99 SI LA FORMULATION MECANIQUE N'EST PAS
C                      DISPONIBLE POUR LE MODELE CONSIDERE OU S'IL Y
C                      A INCOMPATIBILITE ENTRE MFR ET IFOUR
C ====================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
C
      DIMENSION SIGG(*),VARR(*),XMAT(*),XCAR(*),TRAC(*),R(*),S(*),
     1XMU(*),QSI(*),SIGGD(*),W(*)
      DIMENSION DD(NSTRS,NSTRS),BB(NSTRS,NNVARI),F(NCOURB,2)
      DIMENSION DDINV(NSTRS,NSTRS)
C
C =================================
C REMPLISSAGE DES MATRICES DD(HOOK)
C =================================
C
      CALL ELAST1(1,IFOURB,VARR,NNVARI,XMAT,NCOMAT,YUNGV,XNUV,
     1          XCAR,ICARA,MFR,NSTRS,DD,DDINV,KERRE,INDIC,ITHHER)
C
C
C (SI LA FORMULATION MECANIQUE EST NON DISPONIBLE OU S'IL Y A
C  INCOMPATIBILITE ENTRE IFOUR ET MFR, KERRE EST NON NUL:)
      IF (KERRE.NE.0) GOTO 100
C
C =======================
C CALCUL DE LA MATRICE BB
C =======================
      CALL ZERO(BB,NSTRS,NNVARI)
      XD=VARR(3)
      DO 10 I=1,NSTRS
   10                BB(I,2)=-SIGG(I)/(1.D0-XD)
C
C =========================
C CALCUL DES VECT. R ET XMU
C =========================
C
      XNU=XMAT(2)
      EPSD=XMAT(6)
      DC=XMAT(7)
      EPSR=XMAT(8)
      P=VARR(1)
C --------------------
C DEVIATEUR DES CONTR.
C --------------------
      TRACE=SIGG(1)+SIGG(2)+SIGG(3)
      DO 30 I=1,NSTRS
                      A=0.D0
                      IF (I.LE.3) A=1.D0/3.D0
                      SIGGD(I)=SIGG(I)-A*TRACE
   30 CONTINUE
C
C -----------------------------------
C CONTRAINTE EQUIVALENTE DE VON MISES
C -----------------------------------
      SIGEQI=PROCON(SIGGD,SIGGD,NSTRS)
      SIGEQI=SQRT((3.D0/2.D0)*SIGEQI)
      DO 50 I=1,NSTRS
                A=1.D0
                IF (I.GT.3) A=2.D0
                R(I)=((3.D0*SIGGD(I))/(2.D0*(1.D0-XD)*SIGEQI))*A
                XMU(I)=R(I)
   50 CONTINUE
C
C =================
C CALCUL DU VECT. S
C =================
C
      S(1)=1.D0
C
      TRIAXI=TRACE/(3.D0*SIGEQI)
      TRIAXI=TRIAXI*TRIAXI
      TRIAXI=3.D0*(1.D0-2.D0*XNU)*TRIAXI+(2.D0*(1.D0+XNU))/3.D0
C
C --------------------------------------
C CALCUL DE LA CONSTANTE MATERIAU S0PRIM
C --------------------------------------
C
C ..........................................................
C CALCUL DES COEFF. K ET M A PARTIR DE LA COURBE DE TRACTION
C XKM ATTEND LA COURBE DE TRACTION (SIGMA,EPSILON_PLASTIQUE)
C OR TRAC CONTIENT (EPSILON_PLASTIQUE,SIGMA)
C ..........................................................
      DO 70 I=1,NCOURB
                       XTRAC=TRAC(2*I-1)
                       TRAC(2*I-1)=TRAC(2*I)
                       TRAC(2*I)=XTRAC
   70 CONTINUE
      CALL XKM(NCOURB,TRAC,XK,XM,W,F)
      DO 80 I=1,NCOURB
                       XTRAC=TRAC(2*I-1)
                       TRAC(2*I-1)=TRAC(2*I)
                       TRAC(2*I)=XTRAC
   80 CONTINUE
      SIGY=TRAC(2)
      XM2=(2.D0+XM)/XM
      XM1=(1.D0+XM)/XM
      S0PRI1=(EPSR**XM2)-(EPSD**XM2)
      S0PRI1=(XK*XK*S0PRI1)/XM2
      S0PRI2=(EPSR**XM1)-(EPSD**XM1)
      S0PRI2=(2.D0*SIGY*XK*S0PRI2)/XM1
      S0PRI3=SIGY*SIGY*(EPSR-EPSD)
      S0PRIM=(S0PRI1+S0PRI2+S0PRI3)/DC
C
C -----------------------------------------------------------
C "DERTRA" CALCULE LA DERIVEE DE LA COURBE DE TRACTION EN VIR
C -----------------------------------------------------------
      VIR=VARR(2)
      N2COUR=2*NCOURB
      CALL DERTRA(N2COUR,TRAC,VIR,Y,RPRIM,RINF,RSUP)
C
C
      S(2)=((Y**2.D0)*TRIAXI)/(S0PRIM*(1.D0-XD))
C
C ---------------------------------------------
C SI LE SEUIL D'ENDOMMAGEMENT N'EST PAS ATTEINT
C ---------------------------------------------
      IF (INDIC.EQ.0) S(2)=0.D0
C
C ===================
C CALCUL DU VECT. QSI
C ===================
C
      QSI(1)=-RPRIM
      QSI(2)=Y/(1.D0-XD)
C
  100 CONTINUE
      RETURN
      END