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

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