Numérotation des lignes :

C NEWSTA    SOURCE    PV        15/04/13    21:15:15     8474c*****************************************************************************c       sous rouine utilisee dans l'ecoulement plastique de la loic       de GURSONc**********************************************************************        subroutine newsta(wrkgur)c       this subroutine update the material state given the strain increment ecc---    variables      IMPLICIT INTEGER(I-N)        IMPLICIT REAL*8 (A-H,O-Z)        segment wrkgur         real*8 sigbar, sy0,phi0,rho0,g,b,h         real*8 epn,phin,sqrtj2,rho,sig(6)         real*8 e(7),dt         real*8 conv,tol1,tol2        endsegment*        common/prop/sigbar,sy0,phi0,rho0,g,b,h*        common/state/epn,phin,sqrtj2,rho,sig(6)c               e is the strain rate:c                       e(1-6) is the deviatoric partc                       e(7) is the trace/3c               dt ist the time increment*        common /delta/ e(7),dt*       common /toler/ conv,tol1,tol2cc---    executablecc       parametres pour la convergence et les tolerances d'erreur        conv=1.d-4        tol1=1.d-5        tol2=1.d-6cc---    trial stress        dum=(rho0-rho)/rho+e(7)*3.d0*dt+1.d0        rho=rho0/dum        sum=0.d0        do 10 i=1,6             sig(i)=sig(i)+2.d0*g*e(i)*dt             sum=sum+sig(i)*sig(i)             if (i .gt. 3) then               sum=sum+sig(i)*sig(i)             endif   10   continue        sqrtj2=sqrt(1.5d0*sum)        p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))        syn=sy0+h*epncc***    test if phin &lt; phi0 a Prandtl-Reuss model is used        if ((phin .le. phi0) .and. ( p .ge. 0.d0)) then          if (sqrtj2 .gt. syn) then            call prandt(wrkgur)          endif          return        endifc        syach=-1.5d0*p/sigbar*  borner pour eviter un depassement        syach=min(675d0,max(-675d0,syach))        sy=syn**2*(1.d0+phin**2-2.d0*phin*cosh(syach))        yield=sqrtj2*sqrtj2-syCc---    initialize paramter        dphi=0.d0        iter=0cc***    call solve if plastic evolution        if (yield .gt. 0.d0) then          call solve(dphi,iter,wrkgur)cc***    check if increment has been splitted in zero()c               this implies phin=phi0          if ( dphi .ge. 1.d0)  then            p=(1.d0-phi0)*b*(1.d0-rho0/rho)            goto 60          endifcc---    update new state value given new dphi          syn=sy0+h*epn          phin=phin+dphi          p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))          syp=1.d0+phin**2-2.d0*phin*cosh(-1.5d0*p/sigbar)          dep=(sqrtj2-syn*sqrt(syp))/(3*g+h*sqrt(syp))c          if (abs(syp) .ge. tol1*tol1) thenc           a=1.5d0*phin*(1.d0-phin)*sinh(-1.5d0*p/sigbar)/sqrt(syp)cc           if ( a .le. tol2) thenc               dep=(sqrtj2-syn*(1-phin))/(3*g+h*(1-phin))c           elsec                ratio=syn/(2.d0*h)c                dep=sqrt(ratio*ratio+abs(dphi/a)*sigbar/h)-ratioc           endifc            epn=epn+dep            sy=(syn+h*dep)*sqrt(syp)            sqrtj2=sycc---        return stress to yield surface            fac=1.d0+3.d0*g*dep/sqrtj2            do 20 i=1,6              sig(i)=sig(i)/fac   20       continue          else            sqrtj2=0.d0          endif        endifc  60    continuecc***    echo new state        returnc        end

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