Télécharger la plaquette Cast3M


fr en

Télécharger dp_sol.eso

Retour à la liste

Numérotation des lignes : 

dp_sol

  1. C DP_SOL    SOURCE    CB215821  17/11/30    21:16:00     9639           
  2. C DP_SOIL   SOURCE    AF221230  12/10/04    21:15:18     7520
  3.       SUBROUTINE DP_SOL (XMAT,VAR0,VARF,SIG0,SIGF,DEPST,XCAR)
  4. C
  5. C====&===1=========2=========3=========4=========5=========6=========7==
  6. C Commentaires :   Loi Drucker Prager avec ecrouissage non lineaire et
  7. C                  loi non associée
  8. C Traits       : -
  9. C
  10. C Auteur       : A. Frau (Dr - Ing.) - CEA/DEN/DANS/DM2S/SEMT/EMSI
  11. C====&===1=========2=========3=========4=========5=========6=========7==
  12. C
  13. C----DECLARATION GENERALES----------------------------------------------
  14. C
  15.       IMPLICIT INTEGER(I-N)
  16.       IMPLICIT REAL*8(A-H,O-Z)
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. C
  21. C----DECLARATION SEGMENTS-----------------------------------------------
  22. C
  23.       REAL*8 XMAT(*),SIG0(*),SIGF(*),VAR0(*),VARF(*),DEPST(*),XCAR(*)
  24.       REAL*8 XNU,XYOUN,XSY
  25.       REAL*8 XMU,XLAM,XALPHA,XSIGI,XBET1,XGAM1,XDEL1
  26.       REAL*8 XDEEPP(6),XEPP0(6)
  27.       REAL*8 XSIGT(6),XCHIT,XQT,XFHIT,XDSIGT(6),XNDSIGT
  28.       REAL*8 XNUT(6),XPIT(6)
  29.       REAL*8 XTOL1,XVPL1
  30.       REAL*8 XGAMN,XSIGN(6),XEPPN(6),XCHIN,XQN,XFHIN,XDSIGN(6),XNDSIGN
  31.       REAL*8 XGAMN1,XSIGN1(6),XEPPN1(6),XCHIN1,XQN1,XFHIN1
  32.       REAL*8 XDFDG,XDEGAM,XDSIGN1(6),XNDSIGN1,XTSIGN1
  33.       REAL*8 XGAMF,XCHIF,XQF,XEPPF(6),XVALT1,XVALT2
  34. C
  35. C----CONTROLE OPTION CALCUL --------------------------------------------
  36.       IF ((IFOUR.NE.2).AND.(IFOUR.NE.0).AND.(IFOUR.NE.-1)) THEN
  37.         WRITE(IOIMP,10000)
  38. 10000   FORMAT('Option calcul mauvais')
  39.         WRITE(IOIMP,10001)
  40. 10001   FORMAT('  seulement:')
  41.         WRITE(IOIMP,10002)
  42. 10002   FORMAT('    3D MODE PLAN')
  43.         WRITE(IOIMP,10003)
  44. 10003   FORMAT('    2D MODE AXIS')
  45.         WRITE(IOIMP,10004)
  46. 10004   FORMAT('    2D MODE PLAN')
  47.         STOP
  48.       ENDIF
  49. C
  50. C----MISE EN DONNEES----------------------------------------------------
  51. C
  52. C     Module de young
  53.       XYOUN = XMAT(1)
  54. C     Coef Poisson
  55.       XNU  = XMAT(2)
  56. C     alpha
  57.       XALPHA  =  XMAT(5)
  58. C     Sigma_y
  59.       XSY  =  XMAT(6)
  60. C     Sig_inf
  61.       XSIGI  =  XMAT(7)
  62. C     Beta
  63.       XBET1  =  XMAT(8)
  64. C     Gamma
  65.       XGAM1  =  XMAT(9)
  66. C     Delta
  67.       XDEL1  =  XMAT(10)
  68. C
  69. C----VARIABLES INTERNES-------------------------------------------------
  70.  
  71. c      Coef mu elastique
  72.        XMU = (0.5D0)*((XYOUN)/(1.0D0 + XNU))
  73.  
  74. c      Coef lambda elastique
  75.        XLAM = ((XNU)/(1.0D0-((2.0D0)*(XNU))))
  76. C
  77. C----ACTUALISATION DEFORMATION------------------------------------------
  78. C
  79. C
  80. C----PREDICTION ELASTIQUE-----------------------------------------------
  81. C
  82. c      gamma_test = 0 donc le delta_plastique = 0
  83. c      Initialissation Eps_pla_initiales
  84.        DO 11 I1=1,6
  85.          XDEEPP(I1) = 0.0D0
  86.          XEPP0(I1) = VAR0(I1)
  87. 11     CONTINUE
  88. c
  89. c      Calcul de la variable q_test = q_ini -> gamma_test=0
  90.        XCHIT = VAR0(7)
  91. c
  92. c      Calcul de la variable q_test = q_ini -> gamma_test=0
  93.        XQT = ((-1.0D0)*(XSIGI - XSY))*(1.0D0
  94.      &     - EXP(((-1.0D0)*(XBET1))*(XCHIT)))
  95. c
  96. c      Calcul Sigma_test a partir du step initiale
  97.        DO 12 I1=1,6
  98.          IF (I1.LE.3) THEN
  99.            XSIGT(I1) = SIG0(I1) +
  100.      &     ((2.0D0)*(XMU))*(DEPST(I1) + ((XLAM)*(TRACE(DEPST))))
  101.          ELSE
  102.            IF (IFOUR.EQ.2) THEN
  103.              XSIGT(I1) = SIG0(I1) + ((1.0D0)*(XMU))*(DEPST(I1))
  104.            ELSE
  105.              IF (I1.EQ.4) THEN
  106.                XSIGT(I1) = SIG0(I1) + ((1.0D0)*(XMU))*(DEPST(I1))
  107.              ELSE
  108.                XSIGT(I1) = 0.0D0
  109.              ENDIF
  110.            ENDIF
  111.          ENDIF
  112. 12     CONTINUE
  113. C
  114. C      Deviatoire de Sigma_test
  115.        DO 13 I1=1,6
  116.          IF (I1.LE.3) THEN
  117.            XDSIGT(I1) = XSIGT(I1) - ((TRACE(XSIGT)/(3.0D0)))
  118.          ELSE
  119.            IF (IFOUR.EQ.2) THEN
  120.              XDSIGT(I1) = XSIGT(I1)
  121.            ELSE
  122.              IF (I1.EQ.4) THEN
  123.                XDSIGT(I1) = XSIGT(I1)
  124.              ELSE
  125.                XDSIGT(I1) = 0.0D0
  126.              ENDIF
  127.            ENDIF
  128.          ENDIF
  129. 13     CONTINUE
  130. c
  131. c      Calcul de la norme partie deviatoire
  132.        XNDSIGT = 0.0D0
  133.        DO 14 I1=1,6
  134.          IF (I1.LE.3) THEN
  135.            XNDSIGT = XNDSIGT + (XDSIGT(I1)*XDSIGT(I1))
  136.          ELSE
  137.            XNDSIGT = XNDSIGT + ((2.0D0)*(XDSIGT(I1)*XDSIGT(I1)))
  138.          ENDIF
  139. 14     CONTINUE
  140.        XNDSIGT = ((XNDSIGT)**(0.5D0))
  141. C
  142. c      Calcul du tensuer pi
  143.        XPIT(1) = 1.0D0
  144.        XPIT(2) = 1.0D0
  145.        XPIT(3) = 1.0D0
  146.        XPIT(4) = 0.0D0
  147.        XPIT(5) = 0.0D0
  148.        XPIT(6) = 0.0D0
  149.        XPIT = (((3.0D0)**(0.5D0))/(3.0D0))*(XPIT)
  150. C
  151. c      Calcul ((3**0.5)*(trace(sigma_test)))
  152.        XTSIGT = 0.0D0
  153.        DO 23 I1=1,6
  154.          XTSIGT = XTSIGT + (XPIT(I1)*(XSIGT(I1)))
  155. 23     CONTINUE
  156. c
  157. c      Calcul du vecteur nu_test
  158.        IF ((ABS(XNDSIGT/XTSIGT)).LT.(1.0E-8)) THEN
  159.          XNUT(1) = (1.0D0)/(3.0D0)
  160.          XNUT(2) = (1.0D0)/(3.0D0)
  161.          XNUT(3) = (-2.0D0)/(3.0D0)
  162.          XNUT(4) = 0.0D0
  163.          XNUT(5) = 0.0D0
  164.          XNUT(6) = 0.0D0
  165.          XVPL1 = 0.0D0
  166.        ELSE
  167.          XNUT(1) = XDSIGT(1)/XNDSIGT
  168.          XNUT(2) = XDSIGT(2)/XNDSIGT
  169.          XNUT(3) = XDSIGT(3)/XNDSIGT
  170.          XNUT(4) = XDSIGT(4)/XNDSIGT
  171.          XNUT(5) = XDSIGT(5)/XNDSIGT
  172.          XNUT(6) = XDSIGT(6)/XNDSIGT
  173.          XVPL1 = 1.0D0
  174.        ENDIF
  175. c      Calcul du critere fhi_test
  176.        XFHIT = XNDSIGT + (((2.0D0)**(0.5))*(XTSIGT))*(XALPHA)
  177.      &  - (((2.0D0/3.0D0)**(0.5D0))*(XSY - XQT))
  178. C
  179. C-----TEST---------------------------------------------
  180. C
  181.        IF (XFHIT.LT.0) THEN
  182.          GOTO 101
  183. c        fhi_test < 0.0D0  pas de plasticité
  184.        ELSE IF (XFHIT.GE.0) THEN
  185.          GOTO 102
  186. c        fhi_test >= 0.0D0  on est en plasticité
  187.        END IF
  188. C
  189. C-----CALCUL ELASTIQUE---------------------------------------------
  190. C
  191.  
  192. 101   CONTINUE
  193. c
  194. c     si fhi_test < 0  gamma=0
  195.        DO 15 I1=1,6
  196.          SIGF(I1) = XSIGT(I1)
  197. 15     CONTINUE
  198. C
  199. c     pas des variations des variables internes
  200.        DO 16 I1=1,8
  201.          IF (I1.EQ.8) THEN
  202.            VARF(I1) = 0.0D0
  203.          ELSE
  204.            VARF(I1) = VAR0(I1)
  205.          ENDIF
  206. 16     CONTINUE
  207.       RETURN
  208. C
  209. C-----CALCUL PLASTIQUE---------------------------------------------
  210. C
  211. 102   CONTINUE
  212. c
  213. c        si fhi_test >= 0  determination de gamma
  214. c        Nombre max interaction
  215.          NINT1 = 10000
  216. C        Tollerance
  217.          XTOL1 = 1.E-8
  218.       XFHIN = XFHIT
  219.       XQN = XQT
  220.       XGAMN = 0.0D0
  221.       XCHIN = XCHIT
  222. c
  223.       DO 17 I1=1,6
  224.         XEPPN(I1) = XEPP0(I1)
  225.         XSIGN(I1) = XSIGT(I1)
  226. 17    CONTINUE
  227. c
  228.       DO 110 II1=1,NINT1
  229.         IF (II1.EQ.NINT1) THEN
  230. c          PRINT*,'!!!!!   ATTENTION   !!!!'
  231. c          PRINT*,'Nombre des iteractions'
  232. c          PRINT*,'       maximales depasse'
  233. c          STOP
  234.         ENDIF
  235.         XDFDG = (((2.0D0)*(XMU*XVPL1))*(XDEL1))
  236.      &  + ((2.0D0)/(3.0D0))*(((XSIGI - XSY)*(
  237.      &  EXP(((-1.0D0)*(XBET1))*(XCHIN))))*(XBET1)) +
  238.      &  ((4.0D0)*((XALPHA)*(XGAM1)))*
  239.      &  (1.0D0 + ((3.0D0)*(XLAM)))*(XMU)
  240.         XDEGAM = XFHIN/XDFDG
  241.         XGAMN1 = XGAMN + XDEGAM
  242.         XCHIN1 = XCHIN + (((2.0D0/3.0D0)**(0.5D0))*(XDEGAM))
  243.         XQN1 = ((-1.0D0)*(XSIGI - XSY))*(1.0D0
  244.      &       - EXP(((-1.0D0)*(XBET1))*(XCHIN1)))
  245.         IF ((XGAMN1).LT.(0.0D0)) THEN
  246.            WRITE(IOIMP,10005)
  247. 10005      FORMAT('Error....multiplicateur negative')
  248.            STOP
  249.         ENDIF
  250.  
  251.         DO 18 I1=1,6
  252.          IF (I1.LE.3) THEN
  253.             XDEEPP(I1) = ((((XDEL1*XVPL1))*(XNUT(I1)))+
  254.      &       (((XGAM1))*((2.0D0)**(0.5D0)))*(XPIT(I1)))*XDEGAM
  255.          ELSE
  256.            IF (IFOUR.EQ.2) THEN
  257.               XDEEPP(I1) = ((((XDEL1*XVPL1))*(XNUT(I1)))+
  258.      &         (((XGAM1))*((2.0D0)**(0.5D0)))*(XPIT(I1)))*XDEGAM
  259.            ELSE
  260.              IF (I1.EQ.4) THEN
  261.                XDEEPP(I1) = ((((XDEL1*XVPL1))*(XNUT(I1)))+
  262.      &          (((XGAM1))*((2.0D0)**(0.5D0)))*(XPIT(I1)))*XDEGAM
  263.              ELSE
  264.                XDEEPP(I1) = 0.0D0
  265.              ENDIF
  266.            ENDIF
  267.          ENDIF
  268.  
  269.          IF (I1.LE.3) THEN
  270.             XEPPN1(I1) = XEPPN(I1) + XDEEPP(I1)
  271.          ELSE
  272.            IF (IFOUR.EQ.2) THEN
  273.               XEPPN1(I1) = XEPPN(I1) + XDEEPP(I1)
  274.            ELSE
  275.              IF (I1.EQ.4) THEN
  276.                 XEPPN1(I1) = XEPPN(I1) + XDEEPP(I1)
  277.              ELSE
  278.                 XEPPN1(I1) = 0.0D0
  279.              ENDIF
  280.            ENDIF
  281.          ENDIF
  282.  
  283. 18      CONTINUE
  284.  
  285.         DO 21 I1=1,6
  286.         IF (I1.LE.3) THEN
  287.             XSIGN1(I1) = XSIGN(I1) -
  288.      &    ((2.0D0)*(XMU))*(XDEEPP(I1)
  289.      &    + ((XLAM)*(TRACE(XDEEPP))))
  290.           ELSE
  291.            IF (IFOUR.EQ.2) THEN
  292.             XSIGN1(I1) = XSIGN(I1) -
  293.      &    ((1.0D0)*(XMU))*(XDEEPP(I1))
  294.            ELSE
  295.              IF (I1.EQ.4) THEN
  296.                XSIGN1(I1) = XSIGN(I1) -
  297.      &       ((1.0D0)*(XMU))*(XDEEPP(I1))
  298.              ELSE
  299.                XSIGN1(I1) = 0.0D0
  300.              ENDIF
  301.            ENDIF
  302.           ENDIF
  303. 21      CONTINUE
  304.  
  305.         DO 19 I1=1,6
  306.           IF (I1.LE.3) THEN
  307.             XDSIGN1(I1) = XSIGN1(I1)
  308.      &       - ((TRACE(XSIGN1)/(3.0D0)))
  309.           ELSE
  310.            IF (IFOUR.EQ.2) THEN
  311.             XDSIGN1(I1) = XSIGN1(I1)
  312.            ELSE
  313.              IF (I1.EQ.4) THEN
  314.                XDSIGN1(I1) = XSIGN1(I1)
  315.              ELSE
  316.                XDSIGN1(I1) = 0.0D0
  317.              ENDIF
  318.            ENDIF
  319.           ENDIF
  320. 19      CONTINUE
  321.  
  322.         XNDSIGN1 = 0.0D0
  323.         DO 20 I1=1,6
  324.  
  325.          IF (I1.LE.3) THEN
  326.            XNDSIGN1 = XNDSIGN1
  327.      &        + (XDSIGN1(I1)*XDSIGN1(I1))
  328.          ELSE
  329.            XNDSIGN1 = XNDSIGN1
  330.      &        + ((2.0D0)*(XDSIGN1(I1)*XDSIGN1(I1)))
  331.          ENDIF
  332. 20      CONTINUE
  333.         XNDSIGN1 = ((XNDSIGN1)**(0.5D0))
  334.  
  335.         XTSIGN1 = XSIGN1(1) + XSIGN1(2) + XSIGN1(3)
  336.         XTSIGN1 = ((XTSIGN1)/((3.0D0)**(0.5D0)))
  337.         XFHIN1 = XNDSIGN1 +
  338.      &     (((2.0D0)**(0.5))*(XTSIGN1))*(XALPHA)
  339.      &   - (((2.0D0/3.0D0)**(0.5D0))*(XSY - XQN1))
  340. c
  341.         IF (ABS((XFHIN1)/(XFHIT)).LT.(XTOL1)) THEN
  342.           GOTO 103
  343.         ELSE
  344.           XFHIN = XFHIN1
  345.           XGAMN = XGAMN1
  346.           XCHIN = XCHIN1
  347.           XQN = XQN1
  348.           XEPPN(1) = XEPPN1(1)
  349.           XEPPN(2) = XEPPN1(2)
  350.           XEPPN(3) = XEPPN1(3)
  351.           XEPPN(4) = XEPPN1(4)
  352.           XEPPN(5) = XEPPN1(5)
  353.           XEPPN(6) = XEPPN1(6)
  354.           XSIGN(1) = XSIGN1(1)
  355.           XSIGN(2) = XSIGN1(2)
  356.           XSIGN(3) = XSIGN1(3)
  357.           XSIGN(4) = XSIGN1(4)
  358.           XSIGN(5) = XSIGN1(5)
  359.           XSIGN(6) = XSIGN1(6)
  360.         ENDIF
  361. 110   CONTINUE
  362.  
  363.  
  364. 103   CONTINUE
  365.       XCHIF = XCHIN1
  366.       XGAMF = XGAMN1
  367.       SIGF(1) = XSIGN1(1)
  368.       SIGF(2) = XSIGN1(2)
  369.       SIGF(3) = XSIGN1(3)
  370.       SIGF(4) = XSIGN1(4)
  371.       SIGF(5) = XSIGN1(5)
  372.       SIGF(6) = XSIGN1(6)
  373.  
  374.       VARF(1) = XEPPN1(1)
  375.       VARF(2) = XEPPN1(2)
  376.       VARF(3) = XEPPN1(3)
  377.       VARF(4) = XEPPN1(4)
  378.       VARF(5) = XEPPN1(5)
  379.       VARF(6) = XEPPN1(6)
  380.       VARF(7) = XCHIF
  381.       VARF(8) = XGAMF
  382.  
  383.       RETURN
  384.  
  385.       END
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  

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