C CERACA SOURCE CHAT 05/01/12 21:52:12 5004 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 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(IDIM,IDIM),XGLOB(IDIM,IDIM) 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 = DABS(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 2RMAX 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