1 : C ACTI1 SOURCE CHAT 05/01/12 21:16:17 5004 2 : C ACTI1 SOURCE 3 : SUBROUTINE ACTI1(SIG0,SIGF,D,NSTRSS,BETINSA) 4 : C 5 : C =============================================================== 6 : C Un seul critere de traction: DRUCKER PRAGER 3D 7 : C ================================================================ 8 : C CE SOUS-PROGRAMME EST APPELE DANS "BONE3D". 9 : C 10 : IMPLICIT INTEGER(I-N) 11 : IMPLICIT REAL*8(A-H,O-Z) 12 : DIMENSION SIG0(NSTRSS),SIGF(NSTRSS) 13 : DIMENSION DFSIG(6),VEC1(6),DEPSI(6),SIGE(6) 14 : DIMENSION D(6,6),AC(6),A(6,6),AI(6,6),AIM(4,4) 15 : C 16 : SEGMENT BETINSA 17 : REAL*8 RT,RC,YOUN,XNU,GFT,GFC,CAR 18 : REAL*8 DKT,DKC,SEQT,SEQC,ENDT,ENDC 19 : INTEGER IFIS,IPLA,IBB,IGAU 20 : ENDSEGMENT 21 : C 21 : CALL ZERO(SIGE,6,1) 22 : DO 2 I=1,NSTRSS 23 : SIGE(I)=SIGF(I) 24 : 2 CONTINUE 25 : IAPEX=0 26 : PRB=1.D-5 27 : PRB2=1.D-2 28 : ITER=1 29 : ITR = 1500 30 : ITANG=0 31 : IBROY=0 32 : CRIMAX=0.D0 33 : SEQ = 0.D0 34 : SEQ1 = 0.D0 35 : C 36 : 8 CONTINUE 37 : C 38 : DO 4 I=1,NSTRSS 39 : SIGF(I)=SIGE(I) 40 : 4 CONTINUE 41 : DK=DKT 42 : DLAM0=0.D0 42 : CALL DRUTRA(SIGF,SEQTE,BETINSA) 43 : FCRI0 = SEQTE - SEQT 44 : CRIMAX=DABS(100.D0*FCRI0) 45 : 9 CONTINUE 46 : C 47 : C ************ CALCUL DU JACOBIEN INITIAL ******************** 48 : C 48 : CALL DRUTR1(SIGF,DFSIG,BETINSA) 49 : C 50 : DO 25 I=1,NSTRSS 51 : AC(I)=0.D0 52 : DO 25 J=1,NSTRSS 53 : AC(I)=AC(I)+D(I,J)*DFSIG(J) 54 : 25 CONTINUE 55 : FDF=0.D0 56 : DO 30 J=1,NSTRSS 57 : FDF=FDF+DFSIG(J)*AC(J) 58 : 30 CONTINUE 59 : C 59 : CALL ENDAME(1,BETINSA) 59 : CALL FORECR(DK,PAEC,1,SEQ,BETINSA) 60 : C 61 : DJAC0=-(PAEC+FDF) 62 : C 63 : C ************ DEBUT ITERATION INTERNES ************************ 64 : C 65 : 40 CONTINUE 66 : C 67 : C *************** Determination de DK ************************* 68 : C 69 : DLAM1=-FCRI0/DJAC0+DLAM0 70 : DK=DKT+DLAM1 71 : C 72 : C *********** Estimation contrainte quivalente ***************** 73 : C 73 : CALL ENDAME(1,BETINSA) 73 : CALL FORECR(DK,PAEC,1,SEQ,BETINSA) 74 : C 75 : C ************** Determination de DPHI ************************* 76 : C 76 : CALL DRUTR2(SIGE,SEQ,DPHI,DLAM1,VEC1,BETINSA) 77 : C 78 : C **************** Cas de l'apex ******************************* 79 : C 80 : IF (DABS(DPHI).LE.10E-10) THEN 81 : IAPEX=1 82 : DO 50 I=1,NSTRSS 83 : DO 50 J=1,NSTRSS 84 : AI(I,J)=0.D0 85 : 50 CONTINUE 86 : AI(1,1)=1./3. 87 : AI(1,2)=AI(1,1) 88 : AI(1,2)=AI(1,1) 89 : AI(1,3)=AI(1,1) 90 : AI(2,1)=AI(1,1) 91 : AI(2,2)=AI(1,1) 92 : AI(2,3)=AI(1,1) 93 : AI(3,1)=AI(1,1) 94 : AI(3,2)=AI(1,1) 95 : AI(3,3)=AI(1,1) 96 : GOTO 75 97 : ENDIF 98 : C 99 : C ************** Mise a jour des contraintes *************** 100 : C 101 : C ---------------- calcul de la matrice A ------------------ 102 : C 103 : DO 60 I=1,NSTRSS 104 : DO 60 J=1,NSTRSS 105 : A(I,J)=0.D0 106 : 60 CONTINUE 107 : C 108 : DG=YOUN/(1.D0+XNU) 109 : C 110 : A(1,1)=1.D0+2.*(DLAM1*DG)/2.D0/DPHI 111 : A(2,2)=A(1,1) 112 : A(3,3)=A(1,1) 113 : A(1,2)=-(DLAM1*DG)/2.D0/DPHI 114 : A(1,3)=A(1,2) 115 : A(2,1)=A(1,2) 116 : A(2,3)=A(1,2) 117 : A(3,1)=A(1,2) 118 : A(3,2)=A(1,2) 119 : A(4,4)=1.D0+3.*(DLAM1*DG)/2.D0/DPHI 120 : A(5,5)=A(4,4) 121 : A(6,6)=A(4,4) 122 : C 123 : C -------------- invertion de la matrice A ----------------- 124 : C 124 : CALL ZERO(AI,6,6) 124 : CALL ZERO(AIM,6,6) 125 : C 126 : DO 70 I=1,3 127 : DO 70 J=1,3 128 : AIM(I,J)=A(I,J) 129 : 70 CONTINUE 129 : CALL INVMA2(AIM,3,ISING) 130 : IF (ISING.EQ.1) THEN 131 : WRITE(*,*)'MATRICE AIM singuliere ds ACTI1' 132 : ENDIF 133 : DO 72 I=1,3 134 : DO 72 J=1,3 135 : AI(I,J)=AIM(I,J) 136 : 72 CONTINUE 137 : AI(4,4) = 1./A(4,4) 138 : AI(5,5) = 1./A(5,5) 139 : AI(6,6) = 1./A(6,6) 140 : C 141 : C -------------- mise a jour des contraintes ------------ 142 : C 143 : 75 CONTINUE 144 : C 145 : DO 80 I=1,NSTRSS 146 : DEPSI(I)=SIGE(I)-DLAM1*VEC1(I) 147 : 80 CONTINUE 148 : C 149 : DO 90 I=1,NSTRSS 150 : SIGF(I)=0. 151 : DO 90 J=1,NSTRSS 152 : SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 153 : 90 CONTINUE 154 : C 155 : C ******** VERIFICATION DU CRITERE **************** 156 : C 156 : CALL DRUTRA(SIGF,SEQI,BETINSA) 157 : FCRI1 = SEQI - SEQ 158 : C 159 : IF (IBROY.EQ.0.AND.DABS(FCRI1).GE.CRIMAX) THEN 160 : C WRITE(*,*)'****************************************' 161 : C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN' 162 : C WRITE(*,*)'on passe donc a la secante' 163 : C WRITE(*,*)'Dans l element',IBB 164 : C WRITE(*,*)'et au point d intégration',IGAU 165 : C WRITE(*,*)'CRIMAX=',CRIMAX 166 : C WRITE(*,*)'****************************************' 167 : ITER=ITR 168 : ENDIF 169 : C 170 : C ******* Compteur sur la methode de resolution **** 171 : C 172 : IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN 173 : IBROY=1 174 : ITANG=1 175 : ITER=1 176 : IAPEX=0 177 : GOTO 8 178 : ENDIF 179 : C 180 : C ******* non convergence ************************** 181 : C 182 : IF (DABS(FCRI1).GT.PRB.AND.ITER.LT.ITR) THEN 183 : IF (IBROY.EQ.0) THEN 184 : DJAC1=(FCRI0-FCRI1)/(DLAM0-DLAM1) 185 : DLAM0=DLAM1 186 : DJAC0=DJAC1 187 : FCRI0=FCRI1 188 : ITER=ITER+1 189 : GOTO 40 190 : ENDIF 191 : IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN 192 : DLAM0=DLAM1 193 : FCRI0=FCRI1 194 : ITER=ITER+1 195 : IF (ITER.GE.(ITR-5)) THEN 196 : WRITE(*,*)'ITER=',ITER 197 : WRITE(*,*)'FCRI=',FCRI1 198 : ENDIF 199 : GOTO 9 200 : ENDIF 201 : ENDIF 202 : IF (ITER.GE.ITR.AND.DABS(FCRI1).GT.PRB2) THEN 203 : WRITE(*,*)'NON CONVERGENCE INTERNE dans ACTI1' 204 : WRITE(*,*)'FCRI=',FCRI1 205 : WRITE(*,*)'Dans l element',IBB 206 : WRITE(*,*)'et au point d intégration',IGAU 207 : STOP 208 : ENDIF 209 : C 210 : C ******* convergence ************************** 211 : C 212 : SEQT = SEQ 213 : DKT = DK 214 : C 215 : C *********************************************** 216 : C 217 : 100 CONTINUE 218 : C 219 : RETURN 220 : END 221 : 222 : 223 : 224 : 225 :
© Cast3M 2003 - All rights reserved.
Disclaimer