Télécharger gurso2.eso

Retour à la liste

Numérotation des lignes :

  1. C GURSO2 SOURCE FD218221 17/09/07 21:15:02 9528
  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. * Correction SEMI: Mars 2017: le rapport de densite est exprime en deformation logarithmique
  299. *=========================================
  300. * Dans la note 96-566, le rapport de densite est deduit de l'expression lineaire de la variation de volume:
  301. * trace(epsilon)=dV/V0
  302. * avec: V le volume total, V0 le volume total initial et epsilon la deformation totale
  303. * On en deduit: V/V0=1+trace(epsilon)
  304. * Puis le rapport de densite RHO=V0/V=1/(1+trace(epsilon))
  305. * Cette expression n'est pas imposée dans l'article de Needleman et Tvergaard de réference
  306. * "AN ANALYSIS OF DUCTILE RUPTURE MODES AT A CRACK TIP"
  307. * A. NEEDLEMAN, V. TVERGAARD,
  308. * J. Mech. Phys. Solids Vol. 35, No. 2, pp. 151-183, 1987
  309. * Si on exprime la variation de volume par son expression logarithmique:
  310. * trace(depsilon)=dln(V)
  311. * avec depsilon la variation de la deformation totale et dln(V) la variation du logarithme du volume V
  312. * On en deduit: ln(V/V0)=trace(epsilon)
  313. * Puis le rapport de densite RHO=exp(-trace(epsilon))
  314. * Cette expression donne de meilleurs résultats pour les très forts taux de déformation.
  315. *==========================================
  316. * RHOT=1.D0/((1.D0/RHOT)+treps0)
  317. treps00 = LOG((1.D0/RHOT))
  318. treps00 = treps00 + treps0
  319. RHOT = EXP(-1. * treps00)
  320. * Fin Correction SEMI: Mars 2017
  321. IF (RHOT.LT.1.D-10) RHOT=1.D-10
  322. ** IF ((FT.GT.Fmin0).AND.(RHOT.LT.RMIN0)) RHOT=RMIN0
  323. *
  324. *---> Coefficients de nucléation des cavités
  325. *
  326. dif0=VARF(3)-Ytest-SMT
  327. IF ((dif0.LE.1.D-5).AND.(treps0.GE.0.D0)) THEN
  328. BB0=EXP(-.5D0*((VARF(3)-SIGN0)/SNS0)**2.D0)
  329. BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0))
  330. ELSE
  331. BB0=0.D0
  332. ENDIF
  333. *
  334. *---> Cas de la nucléation pilotée en contrainte
  335. *
  336. IF ((FNS0.GT.Fmin0).AND.(recal1.LT.0.D0)) THEN
  337. *
  338. * Si on se trouve dans la "fenetre" où la
  339. * nucléation apparait, on sous incrémente le pas
  340. * de calcul afin de suivre les variations de BB0
  341. *
  342. dif0=Ytest+SMT
  343. IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND.
  344. . (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN
  345. xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0)
  346. nmax1=NINT(xmax1)*nmax0
  347. IF (nmax1.GT.2000) nmax1=2000
  348. IF (nmax1.GT.nmax0) THEN
  349. nmax0=nmax1
  350. iter2=0
  351. recal1=100.D0
  352. GOTO 93
  353. ENDIF
  354. ENDIF
  355. *
  356. ENDIF
  357. *
  358. *---> Fraction de cavité test
  359. * On résoud une équation du second degré
  360. *
  361. AA1=BB0*XK0/((1.D0-F0)*RHOT)
  362. BB1=BB0*XK0-1.D0
  363. CC1=BB0*SMT+1.D0-FT
  364. delta0=BB1*BB1+4.D0*AA1*CC1
  365. *
  366. IF (AA1.GT.1.D-5) THEN
  367. *
  368. * Le terme en contrainte de la nucléation n'est pas négligeable
  369. *
  370. XX1=(BB1+(delta0**.5D0))/(2.D0*AA1)
  371. FT=1.D0-XX1
  372. ELSE
  373. *
  374. * Le terme en contrainte de la nucléation est négligé
  375. *
  376. BB0=0.D0
  377. ENDIF
  378. FT=MAX(FT,0.D0)
  379. *
  380. *---> Contrainte totale et contrainte moyenne test
  381. *
  382. FNST=FNST-BB0*SMT
  383. SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT-
  384. . 1.D0)
  385. FNST=FNST+BB0*SMT
  386. DO 200 I=1,6
  387. A=0.D0
  388. IF (I.LE.3) A=1.D0
  389. SIGT(I)=DEVT(I)+(A*SMT)
  390. 200 CONTINUE
  391. *
  392. *---> Fonction d'endommagement
  393. *
  394. IF (FT.LE.FC0) THEN
  395. FST0=FT
  396. ddel0=1.D0
  397. ELSE
  398. FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0)
  399. ddel0=(FU0-FC0)/(FF0-FC0)
  400. ENDIF
  401. IF (FST0.GT.FU0) FST0=FU0
  402. *
  403. *---> Terme d'endommagement
  404. *
  405. SMTM=MAX(SMT,0.D0)
  406. SMTM=SMT
  407. Gtest=1.D0+(Q0*FST0)*(Q0*FST0)-
  408. . (2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest)))
  409. *
  410. *---> Fonction de charge test PHIT
  411. *
  412. PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest
  413. *
  414. *---> Cas du matériau complètement endommagé
  415. *
  416. 91 IF (FT.GE.FF0) THEN
  417. *
  418. * Les cavités
  419. *
  420. FT=FF0+Fmin0
  421. *
  422. * Les contraintes sont mises a zero
  423. * Les deformations sont plastiques
  424. *
  425. SMT=0.D0
  426. Stest=0.D0
  427. treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/XK0
  428. DO 12 I=1,6
  429. SIGT(I)=0.D0
  430. A=0.D0
  431. IF (I.LE.3) A=1.D0
  432. depel0=SIGT(I)-RSIG0(I)
  433. DEVT(I)=(depel0-(SMT-SMT0)*A)/(2.D0*G0)
  434. depel0=DEVT(I)+(treps1*A/3.D0)
  435. DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  436. IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  437. DEFT(3)=rdeps3+RDEPS(3)-depel0
  438. ENDIF
  439. 12 CONTINUE
  440. *
  441. * Deformation plastique equivalente
  442. *
  443. treps1=DEFT(1)+DEFT(2)+DEFT(3)
  444. DO 14 I=1,6
  445. A=0.D0
  446. IF (I.LE.3) A=1.D0
  447. DEVT(I)=DEFT(I)-treps1*A/3.D0
  448. 14 CONTINUE
  449. IF (treps1.LE.0.D0) THEN
  450. DO 16 I=1,3
  451. DEFT(I)=DEVT(I)
  452. 16 CONTINUE
  453. ENDIF
  454. dlam0=(2.D0*PROCON(DEVT,DEVT,6)/3.D0)**.5D0
  455. EPSPT=EPSPT+dlam0
  456. ** RHOT=RMIN0
  457. *
  458. Gtest=0.D0
  459. PHIT=0.D0
  460. *
  461. ENDIF
  462. *
  463. *
  464. * Vérification du critère de plasticité
  465. *=========================================================
  466. *
  467. *---> Erreur admise
  468. *
  469. ERR0=1.D-7*(ABS(PHIT))
  470. ERR0=MAX(ERR0,1.D-8*Ytest*Ytest)
  471. *
  472. *---> Critère de plasticité
  473. *
  474. PHI0=PHIT
  475. iter0=0
  476. *
  477. 99 IF (PHI0.LE.ERR0) THEN
  478. *
  479. * On est élastique
  480. *=========================================================
  481. *
  482. *--------------------------------------------------------------
  483. * Cas particulier des contraintes planes
  484. *--------------------------------------------------------------
  485. *
  486. IF (IFOUR.EQ.-2) THEN
  487. *
  488. *---> On vérifie qu'on est en contraintes planes
  489. *
  490. iter1=iter1+1
  491. SIG3=ABS(SIGT(3))
  492. SIG30=MAX(Stest*1.D-5,1.D-3)
  493. IF (SIG3.GT.SIG30) THEN
  494. XLAM0=(1.D0-FT)*XK0-(2.D0*G0/3.D0)
  495. RDEPS(3)=-1.D0*(XLAM0/(2.D0*G0+XLAM0))*
  496. . (RDEPS(1)-DEFT(1)+RDEPS(2)-DEFT(2)+def1+def2)
  497. RDEPS(3)=RDEPS(3)+DEFT(3)-def3
  498. IF (iter1.LE.200)THEN
  499. GOTO 98
  500. ELSE
  501. KERRE=460
  502. ENDIF
  503. ELSE
  504. EPSPT1=EPSPT
  505. FT1=FT
  506. SIGMT1=SIGMT
  507. EPSMT1=EPSMT
  508. RHOT1=RHOT
  509. FNST1=FNST
  510. FNET1=FNET
  511. SMT1=SMT
  512. SIGT1(1)=SIGT(1)
  513. SIGT1(2)=SIGT(2)
  514. SIGT1(3)=SIGT(3)
  515. SIGT1(4)=SIGT(4)
  516. SIGT1(5)=SIGT(5)
  517. SIGT1(6)=SIGT(6)
  518. rdeps3=rdeps3+RDEPS(3)
  519. def1=DEFT(1)
  520. def2=DEFT(2)
  521. def3=DEFT(3)
  522. ENDIF
  523. ENDIF
  524. *
  525. *-----------------------------------------------------
  526. * Fin du cas particulier des contraintes planes
  527. *-----------------------------------------------------
  528. *
  529. DO 400 I=1,6
  530. RSIGF(I)=SIGT(I)
  531. RDEFP(I)=DEFT(I)
  532. RDIGT(I)=RSIGF(I)-RSIG0(I)
  533. 400 CONTINUE
  534. VARF(1)=EPSPT
  535. VARF(2)=FT
  536. VARF(3)=MAX(VARF(3),Ytest+SMT)
  537. VARF(4)=MAX(VARF(4),EPSMT)
  538. VARF(5)=RHOT
  539. VARF(6)=FNST
  540. VARF(7)=FNET
  541. *
  542. ELSE
  543. *
  544. * On plastifie
  545. *=========================================================
  546. *
  547. *---> Valeur minimale de Stest au dessous de laquelle
  548. * les contraintes sont considérées comme étant
  549. * purement hydrostatiques
  550. *
  551. Gmin0=2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest))
  552. Gmin0=Gmin0-((Q0*FST0)*(Q0*FST0))
  553. Gmin0=(ABS(Gmin0))**.5D0
  554. Gmin0=Gmin0*Ytest*1.D-6
  555. *
  556. *---> Coefficient de nucléation
  557. *
  558. dif0=VARF(4)-EPSMT
  559. IF (dif0.LE.1.D-5) THEN
  560. DD0=EXP(-.5D0*((VARF(4)-EPSN0)/SNE0)**2.D0)
  561. DD0=DD0*FNE0/(SNE0*((2.D0*XPI)**.5D0))
  562. ELSE
  563. DD0=0.D0
  564. ENDIF
  565. *
  566. *---> Condition de consistance
  567. *
  568. IF (Stest.LT.Gmin0) THEN
  569. tr0=3.D0*Q0*FST0*Ytest*SINH(3.D0*SMTM/(2.D0*Ytest))
  570. DE0=SMT*tr0/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest)
  571. ELSE
  572. tr0=3.D0*Q0*FST0*Ytest*SINH(3.D0*SMTM/(2.D0*Ytest))
  573. tr0=tr0/(2.D0*Stest)
  574. DE0=(SMT*tr0+Stest)/((1.D0-FT+XSRMA*(FNET+FNST))*Ytest)
  575. ENDIF
  576. *
  577. H0=H1*DE0
  578. DS1=XK0*(2.D0*(1.D0-FT+XSRMA*(FNET+FNST))/(1.D0-F0)/RHOT-1.D0)
  579. DF0=((1.D0-FT+XSRMA*(FNET+FNST))*tr0+DD0*DE0+BB0*H0)/
  580. . (1.D0+BB0*DS1)
  581. *
  582. DCOS0=(DS1*DF0/(2.D0*Ytest))+(SMT*H0)/(2.D0*Ytest*Ytest)
  583. DCOS0=DCOS0*3.D0*SINH(3.D0*SMTM/(2.D0*Ytest))
  584. *
  585. DG0=2.D0*Q0*Q0*ddel0*DF0*FST0
  586. DG0=DG0+2.D0*Q0*FST0*DCOS0
  587. DG0=DG0-2.D0*Q0*ddel0*DF0*COSH(3.D0*SMTM/(2.D0*Ytest))
  588. *
  589. IF (Stest.LT.Gmin0) THEN
  590. denom0=Ytest*Ytest*DG0
  591. ELSE
  592. denom0=6.D0*G0*Stest+Ytest*Ytest*DG0
  593. denom0=denom0+2.D0*H0*Gtest*Ytest
  594. ENDIF
  595. *
  596. IF (denom0.GT.1.D-5) THEN
  597. dlam0=PHIT/denom0
  598. ELSE
  599. IF (recal0.LE.0.D0) THEN
  600. recal0=100.D0
  601. nmax0=2000
  602. iter2=0
  603. GOTO 93
  604. ENDIF
  605. ENDIF
  606. *
  607. *
  608. *
  609. * Vérification des hypothèses de calcul
  610. *===============================================================
  611. *
  612. IF (recal0.LE.0.D0) THEN
  613. *
  614. xnum0=10.D0
  615. *
  616. * Correction SEMI: Mars 2017.
  617. * critere de verification de la precision de la linearisation de la loi GTN
  618. * xlim0=1.D-1
  619. xlim0=1.D-2
  620. * Fin Correction SEMI: Mars 2017.
  621. xmax1=0.D0
  622. xmax2=0.D0
  623. xmax3=0.D0
  624. xmax4=0.D0
  625. xmax5=0.D0
  626. xmax6=0.D0
  627. xmax7=0.D0
  628. *
  629. *---> Condition 3.D0*G*dlam0/Stest << 1
  630. *
  631. IF (Stest.GE.Gmin0) THEN
  632. cri01=ABS(3.D0*G0*dlam0/Stest)
  633. IF (cri01.GT.xlim0) THEN
  634. xmax1=xnum0*cri01+5.D0
  635. ENDIF
  636. ENDIF
  637. *
  638. *---> Condition H0*dlam0/Ytest << 1
  639. *
  640. cri01=ABS(H0*dlam0/Ytest)
  641. IF (cri01.GT.xlim0) THEN
  642. xmax2=xnum0*cri01+5.D0
  643. ENDIF
  644. *
  645. *---> Condition DCOS0*dlam0/COSH << 1
  646. *
  647. cri01=ABS(DCOS0*dlam0/COSH(3.D0*SMTM/(2.D0*Ytest)))
  648. IF (cri01.GT.xlim0) THEN
  649. xmax3=xnum0*cri01+5.D0
  650. ENDIF
  651. *
  652. *---> Condition DF0*dlam0/F << 1
  653. *
  654. Fmax0=MAX(FNS0,FNE0)
  655. Fmax0=MAX(Fmax0,F0)
  656. Fmax0=MAX(Fmax0,FT)
  657. Fmax0=MAX(Fmax0,Fmin0)
  658. cri01=ABS(DF0*dlam0/Fmax0)
  659. IF (cri01.GT.xlim0) THEN
  660. xmax4=xnum0*cri01+5.D0
  661. ENDIF
  662. *
  663. *---> Condition DG0*dlam0/G << 1
  664. *
  665. Gmax0=MAX(ABS(Gtest),1.D0)
  666. cri01=ABS(DG0*dlam0/Gmax0)
  667. IF (cri01.GT.xlim0) THEN
  668. xmax5=xnum0*cri01+5.D0
  669. ENDIF
  670. *
  671. *---> Condition DS1*dlam0/SMT << 1
  672. *
  673. SMTmax=MAX(ABS(SMT),(1.D-3*TRAC(1)))
  674. cri01=ABS(DS1*DF0*dlam0/SMTmax)
  675. IF (cri01.GT.xlim0) THEN
  676. xmax6=xnum0*cri01+5.D0
  677. ENDIF
  678. *
  679. *---> Condition tr0*dlam0/tr(EPSP) << 1
  680. *
  681. tremp0=(FT-F0)/(1.D0-F0)/RHOT
  682. IF (ABS(tremp0).GT.1.D-5) THEN
  683. cri01=ABS(tr0*dlam0/tremp0)
  684. IF (cri01.GT.xlim0) THEN
  685. xmax7=xnum0*cri01+5.D0
  686. ENDIF
  687. ENDIF
  688. *
  689. *
  690. IF (xmax1.GT.xmax2) THEN
  691. xmax00=xmax1
  692. ELSE
  693. xmax00=xmax2
  694. ENDIF
  695. IF (xmax3.GT.xmax00) xmax00=xmax3
  696. IF (xmax4.GT.xmax00) xmax00=xmax4
  697. IF (xmax5.GT.xmax00) xmax00=xmax5
  698. IF (xmax6.GT.xmax00) xmax00=xmax6
  699. IF (xmax7.GT.xmax00) xmax00=xmax7
  700. *
  701. nmax00=NINT(xmax00)
  702. recal0=100.D0
  703. IF (nmax00.GT.nmax0) nmax0=nmax00
  704. IF (nmax0.GT.2000) nmax0=2000
  705. *
  706. *---> On recommence le calcul en sous incrémentant le
  707. * pas de temps
  708. *
  709. IF (nmax0.GE.2) THEN
  710. iter2=0
  711. GOTO 93
  712. ENDIF
  713. *
  714. ENDIF
  715. *
  716. * Cas particulier de la nucléation
  717. *============================================================
  718. *
  719. IF (recal1.LT.0.D0) THEN
  720. *
  721. *---> Condition sur la nucléation pilotée en contrainte
  722. *
  723. IF (FNS0.GT.Fmin0) THEN
  724. *
  725. * Si on se trouve dans la fenetre, on sous incremente
  726. *
  727. dif0=Ytest+H0*dlam0+SMT-DS1*DF0*dlam0
  728. IF ((dif0.GE.(SIGN0-SNS0*2.D0)).AND.
  729. . (VARF(3).LE.(SIGN0+SNS0*2.D0))) THEN
  730. xmax1=200.D0*(dif0-VARF(3))/(4.D0*SNS0)
  731. nmax1=NINT(xmax1)*nmax0
  732. IF (nmax1.GT.2000) nmax1=2000
  733. IF (nmax1.GT.nmax0) THEN
  734. nmax0=nmax1
  735. iter2=0
  736. recal1=100.D0
  737. GOTO 93
  738. ENDIF
  739. ENDIF
  740. *
  741. ENDIF
  742. *
  743. *---> Condition sur la nucléation pilotée en deformation
  744. *
  745. IF (FNE0.GT.Fmin0) THEN
  746. *
  747. * Si on se trouve dans la fenetre, on sous incremente
  748. *
  749. dif0=EPSMT+DE0*dlam0
  750. IF ((dif0.GE.(EPSN0-SNE0*2.D0)).AND.
  751. . (VARF(4).LE.(EPSN0+SNE0*2.D0))) THEN
  752. xmax1=200.D0*(dif0-VARF(4))/(4.D0*SNE0)
  753. nmax1=NINT(xmax1)*nmax0
  754. IF (nmax1.GT.2000) nmax1=2000
  755. IF (nmax1.GT.nmax0) THEN
  756. nmax0=nmax1
  757. iter2=0
  758. recal1=100.D0
  759. GOTO 93
  760. ENDIF
  761. ENDIF
  762. *
  763. ENDIF
  764. *
  765. ENDIF
  766. *
  767. *
  768. * Mise à jour des grandeurs après écoulement plastique
  769. *============================================================
  770. *
  771. *---> Variables internes mises à jour
  772. *
  773. IF (Stest.GT.Gmin0) EPSPT=EPSPT+dlam0
  774. EPSMT=EPSMT+DE0*dlam0
  775. FNST=FNST+BB0*(H0-DS1*DF0)*dlam0
  776. FNET=FNET+DD0*DE0*dlam0
  777. FT=FT+DF0*dlam0
  778. IF (FT.LT.Fmin0) FT=0.D0
  779. *
  780. *---> Contraintes mises à jour
  781. *
  782. SMT=(1.D0-FT)*XK0*((1.D0-FT+XSRMA*(FNET+FNST))/RHOT/(1.D0-F0)-
  783. . 1.D0)
  784. DO 601 I=1,6
  785. A=0.D0
  786. IF (I.LE.3) A=1.D0
  787. IF (Stest.LE.Gmin0) THEN
  788. DEVT(I)=0.D0
  789. ELSE
  790. DEVT(I)=DEVT(I)-3.D0*G0*DEVT(I)*dlam0/Stest
  791. ENDIF
  792. SIGT(I)=DEVT(I)+(SMT*A)
  793. 601 CONTINUE
  794. *
  795. *---> Incrément de déformations plastiques mise à jour
  796. *
  797. treps1=(SMT/(1.D0-FT)-SMT0/(1.D0-VAR0(2)))/XK0
  798. DO 610 I=1,6
  799. A=0.D0
  800. IF (I.LE.3) A=1.D0
  801. depel0=(DEVT(I)-RSIG0(I)+(SMT0*A))/(2.D0*G0)
  802. depel0=depel0+(treps1*A/3.D0)
  803. DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  804. IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  805. DEFT(3)=rdeps3+RDEPS(3)-depel0
  806. ENDIF
  807. 610 CONTINUE
  808. *
  809. *---> Contrainte équivalente Stest mise à jour
  810. *
  811. STEST2=3.D0*PROCON(DEVT,DEVT,6)/2.D0
  812. IF (STEST2.LE.1.D-10) STEST2=0.D0
  813. Stest= STEST2 ** 0.5D0
  814. *
  815. *---> Limite d'élasticité lue sur la courbe d'écrouissage
  816. * mise à jour
  817. *
  818. * - Recherche des valeurs de déformations plastiques équivalentes
  819. * encadrant la valeur courante
  820. *
  821. EPSM0=TRAC(2*NCOURB)
  822. IF (EPSMT.GE.EPSM0) THEN
  823. Y1=TRAC(2*NCOURB-1)
  824. Y2=TRAC(2*NCOURB-1)
  825. EPS1=EPSM0
  826. EPS2=2.D0*EPSM0
  827. ELSE
  828. DO 700 I=1,(NCOURB-1)
  829. EPS1=TRAC(2*I)
  830. EPS2=TRAC(2*(I+1))
  831. IF ((EPSMT.GE.EPS1).AND.(EPSMT.LT.EPS2)) THEN
  832. Y1=TRAC(2*I-1)
  833. Y2=TRAC(2*(I+1)-1)
  834. GOTO 97
  835. ENDIF
  836. 700 CONTINUE
  837. ENDIF
  838. *
  839. * - Limite d'élasticité mise à jour
  840. *
  841. 97 H1=(Y2-Y1)/(EPS2-EPS1)
  842. Ytest=(H1*(EPSMT-EPS2))+Y2
  843. *
  844. *---> Mise à jour des coefficients de nucléation
  845. *
  846. IF ((Ytest+SMT).GE.VARF(3)) THEN
  847. BB0=EXP(-.5D0*((Ytest+SMT-SIGN0)/SNS0)**2.D0)
  848. BB0=BB0*FNS0/(SNS0*((2.D0*XPI)**.5D0))
  849. ELSE
  850. BB0=0.D0
  851. ENDIF
  852. IF (EPSMT.GE.VARF(4)) THEN
  853. DD0=EXP(-.5D0*((EPSMT-EPSN0)/SNE0)**2.D0)
  854. DD0=DD0*FNE0/(SNE0*((2.D0*XPI)**.5D0))
  855. ELSE
  856. DD0=0.D0
  857. ENDIF
  858. *
  859. *---> Fonction d'endommagement
  860. *
  861. IF (FT.LE.FC0) THEN
  862. FST0=FT
  863. ddel0=1.D0
  864. ELSE
  865. FST0=FC0+(FU0-FC0)/(FF0-FC0)*(FT-FC0)
  866. ddel0=(FU0-FC0)/(FF0-FC0)
  867. ENDIF
  868. IF (FST0.GT.FU0) FST0=FU0
  869. *
  870. *---> Terme d'endommagement
  871. *
  872. SMTM=MAX(SMT,0.D0)
  873. SMTM=SMT
  874. Gtest=1.D0+(Q0*FST0)*(Q0*FST0)-
  875. . (2.D0*Q0*FST0*COSH(3.D0*SMTM/(2.D0*Ytest)))
  876. *
  877. *---> Fonction de charge test PHIT mise à jour
  878. *
  879. PHIT=(Stest*Stest)-(Ytest*Ytest)*Gtest
  880. PHI0=ABS(PHIT)
  881. *
  882. *---> Cas particulier des matériaux totalement endommagés
  883. *
  884. 92 IF (FT.GE.FF0) THEN
  885. *
  886. FT=FF0+Fmin0
  887. *
  888. SMT=0.D0
  889. treps1=(SMT/(1.D0-FF0)-SMT0/(1.D0-VAR0(2)))/XK0
  890. DO 13 I=1,6
  891. SIGT(I)=0.D0
  892. A=0.D0
  893. IF (I.LE.3) A=1.D0
  894. depel0=SIGT(I)-RSIG0(I)
  895. DEVT(I)=(depel0-((SMT-SMT0)*A))/(2.D0*G0)
  896. depel0=DEVT(I)+(treps1*A/3.D0)
  897. DEFT(I)=RDEPS0(I)*iter2/nmax0-depel0
  898. IF ((IFOUR.EQ.-2).AND.(I.EQ.3)) THEN
  899. DEFT(3)=rdeps3+RDEPS(3)-depel0
  900. ENDIF
  901. 13 CONTINUE
  902. treps1=DEFT(1)+DEFT(2)+DEFT(3)
  903. DO 15 I=1,6
  904. A=0.D0
  905. IF (I.LE.1) A=1.D0
  906. DEVT(I)=DEFT(I)-(DEFT(1)+DEFT(2)+DEFT(3))*A/3.D0
  907. 15 CONTINUE
  908. IF (treps1.LE.0.D0) THEN
  909. DO 17 I=1,3
  910. DEFT(I)=DEVT(I)
  911. 17 CONTINUE
  912. ENDIF
  913. dlam0=(2.D0*PROCON(DEVT,DEVT,6)/3.D0)**.5D0
  914. EPSPT=EPSPT+dlam0
  915. ** RHOT=RMIN0
  916. *
  917. Gtest=0.D0
  918. PHIT=0.D0
  919. PHI0=0.D0
  920. *
  921. ENDIF
  922. *
  923. *---> Test de convergence ou itération suivante
  924. *
  925. iter0=iter0+1
  926. IF (iter0.LT.200) THEN
  927. GOTO 99
  928. ELSE
  929. PHI0=0.D0
  930. KERRE=460
  931. GOTO 99
  932. ENDIF
  933. *
  934. ENDIF
  935. *
  936. *
  937. * Vérification des sous incrémentations
  938. *
  939. IF ((iter2.LT.nmax0).AND.(KERRE.EQ.0)) THEN
  940. GOTO 95
  941. ENDIF
  942. *
  943. * Passage des deformations de cisaillement exprimées
  944. * en déformations aux déformations de cisaillement exprimées
  945. * en GAMA
  946. *
  947. DO 117 I=1,6
  948. A=1.D0
  949. IF (I.GT.3) A=2.D0
  950. RDEFP(I)=RDEFP(I)*A
  951. 117 CONTINUE
  952. *
  953. *
  954. * Passage a l'option de calcul pour les contraintes
  955. *=========================================================
  956. *
  957. IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN
  958. IF (NSTRS .EQ. 6) THEN
  959. *
  960. *---> MASSIF 3D
  961. *
  962. DO 170 I=1,NSTRS
  963. SIGF(I)=RSIGF(I)
  964. DEFP(I)=RDEFP(I)
  965. DSIGT(I)=RDIGT(I)
  966. 170 CONTINUE
  967. ELSE IF ( NSTRS .EQ. 4 ) THEN
  968. *
  969. *---> Calcul axisymétrique ou contraintes planes
  970. *
  971. DO 180 I=1,NSTRS
  972. SIGF(I)=RSIGF(I)
  973. DEFP(I)=RDEFP(I)
  974. DSIGT(I)=RDIGT(I)
  975. 180 CONTINUE
  976. ENDIF
  977. ENDIF
  978. RETURN
  979. *
  980. END
  981.  
  982.  
  983.  
  984.  
  985.  
  986.  
  987.  
  988.  
  989.  
  990.  
  991.  
  992.  
  993.  
  994.  
  995.  
  996.  
  997.  
  998.  

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