User login

You are here

Inaccurate UMAT result in SMP model

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
 

 

 

 

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

Subscribe to Comments for "Inaccurate UMAT result in SMP model"

Recent comments

More comments

Syndicate

Subscribe to Syndicate