Help regarding Johnson Cook VUHARD
Help regarding Johnson Cook VUHARD
(OP)
Dear all,
I am trying to write VUHARD for Johnson Cook plasticity model.
Which is given by:
sigma = [A+ B epsilon^n] [1+C ln epsilon_dot_star] [ 1-T*^m]
This runs and gives a result very different to what the inbuilt JC model (within Abaqus) shows. I have attached the figure showing the difference. Can anyone please point out my mistake.
The mises stress vs time plot is the following:
https://ibb.co/cAmzum
I have written the following VUHARD
C
C User subroutine vuhard
subroutine vuhard (
C Read only -
* nblock,
* nElement, nIntPt, nLayer, nSecPt,
* lAnneal, stepTime, totalTime, dt, cmname,
* nstatev, nfieldv, nprops,
* props, tempOld, tempNew, fieldOld, fieldNew,
* stateOld,
* eqps, eqpsRate,
C Write only -
* yield, dyieldDtemp, dyieldDeqps,
* stateNew )
c
include 'vaba_param.inc'
c
dimension nElement(nblock),
* props(nprops),
* tempOld(nblock),
* fieldOld(nblock,nfieldv),
* stateOld(nblock,nstatev),
* tempNew(nblock),
* fieldNew(nblock,nfieldv),
* eqps(nblock),
* eqpsRate(nblock),
* yield(nblock),
* dyieldDtemp(nblock),
* dyieldDeqps(nblock,2),
* stateNew(nblock,nstatev)
PARAMETER ( zero = 0.0d0, one = 1.0d0 , two = 2.0d0, three = 3.0d0,
1 third = one/three, half = .5d0, twoThirds = two/three,
2 threeHalfs = 1.5d0, eighteen = 18.0d0,
3 four = 4.0d0, fourThirds = four/three)
c
character*80 cmname
A = props(1)
B = props(2)
n = props(3)
C = props(4)
xm = props(5)
T_melt = props(6)
T_trans = props(7)
epsdot0 = props(8)
c
c epsdot0 is 1
c
DO k = 1,nblock
T = tempNew(k)
IF (T .lt. T_trans) THEN
T_star = zero
ELSEIF ( T_trans .ge. T .le. T_melt ) THEN
T_star = (T - T_trans)/(T_melt - T_trans)
ELSE
T_star = one
END IF
strainrate = eqpsRate(k)/epsdot0
strain = eqps(k)
IF (strainrate .gt. zero) THEN
edot_log = LOG( strainrate )
ELSE
edot_log = zero
ENDIF
yield(k) = (A + B*(strain**n))*(one - T_star**xm)
1 * (one + C*edot_log)
IF( strain .le. zero) THEN
dyieldDeqps(k,1) = zero
ELSE
dyieldDeqps(k,1) = (B*n*(strain**(n-one)))*
1 (one - T_star**xm) * (one + C*edot_log)
ENDIF
IF(strainrate .ge. epsdot0) then
dyieldDeqps(k,2) = (A + B*(strain**n))*
1 (one - T_star**xm) * (C/strainrate)
ELSE
dyieldDeqps(k,2) = zero
ENDIF
IF (T .gt. T_trans) THEN
dyieldDtemp(k) = (-xm)*(one/(T - T_trans))*
1 (T_star**(xm))*(A + B*(strain**n))
2 *(one + C*edot_log)
ELSE
dyieldDtemp(k) = zero
ENDIF
END DO
C
RETURN
END
Thanks!!
I am trying to write VUHARD for Johnson Cook plasticity model.
Which is given by:
sigma = [A+ B epsilon^n] [1+C ln epsilon_dot_star] [ 1-T*^m]
This runs and gives a result very different to what the inbuilt JC model (within Abaqus) shows. I have attached the figure showing the difference. Can anyone please point out my mistake.
The mises stress vs time plot is the following:
https://ibb.co/cAmzum
I have written the following VUHARD
C
C User subroutine vuhard
subroutine vuhard (
C Read only -
* nblock,
* nElement, nIntPt, nLayer, nSecPt,
* lAnneal, stepTime, totalTime, dt, cmname,
* nstatev, nfieldv, nprops,
* props, tempOld, tempNew, fieldOld, fieldNew,
* stateOld,
* eqps, eqpsRate,
C Write only -
* yield, dyieldDtemp, dyieldDeqps,
* stateNew )
c
include 'vaba_param.inc'
c
dimension nElement(nblock),
* props(nprops),
* tempOld(nblock),
* fieldOld(nblock,nfieldv),
* stateOld(nblock,nstatev),
* tempNew(nblock),
* fieldNew(nblock,nfieldv),
* eqps(nblock),
* eqpsRate(nblock),
* yield(nblock),
* dyieldDtemp(nblock),
* dyieldDeqps(nblock,2),
* stateNew(nblock,nstatev)
PARAMETER ( zero = 0.0d0, one = 1.0d0 , two = 2.0d0, three = 3.0d0,
1 third = one/three, half = .5d0, twoThirds = two/three,
2 threeHalfs = 1.5d0, eighteen = 18.0d0,
3 four = 4.0d0, fourThirds = four/three)
c
character*80 cmname
A = props(1)
B = props(2)
n = props(3)
C = props(4)
xm = props(5)
T_melt = props(6)
T_trans = props(7)
epsdot0 = props(8)
c
c epsdot0 is 1
c
DO k = 1,nblock
T = tempNew(k)
IF (T .lt. T_trans) THEN
T_star = zero
ELSEIF ( T_trans .ge. T .le. T_melt ) THEN
T_star = (T - T_trans)/(T_melt - T_trans)
ELSE
T_star = one
END IF
strainrate = eqpsRate(k)/epsdot0
strain = eqps(k)
IF (strainrate .gt. zero) THEN
edot_log = LOG( strainrate )
ELSE
edot_log = zero
ENDIF
yield(k) = (A + B*(strain**n))*(one - T_star**xm)
1 * (one + C*edot_log)
IF( strain .le. zero) THEN
dyieldDeqps(k,1) = zero
ELSE
dyieldDeqps(k,1) = (B*n*(strain**(n-one)))*
1 (one - T_star**xm) * (one + C*edot_log)
ENDIF
IF(strainrate .ge. epsdot0) then
dyieldDeqps(k,2) = (A + B*(strain**n))*
1 (one - T_star**xm) * (C/strainrate)
ELSE
dyieldDeqps(k,2) = zero
ENDIF
IF (T .gt. T_trans) THEN
dyieldDtemp(k) = (-xm)*(one/(T - T_trans))*
1 (T_star**(xm))*(A + B*(strain**n))
2 *(one + C*edot_log)
ELSE
dyieldDtemp(k) = zero
ENDIF
END DO
C
RETURN
END
Thanks!!




