C MAZARS SOURCE CHAT 05/01/13 01:38:20 5004 SUBROUTINE MAZARS (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI, 1 NMATT,ISTEP,ICARA,JDIM,IFOUR2) C C C variables en entree C C C WRK0 pointeur sur un segment deformation au pas precedent C C WRK1 pointeur sur un segment increment de deformation C C WRKK2 pointeur sur un segment variables internes au pas precedent C C WRK5 pointeur sur un segment de deformations inelastiques C C XMATER constantes du materiau C C NSTRS nombre de composantes dans les vecteurs des contraintes C et les vecteurs des deformations C C NVARI nombre de variables internes (doit etre egal a 2) C C NMATT nombre de constantes du materiau C C ISTEP flag utilise pour separer les etapes dans un calcul non local C ISTEP=0 -----> calcul local C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils C ISTEP=2 -----> calcul non local etape 2 on continue le calcul C a partir des seuils moyennes C C Modif L.Bode - 14/10/92 C Nouveaux parametres en entree C JDIM Dimension de travail C ( Coques JDIM =2 , Massifs JDIM = IDIM ) C IFOUR2 Type de formulation C ( Coques IFOUR2 = -2 => contraintes planes , C Massifs IFOUR2 = IFOUR) C C variables en sortie C C VARINF variables internes finales C C SIGMAF contraintes finales C C C. LA BORDERIE MARS 1992 C declaration des variables C C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO SEGMENT WRK0 REAL*8 XMAT(NMATT) ENDSEGMENT * SEGMENT WRK1 REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) REAL*8 DEFP(NSTRS),XCAR(ICARA) ENDSEGMENT SEGMENT WRKK2 REAL*8 EPSILI(NSTRS),DSIGT(NSTRS) ENDSEGMENT C Modif L.Bode - 14/10/92 SEGMENT WRK3 REAL*8 EPSILO(NSTRS) ENDSEGMENT C Definition dynamique de EPSILO C Fin modif L.Bode SEGMENT WRK5 REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS) ENDSEGMENT CHARACTER*8 CMATE INTEGER NSTRS,NVARI,NMATT,ISTEP REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3) REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA INTEGER ISTRS,JSTRS,KCAS,IRTD REAL*8 XZERO,UN,DEUX,XPETIT REAL*8 DINI,D,DT,DC,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12) C C C recuperation des variables initiales dans les tableaux C C N=NSTRS CMATE = 'ISOTROPE' YOUN = XMAT(1) XNU = XMAT(2) EPSD0= XMAT(5) ACOM = XMAT(6) BCOM = XMAT(7) ATRA = XMAT(8) BTRA = XMAT(9) BETA = XMAT(10) DINI = VAR0(2) C C calcul des seuils C C C calcul de la deformation totale C SEGINI WRK3 DO 100 ISTRS=1,NSTRS EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS) 100 CONTINUE C C calcul des deformations principales C C C on reecrit les deformations sous forme matricielle C C Modif L.Bode - 14/10/92 C Rajout de IFOUR2 en argument de ENDOCA * print*,'on appelle ENDOCB' CALL ENDOCB (EPSILO,EPS33,2,IFOUR2) * print*,'apres endocb' C Fin modif L.Bode C C et on diagonalise C C Modif L.Bode - 14/10/92 C Pour les elts Coques, on travaille en contraintes planes => JDIM =2 C Pour les elts Massifs JDIM =IDIM * print*,'avant JACOB3' CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33) * print*,'apres JACOB3' C Fin modif L.Bode IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN C C on calcule la matrice de hooke et les contraintes ppales C CMATE = 'ISOTROPE' KCAS=1 C Modif L.Bode - 14/10/92 C IFOUR --> IFOUR2 dans appel DOHMAS * print*,'avant dohmas' CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS,KCAS,DDHOOK,IRTD) * print*,'apres dohmas' C Fin modif L.Bode DO 200 ISTRS=1,3 SIGP(ISTRS)= XZERO DO 210 JSTRS=1,3 SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS) 210 CONTINUE 200 CONTINUE END IF C C on complete la deformation dans le cas des contraintes planes C C Modif L.Bode - 14/10/92 C IFOUR remplace par IFOUR2 IF (IFOUR2.EQ. -2) THEN EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU) END IF C Fin modif L.Bode C C on calcule le epsilontild C EPSTIL=MAX( XZERO , EPSIPP(1) )**2 + 1 MAX( XZERO , EPSIPP(2) )**2 + 2 MAX( XZERO , EPSIPP(3) )**2 EPSTIL=SQRT (EPSTIL) IF (ISTEP .EQ. 0) THEN EPSTIM=EPSTIL VARF(1)=EPSTIL ELSE IF (ISTEP .EQ. 1) THEN VARF(2)=DINI VARF(1)=EPSTIL ELSE IF (ISTEP .EQ. 2) THEN EPSTIM=VAR0(1) VARF(1)=EPSTIM ELSE PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP END IF IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN C C on calcule l'endommagement et les contraintes C C C on calcule ALFAT ALFAC DT et DC puis D C dans le cas ou le seuil initial est depasse C IF ( EPSTIM .GT. EPSD0) THEN C C calcul de l'endommagement C C C on calcule le signe des contraintes elastiques C DO 300 ISTRS=1,3 IF (SIGP(ISTRS).LT. XZERO) THEN SIGPC(ISTRS)=SIGP(ISTRS) SIGPT(ISTRS)=XZERO ELSE SIGPT(ISTRS)=SIGP(ISTRS) SIGPC(ISTRS)=XZERO END IF 300 CONTINUE TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3) TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3) C C on calcule les deformations dues aux contraintes positives C DO 400 ISTRS=1,3 EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN 400 CONTINUE C C on en deduit ALFAT et ALFAC C ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) + 1 MAX(XZERO,EPSIPP(2))*EPSILT(2) + 2 MAX(XZERO,EPSIPP(3))*EPSILT(3) ALFAT = ALFAT/(EPSTIL*EPSTIL) ALFAC = UN - ALFAT C C modification pour la bi ou tricompression C IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+ 1 SIGPC(3)*SIGPC(3) GAMMA=-SQRT(GAMMA)/TRSIGC EPSTIM=EPSTIM*GAMMA END IF C C amelioration de la reponse en cisaillement pour beta > 1. C IF (BETA .GT. UN) THEN IF ( ALFAT .GT. XPETIT ) THEN ALFAT=ALFAT**BETA END IF IF ( ALFAC .GT. XPETIT ) THEN ALFAC=ALFAC**BETA END IF END IF C C on calcule DT et DC puis D C dans le cas ou le seuil initial est depasse C C on est oblige de verifier car on a pu multiplier par gamma C IF (EPSTIM .GT. EPSD0) THEN DT=UN - EPSD0*(UN-ATRA)/EPSTIM - 1 ATRA*EXP(-BTRA*(EPSTIM-EPSD0)) DC=UN - EPSD0*(UN-ACOM)/EPSTIM - 1 ACOM*EXP(-BCOM*(EPSTIM-EPSD0)) ELSE DT=XZERO DC=XZERO END IF D = ALFAT*DT + ALFAC*DC C C on borne la valeur de D a 0.99999 C D=MIN ( D , UN-1.D-05 ) ELSE D=XZERO END IF C C on teste la croissance de D C D=MAX ( D , DINI ) C C on le stocke dans les variables internes finales C VARF(2)= D C C on calcule les contraintes finales C CALL MATVE1 (DDHOOK,EPSILO,NSTRS,NSTRS,SIGF,2) DO 500 ISTRS=1,NSTRS SIGF(ISTRS)=SIGF(ISTRS)*(UN-D) 500 CONTINUE C C et les deformations inelastiques finales C DO 600 ISTRS=1,NSTRS EPINF(ISTRS)=EPSILO(ISTRS)*D 600 CONTINUE END IF SEGSUP WRK3 RETURN END