Télécharger gurso2.eso

Retour à la liste

Numérotation des lignes :

  1. C GURSO2 SOURCE CHAT 05/01/13 00:22:19 5004
  2. SUBROUTINE GURSO2(DEPST,NSTRS,MFR1,IB,IGAU,
  3. 1 DSIGT,NCOMAT,SIG0,VAR0,XMAT,XCAR,NVARI,ICARA,
  4. 2 SIGF,VARF,DEFP,TRAC,KERRE,necou)
  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. * TRAC = COURBE DE TRACTION NON ENDOMMAGEE
  25. *
  26. * SORTIE :
  27. * --------
  28. *
  29. * SIGF(NSTRS) = CONTRAINTES FINALES
  30. * VARF(NVARI) = VARIALES INTERNES A LA FIN DU PAS D'INTEGRATION
  31. * DEFP(NSTRS) = INCREMENT DE DEFORMATION PLASTIQUE A LA FIN
  32. * DU PAS D'INTEGRATION
  33. * ============================================================
  34. * ICI IL FAUT PROGRAMMER EN FORTRAN PUR
  35. *=============================================================
  36. *
  37. IMPLICIT INTEGER(I-N)
  38. IMPLICIT REAL*8(A-H,O-Z)
  39. *
  40. -INC CCREEL
  41. -INC CCOPTIO
  42. *
  43. SEGMENT NECOU
  44. * COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
  45. INTEGER NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
  46. . ITYP,IFOURB,IFLUAG,
  47. . ICINE,ITHER,IFLUPL,ICYCL,IBI,
  48. . JFLUAG,KFLUAG,LFLUAG,
  49. . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
  50. ENDSEGMENT
  51. C COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
  52. C . ITYP,IFOURB,IFLUAG,
  53. C . ICINE,ITHER,IFLUPL,ICYCL,IBI,
  54. C . JFLUAG,KFLUAG,LFLUAG,
  55. C . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
  56. *
  57. DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*),
  58. & VARF(*),DEFP(*),DSIGT(*),TRAC(*)
  59. DIMENSION RDEPS0(6)
  60. DIMENSION RSIG0(6),RDEPS(6),RSIGF(6),RDEFP(6),RDIGT(6)
  61. DIMENSION DEVT(6),SIGT(6),DEFT(6),SIGT1(6)
  62. *
  63. * Adaptation de l'option de calcul vers le 3D massif de SIG0 a RSIG0
  64. *====================================================================
  65. *
  66. IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN
  67. *
  68. *---> 1 formulation massive
  69. *---> 2 formulation quasi incompressible
  70. *---> MASSIF 3D
  71. *
  72. IF (NSTRS .EQ. 6) THEN
  73. DO 110 I=1,NSTRS
  74. RSIG0(I)=SIG0(I)
  75. RDEPS0(I)=DEPST(I)
  76. 110 CONTINUE
  77. ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0)
  78. & .OR.(IFOUR .EQ. -1).OR.(IFOUR.EQ.-2))) THEN
  79. *
  80. *---> Calcul en mode deformations planes ou axisymetrique
  81. * ou contraintes planes
  82. *
  83. DO 115 I=1,NSTRS
  84. RSIG0(I)=SIG0(I)
  85. RDEPS0(I)=DEPST(I)
  86. 115 CONTINUE
  87. RSIG0(5)=0.D0
  88. RSIG0(6)=0.D0
  89. RDEPS0(5)=0.D0
  90. RDEPS0(6)=0.D0
  91. ENDIF
  92. ELSE
  93. KERRE = 99
  94. RETURN
  95. ENDIF
  96. *
  97. *
  98. * Passage des deformations de cisaillement exprimées
  99. * en GAMA aux déformations de cisaillement exprimées
  100. * en déformations
  101. *
  102. DO 116 I=1,6
  103. A=1.D0
  104. IF (I.GT.3) A=2.D0
  105. RDEPS0(I)=RDEPS0(I)/A
  106. RDEPS(I)=0.D0
  107. 116 CONTINUE
  108. *
  109. * Données du materiau
  110. *===========================================================
  111. *
  112. YOUNG=XMAT(1)
  113. XNU=XMAT(2)
  114. Q0=XMAT(6)
  115. FU0=XMAT(7)
  116. FF0=XMAT(8)
  117. FC0=XMAT(9)
  118. FNS0=XMAT(10)
  119. FNE0=XMAT(11)
  120. SNS0=XMAT(12)
  121. SNE0=XMAT(13)
  122. SIGN0=XMAT(14)
  123. EPSN0=XMAT(15)
  124. F0=XMAT(16)
  125. XSRMA=XMAT(17)
  126. Fmin0=1.D-8
  127. IF ((FF0.LE.Fmin0).OR.(FF0.GT.1.D0)) FF0=1.D0
  128. IF ((F0.LE.Fmin0).OR.(F0.GE.1.D0)) F0=0.D0
  129. C RMIN0=(1.D0-FF0)/(1.D0-F0)
  130. *
  131. * Calculs préliminaires
  132. *===================================
  133. *
  134. *---> Module de cisaillement G0 et de pression hydro XK0
  135. *
  136. G0=YOUNG/(2.D0*(1.D0+XNU))
  137. XK0=YOUNG/(3.D0*(1.D0-(2.D0*XNU)))
  138. *
  139. * Début des sous incrémentations
  140. *===================================
  141. *
  142. recal0=-100.D0
  143. iter2=0
  144. nmax0=1
  145. *
  146. *---> Variables internes test
  147. *
  148. * - déformation plastique équivalente
  149. * - fraction volumique de cavités
  150. *
  151. 93 EPSPT=VAR0(1)
  152. FT=VAR0(2)
  153. SIGMT=VAR0(3)
  154. EPSMT=VAR0(4)
  155. RHOT=VAR0(5)
  156. FNST=VAR0(6)
  157. FNET=VAR0(7)
  158. VARF(3)=VAR0(3)
  159. VARF(4)=VAR0(4)
  160. VARF(6)=VAR0(6)
  161. VARF(7)=VAR0(7)
  162. IF (RHOT.LE.1.D-30) THEN
  163. RHOT=1.D0
  164. VAR0(5)=1.D0
  165. ENDIF
  166. IF (FT.LE.Fmin0) FT=0.D0
  167. *
  168. * Contraintes planes
  169. * - on stocke les variables internes au début de
  170. * chaque sous incrémentation
  171. * - on stocke les déformations plastiques au début de
  172. * chaque sous incrémentation
  173. *
  174. IF (IFOUR.EQ.-2) THEN
  175. EPSPT1=EPSPT
  176. FT1=FT
  177. EPSMT1=EPSMT
  178. SIGMT1=SIGMT
  179. RHOT1=RHOT
  180. FNST1=FNST
  181. FNET1=FNET
  182. def1=0.D0
  183. def2=0.D0
  184. def3=0.D0
  185. rdeps3=0.D0
  186. ENDIF
  187. *
  188. * Contraintes planes
  189. * - on stocke les contraintes au début de
  190. * chaque sous incrémentation
  191. *
  192. IF (IFOUR.EQ.-2) THEN
  193. DO 171 I=1,6
  194. SIGT1(I)=RSIG0(I)
  195. 171 CONTINUE
  196. ENDIF
  197. DO 172 I=1,6
  198. SIGT(I)=RSIG0(I)
  199. 172 CONTINUE
  200. SMT0=(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0
  201. SMT=SMT0
  202. IF (IFOUR.EQ.-2) SMT1=SMT0
  203. *
  204. * On divise le chargement pour les sous incrementations
  205. *
  206. 95 DO 118 I=1,6
  207. RDEPS(I)=RDEPS0(I)/nmax0
  208. 118 CONTINUE
  209. iter2=iter2+1
  210. *
  211. * Initialisations
  212. *===================================
  213. *
  214. iter1=0
  215. *
  216. *---> paramètre de recalcul éventuel de la nucléation
  217. * il est ré-initialisé à chaque sous itération
  218. *
  219. recal1=-100.D0
  220. *
  221. *---> Calcul de la trace de l'incrément de déformation
  222. *
  223. IF (IFOUR.EQ.-2) THEN
  224. XLAM0=(1.D0-FT)*XK0-(2.D0*G0/3.D0)
  225. RDEPS(3)=-1.D0*XLAM0/(2.D0*G0+XLAM0)*(RDEPS(1)+RDEPS(2))
  226. ENDIF
  227. 98 treps0=RDEPS(1)+RDEPS(2)+RDEPS(3)
  228. *
  229. *
  230. * Calculs des grandeurs élastiques test
  231. *====================================================
  232. *
  233. *---> Variables internes test: cas contraintes planes
  234. *
  235. * - déformation plastique équivalente
  236. * - fraction volumique de cavités
  237. *
  238. IF (IFOUR.EQ.-2) THEN
  239. EPSPT=EPSPT1
  240. FT=FT1
  241. EPSMT=EPSMT1
  242. SIGMT=SIGMT1
  243. RHOT=RHOT1
  244. FNST=FNST1
  245. FNET=FNET1
  246. SMT=SMT1
  247. ENDIF
  248. *
  249. *---> Calcul des contraintes test élastiques
  250. * et de l'incrément des déformations plastiques test au
  251. * cours du pas
  252. *
  253. DO 101 I=1,6
  254. A=0.D0
  255. IF (I.LE.3) A=1.D0
  256. IF (IFOUR.EQ.-2) SIGT(I)=SIGT1(I)
  257. DEVT(I)=SIGT(I)-(A*SMT)
  258. DEVT(I)=DEVT(I)+(2.D0*G0*(RDEPS(I)-(A*treps0/3.D0)))
  259. DEFT(I)=0.D0
  260. 101 CONTINUE
  261. *
  262. *---> Contrainte équivalente STEST
  263. *
  264. STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0
  265. IF (STEST2.LE.1.D-10) STEST2=0.D0
  266. Stest=STEST2**0.5D0
  267. *
  268. *---> Limite d'élasticité lue sur la courbe d'écrouissage
  269. *
  270. * - Recherche des valeurs de déformations plastiques équivalentes
  271. * encadrant la valeur courante
  272. *
  273. EPSM0=TRAC(2*NCOURB)
  274. IF (EPSMT.GE.EPSM0) THEN
  275. Y1=TRAC(2*NCOURB-1)
  276. Y2=TRAC(2*NCOURB-1)
  277. EPS1=EPSM0
  278. EPS2=2.D0*EPSM0
  279. ELSE
  280. DO 300 I=1,(NCOURB-1)
  281. EPS1=TRAC(2*I)
  282. EPS2=TRAC(2*(I+1))
  283. IF ((EPSMT.GE.EPS1).AND.(EPSMT.LT.EPS2)) THEN
  284. Y1=TRAC(2*I-1)
  285. Y2=TRAC(2*(I+1)-1)
  286. GOTO 96
  287. ENDIF
  288. 300 CONTINUE
  289. ENDIF
  290. *
  291. * - Limite d'élasticité
  292. *
  293. 96 H1=(Y2-Y1)/(EPS2-EPS1)
  294. Ytest=(H1*(EPSMT-EPS2))+Y2
  295. *
  296. *---> Fraction de densité test
  297. *
  298. RHOT=1.D0/((1.D0/RHOT)+treps0)
  299. IF (RHOT.LT.1.D-10) RHOT=1.D-10
  300. ** IF ((FT.GT.Fmin0).AND.(RHOT.LT.RMIN0)) RHOT=RMIN0
  301. *
  302. *---> Coefficients de nucléation des cavités
  303. *
  304. dif0=VARF(3)-Ytest-SMT
  305. IF ((dif0.LE.1.D-5).AND.(treps0.GE.0.D0)) THEN
  306. BB0=EXP(-.5D0*((VARF(3)-SIGN0)/SNS0)**2.D0)
  307. BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0))
  308. ELSE
  309. BB0=0.D0
  310. ENDIF
  311. *
  312. *---> Cas de la nucléation pilotée en contrainte
  313. *
  314. IF ((FNS0.GT.Fmin0).AND.(recal1.LT.0.D0)) THEN
  315. *
  316. * Si on se trouve dans la "fenetre" où la
  317. * nucléation apparait, on sous incrémente le pas
  318. * de calcul afin de suivre les variations de BB0
  319. *
  320. dif0=Ytest+SMT
  321. IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND.
  322. . (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN
  323. xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0)
  324. nmax1=NINT(xmax1)*nmax0
  325. IF (nmax1.GT.2000) nmax1=2000
  326. IF (nmax1.GT.nmax0) THEN
  327. nmax0=nmax1
  328. iter2=0
  329. recal1=100.D0
  330. GOTO 93
  331. ENDIF
  332. ENDIF
  333. *
  334. ENDIF
  335. *
  336. *---> Fraction de cavité test
  337. * On résoud une équation du second degré
  338. *
  339. AA1=BB0*XK0/((1.D0-F0)*RHOT)
  340. BB1=BB0*XK0-1.D0
  341. CC1=BB0*SMT+1.D0-FT
  342. delta0=BB1*BB1+4.D0*AA1*CC1
  343. *
  344. IF (AA1.GT.1.D-5) THEN
  345. *
  346. * Le terme en contrainte de la nucléation n'est pas négligeable
  347. *
  348. XX1=(BB1+(delta0**.5D0))/(2.D0*AA1)
  349. FT=1.D0-XX1
  350. ELSE
  351. *
  352. * Le terme en contrainte de la nucléation est négligé
  353. *
  354. BB0=0.D0
  355. ENDIF
  356. FT=MAX(FT,0.D0)
  357. *
  358. *---> Contrainte totale et contrainte moyenne test
  359. *
  360. FNST=FNST-BB0*SMT
  361. SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT-
  362. . 1.D0)
  363. FNST=FNST+BB0*SMT
  364. DO 200 I=1,6
  365. A=0.D0
  366. IF (I.LE.3) A=1.D0
  367. SIGT(I)=DEVT(I)+(A*SMT)
  368. 200 CONTINUE
  369. *
  370. *---> Fonction d'endommagement
  371. *
  372. IF (FT.LE.FC0) THEN
  373. FST0=FT
  374. ddel0=1.D0
  375. ELSE
  376. FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0)
  377. ddel0=(FU0-FC0)/(FF0-FC0)
  378. ENDIF
  379. IF (FST0.GT.FU0) FST0=FU0
  380. *
  381. *---> Terme d'endommagement
  382. *
  383. SMTM=MAX(SMT,0.D0)
  384. SMTM=SMT
  385. Gtest=1.D0+(Q0*FST0)*(Q0*FST0)-
  386. . (2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest)))
  387. *
  388. *---> Fonction de charge test PHIT
  389. *
  390. PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest
  391. *
  392. *---> Cas du matériau complètement endommagé
  393. *
  394. 91 IF (FT.GE.FF0) THEN
  395. *
  396. * Les cavités
  397. *
  398. FT=FF0+Fmin0
  399. *
  400. * Les contraintes sont mises a zero
  401. * Les deformations sont plastiques
  402. *
  403. SMT=0.D0
  404. Stest=0.D0
  405. treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/XK0
  406. DO 12 I=1,6
  407. SIGT(I)=0.D0
  408. A=0.D0
  409. IF (I.LE.3) A=1.D0
  410. depel0=SIGT(I)-RSIG0(I)
  411. DEVT(I)=(depel0-(SMT-SMT0)*A)/(2.D0*G0)
  412. depel0=DEVT(I)+(treps1*A/3.D0)
  413. DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  414. IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  415. DEFT(3)=rdeps3+RDEPS(3)-depel0
  416. ENDIF
  417. 12 CONTINUE
  418. *
  419. * Deformation plastique equivalente
  420. *
  421. treps1=DEFT(1)+DEFT(2)+DEFT(3)
  422. DO 14 I=1,6
  423. A=0.D0
  424. IF (I.LE.3) A=1.D0
  425. DEVT(I)=DEFT(I)-treps1*A/3.D0
  426. 14 CONTINUE
  427. IF (treps1.LE.0.D0) THEN
  428. DO 16 I=1,3
  429. DEFT(I)=DEVT(I)
  430. 16 CONTINUE
  431. ENDIF
  432. dlam0=(2.D0*PROCON(DEVT,DEVT,6)/3.D0)**.5D0
  433. EPSPT=EPSPT+dlam0
  434. ** RHOT=RMIN0
  435. *
  436. Gtest=0.D0
  437. PHIT=0.D0
  438. *
  439. ENDIF
  440. *
  441. *
  442. * Vérification du critère de plasticité
  443. *=========================================================
  444. *
  445. *---> Erreur admise
  446. *
  447. ERR0=1.D-7*(ABS(PHIT))
  448. ERR0=MAX(ERR0,1.D-8*Ytest*Ytest)
  449. *
  450. *---> Critère de plasticité
  451. *
  452. PHI0=PHIT
  453. iter0=0
  454. *
  455. 99 IF (PHI0.LE.ERR0) THEN
  456. *
  457. * On est élastique
  458. *=========================================================
  459. *
  460. *--------------------------------------------------------------
  461. * Cas particulier des contraintes planes
  462. *--------------------------------------------------------------
  463. *
  464. IF (IFOUR.EQ.-2) THEN
  465. *
  466. *---> On vérifie qu'on est en contraintes planes
  467. *
  468. iter1=iter1+1
  469. SIG3=ABS(SIGT(3))
  470. SIG30=MAX(Stest*1.D-5,1.D-3)
  471. IF (SIG3.GT.SIG30) THEN
  472. XLAM0=(1.D0-FT)*XK0-(2.D0*G0/3.D0)
  473. RDEPS(3)=-1.D0*(XLAM0/(2.D0*G0+XLAM0))*
  474. . (RDEPS(1)-DEFT(1)+RDEPS(2)-DEFT(2)+def1+def2)
  475. RDEPS(3)=RDEPS(3)+DEFT(3)-def3
  476. IF (iter1.LE.200)THEN
  477. GOTO 98
  478. ELSE
  479. KERRE=460
  480. ENDIF
  481. ELSE
  482. EPSPT1=EPSPT
  483. FT1=FT
  484. SIGMT1=SIGMT
  485. EPSMT1=EPSMT
  486. RHOT1=RHOT
  487. FNST1=FNST
  488. FNET1=FNET
  489. SMT1=SMT
  490. SIGT1(1)=SIGT(1)
  491. SIGT1(2)=SIGT(2)
  492. SIGT1(3)=SIGT(3)
  493. SIGT1(4)=SIGT(4)
  494. SIGT1(5)=SIGT(5)
  495. SIGT1(6)=SIGT(6)
  496. rdeps3=rdeps3+RDEPS(3)
  497. def1=DEFT(1)
  498. def2=DEFT(2)
  499. def3=DEFT(3)
  500. ENDIF
  501. ENDIF
  502. *
  503. *-----------------------------------------------------
  504. * Fin du cas particulier des contraintes planes
  505. *-----------------------------------------------------
  506. *
  507. DO 400 I=1,6
  508. RSIGF(I)=SIGT(I)
  509. RDEFP(I)=DEFT(I)
  510. RDIGT(I)=RSIGF(I)-RSIG0(I)
  511. 400 CONTINUE
  512. VARF(1)=EPSPT
  513. VARF(2)=FT
  514. VARF(3)=MAX(VARF(3),Ytest+SMT)
  515. VARF(4)=MAX(VARF(4),EPSMT)
  516. VARF(5)=RHOT
  517. VARF(6)=FNST
  518. VARF(7)=FNET
  519. *
  520. ELSE
  521. *
  522. * On plastifie
  523. *=========================================================
  524. *
  525. *---> Valeur minimale de Stest au dessous de laquelle
  526. * les contraintes sont considérées comme étant
  527. * purement hydrostatiques
  528. *
  529. Gmin0=2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest))
  530. Gmin0=Gmin0-((Q0*FST0)*(Q0*FST0))
  531. Gmin0=(ABS(Gmin0))**.5D0
  532. Gmin0=Gmin0*Ytest*1.D-6
  533. *
  534. *---> Coefficient de nucléation
  535. *
  536. dif0=VARF(4)-EPSMT
  537. IF (dif0.LE.1.D-5) THEN
  538. DD0=EXP(-.5D0*((VARF(4)-EPSN0)/SNE0)**2.D0)
  539. DD0=DD0*FNE0/(SNE0*((2.D0*XPI)**.5D0))
  540. ELSE
  541. DD0=0.D0
  542. ENDIF
  543. *
  544. *---> Condition de consistance
  545. *
  546. IF (Stest.LT.Gmin0) THEN
  547. tr0=3.D0*Q0*FST0*Ytest*SINH(3.D0*SMTM/(2.D0*Ytest))
  548. DE0=SMT*tr0/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest)
  549. ELSE
  550. tr0=3.D0*Q0*FST0*Ytest*SINH(3.D0*SMTM/(2.D0*Ytest))
  551. tr0=tr0/(2.D0*Stest)
  552. DE0=(SMT*tr0+Stest)/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest)
  553. ENDIF
  554. *
  555. H0=H1*DE0
  556. DS1=XK0*(2.D0*(1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT-1.D0)
  557. DF0=((1.D0-FT+XSRMA*(FNET+FNST))*tr0+DD0*DE0+BB0*H0)/
  558. . (1.D0+BB0*DS1)
  559. *
  560. DCOS0=(DS1*DF0/(2.D0*Ytest))+(SMT*H0)/(2.D0*Ytest*Ytest)
  561. DCOS0=DCOS0*3.D0*SINH(3.D0*SMTM/(2.D0*Ytest))
  562. *
  563. DG0=2.D0*Q0*Q0*ddel0*DF0*FST0
  564. DG0=DG0+2.D0*Q0*FST0*DCOS0
  565. DG0=DG0-2.D0*Q0*ddel0*DF0*COSH(3.D0*SMTM/(2.D0*Ytest))
  566. *
  567. IF (Stest.LT.Gmin0) THEN
  568. denom0=Ytest*Ytest*DG0
  569. ELSE
  570. denom0=6.D0*G0*Stest+Ytest*Ytest*DG0
  571. denom0=denom0+2.D0*H0*Gtest*Ytest
  572. ENDIF
  573. *
  574. IF (denom0.GT.1.D-5) THEN
  575. dlam0=PHIT/denom0
  576. ELSE
  577. IF (recal0.LE.0.D0) THEN
  578. recal0=100.D0
  579. nmax0=2000
  580. iter2=0
  581. GOTO 93
  582. ENDIF
  583. ENDIF
  584. *
  585. *
  586. *
  587. * Vérification des hypothèses de calcul
  588. *===============================================================
  589. *
  590. IF (recal0.LE.0.D0) THEN
  591. *
  592. xnum0=10.D0
  593. *
  594. xlim0=1.D-1
  595. xmax1=0.D0
  596. xmax2=0.D0
  597. xmax3=0.D0
  598. xmax4=0.D0
  599. xmax5=0.D0
  600. xmax6=0.D0
  601. xmax7=0.D0
  602. *
  603. *---> Condition 3.D0*G*dlam0/Stest << 1
  604. *
  605. IF (Stest.GE.Gmin0) THEN
  606. cri01=ABS(3.D0*G0*dlam0/Stest)
  607. IF (cri01.GT.xlim0) THEN
  608. xmax1=xnum0*cri01+5.D0
  609. ENDIF
  610. ENDIF
  611. *
  612. *---> Condition H0*dlam0/Ytest << 1
  613. *
  614. cri01=ABS(H0*dlam0/Ytest)
  615. IF (cri01.GT.xlim0) THEN
  616. xmax2=xnum0*cri01+5.D0
  617. ENDIF
  618. *
  619. *---> Condition DCOS0*dlam0/COSH << 1
  620. *
  621. cri01=ABS(DCOS0*dlam0/COSH(3.D0*SMTM/(2.D0*Ytest)))
  622. IF (cri01.GT.xlim0) THEN
  623. xmax3=xnum0*cri01+5.D0
  624. ENDIF
  625. *
  626. *---> Condition DF0*dlam0/F << 1
  627. *
  628. Fmax0=MAX(FNS0,FNE0)
  629. Fmax0=MAX(Fmax0,F0)
  630. Fmax0=MAX(Fmax0,FT)
  631. Fmax0=MAX(Fmax0,Fmin0)
  632. cri01=ABS(DF0*dlam0/Fmax0)
  633. IF (cri01.GT.xlim0) THEN
  634. xmax4=xnum0*cri01+5.D0
  635. ENDIF
  636. *
  637. *---> Condition DG0*dlam0/G << 1
  638. *
  639. Gmax0=MAX(ABS(Gtest),1.D0)
  640. cri01=ABS(DG0*dlam0/Gmax0)
  641. IF (cri01.GT.xlim0) THEN
  642. xmax5=xnum0*cri01+5.D0
  643. ENDIF
  644. *
  645. *---> Condition DS1*dlam0/SMT << 1
  646. *
  647. SMTmax=MAX(ABS(SMT),(1.D-3*TRAC(1)))
  648. cri01=ABS(DS1*DF0*dlam0/SMTmax)
  649. IF (cri01.GT.xlim0) THEN
  650. xmax6=xnum0*cri01+5.D0
  651. ENDIF
  652. *
  653. *---> Condition tr0*dlam0/tr(EPSP) << 1
  654. *
  655. tremp0=(FT-F0)/(1.D0-F0)/RHOT
  656. IF (ABS(tremp0).GT.1.D-5) THEN
  657. cri01=ABS(tr0*dlam0/tremp0)
  658. IF (cri01.GT.xlim0) THEN
  659. xmax7=xnum0*cri01+5.D0
  660. ENDIF
  661. ENDIF
  662. *
  663. *
  664. IF (xmax1.GT.xmax2) THEN
  665. xmax00=xmax1
  666. ELSE
  667. xmax00=xmax2
  668. ENDIF
  669. IF (xmax3.GT.xmax00) xmax00=xmax3
  670. IF (xmax4.GT.xmax00) xmax00=xmax4
  671. IF (xmax5.GT.xmax00) xmax00=xmax5
  672. IF (xmax6.GT.xmax00) xmax00=xmax6
  673. IF (xmax7.GT.xmax00) xmax00=xmax7
  674. *
  675. nmax00=NINT(xmax00)
  676. recal0=100.D0
  677. IF (nmax00.GT.nmax0) nmax0=nmax00
  678. IF (nmax0.GT.2000) nmax0=2000
  679. *
  680. *---> On recommence le calcul en sous incrémentant le
  681. * pas de temps
  682. *
  683. IF (nmax0.GE.2) THEN
  684. iter2=0
  685. GOTO 93
  686. ENDIF
  687. *
  688. ENDIF
  689. *
  690. * Cas particulier de la nucléation
  691. *============================================================
  692. *
  693. IF (recal1.LT.0.D0) THEN
  694. *
  695. *---> Condition sur la nucléation pilotée en contrainte
  696. *
  697. IF (FNS0.GT.Fmin0) THEN
  698. *
  699. * Si on se trouve dans la fenetre, on sous incremente
  700. *
  701. dif0=Ytest+H0*dlam0+SMT-DS1*DF0*dlam0
  702. IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND.
  703. . (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN
  704. xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0)
  705. nmax1=NINT(xmax1)*nmax0
  706. IF (nmax1.GT.2000) nmax1=2000
  707. IF (nmax1.GT.nmax0) THEN
  708. nmax0=nmax1
  709. iter2=0
  710. recal1=100.D0
  711. GOTO 93
  712. ENDIF
  713. ENDIF
  714. *
  715. ENDIF
  716. *
  717. *---> Condition sur la nucléation pilotée en deformation
  718. *
  719. IF (FNE0.GT.Fmin0) THEN
  720. *
  721. * Si on se trouve dans la fenetre, on sous incremente
  722. *
  723. dif0=EPSMT+DE0*dlam0
  724. IF ((dif0.GE.(EPSN0-SNE0*2.D0)).AND.
  725. . (VARF(4).LE.(EPSN0+SNE0*2.D0))) THEN
  726. xmax1=200.D0*(dif0-VARF(4))/(4.D0*SNE0)
  727. nmax1=NINT(xmax1)*nmax0
  728. IF (nmax1.GT.2000) nmax1=2000
  729. IF (nmax1.GT.nmax0) THEN
  730. nmax0=nmax1
  731. iter2=0
  732. recal1=100.D0
  733. GOTO 93
  734. ENDIF
  735. ENDIF
  736. *
  737. ENDIF
  738. *
  739. ENDIF
  740. *
  741. *
  742. * Mise à jour des grandeurs après écoulement plastique
  743. *============================================================
  744. *
  745. *---> Variables internes mises à jour
  746. *
  747. IF (Stest.GT.Gmin0) EPSPT=EPSPT+dlam0
  748. EPSMT=EPSMT+DE0*dlam0
  749. FNST=FNST+BB0*(H0-DS1*DF0)*dlam0
  750. FNET=FNET+DD0*DE0*dlam0
  751. FT=FT+DF0*dlam0
  752. IF (FT.LT.Fmin0) FT=0.D0
  753. *
  754. *---> Contraintes mises à jour
  755. *
  756. SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/RHOT/(1.D0-F0)-
  757. . 1.D0)
  758. DO 601 I=1,6
  759. A=0.D0
  760. IF (I.LE.3) A=1.D0
  761. IF (Stest.LE.Gmin0) THEN
  762. DEVT(I)=0.D0
  763. ELSE
  764. DEVT(I)=DEVT(I)-3.D0*G0*DEVT(I)*dlam0/Stest
  765. ENDIF
  766. SIGT(I)=DEVT(I)+(SMT*A)
  767. 601 CONTINUE
  768. *
  769. *---> Incrément de déformations plastiques mise à jour
  770. *
  771. treps1=(SMT/(1.D0-FT)-SMT0/(1.D0-VAR0(2)))/XK0
  772. DO 610 I=1,6
  773. A=0.D0
  774. IF (I.LE.3) A=1.D0
  775. depel0=(DEVT(I)-RSIG0(I)+(SMT0*A))/(2.D0*G0)
  776. depel0=depel0+(treps1*A/3.D0)
  777. DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  778. IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  779. DEFT(3)=rdeps3+RDEPS(3)-depel0
  780. ENDIF
  781. 610 CONTINUE
  782. *
  783. *---> Contrainte équivalente Stest mise à jour
  784. *
  785. STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0
  786. IF (STEST2.LE.1.D-10) STEST2=0.D0
  787. Stest= STEST2 ** 0.5D0
  788. *
  789. *---> Limite d'élasticité lue sur la courbe d'écrouissage
  790. * mise à jour
  791. *
  792. * - Recherche des valeurs de déformations plastiques équivalentes
  793. * encadrant la valeur courante
  794. *
  795. EPSM0=TRAC(2*NCOURB)
  796. IF (EPSMT.GE.EPSM0) THEN
  797. Y1=TRAC(2*NCOURB-1)
  798. Y2=TRAC(2*NCOURB-1)
  799. EPS1=EPSM0
  800. EPS2=2.D0*EPSM0
  801. ELSE
  802. DO 700 I=1,(NCOURB-1)
  803. EPS1=TRAC(2*I)
  804. EPS2=TRAC(2*(I+1))
  805. IF ((EPSMT.GE.EPS1).AND.(EPSMT.LT.EPS2)) THEN
  806. Y1=TRAC(2*I-1)
  807. Y2=TRAC(2*(I+1)-1)
  808. GOTO 97
  809. ENDIF
  810. 700 CONTINUE
  811. ENDIF
  812. *
  813. * - Limite d'élasticité mise à jour
  814. *
  815. 97 H1=(Y2-Y1)/(EPS2-EPS1)
  816. Ytest=(H1*(EPSMT-EPS2))+Y2
  817. *
  818. *---> Mise à jour des coefficients de nucléation
  819. *
  820. IF ((Ytest+SMT).GE.VARF(3)) THEN
  821. BB0=EXP(-.5D0*((Ytest+SMT-SIGN0)/SNS0)**2.D0)
  822. BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0))
  823. ELSE
  824. BB0=0.D0
  825. ENDIF
  826. IF (EPSMT.GE.VARF(4)) THEN
  827. DD0=EXP(-.5D0*((EPSMT-EPSN0)/SNE0)**2.D0)
  828. DD0=DD0*FNE0/(SNE0*((2.D0*XPI)**.5D0))
  829. ELSE
  830. DD0=0.D0
  831. ENDIF
  832. *
  833. *---> Fonction d'endommagement
  834. *
  835. IF (FT.LE.FC0) THEN
  836. FST0=FT
  837. ddel0=1.D0
  838. ELSE
  839. FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0)
  840. ddel0=(FU0-FC0)/(FF0-FC0)
  841. ENDIF
  842. IF (FST0.GT.FU0) FST0=FU0
  843. *
  844. *---> Terme d'endommagement
  845. *
  846. SMTM=MAX(SMT,0.D0)
  847. SMTM=SMT
  848. Gtest=1.D0+(Q0*FST0)*(Q0*FST0)-
  849. . (2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest)))
  850. *
  851. *---> Fonction de charge test PHIT mise à jour
  852. *
  853. PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest
  854. PHI0=ABS(PHIT)
  855. *
  856. *---> Cas particulier des matériaux totalement endommagés
  857. *
  858. 92 IF (FT.GE.FF0) THEN
  859. *
  860. FT=FF0+Fmin0
  861. *
  862. SMT=0.D0
  863. treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/XK0
  864. DO 13 I=1,6
  865. SIGT(I)=0.D0
  866. A=0.D0
  867. IF (I.LE.3) A=1.D0
  868. depel0=SIGT(I)-RSIG0(I)
  869. DEVT(I)=(depel0-((SMT-SMT0)*A))/(2.D0*G0)
  870. depel0=DEVT(I)+(treps1*A/3.D0)
  871. DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  872. IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  873. DEFT(3)=rdeps3+RDEPS(3)-depel0
  874. ENDIF
  875. 13 CONTINUE
  876. treps1=DEFT(1)+DEFT(2)+DEFT(3)
  877. DO 15 I=1,6
  878. A=0.D0
  879. IF (I.LE.1) A=1.D0
  880. DEVT(I)=DEFT(I)-(DEFT(1)+DEFT(2)+DEFT(3))*A/3.D0
  881. 15 CONTINUE
  882. IF (treps1.LE.0.D0) THEN
  883. DO 17 I=1,3
  884. DEFT(I)=DEVT(I)
  885. 17 CONTINUE
  886. ENDIF
  887. dlam0=(2.D0*PROCON(DEVT,DEVT,6)/3.D0)**.5D0
  888. EPSPT=EPSPT+dlam0
  889. ** RHOT=RMIN0
  890. *
  891. Gtest=0.D0
  892. PHIT=0.D0
  893. PHI0=0.D0
  894. *
  895. ENDIF
  896. *
  897. *---> Test de convergence ou itération suivante
  898. *
  899. iter0=iter0+1
  900. IF (iter0.LT.200) THEN
  901. GOTO 99
  902. ELSE
  903. PHI0=0.D0
  904. KERRE=460
  905. GOTO 99
  906. ENDIF
  907. *
  908. ENDIF
  909. *
  910. *
  911. * Vérification des sous incrémentations
  912. *
  913. IF ((iter2.LT.nmax0).AND.(KERRE.EQ.0)) THEN
  914. GOTO 95
  915. ENDIF
  916. *
  917. * Passage des deformations de cisaillement exprimées
  918. * en déformations aux déformations de cisaillement exprimées
  919. * en GAMA
  920. *
  921. DO 117 I=1,6
  922. A=1.D0
  923. IF (I.GT.3) A=2.D0
  924. RDEFP(I)=RDEFP(I)*A
  925. 117 CONTINUE
  926. *
  927. *
  928. * Passage a l'option de calcul pour les contraintes
  929. *=========================================================
  930. *
  931. IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN
  932. IF (NSTRS .EQ. 6) THEN
  933. *
  934. *---> MASSIF 3D
  935. *
  936. DO 170 I=1,NSTRS
  937. SIGF(I)=RSIGF(I)
  938. DEFP(I)=RDEFP(I)
  939. DSIGT(I)=RDIGT(I)
  940. 170 CONTINUE
  941. ELSE IF ( NSTRS .EQ. 4 ) THEN
  942. *
  943. *---> Calcul axisymétrique ou contraintes planes
  944. *
  945. DO 180 I=1,NSTRS
  946. SIGF(I)=RSIGF(I)
  947. DEFP(I)=RDEFP(I)
  948. DSIGT(I)=RDIGT(I)
  949. 180 CONTINUE
  950. ENDIF
  951. ENDIF
  952. RETURN
  953. *
  954. END
  955.  
  956.  
  957.  
  958.  
  959.  
  960.  
  961.  
  962.  
  963.  
  964.  
  965.  
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  

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