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

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