Télécharger coul1.eso

Retour à la liste

Numérotation des lignes :

  1. C COUL1 SOURCE CB215821 16/04/21 21:16:05 8920
  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 (IDINT(XMAT(IDEBMA+NSTRS+5)).EQ.1) THEN
  86. IDIRN=1
  87. ELSE IF (IDINT(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. <3pan Styme="color: #666666; font-style: italic;">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. GO TO 40
  463. ENDIF
  464. C
  465. * IF ((CRICIP.GT.PRECIS).AND.(ABSTAU.GT.C)) THEN
  466. IF ((CRICIP.GT.PRECIS).AND.(ABSTAU.GT.PRECIS)) THEN
  467. GO TO 20
  468. ENDIF
  469. C
  470. * IF ((ABSTAU.LE.C).AND.(CRINOR.GT.PRECIN)) THEN
  471. IF ((ABSTAU.LE.PRECIS).AND.(CRINOR.GT.PRECIN)) THEN
  472. GO TO 30
  473. ENDIF
  474. C
  475. C-----------CAS DE LA CONVERGENCE A L ITERATION 0 :
  476. C
  477. 10 CONTINUE
  478. C
  479. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  480. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  481. C
  482. C
  483. IF (IMPRI.EQ.1) THEN
  484. WRITE(6,*) 'ETIQUETTE 10 - LES DEUX CRITERES SONT VERIFIES'
  485. ENDIF
  486. C
  487. *- cs deb
  488. SNNP11=SNNP10
  489. *- cs fin
  490. C
  491. IF (INDIC.EQ.1) THEN
  492. C FERMETURE DU JOINT SANS COMPRESSION
  493. * print *,'10 ok1'
  494. S1INT = S1NP10
  495. DEFP11= 0.0D0
  496. S2INT = S2NP10
  497. DEFP21= 0.0D0
  498. DEFPN = -(epin0(3)-DEFPNC)
  499. GO TO 98
  500. C
  501. ELSE
  502. C
  503. IF (IFIN.EQ.0) THEN
  504. C
  505. C
  506. c
  507. defp(1)=0.0d0
  508. defp(2)=0.0d0
  509. C
  510. if ( abs(snnp10).le.presp ) then
  511. * print *,'10 ok2'
  512. c
  513. c deformation plastique du pas (n+1) moins
  514. c deformation plastique du pas (n)
  515. c
  516. defp(3) = epnnp1 - epsn1
  517. c
  518. else
  519. * print *,'10 ok3'
  520. c
  521. c comportement elastoplastique du joint en fermeture
  522. c
  523. *- cs deb
  524. * DEFP(3) = 0.0D0
  525. CALL PLAQP(EC,QT,FNE,SNNP10,EPCO,XLAM3,SNNP11,EPCOF,XLAM3F)
  526. DEFP(3) = EPCOF-EPCO
  527. *- cs fin
  528. c
  529. endif
  530. c
  531. ELSE
  532. * print *,'10 ok4'
  533. c
  534. c cas du passage direct d'un etat ouvert a un etat
  535. c comprime
  536. c
  537. DEFP12 = 0.0D0
  538. DEFP(1) = 0.0d0
  539. DEFP22 = 0.0D0
  540. DEFP(2) = 0.0d0
  541. cc DEFP(3) = 0.0d0
  542. *- cs deb
  543. * DEFP(3) = defpn
  544. CALL PLAQP(EC,QT,FNE,SNNP10,EPCO,XLAM3,SNNP11,EPCOF,XLAM3F)
  545. DEFP(3) = DEFPN+(EPCOF-EPCO)
  546.  
  547. *- cs fin
  548. IFIN = 0
  549. ENDIF
  550. C
  551. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  552. C
  553. c
  554. c calcul du delta epsilon plastique entre
  555. c le pas actuel et le pas precedent
  556. c
  557. delep1 = defp(1)
  558. delep2 = defp(2)
  559. delepn = defp(3)
  560. c
  561. c produit contracte pour le calcul de la defo. plasti. cumulee
  562. c
  563. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  564. *- cs depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  565. depeq = SQRT(pro1)
  566. VARF(1) = EPSEQ + DEPEQ
  567. VARF(4) = XLAM1
  568. VARF(9) = XLAM3F
  569. c
  570. c calcul de la deformation normale due a l'ouverture seule
  571. c
  572. if ( abs(snnp11).le.presp ) then
  573. *- cs deb
  574. epoun = defp(3)-(EPCOF-EPCO)
  575. *- cs fin
  576. VARF(2) = EPOUN
  577. else
  578. * epoun = 0.0d0
  579. epoun = DEFPN
  580. VARF(2) = EPOUN
  581. endif
  582. C
  583. *- cs deb
  584. C deformation normale plastique due a la compression
  585. C
  586. VARF(8)=EPCOF
  587. C
  588. *- cs fin
  589. c
  590. c etat du joint stat = 0 cisaillement pur, levres en contact
  591. c stat = 1 ouverture du joint contrainte nulle
  592. c stat = 2 compression du joint elasticite pente EC
  593. c stat = 3 compression du joint plasticite
  594. c
  595. stat = 2.D0
  596. eptmco = epnnp1-EPCOF
  597. if ( abs(snnp11).le.presp) then
  598. stat = 1.D0
  599. if ( abs(eptmco).le.1.D-8 ) stat = 0.D0
  600. endif
  601. F3=CRIPQ(EC,QT,FNE,SNNP11,XLAM3F)
  602. if (ABS(F3).LE.PRECIS) stat = 3.D0
  603. VARF(3) = STAT
  604. C
  605. SIGF(1) = S1NP10
  606. SIGF(2) = S2NP10
  607. SIGF(3) = SNNP11
  608. C
  609. GOTO 9900
  610. c
  611. ENDIF
  612. C
  613. C----------CRITERES DE L'EFFORT NORMAL ET DE CISAILLEMENT
  614. C----------TOUS LES DEUX VIOLES EN MEME TEMPS
  615. C
  616. 20 CONTINUE
  617. C
  618. C ECOULEMENT SELON LES DEUX CRITERES A L'ITERATION ITER
  619. C
  620. ITER = 1
  621. DELAM1=0.D0
  622. C
  623. 2 CONTINUE
  624. C
  625. IF (ITER.GT.MAXITE) THEN
  626. WRITE(6,*) ' LES DEUX CRITERES SONT VIOLES '
  627. WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
  628. STOP
  629. ENDIF
  630. C
  631. IF (IMPRI.EQ.1) THEN
  632. write(6,*) 'ETIQUETTE 20 - LES DEUX CRITERES SONT VIOLES'
  633. ENDIF
  634. C
  635. C
  636. C CALCUL DES CORRECTEURS PLASTIQUES
  637. C
  638. C
  639. C DERIVEE DU CRITERE F SUIVANT SIGMA
  640. C
  641. DFIDS1 = 0.0D0
  642. DFIDS2 = 0.0D0
  643. DFIDSN = 1.0D0
  644. C
  645. C DERIVEE DU CRITERE G SUIVANT SIGMA
  646. C
  647. DENO = SQRT ( (S1NP10**2) + (S2NP10**2) )
  648. DGIDS1 = S1NP10 / DENO
  649. DGIDS2 = S2NP10 / DENO
  650. DGIDSN = TAN(PHI)
  651. DZGIDS = TAN(ZMU)
  652. C
  653. C MATRICE DE HOOK * DFI/DSIGMA
  654. C
  655. HODF1 = G1 * DFIDS1
  656. HODF2 = G2 * DFIDS2
  657. HODFN = EC * DFIDSN
  658. C
  659. C MATRICE DE HOOK * DGI/DSIGMA
  660. C
  661. HODG1 = G1 * DGIDS1
  662. HODG2 = G2 * DGIDS2
  663. HODGN = EC * DZGIDS
  664. C
  665. C FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  666. C
  667. FIFI = (DFIDS1 * HODF1) + (DFIDS2 * HODF2) + (DFIDSN * HODFN)
  668. C
  669. C FIGI = DFI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
  670. C
  671. FIGI = (DFIDS1 * HODG1) + (DFIDS2 * HODG2) + (DFIDSN * HODGN)
  672. C
  673. C GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  674. C
  675. GIFI = (DGIDS1 * HODF1) + (DGIDS2 * HODF2) + (DGIDSN * HODFN)
  676. C
  677. C GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  678. C
  679. GIGI = (DGIDS1 * HODG1) + (DGIDS2 * HODG2) + (DGIDSN * HODGN)
  680. C
  681. C RESOLUTION DU SYSTEME POUR TROUVER LES CORRECTEURS
  682. C PLASTIQUES DLAM1 ET DLAM2
  683. C
  684. DET = (FIFI * GIGI) - (FIGI * GIFI)
  685. C
  686. IF (ABS(DET).LE.PRECIS) THEN
  687. WRITE(6,*) ' SYSTEME INDETERMINE (F G): ITERATION ',ITER
  688. STOP
  689. ENDIF
  690. C
  691. DLAM1 = (CRINOR * GIGI) - (CRICIS * FIGI)
  692. DLAM1 = DLAM1 / DET
  693. DELAM1 = DELAM1 + DLAM1
  694. RTRAC=0.D0
  695. C
  696. DLAM2 = (CRICIS * FIFI) - (CRINOR * GIFI)
  697. DLAM2 = DLAM2 / DET
  698. C
  699. C CALCUL DU NOUVEL ETAT DE CONTRAINTES
  700. C
  701. S1NP11=S1NP10 - (G1 * ( (DLAM1*DFIDS1) + (DLAM2*DGIDS1)))
  702. S2NP11=S2NP10 - (G2 * ( (DLAM1*DFIDS2) + (DLAM2*DGIDS2)))
  703. SNNP11=SNNP10 - (EC * ( (DLAM1*DFIDSN) + (DLAM2*DZGIDS)))
  704. C
  705. C ALLER EN 52 POUR LE TEST DE LA CONVERGENCE
  706. GOTO 52
  707. C
  708. C
  709. 52 CONTINUE
  710. C
  711. C CONVERGENCE CHECK
  712. C
  713. CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
  714. CRINOR = SNNP11 - RTRAC
  715. C
  716. IF ((ABS(CRICIS).LE.PRECIS).AND.(ABS(CRINOR).LE.PRECIS)) THEN
  717. C
  718. C----------LES DEUX CRITERES SONT VERIFIES
  719. C
  720. C
  721. IF (IMPRI.EQ.1) THEN
  722. write(6,*) 'PASSAGE N. 5 - 2'
  723. write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
  724. write(6,*) 'CONVERGENCE A L ITERATION',ITER
  725. ENDIF
  726. C
  727. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  728. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  729. C
  730. C
  731. SIGF(1) = S1NP11
  732. SIGF(2) = S2NP11
  733. SIGF(3) = SNNP11
  734. C
  735. DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
  736. DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
  737. c
  738. if ( abs(sigf(3)) .le. presp ) then
  739. c
  740. c deformation plastique du pas (n+1) moins
  741. c deformation plastique du pas (n)
  742. c
  743. DEFP(3) = EPNNP1 - EPSN1
  744. else
  745. c
  746. c le joint se deforme toujours elastiquement en fermeture
  747. c
  748. DEFP(3) = 0.0D0
  749. c
  750. endif
  751. C
  752. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  753. C
  754. c
  755. c calcul du delta epsilon plastique entre
  756. c le pas actuel et le pas precedent
  757. c
  758. delep1 = defp(1)
  759. delep2 = defp(2)
  760. delepn = defp(3)
  761. c
  762. c produit contracte pour le calcul de la defo. plasti. cumulee
  763. c
  764. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  765. *- cs depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  766. depeq = SQRT(pro1)
  767. VARF(1) = EPSEQ + DEPEQ
  768. VARF(4) = XLAM1 + DELAM1
  769. c
  770. c calcul de la deformation normale due a l'ouverture seule
  771. c
  772. if ( abs(snnp11).le.presp ) then
  773. epoun = delepn
  774. VARF(2) = EPOUN
  775. else
  776. epoun = 0.0d0
  777. VARF(2) = EPOUN
  778. endif
  779. c
  780. c etat du joint stat = 0 cisaillement pur, levres en contact
  781. c stat = 1 ouverture du joint contrainte nulle
  782. c stat = 2 compression du joint elasticite pente EC
  783. c stat = 3 compression du joint plasticite
  784. c
  785. stat = 2.D0
  786. eptmco = epnnp1-EPCOF
  787. if ( abs(snnp11).le.presp) then
  788. stat = 1.D0
  789. if ( abs(eptmco).le.1.D-8 ) stat = 0.D0
  790. endif
  791. F3=CRIPQ(EC,QT,FNE,SNNP11,XLAM3F)
  792. if (ABS(F3).LE.PRECIS) stat = 3.D0
  793. VARF(3) = STAT
  794. c
  795. GOTO 9900
  796. ELSE
  797. C
  798. S1NP10=S1NP11
  799. S2NP10=S2NP11
  800. SNNP10=SNNP11
  801. CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
  802. CRINOR = SNNP10 - RTRAC
  803. ITER = ITER + 1
  804. C
  805. GO TO 2
  806. ENDIF
  807. C
  808. C
  809. C____________________________________________________________________
  810. C
  811. C
  812. C----------CRITERE DE L'EFFORT NORMAL VIOLE MAIS
  813. C----------CRITERE DE CISAILLEMENT VERIFIE
  814. C
  815. 30 CONTINUE
  816. C
  817. C ECOULEMENT SELON LE CRITERE F A L'ITERATION ITER
  818. C
  819. ITER = 1
  820. DELAM1=0.D0
  821. C
  822. 3 CONTINUE
  823. C
  824. IF (ITER.GT.MAXITE) THEN
  825. WRITE(6,*) ' LE CRITERE NORMAL EST VIOLE '
  826. WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
  827. STOP
  828. ENDIF
  829. C
  830. IF (IMPRI.EQ.1) THEN
  831. write(6,*) 'ETIQUETTE 30 - CRITERE DE L EFFORT NORMAL VIOLE'
  832. write(6,*) 'ITERATION',ITER
  833. ENDIF
  834. C
  835. C
  836. C
  837. C CALCUL DU CORRECTEUR PLASTIQUE
  838. C
  839. C
  840. C DERIVEE DU CRITERE F SUIVANT SIGMA
  841. C
  842. DFIDS1 = 0.0D0
  843. DFIDS2 = 0.0D0
  844. DFIDSN = 1.0D0
  845. C
  846. C MATRICE DE HOOK * DGI/DSIGMA
  847. C
  848. HODF1 = G1 * DFIDS1
  849. HODF2 = G2 * DFIDS2
  850. HODFN = EC * DFIDSN
  851. C
  852. C FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  853. C
  854. FIFI = (DFIDS1 * HODF1) + (DFIDS2 * HODF2) + (DFIDSN * HODFN)
  855. C
  856. C RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
  857. C PLASTIQUE DLAM1
  858. C
  859. IF (ABS(FIFI).LE.PRECIS) THEN
  860. WRITE(6,*) ' SYSTEME INDETERMINE (F) : ITERATION ',ITER
  861. STOP
  862. ENDIF
  863. C
  864. DLAM1 = CRINOR / FIFI
  865. DELAM1 = DELAM1 + DLAM1
  866. RTRAC = 0.D0
  867. C
  868. C CALCUL DU NOUVEL ETAT DE CONTRAINTES
  869. C
  870. S1NP11=S1NP10 - (G1 * (DLAM1*DFIDS1) )
  871. S2NP11=S2NP10 - (G2 * (DLAM1*DFIDS2) )
  872. SNNP11=SNNP10 - (EC * (DLAM1*DFIDSN) )
  873. C
  874. C ALLER EN 53 POUR LE TEST DE LA CONVERGENCE
  875. GOTO 53
  876. C
  877. C
  878. C
  879. 53 CONTINUE
  880. C
  881. C CONVERGENCE CHECK
  882. C
  883. CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
  884. CRINOR = SNNP11 - RTRAC
  885. C
  886. IF (ABS(CRINOR).LE.PRECIS) THEN
  887. C
  888. C----------LES DEUX CRITERES SONT VERIFIES
  889. C
  890. C
  891. * IF (IMPRI.EQ.1) THEN
  892. * write(6,*) 'PASSAGE N. 5 - 3'
  893. * write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
  894. * write(6,*) 'CONVERGENCE A L ITERATION',ITER
  895. * ENDIF
  896. C
  897. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  898. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  899. C
  900. C
  901. SIGF(1) = S1NP11
  902. SIGF(2) = S2NP11
  903. SIGF(3) = SNNP11
  904. C
  905. DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
  906. DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
  907. c
  908. if ( abs(sigf(3)) .le. presp ) then
  909. c
  910. c deformation plastique du pas (n+1) moins
  911. c deformation plastique du pas (n)
  912. c
  913. DEFP(3) = EPNNP1 - EPSN1
  914. else
  915. c
  916. c le joint se deforme toujours elastiquement en fermeture
  917. c
  918. DEFP(3) = 0.0D0
  919. c
  920. endif
  921. C
  922. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  923. C
  924. c
  925. c calcul du delta epsilon plastique entre
  926. c le pas actuel et le pas precedent
  927. c
  928. delep1 = defp(1)
  929. delep2 = defp(2)
  930. delepn = defp(3)
  931. c
  932. c produit contracte pour le calcul de la defo. plasti. cumulee
  933. c
  934. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  935. *- cs depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  936. depeq = SQRT(pro1)
  937. VARF(1) = EPSEQ + DEPEQ
  938. VARF(4) = XLAM1 + DELAM1
  939. c
  940. c calcul de la deformation normale due a l'ouverture seule
  941. c
  942. if ( abs(snnp11).le.presp ) then
  943. epoun = delepn
  944. VARF(2) = EPOUN
  945. else
  946. epoun = 0.0d0
  947. VARF(2) = EPOUN
  948. endif
  949. c
  950. c etat du joint stat = 0 cisaillement pur, levres en contact
  951. c stat = 1 ouverture du joint contrainte nulle
  952. c stat = 2 compression du joint elasticite pente EC
  953. c stat = 3 compression du joint plasticite
  954. c
  955. stat = 2.D0
  956. eptmco = epnnp1-EPCOF
  957. if ( abs(snnp11).le.presp) then
  958. stat = 1.D0
  959. if ( abs(eptmco).le.1.D-8 ) stat = 0.D0
  960. endif
  961. F3=CRIPQ(EC,QT,FNE,SNNP11,XLAM3F)
  962. if (ABS(F3).LE.PRECIS) stat = 3.D0
  963. VARF(3) = STAT
  964. c
  965. GOTO 9900
  966. ELSE
  967. C
  968. S1NP10=S1NP11
  969. S2NP10=S2NP11
  970. SNNP10=SNNP11
  971. CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
  972. CRINOR = SNNP10 - RTRAC
  973. ITER = ITER + 1
  974. C
  975. GO TO 3
  976. ENDIF
  977. C
  978. C
  979. C____________________________________________________________________
  980. C
  981. C
  982. C
  983. C----------CRITERES DE CISAILLEMENT VIOLE MAIS
  984. C----------CRITERE DE L'EFFORT NORMAL VERIFIE
  985. C
  986. 40 CONTINUE
  987. C
  988. C ECOULEMENT SELON LE CRITERE G A L'ITERATION ITER
  989. C
  990. ITER = 1
  991. DELAM3=0.D0
  992. C
  993. 4 CONTINUE
  994. C
  995. IF (ITER.GT.MAXITE) THEN
  996. WRITE(6,*) ' LE CRITERE DE CISAILLEMENT EST VIOLE '
  997. WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
  998. STOP
  999. ENDIF
  1000. C
  1001. IF (IMPRI.EQ.1) THEN
  1002. write(6,*) 'ETIQUETTE 40 - CRITERE DE CISAILLEMENT VIOLE'
  1003. ENDIF
  1004. C
  1005. C
  1006. C LE CRITERE DE PLASTICITE EN COMPRESSION EST-IL VIOLE ?
  1007. C
  1008. F3=CRIPQ(EC,QT,FNE,SNNP10,XLAM3)
  1009. C
  1010. IF (F3.LE.PRECIS) THEN
  1011. C-----------------------------------------------------------------------
  1012. C CALCUL DU CORRECTEUR PLASTIQUE POUR LE CRITERE DE CISAILLEMENT
  1013. C-----------------------------------------------------------------------
  1014. C
  1015. C DERIVEE DU CRITERE DE COULOMB G SUIVANT SIGMA
  1016. C
  1017. DENO = SQRT ( (S1NP10**2) + (S2NP10**2) )
  1018. DGIDS1 = S1NP10 / DENO
  1019. DGIDS2 = S2NP10 / DENO
  1020. DZGIDS = TAN(ZMU)
  1021. DGIDSN = TAN(PHI)
  1022. C
  1023. C MATRICE DE HOOK * DGI/DSIGMA
  1024. C
  1025. HODG1 = G1 * DGIDS1
  1026. HODG2 = G2 * DGIDS2
  1027. HODGN = EC * DZGIDS
  1028. C
  1029. C GIGI = DGI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
  1030. C
  1031. GIGI = (DGIDS1 * HODG1) + (DGIDS2 * HODG2) + (DGIDSN * HODGN)
  1032. C
  1033. C
  1034. C RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
  1035. C PLASTIQUE DLAM2
  1036. C
  1037. IF (ABS(GIGI).LE.PRECIS) THEN
  1038. WRITE(LAM36,*) ' SYSTEME INDETERMINE (G) : ITERATION ',ITER
  1039. STOP
  1040. ENDIF
  1041. C
  1042. C Cisaillement
  1043. DLAM2 = CRICIS / GIGI
  1044. C
  1045. C
  1046. C Plasticité en compression
  1047. DLAM3 = 0.
  1048. DELAM3 = DELAM3 + DLAM3
  1049. XLAM3F=XLAM3+DELAM3
  1050. C
  1051. C CALCUL DU NOUVEL ETAT DE CONTRAINTES
  1052. C
  1053. CC SNNP11=SNNP10 - (EC * (DLAM2*DGIDSN) )
  1054. C
  1055. S1NP11=S1NP10 - (G1 * (DLAM2*DGIDS1) )
  1056. S2NP11=S2NP10 - (G2 * (DLAM2*DGIDS2) )
  1057. SNNP11=SNNP10 - (EC * (DLAM2*DZGIDS) )
  1058. C
  1059. ELSE
  1060. C-----------------------------------------------------------------------
  1061. C CALCUL DES CORRECTEURS PLAST. POUR LE CRITERE DE CISAILLEMENT ET
  1062. C DE PLASTICITE EN COMPRESSION
  1063. C-----------------------------------------------------------------------
  1064.  
  1065. C DERIVEE DU CRITERE DE PLASTICITE F3 SUIVANT SIGMA ET LE
  1066. C MULTIPLICATEUR PLASTIQUE XLAM3
  1067. C
  1068. C
  1069. DFIDS1 = 0.0D0
  1070. DFIDS2 = 0.0D0
  1071. DFIDSN = SNNP10/ABS(SNNP10)
  1072. DFIDLA = -H
  1073. C
  1074. C DERIVEE DU CRITERE DE COULOMB G SUIVANT SIGMA
  1075. C
  1076. DENO = SQRT ( (S1NP10**2) + (S2NP10**2) )
  1077. DGIDS1 = S1NP10 / DENO
  1078. DGIDS2 = S2NP10 / DENO
  1079. DZGIDS = TAN(ZMU)
  1080. DGIDSN = TAN(PHI)
  1081. C
  1082. C MATRICE DE HOOK * DFI/DSIGMA
  1083. C
  1084. HODF1 = G1 * DFIDS1
  1085. HODF2 = G2 * DFIDS2
  1086. HODFN = EC * DFIDSN
  1087. C
  1088. C MATRICE DE HOOK * DGI/DSIGMA
  1089. C
  1090. HODG1 = G1 * DGIDS1
  1091. HODG2 = G2 * DGIDS2
  1092. HODGN = EC * DZGIDS
  1093. C
  1094. C
  1095. C FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA - DFI/DLAM2
  1096. C
  1097. FIFI = (DFIDS1*HODF1)+(DFIDS2*HODF2)+(DFIDSN*HODFN)-DFIDLA
  1098. C
  1099. C FIGI = DFI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
  1100. C
  1101. FIGI = (DFIDS1 * HODG1) + (DFIDS2 * HODG2) + (DFIDSN * HODGN)
  1102. C
  1103. C GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
  1104. C
  1105. GIFI = (DGIDS1 * HODF1) + (DGIDS2 * HODF2) + (DGIDSN * HODFN)
  1106. C
  1107. C GIGI = DGI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
  1108. C
  1109. GIGI = (DGIDS1 * HODG1) + (DGIDS2 * HODG2) + (DGIDSN * HODGN)
  1110. C
  1111. C RESOLUTION DU SYSTEME POUR TROUVER LES CORRECTEURS
  1112. C PLASTIQUES DLAM2 ET DLAM3
  1113. C
  1114. DET = (FIFI * GIGI) - (FIGI * GIFI)
  1115. C
  1116. IF (ABS(DET).LE.PRECIS) THEN
  1117. WRITE(6,*) ' SYSTEME INDETERMINE (F G): ITERATION ',ITER
  1118. STOP
  1119. ENDIF
  1120. C
  1121. C Cisaillement
  1122. DLAM2 = (CRICIS * FIFI) - (F3 * GIFI)
  1123. DLAM2 = DLAM2 / DET
  1124. C
  1125. C Plasticité en compression
  1126. DLAM3 = (F3 * GIGI) - (CRICIS * FIGI)
  1127. DLAM3 = DLAM3 / DET
  1128. DELAM3 = DELAM3 + DLAM3
  1129. XLAM3F=XLAM3+DELAM3
  1130. EPCOF=EPCOF+DLAM3*DFIDSN
  1131. C
  1132. C
  1133. C CALCUL DU NOUVEL ETAT DE CONTRAINTES
  1134. C
  1135. S1NP11=S1NP10 - (G1 * ( (DLAM2*DGIDS1) + (DLAM3*DFIDS1)))
  1136. S2NP11=S2NP10 - (G2 * ( (DLAM2*DGIDS2) + (DLAM3*DFIDS2)))
  1137. SNNP11=SNNP10 - (EC * ( (DLAM2*DZGIDS) + (DLAM3*DFIDSN)))
  1138. C
  1139. ENDIF
  1140. C-----------------------------------------------------------------------
  1141. C FIN DE CALCUL DU OU DES CORRECTEURS PLASTIQUES
  1142. C-----------------------------------------------------------------------
  1143. C
  1144. if (SNNP11.ge.0.0d0) SNNP11 = 0.0d0
  1145. C
  1146. C
  1147. 54 CONTINUE
  1148. C
  1149. C CONVERGENCE CHECK
  1150. C
  1151. CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
  1152. F3=CRIPQ(EC,QT,FNE,SNNP11,XLAM3F)
  1153. C
  1154. IF ((ABS(CRICIS).LE.PRECIS.AND.(F3.LE.PRECIS))) THEN
  1155. C
  1156. IF (INDIC.EQ.1) THEN
  1157. C
  1158. C----------1ER TRAJET : ON PARCOURT ALP * DEPSN
  1159. C
  1160. DEFPN = -(epin0(3)-DEFPNC)
  1161. S1INT = S1NP11
  1162. S2INT = S2NP11
  1163. DEFP11 = EP1NP1 - ( DD(1,1) * S1NP11 )
  1164. DEFP21 = EP2NP1 - ( DD(2,2) * S2NP11 )
  1165. GO TO 98
  1166. ELSE
  1167. IF (IFIN.EQ.0) THEN
  1168. C
  1169. C
  1170. C----------LES DEUX CRITERES SONT VERIFIES
  1171. C
  1172. C
  1173. IF (IMPRI.EQ.1) THEN
  1174. write(6,*) 'CONVERGENCE A L ITERATION',ITER
  1175. ENDIF
  1176. C
  1177. C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
  1178. C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
  1179. C
  1180. C
  1181. SIGF(1) = S1NP11
  1182. SIGF(2) = S2NP11
  1183. SIGF(3) = SNNP11
  1184. C
  1185. DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
  1186. DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
  1187. c
  1188. if ( abs(sigf(3)) .le. presp ) then
  1189. c
  1190. c deformation plastique du pas (n+1) moins
  1191. c deformation plastique du pas (n)
  1192. c
  1193. DEFP(3) = EPNNP1 - EPSN1
  1194. else
  1195. c
  1196. c le joint se deforme plastiquement en fermeture
  1197. c
  1198. C DEFP(3) = 0.0D0
  1199. DEFP(3) = EPCOF-EPCO
  1200. c
  1201. endif
  1202. ELSE
  1203. C
  1204. C----------CAS DU PASSAGE DIRECT D'UN ETAT OUVERT A UN ETAT COMPRIME
  1205. C
  1206. DEFP12 = EP1NP1 - (DD(1,1)*S1NP11)
  1207. DEFP22 = EP2NP1 - (DD(2,2)*S2NP11)
  1208. DEFP(1) = DEFP12 - epin0(1)
  1209. DEFP(2) = DEFP22 - epin0(2)
  1210. * DEFP(3) = -epin0(3)
  1211. DEFP(3) = DEFPN+(EPCOF-EPCO)
  1212. SIGF(1) = S1NP11
  1213. SIGF(2) = S2NP11
  1214. SIGF(3) = SNNP11
  1215. IFIN = 0
  1216. ENDIF
  1217. C
  1218. ENDIF
  1219. C
  1220. C CALCUL DES VARIABLES INTERNES : DEPEQ STAT EPOUN
  1221. C
  1222. c
  1223. c calcul du delta epsilon plastique entre
  1224. c le pas actuel et le pas precedent
  1225. c
  1226. delep1 = defp(1)
  1227. delep2 = defp(2)
  1228. delepn = defp(3)
  1229. c
  1230. c produit contracte pour le calcul de la defo. plasti. cumulee
  1231. c
  1232. pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
  1233. *- cs depeq = sqrt( (2.0D0/3.0D0) * pro1 )
  1234. depeq = SQRT(pro1)
  1235.  
  1236. VARF(1) = EPSEQ + DEPEQ
  1237. VARF(4) = XLAM1
  1238. VARF(8) = EPCOF
  1239. VARF(9) = XLAM3F
  1240. c
  1241. c calcul de la deformation normale due a l'ouverture seule
  1242. c
  1243. if ( abs(snnp11).le.presp ) then
  1244. * epoun = delepn
  1245. epoun = delepn-(EPCOF-EPCO)
  1246. VARF(2) = EPOUN
  1247. else
  1248. * epoun = 0.0d0
  1249. epoun = DEFPN
  1250. VARF(2) = EPOUN
  1251. endif
  1252. c
  1253. c etat du joint
  1254. c stat = 0 levres en contact, avec ou sans cisaillement
  1255. c stat = 1 ouverture du joint pente EC
  1256. c stat = 2 compression du joint elastique pente EC
  1257. c stat = 3 compression du joint plastique
  1258. c
  1259. stat = 2.D0
  1260. eptmco = epnnp1-EPCOF
  1261. if ( abs(snnp11).le.presp) then
  1262. stat = 1.D0
  1263. if ( abs(eptmco).le.1.D-8 ) stat = 0.D0
  1264. endif
  1265. F3=CRIPQ(EC,QT,FNE,SNNP11,XLAM3F)
  1266. if (ABS(F3).LE.PRECIS) stat = 3.D0
  1267. VARF(3) = STAT
  1268. c
  1269. GOTO 9900
  1270. ELSE
  1271. C
  1272. S1NP10=S1NP11
  1273. S2NP10=S2NP11
  1274. SNNP10=SNNP11
  1275. * CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
  1276. ITER = ITER + 1
  1277. C
  1278. GO TO 4
  1279. ENDIF
  1280. C
  1281. C
  1282. C____________________________________________________________________
  1283. C
  1284. 9900 CONTINUE
  1285. C ----------CONTRAINTES A LA FIN DU PAS POUR LES 3 DIRECTIONS
  1286. SIGF1(IDIR1)=SIGF(1)
  1287. SIGF1(IDIR2)=SIGF(2)
  1288. SIGF1(IDIRN)=SIGF(3)
  1289. C
  1290. I1=4
  1291. I2=6
  1292. IF (IDIRN.GE.4) THEN
  1293. I1=1
  1294. I2=3
  1295. ENDIF
  1296. DO 350 I=I1,I2
  1297. SIGF1(I)=XZERO
  1298. IF (VARF(3).EQ.1.D0) GOTO 350
  1299. DO 300 J=I1,I2
  1300. SIGF1(I)=SIGF1(I)+TRAID(I,J)*(EPST00(J)+DEPST0(J))
  1301. 300 CONTINUE
  1302. 350 CONTINUE
  1303. C
  1304. C ----------DEF. PLASTIQUES A LA FIN DU PAS POUR LES 3 DIRECTIONS
  1305. DO I=1,3
  1306. VARF(IDEBLA+I)=EPIN0(I)+DEFP(I)
  1307. END DO
  1308. C
  1309. DEFP1(IDIR1)=VARF(IDEBLA+1)
  1310. DEFP1(IDIR2)=VARF(IDEBLA+2)
  1311. DEFP1(IDIRN)=VARF(IDEBLA+3)
  1312. C
  1313. RETURN
  1314. END
  1315.  
  1316.  
  1317.  
  1318.  
  1319.  
  1320.  
  1321.  
  1322.  

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