Télécharger solve.eso

Retour à la liste

Numérotation des lignes :

solve
  1. C SOLVE SOURCE CHAT 05/01/13 03:21:55 5004
  2. c*****************************************************************************
  3. c sous rouine utilisee dans l'ecoulement plastique de la loi
  4. c de GURSON
  5. c**********************************************************************
  6. subroutine solve(dphik,iter,wrkgur)
  7. c this subroutine solve the non linear equation which solution is
  8. c the phi increment
  9. c
  10. c--- variables
  11. IMPLICIT INTEGER(I-N)
  12. IMPLICIT REAL*8 (A-H,O-Z)
  13. segment wrkgur
  14. real*8 sigbar, sy0,phi0,rho0,g,b,h
  15. real*8 epn,phin,sqrtj2,rho,sig(6)
  16. real*8 e(7),dt
  17. real*8 conv,tol1,tol2
  18. endsegment
  19. * common/prop/sigbar,sy0,phi0,rho0,g,b,h
  20. * common/state/epn,phin,sqrtj2,rho,sig(6)
  21. * common /toler/ conv,tol1,tol2
  22. c
  23. c*** initial data
  24. p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
  25. phimax=1.d0-rho/rho0*(1.d0-phi0)
  26. if ( phimax .lt. 0.d0) phimax=0.d0
  27. c
  28. c*** find the first dphi for which syp>0
  29. call ziro(dphik,wrkgur)
  30. phik=phin+dphik
  31. c
  32. if ( phik .lt. phi0) then
  33. call split(dphik,wrkgur)
  34. return
  35. endif
  36. c
  37. c--- initialize the search
  38. call func(fk,dphik,wrkgur)
  39. iter=1
  40. c
  41. if ( abs(fk) .le. conv) return
  42. c
  43. if ( fk .gt. 0.d0) then
  44. return
  45. endif
  46. c
  47. c*******perform iterations until find fk>0
  48. if ( phimax .gt. 0.d0) then
  49. c
  50. c--- compute the constant used to determine dphi
  51. c=sqrt(abs(phimax-phin)/phimax/h/b)*2.d0*sigbar/3.d0*
  52. & (h*(1-phimax)+3.*g)
  53. c=c/sy0/2.d0
  54. c--- iterates until converges
  55. do 100 i=1,50
  56. c
  57. c check for convergence
  58. if (abs(fk).le. conv) return
  59. c
  60. c check if fk>0 and then exit loop
  61. if (fk .gt. 0.d0) goto 200
  62. c
  63. c check if split is needed
  64. if ( phik .lt. phi0) then
  65. call split(dphik,wrkgur)
  66. return
  67. endif
  68. c
  69. c--- keep track of the previous step
  70. fj=fk
  71. dphij=dphik
  72. c
  73. iter=iter+1
  74. c
  75. c--- find new value of dphik using the interpolation function
  76. d=fk-c*sqrt(1.d0/abs(phimax-phik))
  77. phik=phimax+sign(1.d0,p)*(c/d)**2
  78. dphik=phik-phin
  79. call func(fk,dphik,wrkgur)
  80. c
  81. 100 continue
  82. c
  83. c*******case when phimax is negative
  84. else
  85. c
  86. c--- initialize data for interpolation function
  87. p0=b*(1-rho0/rho/(1-phi0))
  88. c2=(h+3.d0*g)*sqrt(phin*sigbar/1.5d0/h
  89. & /sinh(1.5d0*p0/sigbar))
  90. & /1.d1
  91. c
  92. c--- iterates til converges
  93. do 150 i=1,50
  94. c
  95. c check for convergence
  96. if (abs(fk).le. conv) return
  97. c
  98. c check if fk>0 and then exit loop
  99. if (fk .gt. 0.d0) goto 200
  100. c
  101. c check if split is needed
  102. if ( phik .lt. phi0) then
  103. call split(dphik,wrkgur)
  104. return
  105. endif
  106. c
  107. c--- keep track of the previous step
  108. fj=fk
  109. dphij=dphik
  110. iter=iter+1
  111. c
  112. c--- find new value of dphik using the interpolation function
  113. d=fk-c2*sqrt(1.d0/phik)
  114. phik=sign(1.d0,p)*(c2/d)**2
  115. dphik=phik-phin
  116. call func(fk,dphik,wrkgur)
  117. c
  118. 150 continue
  119. endif
  120. c
  121. c--- the previous algorithm didn't converge
  122. return
  123. 200 continue
  124. c
  125. c*******find zero using Ridders method
  126. c see Numerical recipee in Fortran p351
  127. do 300 i=1,200
  128. iter=iter+1
  129. c
  130. half=(dphik+dphij)/2
  131. call func(fm,half,wrkgur)
  132. s=sqrt(abs(fm**2-fk*fj))
  133. dum=half+(half-dphij)*(sign(1.d0,fj-fk)*fm/s)
  134. c
  135. c--- new interval
  136. call func(fdum,dum,wrkgur)
  137. if (fdum .lt. 0.d0) then
  138. dphij=dum
  139. fj=fdum
  140. else
  141. dphik=dum
  142. fk=fdum
  143. endif
  144. c
  145. if ( phin+dphij .lt. phi0) then
  146. call split(dphik,wrkgur)
  147. return
  148. endif
  149. c
  150. c--- check for convergence
  151. if (abs(fdum) .lt. conv) then
  152. dphik=dum
  153. if ( (phin+dphik) .le. phi0) then
  154. call split(dphik,wrkgur)
  155. return
  156. endif
  157. return
  158. endif
  159. c
  160. 300 continue
  161. return
  162. c
  163. end
  164. c*************************************************************
  165.  
  166.  
  167.  
  168.  
  169.  
  170.  
  171.  
  172.  
  173.  
  174.  

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