Télécharger cfati.eso

Retour à la liste

Numérotation des lignes :

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

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