Télécharger fati.eso

Retour à la liste

Numérotation des lignes :

  1. C FATI SOURCE CHAT 07/12/04 21:15:55 5985
  2. C FATI
  3. SUBROUTINE FATI (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI,
  4. 1 NMATT,ISTEP,ICARA,JDIM,IFOUR2,DT)
  5. C MAZARS SOURCE AM 00/12/13 21:39:37 4045
  6. C SUBROUTINE MAZARS (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI,
  7. C 1 NMATT,ISTEP,ICARA,JDIM,IFOUR2)
  8. *======================================================================
  9. * DB
  10. * FOR THE NONLOCAL CASE CALLED by ECOU40 --> FATTT --> FATI
  11. * New source: one nonlocal damage material model added:
  12. * Fatigue damage
  13. * This routine is inspired based on mazars.eso
  14. *
  15. C Goal: allowing cyclic cumulative damage calculation under simple fatigue loading
  16. C
  17. C Description:
  18. C 1. The threshold is maintiened allowing an endurance limit but not tested
  19. C 2. The equivalent strain is changed replaced by the geometric average of the positive stresses
  20. C induced strain ( no damage growth due to compression induced tension - Poisson effect)
  21. C 3. The damage rate given in [1] is integrated over a cycle with
  22. C dD/dN = f(D) eps_a^(beta+1)/(beta+1)
  23. C 4. Two expression of the f(D) funcion are implemented
  24. C a) LOI = 2. power law as chosen in [2]
  25. C b) LOI = 3. power law including a damage rate acceleration as describe in [1]
  26. C
  27. C References:
  28. C [1] D. Bodin, G. Pijaudier-Cabot, C. de La Roche, J.-M Piau and A. Chabot, (2004)
  29. C A Continuum Damage Approach to Asphalt Concrete Fatigue Modelling, Journal of
  30. C Engineering Mechanics, ASCE, vol. 130 (6), pp. 700-708.
  31. C
  32. C [2] Paas, R. H. J. W., Scheurs, P. J. G., and Brekelmans, W. A. M. (1993). ‘A
  33. C continuum approach to brittle and fatigue damage: Theory and numerical
  34. C procedures.’ Int. J. Solids Struct., 30~4!, 579–599.
  35. C
  36. * DB
  37. *
  38. *======================================================================
  39. C
  40. C
  41. C
  42. C
  43. C variables en entree
  44. C
  45. C
  46. C WRK0 pointeur sur un segment deformation au pas precedent
  47. C
  48. C WRK1 pointeur sur un segment increment de deformation
  49. C
  50. C WRKK2 pointeur sur un segment variables internes au pas precedent
  51. C
  52. C WRK5 pointeur sur un segment de deformations inelastiques
  53. C
  54. C XMATER constantes du materiau
  55. C
  56. C NSTRS nombre de composantes dans les vecteurs des contraintes
  57. C et les vecteurs des deformations
  58. C
  59. C NVARI nombre de variables internes (doit etre egal a 2)
  60. C
  61. C NMATT nombre de constantes du materiau
  62. C
  63. C ISTEP flag utilise pour separer les etapes dans un calcul non local
  64. C ISTEP=0 -----> calcul local
  65. C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
  66. C ISTEP=2 -----> calcul non local etape 2 on continue le calcul
  67. C a partir des seuils moyennes
  68. C
  69. C Modif L.Bode - 14/10/92
  70. C Nouveaux parametres en entree
  71. C JDIM Dimension de travail
  72. C ( Coques JDIM =2 , Massifs JDIM = IDIM )
  73. C IFOUR2 Type de formulation
  74. C ( Coques IFOUR2 = -2 => contraintes planes ,
  75. C Massifs IFOUR2 = IFOUR)
  76. C
  77. C variables en sortie
  78. C
  79. C VARINF variables internes finales
  80. C
  81. C SIGMAF contraintes finales
  82. C
  83. C C. LA BORDERIE MARS 1992
  84. C declaration des variables
  85. C
  86. C
  87. IMPLICIT REAL*8(A-H,O-Z)
  88. -INC CCOPTIO
  89. SEGMENT WRK0
  90. REAL*8 XMAT(NMATT)
  91. ENDSEGMENT
  92. *
  93. SEGMENT WRK1
  94. REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS)
  95. REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI)
  96. REAL*8 DEFP(NSTRS),XCAR(ICARA)
  97. ENDSEGMENT
  98. SEGMENT WRKK2
  99. REAL*8 EPSILI(NSTRS),DSIGT(NSTRS)
  100. ENDSEGMENT
  101. C Modif L.Bode - 14/10/92
  102. SEGMENT WRK3
  103. REAL*8 EPSILO(NSTRS)
  104. ENDSEGMENT
  105. C Definition dynamique de EPSILO
  106. C Fin modif L.Bode
  107. SEGMENT WRK5
  108. REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS)
  109. ENDSEGMENT
  110. C
  111. CHARACTER*8 CMATE
  112. INTEGER NSTRS,NVARI,NMATT,ISTEP
  113. REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3)
  114. REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC
  115. C REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA
  116. REAL*8 YOUN,XNU,EPSD0
  117. C ajout bodin 14 02 2006
  118. REAL*8 BETA,LOI
  119. C pour la loi L2R
  120. REAL*8 C,ALFA
  121. C pour la loi L3R
  122. REAL*8 ALFA1,ALFA2,ALFA3
  123. C fin ajout bodin 14 02 2006
  124. INTEGER ISTRS,JSTRS,KCAS,IRTD
  125. REAL*8 XZERO,UN,DEUX,TROIS
  126. REAL*8 EPSTL0,EPSTM0
  127. REAL*8 DINI,D,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA
  128. REAL*8 dt1
  129. C Fin modif Bodin - 26-09-2005
  130. PARAMETER (XZERO=0.D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0)
  131.  
  132. C
  133. C
  134. C recuperation des variables initiales dans les tableaux
  135. C
  136. C
  137. N=NSTRS
  138. CMATE = 'ISOTROPE'
  139. YOUN = XMAT(1)
  140. XNU = XMAT(2)
  141. EPSD0= XMAT(5)
  142. BETA = XMAT(6)
  143. LOI = XMAT(7)
  144. C = XMAT(8)
  145. ALFA = XMAT(9)
  146. ALFA1 = XMAT(10)
  147. ALFA2 = XMAT(11)
  148. ALFA3 = XMAT(12)
  149. C ajout bodin 17-11-2000
  150. C introduction de la variable EPSTIMI déformation équivalent initiale
  151. EPSTL0 = VAR0(1)
  152. C fin ajout bodin
  153. DINI = VAR0(2)
  154. C
  155. C
  156. C print*,'Je suis dans mazars'
  157. C
  158. C ajout bodin 26 09 2005
  159. C print*,'mazars Variable DT ',DT
  160. dt1 =DT
  161. C print*,'mazars variable dt1',dt1
  162. C fin ajout bodin 26 09 2005
  163.  
  164. C
  165. C calcul des seuils
  166. C
  167. C
  168. C calcul de la deformation totale
  169. C
  170. SEGINI WRK3
  171. DO 100 ISTRS=1,NSTRS
  172. EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS)
  173. 100 CONTINUE
  174. C
  175. C calcul des deformations principales
  176. C
  177. C
  178. C on reecrit les deformations sous forme matricielle
  179. C
  180. C Modif L.Bode - 14/10/92
  181. C Rajout de IFOUR2 en argument de ENDOCA
  182. * print*,'on appelle ENDOCB'
  183. CALL ENDOCB (EPSILO,EPS33,2,IFOUR2)
  184. * print*,'apres endocb'
  185. C Fin modif L.Bode
  186. C
  187. C et on diagonalise
  188. C
  189. C Modif L.Bode - 14/10/92
  190. C Pour les elts Coques, on travaille en contraintes planes => JDIM =2
  191. C Pour les elts Massifs JDIM =IDIM
  192. * print*,'avant JACOB3'
  193. CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33)
  194. * print*,'apres JACOB3'
  195. C Fin modif L.Bode
  196. IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN
  197. C
  198. C on calcule la matrice de hooke et les contraintes ppales
  199. C
  200. C CMATE = 'ISOTROPE'
  201. C KCAS=1
  202. C Modif L.Bode - 14/10/92
  203. C IFOUR --> IFOUR2 dans appel DOHMAS
  204. C* print*,'avant dohmas'
  205. C CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS,KCAS,DDHOOK,IRTD)
  206. C* print*,'apres dohmas'
  207. C Fin modif L.Bode
  208. C DO 200 ISTRS=1,3
  209. C SIGP(ISTRS)= XZERO
  210. C DO 210 JSTRS=1,3
  211. C SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  212. C210 CONTINUE
  213. C200 CONTINUE
  214. END IF
  215. C
  216. C on complete la deformation dans le cas des contraintes planes
  217. C
  218. C Modif L.Bode - 14/10/92
  219. C IFOUR remplace par IFOUR2
  220. IF (IFOUR2.EQ. -2) THEN
  221. EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU)
  222. END IF
  223. C Fin modif L.Bode
  224. C
  225. C on calcule le epsilontild
  226. Cmodif bodin 6 juillet 2001 pour le calcul de epsilon equivalente
  227. C on calcule la matrice de hooke et les contraintes ppales
  228. C
  229. CMATE = 'ISOTROPE'
  230. KCAS=1
  231. C Modif L.Bode - 14/10/92
  232. C IFOUR --> IFOUR2 dans appel DOHMAS
  233. * print*,'avant dohmas'
  234. CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS,KCAS,DDHOOK,IRTD)
  235. * print*,'apres dohmas'
  236. C Fin modif L.Bode
  237.  
  238. DO 200 ISTRS=1,3
  239. SIGP(ISTRS)= XZERO
  240. DO 210 JSTRS=1,3
  241. SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  242. 210 CONTINUE
  243. 200 CONTINUE
  244. C
  245. C on calcule le signe des contraintes elastiques
  246. C
  247. DO 300 ISTRS=1,3
  248. IF (SIGP(ISTRS).LT. XZERO) THEN
  249. SIGPC(ISTRS)=SIGP(ISTRS)
  250. SIGPT(ISTRS)=XZERO
  251. ELSE
  252. SIGPT(ISTRS)=SIGP(ISTRS)
  253. SIGPC(ISTRS)=XZERO
  254. END IF
  255. 300 CONTINUE
  256. TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  257. TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  258. C
  259. C on calcule les deformations dues aux contraintes positives
  260. C
  261. DO 400 ISTRS=1,3
  262. EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  263. 400 CONTINUE
  264. C
  265. EPSTIL=MAX( XZERO , EPSILT(1) )**2 +
  266. 1 MAX( XZERO , EPSILT(2) )**2 +
  267. 2 MAX( XZERO , EPSILT(3) )**2
  268. EPSTIL=SQRT (EPSTIL)
  269. C
  270. C EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
  271. C 1 MAX( XZERO , EPSIPP(2) )**2 +
  272. C 2 MAX( XZERO , EPSIPP(3) )**2
  273. C EPSTIL=SQRT (EPSTIL)
  274. C fin modif bodin 6 juillet 2001 pour le calcul de epsilon equivalente
  275. IF (ISTEP .EQ. 0) THEN
  276. EPSTIM=EPSTIL
  277. EPSTM0=EPSTL0
  278. VARF(1)=EPSTIL
  279. ELSE IF (ISTEP .EQ. 1) THEN
  280. VARF(2)=DINI
  281. VARF(1)=EPSTIL
  282. C ajout bodin 13-03-2001 pour une nouvelle variable interne
  283. C VARx(3) introduite pour garder la mémoire de la déf eq initiale
  284. C marche en non local si le champ de EPSTL0 est non local par exemple identiquement nul
  285. VARF(3)=EPSTL0
  286. C fin ajout bodin 13-03-2001
  287. ELSE IF (ISTEP .EQ. 2) THEN
  288. C ajout bodin 13-03-2001 pour une nouvelle variable interne
  289. C VARx(3) introduite pour garder la mémoire de la déf eq initiale
  290. EPSTM0=VAR0(3)
  291. C fin ajout bodin 13-03-2001
  292. EPSTIM=VAR0(1)
  293. VARF(1)=EPSTIM
  294. ELSE
  295. PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP
  296. END IF
  297. IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN
  298. C
  299. C on calcule l'endommagement et les contraintes
  300. C
  301. C
  302. C on calcule ALFAT ALFAC DT et DC puis D
  303. C dans le cas ou le seuil initial est depasse
  304. C
  305. C
  306. IF ( EPSTIM .GT. EPSD0) THEN
  307. C
  308. C calcul de l'endommagement
  309. C
  310. C
  311. C on calcule le signe des contraintes elastiques
  312. C
  313. C DO 300 ISTRS=1,3
  314. C IF (SIGP(ISTRS).LT. XZERO) THEN
  315. C SIGPC(ISTRS)=SIGP(ISTRS)
  316. C SIGPT(ISTRS)=XZERO
  317. C ELSE
  318. C SIGPT(ISTRS)=SIGP(ISTRS)
  319. C SIGPC(ISTRS)=XZERO
  320. C END IF
  321. C300 CONTINUE
  322. C TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  323. C TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  324. C
  325. C on calcule les deformations dues aux contraintes positives
  326. C
  327. C DO 400 ISTRS=1,3
  328. C EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  329. C400 CONTINUE
  330. C
  331. C on en deduit ALFAT et ALFAC
  332. C
  333. C ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
  334. C 1 MAX(XZERO,EPSIPP(2))*EPSILT(2) +
  335. C 2 MAX(XZERO,EPSIPP(3))*EPSILT(3)
  336. C ALFAT = ALFAT/(EPSTIL*EPSTIL)
  337. C ALFAC = UN - ALFAT
  338. C PRINT*,'alfat = ',ALFAT,' et alfac = ',ALFAC
  339. C PRINT*,'alfa1=',A99,'beta1 =',B99,'C1 = ',C99
  340. C
  341. C modification pour la bi ou tricompression
  342. C
  343. C IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN
  344. C GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
  345. C 1 SIGPC(3)*SIGPC(3)
  346. C GAMMA=-SQRT(GAMMA)/TRSIGC
  347. C EPSTIM=EPSTIM*GAMMA
  348. C END IF
  349. C
  350. C amelioration de la reponse en cisaillement pour beta > 1.
  351. C
  352. C IF (BETA .GT. UN) THEN
  353. C IF ( ALFAT .GT. XPETIT ) THEN
  354. C ALFAT=ALFAT**BETA
  355. C END IF
  356. C IF ( ALFAC .GT. XPETIT ) THEN
  357. C ALFAC=ALFAC**BETA
  358. C END IF
  359. C END IF
  360. C
  361. C
  362. C on calcule DT et DC puis D
  363. C dans le cas ou le seuil initial est depasse
  364. C
  365. C on est oblige de verifier car on a pu multiplier par gamma
  366. C
  367. C modif bodin 26-02-2001
  368. C modif bodin 26-02-2001
  369. C calcul incrémental pour le pas EPSTIM-EPSTIMI de déformation
  370. C équivalente
  371. C
  372. C
  373. C
  374. C
  375. C
  376. C
  377. C PRINT*,'Je suis dans mazars_bod et je suis content !!'
  378. C PRINT*,'La def initiale est',EPSTM0
  379. C PRINT*,'La def finale est',EPSTIM
  380.  
  381. IF (LOI.EQ.DEUX) THEN
  382. D=((DINI**(1-ALFA))+
  383. 1 ((C*(1-ALFA)/(BETA+1))*
  384. 2 ((EPSTIM**(BETA+1))-(EPSD0**(BETA+1)))*(dt1*2)
  385. 3 ))**(1/(1-ALFA))
  386. C
  387. ELSE IF (LOI.EQ.TROIS) THEN
  388. C on traite la loi L3R
  389. C calcul d'un term intermediaire
  390. TERM_A= (((ALFA1*(UN -EXP(-1*((DINI/ALFA2)**ALFA3))))+
  391. 1 ((((EPSTIM**(BETA+1))-(EPSD0**(BETA+1)))
  392. 2 /(BETA+UN))*(dt1*2))
  393. 3 )/ALFA1)
  394. C
  395. IF (TERM_A .LT. UN) THEN
  396. D= ALFA2 *((-1*(LOG(UN - TERM_A)))**(UN/ALFA3))
  397. ELSE
  398. D=UN
  399. END IF
  400. C sinon roblème de définition de LOI
  401. ELSE
  402. WRITE(*,*) 'LOI doit etre egal à 2. ou 3.'
  403. STOP
  404. END IF
  405. C
  406. C on borne la valeur de D a 0.99999
  407. C
  408. D=MIN ( D , UN-1.D-05 )
  409. ELSE
  410. C on n'a pas dépassé le seuil
  411. D=XZERO
  412. END IF
  413. C
  414. C on teste la croissance de D
  415. C
  416. D=MAX ( D , DINI )
  417. C
  418. C on le stocke dans les variables internes finales
  419. C
  420. VARF(2)= D
  421. C
  422. C on calcule les contraintes finales
  423. C
  424. CALL MATVE1 (DDHOOK,EPSILO,NSTRS,NSTRS,SIGF,2)
  425. DO 500 ISTRS=1,NSTRS
  426. SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
  427. 500 CONTINUE
  428. C
  429. C et les deformations inelastiques finales
  430. C
  431. DO 600 ISTRS=1,NSTRS
  432. EPINF(ISTRS)=EPSILO(ISTRS)*D
  433. 600 CONTINUE
  434.  
  435. END IF
  436. SEGSUP WRK3
  437. RETURN
  438. END
  439.  
  440.  

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