C BEHAV3 SOURCE CB215821 26/08/24 21:15:15 12622 SUBROUTINE BEHAV3(SIGR,DSTRN,DEPS1,DEPS2,IPLA, & SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, & BETJEF,VISCO,NECH0,NECH1) C C ================================================================== C C MODELE DE PLASTICITE EN TRAITEMENT PRE ET POST PIC C Deux criteres de traction compression: C Rankine et Von Mises C C ================================================================== C CE SOUS-PROGRAMME EST APPELE DANS "BONE". C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION DFSIG1(4),DFSIG2(4),SIGF(4),DSIGT(4),VEC1(4),VEC2(4) DIMENSION V1(4),AC1(4),AC2(4),SIGE(4),SIGP(4) DIMENSION DLAM1(2),DLAM0(2),SIGR(4),IS1(2) DIMENSION FCRI0(2),FCRI1(2),DEPSI(4),DSTRN(4),SIGRV(4) DIMENSION D2FSI1(4,4),D2FSI2(4,4),DEP(4,4),P1(4,4),P2(4,4) DIMENSION D(4,4),A(4,4),AI(4,4),DJAC0(4,4),DJI(4,4),AH(4,4) C * SEGMENT BETJEF REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF, & TCON,DPSTF1,DPSTF2,TETA,PDT,TP0 INTEGER ICT,ICC,IMOD,IVIS,ITR, & ISIM,IBB,IGAU,IZON ENDSEGMENT SEGMENT VISCO REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV ENDSEGMENT SEGMENT NECH0 REAL*8 DT,DC,ALFG,S0,ENDO ENDSEGMENT SEGMENT NECH1 REAL*8 ENDL ENDSEGMENT * * * COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF, * & TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR, * & ISIM,IBB,IGAU,IZON C MNSTRS=NSTRS C IAPEX=0 Ft=PALF*COLI Rb=COLI PRB=1.D-10 PRB2=1.D-6 ITER=1 IEC1=0 IEC2=0 IBROY=0 ITANG=0 CRIMAX=0.D0 SEQ = 0.D0 SEQ1 = 0.D0 SEQ2 = 0.D0 CALL ZERO(SIGE,4,1) CALL ZERO(V1,4,1) CALL ZERO(P1,4,4) CALL ZERO(P2,4,4) CALL ZERO(D,4,4) C DO 10 I=1,NSTRS SIGE(I)=SIGF(I) 10 CONTINUE C IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN AD=YOUN/(1.D0-XNU*XNU) D(1,1)=AD D(2,2)=D(1,1) D(3,3)=AD*(1.D0-XNU)/2.D0 D(1,2)=AD*XNU D(2,1)=D(1,2) ENDIF C IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU)) D(1,1)=ADD*(1.D0-XNU) D(2,2)=D(1,1) D(3,3)=D(1,1) D(1,2)=ADD*XNU D(2,1)=D(1,2) D(1,3)=D(1,2) D(2,3)=D(1,2) D(3,1)=D(1,2) D(3,2)=D(1,2) D(4,4)=0.5*ADD*(1.D0-2.D0*XNU) ENDIF C C ************ Le point est fissure pour 1ere fois *********** C IF (IFIS.EQ.0) THEN CALL PRINC(SIGF,V1,NSTRS) DEPS1=0.D0 SIG1=Ft IF (V1(1).GT.Ft) THEN IFIS=1 ENDIF ENDIF IF (IPLA.EQ.0) THEN IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF) ENDIF IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF) ENDIF DEPS2=0.D0 SIG2=Rb*AA IF (SEQ.GT.SIG2) THEN IPLA=1 ENDIF ENDIF C 12 CONTINUE C CALL PRINC(SIGE,V1,NSTRS) TETA=V1(4) IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF) ENDIF IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF) ENDIF FCRI0(1)=V1(1)-SIG1 FCRI0(2)=SEQ-SIG2 IF (ABS(FCRI0(1)).GT.ABS(FCRI0(2))) THEN CRIMAX=ABS(100.D0*FCRI0(1)) ELSE CRIMAX=ABS(100.D0*FCRI0(2)) ENDIF IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN IEC1=1 IEC2=1 ENDIF IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN IEC1=1 IEC2=0 ENDIF IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN IEC1=0 IEC2=1 ENDIF C 15 CONTINUE C DO 16 I=1,NSTRS SIGF(I)=SIGE(I) 16 CONTINUE C IF (IEC1.EQ.1.AND.IEC2.EQ.1) THEN DO 101 I=1,NSTRS DO 17 J=1,NSTRS DEP(I,J)=D(I,J) 17 CONTINUE 101 CONTINUE GOTO 100 ENDIF C IF (IEC1.EQ.1.AND.IEC2.EQ.0) THEN CALL BEHAV2(SIGR,DSTRN,DEPS1,DEPS2,IPLA & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, & BETJEF,VISCO,NECH0,NECH1) GOTO 100 ENDIF C IF (IEC1.EQ.0.AND.IEC2.EQ.1) THEN CALL BEHAV1(SIGR,DSTRN,DEPS1,DEPS2,IPLA & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, & BETJEF,VISCO,NECH0,NECH1) GOTO 100 ENDIF C DK1=DEPS1 DK2=DEPS2 DLAM0(1)=0.D0 DLAM0(2)=0.D0 C 18 CONTINUE C C ************ Preparation des criteres ******************** C C ------------------ Von Mises ----------------------------- IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN A1=0.D0 A2=1.D0 ENDIF C ----------------- Drucker Prager ------------------------- IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN A1=(BETA-1.D0)/(2*BETA-1.D0) A2=BETA/(2*BETA-1.D0) ENDIF C ---------------------------------------------------------- C C ---------------- direction de traction ------------------ C C--------------- Determination du parametre d'ecrouissage ---------- C IF(IVIS.LE.2) THEN CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF) ELSE CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1) ENDIF C C------------------------------------------------------------ PI=4.D0*ATAN(1.D0) PHIC=V1(4)*(PI/180.D0) COSA=COS(PHIC) SINA=SIN(PHIC) C IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN DFSIG1(1)=COSA*COSA DFSIG1(2)=SINA*SINA DFSIG1(3)=2.D0*SINA*COSA DO 102 I=1,NSTRS DO 19 J=1,NSTRS P1(I,J)=0.D0 19 CONTINUE 102 CONTINUE P1(1,1)=1.D0/2.D0 P1(1,2)=-1.D0/2.D0 P1(2,1)=-1.D0/2.D0 P1(2,2)=1.D0/2.D0 P1(3,3)=2.D0 ENDIF C IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN DFSIG1(1)=COSA*COSA DFSIG1(2)=SINA*SINA DFSIG1(3)=0.D0 DFSIG1(4)=2.D0*SINA*COSA DO 103 I=1,NSTRS DO 21 J=1,NSTRS P1(I,J)=0.D0 21 CONTINUE 103 CONTINUE P1(1,1)=1.D0/2.D0 P1(1,2)=-1.D0/2.D0 P1(2,1)=-1.D0/2.D0 P1(2,2)=1.D0/2.D0 P1(4,4)=2.D0 ENDIF C DO 104 I=1,NSTRS AC1(I)=0.D0 DO 20 J=1,NSTRS AC1(I)=AC1(I)+D(I,J)*DFSIG1(J) 20 CONTINUE 104 CONTINUE C C ---------------- direction de compression ---------------- C IF(IVIS.LE.2) THEN CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF) ELSE CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) ENDIF C IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN SX=SIGF(1) SY=SIGF(2) SXY=SIGF(3) R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY R=SQRT(0.5*R1) IF (ABS(R).LE.1.D-10) THEN WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3' STOP ENDIF DFSIG2(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2 DFSIG2(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2 DFSIG2(3)=3.D0*SXY/R/A2 P2(1,1)=2.D0/A2/A2 P2(1,2)=-1.D0/A2/A2 P2(1,3)=0.D0/A2/A2 P2(2,1)=-1.D0/A2/A2 P2(2,2)=2.D0/A2/A2 P2(2,3)=0.D0/A2/A2 P2(3,1)=0.D0/A2/A2 P2(3,2)=0.D0/A2/A2 P2(3,3)=6.D0/A2/A2 ENDIF C IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN SX=SIGF(1) SY=SIGF(2) SZ=SIGF(3) SXY=SIGF(4) TX=(SX-SY)*(SX-SY) TY=(SX-SZ)*(SX-SZ) TZ=(SY-SZ)*(SY-SZ) R1=TX+TY+TZ+6*SXY*SXY R=SQRT(0.5*R1) IF (ABS(R).LE.1.D-10) THEN WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3' STOP ENDIF DFSIG2(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2 DFSIG2(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2 DFSIG2(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2 DFSIG2(4)=3.D0*SXY/R/A2 P2(1,1)=2.D0/A2/A2 P2(1,2)=-1.D0/A2/A2 P2(1,3)=-1.D0/A2/A2 P2(1,4)=0.D0/A2/A2 P2(2,1)=-1.D0/A2/A2 P2(2,2)=2.D0/A2/A2 P2(2,3)=-1.D0/A2/A2 P2(2,4)=0.D0/A2/A2 P2(3,1)=-1.D0/A2/A2 P2(3,2)=-1.D0/A2/A2 P2(3,3)=2.D0/A2/A2 P2(3,4)=0.D0/A2/A2 P2(4,1)=0.D0/A2/A2 P2(4,2)=0.D0/A2/A2 P2(4,3)=0.D0/A2/A2 P2(4,4)=6.D0/A2/A2 ENDIF C DO 105 I=1,NSTRS AC2(I)=0.D0 DO 25 J=1,NSTRS AC2(I)=AC2(I)+D(I,J)*DFSIG2(J) 25 CONTINUE 105 CONTINUE C C ************ Preparation du jacobien ******************** C F1DF1=0.D0 DO 30 J=1,NSTRS F1DF1=F1DF1+DFSIG1(J)*AC1(J) 30 CONTINUE F1DF2=0.D0 DO 31 J=1,NSTRS F1DF2=F1DF2+DFSIG1(J)*AC2(J) 31 CONTINUE F2DF2=0.D0 DO 32 J=1,NSTRS F2DF2=F2DF2+DFSIG2(J)*AC2(J) 32 CONTINUE F2DF1=0.D0 DO 33 J=1,NSTRS F2DF1=F2DF1+DFSIG2(J)*AC1(J) 33 CONTINUE C DJAC0(1,1)=-(F1DF1+PAECT) DJAC0(1,2)=-(F1DF2) DJAC0(2,2)=-(F2DF2+PAEC) DJAC0(2,1)=-(F2DF1) C DO 106 I=1,2 DO 43 J=1,2 DJI(I,J)=DJAC0(I,J) 43 CONTINUE 106 CONTINUE C CALL INVMA2(DJI,2,ISING) IF (ISING.EQ.1) THEN WRITE(*,*)'MATRICE DJI singuliere ds BEHAV3' ENDIF C C ************ Debut des iterations internes ******************* C 40 CONTINUE C C *************** Determination de DK et DLAM ****************** C C C DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2) DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2) C IF (DLAM1(1).LE.0.D0) IEC1=1 IF (DLAM1(2).LE.0.D0) IEC2=1 IF (IEC1.EQ.1.OR.IEC2.EQ.1) THEN C WRITE(*,*)'Dans BEHAV3, DLAMDA1 est negatif:',DLAM1(1) C WRITE(*,*)'Dans BEHAV3, DLAMDA2 est negatif:',DLAM1(2) C WRITE(*,*)'A l iteration :',ITER GOTO 15 ENDIF C DK1=DEPS1+DLAM1(1) DK2=DEPS2+DLAM1(2) C IF(IVIS.LE.2) THEN CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF) CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF) ELSE CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1) CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) ENDIF C C ************** Determination de DPHI ********************* C IF (IMOD.EQ.1) THEN AD1=YOUN/(1.D0-XNU) coef=2.D0*SEQ2/(2.D0*SEQ2+DLAM1(2)*AD1) DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1)) TO=SIGF(3) DPHI2=SEQ2 ENDIF IF (IMOD.EQ.2) THEN DG=YOUN/2.D0/(1.D0+XNU) AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU) coef1=(SEQ2+2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG) coef2=(2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG) CAS1=(SIGE(1)+SIGE(2)-0.5*DLAM1(1)*2.D0*AD1) CAS2=SIGE(3)-DLAM1(1)*XNU*AD1 DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2) TO=SIGF(4) DPHI2=SEQ2 ENDIF IF (IMOD.EQ.3) THEN DPHI2=SEQ2-A1/A2*(SIGF(1)+SIGF(2)) AD1=YOUN/(1.D0-XNU) coef=2.D0*DPHI2/(2.D0*DPHI2+DLAM1(2)*AD1) SOM=A1/A2*DLAM1(2)*2.D0*AD1 DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1)-SOM) TO=SIGF(3) ENDIF IF (IMOD.EQ.4) THEN SOM=SIGE(1)+SIGE(2)+SIGE(3) DPHI2=A1*SOM/A2 DPHI2=DPHI2-A1*DLAM1(1)*ADD*(1.D0+XNU)/A2 DPHI2=DPHI2-A1*A1*DLAM1(2)*3.D0*ADD*(1.D0+XNU)/A2/A2 DPHI2=SEQ2-DPHI2 DG=YOUN/2.D0/(1.D0+XNU)/A2/A2 AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU) AD2=YOUN/(1.D0-2.D0*XNU) coef1=(DPHI2+2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG) coef2=(2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG) CAS1=(SIGE(1)+SIGE(2)-DLAM1(1)*AD1-2.D0*A1/A2*DLAM1(2)*AD2) CAS2=SIGE(3)-DLAM1(1)*XNU*AD1-A1/A2*DLAM1(2)*AD2 DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2) TO=SIGF(4) ENDIF C IF (DPHI1.LT.0.D0) THEN C WRITE(*,*)'ATTENTION DPHI1 NEGATIF' C WRITE(*,*)'DPHI1=',DPHI1 C ENDIF C IF (DPHI2.LT.0.D0) THEN C WRITE(*,*)'ATTENTION DPHI2 ds BEHAV3 NEGATIF' C WRITE(*,*)'DPHI2=',DPHI2 C ENDIF C C--------------- Cas de l'apex ----------------------------- C IF (ABS(DPHI1).LE.10E-10.AND. *ABS(DPHI2).GT.10E-10) THEN IAPEX=1 C WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX C WRITE(*,*)'Dans l element',IBB C WRITE(*,*)'et au point d intégration',IGAU DO 107 I=1,NSTRS DO 50 J=1,NSTRS AI(I,J)=0.D0 50 CONTINUE 107 CONTINUE AI(1,1)=0.5 AI(1,2)=0.5 AI(2,1)=AI(1,2) AI(2,2)=AI(1,1) IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0 IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0 GOTO 75 ENDIF C C---------- Cas du critere reduit a un point ---------------- C IF (ABS(DPHI2).LE.10E-10) THEN IAPEX=2 C WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX C WRITE(*,*)'Dans l element',IBB C WRITE(*,*)'et au point d intégration',IGAU DO 108 I=1,NSTRS DO 52 J=1,NSTRS AI(I,J)=0.D0 52 CONTINUE 108 CONTINUE IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN DO 109 I=1,3 DO 53 J=1,3 AI(I,J)=1.D0/3.D0 53 CONTINUE 109 CONTINUE ENDIF GOTO 75 ENDIF C C ************** Mise a jour des contraintes *************** C C ---------------- calcul de la matrice A ------------------ C DO 110 I=1,NSTRS DO 60 J=1,NSTRS A(I,J)=0.D0 60 CONTINUE 110 CONTINUE C DG=YOUN/2.D0/(1.D0+XNU) C IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN DW1=AD*(2.D0-XNU) DW2=AD*(2.D0*XNU-1.D0) A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1 A(1,1)=A(1,1)+(DLAM1(2)*DW1)/2.D0/DPHI2/A2/A2 A(2,2)=A(1,1) A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1 A(1,2)=A(1,2)+(DLAM1(2)*DW2)/2.D0/DPHI2/A2/A2 A(2,1)=A(1,2) A(3,3)=(DLAM1(1)*DG)/DPHI1 A(3,3)=A(3,3)+1.D0+(DLAM1(2)*3.D0*DG)/DPHI2/A2/A2 ENDIF C IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1 A(1,1)=A(1,1)+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2 A(2,2)=A(1,1) A(3,3)=1.D0+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2 A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1 A(1,2)=A(1,2)-(DLAM1(2)*DG)/DPHI2/A2/A2 A(2,1)=A(1,2) A(1,3)=-(DLAM1(2)*DG)/DPHI2/A2/A2 A(3,1)=A(1,3) A(3,2)=A(1,3) A(2,3)=A(1,3) A(4,4)=1.D0+(DLAM1(1)*DG)/DPHI1 A(4,4)=A(4,4)+(DLAM1(2)*DG*3.D0)/DPHI2/A2/A2 ENDIF C C -------------- invertion de la matrice A ----------------- C DO 111 I=1,NSTRS DO 70 J=1,NSTRS AI(I,J)=A(I,J) 70 CONTINUE 111 CONTINUE C CALL INVMA2(AI,NSTRS,ISING) IF (ISING.EQ.1) THEN WRITE(*,*)'MATRICE AI singuliere ds BEHAV3' ENDIF C C -------------- mise a jour des contraintes ------------ C 75 CONTINUE VEC1(1)=AD1 VEC1(2)=AD1 VEC1(3)=0.D0 VEC1(4)=0.D0 IF(IMOD.EQ.3) THEN VEC2(1)=AD1*(1.D0+XNU) VEC2(2)=AD1*(1.D0+XNU) VEC2(3)=0.D0 VEC2(4)=0.D0 ENDIF IF(IMOD.EQ.4) THEN VEC2(1)=AD1*(1.D0+XNU) VEC2(2)=AD1*(1.D0+XNU) VEC2(3)=AD1*(1.D0+XNU) VEC2(4)=0.D0 ENDIF IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN VEC2(1)=0.D0 VEC2(2)=0.D0 VEC2(3)=0.D0 VEC2(4)=0.D0 ENDIF C DO 80 I=1,NSTRS DEPSI(I)=SIGE(I)-0.5*DLAM1(1)*VEC1(I) DEPSI(I)=DEPSI(I)-A1*DLAM1(2)*VEC2(I)/A2 80 CONTINUE C DO 112 I=1,NSTRS SIGF(I)=0.0D+00 DO 90 J=1,NSTRS SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 90 CONTINUE 112 CONTINUE C C ******** Verification des criteres **************** C CALL PRINC(SIGF,V1,NSTRS) IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF) ENDIF IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF) ENDIF FCRI1(1)=V1(1)-SEQ1 IF (IAPEX.EQ.2) FCRI1(1)=0.D0 FCRI1(2)=SEQ-SEQ2 C IF (IBROY.EQ.0.AND.(ABS(FCRI1(1)).GE.CRIMAX.OR * .ABS(FCRI1(2)).GE.CRIMAX)) THEN C WRITE(*,*)'****************************************' C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN' C WRITE(*,*)'on passe donc a la secante' C WRITE(*,*)'Dans l element',IBB C WRITE(*,*)'et au point d intégration',IGAU C WRITE(*,*)'CRIMAX=',CRIMAX C WRITE(*,*)'****************************************' ITER=ITR ENDIF C C ******* Compteur sur la methode de resolution **** C IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN IBROY=1 ITANG=1 ITER=1 IAPEX=0 IEC1=0 IEC2=0 GOTO 12 ENDIF C C ******* non convergence ************************** C IF ((ABS(FCRI1(1)).GT.PRB.OR.ABS(FCRI1(2)).GT.PRB) *.AND.ITER.LT.ITR) THEN IF (IBROY.EQ.0) THEN CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1) DLAM0(1)=DLAM1(1) DLAM0(2)=DLAM1(2) FCRI0(1)=FCRI1(1) FCRI0(2)=FCRI1(2) IF (ITER.GE.(ITR-1)) THEN C WRITE(*,*)'***********************' C WRITE(*,*)'BROYDEN n a pas aboutit' C WRITE(*,*)'ITER=',ITER C WRITE(*,*)'FCRIT=',FCRI0(1) C WRITE(*,*)'FCRIC=',FCRI0(2) C WRITE(*,*)'Dans l element',IBB C WRITE(*,*)'et au point d intégration',IGAU C WRITE(*,*)'***********************' ENDIF ITER=ITER+1 GOTO 40 ENDIF IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN DLAM0(1)=DLAM1(1) DLAM0(2)=DLAM1(2) FCRI0(1)=FCRI1(1) FCRI0(2)=FCRI1(2) CALL PRINC(SIGF,V1,NSTRS) IF (ITER.GE.(ITR-5)) THEN C WRITE(*,*)'ITER=',ITER C WRITE(*,*)'FCRIT=',FCRI0(1) C WRITE(*,*)'FCRIC=',FCRI0(2) ENDIF ITER=ITER+1 GOTO 18 ENDIF ENDIF IF (ITER.GE.ITR.AND.(ABS(FCRI1(1)).GT.PRB2.OR. * ABS(FCRI1(2)).GT.PRB2)) THEN WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3' WRITE(*,*)'Dans l element',IBB WRITE(*,*)'et au point d intégration',IGAU WRITE(*,*)'FCRIT=',FCRI0(1) WRITE(*,*)'FCRIC=',FCRI0(2) WRITE(*,*)'IPLA=',IPLA WRITE(*,*)'IFIS=',IFIS C STOP ENDIF C C **************** Fin des iterations internes ******************* C C C ************* Calcul de la DEP **************** C IF (IAPEX.EQ.0) THEN CALL DERI2(D2FSI1,DPHI1,P1,SIGF,NSTRS,BETJEF) CALL DERI2(D2FSI2,DPHI2,P2,SIGF,NSTRS,BETJEF) CALL LADEP1(SIGF,DEP,PAECT,PAEC,NSTRS,IFOU, & DFSIG1,DFSIG2,D2FSI1,D2FSI2,DLAM1,BETJEF) ELSE IF (IAPEX.EQ.1) THEN DO 114 I=1,NSTRS DO 113 J=1,NSTRS AH(I,J)=0.D0 DO 95 K=1,NSTRS AH(I,J)=AH(I,J)+AI(I,K)*D(K,J) 95 CONTINUE 113 CONTINUE 114 CONTINUE DFSIG1(1)=0.5*SQRT(2.D0) DFSIG1(2)=0.5*SQRT(2.D0) DFSIG1(3)=0.D0 IF (IMOD.EQ.2) DFSIG1(4)=0.D0 CALL CREDEP(AH,DFSIG1,PAEC,NSTRS,DEP,BETJEF) ELSE IF (IAPEX.EQ.2) THEN DO 116 I=1,NSTRS DO 115 J=1,NSTRS DEP(I,J)=0.D0 DO 96 K=1,NSTRS DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J) 96 CONTINUE 115 CONTINUE 116 CONTINUE ENDIF IAPEX=0 C DO 97 I=1,NSTRS C DO 97 J=1,NSTRS C DEP(I,J)=D(I,J) C 97 CONTINUE C C *********************************************** C DEPS1=DK1 DEPS2=DK2 SIG1=SEQ1 SIG2=SEQ2 IF(IVIS.LE.2) THEN CALL INDICA(DK1,DK2,IFIS,IPLA,3,BETJEF) ELSE CALL INDIC1(DK1,DK2,IFIS,IPLA,3,BETJEF,NECH0,NECH1) ENDIF C ********************** CALCUL VISCOPLASTIQUE *********************** IF (IVIS.EQ.1) THEN CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2, & SIGP,SIGRV,DSTRN,D,BETJEF,VISCO) ENDIF C ******************************************************************** 100 CONTINUE IEC1=0 IEC2=0 CRIMAX=0.D0 C RETURN END