Hello:
I have written a UMAT subroutine similar to the one shown in abaqus 6.11 subroutine manual 1.1.40 with a slight differences indicative of the material's rheological model.
my expanded isotropic model looks like:
sigma_x - sigmav_x/3 + vbar*sigma_dot_x + kappa*sigma_dot_x =
lambda*eplv + 2*mu*epl_x + lamdabar*eplv_dot_x +2*mubar*epl_x
where:
sigma = stress
epl = strain
_dot_ are my rate change with time
sigmav = sigmaxx + sigma_yy + sigma_zz ( hydrostatic stress)
eplv = elp_xx + elp_yy + epl_zz ( hydrostatic strain)
lambda, vbar, kappa, mu, lambabar, mubar are all material constants that changes has a function of temp.
I have accounted for time evolution based on midpoint method similar to the example. I have also consider time evolution in my material properties since they are temp dependent.
All material props are defined in the UMAT accordingly.
Some assumption i made were
DELTAsigmav ~ 3*DELTAsigma_x
for material properties:
lambda = lambda(T[i]) + DELTAlambda(T[i]-T[i-1])
...
To test out my umat i have a single element inp file with the following steps similar to experimental approach:
step 1: smp material at 348K is loaded by -5e5MPa
step 2: keep deformation constant and cool material down to 308K
Step 3: remove constraint and allow material to deform at 308K
Step 4: heat material back up to 348K
Based on the experimental results the expand in step 1 elastically due to load, strain remains the same in step 2 due to constrain, contract to a little in step 3 due to low temp as it freezes at the elastic stretch and contract back downward at step 4 since no load applied.
My results are not indicative of this results (see attachment). While I notice similar trends from steps1 through 3, step 4 is completely off.
I am not sure if my material time evolution assumption is wrong? My step 4 is not indicative of the test result as mentioned.
The total strain in step 1 looks to be off by a factor of 2
I will appreciate any help i can get.
Thanks
c 1 2 3 4 5
c23456789012345678901234567890123456789012345678901234567890
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),
2 DDSDDT(NTENS),DRPLDE(NTENS),
3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
DIMENSION DSTRES(6),D(3,3)
DIMENSION SIGold(NTENS),ETHERM(NTENS),DTHERM(NTENS),
1 EELAS(NTENS),ECRP(NTENS),ESLIP(NTENS),DSTRESS(NTENS),
2 DESLIP(NTENS),DECRP(NTENS), SEGold(NTENS),SIGold2(NTENS)
REAL KE, GE, KMe, GMe, GV, KV, alpha
REAL lambda, lambdabar, mu, mubar, rho, kappa
REAL KE0, GE0, KMe0, GMe0, GV0, KV0
REAL dlambda, dlambdabar, dmu, dmubar, drho, dkappa
REAL SIGsq, TrSIG, SIGterm, omegabar, zetabar,thetabar
REAL Eg, ae, GVg, agv, Tg, Psng, CRPLg, acr, CSLIPg,
1 acsl, rtimeg, art
C REAL E,EE,EMe, Psn, CSLIP, CRPL, rtime, TEMPo, TEMPcr
C REAL E0,EE0,EMe0, Psn0, CSLIP0, CRPL0, rtime0
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Props(1) = modulus_glass; Props(2) = mod_gls _const; Props(3) = viscosity_glass
C Props(4) = vis_gls_const; Props(5) = poisson_glass; Props(6) = creep_L_glass
C Props(7) = crpL_gls_const; Props(8)= slip_const_glass; Props(9) = slpC_gls_const
C Props(10) = relax_time_glass; Props(11) = relax_gls_cnst; Props(12)= Initial Temp
C Props(13) = glass_temp; Props(14) = CTE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
Eg = PROPS(1)
ae = PROPS(2)
GVg = PROPS(3)
agv = PROPS(4)
Psng = PROPS(5)
CRPLg = PROPS(6)
acr = PROPS(7)
CSLIPg = PROPS(8)
acsl = PROPS(9)
rtimeg = PROPS(10)
art = PROPS(11)
TEMPo = PROPS(12)
Tg = PROPS(13)
alpha = PROPS(14)
CCCCCCCCCCCCCCCCCCCCCCCCC
IF (TEMP .EQ. 0.0) TEMP = TEMPo
C PROPERTIES BEFORE INCREMENT
IF ( (TEMP - Tg) .GT. 15.0) THEN
TEMPcr = Tg + 15
ELSEIF ( (TEMP - Tg ) .LT. -15.) THEN
TEMPcr = Tg - 15
ELSE
TEMPcr = TEMP
ENDIF
rtime0 = rtimeg * exp(art*(Tg/(TEMPcr) - 1))
C Compute both mouludus Ee and Eme
E0 = Eg * exp(ae*(Tg/(TEMPcr) - 1))
C Need to consider poisson ratio
apsn = 5
Psn0 = Psng * exp(-apsn*(Tg/(TEMPcr) - 1))
IF ( Psn0 .GT. 0.48)THEN
Psn0 = 0.48
ELSE IF ( Psn0 .LT. 0.33)THEN
Psn0 = 0.33
END IF
GV0 = GVg * exp(agv*(Tg/(TEMPcr) - 1))
EE0 = GV0/rtime0
EMe0 = E0 - EE0
KMe0 = (EMe0)/(3*(1-2*Psn0))
GMe0 = (EMe0)/(2*(1+Psn0))
KE0 = (EE0)/(3*(1-2*Psn0))
GE0 = (EE0)/(2*(1+Psn0))
C PROPERTIES AFTER INCREMENT
IF ( ((TEMP+DTEMP) - Tg) .GT. 15.0) THEN
TEMPcr = Tg + 15
ELSEIF ( ((TEMP+DTEMP) - Tg ) .LT. -15.) THEN
TEMPcr = Tg - 15
ELSE
TEMPcr = TEMP+DTEMP
ENDIF
rtime = rtimeg * exp(art*(Tg/(TEMPcr) - 1))
C Compute both mouludus Ee and Eme
E = Eg * exp(ae*(Tg/(TEMPcr) - 1))
C Need to consider poisson ratio
apsn = 5
Psn = Psng * exp(-apsn*(Tg/(TEMPcr) - 1))
IF ( Psn .GT. 0.48)THEN
Psn = 0.48
ELSE IF ( Psn .LT. 0.33)THEN
Psn = 0.33
END IF
GV = GVg * exp(agv*(Tg/(TEMPcr) - 1))
EE = GV/rtime
EMe = E - EE
KMe = (EMe)/(3*(1-2*Psn))
GMe = (EMe)/(2*(1+Psn))
KE = (EE)/(3*(1-2*Psn))
GE = (EE)/(2*(1+Psn))
IF (DTEMP .LT. 0.001) THEN
dlambda = 0.
dlambdabar = 0.
dmu = 0.
dmubar = 0.
drho = 0.
dkappa = 0.
ELSE
dlambda = 0.5*((-2.0*(GE-GEO)) / 3.0)
dlambdabar = 0.5*(2.0* (GV-GVO)*( - (GE-GEO)
1 /(GMe-GMeO) + (KE-KEO)/(KMe-KMeO) / 3.0))
dmu = 0.5*(GE-GEO)
dmubar = 0.5*(GV*( 1.0 + (GE-GEO)/(GMe-GMeO)))
drho =0.5* ((GV-GVO)/(GMe-GMeO))
dkappa = 0.5*((2*(GV-GVO)*(1/(3*(KMe-KMeO)
1 + 1/(2*(GMe-GMeO)))))/3.0)
END IF
lambda = ((-2.0*GE) / 3.0) + dlambda
lambdabar = (2.0*GV*( - GE/GMe + KE/KMe )/3.0)
1 + dlambdabar
mu = GE + dmu
mubar = (GV*( 1.0 + GE/GMe ))+dmubar
rho = (GV/GMe)+drho
kappa = 2*GV*((1/(3*KMe) + 1/(2*GMe)))/3.0 +dkappa
C CALCULATE THERMAL STRAIN
DO K1=1,NDI
ETHERM(K1)=alpha*(TEMP-TEMPo)
DTHERM(K1)=alpha*(DTEMP)
END DO
DO K1=NDI+1,NTENS
ETHERM(K1)=0.0
DTHERM(K1)=0.0
END DO
C EVALUATE NEW STRESS TENSOR
C
EV = 0.
DEV = 0.
SV = 0.
DO K1=1,NDI
EV = EV + (STRAN(K1)-ETHERM(K1))
DEV = DEV + (DSTRAN(K1)-DTHERM(K1))
SV = SV+STRESS(K1)
C write (*,*) "EV,DEV,SV = ", EV,DEV,SV
END DO
C
TERM1 = .5*DTIME + rho + 3.*kappa - .5*DTIME
TERM1I = 1./TERM1
TERM2 = (.5*DTIME*lambda+lambdabar)*TERM1I*DEV
TERM3 = (DTIME*mu+2.*mubar)*TERM1I
C
DO K1=1,NDI
DSTRES(K1) = TERM2+TERM3*(DSTRAN(K1)-DTHERM(K1))
1 +DTIME*TERM1I*(lambda*EV
2 +2.*mu*(STRAN(K1)-ETHERM(K1))
3 - STRESS(K1)+ SV/3)
STRESS(K1) = STRESS(K1) + DSTRES(K1)
END DO
C
TERM2 = (.5*DTIME*mu + mubar)*TERM1I
I1 = NDI
DO K1=1,NSHR
I1 = I1+1
DSTRES(I1) = TERM2*DSTRAN(I1)+
1 DTIME*TERM1I*(mu*STRAN(I1)-STRESS(I1))
STRESS(I1) = STRESS(I1)+DSTRES(I1)
END DO
C
C CREATE NEW JACOBIAN
C
TERM2 = (DTIME*(.5*lambda+mu)+lambdabar+
1 2.*mubar)*TERM1I
TERM3 = (.5*DTIME*lambda+lambdabar)*TERM1I
DO K1=1,NTENS
DO K2=1,NTENS
DDSDDE(K2,K1) = 0.
END DO
END DO
C
DO K1=1,NDI
DDSDDE(K1,K1) = TERM2
END DO
C
DO K1=2,NDI
N2 = K1-1
DO K2=1,N2
DDSDDE(K2,K1) = TERM3
DDSDDE(K1,K2) = TERM3
END DO
END DO
TERM2 = (.5*DTIME*mu+mubar)*TERM1I
I1 = NDI
DO K1=1,NSHR
I1 = I1+1
DDSDDE(I1,I1) = TERM2
END DO
C
C TOTAL CHANGE IN SPECIFIC ENERGY
C
TDE = 0.
DO K1=1,NTENS
TDE = TDE + (STRESS(K1)-.5*DSTRES(K1))
1 *(DSTRAN(K1)-DTHERM(K1))
END DO
C
C CHANGE IN SPECIFIC ELASTIC STRAIN ENERGY
C
TERM1 = lambda + 2.*mu
DO K1=1,NDI
D(K1,K1) = TERM1
END DO
DO K1=2,NDI
N2 = K1-1
DO K2=1,N2
D(K1,K2) = lambda
D(K2,K1) = lambda
END DO
END DO
DEE = 0.
DO K1=1,NDI
TERM1 = 0.
TERM2 = 0.
DO K2=1,NDI
TERM1 = TERM1 + D(K1,K2)
1 *(STRAN(K2)-ETHERM(K2))
TERM2 = TERM2 + D(K1,K2)
1 *(DSTRAN(K2)-DTHERM(K2))
END DO
DEE = DEE + (TERM1+.5*TERM2)
1 * (DSTRAN(K1)-DTHERM(K1))
END DO
I1 = NDI
DO K1=1,NSHR
I1 = I1+1
DEE = DEE + mu*(STRAN(I1)+.5*DSTRAN(I1))*DSTRAN(I1)
END DO
SSE = SSE + DEE
SCD = SCD + TDE - DEE
DO K1=1, NTENS
ETHERM(K1)=ETHERM(K1)+DTHERM(K1)
EELAS(K1)=STRAN(K1)+DSTRAN(K1)-ETHERM(K1)
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=ETHERM(K1)
END DO
RETURN
END
INPUT FILE
*Heading
** Job name: smp_1_ele Model name: Model-1
** Generated by: Abaqus/CAE 6.11-2
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
**INCLUDE,INPUT=tobvisomod.f
** PARTS
**
*Part, name=Part-1
*Node
1, -0.0199999996, -0.0199999996, 0.00399999991
2, -0.0199999996, 0.0199999996, 0.00399999991
3, -0.0199999996, -0.0199999996, 0.
4, -0.0199999996, 0.0199999996, 0.
5, 0.0199999996, -0.0199999996, 0.00399999991
6, 0.0199999996, 0.0199999996, 0.00399999991
7, 0.0199999996, -0.0199999996, 0.
8, 0.0199999996, 0.0199999996, 0.
*Element, type=C3D8T
1, 5, 6, 8, 7, 1, 2, 4, 3
*Nset, nset=_PickedSet5, internal, generate
1, 8, 1
*Elset, elset=_PickedSet5, internal
1,
** Section: smp
*Solid Section, elset=_PickedSet5, material=smpmat
,
*End Part
**
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=smpelem1, part=Part-1
*End Instance
**
*Nset, nset=_PickedSet4, internal, instance=smpelem1, generate
1, 4, 1
*Elset, elset=_PickedSet4, internal, instance=smpelem1
1,
*Nset, nset=_PickedSet4b, internal, instance=smpelem1, generate
5, 8, 1
*Elset, elset=_PickedSet4b, internal, instance=smpelem1
1,
*Elset, elset=__PickedSurf5_S1, internal, instance=smpelem1
1,
*Nset, nset=_PickedSet9, internal, instance=smpelem1, generate
1, 8, 1
*Elset, elset=_PickedSet9, internal, instance=smpelem1
1,
*Surface, type=ELEMENT, name=_PickedSurf5, internal
__PickedSurf5_S1, S1
*Elset, elset=__PickedSurf6_S1, internal, instance=smpelem1
1,
*Surface, type=ELEMENT, name=_PickedSurf6, internal
__PickedSurf6_S1, S1
*End Assembly
**
** MATERIALS
**
*Material, name=Poly
*Conductivity
0.9,
*Density
1.,
*Elastic
1.46e+08, 0.45
*Material, name=smpmat
*Density
1,
*Conductivity
0.17,
*User Material, constants=14, type=MECHANICAL
146e6,38.1,14e9,44.2,0.42,0.003,58.2,0.112
38.7,521,35.4,348,328,11.6e-5
** 4.75e+06, 1.19e+06, 4.76e-06, 1.19e-06, 0.667
*DEPVAR
24
** INITIAL CONDITIONS
**
*INITIAL CONDITION, TYPE=TEMP
_PickedSet9, 348
**
** ----------------------------------------------------------------
**
** STEP: Step-1
**
*Step, name=Step-1,nlgeom,inc=100
*Coupled Temperature-displacement, creep=none, steady state
1., 1., 1e-05, 1.
*** Name: BC-1 Type: Displacement/Rotation
*Boundary
_PickedSet4, 1, 1
_PickedSet4, 2, 2
_PickedSet4, 3, 3
_PickedSet4, 4, 4
_PickedSet4, 5, 5
_PickedSet4, 6, 6
** LOADS
**
** Name: Load-2 Type: Pressure
*Dsload, op=NEW
_PickedSurf6, P, -5e5.
**
** BOUNDARY CONDITIONS
**
** Name: BC-2 Type: Temperature
*Boundary
_PickedSet9, 11, 11, 348
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
CF, RF, U
*Element Output, directions=YES
E, S, SDV, TEMP
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
** ----------------------------------------------------------------
**
** STEP: Step-2
**
*Step, name=Step-2,nlgeom,inc=100
*Coupled Temperature-displacement, creep=none, steady state
0.2, 1., 1e-05, 1.
**
** BOUNDARY CONDITIONS
**
** Name: BC-2 Type: Velocity/Angular velocity
*Boundary, op=NEW, type=VELOCITY
_PickedSet4b, 1, 1
_PickedSet4b, 2, 2
_PickedSet4b, 3, 3
_PickedSet4b, 4, 4
_PickedSet4b, 5, 5
_PickedSet4b, 6, 6
*** Name: BC-1 Type: Displacement/Rotation
*Boundary, op=NEW
_PickedSet4, 1, 1
_PickedSet4, 2, 2
_PickedSet4, 3, 3
_PickedSet4, 4, 4
_PickedSet4, 5, 5
_PickedSet4, 6, 6
** LOADS
**
** Name: Load-2 Type: Pressure
*Dsload, op=NEW
**_PickedSurf6, P, -12e5.
** BOUNDARY CONDITIONS
**
** Name: BC-2 Type: Temperature
*Boundary, op=NEW
_PickedSet9, 11, 11, 308
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
CF, RF, U, NT
*Element Output, directions=YES
E, S, SDV, TEMP
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
** ----------------------------------------------------------------
**
** STEP: Step-3
**
*Step, name=Step-3 ,nlgeom,inc=100
*Coupled Temperature-displacement, creep=none, steady state
0.3, 1., 1e-05, 1
*** Name: BC-2 Type: Velocity
*Boundary, op=NEW
** Name: BC-1 Type: Displacement/Rotation
*Boundary, op=NEW
_PickedSet4, 1, 1
_PickedSet4, 2, 2
_PickedSet4, 3, 3
_PickedSet4, 4, 4
_PickedSet4, 5, 5
_PickedSet4, 6, 6
** LOADS
**
** Name: Load-2 Type: Pressure
*Dsload, op=NEW
**_PickedSurf6, P, 0
** BOUNDARY CONDITIONS
**
** Name: BC-2 Type: Temperature
*Boundary, op=New
_PickedSet9, 11, 11, 308
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
CF, RF, U, NT
*Element Output, directions=YES
E, S, SDV, TEMP
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
** ----------------------------------------------------------------
**
** STEP: Step-4
**
*Step, name=Step-4,nlgeom,inc=100
*Coupled Temperature-displacement, creep=none, steady state
0.2, 1., 1e-05, 1
** Name: BC-1 Type: Displacement/Rotation
*Boundary, op=NEW
_PickedSet4, 1, 1
_PickedSet4, 2, 2
_PickedSet4, 3, 3
_PickedSet4, 4, 4
_PickedSet4, 5, 5
_PickedSet4, 6, 6
*** Name: BC-2 Type: Displacement/Rotation
*Boundary, op=NEW
** LOADS
**
** Name: Load-2 Type: Pressure
*Dsload, op=NEW
**_PickedSurf6, P, 0
** BOUNDARY CONDITIONS
**
** Name: BC-2 Type: Temperature
**Boundary, op=NEW
*Boundary, op=NEW
_PickedSet9, 11, 11, 348
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Node Output
CF, RF, U, NT
*Element Output, directions=YES
E, S, SDV, TEMP
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
SMP effects
Hi ,
I am also working on shape memory polymer. Curious to know if you got a fix for your model.
Did you have to change your umat.
Regards
Shouvik