Télécharger la plaquette Cast3M


fr en

Télécharger pakgon.eso

Retour à la liste

Numérotation des lignes : 

pakgon

  1. C PAKGON    SOURCE    CHAT      05/01/13    02:08:47     5004
  2.       SUBROUTINE PAKGON(XMAT,SIG0,VAR0,SIGF,DEFF,VARF,
  3.      .                  DEPST,KERRE,IB,IGAU,PCT1,
  4.      .                  SUCF0,RSIG0)
  5. *
  6. *=============================================================
  7. *  Modèle d'argile de PAKZAD: surface de gonflement
  8. *=============================================================
  9. *
  10. *  Entrées
  11. *-------------------------------------------------------------
  12. *
  13. *  XMAT   : caractéristiques matériau
  14. *  SIG0   : contraintes initiales TEST
  15. *  VAR0   : variables internes initiales TEST
  16. *  DEPST  : incrément des déformations totales
  17. *  IB     : numéro de l'élément courrant
  18. *  IGAU   : numéro du point de Gauss courrant
  19. *  PCT1   : pression de consolidation "fictive" de l'état saturé
  20. *  SUCF0  : succion finale
  21. *  RSIG0  : contraintes au début du pas
  22. *
  23. *-------------------------------------------------------------
  24. *
  25. *
  26. *  Sorties
  27. *-------------------------------------------------------------
  28. *
  29. *  SIGF   : contraintes finales
  30. *  DEFF   : incrémént de déformation plastique final
  31. *  VARF   : variables internes finales
  32. *  KERRE  : message d'erreur
  33. *
  34. *-------------------------------------------------------------
  35. *  ICI ON PROGRAMME EN FORTRAN PUR
  36. *
  37.       IMPLICIT INTEGER(I-N)
  38.       IMPLICIT REAL*8(A-H,O-Z)
  39. *
  40.       DIMENSION XMAT(*),SIG0(*),VAR0(*),DEPST(*)
  41.       DIMENSION SIGF(*),DEFF(*),VARF(*),RSIG0(*)
  42.       DIMENSION SIGT(6),DEVT(6),DEV(6),DEV0(6),DEFP(6)
  43. *
  44. *
  45. *   Données du materiau
  46. *===========================================================
  47. *
  48.       YOUN0=XMAT(1)
  49.       XNU0=XMAT(2)
  50.       XN0=XMAT(5)
  51.       XKA0=XMAT(6)
  52.       XGA0=XMAT(7)
  53.       XPA0=XMAT(8)
  54.       XPC0=XMAT(9)
  55.       XM0=XMAT(10)
  56.       XBET1=XMAT(11)
  57.       XA0=XMAT(12)
  58.       XPCR0=XMAT(13)
  59.       XS0=XMAT(14)
  60.       XM1=XMAT(15)
  61.       XM2=XMAT(16)
  62.       XM3=XMAT(17)
  63.       XBET2=XMAT(18)
  64.       XTAU0=XMAT(19)
  65.       dif0=ABS(XN0-1.D0)
  66. *
  67. *---> Initialisation du rayon initial de l'ellipse
  68. *
  69.       IF (VAR0(5).LE.1.D-10) VAR0(5)=XA0
  70.       IF (VAR0(6).LE.1.D-10) VAR0(6)=XA0
  71. *
  72. *---> Charge ou décharge ?
  73. *     . NCHAR0 = 0 : décharge
  74. *     . NCHAR0 = 1 : recharge
  75. *
  76.       TRELT=DEPST(1)+DEPST(2)+DEPST(3)
  77.       NCHAR0=0
  78.       IF (TRELT.LE.0.D0) NCHAR0=1
  79. *
  80. *---> Contraintes test, incrément de déformation plastique test,
  81. *     déviateur initial et pression initiale
  82. *
  83.       PRES00=-1.D0/3.D0*(RSIG0(1)+RSIG0(2)+RSIG0(3))
  84.       DO 05 I=1,6
  85.         A=0.D0
  86.         IF (I.LE.3) A=1.D0
  87.         SIGT(I)=SIG0(I)
  88.         DEFP(I)=0.D0
  89.         DEVT(I)=RSIG0(I)+PRES00*A
  90.  05   CONTINUE
  91. *
  92. *---> Cas où XN0 n'est pas égal à 1
  93. *   Les constantes de raideur sont prises au début du pas
  94. *
  95.       IF  (dif0.GT.1.D-5) THEN
  96.         XKA1=XKA0*XPA0*((PRES00/XPA0)**XN0)
  97.         XGA1=XGA0*XPA0*((PRES00/XPA0)**XN0)
  98.       ENDIF
  99. *
  100. *---> Contrainte équivalente initiale
  101. *
  102.       Stest=PROCON(DEVT,DEVT,6)
  103.       IF (Stest.LE.(1.D-10*YOUN0)) Stest=0.D0
  104.       St0=(3.D0*Stest/2.D0)**(.5D0)
  105. *
  106. *---> Calcul . du terme de contrainte équivalente
  107. *            . du centre de la surface
  108. *            . du rayon
  109. *
  110.       PRES1=-1.D0/3.D0*(SIGT(1)+SIGT(2)+SIGT(3))
  111.       IF (NCHAR0.EQ.1) THEN
  112. *
  113. *  On charge
  114. *-----------------------------------------------------
  115. *
  116. *  Centre test et au début du pas
  117. *
  118.         DO 10 I=1,6
  119.           A=0.D0
  120.           IF (I.LE.3) A=1.D0
  121.           DEVT(I)=SIGT(I)+PRES1*A
  122.           DEV(I)=(DEVT(I)-VAR0(14+I))/2.D0
  123.           DEV0(I)=(RSIG0(I)+PRES00*A-VAR0(14+I))/2.D0
  124.  10     CONTINUE
  125.         APT=(PRES1+VAR0(8))/2.D0
  126.         APT0=(PRES00+VAR0(8))/2.D0
  127. *
  128. *  Rayon test et au début du pas
  129. *
  130.         AIS=VAR0(6)
  131.         AIS0=AIS
  132. *
  133. *  Contrainte équivalente du point fixe
  134. *
  135.         Stest0=PROCON(VAR0(15),VAR0(15),6)
  136.         IF (Stest0.LE.(1.D-10*YOUN0)) Stest0=0.D0
  137.         Stest0=(3.D0*Stest0/2.D0)**(.5D0)
  138.       ELSE
  139. *
  140. *  On décharge
  141. *-----------------------------------------------------
  142. *
  143. *  Centre test et au début du pas
  144. *
  145.         DO 11 I=1,6
  146.           A=0.D0
  147.           IF (I.LE.3) A=1.D0
  148.           DEVT(I)=SIGT(I)+PRES1*A
  149.           DEV(I)=(DEVT(I)-VAR0(8+I))/2.D0
  150.           DEV0(I)=(RSIG0(I)+PRES00*A-VAR0(8+I))/2.D0
  151.  11     CONTINUE
  152.         APT=(PRES1+VAR0(7))/2.D0
  153.         APT0=(PRES00+VAR0(7))/2.D0
  154. *
  155. *  Rayon test et au début du pas
  156. *
  157.         AIS=VAR0(5)
  158.         AIS0=AIS
  159. *
  160. *  Contrainte équivalente du point fixe
  161. *
  162.         Stest0=PROCON(VAR0(9),VAR0(9),6)
  163.         IF (Stest0.LE.(1.D-10*YOUN0)) Stest0=0.D0
  164.         Stest0=(3.D0*Stest0/2.D0)**(.5D0)
  165.       ENDIF
  166. *
  167. *---> Terme déviatoire de la distance au centre de la surface test
  168. *
  169.       AJ2=PROCON(DEV,DEV,6)
  170.       IF (AJ2.LE.(1.D-10*YOUN0)) AJ2=0.D0
  171.       AJ2=(3.D0/2.D0*AJ2)**(.5D0)
  172. *
  173. *---> Terme déviatoire de la distance au centre de la surface au
  174. *  début du pas
  175. *
  176.       AJ20=PROCON(DEV0,DEV0,6)
  177.       IF (AJ20.LE.(1.D-10*YOUN0)) AJ20=0.D0
  178.       AJ20=(3.D0/2.D0*AJ20)**(.5D0)
  179. *
  180. *---> Pressions critiques au début du pas
  181. *
  182.       APT10=2.D0*APT0-PRES00
  183.       PC2=(St0*St0+XM0*XM0*PRES00*PRES00)/(XM0*XM0*PRES00)
  184.       PCD1=(Stest0*Stest0+XM0*XM0*APT10*APT10)/(XM0*XM0*APT10)
  185. *
  186. *---> Initialisation de la pression initiale du "CENTRE"
  187. *
  188.       IF (NCHAR0.EQ.1) THEN
  189.         IF (ABS(VAR0(22)).LE.1.D-10) THEN
  190.           AP1=APT
  191.           AP2=0.D0
  192.         ELSE
  193.           AP1=VAR0(22)
  194.           AP2=AP1
  195.         ENDIF
  196.       ELSE
  197.         IF (ABS(VAR0(21)).LE.1.D-10) THEN
  198.           AP1=APT
  199.           AP2=0.D0
  200.         ELSE
  201.           AP1=VAR0(21)
  202.           AP2=AP1
  203.         ENDIF
  204.       ENDIF
  205. *
  206. *---> Terme d'écrouissage BETA2 au début du pas
  207. *     Terme ALP0 non associatif de l'écoulemment
  208. *
  209.       IF (NCHAR0.EQ.1) THEN
  210.        BETA2=(XPCR0/PC2+1.D0)*PRES00*2.D0*VAR0(4)
  211.        BETA2=XBET1*BETA2/(3.D0*PCD1*PC2)
  212.        ALP0=(1.D0-(ABS(APT10-AP1)/XA0))
  213.        ALP0=ALP0*(1.D0-1.5D0*PC2/VAR0(4))
  214.       ELSE
  215.        BETA2=XBET1*(XPCR0/PCD1+1.D0)*PC2/PRES00
  216.        ALP0=(1.D0-(ABS(APT10-AP1)/XA0))
  217.        ALP0=ALP0*(1.5D0*PC2/PCD1-.5D0)
  218.       ENDIF
  219.       BETA2=BETA2*EXP(XTAU0*SUCF0)
  220.       ALP0=-1.D0*ALP0
  221. *
  222. *---> Termes d'accroissement du rayon de l'ellipse
  223. *
  224.       XB1=PRES00*BETA2*(PRES00-APT0)/(2.D0*AIS0)
  225.       XB2=PRES00*BETA2*AJ20/(2.D0*AIS0)
  226. *
  227. *---> Surface de charge de gonflement
  228. *
  229.       PHIT=AJ2*AJ2+XM0*XM0*(PRES1-APT+AIS)*(PRES1-APT-AIS)
  230. *
  231. *
  232. *
  233. *  Vérification du critère de plasticité
  234. *=========================================================
  235. *
  236. *---> Erreur admise
  237. *
  238.       ERR0=1.D-7*ABS(PHIT)
  239.       ERR0=MAX(ERR0,1.D-10*YOUN0)
  240. *
  241. *---> Critère de plasticité
  242. *
  243.       PHI0=PHIT
  244.       iter0=0
  245. *
  246.  99   IF (PHI0.LE.ERR0) THEN
  247. *
  248. *  On est élastique par rapport à la surface de gonflement
  249. *=========================================================
  250. *
  251.       IF (NCHAR0.EQ.1) THEN
  252. *
  253. *  On charge
  254. *
  255.         DO 20 I=1,6
  256.           SIGF(I)=SIGT(I)
  257.           DEFF(I)=DEFP(I)
  258.           VARF(8+I)=DEVT(I)
  259.           VARF(14+I)=VAR0(14+I)
  260.  20     CONTINUE
  261.         VARF(5)=XA0
  262.         VARF(6)=AIS
  263.         VARF(7)=PRES1
  264.         VARF(8)=VAR0(8)
  265.         VARF(21)=0.D0
  266.         VARF(22)=AP2
  267.       ELSE
  268. *
  269. *  On décharge
  270. *
  271.         DO 21 I=1,6
  272.           SIGF(I)=SIGT(I)
  273.           DEFF(I)=DEFP(I)
  274.           VARF(8+I)=VAR0(8+I)
  275.           VARF(14+I)=DEVT(I)
  276.  21     CONTINUE
  277.         VARF(5)=AIS
  278.         VARF(6)=XA0
  279.         VARF(7)=VAR0(7)
  280.         VARF(8)=PRES1
  281.         VARF(21)=AP2
  282.         VARF(22)=0.D0
  283.       ENDIF
  284. *
  285.       ELSE
  286. *
  287. *  On est plastique par rapport à la surface de gonflement
  288. *=========================================================
  289. *
  290. *---> Pression initiale du "CENTRE"
  291. *
  292.       IF (NCHAR0.EQ.1) THEN
  293.         IF (ABS(VAR0(22)).LE.1.D-10) THEN
  294.           AP1=APT
  295.           AP2=AP1
  296.         ENDIF
  297.       ELSE
  298.         IF (ABS(VAR0(21)).LE.1.D-10) THEN
  299.           AP1=APT
  300.           AP2=AP1
  301.         ENDIF
  302.       ENDIF
  303. *
  304. *---> Condition de consistance
  305. *
  306.       IF (dif0.LE.1.D-5) THEN
  307. *
  308. *  Cas où XN0=1
  309. *-------------------------------------------------------------
  310. *
  311.         tr0=2.D0*XM0*XM0*(PRES1-APT)
  312.         DEQ=ABS(2.D0*AJ2-ALP0*tr0)
  313.         DF0=3.D0*XGA0*PRES00*AJ2*DEQ
  314.         DF1=XKA0*PRES1*tr0*tr0/2.D0
  315.         DF2=(XB1*tr0+XB2*DEQ)*2.D0*XM0*XM0*AIS
  316. *
  317.         denom0=DF0+DF1+DF2
  318. *
  319.         dlam0=PHIT/denom0
  320. *
  321. *---> Mise à jour des contraintes
  322. *
  323.         PRES1=PRES1*EXP(-1.D0*XKA0*tr0*dlam0)
  324.         DO 30 I=1,6
  325.           A=0.D0
  326.           IF (I.LE.3) A=1.D0
  327.           IF (AJ2.GT.(YOUN0*1.D-10)) THEN
  328.            DEVT(I)=DEVT(I)-3.D0*XGA0*PRES00*DEV(I)/AJ2*DEQ*dlam0
  329.           ELSE
  330.            DEVT(I)=DEVT(I)-6.D0*XGA0*PRES00*DEV(I)*dlam0
  331.           ENDIF
  332.           SIGT(I)=DEVT(I)-PRES1*A
  333.  30     CONTINUE
  334. *
  335. *---> Mise à jour du rayon de l'ellipse
  336. *
  337.         AIS=AIS+tr0*XB1*dlam0+XB2*DEQ*dlam0
  338. *
  339. *---> Mise à jour du centre de l'ellipse
  340. *
  341.         IF (NCHAR0.EQ.1) THEN
  342.           APT=(PRES1+VAR0(8))/2.D0
  343.           DO 35 I=1,6
  344.             DEV(I)=(DEVT(I)-VAR0(14+I))/2.D0
  345.  35       CONTINUE
  346.         ELSE
  347.           APT=(PRES1+VAR0(7))/2.D0
  348.           DO 36 I=1,6
  349.             DEV(I)=(DEVT(I)-VAR0(8+I))/2.D0
  350.  36       CONTINUE
  351.         ENDIF
  352.         AJ2=PROCON(DEV,DEV,6)
  353.         IF (AJ2.LE.(1.D-10*YOUN0)) AJ2=0.D0
  354.         AJ2=(3.D0/2.D0*AJ2)**(.5D0)
  355. *
  356. *---> Mise à jour des déformations élastiques et plastiques
  357. *
  358.         TRELT=-1.D0*LOG(PRES1/PRES00)/XKA0
  359.         DO 40 I=1,6
  360.           A=0.D0
  361.           IF (I.LE.3) A=1.D0
  362.           DEFP(I)=(DEVT(I)-(RSIG0(I)+PRES00*A))
  363.      .            /(2.D0*XGA0*PRES00)
  364.           DEFP(I)=DEPST(I)-DEFP(I)-TRELT*A/3.D0
  365.  40     CONTINUE
  366. *
  367.       ELSE
  368. *
  369. *   Autres cas
  370. *-------------------------------------------------------------
  371. *
  372.         tr0=2.D0*XM0*XM0*(PRES1-APT)
  373.         DEQ=ABS(2.D0*AJ2-ALP0*tr0)
  374.         DF0=3.D0*XGA1*AJ2*DEQ
  375.         DF1=XKA1*tr0*tr0/2.D0
  376.         DF2=(XB1*tr0+XB2*DEQ)*2.D0*XM0*XM0*AIS
  377. *
  378.         denom0=DF0+DF1+DF2
  379. *
  380.         dlam0=PHIT/denom0
  381. *
  382. *---> Mise à jour des contraintes
  383. *
  384.         PRES1=PRES1-1.D0*XKA1*tr0*dlam0
  385.         DO 31 I=1,6
  386.           A=0.D0
  387.           IF (I.LE.3) A=1.D0
  388.           IF (AJ2.GT.(YOUN0*1.D-10)) THEN
  389.            DEVT(I)=DEVT(I)-3.D0*XGA1*DEV(I)/AJ2*DEQ*dlam0
  390.           ELSE
  391.            DEVT(I)=DEVT(I)-6.D0*XGA1*DEV(I)*dlam0
  392.           ENDIF
  393.           SIGT(I)=DEVT(I)-PRES1*A
  394.  31     CONTINUE
  395. *
  396. *---> Mise à jour du rayon de l'ellipse
  397. *
  398.         AIS=AIS+tr0*XB1*dlam0+XB2*DEQ*dlam0
  399. *
  400. *---> Mise à jour du centre de l'ellipse
  401. *
  402.         IF (NCHAR0.EQ.1) THEN
  403.           APT=(PRES1+VAR0(8))/2.D0
  404.           DO 350 I=1,6
  405.             DEV(I)=(DEVT(I)-VAR0(14+I))/2.D0
  406.  350       CONTINUE
  407.         ELSE
  408.           APT=(PRES1+VAR0(7))/2.D0
  409.           DO 360 I=1,6
  410.             DEV(I)=(DEVT(I)-VAR0(8+I))/2.D0
  411.  360       CONTINUE
  412.         ENDIF
  413.         AJ2=PROCON(DEV,DEV,6)
  414.         IF (AJ2.LE.(1.D-10*YOUN0)) AJ2=0.D0
  415.         AJ2=(3.D0/2.D0*AJ2)**(.5D0)
  416. *
  417. *---> Mise à jour des déformations élastiques et plastiques
  418. *
  419.         TRELT=-1.D0*(PRES1-PRES00)/XKA1
  420.         DO 41 I=1,6
  421.           A=0.D0
  422.           IF (I.LE.3) A=1.D0
  423.           DEFP(I)=(DEVT(I)-(RSIG0(I)+PRES00*A))
  424.      .            /(2.D0*XGA1)
  425.           DEFP(I)=DEPST(I)-DEFP(I)-TRELT*A/3.D0
  426.  41     CONTINUE
  427. *
  428. *   Fin du cas sur XN0
  429. *-------------------------------------------------------------
  430. *
  431.       ENDIF
  432. *
  433. *---> Surface de charge de gonflement
  434. *
  435.       PHIT=AJ2*AJ2+XM0*XM0*(PRES1-APT+AIS)*(PRES1-APT-AIS)
  436.       PHI0=ABS(PHIT)
  437. *
  438. *
  439. *---> Test de convergence ou itération suivante
  440. *
  441.       iter0=iter0+1
  442.       IF (iter0.LT.200) THEN
  443.          GOTO 99
  444.       ELSE
  445.          PHI0=0.D0
  446.          KERRE=460
  447.          GOTO 99
  448.       ENDIF
  449. *
  450. *
  451.       ENDIF
  452. *
  453.       RETURN
  454.       END
  455.  
  456.  
  457.  
  458.  

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