Télécharger ceraca.eso

Retour à la liste

ceraca
C CERACA    SOURCE    BP208322  17/03/01    21:15:27     9325           
      SUBROUTINE CERACA(WRK0,WRK1,WRK5,WTRAV,INPLAS,MFR,DT,NSTRS,NVARI,
     1                  PRECIS,MSOUPA,JECHER,DTT,NSSINC,INV,KERRE,
     2                  ICARA,IFOURB,CMATE,N2EL,N2PTEL,IB,IGAU,EPAIST,
     4                  NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,
     5                  CRIGI)
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
C---------------------------------------------------------------------
C Entree: INPLAS type de materiau
C         MFR indice de la formulation mecanique(seulement massif ou coque
C           pour les materiaux endommageables)
C         DEPST(NSTRS) increment des deformations totales
C         SIG0(NSTRS) contraintes initiales
C         EPIN0(NSTRS) deformations viscoplastiques initiales
C         VAR0(NVARI) variables internes initiales
C         NVARI nombre de variables internes
C         XMAT(NCOMAT) materiau
C         XCAR(ICARA) caracteristiques geometriques
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         NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
C         NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
C         SECT  = SECTION
C         LHOOK = TAILLE DE LA MATRICE DE HOOKE
C         DT pas de temps
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               <> 0 si entrees incoherentes
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER(I-N)
       IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
*
      SEGMENT WRK0
       REAL*8 XMAT(NCXMAT)
      ENDSEGMENT
*
      SEGMENT WRK1
        REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS)
        REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI)
        REAL*8 DEFP(NSTRS),XCAR(ICARA)
      ENDSEGMENT
*
      SEGMENT WRK5
        REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS)
      ENDSEGMENT
*
      SEGMENT WTRAV
        REAL*8 DDAUX(LHOOK,LHOOK),VALMAT(NUMAT)
        REAL*8 VALCAR(NUCAR),DSIGT(NSTRS)
        REAL*8 TXR(IDIM,IDIM),DDHOMU(LHOOK,LHOOK)
        REAL*8 XLOC(3,3),XGLOB(3,3)
        REAL*8 D1HOOK(LHOOK,LHOOK),ROTHOO(LHOOK,LHOOK)
        ENDSEGMENT
 
*
       DIMENSION SIG(6),EPSV(6),VAR(100)
       DIMENSION DSPT(6),XX(8),SIG1(6),EPSV1(6)
       DIMENSION VAR1(100),EVP1(6),VARP1(100)
       DIMENSION EVP2(6),VARP2(100)
       DIMENSION CRIGI(12)
*
       CHARACTER*8 CMATE
C**** debut ajout Eloi
       logical dtlibr
       dtlibr=.TRUE.
C**** fin ajout Eloi
 
 
C        ==================================================
C     RECHERCHE DU NUMERO DE LA VARIABLE INTERNE EGALE A
C     1 SI ON A ENDOMAGEMENT GENERALISE
           IBID2 = NVARI -1
C*************************************
C    On cherche le numero de la propriété du matériau correspondant à ENDG
       IF (MFR.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
C        fluage
        IF (VAR0(IBID2).EQ.1.D0) THEN
            DO 1032 IJ=1,NVARI
            VARF(IJ) = VAR0(IJ)
1032     CONTINUE
         DO 1022 K=1,NSTRS
            EPINF(K)=EPIN0(K)
           SIGF(K) = 0.D0
 1022  CONTINUE
           GO TO 998
            ENDIF
         DO 1010 I =1,NSTRS
          TOTO = ABS(EPIN0(I))
          IF(TOTO.GE.XMAT(IBID1)) THEN
          DO 1021 K=1,NSTRS
            EPINF(K)=EPIN0(K)
          SIGF(K) = 0.D0
 1021  CONTINUE
            DO 1030 IJ=1,(NVARI-1)
            VARF(IJ) = VAR0(IJ)
1030     CONTINUE
          VARF(IBID2) = 1.D0
            GO TO 998
          ENDIF
1010   CONTINUE
C*************************************
 
 
C =========================================================
       KERRE = 0
       IF(MFR.NE.1.AND.MFR.NE.3.OR.IFOURB.EQ.1)THEN
           KERRE = 99
           RETURN
       ENDIF
 
       IF(MFR.EQ.3) THEN
       THICK = XCAR(1)
       ALFA  = XCAR(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,NSTRS,CMATE,VALMAT,
     1   VALCAR,N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
     2   NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,XLOC,
     3           XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
C         PRINT *,'DSIGT =',DSIGT
*
       IF(IRTD.NE.1) THEN
          KERRE=69
          GOTO 998
       ENDIF
*
       IF(MFR.EQ.3) THEN
           DO 10 I=1,NSTRS/2
             SIG0(        I) = SIG0(        I)/THICK
             SIG0(NSTRS/2+I) = SIG0(NSTRS/2+I)*6.0D0/THICK/THICK
              DSIGT(       I) = DSIGT(        I)/THICK
              DSIGT(NSTRS/2+I)= DSIGT(NSTRS/2+I)*6.0D0/THICK/THICK
10         CONTINUE
           IF(IFOURB.EQ.-2) THEN
               SIG0 (2) =0.0
               SIG0 (4) =0.0
               DSIGT(2) =0.0
               DSIGT(4) =0.0
           ENDIF
       ENDIF
C
C------------------------------------------
C CONTROLE DE LA COHERENCE DES ENTREES
C Supprimée dans ce cas
C------------------------------------------
       IF (DT.LE.0.0) 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(MFR.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
*
30     CONTINUE
C
C-----------------------------
       IF (KERRE.NE.0) THEN
         GOTO 999
       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      MFR=1  (MASSIF)
C      MFR=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 (MFR.EQ.1) THEN
           DO 11 I=1,NSTRS
             A=1.D0
             IF (I.GT.3) A=2.D0
             EPIN0(I)=EPIN0(I)/A
 11        CONTINUE
 
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 (MFR.EQ.3) THEN
           IF (IFOURB.EQ.2) THEN
             DO 13 I=1,NSTRS
               A=1.D0
               IF (I.EQ.3) A=2.D0
               IF (I.EQ.6) A=2.D0
               EPIN0(I)=EPIN0(I)/A
 13          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,NSTRS))
       ERRABS = PRECIS*ASIG
       IF (XMAX.GT.ASIG) ERRABS = PRECIS*XMAX
       DO 40 I=1,NSTRS
         SIG(I) = SIG0(I)
         EPSV(I) = EPIN0(I)
         DSPT(I) = DSIGT(I)/DT
40     CONTINUE
C
 
       DO 48 I=1,NVARI
48       VAR(I)=VAR0(I)
C
C ---------------------------------------------------------------------
       INV = 0
       NSSINC = 0
C ---------------------------------------------------------------------
C       DEBUT DES ITERATIONS EN SSINCREMENTS  /FIN SI DTLEFT = 0
C ---------------------------------------------------------------------
  70   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 999
       ENDIF
C
C---------------------------------------------------------------------
C      START OF CALCULATIONS
C_____________________________________________________________________
 
      IF(MFR.NE.3) THEN
 
         CALL INCR11(TAU,SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS,NVARI,MFR)
       ELSE IF(MFR.EQ.3) THEN
         CALL INCR33(TAU,SIG,EPSV,VAR,XMAT,EVP1,VARP1,ALFA,NSTRS,NVARI)
       ENDIF
C
       NITERA  = 0
 
C --------------------------------------------------------------------
C      DEBUT DES ITERATIONS SUR TAU OPTIMAL /FIN SI RA PETIT
C --------------------------------------------------------------------
  80   NITERA  = NITERA + 1
       IF(MFR.EQ.3) GO TO 150
 
           CALL ADVA11(TAU,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
     &                 XMAT,NSTRS,NVARI,MFR)
 
           CALL INCR11(TAU,SIG1,EPSV1,VAR1,EVP2,VARP2,XMAT,NSTRS,NVARI,
     &                 MFR)
           DO 110 I=1,NSTRS
110          EVP2(I) = 0.5*(EVP1(I)+EVP2(I))
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,NSTRS,NVARI,MFR)
 
       GO TO 250
C _____________________________________________________________________
150    CONTINUE
C----------------------------------------------------------------------
C  CALCULATIONS FOR GENERALISED STRESS/STRAIN FORMULATIONS
C----------------------------------------------------------------------
           CALL ADVA22(TAU,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
     &                 XMAT,NSTRS,NVARI)
           CALL INCR33(TAU,SIG1,EPSV1,VAR1,XMAT,EVP2,VARP2,ALFA,NSTRS,
     &                 NVARI)
           DO 160 I=1,NSTRS,1
160          EVP2(I) = 0.5*(EVP1(I)+EVP2(I))
CCC Eloi : on peut effectuer la boucle suivante sur
CCC        les variables internes propres a Ottosen
CCC        qui n'ont pas ete modifiees
              DO 170 I=1,NVARI
                VARP2(I)= 0.5*(VARP1(I) + VARP2(I))
 170          CONTINUE
                   CALL ADVA22(TAU,SIG,EPSV,VAR,SIGF,EPINF,VARF,DSPT,
     &                         EVP2,VARP2,XMAT,NSTRS,NVARI)
250    CONTINUE
 
C ---------------------------------------------------------------------
C   CALCUL DU RAPPORT : ERREUR CALCULEE / ERREUR ADMISE
C ---------------------------------------------------------------------
       DO 260 I=1,NSTRS
260      XX(I) = SIGF(I)-SIG1(I)
       RA = SQRT(PROCON(XX,XX,NSTRS))/(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<RA<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 80
         ELSEIF ( RA.GT.BORNE) THEN
           TAU = TAU/SQRA
           GOTO 80
       ENDIF
C**** debut ajout Eloi
      endif
C**** fin ajout Eloi
C ---------------------------------------------------------------------
C       ici  ra < borne      cas JECHER  :
C ---------------------------------------------------------------------
       IF (JECHER.EQ.1) THEN
           DTT = TAU
           NSSINC = NITERA
           IF ((NSSINC.EQ.1).AND.(RA.EQ.0.0)) GOTO 999
           IF (NITERA.GE.8) GOTO 999
           IF (FAC*SQRA.LT.1.0) THEN
               TAU = TAU*FAC
               GOTO 80
             ELSEIF ((SQRA.LT.RMIN).OR.(SQRA.GT.RMAX)) THEN
               TAU = TAU/SQRA
               GOTO 80
           ENDIF
C ---------------------------------------------------------------------
C         ici  rmin < sqra < rmax    et  nitera < 8
C         pas de mise @ jour des variables
C ---------------------------------------------------------------------
          GOTO 999
       ENDIF
C ----------------------------------------------------------------------
C  FIN D'ITERATIONS  / MISE A JOUR DES VARIABLES
C  ici RA < 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 270 I=1,NSTRS,1
         SIG(I)  = SIGF(I)
         EPSV(I) = EPINF(I)
270    CONTINUE
       DO 280 I=1,NVARI,1
         VAR(I)  = VARF(I)
280    CONTINUE
C*************************************
C      Test pour savoir si on a dépassé la limite de déformation totale
 
         DO 2010 I =1,NSTRS
         TOTO = ABS(EPINF(I))
          IF(TOTO.GE.XMAT(IBID1)) THEN
          VARF(IBID2) = 1.D0
          DO 2021 KI=1,NSTRS,1
          SIGF(KI) = 0.D0
 2021  CONTINUE
            GO TO 998
            ELSE
             VARF(IBID2) = 0.D0
         ENDIF
2010   CONTINUE
C*************************************
C
C --------------------------------------------------------------------
C  TEST DE FIN SS INCREMENTS /  MISE A JOUR DE TAU
C  si   SQRA<1/3        TAU = TAU*3
C  si   1/3<SQRA<RMIN   on vise RA = 1
C  si   RMIN<SQRA<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 70
      ENDIF
C
      IF (ABS(TAU-DTLEFT).GT.(TAU/1000.)) THEN
         WRITE ( IOIMP,* ) ' PROBLEME   TAU > DTLEFT '
         KERRE = 223
      ENDIF
C-----------------------------------------------------------------------
999   CONTINUE
      IF(MFR.EQ.3) THEN
         DO 1000 I=1,NSTRS/2
           SIGF(        I) =SIGF(         I)*THICK
           SIGF(NSTRS/2+I) =SIGF(NSTRS/2+ I)*THICK*THICK/6.0
*           DSIGT(        I)=DSIGT(        I)*THICK
*           DSIGT(NSTRS/2+I)=DSIGT(NSTRS/2+I)*THICK*THICK/6.0
1000     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      MFR=1  (MASSIF)
C      MFR=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 (MFR.EQ.1) THEN
           DO 14 I=1,NSTRS
             A=1.D0
             IF (I.GT.3) A=2.D0
             EPIN0(I)=EPIN0(I)*A
             EPINF(I)=EPINF(I)*A
 14        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 (MFR.EQ.3) THEN
           IF (IFOURB.EQ.2) THEN
             DO 16 I=1,NSTRS
               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
 16          CONTINUE
           ENDIF
         ENDIF
C
C===========================================================
C
 998  RETURN
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales