1 : C ALON2 SOURCE CHAT 05/01/12 21:20:54 5004 2 : SUBROUTINE ALON2(Pat0,XMAT,DEPM,RSIG0,VAR0, 3 : . RSIGF,VARF,RDEFP,KERRE,IB,IGAU) 4 : 5 : * 6 : *============================================================ 7 : * 8 : * Modèle d'argile partiellement saturé d'ALONSO 9 : * 10 : *============================================================ 11 : * 12 : * ENTREES 13 : *----------------------------- 14 : * 15 : * Pat0 : pression atmosphérique 16 : * XMAT : champ des caractéristiques materiau 17 : * DEPM : incrément des déformations plastiques 18 : * RSIG0 : contraintes initiales 19 : * VAR0 : variables internes initiales 20 : * IB : numéro de l'élément 21 : * IGAU : numéro du point de Gauss 22 : * 23 : * SORTIES 24 : *----------------------------- 25 : * 26 : * RSIGF : contraintes finales 27 : * VARF : variables internes finales 28 : * RDEFP : increment de déformation plastique au cours du pas 29 : * KERRE : message d'erreur 30 : * 31 : *============================================================ 32 : * ICI IL FAUT PROGRAMMER EN FORTRAN PUR 33 : *============================================================= 34 : * 35 : IMPLICIT INTEGER(I-N) 36 : IMPLICIT REAL*8(A-H,O-Z) 37 : * 38 : DIMENSION VAR0(*),VARF(*),XMAT(*) 39 : DIMENSION RSIG0(*),RSIGF(*),RDEFP(*) 40 : DIMENSION DEPM(*) 41 : DIMENSION DEPT(6),DEVT(6),SIGT(6),DEFT(6),DEPM0(6) 42 : * 43 : * 44 : * Données du materiau 45 : *=========================================================== 46 : * 47 : YOUN0=XMAT(1) 48 : XNU0=XMAT(2) 49 : XKS0=XMAT(5) 50 : XLAM0=XMAT(6) 51 : XM0=XMAT(7) 52 : XKK0=XMAT(8) 53 : PC00=XMAT(9) 54 : P0=XMAT(10) 55 : XLAM1=XMAT(11) 56 : T0=XMAT(12) 57 : TAU0=XMAT(13) 58 : XG0=XMAT(14) 59 : XK0=XMAT(15) 60 : XE0=XMAT(16) 61 : * 62 : * Quelques initialisations 63 : *============================================================= 64 : * 65 : *---> Terme GAMA0 intervenant dans l'écoulement non associatif 66 : * 67 : GAMA0=XM0*(XM0-9.D0)*(XM0-3.D0)*XLAM1 68 : GAMA0=GAMA0/(9.D0*(6.D0-XM0)*(XLAM1-XK0)) 69 : * 70 : *---> Initialisations des sous incrémentations 71 : * 72 : recal0=-100.D0 73 : iter2=0 74 : nmax0=1 75 : 93 DO 100 I=1,6 76 : DEPM0(I)=DEPM(I)/nmax0 77 : SIGT(I)=RSIG0(I) 78 : DEFT(I)=0.D0 79 : 100 CONTINUE 80 : PRES00=-1.D0*(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0 81 : * 82 : *---> Variables internes initiales 83 : * . pression de consolidation à l'état saturé 84 : * . limite élastique en succion 85 : * . succions initiale et finale 86 : * La succion finale est déstinée à etre calculée au préalable 87 : * dans un opérateur particulier indépendant de ECOULE 88 : * 89 : VOLU0=1.D0+XE0 90 : PC0=VAR0(2) 91 : IF (PC0.LE.1.D-30) PC0=PC00 92 : SLIM0=VAR0(3) 93 : SUCI0=VAR0(4) 94 : * 95 : * MODIF POUR CAS SATURE 96 : * 97 : SUCAUX=VARF(4) 98 : IF(SUCI0.LT.0.D0) THEN 99 : SUCI0=0.D0 100 : ENDIF 101 : IF(SUCAUX.LT.0.D0) THEN 102 : SUCAUX=0.D0 103 : ENDIF 104 : * 105 : PC1=PC0 106 : SLIM1=SLIM0 107 : * 108 : *---> Succion a la fin de la sous incrementation 109 : * 110 : 95 iter2=iter2+1 111 : SUCF0=iter2*(SUCAUX-SUCI0)/nmax0+SUCI0 112 : * 113 : * MODIF POUR CAS SATURE 114 : * 115 : IF(SUCF0.LT.0.D0) THEN 116 : SUCF0=0.D0 117 : ENDIF 118 : * 119 : *---> Pente élasto-plastique (courbe de consolidation) 120 : * 121 : XLAM2=XLAM1*((1.D0-T0)*EXP(-1.D0*TAU0*SUCF0)+T0) 122 : * 123 : *---> Pression de saturation 124 : * 125 : PS0=XKK0*SUCF0 126 : * 127 : * 128 : * Calcul élastique test 129 : *================================================================ 130 : * 131 : *---> Déformations élastiques test et déformation plastique 132 : * 133 : DO 11 I=1,6 134 : A=0.D0 135 : IF (I.LE.3) A=1.D0 136 : DEPT(I)=DEPM0(I) 137 : 11 CONTINUE 138 : * 139 : *---> Contraintes 140 : * 141 : PRES0=-1.D0*(SIGT(1)+SIGT(2)+SIGT(3))/3.D0 142 : TREPS0=DEPT(1)+DEPT(2)+DEPT(3) 143 : PRES1=PRES0*EXP(-1.D0*VOLU0/XK0*TREPS0) 144 : DO 12 I=1,6 145 : A=0.D0 146 : IF (I.LE.3) A=1.D0 147 : DEVT(I)=SIGT(I)+PRES0*A 148 : DEVT(I)=DEVT(I)+2.D0*XG0*(DEPT(I)-TREPS0/3.D0*A) 149 : SIGT(I)=DEVT(I)-PRES1*A 150 : 12 CONTINUE 151 : * 152 : *---> Contrainte équivalente 153 : * 154 : STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0 155 : IF (STEST2.LE.1.D-10) STEST2=0.D0 156 : Stest=(STEST2)**(0.5D0) 157 : * 158 : *---> Pression de consolidation 159 : * . à l'état saturé 160 : * . à la succion finale 161 : * 162 : PCS1=P0*((PC1/P0)**((XLAM1-XK0)/(XLAM2-XK0))) 163 : * 164 : *---> Fonction de charge 165 : * 166 : PHIT=Stest*Stest+XM0*XM0*(PRES1+PS0)*(PRES1-PCS1) 167 : * 168 : * 169 : * 170 : * Vérification du critère de plasticité 171 : *========================================================= 172 : * 173 : *---> Erreur admise 174 : * 175 : ERR0=1.D-7*ABS(PHIT) 176 : ERR0=MAX(ERR0,1.D-10*YOUN0) 177 : * 178 : *---> Critère de plasticité 179 : * 180 : PHI0=PHIT 181 : iter0=0 182 : * 183 : 99 IF (PHI0.LE.ERR0) THEN 184 : * 185 : * On est élastique 186 : *========================================================= 187 : * 188 : * 189 : DO 400 I=1,6 190 : RSIGF(I)=SIGT(I) 191 : RDEFP(I)=DEFT(I) 192 : A=0.D0 193 : IF (I.LE.3) A=1.D0 194 : DEVT(I)=DEFT(I)-(DEFT(1)+DEFT(2)+DEFT(3))/3.D0*A 195 : 400 CONTINUE 196 : EPSPT=PROCON(DEVT,DEVT,6) 197 : IF (EPSPT.LE.1.D-20) EPSPT=0.D0 198 : EPSPT=(2.D0/3.D0*EPSPT)**.5D0 199 : VARF(1)=VAR0(1)+EPSPT 200 : VARF(2)=PC1 201 : VARF(3)=MAX(SLIM1,0.D0) 202 : * 203 : ELSE 204 : * 205 : * On plastifie 206 : *========================================================= 207 : * 208 : *---> Condition de consistance: calcul du paramètre plastique 209 : * 210 : DF1=XM0*XM0*(2.D0*PRES1+PS0-PCS1) 211 : DF2=XM0*XM0*(PRES1+PS0)*(XLAM1-XK0)/(XLAM2-XK0) 212 : ***** DF2=DF2*((PC1/P0)**((XLAM1-XLAM2)/(XLAM1-XK0))) 213 : ***** AM 10/12/03 214 : DF2=DF2*((PC1/P0)**((XLAM1-XLAM2)/(XLAM2-XK0))) 215 : denom0=12.D0*Stest*Stest*XG0*GAMA0 216 : denom0=denom0+DF1*PRES1*VOLU0/XK0*DF1 217 : denom0=denom0+DF2*PC1*VOLU0/(XLAM1-XK0)*DF1 218 : * 219 : dlam0=PHIT/denom0 220 : * 221 : * 222 : *---> Critères de sous incrémentation 223 : * 224 : IF (recal0.LT.0.D0) THEN 225 : * 226 : xmax0=2000.D0 227 : xmax1=0.D0 228 : xmax2=0.D0 229 : xnum0=100.D0 230 : * 231 : cri0=ABS(VOLU0/XK0*DF1*dlam0) 232 : IF (cri0.GT.1.D-2) THEN 233 : xmax1=xnum0*cri0+10.D0 234 : ENDIF 235 : * 236 : cri1=ABS(VOLU0/(XK0-XLAM1)*DF1*dlam0) 237 : IF (cri1.GT.1.D-2) THEN 238 : xmax2=xnum0*cri1+10.D0 239 : ENDIF 240 : * 241 : xmax00=xmax2 242 : IF (xmax1.GE.xmax2) xmax00=xmax1 243 : IF (xmax00.GE.xmax0) xmax00=xmax0 244 : nmax00=NINT(xmax00) 245 : IF (nmax00.GT.nmax0) nmax0=nmax00 246 : recal0=100.D0 247 : * 248 : IF (nmax0.GT.1) THEN 249 : iter2=0 250 : GOTO 93 251 : ENDIF 252 : * 253 : ENDIF 254 : * 255 : *---> Mise à jour des contraintes 256 : * 257 : PRES1=PRES1*EXP(-1.D0*VOLU0/XK0*DF1*dlam0) 258 : DO 10 I=1,6 259 : A=0.D0 260 : IF (I.LE.3) A=1.D0 261 : DEVT(I)=DEVT(I)-6.D0*XG0*DEVT(I)*GAMA0*dlam0 262 : SIGT(I)=DEVT(I)-PRES1*A 263 : 10 CONTINUE 264 : * 265 : *---> Mise à jour de la pression de consolidation 266 : * . dans le domaine saturé 267 : * . à la succion finale 268 : * 269 : PC1=PC1*EXP(VOLU0/(XLAM1-XK0)*DF1*dlam0) 270 : PCS1=P0*((PC1/P0)**((XLAM1-XK0)/(XLAM2-XK0))) 271 : * 272 : *---> Déformation élastiques et plastiques 273 : * 274 : TREPS0=-1.D0*XK0/VOLU0*LOG(PRES1/PRES00) 275 : DO 20 I=1,6 276 : A=0.D0 277 : IF (I.LE.3) A=1.D0 278 : DEPT(I)=(DEVT(I)-RSIG0(I)-PRES00*A)/(2.D0*XG0) 279 : DEPT(I)=DEPT(I)+TREPS0/3.D0*A 280 : DEFT(I)=iter2*DEPM(I)/nmax0-DEPT(I) 281 : 20 CONTINUE 282 : TREPP0=DEFT(1)+DEFT(2)+DEFT(3) 283 : * 284 : *---> Mise à jour de la limite en succion 285 : * 286 : SLIM1=(SLIM0+Pat0)*EXP(VOLU0/(XKS0-XLAM0)*TREPP0)-Pat0 287 : * 288 : *---> Mise à jour de la contrainte équivalente 289 : * 290 : STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0 291 : IF (STEST2.LE.1.D-10) STEST2=0.D0 292 : Stest=(STEST2)**(0.5D0) 293 : * 294 : *---> Mise à jour de la fonction de charge 295 : * 296 : PHIT=Stest*Stest+XM0*XM0*(PRES1+PS0)*(PRES1-PCS1) 297 : PHI0=ABS(PHIT) 298 : * 299 : *---> Test de convergence ou itération suivante 300 : * 301 : iter0=iter0+1 302 : IF (iter0.LT.200) THEN 303 : GOTO 99 304 : ELSE 305 : PHI0=0.D0 306 : KERRE=460 307 : GOTO 99 308 : ENDIF 309 : * 310 : ENDIF 311 : * 312 : * 313 : * Vérification des sous incrémentations 314 : * 315 : IF ((iter2.LT.nmax0).AND.(KERRE.EQ.0)) THEN 316 : GOTO 95 317 : ENDIF 318 : * 319 : * 320 : * 321 : * 322 : RETURN 323 : END 324 : 325 : 326 : 327 : 328 :
© Cast3M 2003 - All rights reserved.
Disclaimer