Télécharger newsta.eso

Retour à la liste

Numérotation des lignes :

  1. C NEWSTA SOURCE PV 15/04/13 21:15:15 8474
  2. c*****************************************************************************
  3. c sous rouine utilisee dans l'ecoulement plastique de la loi
  4. c de GURSON
  5. c**********************************************************************
  6. subroutine newsta(wrkgur)
  7. c this subroutine update the material state given the strain increment e
  8. c
  9. c--- variables
  10. IMPLICIT INTEGER(I-N)
  11. IMPLICIT REAL*8 (A-H,O-Z)
  12. segment wrkgur
  13. real*8 sigbar, sy0,phi0,rho0,g,b,h
  14. real*8 epn,phin,sqrtj2,rho,sig(6)
  15. real*8 e(7),dt
  16. real*8 conv,tol1,tol2
  17. endsegment
  18. * common/prop/sigbar,sy0,phi0,rho0,g,b,h
  19. * common/state/epn,phin,sqrtj2,rho,sig(6)
  20. c e is the strain rate:
  21. c e(1-6) is the deviatoric part
  22. c e(7) is the trace/3
  23. c dt ist the time increment
  24. * common /delta/ e(7),dt
  25. * common /toler/ conv,tol1,tol2
  26. c
  27. c--- executable
  28. c
  29. c parametres pour la convergence et les tolerances d'erreur
  30. conv=1.d-4
  31. tol1=1.d-5
  32. tol2=1.d-6
  33. c
  34. c--- trial stress
  35. dum=(rho0-rho)/rho+e(7)*3.d0*dt+1.d0
  36. rho=rho0/dum
  37. sum=0.d0
  38. do 10 i=1,6
  39. sig(i)=sig(i)+2.d0*g*e(i)*dt
  40. sum=sum+sig(i)*sig(i)
  41. if (i .gt. 3) then
  42. sum=sum+sig(i)*sig(i)
  43. endif
  44. 10 continue
  45. sqrtj2=sqrt(1.5d0*sum)
  46. p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
  47. syn=sy0+h*epn
  48. c
  49. c*** test if phin < phi0 a Prandtl-Reuss model is used
  50. if ((phin .le. phi0) .and. ( p .ge. 0.d0)) then
  51. if (sqrtj2 .gt. syn) then
  52. call prandt(wrkgur)
  53. endif
  54. return
  55. endif
  56. c
  57. syach=-1.5d0*p/sigbar
  58. * borner pour eviter un depassement
  59. syach=min(675d0,max(-675d0,syach))
  60. sy=syn**2*(1.d0+phin**2-2.d0*phin*cosh(syach))
  61. yield=sqrtj2*sqrtj2-sy
  62. C
  63. c--- initialize paramter
  64. dphi=0.d0
  65. iter=0
  66. c
  67. c*** call solve if plastic evolution
  68. if (yield .gt. 0.d0) then
  69. call solve(dphi,iter,wrkgur)
  70. c
  71. c*** check if increment has been splitted in zero()
  72. c this implies phin=phi0
  73. if ( dphi .ge. 1.d0) then
  74. p=(1.d0-phi0)*b*(1.d0-rho0/rho)
  75. goto 60
  76. endif
  77. c
  78. c--- update new state value given new dphi
  79. syn=sy0+h*epn
  80. phin=phin+dphi
  81. p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
  82. syp=1.d0+phin**2-2.d0*phin*cosh(-1.5d0*p/sigbar)
  83. dep=(sqrtj2-syn*sqrt(syp))/(3*g+h*sqrt(syp))
  84. c
  85. if (abs(syp) .ge. tol1*tol1) then
  86. c a=1.5d0*phin*(1.d0-phin)*sinh(-1.5d0*p/sigbar)/sqrt(syp)
  87. c
  88. c if ( a .le. tol2) then
  89. c dep=(sqrtj2-syn*(1-phin))/(3*g+h*(1-phin))
  90. c else
  91. c ratio=syn/(2.d0*h)
  92. c dep=sqrt(ratio*ratio+abs(dphi/a)*sigbar/h)-ratio
  93. c endif
  94. c
  95. epn=epn+dep
  96. sy=(syn+h*dep)*sqrt(syp)
  97. sqrtj2=sy
  98. c
  99. c--- return stress to yield surface
  100. fac=1.d0+3.d0*g*dep/sqrtj2
  101. do 20 i=1,6
  102. sig(i)=sig(i)/fac
  103. 20 continue
  104. else
  105. sqrtj2=0.d0
  106. endif
  107. endif
  108. c
  109. 60 continue
  110. c
  111. c*** echo new state
  112. return
  113. c
  114. end
  115.  
  116.  
  117.  
  118.  
  119.  
  120.  
  121.  
  122.  
  123.  
  124.  
  125.  
  126.  

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