Esope source files

   1 : C ZIRO      SOURCE    CHAT      05/01/13    04:23:01     5004           
   2 : c******************************************************************
   3 : c       sous routine utilisee dans l'ecoulement plastique siuvant *
   4 : c       la loi de gurson                                          *
   5 : c******************************************************************
   6 :         subroutine ziro(dphi,wrkgur)
   7 : c       this subroutine solve for phi1 the solution of
   8 : c       phi=exp(3.|P|./2/sigbar)
   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***    init data
  24 :         phimax=1.d0-rho/rho0*(1.d0-phi0)
  25 :         phi1=phimax
  26 :         p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
  27 :         f=log(phin)+1.5d0*b/sigbar*(1-phin)*(phin-phimax)/(1.d0-phimax)
  28 : c
  29 : c***    material in compression
  30 :         if ( p .gt. 0.d0) then
  31 : c
  32 : c---    test if there is a root inside the intervall
  33 :           if ( f .le. 0.d0) then
  34 :             dphi=0.d0
  35 :             return
  36 :           endif
  37 : c
  38 : c---    test if the answer will be smaller than phi0
  39 :           f=log(phi0)+1.5d0*b/sigbar*(1-phi0)
  40 :      &          *(phi0-phimax)/(1.d0-phimax)
  41 :           if (f .gt. 0.d0) then
  42 :             dphi=phi0*0.99d0-phin
  43 :             return
  44 :           else
  45 :             if (phimax .lt. phi0) phi1=phi0
  46 :           endif
  47 : c
  48 : c---    iterate the tangent approximation
  49 :           do 100 iter=1,50
  50 :           f=log(phi1)+1.5d0*b/sigbar*(1-phi1)*(phi1-phimax)
  51 :      &            /(1.d0-phimax) 
  52 : c
  53 : c---    check for convergence
  54 :           syp=sqrt(abs(1.d0+phi1*phi1-2.d0*phi1*cosh(-1.5d0*b*(1-phi1)*
  55 :      &         (phi1-phimax)/(1.d0-phimax)/sigbar)))
  56 : c
  57 :           if ( syp .le. tol1) then
  58 :             dphi=phi1-phin
  59 :             return
  60 :           endif
  61 : c---
  62 :           fprime=1.d0/phi1+1.5d0*b/sigbar*(1.d0+phimax-2.d0*phi1)/
  63 :      &          (1.d0-phimax)
  64 :           phi1=phi1-f/fprime
  65 : c
  66 :   100     continue
  67 : c
  68 : c***    material in tension
  69 :         else
  70 :           a1=1+phimax
  71 : c
  72 : c---    loop over the process
  73 :           do 200 iter=1,90
  74 : c
  75 :           a2=2.d0/3.d0*sigbar/b*(1.d0-phimax)*log(phi1)+phimax
  76 :           phi1=(a1-sqrt(a1*a1-4.d0*a2))*0.5d0
  77 : c
  78 : c---    test if the search is outside the interval
  79 :           if ( phi1 .le. phin) then
  80 :             dphi=0.d0
  81 :             return
  82 :           endif
  83 : c
  84 : c---    check for convergence 
  85 :           syp=sqrt(abs(1.d0+phi1*phi1-2.d0*phi1*cosh(-1.5d0*b*(1-phi1)*
  86 :      &         (phi1-phimax)/(1.d0-phimax)/sigbar)))
  87 : c
  88 :           if (syp .le. tol1) then 
  89 :             dphi=phi1-phin 
  90 :             return
  91 :           endif  
  92 : c
  93 :   200     continue
  94 : c
  95 :         endif
  96 :         dphi=phi1-phin
  97 : c       the subroutine didnot converge return anyway
  98 :         return
  99 : c
 100 :         end
 101 : 
 102 : 
 103 :  
 104 :