Converting user material subroutines from ABAQUS to CALCULIX can significantly boost flexibility in open-source finite element simulations. Thanks to the ABAQUS UMAT examples shared on GitHub and similar platforms, it is often easier to start from an existing model than to build one from scratch. In this post, we explore how to adapt one of these ABAQUS UMATs for use in CALCULIX, covering the key differences and steps in the conversion process.
Adapting an ABAQUS User Material (UMAT) subroutine for CALCULIX involves several critical steps to ensure compatibility and functionality. Building upon the foundational UMAT structure discussed in the article «Developing User Material Subroutines in CALCULIX: 1 A Comprehensive Guide with UMAT Structure Examples in ABAQUS,» the following key steps outline the process:
Familiarize yourself with the components and workflow of the ABAQUS UMAT, as detailed in the referenced article. This includes input and output variables, as well as the general structure of the subroutine.
Examine the differences in subroutine interfaces between ABAQUS and CALCULIX. While both use Fortran-based UMATs, the specific arguments and their ordering may differ. Ensure that the subroutine header and variable declarations are adjusted to match CALCULIX’s expectations. Adapt the input handling to align with CALCULIX’s data structures.
CALCULIX input variables
! INPUT: ! ! amat material name ! iel element number ! iint integration point number ! ! kode material type (-100-#of constants entered ! under *USER MATERIAL): can be used for materials ! with varying number of constants ! ! elconloc(*) user defined constants defined by the keyword ! card *USER MATERIAL (actual # = ! -kode-100), interpolated for the ! actual temperature t1l ! ! emec(6) Lagrange mechanical strain tensor (component order: ! 11,22,33,12,13,23) at the end of the increment ! (thermal strains are subtracted) ! emec0(6) Lagrange mechanical strain tensor at the start of the ! increment (thermal strains are subtracted) ! beta(6) residual stress tensor (the stress entered under ! the keyword *INITIAL CONDITIONS,TYPE=STRESS) ! ! xokl(3,3) deformation gradient at the start of the increment ! voj Jacobian at the start of the increment ! xkl(3,3) deformation gradient at the end of the increment ! vj Jacobian at the end of the increment ! ! ithermal 0: no thermal effects are taken into account ! >0: thermal effects are taken into account (triggered ! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) ! t1l temperature at the end of the increment ! dtime time length of the increment ! time step time at the end of the current increment ! ttime total time at the start of the current step ! ! icmd not equal to 3: calculate stress and stiffness ! 3: calculate only stress ! ielas 0: no elastic iteration: irreversible effects ! are allowed ! 1: elastic iteration, i.e. no irreversible ! deformation allowed ! ! mi(1) max. # of integration points per element in the ! model ! nstate_ max. # of state variables in the model ! ! xstateini(nstate_,mi(1),# of elements) ! state variables at the start of the increment ! xstate(nstate_,mi(1),# of elements) ! state variables at the end of the increment ! ! stre(6) Piola-Kirchhoff stress of the second kind ! at the start of the increment ! ! iorien number of the local coordinate axis system ! in the integration point at stake (takes the value ! 0 if no local system applies) ! pgauss(3) global coordinates of the integration point ! orab(7,*) description of all local coordinate systems. ! If a local coordinate system applies the global ! tensors can be obtained by premultiplying the local ! tensors with skl(3,3). skl is determined by calling ! the subroutine transformatrix: ! call transformatrix(orab(1,iorien),pgauss,skl)
Retrieve the material properties, the elastic modulus and the Poisson ratio
EMOD = elconloc(1)
ENU=elconloc(2)
Compute the strain increment as the difference between the final total strain and the initial total strain. CALCULIX does not pass the strain increment directly as ABAQUS does—it must be computed manually:
C
C Strain increment
C
DO K1=1,6
DSTRAN(K1)=emec(K1)-emec0(K1)
ENDDO
Retrieve the equivalent plastic strain, stored in a user variable.
EQPLAS = xstateini(1,iint,iel)
Importantly, in CALCULIX the tangent stiffness must be returned in vectorial form (Voigt notation, 6×6 for 3D problems). Update accordingly, ensuring proper ordering of shear components.
stiff(1)=G(1,1)
stiff(2)=G(1,2)
stiff(3)=G(2,2)
stiff(4)=G(1,3)
stiff(5)=G(2,3)
stiff(6)=G(3,3)
stiff(7)=G(1,4)
stiff(8)=G(2,4)
stiff(9)=G(3,4)
stiff(10)=G(4,4)
stiff(11)=G(1,5)
stiff(12)=G(2,5)
stiff(13)=G(3,5)
stiff(14)=G(4,5)
stiff(15)=G(5,5)
stiff(16)=G(1,6)
stiff(17)=G(2,6)
stiff(18)=G(3,6)
stiff(19)=G(4,6)
stiff(20)=G(5,6)
stiff(21)=G(6,6)
Importantly, in CALCULIX the tangent stiffness must be returned in vectorial form. Update accordingly, ensuring proper ordering of shear components.
DO K1=1,3
stre(K1)=FLOW(K1)*SYIELD+SHYDRO
END DO
DO K1=4,6
stre(K1)=FLOW(K1)*SYIELD
END DO
The equivalent plastic strain (often noted as PEEQ in ABAQUS) must be manually stored and updated in the user variable.
EQPLAS=EQPLAS+DEQPL
xstate(1,iint,iel) = EQPLAS
home = $(HOME)/ARPACK to your ARPACK directory
PLAT = SUN4 to PLAT = linux
FC = f77 to FC = gfortran
FFLAGS = subroutine umat_user(amat,iel,iint,kode,elconloc,emec,emec0,
& beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime,
& icmd,ielas,mi,nstate_,xstateini,xstate,stre,stiff,
& iorien,pgauss,orab,pnewdt,ipkon)
!
! calculates stiffness and stresses for a user defined material
! law
!
! icmd=3: calculates stress at mechanical strain
! else: calculates stress at mechanical strain and the stiffness
! matrix
!
! INPUT:
!
! amat material name
! iel element number
! iint integration point number
!
! kode material type (-100-#of constants entered
! under *USER MATERIAL): can be used for materials
! with varying number of constants
!
! elconloc(*) user defined constants defined by the keyword
! card *USER MATERIAL (actual # =
! -kode-100), interpolated for the
! actual temperature t1l
!
! emec(6) Lagrange mechanical strain tensor (component order:
! 11,22,33,12,13,23) at the end of the increment
! (thermal strains are subtracted)
! emec0(6) Lagrange mechanical strain tensor at the start of the
! increment (thermal strains are subtracted)
! beta(6) residual stress tensor (the stress entered under
! the keyword *INITIAL CONDITIONS,TYPE=STRESS)
!
! xokl(3,3) deformation gradient at the start of the increment
! voj Jacobian at the start of the increment
! xkl(3,3) deformation gradient at the end of the increment
! vj Jacobian at the end of the increment
!
! ithermal 0: no thermal effects are taken into account
! >0: thermal effects are taken into account (triggered
! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE)
! t1l temperature at the end of the increment
! dtime time length of the increment
! time step time at the end of the current increment
! ttime total time at the start of the current step
!
! icmd not equal to 3: calculate stress and stiffness
! 3: calculate only stress
! ielas 0: no elastic iteration: irreversible effects
! are allowed
! 1: elastic iteration, i.e. no irreversible
! deformation allowed
!
! mi(1) max. # of integration points per element in the
! model
! nstate_ max. # of state variables in the model
!
! xstateini(nstate_,mi(1),# of elements)
! state variables at the start of the increment
! xstate(nstate_,mi(1),# of elements)
! state variables at the end of the increment
!
! stre(6) Piola-Kirchhoff stress of the second kind
! at the start of the increment
!
! iorien number of the local coordinate axis system
! in the integration point at stake (takes the value
! 0 if no local system applies)
! pgauss(3) global coordinates of the integration point
! orab(7,*) description of all local coordinate systems.
! If a local coordinate system applies the global
! tensors can be obtained by premultiplying the local
! tensors with skl(3,3). skl is determined by calling
! the subroutine transformatrix:
! call transformatrix(orab(1,iorien),pgauss,skl)
!
!
! OUTPUT:
!
! xstate(nstate_,mi(1),# of elements)
! updated state variables at the end of the increment
! stre(6) Piola-Kirchhoff stress of the second kind at the
! end of the increment
! stiff(21): consistent tangent stiffness matrix in the material
! frame of reference at the end of the increment. In
! other words: the derivative of the PK2 stress with
! respect to the Lagrangian strain tensor. The matrix
! is supposed to be symmetric, only the upper half is
! to be given in the same order as for a fully
! anisotropic elastic material (*ELASTIC,TYPE=ANISO).
! Notice that the matrix is an integral part of the
! fourth order material tensor, i.e. the Voigt notation
! is not used.
! pnewdt to be specified by the user if the material
! routine is unable to return the stiffness matrix
! and/or the stress due to divergence within the
! routine. pnewdt is the factor by which the time
! increment is to be multiplied in the next
! trial and should exceed zero but be less than 1.
! Default is -1 indicating that the user routine
! has converged.
! ipkon(*) ipkon(iel) points towards the position in field
! kon prior to the first node of the element's
! topology. If ipkon(iel) is smaller than 0, the
! element is not used.
!
* implicit none
!
character*80 amat
!
integer ithermal(*),icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien,
& i,K1,K2,N
!
real*8 elconloc(*),stiff(21),emec(6),emec0(6),beta(6),stre(6),
& vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*),
& time,ttime
!
real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
!
C ----------------------------------------------------------------
C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC MISES PLASTICITY
C CANNOT BE USED FOR PLANE STRESS
C ----------------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3..) - YIELD AND HARDENING DATA
C CALLS UHARD FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
C ----------------------------------------------------------------
C ----------------------------------------------------------------
C UMAT FOR ISOTROPIC ELASTICITY
C CANNOT BE USED FOR PLANE STRESS
C ----------------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C ----------------------------------------------------------------
C
real*8 EMOD,NU,EBULK3,EG2,EG3,ELAM,G(6,6),DSTRAN(6)
real*8 EQPLAS,SMISES,SYIEL0,HARD,SHYDRO,FLOW(6),SYIELD,RHS
real*8 EFFG,EFFG2,EFFG3,EFFLAM,EFFHRD
real*8 ZERO,ONE,TWO,THREE,SIX,ENUMAX,TOLER,ENU
INTEGER KEWTON
PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,
1 ENUMAX=.4999D0, NEWTON=10, TOLER=1.0D-6)
!
C
C Elastic matrix
C
C
EMOD = elconloc(1)
ENU=elconloc(2)
if (ENU.GT.ENUMAX) THEN
ENU = ENUMAX
ENDIF
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
C
C ELASTIC STIFFNESS
C
DO K1=1,6
DO K2=1,6
G(K1,K2)=0.D0
ENDDO
ENDDO
DO K1=1, 3
DO K2=1, 3
G(K2, K1)=ELAM
END DO
G(K1, K1)=EG2+ELAM
END DO
DO K1=4, 6
G(K1 ,K1)=EG
END DO
C
C Strain increment
C
DO K1=1,6
DSTRAN(K1)=emec(K1)-emec0(K1)
ENDDO
C
C DEFINING Initial STRESS AND STRAIN RELATION
C
stre(1)=stre(1)+G(1,1)*DSTRAN(1)+G(1,2)*DSTRAN(2)
& +G(1,3)*DSTRAN(3)
stre(2)=stre(2)+G(2,1)*DSTRAN(1)+G(2,2)*DSTRAN(2)
& +G(2,3)*DSTRAN(3)
stre(3)=stre(3)+G(3,1)*DSTRAN(1)+G(3,2)*DSTRAN(2)
& +G(3,3)*DSTRAN(3)
stre(4)=stre(4)+G(4,4)*DSTRAN(4)
stre(5)=stre(5)+G(5,5)*DSTRAN(5)
stre(6)=stre(6)+G(6,6)*DSTRAN(6)
C
C NOW THE STATE VARAIBLES WILL BE UPDATED
C
EQPLAS = xstateini(1,iint,iel)
C
C CALCULATE EQUIVALENT VON MISES STRESS
C
SMISES=(stre(1)-stre(2))**2+(stre(2)-stre(3))**2
1 +(stre(3)-stre(1))**2
DO K1=4,6
SMISES=SMISES+SIX*stre(K1)**2
END DO
SMISES=SQRT(SMISES/TWO)
C
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
C
NVALUE=(-kode-100)/2-1
c
CALL UHARD(SYIEL0, HARD, EQPLAS, elconloc(3),NVALUE)
C
C DETERMINE IF ACTIVELY YIELDING
C
IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN
C
C ACTIVELY YIELDING
C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS
C CALCULATE THE FLOW DIRECTION
C
SHYDRO=(stre(1)+stre(2)+stre(3))/THREE
DO K1=1,3
FLOW(K1)=(stre(K1)-SHYDRO)/SMISES
END DO
DO K1=4, 6
FLOW(K1)=stre(K1)/SMISES
END DO
C
C SOLVE FOR EQUIVALENT VON MISES STRESS
C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION
C
SYIELD=SYIEL0
DEQPL=ZERO
DO KEWTON=1, NEWTON
RHS=SMISES-EG3*DEQPL-SYIELD
DEQPL=DEQPL+RHS/(EG3+HARD)
CALL UHARD(SYIELD,HARD,EQPLAS+DEQPL,elconloc(3),NVALUE)
IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10
END DO
C
C WRITE WARNING MESSAGE TO THE .MSG FILE
C
WRITE(*,*) 'PLASTICITY ALGORITHM DID NOT CONVERGE'
10 CONTINUE
C
C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND
C EQUIVALENT PLASTIC STRAIN
C
DO K1=1,3
stre(K1)=FLOW(K1)*SYIELD+SHYDRO
END DO
DO K1=4,6
stre(K1)=FLOW(K1)*SYIELD
END DO
EQPLAS=EQPLAS+DEQPL
C
C CALCULATE PLASTIC DISSIPATION
C
xstate(1,iint,iel) = EQPLAS
c SPD=DEQPL*(SYIEL0+SYIELD)/TWO
C
C FORMULATE THE JACOBIAN (MATERIAL TANGENT)
C FIRST CALCULATE EFFECTIVE MODULI
C
EFFG=EG*SYIELD/SMISES
EFFG2=TWO*EFFG
EFFG3=THREE/TWO*EFFG2
EFFLAM=(EBULK3-EFFG2)/THREE
EFFHRD=EG3*HARD/(EG3+HARD)-EFFG3
DO K1=1, 3
DO K2=1, 3
G(K2, K1)=EFFLAM
END DO
G(K1, K1)=EFFG2+EFFLAM
END DO
DO K1=4, 6
G(K1, K1)=EFFG
END DO
DO K1=1, 6
DO K2=1, 6
G(K2, K1)=G(K2, K1)+EFFHRD*FLOW(K2)*FLOW(K1)
END DO
END DO
ENDIF
C
C
!
if(icmd.ne.3) then
!
! insert here code to calculate the stiffness matrix
!
stiff(1)=G(1,1)
stiff(2)=G(1,2)
stiff(3)=G(2,2)
stiff(4)=G(1,3)
stiff(5)=G(2,3)
stiff(6)=G(3,3)
stiff(7)=G(1,4)
stiff(8)=G(2,4)
stiff(9)=G(3,4)
stiff(10)=G(4,4)
stiff(11)=G(1,5)
stiff(12)=G(2,5)
stiff(13)=G(3,5)
stiff(14)=G(4,5)
stiff(15)=G(5,5)
stiff(16)=G(1,6)
stiff(17)=G(2,6)
stiff(18)=G(3,6)
stiff(19)=G(4,6)
stiff(20)=G(5,6)
stiff(21)=G(6,6)
endif
!
return
end
SUBROUTINE UHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
C
REAL*8 TABLE(2, NVALUE),SYIELD,HARD,EQPLAS
REAL*8 EQPL0,EQPL1,DEQPL,SYIEL0,SYIEL1,DSYIEL
INTEGER NVALUE
C
PARAMETER(ZERO=0.D0)
C
C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
C
SYIELD=TABLE(1, NVALUE)
HARD=ZERO
c
C IF MORE THAN ONE ENTRY, SEARCH TABLE
C
IF(NVALUE.GT.1) THEN
DO K1=1, NVALUE-1
EQPL1=TABLE(2,K1+1)
IF(EQPLAS.LT.EQPL1) THEN
EQPL0=TABLE(2, K1)
C
C CURRENT YIELD STRESS AND HARDENING
C
DEQPL=EQPL1-EQPL0
SYIEL0=TABLE(1, K1)
SYIEL1=TABLE(1, K1+1)
DSYIEL=SYIEL1-SYIEL0
if (DEQPL.NE.ZERO) THEN
HARD=DSYIEL/DEQPL
else
HARD = ZERO
ENDIF
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD
GOTO 10
ENDIF
END DO
10 CONTINUE
ENDIF
RETURN
END