1 : C ACTI3 SOURCE PV 08/09/11 21:15:01 6150 2 : C ACTI3 3 : SUBROUTINE ACTI3(SIG0,SIGF,D,NSTRSS,BETINSA) 4 : C 5 : C ================================================================== 6 : C Deux criteres de traction compression: DRUCKER PRAGER 7 : C 3D 8 : C ================================================================== 9 : C CE SOUS-PROGRAMME EST APPELE DANS "BONE3D". 10 : C 11 : IMPLICIT INTEGER(I-N) 12 : IMPLICIT REAL*8(A-H,O-Z) 13 : DIMENSION SIG0(NSTRSS),SIGF(NSTRSS) 14 : DIMENSION DFSIG1(6),DFSIG2(6),VEC1(6),VEC2(6),SIGE(6) 15 : DIMENSION AC1(6),AC2(6) 16 : DIMENSION DJAC0(2,2),DJI(4,4),A(6,6),AI(6,6),AIM(6,6),D(6,6) 17 : DIMENSION FCRI0(2),FCRI1(2),DEPSI(6),DLAM1(2),DLAM0(2) 18 : C 19 : SEGMENT BETINSA 20 : REAL*8 RT,RC,YOUN,XNU,GFT,GFC,CAR 21 : REAL*8 DKT,DKC,SEQT,SEQC,ENDT,ENDC 22 : INTEGER IFIS,IPLA,IBB,IGAU 23 : ENDSEGMENT 24 : C 25 : IAPEX=0 26 : PRB=1.D-5 27 : PRB2=1.D-2 28 : ITER=1 29 : ITR = 1500 30 : IET=0 31 : IEC=0 32 : IBROY=0 33 : ITANG=0 34 : CRIMAX=0.D0 35 : SEQ = 0.D0 36 : SEQ1 = 0.D0 37 : SEQ2 = 0.D0 37 : CALL ZERO(SIGE,6,1) 38 : C 39 : DO 10 I=1,NSTRSS 40 : SIGE(I)=SIGF(I) 41 : 10 CONTINUE 42 : C 43 : C ************ Le point est fissure pour 1ere fois *********** 44 : C 45 : 12 CONTINUE 45 : CALL DRUTRA(SIGF,SEQTE,BETINSA) 45 : CALL DRUCOM(SIGF,SEQCE,BETINSA) 46 : FCRI0(1) = SEQTE - SEQT 47 : FCRI0(2) = SEQCE - SEQC 48 : CRIMAX=DABS(100.D0*FCRI0(1)) 49 : IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN 50 : IET=1 51 : IEC=1 52 : ENDIF 53 : IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN 54 : IET=1 55 : IEC=0 56 : ENDIF 57 : IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN 58 : IET=0 59 : IEC=1 60 : ENDIF 61 : C 62 : 15 CONTINUE 63 : C 64 : DO 16 I=1,NSTRSS 65 : SIGF(I)=SIGE(I) 66 : 16 CONTINUE 67 : C 68 : IF (IET.EQ.1.AND.IEC.EQ.1) THEN 69 : GOTO 100 70 : ENDIF 71 : C 72 : IF (IET.EQ.1.AND.IEC.EQ.0) THEN 72 : CALL ACTI2(SIG0,SIGF,D,NSTRSS,BETINSA) 73 : GOTO 100 74 : ENDIF 75 : C 76 : IF (IET.EQ.0.AND.IEC.EQ.1) THEN 76 : CALL ACTI1(SIG0,SIGF,D,NSTRSS,BETINSA) 77 : GOTO 100 78 : ENDIF 79 : C 80 : DK1=DKT 81 : DK2=DKC 82 : DLAM0(1)=0.D0 83 : DLAM0(2)=0.D0 84 : C 85 : 18 CONTINUE 86 : C 87 : C ************ CALCUL DU JACOBIEN INITIAL ******************** 88 : C 89 : C ---------------- direction de traction ------------------ 90 : C 90 : CALL DRUTR1(SIGF,DFSIG1,BETINSA) 91 : C 92 : DO 20 I=1,NSTRSS 93 : AC1(I)=0.D0 94 : DO 20 J=1,NSTRSS 95 : AC1(I)=AC1(I)+D(I,J)*DFSIG1(J) 96 : 20 CONTINUE 97 : C 98 : C ---------------- direction de compression ---------------- 99 : C 99 : CALL DRUCO1(SIGF,DFSIG2,BETINSA) 100 : C 101 : DO 25 I=1,NSTRSS 102 : AC2(I)=0.D0 103 : DO 25 J=1,NSTRSS 104 : AC2(I)=AC2(I)+D(I,J)*DFSIG2(J) 105 : 25 CONTINUE 106 : C------------------------------------------------------------ 107 : F1DF1=0.D0 108 : DO 30 J=1,NSTRSS 109 : F1DF1=F1DF1+DFSIG1(J)*AC1(J) 110 : 30 CONTINUE 111 : F1DF2=0.D0 112 : DO 31 J=1,NSTRSS 113 : F1DF2=F1DF2+DFSIG1(J)*AC2(J) 114 : 31 CONTINUE 115 : F2DF2=0.D0 116 : DO 32 J=1,NSTRSS 117 : F2DF2=F2DF2+DFSIG2(J)*AC2(J) 118 : 32 CONTINUE 119 : F2DF1=0.D0 120 : DO 33 J=1,NSTRSS 121 : F2DF1=F2DF1+DFSIG2(J)*AC1(J) 122 : 33 CONTINUE 123 : C 123 : CALL ENDAME(1,BETINSA) 123 : CALL FORECR(DK1,PAECT,1,SEQ,BETINSA) 123 : CALL ENDAME(2,BETINSA) 123 : CALL FORECR(DK2,PAECC,2,SEQ,BETINSA) 124 : C 125 : DJAC0(1,1)=-(F1DF1+PAECT) 126 : DJAC0(1,2)=-(F1DF2) 127 : DJAC0(2,2)=-(F2DF2+PAECC) 128 : DJAC0(2,1)=-(F2DF1) 129 : C 130 : DO 43 I=1,2 131 : DO 43 J=1,2 132 : DJI(I,J)=DJAC0(I,J) 133 : 43 CONTINUE 134 : C 134 : CALL INVMA2(DJI,2,ISING) 135 : IF (ISING.EQ.1) THEN 136 : WRITE(*,*)'MATRICE DJI singuliere ds ACTI3' 137 : ENDIF 138 : C 139 : C ************ DEBUT ITERATION INTERNES ******************* 140 : C 141 : 40 CONTINUE 142 : C 143 : C *************** Determination de DK et DLAM ****************** 144 : C 145 : C 146 : C 147 : DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2) 148 : DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2) 149 : C 150 : IF (DLAM1(1).LE.0.D0) IET=1 151 : IF (DLAM1(2).LE.0.D0) IEC=1 152 : IF (IET.EQ.1.OR.IEC.EQ.1) THEN 153 : C WRITE(*,*)'Dans ACTI3, DLAMDA1 est negatif:',DLAM1(1) 154 : C WRITE(*,*)'Dans ACTI3, DLAMDA2 est negatif:',DLAM1(2) 155 : C WRITE(*,*)'A l iteration :',ITER 156 : GOTO 15 157 : ENDIF 158 : C 159 : DK1=DKT+DLAM1(1) 160 : DK2=DKC+DLAM1(2) 161 : C 161 : CALL ENDAME(1,BETINSA) 161 : CALL FORECR(DK1,PAECT,1,SEQ1,BETINSA) 161 : CALL ENDAME(2,BETINSA) 161 : CALL FORECR(DK2,PAECC,2,SEQ2,BETINSA) 162 : C 163 : C ************** Determination de DPHI1 et 2 ****************** 164 : C 164 : CALL DRUTR2(SIGE,SEQ1,DPHI1,DLAM1(1),VEC1,BETINSA) 164 : CALL DRUCO2(SIGE,SEQ2,DPHI2,DLAM1(2),VEC2,BETINSA) 165 : C 166 : C--------------- Cas de l'apex ----------------------------- 167 : C 168 : IF (DABS(DPHI1).LE.10E-10.AND. 169 : *DABS(DPHI2).GT.10E-10) THEN 170 : IAPEX=1 171 : C WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX 172 : C WRITE(*,*)'Dans l element',IBB 173 : C WRITE(*,*)'et au point d intégration',IGAU 174 : DO 50 I=1,NSTRSS 175 : DO 50 J=1,NSTRSS 176 : AI(I,J)=0.D0 177 : 50 CONTINUE 178 : AI(1,1)=1./3. 179 : AI(1,2)=AI(1,1) 180 : AI(1,2)=AI(1,1) 181 : AI(1,3)=AI(1,1) 182 : AI(2,1)=AI(1,1) 183 : AI(2,2)=AI(1,1) 184 : AI(2,3)=AI(1,1) 185 : AI(3,1)=AI(1,1) 186 : AI(3,2)=AI(1,1) 187 : AI(3,3)=AI(1,1) 188 : GOTO 75 189 : ENDIF 190 : C 191 : C---------- Cas du critere reduit a un point ---------------- 192 : C 193 : IF (DABS(DPHI2).LE.10E-10) THEN 194 : IAPEX=2 195 : C WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX 196 : C WRITE(*,*)'Dans l element',IBB 197 : C WRITE(*,*)'et au point d intégration',IGAU 198 : DO 52 I=1,NSTRSS 199 : DO 52 J=1,NSTRSS 200 : AI(I,J)=0.D0 201 : 52 CONTINUE 202 : AI(1,1)=1./3. 203 : AI(1,2)=AI(1,1) 204 : AI(1,2)=AI(1,1) 205 : AI(1,3)=AI(1,1) 206 : AI(2,1)=AI(1,1) 207 : AI(2,2)=AI(1,1) 208 : AI(2,3)=AI(1,1) 209 : AI(3,1)=AI(1,1) 210 : AI(3,2)=AI(1,1) 211 : AI(3,3)=AI(1,1) 212 : GOTO 75 213 : ENDIF 214 : C 215 : C ************** Mise a jour des contraintes *************** 216 : C 217 : C ---------------- calcul de la matrice A ------------------ 218 : C 219 : DO 60 I=1,NSTRSS 220 : DO 60 J=1,NSTRSS 221 : A(I,J)=0.D0 222 : 60 CONTINUE 223 : C 224 : DG=YOUN/(1.D0+XNU) 225 : C 226 : A(1,1)=1.D0+2.*(DLAM1(1)*DG)/2.D0/DPHI1 227 : A(1,1)=A(1,1)+2.*(DLAM1(2)*DG)/2.D0/DPHI2 228 : A(2,2)=A(1,1) 229 : A(3,3)=A(1,1) 230 : A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1-(DLAM1(2)*DG)/2.D0/DPHI2 231 : A(1,3)=A(1,2) 232 : A(2,1)=A(1,2) 233 : A(2,3)=A(1,2) 234 : A(3,1)=A(1,2) 235 : A(3,2)=A(1,2) 236 : A(4,4)=1.D0+3.*(DLAM1(1)*DG)/2.D0/DPHI1 237 : A(4,4)=A(4,4)+3.*(DLAM1(2)*DG)/2.D0/DPHI2 238 : A(5,5)=A(4,4) 239 : A(6,6)=A(4,4) 240 : C 241 : C -------------- invertion de la matrice A ----------------- 242 : C 242 : CALL ZERO(AI,6,6) 242 : CALL ZERO(AIM,6,6) 243 : C 244 : DO 70 I=1,3 245 : DO 70 J=1,3 246 : AIM(I,J)=A(I,J) 247 : 70 CONTINUE 247 : CALL INVMA2(AIM,3,ISING) 248 : IF (ISING.EQ.1) THEN 249 : WRITE(*,*)'MATRICE AIM singuliere ds ACTI3' 250 : ENDIF 251 : DO 72 I=1,3 252 : DO 72 J=1,3 253 : AI(I,J)=AIM(I,J) 254 : 72 CONTINUE 255 : AI(4,4) = 1./A(4,4) 256 : AI(5,5) = 1./A(5,5) 257 : AI(6,6) = 1./A(6,6) 258 : C 259 : C -------------- mise a jour des contraintes ------------ 260 : C 261 : 75 CONTINUE 262 : C 263 : DO 80 I=1,NSTRSS 264 : DEPSI(I)=SIGE(I)-DLAM1(1)*VEC1(I) 265 : DEPSI(I)=DEPSI(I)-DLAM1(2)*VEC2(I) 266 : 80 CONTINUE 267 : C 268 : DO 90 I=1,NSTRSS 269 : SIGF(I)=0.0D+00 270 : DO 90 J=1,NSTRSS 271 : SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 272 : 90 CONTINUE 273 : C 274 : C ******** Verification des criteres **************** 275 : C 275 : CALL DRUTRA(SIGF,SEQTT,BETINSA) 276 : FCRI1(1) = SEQTT - SEQ1 276 : CALL DRUCOM(SIGF,SEQCC,BETINSA) 277 : FCRI1(2) = SEQCC - SEQ2 278 : C 279 : IF (IAPEX.EQ.2) FCRI1(1)=0.D0 280 : C 281 : IF (IBROY.EQ.0.AND.(DABS(FCRI1(1)).GE.CRIMAX.OR 282 : * .DABS(FCRI1(2)).GE.CRIMAX)) THEN 283 : C WRITE(*,*)'****************************************' 284 : C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN' 285 : C WRITE(*,*)'on passe donc a la secante' 286 : C WRITE(*,*)'Dans l element',IBB 287 : C WRITE(*,*)'et au point d intégration',IGAU 288 : C WRITE(*,*)'CRIMAX=',CRIMAX 289 : C WRITE(*,*)'****************************************' 290 : ITER=ITR 291 : ENDIF 292 : C 293 : C ******* Compteur sur la methode de resolution **** 294 : C 295 : IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN 296 : IBROY=1 297 : ITANG=1 298 : ITER=1 299 : IAPEX=0 300 : IET=0 301 : IEC=0 302 : GOTO 12 303 : ENDIF 304 : C 305 : C ******* non convergence ************************** 306 : C 307 : IF ((DABS(FCRI1(1)).GT.PRB.OR.DABS(FCRI1(2)).GT.PRB) 308 : *.AND.ITER.LT.ITR) THEN 309 : IF (IBROY.EQ.0) THEN 309 : CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1) 310 : DLAM0(1)=DLAM1(1) 311 : DLAM0(2)=DLAM1(2) 312 : FCRI0(1)=FCRI1(1) 313 : FCRI0(2)=FCRI1(2) 314 : IF (ITER.GE.(ITR-1)) THEN 315 : C WRITE(*,*)'***********************' 316 : C WRITE(*,*)'BROYDEN n a pas aboutit' 317 : C WRITE(*,*)'ITER=',ITER 318 : C WRITE(*,*)'FCRIT=',FCRI0(1) 319 : C WRITE(*,*)'FCRIC=',FCRI0(2) 320 : C WRITE(*,*)'Dans l element',IBB 321 : C WRITE(*,*)'et au point d intégration',IGAU 322 : C WRITE(*,*)'***********************' 323 : ENDIF 324 : ITER=ITER+1 325 : GOTO 40 326 : ENDIF 327 : IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN 328 : DLAM0(1)=DLAM1(1) 329 : DLAM0(2)=DLAM1(2) 330 : FCRI0(1)=FCRI1(1) 331 : FCRI0(2)=FCRI1(2) 332 : IF (ITER.GE.(ITR-5)) THEN 333 : C WRITE(*,*)'ITER=',ITER 334 : C WRITE(*,*)'FCRIT=',FCRI0(1) 335 : C WRITE(*,*)'FCRIC=',FCRI0(2) 336 : ENDIF 337 : ITER=ITER+1 338 : GOTO 18 339 : ENDIF 340 : ENDIF 341 : IF (ITER.GE.ITR.AND.(DABS(FCRI1(1)).GT.PRB2.OR. 342 : * DABS(FCRI1(2)).GT.PRB2)) THEN 343 : WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3' 344 : WRITE(*,*)'Dans l element',IBB 345 : WRITE(*,*)'et au point d intégration',IGAU 346 : WRITE(*,*)'FCRIT=',FCRI0(1) 347 : WRITE(*,*)'FCRIC=',FCRI0(2) 348 : WRITE(*,*)'IPLA=',IPLA 349 : WRITE(*,*)'IFIS=',IFIS 350 : C STOP 351 : ENDIF 352 : C 353 : C **************** Fin des iterations internes ******************* 354 : C 355 : DKT=DK1 356 : DKC=DK2 357 : SEQT=SEQ1 358 : SEQC=SEQ2 359 : C 360 : C ******************************************************************** 361 : 100 CONTINUE 362 : IET=0 363 : IEC=0 364 : CRIMAX=0.D0 365 : C 366 : RETURN 367 : END 368 : 369 : 370 : 371 : 372 :
© Cast3M 2003 - All rights reserved.
Disclaimer