Hello. I have been troubleshooting a UHYPER subroutine, and I am unsure how to continue after reviewing the manual and searching online forums. I am hoping somebody can get me pointed in the right direction.

To start, I coded UHYPER with the incompressible form of the Yeoh strain-energy function (SEF) to ensure proper linking and execution:

U = C1*(I1 – 3) + C2*(I1 – 3)^2 + C3*(I1 – 3)^3

where C1,C2,C3 are constants and I1 is the first invariant. (Side note: In this thread, I am dropping ‘bars’ from any variables barred in the Abaqus manual since, if I interpret correctly, barred and unbarred terms are the same for incompressible material at constant temperature.) Indeed, my UHYPER model matched Abaqus’ built-in Yeoh model with a unit-cube 3d element (C3D8RH w/ enhanced hourglass control) displaced 10% in uniaxial extension.

My goal is to implement a modified Yeoh SEF (Amin, Alam, & Okui 2002) with non-integer exponents. For simplicity, let C3 = 0 so that:

U = C1*(I1 – 3) + C2*(I1 – 3)^m
dU/dI1 = C1 + m*C2*(I1 – 3)^(m – 1)

If m < 1, Abaqus returns I1 = NaN. In such cases, (I1 – 3) is in the denominator of the second term in dU/dI1 which initializes dU/dI1 = Infinity. Other initialization values are I1 = 3 and U = 0, both correct. After the initialization step, I1 = NaN which corrupts stress calculations. The simulation actually completes, but no results plot. If m >= 1, all computations appear correct and the visualization module works as expected.

The main question: why is Abaqus computing I1 = NaN when m < 1? Is it an implementation issue, or is something fundamentally wrong with my SEF?

Some additional observations & thoughts:

1. I specify a single step and Abaqus conducts 4 equilibrium iterations when solving.

2. From a phenomenological perspective, I don’t see any reason m < 1 is unacceptable. I have no problems calculating reasonable uniaxial, pure shear, and equi-biaxial responses with the SEF in Excel when m < 1. If somebody thinks differently, please advise.

3. For a fixed uniaxial displacement, shouldn’t I1 be independent of m? I ran simulations with m = 1.2 and m = 2, and both gave the same I1. Of course, stresses did change.

4. Section 4.6.1 of the Abaqus-6.14 manual shows I1 = trace(F*F^T) where F is deformation gradient and ‘T’ indicates transpose. I am unable to figure out how to write any values of the deformation gradient (which are stored in DFGRD0) when using UHYPER subroutine, so I can’t troubleshoot this part of the execution. The mathematical details somewhat elude me, but I can’t think of any reason I1 would calculate as NaN from its defining equation.

5. It looks like UMAT offers more de-bugging capability since it has many additional variables in its arguments, but ideally I want to avoid UMAT.

6. Switching to fully-integrated element does not fix issue.

7. Solver: Abaqus/Standard with static step; nlgeom=ON.

8. I’ve read that numerical precision can be an issue. I specify full precision in CAE.

9. Warning 1 from .dat file: USER SUBROUTINE UHYPER WILL BE USED WITH THE STAVEV ARRAY DIMENSIONED TO ZERO SINCE THE *DEPVAR OPTION IS NOT USED WITH THIS MATERIAL. CONSEQUENTLY, DEFINING STATEV ENTRIES IN SUBROUTINE UHYPER WILL CAUSE CODE EXECUTION ERRORS.

10. Warning 2 from .dat file: THE PARAMETER HOURGLASS = ENHANCED ON THE SECTION CONTROLS OPTION IS RELEVANT FOR THESE ELEMENTS: C3D8R, CAX4R, CGAX4R, CPEG4R, CPE4R, CPS4R, 3D4R, S4R, SC8R AND THEIR HYBRID, THERMAL AND PRESSURE COUNTERPARTS WHEREVER APPLICABLE. IT IS ALSO RELEVANT FOR ALL TYPES OF MODIFIED TRIANGULAR AND TETRAHEDRAL ELEMENTS. THIS WARNING CAN BE IGNORED IF THE FEATURE IS APPLIED TO THESE ELEMENT TYPES ONLY.

*****************
SUBROUTINE
*****************

SUBROUTINE UHYPER(BI1,BI2,AJ,U,UI1,UI2,UI3,TEMP,NOEL,CMNAME,
\\\$ INCMPFLAG,NUMSTATEV,STATEV,NUMFIELDV,
\\\$ FIELDV,FIELDVINC,NUMPROPS,PROPS)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION U(2),UI1(3),UI2(6),UI3(6),STATEV(*),FIELDV(*),
\\\$ FIELDVINC(*),PROPS(*)
C
C10 = 5.0
C20 = 2.0
EMM = 0.5
C
U(1) = C10*(BI1-3.0)+C20*(BI1-3.0)**EMM
U(2) = 0.0
C
UI1(1) = C10+EMM*C20*(BI1-3.0)**(EMM-1.0)
UI1(2) = 0.0
UI1(3) = 0.0
C
UI2(1) = EMM*(EMM-1.0)*C20*(BI1-3.0)**(EMM-2.0)
UI2(2) = 0.0
UI2(3) = 0.0
UI2(4) = 0.0
UI2(5) = 0.0
UI2(6) = 0.0
C
UI3(1) = 0.0
UI3(2) = 0.0
UI3(3) = 0.0
UI3(4) = 0.0
UI3(5) = 0.0
UI3(6) = 0.0
C
WRITE(*,205) 'I1=',BI1,'U=',U(1),'dU/dI1=',UI1(1)
205 FORMAT(' ',A3,F8.4,A6,F8.4,A11,F8.4,F8.4)
C
RETURN
END

******************************************************************
OUTPUTS FROM ‘WRITE’ COMMAND IN SUBROUTINE
******************************************************************
I1= 3.0000 U= 0.0000 dU/dI1=Infinity
I1= 3.0000 U= 0.0000 dU/dI1=Infinity
I1= NaN U= NaN dU/dI1= NaN
I1= NaN U= NaN dU/dI1= NaN
I1= NaN U= NaN dU/dI1= NaN
I1= NaN U= NaN dU/dI1= NaN
I1= NaN U= NaN dU/dI1= NaN
I1= NaN U= NaN dU/dI1= NaN
I1= NaN U= NaN dU/dI1= NaN