In this post, we explore how to integrate a neural network in a user material subroutine (UMAT) for CalculiX. The main idea is that the material’s plastic behavior is learned by a neural network and embedded into the UMAT code, in combination with a classic linear elastic model.
We’ll use a UMAT example with the following characteristics:
- Linear isotropic elasticity
- Von Mises plasticity with isotropic hardening
The neural network in this example takes two normalized inputs:
- The equivalent plastic strain
- The trial Von Mises, a calculated stress state based on the current strain increment.
Normalization is done using mean and standard deviation values defined as:
x(1) = (EQPLAS - inputmu(1)) / inputstd(1) x(2) = (SMISES - inputmu(2)) / inputstd(2)
The example contains a simple feedforward neural network with one hidden layer of 10 neurons using ReLU activation. The trained weights and biases are manually stored in Fortran arrays:
do i = 1, NHIDDEN
do j = 1,2
x1(i) = x1(i) + x(j) * L1(i,j)
end do
x1(i) = max(0.0, x1(i) + b1(i))
end do
Then, the output is computed and denormalized:
y = y * outputstd + outputmu
The following updates are performed after inference:
- Stress tensor
stre(6)is recalculated using the plastic flow direction - The internal variable
EQPLASis updated and stored inxstate
We also calculate the consistent tangent stiffness matrix stiff(21) using effective elastic and hardening moduli derived at the updated point.
- This UMAT is only valid for solid elements (not plane stress elements).
- The neural network must be trained offline, using clean, normalized datasets.
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
INTEGER KEWTON,NVALUE,NHIDDEN
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)
real*8 L1(10,2),L2(10)
real*8 b1(10),b2
real*8 x(2),x1(10),y,inputmu(2),inputstd(2),outputmu,outputstd
DIMENSION xmatrix(10,2)
! Initialize the 4x3 matrix with the specified values
DATA L1 /0.122878104,0.50649923,
& -0.33488736,-0.02137515,
& -0.034247752,0.176506,
& -0.010768869,-0.111633584,
& 0.031866815,-0.059505835,
& -0.19340435,-0.79933566,
& 0.528451,-0.61828786,
& -0.64486074,0.28845328,
& -0.38559982,-0.26867574,
& 0.52625835,-0.20963475/
DATA B1 /-0.052253507,-0.20429772,
& 0.13519116,1.0849684,
& 1.1447347,0.79870075,
& 0.6741626,0.63055074,
& 0.95799935,0.44531068/
DATA L2 /-0.024199534,-0.25845042,
& 0.40025997,-0.20823033,
& -0.23183896,0.3965015,
& -0.21810858,-0.25155863,
& 0.41102442,-0.13165537/
DATA B2 /0.09221251/
DATA inputmu /5.0861180e-01, 5.0296434e+04/
DATA inputstd /2.8188679e-01, 2.8522984e+04/
DATA outputmu /0.21653831/
DATA outputstd /0.12359327/
NHIDDEN = 10
C
C Elastic matrix
C
C
EMOD = elconloc(1)
ENU=MIN(elconloc(2), ENUMAX)
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
x(1) = (EQPLAS - inputmu(1)) / inputstd(1)
x(2) = (SMISES - inputmu(2)) / inputstd(2)
x1 = 0.0
do i = 1, NHIDDEN
do j =1,2
x1(i) = x1(i) + x(j) * L1(i,j)
end do
x1(i) = x1(i) + b1(i)
end do
do i = 1, 10
x1(i) = max(0.0, x1(i))
end do
y = 0
do i = 1, NHIDDEN
y = y + L2(i) * x1(i)
end do
y = y + b2
y = y * outputstd + outputmu
DEQPL = y
CALL UHARD(SYIELD,HARD,EQPLAS+DEQPL,elconloc(3),NVALUE)
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
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)
! write(*,*) 'stiff',(stiff(i),i=1,21)
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