Télécharger la plaquette Cast3M


fr en

Télécharger ccoin0.eso

Retour à la liste

Numérotation des lignes : 

ccoin0

  1. C CCOIN0    SOURCE    OF166741  25/11/04    21:15:09     12349          
  2.  
  3.       SUBROUTINE CCOIN0(wrk52,wrk53,wrk54,wrk2,wrk3,
  4.      &  IB,IGAU,NBPGAU,LTRAC2,IFOUR2,iecou)
  5.  
  6. C-----------------------------------------------------------------------
  7. C               ECOULEMENT PLASTIQUE POUR UN POINT
  8. C               ALGORITHME ORTIZ ET SIMO
  9. C
  10. C      EN ENTREE :
  11. C
  12. C                  SIG0       CONTRAINTES AU DEBUT DU PAS
  13. C                  DEPST      INCREMENT DE DEFORMATIONS TOTALES
  14. C                               ( THERMIQUE ENLEVEE )
  15. C                  VAR0     VARIABLES INTERNES DEDUT DU PAS
  16. C                  VAREX0     VARIABLES EXTERNES DEBUT DU PAS
  17. C                  VAREXF     VARIABLES EXTERNES FIN DU PAS
  18. C                  XMAT       COEFFICIENTS DU MATERIAU
  19. C                  PRECIS     PRECISION POUR ITERATIONS INTERNES
  20. C                  WORK       DES CARACTERISTIQUES
  21. C                  TRAC       COURBE DE TRACTION
  22. C                  MFR1        INDICE DE FORMULATION (iecou)
  23. C                  NSTRSS      NOMBRE DE CONTRAINTES (iecou)
  24. C                  INPLAS     NUMERO DU MODELE DE PLASTICITE
  25. C                  DDAUX = MATRICE DE HOOKE ELASTIQUE
  26. C                  CMATE = NOM DU MATERIAU
  27. C                  VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
  28. C                  VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
  29. C                  N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
  30. C                  N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
  31. C                  IFOUR2  = OPTION DE CALCUL (redondant avec IFOUR, iecou.IFOURB
  32. C                  IB    = NUMERO DE L ELEMENT COURANT
  33. C                  IGAU  = NUMERO DU POINT COURANT
  34. C                  EPAIST= EPAISSEUR
  35. C                  NBPGAU= NBRE DE POINTS DE GAUSS
  36. C                  MELE  = NUMERO DE L ELEMENT FINI
  37. C                  NPINT = NBRE DE POINTS D INTEGRATION
  38. C                  NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES (iecou)
  39. C                  NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES (iecou)
  40. C                  SECT  = SECTION
  41. C                  LHOOK = TAILLE DE LA MATRICE DE HOOKE
  42. C                  TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES
  43. C                  UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE
  44. C
  45. C     EN SORTIE :
  46. C
  47. C                  SIGF        CONTRAINTES A LA FIN DU PAS
  48. C                  VARF      VARIABLES INTERNES FIN DU PAS
  49. C                  DEFP        INCR. DE DEFORMATIONS PLASTIQUES
  50. C                  KERRE       CODE D'ERREUR
  51. C                   = 0  SI TOUT OK
  52. C                   = 99 SI FORMULATION NON DISPONIBLE
  53. C    EN PLASTICITE
  54. C                   = 1  SI DEPS EST NEGATIF
  55. C                   = 2  SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE
  56. C
  57. C         IFOUR : OPTION DE CALCUL
  58. C
  59. C         IFOUR = -3   DEFORMATION PLANE GENERALISEE
  60. C         IFOUR = -2   CONTRAINTES PLANES
  61. C         IFOUR = -1   DEFORMATIONS PLANES
  62. C         IFOUR =  0   AXISYMETRIQUE
  63. C         IFOUR =  1   SERIE DE FOURIER
  64. C         IFOUR =  2   TRIDIMENSIONNEL
  65. C
  66. C         MFR1 : NUMERO DE LA FORMULATION ELEMENT FINI
  67. C
  68. C         MFR1 = 1     MASSIF
  69. C         MFR1 = 3     COQUE
  70. C         MFR1 = 5     COQUE EPAISSE
  71. C         MFR1 = 7     POUTRE
  72. C         MFR1 = 9     COQUE AVEC CISAILLEMENT TRANSVERSE
  73. Ckich     MFR1 = 31    pondération réduite termes diagonaux matrice B,
  74. C  dite MASSIF INCOMPRESSIBLE. Utilisation en contraintes planes a justifier
  75. c doublon historique MFR/MFR1
  76. C
  77. C         INPLAS : NUMERO DU MODELE DE PLASTICITE
  78. C
  79. C         INPLAS = 1      PARFAIT
  80. C         INPLAS = 4      CINEMATIQUE
  81. C         INPLAS = 5      ISOTROPE
  82. C         INPLAS = 7      CHABOCHE1
  83. C         INPLAS = 12     CHABOCHE2
  84. C
  85. C-----------------------------------------------------------------------
  86. C         CONVENTION DE REMPLISSAGE DES MEMOIRES
  87. C-----------------------------------------------------------------------
  88. C
  89. C     XMAT(1)    MODULE D'YOUNG
  90. C     XMAT(2)    COEFFICIENT DE POISSON
  91. C
  92. C     TRAC(1 A 2*LTRAC) COURBE DE TRACTION
  93. C     WORK(  " +1) ALFA1 POUR COQUES PLASTICITE GLOBALE
  94. C     WORK(  " +..) DONNEES POUR CRITERE POUTRES GLOBALES
  95. C
  96. C    MODELE ISOTROPE
  97. C            VAR0(1)  EPS*
  98. C
  99. C    MODELE CINEMATIQUE LINEAIRE
  100. C            VAR0(1)  EPS*
  101. C            VAR0(2) A VAR0(1+NSTRSS)  DEPLACEMENT DE LA SPHERE
  102. C
  103. C    MODELE CHABOCHE
  104. C            XMAT(5) ....      COEFFICIENTS A,C,...
  105. C            VAR0(1)  EPS*
  106. C            VAR0(2) A VAR0(1+NSTRSS) DEPLACEMENT DE LA SPHERE 1
  107. C            VAR0(2+NSTRSS) A VAR0(1+2*NSTRSS) "    "  "   "    2
  108. C
  109. C-----------------------------------------------------------------------
  110. C     20/09/2017 : modif SG critere de convergence trop serre
  111. C                  TEST=PETI*APHI0 + utilisation CCREEL
  112. C                  voir aussi ecoin0.eso, syco12.eso
  113.  
  114.       IMPLICIT INTEGER(I-N)
  115.       IMPLICIT REAL*8(A-H,O-Z)
  116.  
  117. -INC PPARAM
  118. -INC CCOPTIO
  119. -INC CCREEL
  120. -INC DECHE
  121.  
  122. -INC TECOU
  123.  
  124.       SEGMENT WRK2
  125.         REAL*8 TRAC(LTRAC2)
  126.       ENDSEGMENT
  127.  
  128.       SEGMENT WRK3
  129.         REAL*8 WORK(LW),WORK2(LW2)
  130.       ENDSEGMENT
  131.  
  132.       DIMENSION SIG(130),EPS(130)
  133.       DIMENSION S(8),SX(8),DS(8),DSIG(8),SPHER(8),SPHER1(8),SPHER2(8)
  134.       DIMENSION DSPHER1(8),DSPHER2(8),DEPSE(8),DEPSP(8),DDEPSE(8)
  135.       DIMENSION F(8),WF1(8),WF2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8)
  136.       DIMENSION CRIGI(12)
  137.  
  138.       jNPLAS = INPLAS
  139.       MFRl   = iecou.MFR1
  140.       nstrsl = iecou.NSTRSS
  141.       ncara = commat(/2)
  142.       if(ib.eq.1.and.igau.eq.1) then
  143.         do iaca = 1,ncara
  144.          if(commat(iaca).eq.'LIMP') iecou.JELEM = iaca
  145.         enddo
  146.       endif
  147.       icow21 = iecou.JELEM
  148.       if (icow21.gt.0) xlimp = valma0(icow21)
  149.  
  150.       do jj = 1,8
  151.         sx(jj) = 0.d0
  152.       enddo
  153.  
  154. C---------COQUES AVEC CT------------------------------------------------
  155. C         ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES
  156.       IF (MFRl.EQ.9) THEN
  157.         NCONT=6
  158.       ELSE
  159.         NCONT=iecou.NSTRSS
  160.       ENDIF
  161.  
  162.       itracb=0
  163.  
  164. C-----------------------------------------------------------------------
  165.  2222 continue
  166.       KERRE=0
  167.  
  168.       DO I=1,nstrsl
  169.         S(I)=0.D0
  170.         SPHER(I)=0.D0
  171.         SPHER1(I)=0.D0
  172.         SPHER2(I)=0.D0
  173.       ENDDO
  174.       YUNG=XMAT(1)
  175.       XNU=XMAT(2)
  176.  
  177. C---------CARACTERISTIQUES GEOMETRIQUES---------------------------------
  178. C
  179. C         COQUES
  180. C
  181.       IF (MFRl.EQ.3.OR.MFRl.EQ.9) THEN
  182.         EP1   = WORK(1)
  183.         ALFA1 = WORK(2)**2
  184.       ELSE IF (MFRl.EQ.5) THEN
  185.         EP1   = WORK(1)
  186.         ALFA1 = 1.D0
  187.       ELSE
  188.         EP1   = 1.D0
  189.         ALFA1 = 1.D0
  190.       ENDIF
  191.  
  192. C---------COEFFICIENTS CHABOCHE-----------------------------------------
  193.  
  194.       IF(jNPLAS.EQ.7) THEN
  195.         A1=XMAT(5)
  196.         C1=XMAT(6)
  197.         R0=XMAT(7)
  198.         PSI=XMAT(8)
  199.         OME=XMAT(9)
  200.         RM=XMAT(10)
  201.         B=XMAT(11)
  202.         A2=0.D0
  203.         C2=0.D0
  204.       ELSE IF(jNPLAS.EQ.12) THEN
  205.         A1=XMAT(5)
  206.         C1=XMAT(6)
  207.         A2=XMAT(7)
  208.         C2=XMAT(8)
  209.         R0=XMAT(9)
  210.         PSI=XMAT(10)
  211.         OME=XMAT(11)
  212.         RM=XMAT(12)
  213.         B=XMAT(13)
  214.       ELSE IF(jNPLAS.EQ.4) THEN
  215.         SLI0 = XMAT(5)
  216.         HECL = XMAT(6)
  217. c*dbg        if (ib.eq.1 .and. igau.eq.1) then
  218. c*dbg          write(ioimp,*) 'Ecrouissage',jNPLAS,SLI0,HECL
  219. c*dbg        endif
  220.       ELSE
  221.         DO I=1,LTRAC2
  222.           SIG(I)=TRAC(2*I-1)
  223.           EPS(I)=TRAC(2*I)
  224.         ENDDO
  225.       ENDIF
  226.  
  227.       EPST=VAR0(1)
  228. C---------ECROUISSAGE CINEMATIQUE---------------------------------------
  229.  
  230.       IF(jNPLAS.EQ.4.OR.jNPLAS.EQ.7.OR.jNPLAS.EQ.12) THEN
  231.         DO I=1,nstrsl
  232.           SPHER1(I)=VAR0(I+1)
  233.         ENDDO
  234.         IF(jNPLAS.EQ.12) THEN
  235.           DO I=1,nstrsl
  236.             SPHER2(I)=VAR0(nstrsl+1+I)
  237.           ENDDO
  238.         ENDIF
  239.         DO I=1,nstrsl
  240.           SPHER(I)=SPHER1(I)+SPHER2(I)
  241.         ENDDO
  242. c*dbg        if (ib.eq.1 .and.igau.eq.1) then
  243. c*dbg          write(ioimp,*) '  SPHERi',nstrsl
  244. c*dbg          write(ioimp,*) '  ',(spher1(iof),iof=1,nstrsl)
  245. c*dbg          write(ioimp,*) '  ',(spher(iof),iof=1,nstrsl)
  246. c*dbg        endif
  247. C-----------------------------------------------------------------------
  248. C         PREDICTEUR ELASTIQUE
  249. C         S     : CONTRAINTE
  250. C         SPHER : VARIABLE D'ECROUISSAGE CINEMATIQUE
  251. C         SX = S - X
  252. C-----------------------------------------------------------------------
  253.  
  254. *  en elastique non lineaire on annule les contraintes initiales
  255. * mais on cumule les epsilons elastiques
  256.       ELSE IF(jNPLAS.EQ.87) THEN
  257.         EPST=0.D0
  258.         DO I=1,nstrsl
  259.           SIG0(I)=0.D0
  260.           DEPST(I)=DEPST(I) + VAR0(I+1)
  261.         ENDDO
  262.       ENDIF
  263.       CALL CALSIG(DEPST,DDAUX,nstrsl,CMATE,VALMAT,VALCAR,
  264.      &            N2EL,N2PTEL,MFRl,IFOUR2,IB,IGAU,EPAIST,
  265.      &            NBPGAU,MELE,NPINT,iecou.NBGMAT,iecou.NELMAT,
  266.      &            SECT,LHOOK,TXR,
  267.      &            XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
  268.  
  269.       IF(IRTD.NE.1) THEN
  270.          KERRE=69
  271.          GOTO 1000
  272.       ENDIF
  273.  
  274.       IF ((mfr.eq.1.or.mfr.eq.31).and. IFOUR2.eq.-2) then
  275. * en cas de contraintes planes on repasse en 3D
  276.         do ju=1,6
  277.            do iu=1,6
  278.              DDA(iu,ju)=0.d0
  279.            enddo
  280.         enddo
  281.         cte_cp = xnu / (1.d0 - xnu)
  282.         aux= YUNG / ((1.d0+xnu)*(1.d0-2.d0*xnu))
  283.         aux1= aux * xnu
  284.         aux2= aux * (1.d0-xnu)
  285.          gege = Yung / (2.d0 *(1.d0 +xnu))
  286.           DDA(1,1)=AUX2
  287.           DDA(2,1)=AUX1
  288.           DDA(1,2)=AUX1
  289.           DDA(2,2)=AUX2
  290.           DDA(3,3)=aux2
  291.           DDA(1,3)=aux1
  292.           DDA(2,3)=aux1
  293.           DDA(3,1)=aux1
  294.           DDA(3,2)=aux1
  295.           DDA(4,4)=gege
  296.           DDA(5,5)=gege
  297.           DDA(6,6)=GEGE
  298.  
  299.       ELSE IF ((MFR.EQ.3.AND.NPINT.EQ.0).OR.MFR.EQ.9) THEN
  300.         AUX=YUNG/(1.D0-XNU*XNU)
  301.         AUX1=AUX*XNU
  302.         DO J=1,NCONT
  303.           DO I=1,NCONT
  304.             DDAUX(I,J)=0.D0
  305.           ENDDO
  306.         ENDDO
  307. C
  308. C     CAS TRIDIMENSIONNEL  ET FOURIER
  309. C
  310.         IF(IFOUR2.EQ.2.OR.IFOUR2.EQ.1) THEN
  311.  
  312.           GEGE=0.5D0*YUNG/(1.D0+XNU)
  313.           DDAUX(1,1)=AUX
  314.           DDAUX(2,1)=AUX1
  315.           DDAUX(1,2)=AUX1
  316.           DDAUX(2,2)=AUX
  317.           DDAUX(3,3)=GEGE
  318.           DDAUX(4,4)=AUX
  319.           DDAUX(5,4)=AUX1
  320.           DDAUX(4,5)=AUX1
  321.           DDAUX(5,5)=AUX
  322.           DDAUX(6,6)=GEGE
  323. C
  324. C     CAS AXISYMETRIQUE ET DEFORMATIONS PLANES
  325. C
  326.         ELSE IF(IFOUR2.EQ.0.OR.IFOUR2.EQ.-1.OR.IFOUR2.EQ.-3) THEN
  327.  
  328.           DDAUX(1,1)=AUX
  329.           DDAUX(2,1)=AUX1
  330.           DDAUX(1,2)=AUX1
  331.           DDAUX(2,2)=AUX
  332.           DDAUX(3,3)=AUX
  333.           DDAUX(4,3)=AUX1
  334.           DDAUX(3,4)=AUX1
  335.           DDAUX(4,4)=AUX
  336. C
  337. C     CAS CONTRAINTES PLANES
  338. C
  339.         ELSE IF(IFOUR2.EQ.-2) THEN
  340.           DDAUX(1,1)=YUNG
  341.           DDAUX(3,3)=YUNG
  342.         ENDIF
  343.       ENDIF
  344. *
  345.       DO I=1,nstrsl
  346.         S(I)=SIG0(I)+DSIGT(I)
  347.         SIGB(I)=S(I)
  348.         SX(I)=S(I)-SPHER(I)
  349.       ENDDO
  350. C---------CAS DES POUTRES-----------------------------------------------
  351.  
  352.       IF(MFRl.EQ.7) THEN
  353. **      write(6,*) 'work',(work(ic),ic=1,LW)
  354.         DIV(1)=1.D0/WORK(4)
  355.         DIV(2)=1.D0
  356.         DIV(3)=1.D0
  357.         DIV(4)=WORK(10)/WORK(1)
  358.         DIV(5)=WORK(11)/WORK(2)
  359.         DIV(6)=WORK(12)/WORK(3)
  360.         IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(WORK(1)*WORK(4))
  361.         IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(WORK(2)*WORK(4))
  362.         IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(WORK(3)*WORK(4))
  363.         DO I=1,NCONT
  364.           S(I)  = S(I) *DIV(I)
  365.           SX(I) = SX(I)*DIV(I)
  366.           SPHER(I)  = SPHER(I) *DIV(I)
  367.           SPHER1(I) = SPHER1(I)*DIV(I)
  368.           SPHER2(I) = SPHER2(I)*DIV(I)
  369.         END DO
  370.       ENDIF
  371. **         if (ib.eq.1 .AND.igau.eq.1) then
  372. **           write(ioimp,*) '  SIGi',NSTRSS
  373. **           write(ioimp,*) '  S',(s(iof),iof=1,NSTRSS)
  374. **           write(ioimp,*) '  SX',(sx(iof),iof=1,NSTRSS)
  375. **         endif
  376.  
  377. C-----------------------------------------------------------------------
  378. C         CALCUL DE LA LIMITE ELASTIQUE SI
  379. C-----------------------------------------------------------------------
  380.       IF (jNPLAS.EQ.1) THEN
  381.         SI=TRAC(1)
  382. C* Modele a ECROUISSAGE CINEMATIQUE LINEAIRE : la limite d'elasticite
  383. C* est la limite initiale donnee par SIGY
  384.       ELSE IF (jNPLAS.EQ.4) THEN
  385.         SI = SLI0
  386.       ELSE IF (jNPLAS.EQ.5.OR.jNPLAS.EQ.87) THEN
  387.         CALL TRACTI(SI,EPST,SIG,EPS,LTRAC2,2,izz)
  388.         IF (izz.EQ.1) THEN
  389.           KERRE=75
  390.           GOTO 1000
  391.         ENDIF
  392. C* Modele de CHABOCHE (prise en compte ecrouissage)
  393.       ELSE IF (jNPLAS.EQ.7 .OR. jNPLAS.EQ.12) THEN
  394.         RMmRR = (RM - R0) * EXP(-B*EPST)
  395.         SI = RM - RMmRR
  396.       ENDIF
  397. c*dbg      if (ib.eq.1 .AND.igau.eq.1) then
  398. c*dbg        write(ioimp,*) 'Limite d elasticite SI=',SI,TRAC(1),EPST
  399. c*dbg      endif
  400.  
  401. **
  402. * kich : application pression limite trace sigma
  403. **
  404.        yxsxii = 0.d0
  405.        if (icow21.gt.0) then
  406.           ytr = (sx(1) + sx(2) + sx(3))/3.D0
  407.           if (ytr.gt.xlimp) then
  408.             yxsxii = ytr - xlimp
  409.           else
  410.             r_z = ytr + xlimp
  411.             if (r_z.lt.0.D0) yxsxii = r_z
  412.           endif
  413.           if (yxsxii.ne.0.d0) then
  414.             dsigii = dsigt(1) + dsigt(2) + dsigt(3)
  415.             if (dsigii.ne.0.D0) then
  416.              r_z = 3.D0*yxsxii/dsigii
  417.              do jj = 1,3
  418.                s(jj) = s(jj) - (dsigt(jj)*r_z)
  419.                sx(jj) = sx(jj) - (dsigt(jj)*r_z)
  420.              enddo
  421.             else
  422.              do jj = 1,3
  423.                s(jj) = s(jj) - yxsxii
  424.                sx(jj) = sx(jj) - yxsxii
  425.              enddo
  426.             endif
  427.           endif
  428.        endif
  429.  
  430. C-----------------------------------------------------------------------
  431. C         CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
  432. C-----------------------------------------------------------------------
  433. *  attention en contraintes planes on se declare en 3D
  434. * rien besoin de faire dans vonmis0 car IFOUR2 n'est pas utilise
  435. * et les contraintes sont dimensionnees a 6
  436.       SEQ=VONMIS0(SX,nstrsl,MFRl,IFOUR2,EP1,ALFA1)
  437. **           if (ib.eq.1 .AND.igau.eq.1) then
  438. **           write(ioimp,*) 'VONMISES',SEQ,MFRl,IFOUR2,EP1,ALFA1
  439. **           endif
  440.  
  441. C-----------------------------------------------------------------------
  442. C         LE CRITERE EST-IL VERIFIE
  443. C-----------------------------------------------------------------------
  444.       IF (SEQ.LE.SI) THEN
  445.         ITRY = 1
  446.       ELSE
  447.         PETI=1.1D0*0.5D0*PRECIS*SEQ
  448.         CALL EPSPRE(SEQ,SI,PETI,ITRY)
  449.       ENDIF
  450.       IF (ITRY.EQ.1) THEN
  451. * rien a faire on n'a pas plastifie
  452. **        if (ib.eq.1 .and. igau.eq.1) then
  453. **          write(ioimp,*) 'pas de plastification'
  454. **        endif
  455.         IF (MFRl.EQ.7) THEN
  456.          DO I=1,NCONT
  457.             S(I)=S(I)/DIV(I)
  458.           ENDDO
  459.         ENDIF
  460.         DO I=1,NCONT
  461.           SIGF(I)=S(I)
  462.           DEFP(I)=0.D0
  463.         ENDDO
  464.         IF(MFRl.EQ.9) THEN
  465.           DEFP(7)=0.D0
  466.           DEFP(8)=0.D0
  467.           SIGF(7)=S(7)
  468.           SIGF(8)=S(8)
  469.         ENDIF
  470. **             if (ib.eq.1 .and. igau.eq.1) then
  471. **               write(ioimp,*) 'SIGF',(SIGF(iof),iof=1,nstrsl)
  472. **             endif
  473.  
  474.         VARF(1)=VAR0(1)
  475.         IF(jNPLAS.EQ.4.OR.jNPLAS.EQ.7) THEN
  476.           DO I=1,nstrsl
  477.             VARF(I+1)=VAR0(I+1)
  478.           ENDDO
  479.         ELSE IF(jNPLAS.EQ.12) THEN
  480.           DO I=1,2*nstrsl
  481.             VARF(I+1)=VAR0(I+1)
  482.           ENDDO
  483.         ELSE IF(jNPLAS.EQ.87) THEN
  484.           DO I=1,nstrsl
  485.             VARF(I+1)=DEPST(I)
  486.           ENDDO
  487.         ENDIF
  488. c*dbg        if (ib.eq.1 .and. igau.eq.1) then
  489. c*dbg          write(ioimp,*) 'VARF',EPST,'---',(VARF(iof),iof=2,1+nstrsl)
  490. c*dbg        endif
  491.         RETURN
  492.       ENDIF
  493.  
  494. C-----------------------------------------------------------------------
  495. C         ON A PLASTIFIE
  496. C-----------------------------------------------------------------------
  497.       NITER=0
  498.       PHI=SEQ-SI
  499. c*dbg       if (ib.eq.1 .and. igau.eq.1) then
  500. c*dbg         write(ioimp,*) 'niter,phi,si,seq,peti,precis=',
  501. c*dbg      $                 niter,phi,si,seq,peti,precis
  502. c*dbg       endif
  503.         cri0=  si * 1.D-8
  504.         PHI0=PHI
  505.         SI0=SI
  506.         RR=0.D0
  507.  
  508.         DO I=1,NCONT
  509.           DSIG(I)   =0.D0
  510.           DEPSP(I)  =0.D0
  511.           DSPHER1(I)=0.D0
  512.           DSPHER2(I)=0.D0
  513.         ENDDO
  514.  
  515. C-----------------------------------------------------------------------
  516. C         DEBUT DE LA BOUCLE D'ITERATIONS INTERNES
  517. C-----------------------------------------------------------------------
  518.         sx1in=sx(1)
  519.         sx2in=sx(2)
  520.         sx3in=sx(3)
  521.         s1in=s(1)
  522.         s2in=s(2)
  523.  
  524.         iderin=0
  525.  10     CONTINUE
  526.         NITER=NITER+1
  527.  
  528.          phi= seq - si
  529.  
  530. C---------CALCUL DE WF1=DF/D(SIGMA)--------------------------------------
  531.  
  532. C---------ELEMENTS MASSIFS----------------------------------------------
  533.  
  534.         IF(MFRl.EQ.1.OR.MFRl.EQ.31) THEN
  535.  
  536.           F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
  537.           F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
  538.           F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
  539.           DO I=4,nstrsl
  540.             F(I)=SX(I)
  541.           ENDDO
  542.           DO I=1,3
  543.             WF1(I)=1.5D0*F(I)/SEQ
  544.             Z1(I)=WF1(I)
  545.           ENDDO
  546.           DO I=4,NCONT
  547.             WF1(I)=3.D0*F(I)/SEQ
  548.             Z1(I)=1.5D0*F(I)/SEQ
  549.           ENDDO
  550.  
  551. C---------COQUES MINCES-------------------------------------------------
  552.  
  553.         ELSE IF(MFRl.EQ.3.OR.MFRl.EQ.9) THEN
  554.  
  555.           IF(IFOUR2.GE.1) THEN
  556.             WF1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1)
  557.             WF1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1)
  558.             WF1(3)=3.D0*SX(3)/(SEQ*EP1)
  559.             WF1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(5))/(SEQ*EP1*EP1)
  560.             WF1(5)=3.D0*WORK(2)*(2.D0*SX(5)-SX(4))/(SEQ*EP1*EP1)
  561.             WF1(6)=18.D0*WORK(2)*SX(6)/(SEQ*EP1*EP1)
  562.           ELSE
  563.             WF1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1)
  564.             WF1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1)
  565.             WF1(3)=3.D0*WORK(2)*(2.D0*SX(3)-SX(4))/(SEQ*EP1*EP1)
  566.             WF1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(3))/(SEQ*EP1*EP1)
  567.           ENDIF
  568.  
  569. C---------COQUES EPAISSES-----------------------------------------------
  570.         ELSE IF(MFRl.EQ.5) THEN
  571.           F(1)=(2.D0*SX(1)-SX(2))/3.D0
  572.           F(2)=(2.D0*SX(2)-SX(1))/3.D0
  573.           DO I=3,5
  574.             F(I)=SX(I)
  575.           ENDDO
  576.           DO I=1,2
  577.             WF1(I)=1.5D0*F(I)/SEQ
  578.             Z1(I)=WF1(I)
  579.           ENDDO
  580.           DO I=3,5
  581.             WF1(I)=3.D0*F(I)/SEQ
  582.             Z1(I)=1.5D0*F(I)/SEQ
  583.           ENDDO
  584.  
  585. C---------POUTRES-------------------------------------------------------
  586.  
  587.         ELSE IF(MFRl.EQ.7) THEN
  588.  
  589.           DO J=1,NCONT
  590.             DO I=1,NCONT
  591.               DDA(I,J)=0.D0
  592.             ENDDO
  593.           ENDDO
  594.           DDA(1,1)=YUNG
  595.           DDA(4,4)=0.5D0*YUNG/(1.D0+XNU)
  596.           DDA(5,5)=YUNG
  597.           DDA(6,6)=YUNG
  598.           WF1(1)=SX(1)/SEQ
  599.           WF1(2)=0.D0
  600.           WF1(3)=0.D0
  601.           WF1(4)=SX(4)/SEQ
  602.           WF1(5)=SX(5)/SEQ
  603.           WF1(6)=SX(6)/SEQ
  604.         ENDIF
  605.  
  606.         DO I=1,NCONT
  607.           WF2(I)=0.D0
  608.         ENDDO
  609.  
  610.         IF(MFRl.EQ.7) THEN
  611.           DO J=1,NCONT
  612.             XFLO1=WF1(J)
  613.             DO I=1,NCONT
  614.               WF2(I)=WF2(I)+XFLO1*DDA(I,J)
  615.             ENDDO
  616.           ENDDO
  617.  
  618.         ELSE
  619.           IF((mfr.eq.1.or.mfr.eq.31).and. IFOUR2.eq.-2) then
  620.             DO J=1,NCONT
  621.               XFLO1=WF1(J)
  622.               DO I=1,NCONT
  623.                 WF2(I)=WF2(I)+XFLO1*DDA(I,J)
  624.               ENDDO
  625.             ENDDO
  626.           else
  627.             DO J=1,NCONT
  628.               XFLO1=WF1(J)
  629.               DO I=1,NCONT
  630.                 WF2(I)=WF2(I)+XFLO1*DDAUX(I,J)
  631.               ENDDO
  632.             ENDDO
  633.           endif
  634.         ENDIF
  635.         COEF=0.D0
  636.         DO I=1,NCONT
  637.           COEF=COEF+WF1(I)*WF2(I)
  638.         ENDDO
  639.  
  640. C-----------------------------------------------------------------------
  641. C     PLASTICITE PARFAITE, ECROUISSAGE ISOTROPE ET CINEMATIQUE ZIEGLER
  642. C-----------------------------------------------------------------------
  643.  
  644.         IF (jNPLAS.EQ.1.OR.jNPLAS.EQ.4.OR.jNPLAS.EQ.5
  645.      $     .OR.jNPLAS.EQ.87) THEN
  646.  
  647.           IF (jNPLAS.EQ.4) THEN
  648.             RP=0.D0
  649.             C=HECL
  650.           ELSE
  651.           CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC2,1,izz)
  652.  
  653.           IF (izz.EQ.1) THEN
  654.             KERRE=75
  655.             GOTO 1000
  656.           ENDIF
  657.  
  658.           IF (jNPLAS.EQ.1) THEN
  659.             RP=0.D0
  660.             C=0.D0
  661.           ELSE IF(jNPLAS.EQ.5.OR.jNPLAS.EQ.87) THEN
  662.             RP=PENTE
  663.             C=0.D0
  664.           ENDIF
  665.           ENDIF
  666. c*dbg            if (ib.eq.1 .and. igau.eq.1) then
  667. c*dbg              write(ioimp,*) 'RP,C=',RP,C,EPST
  668. c*dbg            endif
  669.  
  670.           DENOM=COEF+C+RP
  671.           DELTA=PHI/DENOM
  672.           DMU=C*DELTA/SEQ
  673.  
  674.           DO I=1,NCONT
  675.             DSIG(I)=-DELTA*WF2(I)
  676.             DSPHER1(I)=DMU*SX(I)
  677.           ENDDO
  678.  
  679. * Cas des contraintes planes en massif
  680.           if ((mfr.eq.1.or.mfr.eq.31).and.IFOUR2.eq.-2) then
  681.  
  682.              bb= abs(dsig(3)+ s(3) )
  683.              r_z = dsig(3) * cte_cp
  684.              sx(3)= sx3in - dsig(3)
  685.              sx(1)= sx1in - r_z
  686.              sx(2)= sx2in - r_z
  687.              SEQ=VONMIS0(SX,nstrsl,MFRl,IFOUR2,EP1,ALFA1)
  688.               s(3)=  - dsig(3)
  689.               s(1)= s1in - r_z
  690.               s(2)= s2in - r_z
  691.               if (bb.gt.cri0) then
  692.                 if (iderin.eq.0) then
  693.                   niter=niter - 1
  694.                 endif
  695.                iderin=iderin+1
  696.                if (iderin.gt.50) then
  697.                  write(ioimp,*) '   probleme dans iterations internes'
  698.                           KERRE=2
  699.                           GO TO 1000
  700.                endif
  701.                go to 10
  702.               endif
  703.               DMU=C*DELTA/SEQ
  704.               DO I=1,NCONT
  705.                 DSPHER1(I)=DMU*SX(I)
  706.               ENDDO
  707.           endif
  708.  
  709.           iderin=0
  710.           DP=DELTA
  711.           DR=RP*DP
  712.  
  713.         ELSE
  714.  
  715. C-----------------------------------------------------------------------
  716. C         MODELE DE CHABOCHE
  717. C-----------------------------------------------------------------------
  718.  
  719. C---------UNIQUEMENT POUR LES ELEMENTS MASSIFS--------------------------
  720.  
  721.            XPRO1=0.D0
  722.            XPRO2=0.D0
  723.            DO I=1,NCONT
  724.             XPRO1=XPRO1+WF1(I)*SPHER1(I)
  725.             XPRO2=XPRO2+WF1(I)*SPHER2(I)
  726.            ENDDO
  727.  
  728.            FIP=1.D0+(PSI-1.D0)*EXP(-OME*EPST)
  729.  
  730.            DENOM=COEF+(A1*C1+A2*C2)*FIP-C1*XPRO1-C2*XPRO2+B*RMmRR
  731.            DELTA=PHI/DENOM
  732.  
  733.            DO I=1,NCONT
  734.             DSIG(I)=-DELTA*WF2(I)
  735.             DSPHER1(I)=(2.D0*A1*FIP*Z1(I)/3.D0-SPHER1(I))*C1*DELTA
  736.             DSPHER2(I)=(2.D0*A2*FIP*Z1(I)/3.D0-SPHER2(I))*C2*DELTA
  737.            ENDDO
  738.  
  739.            DR=B* RMmRR *DELTA
  740.            DP=DELTA
  741.         ENDIF
  742.  
  743.         RR=RR+DR
  744.         EPST=EPST+DP
  745.         EPST=MAX(0.D0,EPST)
  746. c*dbg        if (ib.eq.1 .and. igau.eq.1) then
  747. c*dbg          write(ioimp,*) 'MAJ=',EPST,DP,RR,DR
  748. c*dbg        endif
  749.  
  750.         IF (MFRl.EQ.3.OR.MFRl.EQ.9) THEN
  751.            r_z = EP1*EP1/(6.D0*WORK(2))
  752.            IF (IFOUR2.GE.1) THEN
  753.              DSIG(1)=DSIG(1)*EP1
  754.              DSIG(2)=DSIG(2)*EP1
  755.              DSIG(3)=DSIG(3)*EP1
  756.              DSIG(4)=DSIG(4)*r_z
  757.              DSIG(5)=DSIG(5)*r_z
  758.              DSIG(6)=DSIG(6)*r_z
  759.            ELSE
  760.              DSIG(1)=DSIG(1)*EP1
  761.              DSIG(2)=DSIG(2)*EP1
  762.              DSIG(3)=DSIG(3)*r_z
  763.              DSIG(4)=DSIG(4)*r_z
  764.            ENDIF
  765.         ENDIF
  766. C mise a jour des contraintes
  767.           DO I=1,NCONT
  768.             S(I)=S(I)+DSIG(I)
  769.             SPHER1(I)=SPHER1(I)+DSPHER1(I)
  770.             SPHER2(I)=SPHER2(I)+DSPHER2(I)
  771.             SPHER(I)=SPHER1(I)+SPHER2(I)
  772.             SX(I)=S(I)-SPHER(I)
  773.           ENDDO
  774.           if(IFOUR2.eq.-2.and.(mfr.eq.1.or.mfr.eq.31)) then
  775.             s(3)=0.d0
  776.           endif
  777.  
  778.         if (icow21.gt.0) then
  779.           yxsxii = 0.D0
  780. * kich borne trace sigma
  781.           ytr = (sx(1) + sx(2) + sx(3))/3.D0
  782.           if (ytr.gt.xlimp) then
  783.             yxsxii = ytr - xlimp
  784.           else
  785.             r_z = ytr + xlimp
  786.             if (r_z.lt.0.D0) yxsxii = r_z
  787.           endif
  788.           dsigii = dsigt(1) + dsigt(2) + dsigt(3)
  789.           if (dsigii.ne.0.D0) then
  790.             r_z = 3.D0*yxsxii/dsigii
  791.             do jj = 1,3
  792.               s(jj) = s(jj) - (dsigt(jj)*r_z)
  793.               sx(jj) = sx(jj) - (dsigt(jj)*r_z)
  794.             enddo
  795.           else
  796.             do jj = 1,3
  797.               s(jj) = s(jj) - yxsxii
  798.               sx(jj) = sx(jj) - yxsxii
  799.             enddo
  800.           endif
  801.         endif
  802.  
  803.         SEQ=VONMIS0(SX,nstrsl,MFRl,IFOUR2,EP1,ALFA1)
  804.  
  805. C---------CONTRAINTES PLANES--------------------------------------------
  806.  
  807.         IF(IFOUR2.EQ.-2) THEN
  808.  
  809.           IF(MFRl.EQ.1.OR.MFRl.EQ.31) THEN
  810.             F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
  811.             F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
  812.             F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
  813.             DO I=4,nstrsl
  814.               F(I)=SX(I)
  815.             ENDDO
  816.             DO I=1,3
  817.               WF1(I)=1.5D0*F(I)/SEQ
  818.             ENDDO
  819.             DO I=4,nstrsl
  820.               WF1(I)=3.D0*F(I)/SEQ
  821.             ENDDO
  822.  
  823.           ELSE IF(MFRl.EQ.3.OR.MFRl.EQ.9) THEN
  824.             AUX=EP1*EP1*EP1*EP1
  825.             WF1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1*EP1)
  826.             WF1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1*EP1)
  827.             WF1(3)=18.D0*ALFA1*(2.D0*SX(3)-SX(4))/(SEQ*AUX)
  828.             WF1(4)=18.D0*ALFA1*(2.D0*SX(4)-SX(3))/(SEQ*AUX)
  829.           ENDIF
  830.  
  831.           DO I=1,nstrsl
  832.             DEPSP(I)=DEPSP(I)+DELTA*WF1(I)
  833.           ENDDO
  834.         ENDIF
  835.  
  836. C-----------------------------------------------------------------------
  837. C         TEST
  838. C         CALCUL DE LA NOUVELLE VALEUR DE PHI
  839. C-----------------------------------------------------------------------
  840.         IF (jNPLAS.EQ.5.OR.jNPLAS.EQ.87) THEN
  841.            CALL TRACTI(SI,EPST,SIG,EPS,LTRAC2,2,izz)
  842. C* Modele de CHABOCHE (prise en compte ecrouissage)
  843.         ELSE IF (jNPLAS.EQ.7 .OR. jNPLAS.EQ.12) THEN
  844.           RMmRR = (RM - R0) * EXP(-B*EPST)
  845.           SI = RM - RMmRR
  846.         ELSE IF (jNPLAS.EQ.4) THEN
  847.            SI=SLI0
  848.         ELSE
  849.            SI=RR+SI0
  850.         ENDIF
  851.         PHI=SEQ-SI
  852.  
  853.         PETI=1.D-8
  854.         APHI=ABS(PHI)
  855.         APHI0=ABS(PHI0)
  856.         TEST=max(PETI*APHI0,XZPREC*100.D0*SEQ)
  857. *sg        TEST=PETI*APHI0
  858. *sg          write(ioimp,*) 'niter,phi,phi0,si,seq,rmmrr,test=',
  859. *sg     $       niter,phi,phi0,si,seq,rmmrr,test
  860.         IF(NITER.GT.50) THEN
  861.           if(itracb.eq.0) then
  862.              itracb=1
  863.              go to 2222
  864.           endif
  865. C     write(ioimp,*) ' CCOIN0, APHI, TEST',APHI,TEST
  866. C     write(ioimp,*) ' CCOIN0, KERRE = 2, NITER =',NITER
  867.           KERRE=-2
  868.           call soucis(268)
  869. C         GO TO 1000
  870.         ENDIF
  871.         IF(APHI.LE.TEST.OR.KERRE.LT.0) THEN
  872.           IF (KERRE.LT.0) KERRE = 0
  873.  
  874.           IF (MFRl.EQ.7) THEN
  875.             DO I=1,NCONT
  876.               S(I)=S(I)/DIV(I)
  877.               SPHER1(I)=SPHER1(I)/DIV(I)
  878.               SPHER2(I)=SPHER2(I)/DIV(I)
  879.             ENDDO
  880.           ENDIF
  881.  
  882. C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES-------------------
  883.  
  884.           IF(IFOUR2.NE.-2) THEN
  885.             DO I=1,NCONT
  886.               DS(I)=S(I)-SIGB(I)
  887.             ENDDO
  888.             CALL EPSIG0(DS,DDEPSE,MFRl,IFOUR2,YUNG,XNU,WORK,nstrsl)
  889.             DO I=1,NCONT
  890.               DEPSE(I)=DEPST(I)+DDEPSE(I)
  891.               DEPSP(I)=DEPST(I)-DEPSE(I)
  892.             ENDDO
  893.           ENDIF
  894.  
  895.           DO I=1,nstrsl
  896.             SIGF(I)=S(I)
  897.             DEFP(I)=DEPSP(I)
  898.           ENDDO
  899. c*dbg          if (ib.eq.1 .and. igau.eq.1) then
  900. c*dbg            write(ioimp,*) 'SIGF',(SIGF(iof),iof=1,nstrsl)
  901. c*dbg          endif
  902.  
  903. C---------COQUES AVEC CISAILLEMENT TRANSVERSE---------------------------
  904.  
  905.           IF(MFRl.EQ.9) THEN
  906.             DEFP(7)=0.D0
  907.             DEFP(8)=0.D0
  908.             SIGF(7)=SIGB(7)
  909.             SIGF(8)=SIGB(8)
  910.           ENDIF
  911.  
  912.           VARF(1)=EPST
  913.           IF(jNPLAS.EQ.4.OR.jNPLAS.EQ.7.OR.jNPLAS.EQ.12) THEN
  914.             DO I=1,nstrsl
  915.               VARF(I+1)=SPHER1(I)
  916.             ENDDO
  917.             IF(jNPLAS.EQ.12) THEN
  918.               DO I=1,nstrsl
  919.                 VARF(nstrsl+1+I)=SPHER2(I)
  920.               ENDDO
  921.             ENDIF
  922.           ELSE IF(jNPLAS.EQ.87) THEN
  923.               DO I=1,nstrsl
  924.                 VARF(1+I)=DEPST(I)
  925.               ENDDO
  926.           ENDIF
  927.           KERRE=0
  928. c*dbg          if (ib.eq.1 .and. igau.eq.1) then
  929. c*dbg            write(ioimp,*) 'VARF',EPST,'---',(VARF(iof),iof=2,1+nstrsl)
  930. c*dbg          endif
  931.           RETURN
  932.  
  933.         ELSE
  934.           sx1in=sx(1)
  935.           sx2in=sx(2)
  936.           s1in=s(1)
  937.           s2in=s(2)
  938.           GOTO 10
  939.         ENDIF
  940.  
  941. 1000  CONTINUE
  942. C      RETURN
  943.       END
  944.  
  945.  
  946.  

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