Télécharger la plaquette Cast3M


fr en

Télécharger ouglov.eso

Retour à la liste

Numérotation des lignes : 

ouglov

  1. C OUGLOV    SOURCE    BR232186  16/09/12    12:43:41     9078           
  2.       SUBROUTINE OUGLOV(XMAT,VAR0,VARF,SIG0,SIGF,DEPST,IFOUR)
  3.  
  4. C
  5.  
  6. C====&===1=========2=========3=========4=========5=========6=========7==
  7.  
  8. C Commentaires : Subroutine permettant de mettre en oeuvre le
  9.  
  10. C                modele OUGLOVA pour representer
  11.  
  12. C                le comportement de l'acier corrodé
  13.  
  14. C
  15.  
  16. C Auteurs      : R. Paredes
  17.  
  18. C              : CEA-DEN/DANS/DM2S/SEMT/EMSI
  19.  
  20. C              : Romili.Paredes@cea.fr
  21.  
  22. C
  23.  
  24. C Date         : Aout 2016
  25.  
  26. C====&===1=========2=========3=========4=========5=========6=========7==
  27.  
  28. C
  29.  
  30. C-----DECLARATION GENERALE----------------------------------------------
  31.  
  32. C
  33.  
  34.       IMPLICIT INTEGER(I-N)
  35.       IMPLICIT REAL*8(A-H,O-Z)
  36.  
  37. C
  38.  
  39. C-----DECLARATION DES VARIABLES-----------------------------------------
  40.  
  41. C
  42.  
  43.       DIMENSION XMAT(*),DEPST(*),VAR0(*),SIGF(*),VARF(*)
  44.  
  45.  
  46.  
  47.       REAL*8 UNIT3(3,3), DFEDPHI(3)
  48.  
  49.       REAL*8 EPSR(3), XEPSP(3,3), XEPSF(3,3), XEPSE(3,3), SIG(3,3)
  50.  
  51.       REAL*8 DFDSG(3,3), DPHDSG(3,3), DEPSP(3,3), DVSIG(3,3), EDPH(3,3)
  52.  
  53.  
  54.  
  55. C-----Parametre pour le nombre d iteration locales internes
  56.  
  57.       IMAX = 1000
  58.  
  59. C
  60.  
  61. C-----DEFINITION DE LA MATRICE UNITE D ORDRE 3--------------------------
  62.  
  63. C
  64.  
  65.       DO I = 1,3
  66.  
  67.          DO J = 1,3
  68.  
  69.             IF (I.EQ.J) THEN
  70.  
  71.                UNIT3(I,J) = 1.0D0
  72.  
  73.             ELSE
  74.  
  75.                UNIT3(I,J) = 0.0D0
  76.  
  77.             ENDIF
  78.  
  79.          ENDDO
  80.  
  81.       ENDDO
  82.  
  83. C
  84.  
  85. C-----DONNEES MATERIAUX-------------------------------------------------
  86.  
  87. C
  88.  
  89.       XYG   = XMAT(1)
  90.  
  91.       XNU   = XMAT(2)
  92.  
  93.       XSIGY = XMAT(5)
  94.  
  95.       XK    = XMAT(6)
  96.  
  97.       XM    = XMAT(7)
  98.  
  99.       XTC   = XMAT(8)
  100.  
  101.       XDC   = XMAT(9)
  102.  
  103.  
  104.  
  105.       XLA   = (XNU*XYG)/((1.0D0 - 2.0D0*XNU)*(1.0D0 + XNU))
  106.  
  107.       XMU   = XYG/(2.0D0*(1.0D0 + XNU))
  108.  
  109. C
  110.  
  111. C-----CORROSION---------------------------------------------------------
  112.  
  113. C   
  114.  
  115. C-----Deformation plastique a rupture
  116.  
  117.       IF (XTC.LE.15.0D0) THEN
  118.  
  119.          EPSR(1) = -0.0111D0*XTC + 0.2345D0
  120.  
  121.       ELSE 
  122.  
  123.          EPSR(1) = -0.0006D0*XTC + 0.051D0
  124.  
  125.       ENDIF
  126.  
  127.  
  128.  
  129. C-----Coefficient de contraction      
  130.  
  131.       XNUSTAR = 0.5D0 - (XSIGY/(XYG*EPSR(1)))*(0.5D0 - XNU)
  132.  
  133.  
  134.  
  135.       EPSR(2) = -EPSR(1)*XNUSTAR
  136.  
  137.       EPSR(3) = -EPSR(1)*XNUSTAR
  138.  
  139.  
  140.  
  141. C-----Seuil de rupture et d'endommagement
  142.  
  143.       TEPSR = EPSR(1)**2.0D0 + EPSR(2)**2.0D0 + EPSR(3)**2.0D0
  144.  
  145.       XPR = ((2.0D0/3.0D0)*TEPSR)**(0.5D0)
  146.  
  147.       XPD = 0.8D0*XPR
  148.  
  149. C
  150.  
  151. C-----ETAT INITIAL------------------------------------------------------
  152.  
  153. C 
  154.  
  155.       XD    = VAR0(1)
  156.  
  157.       XR    = VAR0(2)
  158.  
  159.       XP    = VAR0(3)
  160.  
  161.       XZT   = VAR0(4)
  162.  
  163.       XNRUP = VAR0(5)
  164.  
  165.  
  166.  
  167. C-----------------------------------------------------------------------
  168.  
  169. C-----CALCUL TRIDIMENSIONNEL--------------------------------------------
  170.  
  171. C-----------------------------------------------------------------------    
  172.  
  173.       IF (IFOUR.EQ.2) THEN
  174.  
  175.  
  176.  
  177. C--------Deformation plastique      
  178.  
  179.          XEPSP(1,1) = VAR0(6)
  180.  
  181.          XEPSP(2,2) = VAR0(7)
  182.  
  183.          XEPSP(3,3) = VAR0(8)
  184.  
  185.          XEPSP(1,2) = VAR0(9)/2.0D0
  186.  
  187.          XEPSP(1,3) = VAR0(10)/2.0D0
  188.  
  189.          XEPSP(2,3) = VAR0(11)/2.0D0
  190.  
  191.  
  192.  
  193.          XEPSP(2,1) = XEPSP(1,2)
  194.  
  195.          XEPSP(2,3) = XEPSP(2,3)
  196.  
  197.          XEPSP(3,1) = XEPSP(1,3)
  198.  
  199.          XEPSP(3,2) = XEPSP(2,3)
  200.  
  201.  
  202.  
  203. C--------Deformation totale        
  204.  
  205.          XEPSF(1,1) = DEPST(1) + VAR0(12)
  206.  
  207.          XEPSF(2,2) = DEPST(2) + VAR0(13)
  208.  
  209.          XEPSF(3,3) = DEPST(3) + VAR0(14)
  210.  
  211.          XEPSF(1,2) = (DEPST(4) + VAR0(15))/2.0D0
  212.  
  213.          XEPSF(1,3) = (DEPST(5) + VAR0(16))/2.0D0
  214.  
  215.          XEPSF(2,3) = (DEPST(6) + VAR0(17))/2.0D0
  216.  
  217.  
  218.  
  219.          XEPSF(2,1) = XEPSF(1,2)
  220.  
  221.          XEPSF(2,3) = XEPSF(2,3)
  222.  
  223.          XEPSF(3,1) = XEPSF(1,3)
  224.  
  225.          XEPSF(3,2) = XEPSF(2,3)
  226.  
  227.  
  228.  
  229. C--------Deformation elastique        
  230.  
  231.          DO I = 1,3
  232.  
  233.             DO J = 1,3
  234.  
  235.                XEPSE(I,J) = XEPSF(I,J) - XEPSP(I,J)
  236.  
  237.             ENDDO
  238.  
  239.          ENDDO
  240.  
  241.  
  242.  
  243. C--------Prediction elastique
  244.  
  245.          TREPSE = XEPSE(1,1) + XEPSE(2,2) + XEPSE(3,3)
  246.  
  247.          DO I = 1,3
  248.  
  249.             DO J = 1,3
  250.  
  251.                AUX1 = XLA*TREPSE*UNIT3(I,J)
  252.  
  253.                AUX2 = 2.0D0*XMU*XEPSE(I,J)
  254.  
  255.                SIG(I,J) = (1.0D0 - XD)*(AUX1 + AUX2)
  256.  
  257.             ENDDO
  258.  
  259.          ENDDO
  260.  
  261.          CALL J2(SIG,SIGEQ)
  262.  
  263.          CALL DEVIAT(SIG,DVSIG)
  264.  
  265.  
  266.  
  267. C--------Boucle d'endommagement         
  268.  
  269.          DO I = 1,IMAX
  270.  
  271. C-----------Evaluation de la fonction seuil
  272.  
  273.             FP  = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
  274.  
  275.             FP0 = FP
  276.  
  277.  
  278.  
  279. C-----------Boucle de plasticite      
  280.  
  281.             IF (FP.GT.1.0D0) THEN  
  282.  
  283.  
  284.  
  285.                DO J = 1,IMAX
  286.  
  287. C-----------------Calcul des derives               
  288.  
  289.                   DO K = 1,3
  290.  
  291.                      DO L = 1,3
  292.  
  293.                         AUX1 = ((3.0D0/2.0D0)/(SIGEQ*(1.0D0 - XD)))
  294.  
  295.                         DFDSG(K,L) = AUX1*DVSIG(K,L)
  296.  
  297.                      ENDDO
  298.  
  299.                   ENDDO
  300.  
  301.  
  302.  
  303.                   DO K = 1,3
  304.  
  305.                      DO L = 1,3
  306.  
  307.                         DPHDSG(K,L) = DFDSG(K,L)
  308.  
  309.                      ENDDO
  310.  
  311.                   ENDDO
  312.  
  313.  
  314.  
  315. C-----------------Calcul du multiplicateur plastique
  316.  
  317.                   TRDPHDSG = DPHDSG(1,1) + DPHDSG(2,2) + DPHDSG(3,3)
  318.  
  319.                   DO K = 1,3
  320.  
  321.                      DO L = 1,3
  322.  
  323.                         AUX1 = XLA*TRDPHDSG*UNIT3(K,L)
  324.  
  325.                         AUX2 = 2.0D0*XMU*DPHDSG(K,L)
  326.  
  327.                         EDPH(K,L) = (AUX1 + AUX2)
  328.  
  329.                      ENDDO
  330.  
  331.                   ENDDO
  332.  
  333.                   DO K = 1,3
  334.  
  335.                      DFEDPHI(K) = 0.0D0
  336.  
  337.                      DO L = 1,3
  338.  
  339.                         AUX1 = DFDSG(K,L)*EDPH(L,K)
  340.  
  341.                         DFEDPHI(K) = DFEDPHI(K) + AUX1
  342.  
  343.                      ENDDO
  344.  
  345.                   ENDDO
  346.  
  347.                   AUX4  = DFEDPHI(1) + DFEDPHI(2) + DFEDPHI(3)
  348.  
  349.                   AUX5  = XK/XM
  350.  
  351.                   AUX6  = SIGEQ/(XK*(1.0D0 - XD)) 
  352.  
  353.                   AUX7  = XSIGY/XK
  354.  
  355.                   AUX8  = (1.0D0 - XM)           
  356.  
  357.                   XLAMBDA = FP/(AUX4 + AUX5*((AUX6 + AUX7)**AUX8))
  358.  
  359.  
  360.  
  361. C-----------------Mise a jour des variables
  362.  
  363.                   DO K = 1,3
  364.  
  365.                      DO L = 1,3
  366.  
  367.                         SIG(K,L) = SIG(K,L) - XLAMBDA*EDPH(K,L)
  368.  
  369.                      ENDDO
  370.  
  371.                   ENDDO
  372.  
  373.                   CALL J2(SIG,SIGEQ)
  374.  
  375.                   CALL DEVIAT(SIG,DVSIG)
  376.  
  377.  
  378.  
  379.                   DO K = 1,3
  380.  
  381.                      DO L = 1,3
  382.  
  383.                         XEPSP(K,L) = XEPSP(K,L) + XLAMBDA*DPHDSG(K,L)
  384.  
  385.                      ENDDO
  386.  
  387.                   ENDDO
  388.  
  389.                   XP   = XP + XLAMBDA/(1.0D0 - XD);
  390.  
  391.                   XR   = XK*(XP**(1.0D0/XM))
  392.  
  393.  
  394.  
  395. C-----------------Reevaluation de la fonction seuil
  396.  
  397.                   FP   = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
  398.  
  399.  
  400.  
  401.                   CRIT1 = DABS(FP/FP0)
  402.  
  403.                   CRIT2 = DABS(XLAMBDA)
  404.  
  405.  
  406.  
  407.                   IF ((CRIT1.LE.1.0D-8).OR.(CRIT2.LE.1.0D-10)) THEN 
  408.  
  409.                       GOTO 910
  410.  
  411.                   ENDIF
  412.  
  413.  
  414.  
  415.                ENDDO
  416.  
  417.  
  418.  
  419. C-----Sortie de la boucle de plasticité
  420.  
  421.  910  CONTINUE
  422.  
  423.  
  424.  
  425. C--------------Evaluation du seuil d'endommagement
  426.  
  427.                FD = XP - (XPD + XZT)
  428.  
  429.                IF (I.EQ.1) THEN
  430.  
  431.                    FD0 = FD
  432.  
  433.                ENDIF
  434.  
  435.  
  436.  
  437.                CRIT = DABS(FD/FD0)
  438.  
  439.  
  440.  
  441.                IF ((FD.GT.0.0D0).AND.(CRIT.GT.1.0D-5)) THEN
  442.  
  443.  
  444.  
  445. C-----------------Calcul de l'endommagement
  446.  
  447.                   XD  = (XDC/(XPR - XPD))*(XP - XPD)
  448.  
  449.                   XZT = XP - XPD
  450.  
  451.  
  452.  
  453. C-----------------Endommagement critique                
  454.  
  455.                   IF (XD.GT.XDC) THEN
  456.  
  457.                      XD = XDC
  458.  
  459.                   ENDIF
  460.  
  461.  
  462.  
  463.                ELSE
  464.  
  465.                    GOTO 920
  466.  
  467.                ENDIF
  468.  
  469.  
  470.  
  471.                IF (CRIT.LE.1.0D-5) THEN
  472.  
  473.                    GOTO 920
  474.  
  475.                ENDIF
  476.  
  477.  
  478.  
  479.             ELSE
  480.  
  481.                GOTO 920
  482.  
  483.             ENDIF
  484.  
  485.  
  486.  
  487.          ENDDO
  488.  
  489.  
  490.  
  491. C-----Sortie de la boucle d'endommagement
  492.  
  493.  920  CONTINUE
  494.  
  495.  
  496.  
  497. C--------Rupture
  498.  
  499.          IF (XP.GE.XPR) THEN
  500.  
  501.             XNRUP = 1
  502.  
  503.             XP    = XPR
  504.  
  505.             DO I = 1,3
  506.  
  507.                DO J = 1,3
  508.  
  509.                   SIG(I, J) = 1.0D0
  510.  
  511.                ENDDO
  512.  
  513.             ENDDO
  514.  
  515.          ENDIF
  516.  
  517.  
  518.  
  519. C
  520.  
  521. C-----SOCKAGE EN SORTIE-------------------------------------------------
  522.  
  523. C
  524.  
  525. C--------Les variables internes
  526.  
  527.          VARF(1)  = XD
  528.  
  529.          VARF(2)  = XR
  530.  
  531.          VARF(3)  = XP
  532.  
  533.          VARF(4)  = XZT
  534.  
  535.          VARF(5)  = XNRUP
  536.  
  537.  
  538.  
  539.          VARF(6)  = XEPSP(1,1)
  540.  
  541.          VARF(7)  = XEPSP(2,2)
  542.  
  543.          VARF(8)  = XEPSP(3,3)
  544.  
  545.          VARF(9)  = XEPSP(1,2)*2.0D0
  546.  
  547.          VARF(10) = XEPSP(1,3)*2.0D0
  548.  
  549.          VARF(11) = XEPSP(2,3)*2.0D0
  550.  
  551.  
  552.  
  553.          VARF(12) = XEPSF(1,1)
  554.  
  555.          VARF(13) = XEPSF(2,2)
  556.  
  557.          VARF(14) = XEPSF(3,3)
  558.  
  559.          VARF(15) = XEPSF(1,2)*2.0D0
  560.  
  561.          VARF(16) = XEPSF(1,3)*2.0D0
  562.  
  563.          VARF(17) = XEPSF(2,3)*2.0D0
  564.  
  565.  
  566.  
  567. C--------Les contraintes
  568.  
  569.          SIGF(1) = SIG(1,1)
  570.  
  571.          SIGF(2) = SIG(2,2)
  572.  
  573.          SIGF(3) = SIG(3,3)
  574.  
  575.          SIGF(4) = SIG(1,2)
  576.  
  577.          SIGF(5) = SIG(1,3)
  578.  
  579.          SIGF(6) = SIG(2,3)
  580.  
  581.  
  582.  
  583. C-----------------------------------------------------------------------
  584.  
  585. C-----CALCUL DEFORMATION PLANE------------------------------------------
  586.  
  587. C -----------------------------------------------------------------------  
  588.  
  589.       ELSEIF (IFOUR.EQ.-1) THEN
  590.  
  591.  
  592.  
  593. C--------Deformation plastique      
  594.  
  595.          XEPSP(1,1) = VAR0(6)
  596.  
  597.          XEPSP(2,2) = VAR0(7)
  598.  
  599.          XEPSP(3,3) = VAR0(8)
  600.  
  601.          XEPSP(1,2) = VAR0(9)/2.0D0
  602.  
  603.          XEPSP(1,3) = 0.0D0
  604.  
  605.          XEPSP(2,3) = 0.0D0
  606.  
  607.  
  608.  
  609.          XEPSP(2,1) = XEPSP(1,2)
  610.  
  611.          XEPSP(2,3) = XEPSP(2,3)
  612.  
  613.          XEPSP(3,1) = XEPSP(1,3)
  614.  
  615.          XEPSP(3,2) = XEPSP(2,3)
  616.  
  617.  
  618.  
  619. C--------Deformation totale        
  620.  
  621.          XEPSF(1,1) = DEPST(1) + VAR0(10)
  622.  
  623.          XEPSF(2,2) = DEPST(2) + VAR0(11)
  624.  
  625.          XEPSF(3,3) = DEPST(3) + VAR0(12)
  626.  
  627.          XEPSF(1,2) = (DEPST(4) + VAR0(13))/2.0D0
  628.  
  629.          XEPSF(1,3) = 0.0D0
  630.  
  631.          XEPSF(2,3) = 0.0D0
  632.  
  633.  
  634.  
  635.          XEPSF(2,1) = XEPSF(1,2)
  636.  
  637.          XEPSF(2,3) = XEPSF(2,3)
  638.  
  639.          XEPSF(3,1) = XEPSF(1,3)
  640.  
  641.          XEPSF(3,2) = XEPSF(2,3)
  642.  
  643.  
  644.  
  645. C--------Deformation elastique        
  646.  
  647.          DO I = 1,3
  648.  
  649.             DO J = 1,3
  650.  
  651.                XEPSE(I,J) = XEPSF(I,J) - XEPSP(I,J)
  652.  
  653.             ENDDO
  654.  
  655.          ENDDO
  656.  
  657.  
  658.  
  659. C--------Deformation plane     
  660.  
  661.          XEPSE(1,3) = 0.0D0
  662.  
  663.          XEPSE(2,3) = 0.0D0
  664.  
  665.          XEPSE(3,1) = 0.0D0
  666.  
  667.          XEPSE(3,2) = 0.0D0
  668.  
  669.          XEPSE(3,3) = 0.0D0
  670.  
  671.  
  672.  
  673. C--------Prediction elastique
  674.  
  675.          TREPSE = XEPSE(1,1) + XEPSE(2,2) + XEPSE(3,3)
  676.  
  677.          DO I = 1,3
  678.  
  679.             DO J = 1,3
  680.  
  681.                AUX1 = XLA*TREPSE*UNIT3(I,J)
  682.  
  683.                AUX2 = 2.0D0*XMU*XEPSE(I,J)
  684.  
  685.                SIG(I,J) = (1.0D0 - XD)*(AUX1 + AUX2)
  686.  
  687.             ENDDO
  688.  
  689.          ENDDO
  690.  
  691.          CALL J2(SIG,SIGEQ)
  692.  
  693.          CALL DEVIAT(SIG,DVSIG)
  694.  
  695.  
  696.  
  697. C--------Boucle d'endommagement         
  698.  
  699.          DO I = 1,IMAX
  700.  
  701. C-----------Evaluation de la fonction seuil
  702.  
  703.             FP  = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
  704.  
  705.             FP0 = FP
  706.  
  707.  
  708.  
  709. C-----------Boucle de plasticite      
  710.  
  711.             IF (FP.GT.1.0D0) THEN  
  712.  
  713.  
  714.  
  715.                DO J = 1,IMAX
  716.  
  717. C-----------------Calcul des derives               
  718.  
  719.                   DO K = 1,3
  720.  
  721.                      DO L = 1,3
  722.  
  723.                         AUX1 = ((3.0D0/2.0D0)/(SIGEQ*(1.0D0 - XD)))
  724.  
  725.                         DFDSG(K,L) = AUX1*DVSIG(K,L)
  726.  
  727.                      ENDDO
  728.  
  729.                   ENDDO
  730.  
  731.  
  732.  
  733.                   DO K = 1,3
  734.  
  735.                      DO L = 1,3
  736.  
  737.                         DPHDSG(K,L) = DFDSG(K,L)
  738.  
  739.                      ENDDO
  740.  
  741.                   ENDDO
  742.  
  743.  
  744.  
  745. C-----------------Calcul du multiplicateur plastique
  746.  
  747.                   TRDPHDSG = DPHDSG(1,1) + DPHDSG(2,2) + DPHDSG(3,3)
  748.  
  749.                   DO K = 1,3
  750.  
  751.                      DO L = 1,3
  752.  
  753.                         AUX1 = XLA*TRDPHDSG*UNIT3(K,L)
  754.  
  755.                         AUX2 = 2.0D0*XMU*DPHDSG(K,L)
  756.  
  757.                         EDPH(K,L) = (AUX1 + AUX2)
  758.  
  759.                      ENDDO
  760.  
  761.                   ENDDO
  762.  
  763.                   DO K = 1,3
  764.  
  765.                      DFEDPHI(K) = 0.0D0
  766.  
  767.                      DO L = 1,3
  768.  
  769.                         AUX1 = DFDSG(K,L)*EDPH(L,K)
  770.  
  771.                         DFEDPHI(K) = DFEDPHI(K) + AUX1
  772.  
  773.                      ENDDO
  774.  
  775.                   ENDDO
  776.  
  777.                   AUX4  = DFEDPHI(1) + DFEDPHI(2) + DFEDPHI(3)
  778.  
  779.                   AUX5  = XK/XM
  780.  
  781.                   AUX6  = SIGEQ/(XK*(1.0D0 - XD)) 
  782.  
  783.                   AUX7  = XSIGY/XK
  784.  
  785.                   AUX8  = (1.0D0 - XM)           
  786.  
  787.                   XLAMBDA = FP/(AUX4 + AUX5*((AUX6 + AUX7)**AUX8))
  788.  
  789.  
  790.  
  791. C-----------------Mise a jour des variables
  792.  
  793.                   DO K = 1,3
  794.  
  795.                      DO L = 1,3
  796.  
  797.                         SIG(K,L) = SIG(K,L) - XLAMBDA*EDPH(K,L)
  798.  
  799.                      ENDDO
  800.  
  801.                   ENDDO
  802.  
  803.  
  804.  
  805.                   CALL J2(SIG,SIGEQ)
  806.  
  807.  
  808.  
  809.                   CALL DEVIAT(SIG,DVSIG)
  810.  
  811.  
  812.  
  813.                   DO K = 1,3
  814.  
  815.                      DO L = 1,3
  816.  
  817.                         XEPSP(K,L) = XEPSP(K,L) + XLAMBDA*DPHDSG(K,L)
  818.  
  819.                      ENDDO
  820.  
  821.                   ENDDO
  822.  
  823.                   XP   = XP + XLAMBDA/(1.0D0 - XD);
  824.  
  825.                   XR   = XK*(XP**(1.0D0/XM))
  826.  
  827.  
  828.  
  829. C-----------------Reevaluation de la fonction seuil
  830.  
  831.                   FP   = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
  832.  
  833.  
  834.  
  835.                   CRIT1 = DABS(FP/FP0)
  836.  
  837.                   CRIT2 = DABS(XLAMBDA)
  838.  
  839.  
  840.  
  841.                   IF ((CRIT1.LE.1.0D-10).OR.(CRIT2.LE.1.0D-10)) THEN 
  842.  
  843.                       GOTO 930
  844.  
  845.                   ENDIF
  846.  
  847.  
  848.  
  849.                ENDDO
  850.  
  851.  
  852.  
  853. C-----Sortie de la boucle de plasticité
  854.  
  855.  930  CONTINUE
  856.  
  857.  
  858.  
  859. C--------------Evaluation du seuil d'endommagement
  860.  
  861.                FD = XP - (XPD + XZT)
  862.  
  863.                IF (I.EQ.1) THEN
  864.  
  865.                    FD0 = FD
  866.  
  867.                ENDIF
  868.  
  869.  
  870.  
  871.                CRIT = DABS(FD/FD0)
  872.  
  873.  
  874.  
  875.                IF ((FD.GT.0.0D0).AND.(CRIT.GT.1.0D-5)) THEN
  876.  
  877.  
  878.  
  879. C-----------------Calcul de l'endommagement
  880.  
  881.                   XD  = (XDC/(XPR - XPD))*(XP - XPD)
  882.  
  883.                   XZT = XP - XPD
  884.  
  885.  
  886.  
  887. C-----------------Endommagement critique                
  888.  
  889.                   IF (XD.GT.XDC) THEN
  890.  
  891.                      XD = XDC
  892.  
  893.                   ENDIF
  894.  
  895.  
  896.  
  897.                ELSE
  898.  
  899.                    GOTO 940
  900.  
  901.                ENDIF
  902.  
  903.  
  904.  
  905.                IF (CRIT.LE.1.0D-5) THEN
  906.  
  907.                    GOTO 940
  908.  
  909.                ENDIF
  910.  
  911.  
  912.  
  913.             ELSE
  914.  
  915.                GOTO 940
  916.  
  917.             ENDIF
  918.  
  919.  
  920.  
  921.          ENDDO
  922.  
  923.  
  924.  
  925. C-----Sortie de la boucle d'endommagement
  926.  
  927.  940  CONTINUE
  928.  
  929.  
  930.  
  931. C-----Rupture
  932.  
  933.       IF (XP.GE.XPR) THEN
  934.  
  935.          XNRUP = 1
  936.  
  937.          XP    = XPR
  938.  
  939.          DO I = 1,3
  940.  
  941.             DO J = 1,3
  942.  
  943.                SIG(I, J) = 1.0D0
  944.  
  945.             ENDDO
  946.  
  947.          ENDDO
  948.  
  949.       ENDIF
  950.  
  951.  
  952.  
  953. C
  954.  
  955. C-----SOCKAGE EN SORTIE-------------------------------------------------
  956.  
  957. C
  958.  
  959. C--------Les variables internes
  960.  
  961.          VARF(1)  = XD
  962.  
  963.          VARF(2)  = XR
  964.  
  965.          VARF(3)  = XP
  966.  
  967.          VARF(4)  = XZT
  968.  
  969.          VARF(5)  = XNRUP
  970.  
  971.  
  972.  
  973.          VARF(6)  = XEPSP(1,1)
  974.  
  975.          VARF(7)  = XEPSP(2,2)
  976.  
  977.          VARF(8)  = XEPSP(3,3)
  978.  
  979.          VARF(9)  = XEPSP(1,2)*2.0D0
  980.  
  981.  
  982.  
  983.          VARF(10) = XEPSF(1,1)
  984.  
  985.          VARF(11) = XEPSF(2,2)
  986.  
  987.          VARF(12) = XEPSF(3,3)
  988.  
  989.          VARF(13) = XEPSF(1,2)*2.0D0
  990.  
  991.  
  992.  
  993. C--------Les contraintes
  994.  
  995.          SIGF(1) = SIG(1,1)
  996.  
  997.          SIGF(2) = SIG(2,2)
  998.  
  999.          SIGF(3) = SIG(3,3)
  1000.  
  1001.          SIGF(4) = SIG(1,2)
  1002.  
  1003.  
  1004.  
  1005. C-----------------------------------------------------------------------
  1006.  
  1007. C-----CALCUL CONTRAINTES PLANES-----------------------------------------
  1008.  
  1009. C----------------------------------------------------------------------- 
  1010.  
  1011.       ELSEIF (IFOUR.EQ.-2) THEN
  1012.  
  1013.  
  1014.  
  1015. C--------Deformation plastique      
  1016.  
  1017.          XEPSP(1,1) = VAR0(6)
  1018.  
  1019.          XEPSP(2,2) = VAR0(7)
  1020.  
  1021.          XEPSP(3,3) = VAR0(8)
  1022.  
  1023.          XEPSP(1,2) = VAR0(9)/2.0D0
  1024.  
  1025.          XEPSP(1,3) = 0.0D0
  1026.  
  1027.          XEPSP(2,3) = 0.0D0
  1028.  
  1029.  
  1030.  
  1031.          XEPSP(2,1) = XEPSP(1,2)
  1032.  
  1033.          XEPSP(2,3) = XEPSP(2,3)
  1034.  
  1035.          XEPSP(3,1) = XEPSP(1,3)
  1036.  
  1037.          XEPSP(3,2) = XEPSP(2,3)
  1038.  
  1039.  
  1040.  
  1041. C--------Deformation totale        
  1042.  
  1043.          XEPSF(1,1) = DEPST(1) + VAR0(10)
  1044.  
  1045.          XEPSF(2,2) = DEPST(2) + VAR0(11)
  1046.  
  1047.          XEPSF(3,3) = DEPST(3) + VAR0(12)
  1048.  
  1049.          XEPSF(1,2) = (DEPST(4) + VAR0(13))/2.0D0
  1050.  
  1051.          XEPSF(1,3) = 0.0D0
  1052.  
  1053.          XEPSF(2,3) = 0.0D0
  1054.  
  1055.  
  1056.  
  1057.          XEPSF(2,1) = XEPSF(1,2)
  1058.  
  1059.          XEPSF(2,3) = XEPSF(2,3)
  1060.  
  1061.          XEPSF(3,1) = XEPSF(1,3)
  1062.  
  1063.          XEPSF(3,2) = XEPSF(2,3)
  1064.  
  1065.  
  1066.  
  1067. C-------Modification de la déformation 33
  1068.  
  1069.          XEPSF(3,3) = VAR0(14)
  1070.  
  1071.  
  1072.  
  1073. C--------Boucle contrainte plane
  1074.  
  1075.          DO M = 1,IMAX
  1076.  
  1077.  
  1078.  
  1079. C-----------Deformation elastique        
  1080.  
  1081.             DO I = 1,3
  1082.  
  1083.                DO J = 1,3
  1084.  
  1085.                   XEPSE(I,J) = XEPSF(I,J) - XEPSP(I,J)
  1086.  
  1087.                ENDDO
  1088.  
  1089.             ENDDO
  1090.  
  1091.  
  1092.  
  1093. C-----------Prediction elastique
  1094.  
  1095.             TREPSE = XEPSE(1,1) + XEPSE(2,2) + XEPSE(3,3)
  1096.  
  1097.             DO I = 1,3
  1098.  
  1099.                DO J = 1,3
  1100.  
  1101.                   AUX1 = XLA*TREPSE*UNIT3(I,J)
  1102.  
  1103.                   AUX2 = 2.0D0*XMU*XEPSE(I,J)
  1104.  
  1105.                   SIG(I,J) = (1.0D0 - XD)*(AUX1 + AUX2)
  1106.  
  1107.                ENDDO
  1108.  
  1109.             ENDDO
  1110.  
  1111.             CALL J2(SIG,SIGEQ)
  1112.  
  1113.             CALL DEVIAT(SIG,DVSIG)
  1114.  
  1115.  
  1116.  
  1117. C-----------Boucle d'endommagement
  1118.  
  1119.             DO I = 1,IMAX          
  1120.  
  1121. C--------------Evaluation de la fonction seuil
  1122.  
  1123.                FP  = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
  1124.  
  1125.                FP0 = FP
  1126.  
  1127.  
  1128.  
  1129. C--------------Boucle de plasticite      
  1130.  
  1131.                IF (FP.GT.1.0D0) THEN 
  1132.  
  1133.  
  1134.  
  1135.                   DO J = 1,IMAX              
  1136.  
  1137. C--------------------Calcul des derives               
  1138.  
  1139.                      DO K = 1,3
  1140.  
  1141.                         DO L = 1,3
  1142.  
  1143.                            AUX1 = ((3.0D0/2.0D0)/(SIGEQ*(1.0D0 - XD)))
  1144.  
  1145.                            DFDSG(K,L) = AUX1*DVSIG(K,L)
  1146.  
  1147.                         ENDDO
  1148.  
  1149.                      ENDDO
  1150.  
  1151.  
  1152.  
  1153.                      DO K = 1,3
  1154.  
  1155.                         DO L = 1,3
  1156.  
  1157.                            DPHDSG(K,L) = DFDSG(K,L)
  1158.  
  1159.                         ENDDO
  1160.  
  1161.                      ENDDO
  1162.  
  1163.  
  1164.  
  1165. C--------------------Calcul du multiplicateur plastique
  1166.  
  1167.                      TRDPHDSG = DPHDSG(1,1) + DPHDSG(2,2) + DPHDSG(3,3)
  1168.  
  1169.                      DO K = 1,3
  1170.  
  1171.                         DO L = 1,3
  1172.  
  1173.                            AUX1 = XLA*TRDPHDSG*UNIT3(K,L)
  1174.  
  1175.                            AUX2 = 2.0D0*XMU*DPHDSG(K,L)
  1176.  
  1177.                            EDPH(K,L) = (AUX1 + AUX2)
  1178.  
  1179.                         ENDDO
  1180.  
  1181.                      ENDDO
  1182.  
  1183.                      DO K = 1,3
  1184.  
  1185.                         DFEDPHI(K) = 0.0D0
  1186.  
  1187.                         DO L = 1,3
  1188.  
  1189.                            AUX1 = DFDSG(K,L)*EDPH(L,K)
  1190.  
  1191.                            DFEDPHI(K) = DFEDPHI(K) + AUX1
  1192.  
  1193.                         ENDDO
  1194.  
  1195.                      ENDDO
  1196.  
  1197.                      AUX4  = DFEDPHI(1) + DFEDPHI(2) + DFEDPHI(3)
  1198.  
  1199.                      AUX5  = XK/XM
  1200.  
  1201.                      AUX6  = SIGEQ/(XK*(1.0D0 - XD)) 
  1202.  
  1203.                      AUX7  = XSIGY/XK
  1204.  
  1205.                      AUX8  = (1.0D0 - XM)           
  1206.  
  1207.                      XLAMBDA = FP/(AUX4 + AUX5*((AUX6 + AUX7)**AUX8))
  1208.  
  1209.  
  1210.  
  1211. C--------------------Mise a jour des variables
  1212.  
  1213.                      DO K = 1,3
  1214.  
  1215.                         DO L = 1,3
  1216.  
  1217.                            SIG(K,L) = SIG(K,L) - XLAMBDA*EDPH(K,L)
  1218.  
  1219.                         ENDDO
  1220.  
  1221.                      ENDDO
  1222.  
  1223.  
  1224.  
  1225.                      CALL J2(SIG,SIGEQ)
  1226.  
  1227.  
  1228.  
  1229.                      CALL DEVIAT(SIG,DVSIG)
  1230.  
  1231.  
  1232.  
  1233.                      DO K = 1,3
  1234.  
  1235.                         DO L = 1,3
  1236.  
  1237.                            XEPSP(K,L) = XEPSP(K,L) + XLAMBDA*DPHDSG(K,L)
  1238.  
  1239.                         ENDDO
  1240.  
  1241.                      ENDDO
  1242.  
  1243.                      XP   = XP + XLAMBDA/(1.0D0 - XD);
  1244.  
  1245.                      XR   = XK*(XP**(1.0D0/XM))
  1246.  
  1247.  
  1248.  
  1249. C--------------------Reevaluation de la fonction seuil
  1250.  
  1251.                      FP   = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY)
  1252.  
  1253.  
  1254.  
  1255.                      CRIT1 = DABS(FP/FP0)
  1256.  
  1257.                      CRIT2 = DABS(XLAMBDA)
  1258.  
  1259.  
  1260.  
  1261.                      IF ((CRIT1.LE.1.0D-10).OR.(CRIT2.LE.1.0D-10)) THEN 
  1262.  
  1263.                          GOTO 950
  1264.  
  1265.                      ENDIF
  1266.  
  1267.  
  1268.  
  1269.                   ENDDO
  1270.  
  1271.  
  1272.  
  1273. C-----Sortie de la boucle de plasticité
  1274.  
  1275.  950  CONTINUE
  1276.  
  1277.  
  1278.  
  1279. C-----------------Evaluation du seuil d'endommagement
  1280.  
  1281.                   FD = XP - (XPD + XZT)
  1282.  
  1283.                   IF (I.EQ.1) THEN
  1284.  
  1285.                       FD0 = FD
  1286.  
  1287.                   ENDIF
  1288.  
  1289.  
  1290.  
  1291.                   CRIT = DABS(FD/FD0)
  1292.  
  1293.  
  1294.  
  1295.                   IF ((FD.GT.0.0D0).AND.(CRIT.GT.1.0D-5)) THEN
  1296.  
  1297.  
  1298.  
  1299. C--------------------Calcul de l'endommagement
  1300.  
  1301.                      XD  = (XDC/(XPR - XPD))*(XP - XPD)
  1302.  
  1303.                      XZT = XP - XPD
  1304.  
  1305.  
  1306.  
  1307. C--------------------Endommagement critique                
  1308.  
  1309.                      IF (XD.GT.XDC) THEN
  1310.  
  1311.                         XD = XDC
  1312.  
  1313.                      ENDIF
  1314.  
  1315.  
  1316.  
  1317.                   ELSE
  1318.  
  1319.                       GOTO 960
  1320.  
  1321.                   ENDIF
  1322.  
  1323.  
  1324.  
  1325.                   IF (CRIT.LE.1.0D-5) THEN
  1326.  
  1327.                       GOTO 960
  1328.  
  1329.                   ENDIF
  1330.  
  1331.  
  1332.  
  1333.                ELSE
  1334.  
  1335.                   GOTO 960
  1336.  
  1337.                ENDIF
  1338.  
  1339.  
  1340.  
  1341.             ENDDO
  1342.  
  1343.  
  1344.  
  1345. C-----Sortie de la boucle d'endommagement
  1346.  
  1347.  960  CONTINUE
  1348.  
  1349.  
  1350.  
  1351. C-----------Contraintes Planes 
  1352.  
  1353.             CRIT = SIG(3,3)/SIG(2,2)
  1354.  
  1355.  
  1356.  
  1357.             IF (CRIT.GT.1.0D-8) THEN
  1358.  
  1359.                XEPSF(3,3) = XEPSF(3,3) - SIG(3,3)/XYG
  1360.  
  1361.             ELSE
  1362.  
  1363.                GOTO 970
  1364.  
  1365.             ENDIF
  1366.  
  1367.  
  1368.  
  1369.          ENDDO
  1370.  
  1371.  
  1372.  
  1373. C-----Sortie de la boucle de contraintes planes
  1374.  
  1375.  970  CONTINUE
  1376.  
  1377.  
  1378.  
  1379. C-----Rupture
  1380.  
  1381.       IF (XP.GE.XPR) THEN
  1382.  
  1383.          XNRUP = 1
  1384.  
  1385.          XP    = XPR
  1386.  
  1387.          DO I = 1,3
  1388.  
  1389.             DO J = 1,3
  1390.  
  1391.                SIG(I, J) = 1.0D0
  1392.  
  1393.             ENDDO
  1394.  
  1395.          ENDDO
  1396.  
  1397.       ENDIF
  1398.  
  1399.  
  1400.  
  1401. C
  1402.  
  1403. C-----SOCKAGE EN SORTIE-------------------------------------------------
  1404.  
  1405. C
  1406.  
  1407. C--------Les variables internes
  1408.  
  1409.          VARF(1)  = XD
  1410.  
  1411.          VARF(2)  = XR
  1412.  
  1413.          VARF(3)  = XP
  1414.  
  1415.          VARF(4)  = XZT
  1416.  
  1417.          VARF(5)  = XNRUP
  1418.  
  1419.  
  1420.  
  1421.          VARF(6)  = XEPSP(1,1)
  1422.  
  1423.          VARF(7)  = XEPSP(2,2)
  1424.  
  1425.          VARF(8)  = XEPSP(3,3)
  1426.  
  1427.          VARF(9)  = XEPSP(1,2)*2.0D0
  1428.  
  1429.  
  1430.  
  1431.          VARF(10) = XEPSF(1,1)
  1432.  
  1433.          VARF(11) = XEPSF(2,2)
  1434.  
  1435.          VARF(12) = XEPSF(3,3)
  1436.  
  1437.          VARF(13) = XEPSF(1,2)*2.0D0
  1438.  
  1439.  
  1440.  
  1441.          VARF(14) = XEPSF(3,3)
  1442.  
  1443.  
  1444.  
  1445. C--------Les contraintes
  1446.  
  1447.          SIGF(1) = SIG(1,1)
  1448.  
  1449.          SIGF(2) = SIG(2,2)
  1450.  
  1451.          SIGF(3) = SIG(3,3)
  1452.  
  1453.          SIGF(4) = SIG(1,2)
  1454.  
  1455.  
  1456.  
  1457.       ENDIF
  1458.  
  1459.  
  1460.  
  1461. C-----Fin de l integration
  1462.  
  1463.       RETURN
  1464.  
  1465.       END
  1466.  
  1467.  
  1468.  
  1469.  
  1470.  
  1471.  
  1472.  
  1473.  
  1474.  
  1475.  
  1476.  
  1477.  
  1478.  
  1479.  
  1480.  
  1481.  
  1482.  
  1483.  
  1484.  
  1485.  
  1486.  
  1487.  
  1488.  
  1489.  
  1490.  
  1491.  
  1492.  
  1493.  
  1494.  
  1495.  
  1496.  
  1497.  
  1498.  
  1499.  
  1500.  
  1501.  
  1502.  
  1503.  
  1504.  
  1505.  
  1506.  
  1507.  
  1508.  
  1509.  
  1510.  

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