Télécharger la plaquette Cast3M


fr en

Télécharger cicsic.eso

Retour à la liste

Numérotation des lignes : 

cicsic

  1. C CICSIC    SOURCE    OF166741  25/11/04    21:15:22     12349          
  2.  
  3.       SUBROUTINE CICSIC(WRK52,WRK53,WRK54,WRK22,IB,IGAU,NVARI,NBPGAU
  4.      & ,necou,iecou)
  5. C
  6. C     variables en entree
  7. C
  8. C     WRK0,KRK1,WRK5  pointeurs sur des segments de travail
  9. C
  10. C     WRK0:
  11. C            XMAT(NMATT): tableau composantes materiaux
  12. C
  13. C     WRK1:
  14. C            SIG0(NSTRS) : contraintes au debut du pas
  15. C            VAR0(NVARI) : variables internes au debut du pas
  16. C            DEPST(NSTRS): increment de deformation totale
  17. C
  18. C     WRK22:
  19. C            XXE: coordonnees de l'element en double precison
  20. C
  21. C     WTRAV:
  22. C            TXR: cosinus directeurs des axes locaux de l'element massif
  23. C
  24. C     WRK5:
  25. C
  26. C     NSTRS      nombre de composantes dans les vecteurs des contraintes
  27. C                                        et les vecteurs des deformations
  28. C
  29. C     NVARI      nombre de variables internes (doit etre egal a 3)
  30. C
  31. C     NMATT      nombre de constantes du materiau
  32. C
  33. C     variables en sortie
  34. C
  35. C     VARF      variables internes finales dans WRK0
  36. C
  37. C     SIGF      contraintes finales dans WRK0
  38. C
  39.       IMPLICIT INTEGER(I-N)
  40.       IMPLICIT REAL*8(A-H,O-Z)
  41.  
  42. -INC PPARAM
  43. -INC CCOPTIO
  44. -INC CCREEL
  45. -INC DECHE
  46. -INC TECOU
  47.  
  48.       SEGMENT WRK22
  49.         REAL*8 XXE(3,NBNNbi)
  50.       ENDSEGMENT
  51.  
  52.       SEGMENT WRK6
  53.         REAL*8 SIG0S(NSTRS),DEPSTS(NSTRS)
  54.       END SEGMENT
  55.  
  56.       INTEGER NSTRS1,NVARI
  57.       INTEGER KCAS,IRTD,ISTRS
  58.  
  59.       REAL*8 IDAUX(6,6)
  60.       REAL*8 DOR1(3), DOR2(3), DOR3(3)
  61.       REAL*8 MORTH (3,3), IORTH (3,3), AEQ (3,3)
  62.       REAL*8 H1 (6,6), H2(6,6), H3(6,6)
  63.       REAL*8 H01 (6,6), H02(6,6), H03(6,6)
  64.       REAL*8 K01 (6,6), K02(6,6), K03(6,6), KW(6,6)
  65.       REAL*8 KEFF1(6,6), KEFF2(6,6), KEFF3(6,6)
  66.       REAL*8 RIG0 (6,6), CED0 (6,6)
  67.       REAL*8 SORTH(6), SIG0V(6), SIGFV(6),
  68.      &       DEFI(6), DEPSTV(6), DORTH(6)
  69.       REAL*8 VARTMP
  70.       INTEGER HCLO1, HCLO2, HCLO3
  71.  
  72. * scalar damage indicators
  73.       REAL*8 DOM1, DOM2, DOM3
  74.       REAL*8 YGR1, YGR2, YGR3, YW(6)
  75.       REAL*8 YEQ1, YEQ2, YEQ3
  76. *
  77. *    PHENOMENOLOGICAL COEFFICIENTS
  78. *
  79. * deactivation parameter (scalars)
  80.       REAL*8 ETADS
  81. * closure stress
  82. *      REAL*8 SIGCLO (NSTRSS)
  83. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  84. *ON NE LE PREND PAS EN COMPTE!!!!
  85. *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  86.  
  87. * scalar damage law parameters
  88. * Passes comme param. materiau
  89.       REAL*8 G1DC, G1Y0, G1YC, G1P
  90.       REAL*8 G2DC, G2Y0, G2YC, G2P
  91.       REAL*8 G3DC, G3Y0, G3YC, G3P
  92.  
  93.       DATA ETADS /1./
  94.  
  95.       DATA H1
  96.      &    /1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  97.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  98.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  99.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  100.      &     0.0, 0.0, 0.0, 0.0, 0.7, 0.0,
  101.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.7/
  102.       DATA H2
  103.      &    /0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  104.      &     0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
  105.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  106.      &     0.0, 0.0, 0.0, 0.7, 0.0, 0.0,
  107.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  108.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.7/
  109.       DATA H3
  110.      &    /0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  111.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  112.      &     0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
  113.      &     0.0, 0.0, 0.0, 0.7, 0.0, 0.0,
  114.      &     0.0, 0.0, 0.0, 0.0, 0.7, 0.0,
  115.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0/
  116.       DATA AEQ
  117.      &    /1.0, 0.0, 0.0,
  118.      &     0.0, 1.0, 0.0,
  119.      &     0.0, 0.0, 1.0/
  120.  
  121. ***** INITIALIZATION VARIABLES
  122. *
  123.        kerre=0
  124. * Proprietes materiau
  125.        YG1   = XMAT(1)
  126.        YG2   = XMAT(2)
  127.        YG3   = XMAT(3)
  128.        XNU12 = XMAT(4)
  129.        XNU23 = XMAT(5)
  130.        XNU13 = XMAT(6)
  131.        G12   = XMAT(7)
  132.        G23   = XMAT(8)
  133.        G13   = XMAT(9)
  134.  
  135.        G1DC  = XMAT(16)
  136.        G1Y0  = XMAT(17)
  137.        G1YC  = XMAT(18)
  138.        G1P   = XMAT(19)
  139.        G2DC  = XMAT(20)
  140.        G2Y0  = XMAT(21)
  141.        G2YC  = XMAT(22)
  142.        G2P   = XMAT(23)
  143.        G3DC  = XMAT(24)
  144.        G3Y0  = XMAT(25)
  145.        G3YC  = XMAT(26)
  146.        G3P   = XMAT(27)
  147.  
  148.        H1(5,5)= XMAT(29)
  149.        H3(5,5)= XMAT(29)
  150.  
  151.        H1(6,6)= XMAT(30)
  152.        H2(6,6)= XMAT(30)
  153.  
  154.        H2(4,4)= XMAT(28)
  155.        H3(4,4)= XMAT(28)
  156.  
  157.        AEQ(1,2)= XMAT(35)
  158.        AEQ(2,1)= XMAT(35)
  159.  
  160.        AEQ(1,3)= XMAT(37)
  161.        AEQ(3,1)= XMAT(37)
  162.  
  163.        AEQ(2,3)= XMAT(36)
  164.        AEQ(3,2)= XMAT(36)
  165.  
  166. * Pour deboggage
  167. C*      jPARAM=0
  168.  
  169. * Variables internes du modele
  170.       DOM1=VAR0(2)
  171.       DOM2=VAR0(3)
  172.       DOM3=VAR0(4)
  173.  
  174. * Mise a zero des matrices
  175.       CALL ZERO (RIG0,6,6)
  176.       CALL ZERO (CED0,6,6)
  177.  
  178. ****** CALCUL DE LA DEFORMATION INITIALE A PARTIR
  179. ****** DES CONTRAINTES INITIALES
  180.  
  181. **********Debougage
  182. C*      IF (jPARAM.EQ.1) THEN
  183. C*         WRITE(IOIMP,66770) IB,IGAU
  184. C*66770    format(////////2x,'element  ',i6,2x,'point  ',i3//)
  185. C*
  186. C*         WRITE (IOIMP,*) 'Increment des deformations '
  187. C*         WRITE (IOIMP,99999) (DEPST(ILOOP), ILOOP=1,6)
  188. C*
  189. C*         WRITE (IOIMP,*) 'Contraintes au debut de l''iteration'
  190. C*         WRITE (IOIMP,99999) (SIG0(ILOOP), ILOOP=1,6)
  191. C*      ENDIF
  192. **********Debougage
  193.  
  194. * Controle si il faut calculer la matrice de hooke
  195.  
  196.       IF (IB.EQ.1.AND.IGAU.EQ.1) THEN
  197.          GOTO 100
  198.  
  199.       ELSEIF (N2PTEL.EQ.1.AND.N2EL.EQ.1) THEN
  200.          GOTO 200
  201.  
  202.       ELSEIF (N2PTEL.EQ.1.AND.N2EL.NE.1)  THEN
  203.          IF (IGAU.EQ.1) THEN
  204.             GOTO 100
  205.          ELSE
  206.             GOTO 200
  207.          ENDIF
  208.  
  209.       ELSE
  210.          GOTO 100
  211.       ENDIF
  212.  
  213. * Calcul de la matrice de hook pour un materiau orthotrope
  214. 100   CONTINUE
  215. *      WRITE (IOIMP,*) 'Calcul de la matrice de Hooke, CMATE=', CMATE
  216.  
  217.       IPERR=1
  218.       CALL ZERO (DDAUX,LHOOK,LHOOK)
  219.  
  220.       XNU21=(YG2/YG1)*XNU12
  221.       XNU32=(YG3/YG2)*XNU23
  222.       XNU31=(YG3/YG1)*XNU13
  223.       AUX=(1.D0-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
  224.      &    -2.D0*XNU21*XNU32*XNU13)
  225.       AUX1=AUX/YG1
  226.       AUX2=AUX/YG2
  227.       AUX3=AUX/YG3
  228.  
  229.       DDAUX(1,1)=(1.D0-XNU23*XNU32)/AUX1
  230.       DDAUX(1,2)=(XNU21+XNU31*XNU23)/AUX1
  231.       DDAUX(2,1)=DDAUX(1,2)
  232.       DDAUX(1,3)=(XNU31+XNU21*XNU32)/AUX1
  233.       DDAUX(3,1)=DDAUX(1,3)
  234.       DDAUX(2,2)=(1.D0-XNU13*XNU31)/AUX2
  235.       DDAUX(2,3)=(XNU32+XNU12*XNU31)/AUX2
  236.       DDAUX(3,2)=DDAUX(2,3)
  237.       DDAUX(3,3)=(1.D0-XNU12*XNU21)/AUX3
  238.       DDAUX(4,4)=G23
  239.       DDAUX(5,5)=G13
  240.       DDAUX(6,6)=G12
  241.  
  242. * On recopie la matrice de hook pour l'inversee
  243. 200   DO JLOOP=1,LHOOK
  244.          DO ILOOP=1,LHOOK
  245.             IDAUX(ILOOP,JLOOP)= DDAUX(ILOOP,JLOOP)
  246.          ENDDO
  247.       ENDDO
  248.  
  249. * Inversion de la matrice de Hooke
  250. *      WRITE (IOIMP,*) 'Inversion de la matrice de Hooke'
  251.       TPREC= 0.D0
  252.       IPERR=0
  253.       CALL INVALM(IDAUX,LHOOK,LHOOK,IPERR,TPREC)
  254.       IF (IPERR.NE.0) THEN
  255.         WRITE (IOIMP,*) 'ERREUR DANS L''INVERSION DE LA MATRICE DE HOOK'
  256.       ENDIF
  257.  
  258. * Calcul des tenseurs d'endommagement H01, H02, H03
  259.       CALL SICTEN (H1,IDAUX,H01)
  260.       CALL SICTEN (H2,IDAUX,H02)
  261.       CALL SICTEN (H3,IDAUX,H03)
  262.  
  263. C* Calcul des matrices K01,K02,K03
  264.         CALL MULMAT (KW,H01,DDAUX,LHOOK,LHOOK,LHOOK)
  265.         CALL MULMAT (K01,DDAUX,KW,LHOOK,LHOOK,LHOOK)
  266.  
  267.         CALL MULMAT (KW,H02,DDAUX,LHOOK,LHOOK,LHOOK)
  268.         CALL MULMAT (K02,DDAUX,KW,LHOOK,LHOOK,LHOOK)
  269.  
  270.         CALL MULMAT (KW,H03,DDAUX,LHOOK,LHOOK,LHOOK)
  271.         CALL MULMAT (K03,DDAUX,KW,LHOOK,LHOOK,LHOOK)
  272.  
  273. ***** Recuperation deformation initiale dans le repere orth.
  274.  
  275.       DEFI(1)=VAR0(5)
  276.       DEFI(2)=VAR0(6)
  277.       DEFI(3)=VAR0(7)
  278.       DEFI(4)=VAR0(8)
  279.       DEFI(5)=VAR0(9)
  280.       DEFI(6)=VAR0(10)
  281.  
  282. * Calcul de l'increment de deformation dans le repere orthotrope
  283.       DO 451 ILOOP=1,6
  284.          DEPSTV(ILOOP)=DEPST(ILOOP)
  285. 451   CONTINUE
  286.  
  287.       CALL CICROT (wrk52,wrk53,wrk54,1,DEPSTV,DORTH)
  288.  
  289. * Reorganisation dans SICDEF selon la valeur de iAXEP
  290. *  TC  iaxep n'est pas initialise , je le mets a 0 !!!
  291.       iaxep=0
  292.       CALL SICDEF (DORTH,iAXEP,kerre)
  293. * Reorganisation supplementaire pour le contraintes en cisaillement
  294.       VARTMP= DORTH(4)
  295.       DORTH(4)= DORTH(6)
  296.       DORTH(6)= VARTMP
  297.       DORTH(5)= DORTH(5)
  298.  
  299. * On calcule les deformation totales.
  300.       DO 500 ILOOP=1,NSTRSS
  301.         DORTH(ILOOP)=DEFI(ILOOP)+DORTH(ILOOP)
  302. *       On garde le deformations totales dans le var. internes
  303.         VARF(ILOOP+4)=DORTH(ILOOP)
  304. 500   CONTINUE
  305.  
  306. ********************************************************
  307. **           ON A LES DEFORMATIONS TOTALES            **
  308. **  ON PEUT CALCULER LES CONTRAINTES A LA FIN DU PAS  **
  309. ********************************************************
  310. *      WRITE (IOIMP,*) 'Calcul contraintes'
  311.  
  312. * Calcul des noveaux HCLO
  313.  
  314.       IF (DORTH(1).LT.0.) THEN
  315.          HCLO1=1
  316.       ELSE
  317.          HCLO1=0
  318.       ENDIF
  319.       IF (DORTH(2).LT.0.) THEN
  320.          HCLO2=1
  321.       ELSE
  322.          HCLO2=0
  323.       ENDIF
  324.       IF (DORTH(3).LT.0.) THEN
  325.          HCLO3=1
  326.       ELSE
  327.          HCLO3=0
  328.       ENDIF
  329.  
  330. * Calcul des nouvelles KEFF
  331.       DO JLOOP=1,LHOOK
  332.         DO ILOOP=1,LHOOK
  333.            KEFF1(ILOOP,JLOOP)=K01(ILOOP,JLOOP)
  334.            KEFF2(ILOOP,JLOOP)=K02(ILOOP,JLOOP)
  335.            KEFF3(ILOOP,JLOOP)=K03(ILOOP,JLOOP)
  336.         ENDDO
  337.       ENDDO
  338.       KEFF1(1,1)=K01(1,1)*(1.-ETADS*HCLO1)
  339.       KEFF2(2,2)=K02(2,2)*(1.-ETADS*HCLO2)
  340.       KEFF3(3,3)=K03(3,3)*(1.-ETADS*HCLO3)
  341.  
  342. * Calcul des YGR
  343. *      WRITE (IOIMP,*) 'YW1'
  344.       CALL ZERO (YW,6,1)
  345.       CALL MULMAT (YW,KEFF1,DORTH,6,1,6)
  346. *      WRITE (IOIMP,*) (YW(ILOOP), ILOOP=1,6)
  347.       YGR1=0.
  348.       DO 510 ILOOP=1,6
  349.          YGR1=YGR1+YW(ILOOP)*DORTH(ILOOP)
  350. 510   CONTINUE
  351.       YGR1=0.5D0*YGR1
  352. *      WRITE (IOIMP,*) 'YGR1=',YGR1
  353.  
  354.       CALL ZERO (YW,6,1)
  355.       CALL MULMAT (YW,KEFF2,DORTH,6,1,6)
  356. *      WRITE (IOIMP,*) 'YW2'
  357. *      WRITE (IOIMP,*) (YW(ILOOP), ILOOP=1,6)
  358.       YGR2=0.
  359.       DO 520 ILOOP=1,6
  360.          YGR2=YGR2+YW(ILOOP)*DORTH(ILOOP)
  361. 520   CONTINUE
  362.       YGR2=0.5D0*YGR2
  363.  
  364.       CALL ZERO (YW,6,1)
  365.       CALL MULMAT (YW,KEFF3,DORTH,6,1,6)
  366.       YGR3=0.
  367.       DO 530 ILOOP=1,6
  368.          YGR3=YGR3+YW(ILOOP)*DORTH(ILOOP)
  369. 530   CONTINUE
  370.       YGR3=0.5D0*YGR3
  371.  
  372.       YGR1 = MAX(YGR1,0.D0)
  373.       YGR2 = MAX(YGR2,0.D0)
  374.       YGR3 = MAX(YGR3,0.D0)
  375.  
  376. *Calcul des YEQ
  377.       YEQ1= AEQ(1,1)*YGR1+AEQ(1,2)*YGR2+AEQ(1,3)*YGR3
  378.       YEQ2= AEQ(2,1)*YGR1+AEQ(2,2)*YGR2+AEQ(2,3)*YGR3
  379.       YEQ3= AEQ(3,1)*YGR1+AEQ(3,2)*YGR2+AEQ(3,3)*YGR3
  380.  
  381. *Calcul des DOM
  382.  
  383.       DD1= ((SQRT(YEQ1) - G1Y0)/ G1YC )
  384.       DD2= ((SQRT(YEQ2) - G2Y0)/ G2YC )
  385.       DD3= ((SQRT(YEQ3) - G3Y0)/ G3YC )
  386.  
  387.       DD1 = MAX(DD1,0.D0)
  388.       DD2 = MAX(DD2,0.D0)
  389.       DD3 = MAX(DD3,0.D0)
  390.  
  391.       DOM1= G1DC*(1.D0 - EXP(-1.* ( DD1**G1P ) ) )
  392.       DOM2= G2DC*(1.D0 - EXP(-1.* ( DD2**G2P ) ) )
  393.       DOM3= G3DC*(1.D0 - EXP(-1.* ( DD3**G3P ) ) )
  394.  
  395.       DOM1= MAX(DOM1,VAR0(2))
  396.       DOM2= MAX(DOM2,VAR0(3))
  397.       DOM3= MAX(DOM3,VAR0(4))
  398.  
  399. *Calcul de la matrice de rigidite
  400.  
  401.       DO 600 JLOOP=1,LHOOK
  402.         DO ILOOP=1,LHOOK
  403.            RIG0(ILOOP,JLOOP)=DDAUX(ILOOP,JLOOP)
  404.      &                       -DOM1*KEFF1(ILOOP,JLOOP)
  405.      &                       -DOM2*KEFF2(ILOOP,JLOOP)
  406.      &                       -DOM3*KEFF3(ILOOP,JLOOP)
  407.         ENDDO
  408. 600   CONTINUE
  409.  
  410. * Calcul des contraintes a la fin du pas
  411. * SORTH contient les contraintes a la fin du pas
  412. * dans le repere orthotrope
  413.       CALL MULMAT (SORTH,RIG0,DORTH,6,1,6)
  414.  
  415. * ATTENTION:
  416. * On reorganise les contraintes en cisaillement pour
  417. * calculer les contraintes dans le repere global
  418.       SIGFV(1)=SORTH(1)
  419.       SIGFV(2)=SORTH(2)
  420.       SIGFV(3)=SORTH(3)
  421.  
  422.       SIGFV(4)=SORTH(6)
  423.       SIGFV(5)=SORTH(5)
  424.       SIGFV(6)=SORTH(4)
  425.  
  426.       CALL SICCNT (SIGFV,iAXEP,kerre)
  427.  
  428. * On appele SICROT pour trouver le contraintes dans le repere global
  429. * On utilise SORTH pour garder le resultat
  430. * (on peut pas passer SIGF a la subroutine car il est dans un segment)
  431.       CALL CICROT (wrk52,wrk53,wrk54,0,SIGFV,SORTH)
  432.  
  433. * On recopie SORTH in SIGF
  434.       DO 1000 ILOOP=1,6
  435.          SIGF(ILOOP)=SORTH(ILOOP)
  436. 1000  CONTINUE
  437.  
  438.       VARF(2)=  DOM1
  439.       VARF(3)=  DOM2
  440.       VARF(4)=  DOM3
  441.  
  442. 99999    format(2x,'  ',(6(1x,3pe12.5),/))
  443. 99998    format(2x,'  ',(4(1x,3pe12.5),/))
  444. 99997    format(2x,'  ',(1x,3pe12.5),/)
  445. 99996    format(2x,'  ',(3(1x,1pe12.3),/))
  446.  
  447.       RETURN
  448.       END
  449.  
  450.  
  451.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales