Hello. I have been troubleshooting a UHYPER subroutine, and I am unsure how to continue after reviewing the manual and searching many 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:

Some additional observations & thoughts:

1. 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.

2. 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.

3. 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.

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

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

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?**Some additional observations & thoughts:

1. 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.

2. 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.

3. 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.

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

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