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