Developing User Material Subroutines in CALCULIX: 1 A Comprehensive Guide with UMAT Structure Examples in ABAQUS

Deja un comentario / Material modelling

This blog initiates a series of post about creating user material subroutines, with illustrative examples. Finite element software, as ABAQUS or CALCULIX supports advanced material modeling through user-defined subroutines, specifically the User Material (UMAT). These subroutines allow for defining complex material behaviors not included in the material models library. In this article, we provide an overview of the UMAT structure and demonstrate its use with a simple example. The series focuses on CALCULIX and includes the following parts:

  • Explaining the UMAT structure in ABAQUS.
  • Compiling Calculix
  • Transforming an ABAQUS UMAT to a CALCULIX UMAT.
  • Elastoplastic subroutine example in CALCULIX
  • Using a neural network as a user material subroutine

Structure of the UMAT Subroutine in ABAQUS

The UMAT subroutine is coded in Fortran and allows users to define material constitutive equations with a high degree of customization. When implementing a UMAT, users interact with various arrays and variables provided by ABAQUS to define stresses, strains, and other material properties.

Key Components of the UMAT Structure

Input variables

ABAQUS passes several variables to UMAT, including:

  • STRESS: Stress tensor, which is updated by the UMAT
  • STRAIN: Strain tensor for the current increment.
  • STATEV: Array holding state-dependent variables, useful for history-dependent materials.
  • DSTRAN: Increment in strain.
  • TIME: Time array containing the total time and current increment time.
  • DTIME: Current time increment size.
  • PROPS: Material properties, defined in the ABAQUS material properties section.
  • COORDS: Coordinates of the current integration point.

Output variables

The subroutine must output updated values of the following:

  • STRESS: Updated stress tensor at the integration point.
  • DDSDDE: Tangent modulus, which ABAQUS uses for the iterative solution process. This is the Jacobian matrix of the stress-strain relationship.
  • STATEV: Updated state variables, which can hold values from previous steps to incorporate history effects.

Workflow of a UMAT

  1. Read material properties: Using the PROPS array to retrieve properties.
  2. Compute stress increment: Based on the strain increment (DSTRAN) and any history effects (in STATEV).
  3. Update stress and state variables: Store the updated values of stress and any state-dependent variables.
  4. Compute tangent modulus (DDSDDE): Required by ABAQUS to ensure convergence in nonlinear analyses.

Example 1: Simple Linear Elastic UMAT

This example demonstrates a basic linear elastic, isotropic material implemented in a UMAT subroutine.

SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, 
     1       SCD, RPL, DDSDDT, DRPLDE, DRPLDT, STRAN, 
     2       DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF, 
     3       DPRED, CMNAME, NDI, NSHR, NTENS, NSTATV, 
     4       PROPS, NPROPS, COORDS, DROT, PNEWDT, 
     5       CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, 
     6       KSPT, KSTEP, KINC)
     
      INCLUDE 'ABA_PARAM.INC'
      
      DOUBLE PRECISION STRESS(6), DDSDDE(6,6), STATEV(*), 
     1       STRAN(6), DSTRAN(6), PROPS(*), TIME(2)
      INTEGER NDI, NSHR, NTENS, NSTATV, NPROPS
      
      DOUBLE PRECISION YOUNG, POISSON, LAMBDA, MU
      DOUBLE PRECISION DELTA(6,6)
      INTEGER I, J
      
      YOUNG = PROPS(1)
      POISSON = PROPS(2)
      
      LAMBDA = (YOUNG * POISSON) / ((1 + POISSON) * (1 - 2 * POISSON))
      MU = YOUNG / (2 * (1 + POISSON))
      
      DO I = 1, 6
         DO J = 1, 6
            DELTA(I,J) = 0.0
         END DO
         DELTA(I,I) = 1.0
      END DO
      
      ! Stress update using Hooke's law
      STRESS(1) = STRESS(1) + (LAMBDA + 2 * MU) * DSTRAN(1) + LAMBDA * DSTRAN(2) + LAMBDA * DSTRAN(3)
      STRESS(2) = STRESS(2) + LAMBDA * DSTRAN(1) + (LAMBDA + 2 * MU) * DSTRAN(2) + LAMBDA * DSTRAN(3)
      STRESS(3) = STRESS(3) + LAMBDA * DSTRAN(1) + LAMBDA * DSTRAN(2) + (LAMBDA + 2 * MU) * DSTRAN(3)
      STRESS(4) = STRESS(4) + MU * DSTRAN(4)
      STRESS(5) = STRESS(5) + MU * DSTRAN(5)
      STRESS(6) = STRESS(6) + MU * DSTRAN(6)
      
      ! Tangent modulus for isotropic linear elastic material
      DDSDDE(1,1) = LAMBDA + 2 * MU
      DDSDDE(1,2) = LAMBDA
      DDSDDE(1,3) = LAMBDA
      DDSDDE(2,1) = LAMBDA
      DDSDDE(2,2) = LAMBDA + 2 * MU
      DDSDDE(2,3) = LAMBDA
      DDSDDE(3,1) = LAMBDA
      DDSDDE(3,2) = LAMBDA
      DDSDDE(3,3) = LAMBDA + 2 * MU
      DDSDDE(4,4) = MU
      DDSDDE(5,5) = MU
      DDSDDE(6,6) = MU
      
      RETURN
      END

Explanation of the Code

Define Constants: The YOUNG and POISSON values are extracted from PROPS.

      YOUNG = PROPS(1)
      POISSON = PROPS(2)

Calculate Lame Constants: LAMBDA and MU are derived from Young’s modulus and Poisson’s ratio.

      LAMBDA = (YOUNG * POISSON) / ((1 + POISSON) * (1 - 2 * POISSON))
      MU = YOUNG / (2 * (1 + POISSON))

Stress Update: Using Hooke’s law, the stress tensor is updated.

      ! Stress update using Hooke's law
      STRESS(1) = STRESS(1) + (LAMBDA + 2 * MU) * DSTRAN(1) + LAMBDA * DSTRAN(2) + LAMBDA * DSTRAN(3)
      STRESS(2) = STRESS(2) + LAMBDA * DSTRAN(1) + (LAMBDA + 2 * MU) * DSTRAN(2) + LAMBDA * DSTRAN(3)
      STRESS(3) = STRESS(3) + LAMBDA * DSTRAN(1) + LAMBDA * DSTRAN(2) + (LAMBDA + 2 * MU) * DSTRAN(3)
      STRESS(4) = STRESS(4) + MU * DSTRAN(4)
      STRESS(5) = STRESS(5) + MU * DSTRAN(5)
      STRESS(6) = STRESS(6) + MU * DSTRAN(6)

Tangent Modulus: The constant elastic matrix is defined in DDSDDE.

      ! Tangent modulus for isotropic linear elastic material
      DDSDDE(1,1) = LAMBDA + 2 * MU
      DDSDDE(1,2) = LAMBDA
      DDSDDE(1,3) = LAMBDA
      DDSDDE(2,1) = LAMBDA
      DDSDDE(2,2) = LAMBDA + 2 * MU
      DDSDDE(2,3) = LAMBDA
      DDSDDE(3,1) = LAMBDA
      DDSDDE(3,2) = LAMBDA
      DDSDDE(3,3) = LAMBDA + 2 * MU
      DDSDDE(4,4) = MU
      DDSDDE(5,5) = MU
      DDSDDE(6,6) = MU

Running the UMAT in ABAQUS. Define Material Properties in ABAQUS using *USER MATERIAL, CONSTANTS=2 and specify Young’s modulus and Poisson’s ratio

Example 2: Von Mises elastoplastic material with isotropic hardening

This second example corresponds to an implementation of an elastoplastic UMAT based on the von Mises yield criterion and isotropic hardening.

SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD,
     1       SCD, RPL, DDSDDT, DRPLDE, DRPLDT, STRAN,
     2       DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF,
     3       DPRED, CMNAME, NDI, NSHR, NTENS, NSTATV,
     4       PROPS, NPROPS, COORDS, DROT, PNEWDT,
     5       CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER,
     6       KSPT, KSTEP, KINC)

      INCLUDE 'ABA_PARAM.INC'
      
      DOUBLE PRECISION STRESS(6), DDSDDE(6,6), STATEV(*),
     1       STRAN(6), DSTRAN(6), PROPS(*), TIME(2)

      INTEGER NPROPS, NPOINTS, I, J, NSHR, NDI
      DOUBLE PRECISION E, NU, SIGY0, HARD_PT_STRN, HARD_PT_SIGY
      DOUBLE PRECISION LAMBDA, MU, TRIAL_STRESS(6), DEV_STRESS(6)
      DOUBLE PRECISION PLASTIC_INCREMENT, PLASTIC_STRAIN, NORM_TRIAL
      DOUBLE PRECISION I1, J2, SQJ2, TOL

      PARAMETER (TOL = 1.0D-6)

      ! Extract elastic material properties
      E = PROPS(1)
      NU = PROPS(2)

      ! Lame constants
      LAMBDA = (E * NU) / ((1 + NU) * (1 - 2 * NU))
      MU = E / (2 * (1 + NU))

      ! Number of hardening points
      NPOINTS = (NPROPS - 2) / 2

      ! Initialize plastic strain from state variable
      PLASTIC_STRAIN = STATEV(1)

      ! Elastic predictor step - trial stress
      DO I = 1, 6
         TRIAL_STRESS(I) = STRESS(I) + 2.0 * MU * DSTRAN(I)
      END DO
      TRIAL_STRESS(1) = TRIAL_STRESS(1) + LAMBDA * (DSTRAN(1) + DSTRAN(2) + DSTRAN(3))
      TRIAL_STRESS(2) = TRIAL_STRESS(2) + LAMBDA * (DSTRAN(1) + DSTRAN(2) + DSTRAN(3))
      TRIAL_STRESS(3) = TRIAL_STRESS(3) + LAMBDA * (DSTRAN(1) + DSTRAN(2) + DSTRAN(3))

      ! Compute J2 invariant (deviatoric stress norm)
      I1 = (TRIAL_STRESS(1) + TRIAL_STRESS(2) + TRIAL_STRESS(3)) / 3.0
      DO I = 1, 3
         DEV_STRESS(I) = TRIAL_STRESS(I) - I1
      END DO
      J2 = 0.5 * (DEV_STRESS(1)**2 + DEV_STRESS(2)**2 + DEV_STRESS(3)**2 +
     &              2.0 * (TRIAL_STRESS(4)**2 + TRIAL_STRESS(5)**2 + TRIAL_STRESS(6)**2))
      SQJ2 = SQRT(J2)

      ! Determine yield stress from the hardening curve
      IF (PLASTIC_STRAIN .LE. 0.0) THEN
         SIGY0 = PROPS(3) ! Initial yield stress
      ELSE
         DO I = 1, NPOINTS-1
            HARD_PT_STRN = PROPS(2 + 2*I + 1)
            HARD_PT_SIGY = PROPS(2 + 2*I + 2)
            IF (PLASTIC_STRAIN .LE. HARD_PT_STRN) EXIT
         END DO
         ! Linear interpolation between points
         SIGY0 = PROPS(2 + 2*(I-1) + 2) + 
     &        (PLASTIC_STRAIN - PROPS(2 + 2*(I-1) + 1)) * 
     &        (PROPS(2 + 2*I + 2) - PROPS(2 + 2*(I-1) + 2)) /
     &        (PROPS(2 + 2*I + 1) - PROPS(2 + 2*(I-1) + 1))
      END IF

      ! Check yield condition
      IF (SQJ2 .LE. SIGY0) THEN
         ! Elastic response
         STRESS = TRIAL_STRESS
         DDSDDE = 0.0
         RETURN
      ELSE
         ! Plastic response - return mapping
         PLASTIC_INCREMENT = (SQJ2 - SIGY0) / (2 * MU)
         DO I = 1, 6
            STRESS(I) = TRIAL_STRESS(I) - 2 * MU * (TRIAL_STRESS(I) / SQJ2) * PLASTIC_INCREMENT
         END DO
         STATEV(1) = PLASTIC_STRAIN + PLASTIC_INCREMENT

         ! Tangent modulus for plasticity
         DDSDDE(1,1) = LAMBDA + 2 * MU - (2 * MU**2 / (2 * MU))
         DDSDDE(1,2) = LAMBDA - (2 * MU**2 / (2 * MU))
         DDSDDE(2,1) = DDSDDE(1,2)
         DDSDDE(2,2) = DDSDDE(1,1)
         DDSDDE(3,3) = DDSDDE(1,1)
         DDSDDE(4,4) = MU - (MU**2 / (2 * MU))
         DDSDDE(5,5) = DDSDDE(4,4)
         DDSDDE(6,6) = DDSDDE(4,4)
      END IF

      RETURN
      END

The first two entries of PROPS are E (Young’s modulus) and NU (Poisson’s ratio). The remaining entries contain plastic strain and corresponding yield stress values for the hardening curve.

      ! Extract elastic material properties
      E = PROPS(1)
      NU = PROPS(2)

      ! Lame constants
      LAMBDA = (E * NU) / ((1 + NU) * (1 - 2 * NU))
      MU = E / (2 * (1 + NU))

      ! Number of hardening points
      NPOINTS = (NPROPS - 2) / 2

A trial stress is computed assuming an elastic response.

      ! Initialize plastic strain from state variable
      PLASTIC_STRAIN = STATEV(1)

      ! Elastic predictor step - trial stress
      DO I = 1, 6
         TRIAL_STRESS(I) = STRESS(I) + 2.0 * MU * DSTRAN(I)
      END DO
      TRIAL_STRESS(1) = TRIAL_STRESS(1) + LAMBDA * (DSTRAN(1) + DSTRAN(2) + DSTRAN(3))
      TRIAL_STRESS(2) = TRIAL_STRESS(2) + LAMBDA * (DSTRAN(1) + DSTRAN(2) + DSTRAN(3))
      TRIAL_STRESS(3) = TRIAL_STRESS(3) + LAMBDA * (DSTRAN(1) + DSTRAN(2) + DSTRAN(3))

The code iterates through PROPS to identify the two closest data points on the hardening curve. It then interpolates between these points to get the yield stress for the current plastic strain.

      ! Compute J2 invariant (deviatoric stress norm)
      I1 = (TRIAL_STRESS(1) + TRIAL_STRESS(2) + TRIAL_STRESS(3)) / 3.0
      DO I = 1, 3
         DEV_STRESS(I) = TRIAL_STRESS(I) - I1
      END DO
      J2 = 0.5 * (DEV_STRESS(1)**2 + DEV_STRESS(2)**2 + DEV_STRESS(3)**2 +
     &              2.0 * (TRIAL_STRESS(4)**2 + TRIAL_STRESS(5)**2 + TRIAL_STRESS(6)**2))
      SQJ2 = SQRT(J2)

      ! Determine yield stress from the hardening curve
      IF (PLASTIC_STRAIN .LE. 0.0) THEN
         SIGY0 = PROPS(3) ! Initial yield stress
      ELSE
         DO I = 1, NPOINTS-1
            HARD_PT_STRN = PROPS(2 + 2*I + 1)
            HARD_PT_SIGY = PROPS(2 + 2*I + 2)
            IF (PLASTIC_STRAIN .LE. HARD_PT_STRN) EXIT
         END DO
         ! Linear interpolation between points
         SIGY0 = PROPS(2 + 2*(I-1) + 2) + 
     &        (PLASTIC_STRAIN - PROPS(2 + 2*(I-1) + 1)) * 
     &        (PROPS(2 + 2*I + 2) - PROPS(2 + 2*(I-1) + 2)) /
     &        (PROPS(2 + 2*I + 1) - PROPS(2 + 2*(I-1) + 1))
      END IF

If the trial stress satisfies the yield criterion (von Mises), no plasticity occurs, and the trial stress is assigned directly to STRESS. Otherwise, a return mapping algorithm adjusts STRESS.

      ! Check yield condition
      IF (SQJ2 .LE. SIGY0) THEN
         ! Elastic response
         STRESS = TRIAL_STRESS
         DDSDDE = 0.0
         RETURN
      ELSE

The PLASTIC_INCREMENT and STATEV (holding accumulated plastic strain) are updated accordingly.

      ELSE
         ! Plastic response - return mapping
         PLASTIC_INCREMENT = (SQJ2 - SIGY0) / (2 * MU)
         DO I = 1, 6
            STRESS(I) = TRIAL_STRESS(I) - 2 * MU * (TRIAL_STRESS(I) / SQJ2) * PLASTIC_INCREMENT
         END DO
         STATEV(1) = PLASTIC_STRAIN + PLASTIC_INCREMENT

The tangent modulus matrix DDSDDE is updated to reflect the elastoplastic response, facilitating convergence in the FE solution.

         ! Tangent modulus for plasticity
         DDSDDE(1,1) = LAMBDA + 2 * MU - (2 * MU**2 / (2 * MU))
         DDSDDE(1,2) = LAMBDA - (2 * MU**2 / (2 * MU))
         DDSDDE(2,1) = DDSDDE(1,2)
         DDSDDE(2,2) = DDSDDE(1,1)
         DDSDDE(3,3) = DDSDDE(1,1)
         DDSDDE(4,4) = MU - (MU**2 / (2 * MU))
         DDSDDE(5,5) = DDSDDE(4,4)
         DDSDDE(6,6) = DDSDDE(4,4)

Running the UMAT in ABAQUS.

Define Material Properties, specify *USER MATERIAL, CONSTANTS=<total number of constants> where the first two constants are E and NU, and the remaining constants define the hardening curve as pairs of plastic strain and yield stress values.

Conclusion

This article introduced the basics of creating a UMAT in ABAQUS with two initial examples. Subsequent posts in this series will focus on translating UMATs to CALCULIX and exploring innovative techniques, such as neural network-based UMATs.

Deja un comentario