Télécharger psury.eso

Retour à la liste

Numérotation des lignes :

  1. C PSURY SOURCE CB215821 16/04/21 21:18:08 8920
  2. SUBROUTINE PSURY(ENDO,NENDO,NVARI,NSTRS,MFR1,DEPST,XMAT,VAR0,
  3. 1 RAPP,NRAPP,
  4. 1 SIG0,SIGF,VARF,NMATT,DEFP,KERRE)
  5. *
  6. * Entrées
  7. *
  8. * ENDO: courbe de debut d'endommagement
  9. * NENDO: nombre de points sur la courbe ENDO
  10. * RAPP: courbe d'evolution de l'endommagement en fonction de la
  11. * pseudo porosite
  12. * NRAPP: nombre de points de la courbe RAPP
  13. * NSTRS: nombre de composantes des deformations
  14. * MFR1: numero de la formulation
  15. * DEPST: increment de deformation totale
  16. * XMAT: donnees materiau
  17. * SIGF: contraintes plastiquement admissibles issues de ECOINC
  18. * VAR0: variables internes au debut du pas de temps
  19. * VARF: variables internes issues de ECOINC
  20. * SIG0: contraintes au debut du pas de temps
  21. * NMATT: nombre de composantes matériaux
  22. * DEFP: incrément de la déformation plastique
  23. *
  24. * Sorties
  25. *
  26. * SIGF: contraintes finales ( tiennent compte de l'endommagement)
  27. * VARF: variables internes finales
  28. *
  29. *
  30. IMPLICIT INTEGER(I-N)
  31. IMPLICIT REAL*8(A-H,O-Z)
  32. *
  33. -INC CCOPTIO
  34. -INC CCREEL
  35. DIMENSION SIGF(*),XMAT(*),VARF(*),DEPST(*),ENDO(*)
  36. DIMENSION SIG0(*),VAR0(*),DEFP(*),RAPP(*)
  37. DIMENSION RSIGF(6),RSIG0(6)
  38. *
  39. *
  40. *====================================================================
  41. * Adaptation de l'option de calcul vers le 3D massif de SIGF a RSIGF
  42. *====================================================================
  43. *
  44. IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN
  45. *
  46. *---> 1 formulation massive
  47. *---> 2 formulation quasi incompressible
  48. *---> MASSIF 3D
  49. *
  50. IF (NSTRS .EQ. 6) THEN
  51. DO 110 I=1,NSTRS
  52. RSIGF(I)=SIGF(I)
  53. RSIG0(I)=SIG0(I)
  54. 110 CONTINUE
  55. ELSE IF ( NSTRS .EQ. 4 .AND. ((IFOUR .EQ. 0)
  56. & .OR.(IFOUR.EQ.-1).OR.(IFOUR.EQ.-2))) THEN
  57. *
  58. *---> Calcul en mode deformations planes ou axisymetrique
  59. *---> Calcul en mode contraintes planes
  60. *
  61. DO 115 I=1,NSTRS
  62. RSIGF(I)=SIGF(I)
  63. RSIG0(I)=SIG0(I)
  64. 115 CONTINUE
  65. RSIGF(5)=0.D0
  66. RSIG0(5)=0.D0
  67. RSIGF(6)=0.D0
  68. RSIG0(6)=0.D0
  69. ENDIF
  70. ELSE
  71. KERRE = 99
  72. RETURN
  73. ENDIF
  74. *
  75. *========================================
  76. * Calcul de la densité du materiau
  77. *========================================
  78. *
  79. * Si on a oublie de rentrer la densite initiale
  80. *
  81. * IF (IFOUR.EQ.-2) THEN
  82. * DENS0=XMAT(NMATT)
  83. * ELSE
  84. * DENS0=XMAT(3)
  85. * ENDIF
  86. *
  87. * DENS0=XMAT(3)
  88. DENS0=XMAT(NMATT)
  89. IF (DENS0.LE.1.D-10) THEN
  90. KERRE=71
  91. RETURN
  92. ENDIF
  93. *
  94. * Initialisation de la densite en variable interne
  95. *
  96. IF (VAR0(2).LE.1.D-10) THEN
  97. VAR0(2)=DENS0
  98. ENDIF
  99. IF (IFOUR.EQ.-2) THEN
  100. treps0=(1.D0-2.D0*XMAT(2))/(1.D0-XMAT(2))
  101. treps0=treps0*(DEPST(1)+DEPST(2)-DEFP(1)-DEFP(2))
  102. ELSE
  103. treps0=DEPST(1)+DEPST(2)+DEPST(3)
  104. ENDIF
  105. rho=DENS0/((DENS0/VAR0(2))+treps0)
  106. *
  107. *
  108. *========================================
  109. * A t'on déja endommagé ?
  110. *========================================
  111. *
  112. IF (ABS(VAR0(3)).GT.1.D-10) THEN
  113. *
  114. * On a déja endommagé
  115. *----------------------------------------
  116. *
  117. * Calcul de la pseudo porosité du temps n+1
  118. *
  119. alpha=(VAR0(3)-rho)/VAR0(3)
  120. *
  121. * Calcul de la variable d'endommagement
  122. *
  123. D_end1=0.D0
  124. DO 12 I=1,NRAPP-1
  125. DD1=RAPP(2*(I-1)+1)
  126. DD2=RAPP(2*I+1)
  127. AL1=RAPP(2*I)
  128. AL2=RAPP(2*(I+1))
  129. IF ((alpha.GE.AL1).AND.(alpha.LT.AL2)) THEN
  130. D_end1=(DD2-DD1)/(AL2-AL1)*(alpha-AL1)+DD1
  131. ENDIF
  132. 12 CONTINUE
  133. IF ((D_end1.LE.0.D0).OR.(alpha.LE.0.D0)) D_end1=0.D0
  134. IF ((D_end1.GE.1.D0).OR.(alpha.GE.1.D0)) D_end1=1.D0
  135. D_max=VAR0(5)
  136. *
  137. * Calcul de l'ancienne pseudo porosite ( temps n)
  138. *
  139. alp_old=VAR0(4)
  140. *
  141. * Calcul de la fonction d'endommagement g0
  142. *
  143. IF ((alpha.GT.0.D0).AND.(alpha.GT.alp_old)) THEN
  144. g0=1.D0-D_end1
  145. ELSE IF ((alpha.GT.0.D0).AND.(alpha.LE.alp_old)) THEN
  146. g0=1.D0-D_max
  147. ELSE
  148. g0=1.D0
  149. ENDIF
  150. g0=MAX(g0,0.D0)
  151. *
  152. * Calcul des contraintes vérifiant l'endommagement
  153. *
  154. DO 10 I=1,6
  155. RSIGF(I)=RSIGF(I)*g0
  156. 10 CONTINUE
  157. *
  158. * Mise a jour des variables internes
  159. *
  160. *---> Densite
  161. VARF(2)=rho
  162. *---> Densite de debut de fracture
  163. VARF(3)=VAR0(3)
  164. *---> Pseudo porosite
  165. VARF(4)=alpha
  166. *---> Pseudo porosite maximale
  167. VARF(5)=MAX(VAR0(5),D_end1)
  168. *---> Coefficient d'endommagement
  169. VARF(6)=D_end1
  170. *---> Fonction d'endommagement
  171. VARF(7)=g0
  172. *
  173. ELSE
  174. *
  175. * On a pas encore endommagé
  176. *-------------------------------------------
  177. *
  178. * Calcul de P=1/3*trace(RSIGF)
  179. *
  180. P0=(RSIGF(1)+RSIGF(2)+RSIGF(3))/3.D0
  181. *
  182. * Calcul de la contrainte équivalente
  183. *
  184. Y1=(RSIGF(1)*RSIGF(1))+(RSIGF(2)*RSIGF(2))
  185. Y1=Y1+(RSIGF(3)*RSIGF(3))
  186. Y1=Y1-(RSIGF(1)*RSIGF(2))-(RSIGF(2)*RSIGF(3))
  187. Y1=Y1-(RSIGF(3)*RSIGF(1))
  188. Y2=(RSIGF(4)*RSIGF(4))+(RSIGF(5)*RSIGF(5))
  189. Y2=Y2+(RSIGF(6)*RSIGF(6))
  190. Y0=Y1+(3.D0*Y2)
  191. Y0=(Y0)**(0.5D0)
  192. *
  193. * Rapport de triaxialite
  194. *
  195. Rapp0=P0/Y0
  196. *
  197. * Rapport de triaxialite de debut d'endommagement Rapp_th
  198. *
  199. EPSE=VARF(1)
  200. IF (EPSE.LT.ENDO(2)) THEN
  201. Rapp_th=1.D30
  202. pente1=1.D30
  203. Rmax=Rapp_th
  204. epsmin=ENDO(2)
  205. ELSE IF (EPSE.GE.ENDO(2*NENDO)) THEN
  206. Rapp_th=ENDO(2*NENDO-1)
  207. pente1=0.D0
  208. Rmax=Rapp_th
  209. epsmin=ENDO(2*NENDO)
  210. ELSE
  211. DO 100 I=1,(NENDO-1)
  212. IF ((EPSE.GE.ENDO(2*I)).AND.(EPSE.LT.ENDO(2*I+2))) THEN
  213. IND0=I
  214. ENDIF
  215. 100 CONTINUE
  216. epsmin=ENDO(2*IND0)
  217. epsmax=ENDO(2*IND0+2)
  218. dif0=ABS(epsmax-epsmin)
  219. Rmax=ENDO(2*IND0-1)
  220. Rmin=ENDO(2*IND0+1)
  221. IF (dif0.GT.1.D-10) THEN
  222. pente1=(Rmax-Rmin)/(epsmin-epsmax)
  223. Rapp_th=pente1*(EPSE-epsmin)+Rmax
  224. ELSE
  225. pente1=1.D30
  226. Rapp_th=Rmin
  227. ENDIF
  228. ENDIF
  229. *
  230. * Comparaison au rapport de triaxialite d'endommagement
  231. *
  232. IF (Rapp0.GT.Rapp_th) THEN
  233. *
  234. * On endommage
  235. *
  236. * Calcul de la densite de debut d'endommagement : rho_f
  237. *---------------------------------------------------
  238. *
  239. *---> Rapport de triaxialite au debut du pas
  240. *
  241. P_old=(RSIG0(1)+RSIG0(2)+RSIG0(3))/3.D0
  242. Y_old1=(RSIG0(1)*RSIG0(1))+(RSIG0(2)*RSIG0(2))
  243. Y_old1=Y_old1+(RSIG0(3)*RSIG0(3))
  244. Y_old1=Y_old1-(RSIG0(1)*RSIG0(2))-(RSIG0(2)*RSIG0(3))
  245. Y_old1=Y_old1-(RSIG0(3)*RSIG0(1))
  246. Y_old2=(RSIG0(4)*RSIG0(4))+(RSIG0(5)*RSIG0(5))
  247. Y_old2=Y_old2+(RSIG0(6)*RSIG0(6))
  248. Y_old0=Y_old1+(3.D0*Y2)
  249. Y_old=(Y_old0)**(0.5D0)
  250. R_old=P_old/Y_old
  251. *
  252. *---> Courbe d'evolution du rapport P/Y entre le debut et la fin du pas
  253. * On l'approxime par une droite de pente pente0
  254. *
  255. dvar1=VARF(1)-VAR0(1)
  256. * IF (dvar1.GT.1.D-10) THEN
  257. IF (dvar1.GT.abs(varf(1)*xzprec)) THEN
  258. pente0=(Rapp0-R_old)/dvar1
  259. ELSE
  260. pente0=abs(varf(1))/XZPREC
  261. ENDIF
  262. *
  263. *---> intersection avec la courbe de debut d'endommagement
  264. * point (R_int,eps_int)
  265. *
  266. 150 IF (pente1.GE.1.D30) THEN
  267. eps_int=epsmin
  268. ELSE IF (pente0.GE.1.D30) THEN
  269. eps_int=VAR0(1)
  270. ELSE
  271. eps_int=Rapp0-Rmax+(pente1*epsmin)-(pente0*VARF(1))
  272. eps_int=eps_int/(pente1-pente0)
  273. ENDIF
  274. IF ((eps_int.LT.epsmin).AND.(IND0.GE.2)) THEN
  275. IND0=IND0-1
  276. epsmin=ENDO(2*IND0)
  277. epsmax=ENDO(2*IND0+2)
  278. dif0=ABS(epsmin-epsmax)
  279. Rmax=ENDO(2*IND0-1)
  280. Rmin=ENDO(2*IND0+1)
  281. IF (dif0.GT.1.D-10) THEN
  282. pente1=(Rmax-Rmin)/(epsmin-epsmax)
  283. ELSE
  284. pente1=1.D30
  285. ENDIF
  286. GOTO 150
  287. ELSE IF ((eps_int.LT.epsmin).AND.(IND0.EQ.1)) THEN
  288. pente1=1.D30
  289. eps_int=ENDO(2)
  290. epsmin=ENDO(2)
  291. epsmax=ENDO(2)
  292. Rmin=ENDO(1)
  293. Rmax=1.D30
  294. ENDIF
  295. IF (pente1.GE.1.D30) THEN
  296. R_int=pente0*(eps_int-VARF(1))+Rapp0
  297. ELSE
  298. R_int=pente1*(eps_int-epsmin)+Rmax
  299. ENDIF
  300. *
  301. *---> Calcul de rho_f en supposant le module d'ecrouissage
  302. * constant entre le debut et la fin du pas
  303. *
  304. IF (dvar1.GT.1.D-10) THEN
  305. H0=(Y0-Y_old)/dvar1
  306. Y_int=H0*(eps_int-VAR0(1))+Y_old
  307. ELSE
  308. *
  309. *---> Endommagement dans un cas élastique
  310. *
  311. IF ((ABS(Rapp0-R_old)).GT.1.D-20) THEN
  312. alfa0=(R_int-R_old)/(Rapp0-R_old)
  313. Y_int=(alfa0*Y0)+((1.D0-alfa0)*Y_old)
  314. ELSE
  315. Y_int=Y_old
  316. ENDIF
  317. ENDIF
  318. IF (treps0.GT.1.D-10) THEN
  319. XK0=XMAT(1)/(3.D0*(1.D0-2.D0*XMAT(2)))
  320. P_int=R_int*Y_int
  321. tr_eps1=(P_int-P_old)/XK0
  322. ELSE
  323. tr_eps1=0.D0
  324. ENDIF
  325. rho_f=DENS0/((DENS0/VAR0(2))+tr_eps1)
  326. IF ((rho_f.GT.1.D10).OR.(rho_f.LT.0.D0)) THEN
  327. rho_f=rho
  328. ENDIF
  329. *
  330. * Calcul de la pseudo-porosite
  331. *
  332. alpha=(rho_f-rho)/rho_f
  333. *
  334. * Calcul de la variable d'endommagement
  335. *
  336. D_end1=0.D0
  337. DO 13 I=1,NRAPP-1
  338. DD1=RAPP(2*(I-1)+1)
  339. DD2=RAPP(2*I+1)
  340. AL1=RAPP(2*I)
  341. AL2=RAPP(2*(I+1))
  342. IF ((alpha.GE.AL1).AND.(alpha.LT.AL2)) THEN
  343. D_end1=(DD2-DD1)/(AL2-AL1)*(alpha-AL1)+DD1
  344. ENDIF
  345. 13 CONTINUE
  346. IF ((D_end1.LE.0.D0).OR.(alpha.LE.0.D0)) D_end1=0.D0
  347. IF ((D_end1.GE.1.D0).OR.(alpha.GE.1.D0)) D_end1=1.D0
  348. D_max=VAR0(5)
  349. *
  350. * Calcul de l'ancienne pseudo porosite ( temps n)
  351. *
  352. alp_old=VAR0(4)
  353. *
  354. * Calcul de la fonction d'endommagement g0
  355. *
  356. IF ((alpha.GT.0.D0).AND.(alpha.GT.alp_old)) THEN
  357. g0=1.D0-D_end1
  358. ELSE IF ((alpha.GT.0.D0).AND.(alpha.LE.alp_old)) THEN
  359. g0=1.D0-D_max
  360. ELSE
  361. g0=1.D0
  362. ENDIF
  363. g0=MAX(g0,0.D0)
  364. *
  365. * Calcul des contraintes vérifiant l'endommagement
  366. *
  367. DO 15 I=1,6
  368. RSIGF(I)=RSIGF(I)*g0
  369. 15 CONTINUE
  370. *
  371. * Mise a jour des variables internes
  372. *
  373. VARF(2)=rho
  374. VARF(3)=rho_f
  375. VARF(4)=alpha
  376. VARF(5)=MAX(VAR0(5),D_end1)
  377. VARF(6)=D_end1
  378. VARF(7)=g0
  379. IF (D_end1.GT.0.9) THEN
  380. write(*,*) ' Debut endommagement'
  381. write(*,*) 'P_int,P_old=',P_int,P_old
  382. write(*,*) 'R_int,Y_int=',R_int,Y_int
  383. write(*,*) 'r0,rho,tr_eps1=',DENS0,VAR0(2),tr_eps1
  384. write(*,*) 'D_end1,gg0=',D_end1,g0
  385. write(*,*) 'alpha,alp_old=',alpha,alp_old
  386. write(*,*) 'DD1,DD2,AL1,AL2=',DD1,DD2,AL1,AL2
  387. write(*,*) 'rho,rho_f=',rho,rho_f
  388. write(*,*) 'DENS0,VAR0(2)=',DENS0,VAR0(2)
  389. ENDIF
  390. *
  391. ELSE
  392. *
  393. * On n'endommage pas
  394. *
  395. VARF(2)=rho
  396. VARF(3)=0.D0
  397. VARF(4)=0.D0
  398. VARF(5)=0.D0
  399. VARF(6)=0.D0
  400. VARF(7)=1.D0
  401. ENDIF
  402. *
  403. *=======================================
  404. * Fin du calcul d'endommagement
  405. *=======================================
  406. ENDIF
  407. *
  408. *
  409. *=========================================================
  410. * Passage a l'option de calcul pour les contraintes
  411. *=========================================================
  412. *
  413. IF (MFR1 .EQ. 1 .OR. MFR1 .EQ. 31) THEN
  414. IF (NSTRS .EQ. 6) THEN
  415. *
  416. *---> MASSIF 3D
  417. *
  418. DO 170 I=1,NSTRS
  419. SIGF(I)=RSIGF(I)
  420. 170 CONTINUE
  421. ELSE IF ( NSTRS .EQ. 4 ) THEN
  422. *
  423. *---> Calcul axisymétrique ou contraintes planes
  424. *
  425. DO 180 I=1,NSTRS
  426. SIGF(I)=RSIGF(I)
  427. 180 CONTINUE
  428. ENDIF
  429. ENDIF
  430. RETURN
  431. *
  432. END
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  

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