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.
-
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. -
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
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/0B8z7TOlNv4UGMTFjYjY0NDItOGQ2My00NThmLT…
https://drive.google.com/file/d/0B8z7TOlNv4UGZGY0MTNlZTItZWI2OS00ZmY4LW…
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 [at] hotmail.com (guangfengqu[at]hotmail[dot]com)
Links do not work
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
In reply to Links do not work by cjfUNSW
Here is the
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
Input files
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
umat
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