C SYCO12 SOURCE YURI 11/12/09 21:15:01 7199 C SYCO12 SOURCE FANDEUR 11/04/12 21:15:09 6938 SUBROUTINE SYCO12(wrk52,wrk53,wrk54,wrk2,wrk3, & IB,IGAU,NBPGAU,LTRAC,IFOURB,iecou,xecou) 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 = 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*LTRAC) COURBE DE TRACTION C C MODELE ISOTROPE C VAR0(1) EPSE C VAR0(2) VP C 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 * * Commun XECOU: sert de fourre-tout pour les initialisations * de réels * SEGMENT XECOU * COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00 REAL*8 DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00 ENDSEGMENT C COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00 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),W1(8),W2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8) DIMENSION CRIGI(12) NCONT=NSTRSS C C--------PAS DE TEMPS ----------------------- TSYC9 = DT itracb=1 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 EP1=1.D0 ALFAH=1.D0 IF(INPLAS.EQ.153) THEN C write(6,*) 'INPLAS', inplas 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) ENDIF IF(INPLAS.EQ.154) THEN C write(6,*) 'INPLAS', inplas 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) ENDIF DO I=1,LTRAC 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(INPLAS.EQ.153) THEN C deno1 = 1. + ( (VP0/DSYC9)**(1./PSYC9) ) C deno2 = 1. + ( (VP/DSYC9)**(1./PSYC9) ) C else if(INPLAS.EQ.154) then C HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) C deno1 = 1. + (HSYC9 * ((VP0)**(1./PSYC9))) C deno2 = 1. + (HSYC9 * ((VP)**(1./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(INPLAS.EQ.153) THEN VP0 = 1.D-4 VP = VP0 else if(INPLAS.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,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 (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 ENDIF * DO I=1,NSTRSS 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,LTRAC,2,IBI) c IF(IBI.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 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 C C SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) C calcul de la contrainte de Von Mises SOM1=0.D0 DO I=4,NSTRSS 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----------------------------------------------------------------------- C LE CRITERE EST-IL VERIFIE C----------------------------------------------------------------------- C if(itracb.eq.1) then C write(6,*)' en entrée 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 plastifié 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 NITER=NITER+1 phi= seq - si C---------CALCUL DE W1=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,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 if(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 COEF=0.D0 DO I=1,NCONT COEF=COEF+W1(I)*W2(I) ENDDO C----------------------------------------------------------------------- C ECROUISSAGE ISOTROPE C----------------------------------------------------------------------- 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 RP=PENTE C=0.D0 DENOM=COEF+C C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- C---------CALCUL des termes en sup dus a la viscosité---------------- C il s'agit de -qsi1 = -d(phi)/d(°p) = +d(SY)/d(°p) c (où °p = dp/dt cad la vitesse inelastiq equivalente) IF(INPLAS.EQ.153) THEN deno1 = 1. + ( (VP0/DSYC9)**(1./PSYC9) ) deno2 = 1. + ( (VP/DSYC9)**(1./PSYC9) ) qsi1 = (SI/(PSYC9*DSYC9))*((VP/DSYC9)**( (1./PSYC9)-1.)) qsi1 = qsi1 / (deno2 * TSYC9) DENOM = DENOM + qsi1 DENOM = DENOM + (RP * deno2/deno1) C ELSE IF (INPLAS.EQ.154) THEN HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) deno1 = 1. + (HSYC9 * ((VP0)**(1./PSYC9))) deno2 = 1. + (HSYC9 * ((VP)**(1./PSYC9))) qsi1 = (SI*HSYC9/PSYC9) * (VP**((1./PSYC9)-1.)) qsi1 = qsi1 / (deno2 * TSYC9) DENOM = DENOM + qsi1 DENOM = DENOM + (RP * deno2/deno1) C qsi2 = (SI/deno2)*(VP**(1./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*W2(I) DSPHER1(I)=DMU*SX(I) ENDDO * Cas des contraintes planes en massif if(ifourb.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,NSTRSS,MFR1,IFOURB,EP1,ALFAH) C calcul de la contrainte de Von Mises SOM1=0.D0 DO I=4,NSTRSS 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(ioimp,*) ' 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 à 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) 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,NSTRSS,MFR1,IFOURB,EP1,ALFAH) C calcul de la contrainte de Von Mises SOM1=0.D0 DO I=4,NSTRSS 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(IFOURB.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,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 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----------------------------------------------------------------------- CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI) 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 viscosité IF(INPLAS.EQ.153) THEN C materiau SYCO1 deno1 = 1. + ( (VP0/DSYC9)**(1./PSYC9) ) VP = MAX( (VP+(DELTA/TSYC9)) , VP0 ) C deno2 = 1. + ( (VP/DSYC9)**(1./PSYC9) ) SI = SI * (1. + ( (VP/DSYC9)**(1./PSYC9) )) / deno1 C else if (INPLAS.EQ.154) then C materiau SYCO2 HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) deno1 = 1. + (HSYC9 * ((VP0)**(1./PSYC9))) VP = MAX( (VP+(DELTA/TSYC9)) , VP0 ) C deno2 = 1. + (HSYC9 * ((VP)**(1./PSYC9))) SI = SI * (1. + (HSYC9 * (VP**(1./PSYC9)))) / deno1 endif C****************************************************** PHI=SEQ-SI PETI=1.D-7 APHI=ABS(PHI) APHI0=ABS(PHI0) TEST=PETI*APHI0 IF(NITER.GT.200) THEN write(*,*) 'on a bien planté ici,car APHI=',APHI write(*,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI TEST2 = 1.E-6*APHI0 if(APHI.gt.TEST2) then KERRE=2 GO TO 1000 else TEST=TEST2 endif ELSEIF(NITER.eq.101) THEN write(*,*) 'on aurait du planté ici,car APHI=',APHI write(*,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI write(*,*) '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(IFOURB.NE.-2) THEN DO I=1,NCONT DS(I)=S(I)-SIGB(I) ENDDO C CALL EPSIG0(DS,DDEPSE,MFR1,IFOURB,YUNG,XNU,WORK,NSTRSS) C copie d'un bout de EPSIG0 DO I=1,NSTRSS DDEPSE(I)=0.D0 ENDDO 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,NSTRSS 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,NSTRSS 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