C XXCRE1 SOURCE CB215821 17/11/30 21:17:47 9639 SUBROUTINE XXCRE1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD, &DSPT,XMAT,XMAINF,XMASUP,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV, &DDINVp,YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI, &TPOINT,TINF,TSUP,ITEST,ITHHER,VARF,IB,IGAU,KERRE) C----------------------------------------------------------------------- C Objet: Cette subroutine calcule les derivees des variables internes C pour un materiau viscoplastique a endommagement et ecrouissage C isotrope en regime anisotherme de 2 manieres au choix suivant C la valeur de ITEST,ITHHER C----------------------------------------------------------------------- C C----------------------------------------------------------------------- C Entree: TAU pas de temps C TI temperature au debut du pas C TPOINT derivee de la temperature C TINF et TSUP bornes entre lesquelles TI est comprise C SIG(NSTRS) tenseur des contraintes C EPSV(NSTRS) tenseur des deformations C VAR(NVARI) tableau contenant les variables internes C pilotant les equations C VAR(1)=p ; VAR(2)=D11 ; VAR(3)=D22; VAR(4)=D33; C VAR(5)=D12; VAR(6)=D13; VAR(7)=D23; VAR(8)=DMAX C EPSTHD vitesse de dilatation thermique au debut du pas C DSPT(NSTRS,NSTRS) vitesse de deformation totale C XMAT(NCOMAT) tableau des parametres scalaires du materiau C a une temperature T donnee C XMAINF(NCOMAT) tableau des parametres scalaires du materiau C a la temperature TINF C XMASUP(NCOMAT) tableau des parametres scalaires du materiau C a la temperature TSUP C XMAT(1)=YOUNG ;XMAT(2)=XNU ;XMAT(3)=N C XMAT(4)=M ;XMAT(5)=KK; XMAT(6)=ALF2; C XMAT(7)=BET2 ;XMAT(8)=r; XMAT(9)=A ; C XMAT(10)=q ;XMAT(11)=ALPHAT C XMAT(12)=RHO; XMAT(13)=SIGY C MFR indice de la formulation mecanique(seulement massif C ou coque pour les materiaux endommageables C ICARA nombre de caracteristiques geometriques des elements C finis C XCAR(ICARA) tableau des caracteristiques geometriques des C elements finis C IFOURB= -2 EN CONTR. PLANES C -1 EN DEFORM. PLANES C 0 EN AXISYMETRIE C 1 EN SERIE DE FOURIER C 2 EN TRIDIM C INDIC=0, 1 OU -1 pour plasticite avec endommagement C =2 OU -2 pour viscoplasticite avec endommagement C ITEST = 0 pas uniformite des listes de temperatures pour tous les C coefficients non lineaires C d'ou interpolation sur les coefficients du materiau C = 1 uniformite des listes de temperatures pour tous les C coefficients non lineaires C d'ou moyenne ponderee sur la variable elle-meme C ex: D(T)=A(T)*K**B(T) C T=TETA*TINF+(1-TETA)*TSUP C si ITEST=0 on remplace A(T) par (TETA*A(TINF)+(1-TETA))*A(TSUP) C et B(T) par (TETA*B(TINF)+(1-TETA))*B(TSUP) C si ITEST=1 D(T)=TETA*D(TINF)+(1-TETA)*D(TSUP) C ITHHER = 0 pas de chargement thermique et materiau constant C = 1 chargement thermique et materiau constant C = 2 chargement thermique et materiau(T) C------------------------------------------------------------------------ C C------------------------------------------------------------------------ C Sortie: EPSVD(NSTRS) derivee du tenseur des deformations C VARD(NVARI) tableau contenant les derivees des variables C internes C VARF(NVARI) état des variables internes à la fin du sous programme C SIGD(NSTRS) derivee des contraintes C DD(NSTRS,NSTRS) matrice de Hooke au debut du pas C DDV(NSTRS,NSTRS) derivee de DD C DDINV(NSTRS,NSTRS) inverse de DD C------------------------------------------------------------------------ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XMAINF(*),XMASUP(*),XCAR(*) DIMENSION SIG(*),EPSV(*),VAR(*) DIMENSION SIGD(*),EPSVD(*),VARD(*),VARF(NVARI) DIMENSION XMAT(*),DSPT(*),YOUNG(*),XNU(*),EPSTHD(*) DIMENSION DDV(NSTRS,NSTRS),DDINVp(NSTRS,NSTRS) DIMENSION DD(NSTRS,NSTRS),DDINV(NSTRS,NSTRS) DIMENSION VARD1(3),EPSVD1(6),SIGD1(6),VARD2(3),EPSVD2(6) DIMENSION EPSED(6),XX(6),Di(3),VaP(3),AA(3),BB(6) DIMENSION ROT(3,3),ROTT(3,3),CO(3,3),COM(3,3) DIMENSION D1(3,3),D2(3,3),DV1(3,3),DI1(3,3),DI2(3,3) DIMENSION AAA(3,3),BBB(3,3),CCC(3,3),VcP(3,3),FTH(3,3) PARAMETER (UN=1.D0) LOGICAL endonul C C----- le materiau est constant ou pas, on derive normalement les contraintes C IFLAG=0 IF ( (ITEST.EQ.1).AND.(ITHHER.EQ.2) ) THEN C----- si le materiau depend de la temperature et si les listes sont identiques IF (TINF.NE.TSUP) THEN C-------------- Moyenne ponderee sur les vitesses des variables IFLAG=1 TETA=(TI-TSUP)/(TINF-TSUP) CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD1,EPSVD1,VARD1,EPSTHD, & DSPT,XMAINF,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp, & YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC, & TINF,TPOINT,ITHHER,VARF,IB,IGAU,KERRRE) IF(KERRE.NE.0) GO TO 100 DO 10 I=1,NVARI 10 VARD2(I)=TETA*VARD1(I) DO 20 I=1,NSTRS 20 EPSVD2(I)=TETA*EPSVD1(I) CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD1,EPSVD1,VARD1,EPSTHD, & DSPT,XMASUP,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp, & YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC, & TSUP,TPOINT,ITHHER,VARF,IB,IGAU,KERRRE) IF(KERRE.NE.0) GO TO 100 DO 30 I=1,NVARI 30 VARD2(I)=VARD2(I)+(1.D0-TETA)*VARD1(I) DO 40 I=1,NSTRS 40 EPSVD2(I)=EPSVD2(I)+(1.D0-TETA)*EPSVD1(I) ENDIF ENDIF C CALL DERIV1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD, &DSPT,XMAT,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,YOUNG, &NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI,TPOINT, &ITHHER,VARF,IB,IGAU,KERRE) IF(KERRE.NE.0) GO TO 100 C IF (IFLAG.EQ.0) GOTO 100 C C --- Cas ou ITEST = 1 et ITHHER = 2 C DO I=1,NVARI VARD(I)=VARD2(I) ENDDO DO I=1,NSTRS EPSVD(I)=EPSVD2(I) ENDDO C endonul = .FALSE. IF ((VAR(2).EQ.0).AND.(VAR(3).EQ.0).AND.(VAR(4).EQ.0).AND. & (VAR(5).EQ.0).AND.(VAR(6).EQ.0).AND.(VAR(7).EQ.0).AND. & (VAR(8).EQ.0)) THEN endonul = .TRUE. ENDIF C C---------------------------------------------- C TENSEUR DES CONTRAINTES (dans repere general) C---------------------------------------------- C CALL ZERO(CO,3,3) IF (mfr.eq.1) THEN DO I=1,3 CO(I,I) = SIG(I) ENDDO CO(1,2) = SIG(4) IF(NSTRS.NE.6) GOTO 50 CO(1,3) = SIG(5) CO(2,3) = SIG(6) 50 CONTINUE CO(2,1) = CO(1,2) CO(3,1) = CO(1,3) CO(3,2) = CO(2,3) ENDIF C C------------------------------------ C TENSEUR D ENDOMMAGEMENT ET C TENSEUR D ENDOMMAGEMENT DIAGONALISE C------------------------------------ C XD11=VAR(2) XD22=VAR(3) XD33=VAR(4) XD12=VAR(5) XD13=VAR(6) XD23=VAR(7) CALL ZERO(D1,3,3) D1(1,1) = XD11 D1(2,2) = XD22 D1(3,3) = XD33 D1(1,2) = XD12 D1(1,3) = XD13 D1(2,3) = XD23 D1(2,1) = D1(1,2) D1(3,1) = D1(1,3) D1(3,2) = D1(2,3) C CALL JACOB3(D1,3,Di,ROT) C CALL ZERO(D2,3,3) D2(1,1) = Di(1) D2(2,2) = Di(2) D2(3,3) = Di(3) w1 = 1.D0 - Di(1) w2 = 1.D0 - Di(2) w3 = 1.D0 - Di(3) IF (D2(1,1).EQ.1.D0) w1=1.D-6 C ADMAX = MAX(Di(1),Di(2),Di(3)) C C------------------------------------------------------- C TENSEUR DES CONTRAINTES MOTRICES (dans repere general) C------------------------------------------------------- C EQ=XMAT(10) DO I=1,3 DO J=1,3 ROTT(I,J)=ROT(I,J) ENDDO ENDDO CALL INVALM (ROTT, 3 , 3 , KER, 1.D-12) C IF (endonul) THEN DO I=1,3 DO J=1,3 COM(I,J) = CO(I,J) ENDDO ENDDO ELSE C --- Matrice de [(I-D)/det(I-D)]**(q/2) dans le repère principal de D --- CALL ZERO(DI1,3,3) DO I=1,3 XNB=(1.D0-Di(I))/(w1*w2*w3) IF (XNB.GE.0.D0) THEN DI1(I,I)=(XNB)**(EQ/2.D0) ELSE DI1(I,I)=0.D0 WRITE(6,*) 'erreur' ENDIF ENDDO C --- passage dans le repère de général --- CALL MULMAT(AAA,ROT,DI1,3,3,3) CALL MULMAT(DI2,AAA,ROTT,3,3,3) C CALL MULMAT(AAA,DI2,CO,3,3,3) CALL MULMAT(COM,AAA,DI2,3,3,3) ENDIF C C---------------------------------------------------------------------- C C CALCUL DE J1,J2 C J1: CONTRAINTE MOYENNE C J2: CONTRAINTE DE VON MISES C C---------------------------------------------------------------------- C AJ1 = (SIG(1)+SIG(2)+SIG(3))/(3.D0) DO I=1,3 XX(I) = SIG(I)-AJ1 ENDDO DO I=4,NSTRS XX(I) = SIG(I) ENDDO AJ2 = PROCON(XX,XX,NSTRS) AJ2 = SQRT(1.5D0*AJ2) C C---------------------- C FORCE THERMODYNAMIQUE C---------------------- C C --- Etat de triaxialité de la contrainte --- tri = (AJ2)/(SIG(1)+SIG(2)+SIG(3)) C --- Valeurs propres et vecteurs propres du tenseur des contraintes motrices --- CALL JACOB3(COM,3,VaP,VcP) ER=XMAT(8) EA=XMAT(9) DO I=1,3 AA(I) = ((tri*VaP(I)/(EA))+ABS(tri*VaP(I)/(EA)))/(2.D0) AA(I) = (AA(I))**(ER) ENDDO C -- Dans le repère général --- CALL ZERO(FTH,3,3) FTH(1,1) = AA(1)*(VcP(1,1))**2+AA(2)*(VcP(1,2))**2 FTH(1,1) = FTH(1,1)+AA(3)*(VcP(1,3))**2 FTH(2,2) = AA(1)*(VcP(2,1))**2+AA(2)*(VcP(2,2))**2 FTH(2,2) = FTH(2,2)+AA(3)*(VcP(2,3))**2 FTH(3,3) = AA(1)*(VcP(3,1))**2+AA(2)*(VcP(3,2))**2 FTH(3,3) = FTH(3,3)+AA(3)*(VcP(3,3))**2 FTH(1,2) = AA(1)*VcP(1,1)*VcP(2,1)+AA(2)*VcP(1,2)*VcP(2,2) FTH(1,2) = FTH(1,2)+AA(3)*VcP(1,3)*VcP(2,3) FTH(1,3) = AA(1)*VcP(1,1)*VcP(3,1)+AA(2)*VcP(1,2)*VcP(3,2) FTH(1,3) = FTH(1,3)+AA(3)*VcP(1,3)*VcP(3,3) FTH(2,3) = AA(1)*VcP(2,1)*VcP(3,1)+AA(2)*VcP(2,2)*VcP(3,2) FTH(2,3) = FTH(1,3)+AA(3)*VcP(2,3)*VcP(3,3) FTH(3,1) = FTH(1,3) FTH(2,1) = FTH(1,2) FTH(3,2) = FTH(2,3) C BET2=XMAT(7) DV1(1,1) = (BET2*FTH(1,1))+FTH(2,2)+FTH(3,3) DV1(2,2) = (BET2*FTH(2,2))+FTH(1,1)+FTH(3,3) DV1(3,3) = (BET2*FTH(3,3))+FTH(2,2)+FTH(1,1) DV1(1,2) = (BET2-1.D0)*FTH(1,2) DV1(1,3) = (BET2-1.D0)*FTH(1,3) DV1(2,3) = (BET2-1.D0)*FTH(2,3) DV1(2,1) = DV1(1,2) DV1(3,1) = DV1(1,3) DV1(3,2) = DV1(2,3) C C---------------------------------------------------------------------- C C CALCUL DE LA VARIATION DES CONTRAINTES C C---------------------------------------------------------------------- C CALL DERTRA(NYOUNG,YOUNG,TI,YUNG,YUNGV,TINF,TSUP) CALL DERTRA(NXNU,XNU,TI,ENU,ENUV,TINF,TSUP) C C --- Calcul de l inverse de la matrice de Hooke au debut du pas --- CALL ELAST4(2,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV, &XCAR,ICARA,MFR,NSTRS,DDINV,DDV,KERRE,INDIC,ITHHER) IF(KERRE.NE.0) GO TO 100 C C --- Calcul de la matrice de Hooke et de sa derivee au debut du pas --- CALL ELAST4(1,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV, &XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER) IF(KERRE.NE.0) GO TO 100 C C --- Détermination de la variation de la matrive de Hooke inversée --- CALL HOOKP(VAR,VARD,NVARI,XMAT,NCOMAT,MFR,NSTRS,DDINVp) C C --- Calcul de la variation de contrainte --- CALL ZDANUL(SIGD,NSTRS) DO I=1,3 XX(I) = CO(I,I) ENDDO XX(4) = CO(1,2) XX(5) = CO(1,3) XX(6) = CO(2,3) C CALL ZDANUL(EPSED,NSTRS) DO I=1,NSTRS EPSED(I)=DSPT(I)-EPSVD(I)-EPSTHD(I) ENDDO C CALL MULMAT(BB,DDINVp,XX,NSTRS,1,NSTRS) C DO I=1,NSTRS XX(I)=EPSED(I)-BB(I) ENDDO C CALL MULMAT(SIGD,DD,XX,NSTRS,1,NSTRS) C 100 CONTINUE RETURN END