C BEHAV2 SOURCE CB215821 16/04/21 21:15:18 8920 SUBROUTINE BEHAV2(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 Un seul critere de compression: Von Mises ou Drucker-Prager C C ================================================================== C CE SOUS-PROGRAMME EST APPELE DANS "BONE". C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION DFSIG(4),SIGF(4),DSIGT(4),VEC1(4),DEPSI(4),DSTRN(4) DIMENSION SIGR(4),AC(4),SIGE(4),SIGRV(4),SIGP(4) DIMENSION D(4,4),DEP(4,4),D2FSIG(4,4),P2(4,4),DP(4,4) DIMENSION A(4,4),AI(4,4),AH(4,4) * 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 * * C * 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 IAPEX=0 PRB=1.D-10 PRB2=1.D-6 ITER=1 IBROY=0 ITANG=0 Rb=COLI Ft=PALF*COLI IEC2=0 CRIMAX=0.D0 SEQ = 0.D0 SEQ2 = 0.D0 CALL ZERO(SIGE,4,1) ***** CALL ZERO(V1,4,1) 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 comprime pour 1ere fois *********** C 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 IF (IFIS.EQ.0) THEN DEPS1=0.D0 SIG1=Ft ENDIF C 8 CONTINUE C 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 FCRIC0=SEQ-SIG2 CRIMAX=ABS(100.D0*FCRIC0) IF (FCRIC0.LE.0.D0) THEN IEC2=1 ENDIF DK=DEPS2 DLAM0=0.D0 C 12 CONTINUE C DO 15 I=1,NSTRS SIGF(I)=SIGE(I) 15 CONTINUE C IF (IEC2.EQ.1) THEN DO 16 I=1,NSTRS DO 16 J=1,NSTRS DEP(I,J)=D(I,J) 16 CONTINUE GOTO 100 ENDIF C 9 CONTINUE C IF(IVIS.LE.2) THEN CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF) ELSE CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) ENDIF 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 ---------------------------------------------------------- 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 BEHAV2' STOP ENDIF DFSIG(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2 DFSIG(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2 DFSIG(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 BEHAV2' ENDIF DFSIG(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2 DFSIG(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2 DFSIG(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2 DFSIG(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 25 I=1,NSTRS AC(I)=0.D0 DO 25 J=1,NSTRS AC(I)=AC(I)+D(I,J)*DFSIG(J) 25 CONTINUE C C ************ Preparation du jacobien *********************** C FDF=0.D0 DO 30 J=1,NSTRS FDF=FDF+DFSIG(J)*AC(J) 30 CONTINUE C DJAC0=-(FDF+PAEC) C C ************ Debut des iterations internes ******************* C 40 CONTINUE C C *************** Determination de DK et DLAM ****************** C C DLAM1=DLAM0-FCRIC0/DJAC0 C IF (DLAM1.LT.0.D0) IEC2=1 IF (IEC2.EQ.1) THEN C WRITE(*,*)'Dans BEHAV2, DLAMDA est negatif:',DLAM1 C WRITE(*,*)'A l iteration :',ITER GOTO 12 ENDIF C DK=DEPS2+DLAM1 C IF(IVIS.LE.2) THEN CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF) ELSE CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) ENDIF C C ************** Determination de DPHI2 ********************* C IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN DPHI2=SEQ2 ENDIF IF (IMOD.EQ.3) THEN SOM=SIGF(1)+SIGF(2)+SIGF(3) DPHI2=SEQ2-A1*SOM/A2 ENDIF IF (IMOD.EQ.4) THEN SOM=SIGE(1)+SIGE(2)+SIGE(3) DPHI2=A1*SOM/A2-A1*A1*DLAM1*3.D0*ADD*(1.D0+XNU)/A2/A2 DPHI2=SEQ2-DPHI2 ENDIF C IF (DPHI2.LT.0.D0) THEN C WRITE(*,*)'ATTENTION DPHI2 ds BEHAV2 NEGATIF' C WRITE(*,*)'DPHI2=',DPHI2 C ENDIF C C---------- Cas du critere reduit a un point ---------------- C IF (ABS(DPHI2).LE.10E-10) THEN IAPEX=1 C WRITE(*,*)'IAPEX ds BEHAV2 =',IAPEX DO 50 I=1,NSTRS DO 50 J=1,NSTRS AI(I,J)=0.D0 50 CONTINUE IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN DO 53 I=1,3 DO 53 J=1,3 AI(I,J)=1.D0/3.D0 53 CONTINUE ENDIF GOTO 75 ENDIF C C ************** Mise a jour des contraintes *************** C C ---------------- calcul de la matrice A ------------------ C DO 60 I=1,NSTRS DO 60 J=1,NSTRS A(I,J)=0.D0 60 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*DW1)/2.D0/DPHI2/A2/A2 A(2,2)=A(1,1) A(1,2)=(DLAM1*DW2)/2.D0/DPHI2/A2/A2 A(2,1)=A(1,2) A(3,3)=1.D0+(DLAM1*3.D0*DG)/DPHI2/A2/A2 ENDIF IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN A(1,1)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2 A(2,2)=A(1,1) A(3,3)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2 A(1,2)=-(DLAM1*DG)/DPHI2/A2/A2 A(2,1)=A(1,2) A(1,3)=-(DLAM1*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*DG*3.D0)/DPHI2/A2/A2 ENDIF C C -------------- invertion de la matrice A ----------------- C DO 70 I=1,NSTRS DO 70 J=1,NSTRS AI(I,J)=A(I,J) 70 CONTINUE CALL INVMA2(AI,NSTRS,ISING) IF (ISING.EQ.1) THEN WRITE(*,*)'MATRICE AI singuliere ds BEHAV2' WRITE(*,*)'NSTRS=',NSTRS ENDIF C C -------------- mise a jour des contraintes -------------- C 75 CONTINUE IF(IMOD.EQ.3) THEN VEC1(1)=AD*(1.D0+XNU) VEC1(2)=AD*(1.D0+XNU) VEC1(3)=0.D0 VEC1(4)=0.D0 ENDIF IF(IMOD.EQ.4) THEN VEC1(1)=ADD*(1.D0+XNU) VEC1(2)=ADD*(1.D0+XNU) VEC1(3)=ADD*(1.D0+XNU)*2.D0*XNU VEC1(4)=0.D0 ENDIF IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN VEC1(1)=0.D0 VEC1(2)=0.D0 VEC1(3)=0.D0 VEC1(4)=0.D0 ENDIF C DO 80 I=1,NSTRS DEPSI(I)=SIGE(I)-A1*DLAM1*VEC1(I)/A2 80 CONTINUE C DO 90 I=1,NSTRS SIGF(I)=0.0D+00 DO 90 J=1,NSTRS SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 90 CONTINUE C C ************ Verification du critere ******************** C 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 FCRIC1=SEQ-SEQ2 C IF (IBROY.EQ.0.AND.ABS(FCRIC1).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 IAPEX=0 IEC2=0 ITER=1 GOTO 8 ENDIF C C ******* non convergence ************************** C IF (ABS(FCRIC1).GT.PRB.AND.ITER.LT.ITR) THEN IF (IBROY.EQ.0) THEN DJAC0=(FCRIC0-FCRIC1)/(DLAM0-DLAM1) DLAM0=DLAM1 FCRIC0=FCRIC1 ITER=ITER+1 IF (ITER.GE.(ITR-1)) THEN C WRITE(*,*)'***********************' C WRITE(*,*)'BROYDEN n a pas aboutit' C WRITE(*,*)'Dans l element',IBB C WRITE(*,*)'et au point d intégration',IGAU C WRITE(*,*)'ITER=',ITER C WRITE(*,*)'FCRI=',FCRIC1 C WRITE(*,*)'***********************' ENDIF GOTO 40 ENDIF IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN DLAM0=DLAM1 FCRIC0=FCRIC1 ITER=ITER+1 IF (ITER.GE.(ITR-5)) THEN C WRITE(*,*)'ITER=',ITER C WRITE(*,*)'FCRI=',FCRIC1 ENDIF GOTO 9 ENDIF ENDIF IF (ITER.GE.ITR.AND.ABS(FCRIC1).GT.PRB2) THEN WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV2' WRITE(*,*)'FCRI=',FCRIC1 WRITE(*,*)'Dans l element',IBB WRITE(*,*)'et au point d intégration',IGAU WRITE(*,*)'IPLA=',IPLA WRITE(*,*)'IFIS=',IFIS C STOP ENDIF C C ************* Calcul de la DEP **************** C IF (IAPEX.EQ.0) THEN CALL DERI2(D2FSIG,DPHI2,P2,SIGF,NSTRS,BETJEF) CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU, & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF) ELSE IF (IAPEX.EQ.1) THEN C WRITE(*,*)'ELAS2',IAPEX DO 96 I=1,NSTRS DO 96 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 ENDIF IAPEX=0 CRIMAX=0.D0 * DO 97 I=1,NSTRS * DO 97 J=1,NSTRS * DEP(I,J)=D(I,J) * 97 CONTINUE C C *********************************************** C DEPS2=DK SIG2=SEQ2 IF(IVIS.LE.2) THEN CALL INDICA(0.D0,DK,IFIS,IPLA,2,BETJEF) ELSE CALL INDIC1(0.D0,DK,IFIS,IPLA,2,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 IEC2=0 C RETURN END