Télécharger la plaquette Cast3M


fr en

Télécharger syco12.eso

Retour à la liste

Numérotation des lignes : 

syco12

  1. C SYCO12    SOURCE    OF166741  25/11/04    21:16:09     12349          
  2.  
  3.       SUBROUTINE SYCO12(wrk52,wrk53,wrk54,wrk2,
  4.      &  IB,IGAU,NBPGAU,NTRAC,IFOURL,iecou,DTPS)
  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                  TRAC       COURBE DE TRACTION
  21. C                  MFR1        INDICE DE FORMULATION
  22. C                  NSTRSS      NOMBRE DE CONTRAINTES CA2000
  23. C                  INPLAS     NUMERO DU MODELE DE PLASTICITE
  24. C                  DDAUX = MATRICE DE HOOKE ELASTIQUE
  25. C                  CMATE = NOM DU MATERIAU
  26. C                  VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
  27. C                  VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
  28. C                  N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
  29. C                  N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
  30. C                  IFOU  = OPTION DE CALCUL
  31. C                  IB    = NUMERO DE L ELEMENT COURANT
  32. C                  IGAU  = NUMERO DU POINT COURANT
  33. C                  EPAIST= EPAISSEUR
  34. C                  NBPGAU= NBRE DE POINTS DE GAUSS
  35. C                  MELE  = NUMERO DE L ELEMENT FINI
  36. C                  NPINT = NBRE DE POINTS D INTEGRATION
  37. C                  NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
  38. C                  NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
  39. C                  SECT  = SECTION
  40. C                  LHOOK = TAILLE DE LA MATRICE DE HOOKE
  41. C                  TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES
  42. C                  UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE
  43. C
  44. C     EN SORTIE :
  45. C
  46. C                  SIGF        CONTRAINTES A LA FIN DU PAS
  47. C                  VARF      VARIABLES INTERNES FIN DU PAS
  48. C                  DEFP        INCR. DE DEFORMATIONS PLASTIQUES
  49. C                  KERRE       CODE D'ERREUR
  50. C                   = 0  SI TOUT OK
  51. C                   = 99 SI FORMULATION NON DISPONIBLE
  52. C    EN PLASTICITE
  53. C                   = 1  SI DEPS EST NEGATIF
  54. C                   = 2  SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE
  55. C
  56. C         IFOUR : OPTION DE CALCUL
  57. C
  58. C         IFOUR = -3   DEFORMATION PLANE GENERALISEE
  59. C         IFOUR = -2   CONTRAINTES PLANES
  60. C         IFOUR = -1   DEFORMATIONS PLANES
  61. C         IFOUR =  0   AXISYMETRIQUE
  62. C         IFOUR =  1   SERIE DE FOURIER
  63. C         IFOUR =  2   TRIDIMENSIONNEL
  64. C
  65. C         MFR1 : NUMERO DE LA FORMULATION ELEMENT FINI
  66. C
  67. C         MFR1 = 1     MASSIF
  68. C         MFR1 = 63     XFEM
  69. C
  70. C         ISOTROPE
  71. C         INPLAS : NUMERO DU MODELE DE VISCOPLASTICITE
  72. C
  73. C         INPLAS = 153      SYCO1
  74. C         INPLAS = 154      SYCO2
  75. C
  76. C-----------------------------------------------------------------------
  77. C         CONVENTION DE REMPLISSAGE DES MEMOIRES
  78. C-----------------------------------------------------------------------
  79. C
  80. C     XMAT(1)    MODULE D'YOUNG
  81. C     XMAT(2)    COEFFICIENT DE POISSON
  82. C
  83. C     TRAC(1 A 2*NTRAC) COURBE DE TRACTION
  84. C
  85. C    MODELE ISOTROPE
  86. C            VAR0(1)  EPSE
  87. C            VAR0(2)  VP
  88. C
  89. C
  90. C-----------------------------------------------------------------------
  91. C     20/09/2017 : modif SG critere de convergence trop serre
  92. C                  TEST=PETI*APHI0 + utilisation CCREEL
  93. C                  idem ccoin0.eso
  94.  
  95.       IMPLICIT INTEGER(I-N)
  96.       IMPLICIT REAL*8(A-H,O-Z)
  97.  
  98. -INC PPARAM
  99. -INC CCOPTIO
  100. -INC CCREEL
  101. -INC DECHE
  102.  
  103. -INC TECOU
  104.  
  105.       SEGMENT WRK2
  106.         REAL*8 TRAC(NTRAC)
  107.       ENDSEGMENT
  108.  
  109. C-----------------------------------------------------------------------
  110.       DIMENSION SIG(130),EPS(130)
  111.       DIMENSION S(8),SX(8),DS(8),DSIG(8),SPHER(8),SPHER1(8),SPHER2(8)
  112.       DIMENSION DSPHER1(8),DSPHER2(8),DEPSE(8),DEPSP(8),DDEPSE(8)
  113.       DIMENSION F(8),WF1(8),WF2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8)
  114.       DIMENSION CRIGI(12)
  115.  
  116.       KERRE=0
  117.  
  118.       NCONT = iecou.NSTRSS
  119.       MFRl  = iecou.MFR1
  120.       jnplas = INPLAS
  121. C
  122. C--------PAS DE TEMPS -----------------------
  123.       TSYC9 = DTPS
  124.       itracb=1
  125.  
  126. C-----------------------------------------------------------------------
  127.       DO I=1, NCONT
  128.         S(I)=0.D0
  129.         SPHER(I)=0.D0
  130.         SPHER1(I)=0.D0
  131.         SPHER2(I)=0.D0
  132.       ENDDO
  133.       YUNG=XMAT(1)
  134.       XNU=XMAT(2)
  135.  
  136. C---------CARACTERISTIQUES GEOMETRIQUES---------------------------------
  137.       EP1=1.D0
  138.       ALFAC=1.D0
  139.  
  140.       IF (jnplas.EQ.153) THEN
  141. C     WRITE(6,*) 'INPLAS', jnplas
  142. C     WRITE(6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5)
  143. C     WRITE(6,*) 'XMAT', xmat(6), xmat(7), xmat(8)
  144.         PSYC9=XMAT(6)
  145.         DSYC9=XMAT(7)
  146.       ELSE IF (jnplas.EQ.154) THEN
  147. C     WRITE(6,*) 'INPLAS', jnplas
  148. C     WRITE(6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5)
  149. C     WRITE(6,*) 'XMAT', xmat(6), xmat(7), xmat(8),xmat(9),xmat(10)
  150.         PSYC9=XMAT(6)
  151.         ASYC9=XMAT(7)
  152.         BSYC9=XMAT(8)
  153.         CSYC9=XMAT(9)
  154.       ELSE
  155.       ENDIF
  156.  
  157.       DO I=1,NTRAC
  158.         SIG(I)=TRAC(2*I-1)
  159.         EPS(I)=TRAC(2*I)
  160.       ENDDO
  161.  
  162.       EPST=VAR0(1)
  163. C
  164. C     VP=VAR0(2)
  165. C     if(VP.LE.1.d-4) VP = 1.d-4
  166. C        VP0 = 1.D-4
  167. C     IF (jnplas.EQ.153) THEN
  168. C        deno1 = 1.D0 + (  (VP0/DSYC9)**(1.D0/PSYC9) )
  169. C        deno2 = 1.D0 + (  (VP/DSYC9)**(1.D0/PSYC9) )
  170. C     else if (jnplas.EQ.154) then
  171. C        HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9))
  172. C        deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9)))
  173. C        deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9)))
  174. C     endif
  175. C
  176.  
  177. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  178. C-xxxx pour le materiau syco1 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  179. C        DSYC9  parametre D de la loi Symonds Cowper
  180. C        PSYC9  parametre P de la loi Symonds Cowper
  181. C        VP valeur minimum de la vitesse de deformation equivalente
  182.       IF (jnplas.EQ.153) THEN
  183.          VP0 = 1.D-4
  184.         VP = VP0
  185.       else if(jnplas.EQ.154) then
  186. C        VP valeur minimum de la vitesse de deformation equivalente
  187.          VP0 = 1.D-4
  188.          VP = VP0
  189.       endif
  190. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  191.       CALL CALSIG(DEPST,DDAUX,NCONT,CMATE,VALMAT,VALCAR,
  192.      &        N2EL,N2PTEL,MFRl,IFOURL,IB,IGAU,EPAIST,
  193.      &      NBPGAU,MELE,NPINT,iecou.NBGMAT,iecou.NELMAT,SECT,LHOOK,TXR,
  194.      &    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
  195.  
  196.       IF(IRTD.NE.1) THEN
  197.          KERRE=69
  198.          GOTO 1000
  199.       ENDIF
  200.  
  201.       IF (ifourl.eq.-2) then
  202. * en cas de contraintes planes on repasse en 3D
  203.         do ju=1,6
  204.            do iu=1,6
  205.              DDA(iu,ju)=0.d0
  206.            enddo
  207.         enddo
  208.         cte_cp = xnu / (1.d0 - xnu)
  209.         aux= YUNG / ((1.d0+xnu)*(1.d0-2.d0*xnu))
  210.         aux1= aux * xnu
  211.         aux2= aux * (1.d0-xnu)
  212.         gege = Yung / (2.d0 *(1.d0 +xnu))
  213.         DDA(1,1)=AUX2
  214.         DDA(2,1)=AUX1
  215.         DDA(1,2)=AUX1
  216.         DDA(2,2)=AUX2
  217.         DDA(3,3)=aux2
  218.         DDA(1,3)=aux1
  219.         DDA(2,3)=aux1
  220.         DDA(3,1)=aux1
  221.         DDA(3,2)=aux1
  222.  
  223.         DDA(4,4)=gege
  224.         DDA(5,5)=gege
  225.         DDA(6,6)=GEGE
  226.       ENDIF
  227.  
  228.       DO I=1,NCONT
  229.         S(I)=SIG0(I)+DSIGT(I)
  230.         SIGB(I)=S(I)
  231.         SX(I)=S(I)-SPHER(I)
  232.       ENDDO
  233.  
  234. C-----------------------------------------------------------------------
  235. C         CALCUL DE LA LIMITE ELASTIQUE SI
  236. C-----------------------------------------------------------------------
  237. C
  238.         CALL TRACTI(SI,EPST,SIG,EPS,NTRAC,2,izz)
  239. c
  240.         IF (izz.EQ.1) THEN
  241.           KERRE=75
  242.           GOTO 1000
  243.         ENDIF
  244.  
  245. C-----------------------------------------------------------------------
  246. C         CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
  247. C-----------------------------------------------------------------------
  248. *  attention en contraintes planes on se declare en 3D
  249. * rien besoin de faire dans vonmis0 car ifourl n'est pas utilise
  250. * et les contraintes sont dimensionnees a 6
  251. *       if(igau.eq.1) then
  252. *       WRITE(6,*) ' en entree sx ',sx(1),sx(2),sx(3)
  253. *       endif
  254. C
  255. C     SEQ=VONMIS0(SX,NCONT,MFRl,ifourl,EP1,ALFAC)
  256. C calcul de la contrainte de Von Mises
  257.         SOM1=0.D0
  258.         DO I=4,NCONT
  259.           SOM1=SOM1+SX(I)*SX(I)
  260.         ENDDO
  261.  
  262.         SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3)
  263.      &  -SX(2)*SX(3)+3.D0*SOM1
  264.         SEQ=SQRT(ABS(SEQ))
  265.  
  266. C-----------------------------------------------------------------------
  267. C         LE CRITERE EST-IL VERIFIE ?
  268. C-----------------------------------------------------------------------
  269. C     if(itracb.eq.1) then
  270. C      WRITE(6,*)' en entree si seq', si, seq
  271. C     endif
  272. C
  273. c  test de plasticite sur SI statique
  274. C     SI = SI * deno2/deno1
  275. C
  276.       PHI=SEQ-SI
  277.       NITER=0
  278.       PETI=1.1D0*0.5D0*PRECIS*SEQ
  279.       CALL EPSPRE(SEQ,SI,PETI,ITRY)
  280. C
  281.       IF((ITRY.EQ.1).OR.(SEQ.LE.SI)) THEN
  282. * rien a faire on n'a pas plastifie
  283.         DO I=1,NCONT
  284.           SIGF(I)=S(I)
  285.           DEFP(I)=0.D0
  286.         ENDDO
  287.  
  288.         VARF(1)=VAR0(1)
  289.         VARF(2)=VP0
  290.         RETURN
  291.       ENDIF
  292.  
  293. C-----------------------------------------------------------------------
  294. C         ON A PLASTIFIE
  295. C-----------------------------------------------------------------------
  296. C
  297.         cri0=  si * 1.d-8
  298.         PHI0=PHI
  299.         SI0 = SI
  300.         RR=0.D0
  301.  
  302.         DO I=1,NCONT
  303.           DSIG(I)=0.D0
  304.           DEPSP(I)=0.D0
  305.           DSPHER1(I)=0.D0
  306.           DSPHER2(I)=0.D0
  307.         ENDDO
  308.  
  309. C-----------------------------------------------------------------------
  310. C         DEBUT DE LA BOUCLE D'ITERATIONS INTERNES
  311. C-----------------------------------------------------------------------
  312.         sx1in=sx(1)
  313.         sx2in=sx(2)
  314.         sx3in= sx(3)
  315.         s1in=s(1)
  316.         s2in=s(2)
  317.         iderin=0
  318.  
  319. 10      CONTINUE
  320.          NITER=NITER+1
  321.  
  322.          phi= seq - si
  323.  
  324. C---------CALCUL DE WF1=DF/D(SIGMA)--------------------------------------
  325.  
  326. C---------ELEMENTS MASSIFS----------------------------------------------
  327.  
  328.           F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
  329.           F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
  330.           F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
  331.           DO I=4,NCONT
  332.             F(I)=SX(I)
  333.           ENDDO
  334.           DO I=1,3
  335.             WF1(I)=1.5D0*F(I)/SEQ
  336.             Z1(I)=WF1(I)
  337.           ENDDO
  338.           DO I=4,NCONT
  339.             WF1(I)=3.D0*F(I)/SEQ
  340.             Z1(I)=1.5D0*F(I)/SEQ
  341.           ENDDO
  342.  
  343.       if(ifourl.eq.-2) then
  344.          DO I=1,NCONT
  345.             r_z = 0.D0
  346.             DO J=1,NCONT
  347.               r_z=r_z+DDA(I,J)*WF1(J)
  348.             ENDDO
  349.             WF2(I)=r_z
  350.           ENDDO
  351.       else
  352.          DO I=1,NCONT
  353.             r_z = 0.D0
  354.             DO J=1,NCONT
  355.               r_z=r_z+DDAUX(I,J)*WF1(J)
  356.             ENDDO
  357.             WF2(I)=r_z
  358.           ENDDO
  359.       endif
  360.         COEF=0.D0
  361.         DO I=1,NCONT
  362.           COEF=COEF+WF1(I)*WF2(I)
  363.         ENDDO
  364.  
  365. C-----------------------------------------------------------------------
  366. C      ECROUISSAGE ISOTROPE
  367. C-----------------------------------------------------------------------
  368.           CALL TRACTI(PENTE,EPST,SIG,EPS,NTRAC,1,izz)
  369. *          WRITE(6,*) ' pente  epst', pente,epst
  370.           IF (izz.EQ.1) THEN
  371.             KERRE=75
  372.             GOTO 1000
  373.           ENDIF
  374.  
  375.           RP=PENTE
  376.           C=0.D0
  377.  
  378.           DENOM=COEF+C
  379. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  380. C---------CALCUL des termes en sup dus a la viscosite----------------
  381. C    il s'agit de -qsi1 = -d(phi)/d(vp) = +d(SY)/d(vp)
  382. c    (ou  vp = dp/dt  cad la vitesse inelastiq equivalente)
  383.       IF (jnplas.EQ.153) THEN
  384.          deno1 = 1.D0 + (  (VP0/DSYC9)**(1.D0/PSYC9) )
  385.          deno2 = 1.D0 + (  (VP/DSYC9)**(1.D0/PSYC9) )
  386.  
  387.          qsi1 = (SI/(PSYC9*DSYC9))*((VP/DSYC9)**( (1.D0/PSYC9)-1.D0))
  388.          qsi1 = qsi1 / (deno2 * TSYC9)
  389.          DENOM = DENOM + qsi1
  390.          DENOM = DENOM + (RP * deno2/deno1)
  391. C
  392.         ELSE IF (jnplas.EQ.154) THEN
  393.          HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9))
  394.          deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9)))
  395.          deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9)))
  396.  
  397.          qsi1 = (SI*HSYC9/PSYC9) * (VP**((1.D0/PSYC9)-1.D0))
  398.          qsi1 = qsi1 / (deno2 * TSYC9)
  399.          DENOM = DENOM + qsi1
  400.          DENOM = DENOM + (RP * deno2/deno1)
  401. C
  402.          qsi2 = (SI/deno2)*(VP**(1.D0/PSYC9))*((ASYC9-HSYC9)/CSYC9)
  403.          DENOM = DENOM + qsi2
  404. C
  405.         ENDIF
  406. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  407. C23456789012345678901234567890123456789012345678901234567890123456789012
  408.           DELTA=PHI/DENOM
  409.           DMU=C*DELTA/SEQ
  410. C          if(itracb.eq.1) then
  411. C           WRITE(6,*)'pente coef  denom delta dmu  niter si seq'
  412. C            WRITE(6,*)pente ,coef , denom ,delta ,dmu,niter,si,seq
  413. C           endif
  414.  
  415.           DO I=1,NCONT
  416.             DSIG(I)=-DELTA*WF2(I)
  417.             DSPHER1(I)=DMU*SX(I)
  418.           ENDDO
  419.  
  420. * Cas des contraintes planes en massif
  421.           if(ifourl.eq.-2) then
  422.  
  423.              bb= abs(dsig(3)+ s(3) )
  424.  
  425. C             if(itracb.eq.1) then
  426. C               WRITE(6,*) ' iderin bb ' , iderin , bb
  427. C             endif
  428.  
  429.              r_z = dsig(3) * cte_cp
  430.              sx(3)= sx3in - dsig(3)
  431.              sx(1)= sx1in - r_z
  432.              sx(2)= sx2in - r_z
  433. C
  434. C            SEQ=VONMIS0(SX,NCONT,MFRl,ifourl,EP1,ALFAC)
  435. C calcul de la contrainte de Von Mises
  436.         SOM1=0.D0
  437.         DO I=4,NCONT
  438.           SOM1=SOM1+SX(I)*SX(I)
  439.         ENDDO
  440.  
  441.         SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3)
  442.      &  -SX(2)*SX(3)+3.D0*SOM1
  443.         SEQ=SQRT(ABS(SEQ))
  444. C
  445.               s(3)=  - dsig(3)
  446.               s(1)= s1in - r_z
  447.               s(2)= s2in - r_z
  448.  
  449.         cri0=  si * 1.d-8
  450.               if( bb.gt.cri0) then
  451.  
  452.                 if(iderin.eq.0) then
  453.                   niter=niter - 1
  454.                 endif
  455.  
  456. C                if(itracb.eq.1) then
  457. C               WRITE(6,* ) ' cri0 bb seq  delta ',cri0,bb,seq , delta
  458. C               endif
  459.  
  460.                iderin = iderin + 1
  461.  
  462. **               WRITE(6,*) ' iderin ' , iderin
  463.  
  464.                if(iderin.gt.50) then
  465.                WRITE(6,*) 'on iderin superieur a 50'
  466.                  write(6,*) '   probleme dans iterations internes'
  467.                           KERRE=2
  468.                           GO TO 1000
  469.                endif
  470.  
  471. C23456789012345678901234567890123456789012345678901234567890123456789012
  472. C     WRITE(6,*)'cri0 bb niter iderin',cri0,bb,niter,iderin,dsig(3),s(3)
  473.                go to 10
  474.               endif
  475.  
  476.               DMU=C*DELTA/SEQ
  477.               DO I=1,NCONT
  478.                 DSPHER1(I)=DMU*SX(I)
  479.               ENDDO
  480.          endif
  481.  
  482.           iderin=0
  483. C
  484.         DR=RP*DELTA
  485.         RR=RR+DR
  486.         EPST=EPST+DELTA
  487.  
  488. C mise a jour des contraintes
  489.           DO I=1,NCONT
  490.             S(I)=S(I)+DSIG(I)
  491.             SPHER1(I)=SPHER1(I)+DSPHER1(I)
  492.             SPHER2(I)=SPHER2(I)+DSPHER2(I)
  493.             SPHER(I)=SPHER1(I)+SPHER2(I)
  494.             SX(I)=S(I)-SPHER(I)
  495.           ENDDO
  496.           if(ifourl.eq.-2) then
  497.             s(3)=0.d0
  498.           endif
  499. C           if(igau.eq.1) then
  500. C            WRITE(6,*) ' niter',niter
  501. C           WRITE(6,*) ' dsig' , dsig(1),dsig(2),dsig(3),dsig(4)
  502. C           WRITE(6,*) ' s' , s(1),s(2),s(3),s(4)
  503. C           WRITE(6,*) ' sx' , sx(1),sx(2),sx(3),sx(4)
  504. C           WRITE(6,*) ' spher' , spher(1),spher(2),spher(3),spher(4)
  505. C           endif
  506.  
  507. C     SEQ=VONMIS0(SX,NCONT,MFRl,ifourl,EP1,ALFAC)
  508. C calcul de la contrainte de Von Mises
  509.         SOM1=0.D0
  510.         DO I=4,NCONT
  511.           SOM1=SOM1+SX(I)*SX(I)
  512.         ENDDO
  513.  
  514.         SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3)
  515.      &  -SX(2)*SX(3)+3.D0*SOM1
  516.         SEQ=SQRT(ABS(SEQ))
  517. C
  518. C---------CONTRAINTES PLANES--------------------------------------------
  519.  
  520.         IF(ifourl.EQ.-2) THEN
  521.  
  522.             F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
  523.             F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
  524.             F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
  525.             DO I=4,NCONT
  526.               F(I)=SX(I)
  527.             ENDDO
  528.             DO I=1,3
  529.               WF1(I)=1.5D0*F(I)/SEQ
  530.             ENDDO
  531.             DO I=4,NCONT
  532.               WF1(I)=3.D0*F(I)/SEQ
  533.             ENDDO
  534.  
  535.             DO I=1,NCONT
  536.               DEPSP(I)=DEPSP(I)+DELTA*WF1(I)
  537.             ENDDO
  538.         ENDIF
  539.  
  540. C-----------------------------------------------------------------------
  541. C         TEST
  542. C         CALCUL DE LA NOUVELLE VALEUR DE PHI
  543. C-----------------------------------------------------------------------
  544.            CALL TRACTI(SI,EPST,SIG,EPS,NTRAC,2,izz)
  545. C******************************************************
  546. c-xxxx  calcul de la valeur (i+1) de H, de la vitesse de def VP,
  547. c       et de SI en prenant en compte la viscosite
  548.       IF(jnplas.EQ.153) THEN
  549. C materiau  SYCO1
  550.           deno1 = 1.D0 + (  (VP0/DSYC9)**(1.D0/PSYC9) )
  551.               VP = MAX( (VP+(DELTA/TSYC9)) , VP0 )
  552. C         deno2 = 1.D0 + (  (VP/DSYC9)**(1.D0/PSYC9) )
  553.               SI = SI * (1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) )) / deno1
  554. C
  555.            else if (jnplas.EQ.154) then
  556. C materiau SYCO2
  557.               HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9))
  558.           deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9)))
  559.               VP = MAX( (VP+(DELTA/TSYC9)) , VP0 )
  560. C         deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9)))
  561.               SI = SI * (1.D0 + (HSYC9 * (VP**(1.D0/PSYC9)))) / deno1
  562.            endif
  563. C******************************************************
  564.  
  565.         PHI=SEQ-SI
  566.  
  567.         PETI=1.D-7
  568.         APHI=ABS(PHI)
  569.         APHI0=ABS(PHI0)
  570.         TEST=max(PETI*APHI0,XZPREC*100.D0*SEQ)
  571. *sg        TEST=PETI*APHI0
  572.  
  573.         IF(NITER.GT.200) THEN
  574.           WRITE(6,*) 'on a bien plante ici,car APHI=',APHI
  575.           WRITE(6,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI
  576.           TEST2=max(1.D-6*APHI0,XZPREC*100.D0*SEQ)
  577. *sg          TEST2 = 1.D-6*APHI0
  578.           if(APHI.gt.TEST2) then
  579.             KERRE=2
  580.             GO TO 1000
  581.           else
  582.             TEST=TEST2
  583.           endif
  584.         ELSEIF(NITER.eq.101) THEN
  585.           WRITE(6,*) 'on aurait du plante ici,car APHI=',APHI
  586.           WRITE(6,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI
  587.           WRITE(6,*) 'mais on va continuer 50pas de plus'
  588.         ENDIF
  589.  
  590. *         if(igau.eq.1) then
  591. *        WRITE(6,*)'test conver niter phi seq si testphi0',niter,
  592. *     $    phi,seq,si,test,aphi0
  593. *        WRITE(6,*) ' rr si0' , rr, si0
  594. *         endif
  595.  
  596.         IF(APHI.LE.TEST) THEN
  597. C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES-------------------
  598.  
  599.           IF(ifourl.NE.-2) THEN
  600.             DO I=1,NCONT
  601.               DS(I)=S(I)-SIGB(I)
  602.             ENDDO
  603. C copie d'un bout de EPSIG0
  604.             AUX1=1.D0/YUNG
  605.             AUX2=2.D0*(1.D0+XNU)
  606.             AUX3=AUX1*AUX2
  607.             DDEPSE(1)=(DS(1)-XNU*(DS(2)+DS(3)))*AUX1
  608.             DDEPSE(2)=(DS(2)-XNU*(DS(1)+DS(3)))*AUX1
  609.             DDEPSE(3)=(DS(3)-XNU*(DS(1)+DS(2)))*AUX1
  610.             DO I=4,NCONT
  611.               DDEPSE(I)=AUX3*DS(I)
  612.             ENDDO
  613. C fin copie d'un bout de EPSIG0
  614.             DO I=1,NCONT
  615.               DEPSE(I)=DEPST(I)+DDEPSE(I)
  616.               DEPSP(I)=DEPST(I)-DEPSE(I)
  617.             ENDDO
  618.           ENDIF
  619.  
  620.           DO I=1,NCONT
  621.             SIGF(I)=S(I)
  622.             DEFP(I)=DEPSP(I)
  623.           ENDDO
  624.  
  625.           VARF(1)=EPST
  626.           VARF(2)=VP
  627.  
  628.           KERRE=0
  629.           RETURN
  630.  
  631.         ELSE
  632.           sx1in=sx(1)
  633.           sx2in=sx(2)
  634.           s1in=s(1)
  635.           s2in=s(2)
  636.  
  637.           GOTO 10
  638.         ENDIF
  639.  
  640. 1000  CONTINUE
  641.       RETURN
  642.       END
  643.  
  644.  
  645.  

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