C CCOIN0 SOURCE OF166741 25/11/04 21:15:09 12349 SUBROUTINE CCOIN0(wrk52,wrk53,wrk54,wrk2,wrk3, & IB,IGAU,NBPGAU,LTRAC2,IFOUR2,iecou) C----------------------------------------------------------------------- C ECOULEMENT PLASTIQUE POUR UN POINT C ALGORITHME ORTIZ ET SIMO C C EN ENTREE : C C SIG0 CONTRAINTES AU DEBUT DU PAS C DEPST INCREMENT DE DEFORMATIONS TOTALES C ( THERMIQUE ENLEVEE ) C VAR0 VARIABLES INTERNES DEDUT DU PAS C VAREX0 VARIABLES EXTERNES DEBUT DU PAS C VAREXF VARIABLES EXTERNES FIN DU PAS C XMAT COEFFICIENTS DU MATERIAU C PRECIS PRECISION POUR ITERATIONS INTERNES C WORK DES CARACTERISTIQUES C TRAC COURBE DE TRACTION C MFR1 INDICE DE FORMULATION (iecou) C NSTRSS NOMBRE DE CONTRAINTES (iecou) C INPLAS NUMERO DU MODELE DE PLASTICITE C DDAUX = MATRICE DE HOOKE ELASTIQUE 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 IFOUR2 = OPTION DE CALCUL (redondant avec IFOUR, iecou.IFOURB 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 (iecou) C NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES (iecou) C SECT = SECTION C LHOOK = TAILLE DE LA MATRICE DE HOOKE C TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES C UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE C C EN SORTIE : C C SIGF CONTRAINTES A LA FIN DU PAS C VARF VARIABLES INTERNES FIN DU PAS C DEFP INCR. DE DEFORMATIONS PLASTIQUES C KERRE CODE D'ERREUR C = 0 SI TOUT OK C = 99 SI FORMULATION NON DISPONIBLE C EN PLASTICITE C = 1 SI DEPS EST NEGATIF C = 2 SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE C C IFOUR : OPTION DE CALCUL C C IFOUR = -3 DEFORMATION PLANE GENERALISEE C IFOUR = -2 CONTRAINTES PLANES C IFOUR = -1 DEFORMATIONS PLANES C IFOUR = 0 AXISYMETRIQUE C IFOUR = 1 SERIE DE FOURIER C IFOUR = 2 TRIDIMENSIONNEL C C MFR1 : NUMERO DE LA FORMULATION ELEMENT FINI C C MFR1 = 1 MASSIF C MFR1 = 3 COQUE C MFR1 = 5 COQUE EPAISSE C MFR1 = 7 POUTRE C MFR1 = 9 COQUE AVEC CISAILLEMENT TRANSVERSE Ckich MFR1 = 31 pondération réduite termes diagonaux matrice B, C dite MASSIF INCOMPRESSIBLE. Utilisation en contraintes planes a justifier c doublon historique MFR/MFR1 C C INPLAS : NUMERO DU MODELE DE PLASTICITE C C INPLAS = 1 PARFAIT C INPLAS = 4 CINEMATIQUE C INPLAS = 5 ISOTROPE C INPLAS = 7 CHABOCHE1 C INPLAS = 12 CHABOCHE2 C C----------------------------------------------------------------------- C CONVENTION DE REMPLISSAGE DES MEMOIRES C----------------------------------------------------------------------- C C XMAT(1) MODULE D'YOUNG C XMAT(2) COEFFICIENT DE POISSON C C TRAC(1 A 2*LTRAC) COURBE DE TRACTION C WORK( " +1) ALFA1 POUR COQUES PLASTICITE GLOBALE C WORK( " +..) DONNEES POUR CRITERE POUTRES GLOBALES C C MODELE ISOTROPE C VAR0(1) EPS* C C MODELE CINEMATIQUE LINEAIRE C VAR0(1) EPS* C VAR0(2) A VAR0(1+NSTRSS) DEPLACEMENT DE LA SPHERE C C MODELE CHABOCHE C XMAT(5) .... COEFFICIENTS A,C,... C VAR0(1) EPS* C VAR0(2) A VAR0(1+NSTRSS) DEPLACEMENT DE LA SPHERE 1 C VAR0(2+NSTRSS) A VAR0(1+2*NSTRSS) " " " " 2 C C----------------------------------------------------------------------- C 20/09/2017 : modif SG critere de convergence trop serre C TEST=PETI*APHI0 + utilisation CCREEL C voir aussi ecoin0.eso, syco12.eso IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC DECHE -INC TECOU SEGMENT WRK2 REAL*8 TRAC(LTRAC2) ENDSEGMENT SEGMENT WRK3 REAL*8 WORK(LW),WORK2(LW2) ENDSEGMENT DIMENSION SIG(130),EPS(130) DIMENSION S(8),SX(8),DS(8),DSIG(8),SPHER(8),SPHER1(8),SPHER2(8) DIMENSION DSPHER1(8),DSPHER2(8),DEPSE(8),DEPSP(8),DDEPSE(8) DIMENSION F(8),WF1(8),WF2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8) DIMENSION CRIGI(12) jNPLAS = INPLAS MFRl = iecou.MFR1 nstrsl = iecou.NSTRSS ncara = commat(/2) if(ib.eq.1.and.igau.eq.1) then do iaca = 1,ncara if(commat(iaca).eq.'LIMP') iecou.JELEM = iaca enddo endif icow21 = iecou.JELEM if (icow21.gt.0) xlimp = valma0(icow21) do jj = 1,8 sx(jj) = 0.d0 enddo C---------COQUES AVEC CT------------------------------------------------ C ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES IF (MFRl.EQ.9) THEN NCONT=6 ELSE NCONT=iecou.NSTRSS ENDIF itracb=0 C----------------------------------------------------------------------- 2222 continue KERRE=0 DO I=1,nstrsl S(I)=0.D0 SPHER(I)=0.D0 SPHER1(I)=0.D0 SPHER2(I)=0.D0 ENDDO YUNG=XMAT(1) XNU=XMAT(2) C---------CARACTERISTIQUES GEOMETRIQUES--------------------------------- C C COQUES C IF (MFRl.EQ.3.OR.MFRl.EQ.9) THEN EP1 = WORK(1) ALFA1 = WORK(2)**2 ELSE IF (MFRl.EQ.5) THEN EP1 = WORK(1) ALFA1 = 1.D0 ELSE EP1 = 1.D0 ALFA1 = 1.D0 ENDIF C---------COEFFICIENTS CHABOCHE----------------------------------------- IF(jNPLAS.EQ.7) THEN A1=XMAT(5) C1=XMAT(6) R0=XMAT(7) PSI=XMAT(8) OME=XMAT(9) RM=XMAT(10) B=XMAT(11) A2=0.D0 C2=0.D0 ELSE IF(jNPLAS.EQ.12) THEN A1=XMAT(5) C1=XMAT(6) A2=XMAT(7) C2=XMAT(8) R0=XMAT(9) PSI=XMAT(10) OME=XMAT(11) RM=XMAT(12) B=XMAT(13) ELSE IF(jNPLAS.EQ.4) THEN SLI0 = XMAT(5) HECL = XMAT(6) c*dbg if (ib.eq.1 .and. igau.eq.1) then c*dbg write(ioimp,*) 'Ecrouissage',jNPLAS,SLI0,HECL c*dbg endif ELSE DO I=1,LTRAC2 SIG(I)=TRAC(2*I-1) EPS(I)=TRAC(2*I) ENDDO ENDIF EPST=VAR0(1) C---------ECROUISSAGE CINEMATIQUE--------------------------------------- IF(jNPLAS.EQ.4.OR.jNPLAS.EQ.7.OR.jNPLAS.EQ.12) THEN DO I=1,nstrsl SPHER1(I)=VAR0(I+1) ENDDO IF(jNPLAS.EQ.12) THEN DO I=1,nstrsl SPHER2(I)=VAR0(nstrsl+1+I) ENDDO ENDIF DO I=1,nstrsl SPHER(I)=SPHER1(I)+SPHER2(I) ENDDO c*dbg if (ib.eq.1 .and.igau.eq.1) then c*dbg write(ioimp,*) ' SPHERi',nstrsl c*dbg write(ioimp,*) ' ',(spher1(iof),iof=1,nstrsl) c*dbg write(ioimp,*) ' ',(spher(iof),iof=1,nstrsl) c*dbg endif C----------------------------------------------------------------------- C PREDICTEUR ELASTIQUE C S : CONTRAINTE C SPHER : VARIABLE D'ECROUISSAGE CINEMATIQUE C SX = S - X C----------------------------------------------------------------------- * en elastique non lineaire on annule les contraintes initiales * mais on cumule les epsilons elastiques ELSE IF(jNPLAS.EQ.87) THEN EPST=0.D0 DO I=1,nstrsl SIG0(I)=0.D0 DEPST(I)=DEPST(I) + VAR0(I+1) ENDDO ENDIF CALL CALSIG(DEPST,DDAUX,nstrsl,CMATE,VALMAT,VALCAR, & N2EL,N2PTEL,MFRl,IFOUR2,IB,IGAU,EPAIST, & NBPGAU,MELE,NPINT,iecou.NBGMAT,iecou.NELMAT, & SECT,LHOOK,TXR, & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD) IF(IRTD.NE.1) THEN KERRE=69 GOTO 1000 ENDIF IF ((mfr.eq.1.or.mfr.eq.31).and. IFOUR2.eq.-2) then * en cas de contraintes planes on repasse en 3D do ju=1,6 do iu=1,6 DDA(iu,ju)=0.d0 enddo enddo cte_cp = xnu / (1.d0 - xnu) aux= YUNG / ((1.d0+xnu)*(1.d0-2.d0*xnu)) aux1= aux * xnu aux2= aux * (1.d0-xnu) gege = Yung / (2.d0 *(1.d0 +xnu)) DDA(1,1)=AUX2 DDA(2,1)=AUX1 DDA(1,2)=AUX1 DDA(2,2)=AUX2 DDA(3,3)=aux2 DDA(1,3)=aux1 DDA(2,3)=aux1 DDA(3,1)=aux1 DDA(3,2)=aux1 DDA(4,4)=gege DDA(5,5)=gege DDA(6,6)=GEGE ELSE IF ((MFR.EQ.3.AND.NPINT.EQ.0).OR.MFR.EQ.9) THEN AUX=YUNG/(1.D0-XNU*XNU) AUX1=AUX*XNU DO J=1,NCONT DO I=1,NCONT DDAUX(I,J)=0.D0 ENDDO ENDDO C C CAS TRIDIMENSIONNEL ET FOURIER C IF(IFOUR2.EQ.2.OR.IFOUR2.EQ.1) THEN GEGE=0.5D0*YUNG/(1.D0+XNU) DDAUX(1,1)=AUX DDAUX(2,1)=AUX1 DDAUX(1,2)=AUX1 DDAUX(2,2)=AUX DDAUX(3,3)=GEGE DDAUX(4,4)=AUX DDAUX(5,4)=AUX1 DDAUX(4,5)=AUX1 DDAUX(5,5)=AUX DDAUX(6,6)=GEGE C C CAS AXISYMETRIQUE ET DEFORMATIONS PLANES C ELSE IF(IFOUR2.EQ.0.OR.IFOUR2.EQ.-1.OR.IFOUR2.EQ.-3) THEN DDAUX(1,1)=AUX DDAUX(2,1)=AUX1 DDAUX(1,2)=AUX1 DDAUX(2,2)=AUX DDAUX(3,3)=AUX DDAUX(4,3)=AUX1 DDAUX(3,4)=AUX1 DDAUX(4,4)=AUX C C CAS CONTRAINTES PLANES C ELSE IF(IFOUR2.EQ.-2) THEN DDAUX(1,1)=YUNG DDAUX(3,3)=YUNG ENDIF ENDIF * DO I=1,nstrsl S(I)=SIG0(I)+DSIGT(I) SIGB(I)=S(I) SX(I)=S(I)-SPHER(I) ENDDO C---------CAS DES POUTRES----------------------------------------------- IF(MFRl.EQ.7) THEN ** write(6,*) 'work',(work(ic),ic=1,LW) DIV(1)=1.D0/WORK(4) DIV(2)=1.D0 DIV(3)=1.D0 DIV(4)=WORK(10)/WORK(1) DIV(5)=WORK(11)/WORK(2) DIV(6)=WORK(12)/WORK(3) IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(WORK(1)*WORK(4)) IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(WORK(2)*WORK(4)) IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(WORK(3)*WORK(4)) DO I=1,NCONT S(I) = S(I) *DIV(I) SX(I) = SX(I)*DIV(I) SPHER(I) = SPHER(I) *DIV(I) SPHER1(I) = SPHER1(I)*DIV(I) SPHER2(I) = SPHER2(I)*DIV(I) END DO ENDIF ** if (ib.eq.1 .AND.igau.eq.1) then ** write(ioimp,*) ' SIGi',NSTRSS ** write(ioimp,*) ' S',(s(iof),iof=1,NSTRSS) ** write(ioimp,*) ' SX',(sx(iof),iof=1,NSTRSS) ** endif C----------------------------------------------------------------------- C CALCUL DE LA LIMITE ELASTIQUE SI C----------------------------------------------------------------------- IF (jNPLAS.EQ.1) THEN SI=TRAC(1) C* Modele a ECROUISSAGE CINEMATIQUE LINEAIRE : la limite d'elasticite C* est la limite initiale donnee par SIGY ELSE IF (jNPLAS.EQ.4) THEN SI = SLI0 ELSE IF (jNPLAS.EQ.5.OR.jNPLAS.EQ.87) THEN CALL TRACTI(SI,EPST,SIG,EPS,LTRAC2,2,izz) IF (izz.EQ.1) THEN KERRE=75 GOTO 1000 ENDIF C* Modele de CHABOCHE (prise en compte ecrouissage) ELSE IF (jNPLAS.EQ.7 .OR. jNPLAS.EQ.12) THEN RMmRR = (RM - R0) * EXP(-B*EPST) SI = RM - RMmRR ENDIF c*dbg if (ib.eq.1 .AND.igau.eq.1) then c*dbg write(ioimp,*) 'Limite d elasticite SI=',SI,TRAC(1),EPST c*dbg endif ** * kich : application pression limite trace sigma ** yxsxii = 0.d0 if (icow21.gt.0) then ytr = (sx(1) + sx(2) + sx(3))/3.D0 if (ytr.gt.xlimp) then yxsxii = ytr - xlimp else r_z = ytr + xlimp if (r_z.lt.0.D0) yxsxii = r_z endif if (yxsxii.ne.0.d0) then dsigii = dsigt(1) + dsigt(2) + dsigt(3) if (dsigii.ne.0.D0) then r_z = 3.D0*yxsxii/dsigii do jj = 1,3 s(jj) = s(jj) - (dsigt(jj)*r_z) sx(jj) = sx(jj) - (dsigt(jj)*r_z) enddo else do jj = 1,3 s(jj) = s(jj) - yxsxii sx(jj) = sx(jj) - yxsxii enddo endif endif endif C----------------------------------------------------------------------- C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ C----------------------------------------------------------------------- * attention en contraintes planes on se declare en 3D * rien besoin de faire dans vonmis0 car IFOUR2 n'est pas utilise * et les contraintes sont dimensionnees a 6 SEQ=VONMIS0(SX,nstrsl,MFRl,IFOUR2,EP1,ALFA1) ** if (ib.eq.1 .AND.igau.eq.1) then ** write(ioimp,*) 'VONMISES',SEQ,MFRl,IFOUR2,EP1,ALFA1 ** endif C----------------------------------------------------------------------- C LE CRITERE EST-IL VERIFIE C----------------------------------------------------------------------- IF (SEQ.LE.SI) THEN ITRY = 1 ELSE PETI=1.1D0*0.5D0*PRECIS*SEQ CALL EPSPRE(SEQ,SI,PETI,ITRY) ENDIF IF (ITRY.EQ.1) THEN * rien a faire on n'a pas plastifie ** if (ib.eq.1 .and. igau.eq.1) then ** write(ioimp,*) 'pas de plastification' ** endif IF (MFRl.EQ.7) THEN DO I=1,NCONT S(I)=S(I)/DIV(I) ENDDO ENDIF DO I=1,NCONT SIGF(I)=S(I) DEFP(I)=0.D0 ENDDO IF(MFRl.EQ.9) THEN DEFP(7)=0.D0 DEFP(8)=0.D0 SIGF(7)=S(7) SIGF(8)=S(8) ENDIF ** if (ib.eq.1 .and. igau.eq.1) then ** write(ioimp,*) 'SIGF',(SIGF(iof),iof=1,nstrsl) ** endif VARF(1)=VAR0(1) IF(jNPLAS.EQ.4.OR.jNPLAS.EQ.7) THEN DO I=1,nstrsl VARF(I+1)=VAR0(I+1) ENDDO ELSE IF(jNPLAS.EQ.12) THEN DO I=1,2*nstrsl VARF(I+1)=VAR0(I+1) ENDDO ELSE IF(jNPLAS.EQ.87) THEN DO I=1,nstrsl VARF(I+1)=DEPST(I) ENDDO ENDIF c*dbg if (ib.eq.1 .and. igau.eq.1) then c*dbg write(ioimp,*) 'VARF',EPST,'---',(VARF(iof),iof=2,1+nstrsl) c*dbg endif RETURN ENDIF C----------------------------------------------------------------------- C ON A PLASTIFIE C----------------------------------------------------------------------- NITER=0 PHI=SEQ-SI c*dbg if (ib.eq.1 .and. igau.eq.1) then c*dbg write(ioimp,*) 'niter,phi,si,seq,peti,precis=', c*dbg $ niter,phi,si,seq,peti,precis c*dbg endif cri0= si * 1.D-8 PHI0=PHI SI0=SI RR=0.D0 DO I=1,NCONT DSIG(I) =0.D0 DEPSP(I) =0.D0 DSPHER1(I)=0.D0 DSPHER2(I)=0.D0 ENDDO C----------------------------------------------------------------------- C DEBUT DE LA BOUCLE D'ITERATIONS INTERNES C----------------------------------------------------------------------- sx1in=sx(1) sx2in=sx(2) sx3in=sx(3) s1in=s(1) s2in=s(2) iderin=0 10 CONTINUE NITER=NITER+1 phi= seq - si C---------CALCUL DE WF1=DF/D(SIGMA)-------------------------------------- C---------ELEMENTS MASSIFS---------------------------------------------- IF(MFRl.EQ.1.OR.MFRl.EQ.31) THEN F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0 F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0 F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0 DO I=4,nstrsl F(I)=SX(I) ENDDO DO I=1,3 WF1(I)=1.5D0*F(I)/SEQ Z1(I)=WF1(I) ENDDO DO I=4,NCONT WF1(I)=3.D0*F(I)/SEQ Z1(I)=1.5D0*F(I)/SEQ ENDDO C---------COQUES MINCES------------------------------------------------- ELSE IF(MFRl.EQ.3.OR.MFRl.EQ.9) THEN IF(IFOUR2.GE.1) THEN WF1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1) WF1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1) WF1(3)=3.D0*SX(3)/(SEQ*EP1) WF1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(5))/(SEQ*EP1*EP1) WF1(5)=3.D0*WORK(2)*(2.D0*SX(5)-SX(4))/(SEQ*EP1*EP1) WF1(6)=18.D0*WORK(2)*SX(6)/(SEQ*EP1*EP1) ELSE WF1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1) WF1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1) WF1(3)=3.D0*WORK(2)*(2.D0*SX(3)-SX(4))/(SEQ*EP1*EP1) WF1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(3))/(SEQ*EP1*EP1) ENDIF C---------COQUES EPAISSES----------------------------------------------- ELSE IF(MFRl.EQ.5) THEN F(1)=(2.D0*SX(1)-SX(2))/3.D0 F(2)=(2.D0*SX(2)-SX(1))/3.D0 DO I=3,5 F(I)=SX(I) ENDDO DO I=1,2 WF1(I)=1.5D0*F(I)/SEQ Z1(I)=WF1(I) ENDDO DO I=3,5 WF1(I)=3.D0*F(I)/SEQ Z1(I)=1.5D0*F(I)/SEQ ENDDO C---------POUTRES------------------------------------------------------- ELSE IF(MFRl.EQ.7) THEN DO J=1,NCONT DO I=1,NCONT DDA(I,J)=0.D0 ENDDO ENDDO DDA(1,1)=YUNG DDA(4,4)=0.5D0*YUNG/(1.D0+XNU) DDA(5,5)=YUNG DDA(6,6)=YUNG WF1(1)=SX(1)/SEQ WF1(2)=0.D0 WF1(3)=0.D0 WF1(4)=SX(4)/SEQ WF1(5)=SX(5)/SEQ WF1(6)=SX(6)/SEQ ENDIF DO I=1,NCONT WF2(I)=0.D0 ENDDO IF(MFRl.EQ.7) THEN DO J=1,NCONT XFLO1=WF1(J) DO I=1,NCONT WF2(I)=WF2(I)+XFLO1*DDA(I,J) ENDDO ENDDO ELSE IF((mfr.eq.1.or.mfr.eq.31).and. IFOUR2.eq.-2) then DO J=1,NCONT XFLO1=WF1(J) DO I=1,NCONT WF2(I)=WF2(I)+XFLO1*DDA(I,J) ENDDO ENDDO else DO J=1,NCONT XFLO1=WF1(J) DO I=1,NCONT WF2(I)=WF2(I)+XFLO1*DDAUX(I,J) ENDDO ENDDO endif ENDIF COEF=0.D0 DO I=1,NCONT COEF=COEF+WF1(I)*WF2(I) ENDDO C----------------------------------------------------------------------- C PLASTICITE PARFAITE, ECROUISSAGE ISOTROPE ET CINEMATIQUE ZIEGLER C----------------------------------------------------------------------- IF (jNPLAS.EQ.1.OR.jNPLAS.EQ.4.OR.jNPLAS.EQ.5 $ .OR.jNPLAS.EQ.87) THEN IF (jNPLAS.EQ.4) THEN RP=0.D0 C=HECL ELSE CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC2,1,izz) IF (izz.EQ.1) THEN KERRE=75 GOTO 1000 ENDIF IF (jNPLAS.EQ.1) THEN RP=0.D0 C=0.D0 ELSE IF(jNPLAS.EQ.5.OR.jNPLAS.EQ.87) THEN RP=PENTE C=0.D0 ENDIF ENDIF c*dbg if (ib.eq.1 .and. igau.eq.1) then c*dbg write(ioimp,*) 'RP,C=',RP,C,EPST c*dbg endif DENOM=COEF+C+RP DELTA=PHI/DENOM DMU=C*DELTA/SEQ DO I=1,NCONT DSIG(I)=-DELTA*WF2(I) DSPHER1(I)=DMU*SX(I) ENDDO * Cas des contraintes planes en massif if ((mfr.eq.1.or.mfr.eq.31).and.IFOUR2.eq.-2) then bb= abs(dsig(3)+ s(3) ) r_z = dsig(3) * cte_cp sx(3)= sx3in - dsig(3) sx(1)= sx1in - r_z sx(2)= sx2in - r_z SEQ=VONMIS0(SX,nstrsl,MFRl,IFOUR2,EP1,ALFA1) s(3)= - dsig(3) s(1)= s1in - r_z s(2)= s2in - r_z if (bb.gt.cri0) then if (iderin.eq.0) then niter=niter - 1 endif iderin=iderin+1 if (iderin.gt.50) then write(ioimp,*) ' probleme dans iterations internes' KERRE=2 GO TO 1000 endif go to 10 endif DMU=C*DELTA/SEQ DO I=1,NCONT DSPHER1(I)=DMU*SX(I) ENDDO endif iderin=0 DP=DELTA DR=RP*DP ELSE C----------------------------------------------------------------------- C MODELE DE CHABOCHE C----------------------------------------------------------------------- C---------UNIQUEMENT POUR LES ELEMENTS MASSIFS-------------------------- XPRO1=0.D0 XPRO2=0.D0 DO I=1,NCONT XPRO1=XPRO1+WF1(I)*SPHER1(I) XPRO2=XPRO2+WF1(I)*SPHER2(I) ENDDO FIP=1.D0+(PSI-1.D0)*EXP(-OME*EPST) DENOM=COEF+(A1*C1+A2*C2)*FIP-C1*XPRO1-C2*XPRO2+B*RMmRR DELTA=PHI/DENOM DO I=1,NCONT DSIG(I)=-DELTA*WF2(I) DSPHER1(I)=(2.D0*A1*FIP*Z1(I)/3.D0-SPHER1(I))*C1*DELTA DSPHER2(I)=(2.D0*A2*FIP*Z1(I)/3.D0-SPHER2(I))*C2*DELTA ENDDO DR=B* RMmRR *DELTA DP=DELTA ENDIF RR=RR+DR EPST=EPST+DP EPST=MAX(0.D0,EPST) c*dbg if (ib.eq.1 .and. igau.eq.1) then c*dbg write(ioimp,*) 'MAJ=',EPST,DP,RR,DR c*dbg endif IF (MFRl.EQ.3.OR.MFRl.EQ.9) THEN r_z = EP1*EP1/(6.D0*WORK(2)) IF (IFOUR2.GE.1) THEN DSIG(1)=DSIG(1)*EP1 DSIG(2)=DSIG(2)*EP1 DSIG(3)=DSIG(3)*EP1 DSIG(4)=DSIG(4)*r_z DSIG(5)=DSIG(5)*r_z DSIG(6)=DSIG(6)*r_z ELSE DSIG(1)=DSIG(1)*EP1 DSIG(2)=DSIG(2)*EP1 DSIG(3)=DSIG(3)*r_z DSIG(4)=DSIG(4)*r_z ENDIF ENDIF C mise a jour des contraintes DO I=1,NCONT S(I)=S(I)+DSIG(I) SPHER1(I)=SPHER1(I)+DSPHER1(I) SPHER2(I)=SPHER2(I)+DSPHER2(I) SPHER(I)=SPHER1(I)+SPHER2(I) SX(I)=S(I)-SPHER(I) ENDDO if(IFOUR2.eq.-2.and.(mfr.eq.1.or.mfr.eq.31)) then s(3)=0.d0 endif if (icow21.gt.0) then yxsxii = 0.D0 * kich borne trace sigma ytr = (sx(1) + sx(2) + sx(3))/3.D0 if (ytr.gt.xlimp) then yxsxii = ytr - xlimp else r_z = ytr + xlimp if (r_z.lt.0.D0) yxsxii = r_z endif dsigii = dsigt(1) + dsigt(2) + dsigt(3) if (dsigii.ne.0.D0) then r_z = 3.D0*yxsxii/dsigii do jj = 1,3 s(jj) = s(jj) - (dsigt(jj)*r_z) sx(jj) = sx(jj) - (dsigt(jj)*r_z) enddo else do jj = 1,3 s(jj) = s(jj) - yxsxii sx(jj) = sx(jj) - yxsxii enddo endif endif SEQ=VONMIS0(SX,nstrsl,MFRl,IFOUR2,EP1,ALFA1) C---------CONTRAINTES PLANES-------------------------------------------- IF(IFOUR2.EQ.-2) THEN IF(MFRl.EQ.1.OR.MFRl.EQ.31) THEN F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0 F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0 F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0 DO I=4,nstrsl F(I)=SX(I) ENDDO DO I=1,3 WF1(I)=1.5D0*F(I)/SEQ ENDDO DO I=4,nstrsl WF1(I)=3.D0*F(I)/SEQ ENDDO ELSE IF(MFRl.EQ.3.OR.MFRl.EQ.9) THEN AUX=EP1*EP1*EP1*EP1 WF1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1*EP1) WF1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1*EP1) WF1(3)=18.D0*ALFA1*(2.D0*SX(3)-SX(4))/(SEQ*AUX) WF1(4)=18.D0*ALFA1*(2.D0*SX(4)-SX(3))/(SEQ*AUX) ENDIF DO I=1,nstrsl DEPSP(I)=DEPSP(I)+DELTA*WF1(I) ENDDO ENDIF C----------------------------------------------------------------------- C TEST C CALCUL DE LA NOUVELLE VALEUR DE PHI C----------------------------------------------------------------------- IF (jNPLAS.EQ.5.OR.jNPLAS.EQ.87) THEN CALL TRACTI(SI,EPST,SIG,EPS,LTRAC2,2,izz) C* Modele de CHABOCHE (prise en compte ecrouissage) ELSE IF (jNPLAS.EQ.7 .OR. jNPLAS.EQ.12) THEN RMmRR = (RM - R0) * EXP(-B*EPST) SI = RM - RMmRR ELSE IF (jNPLAS.EQ.4) THEN SI=SLI0 ELSE SI=RR+SI0 ENDIF PHI=SEQ-SI PETI=1.D-8 APHI=ABS(PHI) APHI0=ABS(PHI0) TEST=max(PETI*APHI0,XZPREC*100.D0*SEQ) *sg TEST=PETI*APHI0 *sg write(ioimp,*) 'niter,phi,phi0,si,seq,rmmrr,test=', *sg $ niter,phi,phi0,si,seq,rmmrr,test IF(NITER.GT.50) THEN if(itracb.eq.0) then itracb=1 go to 2222 endif C write(ioimp,*) ' CCOIN0, APHI, TEST',APHI,TEST C write(ioimp,*) ' CCOIN0, KERRE = 2, NITER =',NITER KERRE=-2 call soucis(268) C GO TO 1000 ENDIF IF(APHI.LE.TEST.OR.KERRE.LT.0) THEN IF (KERRE.LT.0) KERRE = 0 IF (MFRl.EQ.7) THEN DO I=1,NCONT S(I)=S(I)/DIV(I) SPHER1(I)=SPHER1(I)/DIV(I) SPHER2(I)=SPHER2(I)/DIV(I) ENDDO ENDIF C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES------------------- IF(IFOUR2.NE.-2) THEN DO I=1,NCONT DS(I)=S(I)-SIGB(I) ENDDO CALL EPSIG0(DS,DDEPSE,MFRl,IFOUR2,YUNG,XNU,WORK,nstrsl) DO I=1,NCONT DEPSE(I)=DEPST(I)+DDEPSE(I) DEPSP(I)=DEPST(I)-DEPSE(I) ENDDO ENDIF DO I=1,nstrsl SIGF(I)=S(I) DEFP(I)=DEPSP(I) ENDDO c*dbg if (ib.eq.1 .and. igau.eq.1) then c*dbg write(ioimp,*) 'SIGF',(SIGF(iof),iof=1,nstrsl) c*dbg endif C---------COQUES AVEC CISAILLEMENT TRANSVERSE--------------------------- IF(MFRl.EQ.9) THEN DEFP(7)=0.D0 DEFP(8)=0.D0 SIGF(7)=SIGB(7) SIGF(8)=SIGB(8) ENDIF VARF(1)=EPST IF(jNPLAS.EQ.4.OR.jNPLAS.EQ.7.OR.jNPLAS.EQ.12) THEN DO I=1,nstrsl VARF(I+1)=SPHER1(I) ENDDO IF(jNPLAS.EQ.12) THEN DO I=1,nstrsl VARF(nstrsl+1+I)=SPHER2(I) ENDDO ENDIF ELSE IF(jNPLAS.EQ.87) THEN DO I=1,nstrsl VARF(1+I)=DEPST(I) ENDDO ENDIF KERRE=0 c*dbg if (ib.eq.1 .and. igau.eq.1) then c*dbg write(ioimp,*) 'VARF',EPST,'---',(VARF(iof),iof=2,1+nstrsl) c*dbg endif RETURN ELSE sx1in=sx(1) sx2in=sx(2) s1in=s(1) s2in=s(2) GOTO 10 ENDIF 1000 CONTINUE C RETURN END