C CCOIN0 SOURCE FANDEUR 11/04/12 21:15:09 6938 SUBROUTINE CCOIN0(wrk52,wrk53,wrk54,wrk2,wrk3, & IB,IGAU,NBPGAU,LTRAC,IFOURB,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 C NSTRSS NOMBRE DE CONTRAINTES CA2000 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 IFOU = OPTION DE CALCUL 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 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 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) ALFAH 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----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC CCOPTIO -INC DECHE SEGMENT IECOU INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7, 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16, 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24, 3 icow25,icow26,icow27,icow28,icow29,icow30,icow31, 4 icow32,icow33,NSTRSS,MFR1, NBGMAT,NELMAT,icow38, 5 icow39,icow40,icow41,icow42,icow43,icow44 INTEGER icow45,icow46,icow47,icow48,icow49,icow50, . icow51,icow52,icow53,icow54,icow55,icow56 . icow57,icow58 ENDSEGMENT * SEGMENT WRK2 REAL*8 TRAC(LTRAC) 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),W1(8),W2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8) DIMENSION CRIGI(12) C---------COQUES AVEC CT------------------------------------------------ C ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES IF(MFR1.EQ.9) THEN NCONT=6 ELSE NCONT=NSTRSS ENDIF itracb=0 2222 continue C----------------------------------------------------------------------- DO I=1,NSTRSS S(I)=0.D0 SPHER(I)=0.D0 SPHER1(I)=0.D0 SPHER2(I)=0.D0 ENDDO KERRE=0 YUNG=XMAT(1) XNU=XMAT(2) C---------CARACTERISTIQUES GEOMETRIQUES--------------------------------- C C COQUES C ALFAH=1.D0 IF(MFR1.EQ.3.OR.MFR1.EQ.5.OR.MFR1.EQ.9) THEN EP1=WORK(1) IF(MFR1.NE.5) ALFAH=WORK(2)**2 ENDIF C---------COEFFICIENTS CHABOCHE----------------------------------------- IF(INPLAS.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(INPLAS.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 DO I=1,LTRAC SIG(I)=TRAC(2*I-1) EPS(I)=TRAC(2*I) ENDDO ENDIF EPST=VAR0(1) C---------ECROUISSAGE CINEMATIQUE--------------------------------------- IF(INPLAS.EQ.4.OR.INPLAS.EQ.7.OR.INPLAS.EQ.12) THEN DO I=1,NSTRSS SPHER1(I)=VAR0(I+1) ENDDO IF(INPLAS.EQ.12) THEN DO I=1,NSTRSS SPHER2(I)=VAR0(NSTRSS+1+I) ENDDO ENDIF DO I=1,NSTRSS SPHER(I)=SPHER1(I)+SPHER2(I) ENDDO 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(INPLAS.EQ.87) THEN EPST=0.D0 DO I=1,NSTRSS SIG0(I)=0.D0 DEPST(I)=DEPST(I) + VAR0(I+1) ENDDO ENDIF CALL CALSIG(DEPST,DDAUX,NSTRSS,CMATE,VALMAT,VALCAR, & N2EL,N2PTEL,MFR1,IFOURB,IB,IGAU,EPAIST, & NBPGAU,MELE,NPINT,NBGMAT,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.and. ifourb.eq.-2) then * en cas de contraintes planes on repasse en 3D do iu=1,6 do ju=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 I=1,NCONT DO J=1,NCONT DDAUX(I,J)=0.D0 ENDDO ENDDO C C CAS TRIDIMENSIONNEL ET FOURIER C IF(IFOURB.EQ.2.OR.IFOURB.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(IFOURB.EQ.0.OR.IFOURB.EQ.-1.OR.IFOURB.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(IFOURB.EQ.-2) THEN DDAUX(1,1)=YUNG DDAUX(3,3)=YUNG ENDIF ENDIF * DO I=1,NSTRSS S(I)=SIG0(I)+DSIGT(I) SIGB(I)=S(I) SX(I)=S(I)-SPHER(I) ENDDO * if( igau.eq.1.and.itracb.eq.1) then * write(6,*) ' entre dans cconi0' * write(6,*) ' sig0' , sig0(1),sig0(2),sig0(3),sig0(4) * write(6,*) ' dsigt' , dsigt(1),dsigt(2),dsigt(3),dsigt(4) * write(6,*) ' s' , s(1),s(2),s(3),s(4) * write(6,*) ' sx' , sx(1),sx(2),sx(3),sx(4) * write(6,*) ' spher' , spher(1),spher(2),spher(3),spher(4) * endif C---------CAS DES POUTRES----------------------------------------------- IF(MFR1.EQ.7) THEN 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) ENDDO ENDIF C----------------------------------------------------------------------- C CALCUL DE LA LIMITE ELASTIQUE SI C----------------------------------------------------------------------- IF(INPLAS.EQ.1.OR.INPLAS.EQ.4) THEN SI=TRAC(1) ELSE IF(INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI) IF(IBI.EQ.1) THEN KERRE=75 GOTO 1000 ENDIF C* Modele de CHABOCHE (prise en compte ecrouissage) ELSE IF (INPLAS.EQ.7 .OR. INPLAS.EQ.12) THEN RMmRR = (RM - R0) * EXP(-B*EPST) SI = RM - RMmRR 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 ifourb n'est pas utilisé * et les contraintes sont dimensionnées à 6 * if(igau.eq.1) then * write(6,*) ' en entrée sx ',sx(1),sx(2),sx(3) * endif SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) C----------------------------------------------------------------------- C LE CRITERE EST-IL VERIFIE C----------------------------------------------------------------------- * if(itracb.eq.1) then * write(6,*)' en entrée si seq', si, seq * endif PHI=SEQ-SI NITER=0 PETI=1.1D0*0.5D0*PRECIS*SEQ CALL EPSPRE(SEQ,SI,PETI,ITRY) IF((ITRY.EQ.1).OR.(SEQ.LE.SI)) THEN * rien a faire on n'a pas plastifié IF(MFR1.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(MFR1.EQ.9) THEN DEFP(7)=0.D0 DEFP(8)=0.D0 SIGF(7)=S(7) SIGF(8)=S(8) ENDIF VARF(1)=VAR0(1) IF(INPLAS.EQ.4.OR.INPLAS.EQ.7) THEN DO I=1,NSTRSS VARF(I+1)=VAR0(I+1) ENDDO ELSE IF(INPLAS.EQ.12) THEN DO I=1,2*NSTRSS VARF(I+1)=VAR0(I+1) ENDDO ELSE IF(INPLAS.EQ.87) THEN DO I=1,NSTRSS VARF(I+1)=DEPST(I) ENDDO ENDIF RETURN ENDIF C----------------------------------------------------------------------- C ON A PLASTIFIE C----------------------------------------------------------------------- cri0= si * 1.e-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 NITER=NITER+1 phi= seq - si c* if(itracb.eq.1) then c* write(6,*) ' niter seq si phi ' , niter,seq,si ,phi c* endif C---------CALCUL DE W1=DF/D(SIGMA)-------------------------------------- C---------ELEMENTS MASSIFS---------------------------------------------- IF(MFR1.EQ.1) 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,NSTRSS F(I)=SX(I) ENDDO DO I=1,3 W1(I)=1.5D0*F(I)/SEQ Z1(I)=W1(I) ENDDO DO I=4,NCONT W1(I)=3.D0*F(I)/SEQ Z1(I)=1.5D0*F(I)/SEQ ENDDO c* if(itracb.eq.1) then c* write(6,*)' W1' , w1(1), w1(2), w1(3), w1(4) c* endif C---------COQUES MINCES------------------------------------------------- ELSE IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN IF(IFOURB.GE.1) THEN W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1) W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1) W1(3)=3.D0*SX(3)/(SEQ*EP1) W1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(5))/(SEQ*EP1*EP1) W1(5)=3.D0*WORK(2)*(2.D0*SX(5)-SX(4))/(SEQ*EP1*EP1) W1(6)=18.D0*WORK(2)*SX(6)/(SEQ*EP1*EP1) ELSE W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1) W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1) W1(3)=3.D0*WORK(2)*(2.D0*SX(3)-SX(4))/(SEQ*EP1*EP1) W1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(3))/(SEQ*EP1*EP1) ENDIF C---------COQUES EPAISSES----------------------------------------------- ELSE IF(MFR1.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 W1(I)=1.5D0*F(I)/SEQ Z1(I)=W1(I) ENDDO DO I=3,5 W1(I)=3.D0*F(I)/SEQ Z1(I)=1.5D0*F(I)/SEQ ENDDO C---------POUTRES------------------------------------------------------- ELSE IF(MFR1.EQ.7) THEN DO I=1,NCONT DO J=1,NCONT DDA(I,J)=0.D0 ENDDO ENDDO DDA(1,1)=YUNG DDA(4,4)=0.5*YUNG/(1.+XNU) DDA(5,5)=YUNG DDA(6,6)=YUNG W1(1)=SX(1)/SEQ W1(2)=0.D0 W1(3)=0.D0 W1(4)=SX(4)/SEQ W1(5)=SX(5)/SEQ W1(6)=SX(6)/SEQ ENDIF IF(MFR1.EQ.7) THEN DO I=1,NCONT r_z = 0.D0 DO J=1,NCONT r_z=r_z+DDA(I,J)*W1(J) ENDDO W2(I)=r_z ENDDO ELSE IF(mfr.eq.1.and. ifourb.eq.-2) then DO I=1,NCONT r_z = 0.D0 DO J=1,NCONT r_z=r_z+DDA(I,J)*W1(J) ENDDO W2(I)=r_z ENDDO else DO I=1,NCONT r_z = 0.D0 DO J=1,NCONT r_z=r_z+DDAUX(I,J)*W1(J) ENDDO W2(I)=r_z ENDDO endif ENDIF COEF=0.D0 DO I=1,NCONT COEF=COEF+W1(I)*W2(I) ENDDO C----------------------------------------------------------------------- C PLASTICITE PARFAITE, ECROUISSAGE ISOTROPE ET CINEMATIQUE ZIEGLER C----------------------------------------------------------------------- IF(INPLAS.EQ.1.OR.INPLAS.EQ.4.OR.INPLAS.EQ.5 $ .OR.INPLAS.EQ.87) THEN CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC,1,IBI) * write(6,*) ' pente epst', pente,epst IF(IBI.EQ.1) THEN KERRE=75 GOTO 1000 ENDIF IF(INPLAS.EQ.1) THEN RP=0.D0 C=0.D0 ELSE IF(INPLAS.EQ.4) THEN RP=0.D0 C=PENTE ELSE IF(INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN RP=PENTE C=0.D0 ENDIF DENOM=COEF+C+RP DELTA=PHI/DENOM DMU=C*DELTA/SEQ * if(itracb.eq.1) then * write(6,*)'pente coef denom delta dmu ' * write(6,*)pente ,coef , denom ,delta ,dmu * endif DO I=1,NCONT DSIG(I)=-DELTA*W2(I) DSPHER1(I)=DMU*SX(I) ENDDO * Cas des contraintes planes en massif if(mfr.eq.1.and.ifourb.eq.-2) then * if(itracb.eq.1) then * write(6,*)'niter dsig delta',niter,dsig(1),dsig(2),dsig(3),delta * endif bb= abs(dsig(3)+ s(3) ) * if(itracb.eq.1) then * write(6,*) ' iderin bb ' , iderin , bb * endif r_z = dsig(3) * cte_cp sx(3)= sx3in - dsig(3) sx(1)= sx1in - r_z sx(2)= sx2in - r_z SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) 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 * if(itracb.eq.1) then * write(6,* ) ' cri0 bb seq delta ',cri0,bb,seq , delta * endif iderin=iderin+1 ** write(6,*) ' iderin ' , iderin 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+W1(I)*SPHER1(I) XPRO2=XPRO2+W1(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*W2(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 IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN IF(IFOURB.GE.1) THEN DO I=1,3 DSIG(I)=DSIG(I)*EP1 ENDDO DO I=4,6 DSIG(I)=EP1*EP1*DSIG(I)/(6.D0*WORK(2)) ENDDO ELSE DSIG(1)=DSIG(1)*EP1 DSIG(2)=DSIG(2)*EP1 DSIG(3)=EP1*EP1*DSIG(3)/(6.D0*WORK(2)) DSIG(4)=EP1*EP1*DSIG(4)/(6.D0*WORK(2)) ENDIF ENDIF C mise à 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(ifourb.eq.-2.and.mfr.eq.1) then s(3)=0.d0 endif * if(igau.eq.1) then * write(6,*) ' niter',niter * write(6,*) ' dsig' , dsig(1),dsig(2),dsig(3) * write(6,*) ' s' , s(1),s(2),s(3) * write(6,*) ' sx' , sx(1),sx(2),sx(3) * write(6,*) ' spher' , spher(1),spher(2),spher(3) * endif SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) C---------CONTRAINTES PLANES-------------------------------------------- IF(IFOURB.EQ.-2) THEN IF(MFR1.EQ.1) 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,NSTRSS F(I)=SX(I) ENDDO DO I=1,3 W1(I)=1.5D0*F(I)/SEQ ENDDO DO I=4,NSTRSS W1(I)=3.D0*F(I)/SEQ ENDDO ELSE IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN AUX=EP1*EP1*EP1*EP1 W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1*EP1) W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1*EP1) W1(3)=18.D0*ALFAH*(2.D0*SX(3)-SX(4))/(SEQ*AUX) W1(4)=18.D0*ALFAH*(2.D0*SX(4)-SX(3))/(SEQ*AUX) ENDIF DO I=1,NSTRSS DEPSP(I)=DEPSP(I)+DELTA*W1(I) ENDDO ENDIF C----------------------------------------------------------------------- C TEST C CALCUL DE LA NOUVELLE VALEUR DE PHI C----------------------------------------------------------------------- IF(INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI) C* Modele de CHABOCHE (prise en compte ecrouissage) ELSE IF (INPLAS.EQ.7 .OR. INPLAS.EQ.12) THEN RMmRR = (RM - R0) * EXP(-B*EPST) SI = RM - RMmRR ELSE SI=RR+SI0 ENDIF PHI=SEQ-SI PETI=1.E-7 APHI=ABS(PHI) APHI0=ABS(PHI0) TEST=PETI*APHI0 IF(NITER.GT.50) THEN if(itracb.eq.0) then itracb=1 go to 2222 endif KERRE=2 GO TO 1000 ENDIF * if(igau.eq.1) then * write(6 ,*)'test conver niter phi seq si testphi0',niter, * $ phi,seq,si,test,aphi0 * write(6,*) ' rr si0' , rr, si0 * endif IF(APHI.LE.TEST) THEN IF(MFR1.EQ.7) THEN DO I=1,NCONT S(I)=S(I)/DIV(I) ENDDO ENDIF C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES------------------- IF(IFOURB.NE.-2) THEN DO I=1,NCONT DS(I)=S(I)-SIGB(I) ENDDO CALL EPSIG0(DS,DDEPSE,MFR1,IFOURB,YUNG,XNU,WORK,NSTRSS) DO I=1,NCONT DEPSE(I)=DEPST(I)+DDEPSE(I) DEPSP(I)=DEPST(I)-DEPSE(I) ENDDO ENDIF DO I=1,NSTRSS SIGF(I)=S(I) DEFP(I)=DEPSP(I) ENDDO C---------COQUES AVEC CISAILLEMENT TRANSVERSE--------------------------- IF(MFR1.EQ.9) THEN DEFP(7)=0.D0 DEFP(8)=0.D0 SIGF(7)=SIGB(7) SIGF(8)=SIGB(8) ENDIF VARF(1)=EPST IF(INPLAS.EQ.4.OR.INPLAS.EQ.7.OR.INPLAS.EQ.12) THEN DO I=1,NSTRSS VARF(I+1)=SPHER1(I) ENDDO IF(INPLAS.EQ.12) THEN DO I=1,NSTRSS VARF(NSTRSS+1+I)=SPHER2(I) ENDDO ENDIF ENDIF IF(INPLAS.EQ.87) THEN DO I=1,NSTRSS VARF(1+I)=DEPST(I) ENDDO * write(6,*) ' sortie ecoin0' * write(6,*) 'depst',(depst(I),i=1,nstrs) * write(6,*) 'depsp',(depsp(i),i=1,nstrs) ENDIF KERRE=0 RETURN ELSE sx1in=sx(1) sx2in=sx(2) s1in=s(1) s2in=s(2) * if(igau.eq.1) then * write(6,*) ' fin iter', niter,' nouvel sx',sx(1),sx(2),sx(3) * endif GOTO 10 ENDIF C ENDIF 1000 CONTINUE RETURN END