1 : C DYNAR SOURCE KICH 09/09/11 21:18:03 6505 2 : SUBROUTINE DYNAR(WRK52,WRK53,WRK54,IECOU) 3 : c 4 : C======================================================================= 5 : C 6 : C MODELE VISCOPLASTIQUE VISCOENDOMMAGEABLE POUR LA DYNAMIQUE 7 : C - PROGRAMME PRINCIPAL - 8 : C VERSION 2.0 9 : C 10 : C======================================================================= 11 : C CREATION :F.GATUINGT 12 : C======================================================================= 13 : C=========== Hypopthese: au temps t=0 tout est nul =================== 14 : C======================================================================= 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : -INC CCOPTIO 17 : -INC DECHE 18 : DIMENSION EPSPTDT(6),EPSTE(6),EPSPOPT(6),EPSPO(6),SIGPO(6), 19 : & sigt(6),sigtdt(6),EPST(6),EPSTDT(6),EPSPT(6),EPSTDTE(6), 20 : & S(3,3),EPSIPP(3),ARRAY(3,3),SIGTDTE(6),DFIDSIG(6),EPSPOPTDT(6) 21 : C 22 : PARAMETER (XUnDemi=0.5D0,XZERO=0.D0,UN=1.D0,DEUX=2.D0) 23 : PARAMETER (XPETIT=1.D-12) 24 : PARAMETER (TROIS=3.D0) 25 : C 26 : SEGMENT IECOU 27 : * COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK, 28 : INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7, 29 : C INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM, NYKK, 30 : 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16, 31 : C . NYALF1,NYBET1,NYR, NYA, NYRHO,NSIGY, NNKX, NYKX, IND, 32 : 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24, 33 : C . NSOM, NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT, NCXMAT, 34 : 3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA, 35 : C . LTRAC, MFR, IELE, NHRM, NBNN, NBELEM,ICARA, 36 : 4 icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38, 37 : C . LW2, NDEF, NSTRSS,MFR1, NBGMAT,NELMAT,MSOUPA, 38 : 5 icow39,icow40,icow41,icow42,icow43,icow44 39 : C . NUMAT1,LENDO, NBBB, NNVARI,KERR1, MELEME 40 : INTEGER icow45,icow46,icow47,icow48,icow49,icow50, 41 : . icow51,icow52,icow53,icow54,icow55,icow56 42 : . icow57,icow58 43 : ENDSEGMENT 44 : C 45 : C 46 : C 47 : C 48 : C - Recuperation des donnees materiau 49 : C 50 : Eini=XMAT(1) 51 : XNU=XMAT(2) 52 : D0=XMAT(5) 53 : Q1=XMAT(6) 54 : Q2=XMAT(7) 55 : Q3=XMAT(8) 56 : SIGM0=XMAT(9) 57 : XN=XMAT(10) 58 : XKPLAS = XMAT(12) 59 : XMPLAS = XMAT(11) 60 : XK=XMAT(13) 61 : XMDt=XMAT(14) 62 : XNDt=XMAT(15) 63 : At=XMAT(16) 64 : Bt=XMAT(17) 65 : XMDc=XMAT(18) 66 : XNDc=XMAT(19) 67 : Ac=XMAT(20) 68 : Bc=XMAT(21) 69 : ED0=XMAT(22) 70 : C - Pas de temps 71 : DELTAT=TEMPF-TEMP0 72 : C - Contrainte au temps T 73 : SIGT(1)=SIG0(1) 74 : SIGT(2)=SIG0(2) 75 : SIGT(3)=SIG0(3) 76 : SIGT(4)=SIG0(4) 77 : SIGT(5)=SIG0(5) 78 : SIGT(6)=SIG0(6) 79 : C - Deformation au temps T 80 : EPST(1)=EPST0(1) 81 : EPST(2)=EPST0(2) 82 : EPST(3)=EPST0(3) 83 : EPST(4)=EPST0(4)/2. 84 : EPST(5)=EPST0(5)/2. 85 : EPST(6)=EPST0(6)/2. 86 : C - Increment de deformation 87 : DEPST(1)=DEPST(1) 88 : DEPST(2)=DEPST(2) 89 : DEPST(3)=DEPST(3) 90 : DEPST(4)=DEPST(4)/2. 91 : DEPST(5)=DEPST(5)/2. 92 : DEPST(6)=DEPST(6)/2. 93 : C====================================================================== 94 : C - Subcycling si le pas de temps est trop grand======================= 95 : C====================================================================== 96 : IF (DELTAT.GT.1.D-3) THEN 97 : coef=100 98 : ELSE 99 : coef=max (1.d0 , deltat/ 1.d-5) 100 : ENDIF 101 : nbiter=int(coef) 102 : DO i=1,6 103 : DEPST(i)=DEPST(i)/nbiter 104 : ENDDO 105 : DELTAT=DELTAT/nbiter 106 : l=0 107 : DO WHILE (l.LT.nbiter) 108 : C - Debut des boucles 109 : C======================================================================= 110 : C - Recuperation des variables historiques 111 : C======================================================================= 112 : DT=VAR0(1) 113 : fT=VAR0(2)+D0 114 : EPSPT(1)=VAR0(3) 115 : EPSPT(2)=VAR0(4) 116 : EPSPT(3)=VAR0(5) 117 : EPSPT(4)=VAR0(6) 118 : EPSPT(5)=VAR0(7) 119 : EPSPT(6)=VAR0(8) 120 : SIGMT=VAR0(9)+SIGM0 121 : dtt=VAR0(10) 122 : dct=VAR0(11) 123 : ES=VAR0(12)+ED0 124 : C====================================================================== 125 : C - Deformation au temps t+dt 126 : C====================================================================== 127 : EPSTDT(1)=EPST(1)+DEPST(1) 128 : EPSTDT(2)=EPST(2)+DEPST(2) 129 : EPSTDT(3)=EPST(3)+DEPST(3) 130 : EPSTDT(4)=EPST(4)+DEPST(4) 131 : EPSTDT(5)=EPST(5)+DEPST(5) 132 : EPSTDT(6)=EPST(6)+DEPST(6) 133 : TRACEPSTDT=EPSTDT(1)+EPSTDT(2)+EPSTDT(3) 134 : C====================================================================== 135 : C - Vitesse de Deformation 136 : C====================================================================== 137 : DO i=1,6 138 : EPSPO(i)=DEPST(i)/DELTAT 139 : ENDDO 140 : TRACEPSPO=EPSPO(1)+EPSPO(2)+EPSPO(3) 141 : C======================================================================= 142 : C - Module d'elasticite tangent de la Matrice 143 : C======================================================================= 144 : ETAN=(Eini)/XN*(SIGMt/SIGM0)**(1-XN) 145 : C======================================================================= 146 : C - Coefficients de compressibilite et cisaillement de la matrice sans les pores 147 : C======================================================================= 148 : XKm=Eini/(TROIS*(UN-DEUX*XNU)) 149 : XGm=Eini/(DEUX*(UN+XNU)) 150 : C======================================================================= 151 : C - Coeficients de compressibilite et cisaillement de la matrice avec les pores (Mori-Tanaka) 152 : C======================================================================= 153 : XKporo=4.D0*XKm*XGm*(UN-fT)/(4.D0*XGm+3.D0*XKm*fT) 154 : XGporo=XGm*(UN-fT)/(UN+fT* 155 : & (6.D0*XKm+12.D0*XGm)/(9.D0*XKm+8.D0*XGm)) 156 : C ====================================================================== 157 : C - Coefficients de compressibilite et cisaillement de la matrice endommagee 158 : C ====================================================================== 159 : XKendo=(UN-Dt)*XKporo 160 : XGendo=(UN-Dt)*XGporo 161 : C ====================================================================== 162 : C Calcul des invariants de la contrainte 163 : C ====================================================================== 164 : XI1t=SIGT(1)+SIGT(2)+SIGT(3) 165 : XJ2t=0.5D0*((SIGT(1)-SIGT(2))**2 166 : & +(SIGT(2)-SIGT(3))**2+(SIGT(3)-SIGT(1))**2 167 : & +6.D0*(SIGT(4)**2+SIGT(5)**2+SIGT(6)**2)) 168 : C====================================================================== 169 : C - Deformation elastique au temps t+dt 170 : C====================================================================== 171 : C Calcul du multiplicateur plastique au temps t 172 : Fit=XJ2t/(SIGMT**2)+2*Q1*fT*COSH((Q2*XI1t)/(2*SIGMT)) 173 : & -(1+(Q3*fT)**2) 174 : IF (Fit.LT.0.) THEN 175 : Fit=0. 176 : ENDIF 177 : xlambdapo=(fT/(1-fT))*(Fit/XKPLAS)**XMPLAS 178 : 179 : C - Calcul des deformations plastique au temps t 180 : C Derivees partielles de fi 181 : DFIDSIG(1)=0.5*(4*SIGT(1)-2*SIGT(2)-2*SIGT(3))/SIGMT**2 182 : & +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt 183 : DFIDSIG(2)=0.5*(-2*SIGT(1)+4*SIGT(2)-2*SIGT(3))/SIGMT**2 184 : & +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt 185 : DFIDSIG(3)=0.5*(-2*SIGT(1)-2*SIGT(2)+4*SIGT(3))/SIGMT**2 186 : & +Q1*Q2*ft*SINH(Q2*XI1t/(2*SIGMt))/SIGMt 187 : DFIDSIG(4)=6/(SIGMt**2)*SIGT(4) 188 : DFIDSIG(5)=6/(SIGMt**2)*SIGT(5) 189 : DFIDSIG(6)=6/(SIGMt**2)*SIGT(6) 190 : C 191 : DO i=1,6 192 : EPSPOPT(i)=xlambdapo*DFIDSIG(i) 193 : ENDDO 194 : TRACEPSPOPT=EPSPOPT(1)+EPSPOPT(2)+EPSPOPT(3) 195 : C - Déformation viscoplastique 196 : DO i=1,6 197 : EPSPTDT(i)=EPSPT(i)+EPSPOPT(i)*DELTAT 198 : ENDDO 199 : TRACEPSPTDT=EPSPTDT(1)+EPSPTDT(2)+EPSPTDT(3) 200 : C - Deformation elastique 201 : DO i=1,6 202 : EPSTDTE(i)=EPSTDT(i)-EPSPTDT(i) 203 : ENDDO 204 : TRACEPSTDTE=EPSTDTE(1)+EPSTDTE(2)+EPSTDTE(3) 205 : C====================================================================== 206 : C - Deformation equivalente elastique au sens de Mazars a t 207 : ARRAY(1,1)=EPSTDTE(1) 208 : ARRAY(2,2)=EPSTDTE(2) 209 : ARRAY(3,3)=EPSTDTE(3) 210 : ARRAY(1,2)=EPSTDTE(4) 211 : ARRAY(2,3)=EPSTDTE(5) 212 : ARRAY(1,3)=EPSTDTE(6) 213 : ARRAY(3,1)=ARRAY(1,3) 214 : ARRAY(3,2)=ARRAY(2,3) 215 : ARRAY(2,1)=ARRAY(1,2) 216 : C 216 : CALL JACOB3(ARRAY,3,EPSIPP,S) 217 : EPSE=MAX( XZERO , EPSIPP(1) )**2 + 218 : 1 MAX( XZERO , EPSIPP(2) )**2 + 219 : 2 MAX( XZERO , EPSIPP(3) )**2 220 : EPSE=DSQRT(EPSE) 221 : C - Nouveau seuil de Mazars 222 : ES=max (ES, EPSE) 223 : C ============================================================== 224 : C - Prediction élastique de la contrainte (pas d'évolution de d) 225 : DO i=1,3 226 : SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo* 227 : & ((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i)-1./3.*TRACEPSPOPT)) 228 : ENDDO 229 : DO i=4,6 230 : SIGPO(i)=2*XGendo*((EPSPO(i)) 231 : & -(EPSPOPT(i))) 232 : ENDDO 233 : C 234 : DO i=1,6 235 : SIGTDTE(i)=SIGT(i)+SIGPO(i)*DELTAT 236 : ENDDO 237 : C ============================================================ 238 : C - Calcul des parametres avec la contrainte predite elastique 239 : C - Calcul des invariants à t 240 : XI1tdte=SIGTDTE(1)+SIGTDTE(2)+SIGTDTE(3) 241 : C 242 : XJ2tdte=0.5*((SIGTDTE(1)-SIGTDTE(2))**2 243 : & +(SIGTDTE(2)-SIGTDTE(3))**2+(SIGTDTE(3)-SIGTDTE(1))**2 244 : & +6*(SIGTDTE(4)**2+SIGTDTE(5)**2+SIGTDTE(6)**2)) 245 : C 246 : C Derivees partielles de fi 247 : DFIDSIG(1)=0.5*(4*SIGTDTE(1) 248 : & -2*SIGTDTE(2)-2*SIGTDTE(3))/SIGMT**2 249 : & +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt 250 : DFIDSIG(2)=0.5*(-2*SIGTDTE(1)+4*SIGTDTE(2) 251 : & -2*SIGTDTE(3))/SIGMT**2 252 : & +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt 253 : DFIDSIG(3)=0.5*(-2*SIGTDTE(1)-2*SIGTDTE(2) 254 : & +4*SIGTDTE(3))/SIGMT**2 255 : & +Q1*Q2*ft*SINH(Q2*XI1tdte/(2*SIGMt))/SIGMt 256 : DFIDSIG(4)=6/(SIGMt**2)*SIGTDTE(4) 257 : DFIDSIG(5)=6/(SIGMt**2)*SIGTDTE(5) 258 : DFIDSIG(6)=6/(SIGMt**2)*SIGTDTE(6) 259 : C 260 : TRACEDFIDSIG=DFIDSIG(1)+DFIDSIG(2)+DFIDSIG(3) 261 : C 262 : DFIDSIGM=-2*XJ2tdte/(SIGMt**3)-Q1*Q2*ft*XI1tdte 263 : & *SINH(Q2*XI1tdte/(2*SIGMt))/(SIGMt**2) 264 : C 265 : DFIDft=2*Q1*COSH(Q2*XI1tdte/(2*SIGMt))-2*ft*Q3*Q3 266 : C=================================================================== 267 : C - Calcul de la surface seuil fi à l'aide de la contrainte prédite 268 : Fi=XJ2tdte/(SIGMT**2)+2*Q1*fT*COSH((Q2*XI1tdte)/(2*SIGMT)) 269 : & -(1+(Q3*fT)**2) 270 : C 271 : C====================================================================== 272 : C========= Organigramme du calcul de la contrainte a t+Dt ============= 273 : C====================================================================== 274 : C 275 : C Test sur le signe de la deformation equivalente 276 : 277 : 278 : IF(EPSE.GT.XZERO) THEN 279 : IF(fi.LT.XZERO)THEN 280 : IF(EPSE.LT.ES) THEN 281 : C ===================================================================== 282 : C - Elasticite Seule 283 : dtdt=dt 284 : ftdt=ft 285 : dttdt=dtt 286 : dctdt=dct 287 : sigmtdt=sigmt 288 : DO i=1,6 289 : sigtdt(i)=sigtdte(i) 290 : ENDDO 291 : ELSE 292 : C ===================================================================== 293 : C - Visco-Endommagement Seul 294 : sigmtdt=sigmt 295 : ftdt=ft 295 : CALL VISCOE(XMAT,EPSIPP,EPSE,DTT,DCT,DT 296 : & ,DELTAT,DTTDT,DCTDT,DTDT) 297 : C 298 : Dpo=(Dtdt-Dt)/DELTAT 299 : C 300 : C - Raideur endommagee et ses derivees 301 : XKendo=(1-Dtdt)*XKporo 302 : XGendo=(1-Dtdt)*XGporo 303 : XKpo=-Dpo*XKporo 304 : XGpo=-Dpo*XGporo 305 : C 306 : C - Calcul de la contrainte 307 : DO i=1,3 308 : SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo 309 : & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i) 310 : & -1./3.*TRACEPSPOPT)) 311 : & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo 312 : & *((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i) 313 : & -1./3.*TRACEPSPTDT)) 314 : ENDDO 315 : DO i=4,6 316 : SIGPO(i)=2*XGendo*((EPSPO(i)) 317 : & -(EPSPOPT(i)))+2*XGpo 318 : & *((EPSTDT(i))-(EPSPTDT(i))) 319 : ENDDO 320 : C 321 : DO i=1,6 322 : SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT 323 : ENDDO 324 : ENDIF 325 : ELSE 326 : IF(EPSE.LT.ES) THEN 327 : C ===================================================================== 328 : C - Plasticite sans endommagement 329 : Dtdt=Dt 330 : dttdt=dtt 331 : dctdt=dct 332 : C - Evolution de la contrainte dans la matrice 333 : EPSPOPM=-xlambdapo*DFIDSIGM 334 : SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM 335 : SIGMTDT=SIGMT+SIGPOM*DELTAT 336 : C - Evolution de la porosite 337 : fpo=XK*((UN-ft)*ft*TRACEPSPOPT) 338 : ftdt=ft+fpo*DELTAT 339 : IF (ftdt.LT.XPETIT) THEN 340 : ftdt=XZERO 341 : ENDIF 342 : C - Raideur endommagee et ses derives 343 : XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT) 344 : XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm)) 345 : XKendo=(1-Dtdt)*XKporo 346 : XGendo=(1-Dtdt)*XGporo 347 : XKpo=(1-Dtdt)*(-4*XKm*XGm*(4*XGm+3*XKm) 348 : & /(4*XGm+3*XKm*fTDT)**2)*fpo 349 : XGpo=(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2) 350 : & /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo 351 : C ============================= 352 : C - Correction de la contrainte 353 : DO i=1,3 354 : SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo 355 : & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i) 356 : & -1./3.*TRACEPSPOPT)) 357 : & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo* 358 : & ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i) 359 : & -1./3.*TRACEPSPTDT)) 360 : ENDDO 361 : DO i=4,6 362 : SIGPO(i)=2*XGendo*((EPSPO(i)) 363 : & -(EPSPOPT(i)))+2*XGpo* 364 : & ((EPSTDT(i))-(EPSPTDT(i))) 365 : ENDDO 366 : C 367 : DO i=1,6 368 : SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT 369 : ENDDO 370 : ELSE 371 : C ===================================================================== 372 : C - Plasticite avec endommagement 373 : C - Evolution de l'endommagement 373 : CALL VISCOE(XMAT,EPSIPP,EPSE,DTT,DCT,DT 374 : & ,DELTAT,DTTDT,DCTDT,DTDT) 375 : C 376 : Dpo=(Dtdt-Dt)/DELTAT 377 : IF (Dpo.LE.0.) THEN 378 : Dtdt=Dt 379 : Dpo=0. 380 : ENDIF 381 : C - Evolution de la contrainte dans la matrice 382 : EPSPOPM=-xlambdapo*DFIDSIGM 383 : SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM 384 : SIGMTDT=SIGMT+SIGPOM*DELTAT 385 : C - Evolution de la porosite 386 : fpo=XK*((1-ft)*ft*TRACEPSPOPT) 387 : ftdt=ft+fpo*DELTAT 388 : IF (ftdt.LT.0.) THEN 389 : ftdt=0. 390 : ENDIF 391 : C - Raideur endommagee et ses derives 392 : XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT) 393 : XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm)) 394 : XKendo=(1-dtdt)*XKporo 395 : XGendo=(1-dtdt)*XGporo 396 : XKpo=-Dpo*XKporo+(1-Dtdt)* 397 : & (-4*XKm*XGm*(4*XGm+3*XKm)/(4*XGm+3*XKm*fTDT)**2)*fpo 398 : XGpo=-Dpo*XGporo 399 : & +(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2) 400 : & /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo 401 : C ============================= 402 : C - Correction de la contrainte 403 : DO i=1,3 404 : SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo 405 : & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i) 406 : & -1./3.*TRACEPSPOPT)) 407 : & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo* 408 : & ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i) 409 : & -1./3.*TRACEPSPTDT)) 410 : ENDDO 411 : DO i=4,6 412 : SIGPO(i)=2*XGendo*((EPSPO(i)) 413 : & -(EPSPOPT(i)))+2*XGpo* 414 : & ((EPSTDT(i))-(EPSPTDT(i))) 415 : ENDDO 416 : C 417 : DO i=1,6 418 : SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT 419 : ENDDO 420 : ENDIF 421 : ENDIF 422 : ELSE 423 : IF(fi.LT.XZERO)THEN 424 : C ===================================================================== 425 : C - Elasticite Seule 426 : Dtdt=Dt 427 : ftdt=ft 428 : dttdt=dtt 429 : dctdt=dct 430 : sigmtdt=sigmt 431 : DO i=1,6 432 : sigtdt(i)=sigtdte(i) 433 : ENDDO 434 : ELSE 435 : C ===================================================================== 436 : C - Plasticite sans endommagement 437 : Dtdt=Dt 438 : dttdt=dtt 439 : dctdt=dct 440 : C - Evolution de la contrainte dans la matrice 441 : EPSPOPM=-xlambdapo*DFIDSIGM 442 : SIGPOM=(Etan*Eini/(Eini-Etan))*EPSPOPM 443 : SIGMTDT=SIGMT+SIGPOM*DELTAT 444 : C - Evolution de la porosite 445 : fpo=XK*((1-ft)*ft*TRACEPSPOPT) 446 : ftdt=ft+fpo*DELTAT 447 : IF (ftdt.LT.0.) THEN 448 : ftdt=0. 449 : ENDIF 450 : C - Raideur endommagee et ses derives 451 : XKporo=4*XKm*XGm*(1-fTDT)/(4*XGm+3*XKm*fTDT) 452 : XGporo=XGm*(1-fTDT)/(1+fTDT*(6*XKm+12*XGm)/(9*XKm+8*XGm)) 453 : XKendo=(1-Dtdt)*XKporo 454 : XGendo=(1-Dtdt)*XGporo 455 : XKpo=(1-Dtdt)*(-4*XKm*XGm*(4*XGm+3*XKm) 456 : & /(4*XGm+3*XKm*fTDT)**2)*fpo 457 : XGpo=(1-Dtdt)*(-5*XGm*(32*XGm**2+60*XGm*XKm+27*XKm**2) 458 : & /(4*(2+3*fTDT)*XGm+3*(3+2*fTDT)*XKm)**2)*fpo 459 : C ============================= 460 : C - Correction de la contrainte 461 : DO i=1,3 462 : SIGPO(i)=XKendo*(TRACEPSPO-TRACEPSPOPT)+2*XGendo 463 : & *((EPSPO(i)-1./3.*TRACEPSPO)-(EPSPOPT(i) 464 : & -1./3.*TRACEPSPOPT)) 465 : & +XKpo*(TRACEPSTDT-TRACEPSPTDT)+2*XGpo* 466 : & ((EPSTDT(i)-1./3.*TRACEPSTDT)-(EPSPTDT(i) 467 : & -1./3.*TRACEPSPTDT)) 468 : ENDDO 469 : DO i=4,6 470 : SIGPO(i)=2*XGendo*((EPSPO(i)) 471 : & -(EPSPOPT(i)))+2*XGpo* 472 : & ((EPSTDT(i))-(EPSPTDT(i))) 473 : ENDDO 474 : C 475 : DO i=1,6 476 : SIGTDT(i)=SIGT(i)+SIGPO(i)*DELTAT 477 : ENDDO 478 : ENDIF 479 : ENDIF 480 : C 481 : C====================================================================== 482 : C========= Fin de l'organigramme du calcul de la contrainte a t+dt ==== 483 : C====================================================================== 484 : C 485 : C ===================================================================== 486 : C - Evolution des variables historiques 487 : VAR0(1)=DTDT 488 : VAR0(2)=fTDT-D0 489 : VAR0(3)=EPSPTDT(1) 490 : VAR0(4)=EPSPTDT(2) 491 : VAR0(5)=EPSPTDT(3) 492 : VAR0(6)=EPSPTDT(4) 493 : VAR0(7)=EPSPTDT(5) 494 : VAR0(8)=EPSPTDT(6) 495 : VAR0(9)=SIGMTDT-SIGM0 496 : VAR0(10)=dttdt 497 : VAR0(11)=dctdt 498 : VAR0(12)=ES-ED0 499 : C ===================================================================== 500 : SIGT(1)=SIGTDT(1) 501 : SIGT(2)=SIGTDT(2) 502 : SIGT(3)=SIGTDT(3) 503 : SIGT(4)=SIGTDT(4) 504 : SIGT(5)=SIGTDT(5) 505 : SIGT(6)=SIGTDT(6) 506 : C ===================================================================== 507 : EPST(1)=EPSTDT(1) 508 : EPST(2)=EPSTDT(2) 509 : EPST(3)=EPSTDT(3) 510 : EPST(4)=EPSTDT(4) 511 : EPST(5)=EPSTDT(5) 512 : EPST(6)=EPSTDT(6) 513 : C ===================================================================== 514 : C - Compteur du Subcycling 515 : l=l+1 516 : ENDDO 517 : C ===================================================================== 518 : C ========Fin du Subcycling================= 519 : C ===================================================================== 520 : C - Evolution des variables historiques 521 : VARF(1)=DTDT 522 : VARF(2)=fTDT-D0 523 : VARF(3)=EPSPTDT(1) 524 : VARF(4)=EPSPTDT(2) 525 : VARF(5)=EPSPTDT(3) 526 : VARF(6)=EPSPTDT(4) 527 : VARF(7)=EPSPTDT(5) 528 : VARF(8)=EPSPTDT(6) 529 : VARF(9)=SIGMTDT-SIGM0 530 : VARF(10)=dttdt 531 : VARF(11)=dctdt 532 : VARF(12)=ES-ED0 533 : C====================================================================== 534 : C - Mise a jour des contraintes 535 : SIGF(1)=SIGTDT(1) 536 : SIGF(2)=SIGTDT(2) 537 : SIGF(3)=SIGTDT(3) 538 : SIGF(4)=SIGTDT(4) 539 : SIGF(5)=SIGTDT(5) 540 : SIGF(6)=SIGTDT(6) 541 : C 542 : END 543 : 544 : 545 : 546 : 547 :
© Cast3M 2003 - All rights reserved.
Disclaimer