代码

subroutine vusdfld(
 c Read only -
     *   nblock, nstatev, nfieldv, nprops, ndir, nshr,
     *   jElem, kIntPt, kLayer, kSecPt,
     *   stepTime, totalTime, dt, cmname,
     *   coordMp, direct, T, charLength, props,
     *   stateOld,
c Write only -
     *   stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension props(nprops),
     *          jElem(nblock), coordMp(nblock, *),
     *          direct(nblock, 3, 3), T(nblock,3,3),
     *          charLength(nblock),
     *          stateOld(nblock, nstatev),
     *          stateNew(nblock, nstatev),
     *          field(nblock, nfieldv)
      character*80 cmname
C
C Local arrays from vgetvrm are dimensioned to
C maximum block size (maxblk)
C
      parameter( nrData=6 )
      character*3 cData(maxblk*nrData)
      dimension rData(maxblk*nrData), jData(maxblk*nrData)
      dimension StatusOld(NBLOCK,1),StatusNew(NBLOCK,1)
     
      PARAMETER(ZERO=0.0,ONE=1.0,TWO=2.0,THREE=3.0)
      DIMENSION STRESSNEW(NDIR+NSHR)    
C
C     INPUT CONSTANT PARAMETERS
C
        EMOD=80000
        ENU=0.33
        DC1=5.0992
        DC2=0.7247
        DC3=0.3611
        EA=301.0
        EQ=205.0
        EB=18.0
       
        EBULK3=EMOD/(ONE-TWO*ENU)
        EG2=EMOD/(ONE+ENU)
        EG=EG2/TWO
        EG3=THREE*EG
        ELAM=(EBULK3-EG2)/THREE
      
      DO K0=1,NBLOCK
       
C
C  INITIAL FAILURE FLAGS FROM STATE VARIABLES
C
       DAMAGE=stateOld(K0,2)
       CALL VGETVRM('S',RDATA,JDATA,CDATA,JSTATUS)
         STRESSNEW(1)=RDATA(K0*1)
         STRESSNEW(2)=RDATA(K0*2)
         STRESSNEW(3)=RDATA(K0*3)
         STRESSNEW(4)=RDATA(K0*4)
         IF(NSHR.GT.1) THEN
         STRESSNEW(5)=RDATA(K0*5) 
         STRESSNEW(6)=RDATA(K0*6)
         END IF
      
       CALL VGETVRM('PEEQ',RDATA,JDATA,CDATA,JSTATUS)
       EQPLAS=RDATA(K0)
      
       IF(NSHR.EQ.1) THEN
         POTENIAL=0.50*((STRESSNEW(1)-STRESSNEW(2))**2+
     1   (STRESSNEW(2)-STRESSNEW(3))**2+(STRESSNEW(3)-STRESSNEW(1))**2
     2   +6*STRESSNEW(4)**2)
       ELSE
         POTENIAL=0.50*((STRESSNEW(1)-STRESSNEW(2))**2+
     1   (STRESSNEW(2)-STRESSNEW(3))**2+(STRESSNEW(3)-STRESSNEW(1))**2
     2   +6*(STRESSNEW(4)**2+STRESSNEW(5)**2+STRESSNEW(6)**2))
       END IF
         EFS=SQRT(POTENIAL)
      
         HARD=0.0
       IF(EQPLAS.EQ.0.0) THEN
         SYIELD=EA
       ELSE
         HARD=EQ*EB*EXP(-1*EB*EQPLAS)
         SYIELD=EA+EQ*(1-EXP(-1*EB*EQPLAS))
       END IF
      
         RHS=EFS-SYIELD
C
C  JUDGE IF YIELD
C
       IF(RHS.GT.0.0) THEN
         DEQPLAS=RHS/(EG3+HARD)
        
         STRESSM=(STRESSNEW(1)+STRESSNEW(2)+STRESSNEW(3))/3
       
        HJ2=((STRESSNEW(1)-STRESSNEW(2))**2+(STRESSNEW(2)-STRESSNEW(3))
     1      **2+(STRESSNEW(3)-STRESSNEW(1))**2+6.0*STRESSNEW(4)**2)/6.0
        HJ3=(STRESSNEW(1)-STRESSM)*(STRESSNEW(2)-STRESSM)*
     1    (STRESSNEW(3)-STRESSM)-(STRESSNEW(3)-STRESSM)*STRESSNEW(4)**2
        THETA=ACOS(3.0**1.5*HJ3/(2.0*HJ2**1.50))/3
       
        CR11=COS(THETA)
        CR12=SIN(THETA)
        CR21=-1.0*SIN(THETA)
        CR22=COS(THETA)
       
        SP1=CR11**2*STRESSNEW(1)+CR21**2*STRESSNEW(2)+
     1      2*CR11*CR21*STRESSNEW(4)
        SP2=CR12**2*STRESSNEW(1)+CR22**2*STRESSNEW(2)+
     1      2*CR12*CR22*STRESSNEW(4)
        SP3=STRESSNEW(3)
       
        HYSTRESS=(SP1+SP2+SP3)/3.0
     
        TNH=HYSTRESS/EFS
       
        SP1=MAX(SP1,SP2,SP3)
        SP3=MIN(SP1,SP2,SP3)
       
        DDAMAGE=(((SP1-SP3)/EFS)**DC1+(TNH-1.0/3.0))**DC2*DEQPLAS
     1        /DC3
       
      
       END IF
      
       DAMAGE=DAMAGE+DDAMAGE
       stateNew(K0,2)=DAMAGE
       
       IF(stateNew(K0,2).GE.1) THEN
          stateNew(K0,1)=0.0
       ELSE
          stateNew(K0,1)=1.0
       END IF
        
     
      END DO  
      RETURN
      END