Télécharger cfati.eso

Retour à la liste

Numérotation des lignes :

  1. C CFATI SOURCE PV 17/12/08 21:15:35 9660
  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.  
  91. -INC PPARAM
  92. -INC CCOPTIO
  93. -INC DECHE
  94. *
  95. SEGMENT WRKK2
  96. REAL*8 EPSILI(NSTRS1)
  97. ENDSEGMENT
  98. C Modif L.Bode - 14/10/92
  99. SEGMENT WRK3
  100. REAL*8 EPSILO(NSTRS1)
  101. ENDSEGMENT
  102. C Definition dynamique de EPSILO
  103. C Fin modif L.Bode
  104. INTEGER NSTRS1,NVARI
  105. REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3)
  106. REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC
  107. REAL*8 YOUN,XNU,EPSD0
  108. C ajout bodin 14 02 2006
  109. REAL*8 BETA,LOI
  110. C pour la loi L2R
  111. REAL*8 C,ALFA
  112. C pour la loi L3R
  113. REAL*8 ALFA1,ALFA2,ALFA3
  114. C fin ajout bodin 14 02 2006
  115. INTEGER ISTRS,JSTRS,KCAS,IRTD
  116. REAL*8 XZERO,UN,DEUX,TROIS
  117. REAL*8 EPSTL0,EPSTM0
  118. REAL*8 DINI,D,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA
  119. REAL*8 dt1
  120. C Fin modif Bodin - 26-09-2005
  121. PARAMETER (XZERO=0.D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0)
  122.  
  123. C
  124. C
  125. C recuperation des variables initiales dans les tableaux
  126. C
  127. C
  128. N=NSTRS1
  129. CMATE = 'ISOTROPE'
  130. YOUN = XMAT(1)
  131. XNU = XMAT(2)
  132. EPSD0= XMAT(5)
  133. BETA = XMAT(6)
  134. LOI = XMAT(7)
  135. C = XMAT(8)
  136. ALFA = XMAT(9)
  137. ALFA1 = XMAT(10)
  138. ALFA2 = XMAT(11)
  139. ALFA3 = XMAT(12)
  140. C ajout bodin 17-11-2000
  141. C introduction de la variable EPSTIMI déformation équivalent initiale
  142. EPSTL0 = VAR0(1)
  143. C fin ajout bodin
  144. DINI = VAR0(2)
  145.  
  146. C ajout bodin 26 09 2005
  147. C récupération des temps temp0 et tempf qui sont dans le segment WRK52
  148. C et calcul du pas de temps dt1
  149. C print*,'temps temp0 ',temp0
  150. C print*,'temps tempf ',tempf
  151. dt1=tempf-temp0
  152. C fin ajout bodin 26 09 2005
  153.  
  154.  
  155. C
  156. C calcul des seuils
  157. C
  158. C
  159. C calcul de la deformation totale
  160. C
  161. SEGINI WRK3
  162. DO 100 ISTRS=1,NSTRS1
  163. EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS)
  164. 100 CONTINUE
  165. C
  166. C calcul des deformations principales
  167. C
  168. C
  169. C on reecrit les deformations sous forme matricielle
  170. C
  171. C Modif L.Bode - 14/10/92
  172. C Rajout de IFOUR2 en argument de ENDOCA
  173. * print*,'on appelle ENDOCB'
  174. CALL ENDOCB (EPSILO,EPS33,2,IFOUR2)
  175. * print*,'apres endocb'
  176. C Fin modif L.Bode
  177. C
  178. C et on diagonalise
  179. C
  180. C Modif L.Bode - 14/10/92
  181. C Pour les elts Coques, on travaille en contraintes planes => JDIM =2
  182. C Pour les elts Massifs JDIM =IDIM
  183. * print*,'avant JACOB3'
  184. CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33)
  185. * print*,'apres JACOB3'
  186. C Fin modif L.Bode
  187. IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN
  188. C
  189. C on calcule la matrice de hooke et les contraintes ppales
  190. C
  191. C CMATE = 'ISOTROPE'
  192. C KCAS=1
  193. C Modif L.Bode - 14/10/92
  194. C IFOUR --> IFOUR2 dans appel DOHMAS
  195. C* print*,'avant dohmas'
  196. C CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS,KCAS,DDHOOK,IRTD)
  197. C* print*,'apres dohmas'
  198. C Fin modif L.Bode
  199. C DO 200 ISTRS=1,3
  200. C SIGP(ISTRS)= XZERO
  201. C DO 210 JSTRS=1,3
  202. C SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  203. C210 CONTINUE
  204. C200 CONTINUE
  205. END IF
  206. C
  207. C on complete la deformation dans le cas des contraintes planes
  208. C
  209. C Modif L.Bode - 14/10/92
  210. C IFOUR remplace par IFOUR2
  211. IF (IFOUR2.EQ. -2) THEN
  212. EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU)
  213. END IF
  214. C Fin modif L.Bode
  215. C
  216. C on calcule le epsilontild
  217. Cmodif bodin 6 juillet 2001 pour le calcul de epsilon equivalente
  218. C on calcule la matrice de hooke et les contraintes ppales
  219. C
  220. CMATE = 'ISOTROPE'
  221. KCAS=1
  222. C Modif L.Bode - 14/10/92
  223. C IFOUR --> IFOUR2 dans appel DOHMAS
  224. * print*,'avant dohmas'
  225. CALL DOHMAS(XMAT,CMATE,IFOUR2,LHOOK,KCAS,DDHOOK,IRTD)
  226. * print*,'apres dohmas'
  227. C Fin modif L.Bode
  228.  
  229. DO 200 ISTRS=1,3
  230. SIGP(ISTRS)= XZERO
  231. DO 210 JSTRS=1,3
  232. SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  233. 210 CONTINUE
  234. 200 CONTINUE
  235. C
  236. C on calcule le signe des contraintes elastiques
  237. C
  238. DO 300 ISTRS=1,3
  239. IF (SIGP(ISTRS).LT. XZERO) THEN
  240. SIGPC(ISTRS)=SIGP(ISTRS)
  241. SIGPT(ISTRS)=XZERO
  242. ELSE
  243. SIGPT(ISTRS)=SIGP(ISTRS)
  244. SIGPC(ISTRS)=XZERO
  245. END IF
  246. 300 CONTINUE
  247. TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  248. TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  249. C
  250. C on calcule les deformations dues aux contraintes positives
  251. C
  252. DO 400 ISTRS=1,3
  253. EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  254. 400 CONTINUE
  255. C
  256. EPSTIL=MAX( XZERO , EPSILT(1) )**2 +
  257. 1 MAX( XZERO , EPSILT(2) )**2 +
  258. 2 MAX( XZERO , EPSILT(3) )**2
  259. EPSTIL=SQRT (EPSTIL)
  260. C
  261. C EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
  262. C 1 MAX( XZERO , EPSIPP(2) )**2 +
  263. C 2 MAX( XZERO , EPSIPP(3) )**2
  264. C EPSTIL=SQRT (EPSTIL)
  265. C fin modif bodin 6 juillet 2001 pour le calcul de epsilon equivalente
  266. IF (ISTEP .EQ. 0) THEN
  267. EPSTIM=EPSTIL
  268. EPSTM0=EPSTL0
  269. VARF(1)=EPSTIL
  270. ELSE IF (ISTEP .EQ. 1) THEN
  271. VARF(2)=DINI
  272. VARF(1)=EPSTIL
  273. C ajout bodin 13-03-2001 pour une nouvelle variable interne
  274. C VARx(3) introduite pour garder la mémoire de la déf eq initiale
  275. C marche en non local si le champ de EPSTL0 est non local par exemple identiquement nul
  276. VARF(3)=EPSTL0
  277. C fin ajout bodin 13-03-2001
  278. ELSE IF (ISTEP .EQ. 2) THEN
  279. C ajout bodin 13-03-2001 pour une nouvelle variable interne
  280. C VARx(3) introduite pour garder la mémoire de la déf eq initiale
  281. EPSTM0=VAR0(3)
  282. C fin ajout bodin 13-03-2001
  283. EPSTIM=VAR0(1)
  284. VARF(1)=EPSTIM
  285. ELSE
  286. PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP
  287. END IF
  288. IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN
  289. C
  290. C on calcule l'endommagement et les contraintes
  291. C
  292. C
  293. C on calcule ALFAT ALFAC DT et DC puis D
  294. C dans le cas ou le seuil initial est depasse
  295. C
  296. IF ( EPSTIM .GT. EPSD0) THEN
  297. C
  298. C calcul de l'endommagement
  299. C
  300. C
  301. C on calcule le signe des contraintes elastiques
  302. C
  303. C DO 300 ISTRS=1,3
  304. C IF (SIGP(ISTRS).LT. XZERO) THEN
  305. C SIGPC(ISTRS)=SIGP(ISTRS)
  306. C SIGPT(ISTRS)=XZERO
  307. C ELSE
  308. C SIGPT(ISTRS)=SIGP(ISTRS)
  309. C SIGPC(ISTRS)=XZERO
  310. C END IF
  311. C300 CONTINUE
  312. C TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  313. C TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  314. C
  315. C on calcule les deformations dues aux contraintes positives
  316. C
  317. C DO 400 ISTRS=1,3
  318. C EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  319. C400 CONTINUE
  320. C
  321. C on en deduit ALFAT et ALFAC
  322. C
  323. C ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
  324. C 1 MAX(XZERO,EPSIPP(2))*EPSILT(2) +
  325. C 2 MAX(XZERO,EPSIPP(3))*EPSILT(3)
  326. C ALFAT = ALFAT/(EPSTIL*EPSTIL)
  327. C ALFAC = UN - ALFAT
  328. C PRINT*,'alfat = ',ALFAT,' et alfac = ',ALFAC
  329. C PRINT*,'alfa1=',A99,'beta1 =',B99,'C1 = ',C99
  330. C
  331. C modification pour la bi ou tricompression
  332. C
  333. C IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN
  334. C GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
  335. C 1 SIGPC(3)*SIGPC(3)
  336. C GAMMA=-SQRT(GAMMA)/TRSIGC
  337. C EPSTIM=EPSTIM*GAMMA
  338. C END IF
  339. C
  340. C amelioration de la reponse en cisaillement pour beta > 1.
  341. C
  342. C IF (BETA .GT. UN) THEN
  343. C IF ( ALFAT .GT. XPETIT ) THEN
  344. C ALFAT=ALFAT**BETA
  345. C END IF
  346. C IF ( ALFAC .GT. XPETIT ) THEN
  347. C ALFAC=ALFAC**BETA
  348. C END IF
  349. C END IF
  350.  
  351. C
  352. C on calcule DT et DC puis D
  353. C dans le cas ou le seuil initial est depasse
  354. C
  355. C on est oblige de verifier car on a pu multiplier par gamma
  356. C
  357. C modif bodin 26-02-2001
  358. C calcul incrémental pour le pas EPSTIM-EPSTIMI de déformation
  359. C équivalente
  360. C
  361. C PRINT*,'La def initiale est',EPSTM0
  362. C PRINT*,'La def finale est',EPSTIM
  363. C
  364. IF (LOI.EQ.DEUX) THEN
  365. D=((DINI**(1-ALFA))+
  366. 1 ((C*(1-ALFA)/(BETA+1))*
  367. 2 ((EPSTIM**(BETA+1))-(EPSD0**(BETA+1)))*(dt1*2)
  368. 3 ))**(1/(1-ALFA))
  369. C
  370. ELSE IF (LOI.EQ.TROIS) THEN
  371. C on traite la loi L3R
  372. C calcul d'un term intermediaire
  373. TERM_A= (((ALFA1*(UN -EXP(-1*((DINI/ALFA2)**ALFA3))))+
  374. 1 ((((EPSTIM**(BETA+1))-(EPSD0**(BETA+1)))
  375. 2 /(BETA+UN))*(dt1*2))
  376. 3 )/ALFA1)
  377. C
  378. IF (TERM_A .LT. UN) THEN
  379. D= ALFA2 *((-1*(LOG(UN - TERM_A)))**(UN/ALFA3))
  380. ELSE
  381. D=UN
  382. END IF
  383. C sinon roblème de définition de LOI
  384. ELSE
  385. WRITE(*,*) 'LOI doit etre egal à 2. ou 3.'
  386. STOP
  387. END IF
  388. C
  389. C on borne la valeur de D a 0.99999
  390. C
  391. D=MIN ( D , UN-1.D-05 )
  392. ELSE
  393. C on n'a pas dépassé le seuil
  394. D=XZERO
  395. END IF
  396. C
  397. C on teste la croissance de D
  398. C
  399. D=MAX ( D , DINI )
  400. C
  401. C on le stocke dans les variables internes finales
  402. C
  403. VARF(2)= D
  404. C
  405. C on calcule les contraintes finales
  406. C
  407. * CALL MATVE1 (DDHOOK,EPSILO,NSTRS,NSTRS,SIGF,2)
  408. DO istrs=1,NSTRS
  409. SIGF(istrs)=0.d0
  410. ENDDO
  411. DO jstrs=1,NSTRS
  412. DO istrs=1,nstrs
  413. SIGF(jstrs)=SIGF(jstrs)+DDHOOK(jstrs,istrs)*epsilo(istrs)
  414. ENDDO
  415. ENDDO
  416.  
  417.  
  418. DO 500 ISTRS=1,NSTRS
  419. SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
  420. 500 CONTINUE
  421. C
  422. C et les deformations inelastiques finales
  423. C
  424. DO 600 ISTRS=1,NSTRS
  425. EPINF(ISTRS)=EPSILO(ISTRS)*D
  426. 600 CONTINUE
  427.  
  428. END IF
  429. SEGSUP WRK3
  430. RETURN
  431. END
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  

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