Télécharger pbefe1.eso

Retour à la liste

Numérotation des lignes :

  1. C PBEFE1 SOURCE BR232186 16/12/06 10:42:33 9242
  2. SUBROUTINE PBEFE1 (XMAT,VAR0,VARF,SIG0,SIGF,DEPST,XCAR,IB,XEL,
  3. & ISTEP)
  4. C
  5. C====&===1=========2=========3=========4=========5=========6=========7==
  6. C COMMENTAIRES : SUBROUTINE PERMETTANT DE METTRE EN OEUVRE LE
  7. C MODELE EFEM
  8. C
  9. C TRAITS :
  10. C
  11. C
  12. C
  13. C AUTEUR : F. RICCARDI - CEA/DEN/DANS/DM2S/SEMT/EMSI
  14. C B. RICHARD - CEA/DEN/DANS/DM2S/SEMT/EMSI
  15. C====&===1=========2=========3=========4=========5=========6=========7==
  16. C
  17. C-------------------------------------------------------------------------------------------
  18. C-----DECLARATION GENERALES
  19. C-------------------------------------------------------------------------------------------
  20. IMPLICIT INTEGER(I-N)
  21. IMPLICIT REAL*8(A-H,O-Z)
  22. C-------------------------------------------------------------------------------------------
  23. C-----DECLARATION SEGMENTS
  24. C-------------------------------------------------------------------------------------------
  25.  
  26. C-----REAL 8
  27. REAL*8 EPSIM(2,2),UNIT(2,2),SIGMAM(2,2),XEL(3,3)
  28. REAL*8 XMAT(*),SIG0(4),SIGF(4),VAR0(*),VARF(*),DEPST(4),XCAR(*)
  29. REAL*8 AD(4),XRHO(4),XECR(4),TD(2),XOUV(2),XNIE(2),XIE(2)
  30. REAL*8 XN(2),XT(2),P1(2),P2(2),P3(2)
  31. REAL*8 BNOD1(1,2),BNOD2(1,2),BNOD3(1,2),BGRAD(3,2)
  32. REAL*8 XANG(3),G(3,2),D(3,3),N(3,2)
  33. REAL*8 EPAUX(2,2),EPSVOIGT(3,1),XG(2)
  34. REAL*8 SIGVOI(3,1),M(2,3),T(2),SIGMA(2,2),VEP(2,2)
  35. REAL*8 VAP(2),DP_MAX(2),VP(1,2),FF(3)
  36. REAL*8 VPT(2),SIGAUX(2,2),T_CN0(1,1)
  37. REAL*8 XNT(1,2),XTT(1,2),T_CN(1),T_CT(1)
  38. REAL*8 EPSD(3),EPSENR_V(3),EPSENR(2,2),SIG(2,2),T_CN1(2)
  39. REAL*8 UN1(2,2),UN2(2,2),EPSI_MAT(2,2),XOUVT(1,2)
  40. REAL*8 T_C(2),T_U(2),EPS1(2,2),EPS2(2,2),EPS4(2,2),EPS7(2,2)
  41. REAL*8 EPS1P(2,2),EPS2P(2,2),EPS4P(2,2),EPS7P(2,2),TEND1(2,2)
  42. REAL*8 TEND2(2,2),TEND4(2,2),TEND7(2,2),TEND(2,2),TENDP(2)
  43. REAL*8 K(2,2),XDG(3,2),F(2,2),RES(2),KINV(2,2),XDGN(3)
  44. REAL*8 XNP(2),XNPT(1,2),DUN(2),EPNL(2,2)
  45. C-------------------------------------------------------------------------------------------
  46. C-----DECLARATION PARAMETRES NUMERIQUES
  47. C-------------------------------------------------------------------------------------------
  48. C
  49. C-----NOMBRE MAXIMUM D ITERATIONS POUR LA CONVERGENCE SUR LES CONTRAINTES PLANES
  50. ITMAX1 = 1000
  51. C-----NOMBRE MAXIMUM D ITERATIONS POUR LA CONVERGENCE SUR LA CONTINUITE
  52. ITMAX2 = 1000
  53. C
  54. C-----CHOIX DE L ECROUISSAGE
  55. C NECR = 1 -- ECROUISSAGE DE BARGELLINI
  56. C NECR = 2 -- ECROUISSAGE LINEAIRE
  57. C NECR = 3 -- ECROUISSAGE EXPONENTIEL
  58. NECR = 3
  59. XKAPPA = 70.0D0
  60. C-----TOLERANCE SUR LA CONDITION DE CONTRAINTES PLANES
  61. TOL1 = 1.0D-6
  62. C
  63. C-----TOLERANCE SUR LA CONDITION DE CONTINUITE
  64. TOL2 = 1.0D-6
  65. C
  66. C-------------------------------------------------------------------------------------------
  67. C-----CONSTRUCTION DE LA MATRICE IDENTITE D ORDRE DEUX
  68. C-------------------------------------------------------------------------------------------
  69. C
  70. DO I=1,2
  71. DO J=1,2
  72. IF (I.EQ.J) THEN
  73. UNIT(I,J) = 1.0D0
  74. ELSE
  75. UNIT(I,J) = 0.0D0
  76. ENDIF
  77. ENDDO
  78. ENDDO
  79. C-------------------------------------------------------------------------------------------
  80. C-----RECUPERATION DES PROPRIETES MATERIAUX
  81. C-------------------------------------------------------------------------------------------
  82. C
  83. C-----MODULE D YOUNG
  84. YOUN = XMAT(1)
  85. C
  86. C-----COEFFICIENT DE POISSON
  87. XNU = XMAT(2)
  88. C
  89. C-----COEFFICIENTS DE LAME
  90. XLAM = XNU*YOUN/((1.0D0+XNU)*(1.0D0-(2.0D0*XNU)))
  91. XMU = YOUN/(2.0D0*(1.0D0+XNU))
  92. C
  93. C-----CALCUL DE LA MATRICE DE HOOKE EN CONTRAINTES PLANES
  94. AM = (YOUN)/(1.0D0-(XNU**2.0D0))
  95. C
  96. D(1,1) = AM
  97. D(1,2) = AM*XNU
  98. D(1,3) = 0.0D0
  99. D(2,1) = AM*XNU
  100. D(2,2) = AM
  101. D(2,3) = 0.0D0
  102. D(3,1) = 0.0D0
  103. D(3,2) = 0.0D0
  104. D(3,3) = AM*(1-XNU)/2.0D0
  105. C
  106. C-----PARAMETRES NON LINEAIRES
  107. ALP1 = XMAT(5)
  108. BET1 = XMAT(6)
  109. AD(1) = XMAT(7)
  110. AD(2) = XMAT(8)
  111. AD(3) = XMAT(9)
  112. AD(4) = XMAT(10)
  113. XFT = XMAT(11)
  114. XECR0 = XMAT(12)
  115. XNIE(1) = XMAT(13)
  116. XNIE(2) = XMAT(14)
  117. XUCR = XMAT(15)
  118. XINDE = INT(XMAT(16))
  119. XIE(1) = XMAT(17)
  120. XIE(2) = XMAT(18)
  121. C
  122. C-------------------------------------------------------------------------------------------
  123. C-----RECUPERATION DES VARIABLES INTERNES
  124. C-------------------------------------------------------------------------------------------
  125. C-----DEFORMATION TOTALE
  126. EPSIM(1,1) = VAR0(5)
  127. EPSIM(2,2) = VAR0(6)
  128. EPSIM(1,2) = VAR0(7)
  129. EPSIM(2,1) = VAR0(7)
  130. C-----DEFORMATION TRANSVERSALE
  131. EPSZZ = VAR0(19)
  132. C-----VECTEUR TRACTION DANS LE REPERE GLOBAL
  133. DO I=1,2
  134. T_U(I) = VAR0(10 + I)
  135. ENDDO
  136. C-----VECTEUR DISCONTINUTE DANS LE REPERE GLOBAL
  137. DO I=1,2
  138. XOUV(I) = VAR0(12 + I)
  139. ENDDO
  140. C-----FLAG ACTIVATION DE LA DISCONTINUITE
  141. XNFLA = VAR0(31)
  142. C-----NORMALE A LA FISSURE ET TRANSPOSEE DE LA NORMALE
  143. XN(1) = VAR0(33)
  144. XN(2) = VAR0(34)
  145. XNT(1,1) = XN(1)
  146. XNT(1,2) = XN(2)
  147. C-----ECROUISSAGE POUR LA CONDITION DE CONTINUITE
  148. Z_U = VAR0(36)
  149. C-----OUVERTURE MAXIMALE DANS LE REPERE DE LA DISCONTINUITE
  150. OUVN_MAX = VAR0(37)
  151. C-----OUVERTURE COURANTE DANS LE REPERE DE LA DISCONTINUITE
  152. OUVN = VAR0(38)
  153. OUVT = VAR0(39)
  154. C-----VECTEUR TRACTION DANS LE REPERE DE LA FISSURE
  155. T_UN = VAR0(43)
  156. T_UT = VAR0(44)
  157. C-----NORME VP
  158. XNORM_VP = VAR0(50)
  159. C-----ENERGIE DISSIPEE
  160. XPSI = VAR0(47)
  161. C
  162. VPT(1) = VAR0(48)
  163. VPT(2) = VAR0(49)
  164. C-------------------------------------------------------------------------------------------
  165. C-----PRECALCUL POUR LA GESTION DE LA NORMALE
  166. C-------------------------------------------------------------------------------------------
  167. C
  168. C-----STOCKAGE DES POINTS
  169. DO J=1,2
  170. P1(J) = XEL(J,1)
  171. P2(J) = XEL(J,2)
  172. P3(J) = XEL(J,3)
  173. ENDDO
  174. C-----CALCUL DE LA SURFACE DE L'ELEMENT FINI
  175. XA0 = 0.5D0*((P2(1)*P3(2))-(P3(1)*P2(2))+(P3(1)*P1(2))-
  176. & (P1(1)*P3(2))+(P1(1)*P2(2))-(P1(2)*P2(1)))
  177. XA = ABS(XA0)
  178. VARF(59) = (SQRT(XA))/2.0D0
  179. C
  180. C-----CALCUL DU CENTRE DE GRAVITE DE L ELEMENT
  181. XG(1) = (P1(1)+P2(1)+P3(1))/3.0D0
  182. XG(2) = (P1(2)+P2(2)+P3(2))/3.0D0
  183. C
  184. C-----CALCUL DES VECTEURS GRADIENT DEFINIS POUR CHAQUE NOEUD DE L'ELEMENT
  185. BNOD1(1,1) = -1.0D0*(P3(2)-P2(2))
  186. BNOD1(1,2) = P3(1) - P2(1)
  187. XNORBN1 = SQRT((BNOD1(1,1)**2) + (BNOD1(1,2)**2))
  188. C
  189. BNOD2(1,1) = -1.0D0*(P1(2)-P3(2))
  190. BNOD2(1,2) = P1(1)-P3(1)
  191. XNORBN2 = SQRT((BNOD2(1,1)**2) + (BNOD2(1,2)**2))
  192. C
  193. BNOD3(1,1) = -1.0D0*(P2(2)-P1(2))
  194. BNOD3(1,2) = P2(1)-P1(1)
  195. XNORBN3 = SQRT((BNOD3(1,1)**2) + (BNOD3(1,2)**2))
  196. C
  197. C-----CONSTRUCTIONS DE LA MATRICE GRADIENT
  198. DO J=1,2
  199. BGRAD(1,J) = BNOD1(1,J)
  200. BGRAD(2,J) = BNOD2(1,J)
  201. BGRAD(3,J) = BNOD3(1,J)
  202. ENDDO
  203. C
  204. C-------------------------------------------------------------------------------------------
  205. C-----ON REALISE LE PREDICTEUR POUR VERIFIER LA CONDITION
  206. C-----D ACTIVATION DE LA DISCONTINUITE
  207. C-------------------------------------------------------------------------------------------
  208. C
  209. C-----ACTUALISATION DE LA DEFORMATION TOTALE
  210. EPSIM(1,1) = EPSIM(1,1) + DEPST(1)
  211. EPSIM(2,2) = EPSIM(2,2) + DEPST(2)
  212. EPSIM(1,2) = EPSIM(1,2) + DEPST(4)/(2.0D0)
  213. EPSIM(2,1) = EPSIM(2,1) + DEPST(4)/(2.0D0)
  214. C
  215. C-----REECRITURE SOUS LA FORME DE VOIGT
  216. EPSVOIGT(1,1) = EPSIM(1,1)
  217. EPSVOIGT(2,1) = EPSIM(2,2)
  218. EPSVOIGT(3,1) = 2.0D0*EPSIM(1,2)
  219. C
  220. C-----CAS OU LA DISCONTINUITE N EST PAS ACTIVE
  221. IF (XNFLA.EQ.0.0D0) THEN
  222. C--------ON CALCUL LES CONTRAINTES ELASTIQUES
  223. CALL PRDFE1(D,EPSVOIGT,3,3,1,SIGVOI)
  224. C
  225. C--------ON REECRIT SOUS LA FORME MATRICIELLE
  226. CALL VECMAT(SIGVOI,SIGMA)
  227. C--------ON INITIALISE LA DEFORMATION TRANSVERSALE
  228. EPSZZ = -XNU/YOUN*(SIGVOI(1,1)+SIGVOI(2,1))
  229. C--------ON INITIALISE LA CONTRAINTES TRANSVERSALE - POUR STABILITE NUMERIQUE
  230. SIG_ZZ = XLAM*((EPSIM(1,1)+EPSIM(2,2))+EPSZZ)+(2*XMU*EPSZZ)
  231. C--------INTEGRATION DU NONLOCAL POUR LE CALCUL DE LA NORMALE
  232. C
  233. IF (ISTEP.EQ.0) THEN
  234. C
  235. EPNL1 = EPSIM(1,1)
  236. EPNL2 = EPSIM(2,2)
  237. EPNL3 = EPSIM(1,2)
  238. C
  239. VARF(51) = EPNL1
  240. VARF(52) = EPNL2
  241. VARF(53) = EPNL3
  242. C
  243. C
  244. ELSEIF (ISTEP.EQ.1) THEN
  245. C
  246. VARF(51) = EPSIM(1,1)
  247. VARF(52) = EPSIM(2,2)
  248. VARF(53) = EPSIM(1,2)
  249. C
  250. DO I=1,50
  251. C
  252. VARF(I)=VAR0(I)
  253. C
  254. ENDDO
  255. C
  256. GOTO 2000
  257. C
  258. ELSEIF (ISTEP.EQ.2) THEN
  259. C
  260. EPNL1 =VAR0(51)
  261. EPNL2 =VAR0(52)
  262. EPNL3 =VAR0(53)
  263. C
  264. VARF(51) = EPNL1
  265. VARF(52) = EPNL2
  266. VARF(53) = EPNL3
  267. C
  268. ENDIF
  269. C
  270. EPNL(1,1) = EPNL1
  271. EPNL(2,2) = EPNL2
  272. EPNL(1,2) = EPNL3
  273. EPNL(2,1) = EPNL3
  274. C
  275. C--------CALCUL VALEUR PROPRES ET VECTEURS PROPRES
  276. CALL DIAFE1(EPNL,VAP,VEP)
  277. C
  278. C--------IDENTIFICATION DE LA CONTRAINTE PRINCIPALE MAJEURE ET DE SON INDICE
  279. EPSI_MAX=DMAX1(VAP(1),VAP(2))
  280. C
  281. DO JJ=1,2
  282. IF (VAP(JJ).EQ.EPSI_MAX) THEN
  283. INDI = JJ
  284. ENDIF
  285. ENDDO
  286. C
  287. C--------IDENTIFICATION DU VECTEUR PROPRE
  288. C
  289. DO J=1,2
  290. DP_MAX(J) = VEP(J,INDI)
  291. ENDDO
  292. C
  293. XNDP = SQRT((DP_MAX(1)**2.0D0)+(DP_MAX(2)**2.0D0))
  294. C
  295. XN(1) = XNIE(1)
  296. XN(2) = XNIE(2)
  297. C
  298. CALL TRSPOS(XN,1,2,XNT)
  299. C
  300. C--------MATRICE DE LA NORMALE - SIMPLICITE DE CALCUL
  301. N(1,1) = XN(1)
  302. N(1,2) = 0.0D0
  303. N(2,1) = 0.0D0
  304. N(2,2) = XN(2)
  305. N(3,1) = XN(2)
  306. N(3,2) = XN(1)
  307. C
  308. C--------TRANSPOSE DE LA MATRICE DE LA NORMALE
  309. DO II=1,2
  310. DO JJ=1,3
  311. M(II,JJ) = N(JJ,II)
  312. ENDDO
  313. ENDDO
  314. C
  315. C--------VECTEUR TRACTION T = SIGMA.N POUR L INITIATION DE LA DISCONTINUITE
  316. CALL PRDFE1(M,SIGVOI,2,3,1,T)
  317. CALL PRDFE1(XNT,T,1,2,1,T_CN0)
  318. C
  319. VP(1,1) = 0.0D0;
  320. VP(1,2) = 0.0D0;
  321. C
  322. IF (XINDE.EQ.1) THEN
  323. C
  324. FP1 = (1.0D0*XNIE(1)*P1(1))+(XNIE(2)*P1(2))-
  325. & ((1.0D0*XNIE(1)*XIE(1))+(XNIE(2)*XIE(2)))
  326. C
  327. FP2 = (1.0D0*XNIE(1)*P2(1))+(XNIE(2)*P2(2))-
  328. & ((1.0D0*XNIE(1)*XIE(1))+(XNIE(2)*XIE(2)))
  329. C
  330. FP3 = (1.0D0*XNIE(1)*P3(1))+(XNIE(2)*P3(2))-
  331. & ((1.0D0*XNIE(1)*XIE(1))+(XNIE(2)*XIE(2)))
  332. C
  333. C 2 - VERIFICATION DES SIGNES
  334. CALL SIGEK(FP1,F1)
  335. CALL SIGEK(FP2,F2)
  336. CALL SIGEK(FP3,F3)
  337. C
  338. FF(1) = F1
  339. FF(2) = F2
  340. FF(3) = F3
  341. C
  342. DO II = 1,3
  343. IF (FF(II).GT.0.0D0) THEN
  344. DO J = 1,2
  345. VP(1,J) = VP(1,J) + BGRAD(II,J)
  346. ENDDO
  347. ENDIF
  348. ENDDO
  349. C
  350. C--------CONDITION D INITIATION
  351. IF ((T_CN0(1,1).GT.XFT)) THEN
  352. T_UN = XFT
  353. T_UT = 0
  354. OUVN = 1.0D-7
  355. VAR0(38) = OUVN
  356. OUVT = 1.0D-12
  357. VAR0(39) = OUVT
  358. Z_U = 0.0D0
  359. XNFLA = 1.0D0
  360. ELSE
  361. XNFLA = 0.0D0
  362. ENDIF
  363. ENDIF
  364. C
  365. CALL TRSPOS(VP,2,1,VPT)
  366. XNORM_VP = SQRT((VPT(1)**2)+(VPT(2)**2))
  367. C
  368. ENDIF
  369. C--------FIN DE LA PREDICTION ELASTIQUE NORMALE
  370. C
  371. C--------DEFINITION DE LA MATRICE G
  372. G(1,1) = -1.D0*VPT(1)/(2.0D0*XA)
  373. G(1,2) = 0.0D0/(2.0D0*XA)
  374. G(2,1) = 0.0D0/(2.0D0*Xa)
  375. G(2,2) = -1.D0*VPT(2)/(2.0D0*XA)
  376. G(3,1) = -1.D0*VPT(2)/(2.0D0*XA)
  377. G(3,2) = -1.D0*VPT(1)/(2.0D0*XA)
  378. C
  379. C--------CALCUL DE DERIVE T_CN
  380. CALL PRDFE1(D,G,3,3,2,XDG)
  381. CALL PRDFE1(XDG,XN,3,2,1,XDGN)
  382. CALL VECMAT(XDGN,F)
  383. C
  384. C--------CALCUL DU VECTEUR TANGENT A LA FISSURE
  385. XT(1) = -XN(2)
  386. XT(2) = XN(1)
  387. C
  388. C CALL TRSPOS(XT,1,2,XTT)
  389.  
  390. XTT(1,1) = XT(1)
  391. XTT(1,2) = XT(2)
  392.  
  393. C print*,'XTT=',XTT
  394. C print*,'XT=',XT
  395.  
  396. C stop
  397.  
  398. C-------------------------------------------------------------------------------------------
  399. C-----ON REALISE DE NOUVEAU UNE ETAPE DE PREDICTION ELASTIQUE EN PRENANT EN COMPTE
  400. C-----UNE EVENTUELLE DISCONTINUITE
  401. C-------------------------------------------------------------------------------------------
  402. C
  403. C-----CALCUL DU TERME G X [U]
  404. CALL PRDFE1(G,XOUV,3,2,1,EPSD)
  405. C
  406. C-----CALCUL DU TERME EPSILON + G X [U]
  407. DO II=1,3
  408. EPSENR_V(II) = EPSVOIGT(II,1) + EPSD(II)
  409. ENDDO
  410. C
  411. C-----REECRITURE EN MATRCIEL
  412. CALL VECDEF(EPSENR_V,EPSENR)
  413. C
  414. C-----CALCUL DES CONTRAINTE DANS LE CONTINUUM
  415. DO II=1,2
  416. DO JJ=1,2
  417. SIG(II,JJ) = XLAM*(EPSENR(1,1)+EPSENR(2,2)+
  418. & EPSZZ)*UNIT(II,JJ)+
  419. & 2.0D0*XMU*EPSENR(II,JJ)
  420. ENDDO
  421. ENDDO
  422. C
  423. C-----CALCUL DU VECTEUR TRACTION
  424. CALL PRDFE1(SIG,XN,2,2,1,T_CN1)
  425. CALL PRDFE1(XNT,T_CN1,1,2,1,T_CN)
  426. CALL PRDFE1(XTT,T_CN1,1,2,1,T_CT)
  427.  
  428. C CALL PRODMA(SIG,XN,2,2,1,T_CN1)
  429. C CALL PRODMA(XNT,T_CN1,1,2,1,T_CN)
  430. C CALL PRODMA(XTT,T_CN1,1,2,1,T_CT)
  431.  
  432. C print*,'XNT=',XNT
  433. C print*,'XTT=',XTT
  434.  
  435. C print*,'T_CN=',T_CN(1)
  436. C print*,'T_CT=',T_CT(1)
  437. C STOP
  438. C
  439. C-------------------------------------------------------------------------------------------
  440. C-----BOUCLE SUR LA CONDITION DE CONTRAINTE PLANES - START
  441. C-------------------------------------------------------------------------------------------
  442. C
  443. DO II=1,ITMAX1
  444. C
  445. C--------SI PAS DE DISCONTINUITE ON RENVOIE LA SOLUTION ELASTIQUE
  446. IF (XNFLA.EQ.1.0D0) THEN
  447. C-----------INITIALISATION DU RESIDU SUR LA CONTINUITE
  448. RES(1) = T_CN(1) - T_UN
  449. RES(2) = T_CT(1) - T_UT
  450. C
  451. C-------------------------------------------------------------------------------------------
  452. C-------------BOUCLE SUR LA CONDITION DE CONTINUITE - START
  453. C-------------------------------------------------------------------------------------------
  454. DO III=1,ITMAX2
  455. C
  456. C----------------CALCUL DE LA CORRECTION DE DISCONTINUITE
  457. C
  458. DO JJ = 1,2
  459. DO KK = 1,2
  460. K(JJ,KK) = 1.0D0*YOUN*UNIT(JJ,KK)/
  461. & DSQRT(2.0D0*XA)
  462. ENDDO
  463. ENDDO
  464. K(1,2) = 1.0D0*YOUN / 10.0D0
  465. K(2,1) = 1.0D0*YOUN / 10.0D0
  466. CALL HINV22(K,KINV)
  467. CALL PRDFE1(KINV,RES,2,2,1,DUN)
  468. C
  469. DUNN = DUN(1)
  470. DUNT = DUN(2)
  471. C
  472. C----------------ACTUALISATION DE LA DISCONTINUITE (RELAXATION NUMERIQUE)
  473. OUVN = OUVN + (1.0D-3*DUNN)
  474. OUVT = OUVT + (1.0D-3*DUNT)
  475. C
  476. C----------------CHARGE OU DECHARGE (POUR LA CONDITION DE CONTINUITE)
  477. PHI_U = 1.0D0
  478. C
  479. C
  480. C----------------CONDITION D UNILATERALITE
  481. IF (OUVN.LT.0.0D0) THEN
  482. OUVN = 1.0D-12
  483. GOTO 9
  484. ENDIF
  485. C
  486. C----------------ACTUALISATION DE L OUVERTURE NORMALE MAXIMALE
  487. OUVN_MAX = DMAX1(OUVN,VAR0(37))
  488. C
  489. C----------------SWITCH VERS LA LOI DISCRETE CONSIDEREE
  490. C
  491. C 3. ECROUISSAGE EXPONENTIEL
  492. IF (NECR.EQ.3) THEN
  493. C
  494. IF (PHI_U.GT.0.0D0) THEN
  495. C
  496. IF ((II.EQ.1).AND.(III.EQ.1)) THEN
  497. VMOUV0 = DSQRT((VAR0(38)**2.0D0)+(VAR0(39)**2.0D0))
  498. VMOUV = VMOUV0
  499. ELSE
  500. C
  501. VMOUV = DSQRT((OUVN**2.0D0)+(OUVT**2.0D0))
  502. ENDIF
  503. T_UN = XFT*(OUVN/VMOUV)*EXP(-1.0D0*VMOUV/XKAPPA*XFT)
  504. T_UT = XFT*(OUVT/VMOUV)*EXP(-1.0D0*VMOUV/XKAPPA*XFT)
  505.  
  506. Z_U = T_UN
  507. ELSE
  508. C
  509. T_UN = (Z_U/OUVN_MAX)*OUVN
  510. C
  511. ENDIF
  512. ENDIF
  513. C
  514. C----------------CALCUL DU VECTEUR DISCONTINUITE
  515. DO JJ=1,2
  516. XOUV(JJ) = (OUVN*XN(JJ))+(OUVT*XT(JJ))
  517. ENDDO
  518. C
  519. C----------------CALCUL DU TERME G X [U]
  520. CALL PRDFE1(G,XOUV,3,2,1,EPSD)
  521. C
  522. C----------------CALCUL DU TERME EPSILON + G X [U]
  523. DO JJ=1,3
  524. EPSENR_V(JJ) = EPSVOIGT(JJ,1) + EPSD(JJ)
  525. ENDDO
  526. C
  527. C----------------REECRITURE EN MATRCIEL
  528. CALL VECDEF(EPSENR_V,EPSENR)
  529. C
  530. C----------------CALCUL DES CONTRAINTE DANS LE CONTINUUM
  531. DO JJ=1,2
  532. DO KK=1,2
  533. SIG(JJ,KK) = XLAM*((EPSENR(1,1)+
  534. & EPSENR(2,2))+ EPSZZ)*UNIT(JJ,KK)+
  535. & 2.0D0*XMU*EPSENR(JJ,KK)
  536. ENDDO
  537. ENDDO
  538. C
  539. C----------------CALCUL DU VECTEUR CONTRAINTE
  540. CALL PRDFE1(SIG,XN,2,2,1,T_CN1)
  541. CALL PRDFE1(XNT,T_CN1,1,2,1,T_CN)
  542. CALL PRDFE1(XTT,T_CN1,1,2,1,T_CT)
  543.  
  544. C print*,'T_CN=',T_CN(1)
  545. C print*,'T_CT=',T_CT(1)
  546. C STOP
  547. C
  548. C----------------CALCUL DU RESIDU ACTUALISE
  549. RES(1) = T_CN(1) - T_UN
  550. RES(2) = T_CT(1) - T_UT
  551. C
  552. C
  553. XNORMRES = DSQRT((RES(1)*RES(1))+(RES(2)*RES(2)))
  554. C
  555. C----------------CALCUL DE L ENERGIE DISSIPEE
  556. XPSI = T_UN*OUVN
  557. C
  558. C----------------TEST DE CONVERGENCE
  559. IF ((XNORMRES).LT.TOL2) THEN
  560. C
  561. C-------------------CALCUL DE L OUVERTURE DANS LE REPERE GLOBAL
  562. DO JJ=1,2
  563. XOUV(JJ) = (OUVN*XN(JJ))+(OUVT*XT(JJ))
  564. ENDDO
  565. C
  566. C-------------------CALCUL DES VECTEURS TRACTIONS DANS LE CONTINUUM ET
  567. C-------------------SUR LA DISCONTINUITE
  568. DO JJ=1,2
  569. T_C(JJ) = (T_CN(1)*XN(JJ))+(T_CT(1)*XT(JJ))
  570. T_U(JJ) = (T_UN*XN(JJ))+(T_UT*XT(JJ))
  571. ENDDO
  572. C-------------------BALISE DE SORTIE
  573. GOTO 9
  574. C
  575. ENDIF
  576. C
  577. ENDDO
  578. C
  579. C-------------------------------------------------------------------------------------------
  580. C-------------BOUCLE SUR LA CONDITION DE CONTINUITE - END
  581. C-------------------------------------------------------------------------------------------
  582. ENDIF
  583. C
  584. 9 CONTINUE
  585. C
  586. C
  587. C-----------CALCUL DES CONTRAINTES DE CAUCHY A CONVERGENCE
  588. DO JJ=1,2
  589. XOUV(JJ) = (OUVN*XN(JJ))+(OUVT*XT(JJ))
  590. ENDDO
  591. C
  592. CALL PRDFE1(G,XOUV,3,2,1,EPSD)
  593. C
  594. DO JJ=1,3
  595. EPSENR_V(JJ) = EPSVOIGT(JJ,1) + EPSD(JJ)
  596. ENDDO
  597. C
  598. CALL VECDEF(EPSENR_V,EPSENR)
  599. C
  600. DO JJ=1,2
  601. DO KK=1,2
  602. SIG(JJ,KK) = XLAM*((EPSENR(1,1)+
  603. & EPSENR(2,2))+ EPSZZ)*UNIT(JJ,KK)+
  604. & 2.0D0*XMU*EPSENR(JJ,KK)
  605. ENDDO
  606. ENDDO
  607. C
  608. C-----------CALCUL CONTRAINTES HORS PLAN
  609. SIG_ZZ = XLAM*((EPSENR(1,1)+EPSENR(2,2))+
  610. & EPSZZ)+(2*XMU*EPSZZ)
  611. C
  612. C-----------TEST DE CONVERGENCE SUR LA CONDITION DE CONTRAINTES PLANES
  613. IF (ABS(SIG_ZZ).LT.TOL1) THEN
  614. C
  615. GOTO 99
  616. C
  617. ELSE
  618. C
  619. EPSZZ = EPSZZ - (SIG_ZZ/(XLAM+(2*XMU)))
  620. ENDIF
  621. C
  622. C
  623. ENDDO
  624. C-------------------------------------------------------------------------------------------
  625. C-----BOUCLE SUR LA CONDITION DE CONTRAINTE PLANES - END
  626. C-------------------------------------------------------------------------------------------
  627. C
  628. 99 CONTINUE
  629. C
  630. C-------------------------------------------------------------------------------------------
  631. C-----CONTRAINTES DE CAUCHY EN SORTIE POUR CALCUL DES EFFORTS INTERNES
  632. C-------------------------------------------------------------------------------------------
  633. IF (XNFLA.EQ.0.0D0) THEN
  634. CALL VECMAT(SIGVOI,SIGMAM)
  635. ELSE
  636. DO JJ=1,2
  637. DO KK=1,2
  638. SIGMAM(JJ,KK) = SIG(JJ,KK)
  639. ENDDO
  640. ENDDO
  641. ENDIF
  642. C
  643. SIGF(1) = SIGMAM(1,1)
  644. SIGF(2) = SIGMAM(2,2)
  645. SIGF(3) = 0.0D0
  646. SIGF(4) = SIGMAM(1,2)
  647. C
  648. C
  649. C-------------------------------------------------------------------------------------------
  650. C-----VARIABLES INTERNES FINALES
  651. C-------------------------------------------------------------------------------------------
  652. C-----VECTEUR TRACTION CALCULES PAR LA LOI DISCRETE DANS LE REPERE GLOBAL
  653. DO I=1,2
  654. VARF(10+I) = T_U(I)
  655. ENDDO
  656. C
  657. C-----VECTEUR DISCONTINUITE DANS LE REPERE GLOBAL
  658. DO I=1,2
  659. VARF(12+I) = XOUV(I)
  660. ENDDO
  661. C
  662. C-----DEFORMATIONS TOTALES
  663. VARF(5) = EPSIM(1,1)
  664. VARF(6) = EPSIM(2,2)
  665. VARF(7) = EPSIM(1,2)
  666. C-----DEFORMATION TRANSVERSALE
  667. VARF(19) = EPSZZ
  668. C-----FLAG ACTIVATION DISCONTINUITE
  669. VARF(31) = XNFLA
  670. C-----VECTEUR DE LA NORMALE DANS LE REPERE GLOBAL
  671. VARF(33) = XN(1)
  672. VARF(34) = XN(2)
  673. C-----ECROUISSAGE LIE A LA LOI DISCRETE
  674. VARF(36) = Z_U
  675. C-----OUVERTURE NORMALE MAXIMALE
  676. VARF(37) = OUVN_MAX
  677. C-----OUVERTURE NORMALE COURANTE
  678. VARF(38) = OUVN
  679. VARF(39) = OUVT
  680. C-----VECTEUR TRACTION DISCRET
  681. VARF(43) = T_UN
  682. VARF(44) = T_UT
  683. C-----NORME DE VP
  684. VARF(50) = XNORM_VP
  685. IF (XNFLA.EQ.1.0D0) THEN
  686. VARF(45) = XN(1)
  687. VARF(46) = XN(2)
  688. ENDIF
  689. C
  690. C-----ENERGIE DISSIPEE
  691. C
  692. VARF(47) = XPSI
  693. VARF(48) = VPT(1)
  694. VARF(49) = VPT(2)
  695. C-----CONTRAINTES
  696. VARF(54) = SIGF(1)
  697. VARF(55) = SIGF(2)
  698. VARF(56) = SIGF(3)
  699. VARF(57) = SIGF(4)
  700. 2000 CONTINUE
  701. C
  702. RETURN
  703. END
  704.  
  705.  
  706.  

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