1 : C NEWSTA SOURCE CHAT 05/01/13 01:56:55 5004 2 : c***************************************************************************** 3 : c sous rouine utilisee dans l'ecoulement plastique de la loi 4 : c de GURSON 5 : c********************************************************************** 6 : subroutine newsta(wrkgur) 7 : c this subroutine update the material state given the strain increment e 8 : c 9 : c--- variables 10 : IMPLICIT INTEGER(I-N) 11 : IMPLICIT REAL*8 (A-H,O-Z) 12 : segment wrkgur 13 : real*8 sigbar, sy0,phi0,rho0,g,b,h 14 : real*8 epn,phin,sqrtj2,rho,sig(6) 15 : real*8 e(7),dt 16 : real*8 conv,tol1,tol2 17 : endsegment 18 : * common/prop/sigbar,sy0,phi0,rho0,g,b,h 19 : * common/state/epn,phin,sqrtj2,rho,sig(6) 20 : c e is the strain rate: 21 : c e(1-6) is the deviatoric part 22 : c e(7) is the trace/3 23 : c dt ist the time increment 24 : * common /delta/ e(7),dt 25 : * common /toler/ conv,tol1,tol2 26 : c 27 : c--- executable 28 : c 29 : c parametres pour la convergence et les tolerances d'erreur 30 : conv=1.d-4 31 : tol1=1.d-5 32 : tol2=1.d-6 33 : c 34 : c--- trial stress 35 : dum=(rho0-rho)/rho+e(7)*3.d0*dt+1.d0 36 : rho=rho0/dum 37 : sum=0.d0 38 : do 10 i=1,6 39 : sig(i)=sig(i)+2.d0*g*e(i)*dt 40 : sum=sum+sig(i)*sig(i) 41 : if (i .gt. 3) then 42 : sum=sum+sig(i)*sig(i) 43 : endif 44 : 10 continue 45 : sqrtj2=sqrt(1.5d0*sum) 46 : p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho)) 47 : syn=sy0+h*epn 48 : c 49 : c*** test if phin < phi0 a Prandtl-Reuss model is used 50 : if ((phin .le. phi0) .and. ( p .ge. 0.d0)) then 51 : if (sqrtj2 .gt. syn) then 51 : call prandt(wrkgur) 52 : endif 53 : return 54 : endif 55 : c 56 : sy=syn**2*(1.d0+phin**2-2.d0*phin*cosh(-1.5d0*p/sigbar)) 57 : yield=sqrtj2*sqrtj2-sy 58 : C 59 : c--- initialize paramter 60 : dphi=0.d0 61 : iter=0 62 : c 62 : c*** call solve if plastic evolution 63 : if (yield .gt. 0.d0) then 63 : call solve(dphi,iter,wrkgur) 64 : c 65 : c*** check if increment has been splitted in zero() 66 : c this implies phin=phi0 67 : if ( dphi .ge. 1.d0) then 68 : p=(1.d0-phi0)*b*(1.d0-rho0/rho) 69 : goto 60 70 : endif 71 : c 72 : c--- update new state value given new dphi 73 : syn=sy0+h*epn 74 : phin=phin+dphi 75 : p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho)) 76 : syp=1.d0+phin**2-2.d0*phin*cosh(-1.5d0*p/sigbar) 77 : dep=(sqrtj2-syn*sqrt(syp))/(3*g+h*sqrt(syp)) 78 : c 79 : if (abs(syp) .ge. tol1*tol1) then 80 : c a=1.5d0*phin*(1.d0-phin)*sinh(-1.5d0*p/sigbar)/sqrt(syp) 81 : c 82 : c if ( a .le. tol2) then 83 : c dep=(sqrtj2-syn*(1-phin))/(3*g+h*(1-phin)) 84 : c else 85 : c ratio=syn/(2.d0*h) 86 : c dep=sqrt(ratio*ratio+abs(dphi/a)*sigbar/h)-ratio 87 : c endif 88 : c 89 : epn=epn+dep 90 : sy=(syn+h*dep)*sqrt(syp) 91 : sqrtj2=sy 92 : c 93 : c--- return stress to yield surface 94 : fac=1.d0+3.d0*g*dep/sqrtj2 95 : do 20 i=1,6 96 : sig(i)=sig(i)/fac 97 : 20 continue 98 : else 99 : sqrtj2=0.d0 100 : endif 101 : endif 102 : c 103 : 60 continue 104 : c 105 : c*** echo new state 106 : return 107 : c 108 : end 109 : 110 : 111 : 112 : 113 : 114 : 115 : 116 : 117 : 118 :
© Cast3M 2003 - All rights reserved.
Disclaimer