C ELAST1 SOURCE CHAT 05/01/12 23:29:27 5004 SUBROUTINE ELAST1(ISENS,IFOUR,VAR,NNVARI,XMATT,NCOMAT, & YUNGV,XNUV,XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER) C C---------------------------------------------------------------------------- C Objet: Calcul de la matrice de Hooke endommagee(ISENS=1) C ou son inverse(ISENS=2) et de sa derivee par rapport C a T a la temperature T0 dans le cas de la plasticite avec C ecrouissage et endommagement isotropes(materiau 26) C et dans le cas de la viscoplasticite avec endommagement C isotrope(materiau 29) en regime anisotherme C---------------------------------------------------------------------------- C C---------------------------------------------------------------------------- C Entree: ISENS= 1 DD en sortie est la matrice de Hooke endommagee C d'ou SIG=DD*EPSELAS C = 2 DD en sortie est la matrice de Hooke inverse C endommagee d'ou EPSELAS=DD*SIG C avec SIG contraintes et EPSELAS deformations elastiques C IFOUR= -2 EN CONTR. PLANES C -1 EN DEFORM. PLANES C 0 EN AXISYMETRIE C 1 EN SERIE DE FOURIER C 2 EN TRIDIM C NSTRS nombre de composantes des contraintes ou C des deformations C NNVARI nombre de variables internes pilotant les C equations du modele C VAR(NVARI) tableau des variables internes C XMATT(NCOMAT) tableau des parametres scalaires du materiau C a une temperature T donnee C dont XMATT(1) module d'Young et XMATT(2) coeff de Poisson C YUNGV derivee /T du module d'Young a T C XNUV derivee /T du coeff de Poisson a T C MFR indice de la formulation mecanique(seulement massif C 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 INDIC=0, 1 OU -1 pour plasticite avec endommagement C =2 OU -2 pour viscoplasticite avec endommagement 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: DD(NSTRS,NSTRS) matrice de Hooke endommagee si ISENS=1 C son inverse si ISENS=2 C DDV(NSTRS,NSTRS) derivee de la matrice de Hooke C endommagee si ISENS=1 C son inverse si ISENS=2 C KERRE indice qui regit les erreurs C Il vaut 99 si la formulation mecanique n'est pas disponible ou C s'il y a incompatibilite entre MFR et IFOUR C _______________________________________________________________________________ C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION VAR(*),XCAR(*),XMATT(*) DIMENSION DD(NSTRS,NSTRS),DDV(NSTRS,NSTRS) C CALL ZERO(DD,NSTRS,NSTRS) CALL ZERO(DDV,NSTRS,NSTRS) KERRE=0 C YUNG=XMATT(1) XNU=XMATT(2) XD=VAR(3) * write(6,*) 'elast1 ifour mfr ',ifour,mfr C D1=1.D0-XD C C ==================================================================== C SI LE PT. DE GAUSS EST ROMPU, ON DIV. SON MOD. D'YOUNG PAR 1 MILLION C ==================================================================== IF (INDIC.EQ.2.OR.INDIC.EQ.-2) THEN IF (XD.EQ.1.D0) D1=1.D-6 GOTO 3 ENDIF C DC=XMATT(7) IF (INDIC.EQ.-1.AND.XD.GE.DC) D1=1.D-6 C 3 CONTINUE COEFF1=(D1*YUNG)/(1.D0-XNU*XNU) COEFF2=1.D0/(D1*YUNG) COEFF3=(D1*YUNG)/((1.D0+XNU)*(1.D0-2.D0*XNU)) COEFF4=(1.D0-2.D0*XNU)/2.D0 C IF (ITHHER.EQ.2) THEN COEF1V=YUNGV+(2.D0*XNUV*YUNG*XNU/(1.D0-XNU*XNU)) COEF1V=COEF1V*D1/(1.D0-XNU*XNU) C COEF2V=-YUNGV/(YUNG*YUNG*D1) C DEN=(1.D0+XNU)*(1.D0-2.D0*XNU) COEF3V=(1.D0+(4.D0*XNU))*YUNG*XNUV COEF3V=COEF3V/DEN COEF3V=COEF3V+YUNGV COEF3V=D1*COEF3V/DEN C COEF4V=-XNUV ENDIF C IF (MFR.NE.1) GOTO 20 C C------------ C MASSIFS C------------ IF (IFOUR.EQ.-2) GOTO 10 C C ---------------------------------- C DEFORMATIONS PLANES ET AXISYMETRIE C ---------------------------------- IF (ISENS.EQ.1) THEN C C---------------- C SIG=DD*EPSELAS C---------------- DD(1,1)=COEFF3*(1.D0-XNU) DD(1,2)=COEFF3*XNU DD(1,3)=DD(1,2) DD(2,2)=DD(1,1) DD(2,3)=DD(1,2) DD(3,3)=DD(1,1) DD(2,1)=DD(1,2) DD(3,1)=DD(1,2) DD(3,2)=DD(1,2) C IF (ITHHER.EQ.2) THEN DDV(1,1)=(COEF3V*(1.D0-XNU))-(COEFF3*XNUV) DDV(1,2)=(COEF3V*XNU)+(COEFF3*XNUV) DDV(1,3)=DDV(1,2) DDV(2,2)=DDV(1,1) DDV(2,3)=DDV(1,2) DDV(3,3)=DDV(1,1) DDV(2,1)=DDV(1,2) DDV(3,1)=DDV(1,2) DDV(3,2)=DDV(1,2) ENDIF LA=4 LB=4 C C --------------------------- C LB=6 POUR TRIDIM ET FOURIER C --------------------------- IF (IFOUR.GT.0) LB=6 DO 5 IA=LA,LB DD(IA,IA)=COEFF3*COEFF4 IF (ITHHER.EQ.2) THEN DDV(IA,IA)=(COEF3V*COEFF4)+(COEFF3*COEF4V) ENDIF 5 CONTINUE ELSE C C------------------- C EPSELAS=DD*SIG C------------------- DD(1,1)=COEFF2 DD(1,2)=-XNU*COEFF2 DD(1,3)=DD(1,2) DD(2,2)=DD(1,1) DD(2,3)=DD(1,2) DD(3,3)=DD(1,1) DD(2,1)=DD(1,2) DD(3,1)=DD(1,3) DD(3,2)=DD(2,3) IF (ITHHER.EQ.2) THEN DDV(1,1)=COEF2V DDV(1,2)=-( (XNUV*COEFF2)+(XNU*COEF2V) ) DDV(1,3)=DDV(1,2) DDV(2,2)=DDV(1,1) DDV(2,3)=DDV(1,2) DDV(3,3)=DDV(1,1) DDV(2,1)=DDV(1,2) DDV(3,1)=DDV(1,3) DDV(3,2)=DDV(2,3) ENDIF LA=4 LB=4 C C --------------------------- C LB=6 POUR TRIDIM ET FOURIER C --------------------------- IF (IFOUR.GT.0) LB=6 DO 7 IA=LA,LB DD(IA,IA)=(1.D0+XNU)*DD(1,1)*2.D0 IF (ITHHER.EQ.2) THEN DDV(IA,IA)=2.D0*( (XNUV*COEFF2)+((1.D0+XNU)*COEF2V) ) ENDIF 7 CONTINUE C ENDIF C GOTO 500 C C ------------- C CONTR. PLANES C ------------- 10 IF (ISENS.EQ.1) THEN C C----------------- C SIG=DD*EPSELAS C----------------- DD(1,1)=COEFF1 DD(1,2)=COEFF1*XNU DD(2,2)=COEFF1 DD(4,4)=COEFF1*((1.D0-XNU)/2.D0) DD(2,1)=DD(1,2) IF (ITHHER.EQ.2) THEN DDV(1,1)=COEF1V DDV(1,2)=(COEF1V*XNU)+(COEFF1*XNUV) DDV(2,2)=COEF1V DDV(4,4)=0.5D0*( (COEF1V*(1.D0-XNU))-(COEFF1*XNUV) ) DDV(2,1)=DDV(1,2) ENDIF ELSE C C------------------ C EPSELAS=DD*SIG C------------------ DD(1,1)=COEFF2 DD(1,2)=-COEFF2*XNU DD(2,2)=COEFF2 DD(4,4)=COEFF2*2.D0*(1.D0+XNU) DD(2,1)=DD(1,2) DD(3,1)=DD(1,2) DD(3,2)=DD(1,2) C IF (ITHHER.EQ.2) THEN DDV(1,1)=COEF2V DDV(1,2)=-( (COEF2V*XNU)+(COEFF2*XNUV) ) DDV(2,2)=COEF2V DDV(4,4)=2.D0*( ((1.D0+XNU)*COEF2V)+(XNUV*COEFF2) ) DDV(2,1)=DDV(1,2) DDV(3,1)=DDV(1,2) DDV(3,2)=DDV(1,2) ENDIF ENDIF C GOTO 500 C 20 IF (MFR.NE.5) GOTO 30 C C------------------------------------------ C COQUES EPAISSES TRIDIM C------------------------------------------ C C IF(IFOUR.EQ.2) THEN XK=1.2d0 IF (ISENS.EQ.1) THEN C sig= DD * EPSELAS AUX=COEFF1 AUX1=COEFF1*XNU GEGE=YUNG*D1*0.5D0/(1.D0+XNU) C DD(1,1)=AUX DD(2,1)=AUX1 DD(1,2)=AUX1 DD(2,2)=AUX DD(3,3)=GEGE DD(4,4)=GEGE/XK DD(5,5)=DD(4,4) IF (ITHHER.EQ.2) THEN DDV(1,1)=COEF1V DDV(1,2)=(COEF1V*XNU)+(COEFF1*XNUV) DDV(2,2)=COEF1V DDV(2,1)=DDV(1,2) DDV(3,3)=0.5D0*( (COEF1V*(1.D0-XNU))-(COEFF1*XNUV) ) DDV(4,4)=DDV(3,3) * XK DDV(5,5)=DDV(3,3) * XK ENDIF ELSE C eps = DD* SIG C------------------ DD(1,1)=COEFF2 DD(1,2)=-COEFF2*XNU DD(2,2)=COEFF2 DD(3,3)=COEFF2*2.D0*(1.D0+XNU) DD(2,1)=DD(1,2) DD(4,4)=DD(3,3)/XK DD(5,5)=DD(3,3)/XK C IF (ITHHER.EQ.2) THEN DDV(1,1)=COEF2V DDV(1,2)=-( (COEF2V*XNU)+(COEFF2*XNUV) ) DDV(2,2)=COEF2V DDV(3,3)=2.D0*( ((1.D0+XNU)*COEF2V)+(XNUV*COEFF2) ) DDV(2,1)=DDV(1,2) DDV(4,4)=DDV(3,3)/XK DDV(5,5)=DDV(3,3)/XK ENDIF ENDIF ELSE IRET=0 ENDIF GO TO 500 30 IF(MFR.NE.3) go to 40 IF (ISENS.EQ.1) THEN C C----------------- C SIG=DD*EPSELAS C----------------- DD(1,1)=COEFF1 DD(1,2)=COEFF1*XNU DD(2,2)=COEFF1 DD(4,4)=COEFF1*((1.D0-XNU)/2.D0) C voir petite ruse de hookis pour (3,3) DD(3,3)=COEFF1*1.d-6 DD(2,1)=DD(1,2) IF (ITHHER.EQ.2) THEN DDV(1,1)=COEF1V DDV(1,2)=(COEF1V*XNU)+(COEFF1*XNUV) DDV(2,2)=COEF1V DDV(4,4)=0.5D0*( (COEF1V*(1.D0-XNU))-(COEFF1*XNUV) ) DDV(2,1)=DDV(1,2) ENDIF ELSE C C------------------ C EPSELAS=DD*SIG C------------------ DD(1,1)=COEFF2 DD(1,2)=-COEFF2*XNU DD(2,2)=COEFF2 DD(4,4)=COEFF2*2.D0*(1.D0+XNU) DD(2,1)=DD(1,2) DD(3,1)=DD(1,2) DD(3,2)=DD(1,2) C IF (ITHHER.EQ.2) THEN DDV(1,1)=COEF2V DDV(1,2)=-( (COEF2V*XNU)+(COEFF2*XNUV) ) DDV(2,2)=COEF2V DDV(4,4)=2.D0*( ((1.D0+XNU)*COEF2V)+(XNUV*COEFF2) ) DDV(2,1)=DDV(1,2) DDV(3,1)=DDV(1,2) DDV(3,2)=DDV(1,2) ENDIF ENDIF GO TO 500 C FORMULATION NON DISPONIBLE 40 KERRE=99 500 CONTINUE RETURN END