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    PV090527  24/04/04    21:15:04     11875          
  2.  
  3.       SUBROUTINE CCOIN0(wrk52,wrk53,wrk54,wrk2,wrk3,
  4.      &  IB,IGAU,NBPGAU,LTRAC,IFOURB,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
  23. C                  NSTRSS      NOMBRE DE CONTRAINTES CA2000
  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                  IFOU  = OPTION DE CALCUL
  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
  39. C                  NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
  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) ALFAH 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.        SEGMENT IECOU
  123.        INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
  124.      1  icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
  125.      2  icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
  126.      3  icow25,icow26,icow27,icow28,icow29,icow30,icow31,
  127.      4  icow32,icow33,NSTRSS,MFR1,  NBGMAT,NELMAT,icow38,
  128.      5  icow39,icow40,icow41,icow42,icow43,icow44
  129.        INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
  130.      .         icow51,icow52,icow53,icow54,icow55,icow56
  131.      .         icow57,icow58
  132.       ENDSEGMENT
  133. *
  134.       SEGMENT WRK2
  135.         REAL*8 TRAC(LTRAC)
  136.       ENDSEGMENT
  137. *
  138.       SEGMENT WRK3
  139.         REAL*8 WORK(LW),WORK2(LW2)
  140.       ENDSEGMENT
  141.  
  142.       DIMENSION SIG(130),EPS(130)
  143.       DIMENSION S(8),SX(8),DS(8),DSIG(8),SPHER(8),SPHER1(8),SPHER2(8)
  144.       DIMENSION DSPHER1(8),DSPHER2(8),DEPSE(8),DEPSP(8),DDEPSE(8)
  145.       DIMENSION F(8),W1(8),W2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8)
  146.       DIMENSION CRIGI(12)
  147.  
  148. *      logical lvisu
  149. *      lvisu = ib.eq.125.and.(igau.eq.1.or.igau.eq.5)
  150.  
  151.       ncara = commat(/2)
  152.       if(ib.eq.1.and.igau.eq.1) then
  153.         do iaca = 1,ncara
  154.          if(commat(iaca).eq.'LIMP') icow21 = iaca
  155.         enddo
  156.       endif
  157.       if(icow21.gt.0) xlimp = valma0(icow21)
  158.       do jj = 1,8
  159.        sx(jj) = 0.d0
  160.       enddo
  161. C---------COQUES AVEC CT------------------------------------------------
  162. C         ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES
  163.       IF(MFR1.EQ.9) THEN
  164.         NCONT=6
  165.       ELSE
  166.         NCONT=NSTRSS
  167.       ENDIF
  168.       itracb=0
  169.  
  170. C-----------------------------------------------------------------------
  171.  2222 continue
  172.       KERRE=0
  173.  
  174.       DO I=1,NSTRSS
  175.         S(I)=0.D0
  176.         SPHER(I)=0.D0
  177.         SPHER1(I)=0.D0
  178.         SPHER2(I)=0.D0
  179.       ENDDO
  180.       YUNG=XMAT(1)
  181.       XNU=XMAT(2)
  182.  
  183. C---------CARACTERISTIQUES GEOMETRIQUES---------------------------------
  184. C
  185. C         COQUES
  186. C
  187.       ALFAH=1.D0
  188.       EP1  =1.D0
  189.       IF(MFR1.EQ.3.OR.MFR1.EQ.5.OR.MFR1.EQ.9) THEN
  190.         EP1=WORK(1)
  191.         IF(MFR1.NE.5) ALFAH=WORK(2)**2
  192.       ENDIF
  193.  
  194. C---------COEFFICIENTS CHABOCHE-----------------------------------------
  195.  
  196.       IF(INPLAS.EQ.7) THEN
  197.         A1=XMAT(5)
  198.         C1=XMAT(6)
  199.         R0=XMAT(7)
  200.         PSI=XMAT(8)
  201.         OME=XMAT(9)
  202.         RM=XMAT(10)
  203.         B=XMAT(11)
  204.         A2=0.D0
  205.         C2=0.D0
  206.       ELSE IF(INPLAS.EQ.12) THEN
  207.         A1=XMAT(5)
  208.         C1=XMAT(6)
  209.         A2=XMAT(7)
  210.         C2=XMAT(8)
  211.         R0=XMAT(9)
  212.         PSI=XMAT(10)
  213.         OME=XMAT(11)
  214.         RM=XMAT(12)
  215.         B=XMAT(13)
  216.       ELSE IF(INPLAS.EQ.4) THEN
  217.         SLI0 = XMAT(5)
  218.         HECL = XMAT(6)
  219. c*dbg        if (ib.eq.1 .and. igau.eq.1) then
  220. c*dbg          write(ioimp,*) 'Ecrouissage',INPLAS,SLI0,HECL
  221. c*dbg        endif
  222.       ELSE
  223.         DO I=1,LTRAC
  224.           SIG(I)=TRAC(2*I-1)
  225.           EPS(I)=TRAC(2*I)
  226.         ENDDO
  227.       ENDIF
  228.  
  229.       EPST=VAR0(1)
  230. C---------ECROUISSAGE CINEMATIQUE---------------------------------------
  231.  
  232.       IF(INPLAS.EQ.4.OR.INPLAS.EQ.7.OR.INPLAS.EQ.12) THEN
  233.         DO I=1,NSTRSS
  234.           SPHER1(I)=VAR0(I+1)
  235.         ENDDO
  236.         IF(INPLAS.EQ.12) THEN
  237.           DO I=1,NSTRSS
  238.             SPHER2(I)=VAR0(NSTRSS+1+I)
  239.           ENDDO
  240.         ENDIF
  241.         DO I=1,NSTRSS
  242.           SPHER(I)=SPHER1(I)+SPHER2(I)
  243.         ENDDO
  244. c*dbg        if (ib.eq.1 .and.igau.eq.1) then
  245. c*dbg          write(ioimp,*) '  SPHERi',NSTRSS
  246. c*dbg          write(ioimp,*) '  ',(spher1(iof),iof=1,NSTRSS)
  247. c*dbg          write(ioimp,*) '  ',(spher(iof),iof=1,NSTRSS)
  248. c*dbg        endif
  249. C-----------------------------------------------------------------------
  250. C         PREDICTEUR ELASTIQUE
  251. C         S     : CONTRAINTE
  252. C         SPHER : VARIABLE D'ECROUISSAGE CINEMATIQUE
  253. C         SX = S - X
  254. C-----------------------------------------------------------------------
  255.  
  256. *  en elastique non lineaire on annule les contraintes initiales
  257. * mais on cumule les epsilons elastiques
  258.       ELSE IF(INPLAS.EQ.87) THEN
  259.         EPST=0.D0
  260.         DO I=1,NSTRSS
  261.           SIG0(I)=0.D0
  262.           DEPST(I)=DEPST(I) + VAR0(I+1)
  263.         ENDDO
  264.       ENDIF
  265.       CALL CALSIG(DEPST,DDAUX,NSTRSS,CMATE,VALMAT,VALCAR,
  266.      &        N2EL,N2PTEL,MFR1,IFOURB,IB,IGAU,EPAIST,
  267.      &      NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
  268.      &    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
  269.  
  270.       IF(IRTD.NE.1) THEN
  271.          KERRE=69
  272.          GOTO 1000
  273.       ENDIF
  274.  
  275.       IF ((mfr.eq.1.or.mfr.eq.31).and. ifourb.eq.-2) then
  276. * en cas de contraintes planes on repasse en 3D
  277.         do ju=1,6
  278.            do iu=1,6
  279.              DDA(iu,ju)=0.d0
  280.            enddo
  281.         enddo
  282.         cte_cp = xnu / (1.d0 - xnu)
  283.         aux= YUNG / ((1.d0+xnu)*(1.d0-2.d0*xnu))
  284.         aux1= aux * xnu
  285.         aux2= aux * (1.d0-xnu)
  286.          gege = Yung / (2.d0 *(1.d0 +xnu))
  287.           DDA(1,1)=AUX2
  288.           DDA(2,1)=AUX1
  289.           DDA(1,2)=AUX1
  290.           DDA(2,2)=AUX2
  291.           DDA(3,3)=aux2
  292.           DDA(1,3)=aux1
  293.           DDA(2,3)=aux1
  294.           DDA(3,1)=aux1
  295.           DDA(3,2)=aux1
  296.           DDA(4,4)=gege
  297.           DDA(5,5)=gege
  298.           DDA(6,6)=GEGE
  299.  
  300.       ELSE IF ((MFR.EQ.3.AND.NPINT.EQ.0).OR.MFR.EQ.9) THEN
  301.         AUX=YUNG/(1.D0-XNU*XNU)
  302.         AUX1=AUX*XNU
  303.         DO J=1,NCONT
  304.           DO I=1,NCONT
  305.             DDAUX(I,J)=0.D0
  306.           ENDDO
  307.         ENDDO
  308. C
  309. C     CAS TRIDIMENSIONNEL  ET FOURIER
  310. C
  311.         IF(IFOURB.EQ.2.OR.IFOURB.EQ.1) THEN
  312.  
  313.           GEGE=0.5D0*YUNG/(1.D0+XNU)
  314.           DDAUX(1,1)=AUX
  315.           DDAUX(2,1)=AUX1
  316.           DDAUX(1,2)=AUX1
  317.           DDAUX(2,2)=AUX
  318.           DDAUX(3,3)=GEGE
  319.           DDAUX(4,4)=AUX
  320.           DDAUX(5,4)=AUX1
  321.           DDAUX(4,5)=AUX1
  322.           DDAUX(5,5)=AUX
  323.           DDAUX(6,6)=GEGE
  324. C
  325. C     CAS AXISYMETRIQUE ET DEFORMATIONS PLANES
  326. C
  327.         ELSE IF(IFOURB.EQ.0.OR.IFOURB.EQ.-1.OR.IFOURB.EQ.-3) THEN
  328.  
  329.           DDAUX(1,1)=AUX
  330.           DDAUX(2,1)=AUX1
  331.           DDAUX(1,2)=AUX1
  332.           DDAUX(2,2)=AUX
  333.           DDAUX(3,3)=AUX
  334.           DDAUX(4,3)=AUX1
  335.           DDAUX(3,4)=AUX1
  336.           DDAUX(4,4)=AUX
  337. C
  338. C     CAS CONTRAINTES PLANES
  339. C
  340.         ELSE IF(IFOURB.EQ.-2) THEN
  341.           DDAUX(1,1)=YUNG
  342.           DDAUX(3,3)=YUNG
  343.         ENDIF
  344.       ENDIF
  345. *
  346.       DO I=1,NSTRSS
  347.         S(I)=SIG0(I)+DSIGT(I)
  348.         SIGB(I)=S(I)
  349.         SX(I)=S(I)-SPHER(I)
  350.       ENDDO
  351. C---------CAS DES POUTRES-----------------------------------------------
  352.  
  353.       IF(MFR1.EQ.7) THEN
  354. **      write(6,*) 'work',(work(ic),ic=1,LW)
  355.         DIV(1)=1.D0/WORK(4)
  356.         DIV(2)=1.D0
  357.         DIV(3)=1.D0
  358.         DIV(4)=WORK(10)/WORK(1)
  359.         DIV(5)=WORK(11)/WORK(2)
  360.         DIV(6)=WORK(12)/WORK(3)
  361.         IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(WORK(1)*WORK(4))
  362.         IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(WORK(2)*WORK(4))
  363.         IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(WORK(3)*WORK(4))
  364.         DO I=1,NCONT
  365.           S(I)  = S(I) *DIV(I)
  366.           SX(I) = SX(I)*DIV(I)
  367.           SPHER(I)  = SPHER(I) *DIV(I)
  368.           SPHER1(I) = SPHER1(I)*DIV(I)
  369.           SPHER2(I) = SPHER2(I)*DIV(I)
  370.         END DO
  371.       ENDIF
  372.            if (ib.eq.1 .AND.igau.eq.1) then
  373. **           write(ioimp,*) '  SIGi',NSTRSS
  374. **           write(ioimp,*) '  S',(s(iof),iof=1,NSTRSS)
  375. **           write(ioimp,*) '  SX',(sx(iof),iof=1,NSTRSS)
  376.            endif
  377.  
  378. C-----------------------------------------------------------------------
  379. C         CALCUL DE LA LIMITE ELASTIQUE SI
  380. C-----------------------------------------------------------------------
  381.       IF (INPLAS.EQ.1) THEN
  382.         SI=TRAC(1)
  383. C* Modele a ECROUISSAGE CINEMATIQUE LINEAIRE : la limite d'elasticite 
  384. C* est la limite initiale donnee par SIGY
  385.       ELSE IF (INPLAS.EQ.4) THEN
  386.         SI = SLI0
  387.       ELSE IF (INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN
  388.         CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI)
  389.         IF(IBI.EQ.1) THEN
  390.           KERRE=75
  391.         GOTO 1000
  392.         ENDIF
  393. C* Modele de CHABOCHE (prise en compte ecrouissage)
  394.       ELSE IF (INPLAS.EQ.7 .OR. INPLAS.EQ.12) THEN
  395.         RMmRR = (RM - R0) * EXP(-B*EPST)
  396.         SI = RM - RMmRR
  397.       ENDIF
  398. c*dbg      if (ib.eq.1 .AND.igau.eq.1) then
  399. c*dbg        write(ioimp,*) 'Limite d elasticite SI=',SI,TRAC(1),EPST
  400. c*dbg      endif
  401.  
  402. **
  403. * kich : application pression limite trace sigma
  404. **
  405.        yxsxii = 0.d0
  406.        if (icow21.gt.0) then
  407.           ytr = sx(1) + sx(2) + sx(3)
  408.           ytr = ytr/3.D0
  409.           if (ytr.gt.xlimp) then
  410.             yxsxii = ytr - xlimp
  411.           else if((ytr + xlimp).lt.0) then
  412.             yxsxii = ytr + xlimp
  413.           endif
  414.           if (yxsxii.ne.0) then
  415.             dsigii = dsigt(1) + dsigt(2) + dsigt(3)
  416.             if (dsigii.ne.0.D0) then
  417.              do jj = 1,3
  418.                s(jj) = s(jj) - (dsigt(jj)/dsigii*3.D0*yxsxii)
  419.                sx(jj) = sx(jj) - (dsigt(jj)/dsigii*3.D0*yxsxii)
  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 ifourb n'est pas utilise
  435. * et les contraintes sont dimensionnees a 6
  436.       SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH)
  437.            if (ib.eq.1 .AND.igau.eq.1) then
  438. **           write(ioimp,*) 'VONMISES',SEQ,MFR1,IFOURB,EP1,ALFAH
  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 (MFR1.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(MFR1.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,NSTRSS)
  472.              endif
  473.  
  474.         VARF(1)=VAR0(1)
  475.         IF(INPLAS.EQ.4.OR.INPLAS.EQ.7) THEN
  476.           DO I=1,NSTRSS
  477.             VARF(I+1)=VAR0(I+1)
  478.           ENDDO
  479.         ELSE IF(INPLAS.EQ.12) THEN
  480.           DO I=1,2*NSTRSS
  481.             VARF(I+1)=VAR0(I+1)
  482.           ENDDO
  483.         ELSE IF(INPLAS.EQ.87) THEN
  484.           DO I=1,NSTRSS
  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+NSTRSS)
  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 W1=DF/D(SIGMA)--------------------------------------
  531.  
  532. C---------ELEMENTS MASSIFS----------------------------------------------
  533.  
  534.         IF(MFR1.EQ.1.OR.MFR1.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,NSTRSS
  540.             F(I)=SX(I)
  541.           ENDDO
  542.           DO I=1,3
  543.             W1(I)=1.5D0*F(I)/SEQ
  544.             Z1(I)=W1(I)
  545.           ENDDO
  546.           DO I=4,NCONT
  547.             W1(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(MFR1.EQ.3.OR.MFR1.EQ.9) THEN
  554.  
  555.           IF(IFOURB.GE.1) THEN
  556.             W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1)
  557.             W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1)
  558.             W1(3)=3.D0*SX(3)/(SEQ*EP1)
  559.             W1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(5))/(SEQ*EP1*EP1)
  560.             W1(5)=3.D0*WORK(2)*(2.D0*SX(5)-SX(4))/(SEQ*EP1*EP1)
  561.             W1(6)=18.D0*WORK(2)*SX(6)/(SEQ*EP1*EP1)
  562.           ELSE
  563.             W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1)
  564.             W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1)
  565.             W1(3)=3.D0*WORK(2)*(2.D0*SX(3)-SX(4))/(SEQ*EP1*EP1)
  566.             W1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(3))/(SEQ*EP1*EP1)
  567.           ENDIF
  568.  
  569.  
  570. C---------COQUES EPAISSES-----------------------------------------------
  571.         ELSE IF(MFR1.EQ.5) THEN
  572.           F(1)=(2.D0*SX(1)-SX(2))/3.D0
  573.           F(2)=(2.D0*SX(2)-SX(1))/3.D0
  574.           DO I=3,5
  575.             F(I)=SX(I)
  576.           ENDDO
  577.           DO I=1,2
  578.             W1(I)=1.5D0*F(I)/SEQ
  579.             Z1(I)=W1(I)
  580.           ENDDO
  581.           DO I=3,5
  582.             W1(I)=3.D0*F(I)/SEQ
  583.             Z1(I)=1.5D0*F(I)/SEQ
  584.           ENDDO
  585.  
  586. C---------POUTRES-------------------------------------------------------
  587.  
  588.         ELSE IF(MFR1.EQ.7) THEN
  589.  
  590.           DO J=1,NCONT
  591.             DO I=1,NCONT
  592.               DDA(I,J)=0.D0
  593.             ENDDO
  594.           ENDDO
  595.           DDA(1,1)=YUNG
  596.           DDA(4,4)=0.5D0*YUNG/(1.D0+XNU)
  597.           DDA(5,5)=YUNG
  598.           DDA(6,6)=YUNG
  599.           W1(1)=SX(1)/SEQ
  600.           W1(2)=0.D0
  601.           W1(3)=0.D0
  602.           W1(4)=SX(4)/SEQ
  603.           W1(5)=SX(5)/SEQ
  604.           W1(6)=SX(6)/SEQ
  605.         ENDIF
  606.  
  607.         DO I=1,NCONT
  608.           W2(I)=0.D0
  609.         ENDDO
  610.  
  611.         IF(MFR1.EQ.7) THEN
  612.           DO J=1,NCONT
  613.             XFLO1=W1(J)
  614.             DO I=1,NCONT
  615.               W2(I)=W2(I)+XFLO1*DDA(I,J)
  616.             ENDDO
  617.           ENDDO
  618.  
  619.         ELSE
  620.           IF((mfr.eq.1.or.mfr.eq.31).and. ifourb.eq.-2) then
  621.             DO J=1,NCONT
  622.               XFLO1=W1(J)
  623.               DO I=1,NCONT
  624.                 W2(I)=W2(I)+XFLO1*DDA(I,J)
  625.               ENDDO
  626.             ENDDO
  627.           else
  628.             DO J=1,NCONT
  629.               XFLO1=W1(J)
  630.               DO I=1,NCONT
  631.                 W2(I)=W2(I)+XFLO1*DDAUX(I,J)
  632.               ENDDO
  633.             ENDDO
  634.           endif
  635.         ENDIF
  636.         COEF=0.D0
  637.         DO I=1,NCONT
  638.           COEF=COEF+W1(I)*W2(I)
  639.         ENDDO
  640.  
  641. C-----------------------------------------------------------------------
  642. C     PLASTICITE PARFAITE, ECROUISSAGE ISOTROPE ET CINEMATIQUE ZIEGLER
  643. C-----------------------------------------------------------------------
  644.  
  645.         IF (INPLAS.EQ.1.OR.INPLAS.EQ.4.OR.INPLAS.EQ.5
  646.      $     .OR.INPLAS.EQ.87) THEN
  647.  
  648.           IF (INPLAS.EQ.4) THEN
  649.             RP=0.D0
  650.             C=HECL
  651.           ELSE
  652.           CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC,1,IBI)
  653.  
  654.           IF (IBI.EQ.1) THEN
  655.             KERRE=75
  656.             GOTO 1000
  657.           ENDIF
  658.  
  659.           IF (INPLAS.EQ.1) THEN
  660.             RP=0.D0
  661.             C=0.D0
  662.           ELSE IF(INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN
  663.             RP=PENTE
  664.             C=0.D0
  665.           ENDIF
  666.           ENDIF
  667. c*dbg            if (ib.eq.1 .and. igau.eq.1) then
  668. c*dbg              write(ioimp,*) 'RP,C=',RP,C,EPST
  669. c*dbg            endif
  670.  
  671.           DENOM=COEF+C+RP
  672.           DELTA=PHI/DENOM
  673.           DMU=C*DELTA/SEQ
  674.  
  675.           DO I=1,NCONT
  676.             DSIG(I)=-DELTA*W2(I)
  677.             DSPHER1(I)=DMU*SX(I)
  678.           ENDDO
  679.  
  680. * Cas des contraintes planes en massif
  681.           if ((mfr.eq.1.or.mfr.eq.31).and.ifourb.eq.-2) then
  682.  
  683.              bb= abs(dsig(3)+ s(3) )
  684.              r_z = dsig(3) * cte_cp
  685.              sx(3)= sx3in - dsig(3)
  686.              sx(1)= sx1in - r_z
  687.              sx(2)= sx2in - r_z
  688.              SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH)
  689.               s(3)=  - dsig(3)
  690.               s(1)= s1in - r_z
  691.               s(2)= s2in - r_z
  692.               if (bb.gt.cri0) then
  693.                 if (iderin.eq.0) then
  694.                   niter=niter - 1
  695.                 endif
  696.                iderin=iderin+1
  697.                if (iderin.gt.50) then
  698.                  write(ioimp,*) '   probleme dans iterations internes'
  699.                           KERRE=2
  700.                           GO TO 1000
  701.                endif
  702.                go to 10
  703.               endif
  704.               DMU=C*DELTA/SEQ
  705.               DO I=1,NCONT
  706.                 DSPHER1(I)=DMU*SX(I)
  707.               ENDDO
  708.           endif
  709.  
  710.           iderin=0
  711.           DP=DELTA
  712.           DR=RP*DP
  713.  
  714.         ELSE
  715.  
  716. C-----------------------------------------------------------------------
  717. C         MODELE DE CHABOCHE
  718. C-----------------------------------------------------------------------
  719.  
  720. C---------UNIQUEMENT POUR LES ELEMENTS MASSIFS--------------------------
  721.  
  722.            XPRO1=0.D0
  723.            XPRO2=0.D0
  724.            DO I=1,NCONT
  725.             XPRO1=XPRO1+W1(I)*SPHER1(I)
  726.             XPRO2=XPRO2+W1(I)*SPHER2(I)
  727.            ENDDO
  728.  
  729.            FIP=1.D0+(PSI-1.D0)*EXP(-OME*EPST)
  730.  
  731.            DENOM=COEF+(A1*C1+A2*C2)*FIP-C1*XPRO1-C2*XPRO2+B*RMmRR
  732.            DELTA=PHI/DENOM
  733.  
  734.            DO I=1,NCONT
  735.             DSIG(I)=-DELTA*W2(I)
  736.             DSPHER1(I)=(2.D0*A1*FIP*Z1(I)/3.D0-SPHER1(I))*C1*DELTA
  737.             DSPHER2(I)=(2.D0*A2*FIP*Z1(I)/3.D0-SPHER2(I))*C2*DELTA
  738.            ENDDO
  739.  
  740.            DR=B* RMmRR *DELTA
  741.            DP=DELTA
  742.         ENDIF
  743.  
  744.         RR=RR+DR
  745.         EPST=EPST+DP
  746.         EPST=MAX(0.D0,EPST)
  747. c*dbg        if (ib.eq.1 .and. igau.eq.1) then
  748. c*dbg          write(ioimp,*) 'MAJ=',EPST,DP,RR,DR
  749. c*dbg        endif
  750.  
  751.         IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN
  752.            IF(IFOURB.GE.1) THEN
  753.              DO I=1,3
  754.                DSIG(I)=DSIG(I)*EP1
  755.              ENDDO
  756.              DO I=4,6
  757.                DSIG(I)=EP1*EP1*DSIG(I)/(6.D0*WORK(2))
  758.              ENDDO
  759.            ELSE
  760.              DSIG(1)=DSIG(1)*EP1
  761.              DSIG(2)=DSIG(2)*EP1
  762.              DSIG(3)=EP1*EP1*DSIG(3)/(6.D0*WORK(2))
  763.              DSIG(4)=EP1*EP1*DSIG(4)/(6.D0*WORK(2))
  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(ifourb.eq.-2.and.(mfr.eq.1.or.mfr.eq.31)) then
  775.             s(3)=0.d0
  776.           endif
  777.  
  778.         yxsxii = 0.D0
  779.         if (icow21.gt.0) then
  780. * kich borne trace sigma
  781.           ytr = sx(1) + sx(2) + sx(3)
  782.           ytr = ytr/3.D0
  783.           if (ytr.gt.xlimp) then
  784.             yxsxii = ytr - xlimp
  785.           else if((ytr + xlimp).lt.0) then
  786.             yxsxii = ytr + xlimp
  787.           endif
  788.             dsigii = dsigt(1) + dsigt(2) + dsigt(3)
  789.             if (dsigii.ne.0.D0) then
  790.              do jj = 1,3
  791.                s(jj) = s(jj) - (dsigt(jj)/dsigii*3.D0*yxsxii)
  792.                sx(jj) = sx(jj) - (dsigt(jj)/dsigii*3.D0*yxsxii)
  793.              enddo
  794.             else
  795.              do jj = 1,3
  796.                s(jj) = s(jj) - yxsxii
  797.                sx(jj) = sx(jj) - yxsxii
  798.              enddo
  799.             endif
  800.         endif
  801.  
  802.         SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH)
  803.  
  804. C---------CONTRAINTES PLANES--------------------------------------------
  805.  
  806.         IF(IFOURB.EQ.-2) THEN
  807.  
  808.           IF(MFR1.EQ.1.OR.MFR1.EQ.31) THEN
  809.             F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
  810.             F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
  811.             F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
  812.             DO I=4,NSTRSS
  813.               F(I)=SX(I)
  814.             ENDDO
  815.             DO I=1,3
  816.               W1(I)=1.5D0*F(I)/SEQ
  817.             ENDDO
  818.             DO I=4,NSTRSS
  819.               W1(I)=3.D0*F(I)/SEQ
  820.             ENDDO
  821.  
  822.           ELSE IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN
  823.             AUX=EP1*EP1*EP1*EP1
  824.             W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1*EP1)
  825.             W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1*EP1)
  826.             W1(3)=18.D0*ALFAH*(2.D0*SX(3)-SX(4))/(SEQ*AUX)
  827.             W1(4)=18.D0*ALFAH*(2.D0*SX(4)-SX(3))/(SEQ*AUX)
  828.           ENDIF
  829.  
  830.           DO I=1,NSTRSS
  831.             DEPSP(I)=DEPSP(I)+DELTA*W1(I)
  832.           ENDDO
  833.         ENDIF
  834.  
  835. C-----------------------------------------------------------------------
  836. C         TEST
  837. C         CALCUL DE LA NOUVELLE VALEUR DE PHI
  838. C-----------------------------------------------------------------------
  839.         IF (INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN
  840.            CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI)
  841. C* Modele de CHABOCHE (prise en compte ecrouissage)
  842.         ELSE IF (INPLAS.EQ.7 .OR. INPLAS.EQ.12) THEN
  843.           RMmRR = (RM - R0) * EXP(-B*EPST)
  844.           SI = RM - RMmRR
  845.         ELSE IF (INPLAS.EQ.4) THEN
  846.            SI=SLI0
  847.         ELSE
  848.            SI=RR+SI0
  849.         ENDIF
  850.         PHI=SEQ-SI
  851.  
  852.         PETI=1.D-8
  853.         APHI=ABS(PHI)
  854.         APHI0=ABS(PHI0)
  855.         TEST=max(PETI*APHI0,XZPREC*100.D0*SEQ)
  856. *sg        TEST=PETI*APHI0
  857. *sg          write(ioimp,*) 'niter,phi,phi0,si,seq,rmmrr,test=',
  858. *sg     $       niter,phi,phi0,si,seq,rmmrr,test
  859.         IF(NITER.GT.50) THEN
  860.           if(itracb.eq.0) then
  861.              itracb=1
  862.              go to 2222
  863.           endif
  864. C     write(ioimp,*) ' CCOIN0, APHI, TEST',APHI,TEST
  865. C     write(ioimp,*) ' CCOIN0, KERRE = 2, NITER =',NITER
  866.           KERRE=-2
  867.           call soucis(268)
  868. C         GO TO 1000
  869.         ENDIF
  870.         IF(APHI.LE.TEST.OR.KERRE.LT.0) THEN
  871.           IF (KERRE.LT.0) KERRE = 0
  872.  
  873.           IF(MFR1.EQ.7) THEN
  874.             DO I=1,NCONT
  875.               S(I)=S(I)/DIV(I)
  876.               SPHER1(I)=SPHER1(I)/DIV(I)
  877.               SPHER2(I)=SPHER2(I)/DIV(I)
  878.             ENDDO
  879.           ENDIF
  880.  
  881. C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES-------------------
  882.  
  883.           IF(IFOURB.NE.-2) THEN
  884.             DO I=1,NCONT
  885.               DS(I)=S(I)-SIGB(I)
  886.             ENDDO
  887.             CALL EPSIG0(DS,DDEPSE,MFR1,IFOURB,YUNG,XNU,WORK,NSTRSS)
  888.             DO I=1,NCONT
  889.               DEPSE(I)=DEPST(I)+DDEPSE(I)
  890.               DEPSP(I)=DEPST(I)-DEPSE(I)
  891.             ENDDO
  892.           ENDIF
  893.  
  894.           DO I=1,NSTRSS
  895.             SIGF(I)=S(I)
  896.             DEFP(I)=DEPSP(I)
  897.           ENDDO
  898. c*dbg          if (ib.eq.1 .and. igau.eq.1) then
  899. c*dbg            write(ioimp,*) 'SIGF',(SIGF(iof),iof=1,NSTRSS)
  900. c*dbg          endif
  901.  
  902. C---------COQUES AVEC CISAILLEMENT TRANSVERSE---------------------------
  903.  
  904.           IF(MFR1.EQ.9) THEN
  905.             DEFP(7)=0.D0
  906.             DEFP(8)=0.D0
  907.             SIGF(7)=SIGB(7)
  908.             SIGF(8)=SIGB(8)
  909.           ENDIF
  910.  
  911.           VARF(1)=EPST
  912.           IF(INPLAS.EQ.4.OR.INPLAS.EQ.7.OR.INPLAS.EQ.12) THEN
  913.             DO I=1,NSTRSS
  914.               VARF(I+1)=SPHER1(I)
  915.             ENDDO
  916.             IF(INPLAS.EQ.12) THEN
  917.               DO I=1,NSTRSS
  918.                 VARF(NSTRSS+1+I)=SPHER2(I)
  919.               ENDDO
  920.             ENDIF
  921.           ELSE IF(INPLAS.EQ.87) THEN
  922.               DO I=1,NSTRSS
  923.                 VARF(1+I)=DEPST(I)
  924.               ENDDO
  925.           ENDIF
  926.           KERRE=0
  927. c*dbg          if (ib.eq.1 .and. igau.eq.1) then
  928. c*dbg            write(ioimp,*) 'VARF',EPST,'---',(VARF(iof),iof=2,1+NSTRSS)
  929. c*dbg          endif
  930.           RETURN
  931.  
  932.         ELSE
  933.           sx1in=sx(1)
  934.           sx2in=sx(2)
  935.           s1in=s(1)
  936.           s2in=s(2)
  937.           GOTO 10
  938.         ENDIF
  939.  
  940. 1000  CONTINUE
  941.  
  942. C      RETURN
  943.       END
  944.  
  945.  
  946.  
  947.  
  948.  
  949.  
  950.  
  951.  

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