C SYCO12 SOURCE OF166741 25/11/04 21:16:09 12349 SUBROUTINE SYCO12(wrk52,wrk53,wrk54,wrk2, & IB,IGAU,NBPGAU,NTRAC,IFOURL,iecou,DTPS) 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 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 = 63 XFEM C C ISOTROPE C INPLAS : NUMERO DU MODELE DE VISCOPLASTICITE C C INPLAS = 153 SYCO1 C INPLAS = 154 SYCO2 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*NTRAC) COURBE DE TRACTION C C MODELE ISOTROPE C VAR0(1) EPSE C VAR0(2) VP C C C----------------------------------------------------------------------- C 20/09/2017 : modif SG critere de convergence trop serre C TEST=PETI*APHI0 + utilisation CCREEL C idem ccoin0.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(NTRAC) ENDSEGMENT C----------------------------------------------------------------------- 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) KERRE=0 NCONT = iecou.NSTRSS MFRl = iecou.MFR1 jnplas = INPLAS C C--------PAS DE TEMPS ----------------------- TSYC9 = DTPS itracb=1 C----------------------------------------------------------------------- DO I=1, NCONT 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--------------------------------- EP1=1.D0 ALFAC=1.D0 IF (jnplas.EQ.153) THEN C WRITE(6,*) 'INPLAS', jnplas C WRITE(6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5) C WRITE(6,*) 'XMAT', xmat(6), xmat(7), xmat(8) PSYC9=XMAT(6) DSYC9=XMAT(7) ELSE IF (jnplas.EQ.154) THEN C WRITE(6,*) 'INPLAS', jnplas C WRITE(6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5) C WRITE(6,*) 'XMAT', xmat(6), xmat(7), xmat(8),xmat(9),xmat(10) PSYC9=XMAT(6) ASYC9=XMAT(7) BSYC9=XMAT(8) CSYC9=XMAT(9) ELSE ENDIF DO I=1,NTRAC SIG(I)=TRAC(2*I-1) EPS(I)=TRAC(2*I) ENDDO EPST=VAR0(1) C C VP=VAR0(2) C if(VP.LE.1.d-4) VP = 1.d-4 C VP0 = 1.D-4 C IF (jnplas.EQ.153) THEN C deno1 = 1.D0 + ( (VP0/DSYC9)**(1.D0/PSYC9) ) C deno2 = 1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) ) C else if (jnplas.EQ.154) then C HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) C deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9))) C deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9))) C endif C C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- C-xxxx pour le materiau syco1 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- C DSYC9 parametre D de la loi Symonds Cowper C PSYC9 parametre P de la loi Symonds Cowper C VP valeur minimum de la vitesse de deformation equivalente IF (jnplas.EQ.153) THEN VP0 = 1.D-4 VP = VP0 else if(jnplas.EQ.154) then C VP valeur minimum de la vitesse de deformation equivalente VP0 = 1.D-4 VP = VP0 endif C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- CALL CALSIG(DEPST,DDAUX,NCONT,CMATE,VALMAT,VALCAR, & N2EL,N2PTEL,MFRl,IFOURL,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 (ifourl.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 ENDIF DO I=1,NCONT S(I)=SIG0(I)+DSIGT(I) SIGB(I)=S(I) SX(I)=S(I)-SPHER(I) ENDDO C----------------------------------------------------------------------- C CALCUL DE LA LIMITE ELASTIQUE SI C----------------------------------------------------------------------- C CALL TRACTI(SI,EPST,SIG,EPS,NTRAC,2,izz) c IF (izz.EQ.1) THEN KERRE=75 GOTO 1000 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 ifourl n'est pas utilise * et les contraintes sont dimensionnees a 6 * if(igau.eq.1) then * WRITE(6,*) ' en entree sx ',sx(1),sx(2),sx(3) * endif C C SEQ=VONMIS0(SX,NCONT,MFRl,ifourl,EP1,ALFAC) C calcul de la contrainte de Von Mises SOM1=0.D0 DO I=4,NCONT SOM1=SOM1+SX(I)*SX(I) ENDDO SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3) & -SX(2)*SX(3)+3.D0*SOM1 SEQ=SQRT(ABS(SEQ)) C----------------------------------------------------------------------- C LE CRITERE EST-IL VERIFIE ? C----------------------------------------------------------------------- C if(itracb.eq.1) then C WRITE(6,*)' en entree si seq', si, seq C endif C c test de plasticite sur SI statique C SI = SI * deno2/deno1 C PHI=SEQ-SI NITER=0 PETI=1.1D0*0.5D0*PRECIS*SEQ CALL EPSPRE(SEQ,SI,PETI,ITRY) C IF((ITRY.EQ.1).OR.(SEQ.LE.SI)) THEN * rien a faire on n'a pas plastifie DO I=1,NCONT SIGF(I)=S(I) DEFP(I)=0.D0 ENDDO VARF(1)=VAR0(1) VARF(2)=VP0 RETURN ENDIF C----------------------------------------------------------------------- C ON A PLASTIFIE C----------------------------------------------------------------------- C 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---------------------------------------------- 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,NCONT 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 if(ifourl.eq.-2) then DO I=1,NCONT r_z = 0.D0 DO J=1,NCONT r_z=r_z+DDA(I,J)*WF1(J) ENDDO WF2(I)=r_z ENDDO else DO I=1,NCONT r_z = 0.D0 DO J=1,NCONT r_z=r_z+DDAUX(I,J)*WF1(J) ENDDO WF2(I)=r_z ENDDO endif COEF=0.D0 DO I=1,NCONT COEF=COEF+WF1(I)*WF2(I) ENDDO C----------------------------------------------------------------------- C ECROUISSAGE ISOTROPE C----------------------------------------------------------------------- CALL TRACTI(PENTE,EPST,SIG,EPS,NTRAC,1,izz) * WRITE(6,*) ' pente epst', pente,epst IF (izz.EQ.1) THEN KERRE=75 GOTO 1000 ENDIF RP=PENTE C=0.D0 DENOM=COEF+C C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- C---------CALCUL des termes en sup dus a la viscosite---------------- C il s'agit de -qsi1 = -d(phi)/d(vp) = +d(SY)/d(vp) c (ou vp = dp/dt cad la vitesse inelastiq equivalente) IF (jnplas.EQ.153) THEN deno1 = 1.D0 + ( (VP0/DSYC9)**(1.D0/PSYC9) ) deno2 = 1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) ) qsi1 = (SI/(PSYC9*DSYC9))*((VP/DSYC9)**( (1.D0/PSYC9)-1.D0)) qsi1 = qsi1 / (deno2 * TSYC9) DENOM = DENOM + qsi1 DENOM = DENOM + (RP * deno2/deno1) C ELSE IF (jnplas.EQ.154) THEN HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9))) deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9))) qsi1 = (SI*HSYC9/PSYC9) * (VP**((1.D0/PSYC9)-1.D0)) qsi1 = qsi1 / (deno2 * TSYC9) DENOM = DENOM + qsi1 DENOM = DENOM + (RP * deno2/deno1) C qsi2 = (SI/deno2)*(VP**(1.D0/PSYC9))*((ASYC9-HSYC9)/CSYC9) DENOM = DENOM + qsi2 C ENDIF C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- C23456789012345678901234567890123456789012345678901234567890123456789012 DELTA=PHI/DENOM DMU=C*DELTA/SEQ C if(itracb.eq.1) then C WRITE(6,*)'pente coef denom delta dmu niter si seq' C WRITE(6,*)pente ,coef , denom ,delta ,dmu,niter,si,seq C endif DO I=1,NCONT DSIG(I)=-DELTA*WF2(I) DSPHER1(I)=DMU*SX(I) ENDDO * Cas des contraintes planes en massif if(ifourl.eq.-2) then bb= abs(dsig(3)+ s(3) ) C if(itracb.eq.1) then C WRITE(6,*) ' iderin bb ' , iderin , bb C endif r_z = dsig(3) * cte_cp sx(3)= sx3in - dsig(3) sx(1)= sx1in - r_z sx(2)= sx2in - r_z C C SEQ=VONMIS0(SX,NCONT,MFRl,ifourl,EP1,ALFAC) C calcul de la contrainte de Von Mises SOM1=0.D0 DO I=4,NCONT SOM1=SOM1+SX(I)*SX(I) ENDDO SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3) & -SX(2)*SX(3)+3.D0*SOM1 SEQ=SQRT(ABS(SEQ)) C s(3)= - dsig(3) s(1)= s1in - r_z s(2)= s2in - r_z cri0= si * 1.d-8 if( bb.gt.cri0) then if(iderin.eq.0) then niter=niter - 1 endif C if(itracb.eq.1) then C WRITE(6,* ) ' cri0 bb seq delta ',cri0,bb,seq , delta C endif iderin = iderin + 1 ** WRITE(6,*) ' iderin ' , iderin if(iderin.gt.50) then WRITE(6,*) 'on iderin superieur a 50' write(6,*) ' probleme dans iterations internes' KERRE=2 GO TO 1000 endif C23456789012345678901234567890123456789012345678901234567890123456789012 C WRITE(6,*)'cri0 bb niter iderin',cri0,bb,niter,iderin,dsig(3),s(3) go to 10 endif DMU=C*DELTA/SEQ DO I=1,NCONT DSPHER1(I)=DMU*SX(I) ENDDO endif iderin=0 C DR=RP*DELTA RR=RR+DR EPST=EPST+DELTA 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(ifourl.eq.-2) then s(3)=0.d0 endif C if(igau.eq.1) then C WRITE(6,*) ' niter',niter C WRITE(6,*) ' dsig' , dsig(1),dsig(2),dsig(3),dsig(4) C WRITE(6,*) ' s' , s(1),s(2),s(3),s(4) C WRITE(6,*) ' sx' , sx(1),sx(2),sx(3),sx(4) C WRITE(6,*) ' spher' , spher(1),spher(2),spher(3),spher(4) C endif C SEQ=VONMIS0(SX,NCONT,MFRl,ifourl,EP1,ALFAC) C calcul de la contrainte de Von Mises SOM1=0.D0 DO I=4,NCONT SOM1=SOM1+SX(I)*SX(I) ENDDO SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3) & -SX(2)*SX(3)+3.D0*SOM1 SEQ=SQRT(ABS(SEQ)) C C---------CONTRAINTES PLANES-------------------------------------------- IF(ifourl.EQ.-2) 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,NCONT F(I)=SX(I) ENDDO DO I=1,3 WF1(I)=1.5D0*F(I)/SEQ ENDDO DO I=4,NCONT WF1(I)=3.D0*F(I)/SEQ ENDDO DO I=1,NCONT DEPSP(I)=DEPSP(I)+DELTA*WF1(I) ENDDO ENDIF C----------------------------------------------------------------------- C TEST C CALCUL DE LA NOUVELLE VALEUR DE PHI C----------------------------------------------------------------------- CALL TRACTI(SI,EPST,SIG,EPS,NTRAC,2,izz) C****************************************************** c-xxxx calcul de la valeur (i+1) de H, de la vitesse de def VP, c et de SI en prenant en compte la viscosite IF(jnplas.EQ.153) THEN C materiau SYCO1 deno1 = 1.D0 + ( (VP0/DSYC9)**(1.D0/PSYC9) ) VP = MAX( (VP+(DELTA/TSYC9)) , VP0 ) C deno2 = 1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) ) SI = SI * (1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) )) / deno1 C else if (jnplas.EQ.154) then C materiau SYCO2 HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9))) VP = MAX( (VP+(DELTA/TSYC9)) , VP0 ) C deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9))) SI = SI * (1.D0 + (HSYC9 * (VP**(1.D0/PSYC9)))) / deno1 endif C****************************************************** PHI=SEQ-SI PETI=1.D-7 APHI=ABS(PHI) APHI0=ABS(PHI0) TEST=max(PETI*APHI0,XZPREC*100.D0*SEQ) *sg TEST=PETI*APHI0 IF(NITER.GT.200) THEN WRITE(6,*) 'on a bien plante ici,car APHI=',APHI WRITE(6,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI TEST2=max(1.D-6*APHI0,XZPREC*100.D0*SEQ) *sg TEST2 = 1.D-6*APHI0 if(APHI.gt.TEST2) then KERRE=2 GO TO 1000 else TEST=TEST2 endif ELSEIF(NITER.eq.101) THEN WRITE(6,*) 'on aurait du plante ici,car APHI=',APHI WRITE(6,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI WRITE(6,*) 'mais on va continuer 50pas de plus' 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 C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES------------------- IF(ifourl.NE.-2) THEN DO I=1,NCONT DS(I)=S(I)-SIGB(I) ENDDO C copie d'un bout de EPSIG0 AUX1=1.D0/YUNG AUX2=2.D0*(1.D0+XNU) AUX3=AUX1*AUX2 DDEPSE(1)=(DS(1)-XNU*(DS(2)+DS(3)))*AUX1 DDEPSE(2)=(DS(2)-XNU*(DS(1)+DS(3)))*AUX1 DDEPSE(3)=(DS(3)-XNU*(DS(1)+DS(2)))*AUX1 DO I=4,NCONT DDEPSE(I)=AUX3*DS(I) ENDDO C fin copie d'un bout de EPSIG0 DO I=1,NCONT DEPSE(I)=DEPST(I)+DDEPSE(I) DEPSP(I)=DEPST(I)-DEPSE(I) ENDDO ENDIF DO I=1,NCONT SIGF(I)=S(I) DEFP(I)=DEPSP(I) ENDDO VARF(1)=EPST VARF(2)=VP KERRE=0 RETURN ELSE sx1in=sx(1) sx2in=sx(2) s1in=s(1) s2in=s(2) GOTO 10 ENDIF 1000 CONTINUE RETURN END