1 : C XXCRE1 SOURCE AM 09/12/15 21:15:17 6587 2 : SUBROUTINE XXCRE1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD, 3 : &DSPT,XMAT,XMAINF,XMASUP,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV, 4 : &DDINVp,YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI, 5 : &TPOINT,TINF,TSUP,ITEST,ITHHER,VARF,IB,IGAU,KERRE) 6 : C----------------------------------------------------------------------- 7 : C Objet: Cette subroutine calcule les derivees des variables internes 8 : C pour un materiau viscoplastique a endommagement et ecrouissage 9 : C isotrope en regime anisotherme de 2 manieres au choix suivant 10 : C la valeur de ITEST,ITHHER 11 : C----------------------------------------------------------------------- 12 : C 13 : C----------------------------------------------------------------------- 14 : C Entree: TAU pas de temps 15 : C TI temperature au debut du pas 16 : C TPOINT derivee de la temperature 17 : C TINF et TSUP bornes entre lesquelles TI est comprise 18 : C SIG(NSTRS) tenseur des contraintes 19 : C EPSV(NSTRS) tenseur des deformations 20 : C VAR(NVARI) tableau contenant les variables internes 21 : C pilotant les equations 22 : C VAR(1)=p ; VAR(2)=D11 ; VAR(3)=D22; VAR(4)=D33; 23 : C VAR(5)=D12; VAR(6)=D13; VAR(7)=D23; VAR(8)=DMAX 24 : C EPSTHD vitesse de dilatation thermique au debut du pas 25 : C DSPT(NSTRS,NSTRS) vitesse de deformation totale 26 : C XMAT(NCOMAT) tableau des parametres scalaires du materiau 27 : C a une temperature T donnee 28 : C XMAINF(NCOMAT) tableau des parametres scalaires du materiau 29 : C a la temperature TINF 30 : C XMASUP(NCOMAT) tableau des parametres scalaires du materiau 31 : C a la temperature TSUP 32 : C XMAT(1)=YOUNG ;XMAT(2)=XNU ;XMAT(3)=N 33 : C XMAT(4)=M ;XMAT(5)=KK; XMAT(6)=ALF2; 34 : C XMAT(7)=BET2 ;XMAT(8)=r; XMAT(9)=A ; 35 : C XMAT(10)=q ;XMAT(11)=ALPHAT 36 : C XMAT(12)=RHO; XMAT(13)=SIGY 37 : C MFR indice de la formulation mecanique(seulement massif 38 : C ou coque pour les materiaux endommageables 39 : C ICARA nombre de caracteristiques geometriques des elements 40 : C finis 41 : C XCAR(ICARA) tableau des caracteristiques geometriques des 42 : C elements finis 43 : C IFOURB= -2 EN CONTR. PLANES 44 : C -1 EN DEFORM. PLANES 45 : C 0 EN AXISYMETRIE 46 : C 1 EN SERIE DE FOURIER 47 : C 2 EN TRIDIM 48 : C INDIC=0, 1 OU -1 pour plasticite avec endommagement 49 : C =2 OU -2 pour viscoplasticite avec endommagement 50 : C ITEST = 0 pas uniformite des listes de temperatures pour tous les 51 : C coefficients non lineaires 52 : C d'ou interpolation sur les coefficients du materiau 53 : C = 1 uniformite des listes de temperatures pour tous les 54 : C coefficients non lineaires 55 : C d'ou moyenne ponderee sur la variable elle-meme 56 : C ex: D(T)=A(T)*K**B(T) 57 : C T=TETA*TINF+(1-TETA)*TSUP 58 : C si ITEST=0 on remplace A(T) par (TETA*A(TINF)+(1-TETA))*A(TSUP) 59 : C et B(T) par (TETA*B(TINF)+(1-TETA))*B(TSUP) 60 : C si ITEST=1 D(T)=TETA*D(TINF)+(1-TETA)*D(TSUP) 61 : C ITHHER = 0 pas de chargement thermique et materiau constant 62 : C = 1 chargement thermique et materiau constant 63 : C = 2 chargement thermique et materiau(T) 64 : C------------------------------------------------------------------------ 65 : C 66 : C------------------------------------------------------------------------ 67 : C Sortie: EPSVD(NSTRS) derivee du tenseur des deformations 68 : C VARD(NVARI) tableau contenant les derivees des variables 69 : C internes 70 : C VARF(NVARI) état des variables internes à la fin du sous programme 71 : C SIGD(NSTRS) derivee des contraintes 72 : C DD(NSTRS,NSTRS) matrice de Hooke au debut du pas 73 : C DDV(NSTRS,NSTRS) derivee de DD 74 : C DDINV(NSTRS,NSTRS) inverse de DD 75 : C------------------------------------------------------------------------ 76 : IMPLICIT INTEGER(I-N) 77 : IMPLICIT REAL*8 (A-H,O-Z) 78 : DIMENSION XMAINF(*),XMASUP(*),XCAR(*) 79 : DIMENSION SIG(*),EPSV(*),VAR(*) 80 : DIMENSION SIGD(*),EPSVD(*),VARD(*),VARF(NVARI) 81 : DIMENSION XMAT(*),DSPT(*),YOUNG(*),XNU(*),EPSTHD(*) 82 : DIMENSION DDV(NSTRS,NSTRS),DDINVp(NSTRS,NSTRS) 83 : DIMENSION DD(NSTRS,NSTRS),DDINV(NSTRS,NSTRS) 84 : DIMENSION VARD1(3),EPSVD1(6),SIGD1(6),VARD2(3),EPSVD2(6) 85 : DIMENSION EPSED(6),XX(6),Di(3),VaP(3),AA(3),BB(6) 86 : DIMENSION ROT(3,3),ROTT(3,3),CO(3,3),COM(3,3) 87 : DIMENSION D1(3,3),D2(3,3),DV1(3,3),DI1(3,3),DI2(3,3) 88 : DIMENSION AAA(3,3),BBB(3,3),CCC(3,3),VcP(3,3),FTH(3,3) 89 : PARAMETER (UN=1.D0) 90 : LOGICAL endonul 91 : C 92 : C----- le materiau est constant ou pas, on derive normalement les contraintes 93 : C 94 : IFLAG=0 95 : IF ( (ITEST.EQ.1).AND.(ITHHER.EQ.2) ) THEN 96 : C----- si le materiau depend de la temperature et si les listes sont identiques 97 : IF (TINF.NE.TSUP) THEN 98 : C-------------- Moyenne ponderee sur les vitesses des variables 99 : IFLAG=1 100 : TETA=(TI-TSUP)/(TINF-TSUP) 100 : CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD1,EPSVD1,VARD1,EPSTHD, 101 : & DSPT,XMAINF,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp, 102 : & YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC, 103 : & TINF,TPOINT,ITHHER,VARF,IB,IGAU,KERRRE) 104 : IF(KERRE.NE.0) GO TO 100 105 : 106 : DO 10 I=1,NVARI 107 : 10 VARD2(I)=TETA*VARD1(I) 108 : DO 20 I=1,NSTRS 109 : 20 EPSVD2(I)=TETA*EPSVD1(I) 109 : CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD1,EPSVD1,VARD1,EPSTHD, 110 : & DSPT,XMASUP,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp, 111 : & YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC, 112 : & TSUP,TPOINT,ITHHER,VARF,IB,IGAU,KERRRE) 113 : IF(KERRE.NE.0) GO TO 100 114 : 115 : DO 30 I=1,NVARI 116 : 30 VARD2(I)=VARD2(I)+(1.D0-TETA)*VARD1(I) 117 : DO 40 I=1,NSTRS 118 : 40 EPSVD2(I)=EPSVD2(I)+(1.D0-TETA)*EPSVD1(I) 119 : ENDIF 120 : ENDIF 121 : C 121 : CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD, 122 : &DSPT,XMAT,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,YOUNG, 123 : &NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI,TPOINT, 124 : &ITHHER,VARF,IB,IGAU,KERRE) 125 : IF(KERRE.NE.0) GO TO 100 126 : C 127 : IF (IFLAG.EQ.0) GOTO 100 128 : C 129 : C --- Cas ou ITEST = 1 et ITHHER = 2 130 : C 131 : DO I=1,NVARI 132 : VARD(I)=VARD2(I) 133 : ENDDO 134 : DO I=1,NSTRS 135 : EPSVD(I)=EPSVD2(I) 136 : ENDDO 137 : C 138 : endonul = .FALSE. 139 : IF ((VAR(2).EQ.0).AND.(VAR(3).EQ.0).AND.(VAR(4).EQ.0).AND. 140 : & (VAR(5).EQ.0).AND.(VAR(6).EQ.0).AND.(VAR(7).EQ.0).AND. 141 : & (VAR(8).EQ.0)) THEN 142 : endonul = .TRUE. 143 : ENDIF 144 : C 145 : C---------------------------------------------- 146 : C TENSEUR DES CONTRAINTES (dans repere general) 147 : C---------------------------------------------- 148 : C 148 : CALL ZERO(CO,3,3) 149 : IF (mfr.eq.1) THEN 150 : DO I=1,3 151 : CO(I,I) = SIG(I) 152 : ENDDO 153 : CO(1,2) = SIG(4) 154 : IF(NSTRS.NE.6) GOTO 50 155 : CO(1,3) = SIG(5) 156 : CO(2,3) = SIG(6) 157 : 50 CONTINUE 158 : CO(2,1) = CO(1,2) 159 : CO(3,1) = CO(1,3) 160 : CO(3,2) = CO(2,3) 161 : ENDIF 162 : C 163 : C------------------------------------ 164 : C TENSEUR D ENDOMMAGEMENT ET 165 : C TENSEUR D ENDOMMAGEMENT DIAGONALISE 166 : C------------------------------------ 167 : C 168 : XD11=VAR(2) 169 : XD22=VAR(3) 170 : XD33=VAR(4) 171 : XD12=VAR(5) 172 : XD13=VAR(6) 173 : XD23=VAR(7) 173 : CALL ZERO(D1,3,3) 174 : D1(1,1) = XD11 175 : D1(2,2) = XD22 176 : D1(3,3) = XD33 177 : D1(1,2) = XD12 178 : D1(1,3) = XD13 179 : D1(2,3) = XD23 180 : D1(2,1) = D1(1,2) 181 : D1(3,1) = D1(1,3) 182 : D1(3,2) = D1(2,3) 183 : C 183 : CALL JACOB3(D1,3,Di,ROT) 184 : C 184 : CALL ZERO(D2,3,3) 185 : D2(1,1) = Di(1) 186 : D2(2,2) = Di(2) 187 : D2(3,3) = Di(3) 188 : w1 = 1.D0 - Di(1) 189 : w2 = 1.D0 - Di(2) 190 : w3 = 1.D0 - Di(3) 191 : IF (D2(1,1).EQ.1.D0) w1=1.D-6 192 : C 193 : ADMAX = DMAX1(Di(1),Di(2),Di(3)) 194 : C 195 : C------------------------------------------------------- 196 : C TENSEUR DES CONTRAINTES MOTRICES (dans repere general) 197 : C------------------------------------------------------- 198 : C 199 : EQ=XMAT(10) 200 : DO I=1,3 201 : DO J=1,3 202 : ROTT(I,J)=ROT(I,J) 203 : ENDDO 204 : ENDDO 204 : CALL INVALM (ROTT, 3 , 3 , KER, 1.D-12) 205 : C 206 : IF (endonul) THEN 207 : DO I=1,3 208 : DO J=1,3 209 : COM(I,J) = CO(I,J) 210 : ENDDO 211 : ENDDO 212 : ELSE 213 : C --- Matrice de [(I-D)/det(I-D)]**(q/2) dans le repère principal de D --- 213 : CALL ZERO(DI1,3,3) 214 : DO I=1,3 215 : XNB=(1.D0-Di(I))/(w1*w2*w3) 216 : IF (XNB.GE.0.D0) THEN 217 : DI1(I,I)=(XNB)**(EQ/2.D0) 218 : ELSE 219 : DI1(I,I)=0.D0 220 : WRITE(6,*) 'erreur' 221 : ENDIF 222 : ENDDO 223 : C --- passage dans le repère de général --- 223 : CALL MULMAT(AAA,ROT,DI1,3,3,3) 223 : CALL MULMAT(DI2,AAA,ROTT,3,3,3) 224 : C 224 : CALL MULMAT(AAA,DI2,CO,3,3,3) 224 : CALL MULMAT(COM,AAA,DI2,3,3,3) 225 : ENDIF 226 : C 227 : C---------------------------------------------------------------------- 228 : C 229 : C CALCUL DE J1,J2 230 : C J1: CONTRAINTE MOYENNE 231 : C J2: CONTRAINTE DE VON MISES 232 : C 233 : C---------------------------------------------------------------------- 234 : C 235 : AJ1 = (SIG(1)+SIG(2)+SIG(3))/(3.D0) 236 : DO I=1,3 237 : XX(I) = SIG(I)-AJ1 238 : ENDDO 239 : DO I=4,NSTRS 240 : XX(I) = SIG(I) 241 : ENDDO 242 : AJ2 = PROCON(XX,XX,NSTRS) 243 : AJ2 = SQRT(1.5D0*AJ2) 244 : C 245 : C---------------------- 246 : C FORCE THERMODYNAMIQUE 247 : C---------------------- 248 : C 249 : C --- Etat de triaxialité de la contrainte --- 250 : tri = (AJ2)/(SIG(1)+SIG(2)+SIG(3)) 251 : C --- Valeurs propres et vecteurs propres du tenseur des contraintes motrices --- 251 : CALL JACOB3(COM,3,VaP,VcP) 252 : ER=XMAT(8) 253 : EA=XMAT(9) 254 : DO I=1,3 255 : AA(I) = ((tri*VaP(I)/(EA))+ABS(tri*VaP(I)/(EA)))/(2.D0) 256 : AA(I) = (AA(I))**(ER) 257 : ENDDO 258 : C -- Dans le repère général --- 258 : CALL ZERO(FTH,3,3) 259 : FTH(1,1) = AA(1)*(VcP(1,1))**2+AA(2)*(VcP(1,2))**2 260 : FTH(1,1) = FTH(1,1)+AA(3)*(VcP(1,3))**2 261 : FTH(2,2) = AA(1)*(VcP(2,1))**2+AA(2)*(VcP(2,2))**2 262 : FTH(2,2) = FTH(2,2)+AA(3)*(VcP(2,3))**2 263 : FTH(3,3) = AA(1)*(VcP(3,1))**2+AA(2)*(VcP(3,2))**2 264 : FTH(3,3) = FTH(3,3)+AA(3)*(VcP(3,3))**2 265 : FTH(1,2) = AA(1)*VcP(1,1)*VcP(2,1)+AA(2)*VcP(1,2)*VcP(2,2) 266 : FTH(1,2) = FTH(1,2)+AA(3)*VcP(1,3)*VcP(2,3) 267 : FTH(1,3) = AA(1)*VcP(1,1)*VcP(3,1)+AA(2)*VcP(1,2)*VcP(3,2) 268 : FTH(1,3) = FTH(1,3)+AA(3)*VcP(1,3)*VcP(3,3) 269 : FTH(2,3) = AA(1)*VcP(2,1)*VcP(3,1)+AA(2)*VcP(2,2)*VcP(3,2) 270 : FTH(2,3) = FTH(1,3)+AA(3)*VcP(2,3)*VcP(3,3) 271 : FTH(3,1) = FTH(1,3) 272 : FTH(2,1) = FTH(1,2) 273 : FTH(3,2) = FTH(2,3) 274 : C 275 : BET2=XMAT(7) 276 : DV1(1,1) = (BET2*FTH(1,1))+FTH(2,2)+FTH(3,3) 277 : DV1(2,2) = (BET2*FTH(2,2))+FTH(1,1)+FTH(3,3) 278 : DV1(3,3) = (BET2*FTH(3,3))+FTH(2,2)+FTH(1,1) 279 : DV1(1,2) = (BET2-1.D0)*FTH(1,2) 280 : DV1(1,3) = (BET2-1.D0)*FTH(1,3) 281 : DV1(2,3) = (BET2-1.D0)*FTH(2,3) 282 : DV1(2,1) = DV1(1,2) 283 : DV1(3,1) = DV1(1,3) 284 : DV1(3,2) = DV1(2,3) 285 : C 286 : C---------------------------------------------------------------------- 287 : C 288 : C CALCUL DE LA VARIATION DES CONTRAINTES 289 : C 290 : C---------------------------------------------------------------------- 291 : C 291 : CALL DERTRA(NYOUNG,YOUNG,TI,YUNG,YUNGV,TINF,TSUP) 291 : CALL DERTRA(NXNU,XNU,TI,ENU,ENUV,TINF,TSUP) 292 : C 293 : C --- Calcul de l inverse de la matrice de Hooke au debut du pas --- 293 : CALL ELAST4(2,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV, 294 : &XCAR,ICARA,MFR,NSTRS,DDINV,DDV,KERRE,INDIC,ITHHER) 295 : IF(KERRE.NE.0) GO TO 100 296 : C 297 : C --- Calcul de la matrice de Hooke et de sa derivee au debut du pas --- 297 : CALL ELAST4(1,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV, 298 : &XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER) 299 : IF(KERRE.NE.0) GO TO 100 300 : C 301 : C --- Détermination de la variation de la matrive de Hooke inversée --- 301 : CALL HOOKP(VAR,VARD,NVARI,XMAT,NCOMAT,MFR,NSTRS,DDINVp) 302 : C 303 : C --- Calcul de la variation de contrainte --- 303 : CALL ZDANUL(SIGD,NSTRS) 304 : DO I=1,3 305 : XX(I) = CO(I,I) 306 : ENDDO 307 : XX(4) = CO(1,2) 308 : XX(5) = CO(1,3) 309 : XX(6) = CO(2,3) 310 : C 310 : CALL ZDANUL(EPSED,NSTRS) 311 : DO I=1,NSTRS 312 : EPSED(I)=DSPT(I)-EPSVD(I)-EPSTHD(I) 313 : ENDDO 314 : C 314 : CALL MULMAT(BB,DDINVp,XX,NSTRS,1,NSTRS) 315 : C 316 : DO I=1,NSTRS 317 : XX(I)=EPSED(I)-BB(I) 318 : ENDDO 319 : C 319 : CALL MULMAT(SIGD,DD,XX,NSTRS,1,NSTRS) 320 : C 321 : 100 CONTINUE 322 : RETURN 323 : END 324 : 325 : 326 : 327 : 328 : 329 : 330 : 331 :
© Cast3M 2003 - All rights reserved.
Disclaimer