Télécharger la plaquette Cast3M


fr en

Télécharger ricj3.eso

Retour à la liste

Numérotation des lignes : 

ricj3

  1. C RICJ3     SOURCE    PV        21/11/04    21:15:16     11174          
  2. C sub ricjoi3d
  3. C RICJOI3D     SOURCE    CHAT      05/01/12    22:24:09     5004
  4.       SUBROUTINE RICJ3(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI,
  5.      &     DEPST,IFOUR,TETA1,TETA2,XMAT,
  6.      &     NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE)
  7. C-----------------------------------------------------------------------
  8. C Modele RICJOI en version element joint 3D
  9. C-----------------------------------------------------------------------
  10. C Développé par :       Benjamin Richard
  11. C
  12. C Supervisé par :       Frederic Ragueneau
  13. C                               Lucas Hector Adelaide
  14. C                               Christian Cremona
  15. C                               Jean Louis Tailhan
  16. C-----------------------------------------------------------------------
  17.  
  18. C
  19. C-----DECLARATION GENERALE----------------------------------------------
  20. C
  21.       IMPLICIT INTEGER(I-N)
  22.       IMPLICIT REAL*8(A-H,O-Z)
  23. C
  24. C-----APPEL AUX INCLUDES------------------------------------------------
  25. C
  26. -INC CCREEL
  27. C
  28. C-----PARAMETRES
  29. C
  30.       PARAMETER(XDILAT = 80.0D0)
  31. C
  32. C-----DIMENSION---------------------------------------------------------
  33. C
  34.       DIMENSION SIG0(*),EPIN0(*),VAR0(*),DEPST(*),XMAT(*)
  35.       DIMENSION DD(NSTRS,NSTRS)
  36.       DIMENSION SIGF(*),DEFP(*),VARF(*)
  37.       REAL*8 NCOEF,KCOEF
  38.  
  39.       FVARLIM = 2.0D-5
  40. C
  41. C-----COMPTEUR POUR LE NOMBRE DE PASSAGE DAS LA BOUCLE
  42. C
  43.       COUN  = VAR0(17)
  44. C
  45. C-----ETAT PRECEDENT----------------------------------------------------
  46. C
  47.       D     = VAR0(1 )
  48.       EPSF1 = VAR0(2 )
  49.       EPSF2 = VAR0(3 )
  50.       XEC1  = VAR0(4 )
  51.       XEC2  = VAR0(5 )
  52.       EPST1 = VAR0(6 )
  53.       EPST2 = VAR0(7 )
  54.       EPSN  = VAR0(8 )
  55.       CIN1  = VAR0(9 )
  56.       CIN2  = VAR0(10)
  57.       EPSR  = VAR0(11)
  58.       IFLAG = nint(VAR0(12))
  59.       FVAR  = VAR0(13)
  60.       SIGM  = VAR0(14)
  61.       XKS   = VAR0(15)
  62.       XKN   = VAR0(16)
  63.       crit2tol=1.001d0
  64. C
  65. C-----INITIALISATION SPECIFIQUE POUR LES VARIABLES INTERNES DIFF. DE ZERO
  66. C
  67.       IF (COUN.EQ.0.0D0) THEN
  68.          FVAR = 0.2D0
  69.          SIGM = XMAT(13)
  70.          XKS  = XMAT(1)
  71.          XKN  = XMAT(2)
  72.       ENDIF
  73. C
  74. C-----CONSTANTES MATERIAUX----------------------------------------------
  75. C
  76.       XKS0 = XMAT(1)
  77.       XKN0 = XMAT(2)
  78.  
  79.       XAD  = XMAT(5)
  80.       XY0  = XMAT(6)
  81.       XAL  = XMAT(7)
  82.       XGA  = XMAT(8)
  83.       XAA  = XMAT(9)
  84.  
  85.       Q1COE = XMAT(10)
  86.       Q2COE = XMAT(11)
  87.       Q3COE = XMAT(12)
  88.       SYCOE = XMAT(13)
  89.       NCOEF = XMAT(14)
  90.       KCOEF = XMAT(15)
  91.  
  92.       TCCOR = XMAT(16)
  93.       NGONF = nint(XMAT(17))
  94.  
  95.       CALL APP3(TCCOR,YR,15)
  96.  
  97.  
  98. C
  99. C-----ACTUALISATION DE LA DEFORMATION-----------------------------------
  100. C
  101.       EPST1 = EPST1 + DEPST(1)
  102.       EPST2 = EPST2 + DEPST(2)
  103.       EPSN  = EPSN  + DEPST(3)
  104.  
  105.       IF (NGONF.EQ.1) THEN
  106.          GOTO 123
  107.       ENDIF
  108.  
  109.       XSTE  = XKN*EPSN
  110.       XSAD  = XDILAT*0.5D0*(-ABS(XSTE)+XSTE)
  111.  
  112. C
  113. C-----CALCUL DE L ENERGIE-----------------------------------------------
  114. C
  115.       EPSNP = 0.5D0*(ABS(EPSN)+EPSN)
  116.       YN    = 0.5D0*XKN*EPSNP*EPSNP
  117.       YD1   = 0.5D0*XKS*EPST1*EPST1
  118.       YD2   = 0.5D0*XKS*EPST2*EPST2
  119.  
  120. C
  121. C-----CALCUL DU SEUIL D ENDOMMAGEMENT 1---------------------------------
  122. C
  123.       XFD1 = XAL*YN+YD1-(XY0+XEC1+YR)
  124.  
  125. C
  126. C-----TEST SUR LE SEUIL D ENDOMMAGEMENT---------------------------------
  127. C
  128.       IF (XFD1.GT.0.0D0) THEN
  129.          Dn  = 1.0D0-1.0D0/(1+XAD*(XAL*YN+YD1-XY0-YR))
  130.          XEC1 = XAL*YN+YD1-XY0-YR
  131.          IF (Dn.GT.D) THEN
  132.             D = Dn
  133.          ENDIF
  134.       ENDIF
  135.       IF (D.GE.1.0D0) THEN
  136.          D = 0.999D0
  137.       ENDIF
  138.  
  139. C
  140. C-----CALCUL DU SEUIL D ENDOMMAGEMENT 2---------------------------------
  141. C
  142.       XFD2 = XAL*YN+YD2-(XY0+XEC2+YR)
  143.  
  144. C
  145. C-----TEST SUR LE SEUIL D ENDOMMAGEMENT---------------------------------
  146. C
  147.       IF (XFD2.GT.0.0D0) THEN
  148.          Dn  = 1.0D0-1.0D0/(1+XAD*(XAL*YN+YD2-XY0-YR))
  149.          XEC2 = XAL*YN+YD2-XY0-YR
  150.          IF (Dn.GT.D) THEN
  151.            D = Dn
  152.          ENDIF
  153.       ENDIF
  154.       IF (D.GE.1.0D0) THEN
  155.          D = 0.999D0
  156.       ENDIF
  157.  
  158. C
  159. C-----FROTTEMENT 1-------------------------------------------------------
  160. C
  161.         SPI1  = D*XKS*(EPST1-EPSF1)
  162.  
  163.         CRITF1= ABS(SPI1-CIN1)+XSAD
  164.  
  165.            IF (CRITF1.GT.(0.0D0)) THEN
  166.               XF0   = CRITF1
  167.               XDLAM = 1.0D0
  168.  
  169.               DO I=1,500
  170.  
  171.                  SEUIL = (ABS(SPI1-CIN1)+XSAD)/XF0
  172.  
  173.                  IF ((SEUIL.LE.1.0D-5).OR.(XDLAM.LE.1.0D-10)) THEN
  174.                       GOTO 900
  175.                  ELSE
  176.                       XSG   = ABS(SPI1-CIN1)/(SPI1-CIN1)
  177.                       XDLAM = (ABS(SPI1-CIN1)+XSAD)/
  178.      &                    (XKS*D-XSG*XGA*(-XSG+XAA*CIN1))
  179.                     CIN1  = CIN1-XGA*XDLAM*(-XSG+XAA*CIN1)
  180.                       SPI1  = SPI1-D*XKS*XDLAM*XSG
  181.                  ENDIF
  182.  
  183.                  EPSF1 = -SPI1/(D*XKS)+EPST1
  184.  
  185.               ENDDO
  186.            ENDIF
  187.  
  188.  900     CONTINUE
  189.  
  190. C
  191. C-----FROTTEMENT 2-------------------------------------------------------
  192. C
  193.         SPI2  = D*XKS*(EPST2-EPSF2)
  194.  
  195.         CRITF2= ABS(SPI2-CIN2)+XSAD
  196.  
  197.            IF (CRITF2.GT.(0.0D0)) THEN
  198.               XF0   = CRITF2
  199.               XDLAM = 1.0D0
  200.  
  201.               DO I=1,500
  202.  
  203.                  SEUIL = (ABS(SPI2-CIN2)+XSAD)/XF0
  204.  
  205.                  IF ((SEUIL.LE.1.0D-5).OR.(XDLAM.LE.1.0D-10)) THEN
  206.                       GOTO 901
  207.                  ELSE
  208.                       XSG   = ABS(SPI2-CIN2)/(SPI2-CIN2)
  209.                       XDLAM = (ABS(SPI2-CIN2)+XSAD)/
  210.      &                    (XKS*D-XSG*XGA*(-XSG+XAA*CIN2))
  211.                     CIN2  = CIN2-XGA*XDLAM*(-XSG+XAA*CIN2)
  212.                       SPI2  = SPI2-D*XKS*XDLAM*XSG
  213.                  ENDIF
  214.  
  215.                  EPSF2 = -SPI2/(D*XKS)+EPST2
  216.  
  217.               ENDDO
  218.            ENDIF
  219.  
  220.  901     CONTINUE
  221.  
  222. C
  223. C-----ANELASTICITE--------------------------------------------------------
  224. C
  225.       IF (EPSN.GT.0.0D0) THEN
  226.          SIGN = (1-D)*XKN*EPSN
  227.       ELSE
  228.          SIGN =       XKN*(EPSN-EPSR)
  229.       ENDIF
  230.  
  231. C
  232. C-----CONTRAINTE DUE AUX EFFETS INELASTIQUE DE LA ROUILLE-----------------
  233. C
  234.       sigr=sign
  235.       IF (EPSN.LE.0.0D0) THEN
  236.          SIGR = -1.0D0*SIGN
  237.       ENDIF
  238.  
  239. C
  240. C-----CRITERE D ANELASTICITE DE LA ROUILLE--------------------------------
  241. C
  242.         IF (IFLAG.EQ.1) THEN
  243.                 GOTO 20
  244.         ELSE
  245.                 CRIT2 = 2.0D0*Q1COE*FVAR*DCOSH(Q2COE*SIGR/2.0D0/SYCOE)-
  246.      &                  (1.0D0+(Q3COE*FVAR)**2.0D0)
  247.         ENDIF
  248.  
  249. *      print*,'crit2=',crit2
  250.  
  251. C
  252. C SI LE CRITERE EST DEPASSE ON DEMARRE LE RETURN MAPPING
  253. C
  254.       IF ((CRIT2.GT.0.0D0).AND.
  255.      1    (EPSN .LE.0.0D0).AND.
  256.      2    (TCCOR.GT.0.0D0).AND.
  257.      3    (FVAR.GT.FVARLIM).AND.
  258.      4    (IFLAG.EQ.0).AND.
  259.      5    ((EPSN*DEPST(3)).GT.0.0D0)) THEN
  260.          CRIT20 = CRIT2
  261.            DLAM2  = 1.0D0
  262.          GOTO 10
  263.       ELSE IF (FVAR.EQ.FVARLIM) THEN
  264.            IFLAG = 1
  265.            GOTO 20
  266.         ELSE
  267.          GOTO 30
  268.       END IF
  269. C
  270. C-----DEBUT DES ITERATIONS INTERNES
  271. C
  272.  10   DO I=1,10000
  273. C
  274. C EVALUATION DU CRITERE
  275. C
  276.          CRIT2 = 2.0D0*Q1COE*FVAR*DCOSH(Q2COE*SIGR/2.0D0/SYCOE)-
  277.      &           (1.0D0+(Q3COE*FVAR)**2.0D0)
  278. C
  279. C TEST DE CONVERGENCE SUR LA VALEUR DU CRITERE RELATIVE
  280. C
  281.          IF ((ABS(CRIT2/CRIT20).LE.CRIT2TOL).OR.
  282.      &       (DLAM2.LT.0.0D0).OR.
  283.      &       (CRIT2.LT.0.0D0)) THEN
  284.             GOTO 20
  285.          END IF
  286. C
  287. C CALCUL DES DERIVEES DU CRITERES
  288. C
  289. C     dFr/dSIGR : DERIV21
  290. C     dFr/dSIGM : DERIV22
  291. C     dFr/dFVAR : DERIV23
  292. C
  293.          DERIV21 = SINH(Q2COE*SIGR/2.0D0/SIGM)*
  294.      &             Q1COE*Q2COE*FVAR/SIGM
  295.  
  296.          DERIV22 = SINH(Q2COE*SIGR/2.0D0/SIGM)*
  297.      &             (-1.0D0)*Q1COE*Q2COE*FVAR*SIGR/SIGM**2.0D0
  298.  
  299.          DERIV23 = 2.0D0*Q1COE*DCOSH(Q2COE*SIGR/2.0D0/SIGM)-
  300.      &             2.0D0*(Q3COE**2.0D0)*FVAR
  301. C
  302. C CALUL DU MODULE TANGENT COHERENTS ENTRE DSIGM ET DEPSM
  303. C
  304.          ETAN    = XKN/NCOEF*(SYCOE/SIGM)**(NCOEF-1.0D0)
  305.          ESTA    = 1.0D0/ETAN - 1.0D0/XKN
  306.          ESTA    = 1.0D0/ESTA
  307. C
  308. C CALCUL DU MULTIPLICATEUR DE LAGRANGE
  309. C
  310.          DLAM2   = -CRIT2/(XKN*(DERIV21)**2.0D0-
  311.      &             DERIV23*KCOEF*FVAR*(1.0D0-FVAR)*DERIV21-
  312.      &             ESTA/((1.0D0-FVAR)*SIGM)*DERIV21*SIGR*DERIV22)
  313. C
  314. C ACTUALISATION DES VARIABLES FORCES
  315. C
  316.          SIGR    = SIGR + DLAM2 * XKN * DERIV21
  317.          FVAR    = FVAR - DLAM2 * KCOEF * FVAR * (1.0D0-FVAR) * DERIV21
  318.          SIGM    = SIGM - DLAM2 * SIGR * DERIV21 *
  319.      &             ESTA/(1.0D0-FVAR)/SIGM
  320. C
  321. C SECURITE SUR LA VALEUR DE FVAR (CAS OU ELLE DEVIENT TROP FAIBLE)
  322. C
  323.          IF (FVAR.LE.FVARLIM) THEN
  324.               FVAR = FVARLIM
  325.                   MODR = nint(DLAM2 * XKN * DERIV21/DEPST(2) )
  326.          END IF
  327. C
  328. C --------------------------------------------
  329. C RETURN MAPPING : FIN DES ITERATIONS INTERNES
  330. C --------------------------------------------
  331. C
  332.       END DO
  333. C
  334. C BALISE DE SORTIE DU RETURN MAPPING AVANT ACUALISATION
  335. C
  336.  20   CONTINUE
  337. C
  338. C CAS OU LA FRACTION VOLUMIQUE DES PORES EST EGALE A 0.02
  339. C
  340.         IF ((FVAR.EQ.FVARLIM).AND.
  341.      &    (IFLAG.EQ.1)) THEN
  342.                 SIGR = SIGR - MODR * DEPST(3)
  343.         END IF
  344. C
  345. C CALCUL DES QUANTITES RESULTATS (EPSR, VKACT ET XKN)
  346. C
  347.       SIGN   = -SIGR
  348.       EPSR   = -(EPSN - SIGR/XKN)
  349. C
  350. C BALISE DE SORTIE DU RETURN MAPPING APRES ACUALISATION
  351. C
  352.  30   CONTINUE
  353. C
  354. C ACTUALISATION DES PARAMETRES ELASTIQUES
  355. C
  356.         XKNACT = 4.0D0*XKN0*XKS0*(1.0D0-FVAR)/
  357.      &      (4.0D0*XKS0+3.0D0*XKN0*FVAR)
  358.       XKSACT = XKS0*(1.0D0-FVAR)/
  359.      &        (1.0D0+((6.0D0*XKN0+12.0D0*XKS0)/
  360.      &        (9.0D0*XKN0+8.0D0*XKS0))*FVAR)
  361.  
  362.         XKN      = XKNACT
  363.         XKS    = XKSACT
  364.  
  365. C
  366. C-----CALCUL DES CONTRAINTES FINALES------------------------------------
  367. C
  368.       SIGF(1) = XKS*(1-D)*EPST1 + SPI1
  369.       SIGF(2) = XKS*(1-D)*EPST2 + SPI2
  370.  
  371.       IF (EPSN.GT.0.0D0) THEN
  372.          SIGF(3) = XKN*(1-D)*EPSN
  373.       ELSE
  374.          SIGF(3) = XKN*(EPSN-EPSR)
  375.       ENDIF
  376. C
  377. C-----STOCKAGE EN SORTIE DES VARIABLES INTERNES-------------------------
  378. C
  379.       VARF(1 ) = D
  380.       VARF(2 ) = EPSF1
  381.       VARF(3 ) = EPSF2
  382.       VARF(4 ) = XEC1
  383.       VARF(5 ) = XEC2
  384.       VARF(9 ) = CIN1
  385.       VARF(10) = CIN2
  386.       VARF(11) = EPSR
  387.       VARF(12) = IFLAG
  388.       VARF(13) = FVAR
  389.       VARF(14) = SIGM
  390.       VARF(15) = XKS
  391.       VARF(16) = XKN
  392.       VARF(17) = COUN + 1
  393.  
  394.  123  CONTINUE
  395.  
  396.       IF (NGONF.EQ.1) THEN
  397.          SIGF(1) = XMAT(1)*EPST1
  398.          SIGF(2) = XMAT(1)*EPST2
  399.          SIGF(3) = XMAT(2)*EPSN
  400.       ENDIF
  401.  
  402.       VARF(6 ) = EPST1
  403.       VARF(7 ) = EPST2
  404.       VARF(8 ) = EPSN
  405.  
  406.       RETURN
  407.       END
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  

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