User login

Navigation

You are here

Sharing an abaqus UMAT using an elastic viscoplastic model for soil

Here is a UMAT using an  elastic viscoplastic model for soil

Modifed Cam-clay model coupled with viscosity.

 

ABAQUS - UMAT (Free Download for Research Purposes)


A user-defined subroutine in ABAQUS for viscous soils.

Introduction:

The ABAQUS subroutine is for viscous soils, which
usually have the following characteristics.

  1. Rate-sensitivity. The stress-strain behaviour of
    viscous soils is rate-dependent. Specifically, Two important engineering
    parameters, i.e., undrained shear strength and preconsolidation pressure
    , typically increase with loading rate.

  2. Time-dependence, including creep and
    relaxation.

A general 3D elasto-viscoplastic (EVP) constitutive
model has been implemented in the user-defined subroutine in ABAQUS. The
objective of developing the subroutine is to make users easily take
account of soil’s viscous behaviours.

The report of AN ABAQUS USER SUBROUTINE INCORPORATING
ELASTIC-VISCOPLASTIC CONSTITUTIVE MODEL here provides more detailed
information about the subroutine and its application. This report first
reviews the theoretical background the EVP constitutive model. The next
section introduces the numerical solution algorithm adopted in the
subroutine. The following section illustrates how to determine the input
parameters according the given viscous soil. Finally, the appendix of
this report provides a series of verification numerical tests in a
format of ABAQUS Benchmark problem and evaluates the obtained numerical
results in reference to the corresponding theoretical solution.

  • Benchmark problem 1 - CIU compression tests on NC
    soil at various strain rates

  • Benchmark problem 2 - CIU
    compression tests on OC soil at various strain rates

  • Benchmark problem 3 - CIU
    compression tests on NC soil at step-changed strain rates

  • Benchmark problem 4 – CID (Drained) compression tests
    at various strain rates

  • Benchmark problem 5 - ISO compression creep test

  • Benchmark problem 6 – Undrained creep test

  • Benchmark problem 7 – A case study on a retaining
    wall in soft clay.

This user-defined subroutine is now free to download
for research purposes.

UMAT Code

https://docs.google.com/file/d/0B8z7TOlNv4UGMzJjOTVjOTUtOTJjZi00NzAyLTg2MjktMDlmMDA0Njg3ODE4/edit?pli=1

Memo: AN ABAQUS USER SUBROUTINE INCORPORATING ELASTIC-VISCOPLASTIC CONSTITUTIVE MODEL

https://drive.google.com/file/d/0B8z7TOlNv4UGYzEwNDI4NmMtY2I5OS00NTIyLTk...

   Benchmark problems

 

 https://drive.google.com/file/d/0B8z7TOlNv4UGMTFjYjY0NDItOGQ2My00NThmLTl...

 

 https://drive.google.com/file/d/0B8z7TOlNv4UGZGY0MTNlZTItZWI2OS00ZmY4LWE...

https://drive.google.com/file/d/0B8z7TOlNv4UGODU2MjE4ZGYtZDFhMC00Njk3LWE...

https://drive.google.com/file/d/0B8z7TOlNv4UGN2VmYTkyY2UtZTQxZC00MmJjLTg...

https://drive.google.com/file/d/0B8z7TOlNv4UGYjE4NzI2NzctOGFjNC00YzkzLTk...

 

 

For inquiries, please contact:

Dr. G. Qu


E-mail:
guangfengqu@hotmail.com

 

 

Comments

Dear Dr. G. Qu,

I am very interested in how to include viscoplastic in umat. The links you provided in this post do not work, could you please fix this? Thank you so much.

Best regards,

Jingfen

Here is the UMAT

----------------------

!  UMAT VISCO.f90
!

!****************************************************************************
!
!  PROGRAM:    UMAT FOR ABAQUS 6.7.2
!  OBJECTIVE:    TO IMPLEMENT AN OVERSTRESS TYPE OF ELASTO VISCOPLASTICITY MODEL                    APRIL/2008
!  BY:        GUANGFENG QU                    APRIL/21 008
!
!****************************************************************************
C    Changes:
!    2008 DEC 09  ADD an option to allow set an elastic state for GEOstatic step to avoid inconvergence in case of initial stress > PCSY, using a switch variavle IGEOSTEP,PROPS(18)
!                    PROPS(18)=0  default value, no special measure is taken.
!                     PROPS(18)=1  the special measure is invoken.
                    ! UMAT200812b.for
!    2008 Dec 01  Correct a bug in the rewritten code in 2008 Nov 28 ! UMAT200812.for
!    2008 Nov 28  Rewrite the code of Iyield Zones. The detailed flow chart was recorded in the note. UMAT200811.for
!    2008 Nov 27  DILATANCY    = 0.2
!    2008 Nov 25. Add an option to use "PARALLEL APPROACH" in the O/C zone to determine Overstress, PCDY, and plastic potential.
!            Set TALOC as "-1" to activate this option. IF TALOC > 0, then take the O/C limit approach
!                1. add the subroutine KOCPARALLEL for determine  PCDY using "PARALLEL APPROACH"
!                2. Main program: a) Plastic potential(dilation based on M,i.e. TALNC); b) add option to determine Iyield based on "PARALLEL APPROACH"

!    2008 NOV.  1. Add the subroutine (PCSY0PROFILE) to assign a profile of static yield surface (5 points max.).
!            INPUT: Write "-1" for the CAE interaction input of Static yield stress
!            Main Pragram: Rewrite codes to call this Sub.
!        2. Correct the part obtaining Overstress. Considering "The CASE Where the stress lies ahead the centre of static surface and above the Critical State Line" in the main program
!        3. Add a check when the mean stress became negative in the main program

C    NOTE:
C    1.    UMAT FOR STRESS-DEPENDENT ELASTICITY AND DP/CAP PLASTICITY
C    2.    CANNOT BE USED FOR PLANE STRESS
C    3.   IN THIS SUBROUTINE, COMPRESSION STRESS AND STRAIN IS DEFINED AS POSITIVE
C        BE CONSISTENT WITH THE COMMON DEFINITION IN SOIL MECHANICS, THUS, A TRANSLATION
C        AT THE BEGINING AND THE END IS PROVIDED TO SERVE THIS OBJECTIVE.

    SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)

C    LOCAL ARRAY
      DIMENSION    DE(NTENS,NTENS),DVPR(NTENS),C(NTENS,NTENS),
     1            DC(NTENS,NTENS),H(NTENS,NTENS),DEP(NTENS,NTENS),
     2            DEP1(NTENS,NTENS),DVPRDSTR(NTENS,NTENS)
!      DIMENSION DE(4,4),DVPR(4),C(4,4),DC(4,4),H(4,4),DEP(4,4),DEP1(4,4)
C
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=30.D0,TOLER=1.0D-6)
C SMEANMIN:    THE LOW BOUND OF MEAN STRESS TO CALCULATE 
C            THE STRESS-DEPENDENT ELASTIC MODULUS
C

    COMMON /KPROPS1/                        !ALL IN PROPS()
     1    ENU, RKAP,RLAM,    VOIDS, TALOC, TALNC, RCNC, ASPCT, PCSY0,        ! 1-9
     2    GAMA,RSEN,YITA, THETA                !10-13
    COMMON /KPROPS2/ NLGEOM,ITRATN,ITRMAX !14-17
    COMMON /TEMPPOSTRESS/POSTRESS
    ! temperate data only for out put the accurate stress position

C
C...    TRANFER STRESS AND STRAIN TO CONVENTIONAL DEFINITION IN SOIL MECHANICS.
C            (COMPRESSION IS POSITIVE)
    DO N=1,NTENS
    STRESS(N) = -STRESS(N)
    STRAN(N)  = -STRAN(N)
    DSTRAN(N) = -DSTRAN(N)
    ENDDO

C ----------------------------------------------------------------
C    INITIALIZE THE PARAMETERS IN COMMEN ZONE   
C ----------------------------------------------------------------
C....THE PROPS                                                    
C      ENU  =PROPS(1) POISSON'S RATIO                             
C      RKAP =PROPS(2) Slope of swelling curve.                    
C      RLAM =PROPS(3) SLOPE OF COMPRESSION CURVE                  
C      VOIDS=PROPS(4) voids ratio                                 
C    TALOC=PROPS(5)= tan(alpha) slope of Drucker Prager failure   
C             envelope in the O/C stress range.                   
C    TALNC=PROPS(6) = tan(alpha) slope of Drucker Prager failure  
C             envelope in the N/C stress range.                   
C    RCNC =PROPS(7) = intercept of Drucker Prager failure envelope
C             with SR2J2 axis in N/C stress range.                
C    ASPCT=PROPS(8)  = Aspect ratio of elliptical yield surface (X/Y)
C                        OR = LATERAL RADIUS / VERTICAL RADIUS .  
                                                                  
C....Viscoplastic Parameters:                                     
C    PCSY0=PROPS(9)    = Position of "INITIAL" static yield surface          
C    RAMA=PROPS(10)   = fluidity constant                         
C    RSEN=PROPS(11)   = strain-rate SENSITIVITY PAR.              

C....Viscoplastic VARIABLES:                                 
C    D      = measure of overstress.
C    DVPR() = current incremental viscoplastic strain rates
C    DT     = time increment
C    YRAT   = current yield ratio       

    ENU        =    PROPS(1)   
    RKAP    =    PROPS(2)
       
    RLAM    =    PROPS(3)
    VOIDS    =    PROPS(4)

    TALOC    =    PROPS(5)  !IF the input value is negative, the O/C Parallel Approach will be invoked; otherwise, use the OC limit approach
    TALNC    =    PROPS(6) 
    RCNC    =    PROPS(7) 
                             
    ASPCT    =    PROPS(8) 
                             
    PCSY0    =    PROPS(9)  !IF the input value is negative, the PCSY0 will be deterimed by a subrountine
    GAMA    =    PROPS(10)
    RSEN     =    PROPS(11)
                             
    YITA    =    PROPS(12)  !TIME STEP CONTROL, DEFAULT=0.1
    THETA    =    PROPS(13)
   
    NLGEOM    =    PROPS(14)

    ITCRL    =    PROPS(15)   !TIME INCREMENTAL INCREASE RATIO CONTROL ? 1=YES 0=NO
    RATIOT    =    PROPS(16)   !TIME RATIO. DEFAULT=1.02
    TLIMIT  =    PROPS(17)   !TIME LIMIT, AFTER THIS, TIME INCREMENT CONTROL IS APPLIED.

    IGEOSTEP=    PROPS(18)   !Swith to address the geostatic step problem. =0 default;;  =1 using the subroutine to deal with Geostatic step problem. 

C ----------------------------------------------------------------
C

C    LOCAL VARIBLES:
C
C....Stress State:
C    SMEAN  = Mean effective stress
C    SR2J2  = Square root of 2xJ2.
C

       
C
C------DEFINING PCSY
C

    IF ( IGEOSTEP.EQ.ONE. AND. TIME(2).LT.(ONE-TOLER) ) THEN
    !WHEN SPECIAL MEASURE IS TAKEN IN THE GEOSTEP WHERE TIME(2)<=1. A LARGE VALUE IS SET FOR PCSY.
        PCSY=1.E10
    ELSE
    !DEFAULT PROCESS TO DEFINING PCSY FOR A) ALL STEPS IF IGEOSTEP=0 or B) THE STEPS AFTER GEOSTEP IF IGEOSTEP=1
        !!CHECK:IS THIS THE FIRST TIME TO REACH PLASTIC STATE?
        IF (STATEV(1). GT. (TOLER+ZERO).
     &        AND. ABS(STATEV(1)-1.E10).GT.ONE ) THEN
            !! THE PCSY VALUE HAS BEEN ASSIGNED IN PREVIOUS STEPS, EXCLUDING THE SPECIAL MEASURE (PCSY=1.E10) IN THE GEOSTATIC STEP
            PCSY    =    STATEV(1)
        ELSE   
            !!THE FIRST TIME TO REACH PLASTIC STATE or RESUME THE INITIAL VALUE AFTER THE GEOSTATIC STEP
            IF ( PCSY0 . LT. 0.0) THEN
            !! *PCSY0 PROFILE* will be determined by a function, DEFINING THE VARIATION OF PCSY WITH DEPTH IN Y FOR 2D PROBLEM.
                CALL PCSY0PROFILE(PCSYINI,NDI,NSHR,COORDS)
                PCSY    =    PCSYINI
            ELSE
                PCSY    =    PCSY0
            ENDIF
        ENDIF
    ENDIF

C
C------ELASTIC MATRIX
C

C
C    LOCAL VARIABLES
C    Elastic Constants:
C    ENU    = Poisson's ratio
C    RKAP   = Slope of swelling curve.
C   
C    VOIDS  = voids ratio
C    DE     = Elastic stress strain matrix

    SMEAN = (STRESS(1)+STRESS(2)+STRESS(3)) / THREE   
!....   Check the validation of SMEAN,
    IF ( SMEAN.LE.0.) THEN
           WRITE(7,*)'*******WARNING! SMEAN IS NEGATIVE OR ZERO******'
           WRITE(7,*)'SMEAN,SR2J2,COORDS(1-2)'
           WRITE(7,*) SMEAN,SR2J2,COORDS(1),COORDS(2)   
    ENDIF       
C      EMOD = THREE*(ONE - TWO*ENU)*(ONE + VOIDS)*SMEAN/RKAP
      EMOD = THREE*(ONE - TWO*ENU)*(ONE + VOIDS)*PCSY/RKAP
!debug Dec 08
!         WRITE(7,*)'*******Time=*PCSY*STATEV(1) Coord****'
!          WRITE(7,*)i,Time(1), time(2),PCSY,STATEV(1),COORDS(1),COORDS(2)

    ! Debug by Qu 2008 Nov.
C    IF (PCSY.LE.50.0) THEN
C         WRITE(7,*)'*******WARNING! PCSY IS NEGATIVE OR ZERO******'
C           WRITE(7,*)'SMEAN,SR2J2,COORDS(1-2)'
C           WRITE(7,*) SMEAN,SR2J2,COORDS(1),COORDS(2)
C    ENDIF
            

      EMIN = THREE*(ONE - TWO*ENU)*(ONE + VOIDS)*SMEANMIN/RKAP
      IF(EMOD.LT.EMIN)EMOD=EMIN
      G  = EMOD/(2.*(1.0+ENU))
      AL = 1./EMOD
      DL = -ENU/EMOD
C      !EMOD    !E
       !ENU     !¦Í
      EBULK3=EMOD/(ONE-TWO*ENU)  !3K
      EG2=EMOD/(ONE+ENU)   !2G
      EG=EG2/TWO          !G
      EG3=THREE*EG  !3G
      ELAM=(EBULK3-EG2)/THREE  !¦Ë
C
C    Zero [D] ELASTIC matrix, [C] , [DE] AND STRAIN-RATE VECTOR, DVPR()
C
      DO 10 I = 1,NTENS
    DVPR(I) = 0.
       DO 10 J = 1,NTENS
        DE(I,J) = 0.

10    CONTINUE

       
      DE(1,1) = AL
      DE(2,1) = DL
      DE(3,1) = DL
      DE(1,2) = DL
      DE(2,2) = AL
      DE(3,2) = DL
      DE(1,3) = DL
      DE(2,3) = DL
      DE(3,3) = AL
      DE(4,4) = 1.0/G
    IF (NSHR.GT.1)THEN
    DE(5,5) = 1.0/G
    DE(6,6) = 1.0/G
    ENDIF
C    DE : The inverse of the elastic  matrix
C
C         1/E         -v/E     -v/E  
C         -v/E        1/E         -v/E
C   DE=   -v/E         -v/E     1/E       
C                                          1/G   
C                                            1/G 
C                                                1/G

      DO K1=1,NTENS
        DO K2=1,NTENS
           DDSDDE(K1,K2)=ZERO
        END DO
      END DO
      DO K1=1,NDI
        DO K2=1,NDI
          DDSDDE(K2,K1)=ELAM
        END DO
        DDSDDE(K1,K1)=EG2+ELAM
      END DO
      DO K1=NDI+1,NTENS
        DDSDDE(K1,K1)=EG
      END DO
C         ¦Ë +2G            ¦Ë         ¦Ë  
C         ¦Ë                 ¦Ë +2G     ¦Ë        
C   J=    ¦Ë                 ¦Ë         ¦Ë +2G      
C                                          G   
C                                            G 
C                                                G

 

C===============================
C------VISCOPLASTIC STRAIN RATE
C===============================

C----1. obtain the yield status of current stress point   
C     IYIELD:=0 ELASTIC
C            =1 OC  YIELDED
C            =3 CAP YIELDED
C            =2 CRITICAL STATE
C    2. DETERMINE THE VP RATE AT STEP 'N' USING DVPR

    CALL KVPRATE
    1  (STRESS,NDI,NSHR,NTENS,STATEV,NSTATV, PROPS,NPROPS,
    2   PCSY, IYIELD, PCDY, OVERSTRESS, DVPR)

    IF (IYIELD.LT.ONE) THEN
!        GO TO 8686 !END SUBROUTINE, ELASTIC DOMAIN  ! if elastic, update stress using elastic modulus
    ENDIF
C========================
C   CALCULATE DERIVATION : D VP-RATE / D STRESS
C========================

    CALL KDVPRDSTR
    1  (STRESS,NDI,NSHR,NTENS,STATEV,NSTATV, PROPS,NPROPS,
    2  DVPR, DVPRDSTR)

     CONTINUE

    CALL KJACOBIAN
    1  (NDI,NSHR,NTENS,DTIME,IYIELD,DE,DVPRDSTR,DDSDDE)

C===============================
C------CALCULATE STRESS
C===============================
      DO K1=1,NTENS
        DO K2=1,NTENS
          STRESS(K2)=STRESS(K2)
    1               + DDSDDE(K2,K1)* (DSTRAN(K1)-DVPR(K1) * DTIME)
        END DO
      END DO

    CONTINUE

C===============================
C------UPDATE STATE VARIABLES
C===============================

C    STATEV(1)   PCSY --- VOL-VP-STRAIN   
    DVOLVP = 0.D0
    DO N=1,NDI
        DVOLVP = DVPR(N)*DTIME + DVOLVP
    ENDDO
    PCSYF = EXP(   ( (1.+VOIDS)/(RLAM-RKAP) ) * DVOLVP + LOG(PCSY)  )
     
    STATEV(1) = PCSYF

    STATEV(2) = OVERSTRESS
    STATEV(3) = DVOLVP
   
           
    STATEV(4) = STATEV(4) + ONE/THREE/THREE*
    1        (1.+VOIDS)/(RLAM-RKAP)*DVOLVP*RLAM   
C      THIS IS FOR THE STRAIN 11 FROM THE CHANGE OF VOID RATIO   

C    STATEV(5) = STATEV(5) + DSTRAN(1)   
C      THIS IS FOR THE STRAIN 11 FROM THE CHANGE OF VOID RATIO

    STATEV(5) = POSTRESS

    STATEV(6) = IYIELD 
C     IYIELD:=0 ELASTIC
C            =1 OC  YIELDED
C            =3 CAP YIELDED
C            =2 CRITICAL STATE
8686  CONTINUE

C+===================
C  time incremental control, only for creep tests
C=====================
    IF (ITCRL.EQ.1)THEN
        IF (TIME(2).GT.TLIMIT)THEN
            PNEWDT = RATIOT
        ENDIF
    ENDIF
C
C==== CONTROL MAGNITUDE OF TIME STEP - NUMERICAL STABILITY
C        THE PAR. PNEWDT=NEW DT / OLD DT, i.e., NEW DT = PNEWDT* (OLD DT)

C    CALL KTIMESTEP(YITA,NDI,NSHR,NTENS,DSTRAN,DVPR,DTIME,PNEWDT)

!DEBUG BY QU

!    WRITE(7,*)'------------------------------'
!    WRITE(7,*)'------------------------------'
!    WRITE(7,*)'SMEAN,SR2J2,PCSY,PCDY,IYIELD,PCSYF'
!    WRITE(7,*)SMEAN,SR2J2,PCSY,PCDY,IYIELD,PCSYF
!    WRITE(7,*)'--------------------'
!    WRITE(7,*)'EXP((1.+VOIDS)/(RLAM-RKAP) ) * DVOLVP + LOG(PCSY)'
!    WRITE(7,*) VOIDS,RLAM-RKAP,DVOLVP,dtime,time
!    WRITE(7,*)'DDSDDE(1,1),DDSDDE(2,1),DDSDDE(4,4)'
!    WRITE(7,*)DDSDDE(1,1),DDSDDE(2,1),DDSDDE(4,4)
!    WRITE(7,*)'--------------------'
!    WRITE(7,*) 'STRAN(1), DSTRAN(1), DVPR(1),DVOLVP'
!    WRITE(7,*) STRAN(1), DSTRAN(1), DVPR(1), DVOLVP
!   
!    WRITE(7,*) 'PNEWDT, RATOT, TLIMIT,ITCRL'
!    WRITE(7,*) PNEWDT, RATOT, TLIMIT, ITCRL
!
!    WRITE(7,*) 'DVOLVP, DVPR(1) ' ,DVOLVP, DVPR(1)
C===========================================================================
C===    TRANFER STRESS AND STRAIN TO CONVENTIONAL DEFINITION IN SOIL MECHANICS.
C            (COMPRESSION IS POSITIVE)
C============================================================================
    DO N=1,NTENS
    STRESS(N) = -STRESS(N)
    STRAN(N)  = -STRAN(N)
    DSTRAN(N) = -DSTRAN(N)
    ENDDO

      RETURN
      END

!****************************************************************************
!
!  OBJECTIVE:    OBTAIN 1. IYIELD AND 2. VP RATE
!  BY:        GUANGFENG QU                    APRIL/2008
!
!****************************************************************************

    SUBROUTINE KVPRATE
    1  (STRESS,NDI,NSHR,NTENS,STATEV,NSTATV, PROPS,NPROPS,
    2   PCSY, IYIELD, PCDY, OVERSTRESS, DVPR)

C
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 PROPS(NPROPS), DVPR(NTENS)

C
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=10.D0,TOLER=1.0D-6,TOLMLINE=1.0D-3)
C SMEANMIN:    THE LOW BOUND OF MEAN STRESS TO CALCULATE 
C            THE STRESS-DEPENDENT ELASTIC MODULUS
C TOLMLINE:   THE ASSUMED THICKNESS RATIO (ONE-SIDE) AROUND THE M-LINE TO DEFINE THE ZONE OF IYIELD=2

    COMMON /KPROPS1/                        !ALL IN PROPS()
     1    ENU, RKAP,RLAM,    VOIDS, TALOC, TALNC, RCNC, ASPCT, PCSY0,        ! 1-9
     2    GAMA,RSEN,YITA, THETA                !10-13
    COMMON /KPROPS2/ NLGEOM,ITRATN,ITRMAX !14-17
    COMMON /TEMPPOSTRESS/POSTRESS
C ----------------------------------------------------------------

C    LOCAL VARIBLES:
C
C....Stress State:
C    SMEAN  = Mean effective stress
C    SR2J2  = Square root of 2xJ2.
C    OVERSTRESS    = measure of overstress.
C    IYIELD = yield status 0 = elastic
C                          1 = o/c yield
C                          2 = critical state
C                          3 = n/c yield

C....    CALCULATE SMEAN, AND SR2J2
    SIG11    = STRESS(1)
    SIG22    = STRESS(2)
    SIG33    = STRESS(3)
    SMEAN    =    (SIG11+SIG22+SIG33) / THREE
    !    FOR 3D AND PLANE STRAIN / AXISYMETRIC
    IF ( NSHR.EQ.1) THEN
        SR2J2=SQRT( (ONE/THREE*((SIG11-SIG22)**TWO + (SIG22-SIG33)**TWO
     +              + (SIG33-SIG11)**TWO) + TWO*STRESS(4)**TWO) )
      ELSE
      IF ( NSHR.EQ.3) THEN
        SR2J2=SQRT( ONE/THREE* ( (SIG11-SIG22)**TWO + (SIG22-SIG33)**TWO
     +              + (SIG33-SIG11)**TWO ) 
     +  +    TWO* (STRESS(4)**TWO + STRESS(5)**TWO +STRESS(6)**TWO)  )
      ENDIF
    ENDIF

C... DETERMINE the LOCATION of the STATIC CAP AND O/C LIMIT LINE

C     CENTRE OF STATIC CAP  (SMEAN-COORD)
C    HEIGHT OF STATIC CAP  (RS2J2-COORD)
      SCHECK = PCSY* (PCSY - ASPCT*RCNC)/(PCSY*(ASPCT*TALNC + ONE))
      R2J2SY = (PCSY-SCHECK)/ASPCT
C
C... Determine whether stress state is ELASTIC, CRITICAL STATE, CAP or O/C ,
C    Also set yield ratio for output

C    POSTRESS indicates more accurate the position of stress:: as following
C                 || |       ID=1.5                                                              
C                    |       |    /                                                             
C                    |       |   -ID=2                                                               
C                    | ID=1.0|  /                                                             
C                    |       | /                                                              
C                    |       |/                                                               
C                    |       /. . .     ID=3                                                           
C                    |      /|     .                                                            
C                    |     / |      .                                                          
C                    |    /  |       .                                                           
C                    |   /   |        .                                                        
C                    |  /    | ID=0.5  .                                                       
C                    | /ID=0 |          .-Apparent yield surface                                                                                                                                                    
C                    |/______|___________.____________________ 
C     Start with the elastic assumption
      IYIELD = 0
    POSTRESS = -1
      OVERSTRESS=0
      YRAT  =  0
    R2J2M =  SMEAN*TALNC + RCNC  !The sqrt(2J2) on the M-line corresponding to mean stress.

      IF(SMEAN.LT.SCHECK)THEN
      ! Only Two cases for this condition: Iyield = 0 or 1
        ! TWO options:  *O/C Limit Approach* or **O/C Parallel Approach**
      !                -------------------      ---------------------
         
      IF ( TALOC .GT.0) THEN
         !! O/C Limit Approach
C          R2J2F is the value of sqrt(2J2) of point on O/C LIMIT-line corresponding to
C          the current mean stress.   
           R2J2F = R2J2SY - (SCHECK-SMEAN)*TALOC
          
           IF(SR2J2.GE.R2J2F) THEN
            !!!When Stress is above OR ON the Static OC Limit Line
             IYIELD = 1
             POSTRESS = 1.
C        .... GET THE "PCDY"
            CALL KOCLIMIT(SMEAN,SR2J2,TALNC,TALOC,RCNC,ASPCT,PCDY)
              YRAT   = PCDY/PCSY
            OVERSTRESS  = PCDY-PCSY
               
         ELSE
         !!! When Stress is below  the Static OC Line, Iyield=0 Elastic state 
         POSTRESS = 0.
         ENDIF
        ELSE  !!TALOC .<=.0 i.e. O/C Parallel Approach
           IF(SR2J2.GE.R2J2M) THEN
            !!!When Stress is above OR ON the the Critical State Line
             IYIELD = 1
C         .... GET THE "PCDY"
            POSTRESS = 1.
            CALL KOCPARALLEL
     1          (SMEAN,SR2J2,TALNC,TALOC,RCNC,ASPCT,SCHECK,R2J2SY,PCDY) 
              YRAT   = PCDY/PCSY
            OVERSTRESS  = PCDY-PCSY   
         ELSE
            POSTRESS = 0
            !!! When Stress is below  the Static OC Line, Iyield=0 Elastic state 
         ENDIF
      ENDIF ! over for the two options when SMEAN.LT.SCHECK
    ELSE
   
    ! i.e. if SMEAN.>=.SCHECK
      IF (ABS( (SR2J2-R2J2M)/R2J2M).LT.TOLMLINE) THEN
          !! When Stress on or very close to the M line.
          IYIELD = 2
          POSTRESS = 2.
          PCDY   = SMEAN + R2J2M * ASPCT
            YRAT   = PCDY/PCSY
          OVERSTRESS  = PCDY-PCSY
      ELSE
        !!The stress is not on M-line.
          IF ( SR2J2.GE.R2J2M)THEN
          IYIELD = 1
          POSTRESS = 1.5
            !!! TWO options:  *O/C Limit Approach* or **O/C Parallel Approach**
          !!!                -------------------      ---------------------
            IF ( TALOC .GT.0) THEN
            !!!! O/C Limit Approach
             CALL KOCLIMIT(SMEAN,SR2J2,TALNC,TALOC,RCNC,ASPCT,PCDY) 
             YRAT   = PCDY/PCSY
           OVERSTRESS  = PCDY-PCSY
          ELSE
            !!!! O/C Parallel Approach
            CALL KOCPARALLEL
     1          (SMEAN,SR2J2,TALNC,TALOC,RCNC,ASPCT,SCHECK,R2J2SY,PCDY) 
             YRAT   = PCDY/PCSY
           OVERSTRESS  = PCDY-PCSY
          ENDIF
        ELSE !!!SR2J2.<.R2J2M
            !!! The stress lies Below the Critical State Line
            !!! Determine the position of the "SEMI-DYNAMIC" elliptical cap passing
            !!! through current state of stresses.
            CALL KCAPPC(SMEAN,SR2J2,TALNC,RCNC,ASPCT,PCDY)
            YRAT   = PCDY/PCSY
          OVERSTRESS  = PCDY-PCSY     
            IF(YRAT.GE.1.0D+00) THEN
            !!!! CAP PLASTIC ZONE
            IYIELD=3
            POSTRESS = 3.
          ELSE
            !!!!IN ELASTIC ZONE, INSIDE THE CAP.
            IYIELD = 0
            POSTRESS = 0.5
            OVERSTRESS = 0
            PCDY  =  0
          ENDIF !!!! OVER FOR THE CONDITION OF "YRAT.GE.1.0D+00"
        ENDIF !!! OVER FOR THE CASES CONCERNING SR2J2.V.S.R2J2M
      ENDIF !! OVER FOR THE CASE  WHEN The stress is not on M-line.
   
    ENDIF ! OVER FOR THE CONDITION OF "SMEAN.LT.SCHECK"

C.....OVER

    IF(IYIELD.EQ.0)THEN
C....    ELASTIC DOMAIN , STEP OUT FROM THIS SUBROUTINE
        GO TO 8686
    ENDIF
C========================================
C    OBTAIN THE VISCOPLASTIC STRAIN-RATE
C========================================

C.....MAGNITUDE OF THE VP-STRAIN-RATE
    VPHIF = GAMA*( (PCDY/PCSY)**(ONE/RSEN) - ONE )   

C.....THE PLASTIC POTENTIAL FOR VP-STRAIN-RATE
         
    IF(IYIELD.EQ.3)THEN                                                       
                                                                                  
                                                                                      
        SCHECKDY= PCDY* (PCDY - ASPCT*RCNC)/(PCDY*(ASPCT*TALNC + ONE)) !L
        A1T=TWO/THREE*(SMEAN-SCHECKDY)    ! DF/DSII PRICIPLE                                       
        A2T=TWO*ASPCT**TWO                    ! DF/DEIJ DEVIATERIC
       
C        NORMALIZED BY df/unit stress vecto                                
        AN = SQRT( 
     +    ( TWO*(SMEAN-SCHECKDY) )**TWO + (TWO*SR2J2*ASPCT**TWO)**TWO  )             

C        FOR NORMAL STRAIN-RATE 11,22,33                                       
           DO N=1,NDI                                      
              DVPR(N) = VPHIF* (A1T + (STRESS(N)-SMEAN)*A2T)/AN       
        ENDDO                                                             
C        FOR SHEAR STRAIN-RATE 12,13,23                                      
        DO N=1,NSHR                                                       
            DVPR(N+NDI) = VPHIF*( A2T*TWO*STRESS(N+NDI) ) / AN              
        ENDDO                                                             
      ENDIF                                                                            

C.... Evaluate the viscoplastic strain rate vector using
C     associated OR unassociated flow

      IF(IYIELD.EQ.2.OR.IYIELD.EQ.1)THEN
C   .....O/C ZONE
C           A) O/C LIMIT APPROACH: Dilation angle is from TALOC 
C        or B) O/C PARALLEL APPROACH: Dilation angle is from TALNC , i.e. M line  
C         DILATANCY is used to ajust the dilation angle based on DILATANCY*TALOC or DILATANCY*TALNC ( unassociated)
C         IF DILATANCY =1, USING AN ASSOCIATED FLOW
C       FACT is only used to further ajustment in debug. Not use so far.
C       TAL_ONC_C: the variable replaces TALOC in the origin code.    =TALOC for O/C LIMIT APPROACH;
C                                                                    =TALNC for Parallel APPROACH;
C                                                                    =0 For IYIELD.EQ.2.i.e.Critical state
        DILATANCY    = 0.2     !1.0D+00
        FACT    = 1.0D+00
          IF ( TALOC .GT. 0.0 ) THEN
          !O/C LIMIT APPROACH:
          TAL_ONC_C = TALOC
          ELSE
        !Parallel APPROACH:
          TAL_ONC_C = TALNC
        ENDIF
        !Critical state
        IF(IYIELD.EQ.2.) TAL_ONC_C = 0.0

               
        DEVPCOF=SQRT(1+(DILATANCY*TAL_ONC_C)**2)
C        FOR NORMAL STRAIN-RA1,22,33                                       
           DO N=1,NDI                                      
            DVPR(N) = VPHIF*FACT*((STRESS(N)-SMEAN)/(DEVPCOF*SR2J2) -
     &            (DILATANCY)*TAL_ONC_C/(DEVPCOF*THREE) )
        ENDDO
                                                                    
C        FOR SHEAR STRAIN-RAT12,13,23                                      
        DO N=1,NSHR                                                       
        DVPR(N+NDI)=VPHIF*FACT*(TWO*STRESS(N+NDI)/(DEVPCOF*SR2J2) )              
        ENDDO
                                                                     
      ENDIF

8686  CONTINUE   
      RETURN

    END
      SUBROUTINE KTIMESTEP(YITA,NDI,NSHR,NTENS,DSTRAN,DVPR,DTIME,PNEWDT)
C
      INCLUDE 'ABA_PARAM.INC'

      DIMENSION DSTRAN(NTENS), DVPR(NTENS),PDSTRAN(NDI),PDVPR(NDI)

C

      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=10.D0,TOLER=1.0D-6)
     
 
C.... GET THE MAGNITUDES OF STRAIN AND VP STRAIN-RATE
C     ( SEE THE DETAILS IN THE REPORT BY G.QU)

C        GET PRINCIPLE VALUES   
    CALL    SPRINC(DSTRAN,PDSTRAN,2,NDI,NSHR)
    CALL    SPRINC(DVPR,PDVPR,2,NDI,NSHR)
   

    DSTRANMAG = SQRT(PDSTRAN(1)**TWO+PDSTRAN(2)**TWO+PDSTRAN(3)**TWO)
    DSDVPRMAG = SQRT(PDVPR(1)**TWO+PDVPR(2)**TWO+PDVPR(3)**TWO) 

C.....THE MAXIMUM ALLOWABLE TIME STEP
    DTMAX = YITA*DSTRANMAG/DSDVPRMAG

C.... MAKE THE NEXT TIME STEP >= DTMAX

    IF ( DTMAX.LT.DTIME) PNEWDT = DTMAX/DTIME

C.....IN CASE OF THE FIRST TIME RUN, DTMAX=0.0, THEN SET PNEWDT=1
    IF (DTMAX.LT.TOLER) PNEWDT = ONE 

C. DEBUG BY QU

    WRITE(7,*) 'DTIME, DTMAX,PNEWDT',DTIME, DTMAX,PNEWDT
        RETURN
        END

      SUBROUTINE SPRINC(DSTRAN,PDSTRAN,NIT,NDI,NSHR)
C
      INCLUDE 'ABA_PARAM.INC'

      DIMENSION DSTRAN(NDI+NSHR), PDSTRAN(NDI)

C

      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0)

   
   
   
    DO I=1,3
      PDSTRAN(I)= DSTRAN(I)
    ENDDO

        RETURN
        END

!****************************************************************************
!
!  OBJECTIVE:    OBTAIN deviation of vp strain-rate over deviation of stress
!  BY:        GUANGFENG QU                    APRIL/2008
!
!****************************************************************************

    SUBROUTINE KDVPRDSTR
    1  (STRESS,NDI,NSHR,NTENS,STATEV,NSTATV, PROPS,NPROPS,
    2  DVPR0, DVPRDSTR)

      INCLUDE 'ABA_PARAM.INC'
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 PROPS(NPROPS), DVPR0(NTENS), STRESST(NTENS),DVPRT(NTENS),
     2 DVPRDSTR(NTENS,NTENS)

      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=10.D0,TOLER=1.0D-6)

      PARAMETER(DELTSTR=1.D-5)
   
    COMMON /KPROPS1/                        !ALL IN PROPS()
     1    ENU, RKAP,RLAM,    VOIDS, TALOC, TALNC, RCNC, ASPCT, PCSY0,        ! 1-9
     2    GAMA,RSEN,YITA, THETA                !10-13
    COMMON /KPROPS2/ NLGEOM,ITRATN,ITRMAX !14-17

C ----------------------------------------------------------------------------------
C ----------------------------------------------------------------------------------

    DO I = 1, NTENS ! FOR STRESS 11,22,33, 12,13,23

        DO K = 1,NTENS
            STRESST(K) = STRESS(K)
        ENDDO
        STRESST(I) = STRESS(I) + DELTSTR
        CALL KVPRATE
    1    (STRESST,NDI,NSHR,NTENS,STATEV,NSTATV, PROPS,NPROPS,
    2    PCSY0, IYIELD, PCDY, OVERSTRESS, DVPRT)

        DO J=1, NTENS     !FOR STRAIN-RATE 11,22,33,12,13,23
            DVPRDSTR(J,I) = ( DVPRT(J) - DVPR0(J) ) / DELTSTR
        ENDDO
    ENDDO
    continue

C========================================
C     CORRECTION FOR UNSYMMETRIC MATRIX TO A SYMMETRIC ONE // TAKE THE LOWER TRIANGLE
C========================================
    DO I=1, NTENS     
        DO J=I+1,NTENS
    ! 1. TAKE THE AVERAGE OF THE TWO DIAGNAL COMPONENTS
    !! NO! THIS MAKES THE NUMERICAL ERROR MORE SIGNIFICANT
            !DVPRDSTR(I,J) = ( DVPRDSTR(I,J) + DVPRDSTR(J,I) ) / TWO
            !DVPRDSTR(J,I) = DVPRDSTR(I,J)
   
    ! 2. TAKE THE LOWER TRIANGLE
            DVPRDSTR(I,J) = DVPRDSTR(J,I)
        ENDDO
    ENDDO
       
    RETURN
    END

!****************************************************************************
!
!  OBJECTIVE:    OBTAIN deviation of vp strain-rate over deviation of stress
!  BY:        GUANGFENG QU                    APRIL/2008
!
!****************************************************************************

    SUBROUTINE KJACOBIAN
    1  (NDI,NSHR,NTENS,DTIME,IYIELD,DE,DVPRDSTR,DDSDDE)
    
C     STRESS,NDI,NSHR,NTENS,STATEV,NSTATV, PROPS,NPROPS,
C    2  DVPR0, DVPRDSTR)

      INCLUDE 'ABA_PARAM.INC'

      DIMENSION
     1 DDSDDE(NTENS,NTENS),
     2 DE(NTENS,NTENS),DVPRDSTR(NTENS,NTENS)

      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=10.D0,TOLER=1.0D-6)

    COMMON /KPROPS1/                        !ALL IN PROPS()
     1    ENU, RKAP,RLAM,    VOIDS, TALOC, TALNC, RCNC, ASPCT, PCSY0,        ! 1-9
     2    GAMA,RSEN,YITA, THETA                !10-13
    COMMON /KPROPS2/ NLGEOM,ITRATN,ITRMAX !14-17   
 
C ----------------------------------------------------------------------------------\
 
C===================   
C.... ELASTIC DOMAIN OR THETA=0 (EXPLICIT ALGORITHM)
C===================

C        1        SET DDSDDE AS ELASTIC MATRIX, SEE THE INITIAL VALUE IN THE MAIN SUBROUTINE
C        2        THEN STEP OUT FROM THIS SUBROUTINE
    IF (IYIELD.EQ.0.OR.THETA.LT.TOLER) THEN

    GOTO 9999 !STEP OUT FROM THIS SUBROUTINE
    ENDIF

C==================
C.... EVP DOMAIN
C==================
   
   
    IF (IYIELD.GT.TOLER) THEN
   

      DO 20 I=1,NTENS
       DO 20 J=1,NTENS
        DDSDDE(I,J) = DE(I,J) + THETA*DTIME*DVPRDSTR(I,J)
20    CONTINUE

    ! GET THE INVERSE OF DDSDDE
      DO 200 N=1,NTENS
      D=DDSDDE(N,N)
      DO 100 J=1,NTENS
  100 DDSDDE(N,J)=-DDSDDE(N,J)/D
      DO 150 I=1,NTENS
      IF(N-I) 110,150,110
  110 DO 140 J=1,NTENS
      IF(N-J) 120,140,120
  120 DDSDDE(I,J)=DDSDDE(I,J)+DDSDDE(I,N)*DDSDDE(N,J)
  140 CONTINUE
  150 DDSDDE(I,N)=DDSDDE(I,N)/D
      DDSDDE(N,N)=ONE/D
  200 continue   
    ENDIF !FOR EVP DOMAIN

9999  CONTINUE
    RETURN
    END

      SUBROUTINE KOCLIMIT(SMEAN,SR2J2,TALNC,TALOC,RCNC,ASPCT,PCDY)
C
      INCLUDE 'ABA_PARAM.INC'
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=10.D0,TOLER=1.0D-6)
     
C OBJECTIVE: TO OBTAIN "PCDY" ACCORESPONDING TO
C            THE CURRENT STRESS STATE IN O/C DOMAIN (IYEILD=0)
C    BY G. QU APRIL/2008

C    STEP1   
C    TWO EQUATIONS:
C        Y-Y0 = KOC * (X-X0)      1   THE LIMIT LINE WITH GIVEN STRESS POINT
C        Y    = KNC*X + CNC       2   THE CRITICAL STATE LINE OR M LINE
C    STEP2    SOLVING THESE PROVIDED THE INTERCECTION POINT ON THE M LINE
C            X=(Y0-CNC-KNC*X0)/(KNC-KOC)
C    STEP3    BASED ON ASPECT RATIO, PCDY CAN BE OBTAINED

    SCHECK = ( SR2J2-RCNC-TALOC*SMEAN ) / (TALNC-TALOC) !X
    SCHECK2J2 = SCHECK*TALNC + RCNC

    PCDY = SCHECK2J2*ASPCT + SCHECK
        RETURN
        END

      SUBROUTINE KCAPPC(SMEAN,SR2J2,TALNC,RCNC,ASPCT,PCDY)
C
      INCLUDE 'ABA_PARAM.INC'

C
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=10.D0,TOLER=1.0D-6)
     
         AA = 0.0d0
         BB = 0.0d0
         CC = 0.0d0
         DD = 0.0D0
         AA = (1.0D0-ASPCT*TALNC)/(ASPCT*TALNC+1.0D0)
c
         BB = (-2.0D0*SMEAN - 2.0D0*ASPCT*RCNC)/(ASPCT*TALNC+1.0D0)
         CC = SMEAN*SMEAN + 2.0D0*SMEAN*ASPCT*RCNC/(ASPCT*TALNC+1.0D0)
     $       + ASPCT*ASPCT*SR2J2*SR2J2
c
       IF(AA.GT.0..OR.AA.LT.0.)THEN
          DD = BB*BB-4.0D0*AA*CC
          IF(DD.LT.0.0D0)THEN
           WRITE(7,*)'*******WARNING OVERSTRESS IS LESS THAN ZERO******'
           WRITE(7,*)'SMEAN,SR2J2,DD,  CHECK SUBROUTINE KCAPPC'
           WRITE(7,*) SMEAN,SR2J2,DD
         ENDIF
         PCDY = (-BB - SQRT(ABS(DD)))/(2.0D0*AA)
       ELSE
        PCDY = -CC/BB
       ENDIF
        RETURN
        END

!****************************************************************************
!
!  OBJECTIVE:    Determine PCSY0 according to the PCSYO profile defined by 5 points
!        MAXIUM NUMBER OF POINTS =5 controled by   "HN(5), DVPR(5)"
!  BY:        GUANGFENG QU                    NOV.12/2008
!
!****************************************************************************
      SUBROUTINE PCSY0PROFILE
     1  (PCSY0I,NDI,NSHR,COORDS)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION COORDS(3),
     1 HN(5), PN(5)

C ----------------------------------------------------------------
C -----INPUT DATA: START------------------------------------------

!    NPOINTS    = The total Number of reference points (Max=5)
!    HN(n)    =    Vertical Coord of the nth reference point
!    PN(n)    =    Yield Stress    of the nth reference point
!                                **********************
! NOTE:  the HN() "MUST" be in sort ****ASECNDING*******
!                                **********************

    NPOINTS    =5
    HN(1)    =    140
    HN(2)    =    173
    HN(3)    =    177
    HN(4)    =    180
    HN(5)    =    185

    PN(1)    =    375.4
    PN(2)    =    98.8
    PN(3)    =    98.8 
    PN(4)    =    110.0
    PN(5)    =    110.0

!        H        Pc Profile 3    s(d)y(Profile 1)
!1        140     438                 375.4 
!2        168     115                 98.8 
!3        172     115                 98.8 
!4        175     184                 158.0
!5              180      184                158.0
C -----INPUT DATA: END------------------------------------------
C ----------------------------------------------------------------

C....    OBTAIN VERTICAL COORD for the current Gauss point

    IF ( NSHR.EQ.1) THEN            !    FOR PLANE STRAIN / AXISYMETRIC
        H    =    COORDS(2)
      ELSEIF ( NSHR.EQ.3) THEN        !    FOR 3D
          H    =    COORDS(3)
      ELSE
      WRITE(7,*)'***WARNING: Error in PSCY0FUN subroutine -H******'
    ENDIF

C... DETERMINE the yield stress according to the VERTICAL COORD (H) for the current Gauss point
    PC=0.0
    DO I=2, NPOINTS
      IF ( H.GE.HN(I-1).AND.H.LE.HN(I))Then
        PC    =    PN(I-1)  + 
     &        (H-HN(I-1))*(  (PN(I)-PN(I-1)) / (HN(I)-HN(I-1))  )
      GOTO 1000
      ENDIF
    ENDDO

1000    CONTINUE

    IF (PC.EQ.0.0)THEN
      WRITE(7,*)'***WARNING: Error in PSCY0FUN subroutine -PC******'
    ENDIF
   
    PCSY0I    = PC
   
8686  CONTINUE   
      RETURN
    END

      SUBROUTINE KOCPARALLEL
     1    (SMEAN,SR2J2,TALNC,TALOC,RCNC,ASPCT,SCHECK,R2J2SY,PCDY)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      INCLUDE 'ABA_PARAM.INC'
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
     1          SMEANMIN=10.D0,TOLER=1.0D-6)
     
C OBJECTIVE: TO OBTAIN "PCDY" CORESPONDING TO THE CURRENT STRESS STATE IN O/C DOMAIN (IYEILD=1)
C        ***************************           
!        ****Parallel Approach******
!        ***************************
!        The dynamic yield surface in O/C zone includes an inclinded line with M slope;+ a horizontal line
!              a)an inclinded line with M slope:
!             When stress is on the left of the centr of static surface, OVERSTRESS is defined as the vertical disp from the current stress to M line.
!              b)a horizontal line
!             When stress is on the right of the centr of static surface, OVERSTRESS is the vertical disp from the current stress to the top of static surface.
C    BY G. QU Nov./2008
C    STEP1
C      Location of the current stress:  left or right from the static yield surface centre.
C    STEP2    Obtain overtress           
C    STEP3    BASED ON ASPECT RATIO, PCDY CAN BE OBTAINED

!        SCHECK: Mean stress of the static yield surface centre
!        R2J2SY: 2J2 of the static yield surface top

    IF ( SCHECK.LT.SMEAN)THEN
    ! ON THE RIGHT SIDE
    OVERSTRESS = SR2J2 - R2J2SY
    ELSE
    ! ON THE LEFT SIDE
    OVERSTRESS = SR2J2 - (SMEAN*TALNC+RCNC)
    ENDIF

    PCDY = SCHECK + (R2J2SY + OVERSTRESS)*ASPCT

!    DEBUG
    IF (OVERSTRESS .LE.0)THEN
       WRITE(7,*)'***WARNING! OVERSTRESS IS NEGATIVE*SUB-KOCPARALLEL***'
           WRITE(7,*)'SMEAN,SR2J2,OVERSTRESS'   
    ENDIF
      RETURN
      END    

Dr. G. Qu;

I am interested in running the benchmark examplesbut I could not find the input files, can you provide me with the input files of some of examples?!

Thanks

Amir

Hi

I have a problem with UMAT subroutine usinage.

I wrote a UMAT subroutine for my material and I want to implement it in Abaqus but I don't know how can I do this. How can I introduce this code to Abaqus and run the simulation. I'll be glad if you help me to do this.

Sincerely,

Mehran Ashrafi

Subscribe to Comments for "Sharing an abaqus UMAT using an  elastic viscoplastic model for soil"

Recent comments

More comments

Syndicate

Subscribe to Syndicate