Esope source files

   1 : C SPLIT     SOURCE    PV        08/09/23    21:15:07     6164           
   2 : c*****************************************************************************
   3 : c       sous rouine utilisee dans l'ecoulement plastique de la loi
   4 : c       de GURSON 
   5 : c*****************************************************************************
   6 :         subroutine split(dphi,wrkgur)
   7 : c
   8 : c       this subroutine split up the strain increment in two if a normal
   9 : c       evolution reaches a state where phin < phi0     
  10 : c       xi*de is the first strain increment that makes reaches phin=phi0
  11 : c
  12 : c       rho and sqrtj2 have been updated in newsta
  13 : c
  14 : c---    variables
  15 :       IMPLICIT INTEGER(I-N)
  16 :         IMPLICIT REAL*8 (A-H,O-Z)
  17 :         segment wrkgur
  18 :          real*8 sigbar, sy0,phi0,rho0,g,b,h
  19 :          real*8 epn,phin,sqrtj2,rho,sig(6)
  20 :          real*8 e(7),dt
  21 :          real*8 conv,tol1,tol2
  22 :         endsegment
  23 : *        common/prop/sigbar,sy0,phi0,rho0,g,b,h
  24 : *        common/state/epn,phin,sqrtj2,rho,sig(6)
  25 : c               e is the strain rate:
  26 : c                       e(1-6) is the deviatoric part
  27 : c                       e(7) is the trace/3
  28 : c               dt ist the time increment
  29 : *        common /delta/ e(7),dt
  30 : *       common /toler/ conv,tol1,tol2
  31 : c
  32 : c***    set the bounds
  33 :         p=(1.d0-phi0)*b*(1.d0-rho0/rho)
  34 :         syp=1.d0+phi0*phi0-2.d0*phi0*cosh(1.5d0*p/sigbar)
  35 : c
  36 :         if ( syp .le. 0.d0) then
  37 : c         xi2 is computed such that syp(xi2)=0
  38 :           xi2=1.d0+(2.d0/3.d0*sigbar*log(phi0)+(1.d0-phi0)*b*
  39 :      &       (1.d0-rho0/rho))
  40 :      &       /(1.d0-phi0)/b/e(7)/dt/3.d0
  41 :         else
  42 :           xi2=1.d0
  43 :         endif
  44 :         xi1=0.d0
  45 : c
  45 :         call funk(xi1,f1,wrkgur)
  45 :         call funk(xi2,f2,wrkgur)
  46 : c
  47 :         if ( abs(f2) .le. conv) then
  48 :           xi=xi2
  49 :           goto 200
  50 :         endif
  51 :         if ( abs(f1) .le. conv) then
  52 :           xi=xi1
  53 :           goto 200
  54 :         endif
  55 : c
  56 :         if ( (f2*f1) .gt. 0.d0) then
  57 :           write (6,*)  'f2*f1>0 in split'
  58 :         endif
  59 : c
  60 : c***    start iterative process (Ridder's method)       
  61 :         do 100 j=1,99
  62 : c
  63 :           xi3=(xi1+xi2)*0.5d0
  63 :           call funk (xi3,f3,wrkgur)
  64 : c
  65 :           dum=xi3+(xi3-xi1)*sign(1.d0,f1-f2)*f3/sqrt(abs(f3**2-f2*f1))
  65 :           call funk(dum,fdum,wrkgur)
  66 : c
  67 :           if (fdum .lt. 0.d0) then
  68 :             xi2=dum
  69 :             f2=fdum
  70 :           else
  71 :             xi1=dum
  72 :             f1=fdum 
  73 :           endif
  74 : c
  75 :           if ( abs(fdum) .le. conv) then
  76 :             xi=dum
  77 :             goto 200
  78 :           endif
  79 : c
  80 :   100   continue
  81 : c
  82 : c       split did not converge!'
  83 :         xi=xi1
  84 : c
  85 :   200   continue
  86 : c
  87 : c***    update new state for xi*de strain increment
  88 :         dphi=phi0-phin
  89 :         phin=phi0
  90 :         p=(1.d0-phi0)*b*(1.d0-rho0/rho)+(1-xi)*
  91 :      &     (1.d0-phi0)*b*e(7)*dt*3.d0
  92 :         syn=sy0+h*epn
  93 :         syp=1.d0+phi0**2-2.d0*phi0*cosh(-1.5d0*p/sigbar)
  94 : c
  95 :         if (abs(syp) .gt. tol1*tol1) then
  96 :           a=1.5d0*phi0*(1.d0-phi0)*sinh(-1.5d0*p/sigbar)
  97 :      &          /sqrt(syp)
  98 : c
  99 :           ratio=syn/(2.d0*h)
 100 :           dep=sqrt(ratio**2+abs(dphi/a)*sigbar/h)-ratio
 101 : c
 102 :           epn=epn+dep
 103 :           sy=(syn+h*dep)*sqrt(syp)
 104 :           sqrtj2=sy
 105 : c
 106 :           fac=1.d0+3.d0*g*dep/sqrtj2
 107 : c
 108 :           do 250 i=1,6
 109 :           sig(i)=(sig(i)-(1.d0-xi)*2.d0*g*e(i)*dt)/fac
 110 :   250     continue
 111 : c
 112 :         else
 113 :           sqrtj2=0.d0
 114 :         endif
 115 : c
 116 : c***    evolution under (1-xi)*de
 117 :         sum=0.d0
 118 : c
 119 :         do 400 i=1,6
 120 :           e(i)=(1.d0-xi)*e(i)
 121 :           sig(i)=sig(i)+2*g*e(i)*dt
 122 :           sum=sum+sig(i)*sig(i)
 123 :           if (i .gt. 3) then
 124 :             sum=sum+sig(i)*sig(i)
 125 :           endif
 126 :   400   continue
 127 :         sqrtj2=sqrt(1.5d0*sum)
 128 : c
 129 :         syn=sy0+h*epn
 130 : c
 131 : c---    trial state with respect to yield surface
 132 :         if ( sqrtj2 .gt. syn )then
 132 :           call prandt(wrkgur)
 133 :         endif
 134 : c
 135 :         dphi=1.0d0
 136 :         return
 137 :         end
 138 : 
 139 : 
 140 : 
 141 : 
 142 :  
 143 :  
 144 :