Télécharger la plaquette Cast3M


fr en

Télécharger psury.eso

Retour à la liste

Numérotation des lignes : 

psury

  1. C PSURY     SOURCE    CB215821  20/07/29    21:16:07     10668          
  2.       SUBROUTINE PSURY(ENDO,NENDO,NVARI,NSTRS,MFR1,DEPST,XMAT,VAR0,
  3.      1                 RAPP,NRAPP,
  4.      1                 SIG0,SIGF,VARF,NMATT,DEFP,KERRE)
  5. *
  6. *   Entrées
  7. *
  8. *   ENDO: courbe de debut d'endommagement
  9. *   NENDO: nombre de points sur la courbe ENDO
  10. *   RAPP: courbe d'evolution de l'endommagement en fonction de la
  11. *         pseudo porosite
  12. *   NRAPP: nombre de points de la courbe RAPP
  13. *   NSTRS: nombre de composantes des deformations
  14. *   MFR1: numero de la formulation
  15. *   DEPST: increment de deformation totale
  16. *   XMAT: donnees materiau
  17. *   SIGF: contraintes plastiquement admissibles issues de ECOINC
  18. *   VAR0: variables internes au debut du pas de temps
  19. *   VARF: variables internes issues de ECOINC
  20. *   SIG0: contraintes au debut du pas de temps
  21. *   NMATT: nombre de composantes matériaux
  22. *   DEFP: incrément de la déformation plastique
  23. *
  24. *   Sorties
  25. *
  26. *   SIGF: contraintes finales ( tiennent compte de l'endommagement)
  27. *   VARF: variables internes finales
  28. *
  29. *
  30.       IMPLICIT INTEGER(I-N)
  31.       IMPLICIT REAL*8(A-H,O-Z)
  32. *
  33.  
  34. -INC PPARAM
  35. -INC CCOPTIO
  36. -INC CCREEL
  37.       DIMENSION SIGF(*),XMAT(*),VARF(*),DEPST(*),ENDO(*)
  38.       DIMENSION SIG0(*),VAR0(*),DEFP(*),RAPP(*)
  39.       DIMENSION RSIGF(6),RSIG0(6)
  40. *
  41. *
  42. *====================================================================
  43. *  Adaptation de l'option de calcul vers le 3D massif de SIGF a RSIGF
  44. *====================================================================
  45. *
  46.       IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
  47. *
  48. *---> 1 formulation massive
  49. *---> 2 formulation quasi incompressible
  50. *---> MASSIF 3D
  51. *
  52.         IF (NSTRS .EQ. 6) THEN
  53.           DO 110 I=1,NSTRS
  54.             RSIGF(I)=SIGF(I)
  55.             RSIG0(I)=SIG0(I)
  56.  110      CONTINUE
  57.         ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0)
  58.      &           .OR.(IFOUR.EQ.-1).OR.(IFOUR.EQ.-2))) THEN
  59. *
  60. *---> Calcul en mode deformations planes ou axisymetrique
  61. *---> Calcul en mode contraintes planes
  62. *
  63.           DO 115 I=1,NSTRS
  64.             RSIGF(I)=SIGF(I)
  65.             RSIG0(I)=SIG0(I)
  66.  115      CONTINUE
  67.           RSIGF(5)=0.D0
  68.           RSIG0(5)=0.D0
  69.           RSIGF(6)=0.D0
  70.           RSIG0(6)=0.D0
  71.         ENDIF
  72.       ELSE
  73.         KERRE = 99
  74.         RETURN
  75.       ENDIF
  76. *
  77. *========================================
  78. *   Calcul de la densité du materiau
  79. *========================================
  80. *
  81. *  Si on a oublie de rentrer la densite initiale
  82. *
  83. *      IF (IFOUR.EQ.-2) THEN
  84. *        DENS0=XMAT(NMATT-2)
  85. *      ELSE
  86. *        DENS0=XMAT(3)
  87. *      ENDIF
  88. *
  89. *      DENS0=XMAT(3)
  90.         DENS0=XMAT(NMATT-2)
  91.       IF (DENS0.LE.1.D-10) THEN
  92.          KERRE=71
  93.          RETURN
  94.       ENDIF
  95. *
  96. *  Initialisation de la densite en variable interne
  97. *
  98.       IF (VAR0(2).LE.1.D-10) THEN
  99.          VAR0(2)=DENS0
  100.       ENDIF
  101.       IF (IFOUR.EQ.-2) THEN
  102.          treps0=(1.D0-2.D0*XMAT(2))/(1.D0-XMAT(2))
  103.          treps0=treps0*(DEPST(1)+DEPST(2)-DEFP(1)-DEFP(2))
  104.       ELSE
  105.          treps0=DEPST(1)+DEPST(2)+DEPST(3)
  106.       ENDIF
  107.       rho=DENS0/((DENS0/VAR0(2))+treps0)
  108. *
  109. *
  110. *========================================
  111. *   A t'on déja endommagé ?
  112. *========================================
  113. *
  114.       IF (ABS(VAR0(3)).GT.1.D-10) THEN
  115. *
  116. *  On a déja endommagé
  117. *----------------------------------------
  118. *
  119. *   Calcul de la pseudo porosité du temps n+1
  120. *
  121.          alpha=(VAR0(3)-rho)/VAR0(3)
  122. *
  123. *   Calcul de la variable d'endommagement
  124. *
  125.          D_end1=0.D0
  126.          DO 12 I=1,NRAPP-1
  127.            DD1=RAPP(2*(I-1)+1)
  128.            DD2=RAPP(2*I+1)
  129.            AL1=RAPP(2*I)
  130.            AL2=RAPP(2*(I+1))
  131.            IF ((alpha.GE.AL1).AND.(alpha.LT.AL2)) THEN
  132.               D_end1=(DD2-DD1)/(AL2-AL1)*(alpha-AL1)+DD1
  133.            ENDIF
  134.  12      CONTINUE
  135.          IF ((D_end1.LE.0.D0).OR.(alpha.LE.0.D0)) D_end1=0.D0
  136.          IF ((D_end1.GE.1.D0).OR.(alpha.GE.1.D0)) D_end1=1.D0
  137.          D_max=VAR0(5)
  138. *
  139. *   Calcul de l'ancienne pseudo porosite ( temps n)
  140. *
  141.          alp_old=VAR0(4)
  142. *
  143. *   Calcul de la fonction d'endommagement g0
  144. *
  145.          IF ((alpha.GT.0.D0).AND.(alpha.GT.alp_old)) THEN
  146.             g0=1.D0-D_end1
  147.          ELSE IF ((alpha.GT.0.D0).AND.(alpha.LE.alp_old)) THEN
  148.             g0=1.D0-D_max
  149.          ELSE
  150.             g0=1.D0
  151.          ENDIF
  152.          g0=MAX(g0,0.D0)
  153. *
  154. *   Calcul des contraintes vérifiant l'endommagement
  155. *
  156.          DO 10 I=1,6
  157.             RSIGF(I)=RSIGF(I)*g0
  158.  10      CONTINUE
  159. *
  160. *   Mise a jour des variables internes
  161. *
  162. *---> Densite
  163.          VARF(2)=rho
  164. *---> Densite de debut de fracture
  165.          VARF(3)=VAR0(3)
  166. *---> Pseudo porosite
  167.          VARF(4)=alpha
  168. *---> Pseudo porosite maximale
  169.          VARF(5)=MAX(VAR0(5),D_end1)
  170. *---> Coefficient d'endommagement
  171.          VARF(6)=D_end1
  172. *---> Fonction d'endommagement
  173.          VARF(7)=g0
  174. *
  175.       ELSE
  176. *
  177. *   On a pas encore endommagé
  178. *-------------------------------------------
  179. *
  180. *   Calcul de P=1/3*trace(RSIGF)
  181. *
  182.          P0=(RSIGF(1)+RSIGF(2)+RSIGF(3))/3.D0
  183. *
  184. *   Calcul de la contrainte équivalente
  185. *
  186.          Y1=(RSIGF(1)*RSIGF(1))+(RSIGF(2)*RSIGF(2))
  187.          Y1=Y1+(RSIGF(3)*RSIGF(3))
  188.          Y1=Y1-(RSIGF(1)*RSIGF(2))-(RSIGF(2)*RSIGF(3))
  189.          Y1=Y1-(RSIGF(3)*RSIGF(1))
  190.          Y2=(RSIGF(4)*RSIGF(4))+(RSIGF(5)*RSIGF(5))
  191.          Y2=Y2+(RSIGF(6)*RSIGF(6))
  192.          Y0=Y1+(3.D0*Y2)
  193.          Y0=(Y0)**(0.5D0)
  194. *
  195. *   Rapport de triaxialite
  196. *
  197.          Rapp0=P0/Y0
  198. *
  199. *   Rapport de triaxialite de debut d'endommagement Rapp_th
  200. *
  201.          EPSE=VARF(1)
  202.          IF (EPSE.LT.ENDO(2)) THEN
  203.             Rapp_th=1.D30
  204.             pente1=1.D30
  205.             Rmax=Rapp_th
  206.             epsmin=ENDO(2)
  207.          ELSE IF (EPSE.GE.ENDO(2*NENDO)) THEN
  208.             Rapp_th=ENDO(2*NENDO-1)
  209.             pente1=0.D0
  210.             Rmax=Rapp_th
  211.             epsmin=ENDO(2*NENDO)
  212.          ELSE
  213.             DO 100 I=1,(NENDO-1)
  214.             IF ((EPSE.GE.ENDO(2*I)).AND.(EPSE.LT.ENDO(2*I+2))) THEN
  215.                IND0=I
  216.             ENDIF
  217.  100        CONTINUE
  218.             epsmin=ENDO(2*IND0)
  219.             epsmax=ENDO(2*IND0+2)
  220.             dif0=ABS(epsmax-epsmin)
  221.             Rmax=ENDO(2*IND0-1)
  222.             Rmin=ENDO(2*IND0+1)
  223.             IF (dif0.GT.1.D-10) THEN
  224.              pente1=(Rmax-Rmin)/(epsmin-epsmax)
  225.              Rapp_th=pente1*(EPSE-epsmin)+Rmax
  226.             ELSE
  227.              pente1=1.D30
  228.              Rapp_th=Rmin
  229.             ENDIF
  230.          ENDIF
  231. *
  232. *   Comparaison au rapport de triaxialite d'endommagement
  233. *
  234.          IF (Rapp0.GT.Rapp_th) THEN
  235. *
  236. *   On endommage
  237. *
  238. *   Calcul de la densite de debut d'endommagement : rho_f
  239. *---------------------------------------------------
  240. *
  241. *---> Rapport de triaxialite au debut du pas
  242. *
  243.              P_old=(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0
  244.              Y_old1=(RSIG0(1)*RSIG0(1))+(RSIG0(2)*RSIG0(2))
  245.              Y_old1=Y_old1+(RSIG0(3)*RSIG0(3))
  246.              Y_old1=Y_old1-(RSIG0(1)*RSIG0(2))-(RSIG0(2)*RSIG0(3))
  247.              Y_old1=Y_old1-(RSIG0(3)*RSIG0(1))
  248.              Y_old2=(RSIG0(4)*RSIG0(4))+(RSIG0(5)*RSIG0(5))
  249.              Y_old2=Y_old2+(RSIG0(6)*RSIG0(6))
  250.              Y_old0=Y_old1+(3.D0*Y2)
  251.              Y_old=(Y_old0)**(0.5D0)
  252.              R_old=P_old/Y_old
  253. *
  254. *---> Courbe d'evolution du rapport P/Y entre le debut et la fin du pas
  255. *     On l'approxime par une droite de pente pente0
  256. *
  257.              dvar1=VARF(1)-VAR0(1)
  258. *            IF (dvar1.GT.1.D-10) THEN
  259.              IF (dvar1.GT.abs(varf(1)*xzprec)) THEN
  260.                 pente0=(Rapp0-R_old)/dvar1
  261.              ELSE
  262.                 pente0=abs(varf(1))/XZPREC
  263.              ENDIF
  264. *
  265. *---> intersection avec la courbe de debut d'endommagement
  266. *     point (R_int,eps_int)
  267. *
  268.  150         IF (pente1.GE.1.D30) THEN
  269.               eps_int=epsmin
  270.              ELSE IF (pente0.GE.1.D30) THEN
  271.               eps_int=VAR0(1)
  272.              ELSE
  273.               eps_int=Rapp0-Rmax+(pente1*epsmin)-(pente0*VARF(1))
  274.               eps_int=eps_int/(pente1-pente0)
  275.              ENDIF
  276.              IF ((eps_int.LT.epsmin).AND.(IND0.GE.2)) THEN
  277.                 IND0=IND0-1
  278.                 epsmin=ENDO(2*IND0)
  279.                 epsmax=ENDO(2*IND0+2)
  280.                 dif0=ABS(epsmin-epsmax)
  281.                 Rmax=ENDO(2*IND0-1)
  282.                 Rmin=ENDO(2*IND0+1)
  283.                 IF (dif0.GT.1.D-10) THEN
  284.                  pente1=(Rmax-Rmin)/(epsmin-epsmax)
  285.                 ELSE
  286.                  pente1=1.D30
  287.                 ENDIF
  288.                 GOTO 150
  289.              ELSE IF ((eps_int.LT.epsmin).AND.(IND0.EQ.1)) THEN
  290.                 pente1=1.D30
  291.                 eps_int=ENDO(2)
  292.                 epsmin=ENDO(2)
  293.                 epsmax=ENDO(2)
  294.                 Rmin=ENDO(1)
  295.                 Rmax=1.D30
  296.              ENDIF
  297.              IF (pente1.GE.1.D30) THEN
  298.                 R_int=pente0*(eps_int-VARF(1))+Rapp0
  299.              ELSE
  300.                 R_int=pente1*(eps_int-epsmin)+Rmax
  301.              ENDIF
  302. *
  303. *---> Calcul de rho_f en supposant le module d'ecrouissage
  304. *     constant entre le debut et la fin du pas
  305. *
  306.              IF (dvar1.GT.1.D-10) THEN
  307.                 H0=(Y0-Y_old)/dvar1
  308.                 Y_int=H0*(eps_int-VAR0(1))+Y_old
  309.              ELSE
  310. *
  311. *---> Endommagement dans un cas élastique
  312. *
  313.                 IF ((ABS(Rapp0-R_old)).GT.1.D-20) THEN
  314.                    alfa0=(R_int-R_old)/(Rapp0-R_old)
  315.                    Y_int=(alfa0*Y0)+((1.D0-alfa0)*Y_old)
  316.                 ELSE
  317.                    Y_int=Y_old
  318.                 ENDIF
  319.              ENDIF
  320.              IF (treps0.GT.1.D-10) THEN
  321.                 XK0=XMAT(1)/(3.D0*(1.D0-2.D0*XMAT(2)))
  322.                 P_int=R_int*Y_int
  323.                 tr_eps1=(P_int-P_old)/XK0
  324.              ELSE
  325.                 tr_eps1=0.D0
  326.              ENDIF
  327.              rho_f=DENS0/((DENS0/VAR0(2))+tr_eps1)
  328.          IF ((rho_f.GT.1.D10).OR.(rho_f.LT.0.D0)) THEN
  329.             rho_f=rho
  330.          ENDIF
  331. *
  332. *    Calcul de la pseudo-porosite
  333. *
  334.          alpha=(rho_f-rho)/rho_f
  335. *
  336. *   Calcul de la variable d'endommagement
  337. *
  338.          D_end1=0.D0
  339.          DO 13 I=1,NRAPP-1
  340.            DD1=RAPP(2*(I-1)+1)
  341.            DD2=RAPP(2*I+1)
  342.            AL1=RAPP(2*I)
  343.            AL2=RAPP(2*(I+1))
  344.            IF ((alpha.GE.AL1).AND.(alpha.LT.AL2)) THEN
  345.               D_end1=(DD2-DD1)/(AL2-AL1)*(alpha-AL1)+DD1
  346.            ENDIF
  347.  13      CONTINUE
  348.          IF ((D_end1.LE.0.D0).OR.(alpha.LE.0.D0)) D_end1=0.D0
  349.          IF ((D_end1.GE.1.D0).OR.(alpha.GE.1.D0)) D_end1=1.D0
  350.          D_max=VAR0(5)
  351. *
  352. *   Calcul de l'ancienne pseudo porosite ( temps n)
  353. *
  354.          alp_old=VAR0(4)
  355. *
  356. *   Calcul de la fonction d'endommagement g0
  357. *
  358.          IF ((alpha.GT.0.D0).AND.(alpha.GT.alp_old)) THEN
  359.             g0=1.D0-D_end1
  360.          ELSE IF ((alpha.GT.0.D0).AND.(alpha.LE.alp_old)) THEN
  361.             g0=1.D0-D_max
  362.          ELSE
  363.             g0=1.D0
  364.          ENDIF
  365.          g0=MAX(g0,0.D0)
  366. *
  367. *   Calcul des contraintes vérifiant l'endommagement
  368. *
  369.          DO 15 I=1,6
  370.             RSIGF(I)=RSIGF(I)*g0
  371.  15      CONTINUE
  372. *
  373. *    Mise a jour des variables internes
  374. *
  375.              VARF(2)=rho
  376.              VARF(3)=rho_f
  377.              VARF(4)=alpha
  378.              VARF(5)=MAX(VAR0(5),D_end1)
  379.              VARF(6)=D_end1
  380.              VARF(7)=g0
  381.          IF (D_end1.GT.0.9) THEN
  382.             write(*,*) ' Debut endommagement'
  383.             write(*,*) 'P_int,P_old=',P_int,P_old
  384.             write(*,*) 'R_int,Y_int=',R_int,Y_int
  385.             write(*,*) 'r0,rho,tr_eps1=',DENS0,VAR0(2),tr_eps1
  386.             write(*,*) 'D_end1,gg0=',D_end1,g0
  387.             write(*,*) 'alpha,alp_old=',alpha,alp_old
  388.             write(*,*) 'DD1,DD2,AL1,AL2=',DD1,DD2,AL1,AL2
  389.             write(*,*) 'rho,rho_f=',rho,rho_f
  390.             write(*,*) 'DENS0,VAR0(2)=',DENS0,VAR0(2)
  391.          ENDIF
  392. *
  393.          ELSE
  394. *
  395. *   On n'endommage pas
  396. *
  397.              VARF(2)=rho
  398.              VARF(3)=0.D0
  399.              VARF(4)=0.D0
  400.              VARF(5)=0.D0
  401.              VARF(6)=0.D0
  402.              VARF(7)=1.D0
  403.          ENDIF
  404. *
  405. *=======================================
  406. *  Fin du calcul d'endommagement
  407. *=======================================
  408.       ENDIF
  409. *
  410. *
  411. *=========================================================
  412. *  Passage a l'option de calcul pour les contraintes
  413. *=========================================================
  414. *
  415.       IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
  416.         IF (NSTRS .EQ. 6) THEN
  417. *
  418. *---> MASSIF 3D
  419. *
  420.           DO 170 I=1,NSTRS
  421.             SIGF(I)=RSIGF(I)
  422.  170      CONTINUE
  423.         ELSE IF ( NSTRS .EQ. 4 ) THEN
  424. *
  425. *---> Calcul axisymétrique ou contraintes planes
  426. *
  427.           DO 180 I=1,NSTRS
  428.             SIGF(I)=RSIGF(I)
  429.  180      CONTINUE
  430.         ENDIF
  431.       ENDIF
  432.       RETURN
  433. *
  434.       END
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  

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