1 : C DESMOR SOURCE PV 11/03/07 21:16:19 6885 2 : C DESMORAT (V2) 3 : SUBROUTINE DESMOR(WRK52,WRK53,WRK54,NVARI,Iecou) 4 : C 5 : C variables en entree 6 : C 7 : C WRK52, WRK53, WRK54 pointeur actif a l'entree 8 : C NVARI nombre de variables internes (doit etre egal a ?) 9 : C 10 : C ISTEP flag utilise pour separer les etapes dans un calcul non local 11 : C ISTEP=0 -----> calcul local 12 : C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils 13 : C ISTEP=2 -----> calcul non local etape 2 on continue le calcul 14 : C a partir des seuils moyennes 15 : C 16 : C variables en sortie 17 : C 18 : C VARINF variables internes finales 19 : C 20 : C SIGMAF contraintes finales 21 : C 22 : C 23 : IMPLICIT REAL*8(A-H,O-Z) 24 : C 25 : -INC CCOPTIO 26 : -INC DECHE 27 : C 28 : C 29 : DIMENSION XD0(3,3),DF(3,3),POSEPSI2(3,3),EPS33(3,3) 30 : DIMENSION SIGTILDE33(3,3),SIGDEV33(3,3),SIGD33(3,3),SIGTILDE(6) 31 : DIMENSION UNMOINSDVP(3),HM1(3,3),DFD(3,3),UNMOINSD(3,3),DFVP(3) 32 : DIMENSION ROT(3,3),ROTT(3,3),VECP33T(3,3),TRAV(3,3),VECP33(3,3) 33 : DIMENSION EPSIPOS(3),HM1DIAG(3),EPS(2,2),DELTAF22(2,2),EPSIPP(3) 34 : DIMENSION VPDF(2),VECPDF(2,2),EPSTF(6) 35 : C 36 : PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12) 37 : C 38 : SEGMENT IECOU 39 : * COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK, 40 : INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7, 41 : C INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM, NYKK, 42 : 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16, 43 : C . NYALF1,NYBET1,NYR, NYA, NYRHO,NSIGY, NNKX, NYKX, IND, 44 : 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24, 45 : C . NSOM, NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT, NCXMAT, 46 : 3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA, 47 : C . LTRAC, MFR, IELE, NHRM, NBNN, NBELEM,ICARA, 48 : 4 icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38, 49 : C . LW2, NDEF, NSTRSS,MFR1, NBGMAT,NELMAT,MSOUPA, 50 : 5 icow39,icow40,icow41,icow42,icow43,icow44 51 : C . NUMAT1,LENDO, NBBB, NNVARI,KERR1, MELEME 52 : INTEGER icow45,icow46,icow47,icow48,icow49,icow50, 53 : . icow51,icow52,icow53,icow54,icow55,icow56 54 : . icow57,icow58 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 Debut de l'intégration numerique 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 : PRINT*, 'IFAILURE=', IFAILURE 239 : C mise a jour de XD0 pour traiter IFAILURE=2 dans le meme pas 240 : DO I=1,3 241 : DO J=1,3 242 : XD0(I,J)=DF(I,J) 243 : ENDDO 244 : ENDDO 245 : TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) 246 : TRADF=XD0(1,1)+XD0(2,2)+XD0(3,3) 246 : CALL JACOB3(XD0,3,DFVP,ROT) 247 : ENDIF 248 : ENDIF 249 : C 250 : IF(IFAILURE.EQ.2) THEN 251 : C deux directions d'endom sont fixee a Dc 251 : CALL TRSPOD(ROT,3,3,ROTT) 252 : C 253 : DF33=TRADF-DEUX*Dc 254 : C print*, 'DF33=',DF33 255 : IF(DF33.LT.XZERO) THEN 256 : DF33=XPETIT 257 : ENDIF 258 : C On corrige l'endommagement 258 : CALL ZERO(DF,3,3) 259 : DF(1,1)=Dc 260 : DF(2,2)=Dc 261 : DF(3,3)=DF33 262 : C 262 : CALL MATMAT(DF,ROTT,3,3,3,TRAV) 262 : CALL MATMAT(ROT,TRAV,3,3,3,DF) 263 : C 264 : C On verifie la troisieme direction d'endommagement 265 : IF(DF33.GE.Dc)THEN 266 : IFAILURE=3 267 : C mise a jour de XD0 pour traiter IFAILURE=3 dans le même pas 268 : DO I=1,3 269 : DO J=1,3 270 : XD0(I,J)=DF(I,J) 271 : ENDDO 272 : ENDDO 273 : TRAD0=XD0(1,1)+XD0(2,2)+XD0(3,3) 273 : CALL JACOB3(XD0,3,DFVP,ROT) 274 : ENDIF 275 : ENDIF 276 : C 277 : IF(IFAILURE.EQ.3)THEN 278 : C On est totalement rompu, endommagement isotrope egal a DC 278 : CALL TRSPOD(ROT,3,3,ROTT) 279 : C On corrige l'endommagement 279 : CALL ZERO(DF,3,3) 280 : DF(1,1)=Dc 281 : DF(2,2)=Dc 282 : DF(3,3)=Dc 283 : C 283 : CALL MATMAT(DF,ROTT,3,3,3,TRAV) 283 : CALL MATMAT(ROT,TRAV,3,3,3,DF) 284 : ENDIF 285 : C 286 : C Fin de la correction d'endommagement 287 : C 288 : C 289 : DH=(DF(1,1)+DF(2,2)+DF(3,3))/3.D0 290 : C 291 : ENDIF 292 : C 293 : C Fin des differents cas (elasticite ou endommagement) 294 : C 295 : C Calcul de la contrainte finale 296 : C 297 : DO I=1,3 298 : DO J=1,3 299 : UNMOINSD(I,J)=-UN*DF(I,J) 300 : ENDDO 301 : ENDDO 302 : UNMOINSD(1,1)=UN-DF(1,1) 303 : UNMOINSD(2,2)=UN-DF(2,2) 304 : UNMOINSD(3,3)=UN-DF(3,3) 304 : CALL JACOB3(UNMOINSD,3,UNMOINSDVP,VECP33) 305 : HM1DIAG(1)=DSQRT(UNMOINSDVP(1)) 306 : HM1DIAG(2)=DSQRT(UNMOINSDVP(2)) 307 : HM1DIAG(3)=DSQRT(UNMOINSDVP(3)) 307 : CALL PRODT2(HM1,HM1DIAG,VECP33,3) 308 : C 309 : C contrainte effective 309 : CALL MATVE1(DDHOOK,EPSTF,NSTRS,NSTRS,SIGTILDE,2) 309 : CALL ENDOCB(SIGTILDE,SIGTILDE33,1,IFOUR) 309 : CALL MATMAT(HM1,SIGTILDE33,3,3,3,TRAV) 309 : CALL MATMAT(TRAV,HM1,3,3,3,SIGD33) 310 : XTRAV=XZERO 311 : DO I=1,3 312 : DO J=1,3 313 : XTRAV=XTRAV+UNMOINSD(I,J)*SIGTILDE33(I,J) 314 : ENDDO 315 : ENDDO 316 : DO I=1,3 317 : DO J=1,3 318 : SIGD33(I,J)=SIGD33(I,J)-(XTRAV/(3.D0*(UN-DH)))*UNMOINSD(I,J) 319 : ENDDO 320 : ENDDO 321 : TRASIGTILDE=SIGTILDE(1)+SIGTILDE(2)+SIGTILDE(3) 322 : C 323 : IF(TRASIGTILDE.GE.XZERO) THEN 324 : COEF=UN-etaT*DH 325 : IF(COEF.LE.XZERO) THEN 326 : COEF=UN-Dc 327 : ENDIF 328 : TRASIG=COEF*TRASIGTILDE 329 : ELSE 330 : TRASIG=(UN-etaC*DH)*TRASIGTILDE 331 : ENDIF 332 : C 333 : SIGF(1)=SIGD33(1,1)+TRASIG/3.d0 334 : SIGF(2)=SIGD33(2,2)+TRASIG/3.d0 335 : SIGF(3)=SIGD33(3,3)+TRASIG/3.d0 336 : SIGF(4)=SIGD33(1,2) 337 : SIGF(5)=SIGD33(1,3) 338 : SIGF(6)=SIGD33(2,3) 339 : C 340 : C 341 : C stockage des variables du probleme 342 : C 343 : VARF(2)=DF(1,1) 344 : VARF(3)=DF(2,2) 345 : VARF(4)=DF(3,3) 346 : VARF(5)=DF(1,2) 347 : VARF(6)=DF(1,3) 348 : VARF(7)=DF(2,3) 349 : VARF(8)=ROT(1,1) 350 : VARF(9)=ROT(1,2) 351 : VARF(10)=ROT(1,3) 352 : VARF(11)=ROT(2,1) 353 : VARF(12)=ROT(2,2) 354 : VARF(13)=ROT(2,3) 355 : VARF(14)=ROT(3,1) 356 : VARF(15)=ROT(3,2) 357 : VARF(16)=ROT(3,3) 358 : C 359 : C flag pour l'endommagement maximum 360 : VARF(17)=IFAILURE 361 : C 362 : C Deformations totales finales 363 : DO ISTRS=1,NSTRS 364 : VARF(17+ISTRS)=EPSTF(ISTRS) 365 : ENDDO 366 : C 367 : 2000 CONTINUE 368 : C 369 : C 370 : RETURN 371 : END 372 : C 373 : 374 : 375 : 376 : 377 :
© Cast3M 2003 - All rights reserved.
Disclaimer