Télécharger coul1.eso

Retour à la liste

Numérotation des lignes :

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

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