C ACTI3 SOURCE CB215821 16/04/21 21:15:04 8920 C ACTI3 SUBROUTINE ACTI3(SIG0,SIGF,D,NSTRSS,BETINSA) C C ================================================================== C Deux criteres de traction compression: DRUCKER PRAGER C 3D C ================================================================== C CE SOUS-PROGRAMME EST APPELE DANS "BONE3D". C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION SIG0(NSTRSS),SIGF(NSTRSS) DIMENSION DFSIG1(6),DFSIG2(6),VEC1(6),VEC2(6),SIGE(6) DIMENSION AC1(6),AC2(6) DIMENSION DJAC0(2,2),DJI(4,4),A(6,6),AI(6,6),AIM(6,6),D(6,6) DIMENSION FCRI0(2),FCRI1(2),DEPSI(6),DLAM1(2),DLAM0(2) C SEGMENT BETINSA REAL*8 RT,RC,YOUN,XNU,GFT,GFC,CAR REAL*8 DKT,DKC,SEQT,SEQC,ENDT,ENDC INTEGER IFIS,IPLA,IBB,IGAU ENDSEGMENT C IAPEX=0 PRB=1.D-5 PRB2=1.D-2 ITER=1 ITR = 1500 IET=0 IEC=0 IBROY=0 ITANG=0 CRIMAX=0.D0 SEQ = 0.D0 SEQ1 = 0.D0 SEQ2 = 0.D0 CALL ZERO(SIGE,6,1) C DO 10 I=1,NSTRSS SIGE(I)=SIGF(I) 10 CONTINUE C C ************ Le point est fissure pour 1ere fois *********** C 12 CONTINUE CALL DRUTRA(SIGF,SEQTE,BETINSA) CALL DRUCOM(SIGF,SEQCE,BETINSA) FCRI0(1) = SEQTE - SEQT FCRI0(2) = SEQCE - SEQC CRIMAX=ABS(100.D0*FCRI0(1)) IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN IET=1 IEC=1 ENDIF IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN IET=1 IEC=0 ENDIF IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN IET=0 IEC=1 ENDIF C 15 CONTINUE C DO 16 I=1,NSTRSS SIGF(I)=SIGE(I) 16 CONTINUE C IF (IET.EQ.1.AND.IEC.EQ.1) THEN GOTO 100 ENDIF C IF (IET.EQ.1.AND.IEC.EQ.0) THEN CALL ACTI2(SIG0,SIGF,D,NSTRSS,BETINSA) GOTO 100 ENDIF C IF (IET.EQ.0.AND.IEC.EQ.1) THEN CALL ACTI1(SIG0,SIGF,D,NSTRSS,BETINSA) GOTO 100 ENDIF C DK1=DKT DK2=DKC DLAM0(1)=0.D0 DLAM0(2)=0.D0 C 18 CONTINUE C C ************ CALCUL DU JACOBIEN INITIAL ******************** C C ---------------- direction de traction ------------------ C CALL DRUTR1(SIGF,DFSIG1,BETINSA) C DO 20 I=1,NSTRSS AC1(I)=0.D0 DO 20 J=1,NSTRSS AC1(I)=AC1(I)+D(I,J)*DFSIG1(J) 20 CONTINUE C C ---------------- direction de compression ---------------- C CALL DRUCO1(SIGF,DFSIG2,BETINSA) C DO 25 I=1,NSTRSS AC2(I)=0.D0 DO 25 J=1,NSTRSS AC2(I)=AC2(I)+D(I,J)*DFSIG2(J) 25 CONTINUE C------------------------------------------------------------ F1DF1=0.D0 DO 30 J=1,NSTRSS F1DF1=F1DF1+DFSIG1(J)*AC1(J) 30 CONTINUE F1DF2=0.D0 DO 31 J=1,NSTRSS F1DF2=F1DF2+DFSIG1(J)*AC2(J) 31 CONTINUE F2DF2=0.D0 DO 32 J=1,NSTRSS F2DF2=F2DF2+DFSIG2(J)*AC2(J) 32 CONTINUE F2DF1=0.D0 DO 33 J=1,NSTRSS F2DF1=F2DF1+DFSIG2(J)*AC1(J) 33 CONTINUE C CALL ENDAME(1,BETINSA) CALL FORECR(DK1,PAECT,1,SEQ,BETINSA) CALL ENDAME(2,BETINSA) CALL FORECR(DK2,PAECC,2,SEQ,BETINSA) C DJAC0(1,1)=-(F1DF1+PAECT) DJAC0(1,2)=-(F1DF2) DJAC0(2,2)=-(F2DF2+PAECC) DJAC0(2,1)=-(F2DF1) C DO 43 I=1,2 DO 43 J=1,2 DJI(I,J)=DJAC0(I,J) 43 CONTINUE C CALL INVMA2(DJI,2,ISING) IF (ISING.EQ.1) THEN WRITE(*,*)'MATRICE DJI singuliere ds ACTI3' ENDIF C C ************ DEBUT ITERATION 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) IET=1 IF (DLAM1(2).LE.0.D0) IEC=1 IF (IET.EQ.1.OR.IEC.EQ.1) THEN C WRITE(*,*)'Dans ACTI3, DLAMDA1 est negatif:',DLAM1(1) C WRITE(*,*)'Dans ACTI3, DLAMDA2 est negatif:',DLAM1(2) C WRITE(*,*)'A l iteration :',ITER GOTO 15 ENDIF C DK1=DKT+DLAM1(1) DK2=DKC+DLAM1(2) C CALL ENDAME(1,BETINSA) CALL FORECR(DK1,PAECT,1,SEQ1,BETINSA) CALL ENDAME(2,BETINSA) CALL FORECR(DK2,PAECC,2,SEQ2,BETINSA) C C ************** Determination de DPHI1 et 2 ****************** C CALL DRUTR2(SIGE,SEQ1,DPHI1,DLAM1(1),VEC1,BETINSA) CALL DRUCO2(SIGE,SEQ2,DPHI2,DLAM1(2),VEC2,BETINSA) 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 ACTI3 =',IAPEX C WRITE(*,*)'Dans l element',IBB C WRITE(*,*)'et au point d intégration',IGAU DO 50 I=1,NSTRSS DO 50 J=1,NSTRSS AI(I,J)=0.D0 50 CONTINUE AI(1,1)=1./3. AI(1,2)=AI(1,1) AI(1,2)=AI(1,1) AI(1,3)=AI(1,1) AI(2,1)=AI(1,1) AI(2,2)=AI(1,1) AI(2,3)=AI(1,1) AI(3,1)=AI(1,1) AI(3,2)=AI(1,1) AI(3,3)=AI(1,1) 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 ACTI3 =',IAPEX C WRITE(*,*)'Dans l element',IBB C WRITE(*,*)'et au point d intégration',IGAU DO 52 I=1,NSTRSS DO 52 J=1,NSTRSS AI(I,J)=0.D0 52 CONTINUE AI(1,1)=1./3. AI(1,2)=AI(1,1) AI(1,2)=AI(1,1) AI(1,3)=AI(1,1) AI(2,1)=AI(1,1) AI(2,2)=AI(1,1) AI(2,3)=AI(1,1) AI(3,1)=AI(1,1) AI(3,2)=AI(1,1) AI(3,3)=AI(1,1) GOTO 75 ENDIF C C ************** Mise a jour des contraintes *************** C C ---------------- calcul de la matrice A ------------------ C DO 60 I=1,NSTRSS DO 60 J=1,NSTRSS A(I,J)=0.D0 60 CONTINUE C DG=YOUN/(1.D0+XNU) C A(1,1)=1.D0+2.*(DLAM1(1)*DG)/2.D0/DPHI1 A(1,1)=A(1,1)+2.*(DLAM1(2)*DG)/2.D0/DPHI2 A(2,2)=A(1,1) A(3,3)=A(1,1) A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1-(DLAM1(2)*DG)/2.D0/DPHI2 A(1,3)=A(1,2) A(2,1)=A(1,2) A(2,3)=A(1,2) A(3,1)=A(1,2) A(3,2)=A(1,2) A(4,4)=1.D0+3.*(DLAM1(1)*DG)/2.D0/DPHI1 A(4,4)=A(4,4)+3.*(DLAM1(2)*DG)/2.D0/DPHI2 A(5,5)=A(4,4) A(6,6)=A(4,4) C C -------------- invertion de la matrice A ----------------- C CALL ZERO(AI,6,6) CALL ZERO(AIM,6,6) C DO 70 I=1,3 DO 70 J=1,3 AIM(I,J)=A(I,J) 70 CONTINUE CALL INVMA2(AIM,3,ISING) IF (ISING.EQ.1) THEN WRITE(*,*)'MATRICE AIM singuliere ds ACTI3' ENDIF DO 72 I=1,3 DO 72 J=1,3 AI(I,J)=AIM(I,J) 72 CONTINUE AI(4,4) = 1./A(4,4) AI(5,5) = 1./A(5,5) AI(6,6) = 1./A(6,6) C C -------------- mise a jour des contraintes ------------ C 75 CONTINUE C DO 80 I=1,NSTRSS DEPSI(I)=SIGE(I)-DLAM1(1)*VEC1(I) DEPSI(I)=DEPSI(I)-DLAM1(2)*VEC2(I) 80 CONTINUE C DO 90 I=1,NSTRSS SIGF(I)=0.0D+00 DO 90 J=1,NSTRSS SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 90 CONTINUE C C ******** Verification des criteres **************** C CALL DRUTRA(SIGF,SEQTT,BETINSA) FCRI1(1) = SEQTT - SEQ1 CALL DRUCOM(SIGF,SEQCC,BETINSA) FCRI1(2) = SEQCC - SEQ2 C IF (IAPEX.EQ.2) FCRI1(1)=0.D0 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 IET=0 IEC=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) 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 DKT=DK1 DKC=DK2 SEQT=SEQ1 SEQC=SEQ2 C C ******************************************************************** 100 CONTINUE IET=0 IEC=0 CRIMAX=0.D0 C RETURN END