1 : C PSURY SOURCE CHAT 05/01/13 02:38:25 5004 2 : SUBROUTINE PSURY(ENDO,NENDO,NVARI,NSTRS,MFR1,DEPST,XMAT,VAR0, 3 : 1 RAPP,NRAPP, 4 : 1 SIG0,SIGF,VARF,NMATT,DEFP,KERRE) 5 : * 6 : * Entrées 7 : * 8 : * ENDO: courbe de debut d'endommagement 9 : * NENDO: nombre de points sur la courbe ENDO 10 : * RAPP: courbe d'evolution de l'endommagement en fonction de la 11 : * pseudo porosite 12 : * NRAPP: nombre de points de la courbe RAPP 13 : * NSTRS: nombre de composantes des deformations 14 : * MFR1: numero de la formulation 15 : * DEPST: increment de deformation totale 16 : * XMAT: donnees materiau 17 : * SIGF: contraintes plastiquement admissibles issues de ECOINC 18 : * VAR0: variables internes au debut du pas de temps 19 : * VARF: variables internes issues de ECOINC 20 : * SIG0: contraintes au debut du pas de temps 21 : * NMATT: nombre de composantes matériaux 22 : * DEFP: incrément de la déformation plastique 23 : * 24 : * Sorties 25 : * 26 : * SIGF: contraintes finales ( tiennent compte de l'endommagement) 27 : * VARF: variables internes finales 28 : * 29 : * 30 : IMPLICIT INTEGER(I-N) 31 : IMPLICIT REAL*8(A-H,O-Z) 32 : * 33 : -INC CCOPTIO 34 : DIMENSION SIGF(*),XMAT(*),VARF(*),DEPST(*),ENDO(*) 35 : DIMENSION SIG0(*),VAR0(*),DEFP(*),RAPP(*) 36 : DIMENSION RSIGF(6),RSIG0(6) 37 : * 38 : * 39 : *==================================================================== 40 : * Adaptation de l'option de calcul vers le 3D massif de SIGF a RSIGF 41 : *==================================================================== 42 : * 43 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 44 : * 45 : *---> 1 formulation massive 46 : *---> 2 formulation quasi incompressible 47 : *---> MASSIF 3D 48 : * 49 : IF (NSTRS .EQ. 6) THEN 50 : DO 110 I=1,NSTRS 51 : RSIGF(I)=SIGF(I) 52 : RSIG0(I)=SIG0(I) 53 : 110 CONTINUE 54 : ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0) 55 : & .OR.(IFOUR.EQ.-1).OR.(IFOUR.EQ.-2))) THEN 56 : * 57 : *---> Calcul en mode deformations planes ou axisymetrique 58 : *---> Calcul en mode contraintes planes 59 : * 60 : DO 115 I=1,NSTRS 61 : RSIGF(I)=SIGF(I) 62 : RSIG0(I)=SIG0(I) 63 : 115 CONTINUE 64 : RSIGF(5)=0.D0 65 : RSIG0(5)=0.D0 66 : RSIGF(6)=0.D0 67 : RSIG0(6)=0.D0 68 : ENDIF 69 : ELSE 70 : KERRE = 99 71 : RETURN 72 : ENDIF 73 : * 74 : *======================================== 75 : * Calcul de la densité du materiau 76 : *======================================== 77 : * 78 : * Si on a oublie de rentrer la densite initiale 79 : * 80 : * IF (IFOUR.EQ.-2) THEN 81 : * DENS0=XMAT(NMATT) 82 : * ELSE 83 : * DENS0=XMAT(3) 84 : * ENDIF 85 : * 86 : * DENS0=XMAT(3) 87 : DENS0=XMAT(NMATT) 88 : IF (DENS0.LE.1.D-10) THEN 89 : KERRE=71 90 : RETURN 91 : ENDIF 92 : * 93 : * Initialisation de la densite en variable interne 94 : * 95 : IF (VAR0(2).LE.1.D-10) THEN 96 : VAR0(2)=DENS0 97 : ENDIF 98 : IF (IFOUR.EQ.-2) THEN 99 : treps0=(1.D0-2.D0*XMAT(2))/(1.D0-XMAT(2)) 100 : treps0=treps0*(DEPST(1)+DEPST(2)-DEFP(1)-DEFP(2)) 101 : ELSE 102 : treps0=DEPST(1)+DEPST(2)+DEPST(3) 103 : ENDIF 104 : rho=DENS0/((DENS0/VAR0(2))+treps0) 105 : * 106 : * 107 : *======================================== 108 : * A t'on déja endommagé ? 109 : *======================================== 110 : * 111 : IF (ABS(VAR0(3)).GT.1.D-10) THEN 112 : * 113 : * On a déja endommagé 114 : *---------------------------------------- 115 : * 116 : * Calcul de la pseudo porosité du temps n+1 117 : * 118 : alpha=(VAR0(3)-rho)/VAR0(3) 119 : * 120 : * Calcul de la variable d'endommagement 121 : * 122 : D_end1=0.D0 123 : DO 12 I=1,NRAPP-1 124 : DD1=RAPP(2*(I-1)+1) 125 : DD2=RAPP(2*I+1) 126 : AL1=RAPP(2*I) 127 : AL2=RAPP(2*(I+1)) 128 : IF ((alpha.GE.AL1).AND.(alpha.LT.AL2)) THEN 129 : D_end1=(DD2-DD1)/(AL2-AL1)*(alpha-AL1)+DD1 130 : ENDIF 131 : 12 CONTINUE 132 : IF ((D_end1.LE.0.D0).OR.(alpha.LE.0.D0)) D_end1=0.D0 133 : IF ((D_end1.GE.1.D0).OR.(alpha.GE.1.D0)) D_end1=1.D0 134 : D_max=VAR0(5) 135 : * 136 : * Calcul de l'ancienne pseudo porosite ( temps n) 137 : * 138 : alp_old=VAR0(4) 139 : * 140 : * Calcul de la fonction d'endommagement g0 141 : * 142 : IF ((alpha.GT.0.D0).AND.(alpha.GT.alp_old)) THEN 143 : g0=1.D0-D_end1 144 : ELSE IF ((alpha.GT.0.D0).AND.(alpha.LE.alp_old)) THEN 145 : g0=1.D0-D_max 146 : ELSE 147 : g0=1.D0 148 : ENDIF 149 : g0=MAX(g0,0.D0) 150 : * 151 : * Calcul des contraintes vérifiant l'endommagement 152 : * 153 : DO 10 I=1,6 154 : RSIGF(I)=RSIGF(I)*g0 155 : 10 CONTINUE 156 : * 157 : * Mise a jour des variables internes 158 : * 159 : *---> Densite 160 : VARF(2)=rho 161 : *---> Densite de debut de fracture 162 : VARF(3)=VAR0(3) 163 : *---> Pseudo porosite 164 : VARF(4)=alpha 165 : *---> Pseudo porosite maximale 166 : VARF(5)=MAX(VAR0(5),D_end1) 167 : *---> Coefficient d'endommagement 168 : VARF(6)=D_end1 169 : *---> Fonction d'endommagement 170 : VARF(7)=g0 171 : * 172 : ELSE 173 : * 174 : * On a pas encore endommagé 175 : *------------------------------------------- 176 : * 177 : * Calcul de P=1/3*trace(RSIGF) 178 : * 179 : P0=(RSIGF(1)+RSIGF(2)+RSIGF(3))/3.D0 180 : * 181 : * Calcul de la contrainte équivalente 182 : * 183 : Y1=(RSIGF(1)*RSIGF(1))+(RSIGF(2)*RSIGF(2)) 184 : Y1=Y1+(RSIGF(3)*RSIGF(3)) 185 : Y1=Y1-(RSIGF(1)*RSIGF(2))-(RSIGF(2)*RSIGF(3)) 186 : Y1=Y1-(RSIGF(3)*RSIGF(1)) 187 : Y2=(RSIGF(4)*RSIGF(4))+(RSIGF(5)*RSIGF(5)) 188 : Y2=Y2+(RSIGF(6)*RSIGF(6)) 189 : Y0=Y1+(3.D0*Y2) 190 : Y0=(Y0)**(0.5D0) 191 : * 192 : * Rapport de triaxialite 193 : * 194 : Rapp0=P0/Y0 195 : * 196 : * Rapport de triaxialite de debut d'endommagement Rapp_th 197 : * 198 : EPSE=VARF(1) 199 : IF (EPSE.LT.ENDO(2)) THEN 200 : Rapp_th=1.D30 201 : pente1=1.D30 202 : Rmax=Rapp_th 203 : epsmin=ENDO(2) 204 : ELSE IF (EPSE.GE.ENDO(2*NENDO)) THEN 205 : Rapp_th=ENDO(2*NENDO-1) 206 : pente1=0.D0 207 : Rmax=Rapp_th 208 : epsmin=ENDO(2*NENDO) 209 : ELSE 210 : DO 100 I=1,(NENDO-1) 211 : IF ((EPSE.GE.ENDO(2*I)).AND.(EPSE.LT.ENDO(2*I+2))) THEN 212 : IND0=I 213 : ENDIF 214 : 100 CONTINUE 215 : epsmin=ENDO(2*IND0) 216 : epsmax=ENDO(2*IND0+2) 217 : dif0=ABS(epsmax-epsmin) 218 : Rmax=ENDO(2*IND0-1) 219 : Rmin=ENDO(2*IND0+1) 220 : IF (dif0.GT.1.D-10) THEN 221 : pente1=(Rmax-Rmin)/(epsmin-epsmax) 222 : Rapp_th=pente1*(EPSE-epsmin)+Rmax 223 : ELSE 224 : pente1=1.D30 225 : Rapp_th=Rmin 226 : ENDIF 227 : ENDIF 228 : * 229 : * Comparaison au rapport de triaxialite d'endommagement 230 : * 231 : IF (Rapp0.GT.Rapp_th) THEN 232 : * 233 : * On endommage 234 : * 235 : * Calcul de la densite de debut d'endommagement : rho_f 236 : *--------------------------------------------------- 237 : * 238 : *---> Rapport de triaxialite au debut du pas 239 : * 240 : P_old=(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0 241 : Y_old1=(RSIG0(1)*RSIG0(1))+(RSIG0(2)*RSIG0(2)) 242 : Y_old1=Y_old1+(RSIG0(3)*RSIG0(3)) 243 : Y_old1=Y_old1-(RSIG0(1)*RSIG0(2))-(RSIG0(2)*RSIG0(3)) 244 : Y_old1=Y_old1-(RSIG0(3)*RSIG0(1)) 245 : Y_old2=(RSIG0(4)*RSIG0(4))+(RSIG0(5)*RSIG0(5)) 246 : Y_old2=Y_old2+(RSIG0(6)*RSIG0(6)) 247 : Y_old0=Y_old1+(3.D0*Y2) 248 : Y_old=(Y_old0)**(0.5D0) 249 : R_old=P_old/Y_old 250 : * 251 : *---> Courbe d'evolution du rapport P/Y entre le debut et la fin du pas 252 : * On l'approxime par une droite de pente pente0 253 : * 254 : dvar1=VARF(1)-VAR0(1) 255 : IF (dvar1.GT.1.D-10) THEN 256 : pente0=(Rapp0-R_old)/dvar1 257 : ELSE 258 : pente0=1.D30 259 : ENDIF 260 : * 261 : *---> intersection avec la courbe de debut d'endommagement 262 : * point (R_int,eps_int) 263 : * 264 : 150 IF (pente1.GE.1.D30) THEN 265 : eps_int=epsmin 266 : ELSE IF (pente0.GE.1.D30) THEN 267 : eps_int=VAR0(1) 268 : ELSE 269 : eps_int=Rapp0-Rmax+(pente1*epsmin)-(pente0*VARF(1)) 270 : eps_int=eps_int/(pente1-pente0) 271 : ENDIF 272 : IF ((eps_int.LT.epsmin).AND.(IND0.GE.2)) THEN 273 : IND0=IND0-1 274 : epsmin=ENDO(2*IND0) 275 : epsmax=ENDO(2*IND0+2) 276 : dif0=ABS(epsmin-epsmax) 277 : Rmax=ENDO(2*IND0-1) 278 : Rmin=ENDO(2*IND0+1) 279 : IF (dif0.GT.1.D-10) THEN 280 : pente1=(Rmax-Rmin)/(epsmin-epsmax) 281 : ELSE 282 : pente1=1.D30 283 : ENDIF 284 : GOTO 150 285 : ELSE IF ((eps_int.LT.epsmin).AND.(IND0.EQ.1)) THEN 286 : pente1=1.D30 287 : eps_int=ENDO(2) 288 : epsmin=ENDO(2) 289 : epsmax=ENDO(2) 290 : Rmin=ENDO(1) 291 : Rmax=1.D30 292 : ENDIF 293 : IF (pente1.GE.1.D30) THEN 294 : R_int=pente0*(eps_int-VARF(1))+Rapp0 295 : ELSE 296 : R_int=pente1*(eps_int-epsmin)+Rmax 297 : ENDIF 298 : * 299 : *---> Calcul de rho_f en supposant le module d'ecrouissage 300 : * constant entre le debut et la fin du pas 301 : * 302 : IF (dvar1.GT.1.D-10) THEN 303 : H0=(Y0-Y_old)/dvar1 304 : Y_int=H0*(eps_int-VAR0(1))+Y_old 305 : ELSE 306 : * 307 : *---> Endommagement dans un cas élastique 308 : * 309 : IF ((DABS(Rapp0-R_old)).GT.1.D-20) THEN 310 : alfa0=(R_int-R_old)/(Rapp0-R_old) 311 : Y_int=(alfa0*Y0)+((1.D0-alfa0)*Y_old) 312 : ELSE 313 : Y_int=Y_old 314 : ENDIF 315 : ENDIF 316 : IF (treps0.GT.1.D-10) THEN 317 : XK0=XMAT(1)/(3.D0*(1.D0-2.D0*XMAT(2))) 318 : P_int=R_int*Y_int 319 : tr_eps1=(P_int-P_old)/XK0 320 : ELSE 321 : tr_eps1=0.D0 322 : ENDIF 323 : rho_f=DENS0/((DENS0/VAR0(2))+tr_eps1) 324 : IF ((rho_f.GT.1.D10).OR.(rho_f.LT.0.D0)) THEN 325 : rho_f=rho 326 : ENDIF 327 : * 328 : * Calcul de la pseudo-porosite 329 : * 330 : alpha=(rho_f-rho)/rho_f 331 : * 332 : * Calcul de la variable d'endommagement 333 : * 334 : D_end1=0.D0 335 : DO 13 I=1,NRAPP-1 336 : DD1=RAPP(2*(I-1)+1) 337 : DD2=RAPP(2*I+1) 338 : AL1=RAPP(2*I) 339 : AL2=RAPP(2*(I+1)) 340 : IF ((alpha.GE.AL1).AND.(alpha.LT.AL2)) THEN 341 : D_end1=(DD2-DD1)/(AL2-AL1)*(alpha-AL1)+DD1 342 : ENDIF 343 : 13 CONTINUE 344 : IF ((D_end1.LE.0.D0).OR.(alpha.LE.0.D0)) D_end1=0.D0 345 : IF ((D_end1.GE.1.D0).OR.(alpha.GE.1.D0)) D_end1=1.D0 346 : D_max=VAR0(5) 347 : * 348 : * Calcul de l'ancienne pseudo porosite ( temps n) 349 : * 350 : alp_old=VAR0(4) 351 : * 352 : * Calcul de la fonction d'endommagement g0 353 : * 354 : IF ((alpha.GT.0.D0).AND.(alpha.GT.alp_old)) THEN 355 : g0=1.D0-D_end1 356 : ELSE IF ((alpha.GT.0.D0).AND.(alpha.LE.alp_old)) THEN 357 : g0=1.D0-D_max 358 : ELSE 359 : g0=1.D0 360 : ENDIF 361 : g0=MAX(g0,0.D0) 362 : * 363 : * Calcul des contraintes vérifiant l'endommagement 364 : * 365 : DO 15 I=1,6 366 : RSIGF(I)=RSIGF(I)*g0 367 : 15 CONTINUE 368 : * 369 : * Mise a jour des variables internes 370 : * 371 : VARF(2)=rho 372 : VARF(3)=rho_f 373 : VARF(4)=alpha 374 : VARF(5)=MAX(VAR0(5),D_end1) 375 : VARF(6)=D_end1 376 : VARF(7)=g0 377 : IF (D_end1.GT.0.9) THEN 378 : write(*,*) ' Debut endommagement' 379 : write(*,*) 'P_int,P_old=',P_int,P_old 380 : write(*,*) 'R_int,Y_int=',R_int,Y_int 381 : write(*,*) 'r0,rho,tr_eps1=',DENS0,VAR0(2),tr_eps1 382 : write(*,*) 'D_end1,gg0=',D_end1,g0 383 : write(*,*) 'alpha,alp_old=',alpha,alp_old 384 : write(*,*) 'DD1,DD2,AL1,AL2=',DD1,DD2,AL1,AL2 385 : write(*,*) 'rho,rho_f=',rho,rho_f 386 : write(*,*) 'DENS0,VAR0(2)=',DENS0,VAR0(2) 387 : ENDIF 388 : * 389 : ELSE 390 : * 391 : * On n'endommage pas 392 : * 393 : VARF(2)=rho 394 : VARF(3)=0.D0 395 : VARF(4)=0.D0 396 : VARF(5)=0.D0 397 : VARF(6)=0.D0 398 : VARF(7)=1.D0 399 : ENDIF 400 : * 401 : *======================================= 402 : * Fin du calcul d'endommagement 403 : *======================================= 404 : ENDIF 405 : * 406 : * 407 : *========================================================= 408 : * Passage a l'option de calcul pour les contraintes 409 : *========================================================= 410 : * 411 : IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN 412 : IF (NSTRS .EQ. 6) THEN 413 : * 414 : *---> MASSIF 3D 415 : * 416 : DO 170 I=1,NSTRS 417 : SIGF(I)=RSIGF(I) 418 : 170 CONTINUE 419 : ELSE IF ( NSTRS .EQ. 4 ) THEN 420 : * 421 : *---> Calcul axisymétrique ou contraintes planes 422 : * 423 : DO 180 I=1,NSTRS 424 : SIGF(I)=RSIGF(I) 425 : 180 CONTINUE 426 : ENDIF 427 : ENDIF 428 : RETURN 429 : * 430 : END 431 : 432 : 433 : 434 : 435 : 436 : 437 :
© Cast3M 2003 - All rights reserved.
Disclaimer