1 : C DOHOT1 SOURCE FANDEUR 09/09/23 21:15:06 6374 2 : 3 : SUBROUTINE DOHOT1(IVAMAT,NMATT,IVACON,NSTRS,IVARI,NVART,TRAC, 4 : & LTRAC,IGAU,IB,MATE,MAPL,XPREC,DTPS,IFOU,LHOOK, 5 : & DDHOOK,IRET) 6 : *_______________________________________________________________________ 7 : * 8 : * entrees : 9 : * -------- 10 : * 11 : * ivamat pointeur sur un segment mptval contenant les materiaux 12 : * nmatt nombre de composantes du materiau 13 : * ivacon pointeur sur un segment mptval contenant les contraintes 14 : * nstrs nombre de composantes de contraintes 15 : * ivari pointeur sur un segment mptval contenant les variables 16 : * internes 17 : * nvart nombre de composantes de variables internes 18 : * trac courbe de traction 19 : * ltrac taille du tableau trac 20 : * igau numero du point de gauss ou l'on veut le mat. de hooke 21 : * ib numero de l'element ou l'on veut le mat. de hooke 22 : * mate numeror du materiau lineaire 23 : * mapl numero du materiau plastique 24 : * xprec flottant precision 25 : * dtps flottant pas de temps pour les modeles visqueux 26 : * ifou option ifour de ccoptio 27 : * lhook taille de la matrice de hooke 28 : * 29 : * sorties : 30 : * -------- 31 : * 32 : * ddhook matrice de hooke tangente pour l'element ib 33 : * iret = 1 si ok 34 : * = < 0 : no de l erreur a appeler dans ktanga 35 : * = 0 : autre erreur 36 : * 37 : * passage aux nouveaux chamelem par jm campenon le 05/91 38 : *_______________________________________________________________________ 39 : * 40 : IMPLICIT INTEGER(I-N) 41 : IMPLICIT REAL*8(A-H,O-Z) 42 : * 43 : -INC CCOPTIO 44 : -INC SMCHAML 45 : * 46 : SEGMENT MPTVAL 47 : INTEGER IPOS(NS) ,NSOF(NS) 48 : INTEGER IVAL(NCOSOU) 49 : CHARACTER*16 TYVAL(NCOSOU) 50 : ENDSEGMENT 51 : * 52 : SEGMENT WRK5 53 : INTEGER NTRAT,NTRAC 54 : ENDSEGMENT 55 : * 56 : SEGMENT WTRAV 57 : REAL*8 TXR(IDIM,IDIM) 58 : ENDSEGMENT 59 : 60 : SEGMENT WRKMAT 61 : REAL*8 VALMAT(NMATT) 62 : REAL*8 D1HO(LHOOK,LHOOK),ROTH(LHOOK,LHOOK) 63 : REAL*8 XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM) 64 : ENDSEGMENT 65 : * 66 : PARAMETER(XZER=0.D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0,QUATRE=4.D0) 67 : PARAMETER(CINQ=5.D0,SIX=6.D0,HUIT=8.D0,XNEUF=9.D0,DOUZE=12.D0) 68 : PARAMETER(UNTIER=.33333 33333 33333D0,TRDEMI=1.5D0) 69 : PARAMETER(DETIER=.66666 66666 66666D0,UNDEMI=0.5D0) 70 : PARAMETER(UNQU=.25D0) 71 : * 72 : DIMENSION DDHOOK(LHOOK,*) 73 : DIMENSION TRAC(*) 74 : * 75 : * arguments pour l appel a DOHMAS et DOHMAO 76 : DIMENSION VELA(2) 77 : CHARACTER*8 CMATE 78 : * 79 : DIMENSION XSTRS(8),DEVIAT(8),XMAT(30) 80 : DIMENSION DDF(8),DDG(8),DHDF(8),DHDG(8) 81 : DIMENSION XVAR(14),XKTA(3,3) 82 : DIMENSION PMATRIC(6,6),GMATRIC(6,6),GDMATR(6,6),GINV(6,6),COEF3(6) 83 : * matrice P pour le modele viscoplastique parfait 84 : * la matrice sert a passer de la contrainte a la contrainte deviatorique 85 : DATA PMATRIC /2.D0,-1.D0,-1.D0,0.D0,0.D0,0.D0, 86 : & -1.D0,2.D0,-1.D0,0.D0,0.D0,0.D0, 87 : & -1.D0,-1.D0,2.D0,0.D0,0.D0,0.D0, 88 : & 0.D0, 0.D0,0.D0,6.D0,0.D0,0.D0, 89 : & 0.D0, 0.D0,0.D0,0.D0,6.D0,0.D0, 90 : & 0.D0, 0.D0,0.D0,0.D0,0.D0,6.D0/ 91 : * 92 : DATA COEF3 /1.D0,1.D0,1.D0,2.D0,2.D0,2.D0/ 93 : * 94 : *----------------------------------------------------------------------- 95 : * formulation massive (MFR=1) 96 : * Pour autres formulations (ex : coques) voir dans dohota.eso 97 : *----------------------------------------------------------------------- 98 : IF ((MATE.EQ.1).AND.(MAPL.EQ.54)) THEN 99 : WRK5 = IRET 100 : ELSE IF ((MATE.EQ.2).OR.(MATE.EQ.3).OR.(MATE.EQ.4)) THEN 101 : WTRAV = IRET 102 : ENDIF 103 : * 104 : IRET=0 104 : CALL ZERO(DDHOOK,LHOOK,LHOOK) 104 : CALL ZERO(XMAT,30,1) 105 : * 106 : *----------------------------------------------------------------------- 107 : * MATERIAU ISOTROPE 108 : *----------------------------------------------------------------------- 109 : IF (MATE.EQ.1) THEN 110 : *----------------------------------------------------------------------- 111 : C Initialisations 112 : XXHH =0.D0 113 : SIGY =0.D0 114 : ALP =0.D0 115 : C Initialisation de ILOGEL : sortie de PENT1 116 : C =1 on doit/peut calculer DDHOOK tangent ; =0 DDHOOK elastique 117 : ILOGEL = 0 118 : * 119 : C Contraintes 120 : MPTVAL=IVACON 121 : DO 50 ICOMP=1,NSTRS 122 : MELVAL=IVAL(ICOMP) 123 : IGMN=MIN(IGAU,VELCHE(/1)) 124 : IBMN=MIN(IB ,VELCHE(/2)) 125 : XSTRS(ICOMP)=VELCHE(IGMN,IBMN) 126 : 50 CONTINUE 127 : C 128 : C Deformation plastique 129 : C 130 : MPTVAL=IVARI 131 : MELVAL=IVAL(1) 132 : IGMN=MIN(IGAU,VELCHE(/1)) 133 : IBMN=MIN(IB ,VELCHE(/2)) 134 : EPS =VELCHE(IGMN,IBMN) 135 : C si autres variables internes 136 : IF (NVART.EQ.NSTRS+1) THEN 137 : DO 60 ICOMP=1,NSTRS 138 : MELVAL=IVAL(ICOMP+1) 139 : IGMN=MIN(IGAU,VELCHE(/1)) 140 : IBMN=MIN(IB ,VELCHE(/2)) 141 : XSTRS(ICOMP)= XSTRS(ICOMP)-VELCHE(IGMN,IBMN) 142 : 60 CONTINUE 143 : ENDIF 144 : C 145 : C Recuperation des caracteristiques elastiques isotropes 146 : C 147 : MPTVAL=IVAMAT 148 : C Module d Young 149 : MELVAL=IVAL(1) 150 : IGMN=MIN(IGAU,VELCHE(/1)) 151 : IBMN=MIN(IB ,VELCHE(/2)) 152 : YOU =VELCHE(IGMN,IBMN) 153 : C Coefficient de Poisson 154 : MELVAL=IVAL(2) 155 : IGMN=MIN(IGAU,VELCHE(/1)) 156 : IBMN=MIN(IB ,VELCHE(/2)) 157 : XNU =VELCHE(IGMN,IBMN) 158 : * 159 : *----------------------------------------------------------------------- 160 : * Recuperation des caracteristiques 161 : * Calcul eventuel de la contrainte equivalente 162 : *----------------------------------------------------------------------- 163 : C Plastique parfait 164 : IF (MAPL.EQ.1) THEN 165 : * 166 : MPTVAL=IVAMAT 167 : MELVAL=IVAL(3) 168 : IGMN=MIN(IGAU,VELCHE(/1)) 169 : IBMN=MIN(IB ,VELCHE(/2)) 170 : SIGY =VELCHE(IGMN,IBMN) 171 : * 172 : C Calcul de la contrainte equivalente 173 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 174 : & + XSTRS(3)*XSTRS(3) 175 : YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3) 176 : & + XSTRS(3)*XSTRS(1) 177 : SIGET2 = XXXX - YYYY 178 : IF (NSTRS.GE.4) THEN 179 : DO 61 I=4,NSTRS 180 : SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I) 181 : 61 CONTINUE 182 : ENDIF 183 : SIGET = SQRT(SIGET2) 184 : * 185 : C A t on plastifie ? 185 : CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC, 186 : & YOUTA,ILOGEL) 187 : IF (ILOGEL.EQ.-1) THEN 188 : IRET=-275 189 : RETURN 190 : ENDIF 191 : * 192 : C Drucker simple : ALFA * Tr(S) + Seq = K 193 : ELSE IF (MAPL.EQ.3) THEN 194 : * 195 : MPTVAL=IVAMAT 196 : MELVAL=IVAL(3) 197 : IGMN=MIN(IGAU,VELCHE(/1)) 198 : IBMN=MIN(IB ,VELCHE(/2)) 199 : XLTR =VELCHE(IGMN,IBMN) 200 : MELVAL=IVAL(4) 201 : IGMN=MIN(IGAU,VELCHE(/1)) 202 : IBMN=MIN(IB ,VELCHE(/2)) 203 : XLCS =VELCHE(IGMN,IBMN) 204 : * 205 : DEN = ABS(XLCS) + XLTR 206 : IF(DEN.EQ.0.D0) THEN 207 : IRET=-366 208 : RETURN 209 : ENDIF 210 : C ALP = ALFA ; XXHH = BETA = 1 ; SIGY = K 211 : ALP = (ABS(XLCS) - XLTR)/DEN 212 : XXHH = UN 213 : SIGY = DEUX*ABS(XLCS)*XLTR/DEN 214 : * 215 : * XMAT(1) = 2.0D0*ABS(ALP)*SIGY/DEN 216 : * XMAT(2) = (ABS(ALP) - SIGY)/DEN 217 : * XMAT(3) = 1.D0 218 : XMAT(1) = ALP 219 : XMAT(2) = XXHH 220 : XMAT(3) = SIGY 221 : XMAT(4)=XMAT(1) 222 : XMAT(5)=XMAT(2) 223 : XMAT(6)=XMAT(1) 224 : XMAT(7)=XMAT(2) 225 : XMAT(8)=XMAT(3) 226 : XMAT(9)=0.D0 227 : * 228 : C Calcul de la contrainte equivalente 229 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 230 : & + XSTRS(3)*XSTRS(3) 231 : YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3) 232 : & + XSTRS(3)*XSTRS(1) 233 : SIGET2 = XXXX - YYYY 234 : IF (NSTRS.GE.4) THEN 235 : DO 62 I=4,NSTRS 236 : SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I) 237 : 62 CONTINUE 238 : ENDIF 239 : SIGET = SQRT(SIGET2) 240 : * 241 : C A t on plastifie ? 242 : SIDD = ALP * (XSTRS(1)+XSTRS(2)+XSTRS(3)) + SIGET 242 : CALL PENT1(EPS,SIDD,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC, 243 : & YOUTA,ILOGEL) 244 : IF (ILOGEL.EQ.-1) THEN 245 : IRET=-275 246 : RETURN 247 : ENDIF 248 : * 249 : C Cinematique 250 : ELSE IF (MAPL.EQ.4) THEN 251 : * 252 : MPTVAL=IVAMAT 253 : MELVAL=IVAL(3) 254 : IGMN=MIN(IGAU,VELCHE(/1)) 255 : IBMN=MIN(IB ,VELCHE(/2)) 256 : SIGY =VELCHE(IGMN,IBMN) 257 : * 258 : MELVAL=IVAL(4) 259 : IGMN=MIN(IGAU,VELCHE(/1)) 260 : IBMN=MIN(IB ,VELCHE(/2)) 261 : XXHH=VELCHE(IGMN,IBMN) 262 : * 263 : C Calcul de la contrainte equivalente 264 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 265 : & + XSTRS(3)*XSTRS(3) 266 : YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3) 267 : & + XSTRS(3)*XSTRS(1) 268 : SIGET2 = XXXX - YYYY 269 : IF (NSTRS.GE.4) THEN 270 : DO 63 I=4,NSTRS 271 : SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I) 272 : 63 CONTINUE 273 : ENDIF 274 : SIGET = SQRT(SIGET2) 275 : * 276 : C A t on plastifie ? 276 : CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC, 277 : & YOUTA,ILOGEL) 278 : IF (ILOGEL.EQ.-1) THEN 279 : IRET=-275 280 : RETURN 281 : ENDIF 282 : * 283 : C Plastique isotrope 284 : ELSE IF (MAPL.EQ.5) THEN 285 : * 286 : C Calcul de la contrainte equivalente 287 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 288 : & + XSTRS(3)*XSTRS(3) 289 : YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3) 290 : & + XSTRS(3)*XSTRS(1) 291 : SIGET2 = XXXX - YYYY 292 : IF (NSTRS.GE.4) THEN 293 : DO 64 I=4,NSTRS 294 : SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I) 295 : 64 CONTINUE 296 : ENDIF 297 : SIGET = SQRT(SIGET2) 298 : * 299 : C A t on plastifie ? 299 : CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC, 300 : & YOUTA,ILOGEL) 301 : IF (ILOGEL.EQ.-1) THEN 302 : IRET=-275 303 : RETURN 304 : ENDIF 305 : * 306 : C Drucker Prager 307 : ELSE IF (MAPL.EQ.15) THEN 308 : * 309 : MPTVAL=IVAMAT 310 : DO 1106 IC=3,NMATT 311 : MELVAL=IVAL(IC) 312 : XMAT(IC-2)=VELCHE(IGMN,IBMN) 313 : 1106 CONTINUE 314 : * 315 : ** ALP = ALFA ou ETA ; XXHH = BETA ou MU ; SIGY = K ou KL+H*EPS 316 : IF(EPS.EQ.0.D0) THEN 317 : ALP=XMAT(6) 318 : XXHH=XMAT(7) 319 : SIGY=XMAT(8) 320 : ELSE 321 : ALP=XMAT(1) 322 : XXHH=XMAT(2) 323 : SIGY=XMAT(3)+XMAT(9)*EPS 324 : ENDIF 325 : * 326 : C Calcul de la contrainte equivalente 327 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 328 : & + XSTRS(3)*XSTRS(3) 329 : YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3) 330 : & + XSTRS(3)*XSTRS(1) 331 : SIGET2 = XXXX - YYYY 332 : IF (NSTRS.GE.4) THEN 333 : DO 66 I=4,NSTRS 334 : SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I) 335 : 66 CONTINUE 336 : ENDIF 337 : SIGET = SQRT(SIGET2) 338 : * 339 : C A t on plastifie ? 340 : SIDD = ALP * (XSTRS(1)+XSTRS(2)+XSTRS(3)) + XXHH * SIGET 340 : CALL PENT1(EPS,SIDD,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC, 341 : & YOUTA,ILOGEL) 342 : IF (ILOGEL.EQ.-1) THEN 343 : IRET=-275 344 : RETURN 345 : ENDIF 346 : * 347 : C Viscoplastique parfait 348 : ELSE IF (MAPL.EQ.43) THEN 349 : * 350 : C Recuperation des caracteristiques 351 : MPTVAL=IVAMAT 352 : MELVAL = IVAL(3) 353 : IGMN=MIN(IGAU,VELCHE(/1)) 354 : IBMN=MIN(IB ,VELCHE(/2)) 355 : SIGY = VELCHE(IGMN,IBMN) 356 : MELVAL = IVAL(5) 357 : IGMN=MIN(IGAU,VELCHE(/1)) 358 : IBMN=MIN(IB ,VELCHE(/2)) 359 : CK = VELCHE(IGMN,IBMN) 360 : MELVAL = IVAL(4) 361 : IGMN=MIN(IGAU,VELCHE(/1)) 362 : IBMN=MIN(IB ,VELCHE(/2)) 363 : CN = VELCHE(IGMN,IBMN) 364 : ILOGEL = 1 365 : * 366 : C Betocy 367 : ELSE IF (MAPL.EQ.54) THEN 368 : C 369 : MPTVAL=IVAMAT 370 : DO 1107 IC=1,NMATT-2 371 : MELVAL=IVAL(IC) 372 : IGMN=MIN(IGAU,VELCHE(/1)) 373 : IBMN=MIN(IB ,VELCHE(/2)) 374 : XMAT(IC)=VELCHE(IGMN,IBMN) 375 : 1107 CONTINUE 376 : C Lecture Variables internes 377 : MPTVAL=IVARI 378 : DO 1108 IC=1,NVART 379 : MELVAL=IVAL(IC) 380 : IGMN=MIN(IGAU,VELCHE(/1)) 381 : IBMN=MIN(IB ,VELCHE(/2)) 382 : XVAR(IC)=VELCHE(IGMN,IBMN) 383 : 1108 CONTINUE 384 : ILOGEL = 1 385 : * 386 : C Rotating Crack 387 : ELSE IF (MAPL.EQ.55) THEN 388 : C Lecture materiau 389 : MPTVAL=IVAMAT 390 : DO 1109 IC=1,NMATT 391 : MELVAL=IVAL(IC) 392 : IGMN=MIN(IGAU,VELCHE(/1)) 393 : IBMN=MIN(IB ,VELCHE(/2)) 394 : XMAT(IC)=VELCHE(IGMN,IBMN) 395 : 1109 CONTINUE 396 : C Lecture Variables internes 397 : MPTVAL=IVARI 398 : DO 1110 IC=1,NVART 399 : MELVAL=IVAL(IC) 400 : IGMN=MIN(IGAU,VELCHE(/1)) 401 : IBMN=MIN(IB ,VELCHE(/2)) 402 : XVAR(IC)=VELCHE(IGMN,IBMN) 403 : 1110 CONTINUE 404 : ILOGEL = 1 405 : * 406 : * BCN (start) 407 : ELSE IF ((MAPL.EQ.111).OR.(MAPL.EQ.112).OR.(MAPL.EQ.113)) THEN 408 : IF (MAPL.EQ.112.OR.MAPL.EQ.113) THEN 409 : C J2 (MAPL=112) or Rounded Hyperbolic Mohr-Coulomb (MAPL=113) 410 : C get auxiliary variable for CTM 411 : MPTVAL=IVARI 412 : MELVAL=IVAL(2) 413 : IGMN=MIN(IGAU,VELCHE(/1)) 414 : IBMN=MIN(IB ,VELCHE(/2)) 415 : AUX1 =VELCHE(IGMN,IBMN) 416 : ELSE IF (MAPL.EQ.111) THEN 417 : C MRS-Lade (MAPL=111) 418 : C get auxiliary variable for CTM 419 : MPTVAL=IVARI 420 : MELVAL=IVAL(1) 421 : IGMN=MIN(IGAU,VELCHE(/1)) 422 : IBMN=MIN(IB ,VELCHE(/2)) 423 : EPSCON =VELCHE(IGMN,IBMN) 424 : MPTVAL=IVARI 425 : MELVAL=IVAL(2) 426 : IGMN=MIN(IGAU,VELCHE(/1)) 427 : IBMN=MIN(IB ,VELCHE(/2)) 428 : EPSCAP =VELCHE(IGMN,IBMN) 429 : MPTVAL=IVARI 430 : MELVAL=IVAL(3) 431 : IGMN=MIN(IGAU,VELCHE(/1)) 432 : IBMN=MIN(IB ,VELCHE(/2)) 433 : AUX1CON =VELCHE(IGMN,IBMN) 434 : MPTVAL=IVARI 435 : MELVAL=IVAL(4) 436 : IGMN=MIN(IGAU,VELCHE(/1)) 437 : IBMN=MIN(IB ,VELCHE(/2)) 438 : AUX1CAP =VELCHE(IGMN,IBMN) 439 : ENDIF 440 : C The three models: read material constants in XMAT 441 : XMAT(1)=YOU 442 : XMAT(2)=XNU 443 : XMAT(3)=0.D0 444 : XMAT(4)=0.D0 445 : MPTVAL=IVAMAT 446 : DO IC=3,NMATT 447 : MELVAL=IVAL(IC) 448 : IGMN=MIN(IGAU,VELCHE(/1)) 449 : IBMN=MIN(IB ,VELCHE(/2)) 450 : XMAT(IC+2)=VELCHE(IGMN,IBMN) 451 : ENDDO 452 : ILOGEL = 1 453 : ENDIF 454 : * BCN (end) 455 : * 456 : *----------------------------------------------------------------------- 457 : * Calcul de la matrice elastique 458 : *----------------------------------------------------------------------- 459 : VELA(1) = YOU 460 : VELA(2) = XNU 461 : * MATE = 1 --> CMATE = ISOTROPE 462 : CMATE = 'ISOTROPE' 462 : CALL DOHMAS(VELA,CMATE,IFOU,LHOOK,1,DDHOOK,IRETOU) 463 : IF (IRETOU.EQ.0) THEN 464 : IRET = -81 465 : RETURN 466 : ENDIF 467 : IRET = 1 468 : * 469 : * ILOGEL = 0 on a fini, =1 on calcule la matrice tangente 470 : IF (ILOGEL.EQ.0) RETURN 471 : * series de Fourier, on rend la matrice elastique 472 : IF (IFOU.EQ.1) RETURN 473 : * 474 : *----------------------------------------------------------------------- 475 : * Calcul de la matrice tangente 476 : *----------------------------------------------------------------------- 477 : c MODELE MRS_LADE 478 : IF (MAPL.EQ.111) THEN 479 : * 480 : * notice de l operateur MODE: "option de calcul : plan defo,axis" 481 : * sinon, on rend la matrice elastique 482 : IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN 483 : nescri =0 484 : nues =6 485 : C mrs-lade requiere siempre derivacion numerica 486 : nnumer=3 487 : deltax=2.D0**(int(log10(1.D-6)/log10(2.D0))) 487 : call mtc_MRSMAC(DDHOOK,LHOOK,XSTRS,NSTRS,EPSCON,EPSCAP, 488 : & AUX1CON,AUX1CAP,XMAT,nescri,nues, 489 : & nnumer,deltax,kerre) 490 : * kerre = 1 : formulacion no disponible 491 : IF (KERRE.EQ.1) THEN 492 : IRET = -328 493 : ENDIF 494 : ENDIF 495 : * 496 : c MODELE J2 497 : ELSE IF (MAPL.EQ.112) THEN 498 : * notice de l operateur MODE: "option de calcul : plan defo,axis" 499 : * sinon, on rend la matrice elastique 500 : IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN 501 : nescri =0 502 : nues =6 502 : call mtc_j2(DDHOOK,LHOOK,XSTRS,NSTRS,EPS,AUX1,XMAT, 503 : & nescri,nues,kerre) 504 : * kerre vaut toujours 0 505 : ENDIF 506 : * 507 : c MODELE RH_COULOMB 508 : ELSE IF (MAPL.EQ.113) THEN 509 : * notice de l operateur MODE: "option de calcul : plan defo,axis" 510 : * sinon, on rend la matrice elastique 511 : IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN 512 : nescri =0 513 : nues =6 513 : call mtc_rhmc(DDHOOK,LHOOK,XSTRS,NSTRS,AUX1,XMAT, 514 : & nescri,nues,kerre) 515 : * kerre vaut toujours 0 516 : ENDIF 517 : * 518 : C Viscoplastique parfait 519 : ELSE IF (MAPL.EQ.43) THEN 520 : * 521 : * tout sauf contraintes planes 522 : * sinon on rend la matrice elastique 523 : IF ((IFOU.EQ.-1).OR.(IFOU.EQ.-3).OR.(IFOU.EQ.0).OR. 524 : & (IFOU.EQ.2).OR.(IFOU.EQ.3).OR.(IFOU.EQ.7).OR.(IFOU.EQ.9) 525 : & .OR.(IFOU.EQ.11).OR.(IFOU.EQ.12).OR.(IFOU.EQ.14).OR. 526 : & (IFOU.EQ.15)) THEN 527 : * SIGMA EQUIVALENT: (3/2(s':s'))**0.5D0 528 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 529 : & + XSTRS(3)*XSTRS(3) 530 : YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3) 531 : & + XSTRS(3)*XSTRS(1) 532 : SIGET2 = XXXX - YYYY 533 : IF (NSTRS.GE.4) THEN 534 : DO 67 I=4,NSTRS 535 : SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I) 536 : 67 CONTINUE 537 : ENDIF 538 : SIGEQ = SQRT(SIGET2) 539 : * 540 : XMOY=(XSTRS(1)+XSTRS(2)+XSTRS(3))*UNTIER 541 : DO 200 IA=1,3 542 : DEVIAT(IA)=XSTRS(IA)- XMOY 543 : 200 CONTINUE 544 : IF (NSTRS.GE.4) THEN 545 : DO 300 IA=4,NSTRS 546 : DEVIAT(IA)=XSTRS(IA) 547 : 300 CONTINUE 548 : ENDIF 549 : * 550 : SYN = ( SIGEQ - SIGY ) / CK 551 : IF (SYN .GT. XPREC) THEN 552 : * on ne travaille que si l'evolution est plastique 553 : * 554 : AUX6 = UNDEMI * YOU / (UN+XNU) 555 : * 556 : * ON CALCULE LA MATRICE G 557 : COEF1 = SYN ** CN / SIGEQ*DTPS*2.D0*AUX6/3.D0 558 : COEF2=SYN**(CN-1.D0)*(SYN/SIGEQ+CN/CK)/SIGET2*DTPS 559 : & *2.D0*AUX6 560 : DO 610 I=1,NSTRS 561 : DO 600 J=1,NSTRS 562 : GMATRIC(I,J)=COEF1*PMATRIC(I,J)+ 563 : & COEF2*DEVIAT(I)*DEVIAT(J)*COEF3(I) 564 : 600 CONTINUE 565 : GMATRIC(I,I)=GMATRIC(I,I)+1.D0 566 : 610 CONTINUE 567 : 567 : CALL INVHOK(GMATRIC,GINV,NSTRS) 568 : * 569 : * on calcule InvG . D 570 : DO 650 I=1,NSTRS 571 : DO 640 J=I,NSTRS 572 : GDMATR(I,J)=0.D0 573 : DO 630 K=1,NSTRS 574 : GDMATR(I,J)=GDMATR(I,J)+GINV(I,K)*DDHOOK(K,J) 575 : 630 CONTINUE 576 : GDMATR(J,I)=GDMATR(I,J) 577 : 640 CONTINUE 578 : 650 CONTINUE 579 : * on remet tout dans DHOOK 580 : DO 670 I=1,NSTRS 581 : DO 660 J=1,NSTRS 582 : DDHOOK(I,J)=GDMATR(I,J) 583 : 660 CONTINUE 584 : 670 CONTINUE 585 : ENDIF 586 : ENDIF 587 : * 588 : C Betocy 589 : ELSE IF (MAPL.EQ.54) THEN 590 : * 591 : * notice de l operateur MODE : "en 2D contraintes planes" 592 : * sinon on rend la matrice elastique 593 : IF (IFOU.EQ.-2) THEN 593 : CALL BETOKT(XMAT,TRAC,XVAR,XSTRS,NTRAT,NTRAC,DDHOOK,KERRE) 594 : * kerre toujours nul 595 : ENDIF 596 : * 597 : C Rotating Crack 598 : ELSE IF (MAPL.EQ.55) THEN 599 : * 600 : * notice de MATE : "la direction de fissuration change a chaque pas" 601 : * possible en 1D ? dans ROTKTA DDHOOK(4,4), en 1D on va jusque 3 602 : * sinon on rend la matrice elastique 603 : IF (IFOU.EQ.-2) THEN 603 : CALL ROTKTA(XMAT,XVAR,XSTRS,DDHOOK,KERRE) 604 : * kerre = 2 "impossible d inverser" 605 : IF (KERRE.EQ.2) IRET = 0 606 : ENDIF 607 : * 608 : C Plastique parfait, cinematique, isotrope 609 : ELSE IF ((MAPL.EQ.1).OR.(MAPL.EQ.4).OR.(MAPL.EQ.5)) THEN 610 : * 611 : IF (IFOU.EQ.-2) THEN 612 : * 2D contraintes planes 613 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 614 : YYYY = XSTRS(1)*XSTRS(2) 615 : ZZZZ = XSTRS(4)*XSTRS(4) 616 : AUX4=UNTIER/(UN-XNU*XNU) 617 : U=(DEUX*XSTRS(1)-XSTRS(2)+XNU*(DEUX*XSTRS(2)- 618 : & XSTRS(1)))*AUX4 619 : V=(DEUX*XSTRS(2)-XSTRS(1)+XNU*(DEUX*XSTRS(1)- 620 : & XSTRS(2)))*AUX4 621 : W= XSTRS(4)/(UN+XNU) 622 : A= CINQ*XXXX - HUIT*YYYY 623 : B= CINQ*YYYY - DEUX*XXXX 624 : C= DETIER*(XXXX-YYYY)+DEUX*ZZZZ 625 : D= (A+DEUX*XNU*B)*AUX4 + SIX*ZZZZ/(UN+XNU) 626 : XNUM=TROIS*YOU*(YOU-YOUTA) 627 : DENOM=DEUX*YOUTA*C + (YOU - YOUTA)*D 628 : BETA=XNUM/DENOM 629 : * 630 : DDHOOK(1,1)=DDHOOK(1,1)-BETA*U*U 631 : DDHOOK(2,1)=DDHOOK(2,1)-BETA*V*U 632 : DDHOOK(4,1)=DDHOOK(4,1)-BETA*W*U 633 : DDHOOK(1,2)=DDHOOK(1,2)-BETA*U*V 634 : DDHOOK(2,2)=DDHOOK(2,2)-BETA*V*V 635 : DDHOOK(4,2)=DDHOOK(4,2)-BETA*W*V 636 : DDHOOK(1,4)=DDHOOK(1,4)-BETA*U*W 637 : DDHOOK(2,4)=DDHOOK(2,4)-BETA*V*W 638 : DDHOOK(4,4)=DDHOOK(4,4)-BETA*W*W 639 : * 640 : ELSE IF ((IFOU.EQ.4).OR.(IFOU.EQ.8)) THEN 641 : * DYCZ,GYCZ 642 : XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2) 643 : YYYY = XSTRS(1)*XSTRS(2) 644 : AUX4=UNTIER/(UN-XNU*XNU) 645 : U=(DEUX*XSTRS(1)-XSTRS(2)+XNU*(DEUX*XSTRS(2)- 646 : & XSTRS(1)))*AUX4 647 : V=(DEUX*XSTRS(2)-XSTRS(1)+XNU*(DEUX*XSTRS(1)- 648 : & XSTRS(2)))*AUX4 649 : A= CINQ*XXXX - HUIT*YYYY 650 : B= CINQ*YYYY - DEUX*XXXX 651 : C= DETIER*(XXXX-YYYY) 652 : D= (A+DEUX*XNU*B)*AUX4 653 : XNUM=TROIS*YOU*(YOU-YOUTA) 654 : DENOM=DEUX*YOUTA*C + (YOU - YOUTA)*D 655 : BETA=XNUM/DENOM 656 : * 657 : DDHOOK(1,1)=DDHOOK(1,1)-BETA*U*U 658 : DDHOOK(2,1)=DDHOOK(2,1)-BETA*V*U 659 : DDHOOK(1,2)=DDHOOK(1,2)-BETA*U*V 660 : DDHOOK(2,2)=DDHOOK(2,2)-BETA*V*V 661 : * 662 : ELSE IF ((IFOU.EQ.5).OR.(IFOU.EQ.10).OR.(IFOU.EQ.13) 663 : & .OR.(IFOU.EQ.6)) THEN 664 : * CYDZ,CYGZ,AXCZ,CYCZ 665 : * comment adapter les formules 2D contraintes planes au 1D CY.. ? 666 : * pour le moment, on rend la matrice elastique 667 : * 668 : ELSE 669 : * tout sauf contraintes planes (1D aussi) 670 : * deviateur 671 : XMOY=(XSTRS(1)+XSTRS(2)+XSTRS(3))*UNTIER 672 : DO 201 IA=1,3 673 : DEVIAT(IA)=XSTRS(IA)- XMOY 674 : 201 CONTINUE 675 : IF (NSTRS.GE.4) THEN 676 : DO 301 IA=4,NSTRS 677 : DEVIAT(IA)=XSTRS(IA) 678 : 301 CONTINUE 679 : ENDIF 680 : * 681 : XNUM = XNEUF * YOU * ( YOU - YOUTA) 682 : XXXX = TROIS * ( YOU - YOUTA ) + 683 : & DEUX * ( UN + XNU ) * YOUTA 684 : DENOM= DEUX * SIGET2 * ( UN + XNU ) * XXXX 685 : BETA = XNUM / DENOM 686 : * 687 : DO 400 ICOMP1=1,NSTRS 688 : DO 400 ICOMP2=1,NSTRS 689 : DDHOOK(ICOMP1,ICOMP2)= DDHOOK(ICOMP1,ICOMP2) - 690 : & BETA*DEVIAT(ICOMP1)*DEVIAT(ICOMP2) 691 : 400 CONTINUE 692 : * 693 : ENDIF 694 : * 695 : C Drucker simple, Drucker Prager 696 : ELSE IF ((MAPL.EQ.3).OR.(MAPL.EQ.15)) THEN 697 : * 698 : * si SIGET = 0 (avec tr(S) dans le critere, c est possible) 699 : * on rend la matrice elastique (sinon on divise par SIGET=0) 700 : IF (SIGET.EQ.0.) RETURN 700 : CALL ZERO(DDF,NSTRS,1) 700 : CALL ZERO(DDG,NSTRS,1) 701 : * 702 : IF ((IFOU.EQ.-2).OR.(IFOU.EQ.4).OR.(IFOU.EQ.8)) THEN 703 : * 2D contraintes planes et 1D DYCZ GYCZ 704 : SIGM = (XSTRS(1)-XSTRS(2)*0.5D0)/SIGET 705 : DDF(1)=XMAT(1) + XMAT(2)*SIGM 706 : DDG(1)=XMAT(4) + XMAT(5)*SIGM 707 : SIGM = (XSTRS(2)-XSTRS(1)*0.5D0)/SIGET 708 : DDF(2)=XMAT(1) + XMAT(2)*SIGM 709 : DDG(2)=XMAT(4) + XMAT(5)*SIGM 710 : IF (NSTRS.GE.4) THEN 711 : DDF(4)=6.D0*XMAT(2)*XSTRS(4)/SIGET 712 : DDG(4)=6.D0*XMAT(5)*XSTRS(4)/SIGET 713 : ENDIF 714 : * 715 : ELSE IF ((IFOU.EQ.5).OR.(IFOU.EQ.10).OR.(IFOU.EQ.13)) THEN 716 : * 1D CYDZ CYGZ et AXCZ 717 : SIGM = (XSTRS(1)-XSTRS(3)*0.5D0)/SIGET 718 : DDF(1)=XMAT(1) + XMAT(2)*SIGM 719 : DDG(1)=XMAT(4) + XMAT(5)*SIGM 720 : SIGM = (XSTRS(3)-XSTRS(1)*0.5D0)/SIGET 721 : DDF(3)=XMAT(1) + XMAT(2)*SIGM 722 : DDG(3)=XMAT(4) + XMAT(5)*SIGM 723 : * 724 : ELSE IF (IFOU.EQ.6) THEN 725 : * 1D CYCZ 726 : SIGM = (XSTRS(1)-XSTRS(2)*0.5D0)/SIGET 727 : DDF(1)=XMAT(1) + XMAT(2)*SIGM 728 : DDG(1)=XMAT(4) + XMAT(5)*SIGM 729 : * 730 : ELSE 731 : SIGM = (XSTRS(1)-XSTRS(2)*0.5D0 -XSTRS(3)*0.5D0)/SIGET 732 : DDF(1)=XMAT(1) + XMAT(2)*SIGM 733 : DDG(1)=XMAT(4) + XMAT(5)*SIGM 734 : SIGM = (XSTRS(2)-XSTRS(1)*0.5D0 -XSTRS(3)*0.5D0)/SIGET 735 : DDF(2)=XMAT(1) + XMAT(2)*SIGM 736 : DDG(2)=XMAT(4) + XMAT(5)*SIGM 737 : SIGM = (XSTRS(3)-XSTRS(1)*0.5D0 -XSTRS(2)*0.5D0)/SIGET 738 : DDF(3)=XMAT(1) + XMAT(2)*SIGM 739 : DDG(3)=XMAT(4) + XMAT(5)*SIGM 740 : IF (NSTRS.GE.4) THEN 741 : DO 421 I=4,NSTRS 742 : DDF(I)=6.D0*XMAT(2)*XSTRS(I)/SIGET 743 : DDG(I)=6.D0*XMAT(5)*XSTRS(I)/SIGET 744 : 421 CONTINUE 745 : ENDIF 746 : * 747 : ENDIF 748 : * 749 : DO 422 I=1,NSTRS 750 : DHDF(I)=0.D0 751 : DHDG(I)=0.D0 752 : DO 422 J=1,NSTRS 753 : DHDF(I)=DHDF(I)+DDHOOK(I,J)*DDF(J) 754 : DHDG(I)=DHDG(I)+DDHOOK(I,J)*DDG(J) 755 : 422 CONTINUE 756 : * 757 : DENOM = XMAT(9)* SQRT(2.D0*XMAT(4)*XMAT(4)+ 758 : & XMAT(5)*XMAT(5)) 759 : DO 423 I=1,NSTRS 760 : DENOM = DENOM + DDF(I)*DHDG(I) 761 : 423 CONTINUE 762 : DO 320 ICOMP1=1,NSTRS 763 : DO 320 ICOMP2=1,NSTRS 764 : DDHOOK(ICOMP1,ICOMP2)=DDHOOK(ICOMP1,ICOMP2) 765 : & - DHDG(ICOMP1)*DHDF(ICOMP2)/DENOM 766 : 320 CONTINUE 767 : * 768 : ELSE 769 : * on rend la matrice elastique pour les autres modeles 770 : ENDIF 771 : * 772 : *----------------------------------------------------------------------- 773 : * MATERIAU UNIDIRECTIONNEL (MATE=4) 774 : * Acier unidirectionnel (MAPL=40) 775 : *----------------------------------------------------------------------- 776 : ELSE IF ((MATE.EQ.4).AND.(MAPL.EQ.40).AND.(IFOU.LE.0)) THEN 777 : *----------------------------------------------------------------------- 778 : DO 530 IC=1,3 779 : MPTVAL=IVAMAT 780 : MELVAL=IVAL(IC) 781 : IGMN=MIN(IGAU,VELCHE(/1)) 782 : IBMN=MIN(IB ,VELCHE(/2)) 783 : XMAT(IC)=VELCHE(IGMN,IBMN) 784 : 530 CONTINUE 785 : V1X=TXR(1,1)*XMAT(2)+TXR(1,2)*XMAT(3) 786 : V1Y=TXR(2,1)*XMAT(2)+TXR(2,2)*XMAT(3) 787 : XMAT(2)=V1X 788 : XMAT(3)=V1Y 789 : C Lecture Variables internes 790 : MPTVAL=IVARI 791 : MELVAL=IVAL(4) 792 : IGMN=MIN(IGAU,VELCHE(/1)) 793 : IBMN=MIN(IB ,VELCHE(/2)) 794 : XVAR(1)=VELCHE(IGMN,IBMN) 795 : * 795 : CALL UNIAKT(XMAT,XVAR,IFOU,DDHOOK,KERRE) 796 : * kerre toujours nul 797 : IRET = 1 798 : * 799 : *----------------------------------------------------------------------- 800 : * AUTRE MATERIAU : 801 : * UNIDIRECTIONNEL (sauf MAPL=40 et IFOU = -3,-2,-1,0) 802 : * ORTHOTROPE 803 : * ANISOTROPE 804 : *----------------------------------------------------------------------- 805 : ELSE 806 : *----------------------------------------------------------------------- 807 : IF (MATE.EQ.2) THEN 808 : CMATE = 'ORTHOTRO' 809 : ELSE IF (MATE.EQ.3) THEN 810 : CMATE = 'ANISOTRO' 811 : ELSE IF (MATE.EQ.4) THEN 812 : CMATE = 'UNIDIREC' 813 : ELSE 814 : C cet operateur ne fonctionne pas pour le materiau %m1:8 815 : IRET = -834 816 : RETURN 817 : ENDIF 818 : SEGINI,WRKMAT 819 : MPTVAL = IVAMAT 820 : SEGACT,MPTVAL 821 : DO 9004 IM=1,NMATT 822 : IF (IVAL(IM).NE.0) THEN 823 : MELVAL=IVAL(IM) 824 : IBMN=MIN(IB ,VELCHE(/2)) 825 : IGMN=MIN(IGAU,VELCHE(/1)) 826 : if (ibmn.gt.0.and.igmn.gt.0) then 827 : VALMAT(IM)= VELCHE(IGMN,IBMN) 828 : C* else 829 : C* VALMAT(IM)=0.D0 830 : endif 831 : C* ELSE 832 : C* VALMAT(IM)=0.D0 833 : ENDIF 834 : 9004 CONTINUE 834 : CALL DOHMAO(VALMAT,CMATE,IFOU,IDIM,TXR,XLOC,XGLOB,D1HO, 835 : & ROTH,DDHOOK,LHOOK,1,IRETOU) 836 : SEGSUP,WRKMAT 837 : IF (IRETOU.EQ.0) THEN 838 : IRET = -81 839 : ELSE IF (IRETOU.EQ.1) THEN 840 : IRET = 1 841 : ELSE 842 : IRET = 0 843 : ENDIF 844 : ENDIF 845 : * 846 : RETURN 847 : END 848 : 849 :
© Cast3M 2003 - All rights reserved.
Disclaimer