1 : C MODVMS SOURCE PV 08/09/11 21:16:15 6150 2 : C MODVONMISES 3 : SUBROUTINE MODVONMISES (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI, 4 : 1 NMATT,ISTEP,ICARA,JDIM,IFOUR2) 5 : C 6 : * BCN 7 : C New source: Modified Von Mises nonlocal damage model added (inspired in MAZARS) 8 : C Goal: satisfactory simulation of the single-edge notched beam 9 : C 10 : C Modifications: 11 : C 1. The definition of the equivalent strain (epsilontild) 12 : C is changed [1] 13 : C 2. The distinction between damage in traction (DT) and damage 14 : C in compression (DC) is suppressed. Only one damage (D) is used 15 : C 3. The relationship between the nonlocal equivalent strain and the 16 : C damage is changed : Two choices for damage vs. nonlocal equivalent 17 : C strain. [2] 18 : C a) LOI = 1. Exponential law. 19 : C b) LOI = 0. Polynomial law. 20 : C 21 : C References: 22 : C [1] Peerlings, R.H.J., de Borst, R., Brekelmans, W.A.M., and Geers, M.D. 23 : C (1998), Gradient-enhanced damage modelling of concrete fracture, 24 : C Mechanics of Cohesive-Frictional Materials, 3, 323-342. 25 : C [2] Rodriguez-Ferran A., Huerta A. (2000), Error estimation and 26 : C adaptivity for nonlocal damage models. International Journal of 27 : C Solids and Structures, 37, 7501-7528. 28 : * BCN 29 : C 30 : C variables en entree 31 : C 32 : C WRK0 pointeur sur un segment deformation au pas precedent 33 : C 34 : C WRK1 pointeur sur un segment increment de deformation 35 : C 36 : C WRKK2 pointeur sur un segment variables internes au pas precedent 37 : C 38 : C WRK5 pointeur sur un segment de deformations inelastiques 39 : C 40 : C XMATER constantes du materiau 41 : C 42 : C NSTRS nombre de composantes dans les vecteurs des contraintes 43 : C et les vecteurs des deformations 44 : C 45 : C NVARI nombre de variables internes (doit etre egal a 2) 46 : C 47 : C NMATT nombre de constantes du materiau 48 : C 49 : C ISTEP flag utilise pour separer les etapes dans un calcul non local 50 : C ISTEP=0 -----> calcul local 51 : C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils 52 : C ISTEP=2 -----> calcul non local etape 2 on continue le calcul 53 : C a partir des seuils moyennes 54 : C 55 : C Modif L.Bode - 14/10/92 56 : C Nouveaux parametres en entree 57 : C JDIM Dimension de travail 58 : C ( Coques JDIM =2 , Massifs JDIM = IDIM ) 59 : C IFOUR2 Type de formulation 60 : C ( Coques IFOUR2 = -2 => contraintes planes , 61 : C Massifs IFOUR2 = IFOUR) 62 : C 63 : C variables en sortie 64 : C 65 : C VARINF variables internes finales 66 : C 67 : C SIGMAF contraintes finales 68 : C 69 : C C. LA BORDERIE MARS 1992 70 : C declaration des variables 71 : C 72 : C 73 : IMPLICIT INTEGER(I-N) 74 : IMPLICIT REAL*8(A-H,O-Z) 75 : -INC CCOPTIO 76 : SEGMENT WRK0 77 : REAL*8 XMAT(NMATT) 78 : ENDSEGMENT 79 : * 80 : SEGMENT WRK1 81 : REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) 82 : REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) 83 : REAL*8 DEFP(NSTRS),XCAR(ICARA) 84 : ENDSEGMENT 85 : SEGMENT WRKK2 86 : REAL*8 EPSILI(NSTRS),DSIGT(NSTRS) 87 : ENDSEGMENT 88 : C Modif L.Bode - 14/10/92 89 : SEGMENT WRK3 90 : REAL*8 EPSILO(NSTRS) 91 : ENDSEGMENT 92 : C Definition dynamique de EPSILO 93 : C Fin modif L.Bode 94 : SEGMENT WRK5 95 : REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS) 96 : ENDSEGMENT 97 : CHARACTER*8 CMATE 98 : INTEGER NSTRS,NVARI,NMATT,ISTEP 99 : REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3) 100 : * BCN 101 : C REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC 102 : C REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA 103 : REAL*8 SIGP(3) 104 : REAL*8 YOUN,XNU,EPSD0,B1,B2,RATIO_CT,LOI 105 : * BCN 106 : INTEGER ISTRS,JSTRS,KCAS,IRTD 107 : REAL*8 XZERO,UN,DEUX,XPETIT 108 : * BCN 109 : C REAL*8 DINI,D,DT,DC,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA 110 : REAL*8 DINI,D,EPSTIL,EPSTIM 111 : REAL*8 TINV1,TINV2,TINV3,AUX1,AUX2,AUX3,AUX4,AUX5 112 : * BCN 113 : PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12) 114 : 115 : C 116 : C 117 : C recuperation des variables initiales dans les tableaux 118 : C 119 : C 120 : N=NSTRS 121 : CMATE = 'ISOTROPE' 122 : YOUN = XMAT(1) 123 : XNU = XMAT(2) 124 : EPSD0= XMAT(5) 125 : * BCN 126 : C The material parameters change due to the modification of the model 127 : C EPSD0: damage threshold [2] 128 : C B1 and B2: material parameters in equivalent strain-damage law [2] 129 : C RATIO_CT: ratio of compressive to tensile strength [1,2] 130 : C LOI: Choice of equivalent strain-damage law [2] 131 : C ACOM = XMAT(6) 132 : C BCOM = XMAT(7) 133 : C ATRA = XMAT(8) 134 : C BTRA = XMAT(9) 135 : C BETA = XMAT(10) 136 : B1 = XMAT(6) 137 : B2 = XMAT(7) 138 : RATIO_CT = XMAT(8) 139 : LOI = XMAT(9) 140 : * BCN 141 : DINI = VAR0(2) 142 : C 143 : C calcul des seuils 144 : C 145 : C 146 : C calcul de la deformation totale 147 : C 148 : SEGINI WRK3 149 : DO 100 ISTRS=1,NSTRS 150 : EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS) 151 : 100 CONTINUE 152 : C 153 : C calcul des deformations principales 154 : C 155 : C 156 : C on reecrit les deformations sous forme matricielle 157 : C 158 : C Modif L.Bode - 14/10/92 159 : C Rajout de IFOUR2 en argument de ENDOCA 160 : * print*,'on appelle ENDOCB' 160 : CALL ENDOCB (EPSILO,EPS33,2,IFOUR2) 161 : * print*,'apres endocb' 162 : C Fin modif L.Bode 163 : C 164 : C et on diagonalise 165 : C 166 : C Modif L.Bode - 14/10/92 167 : C Pour les elts Coques, on travaille en contraintes planes => JDIM =2 168 : C Pour les elts Massifs JDIM =IDIM 169 : * print*,'avant JACOB3' 169 : CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33) 170 : * print*,'apres JACOB3' 171 : C Fin modif L.Bode 172 : IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN 173 : C 174 : C on calcule la matrice de hooke et les contraintes ppales 175 : C 176 : CMATE = 'ISOTROPE' 177 : KCAS=1 178 : C Modif L.Bode - 14/10/92 179 : C IFOUR --> IFOUR2 dans appel DOHMAS 180 : * print*,'avant dohmas' 180 : CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS,KCAS,DDHOOK,IRTD) 181 : * print*,'apres dohmas' 182 : C Fin modif L.Bode 183 : DO 200 ISTRS=1,3 184 : SIGP(ISTRS)= XZERO 185 : DO 210 JSTRS=1,3 186 : SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS) 187 : 210 CONTINUE 188 : 200 CONTINUE 189 : END IF 190 : C 191 : C on complete la deformation dans le cas des contraintes planes 192 : C 193 : C Modif L.Bode - 14/10/92 194 : C IFOUR remplace par IFOUR2 195 : IF (IFOUR2.EQ. -2) THEN 196 : EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU) 197 : END IF 198 : C Fin modif L.Bode 199 : C 200 : C on calcule le epsilontild 201 : C 202 : * BCN 203 : C The original definition of epsilontild is replaced. 204 : C TINV1: first strain tensor invariant 205 : C TINV2: second (deviatoric) strain tensor invariant 206 : C 207 : C IMPORTANT! There is a discrepancy by a factor of 6 in the 208 : C definitions of the second invariant in [1] and [2] 209 : C The definition of Askes and Sluys is used here. 210 : C 211 : C EPSTIL=MAX( XZERO , EPSIPP(1) )**2 + 212 : C 1 MAX( XZERO , EPSIPP(2) )**2 + 213 : C 2 MAX( XZERO , EPSIPP(3) )**2 214 : C EPSTIL=SQRT (EPSTIL) 215 : C 216 : TINV1 = EPSIPP(1)+EPSIPP(2)+EPSIPP(3) 217 : TINV2 = (EPSIPP(1)- EPSIPP(2))*(EPSIPP(1)- EPSIPP(2)) + 218 : . (EPSIPP(1)- EPSIPP(3))*(EPSIPP(1)- EPSIPP(3)) + 219 : . (EPSIPP(2)- EPSIPP(3))*(EPSIPP(2)- EPSIPP(3)) 220 : TINV2 = TINV2/6.d0 221 : C 222 : AUX1 = RATIO_CT-UN 223 : AUX2 = UN+XNU 224 : AUX3 = UN-DEUX*XNU 225 : AUX4 = AUX1*TINV1/AUX3 226 : AUX5 = DSQRT((AUX4*AUX4)+((12.d0*RATIO_CT*TINV2)/(AUX2*AUX2))) 227 : C 228 : EPSTIL = (AUX1*TINV1)/(DEUX*RATIO_CT*AUX3) + 229 : . AUX5/(DEUX*RATIO_CT) 230 : * BCN 231 : IF (ISTEP .EQ. 0) THEN 232 : EPSTIM=EPSTIL 233 : VARF(1)=EPSTIL 234 : ELSE IF (ISTEP .EQ. 1) THEN 235 : VARF(2)=DINI 236 : VARF(1)=EPSTIL 237 : ELSE IF (ISTEP .EQ. 2) THEN 238 : EPSTIM=VAR0(1) 239 : VARF(1)=EPSTIM 240 : ELSE 241 : PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP 242 : END IF 243 : IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN 244 : C 245 : C on calcule l'endommagement et les contraintes 246 : C 247 : C 248 : C on calcule ALFAT ALFAC DT et DC puis D 249 : C dans le cas ou le seuil initial est depasse 250 : C 251 : IF ( EPSTIM .GT. EPSD0) THEN 252 : C 253 : C calcul de l'endommagement 254 : * BCN 255 : C The equivalent strain-damage law is modified 256 : C 257 : C on calcule le signe des contraintes elastiques 258 : C 259 : C DO 300 ISTRS=1,3 260 : C IF (SIGP(ISTRS).LT. XZERO) THEN 261 : C SIGPC(ISTRS)=SIGP(ISTRS) 262 : C SIGPT(ISTRS)=XZERO 263 : C ELSE 264 : C SIGPT(ISTRS)=SIGP(ISTRS) 265 : C SIGPC(ISTRS)=XZERO 266 : C END IF 267 : C300 CONTINUE 268 : C TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3) 269 : C TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3) 270 : C 271 : C on calcule les deformations dues aux contraintes positives 272 : C 273 : C DO 400 ISTRS=1,3 274 : C EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN 275 : C400 CONTINUE 276 : C 277 : C on en deduit ALFAT et ALFAC 278 : C 279 : C ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) + 280 : C 1 MAX(XZERO,EPSIPP(2))*EPSILT(2) + 281 : C 2 MAX(XZERO,EPSIPP(3))*EPSILT(3) 282 : C ALFAT = ALFAT/(EPSTIL*EPSTIL) 283 : C ALFAC = UN - ALFAT 284 : C 285 : C modification pour la bi ou tricompression 286 : C 287 : C IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN 288 : C GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+ 289 : C 1 SIGPC(3)*SIGPC(3) 290 : C GAMMA=-SQRT(GAMMA)/TRSIGC 291 : C EPSTIM=EPSTIM*GAMMA 292 : C END IF 293 : C 294 : C amelioration de la reponse en cisaillement pour beta > 1. 295 : C 296 : C IF (BETA .GT. UN) THEN 297 : C IF ( ALFAT .GT. XPETIT ) THEN 298 : C ALFAT=ALFAT**BETA 299 : C END IF 300 : C IF ( ALFAC .GT. XPETIT ) THEN 301 : C ALFAC=ALFAC**BETA 302 : C END IF 303 : C END IF 304 : C 305 : C 306 : C on calcule DT et DC puis D 307 : C dans le cas ou le seuil initial est depasse 308 : C 309 : C on est oblige de verifier car on a pu multiplier par gamma 310 : C 311 : C IF (EPSTIM .GT. EPSD0) THEN 312 : C DT=UN - EPSD0*(UN-ATRA)/EPSTIM - 313 : C 1 ATRA*EXP(-BTRA*(EPSTIM-EPSD0)) 314 : C DC=UN - EPSD0*(UN-ACOM)/EPSTIM - 315 : C 1 ACOM*EXP(-BCOM*(EPSTIM-EPSD0)) 316 : C ELSE 317 : C DT=XZERO 318 : C DC=XZERO 319 : C END IF 320 : C D = ALFAT*DT + ALFAC*DC 321 : C 322 : ***tc Y0 n'etait pas initialisé on suppose que c'est EPSD0 323 : Y0 = epsd0 324 : C Two choices for damage vs. internal variable: 325 : IF (LOI.EQ.XZERO) THEN 326 : D = UN - (UN/(UN+B1*(EPSTIM-Y0)+ 327 : . B2*(EPSTIM-Y0)*(EPSTIM-Y0))) 328 : ELSE IF (LOI.EQ.UN) THEN 329 : D = UN - (Y0*(UN-B2)/EPSTIM)- B2*EXP(B1*(Y0-EPSTIM)) 330 : ELSE 331 : WRITE(*,*) 'LOI SHOULD BE EQUAL TO 1 OR 0.(SEE MATE)' 332 : STOP 333 : ENDIF 334 : * BCN 335 : C 336 : C on borne la valeur de D a 0.99999 337 : C 338 : D=MIN ( D , UN-1.D-05 ) 339 : ELSE 340 : D=XZERO 341 : END IF 342 : C 343 : C on teste la croissance de D 344 : C 345 : D=MAX ( D , DINI ) 346 : C 347 : C on le stocke dans les variables internes finales 348 : C 349 : VARF(2)= D 350 : C 351 : C on calcule les contraintes finales 352 : C 352 : CALL MATVE1 (DDHOOK,EPSILO,NSTRS,NSTRS,SIGF,2) 353 : DO 500 ISTRS=1,NSTRS 354 : SIGF(ISTRS)=SIGF(ISTRS)*(UN-D) 355 : 500 CONTINUE 356 : C 357 : C et les deformations inelastiques finales 358 : C 359 : DO 600 ISTRS=1,NSTRS 360 : EPINF(ISTRS)=EPSILO(ISTRS)*D 361 : 600 CONTINUE 362 : 363 : END IF 364 : SEGSUP WRK3 365 : RETURN 366 : END 367 : 368 : 369 : 370 : 371 : 372 :
© Cast3M 2003 - All rights reserved.
Disclaimer