1 : C SICTENS SOURCE KICH 09/09/11 21:18:20 6505 2 : SUBROUTINE SICTENS(WRK52,WRK53,WRK54,WRK22,IB,IGAU,NVARI,NBPGAU 3 : & ,necou,iecou) 4 : C Modèle ONERA pseudo-tensoriel pour le SICf/SiC tissé 5 : C variables en entree 6 : C 7 : C WRK0,KRK1,WRK5 pointeurs sur des segments de travail 8 : C 9 : C 10 : C WRK0: 11 : C XMAT(NMATT): tableau composantes materiaux 12 : C 13 : C WRK1: 14 : C SIG0(NSTRS) : contraintes au debut du pas 15 : C VAR0(NVARI) : variables internes au debut du pas 16 : C DEPST(NSTRS): increment de deformation totale 17 : C 18 : C WRK22: 19 : C XXE: coordonnees de l'element en double precison 20 : C 21 : C WTRAV: 22 : C TXR: cosinus directeurs des axes locaux de l'element massif 23 : C 24 : C 25 : C WRK5: 26 : C 27 : C 28 : C 29 : C NSTRS nombre de composantes dans les vecteurs des contraintes 30 : C et les vecteurs des deformations 31 : C 32 : C NVARI nombre de variables internes (doit etre egal a 3) 33 : C 34 : C NMATT nombre de constantes du materiau 35 : C 36 : C 37 : C 38 : C variables en sortie 39 : C 40 : C VARF variables internes finales dans WRK0 41 : C 42 : C SIGF contraintes finales dans WRK0 43 : C 44 : IMPLICIT INTEGER(I-N) 45 : IMPLICIT REAL*8(A-H,O-Z) 46 : -INC CCOPTIO 47 : -INC CCREEL 48 : -INC DECHE 49 : * 50 : * Commun NECOU utilisé dans ECOINC 51 : * 52 : SEGMENT NECOU 53 : 54 : INTEGER NCOURB,IPLAST,IT,IMAPLA,ISOTRO, 55 : . ITYP,IFOURB,IFLUAG, 56 : . ICINE,ITHER,IFLUPL,ICYCL,IBI, 57 : . JFLUAG,KFLUAG,LFLUAG, 58 : . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF 59 : ENDSEGMENT 60 : 61 : * 62 : * Commun IECOU: sert de fourre-tout pour les initialisations 63 : * d'entiers 64 : * 65 : SEGMENT IECOU 66 : * 67 : INTEGER NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK, 68 : . NYALF1,NYBET1,NYR,NYA,NYRHO,NSIGY,NNKX,NYKX,IND, 69 : . NSOM,NINV,NINCMA,NCOMP,JELEM,LEGAUS,INAT,NCXMAT, 70 : . LTRAC,MFRbi,IELE,NHRM,NBNNBI,NBELMB,ICARA, 71 : . LW2bi,NDEF,NSTRSS,MFR1,NBGMAT,NELMAT,MSOUPA, 72 : . NUMAT1,LENDO,NBBB,NNVARI,KERR1,MELEME 73 : INTEGER icow45,icow46,icow47,icow48,icow49,icow50, 74 : . icow51,icow52,icow53,icow54,icow55,icow56 75 : . icow57,icow58 76 : ENDSEGMENT 77 : 78 : SEGMENT WRK22 79 : REAL*8 XXE(3,NBNNbi) 80 : ENDSEGMENT 81 : 82 : DIMENSION ZIDAUX(6,6) 83 : DIMENSION H1 (6,6), H2(6,6), H3(6,6), H4(6,6), H5(6,6) 84 : DIMENSION H01 (6,6), H02(6,6), H03(6,6), H04(6,6), H05(6,6) 85 : DIMENSION SEFF(6,6),CEFF(6,6) 86 : * 87 : * DORTH : DEFO DAN S LE REPÈRE OrTHO 88 : * sigfv : CONTRAINTE RÉORGANISÉES ORDRE CASTEM 89 : * DEFI : DÉFORMATIONS À LA FIN DU PAS PRÉCEDENT 90 : * depst : INCREMENT DE DÉFO DANS LE REPÈRE GÉNÉ 91 : * depstV : INCREMENT DE DÉFO DANS LE REPÈRE GÉNÉ 92 : * epsel : DEFO ELASTIQUE CALCULÉE POUR L'OBTENTION DES CONTRAINTES 93 : DIMENSION SORTH(6), SIGFV(6),EPSIN(6),EPSINV(6), 94 : & DEFI(6), DEPSTV(6), DORTH(6), EPSEL(6),EPSTFV(6), 95 : & SORTH1(6),SORTH2(6) 96 : * defoth : DEFO DANS LE REPERE ORTHO SOus FORME (3X3) 97 : * pOUR LE CALCUL DES DEFO PROPRES PUIS PARTIE >0 DES DÉFO PROPRES 98 : * xprop : DIRECTIONS PROPRES/REPERE ortho 99 : * xpinv : MATRICE INVERSE DE CHANGEMENT DE REPERE XPROP 100 : * xPINVT : MATRICE TRANSPOSÉE DE XPINV 101 : * defotht : MATRICE TEMPORAIRE POUR LE CHANGEMENT DE REPÈRE 102 : * epsclo : DEFO DE FERMETURE (CLOSURE) 103 : * EPSBAR : DEFO - DEFO CLOSURE 104 : DIMENSION DEFOTH(3,3),PROP(3),XPROP(3,3),XPINV(3,3), 105 : & DEFOTHT(3,3),XPINVT(3,3),EPSCLO(6), 106 : & EPSBAR(6),EPSIM(6),SIGPLU(6) 107 : DIMENSION ZK1(6,6),ZK2(6,6),ZK3(6,6),ZK4(6,6),ZK5(6,6) 108 : & ,TERME(6,6) 109 : DIMENSION EPSS1(6),EPSR1(6) 110 : * 111 : * DRIVING FORCES NORMALES ET TRANSVERSALES 112 : * dom : VARIABLES D'ENDOMAGEMENT 113 : DIMENSION Y(5),DOM(5) 114 : * DEPSIF (DELTA EPSILON I f : 115 : * LA FERMETURE COMMENCE POUR epsbar = DELTA EPSILON 116 : * finit POUR epsbar = DELTA EPSILON 117 : * ZNUACT INDICE D'ACTIVATION =1 : ACTIVE 118 : * =0 : PASSIVE 119 : DIMENSION DEPSIF(5),ZNUACT(5) 120 : * 121 : * Pseudo contraintes pour le calcul des 122 : * incréments de defo residuelles EPSR1 et stockées EPS1 123 : DIMENSION SOMS(6,6),SOMR(6,6), 124 : $ SISIS(6),SISIR(6) 125 : DIMENSION Y0(5),AIF(5),PY(5), 126 : $ YC(5),DC(5),DEPSIF0(5) 127 : DATA ( ( H1(I,J), J=1,6), I=1,6) 128 : & /1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 129 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 130 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 131 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 132 : & 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 133 : & 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/ 134 : DATA ( ( H2(I,J), J=1,6), I=1,6) 135 : & /0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 136 : & 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 137 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 138 : & 0.0, 0.0, 0.0, 1., 0.0, 0.0, 139 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 140 : & 0.0, 0.0, 0.0, 0.0, 0.0, 1./ 141 : * 142 : DATA ( ( H3(I,J), J=1,6), I=1,6) 143 : & /0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 144 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 145 : & 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 146 : & 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 147 : & 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 148 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/ 149 : * + 45° 150 : DATA ( ( H4(I,J), J=1,6), I=1,6) 151 : & /1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 152 : & 0.0, 1.0, 0.0, 0.0, 0.0, 0.5, 153 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 154 : & 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 155 : & 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 156 : & 0.5, 0.5, 0.0, 0.0, 0.0, 1.0/ 157 : * -45° 158 : DATA ( ( H5(I,J), J=1,6), I=1,6) 159 : & /1.0, 0.0, 0.0, 0.0, 0.0, -0.5, 160 : & 0.0, 1.0, 0.0, 0.0, 0.0, -0.5, 161 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 162 : & 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 163 : & 0.0, 0.0, 0.0, -1.0, 1.0, 0.0, 164 : & -0.5, -0.5, 0.0, 0.0, 0.0, 1.0/ 165 : 166 : * 167 : * Attention on rentre les racines carrée des valeurs 168 : * Y0 et YC telles que définies formules (2.21) et (2.24) 169 : * dans SEMT/LM2S/RT/05-034 170 : * PHENOMENOLOGICAL COEFFICIENTS 171 : * Coefficients pour les matrices H 172 : * 173 : * H1N, H1HP, H1P 174 : * H2N, H2HP, H2P 175 : * H3N, H3HP, H3P 176 : * H4N, H4HP, H4P 177 : * 178 : * scalar damage indicators 179 : * DOM1,DOM2,DOM3,DOM4,DOM5 180 : * YGR1, YGR2, YGR3, YW(6) 181 : * YEQ1, YEQ2, YEQ3 182 : * 183 : ***** INITIALIZATION VARIABLES 184 : * 185 : PI = 4.D0*DATAN(1.D0) 186 : * Proprietes materiau 187 : YG1 = XMAT(1) 188 : YG2 = XMAT(2) 189 : YG3 = XMAT(3) 190 : XNU12 = XMAT(4) 191 : XNU23 = XMAT(5) 192 : XNU13 = XMAT(6) 193 : G12 = XMAT(7) 194 : G23 = XMAT(8) 195 : G13 = XMAT(9) 196 : * Nombre de IIPARAM obligatoires pour le modèle 15 + ALP1 ALP2 ALP3 et RHO =19 197 : NOBL = 19 198 : * On extrait les valeurs facultatives du modèle 199 : H1N = XMAT(NOBL + 1) 200 : IF(IVALMA(NOBL+1).EQ.0.) H1N = 1. 201 : H1HP = XMAT(NOBL + 2) 202 : IF(IVALMA(NOBL+2).EQ.0.) H1HP = 0.7 203 : H1P = XMAT(NOBL + 3) 204 : IF(IVALMA(NOBL+3).EQ.0.) H1P = 0.45 205 : H2N = XMAT(NOBL + 4) 206 : IF(IVALMA(NOBL+4).EQ.0.) H2N = 1. 207 : H2HP = XMAT(NOBL + 5) 208 : IF(IVALMA(NOBL+5).EQ.0.) H2HP = 0.7 209 : H2P = XMAT(NOBL + 6) 210 : IF(IVALMA(NOBL+6).EQ.0.) H2P = 0.45 211 : H3N = XMAT(NOBL+7) 212 : IF(IVALMA(NOBL+7).EQ.0.) H3N = 1. 213 : H3P = XMAT(NOBL + 8) 214 : IF(IVALMA(NOBL+8).EQ.0.) H3P = 0.7 215 : H4N = XMAT(NOBL + 9) 216 : IF(IVALMA(NOBL+9).EQ.0.) H4N = 1. 217 : H4HP = XMAT(NOBL + 10) 218 : IF(IVALMA(NOBL+10).EQ.0.) H4HP = 0.7 219 : H4P = XMAT(NOBL + 11) 220 : IF(IVALMA(NOBL+11).EQ.0.) H4P = 1.2 221 : Y0(1) = XMAT(NOBL + 12) 222 : IF(IVALMA(NOBL+12).EQ.0.) Y0(1) =173.205 223 : Y0(2) = XMAT(NOBL + 13) 224 : IF(IVALMA(NOBL+13).EQ.0.) Y0(2) =173.205 225 : Y0(3) = XMAT(NOBL + 14) 226 : IF(IVALMA(NOBL+14).EQ.0.) Y0(3) =173.205 227 : Y0(4) = XMAT(NOBL + 15) 228 : IF(IVALMA(NOBL+15).EQ.0.) Y0(4) =100.00 229 : Y0(5) = XMAT(NOBL + 16) 230 : IF(IVALMA(NOBL+16).EQ.0.) Y0(5) =100.00 231 : YC(1) = XMAT(NOBL + 17) 232 : IF(IVALMA(NOBL+17).EQ.0.) YC(1) =1870.83 233 : YC(2) = XMAT(NOBL + 18) 234 : IF(IVALMA(NOBL+18).EQ.0.) YC(2) =1870.83 235 : YC(3) = XMAT(NOBL + 19) 236 : IF(IVALMA(NOBL+19).EQ.0.) YC(3) =1870.83 237 : YC(4) = XMAT(NOBL + 20) 238 : IF(IVALMA(NOBL+20).EQ.0.) YC(4) =3464.1 239 : YC(5) = XMAT(NOBL + 21) 240 : IF(IVALMA(NOBL+21).EQ.0.) YC(5) =3464.1 241 : DC(1) = XMAT(NOBL + 22) 242 : IF(IVALMA(NOBL+22).EQ.0.) DC(1) =4. 243 : DC(2) = XMAT(NOBL + 23) 244 : IF(IVALMA(NOBL+23).EQ.0.) DC(2) =4. 245 : DC(3) = XMAT(NOBL + 24) 246 : IF(IVALMA(NOBL+24).EQ.0.) DC(3) =4. 247 : DC(4) = XMAT(NOBL + 25) 248 : IF(IVALMA(NOBL+25).EQ.0.) DC(4) =4. 249 : DC(5) = XMAT(NOBL + 26) 250 : IF(IVALMA(NOBL+26).EQ.0.) DC(5) =4. 251 : PY(1) = XMAT(NOBL + 27) 252 : IF(IVALMA(NOBL+27).EQ.0.) PY(1) =1. 253 : PY(2) = XMAT(NOBL + 28) 254 : IF(IVALMA(NOBL+28).EQ.0.) PY(2) =1. 255 : PY(3) = XMAT(NOBL + 29) 256 : IF(IVALMA(NOBL+29).EQ.0.) PY(3) =1. 257 : PY(4) = XMAT(NOBL + 30) 258 : IF(IVALMA(NOBL+30).EQ.0.) PY(4) =1.2 259 : PY(5) = XMAT(NOBL + 31) 260 : IF(IVALMA(NOBL+31).EQ.0.) PY(5) =1.2 261 : B1 = XMAT(NOBL + 32) 262 : IF(IVALMA(NOBL+32).EQ.0.) B1 =1. 263 : B2 = XMAT(NOBL + 33) 264 : IF(IVALMA(NOBL+33).EQ.0.) B2 =1. 265 : B3 = XMAT(NOBL + 34) 266 : IF(IVALMA(NOBL+34).EQ.0.) B3 =1. 267 : DELTAL = XMAT(NOBL + 35) 268 : IF(IVALMA(NOBL+35).EQ.0.) DELTAL =0. 269 : TZERO = XMAT(NOBL + 36) 270 : IF(IVALMA(NOBL+36).EQ.0.) TZERO =0. 271 : DEPSIF0(1) = XMAT(NOBL + 37) 272 : IF(IVALMA(NOBL+37).EQ.0.) DEPSIF0(1) =0.0003 273 : DEPSIF0(2) = XMAT(NOBL + 38) 274 : IF(IVALMA(NOBL+38).EQ.0.) DEPSIF0(2) =0.0003 275 : DEPSIF0(3) = XMAT(NOBL + 39) 276 : IF(IVALMA(NOBL+39).EQ.0.) DEPSIF0(3) =0.0003 277 : DEPSIF0(4) = XMAT(NOBL + 40) 278 : IF(IVALMA(NOBL+40).EQ.0.) DEPSIF0(4) =0.0005 279 : DEPSIF0(5) = XMAT(NOBL + 41) 280 : IF(IVALMA(NOBL+41).EQ.0.) DEPSIF0(5) =0.0005 281 : AIF(1) = XMAT(NOBL + 42) 282 : IF(IVALMA(NOBL+42).EQ.0.) AIF(1) =0.5 283 : AIF(2) = XMAT(NOBL + 43) 284 : IF(IVALMA(NOBL+43).EQ.0.) AIF(2) =0.5 285 : AIF(3) = XMAT(NOBL + 44) 286 : IF(IVALMA(NOBL+44).EQ.0.) AIF(3) =0.5 287 : AIF(4) = XMAT(NOBL + 45) 288 : IF(IVALMA(NOBL+45).EQ.0.) AIF(4) =1. 289 : AIF(5) = XMAT(NOBL + 46) 290 : IF(IVALMA(NOBL+46).EQ.0.) AIF(5) =1. 291 : ETA1 = XMAT(NOBL + 47) 292 : IF(IVALMA(NOBL+47).EQ.0.) ETA1 =0.1 293 : ETA2 = XMAT(NOBL + 48) 294 : IF(IVALMA(NOBL+48).EQ.0.) ETA2 =0.1 295 : ETA3 = XMAT(NOBL + 49) 296 : IF(IVALMA(NOBL+49).EQ.0.) ETA3 =0.0 297 : ETA4 = XMAT(NOBL + 50) 298 : IF(IVALMA(NOBL+50).EQ.0.) ETA4 =0.1 299 : ETA5 = XMAT(NOBL + 51) 300 : IF(IVALMA(NOBL+51).EQ.0.) ETA5 =0.1 301 : 302 : * Pour debuggage 303 : IIIPARAM=0 304 : * Variables internes du modele 305 : DOM1=VAR0(2) 306 : DOM2=VAR0(3) 307 : DOM3=VAR0(4) 308 : DOM4=VAR0(5) 309 : DOM5=VAR0(6) 310 : * 311 : ****** CALCUL DE LA DEFORMATION INITIALE A PARTIR 312 : ****** DES CONTRAINTES INITIALES 313 : 314 : * 315 : IF (IIIPARAM.EQ.1) THEN 316 : WRITE(IOIMP,66770) IB,IGAU 317 : 66770 format(////////2x,'element ',i6,2x,'point ',i3//) 318 : 319 : WRITE (IOIMP,*) 'Increment des deformations ' 320 : WRITE (IOIMP,99999) (DEPST(ILOOP), ILOOP=1,6) 321 : 322 : WRITE (IOIMP,*) 'Contraintes au debut de l''iteration' 323 : WRITE (IOIMP,99999) (SIG0(ILOOP), ILOOP=1,6) 324 : ENDIF 325 : 326 : * Controle si il faut calculer la matrice de hook 327 : 328 : IF (IB.EQ.1.AND.IGAU.EQ.1) THEN 329 : GOTO 100 330 : ELSEIF (N2PTEL.EQ.1.AND.N2EL.EQ.1) THEN 331 : GOTO 200 332 : ELSEIF (N2PTEL.EQ.1.AND.N2EL.NE.1) THEN 333 : IF (IGAU.EQ.1) THEN 334 : GOTO 100 335 : ELSE 336 : GOTO 200 337 : ENDIF 338 : ELSE 339 : GOTO 100 340 : ENDIF 341 : 342 : 100 CONTINUE 343 : IPERR=1 344 : * Calcul de la matrice de hook pour un materiau orthotrope: C0 344 : CALL ZERO (DDAUX,LHOOK,LHOOK) 345 : XNU21=(YG2/YG1)*XNU12 346 : XNU32=(YG3/YG2)*XNU23 347 : XNU31=(YG3/YG1)*XNU13 348 : AUX=(1.D0-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 349 : & -2.D0*XNU21*XNU32*XNU13) 350 : AUX1=AUX/YG1 351 : AUX2=AUX/YG2 352 : AUX3=AUX/YG3 353 : DDAUX(1,1)=(1.D0-XNU23*XNU32)/AUX1 354 : DDAUX(1,2)=(XNU21+XNU31*XNU23)/AUX1 355 : DDAUX(2,1)=DDAUX(1,2) 356 : DDAUX(1,3)=(XNU31+XNU21*XNU32)/AUX1 357 : DDAUX(3,1)=DDAUX(1,3) 358 : DDAUX(2,2)=(1.D0-XNU13*XNU31)/AUX2 359 : DDAUX(2,3)=(XNU32+XNU12*XNU31)/AUX2 360 : DDAUX(3,2)=DDAUX(2,3) 361 : DDAUX(3,3)=(1.D0-XNU12*XNU21)/AUX3 362 : DDAUX(4,4)=G23 363 : DDAUX(5,5)=G13 364 : DDAUX(6,6)=G12 365 : * On recopie la matrice de hook pour l'inverser 366 : 200 DO 205 ILOOP=1,LHOOK 367 : DO 205 JLOOP=1,LHOOK 368 : ZIDAUX(ILOOP,JLOOP)= DDAUX(ILOOP,JLOOP) 369 : 205 CONTINUE 370 : * Inversion de la matrice de hook 371 : TPREC= 0.D1 372 : IPERR=0 372 : CALL INVALM(ZIDAUX,LHOOK,LHOOK,IPERR,TPREC) 373 : IF (IPERR.NE.0) THEN 374 : WRITE (IOIMP,*) 'ERREUR DANS L''INVERSION DE LA MATRICE DE HOOK' 375 : ENDIF 376 : * Calcul des tenseurs d'endommagement H01, H02, H03 377 : * selon les equations 4.2.3 à 4.2.5 377 : CALL ZERO (H01,LHOOK,LHOOK) 377 : CALL ZERO (H02,LHOOK,LHOOK) 377 : CALL ZERO (H03,LHOOK,LHOOK) 377 : CALL ZERO (H04,LHOOK,LHOOK) 377 : CALL ZERO (H05,LHOOK,LHOOK) 378 : 379 : H01(1,1)=H1(1,1)*H1N*ZIDAUX(1,1) 380 : H01(5,5)=H1(5,5)*H1HP*ZIDAUX(5,5) 381 : H01(6,6)=H1(6,6)*H1P*ZIDAUX(6,6) 382 : H02(2,2)=H2(2,2)*H2N*ZIDAUX(2,2) 383 : H02(4,4)=H2(4,4)*H2HP*ZIDAUX(4,4) 384 : H02(6,6)=H2(6,6)*H2P*ZIDAUX(6,6) 385 : H03(3,3)=H3(3,3)*H3N*ZIDAUX(3,3) 386 : H03(4,4)=H3(4,4)*H3P*ZIDAUX(4,4) 387 : H03(5,5)=H3(5,5)*H3P*ZIDAUX(5,5) 388 : * direction 4 : + 45 5 :-45 389 : * 1 ere ligne 390 : H04(1,1)=H4(1,1)*H4N*ZIDAUX(1,1) 391 : H05(1,1)=H5(1,1)*H4N*ZIDAUX(1,1) 392 : H04(1,6)=H4(1,6)*H4P*ZIDAUX(6,6) 393 : H05(1,6)=H5(1,6)*H4P*ZIDAUX(6,6) 394 : * 2 eme ligne 395 : H04(2,2)=H4(2,2)*H4N*ZIDAUX(1,1) 396 : H05(2,2)=H5(2,2)*H4N*ZIDAUX(1,1) 397 : H04(2,6)=H4(2,6)*H4P*ZIDAUX(6,6) 398 : H05(2,6)=H5(2,6)*H4P*ZIDAUX(6,6) 399 : * 4 eme ligne 400 : H04(4,4)=H4(4,4)*H4HP*ZIDAUX(5,5) 401 : H05(4,4)=H5(4,4)*H4HP*ZIDAUX(5,5) 402 : H04(4,5)=H4(4,5)*H4HP*ZIDAUX(5,5) 403 : H05(4,5)=H5(4,5)*H4HP*ZIDAUX(5,5) 404 : * 5 eme ligne 405 : H04(5,4)=H4(5,4)*H4HP*ZIDAUX(5,5) 406 : H05(5,4)=H5(5,4)*H4HP*ZIDAUX(5,5) 407 : H04(5,5)=H4(5,5)*H4HP*ZIDAUX(5,5) 408 : H05(5,5)=H5(5,5)*H4HP*ZIDAUX(5,5) 409 : * 6 eme ligne 410 : H04(6,1)=H4(6,1)*H4P*ZIDAUX(6,6) 411 : H05(6,1)=H5(6,1)*H4P*ZIDAUX(6,6) 412 : H04(6,2)=H4(6,2)*H4P*ZIDAUX(6,6) 413 : H05(6,2)=H5(6,2)*H4P*ZIDAUX(6,6) 414 : H04(6,6)=H4(6,6)*H4P*ZIDAUX(6,6) 415 : H05(6,6)=H5(6,6)*H4P*ZIDAUX(6,6) 416 : * 417 : IF (IIIPARAM.EQ.1) THEN 418 : WRITE (IOIMP,*) 'Matrice H01' 419 : DO 250 ILOOP= 1,LHOOK 420 : WRITE (IOIMP, 99999 ) (H01(ILOOP,J), J=1,LHOOK) 421 : 250 CONTINUE 422 : WRITE (IOIMP,*) 'Matrice H02' 423 : DO 251 ILOOP= 1,LHOOK 424 : WRITE (IOIMP, 99999 ) (H02(ILOOP,J), J=1,LHOOK) 425 : 251 CONTINUE 426 : WRITE (IOIMP,*) 'Matrice H03' 427 : DO 252 ILOOP= 1,LHOOK 428 : WRITE (IOIMP, 99999 ) (H03(ILOOP,J), J=1,LHOOK) 429 : 252 CONTINUE 430 : ENDIF 431 : ***** Recuperation deformation initiale dans le repere orth. 432 : ***** EN plus les VI sont rangées dans l'ordre de l'ONERA cad 433 : ****** 12 ==> 6, 13 ==> 5 et 23 ==> 4 434 : DEFI(1)=VAR0(7) 435 : DEFI(2)=VAR0(8) 436 : DEFI(3)=VAR0(9) 437 : DEFI(4)=VAR0(10) 438 : DEFI(5)=VAR0(11) 439 : DEFI(6)=VAR0(12) 440 : ** repère global 441 : IF (IIIPARAM.EQ.1) THEN 442 : WRITE (IOIMP,*) 'DEFORMATION RECUPEREE' 443 : WRITE (IOIMP,99999) (DEFI(ILOOP), ILOOP=1,6) 444 : ENDIF 445 : * Calcul de la nouvelle défo repère glob ordre cast 446 : DO 451 ILOOP=1,6 447 : EPSTFV(ILOOP)=EPST0(ILOOP)+DEPST(ILOOP) 448 : * EPSTF(ILOOP)=EPSTFV(ILOOP) 449 : 451 CONTINUE 450 : * WRITE (IOIMP,*) 'DEF IMPOSEE glob cast3m' 451 : * WRITE (IOIMP,99999) (EPSTFV(ILOOP), ILOOP=1,6) 452 : * Par contre les depstv cad après les DORTH sont rangés dans l'ordre castem 453 : * : 12 ==> 4 13 ==> 5 23 ==>6 - 453 : CALL CICROT (wrk52,wrk53,wrk54,1,EPSTFV,DORTH) 454 : * 455 : * On les réorganise donc pour etre dans l'ordre ONERA 456 : VARTMP= DORTH(4) 457 : DORTH(4)= DORTH(6) 458 : DORTH(6)= VARTMP 459 : * 460 : * Calcul de la nouvelle deformation totale dans le repère ortho 461 : * Utilisées pour le calcul des valeurs propres 462 : DEFOTH(1,1) = DORTH(1) 463 : DEFOTH(2,2) = DORTH(2) 464 : DEFOTH(3,3) = DORTH(3) 465 : DEFOTH(1,2) = DORTH(6) 466 : DEFOTH(1,3) = DORTH(5) 467 : DEFOTH(2,1) = DEFOTH(1,2) 468 : DEFOTH(2,3) = DORTH(4) 469 : DEFOTH(3,1) = DEFOTH(1,3) 470 : DEFOTH(3,2) = DEFOTH(2,3) 471 : * 471 : CALL JACOBA(DEFOTH,3,PROP,XPROP) 472 : C 473 : C partie positive des valeurs propres 474 : IF(PROP(1).LT.0.) PROP(1)=0.D0 475 : IF(PROP(2).LT.0.) PROP(2)=0.D0 476 : IF(PROP(3).LT.0.) PROP(3)=0.D0 477 : C Je remet ma matrice de valeurs propres corrigée axes ortho 478 : DEFOTH(1,1) = PROP(1) 479 : DEFOTH(2,2) = PROP(2) 480 : DEFOTH(3,3) = PROP(3) 481 : DEFOTH(1,2) = 0.D0 482 : DEFOTH(1,3) = 0.D0 483 : DEFOTH(2,1) = 0.D0 484 : DEFOTH(2,3) = 0.D0 485 : DEFOTH(3,1) = 0.D0 486 : DEFOTH(3,2) = 0.D0 487 : C XPROP(3,3) est la matrice donnant les coordonnées des vecteurs propre 488 : C par rapport aux axes d'orthotropie 489 : C 1 er indice direction 2° n° vecteur propre 490 : C j'ai besoin de la matrice donnant les axes ortho dans le repère des vecteurs propres 491 : C PRINT *,'AVANT INVER3' 491 : CALL INVER3(XPROP,XPINV) 492 : C PRINT *,'APRES INVER3' 492 : CALL TRSPOD(XPINV,3,3,XPINVT) 492 : CALL MULMAT(DEFOTHT,DEFOTH,XPINV,3,3,3) 492 : CALL MULMAT(DEFOTH,XPINVT,DEFOTHT,3,3,3) 493 : C 494 : C en principe on a la matrice de defo efficasse dans les axes ortho defoth 495 : C calcul des driving forces dans le cas scalaire 496 : C mode normal 497 : C calcul des driving forces dans le cas pseudo-tensoriel 498 : 499 : Y1 =0.5*(DDAUX(1,1)*DEFOTH(1,1)*DEFOTH(1,1)+b1* 500 : & DDAUX(6,6)*DEFOTH(1,2)*DEFOTH(1,2)+b2* 501 : & DDAUX(5,5)*DEFOTH(1,3)*DEFOTH(1,3)) 502 : Y2 =0.5*(DDAUX(2,2)*DEFOTH(2,2)*DEFOTH(2,2)+b1* 503 : & DDAUX(6,6)*DEFOTH(1,2)*DEFOTH(1,2)+b2* 504 : & DDAUX(4,4)*DEFOTH(2,3)*DEFOTH(2,3)) 505 : Y6 =0.25*(DDAUX(1,1)*DEFOTH(1,1)*DEFOTH(1,2)+ 506 : & DDAUX(2,2)*DEFOTH(2,2)*DEFOTH(1,2)+b1* 507 : & DDAUX(6,6)*DEFOTH(1,2)* 508 : & (DEFOTH(1,1)+DEFOTH(2,2)) ) 509 : C driving forces utilisées pour le calcul des dommages 510 : Y1 = Y1 - (ABS(Y6)) 511 : Y2 = Y2 - (ABS(Y6)) 512 : IF(Y6.GT.0.) THEN 513 : Y(4)=DSQRT(Y6) 514 : Y(5)= 0. 515 : ELSE 516 : Y(4)=0. 517 : Y(5)=DSQRT(-Y6) 518 : ENDIF 519 : Y3 = 0.5*(DDAUX(3,3)*DEFOTH(3,3)*DEFOTH(3,3)+b3* 520 : & DDAUX(4,4)*DEFOTH(2,3)*DEFOTH(2,3)+b3* 521 : & DDAUX(5,5)*DEFOTH(1,3)*DEFOTH(1,3)) 522 : C 523 : 524 : IF(Y1.LT.0.) THEN 525 : Y(1)=0. 526 : ELSE 527 : Y(1)=DSQRT(Y1) 528 : ENDIF 529 : IF(Y2.LT.0.) THEN 530 : Y(2)=0. 531 : ELSE 532 : Y(2)=DSQRT(Y2) 533 : ENDIF 534 : IF(Y3.LT.0.) THEN 535 : Y(3)=0. 536 : ELSE 537 : Y(3)=DSQRT(Y3) 538 : ENDIF 539 : C Calcul des variables d'endomagenement 540 : Do 722 IDRF=1,5 541 : FAFA =Y(IDRF)-Y0(IDRF) 542 : IF (FAFA.LT.0.) FAFA=0. 543 : DOM(IDRF)= DC(IDRF)*(1.D0 - DEXP(-1.* 544 : & ((FAFA/YC(IDRF))**PY(IDRF)))) 545 : 722 CONTINUE 546 : DOM1= MAX(DOM(1),VAR0(2)) 547 : DOM2= MAX(DOM(2),VAR0(3)) 548 : DOM3= MAX(DOM(3),VAR0(4)) 549 : DOM4= MAX(DOM(4),VAR0(5)) 550 : DOM5= MAX(DOM(5),VAR0(6)) 551 : C dans WRK52 (DECHE ) on a tur0(1) et turf(1) 552 : 553 : IF (IIIPARAM.EQ.1) THEN 554 : WRITE (IOIMP,*) 'Increment des deformations orth. ' 555 : WRITE (IOIMP,99999) (DORTH(ILOOP), ILOOP=1,6) 556 : ENDIF 557 : * 558 : DO 500 ILOOP=1,NSTRSS 559 : * On garde les déformations totales dans les var. internes 560 : VARF(ILOOP+6)=DORTH(ILOOP) 561 : 500 CONTINUE 562 : * 563 : * On peut calculer EPSILON CLosure 564 : * EN fonction de deltal et tzero 565 : * et de la temperature TURE0(1)ou TUREF(1)du segment WRK52 de l'include DECHE 566 : * 567 : DO 701 IDOU=1,6 568 : EPSCLO(IDOU) = DELTAL*(TURE0(1)-TZERO) 569 : EPSBAR(IDOU) = DORTH(IDOU)- EPSCLO(IDOU) 570 : EPSIM(IDOU) = (DORTH(IDOU)+VAR0(IDOU+6))/2. 571 : 701 CONTINUE 572 : * On redefini les epsbar pour les direction 45 et -45 ° (cf 4.2.8) 573 : EPSBAR(4) = (EPSBAR(1)+EPSBAR(2)+(2.D0*EPSBAR(6)))*0.5 574 : EPSBAR(5) = (EPSBAR(1)+EPSBAR(2)-(2.D0*EPSBAR(6)))*0.5 575 : 576 : * Calcul des valeurs délimitant les déformations pour 577 : * lesquelles les fissures commencent et finissent à/de se fermer 578 : * en vue du calcul de l'index de désactivation 579 : DEPSIF(1)= (1.D0 + (AIF(1)*DOM1))*DEPSIF0(1) 580 : DEPSIF(2)= (1.D0 + (AIF(2)*DOM2))*DEPSIF0(2) 581 : DEPSIF(3)= (1.D0 + (AIF(3)*DOM3))*DEPSIF0(3) 582 : DEPSIF(4)= (1.D0 + (AIF(4)*DOM4))*DEPSIF0(4) 583 : DEPSIF(5)= (1.D0 + (AIF(5)*DOM5))*DEPSIF0(5) 584 : * 585 : * Calcul de l'index d'activation 1 si actif 0 si passif 586 : * 587 : DO 703 IDO2=1,5 588 : IF(EPSBAR(IDO2).GE.DEPSIF(IDO2)) THEN 589 : ZNUACT(IDO2)= 1.D0 590 : ELSEIF(EPSBAR(IDO2).GT.(-1.*(DEPSIF(IDO2)))) THEN 591 : ZNUACT(IDO2)= 0.5D0*(1.D0-(DCOS(PI/2.D0*(EPSBAR(IDO2) 592 : & +DEPSIF(IDO2))/DEPSIF(IDO2)))) 593 : ELSE 594 : ZNUACT(IDO2)= 0.D0 595 : ENDIF 596 : 703 CONTINUE 597 : * 598 : * Calcul des déformations résiduelles et stockées 599 : * 600 : * Calcul des variations du dommage 601 : DELTDO1 = DOM1 - VAR0(2) 602 : DELTDO2 = DOM2 - VAR0(3) 603 : DELTDO3 = DOM3 - VAR0(4) 604 : DELTDO4 = DOM4 - VAR0(5) 605 : DELTDO5 = DOM5 - VAR0(6) 606 : * Calcul des variations d'indice d'activation 607 : DELTNU1 = ZNUACT(1) - VAR0(13) 608 : DELTNU2 = ZNUACT(2) - VAR0(14) 609 : DELTNU3 = ZNUACT(3) - VAR0(15) 610 : DELTNU4 = ZNUACT(4) - VAR0(16) 611 : DELTNU5 = ZNUACT(5) - VAR0(17) 612 : * 613 : * Calcul de Seff =somme sur i de nui di hi0 614 : * 615 : DO 744 INDJ3=1,6 616 : Do 745 INDJ4=1,6 617 : SEFF(INDJ3,INDJ4) = ZIDAUX(INDJ3,INDJ4)+ 618 : & ZNUACT(1)*DOM1*H01(INDJ3,INDJ4)+ 619 : & ZNUACT(2)*DOM2*H02(INDJ3,INDJ4)+ 620 : & ZNUACT(3)*DOM3*H03(INDJ3,INDJ4)+ 621 : & ZNUACT(4)*DOM4*H04(INDJ3,INDJ4)+ 622 : & ZNUACT(5)*DOM5*H05(INDJ3,INDJ4) 623 : CEFF(INDJ3,INDJ4) =SEFF(INDJ3,INDJ4) 624 : 625 : 745 CONTINUE 626 : 744 CONTINUE 627 : * inversion de la matrice de rigidité 628 : TPREC= 0.D1 629 : IPERR=0 629 : CALL INVALM(CEFF,LHOOK,LHOOK,IPERR,TPREC) 630 : IF (IPERR.NE.0) THEN 631 : WRITE (IOIMP,*) 'ERREUR DANS L''INVERSION DE LA MATRICE DE HOOK' 632 : ENDIF 632 : CALL MULMAT (TERME,H01,CEFF,6,6,6) 632 : CALL MULMAT (ZK1,CEFF,TERME,6,6,6) 632 : CALL MULMAT (TERME,H02,CEFF,6,6,6) 632 : CALL MULMAT (ZK2,CEFF,TERME,6,6,6) 632 : CALL MULMAT (TERME,H03,CEFF,6,6,6) 632 : CALL MULMAT (ZK3,CEFF,TERME,6,6,6) 632 : CALL MULMAT (TERME,H04,CEFF,6,6,6) 632 : CALL MULMAT (ZK4,CEFF,TERME,6,6,6) 632 : CALL MULMAT (TERME,H05,CEFF,6,6,6) 632 : CALL MULMAT (ZK5,CEFF,TERME,6,6,6) 633 : * 634 : * Déformation stockées (S) + résiduelle (R) à mettre dans les VI 635 : * 636 : DO 741 INDI=1,6 637 : SISIS(INDI)=0.D0 638 : SISIR(INDI)=0.D0 639 : DO 742 INDJ=1,6 640 : SOMS(INDI,INDJ)= DELTNU1*DOM1*ZK1(INDI,INDJ)+ 641 : $ DELTNU2*DOM2*ZK2(INDI,INDJ)+ 642 : $ DELTNU3*DOM3*ZK3(INDI,INDJ)+ 643 : $ DELTNU4*DOM4*ZK4(INDI,INDJ)+ 644 : $ DELTNU5*DOM5*ZK5(INDI,INDJ) 645 : SISIS(INDI)=SISIS(INDI)-SOMS(INDI,INDJ)*EPSIM(INDJ) 646 : SOMR(INDI,INDJ)= ETA1*ZNUACT(1)*DELTDO1*ZK1(INDI,INDJ)+ 647 : $ ETA2*ZNUACT(2)*DELTDO2*ZK2(INDI,INDJ)+ 648 : $ ETA3*ZNUACT(3)*DELTDO3*ZK3(INDI,INDJ)+ 649 : $ ETA4*ZNUACT(4)*DELTDO4*ZK4(INDI,INDJ)+ 650 : $ ETA5*ZNUACT(5)*DELTDO5*ZK5(INDI,INDJ) 651 : SISIR(INDI)=SISIR(INDI)+SOMR(INDI,INDJ)*EPSIM(INDJ) 652 : 653 : 742 CONTINUE 654 : 741 CONTINUE 655 : C 655 : CALL MULMAT (EPSR1,ZIDAUX,SISIR,6,1,6) 655 : CALL MULMAT (EPSS1,ZIDAUX,SISIS,6,1,6) 656 : DO 743 INDI=1,6 657 : EPSIN(INDI)=VAR0(17+INDI)+EPSR1(INDI)+EPSS1(INDI) 658 : EPSEL(INDI) = DORTH(INDI)-EPSIN(INDI) 659 : 743 CONTINUE 660 : * 661 : * SORTH contient les contraintes a la fin du pas 662 : * dans le repère orthotrope 662 : CALL MULMAT (SORTH1,CEFF,DORTH,6,1,6) 662 : CALL MULMAT (SORTH2,DDAUX,EPSIN,6,1,6) 663 : DO 749 INDIA=1,6 664 : SORTH(INDIA) = SORTH1(INDIA) - SORTH2(INDIA) 665 : 749 CONTINUE 666 : IF (IIIPARAM.EQ.1) THEN 667 : WRITE (IOIMP,*) 'Contraintes calculees (repere orth.)' 668 : WRITE (IOIMP,99999) (SORTH(ILOOP), ILOOP=1,6) 669 : ENDIF 670 : * On réorganise les contraintes en cisaillement pour 671 : * calculer les contraintes dans le repere global (ON REMET DANS L'ORDRE castem) 672 : SIGFV(1)=SORTH(1) 673 : SIGFV(2)=SORTH(2) 674 : SIGFV(3)=SORTH(3) 675 : SIGFV(4)=SORTH(6) 676 : SIGFV(5)=SORTH(5) 677 : SIGFV(6)=SORTH(4) 678 : EPSINV(1)=EPSIN(1) 679 : EPSINV(2)=EPSIN(2) 680 : EPSINV(3)=EPSIN(3) 681 : EPSINV(4)=EPSIN(6) 682 : EPSINV(5)=EPSIN(5) 683 : EPSINV(6)=EPSIN(4) 684 : VARF(1)= 0.D0 685 : DO 186 INDI = 1,6 686 : VARF(17+INDI)= EPSIN(INDI) 687 : VARF(1)= VARF(1)+ EPSIN(INDI) 688 : 186 CONTINUE 689 : * On appele SICROT pour trouver les contraintes dans le repère global 690 : * On utilise SORTH pour garder le resultat 691 : * (on peut pas passer SIGF à la subroutine car il est dans un segment) 691 : CALL CICROT (wrk52,wrk53,wrk54,0,SIGFV,SORTH) 692 : 693 : * On recopie SORTH in SIGF 694 : DO 1000 ILOOP=1,6 695 : SIGF(ILOOP)=SORTH(ILOOP) 696 : 1000 CONTINUE 697 : IF (IIIPARAM.EQ.1) THEN 698 : WRITE (IOIMP,*) 'SIG REP GLOB cast' 699 : WRITE (IOIMP,99999) (SIGF(ILOOP), ILOOP=1,6) 700 : ENDIF 701 : VARF(2)= DOM1 702 : VARF(3)= DOM2 703 : VARF(4)= DOM3 704 : VARF(5)= DOM4 705 : VARF(6)= DOM5 706 : VARF(13)= ZNUACT(1) 707 : VARF(14)= ZNUACT(2) 708 : VARF(15)= ZNUACT(3) 709 : VARF(16)= ZNUACT(4) 710 : VARF(17)= ZNUACT(5) 711 : 99999 format(2x,' ',(6(1x,3pe12.5),/)) 712 : 99998 format(2x,' ',(4(1x,3pe12.5),/)) 713 : 99997 format(2x,' ',(1x,3pe12.5),/) 714 : 99996 format(2x,' ',(3(1x,1pe12.3),/)) 715 : RETURN 716 : END 717 : 718 : 719 : 720 : 721 : 722 : 723 : 724 : 725 : 726 : 727 : 728 : 729 : 730 : 731 : 732 : 733 : 734 : 735 : 736 : 737 : 738 :
© Cast3M 2003 - All rights reserved.
Disclaimer