Télécharger prvik2.eso

Retour à la liste

Numérotation des lignes :

  1. C PRVIK2 SOURCE BP208322 17/03/01 21:18:02 9325
  2. SUBROUTINE PRVIK2(SIG0,NSTRS,DEPST,VAR0,XMAT,NCOMAT,XCAR,ICARA,
  3. 1 NVARI,SIGF,VARF,DEFP,MFR1,DDAUX,CMATE,VALMAT,VALCAR,N2EL,
  4. 2 N2PTEL,NBPGAU,IFOU,IB,IGAU,EPAIST,MELE,NPINT,NBGMAT,
  5. 3 NELMAT,SECT,LHOOK,TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,
  6. 4 CRIGI,DSIGT,KERRE,DT)
  7. comm ******************************************************************
  8. comm * execute le calcul annonce au niveau du p_gauss *
  9. comm * loi viscoplastique : le comportement anelastique est decrit *
  10. comm * par un mecanisme "plastique" (ecrouissage cinematique) *
  11. comm * et un mecanisme "visqueux" (Maxwell) *
  12. comm * en parallele, ce au-dela du seuil de plasticite *
  13. comm * *
  14. comm * il faut donc gerer 2 jeux de variables cachees tensorielles *
  15. comm * *
  16. comm * cette loi est derivee des travaux deja realises a propos du *
  17. comm * comportement mecanique du polyethylene *
  18. comm * (these Kichenin, Ecole Polytechnique, 1992) *
  19. comm * algorithme Nguyen Quoc Son, mise au point J Frelat *
  20. comm * $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$*
  21. comm * RECOMMANDATION : *
  22. comm * la chasse au bugs est bienvenue MAIS PAS *
  23. comm * les ameliorations ou nouveaux developpements dans *
  24. comm * l inspiration de VISK2 (rapport SEMT/LAMS 98 Kichenin, Charras)*
  25. comm * Il est preferable d une part de creer une nouvelle loi, d autre*
  26. comm * part de RECOPIER dans un AUTRE fichier les lignes voulues *
  27. comm * MERCI. JK *
  28. comm ******************************************************************
  29. IMPLICIT INTEGER(I-N)
  30. IMPLICIT REAL*8(A-H,O-Z)
  31. -INC CCOPTIO
  32. DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*),SIGF(*),
  33. & VARF(*),DEFP(*)
  34. DIMENSION RSIG0(6),RDSIGT(6),RDEFP(6),RSIGF(6),DSIGT(*)
  35. DIMENSION VALMAT(*),VALCAR(*)
  36. DIMENSION TXR(IDIM,*),CRIGI(12)
  37. DIMENSION DDAUX(LHOOK,*),DDHOMU(LHOOK,*)
  38. DIMENSION XLOC(3,3),XGLOB(3,3)
  39. DIMENSION D1HOOK(LHOOK,*),ROTHOO(LHOOK,*)
  40. CHARACTER*8 CMATE
  41. real*8 lambdav
  42. logical plast
  43. dimension sigpca(6),sigppa(6),sig2(6)
  44. c..
  45. c.. VOIR PB IFOUR IFOURB IFOU
  46. IF (IFOUR.NE.-2.AND.IFOUR.NE.-1.AND.IFOUR.NE.-3.AND.IFOUR.NE.0)
  47. & THEN
  48. iuv1 = 7
  49. ELSE
  50. iuv1 = 5
  51. ENDIF
  52. c..
  53. if((1.+xmat(2)).eq.0.or.(1-2.*xmat(2)).eq.0) then
  54. interr(5)=1
  55. kerre = 82
  56. return
  57. endif
  58. * suppose fondalement isotropie du comportement
  59. x2mu= xmat(1)/(1.+xmat(2))
  60. x2muv=xmat(8)/(1.+xmat(2))
  61. seuil=xmat(5)/1.73205080757
  62. ecrou=xmat(6)
  63. if (x2mu.le.0.or.seuil.lt.0.or.ecrou.lt.0) then
  64. interr(5) = 4
  65. kerre = 82
  66. return
  67. endif
  68. c.. initialise ecrouissage mecanisme plastique
  69. c.. l indice 1 est la deformation plastique equivalente
  70. varf(1) = var0(1)
  71. varf(2) = var0(2)
  72. varf(3) = var0(3)
  73. varf(4) = var0(4)
  74. varf(5) = var0(5)
  75. if (ifour.gt.0) then
  76. varf(6) = var0(6)
  77. varf(7) = var0(7)
  78. endif
  79. c mecanisme visqueux
  80. svxx=var0(iuv1+1)
  81. svyy=var0(iuv1+2)
  82. svzz=var0(iuv1+3)
  83. svxy=var0(iuv1+4)
  84. if (ifour.gt.0) then
  85. svxz=var0(iuv1+5)
  86. svyz=var0(iuv1+6)
  87. endif
  88. c..
  89. c
  90. C CALCUL DE L'INCREMENT DE CONTRAINTES "elastique"
  91. C
  92. CALL CALSIG(DEPST,DDAUX,NSTRS,CMATE,VALMAT,VALCAR,
  93. 1 N2EL,N2PTEL,MFR1,IFOU,IB,IGAU,EPAIST,NBPGAU,
  94. 2 MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,XLOC,
  95. 3 XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
  96. *
  97. IF(IRTD.NE.1) THEN
  98. KERRE=69
  99. GOTO 1000
  100. ENDIF
  101.  
  102.  
  103. plast = .false.
  104. c.. composante plastique + increment elastique : cinematique admissible
  105. sigpca(1)= sig0(1) - svxx + DSIGT(1)
  106. sigpca(2)= sig0(2) - svyy + DSIGT(2)
  107. sigpca(3)= sig0(3) - svzz + DSIGT(3)
  108. sigpca(4)= sig0(4) - svxy + DSIGT(4)
  109. if (ifour.gt.0)then
  110. sigpca(5)= sig0(5) - svxz + DSIGT(5)
  111. sigpca(6)= sig0(6) - svyz + DSIGT(6)
  112. endif
  113. trsigp=0.33333*(sigpca(1)+sigpca(2)+sigpca(3))
  114. c..
  115. c.. evalue cette composante (deviateur) dans le convexe definissant l'ecrouissage
  116. sig2(1)= sigpca(1)-trsigp-varf(2)
  117. sig2(2)= sigpca(2)-trsigp-varf(3)
  118. sig2(3)= sigpca(3)-trsigp-varf(4)
  119. sig2(4)= sigpca(4)-varf(5)
  120. if (ifour.gt.0)then
  121. sig2(5)= sigpca(5)-varf(6)
  122. sig2(6)= sigpca(6)-varf(7)
  123. endif
  124. c.. calcule la norme de Von Mises
  125. fsig2 = sig2(1)*sig2(1) + sig2(2)*sig2(2) + sig2(3)*sig2(3)
  126. fsig2 = fsig2 + 2.d0*sig2(4)*sig2(4)
  127. if(ifour.gt.0)then
  128. fsig2 = fsig2 + 2.d0*sig2(5)*sig2(5) + 2.d0*sig2(6)*sig2(6)
  129. endif
  130. fsig2 = SQRT(0.5d0*fsig2)
  131.  
  132.  
  133. c.. tests
  134. c a l interieur du convexe ou cas pathologique
  135. if ((fsig2-seuil).le.0.or.fsig2.eq.0) then
  136. plast=.false.
  137. coef = 0.
  138. defp(1) = 0.
  139. defp(2) = 0.
  140. defp(3) = 0.
  141. defp(4) = 0.
  142. if (ifour.gt.0) then
  143. defp(5) = 0.
  144. defp(6) = 0.
  145. endif
  146. goto 100
  147. endif
  148. c.. sinon le seuil est franchi
  149. plast = .true.
  150. c.. et il est necessaire de faire un calcul d ecoulement
  151. coef=(fsig2-seuil)/(x2mu+ecrou)/fsig2
  152. c.. increment de la deformation plastique
  153. defp(1) = coef * sig2(1)
  154. defp(2) = coef * sig2(2)
  155. defp(3) = coef * sig2(3)
  156. defp(4) = coef * sig2(4) * 2.
  157. if (ifour.gt.0) then
  158. defp(5) = coef * sig2(5) * 2.
  159. defp(6) = coef * sig2(6) * 2.
  160. endif
  161. 100 continue
  162. c.. increment de la variable d ecrouissage
  163. varf(1) = varf(1) + coef*fsig2*1.414213562
  164. varf(2) = varf(2) + ecrou*defp(1)
  165. varf(3) = varf(3) + ecrou*defp(2)
  166. varf(4) = varf(4) + ecrou*defp(3)
  167. varf(5) = varf(5) + ecrou*defp(4) / 2.
  168. if (ifour.gt.0) then
  169. varf(6) = varf(6) + ecrou*defp(5) / 2.
  170. varf(7) = varf(7) + ecrou*defp(6) / 2.
  171. endif
  172. c.. contrainte plastiquement admissible
  173. sigppa(1) = sigpca(1) - x2mu*defp(1)
  174. sigppa(2) = sigpca(2) - x2mu*defp(2)
  175. sigppa(3) = sigpca(3) - x2mu*defp(3)
  176. sigppa(4) = sigpca(4) - x2mu*defp(4) / 2.
  177. if(ifour.gt.0) then
  178. sigppa(5) = sigpca(5) - x2mu*defp(5) / 2.
  179. sigppa(6) = sigpca(6) - x2mu*defp(6) / 2.
  180. endif
  181.  
  182. fsigpa = sigppa(1)*sigppa(1) + sigppa(2)*sigppa(2)
  183. & + sigppa(3)*sigppa(3) + 2.d0*sigppa(4)*sigppa(4)
  184. if (ifour.gt.0) then
  185. fsigpa = fsigpa + 2.d0*sigppa(5)*sigppa(5)
  186. & + 2.d0*sigppa(6)*sigppa(6)
  187. endif
  188. fsigpa = SQRT(0.5d0*fsigpa)
  189. c..
  190.  
  191. 300 CONTINUE
  192.  
  193. c..
  194. c.. traite le mecanisme visqueux
  195. fsv1 = 0.
  196. comm *................................................................*
  197. comm * viscosite non-lineaire ... si devis est le parametre cache *
  198. comm * increment de deformation-visqueuse, et sigvis le parametre *
  199. comm * contrainte-visqueuse initial, la loi de comportement s ecrit *
  200. comm * eta*dev = (sigvis + loiv*(du-dev))*|sigvis+loiv*(du-dev)|**m *
  201. comm * (m = n-1) le mecanisme visqueux etant somme toute secondaire *
  202. comm * l approche bricolee consiste a evaluer *
  203. comm * fsv1 = |sigvis+loiv*du|**m, dev=0 *
  204. comm * puis ecrire sigv-final = (sigvis + loiv*(du-dev)) * coefv *
  205. comm * le coefv etant evalue a partir de fsv1 et valant *
  206. comm * coefv = 1./(1+x2muv/eta) quand m=0 *
  207. comm *................................................................*
  208. c.. on suppose le comportement visqueux homogene et isotrope
  209. lambdav=xmat(8)*xmat(2)/(1+xmat(2))/(1-2.*xmat(2))
  210. x2muv=xmat(8)/(1.+xmat(2))
  211. if(x2muv.lt.0) then
  212. interr(5)=5
  213. kerre = 82
  214. return
  215. endif
  216. m = nint(xmat(9)) - 1
  217. c.. coefficient applique a la branche visqueuse : mecanisme actif
  218. c.. au-dela du seuil, inactif sinon
  219. if (x2mu.gt.0) then
  220. xambda2 = x2muv/x2mu
  221. endif
  222. if (xambda2.le.0) xambda2 = 0.
  223. c.. viscosite
  224. if (xmat(7).le.0.) then
  225. c.. mecanisme visqueux inactif,
  226. varf(iuv1+1)=0.
  227. varf(iuv1+2)=0.
  228. varf(iuv1+3)=0.
  229. varf(iuv1+4)=0.
  230. if (ifour.gt.0)then
  231. varf(iuv1+5)=0.
  232. varf(iuv1+6)=0.
  233. endif
  234. goto 400
  235. endif
  236. c.. le mecanisme n'est efficace qu a partir d un pas de temps assez grand
  237. c.. et si l on est dans le domaine plastique
  238. if (dt.gt.1.e-15) then
  239. eta=xmat(7)/dt
  240. x2muv2=x2muv+eta
  241. coedev=1./x2muv2
  242. coesph=lambdav/x2muv2/(3.*lambdav+x2muv2)
  243. c..
  244. c.. incremente "elastiquement" la contrainte-visqueuse
  245. svxx=svxx+x2muv * defp(1)
  246. svyy=svyy+x2muv * defp(2)
  247. svzz=svzz+x2muv * defp(3)
  248. svxy=svxy+x2muv * defp(4)
  249. if(ifour.gt.0)then
  250. svxz=svxz+x2muv * defp(5)
  251. svyz=svyz+x2muv * defp(6)
  252. endif
  253. c..
  254. c.. deduit l'increment de deformation-visqueuse
  255. trsv=svxx+svyy+svzz
  256. svxx=eta*(coedev*svxx-coesph*trsv)
  257. svyy=eta*(coedev*svyy-coesph*trsv)
  258. svzz=eta*(coedev*svzz-coesph*trsv)
  259. evxy=eta*coedev*svxy
  260. if(ifour.gt.0)then
  261. svxz=eta*coedev*svxz
  262. svyz=eta*coedev*svyz
  263. endif
  264. c..
  265. fsv1 = svxx*svxx + svyy*svyy + svzz*svzz + 2.*svxy*svxy
  266. if(ifour.gt.0) then
  267. fsv1= fsv1 + 2.*svxz*svxz+2.*svyz*svyz
  268. fsv1 = fsv1/eta
  269. endif
  270. c
  271. c.. coef multiplicateur
  272. if (fsv1.gt.0.and.m.gt.0) fsv1 = (SQRT(0.5*fsv1))**m
  273. if (fsv1.gt.0.and.m.lt.0) fsv1 = (1./SQRT(0.5*fsv1))**(-1.*m)
  274. if (fsv1.gt.0.and.m.eq.0) fsv1 = 1.
  275. * if(igau.eq.4) write(6,*) ' fsv1 ', fsv1
  276. if (fsv1.eq.0) then
  277. c.. composante visqueuse nulle
  278. varf(iuv1+1) = 0.
  279. varf(iuv1+2) = 0.
  280. varf(iuv1+3) = 0.
  281. varf(iuv1+4) = 0.
  282. if (ifour.gt.0) then
  283. varf(iuv1+5) = 0.
  284. varf(iuv1+6) = 0.
  285. endif
  286. goto 400
  287. endif
  288. c.. coeff reducteur
  289. coefv = fsv1/(1. + (m+1)*x2muv*fsv1/eta)
  290. varf(iuv1+1)=svxx*coefv
  291. varf(iuv1+2)=svyy*coefv
  292. varf(iuv1+3)=svzz*coefv
  293. varf(iuv1+4)=svxy*coefv
  294. if (ifour.gt.0)then
  295. varf(iuv1+5)=svxz*coefv
  296. varf(iuv1+6)=svyz*coefv
  297. endif
  298.  
  299. c.. quand le pas de temps est petit la variable cachee est inchangee
  300. else
  301. varf(iuv1+1)=var0(iuv1+1)
  302. varf(iuv1+2)=var0(iuv1+2)
  303. varf(iuv1+3)=var0(iuv1+3)
  304. varf(iuv1+4)=var0(iuv1+4)
  305. if (ifour.gt.0)then
  306. varf(iuv1+5)=var0(iuv1+5)
  307. varf(iuv1+6)=var0(iuv1+6)
  308. endif
  309. endif
  310. comm *................................................................*
  311. comm * calcul de la contrainte finale *
  312. comm *................................................................*
  313. 400 continue
  314. sigf(1) = sigppa(1) + varf(iuv1+1)
  315. sigf(2) = sigppa(2) + varf(iuv1+2)
  316. sigf(3) = sigppa(3) + varf(iuv1+3)
  317. sigf(4) = sigppa(4) + varf(iuv1+4)
  318. if(ifour.gt.0)then
  319. sigf(5) = sigppa(5) + varf(iuv1+5)
  320. sigf(6) = sigppa(6) + varf(iuv1+6)
  321. endif
  322. return
  323.  
  324. 1000 continue
  325. end
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  

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