Télécharger presto.eso

Retour à la liste

Numérotation des lignes :

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

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