Télécharger codvms.eso

Retour à la liste

Numérotation des lignes :

  1. C CODVMS SOURCE BP208322 17/03/01 21:16:00 9325
  2. SUBROUTINE CODVMS(WRK52,WRK53,WRK54,WRKK2,NSTRS1,NVARI,
  3. $ icara,jdim,ifour2)
  4. C MODVMS SOURCE CHAT 03/01/24 21:27:39 4558
  5. C MODVONMISES
  6. c SUBROUTINE MODVONMISES (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI,
  7. c 1 NMATT,ISTEP,ICARA,JDIM,IFOUR2)
  8. C
  9. * BCN
  10. C New source: Modified Von Mises nonlocal damage model added (inspired in MAZARS)
  11. C Goal: satisfactory simulation of the single-edge notched beam
  12. C
  13. C Modifications:
  14. C 1. The definition of the equivalent strain (epsilontild)
  15. C is changed [1]
  16. C 2. The distinction between damage in traction (DT) and damage
  17. C in compression (DC) is suppressed. Only one damage (D) is used
  18. C 3. The relationship between the nonlocal equivalent strain and the
  19. C damage is changed : Two choices for damage vs. nonlocal equivalent
  20. C strain. [2]
  21. C a) LOI = 1. Exponential law.
  22. C b) LOI = 0. Polynomial law.
  23. C
  24. C References:
  25. C [1] Peerlings, R.H.J., de Borst, R., Brekelmans, W.A.M., and Geers, M.D.
  26. C (1998), Gradient-enhanced damage modelling of concrete fracture,
  27. C Mechanics of Cohesive-Frictional Materials, 3, 323-342.
  28. C [2] Rodriguez-Ferran A., Huerta A. (2000), Error estimation and
  29. C adaptivity for nonlocal damage models. International Journal of
  30. C Solids and Structures, 37, 7501-7528.
  31. * BCN
  32. C
  33. C variables en entree
  34. C
  35. C
  36. C WRKK2 pointeur sur un segment variables internes au pas precedent
  37. C
  38. C XMATER constantes du materiau
  39. C
  40. C NSTRS1 nombre de composantes dans les vecteurs des contraintes
  41. C et les vecteurs des deformations
  42. C
  43. C NVARI nombre de variables internes (doit etre egal a 2)
  44. C
  45. C NMATT nombre de constantes du materiau
  46. C
  47. C ISTEP flag utilise pour separer les etapes dans un calcul non local
  48. C ISTEP=0 -----> calcul local
  49. C ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
  50. C ISTEP=2 -----> calcul non local etape 2 on continue le calcul
  51. C a partir des seuils moyennes
  52. C
  53. C Modif L.Bode - 14/10/92
  54. C Nouveaux parametres en entree
  55. C JDIM Dimension de travail
  56. C ( Coques JDIM =2 , Massifs JDIM = IDIM )
  57. C IFOUR2 Type de formulation
  58. C ( Coques IFOUR2 = -2 => contraintes planes ,
  59. C Massifs IFOUR2 = IFOUR)
  60. C
  61. C variables en sortie
  62. C
  63. C VARINF variables internes finales
  64. C
  65. C SIGMAF contraintes finales
  66. C
  67. C C. LA BORDERIE MARS 1992
  68. C declaration des variables
  69. C
  70. C
  71. IMPLICIT INTEGER(I-N)
  72. IMPLICIT REAL*8(A-H,O-Z)
  73. -INC CCOPTIO
  74. -INC DECHE
  75. SEGMENT WRKK2
  76. REAL*8 EPSILI(nstrs1)
  77. ENDSEGMENT
  78. C Modif L.Bode - 14/10/92
  79. SEGMENT WRK3
  80. REAL*8 EPSILO(nstrs1)
  81. ENDSEGMENT
  82. C Definition dynamique de EPSILO
  83. C Fin modif L.Bode
  84. c INTEGER nstrs1,NVARI,NMATT,ISTEP
  85. REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3)
  86. * BCN
  87. C REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC
  88. C REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA
  89. REAL*8 SIGP(3)
  90. REAL*8 YOUN,XNU,EPSD0,B1,B2,RATIO_CT,LOI
  91. * BCN
  92. c INTEGER ISTRS,JSTRS,KCAS,IRTD
  93. c REAL*8 XZERO,UN,DEUX,XPETIT
  94. * BCN
  95. C REAL*8 DINI,D,DT,DC,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA
  96. * REAL*8 DINI,D,EPSTIL,EPSTIM
  97. * REAL*8 TINV1,TINV2,TINV3,AUX1,AUX2,AUX3,AUX4,AUX5
  98. * BCN
  99. PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  100.  
  101. C
  102. C
  103. C recuperation des variables initiales dans les tableaux
  104. C
  105. C
  106. N=nstrs1
  107. CMATE = 'ISOTROPE'
  108. YOUN = XMAT(1)
  109. XNU = XMAT(2)
  110. EPSD0= XMAT(5)
  111. * BCN
  112. C The material parameters change due to the modification of the model
  113. C EPSD0: damage threshold [2]
  114. C B1 and B2: material parameters in equivalent strain-damage law [2]
  115. C RATIO_CT: ratio of compressive to tensile strength [1,2]
  116. C LOI: Choice of equivalent strain-damage law [2]
  117. C ACOM = XMAT(6)
  118. C BCOM = XMAT(7)
  119. C ATRA = XMAT(8)
  120. C BTRA = XMAT(9)
  121. C BETA = XMAT(10)
  122. B1 = XMAT(6)
  123. B2 = XMAT(7)
  124. RATIO_CT = XMAT(8)
  125. LOI = XMAT(9)
  126. * BCN
  127. DINI = VAR0(2)
  128. C
  129. C calcul des seuils
  130. C
  131. C
  132. C calcul de la deformation totale
  133. C
  134. SEGINI WRK3
  135. DO 100 ISTRS=1,nstrs1
  136. EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS)
  137. 100 CONTINUE
  138. C
  139. C calcul des deformations principales
  140. C
  141. C
  142. C on reecrit les deformations sous forme matricielle
  143. C
  144. C Modif L.Bode - 14/10/92
  145. C Rajout de IFOUR2 en argument de ENDOCA
  146. * print*,'on appelle ENDOCB'
  147. CALL ENDOCB (EPSILO,EPS33,2,IFOUR2)
  148. * print*,'apres endocb'
  149. C Fin modif L.Bode
  150. C
  151. C et on diagonalise
  152. C
  153. C Modif L.Bode - 14/10/92
  154. C Pour les elts Coques, on travaille en contraintes planes => JDIM =2
  155. C Pour les elts Massifs JDIM =IDIM
  156. * print*,'avant JACOB3'
  157. CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33)
  158. * print*,'apres JACOB3'
  159. C Fin modif L.Bode
  160. IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN
  161. C
  162. C on calcule la matrice de hooke et les contraintes ppales
  163. C
  164. CMATE = 'ISOTROPE'
  165. KCAS=1
  166. C Modif L.Bode - 14/10/92
  167. C IFOUR --> IFOUR2 dans appel DOHMAS
  168. * print*,'avant dohmas'
  169. CALL DOHMAS(XMAT,CMATE,IFOUR2,nstrs1,KCAS,DDHOOK,IRTD)
  170. * print*,'apres dohmas'
  171. C Fin modif L.Bode
  172. DO 200 ISTRS=1,3
  173. SIGP(ISTRS)= XZERO
  174. DO 210 JSTRS=1,3
  175. SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  176. 210 CONTINUE
  177. 200 CONTINUE
  178. END IF
  179. C
  180. C on complete la deformation dans le cas des contraintes planes
  181. C
  182. C Modif L.Bode - 14/10/92
  183. C IFOUR remplace par IFOUR2
  184. IF (IFOUR2.EQ. -2) THEN
  185. EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU)
  186. END IF
  187. C Fin modif L.Bode
  188. C
  189. C on calcule le epsilontild
  190. C
  191. * BCN
  192. C The original definition of epsilontild is replaced.
  193. C TINV1: first strain tensor invariant
  194. C TINV2: second (deviatoric) strain tensor invariant
  195. C
  196. C IMPORTANT! There is a discrepancy by a factor of 6 in the
  197. C definitions of the second invariant in [1] and [2]
  198. C The definition of Askes and Sluys is used here.
  199. C
  200. C EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
  201. C 1 MAX( XZERO , EPSIPP(2) )**2 +
  202. C 2 MAX( XZERO , EPSIPP(3) )**2
  203. C EPSTIL=SQRT (EPSTIL)
  204. C
  205. TINV1 = EPSIPP(1)+EPSIPP(2)+EPSIPP(3)
  206. TINV2 = (EPSIPP(1)- EPSIPP(2))*(EPSIPP(1)- EPSIPP(2)) +
  207. . (EPSIPP(1)- EPSIPP(3))*(EPSIPP(1)- EPSIPP(3)) +
  208. . (EPSIPP(2)- EPSIPP(3))*(EPSIPP(2)- EPSIPP(3))
  209. TINV2 = TINV2/6.d0
  210. C
  211. AUX1 = RATIO_CT-UN
  212. AUX2 = UN+XNU
  213. AUX3 = UN-DEUX*XNU
  214. AUX4 = AUX1*TINV1/AUX3
  215. AUX5 = SQRT((AUX4*AUX4)+((12.d0*RATIO_CT*TINV2)/(AUX2*AUX2)))
  216. C
  217. EPSTIL = (AUX1*TINV1)/(DEUX*RATIO_CT*AUX3) +
  218. . AUX5/(DEUX*RATIO_CT)
  219. * BCN
  220. IF (ISTEP .EQ. 0) THEN
  221. EPSTIM=EPSTIL
  222. VARF(1)=EPSTIL
  223. ELSE IF (ISTEP .EQ. 1) THEN
  224. VARF(2)=DINI
  225. VARF(1)=EPSTIL
  226. ELSE IF (ISTEP .EQ. 2) THEN
  227. EPSTIM=VAR0(1)
  228. VARF(1)=EPSTIM
  229. ELSE
  230. PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP
  231. END IF
  232. IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN
  233. C
  234. C on calcule l'endommagement et les contraintes
  235. C
  236. C
  237. C on calcule ALFAT ALFAC DT et DC puis D
  238. C dans le cas ou le seuil initial est depasse
  239. C
  240. IF ( EPSTIM .GT. EPSD0) THEN
  241. C
  242. C calcul de l'endommagement
  243. * BCN
  244. C The equivalent strain-damage law is modified
  245. C
  246. C on calcule le signe des contraintes elastiques
  247. C
  248. C DO 300 ISTRS=1,3
  249. C IF (SIGP(ISTRS).LT. XZERO) THEN
  250. C SIGPC(ISTRS)=SIGP(ISTRS)
  251. C SIGPT(ISTRS)=XZERO
  252. C ELSE
  253. C SIGPT(ISTRS)=SIGP(ISTRS)
  254. C SIGPC(ISTRS)=XZERO
  255. C END IF
  256. C300 CONTINUE
  257. C TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  258. C TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  259. C
  260. C on calcule les deformations dues aux contraintes positives
  261. C
  262. C DO 400 ISTRS=1,3
  263. C EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  264. C400 CONTINUE
  265. C
  266. C on en deduit ALFAT et ALFAC
  267. C
  268. C ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
  269. C 1 MAX(XZERO,EPSIPP(2))*EPSILT(2) +
  270. C 2 MAX(XZERO,EPSIPP(3))*EPSILT(3)
  271. C ALFAT = ALFAT/(EPSTIL*EPSTIL)
  272. C ALFAC = UN - ALFAT
  273. C
  274. C modification pour la bi ou tricompression
  275. C
  276. C IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN
  277. C GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
  278. C 1 SIGPC(3)*SIGPC(3)
  279. C GAMMA=-SQRT(GAMMA)/TRSIGC
  280. C EPSTIM=EPSTIM*GAMMA
  281. C END IF
  282. C
  283. C amelioration de la reponse en cisaillement pour beta > 1.
  284. C
  285. C IF (BETA .GT. UN) THEN
  286. C IF ( ALFAT .GT. XPETIT ) THEN
  287. C ALFAT=ALFAT**BETA
  288. C END IF
  289. C IF ( ALFAC .GT. XPETIT ) THEN
  290. C ALFAC=ALFAC**BETA
  291. C END IF
  292. C END IF
  293. C
  294. C
  295. C on calcule DT et DC puis D
  296. C dans le cas ou le seuil initial est depasse
  297. C
  298. C on est oblige de verifier car on a pu multiplier par gamma
  299. C
  300. C IF (EPSTIM .GT. EPSD0) THEN
  301. C DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
  302. C 1 ATRA*EXP(-BTRA*(EPSTIM-EPSD0))
  303. C DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
  304. C 1 ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
  305. C ELSE
  306. C DT=XZERO
  307. C DC=XZERO
  308. C END IF
  309. C D = ALFAT*DT + ALFAC*DC
  310. C
  311. C Two choices for damage vs. internal variable:
  312. ***tc Y0 n'etait pas initialisé on suppose que c'est EPSD0
  313. Y0 = epsd0
  314. IF (LOI.EQ.XZERO) THEN
  315. D = UN - (UN/(UN+B1*(EPSTIM-Y0)+
  316. . B2*(EPSTIM-Y0)*(EPSTIM-Y0)))
  317. ELSE IF (LOI.EQ.UN) THEN
  318. D = UN - (Y0*(UN-B2)/EPSTIM)- B2*EXP(B1*(Y0-EPSTIM))
  319. ELSE
  320. WRITE(*,*) 'LOI SHOULD BE EQUAL TO 1 OR 0.(SEE MATE)'
  321. STOP
  322. ENDIF
  323. * BCN
  324. C
  325. C on borne la valeur de D a 0.99999
  326. C
  327. D=MIN ( D , UN-1.D-05 )
  328. ELSE
  329. D=XZERO
  330. END IF
  331. C
  332. C on teste la croissance de D
  333. C
  334. D=MAX ( D , DINI )
  335. C
  336. C on le stocke dans les variables internes finales
  337. C
  338. VARF(2)= D
  339. C
  340. C on calcule les contraintes finales
  341. C
  342. CALL MATVE1 (DDHOOK,EPSILO,nstrs1,nstrs1,SIGF,2)
  343. DO 500 ISTRS=1,nstrs1
  344. SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
  345. 500 CONTINUE
  346. C
  347. C et les deformations inelastiques finales
  348. C
  349. DO 600 ISTRS=1,nstrs1
  350. EPINF(ISTRS)=EPSILO(ISTRS)*D
  351. 600 CONTINUE
  352.  
  353. END IF
  354. SEGSUP WRK3
  355. RETURN
  356. END
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  

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