Télécharger la plaquette Cast3M


fr en

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