Télécharger coul2.eso

Retour à la liste

Numérotation des lignes :

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

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