Télécharger la plaquette Cast3M


fr en

Télécharger alon2.eso

Retour à la liste

Numérotation des lignes : 

alon2

  1. C ALON2     SOURCE    CHAT      05/01/12    21:20:54     5004
  2.       SUBROUTINE ALON2(Pat0,XMAT,DEPM,RSIG0,VAR0,
  3.      .           RSIGF,VARF,RDEFP,KERRE,IB,IGAU)
  4.  
  5. *
  6. *============================================================
  7. *
  8. *  Modèle d'argile partiellement saturé d'ALONSO
  9. *
  10. *============================================================
  11. *
  12. *  ENTREES
  13. *-----------------------------
  14. *
  15. * Pat0  : pression atmosphérique
  16. * XMAT  : champ des caractéristiques materiau
  17. * DEPM  : incrément des déformations plastiques
  18. * RSIG0 : contraintes initiales
  19. * VAR0  : variables internes initiales
  20. * IB    : numéro de l'élément
  21. * IGAU  : numéro du point de Gauss
  22. *
  23. * SORTIES
  24. *-----------------------------
  25. *
  26. * RSIGF : contraintes finales
  27. * VARF  : variables internes finales
  28. * RDEFP : increment de déformation plastique au cours du pas
  29. * KERRE : message d'erreur
  30. *
  31. *============================================================
  32. *       ICI IL FAUT PROGRAMMER EN FORTRAN PUR
  33. *=============================================================
  34. *
  35.       IMPLICIT INTEGER(I-N)
  36.       IMPLICIT REAL*8(A-H,O-Z)
  37. *
  38.       DIMENSION VAR0(*),VARF(*),XMAT(*)
  39.       DIMENSION RSIG0(*),RSIGF(*),RDEFP(*)
  40.       DIMENSION DEPM(*)
  41.       DIMENSION DEPT(6),DEVT(6),SIGT(6),DEFT(6),DEPM0(6)
  42. *
  43. *
  44. *   Données du materiau
  45. *===========================================================
  46. *
  47.       YOUN0=XMAT(1)
  48.       XNU0=XMAT(2)
  49.       XKS0=XMAT(5)
  50.       XLAM0=XMAT(6)
  51.       XM0=XMAT(7)
  52.       XKK0=XMAT(8)
  53.       PC00=XMAT(9)
  54.       P0=XMAT(10)
  55.       XLAM1=XMAT(11)
  56.       T0=XMAT(12)
  57.       TAU0=XMAT(13)
  58.       XG0=XMAT(14)
  59.       XK0=XMAT(15)
  60.       XE0=XMAT(16)
  61. *
  62. *   Quelques initialisations
  63. *=============================================================
  64. *
  65. *---> Terme GAMA0 intervenant dans l'écoulement non associatif
  66. *
  67.       GAMA0=XM0*(XM0-9.D0)*(XM0-3.D0)*XLAM1
  68.       GAMA0=GAMA0/(9.D0*(6.D0-XM0)*(XLAM1-XK0))
  69. *
  70. *---> Initialisations des sous incrémentations
  71. *
  72.       recal0=-100.D0
  73.       iter2=0
  74.       nmax0=1
  75.  93   DO 100 I=1,6
  76.         DEPM0(I)=DEPM(I)/nmax0
  77.         SIGT(I)=RSIG0(I)
  78.         DEFT(I)=0.D0
  79.  100  CONTINUE
  80.       PRES00=-1.D0*(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0
  81. *
  82. *---> Variables internes initiales
  83. *     . pression de consolidation à l'état saturé
  84. *     . limite élastique en succion
  85. *     . succions initiale et finale
  86. *  La succion finale est déstinée à etre calculée au préalable
  87. *  dans un opérateur particulier indépendant de ECOULE
  88. *
  89.       VOLU0=1.D0+XE0
  90.       PC0=VAR0(2)
  91.       IF (PC0.LE.1.D-30) PC0=PC00
  92.       SLIM0=VAR0(3)
  93.       SUCI0=VAR0(4)
  94. *
  95. *     MODIF POUR CAS SATURE
  96. *
  97.       SUCAUX=VARF(4)
  98.       IF(SUCI0.LT.0.D0) THEN
  99.         SUCI0=0.D0
  100.       ENDIF
  101.       IF(SUCAUX.LT.0.D0) THEN
  102.         SUCAUX=0.D0
  103.       ENDIF
  104. *
  105.       PC1=PC0
  106.       SLIM1=SLIM0
  107. *
  108. *---> Succion a la fin de la sous incrementation
  109. *
  110.  95   iter2=iter2+1
  111.       SUCF0=iter2*(SUCAUX-SUCI0)/nmax0+SUCI0
  112. *
  113. *     MODIF POUR CAS SATURE
  114. *
  115.       IF(SUCF0.LT.0.D0) THEN
  116.          SUCF0=0.D0
  117.       ENDIF
  118. *
  119. *---> Pente élasto-plastique (courbe de consolidation)
  120. *
  121.       XLAM2=XLAM1*((1.D0-T0)*EXP(-1.D0*TAU0*SUCF0)+T0)
  122. *
  123. *---> Pression de saturation
  124. *
  125.       PS0=XKK0*SUCF0
  126. *
  127. *
  128. *  Calcul élastique test
  129. *================================================================
  130. *
  131. *---> Déformations élastiques test et déformation plastique
  132. *
  133.          DO 11 I=1,6
  134.            A=0.D0
  135.            IF (I.LE.3) A=1.D0
  136.            DEPT(I)=DEPM0(I)
  137.  11      CONTINUE
  138. *
  139. *---> Contraintes
  140. *
  141.          PRES0=-1.D0*(SIGT(1)+SIGT(2)+SIGT(3))/3.D0
  142.          TREPS0=DEPT(1)+DEPT(2)+DEPT(3)
  143.          PRES1=PRES0*EXP(-1.D0*VOLU0/XK0*TREPS0)
  144.          DO 12 I=1,6
  145.            A=0.D0
  146.            IF (I.LE.3) A=1.D0
  147.            DEVT(I)=SIGT(I)+PRES0*A
  148.            DEVT(I)=DEVT(I)+2.D0*XG0*(DEPT(I)-TREPS0/3.D0*A)
  149.            SIGT(I)=DEVT(I)-PRES1*A
  150.  12      CONTINUE
  151. *
  152. *---> Contrainte équivalente
  153. *
  154.          STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0
  155.          IF (STEST2.LE.1.D-10) STEST2=0.D0
  156.          Stest=(STEST2)**(0.5D0)
  157. *
  158. *---> Pression de consolidation
  159. *      . à l'état saturé
  160. *      . à la succion finale
  161. *
  162.          PCS1=P0*((PC1/P0)**((XLAM1-XK0)/(XLAM2-XK0)))
  163. *
  164. *---> Fonction de charge
  165. *
  166.         PHIT=Stest*Stest+XM0*XM0*(PRES1+PS0)*(PRES1-PCS1)
  167. *
  168. *
  169. *
  170. *  Vérification du critère de plasticité
  171. *=========================================================
  172. *
  173. *---> Erreur admise
  174. *
  175.       ERR0=1.D-7*ABS(PHIT)
  176.       ERR0=MAX(ERR0,1.D-10*YOUN0)
  177. *
  178. *---> Critère de plasticité
  179. *
  180.       PHI0=PHIT
  181.       iter0=0
  182. *
  183.  99   IF (PHI0.LE.ERR0) THEN
  184. *
  185. *  On est élastique
  186. *=========================================================
  187. *
  188. *
  189.       DO 400 I=1,6
  190.          RSIGF(I)=SIGT(I)
  191.          RDEFP(I)=DEFT(I)
  192.          A=0.D0
  193.          IF (I.LE.3) A=1.D0
  194.          DEVT(I)=DEFT(I)-(DEFT(1)+DEFT(2)+DEFT(3))/3.D0*A
  195.  400  CONTINUE
  196.       EPSPT=PROCON(DEVT,DEVT,6)
  197.       IF (EPSPT.LE.1.D-20) EPSPT=0.D0
  198.       EPSPT=(2.D0/3.D0*EPSPT)**.5D0
  199.       VARF(1)=VAR0(1)+EPSPT
  200.       VARF(2)=PC1
  201.       VARF(3)=MAX(SLIM1,0.D0)
  202. *
  203.       ELSE
  204. *
  205. *  On plastifie
  206. *=========================================================
  207. *
  208. *---> Condition de consistance: calcul du paramètre plastique
  209. *
  210.       DF1=XM0*XM0*(2.D0*PRES1+PS0-PCS1)
  211.       DF2=XM0*XM0*(PRES1+PS0)*(XLAM1-XK0)/(XLAM2-XK0)
  212. ***** DF2=DF2*((PC1/P0)**((XLAM1-XLAM2)/(XLAM1-XK0)))
  213. ***** AM 10/12/03
  214.       DF2=DF2*((PC1/P0)**((XLAM1-XLAM2)/(XLAM2-XK0)))
  215.       denom0=12.D0*Stest*Stest*XG0*GAMA0
  216.       denom0=denom0+DF1*PRES1*VOLU0/XK0*DF1
  217.       denom0=denom0+DF2*PC1*VOLU0/(XLAM1-XK0)*DF1
  218. *
  219.       dlam0=PHIT/denom0
  220. *
  221. *
  222. *---> Critères de sous incrémentation
  223. *
  224.       IF (recal0.LT.0.D0) THEN
  225. *
  226.          xmax0=2000.D0
  227.          xmax1=0.D0
  228.          xmax2=0.D0
  229.          xnum0=100.D0
  230. *
  231.          cri0=ABS(VOLU0/XK0*DF1*dlam0)
  232.          IF (cri0.GT.1.D-2) THEN
  233.            xmax1=xnum0*cri0+10.D0
  234.          ENDIF
  235. *
  236.          cri1=ABS(VOLU0/(XK0-XLAM1)*DF1*dlam0)
  237.          IF (cri1.GT.1.D-2) THEN
  238.            xmax2=xnum0*cri1+10.D0
  239.          ENDIF
  240. *
  241.          xmax00=xmax2
  242.          IF (xmax1.GE.xmax2) xmax00=xmax1
  243.          IF (xmax00.GE.xmax0) xmax00=xmax0
  244.          nmax00=NINT(xmax00)
  245.          IF (nmax00.GT.nmax0) nmax0=nmax00
  246.          recal0=100.D0
  247. *
  248.          IF (nmax0.GT.1) THEN
  249.            iter2=0
  250.            GOTO 93
  251.          ENDIF
  252. *
  253.       ENDIF
  254. *
  255. *---> Mise à jour des contraintes
  256. *
  257.       PRES1=PRES1*EXP(-1.D0*VOLU0/XK0*DF1*dlam0)
  258.       DO 10 I=1,6
  259.         A=0.D0
  260.         IF (I.LE.3) A=1.D0
  261.         DEVT(I)=DEVT(I)-6.D0*XG0*DEVT(I)*GAMA0*dlam0
  262.         SIGT(I)=DEVT(I)-PRES1*A
  263.  10   CONTINUE
  264. *
  265. *---> Mise à jour de la pression de consolidation
  266. *      . dans le domaine saturé
  267. *      . à la succion finale
  268. *
  269.       PC1=PC1*EXP(VOLU0/(XLAM1-XK0)*DF1*dlam0)
  270.       PCS1=P0*((PC1/P0)**((XLAM1-XK0)/(XLAM2-XK0)))
  271. *
  272. *---> Déformation élastiques et plastiques
  273. *
  274.       TREPS0=-1.D0*XK0/VOLU0*LOG(PRES1/PRES00)
  275.       DO 20 I=1,6
  276.         A=0.D0
  277.         IF (I.LE.3) A=1.D0
  278.         DEPT(I)=(DEVT(I)-RSIG0(I)-PRES00*A)/(2.D0*XG0)
  279.         DEPT(I)=DEPT(I)+TREPS0/3.D0*A
  280.         DEFT(I)=iter2*DEPM(I)/nmax0-DEPT(I)
  281.  20   CONTINUE
  282.       TREPP0=DEFT(1)+DEFT(2)+DEFT(3)
  283. *
  284. *---> Mise à jour de la limite en succion
  285. *
  286.       SLIM1=(SLIM0+Pat0)*EXP(VOLU0/(XKS0-XLAM0)*TREPP0)-Pat0
  287. *
  288. *---> Mise à jour de la contrainte équivalente
  289. *
  290.       STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0
  291.       IF (STEST2.LE.1.D-10) STEST2=0.D0
  292.       Stest=(STEST2)**(0.5D0)
  293. *
  294. *---> Mise à jour de la fonction de charge
  295. *
  296.       PHIT=Stest*Stest+XM0*XM0*(PRES1+PS0)*(PRES1-PCS1)
  297.       PHI0=ABS(PHIT)
  298. *
  299. *---> Test de convergence ou itération suivante
  300. *
  301.       iter0=iter0+1
  302.       IF (iter0.LT.200) THEN
  303.          GOTO 99
  304.       ELSE
  305.          PHI0=0.D0
  306.          KERRE=460
  307.          GOTO 99
  308.       ENDIF
  309. *
  310.       ENDIF
  311. *
  312. *
  313. *   Vérification des sous incrémentations
  314. *
  315.       IF ((iter2.LT.nmax0).AND.(KERRE.EQ.0)) THEN
  316.          GOTO 95
  317.       ENDIF
  318. *
  319. *
  320. *
  321. *
  322.       RETURN
  323.       END
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  

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