1 : C DESMNL SOURCE PV 11/03/07 21:16:17 6885 2 : C DESMORAT 3 : SUBROUTINE DESMNL(WRK0,WRK1,WRK5,NSTRS,NVARI,LHOOK,ISTEP, 4 : & NMATT,ICARA,KERRE,MFR) 5 : C 6 : C variables en entree 7 : C 8 : C WRK0, WRK1 pointeur actif a l'entree 9 : C NVARI nombre de variables internes (doit etre egal a ?) 10 : C NMATT nombre de constantes du materiau 11 : C 12 : C ISTEP flag utilise pour separer les etapes dans un calcul non local 13 : C ISTEP=0 -----> calcul local 14 : C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils 15 : C ISTEP=2 -----> calcul non local etape 2 on continue le calcul 16 : C a partir des seuils moyennes 17 : C 18 : C variables en sortie 19 : C 20 : C VARINF variables internes finales 21 : C 22 : C SIGMAF contraintes finales 23 : C 24 : C 25 : IMPLICIT REAL*8(A-H,O-Z) 26 : C 27 : -INC CCOPTIO 28 : -INC SMINTE 29 : -INC SMCOORD 30 : -INC SMELEME 31 : C 32 : C 33 : CHARACTER*8 CMATE 34 : REAL*8 XD0(3,3),DF(3,3),POSEPSI2(3,3),EPSIN(6),EPS33(3,3) 35 : REAL*8 SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6) 36 : REAL*8 UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3) 37 : REAL*8 ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3) 38 : REAL*8 EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3) 39 : REAL*8 VPDF(2),VECPDF(2,2),EPSTF(6),DHOOK(6,6) 40 : C 41 : PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12) 42 : C 43 : SEGMENT WRK0 44 : REAL*8 XMAT(NMATT) 45 : ENDSEGMENT 46 : C 47 : SEGMENT WRK1 48 : REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) 49 : REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) 50 : REAL*8 DEFP(NSTRS),XCAR(ICARA) 51 : ENDSEGMENT 52 : C 53 : SEGMENT WRK5 54 : REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS) 55 : ENDSEGMENT 56 : C 57 : C 58 : C Recuperation des variables initiales dans les tableaux 59 : C 60 : CMATE= 'ISOTROPE' 61 : YOUN = XMAT(1) 62 : XNU = XMAT(2) 63 : XK0 = XMAT(5) 64 : AA = XMAT(6) 65 : A = XMAT(7) 66 : etaC = XMAT(8) 67 : etaT = XMAT(9) 68 : Xk = XMAT(10) 69 : DC=XMAT(11) 70 : C 71 : C Recuperation de l'endommagement initial 72 : C 73 : XD0(1,1)=VAR0(2) 74 : XD0(2,2)=VAR0(3) 75 : XD0(3,3)=VAR0(4) 76 : XD0(1,2)=VAR0(5) 77 : XD0(1,3)=VAR0(6) 78 : XD0(2,3)=VAR0(7) 79 : XD0(2,1)=XD0(1,2) 80 : XD0(3,1)=XD0(1,3) 81 : XD0(3,2)=XD0(2,3) 82 : C 83 : C Recuperation de la base d'endommagement 84 : C 85 : ROT(1,1)=VAR0(8) 86 : ROT(1,2)=VAR0(9) 87 : ROT(1,3)=VAR0(10) 88 : ROT(2,1)=VAR0(11) 89 : ROT(2,2)=VAR0(12) 90 : ROT(2,3)=VAR0(13) 91 : ROT(3,1)=VAR0(14) 92 : ROT(3,2)=VAR0(15) 93 : ROT(3,3)=VAR0(16) 94 : C 95 : C Recuperation de l'etat d'endommagement critique 96 : C 97 : IFAILURE=nint(VAR0(17)) 98 : C 99 : C Matrice de raideur elastique 99 : CALL ZERO(DDHOOK,6,6) 100 : AUX0=YOUN/((UN+XNU)*(UN-XNU-XNU)) 101 : AUX=AUX0*(UN-XNU) 102 : AUX1=AUX0*XNU 103 : GEGE=0.5D0*YOUN/(UN+XNU) 104 : DDHOOK(1,1)=AUX 105 : DDHOOK(2,1)=AUX1 106 : DDHOOK(3,1)=AUX1 107 : DDHOOK(1,2)=AUX1 108 : DDHOOK(2,2)=AUX 109 : DDHOOK(3,2)=AUX1 110 : DDHOOK(1,3)=AUX1 111 : DDHOOK(2,3)=AUX1 112 : DDHOOK(3,3)=AUX 113 : DDHOOK(4,4)=GEGE 114 : DDHOOK(5,5)=GEGE 115 : DDHOOK(6,6)=GEGE 116 : C 117 : C 118 : C Deformation totale 119 : C 120 : DO ISTRS=1,NSTRS 121 : EPSTF(ISTRS)=VAR0(17+ISTRS)+DEPST(ISTRS) 122 : ENDDO 123 : C - Trace de la deformation 124 : TRAEPS=EPSTF(1)+EPSTF(2)+EPSTF(3) 125 : C 126 : C calcul des deformations principales 126 : CALL ENDOCB(EPSTF,EPS33,2,IFOUR) 126 : CALL JACOB3(EPS33,3,EPSIPP,VECP33) 127 : EPSIPOS(1)= MAX( XZERO , EPSIPP(1) )**2 128 : EPSIPOS(2)= MAX( XZERO , EPSIPP(2) )**2 129 : EPSIPOS(3)= MAX( XZERO , EPSIPP(3) )**2 130 : EPSTILNLOC=SQRT(EPSIPOS(1)+EPSIPOS(2)+EPSIPOS(3)) 131 : C 132 : IF (ISTEP.EQ.0) THEN 133 : EPSTIL=EPSTILNLOC 134 : VARF(1)=EPSTIL 135 : ELSE IF (ISTEP.EQ.1) THEN 136 : VARF(1)=EPSTILNLOC 137 : DO I=2,23 138 : VARF(I)=VAR0(I) 139 : ENDDO 140 : GOTO 2000 141 : ELSE IF (ISTEP.EQ.2) THEN 142 : c Over nonlocal, coef 2 143 : C EPSTIL=DEUX*VAR0(1)-EPSTILNLOC 144 : EPSTIL=VAR0(1) 145 : ENDIF 146 : C 147 : C********************************************************** 148 : C Début de l'intégration numérique du comportement 149 : C********************************************************** 150 : C 151 : C calcul des seuils 152 : C 153 : C partie positive de la matrice des deformations 153 : CALL PRODT2(POSEPSI2,EPSIPOS,VECP33,3) 154 : C 155 : C calcul de la trace de D au pas n 156 : TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) 157 : C 158 : C calcul du seuil 159 : XKAPPA=A*DTAN(TRAD0/(A*AA)+DATAN(XK0/A)) 160 : FCRIT=EPSTIL+Xk*TRAEPS-XKAPPA 161 : C 162 : C TEST sur le seuil 163 : IF(FCRIT.LT.XZERO) THEN 164 : C on est elastique 165 : C endommagement et deformations permanentes inchanges 166 : DO I=1,3 167 : DO J=1,3 168 : DF(I,J)=XD0(I,J) 169 : ENDDO 170 : ENDDO 171 : DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0 172 : ELSE 173 : C 174 : C debut de la correction d'endom. 175 : C 176 : C Calcul de Kappa^(-1)=Tr Dn+1 177 : C 178 : TRADF=A*AA*(ATAN((EPSTIL+Xk*TRAEPS)/A)-ATAN(XK0/A)) 179 : C 180 : C test sur l'eventuelle presence d'une direction d'endommagement fixee 181 : C 182 : IF(IFAILURE.EQ.0)THEN 183 : C actualisation de l'endommagement, Dc pas encore atteint (cas "normal") 184 : IF (EPSTIL.EQ.XZERO) THEN 185 : DLAMBDA=XZERO 186 : ELSE 187 : DLAMBDA=(TRADF-TRAD0)/(EPSTIL**2) 188 : ENDIF 189 : IF(DLAMBDA.LT.XZERO) THEN 190 : DLAMBDA=XZERO 191 : ENDIF 191 : CALL ZERO(DF,3,3) 192 : DO I=1,3 193 : DO J=1,3 194 : DF(I,J)=XD0(I,J)+DLAMBDA*POSEPSI2(I,J) 195 : ENDDO 196 : ENDDO 197 : C On borne l'endommagement D(1) a Dc si D>Dc et Fixed-crack 197 : CALL JACOD3(DF,3,DFVP) 198 : IF(DFVP(1).GE.DC) THEN 199 : IFAILURE=1 200 : C mise a jour de XD0 pour traiter IFAILURE=1 dans le meme pas 201 : DO I=1,3 202 : DO J=1,3 203 : XD0(I,J)=DF(I,J) 204 : ENDDO 205 : ENDDO 206 : TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) 207 : TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3) 207 : CALL JACOB3(XD0,3,DFVP,ROT) 208 : ENDIF 209 : ENDIF 210 : C 211 : IF(IFAILURE.EQ.1)THEN 212 : C une direction d'endom est fixee a Dc 213 : C 214 : C On calcule le multiplicateur d'endommagement dans la bonne base 215 : C On met les deformations positives dans le repère principal 216 : C de l'endommagement 216 : CALL JACOD3(XD0,3,DFVP) 216 : CALL TRSPOD(ROT,3,3,ROTT) 216 : CALL MATMAT(POSEPSI2,ROT,3,3,3,TRAV) 216 : CALL MATMAT(ROTT,TRAV,3,3,3,POSEPSI2) 217 : C 218 : IF(ABS(POSEPSI2(2,2)+POSEPSI2(3,3)).GT.XPETIT) THEN 219 : DLAMBDA=(TRADF-TRAD0)/(POSEPSI2(2,2)+POSEPSI2(3,3)) 220 : ELSE 221 : DLAMBDA=XZERO 222 : ENDIF 223 : IF(DLAMBDA.LT.XZERO) THEN 224 : DLAMBDA=XZERO 225 : ENDIF 226 : C 227 : C On corrige l'endommagement 227 : CALL ZERO(DF,3,3) 228 : DF(1,1)=Dc 229 : DF(2,2)=DFVP(2)+DLAMBDA*POSEPSI2(2,2) 230 : DF(3,3)=DFVP(3)+DLAMBDA*POSEPSI2(3,3) 231 : DF(2,3)=DLAMBDA*POSEPSI2(2,3) 232 : DF(3,2)=DF(2,3) 233 : C 233 : CALL MATMAT(DF,ROTT,3,3,3,TRAV) 233 : CALL MATMAT(ROT,TRAV,3,3,3,DF) 234 : C 235 : C On verifie la deuxieme direction d'endommagement 235 : CALL JACOD3(DF,3,DFVP) 236 : IF(DFVP(2).GE.Dc)THEN 237 : IFAILURE=2 238 : C mise a jour de XD0 pour traiter IFAILURE=2 dans le meme pas 239 : DO I=1,3 240 : DO J=1,3 241 : XD0(I,J)=DF(I,J) 242 : ENDDO 243 : ENDDO 244 : TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) 245 : TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3) 245 : CALL JACOB3(XD0,3,DFVP,ROT) 246 : ENDIF 247 : ENDIF 248 : C 249 : IF(IFAILURE.EQ.2) THEN 250 : C deux directions d'endom sont fixee a Dc 250 : CALL TRSPOD(ROT,3,3,ROTT) 251 : C 252 : DF33=TRADF-DEUX*Dc 253 : C print*, 'DF33=',DF33 254 : IF(DF33.LT.XZERO) THEN 255 : DF33=XPETIT 256 : ENDIF 257 : C On corrige l'endommagement 257 : CALL ZERO(DF,3,3) 258 : DF(1,1)=Dc 259 : DF(2,2)=Dc 260 : DF(3,3)=DF33 261 : C 261 : CALL MATMAT(DF,ROTT,3,3,3,TRAV) 261 : CALL MATMAT(ROT,TRAV,3,3,3,DF) 262 : C 263 : C On verifie la troisieme direction d'endommagement 264 : IF(DF33.GE.Dc)THEN 265 : IFAILURE=3 266 : C mise a jour de XD0 pour traiter IFAILURE=3 dans le même pas 267 : DO I=1,3 268 : DO J=1,3 269 : XD0(I,J)=DF(I,J) 270 : ENDDO 271 : ENDDO 272 : TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) 272 : CALL JACOB3(XD0,3,DFVP,ROT) 273 : ENDIF 274 : ENDIF 275 : C 276 : IF(IFAILURE.EQ.3)THEN 277 : C On est totalement rompu, endommagement isotrope egal a DC 277 : CALL TRSPOD(ROT,3,3,ROTT) 278 : C On corrige l'endommagement 278 : CALL ZERO(DF,3,3) 279 : DF(1,1)=Dc 280 : DF(2,2)=Dc 281 : DF(3,3)=Dc 282 : C 282 : CALL MATMAT(DF,ROTT,3,3,3,TRAV) 282 : CALL MATMAT(ROT,TRAV,3,3,3,DF) 283 : ENDIF 284 : C 285 : C Fin de la correction d'endommagement 286 : C 287 : C 288 : DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0 289 : C 290 : ENDIF 291 : C 292 : C Fin des differents cas (elasticite ou endommagement) 293 : C 294 : C Calcul de la contrainte finale 295 : C 296 : DO I=1,3 297 : DO J=1,3 298 : UNMOINSD(I,J)=-UN*DF(I,J) 299 : ENDDO 300 : ENDDO 301 : UNMOINSD(1,1)=UN-DF(1,1) 302 : UNMOINSD(2,2)=UN-DF(2,2) 303 : UNMOINSD(3,3)=UN-DF(3,3) 303 : CALL JACOB3(UNMOINSD,3,UNMOINSDVP,VECP33) 304 : HM1DIAG(1)=DSQRT(UNMOINSDVP(1)) 305 : HM1DIAG(2)=DSQRT(UNMOINSDVP(2)) 306 : HM1DIAG(3)=DSQRT(UNMOINSDVP(3)) 306 : CALL PRODT2(HM1,HM1DIAG,VECP33,3) 307 : C 308 : C contrainte effective 308 : CALL MATVE1(DDHOOK,EPSTF,NSTRS,NSTRS,SIGTILDE,2) 308 : CALL ENDOCB(SIGTILDE,SIGTILDE33,1,IFOUR) 308 : CALL MATMAT(HM1,SIGTILDE33,3,3,3,TRAV) 308 : CALL MATMAT(TRAV,HM1,3,3,3,SIGD33) 309 : XTRAV=XZERO 310 : DO I=1,3 311 : DO J=1,3 312 : XTRAV=XTRAV+UNMOINSD(I,J)*SIGTILDE33(I,J) 313 : ENDDO 314 : ENDDO 315 : DO I=1,3 316 : DO J=1,3 317 : SIGD33(I,J)=SIGD33(I,J)-(XTRAV/(3.D0*(UN-DH)))*UNMOINSD(I,J) 318 : ENDDO 319 : ENDDO 320 : TRASIGTILDE=SIGTILDE(1)+SIGTILDE(2)+SIGTILDE(3) 321 : C 322 : IF(TRASIGTILDE.GE.XZERO) THEN 323 : COEF=UN-etaT*DH 324 : IF(COEF.LE.XZERO) THEN 325 : COEF=UN-Dc 326 : ENDIF 327 : TRASIG=COEF*TRASIGTILDE 328 : ELSE 329 : TRASIG=(UN-etaC*DH)*TRASIGTILDE 330 : ENDIF 331 : C 332 : SIGF(1)=SIGD33(1,1)+TRASIG/3.d0 333 : SIGF(2)=SIGD33(2,2)+TRASIG/3.d0 334 : SIGF(3)=SIGD33(3,3)+TRASIG/3.d0 335 : SIGF(4)=SIGD33(1,2) 336 : SIGF(5)=SIGD33(1,3) 337 : SIGF(6)=SIGD33(2,3) 338 : C 339 : C 340 : C stockage des variables du probleme 341 : C 342 : VARF(2)=DF(1,1) 343 : VARF(3)=DF(2,2) 344 : VARF(4)=DF(3,3) 345 : VARF(5)=DF(1,2) 346 : VARF(6)=DF(1,3) 347 : VARF(7)=DF(2,3) 348 : VARF(8)=ROT(1,1) 349 : VARF(9)=ROT(1,2) 350 : VARF(10)=ROT(1,3) 351 : VARF(11)=ROT(2,1) 352 : VARF(12)=ROT(2,2) 353 : VARF(13)=ROT(2,3) 354 : VARF(14)=ROT(3,1) 355 : VARF(15)=ROT(3,2) 356 : VARF(16)=ROT(3,3) 357 : C 358 : C flag pour l'endommagement maximum 359 : VARF(17)=IFAILURE 360 : C 361 : C Deformations totales finales 362 : DO ISTRS=1,NSTRS 363 : VARF(17+ISTRS)=EPSTF(ISTRS) 364 : ENDDO 365 : C 366 : C 367 : 2000 CONTINUE 368 : C 369 : RETURN 370 : END 371 : C 372 : 373 : 374 :
© Cast3M 2003 - All rights reserved.
Disclaimer