Télécharger alon2.eso

Retour à la liste

Numérotation des lignes :

  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