Télécharger modvms.eso

Retour à la liste

Numérotation des lignes :

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

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