Télécharger la plaquette Cast3M


fr en

Télécharger presto.eso

Retour à la liste

Numérotation des lignes : 

presto

  1. C PRESTO    SOURCE    CHAT      05/01/13    02:25:58     5004
  2.       SUBROUTINE PRESTO(DEPST,NSTRS,MFR1,IB,IGAU,
  3.      1     DSIGT,NCOMAT,SIG0,VAR0,XMAT,XCAR,NVARI,ICARA,
  4.      2     SIGF,VARF,DEFP,TETA1,TETA2,KERRE,DT0)
  5. *
  6. *_________________________________________________________________
  7. *
  8. *
  9. *    ENTREES :
  10. *    ---------
  11. *
  12. * DEPST = INCREMENT DE DEFORMATIONS TOTALES
  13. * NSTRS = NBRE DE COMPOSANTES DES DEFORMATIONS
  14. * NCOMAT= NBRE DE CARACTERISTIQUES MECANIQUES DU MATERIAU
  15. * MFR1  = NUMERO DE LA FORMULATION
  16. * IB    = NUMERO DE L ELEMENT COURANT
  17. * IGAU  = NUMERO DU POINT COURANT
  18. * NVARI = NBRE DE VARIABLES INTERNES
  19. * SIG0(NSTRS) = CONTRAINTES AU DEBUT DU PAS D'INTEGRATION
  20. * VAR0(NVARI) = VARIABLES INTERNES AU DEBUT DU PAS DE TEMPS
  21. * XMAT(NCOMAT) = CARACTERISTIQUES MECANIQUES DU MATERIAU
  22. * ICARA = NBRE DE CARACTERISTIQUES GEOMETRIQUES DES ELEMENTS FINIS
  23. * XCARA(ICARA) = CARACTERISTIQUES GEOMETRIQUES DES ELEMENTS FINIS
  24. * TETA1 = TEMPERATURE AU DEBUT DU PAS DE TEMPS
  25. * TETA2 = TEMPERATURE A LA FIN DU PAS DE TEMPS
  26. * DT0   = PAS DE TEMPS DU CALCUL
  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.       DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*),
  47.      &     VARF(*),DEFP(*),DSIGT(*)
  48.       DIMENSION RSIG0(6),RDEPS(6),RSIGF(6),RDEFP(6),RDIGT(6)
  49.       DIMENSION DEVT(6),ELT(6),SIGT(6),ST(6),DEFT(6)
  50.       DIMENSION DIJ(6),CIJ(6),CIJ2(6)
  51. *
  52. *  Adaptation de l'option de calcul vers le 3D massif de SIG0 a RSIG0
  53. *====================================================================
  54. *
  55.       IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
  56. *
  57. *---> 1 formulation massive
  58. *---> 2 formulation quasi incompressible
  59. *---> MASSIF 3D
  60. *
  61.          IF (NSTRS .EQ. 6) THEN
  62.             DO 110 I=1,NSTRS
  63.                RSIG0(I)=SIG0(I)
  64.                VARF(I+1)=VAR0(I+1)+DEPST(I)
  65.                RDEPS(I)=DEPST(I)
  66.  110        CONTINUE
  67.          ELSE IF (NSTRS.EQ.4.AND.((IFOUR.EQ.0)
  68.      &           .OR.(IFOUR.EQ.-1).OR.(IFOUR.EQ.-2))) THEN
  69. *
  70. *---> Calcul en mode deformations planes ou axisymetrique
  71. *     ou en contraintes planes
  72. *
  73.             DO 115 I=1,NSTRS
  74.                RSIG0(I)=SIG0(I)
  75.                VARF(I+1)=VAR0(I+1)+DEPST(I)
  76.                RDEPS(I)=DEPST(I)
  77.  115        CONTINUE
  78.             RSIG0(5)=0.D0
  79.             RSIG0(6)=0.D0
  80.             VARF(6)=VAR0(6)
  81.             VARF(7)=VAR0(7)
  82.             RDEPS(5)=0.D0
  83.             RDEPS(6)=0.D0
  84.          ENDIF
  85.       ELSE
  86.          KERRE = 99
  87.          GOTO 98
  88.       ENDIF
  89. *
  90. *   Passage des deformations de cisaillement exprimées
  91. *   en GAMA aux déformations de cisaillement exprimées
  92. *   en déformations
  93. *
  94.       DO 116 I=1,6
  95.          A=1.D0
  96.          IF (I.GT.3) A=2.D0
  97.          RDEPS(I)=RDEPS(I)/A
  98.          VAR0(I+1)=VAR0(I+1)/A
  99.          VARF(I+1)=VARF(I+1)/A
  100.          VAR0(I+8)=VAR0(I+8)/A
  101.          VARF(I+8)=VARF(I+8)/A
  102.  116  CONTINUE
  103. *
  104. *   Données du materiau
  105. *===========================================================
  106. *
  107. *
  108.       YOUNG=XMAT(1)
  109.       XNU=XMAT(2)
  110.       G00=YOUNG/(2.D0*(1.D0+XNU))
  111.       IF (IFOUR.EQ.-2) THEN
  112.          RHO0=XMAT(21)
  113.       ELSE
  114.          RHO0=XMAT(20)
  115.       ENDIF
  116.       TAU0=XMAT(5)
  117.       P0=XMAT(6)
  118.       S0=XMAT(7)
  119.       SINF0=XMAT(8)
  120.       XK00=XMAT(9)
  121.       G0=XMAT(10)
  122.       Y0=XMAT(11)
  123.       YINF0=XMAT(12)
  124.       Y1=XMAT(13)
  125.       Y2=XMAT(14)
  126.       BETA0=XMAT(15)
  127.       Gp1=XMAT(16)
  128.       Gt1=XMAT(17)
  129.       XMU0=XMAT(18)
  130.       TM0=XMAT(19)
  131. *
  132. *  Calculs préliminaires
  133. *===================================
  134. *
  135. *---> Nombre d'iterations internes maximales
  136. *
  137.       NMAX=200
  138.       iter00=0
  139. *
  140. *---> Trace de l'incrément de déformation totale
  141. *
  142.       IF (IFOUR.EQ.-2) THEN
  143.          VARF(4)=VAR0(11)-(XNU/(1.D0-XNU)
  144.      .        *(VARF(2)-VAR0(9)+VARF(3)-VAR0(10)))
  145.          RDEPS(3)=VARF(4)-VAR0(4)
  146.       ENDIF
  147.  96   treps1=RDEPS(1)+RDEPS(2)+RDEPS(3)
  148. *
  149. *---> Vitesse de deformation equivalente
  150. *
  151. *      DT0=5.D-7
  152.       DEPS2=RDEPS(1)*RDEPS(1)+RDEPS(2)*RDEPS(2)+RDEPS(3)*RDEPS(3)
  153.       DEPS2=DEPS2-RDEPS(1)*RDEPS(2)-RDEPS(2)*RDEPS(3)
  154.       DEPS2=DEPS2-RDEPS(3)*RDEPS(1)
  155.       DEPS3=RDEPS(4)*RDEPS(4)+RDEPS(5)*RDEPS(5)+RDEPS(6)*RDEPS(6)
  156.       DEPS2=DEPS2+3.D0*DEPS3
  157.       DEPS0=(DEPS2**(0.5D0))*2.D0/3.D0
  158.       DEPS0=DEPS0+(treps1/3.D0)
  159.       IF (ABS(DT0).LT.1.D-15) THEN
  160.          DEPS=0.D0
  161.       ELSE
  162.          DEPS=DEPS0/DT0
  163.       END IF
  164. *
  165. *---> Trace de la déformation totale
  166. *
  167.       treps0=VARF(2)+VARF(3)+VARF(4)
  168. *
  169. *---> Calcul de la compression a la fin du pas de temps
  170. *
  171.       RHO1=RHO0/(1.D0+treps0)
  172.       ETAnew=RHO1/RHO0
  173. *
  174. *---> Temperature de fusion a la fin (TM2)
  175. *     du pas de temps
  176. *
  177.       TM2=TM0*EXP(2.D0*XMU0*(1.D0-(1.D0/ETAnew)))
  178.       TM2=TM2/(ETAnew**(2.D0/3.D0))
  179. *
  180. *---> On divise le coefficient variable avec la température
  181. *     par la température de fusion
  182. *
  183.       XK0=XK00/TM2
  184. *
  185. *  Calculs des grandeurs élastiques test
  186. *====================================================
  187. *
  188. *---> Déformations
  189. *      -élastiques totales test ELT
  190. *      -trace des déformations élastiques totales TRELT
  191. *      -déviatoires élastiques test DEVT
  192. *      -incrément de déformation plastique DEFT
  193. *      -déformation plastique équivalente test EPSPT
  194. *
  195.       DO 100 I=1,6
  196.          ELT(I)=VARF(I+1)-VAR0(8+I)
  197.  100  CONTINUE
  198.       TRELT=ELT(1)+ELT(2)+ELT(3)
  199.       DO 200 I=1,6
  200.          A=0.D0
  201.          IF (I.LE.3) A=1.D0
  202.          DEVT(I)=ELT(I)-(1.D0*TRELT*A/3.D0)
  203.          DEFT(I)=0.D0
  204.  200  CONTINUE
  205.       EPSPT=VAR0(1)
  206. *
  207. *---> Coefficient proportionnel entre le module de
  208. *     cisaillement et le module de pression hydrostatique
  209. *
  210.       E1=2.D0*(1.D0+XNU)/(3.D0*(1.D0-(2.D0*XNU)))
  211. *
  212. *---> Données matériaux Gtest et XKtest
  213. *
  214.       pente1=E1*Gp1*TRELT/(ETAnew**(1.D0/3.D0))
  215.       Gtest=(G00+Gt1)/(1.D0+pente1)
  216.       XKtest=E1*Gtest
  217. *
  218. *---> Si on dépasse la température de fusion,
  219. *     le solide devient liquide
  220. *        -les caractéritiques matériaux sont constantes
  221. *        -G=0, donc le déviateur élastique est nul
  222. *        -la limite d'élasticité est constante et nulle donc:
  223. *        -la trace de la déformation plastique est nulle
  224. *
  225.       IF (TETA2.GT.TM2) THEN
  226.          Gtest=0.D0
  227.          DEVT(1)=0.D0
  228.          DEVT(2)=0.D0
  229.          DEVT(3)=0.D0
  230.          DEVT(4)=0.D0
  231.          DEVT(5)=0.D0
  232.          DEVT(6)=0.D0
  233.          TREPT=RDEPS(1)+RDEPS(2)+RDEPS(3)
  234.          TRELT=TREPT
  235.          DEFT(1)=RDEPS(1)-(1.D0*TREPT/3.D0)
  236.          DEFT(2)=RDEPS(2)-(1.D0*TREPT/3.D0)
  237.          DEFT(3)=RDEPS(3)-(1.D0*TREPT/3.D0)
  238.          DEFT(4)=RDEPS(4)
  239.          DEFT(5)=RDEPS(5)
  240.          DEFT(6)=RDEPS(6)
  241.          XKtest=VAR0(8)
  242.          DEPS0=PROCON(DEFT,DEFT,6)
  243.          DEPS0=(2.D0*DEPS0/3.D0)**.5D0
  244.          EPSPT=EPSPT+DEPS0
  245.       ENDIF
  246. *
  247. *---> Contraintes test
  248. *
  249.       Ptest=-1.D0*XKtest*TRELT
  250.       DO 300 I=1,6
  251.          ST(I)=2.D0*Gtest*DEVT(I)
  252.          A=0.D0
  253.          IF (I.LE.3) A=1.D0
  254.          SIGT(I)=ST(I)-(Ptest*A)
  255.  300  CONTINUE
  256. *
  257. *---> Contrainte équivalente test STEST
  258. *
  259.       STEST2=(SIGT(1)*SIGT(1))+(SIGT(2)*SIGT(2))
  260.       STEST2=STEST2+(SIGT(3)*SIGT(3))
  261.       STEST2=STEST2-(SIGT(1)*SIGT(2))-(SIGT(2)*SIGT(3))
  262.       STEST2=STEST2-(SIGT(3)*SIGT(1))
  263.       STEST3=(SIGT(4)*SIGT(4))+(SIGT(5)*SIGT(5))
  264.       STEST3=STEST3+(SIGT(6)*SIGT(6))
  265.       STEST2=STEST2+(3.D0*STEST3)
  266.       IF (STEST2.LT.0.D0) STEST2=0.D0
  267.       STEST=STEST2**(0.5D0)
  268. *
  269. *---> Calcul de la pulsation de Debye: W0 et de DTETA0
  270. *
  271.       W0=(Gtest/RHO1)**(.5D0)
  272. *
  273.       DTETA0=(1.D0/6.D0)*((4.D0/XPI)**(.5D0))*W0
  274. *
  275. *---> Calcul de la contrainte de saturation YSAT0
  276. *
  277.       IF (ABS(DT0).LT.1.D-15) THEN
  278.          YS1=SINF0
  279.       ELSE
  280.          FONC0=XK0*LOG(G0*DTETA0/DEPS)
  281.          ERF0=ERGAUS(FONC0)
  282.          YS1=S0-((S0-SINF0)*ERF0)
  283.       ENDIF
  284.       YS2=S0*((DEPS/(G0*DTETA0))**BETA0)
  285.       YSAT0=MAX(YS1,YS2)
  286. *
  287. *---> Calcul de la limite élastique YLIM0
  288. *
  289.       IF (ABS(DT0).LT.1.D-15) THEN
  290.          YLIM1=YINF0
  291.       ELSE
  292.          FONC0=XK0*LOG(G0*DTETA0/DEPS)
  293.          ERF0=ERGAUS(FONC0)
  294.          YLIM1=Y0-((Y0-YINF0)*ERF0)
  295.       ENDIF
  296.       YLIM2=Y1*((DEPS/(G0*DTETA0))**Y2)
  297.       YLIM0=MAX(YLIM1,MIN(YLIM2,YS2))
  298. *
  299. *---> Limite d'élasticité test : CAS CUBIQUE CENTRE C.C.
  300. *
  301.       DYSL0=ABS(YSAT0-YLIM0)
  302.       DYSL1=ABS(S0-YLIM0)
  303.       IF (P0.EQ.0.D0) THEN
  304. *
  305.          IF (DYSL0.LT.1.D-15) THEN
  306.             Yel0=YSAT0
  307.          ELSE
  308.             Yel0=-1.D0*TAU0*EPSPT
  309.             Yel0=EXP(Yel0/(YSAT0-YLIM0))
  310.             Yel0=YSAT0-((YSAT0-YLIM0)*Yel0)
  311.          ENDIF
  312.          Ytest=Yel0*Gtest
  313.          IF (Ytest.LE.0.D0) Ytest=0.D0
  314. *
  315.       ELSE
  316. *
  317. *---> Limite d'élasticité test : AUTRES CAS
  318. *
  319.          IF (DYSL0.LT.1.D-15) THEN
  320.             Yel0=YLIM0
  321.          ELSE IF (DYSL1.LT.1.D-15) THEN
  322.             Yel0=YLIM0
  323.          ELSE
  324.             Yel1=-1.D0*TAU0*EPSPT*P0
  325.             denom0=P0*(YSAT0-YLIM0)/(S0-YLIM0)
  326.             denom1=EXP(denom0)-1.D0
  327.             denom1=denom1*(S0-YLIM0)
  328.             Yel1=EXP(Yel1/denom1)
  329.             denom2=1.D0-(EXP(-1.D0*denom0))
  330.             Yel0=LOG(1.D0-(denom2*Yel1))
  331.             coeff0=(S0-YLIM0)/P0
  332.             Yel0=YSAT0+(coeff0*Yel0)
  333.          ENDIF
  334.          Ytest=Yel0*Gtest
  335.          IF (Ytest.LE.0.D0) Ytest=0.D0
  336. *
  337.       ENDIF
  338. *
  339. *---> Fonction de charge test PHIT
  340. *
  341.       PHIT=STEST-Ytest
  342.       IF (TETA2.GT.TM2) PHIT=0.D0
  343. *
  344. *
  345. *  Vérification du critère de plasticité
  346. *=========================================================
  347. *
  348. *---> Erreur admise
  349. *
  350.       ERR0=1.D-7*ABS(PHIT)
  351. *
  352. *---> Critère de plasticité
  353. *
  354.       PHI0=PHIT
  355.       iter0=0
  356. *
  357. *
  358.  97   IF (PHI0.LE.ERR0) THEN
  359. *
  360. *  On est élastique
  361. *=========================================================
  362. *
  363. *--------------------------------------------------------------
  364. *     Cas particulier des contraintes planes
  365. *--------------------------------------------------------------
  366. *
  367.          IF (IFOUR.EQ.-2) THEN
  368. *
  369. *---> On vérifie qu'on est en contraintes planes
  370. *
  371.             iter00=iter00+1
  372.             SIG3=ABS(SIGT(3))
  373.             IF (SIG3.GT.1.D-2) THEN
  374.                ELT(3)=-1.D0*XNU*(ELT(1)+ELT(2))/(1.D0-XNU)
  375.                VARF(4)=VAR0(11)+ELT(3)+DEFT(3)
  376.                RDEPS(3)=VARF(4)-VAR0(4)
  377.                IF (iter00.LE.NMAX) THEN
  378.                   GOTO 96
  379.                ELSE
  380.                   KERRE=460
  381.                   GOTO 98
  382.                ENDIF
  383.             ENDIF
  384.          ENDIF
  385. *
  386. *-----------------------------------------------------
  387. *    Fin du cas particulier des contraintes planes
  388. *-----------------------------------------------------
  389. *
  390.          DO 400 I=1,6
  391.             RSIGF(I)=SIGT(I)
  392.             RDEFP(I)=DEFT(I)
  393.             RDIGT(I)=RSIGF(I)-RSIG0(I)
  394.             VARF(I+8)=DEFT(I)+VAR0(I+8)
  395.  400     CONTINUE
  396.          VARF(1)=EPSPT
  397.          VARF(8)=XKtest
  398.          DEPST(3)=VARF(4)-VAR0(4)
  399. *
  400. *
  401.       ELSE
  402. *
  403. *  On plastifie
  404. *=========================================================
  405. *
  406. *---> Coefficient donnant la trace de l'incrément de
  407. *     déformation plastique
  408. *
  409.          tr0=Yel0*Gp1/(ETAnew**(1.D0/3.D0))
  410. *
  411. *---> Coefficient donnant la variation de pression due
  412. *     à l'incrément de déformation plastique
  413. *
  414.          coeff1=E1*Gp1*TRELT/(ETAnew**(1.D0/3.D0))
  415.          coeff1=XKtest*tr0/(1.D0+coeff1)
  416. *
  417. *---> Module d'écrouissage H = - d PHI / d EPSeq
  418. *
  419.          IF (P0.EQ.0.D0) THEN
  420.             IF (DYSL0.LT.1.D-15) THEN
  421.                H0=0.D0
  422.             ELSE
  423.                H0=TAU0*EXP(-1.D0*TAU0*EPSPT/(YSAT0-YLIM0))
  424.                H0=H0*Gtest
  425.             ENDIF
  426.          ELSE
  427.             IF ((DYSL0.LT.1.D-15).OR.(DYSL1.LT.1.D-15)) THEN
  428.                H0=0.D0
  429.             ELSE
  430.                H0=coeff0*denom2*P0*TAU0
  431.                fac0=denom2*EXP(-1.D0*P0*TAU0*EPSPT/denom1)
  432.                H0=H0/(denom1*(1.D0-fac0))
  433.                H0=H0*Gtest
  434.             ENDIF
  435.          ENDIF
  436. *
  437. *---> Matrice de la direction d'écoulement DIJ = d PHI / d SIG
  438. *
  439.          DO 500 I=1,6
  440.             A=0.D0
  441.             IF (I.LE.3) A=1.D0
  442.             DIJ(I)=3.D0*ST(I)/(STEST*2.D0)
  443.             DIJ(I)=DIJ(I)+(A*tr0/3.D0)
  444.  500     CONTINUE
  445. *
  446. *---> Matrice de la variation de contrainte due à l'incrément de
  447. *     déformation plastique CIJ = - ( Hooke . d PHI / d SIG)
  448. *
  449.          DO 600 I=1,6
  450.             A=0.D0
  451.             IF (I.LE.3) A=1.D0
  452.             CIJ(I)=3.D0*Gtest*ST(I)/STEST
  453.             CIJ(I)=CIJ(I)-(2.D0*Gp1*coeff1*DEVT(I)/(ETAnew**(1.D0/3.D0))
  454.      $           )
  455.             CIJ(I)=CIJ(I)+(coeff1*A)
  456.             CIJ2(I)=3.D0*Gp1*coeff1*ST(I)/((ETAnew**(1.D0/3.D0))*STEST)
  457.  600     CONTINUE
  458. *
  459. *---> Calcul du produit DIJ : CIJ
  460. *
  461.          PROD0=PROCON(DIJ,CIJ,6)
  462. *
  463. *---> Calcul de DIJ : CIJ2
  464. *
  465.          PROD1=PROCON(DIJ,CIJ2,6)
  466.          PROD11=ABS(PROD1)
  467. *
  468. *---> Calcul de l'incrément du paramètre plastique
  469. *     On résoud une équation du second degré
  470. *     Si le terme (PROD0+H0) est positif, dlam0 est
  471. *  clairement défini par la solution positive de l'équation
  472. *     Sinon, on prend la solution qui est de meme signe que
  473. *  PHIT (une solution est nécessairement
  474. *  de signe opposé a PHIT). Si aucune solution ne convient
  475. *  on prend celle qui a le plus petit module afin d'éviter de
  476. *  diverger.
  477. *
  478.          delta0=(PROD0+H0)*(PROD0+H0)
  479.          delta0=delta0+(4.D0*PHIT*PROD1)
  480.          IF ((delta0.LT.0.D0).OR.(PROD11.LT.1.D-15)) THEN
  481.             dlam0=PHIT/(PROD0+H0)
  482.          ELSE
  483.             IF ((PROD0+H0).GE.0.D0) THEN
  484.                dlam0=(delta0**(.5D0))/(2.D0*PROD1)
  485.                dlam0=dlam0-((PROD0+H0)/(2.D0*PROD1))
  486.             ELSE
  487.                dlam0=(delta0**(.5D0))/(2.D0*PROD1)
  488.                dlam0=dlam0-((PROD0+H0)/(2.D0*PROD1))
  489.                IF ((dlam0*PHIT).LE.0.D0) THEN
  490.                   dlam1=-1.D0*(delta0**(.5D0))/(2.D0*PROD1)
  491.                   dlam1=dlam1-((PROD0+H0)/(2.D0*PROD1))
  492.                   IF (ABS(dlam1).LT.ABS(dlam0)) dlam0=dlam1
  493.                ENDIF
  494.             ENDIF
  495.          ENDIF
  496. *
  497. *
  498. *  Mise à jour des grandeurs après écoulement plastique
  499. *============================================================
  500. *
  501. *---> Déformation plastique équivalente
  502. *
  503.          EPSPT=EPSPT+dlam0
  504. *
  505. *---> Pression mise à jour
  506. *
  507.          Ptest=Ptest+(coeff1*dlam0)
  508. *
  509. *---> Contraintes
  510. *
  511.          DO 800 I=1,6
  512.             A=0.D0
  513.             IF (I.LE.3) A=1.D0
  514.             SIGT(I)=SIGT(I)-(CIJ(I)*dlam0)-(CIJ2(I)*dlam0*dlam0)
  515. *         SIGT(I)=SIGT(I)-(CIJ(I)*dlam0)
  516.             ST(I)=SIGT(I)+(Ptest*A)
  517.  800     CONTINUE
  518. *
  519. *---> Contrainte équivalente mise à jour STEST
  520. *
  521.          STEST2=(SIGT(1)*SIGT(1))+(SIGT(2)*SIGT(2))
  522.          STEST2=STEST2+(SIGT(3)*SIGT(3))
  523.          STEST2=STEST2-(SIGT(1)*SIGT(2))-(SIGT(2)*SIGT(3))
  524.          STEST2=STEST2-(SIGT(3)*SIGT(1))
  525.          STEST3=(SIGT(4)*SIGT(4))+(SIGT(5)*SIGT(5))
  526.          STEST3=STEST3+(SIGT(6)*SIGT(6))
  527.          STEST2=STEST2+(3.D0*STEST3)
  528.          IF (STEST2.LE.0.D0) STEST2=0.D0
  529.          STEST=STEST2**(0.5D0)
  530. *
  531. *---> Calcul des données matériaux mises à jour
  532. *
  533.          Gtest=G00+(Gp1*Ptest/(ETAnew**(1.D0/3.D0)))+Gt1
  534.          XKtest=E1*Gtest
  535. *
  536. *
  537. *  Mise à jour des déformations après écoulement plastique
  538. *============================================================
  539. *
  540. *---> Trace des déformations élastiques mises à jour
  541. *
  542.          TRELT=-1.D0*Ptest/XKtest
  543. *
  544. *---> Déviateur des déformations élastiques
  545. *
  546.          DO 900 I=1,6
  547.             A=0.D0
  548.             IF (I.LE.3) A=1.D0
  549.             DEVT(I)=ST(I)/(2.D0*Gtest)
  550.             IF (IFOUR.EQ.-2) ELT(I)=DEVT(I)+(A*TRELT/3.D0)
  551.  900     CONTINUE
  552. *
  553. *---> Incrément de déformation plastique au cours du pas
  554. *
  555.          DO 700 I=1,6
  556.             A=0.D0
  557.             IF (I.LE.3) A=1.D0
  558.             DEFT(I)=VARF(I+1)-VAR0(I+8)
  559.             DEFT(I)=DEFT(I)-DEVT(I)-(1.D0*TRELT*A/3.D0)
  560.  700     CONTINUE
  561. *
  562. *---> Mise à jour de la pulsation de Debye: W0 et de DTETA0
  563. *
  564.          W0=(Gtest/RHO1)**(.5D0)
  565. *
  566.          DTETA0=(1.D0/6.D0)*((4.D0/XPI)**(.5D0))*W0
  567. *
  568. *---> Mise à jour de la contrainte de saturation YSAT0
  569. *
  570.          IF (ABS(DT0).LT.1.D-15) THEN
  571.             YS1=SINF0
  572.          ELSE
  573.             ERF0=ERGAUS(FONC0)
  574.             YS1=S0-((S0-SINF0)*ERF0)
  575.          ENDIF
  576.          YS2=S0*((DEPS/(G0*DTETA0))**BETA0)
  577.          YSAT0=MAX(YS1,YS2)
  578. *
  579. *---> Mise à jour de la limite élastique YLIM0
  580. *
  581.          IF (ABS(DT0).LT.1.D-15) THEN
  582.             YLIM1=YINF0
  583.          ELSE
  584.             FONC0=XK0*LOG(G0*DTETA0/DEPS)
  585.             ERF0=ERGAUS(FONC0)
  586.             YLIM1=Y0-((Y0-YINF0)*ERF0)
  587.          ENDIF
  588.          YLIM2=Y1*((DEPS/(G0*DTETA0))**Y2)
  589.          YLIM0=MAX(YLIM1,MIN(YLIM2,YS2))
  590. *
  591. *---> limite d'élasticité mise à jour
  592. *
  593.          DYSL0=ABS(YSAT0-YLIM0)
  594.          DYSL1=ABS(S0-YLIM0)
  595.          IF (P0.EQ.0.D0) THEN
  596.             IF (DYSL0.LT.1.D-15) THEN
  597.                Yel0=YSAT0
  598.             ELSE
  599.                Yel0=-1.D0*TAU0*EPSPT
  600.                Yel0=EXP(Yel0/(YSAT0-YLIM0))
  601.                Yel0=YSAT0-((YSAT0-YLIM0)*Yel0)
  602.             ENDIF
  603.             Ytest=Yel0*Gtest
  604.             IF (Ytest.LE.0.D0) Ytest=0.D0
  605.          ELSE
  606.             IF (DYSL0.LT.1.D-15) THEN
  607.                Yel0=YLIM0
  608.             ELSE IF (DYSL1.LT.1.D-15) THEN
  609.                Yel0=YLIM0
  610.             ELSE
  611.                Yel1=-1.D0*TAU0*EPSPT*P0
  612.                denom0=P0*(YSAT0-YLIM0)/(S0-YLIM0)
  613.                denom1=EXP(denom0)-1.D0
  614.                denom1=denom1*(S0-YLIM0)
  615.                Yel1=EXP(Yel1/denom1)
  616.                denom2=1.D0-(EXP(-1.D0*denom0))
  617.                Yel0=LOG(1.D0-(denom2*Yel1))
  618.                coeff0=(S0-YLIM0)/P0
  619.                Yel0=YSAT0+(coeff0*Yel0)
  620.             ENDIF
  621.             Ytest=Yel0*Gtest
  622.          ENDIF
  623. *
  624. *---> Nouvelle fonction de charge
  625. *
  626.          PHIT=STEST-Ytest
  627.          PHI0=ABS(PHIT)
  628. *
  629. *---> Test de convergence ou itération suivante
  630. *
  631.          iter0=iter0+1
  632. *
  633.          IF (iter0.LT.NMAX) THEN
  634.             GOTO 97
  635.          ELSE
  636.             KERRE=460
  637.             GOTO 98
  638.          ENDIF
  639. *
  640.       ENDIF
  641. *
  642. *   Passage des deformations de cisaillement exprimées
  643. *   en déformations aux déformations de cisaillement exprimées
  644. *   en GAMA
  645. *
  646.       DO 117 I=1,6
  647.          A=1.D0
  648.          IF (I.GT.3) A=2.D0
  649.          RDEFP(I)=RDEFP(I)*A
  650.          VAR0(I+1)=VAR0(I+1)*A
  651.          VARF(I+1)=VARF(I+1)*A
  652.          VAR0(I+8)=VAR0(I+8)*A
  653.          VARF(I+8)=VARF(I+8)*A
  654.  117  CONTINUE
  655. *
  656. *  Passage a l'option de calcul pour les contraintes
  657. *=========================================================
  658. *
  659.       IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31)  THEN
  660.          IF (NSTRS .EQ. 6) THEN
  661. *
  662. *---> MASSIF 3D
  663. *
  664.             DO 170 I=1,NSTRS
  665.                SIGF(I)=RSIGF(I)
  666.                DEFP(I)=RDEFP(I)
  667.                DSIGT(I)=RDIGT(I)
  668.  170        CONTINUE
  669.          ELSE IF ( NSTRS .EQ. 4 ) THEN
  670. *
  671. *---> Calcul axisymétrique ou contraintes planes
  672. *
  673.             DO 180 I=1,NSTRS
  674.                SIGF(I)=RSIGF(I)
  675.                DEFP(I)=RDEFP(I)
  676.                DSIGT(I)=RDIGT(I)
  677.  180        CONTINUE
  678.          ENDIF
  679.       ENDIF
  680. *
  681. *
  682.  98   CONTINUE
  683.       RETURN
  684. *
  685.       END
  686.  
  687.  
  688.  
  689.  
  690.  
  691.  
  692.  
  693.  
  694.  
  695.  
  696.  
  697.  
  698.  
  699.  

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