Skip to main content

Abaqus UMAT implementation

Submitted by Mandar Kulkarni on

Hi,

I have a thermo-mechanical constitutive model for polymer and want to implement it in ABAQUS through two user subroutines namely UMAT and UEXPAN but I dont know how to use them together. Also, in UEXPAN, I need to pass in the strain tensor as information but I dont see 'STRAN' variable in the list of variables passed in to UEXPAN in Abaqus documentation. Does anyone know answer to this?

Thanks,

Mandar

Hi Frank,

Thanks for your help. I have a very specific question now. I was studying the .pdf file mentioned in your reply, particularly UMAT for Neo-Hookean material. On slide number L6.46 they define material jacobian for this hyperelastic material C_ijkl. I read through ABAQUS Theory manual section 4.6.1 but I could not understand where the last two terms in the expression for C_ijkl come from. I see that they are related to the DDSDDE definition on slide number L6.50 by terms TRBBAR (which is trace of left Cauchy deformation tensor) and EK which corresponds to bulk modulus of the material. But the theory manual defines material jacobian by deviding the strain into deviatoric and volumetric parts (Eq. 4.6.1-12 in ABAQUS Theory manual). So my question is how both deviatoric and volumetric behaviors can be included in one entry of DDSDDE array by just adding the terms. This might be a silly question and the anwser might just be the rearrangement of the terms. But I am finding it difficult to see it. I hope you can help me on this issue.

Thanks for your help and time in advance.

Mandar

Tue, 02/26/2013 - 23:31 Permalink

Hi, Mandar

I have the same question. According to  ABAQUS Theory manual section 4.6.1,  the material jacobian C_ijkl are consisted of deviatoric part C-s and volumetric part K, see the equation (4.6.1-12). However, based on these, i can not figure out the second to last term in the expression for C_ijkl of Neo-Hookean hyperelastic model on slide number L6.46.

So could you kindly tell me how you solve your problem?

Thank you.

zqctate

Thu, 05/16/2013 - 01:53 Permalink

I used to think that UEXPAN is a good way to implement thermal expasion cases. But after i tried. I would recommend you to integrate them into one UMAT subroutine. It could make your life easier.

You just need to modify your constitutive model by replacing your elactic strain with the total strain deducted by thermal strain.

 

Thu, 02/28/2013 - 19:24 Permalink

I just finished coding UMAT for hyperelastic material and tested it on one element for uniaxial and simple shear boundary conditions. It works fine for one element. But, if I try to run a simulation for more than one element, ABAQUS exits with an error. The message says that the error can be found in message file if the file exists and obviously the .msg file does not exist. Does any one have solution to this problem? or where should I look for the kind of error I am encountering with?

Thanks in advance.

Mandar

Sun, 03/10/2013 - 01:52 Permalink

Hi all,

Just a follow up question on the leading thread. My UMAT for hyperelastic material works fine if I use automatic time stepping provided by ABAQUS but if I use fixed time stepping it exits with an error message which says that FIXED TIME STEP IS TOO LARGE. I tried using fixed time stepping with time increment even smaller than the minimum increment size in autimatic time stepping but without any luck. Does anyone know why this occurs? I am just curious why fixed time increment wouldn't work.

Thanks for your help.

Mandar

Fri, 04/05/2013 - 20:50 Permalink

Please I am currently having problem developing a UMAT for transformation-induced swelling material behavior.

I need help.

 

ABBA

Tue, 06/11/2013 - 18:20 Permalink

Hello All,

I am new to creep modeling using Abaqus. I have been following a tutorial on creep modeling of pipe and modified the inputs based on my problem. I came across this error message while running the job with particular time step."Time increment required is less than the minimum specified"
When i reduced the step time, again a error was shown "Too many attempts made for this increment"

Is there any relation between the step time to be followed and the mesh size while doing creep analysis??
I had used 2 Steps for the analysis, (time in hours)

Step1. Pressure loading
 Time Period = 1
Incrementations:
 Max no: of increments =100
 Initial =1
 Minimum = 1e-5
 Maximum =1

Step 2: Creep
 Time Period = 70
Incrementations:
 Max no: of increments =100
 Initial =1
 Minimum = 1
 Maximum =70

Mesh Size is .001
Element type is C3D20R.

Please let me know your valuable suggestions.

Thanks for the help
Sat, 06/15/2013 - 07:01 Permalink

Here is a creep umat:

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,*

RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,*

TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,*

NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,*

DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)

C INCLUDE 'ABA_PARAM.INC'

C CHARACTER*80 MATERL

DIMENSION STRESS(NTENS),STATEV(NSTATV),

DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),

PROPS(NPROPS),COORDS(3),DROT(3,3),

DFGRD0(3,3),DFGRD1(3,3)

DOUBLE PRECISION J0,J,J2,KSI,NN,N,V,MC

C ELASTIC PROPERTIES

EMOD1=PROPS(1)

EMOD2=PROPS(2)

ENU=PROPS(3)

EG=PROPS(4)

THETA=PROPS(5)

KSI=PROPS(6)

ESI=PROPS(7)

ENTA=PROPS(8)

N=PROPS(9)

E0=PROPS(10)

SIG0=PROPS(11)

V=PROPS(12)

M=PROPS(13)

T0=PROPS(14)

EBULK=1/(1.0-ENU**2*EMOD2/EMOD1)

C

C ELASTIC STIFFNESS

C

IF (KSTEP.EQ.1) THEN

DO K1=1,NTENSDO K2=1,NTENSDDSDDE(K2,K1)=0

END DO

ENDDO

Q11=EMOD1*EBULK

Q12=ENU*EMOD2*EBULK

Q21=ENU*EMOD2*EBULK

Q22=EMOD2*EBULK

Q66=EG

DDSDDE(1,1)=Q11*COS(THETA)**4.0

+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0 +Q22*SIN(THETA)**4.0

DDSDDE(1,2)=(Q11+Q22-

4*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0+Q12*(SIN(THETA)**4.0+COS(THETA)*

*4.0)

DDSDDE(2,1)=DDSDDE(1,2)

DDSDDE(2,2)=Q11*SIN(THETA)**4.0

+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0 +Q22*COS(THETA)**4.0

DDSDDE(NTENS,1)=(Q11-Q12-2*Q66)*SIN(THETA)*COS(THETA)**3.0+ (Q12-

Q22+2*Q66)*SIN(THETA)**3.0*COS(THETA)

DDSDDE(1,NTENS)=DDSDDE(NTENS,1)

DDSDDE(NTENS,2)=(Q11-Q12-2*Q66)*SIN(THETA)**3.0*COS(THETA)+ (Q12-

Q22+2*Q66)*SIN(THETA)*COS(THETA)**3.0

DDSDDE(2,NTENS)=DDSDDE(NTENS,2)

DDSDDE(NTENS,NTENS)=(Q11+Q22-2*Q12-2*Q66)*SIN(THETA)**2*

*COS(THETA)**2.0+ Q66*(SIN(THETA)**4+COS(THETA)**4.0)

C

C1=DDSDDE(1,1)

C2=DDSDDE(1,2)

C3=DDSDDE(1,NTENS)

C4=DDSDDE(2,1)

C5=DDSDDE(2,2)

C6=DDSDDE(2,NTENS)

C7=DDSDDE(NTENS,1)

C8=DDSDDE(NTENS,2)

C9=DDSDDE(NTENS,NTENS)

C

C CALCULATE STRESS FROM ELASTIC STRAINS

C

DO K1=1,NTENSDO

K2=1,NTENSSTRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

END DO

END DO

C

STATEV(1)=STATEV(1)+DSTRAN(1)

STATEV(2)=STATEV(2)+DSTRAN(2)

STATEV(3)=STATEV(3)+DSTRAN(3)

A1=STATEV(1)

A2=STATEV(2)

A3=STATEV(3)

D1=C1

D2=C5

D3=C2C

END IF

C

C CALCULATE CREEP STRESS

IF (KSTEP.EQ.2) THEN

S11=STRESS(1)-(STRESS(1)+STRESS(2)+STRESS(3))/3.D0S22=STRESS(2)-

(STRESS(1)+STRESS(2)+STRESS(3))/3.D0

S22=STRESS(2)-(STRESS(1)+STRESS(2)+STRESS(3))/3.D0S33=STRESS(3)-

(STRESS(1)+STRESS(2)+STRESS(3))/3.D0

S12=STRESS(NTENS)D11=COS(THETA)**2.D0D22=SIN(THETA)**2.D0

D12=COS(THETA)*SIN(THETA)

J0=D11*S11+D22*S22+2.D0*D12*S12

J=D11*(S11**2.D0+S12**2.D0)+D22*(S12**2.D0+S22**2.D0)+D12*S12*(S11+S22)

J2=.5*(S11**2.D0+S22**2.D0+S33**2.D0)+S12**2Q1=JJ0**

2Q2=J0**2TEMP1=(ESI4.D0*ENTA)*(STRESS(1)+STRESS(2))**2/9.D0

TEMP2=(KSI)*Q1

TEMP3=(ESI-ENTA)*Q2TEMP4=3.0*(J2-TEMP2-TEMP3+TEMP1)

PHI=SQRT(TEMP4)/(E0)

TA11=2.D0*S11*D11+2.D0*D12*S12

TA22=2.D0*S22*D22+2.D0*D12*S12

TA12=S12*(D11+D22)+D12*(S11+S22)

TB11=2.D0*J0*D11TB22=2.D0*J0*D22

TB12=2.D0*J0*D12

TONE11=TA11-TB11

TONE22=TA22-TB22

TONE12=TA12-TB12

TTWO11=2.D0*J0*(D11-1.D0/3.D0)

TTWO22=2.D0*J0*(D22-1.D0/3.D0)

TTWO12=2.D0*J0*D12

TAO11=S11-KSI*TONE11-(ESIENTA)*

TTWO11+2.D0*(ESI4.D0*ENTA)*(STRESS(1)+STRESS(2))/9.D0

TAO22=S22-KSI*TONE22-(ESIENTA)*

TTWO22+2.D0*(ESI4.D0*ENTA)*(STRESS(1)+STRESS(2))/9.D0

TAO12=S12-KSI*TONE12-(ESI-ENTA)*TTWO12

RSTRAN11=3.D0*PHI**(N-1.D0)*TAO11/(2.D0*E0*SIG0*100)

RSTRAN22=3.D0*PHI**(N-1.D0)*TAO22/(2.D0*E0*SIG0*100)

RSTRAN12=3.D0*PHI**(N-1.D0)*TAO12/(2.D0*E0*SIG0*100)DC

RSTRAN11=RSTRAN11*DTIMEDC

RSTRAN22=RSTRAN22*DTIMEDC

RSTRAN12=RSTRAN12*DTIME

STATEV(4)=DCRSTRAN11+STATEV(4)

STATEV(5)=DCRSTRAN22+STATEV(5)

STATEV(6)=DCRSTRAN12+STATEV(6)

C

C DAMAGE calculation

C

I=STRESS(1)+STRESS(2)+STRESS(3)

I0=D11*STRESS(1)+2*D12*STRESS(NTENS)+D22*STRESS(2)

RE1=J2+.25*J0**2.0-J

IF (RE1.LE.0) THEN

RE1=0-RE1

END IF

NN=.5*(I-I0)+SQRT(RE1)

IF (NN.LE.0) THEN

NN=0.0

END IF

IF (Q1.LE.0) THEN

Q1=J0**2-J

END IF

S=SQRT(Q1)

DELTA=(NN+0.35*S)/(E0)

STATEV(11)=(1-DELTA**V*TIME(2)/T0)**(1/(1+M))

RDASTRAN11=3.0*PHI**(N-1.0)*TAO11/(2.0*E0*SIG0*100/(STATEV(11)**N)

DDAMSTRAN11=RDASTRAN11*DTIME

DDAMSTRAN22=RDASTRAN22*DTIME

DDAMSTRAN12=RDASTRAN12*DTIME

STATEV(12)=DDAMSTRAN11+STATEV(12)

STATEV(13)=DDAMSTRAN22+STATEV(13)

STATEV(14)=DDAMSTRAN12+STATEV(14)

C

C CALCUALTE UPDATED STRESS

C

DS1=C1*(DSTRAN(1)-DCRSTRAN11)

DSTRESS11=C1*(DSTRAN(1)-DCRSTRAN11)+C2*(DSTRAN(2)-DCRSTRAN22)+

C3*(DSTRAN(NTENS)-DCRSTRAN12)

DSTRESS22=C4*(DSTRAN(1)-DCRSTRAN11)+C5*(DSTRAN(2)-DCRSTRAN22)+

C6*(DSTRAN(NTENS)-DCRSTRAN12)

DSTRESS12=C7*(DSTRAN(1)-DCRSTRAN11)+C8*(DSTRAN(2)-DCRSTRAN22)+

C9*(DSTRAN(NTENS)-DCRSTRAN12)

STRESS(1)=DSTRESS11+STRESS(1)

STRESS(2)=DSTRESS22+STRESS(2)

STRESS(NTENS)=DSTRESS12+STRESS(NTENS)

C

C CALCULATE UPDATED JACOBIAN

DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K2,K1)=0

ENDDO

ENDDO

Q11=EMOD1*EBULK

Q12=ENU*EMOD2*EBULK

Q21=ENU*EMOD2*EBULK

Q22=EMOD2*EBULK

Q66=EG

DDSDDE(1,1)=Q11*COS(THETA)**4.0+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(TH

ETA)**2.0* +Q22*SIN(THETA)**4.0

DDSDDE(1,2)=(Q11+Q22-4*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0*

+Q12*(SIN(THETA)**4.0+COS(THETA)**4.0)

DDSDDE(2,1)=DDSDDE(1,2)

DDSDDE(2,2)=Q11*SIN(THETA)**4.0

+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0* +Q22*COS(THETA)**4.0

DDSDDE(NTENS,1)=(Q11-Q12-2*Q66)*SIN(THETA)*COS(THETA)**3.0+* (Q12-

Q22+2*Q66)*SIN(THETA)**3.0*COS(THETA)

DDSDDE(1,NTENS)=DDSDDE(NTENS,1)

DDSDDE(NTENS,2)=(Q11-Q12-2*Q66)*SIN(THETA)**3.0*COS(THETA)* (Q12-

Q22+2*Q66)*SIN(THETA)*COS(THETA)**3.0

DDSDDE(2,NTENS)=DDSDDE(NTENS,2)

DDSDDE(NTENS,NTENS)=(Q11+Q22-2*Q12-2*Q66)*SIN(THETA)**2*

*COS(THETA)**2.0+* Q66*(SIN(THETA)**4+COS(THETA)**4.0)

RETURN

END

INPUT FILE FOR PRESSURE VESSEL

*Heading

** Job name: Job-1 Model name: Model-1

*Node

1, 0., 20., 0

.……

12445, 19.9938431, -20., 0.496222973

*Element, type=S8R,ELSET=CREEP

1, 1, 291, 399, 4, 4171, 4172, 4173, 4174

*Nset, nset=B1

3, 164, ……., 12445

Surface, type=ELEMENT, name=vessel

creep, SNEG

*shell SECTION,ELSET=CREEP,MATERIAL=com1

0.4

*TRANSVERSE SHEAR STIFFNESS

2230, 2230, 0

** MATERIALS-1

*MATERIAL,NAME=COM1

*USER MATERIAL,CONSTANTS=16,TYPE=MECHANICAL

44e4, 7.3e3, .284, 2.7e3,1.047, 0.57, 0.64, 0.1,

6.5, 46., 50., 10.6, 7.125, 12,2.7e3,2.7e3

*DEPVAR

14

*INITIALCONDITIONS,TYPE=SOLUTION

CREEP,0.0,0.0,0.0,0.0,0.0,0.0,0.0

0.0,0.0,0.0,0.0,0.0,0.0,0.0

*BoundaryB1, 2, 2

*STEP,INC=30

PRESCRIBED TENSILE STRESS

*STATIC

1.E-7,1.E-6

*Dsload

vessel, P, 0.5

*EL PRINT,FREQ=1

S,

SDV1,

SDV4

*END STEP

*STEP,INC=100000

CREEP STEP

*VISCO,CETOL=1E-4

1.E-6,10

** OUTPUT REQUESTS

*PRINT,FREQ=1,RESIDUAL=YES

*EL PRINT,FREQ=1

SDV

*OUTPUT,FIELD,FREQ =1

*ELEMENT OUTPUT

S,E,CE,

SDV

NODE OUTPUT

U

*END STEP

Fri, 07/05/2013 - 19:55 Permalink

Hi all

 I have written an UMAT  to simulate pseudoelasticity effect ( With abaqus CAE) but I am against difficulties.

 

 Does somebody help me to tell me steps to follow i.e :

To create right steps and loadings

Fri, 07/05/2013 - 18:22 Permalink