Sources Esope | Cast3M

ccerac
C CCERAC    SOURCE    PV        17/12/08    21:15:14     9660           
      SUBROUTINE CCERAC(wrk52,wrk53,wrk54,NVARI,
     1           NSSINC,INV,IFOURB,IB,IGAU,NBPGAU,iecou,xecou)
C CERACA    SOURCE    AM        00/12/13    21:15:32     4045
c      SUBROUTINE CERACA(WRK0,WR00,WRK1,WRK5,WRK7,WRK8,WRK9,WTRAV,
c     1 INPLAS,MFR,DT,NSTRS,NVARI,NCOMAT,PRECIS,MSOUPA,JECHER,DTT,
c     2 NSSINC,INV,KERRE,ICARA,IFOURB,NYOG,NYNU,NYALFA,NYSMAX,
c     3 NYN,NYM,NYKK,NYALF1,NYBET1,NYR,NYA,NYKX,NNKX,NYRHO,NSIGY,T0,TF,
c     5 TREF,TLIFE,ITHHER,NCOURB,CMATE,N2EL,N2PTEL,IB,IGAU,EPAIST,
c     7 NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,CRIGI,
c     8 KERREU1)
C
C---------------------------------------------------------------------
C Objet: Calculer au cours d'un pas de temps DT, l'evolution des
C        variables internes a l'aide d'un schema Runge-Kutta 1.2
C ---------------------------------------------------------------------
C MFR1 &lt;- MFR, XCARB &lt;- XCAR, NSTRS1 &lt;- NSTRS,
C
C---------------------------------------------------------------------
C Entree: INPLAS type de materiau
C         MFR1 indice de la formulation mecanique(seulement massif ou coque
C           pour les materiaux endommageables)
C         DEPST(NSTRS1) increment des deformations totales
C         SIG0(NSTRS1) contraintes initiales
C         EPIN0(NSTRS1) deformations viscoplastiques initiales
C         VAR0(NVARI) variables internes initiales
C         NVARI nombre de variables internes
C         SIGY(NSIGY) courbe de la limite elastique en fonction de T°C
C         XMAT(NCOMAT) materiau
C         XCARB(ICARA) caracteristiques geometriques
C         YSMAX(NYSMAX) intervient ds. le test de convergence des iter.
C         PRECIS precision relative sur SIGMA
C         MSOUPA nombre maximal de sous pas autorises
C         JECHER = 0 avancer
C                = 1 rechercher sortie avec DTT
C         IFOURB = -2 EN CONTR.PLANES
C                  -1 EN DEFORM. PLANES
C                   0 EN AXISYMETRIE
C                   1 EN SERIE DE FOURIER
C                   2 EN TRIDIM
C         CMATE = NOM DU MATERIAU
C         VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
C         VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
C         N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
C         N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
C         IB    = NUMERO DE L ELEMENT COURANT
C         IGAU  = NUMERO DU POINT COURANT
C         EPAIST= EPAISSEUR
C         NBPGAU= NBRE DE POINTS DE GAUSS
C         MELE  = NUMERO DE L ELEMENT FINI
C         NPINT = NBRE DE POINTS D INTEGRATION
C         SECT  = SECTION
C         LHOOK = TAILLE DE LA MATRICE DE HOOKE
C         DT pas de temps
C         ITHHER = 0 pas de chargement thermique et materiau constant
C                = 1 chargement thermique et materiau constant
C                = 2 chargement thermique et materiau(T)
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
C Sortie: SIGF(NSTRS) contraintes finales
C         EPINF(NSTRS) deformations viscoplastiques finales
C         VARF(NVARI) variables internes finales
C         DTT sous-increment de temps optimal (si JECHER=1)
C         NSSINC nombre de sous-increments si JECHER=0
C         INV = 1 si inversion
C               0 sinon
C         KERRE = 0 si tout OK
C               &lt;> 0 si entrees incoherentes
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER(I-N)
       IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC DECHE
 
      SEGMENT IECOU
*      COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
       INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
C      INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM,  NYKK,
     1  icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
C    .  NYALF1,NYBET1,NYR,  NYA,   NYRHO,NSIGY,  NNKX,  NYKX,  IND,
     2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
C    .  NSOM,  NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT,  NCXMAT,
     3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
C    .  LTRAC, MFR,  IELE,  NHRM,   NBNN, NBELEM,ICARA,
     4  icow32,icow33,NSTRS1,MFR1,NBGMAT,NELMAT,MSOUPA,
C    .  LW2,   NDEF,  NSTRSS,MFR1,  NBGMAT,NELMAT,MSOUPA,
     5  icow39,icow40,icow41,icow42,icow43,icow44
C    .  NUMAT1,LENDO, NBBB,  NNVARI,KERR1, MELEME
        INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
     .          icow51,icow52,icow53,icow54,icow55,icow56
     .          icow57,icow58
      ENDSEGMENT
 
      SEGMENT XECOU
*      COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00
      REAL*8 xcow1, xcow2,xcow3,DTT,DT,xcow6, xcow7
C     REAL*8 DTOPTI,TSOM, TCAR, DTT,  DT,   TREFA,TEMP00
      ENDSEGMENT
 
 
 
*
*
       DIMENSION SIG(8),EPSV(8),VAR(100)
       DIMENSION DSPT(8),XX(8),SIG1(8),EPSV1(8)
       DIMENSION VAR1(100),EVP1(8),VARP1(100)
       DIMENSION EVP2(8),VARP2(100)
       DIMENSION CRIGI(12)
C**** debut ajout Eloi
       logical dtlibr
C**** fin ajout Eloi
        NCOMAT = NMATT
C**** debut ajout Eloi
        dtlibr=.TRUE.
C**** fin ajout Eloi
C        ==================================================
C RECHERCHE DU NUMERO DE LA VARIABLE INTERNE EGALE A 1
C SI ON A ENDOMMAGEMENT GENERALISE
C
      IBID2 = NVARI-1
C*************************************
C    On cherche le numero de la propriété du matériau correspondant à ENDG
C
      IF (MFR1.EQ.1.AND.IFOMOD.EQ.2) THEN
        IBID1 = 20
      ELSE
        IBID1 = 15
      ENDIF
C Test pour savoir si on a dépassé la limite de déformation de fluage
C
      IF (VAR0(IBID2).EQ.1.D0) THEN
        DO 10 IJ=1,NVARI
          VARF(IJ) = VAR0(IJ)
  10    CONTINUE
        DO 20 K=1,NSTRS1
          EPINF(K) = EPIN0(K)
          SIGF(K) = 0.D0
  20    CONTINUE
        GOTO 1000
      ENDIF
      DO 30 I=1,NSTRS1
        TOTO = ABS(EPIN0(I))
        IF (TOTO.GE.XMAT(IBID1)) THEN
          DO 40 K=1,NSTRS1
            EPINF(K) = EPIN0(K)
            SIGF(K) = 0.D0
  40      CONTINUE
          DO 50 IJ=1,(NVARI-1)
            VARF(IJ) = VAR0(IJ)
  50      CONTINUE
          VARF(IBID2) = 1.D0
          GOTO 1000
        ENDIF
  30  CONTINUE
C*************************************
 
 
C =========================================================
      KERRE = 0
      IF (MFR1.NE.1.AND.MFR1.NE.3.OR.IFOURB.EQ.1) THEN
        KERRE = 99
        RETURN
      ENDIF
 
      IF (MFR1.EQ.3) THEN
        THICK = XCARB(1)
        ALFA  = XCARB(2)
      ENDIF
      BORNE = 2.0
      RMAX  = 1.3
      RMIN  = 0.7
      DIV   = 7.0
      FAC   = 3.0
C
C   CALCUL DES INCREMENTS DE CONTRAINTES ELASTIQUES
C
      CALL CALSIG(DEPST,DDAUX,NSTRS1,CMATE,VALMAT,VALCAR,N2EL,N2PTEL,
     1            MFR1,IFOURB,IB,IGAU,EPAIST,NBPGAU,MELE,NPINT,NBGMAT,
     2            NELMAT,SECT,LHOOK,TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,
     3            CRIGI,DSIGT,IRTD)
C         PRINT *,'DSIGT =',DSIGT
*
      IF (IRTD.NE.1) THEN
        KERRE=69
        GOTO 1000
      ENDIF
*
      IF (MFR1.EQ.3) THEN
        DO 60 I=1,NSTRS1/2
          SIG0(        I) = SIG0(        I)/THICK
          SIG0(NSTRS1/2+I) = SIG0(NSTRS1/2+I)*6.0D0/THICK/THICK
          DSIGT(       I) = DSIGT(        I)/THICK
          DSIGT(NSTRS1/2+I)= DSIGT(NSTRS1/2+I)*6.0D0/THICK/THICK
  60    CONTINUE
        IF (IFOURB.EQ.-2) THEN
          SIG0 (2) = 0.D0
          SIG0 (4) = 0.D0
          DSIGT(2) = 0.D0
          DSIGT(4) = 0.D0
        ENDIF
      ENDIF
C
C------------------------------------------
C CONTROLE DE LA COHERENCE DES ENTREES
C Supprimée dans ce cas
C------------------------------------------
      IF (DT.LE.0.D0) KERRE = 414
      MOTERR(1:8) = ' CONST  '
C      IF (INCOMAT.NE. 15)KERRE = 146
C       IF (IFOURB.EQ.-2.AND.NCOMAT.NE. 9)KERRE = 146
      IF (IFOURB.NE.2.AND.IFOURB.NE.-2.AND.IFOURB.NE.0
     &              .AND.IFOURB.NE.-1) THEN
        KERRE = 194
        MOTERR(1:8) = 'FLUAGE'
      ENDIF
      IF (MFR1.EQ.1.AND.IFOMOD.EQ.2) THEN
        IBID = 15
      ELSE
        IBID = 10
      ENDIF
      XMAX=XMAT(IBID)
*
*  TEST SUR XMAX          MILL  8/3/91
*
      IF (XMAX.EQ.0.D0) THEN
        XMAX=XMAT(1)*1.D-3
      ENDIF
C
C-----------------------------
      IF (KERRE.NE.0) THEN
        GOTO 1100
      ENDIF
C
C===========================================================
C  A PARTIR DE MAINTENANT, LES DEFORMATIONS
C  DE CISAILLEMENT NE SONT PLUS
C  DEFINIES PAR DES GAMA.
C  ON DIVISE DONC LES TERMES DE CISAILLEMENT PAR 2.
C  SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR CERACARO:
C      MFR1=1  (MASSIF)
C      MFR1=3  (COQUES MINCES)
C
C
C  Cas de la formulation massive
C  Les termes de cisaillement apparaissent
C  au delà de la troisieme composante
C
      IF (MFR1.EQ.1) THEN
        DO 70 I=1,NSTRS1
          A=1.D0
          IF (I.GT.3) A=2.D0
          EPIN0(I)=EPIN0(I)/A
  70    CONTINUE
C
C  Cas des coques minces
C  Les termes de cisaillement apparaissent
C  pour la troisieme et la sixieme composante
C  uniquement dans les cas de calculs
C  tridimensionnels ou d'analyse de Fourier
C
      ELSE IF (MFR1.EQ.3) THEN
        IF (IFOURB.EQ.2) THEN
          DO 80 I=1,NSTRS1
            A=1.D0
            IF (I.EQ.3) A=2.D0
            IF (I.EQ.6) A=2.D0
            EPIN0(I)=EPIN0(I)/A
  80      CONTINUE
        ENDIF
      ENDIF
C
C===========================================================
C
C ----------------
C  INITIALISATION
C ----------------
C**** debut ajout Eloi
      ITERO = 0
1200  CONTINUE
 
      itero = 1 + itero
      if (itero.ne.1) THEN
c            write(6,*) 'itero ib igau', itero,ib,igau
        dtlibr = .true.
        precis = precis * 7.
c       write(6,*)  ' precision modifiée ', precis
        if (itero.gt.3) then
****      kerre = 460
          kerre = 268
          return
        endif
      endif
C**** fin ajout Eloi
      DTLEFT = DT
      TAU    = DTLEFT
      ASIG = SQRT(PROCON(SIG0,SIG0,NSTRS1))
      ERRABS = PRECIS*ASIG
      IF (XMAX.GT.ASIG) ERRABS = PRECIS*XMAX
      DO 90 I=1,NSTRS1
        SIG(I) = SIG0(I)
        EPSV(I) = EPIN0(I)
        DSPT(I) = DSIGT(I)/DT
  90  CONTINUE
C
      DO 100 I=1,NVARI
        VAR(I)=VAR0(I)
 100  CONTINUE
C
C ---------------------------------------------------------------------
      INV = 0
      NSSINC = 0
C ---------------------------------------------------------------------
C       DEBUT DES ITERATIONS EN SSINCREMENTS  /FIN SI DTLEFT = 0
C ---------------------------------------------------------------------
1400  CONTINUE
      NSSINC = NSSINC + 1
      IF (NSSINC.GT.MSOUPA) THEN
C**** debut ajout Eloi
        DTLIBR=.FALSE.
        GOTO 1200
C**** fin ajout Eloi
C***       KERRE = 460
C***       GOTO 1100
      ENDIF
C
C---------------------------------------------------------------------
C      START OF CALCULATIONS
C_____________________________________________________________________
 
      IF(MFR1.NE.3) THEN
        CALL INCR11(TAU,SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS1,NVARI,MFR1)
      ELSE IF(MFR1.EQ.3) THEN
        CALL INCR33(TAU,SIG,EPSV,VAR,XMAT,EVP1,VARP1,ALFA,NSTRS1,NVARI)
      ENDIF
C
      NITERA  = 0
 
C --------------------------------------------------------------------
C      DEBUT DES ITERATIONS SUR TAU OPTIMAL /FIN SI RA PETIT
C --------------------------------------------------------------------
 
1300  CONTINUE
      NITERA  = NITERA + 1
      IF (MFR1.NE.3) THEN
 
        CALL ADVA11(TAU,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
     &              XMAT,NSTRS1,NVARI,MFR1)
        CALL INCR11(TAU,SIG1,EPSV1,VAR1,EVP2,VARP2,XMAT,NSTRS1,NVARI,
     &              MFR1)
        DO 110 I=1,NSTRS1
          EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
 110    CONTINUE
CCC Eloi : on peut effectuer la boucle suivante sur
CCC        les variables internes propres a Ottosen
CCC        qui n'ont pas ete modifiees
        DO 120 I=1,NVARI
          VARP2(I)= 0.5*(VARP1(I) + VARP2(I))
 120    CONTINUE
        CALL ADVA11(TAU,SIG,EPSV,VAR,SIGF,EPINF,VARF,DSPT,EVP2,VARP2,
     &              XMAT,NSTRS1,NVARI,MFR1)
      ELSE
C----------------------------------------------------------------------
C  CALCULATIONS FOR GENERALISED STRESS/STRAIN FORMULATIONS
C----------------------------------------------------------------------
        CALL ADVA22(TAU,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
     &              XMAT,NSTRS1,NVARI)
        CALL INCR33(TAU,SIG1,EPSV1,VAR1,XMAT,EVP2,VARP2,ALFA,NSTRS1,
     &              NVARI)
        DO 140 I=1,NSTRS1,1
          EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
 140    CONTINUE
CCC Eloi : on peut effectuer la boucle suivante sur
CCC        les variables internes propres a Ottosen
CCC        qui n'ont pas ete modifiees
        DO 150 I=1,NVARI
          VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
 150    CONTINUE
        CALL ADVA22(TAU,SIG,EPSV,VAR,SIGF,EPINF,VARF,DSPT,EVP2,VARP2,
     &              XMAT,NSTRS1,NVARI)
      ENDIF
 
C ---------------------------------------------------------------------
C   CALCUL DU RAPPORT : ERREUR CALCULEE / ERREUR ADMISE
C ---------------------------------------------------------------------
      DO 170 I=1,NSTRS1
        XX(I) = SIGF(I)-SIG1(I)
 170  CONTINUE
      RA = SQRT(PROCON(XX,XX,NSTRS1))/(ERRABS)
      SQRA = SQRT(RA)
 
C ---------------------------------------------------------------------
C  TEST DE FIN D'ITERATIONS /  MISE A JOUR DE TAU  /OPTION JECHER
C  DIV =7   BORNE = 2
C  SI   SQRA>7   TAU = TAU/7    ET NOUVEL ESSAI
C  SI 2&lt;RA&lt;7*7   ON VISE RA = 1  ET NOUVEL ESSAI
C ------------------------------------------------------------------
C**** debut ajout Eloi
      IF (dtlibr) Then
C**** fin ajout Eloi
        IF (RA.GT.DIV*DIV ) THEN
          TAU = TAU/DIV
          GOTO 1300
        ELSEIF (RA.GT.BORNE) THEN
          TAU = TAU/SQRA
          GOTO 1300
        ENDIF
C**** debut ajout Eloi
      endif
C**** fin ajout Eloi
 
C ---------------------------------------------------------------------
C       ici  ra &lt; borne      cas JECHER  :
C ---------------------------------------------------------------------
      IF (JECHER.EQ.1) THEN
        DTT = TAU
        NSSINC = NITERA
        IF ((NSSINC.EQ.1).AND.(RA.EQ.0.0)) GOTO 1100
        IF (NITERA.GE.8) GOTO 1100
        IF (FAC*SQRA.LT.1.0) THEN
          TAU = TAU*FAC
          GOTO 1300
        ELSEIF ((SQRA.LT.RMIN).OR.(SQRA.GT.RMAX)) THEN
          TAU = TAU/SQRA
          GOTO 1300
        ENDIF
C ---------------------------------------------------------------------
C         ici  rmin &lt; sqra &lt; rmax    et  nitera &lt; 8
C         pas de mise @ jour des variables
C ---------------------------------------------------------------------
        GOTO 1100
      ENDIF
C ----------------------------------------------------------------------
C  FIN D'ITERATIONS  / MISE A JOUR DES VARIABLES
C  ici RA &lt; BORNE
C  fin des boucles sur tau optimal
C  on avance en temps
C  mise @ jour de SIG etc...
C -------------------------------------------------------------------
C --
C -- ajout de ivtest = 1 par chloe
C --
      IVTEST = 1
      INV = INV + IVTEST
      DO 180 I=1,NSTRS1
        SIG(I)  = SIGF(I)
        EPSV(I) = EPINF(I)
 180  CONTINUE
      DO 190 I=1,NVARI
        VAR(I)  = VARF(I)
 190  CONTINUE
C*************************************
C      Test pour savoir si on a dépassé la limite de déformation totale
C
      DO 200 I =1,NSTRS1
        TOTO = ABS(EPINF(I))
        IF(TOTO.GE.XMAT(IBID1)) THEN
          VARF(IBID2) = 1.D0
          DO 210 KI=1,NSTRS1,1
            SIGF(KI) = 0.D0
 210      CONTINUE
          GOTO 1000
        ELSE
          VARF(IBID2) = 0.D0
        ENDIF
 200  CONTINUE
C*************************************
C
C --------------------------------------------------------------------
C  TEST DE FIN SS INCREMENTS /  MISE A JOUR DE TAU
C  si   SQRA&lt;1/3        TAU = TAU*3
C  si   1/3&lt;SQRA&lt;RMIN   on vise RA = 1
C  si   RMIN&lt;SQRA&lt;RMAX  TAU inchangé
C  si   SQRA>RMAX       on vise RA = 1
C    fin des boucles en ss increments  si tau = dtleft
C --------------------------------------------------------------------
C
      IF (TAU.LT.DTLEFT) THEN
        DTLEFT = DTLEFT - TAU
        IF (FAC*SQRA.LT.1.D0) THEN
          TAU=TAU*FAC
        ELSEIF ( (SQRA.LT.RMIN).OR.(SQRA.GT.RMAX) ) THEN
          TAU=TAU/SQRA
        ENDIF
        IF (TAU.GT.DTLEFT) TAU = DTLEFT
        GOTO 1400
      ENDIF
C
      IF (ABS(TAU-DTLEFT).GT.(TAU/1000.)) THEN
         WRITE ( IOIMP,* ) ' PROBLEME   TAU > DTLEFT '
         KERRE = 223
      ENDIF
C-----------------------------------------------------------------------
1100  CONTINUE
      IF (MFR1.EQ.3) THEN
        DO 220 I=1,NSTRS1/2
          SIGF(        I) =SIGF(         I)*THICK
          SIGF(NSTRS1/2+I) =SIGF(NSTRS1/2+ I)*THICK*THICK/6.0
*          DSIGT(        I)=DSIGT(        I)*THICK
*          DSIGT(NSTRS1/2+I)=DSIGT(NSTRS1/2+I)*THICK*THICK/6.0
 220    CONTINUE
      ENDIF
 
C
C===========================================================
C  RETOUR A LA DEFINITION NORMALE DES DEFORMATIONS
C  A SAVOIR: LES DEFORMATIONS DE CISAILLEMENT SONT
C  DEFINIES PAR DES GAMA.
C  SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR Ceracaro:
C      MFR1=1  (MASSIF)
C      MFR1=3  (COQUES MINCES)
C
C  Cas de la formulation massive
C  Les termes de cisaillement apparaissent
C  au delà de la troisieme composante
C
      IF (MFR1.EQ.1) THEN
        DO 230 I=1,NSTRS1
          A=1.D0
          IF (I.GT.3) A=2.D0
          EPIN0(I)=EPIN0(I)*A
          EPINF(I)=EPINF(I)*A
 230    CONTINUE
C
C  Cas des coques minces
C  Les termes de cisaillement apparaissent
C  pour la troisieme et la sixieme composante
C  uniquement dans les cas de calculs
C  tridimensionnels ou d'analyse de Fourier
C
      ELSE IF (MFR1.EQ.3) THEN
        IF (IFOURB.EQ.2) THEN
          DO 240 I=1,NSTRS1
            A=1.D0
            IF (I.EQ.3) A=2.D0
            IF (I.EQ.6) A=2.D0
            EPIN0(I)=EPIN0(I)*A
            EPINF(I)=EPINF(I)*A
 240       CONTINUE
        ENDIF
      ENDIF
C
C===========================================================
C
1000  RETURN
      END