C CMICRO SOURCE OF166741 25/11/04 21:15:29 12349 SUBROUTINE CMICRO (WRK52,WRK53,WRK54,NVARI,iecou) * * modele d'endommagement microplan couple a la plasticite * C. La Borderie + S. Fichant Oct. 95 * routines utilisees: * micro1: plasticite nadai * IDECAL=3 DANS LE CAS ISO IDECAL=8 DANS LE CAS ANISO * jacob3: diagonalisation: * attention jacob3 modifie la matrice a diagonaliser!! * prodt et prodt2 * attention prodt2 ne fonctionne qu'avec la matrice des V. P. !! * IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) INTEGER NVARI INTEGER ISTRS,I,J REAL*8 YOUNG,XNU,EPSD0,BT,DEUXMU,ALFA REAL*8 DOM(6),SIGAN(6),TRSIG,DEF33(3,3),EPSIPP(3),VECP(3,3) REAL*8 DEFPT(6),DEFTOT(6),EPSITOP(3),EPSE1,EPSE REAL*8 VECPT(3,3),DOM33(3,3),DEFRPDE(3,3),DOMRPDE(3,3) REAL*8 D1,D2,D3,DOM3(3) REAL*8 SIGPP(3),SIGPM(3),SIG33(3,3),SIG33P(3,3),SIG33M(3,3) REAL*8 S33PRD(3,3),S33MRD(3,3) REAL*8 LAMBDAP(6), LAMBDAM(6) LOGICAL COMP INTEGER IDECAL -INC PPARAM -INC CCOPTIO -INC DECHE -INC TECOU SEGMENT WRKK1 REAL*8 DEFELA(NSTRS1) REAL*8 DDEFEL(NSTRS1) ENDSEGMENT * print*,'dans microp' * * on recupere les variable materielles * YOUNG=XMAT(1) XNU=XMAT(2) EPSD0=XMAT(5) BT=XMAT(6) ALFA=XMAT(9) DEUXMU=YOUNG/(1.D0+XNU) NSTRS1 = iecou.NSTRSS * * recuperation des variables internes d'endommagement * DO ISTRS=1,NSTRS1 DOM(ISTRS)=VAR0(2+ISTRS) END DO * * on ecoule plastiquement sur la contrainte effective * IDECAL=8 * print*,'appel a micro1' CALL CMICR1(wrk52,wrk53,wrk54,nstrs1,NVARI,IDECAL, & .false.,DEFPT,EPSE1,EPSE) IF (KERRE .NE. 0) THEN print*,'on n''a pas converge dans micro1' CALL CMICR1(wrk52,wrk53,wrk54,nstrs1,NVARI,IDECAL, & .true.,DEFPT,EPSE1,EPSE) RETURN ENDIF * * on ecoule en endommagement sur les deformations elastiques * * print*,'apres micro1' SEGINI WRKK1 * * calcul de l'increment deformations elastiques DDEFEL * 1) on calcule l'increment de deformations totales avec * l'increment de contraintes elastique DSIGT * 2) on retranche l'increment de deformations plastiques DEFP * TRSIG= DSIGT(1)+DSIGT(2)+DSIGT(3) DO ISTRS=1,3 DDEFEL(ISTRS)=( (1.D0+XNU)*DSIGT(ISTRS)-XNU*TRSIG)/YOUNG 1 - DEFP(ISTRS) END DO DO ISTRS=4,NSTRS1 DDEFEL(ISTRS)= (1.D0+XNU)*DSIGT(ISTRS)/YOUNG 1 - 0.5d0*DEFP(ISTRS) END DO * * on diagonalise l'increment de deformations elastiques: * 1) on met sous forme 3x3 avec endoca * 2) on diagonalise avec jacob3 * * print*,'incr de def el dans rpg' * print*,ddefel * print*,'avant endoca ddefl' CALL ENDOCA (DDEFEL,DEF33,1) * print*,'apres endoca ddefl' CALL JACOB3 (DEF33,IDIM,EPSIPP,VECP) * print*,'apres jacob3 def33' * print*,'incr deformations principales' * print*,epsipp * * calcul des deformations elastiques * * print*,'-----sigf------' * print*,sigf(1),sigf(2),sigf(3) TRSIG=SIGF(1)+SIGF(2)+SIGF(3) DO ISTRS=1,3 DEFELA(ISTRS)=( (1.D0+XNU)*SIGF(ISTRS)-XNU*TRSIG)/YOUNG END DO * print*,'-----defela------' * print*,defela(1),defela(2),defela(3) DO ISTRS=4,NSTRS1 DEFELA(ISTRS)= (1.D0+XNU)*SIGF(ISTRS)/YOUNG END DO * print*,'deformations elastiques dans rpg' * print*,defela * * on met les deformations sous forme de matrice 3x3 * puis on ecrit la matrice dans le repere de depsilon:DEFRPDE * CALL ENDOCA (DEFELA,DEF33,1) * print*,'deformations elastiques dans rpg(3x3)' * print*,def33 CALL PRODT (DEFRPDE,DEF33,VECP,3,3) * print*,'def elast dans rpddeps' * print*,defrpde * * on calcule l'endommagement resultant * IF (DEFRPDE(1,1) .GT. EPSD0) THEN D1=1.D0-EPSD0/DEFRPDE(1,1)*EXP(BT*( EPSD0 - DEFRPDE(1,1))) ELSE D1=0.D0 END IF IF (DEFRPDE(2,2) .GT. EPSD0) THEN D2=1.D0-EPSD0/DEFRPDE(2,2)*EXP(BT*( EPSD0 - DEFRPDE(2,2))) ELSE D2=0.D0 END IF IF (DEFRPDE(3,3) .GT. EPSD0) THEN D3=1.D0-EPSD0/DEFRPDE(3,3)*EXP(BT*( EPSD0 - DEFRPDE(3,3))) ELSE D3=0.D0 END IF * print*,defrpde(1,1),epsd0,D1 * print*,defrpde(2,2),epsd0,D2 * print*,defrpde(3,3),epsd0,D3 * * on met l'endommagement initial dans le meme repere * CALL ENDOCA(DOM,DOM33,1) * print*,'prodt DOMRPDE' CALL PRODT (DOMRPDE,DOM33,VECP,3,3) * print*,'endommagement initial dans rpddeps' * print*,domrpde * * et on en deduit l'increment d'endommagement dans RPDE * IF(d1.gt.domrpde(1,1))then domrpde(1,1)=d1 endif IF(d2.gt.domrpde(2,2))then domrpde(2,2)=d2 endif IF(d3.gt.domrpde(3,3))then domrpde(3,3)=d3 endif * print*,'endommagement final dans rpddeps' * print*,domrpde * * on remet D dans RPG * DO J=1,3 DO I=1,3 VECPT(I,J)=VECP(J,I) END DO END DO call prodt(DOM33,DOMRPDE,VECPT,3,3) * print*,'endommagement final dans rpg' * print*,dom33 * * on separe les contraintes effectives en + et - dans rpsigma * CALL ENDOCA (SIGF,SIG33,1) CALL JACOB3 (SIG33,3,SIGPP,VECP) * print*,'contraintes ppales' * print*,sigpp DO I=1,3 IF (SIGPP(I) .LT. 0.D0)THEN SIGPM(I)=SIGPP(I) SIGPP(I)=0.D0 ELSE SIGPM(I)=0.D0 END IF END DO CALL PRODT2(SIG33P,SIGPP,VECP,3) CALL PRODT2(SIG33M,SIGPM,VECP,3) * print*,'contraintes dans rpg' * print*,sig33p * print*,sig33m * * on met le tout dans le repere ppal d'endo * attention jacob3 modifie la matrice fournie * --> on passe une copie DO J=1,3 DO I=1,3 SIG33(I,J)=DOM33(I,J) END DO END DO CALL JACOB3(SIG33,3,DOM3,VECP) * print*,'endom dans rpd' * print*,dom3 ************** * A REVOIR EN ATENDANT MIEUX ON BORNE LES VALEURS PROPRES DE D ENTRE 0. ET 1.* ************** DO I=1,3 DOM3(I)=MAX(DOM3(I),0.D0) DOM3(I)=MIN(DOM3(I),1.D0-1.d-6) END DO CALL PRODT2 (DOM33,DOM3,VECP,3) ************** * FIN A REVOIR ************** * print*,'dom3=',dom3 * print*,'dom33=',dom33 CALL PRODT (S33PRD,SIG33P,VECP,3,3) CALL PRODT (S33MRD,SIG33M,VECP,3,3) * print*,'s33prd',s33prd(3,3) * print*,'s33mrd',s33mrd(3,3) * print*,'callambdap' COMP=.FALSE. CALL CLMBDA(DOM3(1),DOM3(2),DOM3(3),LAMBDAP,YOUNG,XNU,ALFA, 1 COMP) COMP=.TRUE. CALL CLMBDA(DOM3(1),DOM3(2),DOM3(3),LAMBDAM,YOUNG,XNU,ALFA, 1 COMP) * print*,'apres callambdap' SIG33P(1,1)=LAMBDAP(1)*S33PRD(1,1)+LAMBDAM(1)*S33MRD(1,1) SIG33P(1,2)=LAMBDAP(4)*S33PRD(1,2)+LAMBDAM(4)*S33MRD(1,2) SIG33P(1,3)=LAMBDAP(5)*S33PRD(1,3)+LAMBDAM(5)*S33MRD(1,3) SIG33P(2,1)=SIG33P(1,2) SIG33P(2,2)=LAMBDAP(2)*S33PRD(2,2)+LAMBDAM(2)*S33MRD(2,2) SIG33P(2,3)=LAMBDAP(6)*S33PRD(2,3)+LAMBDAM(6)*S33MRD(2,3) SIG33P(3,1)=SIG33P(1,3) SIG33P(3,2)=SIG33P(2,3) SIG33P(3,3)=LAMBDAP(3)*S33PRD(3,3)+LAMBDAM(3)*S33MRD(3,3) * print*,'contraintes dans RPD' * print*,sig33p * print*,'sig33p(1,2)',sig33p(1,2) * print*,'dsigt',dsigt(1),dsigt(2),dsigt(3) * print*,'dsigt',dsigt(4),dsigt(5),dsigt(6) * print*,'sigf',sigf(1),sigf(2),sigf(3) * print*,'sigf',sigf(4),sigf(5),sigf(6) * * on remet le tout dans le repere global * * print*,'matrice des vect proprs' * print*,vecp DO J=1,3 DO I=1,3 VECPT(I,J)=VECP(J,I) END DO END DO CALL PRODT (SIG33,SIG33P,VECPT,3,3) * print*,'contraintes dans rpg' * print*,sig33 * * on rend les contraintes et les variables internes finales * SIGAN(1)=SIGF(1)-SIG33(1,1) SIGF(1)=SIG33(1,1) VARF(3)=MAX(DOM33(1,1),0.d0) SIGAN(2)=SIGF(2)-SIG33(2,2) SIGF(2)=SIG33(2,2) VARF(4)=MAX(DOM33(2,2),0.d0) SIGAN(3)=SIGF(3)-SIG33(3,3) SIGF(3)=SIG33(3,3) VARF(5)=MAX(DOM33(3,3),0.d0) SIGAN(4)=SIGF(4)-SIG33(1,2) SIGF(4)=SIG33(1,2) VARF(6)=MAX(DOM33(1,2),0.d0) VARF(7)=MAX(DOM33(1,3),0.d0) VARF(8)=MAX(DOM33(2,3),0.d0) IF(IFOUR.GE.1.OR.IFOUR.LE.-3) THEN SIGAN(5)=SIGF(5)-SIG33(1,3) SIGF(5)=SIG33(1,3) SIGAN(6)=SIGF(6)-SIG33(2,3) SIGF(6)=SIG33(2,3) ELSE SIGAN(5)=0.D0 SIGAN(6)=0.D0 END IF DO ISTRS=1,6 VARF(ISTRS+8)=SIGAN(ISTRS) END DO * print*,'sigf',sigf(1),sigf(2),sigf(3) * print*,'sigf',sigf(4),sigf(5),sigf(6) * print*,'sigan',sigan(1),sigan(2),sigan(3) * print*,'sigan',sigan(4),sigan(5),sigan(6) SEGSUP WRKK1 RETURN END