Hi,
I intend to simulate a composite part subjected to cool down, then mechanical load in a subsequent step.
I developed a UMAT code for this purpose, but the results are different from what ABAQUS calculates using its built in material model. If the thermal simulation is turned off, the results are identical, thus the rest of the formulation in UMAT apart from the thermal part is correct.
Why is there such a significant difference between results?
I do not open a zip. Post individual text files.
Subscribe to and seek assistance from the
ABAQUS mailing list
[url] https://groups.yahoo.com/group/ABAQUS [/url]
Search the archive of the list before posting in it.
The list does not accept attachments.
For an intro to subroutines get the file
[url] http://imechanica.org/files/Writing [/url] User Subroutines with ABAQUS.pdf
Good luck
Frank
I attach the UMAT in text file, the input files are too large.
Yes, I am aware of that document and the thermal part of my code is based on that.
Are the discrepancies significant ?
I advise you to post in the list. It does not accept attachments.
Good luck
Frank
Yes they are significant. Although the picture size got reduced, in my posted example, at the end of thermal step there is a 6% difference (34.6 vs 36.8 MPa), by the end of the mechanical step there is a 9% difference (368 vs 339 MPa), tried with other models too. It doesnt seem like a difference because of rounding or similar numerical issue, rather like a problem with the constitutive equation or simply because ABAQUS does not use this formulation, even if the documentation says so.
Im waiting for approval for the Yahoo mailing list.
The code:
[FONT=courier new]C********************************************************************
C*****************************************************************
C UMAT for elastic orthotropic material behaviour with thermal residual
C stress
C*****************************************************************
C********************************************************************
C
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,JSTEP,KINC)
C
INCLUDE ABA_PARAM.INC
C
CHARACTER*80 CMNAME
C
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),
4 JSTEP(4), Q(6,6), DTHERM(6)
C
C********************************************************************
C Initialisation
C********************************************************************
C
E1 = PROPS(1)
E2 = PROPS(2)
E3 = PROPS(3)
NU12 = PROPS(4)
NU13 = PROPS(5)
NU23 = PROPS(6)
G12 = PROPS(7)
G13 = PROPS(8)
G23 = PROPS(9)
ALPHA11 = PROPS(10)
ALPHA22 = PROPS(11)
ALPHA33 = PROPS(12)
C
NU21 = (E2/E1)*NU12
NU32 = (E3/E2)*NU23
NU31 = (E3/E1)*NU13
C
C********************************************************************
C Calculate Q matrix
C********************************************************************
C
QQ = 1.0D0/(1.0D0-
1 NU12*NU21-NU23*NU32-NU31*NU13-2.0D0*NU21*NU32*NU13)
C
Q(1,1) = E1*(1.0D0-NU23*NU32)*QQ
Q(2,2) = E2*(1.0D0-NU13*NU31)*QQ
Q(3,3) = E3*(1.0D0-NU12*NU21)*QQ
Q(4,4) = G12
Q(5,5) = G13
Q(6,6) = G23
Q(1,2) = E1*(NU21+NU31*NU23)*QQ
Q(1,3) = E1*(NU31+NU21*NU32)*QQ
Q(2,3) = E2*(NU32+NU12*NU31)*QQ
Q(1,4) = 0.0D0
Q(1,5) = 0.0D0
Q(1,6) = 0.0D0
Q(2,4) = 0.0D0
Q(2,5) = 0.0D0
Q(2,6) = 0.0D0
Q(3,4) = 0.0D0
Q(3,5) = 0.0D0
Q(3,6) = 0.0D0
Q(4,5) = 0.0D0
Q(4,6) = 0.0D0
Q(5,6) = 0.0D0
C
C********************************************************************
C Calculate thermal strain
C********************************************************************
C
DTHERM(1) = ALPHA11*DTEMP
DTHERM(2) = ALPHA22*DTEMP
DTHERM(3) = ALPHA33*DTEMP
DTHERM(4) = 0.0D0
DTHERM(5) = 0.0D0
DTHERM(6) = 0.0D0
C
C********************************************************************
C Calculate stress matrix
C********************************************************************
C
DO K1=1,NTENS
DO K2=1,NTENS
STRESS(K2) = STRESS(K2)+Q(K2,K1)*(DSTRAN(K2)-DTHERM(K2))
END DO
END DO
C
C********************************************************************
C Calculate stiffness matrix
C********************************************************************
C
DO K1=1,NTENS
DO K2=1,NTENS
DDSDDE(K2,K1) = 0.0D0
END DO
END DO
C
DO K1=1,NTENS
DO K2=1,NTENS
DDSDDE(K2,K1) = Q(K2,K1)
END DO
END DO
C
RETURN
END [/FONT]
Why are you not taking care of DROT.
I think DDSDDE should be rotated Q.
Please help me in this as I am stuck to this point.