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

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