Télécharger pakgon.eso

Retour à la liste

Numérotation des lignes :

  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