C CMAZAR SOURCE MB234859 26/01/26 21:15:08 12457 SUBROUTINE CMAZAR (WRK52,WRK53,WRK54,WRKK2,NSTRS1,NVARI, 1 ICARA,JDIM,IFOUR2) C MAZARS SOURCE AM 98/12/23 21:38:30 3409 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 NSTRS1 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 -INC DECHE SEGMENT WRKK2 REAL*8 EPSILI(NSTRS1) REAL*8 EPSILO(NSTRS1) ENDSEGMENT INTEGER NSTRS1,NVARI 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,BETA,KSEUIL,EPS_C_CHAPE INTEGER ISTRS,JSTRS,KCAS,IRTD REAL*8 XZERO,UN,DEUX,XPETIT,PT_DROITE,EPS_C_BARRE REAL*8 DINI,D,DT,DC,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA REAL*8 FTRA,FCOM,GFT,GFC,LCAR REAL*8 EPC2,XK,EPCU,K1,K2,KMAX,KPIC REAL*8 EPSD0,ACOM,BCOM,ATRA,BTRA PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12) PARAMETER (YPETIT=1.D0, RA2=SQRT(DEUX)) C C recuperation des variables initiales dans les tableaux C YOUN = XMAT(1) XNU = XMAT(2) BETA = XMAT(5) DINI = VAR0(2) KSEUIL = XZERO EPS_C_CHAPE = -1.D0*UN * Numéro de la loi jnplas = wrk53.INPLAS * Mazars origine appel : 'MAZARS' IF (jnplas .EQ. 30) THEN ACOM = XMAT(6) BCOM = XMAT(7) ATRA = XMAT(8) BTRA = XMAT(9) EPSD0= XMAT(10) ELSE FTRA = XMAT(6) FCOM = XMAT(7) KPIC = XMAT(8) EPSD0 = FTRA / YOUN * Mazars RTC appel : 'MAZARS_RTC' -> 193 IF (jnplas .EQ. 193) THEN GFC = XMAT(9) LCAR = XMAT(10) GFT = XMAT(11) EPC2 = XMAT(12) XK = (1.05D0*YOUN*KPIC)/FCOM EPCU = (DEUX*GFC)/(LCAR*FCOM) 1 -(EPC2-KPIC) K1 = FCOM/(EPCU-EPC2) K2 = FCOM + (K1*EPC2) BTRA = LCAR * FTRA / 1 (GFT - (LCAR * YOUN * (EPSD0**DEUX) / DEUX)) * Mazars initial appel : 'MAZARS_INI' -> 194 ELSE IF (jnplas .EQ. 194) THEN SRES = XMAT(9) TAU = XMAT(10) BCOM = UN/(KPIC * XNU * RA2) ACOM = (FCOM/(KPIC*YOUN)-(BCOM*EPSD0)) / 1 ((EXP (BCOM*EPSD0 - UN))- (BCOM*EPSD0)) ATRA = UN - (SRES/(YOUN*EPSD0)) BTRA = (UN + TAU ) / EPSD0 * Mazars linéaire appel : 'MAZARS_LIN' -> 195 ELSE IF (jnplas .EQ. 195) THEN KMAX = XMAT(9) BCOM = UN/(KPIC * XNU * RA2) ACOM = (FCOM/(KPIC*YOUN)-(BCOM*EPSD0)) / 1 ((EXP (BCOM*EPSD0 - UN))- (BCOM*EPSD0)) * Mazars RT appel : 'MAZARS_RT' -> 196 ELSE IF (jnplas .EQ. 196) THEN GFT = XMAT(9) LCAR = XMAT(10) BCOM = UN / (KPIC * XNU * RA2) ACOM = (FCOM/(KPIC*YOUN)-(BCOM*EPSD0)) / 1 ((EXP (BCOM*EPSD0 - UN))- (BCOM*EPSD0)) BTRA = LCAR * FTRA / 1 (GFT - (LCAR * YOUN * (EPSD0**DEUX) / DEUX)) END IF END IF C C calcul des seuils C C calcul de la deformation totale C DO 100 ISTRS=1,NSTRS1 EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS) 100 CONTINUE C C calcul des deformations principales 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,NSTRS1,KCAS,DDHOOK,IRTD) * print*,'apres dohmas' C Fin modif L.Bode DO 200 ISTRS=1,3 r_z = XZERO DO 210 JSTRS=1,3 r_z = r_z + DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS) 210 CONTINUE SIGP(ISTRS)= r_z 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) )**DEUX + 1 MAX( XZERO , EPSIPP(2) )**DEUX + 2 MAX( XZERO , EPSIPP(3) )**DEUX EPSTIL=SQRT (EPSTIL) epstil=max(xpetit,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 on calcule ALFAT ALFAC DT et DC puis D C dans le cas ou le seuil initial est depasse C IF (jnplas .EQ. 193) THEN EPS_C_CHAPE = EPSTIM/(RA2*XNU) KSEUIL = (0.05D0*XK*KPIC/(1.05D0+XK*(XK - DEUX))) ENDIF IF (( EPSTIM .GT. EPSD0).OR.(EPS_C_CHAPE .GT. KSEUIL )) THEN C C calcul de l'endommagement 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. -YPETIT .AND. TRSIGT.LT.YPETIT) 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 trois lois d'évolution: C ATRA > 0 : loi de mazars classique C -10 BTRA est alors C la déformation pour laquelle la contrainte s'annule C C C Loi originellement implémentée IF (jnplas .EQ. 30) THEN IF (EPSTIM .GT. EPSD0) THEN IF (ATRA .GT. XZERO) THEN DT=UN - EPSD0*(UN-ATRA)/EPSTIM - 1 ATRA*EXP(BTRA*(EPSD0-EPSTIM)) ELSE IF (ATRA . GT. -10.D0) THEN DT=UN - epsd0/epstim*EXP(BTRA*(EPSD0-EPSTIM)) ELSE IF(EPSTIM .LT. BTRA) THEN DT=UN - EPSD0*(BTRA - EPSTIM)/EPSTIM/ 1 (BTRA - EPSD0) ELSE DT= UN ENDIF END IF DC=UN - EPSD0*(UN-ACOM)/EPSTIM - 1 ACOM*EXP(-BCOM*(EPSTIM-EPSD0)) ELSE DT=XZERO DC=XZERO END IF ELSE c Nouvelle implémentation : M de Carlan c Partie traction IF (jnplas .EQ. 193.OR.jnplas.EQ.196) THEN IF (EPSTIM .GT. EPSD0) THEN DT=UN - epsd0/epstim*EXP(BTRA*(EPSD0-EPSTIM)) ELSE DT=XZERO END IF ELSE IF (jnplas .EQ. 194) THEN IF (EPSTIM .GT. EPSD0) THEN DT=UN - EPSD0*(UN-ATRA)/EPSTIM - 1 ATRA*EXP(BTRA*(EPSD0-EPSTIM)) ELSE DT=XZERO END IF ELSE IF (EPSTIM .GT. EPSD0) THEN IF(EPSTIM .LT. KMAX) THEN DT=UN - EPSD0*(KMAX - EPSTIM)/EPSTIM/ 1 (KMAX - EPSD0) ELSE DT= UN ENDIF ELSE DT= XZERO END IF END IF c Partie compression IF (jnplas .EQ. 194.OR.jnplas.EQ.195 1 .OR.jnplas.EQ.196) THEN IF (EPSTIM .GT. EPSD0) THEN DC=UN - EPSD0*(UN-ACOM)/EPSTIM - 1 ACOM*EXP(-BCOM*(EPSTIM-EPSD0)) ELSE DC= XZERO END IF ELSE IF (EPS_C_CHAPE .GT. KSEUIL ) THEN EPS_C_BARRE = EPS_C_CHAPE / KPIC PT_DROITE = XZERO IF (EPS_C_CHAPE .LT. KPIC) THEN FACT = (XK*EPS_C_BARRE-EPS_C_BARRE*EPS_C_BARRE) 1 *FCOM DENOM=((UN+(XK-DEUX)*EPS_C_BARRE)*YOUN 1 *EPS_C_CHAPE) PT_DROITE = FACT / DENOM ELSE IF (EPS_C_CHAPE .LT. EPC2) THEN PT_DROITE = FCOM / (YOUN * EPS_C_CHAPE) ELSE IF (EPS_C_CHAPE .LT. EPCU) THEN PT_DROITE = K2 / (YOUN * EPS_C_CHAPE) - 1 K1 / YOUN ELSE PT_DROITE = XZERO END IF DC = UN - PT_DROITE ELSE DC = XZERO END IF END IF END IF DC = MAX(DC,XZERO) DT = MAX(DT,XZERO) ALFAT = MAX(ALFAT, XZERO) ALFAC = MAX(ALFAC, XZERO) D = ALFAT*DT + ALFAC*DC C C on borne la valeur de D a 0.99999999 pour avoir une raideur residuelle evitant des points flottants C *pv D=MIN ( D , UN-1.D-20 ) D=MIN ( D , UN-1.D-8 ) 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,NSTRS1,NSTRS1,SIGF,2) DO 500 ISTRS=1,NSTRS1 SIGF(ISTRS)=SIGF(ISTRS)*(UN-D) 500 CONTINUE C C et les deformations inelastiques finales C DO 600 ISTRS=1,NSTRS1 EPINF(ISTRS)=EPSILO(ISTRS)*D 600 CONTINUE ENDIF RETURN END