C CICSIC SOURCE OF166741 25/11/04 21:15:22 12349 SUBROUTINE CICSIC(WRK52,WRK53,WRK54,WRK22,IB,IGAU,NVARI,NBPGAU & ,necou,iecou) C C variables en entree C C WRK0,KRK1,WRK5 pointeurs sur des segments de travail C C WRK0: C XMAT(NMATT): tableau composantes materiaux C C WRK1: C SIG0(NSTRS) : contraintes au debut du pas C VAR0(NVARI) : variables internes au debut du pas C DEPST(NSTRS): increment de deformation totale C C WRK22: C XXE: coordonnees de l'element en double precison C C WTRAV: C TXR: cosinus directeurs des axes locaux de l'element massif C C WRK5: 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 3) C C NMATT nombre de constantes du materiau C C variables en sortie C C VARF variables internes finales dans WRK0 C C SIGF contraintes finales dans WRK0 C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC DECHE -INC TECOU SEGMENT WRK22 REAL*8 XXE(3,NBNNbi) ENDSEGMENT SEGMENT WRK6 REAL*8 SIG0S(NSTRS),DEPSTS(NSTRS) END SEGMENT INTEGER NSTRS1,NVARI INTEGER KCAS,IRTD,ISTRS REAL*8 IDAUX(6,6) REAL*8 DOR1(3), DOR2(3), DOR3(3) REAL*8 MORTH (3,3), IORTH (3,3), AEQ (3,3) REAL*8 H1 (6,6), H2(6,6), H3(6,6) REAL*8 H01 (6,6), H02(6,6), H03(6,6) REAL*8 K01 (6,6), K02(6,6), K03(6,6), KW(6,6) REAL*8 KEFF1(6,6), KEFF2(6,6), KEFF3(6,6) REAL*8 RIG0 (6,6), CED0 (6,6) REAL*8 SORTH(6), SIG0V(6), SIGFV(6), & DEFI(6), DEPSTV(6), DORTH(6) REAL*8 VARTMP INTEGER HCLO1, HCLO2, HCLO3 * scalar damage indicators REAL*8 DOM1, DOM2, DOM3 REAL*8 YGR1, YGR2, YGR3, YW(6) REAL*8 YEQ1, YEQ2, YEQ3 * * PHENOMENOLOGICAL COEFFICIENTS * * deactivation parameter (scalars) REAL*8 ETADS * closure stress * REAL*8 SIGCLO (NSTRSS) *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! *ON NE LE PREND PAS EN COMPTE!!!! *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! * scalar damage law parameters * Passes comme param. materiau REAL*8 G1DC, G1Y0, G1YC, G1P REAL*8 G2DC, G2Y0, G2YC, G2P REAL*8 G3DC, G3Y0, G3YC, G3P DATA ETADS /1./ DATA H1 & /1.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.7, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.7/ DATA H2 & /0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.7, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.7/ DATA H3 & /0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.7, 0.0, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.7, 0.0, & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/ DATA AEQ & /1.0, 0.0, 0.0, & 0.0, 1.0, 0.0, & 0.0, 0.0, 1.0/ ***** INITIALIZATION VARIABLES * kerre=0 * Proprietes materiau YG1 = XMAT(1) YG2 = XMAT(2) YG3 = XMAT(3) XNU12 = XMAT(4) XNU23 = XMAT(5) XNU13 = XMAT(6) G12 = XMAT(7) G23 = XMAT(8) G13 = XMAT(9) G1DC = XMAT(16) G1Y0 = XMAT(17) G1YC = XMAT(18) G1P = XMAT(19) G2DC = XMAT(20) G2Y0 = XMAT(21) G2YC = XMAT(22) G2P = XMAT(23) G3DC = XMAT(24) G3Y0 = XMAT(25) G3YC = XMAT(26) G3P = XMAT(27) H1(5,5)= XMAT(29) H3(5,5)= XMAT(29) H1(6,6)= XMAT(30) H2(6,6)= XMAT(30) H2(4,4)= XMAT(28) H3(4,4)= XMAT(28) AEQ(1,2)= XMAT(35) AEQ(2,1)= XMAT(35) AEQ(1,3)= XMAT(37) AEQ(3,1)= XMAT(37) AEQ(2,3)= XMAT(36) AEQ(3,2)= XMAT(36) * Pour deboggage C* jPARAM=0 * Variables internes du modele DOM1=VAR0(2) DOM2=VAR0(3) DOM3=VAR0(4) * Mise a zero des matrices CALL ZERO (RIG0,6,6) CALL ZERO (CED0,6,6) ****** CALCUL DE LA DEFORMATION INITIALE A PARTIR ****** DES CONTRAINTES INITIALES **********Debougage C* IF (jPARAM.EQ.1) THEN C* WRITE(IOIMP,66770) IB,IGAU C*66770 format(////////2x,'element ',i6,2x,'point ',i3//) C* C* WRITE (IOIMP,*) 'Increment des deformations ' C* WRITE (IOIMP,99999) (DEPST(ILOOP), ILOOP=1,6) C* C* WRITE (IOIMP,*) 'Contraintes au debut de l''iteration' C* WRITE (IOIMP,99999) (SIG0(ILOOP), ILOOP=1,6) C* ENDIF **********Debougage * Controle si il faut calculer la matrice de hooke IF (IB.EQ.1.AND.IGAU.EQ.1) THEN GOTO 100 ELSEIF (N2PTEL.EQ.1.AND.N2EL.EQ.1) THEN GOTO 200 ELSEIF (N2PTEL.EQ.1.AND.N2EL.NE.1) THEN IF (IGAU.EQ.1) THEN GOTO 100 ELSE GOTO 200 ENDIF ELSE GOTO 100 ENDIF * Calcul de la matrice de hook pour un materiau orthotrope 100 CONTINUE * WRITE (IOIMP,*) 'Calcul de la matrice de Hooke, CMATE=', CMATE IPERR=1 CALL ZERO (DDAUX,LHOOK,LHOOK) XNU21=(YG2/YG1)*XNU12 XNU32=(YG3/YG2)*XNU23 XNU31=(YG3/YG1)*XNU13 AUX=(1.D0-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 & -2.D0*XNU21*XNU32*XNU13) AUX1=AUX/YG1 AUX2=AUX/YG2 AUX3=AUX/YG3 DDAUX(1,1)=(1.D0-XNU23*XNU32)/AUX1 DDAUX(1,2)=(XNU21+XNU31*XNU23)/AUX1 DDAUX(2,1)=DDAUX(1,2) DDAUX(1,3)=(XNU31+XNU21*XNU32)/AUX1 DDAUX(3,1)=DDAUX(1,3) DDAUX(2,2)=(1.D0-XNU13*XNU31)/AUX2 DDAUX(2,3)=(XNU32+XNU12*XNU31)/AUX2 DDAUX(3,2)=DDAUX(2,3) DDAUX(3,3)=(1.D0-XNU12*XNU21)/AUX3 DDAUX(4,4)=G23 DDAUX(5,5)=G13 DDAUX(6,6)=G12 * On recopie la matrice de hook pour l'inversee 200 DO JLOOP=1,LHOOK DO ILOOP=1,LHOOK IDAUX(ILOOP,JLOOP)= DDAUX(ILOOP,JLOOP) ENDDO ENDDO * Inversion de la matrice de Hooke * WRITE (IOIMP,*) 'Inversion de la matrice de Hooke' TPREC= 0.D0 IPERR=0 CALL INVALM(IDAUX,LHOOK,LHOOK,IPERR,TPREC) IF (IPERR.NE.0) THEN WRITE (IOIMP,*) 'ERREUR DANS L''INVERSION DE LA MATRICE DE HOOK' ENDIF * Calcul des tenseurs d'endommagement H01, H02, H03 CALL SICTEN (H1,IDAUX,H01) CALL SICTEN (H2,IDAUX,H02) CALL SICTEN (H3,IDAUX,H03) C* Calcul des matrices K01,K02,K03 CALL MULMAT (KW,H01,DDAUX,LHOOK,LHOOK,LHOOK) CALL MULMAT (K01,DDAUX,KW,LHOOK,LHOOK,LHOOK) CALL MULMAT (KW,H02,DDAUX,LHOOK,LHOOK,LHOOK) CALL MULMAT (K02,DDAUX,KW,LHOOK,LHOOK,LHOOK) CALL MULMAT (KW,H03,DDAUX,LHOOK,LHOOK,LHOOK) CALL MULMAT (K03,DDAUX,KW,LHOOK,LHOOK,LHOOK) ***** Recuperation deformation initiale dans le repere orth. DEFI(1)=VAR0(5) DEFI(2)=VAR0(6) DEFI(3)=VAR0(7) DEFI(4)=VAR0(8) DEFI(5)=VAR0(9) DEFI(6)=VAR0(10) * Calcul de l'increment de deformation dans le repere orthotrope DO 451 ILOOP=1,6 DEPSTV(ILOOP)=DEPST(ILOOP) 451 CONTINUE CALL CICROT (wrk52,wrk53,wrk54,1,DEPSTV,DORTH) * Reorganisation dans SICDEF selon la valeur de iAXEP * TC iaxep n'est pas initialise , je le mets a 0 !!! iaxep=0 CALL SICDEF (DORTH,iAXEP,kerre) * Reorganisation supplementaire pour le contraintes en cisaillement VARTMP= DORTH(4) DORTH(4)= DORTH(6) DORTH(6)= VARTMP DORTH(5)= DORTH(5) * On calcule les deformation totales. DO 500 ILOOP=1,NSTRSS DORTH(ILOOP)=DEFI(ILOOP)+DORTH(ILOOP) * On garde le deformations totales dans le var. internes VARF(ILOOP+4)=DORTH(ILOOP) 500 CONTINUE ******************************************************** ** ON A LES DEFORMATIONS TOTALES ** ** ON PEUT CALCULER LES CONTRAINTES A LA FIN DU PAS ** ******************************************************** * WRITE (IOIMP,*) 'Calcul contraintes' * Calcul des noveaux HCLO IF (DORTH(1).LT.0.) THEN HCLO1=1 ELSE HCLO1=0 ENDIF IF (DORTH(2).LT.0.) THEN HCLO2=1 ELSE HCLO2=0 ENDIF IF (DORTH(3).LT.0.) THEN HCLO3=1 ELSE HCLO3=0 ENDIF * Calcul des nouvelles KEFF DO JLOOP=1,LHOOK DO ILOOP=1,LHOOK KEFF1(ILOOP,JLOOP)=K01(ILOOP,JLOOP) KEFF2(ILOOP,JLOOP)=K02(ILOOP,JLOOP) KEFF3(ILOOP,JLOOP)=K03(ILOOP,JLOOP) ENDDO ENDDO KEFF1(1,1)=K01(1,1)*(1.-ETADS*HCLO1) KEFF2(2,2)=K02(2,2)*(1.-ETADS*HCLO2) KEFF3(3,3)=K03(3,3)*(1.-ETADS*HCLO3) * Calcul des YGR * WRITE (IOIMP,*) 'YW1' CALL ZERO (YW,6,1) CALL MULMAT (YW,KEFF1,DORTH,6,1,6) * WRITE (IOIMP,*) (YW(ILOOP), ILOOP=1,6) YGR1=0. DO 510 ILOOP=1,6 YGR1=YGR1+YW(ILOOP)*DORTH(ILOOP) 510 CONTINUE YGR1=0.5D0*YGR1 * WRITE (IOIMP,*) 'YGR1=',YGR1 CALL ZERO (YW,6,1) CALL MULMAT (YW,KEFF2,DORTH,6,1,6) * WRITE (IOIMP,*) 'YW2' * WRITE (IOIMP,*) (YW(ILOOP), ILOOP=1,6) YGR2=0. DO 520 ILOOP=1,6 YGR2=YGR2+YW(ILOOP)*DORTH(ILOOP) 520 CONTINUE YGR2=0.5D0*YGR2 CALL ZERO (YW,6,1) CALL MULMAT (YW,KEFF3,DORTH,6,1,6) YGR3=0. DO 530 ILOOP=1,6 YGR3=YGR3+YW(ILOOP)*DORTH(ILOOP) 530 CONTINUE YGR3=0.5D0*YGR3 YGR1 = MAX(YGR1,0.D0) YGR2 = MAX(YGR2,0.D0) YGR3 = MAX(YGR3,0.D0) *Calcul des YEQ YEQ1= AEQ(1,1)*YGR1+AEQ(1,2)*YGR2+AEQ(1,3)*YGR3 YEQ2= AEQ(2,1)*YGR1+AEQ(2,2)*YGR2+AEQ(2,3)*YGR3 YEQ3= AEQ(3,1)*YGR1+AEQ(3,2)*YGR2+AEQ(3,3)*YGR3 *Calcul des DOM DD1= ((SQRT(YEQ1) - G1Y0)/ G1YC ) DD2= ((SQRT(YEQ2) - G2Y0)/ G2YC ) DD3= ((SQRT(YEQ3) - G3Y0)/ G3YC ) DD1 = MAX(DD1,0.D0) DD2 = MAX(DD2,0.D0) DD3 = MAX(DD3,0.D0) DOM1= G1DC*(1.D0 - EXP(-1.* ( DD1**G1P ) ) ) DOM2= G2DC*(1.D0 - EXP(-1.* ( DD2**G2P ) ) ) DOM3= G3DC*(1.D0 - EXP(-1.* ( DD3**G3P ) ) ) DOM1= MAX(DOM1,VAR0(2)) DOM2= MAX(DOM2,VAR0(3)) DOM3= MAX(DOM3,VAR0(4)) *Calcul de la matrice de rigidite DO 600 JLOOP=1,LHOOK DO ILOOP=1,LHOOK RIG0(ILOOP,JLOOP)=DDAUX(ILOOP,JLOOP) & -DOM1*KEFF1(ILOOP,JLOOP) & -DOM2*KEFF2(ILOOP,JLOOP) & -DOM3*KEFF3(ILOOP,JLOOP) ENDDO 600 CONTINUE * Calcul des contraintes a la fin du pas * SORTH contient les contraintes a la fin du pas * dans le repere orthotrope CALL MULMAT (SORTH,RIG0,DORTH,6,1,6) * ATTENTION: * On reorganise les contraintes en cisaillement pour * calculer les contraintes dans le repere global SIGFV(1)=SORTH(1) SIGFV(2)=SORTH(2) SIGFV(3)=SORTH(3) SIGFV(4)=SORTH(6) SIGFV(5)=SORTH(5) SIGFV(6)=SORTH(4) CALL SICCNT (SIGFV,iAXEP,kerre) * On appele SICROT pour trouver le contraintes dans le repere global * On utilise SORTH pour garder le resultat * (on peut pas passer SIGF a la subroutine car il est dans un segment) CALL CICROT (wrk52,wrk53,wrk54,0,SIGFV,SORTH) * On recopie SORTH in SIGF DO 1000 ILOOP=1,6 SIGF(ILOOP)=SORTH(ILOOP) 1000 CONTINUE VARF(2)= DOM1 VARF(3)= DOM2 VARF(4)= DOM3 99999 format(2x,' ',(6(1x,3pe12.5),/)) 99998 format(2x,' ',(4(1x,3pe12.5),/)) 99997 format(2x,' ',(1x,3pe12.5),/) 99996 format(2x,' ',(3(1x,1pe12.3),/)) RETURN END