Télécharger coul3.eso

Retour à la liste

Numérotation des lignes :

  1. C COUL3 SOURCE CHAT 05/01/12 22:24:13 5004
  2. SUBROUTINE COUL3(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI,
  3. . DEPST,IFOUR,XMAT,NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE)
  4. C--------------------------------------------------------------
  5. C
  6. C PLASTICITE MODELE 3D MOHR COULOMB ELEMENTS JOINTS
  7. C
  8. C ENTREES
  9. C IB = NUMERO DE L'ELEMENT
  10. C IGAUSS = NUMERO DU POINT DE GAUSS
  11. C NSTRS = NOMBRE DE COMPOSANTES DE CONTRAINTES
  12. C SIG0 = CONTRAINTES INITIALES (AU PAS PRECEDENT)
  13. C EPIN0 = DEFORMATIONS INITIALES (AU PAS PRECEDENT)
  14. C VAR0(NVARI) = VARIABLES INTERNES DEBUT
  15. C VAR0(1) = DEFORMATION PLASTIQUE EQUIVALENTE
  16. C VAR0(2) = EPOUN (DEFO. PLASTIQUE DUE A L'OUVERTURE SEULE)
  17. C VAR0(3) = STAT (ETAT DU JOINT)
  18. C VAR0(4) = LAM1 ( MULTIPLICATEUR PLASTIQUE EN TRACTION )
  19. C DSIGT(NSTRS) = INCREMENT DE CONTRAINTES TOTALES
  20. C XMAT(NCOMAT) = COMPOSANTES DE MATERIAU
  21. C NMATT = NOMBRE DE COMPOSANTES DE MATERIAU
  22. C IVAL(NCOMAT) = INDICE DES COMPOSANTES DE MATERIAU
  23. C DD(NSTRS,NSTRS) = MATRICE DE HOOK
  24. C
  25. C SORTIES
  26. C SIGF(NSTRS) = CONTRAINTES FINALES
  27. C DEFP(NSTRS) = INCREMENTS DE DEFORMATIONS INELASTIQUES FINALES
  28. C VARF(NVARI) = VARIABLES INTERNES FINALES
  29. C KERRE = INDICE D'ERREUR
  30. C = 60 PAS DE CONVERGENCE APRES 10 ITERATIONS
  31. C INTERNES (VARIABLE ISSU DANS COULM)
  32. C = 61
  33. C
  34. C-----------------------------------------------------------------------
  35. C
  36. IMPLICIT INTEGER(I-N)
  37. IMPLICIT REAL*8(A-H,O-Z)
  38. C Include contenant quelques constantes dont XPI :
  39. -INC CCREEL
  40. C
  41. DIMENSION SIG0(*),EPIN0(*),VAR0(*),DEPST(*),XMAT(*)
  42. DIMENSION DD(NSTRS,NSTRS)
  43. DIMENSION SIGF(*),DEFP(*),VARF(*)
  44. C
  45. C
  46. C ----------PARAMETRES POUR LE TEST DE NULLITE DES CONTRAINTES ET
  47. C ----------DEFORMATIONS
  48. C
  49. PRESP = +1.D-5
  50. PRESM = -1.D-5
  51. PREEP = +1.D-10
  52. PREEM = -1.D-10
  53. C
  54. C
  55. C ----------PARAMETRE D'IMPPRESSION
  56. C
  57. IMPRI = 0
  58. C
  59. C ----------NOMBRE MAXIMUM D'ITERATIONS INTERNES AUTORISEES
  60. C
  61. MAXITE = 1000
  62. C
  63. C
  64. C------CHARGEMENT DU EPSILON DE PRECISION DES ITERATIONS INTERNES
  65. C
  66. PRECIS = 1.D-6
  67. PRECIN = -1.D-6
  68. C
  69. C------CHARGEMENT DES CARACTERISTIQUES MECANIQUES
  70. C
  71. G1 = XMAT(1)
  72. G2 = XMAT(1)
  73. EC = XMAT(2)
  74. EF = XMAT(5)
  75. ECN = XMAT(6)
  76. C = XMAT(7)
  77. PHI = XMAT(8)
  78. PHI = (PHI*XPI)/180.D0
  79. ccc ZMU = XMAT(8)
  80. ZMU = 0.0D0
  81. ZMU = (ZMU*XPI)/180.D0
  82. FTRC = XMAT(4)
  83. C AM 28/7/95
  84. C FTRC EST EN 4 ET NON EN 9 A CAUSE DE RHO ET ALPH
  85. C VOIR COMMENTAIRE DANS ECOUL2
  86. C
  87. C LA LIMITE EN TRACTION EST MISE BRUTALEMENT A 0
  88. C DES QU'ELLE EST DEPASSEE.
  89. C
  90. C
  91. C------CHARGEMENT DE LA MATRICE DE HOOK INVERSE
  92. C------POUR LE CALCUL DE L'INCREMENT DE DEFORMATION
  93. C------A PARTIR DE L'INCREMENT DE CONTRAINTE
  94. C------ DEPSN = D**-1 DSIGT
  95. C ( D**-1 = DD )
  96. C
  97. DD(1,1) = 1.0D0/G1
  98. DD(1,2) = 0.0D0
  99. DD(1,3) = 0.0D0
  100. DD(2,1) = 0.0D0
  101. DD(2,2) = 1.0D0/G2
  102. DD(2,3) = 0.0D0
  103. DD(3,1) = 0.0D0
  104. DD(3,2) = 0.0D0
  105. DD(3,3) = 1.0D0/EC
  106. C
  107. C contrainte en fin de pente ec
  108. C
  109. sigsei = -ecn * ec
  110. c
  111. C
  112. C-----------INCREMENT DE DEFORMATION
  113. C
  114. DEPS1 = DEPST(1)
  115. DEPS2 = DEPST(2)
  116. DEPSN = DEPST(3)
  117. C
  118. C TRAITEMENT DU PASSAGE DIRECT ENTRE UN ETAT OUVERT
  119. C (I.E. EXISTENCE D'UN JEU AU TEMPS N) A UN ETAT
  120. C COMPRIME AU TEMPS N+1
  121. C
  122. DEFPN=0.0D0
  123. IFIN =0
  124. STOTA1 = 0.0D0
  125. STOTA2 = 0.0D0
  126. DEFP11 = 0.0D0
  127. DEFP12 = 0.0D0
  128. DEFP21 = 0.0D0
  129. DEFP22 = 0.0D0
  130. ALP = 1.0D0
  131. S1INT = 0.0D0
  132. S2INT = 0.0D0
  133. ZJEU = EPIN0(3)
  134. C
  135. C AM CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
  136. C CAS DU RAPPROCHEMENT DES LEVRES
  137. C IF (DEPSN.LT.-1.D-4) THEN
  138. IF (DEPSN.LT.-1.D-4.AND.SIG0(3).GT.PRESP) THEN
  139. C IL EXISTE UN JEU
  140. IF (ZJEU.GT.1.D-3) THEN
  141. C PASSAGE DIRECT D'UN ETAT OUVERT
  142. C A UN ETAT COMPRIME
  143. cc IF (ABS(DEPSN).GT.ZJEU) THEN
  144. valdif=(ABS(DEPSN) - ZJEU)/ABS(DEPSN)
  145. ccc IF (valdif.gt.3.05D-6) THEN
  146. IF (valdif.gt.1.D-3) THEN
  147. C
  148. ALP = ZJEU/(ABS(DEPSN))
  149. STOTA1= DEPS1
  150. DEPS1 = ALP * DEPS1
  151. STOTA2= DEPS2
  152. DEPS2 = ALP * DEPS2
  153. STOCK = DEPSN
  154. DEPSN = ALP * DEPSN
  155. INDIC = 1
  156. GO TO 99
  157. ENDIF
  158. C
  159. ENDIF
  160. C
  161. ENDIF
  162. C
  163. 98 CONTINUE
  164. C
  165. C IL RESTE A PARCOURIR LE TRAJET
  166. C (1-ALP) * DEPSN
  167. C
  168. IF (INDIC.EQ.1) THEN
  169. DEPS1=(1.0D0-ALP)*STOTA1
  170. DEPS2=(1.0D0-ALP)*STOTA2
  171. DEPSN=(1.0D0-ALP)*STOCK
  172. IFIN =1
  173. ENDIF
  174. C
  175. 99 CONTINUE
  176. C
  177. C
  178. C ----CAS GENERAL SANS TRANSITION DIRECTE D'UN ETAT OUVERT A UN ETAT
  179. C ----A UN ETAT COMPRIME
  180. C
  181. C
  182. C
  183. C
  184. IF ((INDIC.EQ.1).OR.(IFIN.EQ.0)) THEN
  185. C
  186. C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS
  187. C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT
  188. C
  189. C------CHARGEMENT DES VALEURS AU PAS PRECEDENT
  190. C
  191. C
  192. C----------DEFORMATIONS
  193. C
  194. C DEFORMATIONS INELASTIQUES
  195. C
  196. EPSS1 = EPIN0(1)
  197. EPSS2 = EPIN0(2)
  198. EPSN1 = EPIN0(3)
  199. C
  200. C DEFORMATIONS ELASTIQUES
  201. C
  202. DEFS10 = DD(1,1)*SIG0(1)
  203. DEFS20 = DD(2,2)*SIG0(2)
  204. c
  205. if (sig0(3).lt.sigsei) then
  206. DEFN0 = ( -1.0d0*ecn ) + ( (sig0(3) - sigsei) / ef )
  207. else if ( (sig0(3).ge.sigsei) .and. (sig0(3).le.presp) ) then
  208. DEFN0 = DD(3,3)*SIG0(3)
  209. else if (sig0(3).gt.presp) then
  210. C
  211. C AM CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
  212. C
  213. C DEFN0 = 0.0D0
  214. DEFN0 = DD(3,3)*SIG0(3)
  215. endif
  216. C
  217. C DEFORMATIONS TOTALES AU PAS PRECEDENT N
  218. C
  219. EPS1N = EPSS1 + DEFS10
  220. EPS2N = EPSS2 + DEFS20
  221. EPSNN = EPSN1 + DEFN0
  222. C
  223. ENDIF
  224. C
  225. C
  226. C
  227. C
  228. C ----CAS DU 2ND PASSAGE : D'UN ETAT FERME ON PASSE A UN ETAT COMPRIME
  229. C
  230. C
  231. C
  232. C
  233. IF (IFIN.EQ.1) THEN
  234. C
  235. C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS
  236. C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT
  237. C
  238. C------CHARGEMENT DES VALEURS AU PAS PRECEDENT
  239. C
  240. C----------DEFORMATIONS
  241. C
  242. C DEFORMATIONS INELASTIQUES
  243. C
  244. EPSS1 = DEFP11
  245. EPSS2 = DEFP21
  246. EPSN1 = 0.0D0
  247. C
  248. C DEFORMATIONS ELASTIQUES
  249. C
  250. DEFS10 = DD(1,1)*S1INT
  251. DEFS20 = DD(2,2)*S2INT
  252. DEFN0 = 0.0D0
  253. C
  254. C DEFORMATIONS TOTALES AU PAS PRECEDENT N
  255. C
  256. EPS1N = EPSS1 + DEFS10
  257. EPS2N = EPSS2 + DEFS20
  258. EPSNN = EPSN1 + DEFN0
  259. C
  260. ENDIF
  261. C
  262. C ESTIMATION DE LA DEFO. TOTALE AU PAS N+1
  263. C
  264. EP1NP1 = EPS1N + DEPS1
  265. EP2NP1 = EPS2N + DEPS2
  266. EPNNP1 = EPSNN + DEPSN
  267. C
  268. C----------EPSILON ELASTIQUE (1 2 N) AU PAS N+1 A L'ITERATION 0
  269. C
  270. EE1P10 = DEFS10 + DEPS1
  271. EE2P10 = DEFS20 + DEPS2
  272. EENP10 = DEFN0 + DEPSN
  273. C
  274. C----------SIGMA (1 2 N) AU PAS N+1 A L'ITERATION 0
  275. C
  276. S1NP10 = G1 * EE1P10
  277. S2NP10 = G2 * EE2P10
  278. c
  279. x6 = -1.0d0 * ecn
  280. c
  281. C AM CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
  282. C
  283. C if (depsn.lt.preem) then
  284. IF (DEPSN.LT.PREEM.AND.SIG0(3).LE.PRESP) THEN
  285. c
  286. c -----------fermeture du joint
  287. c
  288. if (epnnp1.le.preep) then
  289. c
  290. c --------------fermeture du joint avec compression des levres
  291. c
  292. if ( (eenp10.ge.x6) .and. (eenp10.le.preep) ) then
  293. c
  294. c ----------------compression : pente EC
  295. c
  296. SNNP10 = XMAT(2) * EENP10
  297. endif
  298. c
  299. if (eenp10.lt.x6) then
  300. c
  301. c ----------------compression : pente Ef
  302. c
  303. snnp10= ( xmat(2) * x6 ) +
  304. . ( (eenp10 + xmat(6)) * xmat(5) )
  305. endif
  306. c
  307. endif
  308. c
  309. if ((epnnp1.gt.preep).or.(indic.eq.1)) then
  310. c
  311. c --------------rapprochement des levres sans compression
  312. c
  313. SNNP10 = 0.0D0
  314. c
  315. endif
  316. c
  317. ELSE
  318. c
  319. c -----------ouverture du joint OU FERMETURE SI JOINT PAS CASSE
  320. c
  321. SNNP10 = XMAT(2) * EENP10
  322. c
  323. endif
  324. C
  325. C-----------VARIABLES INTERNES
  326. C
  327. C DEFORMATION PLASTIQUE EQUIVALENTE
  328. C
  329. EPSEQ=VAR0(1)
  330. C
  331. C DEFORMATION PLASTIQUE DUE A L'OUVERTURE SEULE
  332. C
  333. EPOUN=VAR0(2)
  334. C
  335. C ETAT DU JOINT
  336. C
  337. C STAT = 0 LEVRES EN CONTACT, AVEC OU SANS CISAILLEMENT
  338. C STAT = 1 JOINT EN OUVERTURE PENTE EC
  339. C STAT = 2 JOINT EN FERMETURE PENTE EC
  340. C STAT = 3 JOINT EN FERMETURE PENTE EF
  341. C
  342. C
  343. C
  344. STAT =VAR0(3)
  345. C
  346. C MULTIPLICATEUR PLASTIQUE EN TRACTION
  347. C
  348. XLAM1=VAR0(4)
  349. C
  350. C MISE A 0 DE LA LIMITE EN TRACTION SI DEJA CASSE OU SI JEU
  351. C
  352. RTRAC=FTRC
  353. IF(XLAM1.GT.0.D0) RTRAC=0.D0
  354.  
  355. C
  356. C
  357. C-----------ETUDE DE L'ECOULEMENT SELON LES DIFFERENTS CRITERES
  358. C
  359. CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
  360. CRINOR = SNNP10 - RTRAC
  361. CRICIP = CALN(S1NP10,S2NP10,SNNP10,PHI,C)
  362. ABSTAU = SQRT( (S1NP10 ** 2) + (S2NP10 ** 2) )
  363. C
  364. C----------AIGUILLAGE SUIVANT LA ZONE 0, 1, 2, 3
  365. C
  366. IF ((CRICIS.LE.PRECIS).AND.(CRINOR.LE.PRECIS)) THEN
  367. GO TO 10
  368. ENDIF
  369. C
  370. IF ((CRICIS.GT.PRECIS).AND.(CRICIP.LE.PRECIS)) THEN
  371. GO TO 40
  372. ENDIF
  373. C
  374. IF ((CRICIP.GT.PRECIS).AND.(ABSTAU.GT.C)) THEN
  375. GO TO 20
  376. ENDIF
  377. C
  378. IF ((ABSTAU.LE.C).AND.(CRINOR.GT.PRECIN)) THEN
  379. GO TO 30
  380. ENDIF
  381. C
  382. C-----------CAS DE LA CONVERGENCE A L ITERATION 0 :
  383. C
  384. 10 CONTINUE
  385. C
  386. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  387. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  388. C
  389. C
  390. IF (IMPRI.EQ.1) THEN
  391. WRITE(6,*) 'PASSAGE N. 1'
  392. WRITE(6,*) 'LES DEUX CRITERES SONT VERIFIES'
  393. WRITE(6,*) 'CONVERGENCE A L ITERATION 0'
  394. ENDIF
  395. C
  396. IF (INDIC.EQ.1) THEN
  397. C FERMETURE DU JOINT SANS COMPRESSION
  398. S1INT = S1NP10
  399. DEFP11= 0.0D0
  400. S2INT = S2NP10
  401. DEFP21= 0.0D0
  402. DEFPN = -epin0(3)
  403. GO TO 98
  404. C
  405. ELSE
  406. C
  407. IF (IFIN.EQ.0) THEN
  408. C
  409. C
  410. c
  411. defp(1)=0.0d0
  412. defp(2)=0.0d0
  413. C
  414. if ( abs(snnp10).le.presp ) then
  415. c
  416. c deformation plastique du pas (n+1) moins
  417. c deformation plastique du pas (n)
  418. c
  419. defp(3) = epnnp1 - epsn1
  420. c
  421. else
  422. c
  423. c le joint se deforme toujours elastiquement en fermeture
  424. c
  425. DEFP(3) = 0.0D0
  426. c
  427. endif
  428. c
  429. ELSE
  430. c
  431. c cas du passage direct d'un etat ouvert a un etat
  432. c comprime
  433. c
  434. DEFP12 = 0.0D0
  435. DEFP(1) = 0.0d0
  436. DEFP22 = 0.0D0
  437. DEFP(2) = 0.0d0
  438. cc DEFP(3) = 0.0d0
  439. DEFP(3) = defpn
  440. IFIN = 0
  441. ENDIF
  442. C
  443. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  444. C
  445. c
  446. c calcul du delta epsilon plastique entre
  447. c le pas actuel et le pas precedent
  448. c
  449. delep1 = defp(1)
  450. delep2 = defp(2)
  451. delepn = defp(3)
  452. c
  453. c produit contracte pour le calcul de la defo. plasti. cumulee
  454. c
  455. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  456. depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  457. VARF(1) = EPSEQ + DEPEQ
  458. VARF(4) = XLAM1
  459. c
  460. c calcul de la deformation normale due a l'ouverture seule
  461. c
  462. if ( abs(snnp10).le.presp ) then
  463. epoun = defp(3)
  464. VARF(2) = EPOUN
  465. else
  466. epoun = 0.0d0
  467. VARF(2) = EPOUN
  468. endif
  469. c
  470. c etat du joint stat = 0 cisaillement pur, levres en contact
  471. c stat = 1 ouverture du joint pente EC
  472. c stat = 2 compression du joint pente EC
  473. c stat = 3 compression du joint pente EF
  474. c
  475. if ( abs(snnp10).le.presp) then
  476. stat = 1.D0
  477. if ( abs(epnnp1).le.1.D-3 ) then
  478. stat = 0.D0
  479. endif
  480. VARF(3) = STAT
  481. endif
  482. if ( (snnp10.gt.sigsei).and.(snnp10.lt.presm) ) then
  483. stat = 2.D0
  484. VARF(3) = STAT
  485. endif
  486. if (snnp10.le.sigsei) then
  487. stat = 3.D0
  488. VARF(3) = STAT
  489. endif
  490. C
  491. SIGF(1) = S1NP10
  492. SIGF(2) = S2NP10
  493. SIGF(3) = SNNP10
  494. C
  495. RETURN
  496. c
  497. ENDIF
  498. C
  499. C----------CRITERES DE L'EFFORT NORMAL ET DE CISAILLEMENT
  500. C----------TOUS LES DEUX VIOLES EN MEME TEMPS
  501. C
  502. 20 CONTINUE
  503. C
  504. C ECOULEMENT SELON LES DEUX CRITERES A L'ITERATION ITER
  505. C
  506. ITER = 1
  507. DELAM1=0.D0
  508. C
  509. 2 CONTINUE
  510. C
  511. IF (ITER.GT.MAXITE) THEN
  512. WRITE(6,*) ' LES DEUX CRITERES SONT VIOLES '
  513. WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
  514. STOP
  515. ENDIF
  516. C
  517. IF (IMPRI.EQ.1) THEN
  518. write(6,*) 'PASSAGE N. 2'
  519. write(6,*) 'LES DEUX CRITERES SONT VIOLES'
  520. write(6,*) 'ITERATION',ITER
  521. ENDIF
  522. C
  523. C
  524. C CALCUL DES CORRECTEURS PLASTIQUES
  525. C
  526. C
  527. C DERIVEE DU CRITERE F SUIVANT SIGMA
  528. C
  529. DFIDS1 = 0.0D0
  530. DFIDS2 = 0.0D0
  531. DFIDSN = 1.0D0
  532. C
  533. C DERIVEE DU CRITERE G SUIVANT SIGMA
  534. C
  535. DENO = SQRT ( (S1NP10**2) + (S2NP10**2) )
  536. DGIDS1 = S1NP10 / DENO
  537. DGIDS2 = S2NP10 / DENO
  538. DGIDSN = TAN(PHI)
  539. DZGIDS = TAN(ZMU)
  540. C
  541. C MATRICE DE HOOK * DFI/DSIGMA
  542. C
  543. HODF1 = G1 * DFIDS1
  544. HODF2 = G2 * DFIDS2
  545. HODFN = EC * DFIDSN
  546. C
  547. C MATRICE DE HOOK * DGI/DSIGMA
  548. C
  549. HODG1 = G1 * DGIDS1
  550. HODG2 = G2 * DGIDS2
  551. HODGN = EC * DZGIDS
  552. C
  553. C FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  554. C
  555. FIFI = (DFIDS1 * HODF1) + (DFIDS2 * HODF2) + (DFIDSN * HODFN)
  556. C
  557. C FIGI = DFI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
  558. C
  559. FIGI = (DFIDS1 * HODG1) + (DFIDS2 * HODG2) + (DFIDSN * HODGN)
  560. C
  561. C GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  562. C
  563. GIFI = (DGIDS1 * HODF1) + (DGIDS2 * HODF2) + (DGIDSN * HODFN)
  564. C
  565. C GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  566. C
  567. GIGI = (DGIDS1 * HODG1) + (DGIDS2 * HODG2) + (DGIDSN * HODGN)
  568. C
  569. C RESOLUTION DU SYSTEME POUR TROUVER LES CORRECTEURS
  570. C PLASTIQUES DLAM1 ET DLAM2
  571. C
  572. DET = (FIFI * GIGI) - (FIGI * GIFI)
  573. C
  574. IF (ABS(DET).LE.PRECIS) THEN
  575. WRITE(6,*) ' SYSTEME INDETERMINE (F G): ITERATION ',ITER
  576. STOP
  577. ENDIF
  578. C
  579. DLAM1 = (CRINOR * GIGI) - (CRICIS * FIGI)
  580. DLAM1 = DLAM1 / DET
  581. DELAM1 = DELAM1 + DLAM1
  582. RTRAC=0.D0
  583. C
  584. DLAM2 = (CRICIS * FIFI) - (CRINOR * GIFI)
  585. DLAM2 = DLAM2 / DET
  586. C
  587. C CALCUL DU NOUVEL ETAT DE CONTRAINTES
  588. C
  589. S1NP11=S1NP10 - (G1 * ( (DLAM1*DFIDS1) + (DLAM2*DGIDS1)))
  590. S2NP11=S2NP10 - (G2 * ( (DLAM1*DFIDS2) + (DLAM2*DGIDS2)))
  591. SNNP11=SNNP10 - (EC * ( (DLAM1*DFIDSN) + (DLAM2*DZGIDS)))
  592. C
  593. C ALLER EN 52 POUR LE TEST DE LA CONVERGENCE
  594. C
  595. C
  596. 52 CONTINUE
  597. C
  598. C CONVERGENCE CHECK
  599. C
  600. CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
  601. CRINOR = SNNP11 - RTRAC
  602. C
  603. IF ((ABS(CRICIS).LE.PRECIS).AND.(ABS(CRINOR).LE.PRECIS)) THEN
  604. C
  605. C----------LES DEUX CRITERES SONT VERIFIES
  606. C
  607. C
  608. IF (IMPRI.EQ.1) THEN
  609. write(6,*) 'PASSAGE N. 5 - 2'
  610. write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
  611. write(6,*) 'CONVERGENCE A L ITERATION',ITER
  612. ENDIF
  613. C
  614. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  615. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  616. C
  617. C
  618. SIGF(1) = S1NP11
  619. SIGF(2) = S2NP11
  620. SIGF(3) = SNNP11
  621. C
  622. DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
  623. DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
  624. c
  625. if ( abs(sigf(3)) .le. presp ) then
  626. c
  627. c deformation plastique du pas (n+1) moins
  628. c deformation plastique du pas (n)
  629. c
  630. DEFP(3) = EPNNP1 - EPSN1
  631. else
  632. c
  633. c le joint se deforme toujours elastiquement en fermeture
  634. c
  635. DEFP(3) = 0.0D0
  636. c
  637. endif
  638. C
  639. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  640. C
  641. c
  642. c calcul du delta epsilon plastique entre
  643. c le pas actuel et le pas precedent
  644. c
  645. delep1 = defp(1)
  646. delep2 = defp(2)
  647. delepn = defp(3)
  648. c
  649. c produit contracte pour le calcul de la defo. plasti. cumulee
  650. c
  651. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  652. depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  653. VARF(1) = EPSEQ + DEPEQ
  654. VARF(4) = XLAM1 + DELAM1
  655. c
  656. c calcul de la deformation normale due a l'ouverture seule
  657. c
  658. if ( abs(snnp11).le.presp ) then
  659. epoun = delepn
  660. VARF(2) = EPOUN
  661. else
  662. epoun = 0.0d0
  663. VARF(2) = EPOUN
  664. endif
  665. c
  666. c etat du joint
  667. c stat = 0 levres en contact, avec ou sans cisaillement
  668. c stat = 1 ouverture du joint pente EC
  669. c stat = 2 compression du joint pente EC
  670. c stat = 3 compression du joint pente EF
  671. c
  672. if ( abs(snnp11).le.presp ) then
  673. stat = 1.D0
  674. if ( abs(epnnp1).le.1.D-8 ) then
  675. stat = 0.D0
  676. endif
  677. VARF(3) = STAT
  678. endif
  679. if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
  680. stat = 2.D0
  681. VARF(3) = STAT
  682. endif
  683. if (snnp11.le.sigsei) then
  684. stat = 3.D0
  685. VARF(3) = STAT
  686. endif
  687. c
  688. RETURN
  689. ELSE
  690. C
  691. S1NP10=S1NP11
  692. S2NP10=S2NP11
  693. SNNP10=SNNP11
  694. CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
  695. CRINOR = SNNP10 - RTRAC
  696. ITER = ITER + 1
  697. C
  698. GO TO 2
  699. ENDIF
  700. C
  701. C
  702. C____________________________________________________________________
  703. C
  704. C
  705. C----------CRITERE DE L'EFFORT NORMAL VIOLE MAIS
  706. C----------CRITERE DE CISAILLEMENT VERIFIE
  707. C
  708. 30 CONTINUE
  709. C
  710. C ECOULEMENT SELON LE CRITERE F A L'ITERATION ITER
  711. C
  712. ITER = 1
  713. DELAM1=0.D0
  714. C
  715. 3 CONTINUE
  716. C
  717. IF (ITER.GT.MAXITE) THEN
  718. WRITE(6,*) ' LE CRITERE NORMAL EST VIOLE '
  719. WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
  720. STOP
  721. ENDIF
  722. C
  723. IF (IMPRI.EQ.1) THEN
  724. write(6,*) 'PASSAGE N. 3'
  725. write(6,*) 'CRITERE DE L EFFORT NORMAL VIOLE'
  726. write(6,*) 'ITERATION',ITER
  727. ENDIF
  728. C
  729. C
  730. C
  731. C CALCUL DU CORRECTEUR PLASTIQUE
  732. C
  733. C
  734. C DERIVEE DU CRITERE F SUIVANT SIGMA
  735. C
  736. DFIDS1 = 0.0D0
  737. DFIDS2 = 0.0D0
  738. DFIDSN = 1.0D0
  739. C
  740. C MATRICE DE HOOK * DGI/DSIGMA
  741. C
  742. HODF1 = G1 * DFIDS1
  743. HODF2 = G2 * DFIDS2
  744. HODFN = EC * DFIDSN
  745. C
  746. C FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  747. C
  748. FIFI = (DFIDS1 * HODF1) + (DFIDS2 * HODF2) + (DFIDSN * HODFN)
  749. C
  750. C RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
  751. C PLASTIQUE DLAM1
  752. C
  753. IF (ABS(FIFI).LE.PRECIS) THEN
  754. WRITE(6,*) ' SYSTEME INDETERMINE (F) : ITERATION ',ITER
  755. STOP
  756. ENDIF
  757. C
  758. DLAM1 = CRINOR / FIFI
  759. DELAM1 = DELAM1 + DLAM1
  760. RTRAC = 0.D0
  761. C
  762. C CALCUL DU NOUVEL ETAT DE CONTRAINTES
  763. C
  764. S1NP11=S1NP10 - (XMAT(1) * (DLAM1*DFIDS1) )
  765. S2NP11=S2NP10 - (XMAT(1) * (DLAM1*DFIDS2) )
  766. SNNP11=SNNP10 - ( EC * (DLAM1*DFIDSN) )
  767. C
  768. C ALLER EN 53 POUR LE TEST DE LA CONVERGENCE
  769. C
  770. C
  771. C
  772. 53 CONTINUE
  773. C
  774. C CONVERGENCE CHECK
  775. C
  776. CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
  777. CRINOR = SNNP11 - RTRAC
  778. C
  779. IF (ABS(CRINOR).LE.PRECIS) THEN
  780. C
  781. C----------LES DEUX CRITERES SONT VERIFIES
  782. C
  783. C
  784. IF (IMPRI.EQ.1) THEN
  785. write(6,*) 'PASSAGE N. 5 - 3'
  786. write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
  787. write(6,*) 'CONVERGENCE A L ITERATION',ITER
  788. ENDIF
  789. C
  790. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  791. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  792. C
  793. C
  794. SIGF(1) = S1NP11
  795. SIGF(2) = S2NP11
  796. SIGF(3) = SNNP11
  797. C
  798. DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
  799. DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
  800. c
  801. if ( abs(sigf(3)) .le. presp ) then
  802. c
  803. c deformation plastique du pas (n+1) moins
  804. c deformation plastique du pas (n)
  805. c
  806. DEFP(3) = EPNNP1 - EPSN1
  807. else
  808. c
  809. c le joint se deforme toujours elastiquement en fermeture
  810. c
  811. DEFP(3) = 0.0D0
  812. c
  813. endif
  814. C
  815. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  816. C
  817. c
  818. c calcul du delta epsilon plastique entre
  819. c le pas actuel et le pas precedent
  820. c
  821. delep1 = defp(1)
  822. delep2 = defp(2)
  823. delepn = defp(3)
  824. c
  825. c produit contracte pour le calcul de la defo. plasti. cumulee
  826. c
  827. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  828. depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  829. VARF(1) = EPSEQ + DEPEQ
  830. VARF(4) = XLAM1 + DELAM1
  831. c
  832. c calcul de la deformation normale due a l'ouverture seule
  833. c
  834. if ( abs(snnp11).le.presp ) then
  835. epoun = delepn
  836. VARF(2) = EPOUN
  837. else
  838. epoun = 0.0d0
  839. VARF(2) = EPOUN
  840. endif
  841. c
  842. c etat du joint
  843. c stat = 0 levres en contact, avec ou sans cisaillement
  844. c stat = 1 ouverture du joint pente EC
  845. c stat = 2 compression du joint pente EC
  846. c stat = 3 compression du joint pente EF
  847. c
  848. if ( abs(snnp11).le.presp ) then
  849. stat = 1.D0
  850. if ( abs(epnnp1).le.1.D-8 ) then
  851. stat = 0.D0
  852. endif
  853. VARF(3) = STAT
  854. endif
  855. if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
  856. stat = 2.D0
  857. VARF(3) = STAT
  858. endif
  859. if (snnp11.le.sigsei) then
  860. stat = 3.D0
  861. VARF(3) = STAT
  862. endif
  863. c
  864. RETURN
  865. ELSE
  866. C
  867. S1NP10=S1NP11
  868. S2NP10=S2NP11
  869. SNNP10=SNNP11
  870. CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
  871. CRINOR = SNNP10 - RTRAC
  872. ITER = ITER + 1
  873. C
  874. GO TO 3
  875. ENDIF
  876. C
  877. C
  878. C____________________________________________________________________
  879. C
  880. C
  881. C
  882. C----------CRITERES DE CISAILLEMENT VIOLE MAIS
  883. C----------CRITERE DE L'EFFORT NORMAL VERIFIE
  884. C
  885. 40 CONTINUE
  886. C
  887. C ECOULEMENT SELON LE CRITERE G A L'ITERATION ITER
  888. C
  889. ITER = 1
  890. C
  891. 4 CONTINUE
  892. C
  893. IF (ITER.GT.MAXITE) THEN
  894. WRITE(6,*) ' LE CRITERE DE CISAILLEMENT EST VIOLE '
  895. WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
  896. STOP
  897. ENDIF
  898. C
  899. IF (IMPRI.EQ.1) THEN
  900. write(6,*) 'PASSAGE N. 4'
  901. write(6,*) 'CRITERE DE CISAILLEMENT VIOLE'
  902. write(6,*) 'ITERATION',ITER
  903. ENDIF
  904. C
  905. C
  906. C CALCUL DU CORRECTEUR PLASTIQUE
  907. C
  908. C
  909. C DERIVEE DU CRITERE G SUIVANT SIGMA
  910. C
  911. DENO = SQRT ( (S1NP10**2) + (S2NP10**2) )
  912. DGIDS1 = S1NP10 / DENO
  913. DGIDS2 = S2NP10 / DENO
  914. DZGIDS = TAN(ZMU)
  915. DGIDSN = TAN(PHI)
  916. C
  917. C MATRICE DE HOOK * DGI/DSIGMA
  918. C
  919. HODG1 = G1 * DGIDS1
  920. HODG2 = G2 * DGIDS2
  921. HODGN = EC * DZGIDS
  922. C
  923. C GIGI = DGI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
  924. C
  925. GIGI = (DGIDS1 * HODG1) + (DGIDS2 * HODG2) + (DGIDSN * HODGN)
  926. C
  927. C RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
  928. C PLASTIQUE DLAM2
  929. C
  930. IF (ABS(GIGI).LE.PRECIS) THEN
  931. WRITE(6,*) ' SYSTEME INDETERMINE (G) : ITERATION ',ITER
  932. STOP
  933. ENDIF
  934. C
  935. DLAM2 = CRICIS / GIGI
  936. C
  937. C CALCUL DU NOUVEL ETAT DE CONTRAINTES
  938. C
  939. cc SNNP11=SNNP10 - (XMAT(2) * (DLAM2*DGIDSN) )
  940. C
  941. S1NP11=S1NP10 - (G1 * (DLAM2*DGIDS1) )
  942. S2NP11=S2NP10 - (G2 * (DLAM2*DGIDS2) )
  943. SNNP11=SNNP10 - (EC * (DLAM2*DZGIDS) )
  944. c
  945. if (SNNP11.ge.0.0d0) then
  946. SNNP11 = 0.0d0
  947. endif
  948. c
  949. cccssszzzeeerrr if (SNNP11.GE.0.0d0) THEN
  950. cccssszzzeeerrr SNNP11 = 0.0d0
  951. cccssszzzeeerrr ENDIF
  952. cccssszzzeeerrrc
  953. cccssszzzeeerrr if (abs(zmu).le.presp) then
  954. cccssszzzeeerrr differ = abs( depsn ) - epsn1
  955. cccssszzzeeerrr differ = abs(differ)/abs( depsn )
  956. cccssszzzeeerrr if (depsn.lt.preem) then
  957. cccssszzzeeerrr if (epnnp1.gt.preem) then
  958. cccssszzzeeerrr SNNP11 = 0.0d0
  959. cccssszzzeeerrr endif
  960. cccssszzzeeerrr else
  961. cccssszzzeeerrrc
  962. cccssszzzeeerrr SNNP11 = 0.0d0
  963. cccssszzzeeerrr endif
  964. cccssszzzeeerrr endif
  965. C
  966. C ALLER EN 54 POUR LE TEST DE LA CONVERGENCE
  967. C
  968. C
  969. C
  970. 54 CONTINUE
  971. C
  972. C CONVERGENCE CHECK
  973. C
  974. CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
  975. CRINOR = SNNP11 - RTRAC
  976. C
  977. IF (ABS(CRICIS).LE.PRECIS) THEN
  978. C
  979. IF (INDIC.EQ.1) THEN
  980. C
  981. C----------1ER TRAJET : ON PARCOURT ALP * DEPSN
  982. C
  983. DEFPN = -epin0(3)
  984. S1INT = S1NP11
  985. S2INT = S2NP11
  986. DEFP11 = EP1NP1 - ( DD(1,1) * S1NP11 )
  987. DEFP21 = EP2NP1 - ( DD(2,2) * S2NP11 )
  988. GO TO 98
  989. ELSE
  990. IF (IFIN.EQ.0) THEN
  991. C
  992. C
  993. C----------LES DEUX CRITERES SONT VERIFIES
  994. C
  995. C
  996. IF (IMPRI.EQ.1) THEN
  997. write(6,*) 'PASSAGE N. 5 - 4'
  998. write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
  999. write(6,*) 'CONVERGENCE A L ITERATION',ITER
  1000. ENDIF
  1001. C
  1002. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  1003. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  1004. C
  1005. C
  1006. SIGF(1) = S1NP11
  1007. SIGF(2) = S2NP11
  1008. SIGF(3) = SNNP11
  1009. C
  1010. DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
  1011. DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
  1012. c
  1013. if ( abs(sigf(3)) .le. presp ) then
  1014. c
  1015. c deformation plastique du pas (n+1) moins
  1016. c deformation plastique du pas (n)
  1017. c
  1018. DEFP(3) = EPNNP1 - EPSN1
  1019. else
  1020. c
  1021. c le joint se deforme toujours elastiquement en fermeture
  1022. c
  1023. DEFP(3) = 0.0D0
  1024. c
  1025. endif
  1026. ELSE
  1027. C
  1028. C----------CAS DU PASSAGE DIRECT D'UN ETAT OUVERT A UN ETAT COMPRIME
  1029. C
  1030. DEFP12 = EP1NP1 - (DD(1,1)*S1NP11)
  1031. DEFP22 = EP2NP1 - (DD(2,2)*S2NP11)
  1032. DEFP(1) = DEFP12 - epin0(1)
  1033. DEFP(2) = DEFP22 - epin0(2)
  1034. DEFP(3) = -epin0(3)
  1035. SIGF(1) = S1NP11
  1036. SIGF(2) = S2NP11
  1037. SIGF(3) = SNNP11
  1038. IFIN = 0
  1039. ENDIF
  1040. C
  1041. ENDIF
  1042. C
  1043. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  1044. C
  1045. c
  1046. c calcul du delta epsilon plastique entre
  1047. c le pas actuel et le pas precedent
  1048. c
  1049. delep1 = defp(1)
  1050. delep2 = defp(2)
  1051. delepn = defp(3)
  1052. c
  1053. c produit contracte pour le calcul de la defo. plasti. cumulee
  1054. c
  1055. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  1056. depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  1057.  
  1058. VARF(1) = EPSEQ + DEPEQ
  1059. VARF(4) = XLAM1
  1060. c
  1061. c calcul de la deformation normale due a l'ouverture seule
  1062. c
  1063. if ( abs(snnp11).le.presp ) then
  1064. epoun = delepn
  1065. VARF(2) = EPOUN
  1066. else
  1067. epoun = 0.0d0
  1068. VARF(2) = EPOUN
  1069. endif
  1070. c
  1071. c etat du joint
  1072. c stat = 0 levres en contact, avec ou sans cisaillement
  1073. c stat = 1 ouverture du joint pente EC
  1074. c stat = 2 compression du joint pente EC
  1075. c stat = 3 compression du joint pente EF
  1076. c
  1077. if ( abs(snnp11).le.presp ) then
  1078. stat = 1.D0
  1079. if ( abs(epnnp1).le.1.D-8 ) then
  1080. stat = 0.D0
  1081. endif
  1082. VARF(3) = STAT
  1083. endif
  1084. if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
  1085. stat = 2.D0
  1086. VARF(3) = STAT
  1087. endif
  1088. if (snnp11.le.sigsei) then
  1089. stat = 3.D0
  1090. VARF(3) = STAT
  1091. endif
  1092. c
  1093. RETURN
  1094. ELSE
  1095. C
  1096. S1NP10=S1NP11
  1097. S2NP10=S2NP11
  1098. SNNP10=SNNP11
  1099. CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
  1100. CRINOR = SNNP10 - RTRAC
  1101. ITER = ITER + 1
  1102. C
  1103. GO TO 4
  1104. ENDIF
  1105. C
  1106. C
  1107. C____________________________________________________________________
  1108. C
  1109. RETURN
  1110. END
  1111.  
  1112.  
  1113.  
  1114.  
  1115.  
  1116.  
  1117.  
  1118.  
  1119.  

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