Hi AsterO'dactyle, thank you for your reply. I tried incorporating your comments into my UMAT file, but I still received the same FPE error. Here's the UMAT file and values inputted in the AsterStudy module for further comments. Thank you.
Gmax = 131323457.5
gamma_ref = 0.01
nu = 0.45
xi = 0.03
===
`MODULE UMAT_Module
IMPLICIT NONE
CONTAINS
SUBROUTINE Compute_GammaEq(STRAN, gamma_eq)
IMPLICIT NONE
REAL(8), DIMENSION(6), INTENT(IN) :: STRAN
REAL(8), INTENT(OUT) :: gamma_eq
REAL(8) :: mean_strain
REAL(8), DIMENSION(6) :: dev
REAL(8) :: g_eq2
! Compute deviatoric part of strain
mean_strain = (STRAN(1) + STRAN(2) + STRAN(3)) / 3.0D0
dev(1:3) = STRAN(1:3) - mean_strain
dev(4:6) = STRAN(4:6)
! Compute squared gamma_eq safely
g_eq2 = 2.0D0 * (dev(1)**2 + dev(2)**2 + dev(3)**2 + 2.0D0*(dev(4)**2 + dev(5)**2 + dev(6)**2))
gamma_eq = SQRT(MAX(1.0D-20, g_eq2)) ! prevent SQRT(negative or 0)
END SUBROUTINE Compute_GammaEq
SUBROUTINE UMAT(STRAN, DSTRAN, STRESS, DDSDDE, STATEV, PROPS, TIME, DTIME) BIND(C, NAME="UMAT_")
IMPLICIT NONE
REAL(8), DIMENSION(6), INTENT(IN) :: STRAN, DSTRAN
REAL(8), DIMENSION(6), INTENT(INOUT) :: STRESS
REAL(8), DIMENSION(6,6), INTENT(OUT) :: DDSDDE
REAL(8), DIMENSION(), INTENT(INOUT) :: STATEV
REAL(8), DIMENSION(), INTENT(IN) :: PROPS
REAL(8), DIMENSION(2), INTENT(IN) :: TIME
REAL(8), INTENT(IN) :: DTIME
! Material parameters
REAL(8) :: Gmax, gamma_ref, nu, xi, eta
REAL(8) :: gamma_eq, G, K, lam, dt, strain_rate, diag_val
REAL(8) :: vol_strain
INTEGER :: i, j
! Clamp and protect material parameters
Gmax = MAX(PROPS(1), 1.0D-6)
gamma_ref = MAX(PROPS(2), 1.0D-6)
nu = MIN(MAX(PROPS(3), 1.0D-4), 0.48D0) ! prevent incompressibility
xi = MAX(PROPS(4), 0.0D0)
dt = MAX(DTIME, 1.0D-9)
eta = 2.0D0 * xi * Gmax
! Elastic bulk modulus (K) and Lamé’s parameter (λ)
K = 2.0D0 * Gmax * (1.0D0 + nu) / MAX(3.0D0 * (1.0D0 - 2.0D0 * nu), 1.0D-6)
lam = K - 2.0D0 * Gmax / 3.0D0
! Compute equivalent shear strain
CALL Compute_GammaEq(STRAN, gamma_eq)
gamma_eq = MAX(gamma_eq, 1.0D-12) ! ensure strictly positive
! G degradation (Hardin-Drnevich type)
G = Gmax / (1.0D0 + gamma_eq / gamma_ref)
G = MAX(MIN(G, Gmax), 1.0D-6) ! clamp G to safe bounds
! Store in state variables
STATEV(1) = gamma_eq
STATEV(2) = G
! Compute volumetric strain
vol_strain = STRAN(1) + STRAN(2) + STRAN(3)
! Stress = elastic + damping
STRESS(1:3) = lam * vol_strain + 2.0D0 * G * STRAN(1:3)
STRESS(4:6) = 2.0D0 * G * STRAN(4:6)
DO i = 1, 6
strain_rate = DSTRAN(i) / dt
strain_rate = MAX(MIN(strain_rate, 1.0D6), -1.0D6)
STRESS(i) = STRESS(i) + 2.0D0 * eta * strain_rate
END DO
! Tangent stiffness
DDSDDE(:,:) = 0.0D0
DO i = 1, 6
diag_val = 2.0D0 * G + 2.0D0 * eta / dt
DDSDDE(i,i) = MAX(MIN(diag_val, 1.0D9), 1.0D-6)
END DO
! Add Lamé contributions to normal-normal entries
DO i = 1, 3
DDSDDE(i,i) = DDSDDE(i,i) + lam
DO j = 1, 3
DDSDDE(i,j) = DDSDDE(i,j) + lam
END DO
END DO
END SUBROUTINE UMAT
END MODULE UMAT_Module
`
