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

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