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