1 : C PAKGON SOURCE CHAT 05/01/13 02:08:47 5004 2 : SUBROUTINE PAKGON(XMAT,SIG0,VAR0,SIGF,DEFF,VARF, 3 : . DEPST,KERRE,IB,IGAU,PCT1, 4 : . SUCF0,RSIG0) 5 : * 6 : *============================================================= 7 : * Modèle d'argile de PAKZAD: surface de gonflement 8 : *============================================================= 9 : * 10 : * Entrées 11 : *------------------------------------------------------------- 12 : * 13 : * XMAT : caractéristiques matériau 14 : * SIG0 : contraintes initiales TEST 15 : * VAR0 : variables internes initiales TEST 16 : * DEPST : incrément des déformations totales 17 : * IB : numéro de l'élément courrant 18 : * IGAU : numéro du point de Gauss courrant 19 : * PCT1 : pression de consolidation "fictive" de l'état saturé 20 : * SUCF0 : succion finale 21 : * RSIG0 : contraintes au début du pas 22 : * 23 : *------------------------------------------------------------- 24 : * 25 : * 26 : * Sorties 27 : *------------------------------------------------------------- 28 : * 29 : * SIGF : contraintes finales 30 : * DEFF : incrémént de déformation plastique final 31 : * VARF : variables internes finales 32 : * KERRE : message d'erreur 33 : * 34 : *------------------------------------------------------------- 35 : * ICI ON PROGRAMME EN FORTRAN PUR 36 : * 37 : IMPLICIT INTEGER(I-N) 38 : IMPLICIT REAL*8(A-H,O-Z) 39 : * 40 : DIMENSION XMAT(*),SIG0(*),VAR0(*),DEPST(*) 41 : DIMENSION SIGF(*),DEFF(*),VARF(*),RSIG0(*) 42 : DIMENSION SIGT(6),DEVT(6),DEV(6),DEV0(6),DEFP(6) 43 : * 44 : * 45 : * Données du materiau 46 : *=========================================================== 47 : * 48 : YOUN0=XMAT(1) 49 : XNU0=XMAT(2) 50 : XN0=XMAT(5) 51 : XKA0=XMAT(6) 52 : XGA0=XMAT(7) 53 : XPA0=XMAT(8) 54 : XPC0=XMAT(9) 55 : XM0=XMAT(10) 56 : XBET1=XMAT(11) 57 : XA0=XMAT(12) 58 : XPCR0=XMAT(13) 59 : XS0=XMAT(14) 60 : XM1=XMAT(15) 61 : XM2=XMAT(16) 62 : XM3=XMAT(17) 63 : XBET2=XMAT(18) 64 : XTAU0=XMAT(19) 65 : dif0=ABS(XN0-1.D0) 66 : * 67 : *---> Initialisation du rayon initial de l'ellipse 68 : * 69 : IF (VAR0(5).LE.1.D-10) VAR0(5)=XA0 70 : IF (VAR0(6).LE.1.D-10) VAR0(6)=XA0 71 : * 72 : *---> Charge ou décharge ? 73 : * . NCHAR0 = 0 : décharge 74 : * . NCHAR0 = 1 : recharge 75 : * 76 : TRELT=DEPST(1)+DEPST(2)+DEPST(3) 77 : NCHAR0=0 78 : IF (TRELT.LE.0.D0) NCHAR0=1 79 : * 80 : *---> Contraintes test, incrément de déformation plastique test, 81 : * déviateur initial et pression initiale 82 : * 83 : PRES00=-1.D0/3.D0*(RSIG0(1)+RSIG0(2)+RSIG0(3)) 84 : DO 05 I=1,6 85 : A=0.D0 86 : IF (I.LE.3) A=1.D0 87 : SIGT(I)=SIG0(I) 88 : DEFP(I)=0.D0 89 : DEVT(I)=RSIG0(I)+PRES00*A 90 : 05 CONTINUE 91 : * 92 : *---> Cas où XN0 n'est pas égal à 1 93 : * Les constantes de raideur sont prises au début du pas 94 : * 95 : IF (dif0.GT.1.D-5) THEN 96 : XKA1=XKA0*XPA0*((PRES00/XPA0)**XN0) 97 : XGA1=XGA0*XPA0*((PRES00/XPA0)**XN0) 98 : ENDIF 99 : * 100 : *---> Contrainte équivalente initiale 101 : * 102 : Stest=PROCON(DEVT,DEVT,6) 103 : IF (Stest.LE.(1.D-10*YOUN0)) Stest=0.D0 104 : St0=(3.D0*Stest/2.D0)**(.5D0) 105 : * 106 : *---> Calcul . du terme de contrainte équivalente 107 : * . du centre de la surface 108 : * . du rayon 109 : * 110 : PRES1=-1.D0/3.D0*(SIGT(1)+SIGT(2)+SIGT(3)) 111 : IF (NCHAR0.EQ.1) THEN 112 : * 113 : * On charge 114 : *----------------------------------------------------- 115 : * 116 : * Centre test et au début du pas 117 : * 118 : DO 10 I=1,6 119 : A=0.D0 120 : IF (I.LE.3) A=1.D0 121 : DEVT(I)=SIGT(I)+PRES1*A 122 : DEV(I)=(DEVT(I)-VAR0(14+I))/2.D0 123 : DEV0(I)=(RSIG0(I)+PRES00*A-VAR0(14+I))/2.D0 124 : 10 CONTINUE 125 : APT=(PRES1+VAR0(8))/2.D0 126 : APT0=(PRES00+VAR0(8))/2.D0 127 : * 128 : * Rayon test et au début du pas 129 : * 130 : AIS=VAR0(6) 131 : AIS0=AIS 132 : * 133 : * Contrainte équivalente du point fixe 134 : * 135 : Stest0=PROCON(VAR0(15),VAR0(15),6) 136 : IF (Stest0.LE.(1.D-10*YOUN0)) Stest0=0.D0 137 : Stest0=(3.D0*Stest0/2.D0)**(.5D0) 138 : ELSE 139 : * 140 : * On décharge 141 : *----------------------------------------------------- 142 : * 143 : * Centre test et au début du pas 144 : * 145 : DO 11 I=1,6 146 : A=0.D0 147 : IF (I.LE.3) A=1.D0 148 : DEVT(I)=SIGT(I)+PRES1*A 149 : DEV(I)=(DEVT(I)-VAR0(8+I))/2.D0 150 : DEV0(I)=(RSIG0(I)+PRES00*A-VAR0(8+I))/2.D0 151 : 11 CONTINUE 152 : APT=(PRES1+VAR0(7))/2.D0 153 : APT0=(PRES00+VAR0(7))/2.D0 154 : * 155 : * Rayon test et au début du pas 156 : * 157 : AIS=VAR0(5) 158 : AIS0=AIS 159 : * 160 : * Contrainte équivalente du point fixe 161 : * 162 : Stest0=PROCON(VAR0(9),VAR0(9),6) 163 : IF (Stest0.LE.(1.D-10*YOUN0)) Stest0=0.D0 164 : Stest0=(3.D0*Stest0/2.D0)**(.5D0) 165 : ENDIF 166 : * 167 : *---> Terme déviatoire de la distance au centre de la surface test 168 : * 169 : AJ2=PROCON(DEV,DEV,6) 170 : IF (AJ2.LE.(1.D-10*YOUN0)) AJ2=0.D0 171 : AJ2=(3.D0/2.D0*AJ2)**(.5D0) 172 : * 173 : *---> Terme déviatoire de la distance au centre de la surface au 174 : * début du pas 175 : * 176 : AJ20=PROCON(DEV0,DEV0,6) 177 : IF (AJ20.LE.(1.D-10*YOUN0)) AJ20=0.D0 178 : AJ20=(3.D0/2.D0*AJ20)**(.5D0) 179 : * 180 : *---> Pressions critiques au début du pas 181 : * 182 : APT10=2.D0*APT0-PRES00 183 : PC2=(St0*St0+XM0*XM0*PRES00*PRES00)/(XM0*XM0*PRES00) 184 : PCD1=(Stest0*Stest0+XM0*XM0*APT10*APT10)/(XM0*XM0*APT10) 185 : * 186 : *---> Initialisation de la pression initiale du "CENTRE" 187 : * 188 : IF (NCHAR0.EQ.1) THEN 189 : IF (ABS(VAR0(22)).LE.1.D-10) THEN 190 : AP1=APT 191 : AP2=0.D0 192 : ELSE 193 : AP1=VAR0(22) 194 : AP2=AP1 195 : ENDIF 196 : ELSE 197 : IF (ABS(VAR0(21)).LE.1.D-10) THEN 198 : AP1=APT 199 : AP2=0.D0 200 : ELSE 201 : AP1=VAR0(21) 202 : AP2=AP1 203 : ENDIF 204 : ENDIF 205 : * 206 : *---> Terme d'écrouissage BETA2 au début du pas 207 : * Terme ALP0 non associatif de l'écoulemment 208 : * 209 : IF (NCHAR0.EQ.1) THEN 210 : BETA2=(XPCR0/PC2+1.D0)*PRES00*2.D0*VAR0(4) 211 : BETA2=XBET1*BETA2/(3.D0*PCD1*PC2) 212 : ALP0=(1.D0-(ABS(APT10-AP1)/XA0)) 213 : ALP0=ALP0*(1.D0-1.5D0*PC2/VAR0(4)) 214 : ELSE 215 : BETA2=XBET1*(XPCR0/PCD1+1.D0)*PC2/PRES00 216 : ALP0=(1.D0-(ABS(APT10-AP1)/XA0)) 217 : ALP0=ALP0*(1.5D0*PC2/PCD1-.5D0) 218 : ENDIF 219 : BETA2=BETA2*EXP(XTAU0*SUCF0) 220 : ALP0=-1.D0*ALP0 221 : * 222 : *---> Termes d'accroissement du rayon de l'ellipse 223 : * 224 : XB1=PRES00*BETA2*(PRES00-APT0)/(2.D0*AIS0) 225 : XB2=PRES00*BETA2*AJ20/(2.D0*AIS0) 226 : * 227 : *---> Surface de charge de gonflement 228 : * 229 : PHIT=AJ2*AJ2+XM0*XM0*(PRES1-APT+AIS)*(PRES1-APT-AIS) 230 : * 231 : * 232 : * 233 : * Vérification du critère de plasticité 234 : *========================================================= 235 : * 236 : *---> Erreur admise 237 : * 238 : ERR0=1.D-7*ABS(PHIT) 239 : ERR0=MAX(ERR0,1.D-10*YOUN0) 240 : * 241 : *---> Critère de plasticité 242 : * 243 : PHI0=PHIT 244 : iter0=0 245 : * 246 : 99 IF (PHI0.LE.ERR0) THEN 247 : * 248 : * On est élastique par rapport à la surface de gonflement 249 : *========================================================= 250 : * 251 : IF (NCHAR0.EQ.1) THEN 252 : * 253 : * On charge 254 : * 255 : DO 20 I=1,6 256 : SIGF(I)=SIGT(I) 257 : DEFF(I)=DEFP(I) 258 : VARF(8+I)=DEVT(I) 259 : VARF(14+I)=VAR0(14+I) 260 : 20 CONTINUE 261 : VARF(5)=XA0 262 : VARF(6)=AIS 263 : VARF(7)=PRES1 264 : VARF(8)=VAR0(8) 265 : VARF(21)=0.D0 266 : VARF(22)=AP2 267 : ELSE 268 : * 269 : * On décharge 270 : * 271 : DO 21 I=1,6 272 : SIGF(I)=SIGT(I) 273 : DEFF(I)=DEFP(I) 274 : VARF(8+I)=VAR0(8+I) 275 : VARF(14+I)=DEVT(I) 276 : 21 CONTINUE 277 : VARF(5)=AIS 278 : VARF(6)=XA0 279 : VARF(7)=VAR0(7) 280 : VARF(8)=PRES1 281 : VARF(21)=AP2 282 : VARF(22)=0.D0 283 : ENDIF 284 : * 285 : ELSE 286 : * 287 : * On est plastique par rapport à la surface de gonflement 288 : *========================================================= 289 : * 290 : *---> Pression initiale du "CENTRE" 291 : * 292 : IF (NCHAR0.EQ.1) THEN 293 : IF (ABS(VAR0(22)).LE.1.D-10) THEN 294 : AP1=APT 295 : AP2=AP1 296 : ENDIF 297 : ELSE 298 : IF (ABS(VAR0(21)).LE.1.D-10) THEN 299 : AP1=APT 300 : AP2=AP1 301 : ENDIF 302 : ENDIF 303 : * 304 : *---> Condition de consistance 305 : * 306 : IF (dif0.LE.1.D-5) THEN 307 : * 308 : * Cas où XN0=1 309 : *------------------------------------------------------------- 310 : * 311 : tr0=2.D0*XM0*XM0*(PRES1-APT) 312 : DEQ=ABS(2.D0*AJ2-ALP0*tr0) 313 : DF0=3.D0*XGA0*PRES00*AJ2*DEQ 314 : DF1=XKA0*PRES1*tr0*tr0/2.D0 315 : DF2=(XB1*tr0+XB2*DEQ)*2.D0*XM0*XM0*AIS 316 : * 317 : denom0=DF0+DF1+DF2 318 : * 319 : dlam0=PHIT/denom0 320 : * 321 : *---> Mise à jour des contraintes 322 : * 323 : PRES1=PRES1*EXP(-1.D0*XKA0*tr0*dlam0) 324 : DO 30 I=1,6 325 : A=0.D0 326 : IF (I.LE.3) A=1.D0 327 : IF (AJ2.GT.(YOUN0*1.D-10)) THEN 328 : DEVT(I)=DEVT(I)-3.D0*XGA0*PRES00*DEV(I)/AJ2*DEQ*dlam0 329 : ELSE 330 : DEVT(I)=DEVT(I)-6.D0*XGA0*PRES00*DEV(I)*dlam0 331 : ENDIF 332 : SIGT(I)=DEVT(I)-PRES1*A 333 : 30 CONTINUE 334 : * 335 : *---> Mise à jour du rayon de l'ellipse 336 : * 337 : AIS=AIS+tr0*XB1*dlam0+XB2*DEQ*dlam0 338 : * 339 : *---> Mise à jour du centre de l'ellipse 340 : * 341 : IF (NCHAR0.EQ.1) THEN 342 : APT=(PRES1+VAR0(8))/2.D0 343 : DO 35 I=1,6 344 : DEV(I)=(DEVT(I)-VAR0(14+I))/2.D0 345 : 35 CONTINUE 346 : ELSE 347 : APT=(PRES1+VAR0(7))/2.D0 348 : DO 36 I=1,6 349 : DEV(I)=(DEVT(I)-VAR0(8+I))/2.D0 350 : 36 CONTINUE 351 : ENDIF 352 : AJ2=PROCON(DEV,DEV,6) 353 : IF (AJ2.LE.(1.D-10*YOUN0)) AJ2=0.D0 354 : AJ2=(3.D0/2.D0*AJ2)**(.5D0) 355 : * 356 : *---> Mise à jour des déformations élastiques et plastiques 357 : * 358 : TRELT=-1.D0*LOG(PRES1/PRES00)/XKA0 359 : DO 40 I=1,6 360 : A=0.D0 361 : IF (I.LE.3) A=1.D0 362 : DEFP(I)=(DEVT(I)-(RSIG0(I)+PRES00*A)) 363 : . /(2.D0*XGA0*PRES00) 364 : DEFP(I)=DEPST(I)-DEFP(I)-TRELT*A/3.D0 365 : 40 CONTINUE 366 : * 367 : ELSE 368 : * 369 : * Autres cas 370 : *------------------------------------------------------------- 371 : * 372 : tr0=2.D0*XM0*XM0*(PRES1-APT) 373 : DEQ=ABS(2.D0*AJ2-ALP0*tr0) 374 : DF0=3.D0*XGA1*AJ2*DEQ 375 : DF1=XKA1*tr0*tr0/2.D0 376 : DF2=(XB1*tr0+XB2*DEQ)*2.D0*XM0*XM0*AIS 377 : * 378 : denom0=DF0+DF1+DF2 379 : * 380 : dlam0=PHIT/denom0 381 : * 382 : *---> Mise à jour des contraintes 383 : * 384 : PRES1=PRES1-1.D0*XKA1*tr0*dlam0 385 : DO 31 I=1,6 386 : A=0.D0 387 : IF (I.LE.3) A=1.D0 388 : IF (AJ2.GT.(YOUN0*1.D-10)) THEN 389 : DEVT(I)=DEVT(I)-3.D0*XGA1*DEV(I)/AJ2*DEQ*dlam0 390 : ELSE 391 : DEVT(I)=DEVT(I)-6.D0*XGA1*DEV(I)*dlam0 392 : ENDIF 393 : SIGT(I)=DEVT(I)-PRES1*A 394 : 31 CONTINUE 395 : * 396 : *---> Mise à jour du rayon de l'ellipse 397 : * 398 : AIS=AIS+tr0*XB1*dlam0+XB2*DEQ*dlam0 399 : * 400 : *---> Mise à jour du centre de l'ellipse 401 : * 402 : IF (NCHAR0.EQ.1) THEN 403 : APT=(PRES1+VAR0(8))/2.D0 404 : DO 350 I=1,6 405 : DEV(I)=(DEVT(I)-VAR0(14+I))/2.D0 406 : 350 CONTINUE 407 : ELSE 408 : APT=(PRES1+VAR0(7))/2.D0 409 : DO 360 I=1,6 410 : DEV(I)=(DEVT(I)-VAR0(8+I))/2.D0 411 : 360 CONTINUE 412 : ENDIF 413 : AJ2=PROCON(DEV,DEV,6) 414 : IF (AJ2.LE.(1.D-10*YOUN0)) AJ2=0.D0 415 : AJ2=(3.D0/2.D0*AJ2)**(.5D0) 416 : * 417 : *---> Mise à jour des déformations élastiques et plastiques 418 : * 419 : TRELT=-1.D0*(PRES1-PRES00)/XKA1 420 : DO 41 I=1,6 421 : A=0.D0 422 : IF (I.LE.3) A=1.D0 423 : DEFP(I)=(DEVT(I)-(RSIG0(I)+PRES00*A)) 424 : . /(2.D0*XGA1) 425 : DEFP(I)=DEPST(I)-DEFP(I)-TRELT*A/3.D0 426 : 41 CONTINUE 427 : * 428 : * Fin du cas sur XN0 429 : *------------------------------------------------------------- 430 : * 431 : ENDIF 432 : * 433 : *---> Surface de charge de gonflement 434 : * 435 : PHIT=AJ2*AJ2+XM0*XM0*(PRES1-APT+AIS)*(PRES1-APT-AIS) 436 : PHI0=ABS(PHIT) 437 : * 438 : * 439 : *---> Test de convergence ou itération suivante 440 : * 441 : iter0=iter0+1 442 : IF (iter0.LT.200) THEN 443 : GOTO 99 444 : ELSE 445 : PHI0=0.D0 446 : KERRE=460 447 : GOTO 99 448 : ENDIF 449 : * 450 : * 451 : ENDIF 452 : * 453 : RETURN 454 : END 455 : 456 : 457 :
© Cast3M 2003 - All rights reserved.
Disclaimer