Télécharger la plaquette Cast3M


fr en

Télécharger gurso2.eso

Retour à la liste

Numérotation des lignes : 

gurso2

  1. C GURSO2    SOURCE    OF166741  25/11/04    21:15:55     12349          
  2.  
  3.       SUBROUTINE GURSO2(DEPST,NSTRS,MFR1,IB,IGAU,
  4.      1     DSIGT,NCOMAT,SIG0,VAR0,XMAT,XCAR,NVARI,ICARA,
  5.      2     SIGF,VARF,DEFP,TRAC,KERRE,NCOURB)
  6.  
  7. *_________________________________________________________________
  8. *
  9. *
  10. *    ENTREES :
  11. *    ---------
  12. *
  13. * DEPST = INCREMENT DE DEFORMATIONS TOTALES
  14. * NSTRS = NBRE DE COMPOSANTES DES DEFORMATIONS
  15. * NCOMAT= NBRE DE CARACTERISTIQUES MECANIQUES DU MATERIAU
  16. * MFR1  = NUMERO DE LA FORMULATION
  17. * IB    = NUMERO DE L ELEMENT COURANT
  18. * IGAU  = NUMERO DU POINT COURANT
  19. * NVARI = NBRE DE VARIABLES INTERNES
  20. * SIG0(NSTRS) = CONTRAINTES AU DEBUT DU PAS D'INTEGRATION
  21. * VAR0(NVARI) = VARIABLES INTERNES AU DEBUT DU PAS DE TEMPS
  22. * XMAT(NCOMAT) = CARACTERISTIQUES MECANIQUES DU MATERIAU
  23. * ICARA = NBRE DE CARACTERISTIQUES GEOMETRIQUES DES ELEMENTS FINIS
  24. * XCARA(ICARA) = CARACTERISTIQUES GEOMETRIQUES DES ELEMENTS FINIS
  25. * TRAC = COURBE DE TRACTION NON ENDOMMAGEE
  26. * NCOURB = NOMBRE DE POINTS DE LA COURBE TRAC
  27. *
  28. *    SORTIE :
  29. *    --------
  30. *
  31. * SIGF(NSTRS) = CONTRAINTES FINALES
  32. * VARF(NVARI) = VARIALES INTERNES A LA FIN DU PAS D'INTEGRATION
  33. * DEFP(NSTRS) = INCREMENT DE DEFORMATION PLASTIQUE A LA FIN
  34. *               DU PAS D'INTEGRATION
  35. * ============================================================
  36. *       ICI IL FAUT PROGRAMMER EN FORTRAN PUR
  37. *=============================================================
  38. *
  39.       IMPLICIT INTEGER(I-N)
  40.       IMPLICIT REAL*8(A-H,O-Z)
  41. *
  42. -INC CCREEL
  43.  
  44. -INC PPARAM
  45. -INC CCOPTIO
  46.  
  47.       DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*),
  48.      &          VARF(*),DEFP(*),DSIGT(*),TRAC(*)
  49.  
  50.       DIMENSION RDEPS0(6)
  51.       DIMENSION RSIG0(6),RDEPS(6),RSIGF(6),RDEFP(6),RDIGT(6)
  52.       DIMENSION DEVT(6),SIGT(6),DEFT(6),SIGT1(6)
  53. *
  54. *  Adaptation de l'option de calcul vers le 3D massif de SIG0 a RSIG0
  55. *====================================================================
  56. *
  57.       IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
  58. *
  59. *---> 1 formulation massive
  60. *---> 2 formulation quasi incompressible
  61. *---> MASSIF 3D
  62. *
  63.          IF (NSTRS .EQ. 6) THEN
  64.             DO 110 I=1,NSTRS
  65.                RSIG0(I)=SIG0(I)
  66.                RDEPS0(I)=DEPST(I)
  67.  110        CONTINUE
  68.          ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0)
  69.      &           .OR.(IFOUR .EQ. -1).OR.(IFOUR.EQ.-2))) THEN
  70. *
  71. *---> Calcul en mode deformations planes ou axisymetrique
  72. *     ou contraintes planes
  73. *
  74.             DO 115 I=1,NSTRS
  75.                RSIG0(I)=SIG0(I)
  76.                RDEPS0(I)=DEPST(I)
  77.  115        CONTINUE
  78.             RSIG0(5)=0.D0
  79.             RSIG0(6)=0.D0
  80.             RDEPS0(5)=0.D0
  81.             RDEPS0(6)=0.D0
  82.          ENDIF
  83.       ELSE
  84.          KERRE = 99
  85.          RETURN
  86.       ENDIF
  87. *
  88. *
  89. *   Passage des deformations de cisaillement exprimées
  90. *   en GAMA aux déformations de cisaillement exprimées
  91. *   en déformations
  92. *
  93.       DO 116 I=1,6
  94.          IF (I.GT.3) RDEPS0(I)=RDEPS0(I)*0.5D0
  95.          RDEPS(I)=0.D0
  96.  116  CONTINUE
  97. *
  98. *   Données du materiau
  99. *===========================================================
  100. *
  101. C Ordre des composants dans XMAT (different de celui de VALMAT !) :
  102. C XMAT(1)  : YOUNG
  103. C XMAT(2)  : NU
  104. C XMAT(3)  : RHO
  105. C XMAT(4)  : ALPH
  106. C XMAT(5)  : ECRO
  107. C XMAT(6)  : Q
  108. C XMAT(7)  : FU
  109. C XMAT(8)  : FF
  110. C XMAT(9)  : FC
  111. C XMAT(10) : FNS0
  112. C XMAT(11) : FNE0
  113. C XMAT(12) : SNS
  114. C XMAT(13) : SNE
  115. C XMAT(14) : SIGN
  116. C XMAT(15) : EPSN
  117. C XMAT(16) : F0
  118. C XMAT(17) : VISQ
  119. C XMAT(18) : SRMA
  120. C XMAT(19) : Q2
  121. C XMAT(20) : Q3
  122. C XMAT(21) : TREF
  123. C XMAT(22) : TALP
  124.  
  125.       YOUNG=XMAT(1)
  126.       XNU=XMAT(2)
  127.       Q1=XMAT(6)
  128.       FU0=XMAT(7)
  129.       FF0=XMAT(8)
  130.       FC0=XMAT(9)
  131.       FNS0=XMAT(10)
  132.       FNE0=XMAT(11)
  133.       SNS0=XMAT(12)
  134.       SNE0=XMAT(13)
  135.       SIGN0=XMAT(14)
  136.       EPSN0=XMAT(15)
  137.       F0=XMAT(16)
  138.       XSRMA=XMAT(18)
  139.       Q2 = XMAT(19)
  140.       Q3 = XMAT(20)
  141.       Fmin0=1.D-8
  142.       IF ((FF0.LE.Fmin0).OR.(FF0.GT.1.D0)) FF0=1.D0
  143.       IF ((F0.LE.Fmin0).OR.(F0.GE.1.D0)) F0=0.D0
  144. C      RMIN0=(1.D0-FF0)/(1.D0-F0)
  145. *
  146. *  Calculs préliminaires
  147. *===================================
  148. *
  149. *---> Module de cisaillement G0 et de pression hydro XK0
  150. *
  151.       G0=YOUNG/(2.D0*(1.D0+XNU))
  152.       XK0=YOUNG/(3.D0*(1.D0-(2.D0*XNU)))
  153. *
  154. *  Début des sous incrémentations
  155. *===================================
  156. *
  157.       recal0=-100.D0
  158.       iter2=0
  159.       nmax0=1
  160. *
  161. *---> Variables internes test
  162. *
  163. *   - déformation plastique équivalente
  164. *   - fraction volumique de cavités
  165. *
  166.  93   EPSPT=VAR0(1)
  167.       FT=VAR0(2)
  168.       SIGMT=VAR0(3)
  169.       EPSMT=VAR0(4)
  170.       RHOT=VAR0(5)
  171.       FNST=VAR0(6)
  172.       FNET=VAR0(7)
  173.       VARF(3)=VAR0(3)
  174.       VARF(4)=VAR0(4)
  175.       VARF(6)=VAR0(6)
  176.       VARF(7)=VAR0(7)
  177.       IF (RHOT.LE.1.D-30) THEN
  178.          RHOT=1.D0
  179.          VAR0(5)=1.D0
  180.       ENDIF
  181.       IF (FT.LE.Fmin0) FT=0.D0
  182. *
  183. *  Contraintes planes
  184. *    - on stocke les variables internes au début de
  185. *    chaque sous incrémentation
  186. *    - on stocke les déformations plastiques au début de
  187. *    chaque sous incrémentation
  188. *
  189.       IF (IFOUR.EQ.-2) THEN
  190.          EPSPT1=EPSPT
  191.          FT1=FT
  192.          EPSMT1=EPSMT
  193.          SIGMT1=SIGMT
  194.          RHOT1=RHOT
  195.          FNST1=FNST
  196.          FNET1=FNET
  197.          def1=0.D0
  198.          def2=0.D0
  199.          def3=0.D0
  200.          rdeps3=0.D0
  201.       ENDIF
  202. *
  203. *  Contraintes planes
  204. *    - on stocke les contraintes au début de
  205. *    chaque sous incrémentation
  206. *
  207.       IF (IFOUR.EQ.-2) THEN
  208.          DO 171 I=1,6
  209.             SIGT1(I)=RSIG0(I)
  210.  171     CONTINUE
  211.       ENDIF
  212.       DO 172 I=1,6
  213.          SIGT(I)=RSIG0(I)
  214.  172  CONTINUE
  215.       SMT0=(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0
  216.       SMT=SMT0
  217.       IF (IFOUR.EQ.-2) SMT1=SMT0
  218. *
  219. *  On divise le chargement pour les sous incrementations
  220. *
  221.  95   DO 118 I=1,6
  222.          RDEPS(I)=RDEPS0(I)/nmax0
  223.  118  CONTINUE
  224.       iter2=iter2+1
  225. *
  226. *  Initialisations
  227. *===================================
  228. *
  229.       iter1=0
  230. *
  231. *---> paramètre de recalcul éventuel de la nucléation
  232. *     il est ré-initialisé à chaque sous itération
  233. *
  234.       recal1=-100.D0
  235. *
  236. *---> Calcul de la trace de l'incrément de déformation
  237. *
  238.       IF (IFOUR.EQ.-2) THEN
  239.          XLAM0=(1.D0-FT)*XK0-(2.D0*G0/3.D0)
  240.          RDEPS(3)=-1.D0*XLAM0/(2.D0*G0+XLAM0)*(RDEPS(1)+RDEPS(2))
  241.       ENDIF
  242.  98   treps0=RDEPS(1)+RDEPS(2)+RDEPS(3)
  243. *
  244. *
  245. *  Calculs des grandeurs élastiques test
  246. *====================================================
  247. *
  248. *---> Variables internes test: cas contraintes planes
  249. *
  250. *   - déformation plastique équivalente
  251. *   - fraction volumique de cavités
  252. *
  253.       IF (IFOUR.EQ.-2) THEN
  254.          EPSPT=EPSPT1
  255.          FT=FT1
  256.          EPSMT=EPSMT1
  257.          SIGMT=SIGMT1
  258.          RHOT=RHOT1
  259.          FNST=FNST1
  260.          FNET=FNET1
  261.          SMT=SMT1
  262.       ENDIF
  263. *
  264. *---> Calcul des contraintes test élastiques
  265. *     et de l'incrément des déformations plastiques test au
  266. *     cours du pas
  267. *
  268.       DO 101 I=1,6
  269.          A=0.D0
  270.          IF (I.LE.3) A=1.D0
  271.          IF (IFOUR.EQ.-2) SIGT(I)=SIGT1(I)
  272.          DEVT(I)=SIGT(I)-(A*SMT)
  273.          DEVT(I)=DEVT(I)+(2.D0*G0*(RDEPS(I)-(A*treps0/3.D0)))
  274.          DEFT(I)=0.D0
  275.  101  CONTINUE
  276. *
  277. *---> Contrainte équivalente STEST
  278. *
  279.       STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0
  280.       IF (STEST2.LE.1.D-10) STEST2=0.D0
  281.       Stest=STEST2**0.5D0
  282. *
  283. *---> Limite d'élasticité lue sur la courbe d'écrouissage
  284. *
  285. *   - Recherche des valeurs de déformations plastiques équivalentes
  286. *     encadrant la valeur courante
  287. *
  288.       EPSM0=TRAC(2*NCOURB)
  289.       IF (EPSMT.GE.EPSM0) THEN
  290.          Y1=TRAC(2*NCOURB-1)
  291.          Y2=TRAC(2*NCOURB-1)
  292.          EPS1=EPSM0
  293.          EPS2=2.D0*EPSM0
  294.       ELSE
  295.          DO 300 I=1,(NCOURB-1)
  296.             EPS1=TRAC(2*I)
  297.             EPS2=TRAC(2*(I+1))
  298.             IF ((EPSMT.GE.EPS1).AND.(EPSMT.LT.EPS2)) THEN
  299.                Y1=TRAC(2*I-1)
  300.                Y2=TRAC(2*(I+1)-1)
  301.                GOTO 96
  302.             ENDIF
  303.  300     CONTINUE
  304.       ENDIF
  305. *
  306. *   - Limite d'élasticité
  307. *
  308.  96   H1=(Y2-Y1)/(EPS2-EPS1)
  309.       Ytest=(H1*(EPSMT-EPS2))+Y2
  310. *
  311. *---> Fraction de densité test
  312. *
  313. * Correction SEMI: Mars 2017: le rapport de densite est exprime en deformation logarithmique 
  314. *=========================================
  315. * Dans la note 96-566, le rapport de densite est deduit de l'expression lineaire de la variation de volume:  
  316. *       trace(epsilon)=dV/V0 
  317. *    avec: V le volume total, V0 le volume total initial et epsilon la deformation totale 
  318. *    On en deduit: V/V0=1+trace(epsilon) 
  319. *    Puis le rapport de densite RHO=V0/V=1/(1+trace(epsilon)) 
  320. *    Cette expression n'est pas imposée dans l'article de Needleman et Tvergaard de réference   
  321. *      "AN ANALYSIS OF DUCTILE RUPTURE MODES AT A CRACK TIP" 
  322. *       A. NEEDLEMAN, V. TVERGAARD, 
  323. *       J. Mech. Phys. Solids Vol. 35, No. 2, pp. 151-183, 1987 
  324. * Si on exprime la variation de volume par son expression logarithmique: 
  325. *       trace(depsilon)=dln(V)
  326. *    avec depsilon la variation de la deformation totale et dln(V) la variation du logarithme du volume V
  327. *    On en deduit: ln(V/V0)=trace(epsilon) 
  328. *    Puis le rapport de densite RHO=exp(-trace(epsilon)) 
  329. *    Cette expression donne de meilleurs résultats pour les très forts taux de déformation. 
  330. *========================================== 
  331. *      RHOT=1.D0/((1.D0/RHOT)+treps0)
  332.       treps00 = LOG((1.D0/RHOT))
  333.       treps00 = treps00 + treps0 
  334.       RHOT = EXP(-1. * treps00)
  335. * Fin Correction SEMI: Mars 2017
  336.       IF (RHOT.LT.1.D-10) RHOT=1.D-10
  337. **      IF ((FT.GT.Fmin0).AND.(RHOT.LT.RMIN0)) RHOT=RMIN0
  338. *
  339. *---> Coefficients de nucléation des cavités
  340. *
  341.       dif0=VARF(3)-Ytest-SMT
  342.       IF ((dif0.LE.1.D-5).AND.(treps0.GE.0.D0)) THEN
  343.          BB0=EXP(-.5D0*((VARF(3)-SIGN0)/SNS0)**2.D0)
  344.          BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0))
  345.       ELSE
  346.          BB0=0.D0
  347.       ENDIF
  348. *
  349. *---> Cas de la nucléation pilotée en contrainte
  350. *
  351.       IF ((FNS0.GT.Fmin0).AND.(recal1.LT.0.D0)) THEN
  352. *
  353. *  Si on se trouve dans la "fenetre" où la
  354. *  nucléation apparait, on sous incrémente le pas
  355. *  de calcul afin de suivre les variations de BB0
  356. *
  357.          dif0=Ytest+SMT
  358.          IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND.
  359.      .        (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN
  360.             xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0)
  361.             nmax1=NINT(xmax1)*nmax0
  362.             IF (nmax1.GT.2000) nmax1=2000
  363.             IF (nmax1.GT.nmax0) THEN
  364.                nmax0=nmax1
  365.                iter2=0
  366.                recal1=100.D0
  367.                GOTO 93
  368.             ENDIF
  369.          ENDIF
  370. *
  371.       ENDIF
  372. *
  373. *---> Fraction de cavité test
  374. *     On résoud une équation du second degré
  375. *
  376.       AA1=BB0*XK0/((1.D0-F0)*RHOT)
  377.       BB1=BB0*XK0-1.D0
  378.       CC1=BB0*SMT+1.D0-FT
  379.       delta0=BB1*BB1+4.D0*AA1*CC1
  380. *
  381.       IF (AA1.GT.1.D-5) THEN
  382. *
  383. *   Le terme en contrainte de la nucléation n'est pas négligeable
  384. *
  385.          XX1=(BB1+(delta0**.5D0))/(2.D0*AA1)
  386.          FT=1.D0-XX1
  387.       ELSE
  388. *
  389. *   Le terme en contrainte de la nucléation est négligé
  390. *
  391.          BB0=0.D0
  392.       ENDIF
  393.       FT=MAX(FT,0.D0)
  394. *
  395. *---> Contrainte totale et contrainte moyenne test
  396. *
  397.       FNST=FNST-BB0*SMT
  398.       SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT-
  399.      . 1.D0)
  400.       FNST=FNST+BB0*SMT
  401.       DO 200 I=1,6
  402.          A=0.D0
  403.          IF (I.LE.3) A=1.D0
  404.          SIGT(I)=DEVT(I)+(A*SMT)
  405.  200  CONTINUE
  406. *
  407. *---> Fonction d'endommagement
  408. *
  409.       IF (FT.LE.FC0) THEN
  410.          FST0=FT
  411.          ddel0=1.D0
  412.       ELSE
  413.          FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0)
  414.          ddel0=(FU0-FC0)/(FF0-FC0)
  415.       ENDIF
  416.       IF (FST0.GT.FU0) FST0=FU0
  417. *
  418. *---> Terme d'endommagement
  419. *
  420.       SMTM=MAX(SMT,0.D0)
  421.       SMTM=SMT
  422.       Gtest=1.D0+(Q3*FST0*FST0)-
  423.      .     (2.D0*Q1*FST0*COSH(3.D0*Q2*SMTM/(2.D0*Ytest)))
  424. *
  425. *---> Fonction de charge test PHIT
  426. *
  427.       PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest
  428. *
  429. *---> Cas du matériau complètement endommagé
  430. *
  431.  91   IF (FT.GE.FF0) THEN
  432. *
  433. *  Les cavités
  434. *
  435.          FT=FF0+Fmin0
  436. *
  437. *  Les contraintes sont mises a zero
  438. *  Les deformations sont plastiques
  439. *
  440.          SMT=0.D0
  441.          Stest=0.D0
  442.          treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/(3.D0*XK0)
  443.          DO 12 I=1,6
  444.             SIGT(I)=0.D0
  445.             A=0.D0
  446.             IF (I.LE.3) A=1.D0
  447.             depel0=SIGT(I)-RSIG0(I)
  448.             DEVT(I)=(depel0-(SMT-SMT0)*A)/(2.D0*G0)
  449.             depel0=DEVT(I)+(treps1*A)
  450.             DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  451.             IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  452.                DEFT(3)=rdeps3+RDEPS(3)-depel0
  453.             ENDIF
  454.  12      CONTINUE
  455. *
  456. *  Deformation plastique equivalente
  457. *
  458.          treps1=(DEFT(1)+DEFT(2)+DEFT(3))/3.D0
  459.          DO 14 I=1,6
  460.             A=0.D0
  461.             IF (I.LE.3) A=1.D0
  462.             DEVT(I)=DEFT(I)-treps1*A
  463.  14      CONTINUE
  464.          IF (treps1.LE.0.D0) THEN
  465.             DO 16 I=1,3
  466.                DEFT(I)=DEVT(I)
  467.  16         CONTINUE
  468.          ENDIF
  469.          dlam0=SQRT((2.D0/3.D0)*PROCON(DEVT,DEVT,6))
  470.          EPSPT=EPSPT+dlam0
  471. **        RHOT=RMIN0
  472. *
  473.          Gtest=0.D0
  474.          PHIT=0.D0
  475. *
  476.       ENDIF
  477. *
  478. *  Vérification du critère de plasticité
  479. *=========================================================
  480. *
  481. *---> Erreur admise
  482. *
  483.       ERR0=1.D-7*(ABS(PHIT))
  484.       ERR0=MAX(ERR0,1.D-8*Ytest*Ytest)
  485. *
  486. *---> Critère de plasticité
  487. *
  488.       PHI0=PHIT
  489.       iter0=0
  490. *
  491.  99   IF (PHI0.LE.ERR0) THEN
  492. *
  493. *  On est élastique
  494. *=========================================================
  495. *
  496. *--------------------------------------------------------------
  497. *     Cas particulier des contraintes planes
  498. *--------------------------------------------------------------
  499. *
  500.          IF (IFOUR.EQ.-2) THEN
  501. *
  502. *---> On vérifie qu'on est en contraintes planes
  503. *
  504.             iter1=iter1+1
  505.             SIG3=ABS(SIGT(3))
  506.             SIG30=MAX(Stest*1.D-5,1.D-3)
  507.             IF (SIG3.GT.SIG30) THEN
  508.                XLAM0=(1.D0-FT)*XK0-(2.D0*G0/3.D0)
  509.                RDEPS(3)=-1.D0*(XLAM0/(2.D0*G0+XLAM0))*
  510.      .              (RDEPS(1)-DEFT(1)+RDEPS(2)-DEFT(2)+def1+def2)
  511.                RDEPS(3)=RDEPS(3)+DEFT(3)-def3
  512.                IF (iter1.LE.200)THEN
  513.                   GOTO 98
  514.                ELSE
  515.                   KERRE=460
  516.                ENDIF
  517.             ELSE
  518.                EPSPT1=EPSPT
  519.                FT1=FT
  520.                SIGMT1=SIGMT
  521.                EPSMT1=EPSMT
  522.                RHOT1=RHOT
  523.                FNST1=FNST
  524.                FNET1=FNET
  525.                SMT1=SMT
  526.                SIGT1(1)=SIGT(1)
  527.                SIGT1(2)=SIGT(2)
  528.                SIGT1(3)=SIGT(3)
  529.                SIGT1(4)=SIGT(4)
  530.                SIGT1(5)=SIGT(5)
  531.                SIGT1(6)=SIGT(6)
  532.                rdeps3=rdeps3+RDEPS(3)
  533.                def1=DEFT(1)
  534.                def2=DEFT(2)
  535.                def3=DEFT(3)
  536.             ENDIF
  537.          ENDIF
  538. *
  539. *-----------------------------------------------------
  540. *    Fin du cas particulier des contraintes planes
  541. *-----------------------------------------------------
  542. *
  543.          DO 400 I=1,6
  544.             RSIGF(I)=SIGT(I)
  545.             RDEFP(I)=DEFT(I)
  546.             RDIGT(I)=RSIGF(I)-RSIG0(I)
  547.  400     CONTINUE
  548.          VARF(1)=EPSPT
  549.          VARF(2)=FT
  550.          VARF(3)=MAX(VARF(3),Ytest+SMT)
  551.          VARF(4)=MAX(VARF(4),EPSMT)
  552.          VARF(5)=RHOT
  553.          VARF(6)=FNST
  554.          VARF(7)=FNET
  555. *
  556.       ELSE
  557. *
  558. *  On plastifie
  559. *=========================================================
  560. *
  561. *---> Valeur minimale de Stest au dessous de laquelle
  562. *     les contraintes sont considérées comme étant
  563. *     purement hydrostatiques
  564. *
  565.          Gmin0=2.D0*Q1*FST0*COSH(3.D0*Q2*SMTM/(2.D0*Ytest))
  566.          Gmin0=Gmin0-((Q3*FST0*FST0))
  567.          Gmin0=(ABS(Gmin0))**.5D0
  568.          Gmin0=Gmin0*Ytest*1.D-6
  569. *
  570. *---> Coefficient de nucléation
  571. *
  572.          dif0=VARF(4)-EPSMT
  573.          IF (dif0.LE.1.D-5) THEN
  574.             DD0=EXP(-.5D0*((VARF(4)-EPSN0)/SNE0)**2.D0)
  575.             DD0=DD0*FNE0/(SNE0*((2.D0*XPI)**.5D0))
  576.          ELSE
  577.             DD0=0.D0
  578.          ENDIF
  579. *
  580. *---> Condition de consistance
  581. *
  582.          tr0=3.D0*Q1*Q2*FST0*Ytest*SINH(3.D0*Q2*SMTM/(2.D0*Ytest))
  583.          IF (Stest.LT.Gmin0) THEN
  584.             DE0=SMT*tr0/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest)
  585.          ELSE
  586.             tr0=tr0/(2.D0*Stest)
  587.             DE0=(SMT*tr0+Stest)/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest)
  588.          ENDIF
  589. *
  590.          H0=H1*DE0
  591.          DS1=XK0*(2.D0*(1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT-1.D0)
  592.          DF0=((1.D0-FT+XSRMA*(FNET+FNST))*tr0+DD0*DE0+BB0*H0)/
  593.      .    (1.D0+BB0*DS1)
  594. *
  595.          DCOS0=(DS1*DF0/(2.D0*Ytest))+(SMT*H0)/(2.D0*Ytest*Ytest)
  596.          DCOS0=DCOS0*3.D0*Q2*SINH(3.D0*Q2*SMTM/(2.D0*Ytest))
  597. *
  598.          DG0=2.D0*Q3*ddel0*DF0*FST0
  599.          DG0=DG0+2.D0*Q1*FST0*DCOS0
  600.          DG0=DG0-2.D0*Q1*ddel0*DF0*COSH(3.D0*Q2*SMTM/(2.D0*Ytest))
  601. *
  602.          IF (Stest.LT.Gmin0) THEN
  603.             denom0=Ytest*Ytest*DG0
  604.          ELSE
  605.             denom0=6.D0*G0*Stest+Ytest*Ytest*DG0
  606.             denom0=denom0+2.D0*H0*Gtest*Ytest
  607.          ENDIF
  608. *
  609.          IF (denom0.GT.1.D-5) THEN
  610.             dlam0=PHIT/denom0
  611.          ELSE
  612.             IF (recal0.LE.0.D0) THEN
  613.                recal0=100.D0
  614.                nmax0=2000
  615.                iter2=0
  616.                GOTO 93
  617.             ENDIF
  618.          ENDIF
  619. *
  620. *
  621. *
  622. *     Vérification des hypothèses de calcul
  623. *===============================================================
  624. *
  625.          IF (recal0.LE.0.D0) THEN
  626. *
  627.             xnum0=10.D0
  628. *
  629. * Correction SEMI: Mars 2017. 
  630. * critere de verification de la precision de la linearisation de la loi GTN
  631. *            xlim0=1.D-1
  632.             xlim0=1.D-2
  633. * Fin Correction SEMI: Mars 2017. 
  634.             xmax1=0.D0
  635.             xmax2=0.D0
  636.             xmax3=0.D0
  637.             xmax4=0.D0
  638.             xmax5=0.D0
  639.             xmax6=0.D0
  640.             xmax7=0.D0
  641. *
  642. *---> Condition 3.D0*G*dlam0/Stest << 1
  643. *
  644.             IF (Stest.GE.Gmin0) THEN
  645.                cri01=ABS(3.D0*G0*dlam0/Stest)
  646.                IF (cri01.GT.xlim0) THEN
  647.                   xmax1=xnum0*cri01+5.D0
  648.                ENDIF
  649.             ENDIF
  650. *
  651. *---> Condition H0*dlam0/Ytest << 1
  652. *
  653.             cri01=ABS(H0*dlam0/Ytest)
  654.             IF (cri01.GT.xlim0) THEN
  655.                xmax2=xnum0*cri01+5.D0
  656.             ENDIF
  657. *
  658. *---> Condition DCOS0*dlam0/COSH << 1
  659. *
  660.             cri01=ABS(DCOS0*dlam0/COSH(3.D0*Q2*SMTM/(2.D0*Ytest)))
  661.             IF (cri01.GT.xlim0) THEN
  662.                xmax3=xnum0*cri01+5.D0
  663.             ENDIF
  664. *
  665. *---> Condition DF0*dlam0/F << 1
  666. *
  667.             Fmax0=MAX(FNS0,FNE0)
  668.             Fmax0=MAX(Fmax0,F0)
  669.             Fmax0=MAX(Fmax0,FT)
  670.             Fmax0=MAX(Fmax0,Fmin0)
  671.             cri01=ABS(DF0*dlam0/Fmax0)
  672.             IF (cri01.GT.xlim0) THEN
  673.                xmax4=xnum0*cri01+5.D0
  674.             ENDIF
  675. *
  676. *---> Condition DG0*dlam0/G << 1
  677. *
  678.             Gmax0=MAX(ABS(Gtest),1.D0)
  679.             cri01=ABS(DG0*dlam0/Gmax0)
  680.             IF (cri01.GT.xlim0) THEN
  681.                xmax5=xnum0*cri01+5.D0
  682.             ENDIF
  683. *
  684. *---> Condition DS1*dlam0/SMT << 1
  685. *
  686.             SMTmax=MAX(ABS(SMT),(1.D-3*TRAC(1)))
  687.             cri01=ABS(DS1*DF0*dlam0/SMTmax)
  688.             IF (cri01.GT.xlim0) THEN
  689.                xmax6=xnum0*cri01+5.D0
  690.             ENDIF
  691. *
  692. *---> Condition tr0*dlam0/tr(EPSP) << 1
  693. *
  694.             tremp0=(FT-F0)/(1.D0-F0)/RHOT
  695.             IF (ABS(tremp0).GT.1.D-5) THEN
  696.                cri01=ABS(tr0*dlam0/tremp0)
  697.                IF (cri01.GT.xlim0) THEN
  698.                   xmax7=xnum0*cri01+5.D0
  699.                ENDIF
  700.             ENDIF
  701. *
  702. *
  703.             IF (xmax1.GT.xmax2) THEN
  704.                xmax00=xmax1
  705.             ELSE
  706.                xmax00=xmax2
  707.             ENDIF
  708.             IF (xmax3.GT.xmax00) xmax00=xmax3
  709.             IF (xmax4.GT.xmax00) xmax00=xmax4
  710.             IF (xmax5.GT.xmax00) xmax00=xmax5
  711.             IF (xmax6.GT.xmax00) xmax00=xmax6
  712.             IF (xmax7.GT.xmax00) xmax00=xmax7
  713. *
  714.             nmax00=NINT(xmax00)
  715.             recal0=100.D0
  716.             IF (nmax00.GT.nmax0) nmax0=nmax00
  717.             IF (nmax0.GT.2000) nmax0=2000
  718. *
  719. *---> On recommence le calcul en sous incrémentant le
  720. *     pas de temps
  721. *
  722.             IF (nmax0.GE.2) THEN
  723.                iter2=0
  724.                GOTO 93
  725.             ENDIF
  726. *
  727.          ENDIF
  728. *
  729. *  Cas particulier de la nucléation
  730. *============================================================
  731. *
  732.          IF (recal1.LT.0.D0) THEN
  733. *
  734. *---> Condition sur la nucléation pilotée en contrainte
  735. *
  736.             IF (FNS0.GT.Fmin0) THEN
  737. *
  738. *  Si on se trouve dans la fenetre, on sous incremente
  739. *
  740.                dif0=Ytest+H0*dlam0+SMT-DS1*DF0*dlam0
  741.                IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND.
  742.      .              (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN
  743.                   xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0)
  744.                   nmax1=NINT(xmax1)*nmax0
  745.                   IF (nmax1.GT.2000) nmax1=2000
  746.                   IF (nmax1.GT.nmax0) THEN
  747.                      nmax0=nmax1
  748.                      iter2=0
  749.                      recal1=100.D0
  750.                      GOTO 93
  751.                   ENDIF
  752.                ENDIF
  753. *
  754.             ENDIF
  755. *
  756. *---> Condition sur la nucléation pilotée en deformation
  757. *
  758.             IF (FNE0.GT.Fmin0) THEN
  759. *
  760. *  Si on se trouve dans la fenetre, on sous incremente
  761. *
  762.                dif0=EPSMT+DE0*dlam0
  763.                IF ((dif0.GE.(EPSN0-SNE0*2.D0)).AND.
  764.      .              (VARF(4).LE.(EPSN0+SNE0*2.D0))) THEN
  765.                   xmax1=200.D0*(dif0-VARF(4))/(4.D0*SNE0)
  766.                   nmax1=NINT(xmax1)*nmax0
  767.                   IF (nmax1.GT.2000) nmax1=2000
  768.                   IF (nmax1.GT.nmax0) THEN
  769.                      nmax0=nmax1
  770.                      iter2=0
  771.                      recal1=100.D0
  772.                      GOTO 93
  773.                   ENDIF
  774.                ENDIF
  775. *
  776.             ENDIF
  777. *
  778.          ENDIF
  779. *
  780. *  Mise à jour des grandeurs après écoulement plastique
  781. *============================================================
  782. *
  783. *---> Variables internes mises à jour
  784. *
  785.          IF (Stest.GT.Gmin0) EPSPT=EPSPT+dlam0
  786.          EPSMT=EPSMT+DE0*dlam0
  787.          FNST=FNST+BB0*(H0-DS1*DF0)*dlam0
  788.          FNET=FNET+DD0*DE0*dlam0
  789.          FT=FT+DF0*dlam0
  790.          IF (FT.LT.Fmin0) FT=0.D0
  791. *
  792. *---> Contraintes mises à jour
  793. *
  794.          SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/RHOT/(1.D0-F0)-
  795.      .    1.D0)
  796.          DO 601 I=1,6
  797.             A=0.D0
  798.             IF (I.LE.3) A=1.D0
  799.             IF (Stest.LE.Gmin0) THEN
  800.                DEVT(I)=0.D0
  801.             ELSE
  802.                DEVT(I)=DEVT(I)-3.D0*G0*DEVT(I)*dlam0/Stest
  803.             ENDIF
  804.             SIGT(I)=DEVT(I)+(SMT*A)
  805.  601     CONTINUE
  806. *
  807. *---> Incrément de déformations plastiques mise à jour
  808. *
  809.          treps1=(SMT/(1.D0-FT)-SMT0/(1.D0-VAR0(2)))/(3.D0*XK0)
  810.          DO 610 I=1,6
  811.             A=0.D0
  812.             IF (I.LE.3) A=1.D0
  813.             depel0=(DEVT(I)-RSIG0(I)+(SMT0*A))/(2.D0*G0)
  814.             depel0=depel0+(treps1*A)
  815.             DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  816.             IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  817.                DEFT(3)=rdeps3+RDEPS(3)-depel0
  818.             ENDIF
  819.  610     CONTINUE
  820. *
  821. *---> Contrainte équivalente Stest mise à jour
  822. *
  823.          STEST2=1.5D0*PROCON(DEVT,DEVT,6)
  824.          IF (STEST2.LE.1.D-10) STEST2=0.D0
  825.          Stest= SQRT(STEST2)
  826. *
  827. *---> Limite d'élasticité lue sur la courbe d'écrouissage
  828. *     mise à jour
  829. *
  830. *   - Recherche des valeurs de déformations plastiques équivalentes
  831. *     encadrant la valeur courante
  832. *
  833.          EPSM0=TRAC(2*NCOURB)
  834.          IF (EPSMT.GE.EPSM0) THEN
  835.             Y1=TRAC(2*NCOURB-1)
  836.             Y2=TRAC(2*NCOURB-1)
  837.             EPS1=EPSM0
  838.             EPS2=2.D0*EPSM0
  839.          ELSE
  840.             DO 700 I=1,(NCOURB-1)
  841.                EPS1=TRAC(2*I)
  842.                EPS2=TRAC(2*(I+1))
  843.                IF ((EPSMT.GE.EPS1).AND.(EPSMT.LT.EPS2)) THEN
  844.                   Y1=TRAC(2*I-1)
  845.                   Y2=TRAC(2*(I+1)-1)
  846.                   GOTO 97
  847.                ENDIF
  848.  700        CONTINUE
  849.          ENDIF
  850. *
  851. *   - Limite d'élasticité mise à jour
  852. *
  853.  97      H1=(Y2-Y1)/(EPS2-EPS1)
  854.          Ytest=(H1*(EPSMT-EPS2))+Y2
  855. *
  856. *---> Mise à jour des coefficients de nucléation
  857. *
  858.          IF ((Ytest+SMT).GE.VARF(3)) THEN
  859.             BB0=EXP(-.5D0*((Ytest+SMT-SIGN0)/SNS0)**2.D0)
  860.             BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0))
  861.          ELSE
  862.             BB0=0.D0
  863.          ENDIF
  864.          IF (EPSMT.GE.VARF(4)) THEN
  865.             DD0=EXP(-.5D0*((EPSMT-EPSN0)/SNE0)**2.D0)
  866.             DD0=DD0*FNE0/(SNE0*(SQRT(2.D0*XPI)))
  867.          ELSE
  868.             DD0=0.D0
  869.          ENDIF
  870. *
  871. *---> Fonction d'endommagement
  872. *
  873.          IF (FT.LE.FC0) THEN
  874.             FST0=FT
  875.             ddel0=1.D0
  876.          ELSE
  877.             FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0)
  878.             ddel0=(FU0-FC0)/(FF0-FC0)
  879.          ENDIF
  880.          IF (FST0.GT.FU0) FST0=FU0
  881. *
  882. *---> Terme d'endommagement
  883. *
  884.          SMTM=MAX(SMT,0.D0)
  885.          SMTM=SMT
  886.          Gtest=1.D0+(Q3*FST0*FST0)-
  887.      .        (2.D0*Q1*FST0*COSH(3.D0*Q2*SMTM/(2.D0*Ytest)))
  888. *
  889. *---> Fonction de charge test PHIT mise à jour
  890. *
  891.          PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest
  892.          PHI0=ABS(PHIT)
  893. *
  894. *---> Cas particulier des matériaux totalement endommagés
  895. *
  896.  92      IF (FT.GE.FF0) THEN
  897. *
  898.             FT=FF0+Fmin0
  899. *
  900.             SMT=0.D0
  901.             treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/(3.D0*XK0)
  902.             DO 13 I=1,6
  903.                SIGT(I)=0.D0
  904.                A=0.D0
  905.                IF (I.LE.3) A=1.D0
  906.                depel0=SIGT(I)-RSIG0(I)
  907.                DEVT(I)=(depel0-((SMT-SMT0)*A))/(2.D0*G0)
  908.                depel0=DEVT(I)+(treps1*A)
  909.                DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  910.                IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  911.                   DEFT(3)=rdeps3+RDEPS(3)-depel0
  912.                ENDIF
  913.  13         CONTINUE
  914.             treps1=(DEFT(1)+DEFT(2)+DEFT(3))/3.D0
  915.             DO 15 I=1,6
  916.                DEVT(I)=DEFT(I)
  917.                IF (I.LE.1) DEVT(I)=DEVT(I)-treps1
  918.  15         CONTINUE
  919.             IF (treps1.LE.0.D0) THEN
  920.                DO 17 I=1,3
  921.                   DEFT(I)=DEVT(I)
  922.  17            CONTINUE
  923.             ENDIF
  924.             dlam0=SQRT((2.D0/3.D0)*PROCON(DEVT,DEVT,6))
  925.             EPSPT=EPSPT+dlam0
  926. **        RHOT=RMIN0
  927. *
  928.             Gtest=0.D0
  929.             PHIT=0.D0
  930.             PHI0=0.D0
  931. *
  932.          ENDIF
  933. *
  934. *---> Test de convergence ou itération suivante
  935. *
  936.          iter0=iter0+1
  937.          IF (iter0.LT.200) THEN
  938.             GOTO 99
  939.          ELSE
  940.             PHI0=0.D0
  941.             KERRE=460
  942.             GOTO 99
  943.          ENDIF
  944. *
  945.       ENDIF
  946. *
  947. *
  948. *   Vérification des sous incrémentations
  949. *
  950.       IF ((iter2.LT.nmax0).AND.(KERRE.EQ.0)) THEN
  951.          GOTO 95
  952.       ENDIF
  953. *
  954. *   Passage des deformations de cisaillement exprimées
  955. *   en déformations aux déformations de cisaillement exprimées
  956. *   en GAMA
  957. *
  958.       DO 117 I=4,6
  959.         RDEFP(I)=RDEFP(I)*2.D0
  960.  117  CONTINUE
  961. *
  962. *
  963. *  Passage a l'option de calcul pour les contraintes
  964. *=========================================================
  965. *
  966.       IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
  967.          IF (NSTRS .EQ. 6) THEN
  968. *
  969. *---> MASSIF 3D
  970. *
  971.             DO 170 I=1,NSTRS
  972.                SIGF(I)=RSIGF(I)
  973.                DEFP(I)=RDEFP(I)
  974.                DSIGT(I)=RDIGT(I)
  975.  170        CONTINUE
  976.          ELSE IF ( NSTRS .EQ. 4 ) THEN
  977. *
  978. *---> Calcul axisymétrique ou contraintes planes
  979. *
  980.             DO 180 I=1,NSTRS
  981.                SIGF(I)=RSIGF(I)
  982.                DEFP(I)=RDEFP(I)
  983.                DSIGT(I)=RDIGT(I)
  984.  180        CONTINUE
  985.          ENDIF
  986.       ENDIF
  987.  
  988.       RETURN
  989.       END
  990.  
  991.  
  992.  

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