Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

55, 3, pp. 787-800, Warsaw 2017
DOI: 10.15632/jtam-pl.55.3.787

FINITE ELEMENT IMPLEMENTATION OF SLIGHTLY COMPRESSIBLE

AND INCOMPRESSIBLE FIRST INVARIANT-BASED HYPERELASTICITY:

THEORY, CODING, EXEMPLARY PROBLEMS

Cyprian Suchocki

Warsaw University of Technology, Department of Mechanics and Armament Technology, Warsaw, Poland

e-mail: c.suchocki@imik.wip.pw.edu.pl

The present study is concernedwith the finite element (FE) implementation of slightly com-
pressible hyperelasticmaterialmodels. A class of constitutive equations is consideredwhere
the isochoric potential functions are based on the first invariant of the right Cauchy-Green
(C-G) deformation tensor. Special attention is paid to the most recently developed model
formulations.The incremental formof hyperelasticity and its numerical implementation into
both commercial and non-commercial FE software are discussed. A Fortran 77UMAT code
is attached which allows for a simple implementation of arbitrary first invariant-based con-
stitutive models into Abaqus and Salome-Meca FE packages. Several exemplary problems
are considered.

Keywords: hyperelasticity, finite element method, UMAT, elasticity tensor

1. Introduction

The hyperelastic constitutive equations are nowadays available in every advanced FE program.
However, thematerial libraries of FE software usually include only a number of standard hyper-
elastic models such as: neo-Hooke, Mooney-Rivlin, Ogden or Yeoh. Less celebrated or newly
developed constitutive models can be implemented into a FE program by taking advantage of
a proper user subroutine. The FE package Abaqus provides three user subroutines which allow
one to define a custom hyperelastic model, i.e. UHYPER (for isotropic hyperelastic materials),
UANISOHYPER (for anisotropic hyperelastic materials) and UMAT (a general purpose subro-
utine which can be utilized for implementing any kind of constitutive equation), cf. Hibbit et
al. (2008). Due to the method of FE implementation used for slightly compressible hyperelasti-
city in Abaqus, it is not recommended to utilize the subroutine UHYPER for all kinds of finite
elements (cf. Jemioło, 2002). Thus, in the case of slightly compressible hyperelastic materials,
i.e. thematerials with decoupled volumetric and isochoric responses, the subroutineUMATmi-
ght be preferred. Both UHYPER and UANISOHYPER subroutines can be utilized to define
nonlinear viscoelastic models based on the viscoelasticity theory used byAbaqus. Alternatively,
a proper option allows one to simulate the Mullins effect in a hyperelastic material defined by
the aforementioned subroutines1. On the other hand, the subroutine UMAT is a much more
powerful tool which enables one to define an arbitrary constitutive theory, including those based
on hyperelasticity such as nonlinear viscoelasticity (e.g. Suchocki 2013) or growth models (e.g.
Young et al., 2010), so that the user is not limited by the built-in options of Abaqus.
The subroutine UMAT is a Fortran 77 code which is called during every iteration of the

Newton-Raphson numerical procedure to calculate components of the stress tensor and thema-
terial Jacobian which is also reffered to as tangent modulus or (in the case of elastic materials)

1Thenonlinear viscoelasticity and theMullins effectmust be used separately asAbaqusdoes not allow
for combining these behaviors.



788 C. Suchocki

the elasticity tensor, cf. Hibbit et al. (2008). The material Jacobian may be defined either in
an approximate or (if possible) an analytical form, which is usually very difficult to determine.
The approximate material Jacobians always worsen the rate of convergence of the numerical
calculations. It was demonstrated by Stein and Sagar (2008) that for the neo-Hooke hypere-
lastic model, the quadratic rate of convergence2 is obtained only when the analytical material
Jacobian is used. The utilization of the approximate material Jacobians resulted in worsening
the convergence rate and, in the case of some of the considered problems and finite element ty-
pes, it caused lack of convergence. Thus, it is always recommended to use an analytical material
Jacobian whenever it is available.
In this study, theFE implementation of slightly compressible isotropic hyperelastic constitu-

tive models that are not included in any of the commercial and non-commercial CAE packages
is discussed. The stored energy functions that are based on the first invariant of the isochoric
right C-G tensor are considered. The focus is on the recently developed models for polymeric
materials (Gent, 1996; Jemioło, 2002; Lopez-Pamies, 2010, da Silva Soares, 2008; Khajehsaeid
et al., 2013) and on some model formulations used in soft tissue biomechanics (Demiray, 1972;
Demiray et al., 1988). The general framework for deriving an analytical material Jacobian is
presented. A subroutineUMAT is attached allowing for using the newly developed exponential-
logarithmic model (Khajehsaeid et al., 2013) in both Abaqus and Salome-Meca FE packages.
The code structure is universal so that any other first invariant-based slightly compressible or in-
compressible hyperelastic model can be easily implemented by simply changing the expressions
for the stored energy derivatives. A number of exemplary problems were solved for selected
energy potentials. The presented UMAT code can be upgraded to define nonlinear viscoelastic,
elastoplastic, viscoplastic or other behavior using arbitrary constitutive theory.

2. Slightly compressible hyperelastic materials

In the following derivations, the multiplicative split of the deformation gradient tensor into the
volumetric and isochoric component is utilized (e.g. Jemioło, 2016), i.e.

F=FvolF Fvol = J
1
31 F= J−

1
3F C=F

T
F= J−

2
3C (2.1)

where J = detF and C is the isochoric right C-G tensor with the following set of algebraic
invariants

Ī1 = trC Ī2 =
1

2

(

(trC)2− trC
2
)

Ī3 =detC=1 (2.2)

In the case of slightly compressible hyperelastic materials, the stored energy function is consi-
dered to be the sum of the volumetric contribution U and the isochoric partW, thus

W(C)=U(J)+W(Ī1, Ī2) S=2
∂W

∂C

∣
∣
∣
∣
C=CT

(2.3)

where the most general form of the constitutive equation is given by Eq. (2.3)2
3. After substi-

tuting Eq. (2.3)1 into Eq. (2.3)2, the decoupled form of the constitutive equation is found

S= JpC−1+J−
2
3 DEV

[

S
]

p=
∂U

∂J
S=2

∂W

∂C

∣
∣
∣
∣
C=C

T
(2.4)

with DEV[•] = [•]− 1
3

(

[•] ·C
)

C
−1
being a deviator in the reference configuration.

2The quadratic convergencemeans that the error at the current iteration is proportional to the square
of the error from the previous iteration.
3The adopted notation emphasizes the fact that symmetrization is carried out after calculating a

derivative.



Finite element implementation of slightly compressible and incompressible... 789

3. Material Jacobian tensor

Taking a directional derivative of Eq. (2.4)1 with respect to C, an incremental constitutive
relation is found, see e.g. Jemioło and Gajewski (2014)

∆S=C ·
1

2
∆C C=2

∂S

∂C

∣
∣
∣
∣
C=CT

=4
∂2W

∂C⊗∂C

∣
∣
∣
∣
C=CT

C=Cvol+Ciso (3.1)

AssumingU =U(J) andW =W(Ī1), the expressions for the volumetric and the isochoric parts
of the elasticity tensor can be derived

C
vol = J

∂U

∂J

(

C
−1⊗C−1−2IC−1

)

+J2
∂2U

∂J2
C
−1⊗C−1

C
iso =−

4

3
J−

2
3
∂W

∂Ī1

[

1⊗C−1+C−1⊗1− I1
(

IC−1 +
1

3
C
−1⊗C−1

)]

+J−
4
3C
W

C
W
=4
∂2W

∂Ī21

[

1⊗1−
1

3
I1(1⊗C

−1+C−1⊗1)+
1

9
I21C

−1⊗C−1
]

(3.2)

where

IC−1 =
1

2

[

(C−1⊗C−1)
(2,3)
T +(C−1⊗C−1)

(2,4)
T

]

=
1

2
(C−1
IK
C−1
JL
+C−1
IL
C−1
JK
)EI⊗EJ⊗EK⊗EL

is the fourth order identity tensor in the reference configuration with the Cartesian base {EK}
(K =1,2,3)4, see e.g. Suchocki (2011).
The incremental constitutive law given byEq. (3.1)1 can be transformed into a form relating

the incremental Oldroyd (convected) rate of the Kirchhoff stress to the increment of the strain
rate tensor, i.e.

Lvτ =∆τ −∆Lτ −τ∆L
T =Cτc ·∆D (3.3)

where∆L=∆FF−1 is the increment of the velocity gradient,whereasCτc is thepushed-forward
form of the material Jacobian

C
τc =FiPFjQFkRFlSCPQRSei⊗ej ⊗ek⊗el (3.4)

with {ek} (k = 1,2,3) being the Cartesian base in the current configuration. The elasticity
tensor takes the following form

C
τc =
4

3

∂W

∂Ī1

[

Ī1
(

I−
1

3
1⊗1

)

−
(
1⊗ dev(B)+ dev(B)⊗1

)
]

+4
∂2W

∂Ī21
dev(B)⊗ dev(B)

+J
[(∂U

∂J
+J
∂2U

∂J2

)

1⊗1−2
∂U

∂J
I
]

(3.5)

where

I=1✸1=
1

2

[

(1⊗1)
(2,3)
T +(1⊗1)

(2,4)
T

]

=
1

2
(δikδjl+ δilδjk)ei⊗ej ⊗ek⊗el

and dev[•] = [•]− 1
3
([•] ·1)1.

4The followingnotation is used: [•]
(µ,ν)
T =

(
[•]ijklei⊗ ej

︸︷︷︸

µ

⊗ek⊗ el
︸︷︷︸

ν

)(µ,ν)
T = [•]ijklei⊗ el

︸︷︷︸

µ

⊗ek⊗ ej
︸︷︷︸

ν

.



790 C. Suchocki

The FE software Abaqus utilizes the incremental constitutive equation written in terms of
the incremental Zaremba-Jaumann rate of the Kirchhoff stress (cf. Hibbit et al. 2008), i.e.

τ
∇ =∆τ −∆Wτ −τ∆WT = JCMJ ·∆D (3.6)

where, respectively

C
MJ =

1

J

(
C
τc+1✸τ +τ✸1

)
τ = Jp1+2

∂W

∂Ī1
dev(B) (3.7)

and

1✸τ =
1

2

[

(1⊗τ)
(2,3)
T +(1⊗τ)

(2,4)
T

]

τ✸1=
1

2

[

(τ ⊗1)
(2,3)
T +(τ ⊗1)

(2,4)
T

]

and

∆W=
1

2
(∆L−∆LT) ∆D=

1

2
(∆L+∆LT) (3.8)

The fourth order tensor CMJ is the material Jacobian which should be coded in the subroutine
UMAT. For the considered class of hyperelastic materials, it takes the form

C
MJ =

2

J

∂W

∂Ī1

[

1✸dev(B)+ dev(B)✸1+
2

3
Ī1
(

I−
1

3
1⊗1

)

−
2

3

(
1⊗ dev(B)+ dev(B)⊗1

)
]

+
4

J

∂2W

∂Ī21
dev(B)⊗ dev(B)+

(∂U

∂J
+J
∂2U

∂J2

)

1⊗1

(3.9)

4. Finite element implementation

4.1. General

In Fig. 1, the interaction of the subroutineUMATwith theAbaqus package is illustrated for
the Newton-Raphson iterative procedure during a single time increment (cf. Hibbit et al. 2008).

Fig. 1. Flow chart for the interaction of Abaqus and UMAT

The subroutine UMAT calculates the components of Cauchy stress and material Jacobian
for each Gauss integration point. These quantities are subsequently used byAbaqus to form up
the element stiffness matrix. Finally, the global stiffness matrix is assembled by Abaqus using
the element stiffnessmatrices. The user subroutines used in other FE packages to define custom
constitutive equations are integrated with the remainder of the program in a similar way and
play the same role.



Finite element implementation of slightly compressible and incompressible... 791

4.2. Variables and dimensions

In the following table, the meaning of the variables used in the Fortran 77 code has been
explained. The dimensions of array variables have been specified in proper indices. The lengthy
definitions of the auxiliary variables have been skipped.

Number of direct stress components NDI
Number of shear stress components NSHR
Array of material constants PROPS(I)
Deformation gradient tensor F3×3 DFGRD1(I,J)
Jacobian determinant DET
Isochoric deformation gradient matrixF3×3 DISTGR(I,J)
Isochoric Left C-G deformation tensor matrixB6×1 BBAR(I)
Trace of B divided by 3 TRBBAR
First partial derivative ∂JU DUDJ
Second partial derivative ∂2

J2
U DDUDDJ

First partial derivative ∂Ī1W DWDI1
Second partial derivative ∂2

Ī2
1
W DDWDDI1

Cauchy stress tensor matrix σ6×1 STRESS(I)

Material Jacobian matrix CMJ6×6 DDSDDE(I,J)
Auxiliary variables EK, PR, SCALE,

TERM1, TERM2,
TERM3

According to the rule adopted in Abaqus, the column matrix components 1,2, . . . ,6 corre-
spond to the scalar components of the second order tensor: 11, 22, 33, 12, 13, 23, respectively.

4.3. User subroutine UMAT

Algorithm for the implementation in ABAQUS

Input data:F3×3 (DFGRD1), NDI, NSHR

1. Calculate Jacobian determinant J (DET)

J =detF3×3

2. Calculate isochoric deformation gradientF3×3 (DISTGR)

F3×3 = J
−
1
3F3×3

3. Calculate left C-G deformation tensor B6×1 (BBAR)

B3×3 =F3×3F
T
3×3 B6×1 =

{
B11 B22 B33 B12 B13 B23

}T

4. Calculate Cauchy stress matrixσ6×1 (STRESS)

5. Calculate Material Jacobian matrix CMJ6×6 (DDSDDE).



792 C. Suchocki

4.4. Coding in Fortran 77

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,

2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,

3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,

4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)

!

INCLUDE ’ABA PARAM.INC’

!

CHARACTER*8 MATERL

DIMENSION STRESS(NTENS),STATEV(NSTATV),

1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

2 STRAN(NTENS),DSTRAN(NTENS),DFGRD0(3,3),DFGRD1(3,3),

3 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3)

!

! LOCAL ARRAYS

! ----------------------------------------------------------------

! BBAR - DEVIATORIC RIGHT CAUCHY-GREEN TENSOR

! DISTGR - DEVIATORIC DEFORMATION GRADIENT (DISTORTION TENSOR)

! ----------------------------------------------------------------

!

REAL*8 BBAR,DISTGR

DIMENSION BBAR(6),DISTGR(3,3)

!

PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0)

!

! ----------------------------------------------------------------

! UMAT FOR COMPRESSIBLE EXPONENTIAL-LOGARITHMIC HYPERELASTICITY

!

! WARSAW UNIVERSITY OF TECHNOLOGY

! CYPRIAN SUCHOCKI, JULY 2015

!

! CANNOT BE USED FOR PLANE STRESS

! ----------------------------------------------------------------

! PROPS(1) - A

! PROPS(2) - A1

! PROPS(3) - B

! PROPS(4) - D1

! ----------------------------------------------------------------

REAL*8 A, A1, B, D1, TERM1, TERM2, TERM3, DUDJ, DDUDDJ,

1 DWDI1, DDWDDI1, TRBBAR, DET, SCALE

!

! ELASTIC PROPERTIES

!

A=0.195

A1=0.018 ! originally a

B=0.22

D1=0.000000033

!

! JACOBIAN AND DISTORTION TENSOR

!

DET=DFGRD1(1, 1)*DFGRD1(2, 2)*DFGRD1(3, 3)

1 -DFGRD1(1, 2)*DFGRD1(2, 1)*DFGRD1(3, 3)

IF(NSHR.EQ.3) THEN

DET=DET+DFGRD1(1, 2)*DFGRD1(2, 3)*DFGRD1(3, 1)

1 +DFGRD1(1, 3)*DFGRD1(3, 2)*DFGRD1(2, 1)

2 -DFGRD1(1, 3)*DFGRD1(3,1)*DFGRD1(2, 2)

3 -DFGRD1(2, 3)*DFGRD1(3, 2)*DFGRD1(1, 1)

END IF

SCALE=DET**(-ONE/THREE)

DO K1=1, 3

DO K2=1, 3

DISTGR(K2, K1)=SCALE*DFGRD1(K2, K1)

END DO

END DO

!

! CALCULATE LEFT CAUCHY-GREEN TENSOR

!

BBAR(1)=DISTGR(1, 1)**2+DISTGR(1, 2)**2+DISTGR(1, 3)**2



Finite element implementation of slightly compressible and incompressible... 793

BBAR(2)=DISTGR(2, 1)**2+DISTGR(2, 2)**2+DISTGR(2, 3)**2

BBAR(3)=DISTGR(3, 3)**2+DISTGR(3, 1)**2+DISTGR(3, 2)**2

BBAR(4)=DISTGR(1, 1)*DISTGR(2, 1)+DISTGR(1, 2)*DISTGR(2, 2)

1 +DISTGR(1, 3)*DISTGR(2, 3)

IF(NSHR.EQ.3) THEN

BBAR(5)=DISTGR(1, 1)*DISTGR(3, 1)+DISTGR(1, 2)*DISTGR(3, 2)

1 +DISTGR(1, 3)*DISTGR(3, 3)

BBAR(6)=DISTGR(2, 1)*DISTGR(3, 1)+DISTGR(2, 2)*DISTGR(3, 2)

1 +DISTGR(2, 3)*DISTGR(3, 3)

END IF

!

! CALCULATE THE STRESS

!

TRBBAR=(BBAR(1)+BBAR(2)+BBAR(3))/THREE

DUDJ=2/D1*(DET-ONE)

DDUDDJ=2/D1

DWDI1=A*(exp(A1*(THREE*TRBBAR-THREE))

1 -B*log(THREE*TRBBAR-TWO))

DDWDDI1=A*(A1*exp(A1*(THREE*TRBBAR-THREE))

1 -B*(THREE*TRBBAR-TWO)**(-ONE))

TERM1=-FOUR/(THREE*DET)*DWDI1

TERM2=FOUR/DET*DDWDDI1

TERM3=DET*DDUDDJ

CALL CALCSTRESS(BBAR,TRBBAR,DET,DUDJ,DWDI1,NDI,NSHR,

1 STRESS)

!

! CALCULATE THE STIFFNESS

!

CALL CALCTANGENT(DDSDDE,STRESS,BBAR,TRBBAR,DUDJ,

1 DWDI1,DDWDDI1,TERM1,TERM2,TERM3,DET,NTENS,NSHR)

!

RETURN

END

! ----------------------------------------------------------------

SUBROUTINE CALCSTRESS(BBAR,TRBBAR,DET,DUDJ,DWDI1,NDI,NSHR,

1 STRESS)

REAL*8 BBAR,TRBBAR,DET,DUDJ,DWDI1,STRESS

DIMENSION BBAR(6),STRESS(6)

PARAMETER(TWO=2.D0)

INTEGER NDI,NSHR,K1

DO K1=1,NDI

STRESS(K1)=TWO/DET*DWDI1*( BBAR(K1)-TRBBAR)+DUDJ

END DO

DO K1=NDI+1,NDI+NSHR

STRESS(K1)=TWO/DET*DWDI1*BBAR(K1)

END DO

RETURN

END

! ----------------------------------------------------------------

SUBROUTINE CALCTANGENT(DDSDDE,STRESS,BBAR,TRBBAR,DUDJ,

1 DWDI1,DDWDDI1,TERM1,TERM2,TERM3,DET,NTENS,NSHR)

REAL*8 DDSDDE,STRESS,BBAR,TRBBAR,DUDJ,DWDI1,DDWDDI1,

1 TERM1,TERM2,TERM3,DET

DIMENSION DDSDDE(6,6),STRESS(6),BBAR(6)

INTEGER NTENS,NSHR,K1,K2

PARAMETER(TWO=2.D0, THREE=3.D0, FOUR=4.D0)

DDSDDE(1, 1)=-DUDJ+TERM3+TWO*TERM1*(BBAR(1)-TWO*TRBBAR)+

1 TERM2*(BBAR(1)**TWO+TRBBAR*(-TWO*BBAR(1)+TRBBAR))+



794 C. Suchocki

2 TWO*STRESS(1)

DDSDDE(2, 2)=-DUDJ+TERM3+TWO*TERM1*(BBAR(2)-TWO*TRBBAR)+

1 TERM2*(BBAR(2)**TWO+TRBBAR*(-TWO*BBAR(2)+TRBBAR))+

2 TWO*STRESS(2)

DDSDDE(3, 3)=-DUDJ+TERM3+TWO*TERM1*(BBAR(3)-TWO*TRBBAR)+

1 TERM2*(BBAR(3)**TWO+TRBBAR*(-TWO*BBAR(3)+TRBBAR))+

2 TWO*STRESS(3)

DDSDDE(1, 2)=DUDJ+TERM3+TERM1*(BBAR(1)+BBAR(2)-TRBBAR)+

1 TERM2*(BBAR(1)*BBAR(2)-TRBBAR*(BBAR(1)+BBAR(2))+

2 TRBBAR**TWO)

DDSDDE(1, 3)=DUDJ+TERM3+TERM1*(BBAR(1)+BBAR(3)-TRBBAR)+

1 TERM2*(BBAR(1)*BBAR(3)-TRBBAR*(BBAR(1)+BBAR(3))+

2 TRBBAR**TWO)

DDSDDE(2, 3)=DUDJ+TERM3+TERM1*(BBAR(2)+BBAR(3)-TRBBAR)+

1 TERM2*(BBAR(2)*BBAR(3)-TRBBAR*(BBAR(2)+BBAR(3))

2 +TRBBAR**TWO)

DDSDDE(1, 4)=FOUR/DET*BBAR(4)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(1)-TRBBAR))+STRESS(4)

DDSDDE(2, 4)=FOUR/DET*BBAR(4)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(2)-TRBBAR))+STRESS(4)

DDSDDE(3, 4)=FOUR/DET*BBAR(4)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(3)-TRBBAR))

DDSDDE(4, 4)=-DUDJ+TWO/DET*(TRBBAR*DWDI1+

1 TWO*DDWDDI1*BBAR(4)**TWO)+(STRESS(1)+STRESS(2))/TWO

IF(NSHR.EQ.3) THEN

DDSDDE(1, 5)=FOUR/DET*BBAR(5)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(1)-TRBBAR))+STRESS(5)

DDSDDE(2, 5)=FOUR/DET*BBAR(5)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(2)-TRBBAR))

DDSDDE(3, 5)=FOUR/DET*BBAR(5)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(3)-TRBBAR))+STRESS(5)

DDSDDE(1, 6)=FOUR/DET*BBAR(6)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(1)-TRBBAR))

DDSDDE(2, 6)=FOUR/DET*BBAR(6)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(2)-TRBBAR))+STRESS(6)

DDSDDE(3, 6)=FOUR/DET*BBAR(6)*(-DWDI1/THREE+

1 DDWDDI1*(BBAR(3)-TRBBAR))+STRESS(6)

DDSDDE(5, 5)=-DUDJ+TWO/DET*(TRBBAR*DWDI1+

1 TWO*DDWDDI1*BBAR(5)**TWO)+(STRESS(1)+STRESS(3))/TWO

DDSDDE(6, 6)=-DUDJ+TWO/DET*(TRBBAR*DWDI1+

1 TWO*DDWDDI1*BBAR(6)**TWO)+(STRESS(2)+STRESS(3))/TWO

DDSDDE(4,5)=TERM2*BBAR(4)*BBAR(5)+STRESS(6)/TWO

DDSDDE(4,6)=TERM2*BBAR(4)*BBAR(6)+STRESS(5)/TWO

DDSDDE(5,6)=TERM2*BBAR(5)*BBAR(6)+STRESS(4)/TWO

END IF

DO K1=1, NTENS

DO K2=1, K1-1

DDSDDE(K1, K2)=DDSDDE(K2, K1)

END DO

END DO

RETURN

END

5. Exemplary problems

Anumber of exemplaryFE simulations have been prepared in order to verify the performance of
the developed UMAT code. Seven different types of the isochoric stored energy potentialW(Ī1)
and two types of the volumetric function U(J) have been tested (see Tables 1 and 2). Two
different approaches were used in order to simulate the material near incompressibility, i.e. the
penalty method and the hybrid formulation (e.g. Liu et al. 1994). The results obtained for the
material near incompressibility in the case of homogenous deformations were compared to the
analytical solutions available in the fully incompressible case.



Finite element implementation of slightly compressible and incompressible... 795

Table 1.Material parameter values

Material Constitutive parameters Units

Jemioło (2002) – µ1 =2.228 [MPa]
Lopez-Pamies (2010) µ2 =1.919 [MPa]

α1 =0.6 [-],α2 =−68.73 [–]

Gent (1996) µ=0.27 [MPa]
Jm =85.91 [–]

Khajehsaeid et al. (2013) A=0.195 [MPa]
a=0.018 [–]
b=0.22 [–]

Demiray (1971) c=0.2 [MPa]
β=16 [–]

Demiray et al. (1988) α=10.74E-10 [MPa]
β=7.548E-9 [MPa]
c=1.17 [–]

Da Silva Soares (2008) µ1 =17.999 [MPa]
µ2 =0.17047 [MPa]
a=477.28 [–]

Knowles (1977) µ=264.069 [MPa]
b=54.19 [–]
n=0.2554 [–]

5.1. Simple tension

In the case of uniaxial tension of an incompressible rectangular block (Fig. 2) along the
X1-direction, the deformation is defined by the set of equations

x1 =λ1X1 x2 =λ
−
1
2
1 X2 x3 =λ

−
1
2
1 X3 (5.1)

where the stretch ratio λ1 > 1 and J =1 is assumed. It follows that

I1 =λ
2
1+
2

λ1
W =W(I1) (5.2)

which yields an equation for the axial component of the Lagrange stress

T11 =2
∂W

∂I1

(

λ1−
1

λ21

)

(5.3)

The analytical Eq. (5.3) was used to verify the results of FE calculations. In numerical
simulation, a 15mm×15mm×15mm block was undergoing a uniaxial tension (Fig. 2). In the
first approach, a single C3D85 element was used with thematerial near incompressibility being
enforcedbyusing thepenaltymethod.Thepenalty parameterD1 =33E-9MPa

−1. In the second
approach, a hybrid element C3D8H was utilized. The comparison of the numerical results and
the analytical solution for the incompressiblematerial can be seen inFig. 3. TheFE simulations
were later repeated for the block meshed with 125 elements which produced exactly the same
results.

5Cubic, three-dimensional, 8 nodes.



7
9
6

C
.
S
u
ch
o
ck
i

Table 2.Exemplary isochoric and volumetric stored-energy functions and their derivatives

Material Energy potentialW(Ī1) 1st derivative ∂W/∂Ī1 2nd derivative ∂
2W/∂Ī21

Jemioło (2002) –
∑M
r=1

31−αr
2αr
µr
(
Īαr1 −3

αr
) ∑M

r=1
31−αr
2
µrĪ
αr−1
1

∑M
r=1

31−αr
2
µr(αr −1)Ī

αr−2
1

Lopez-Pamies (2010)

Gent (1996) −
µJm
2
ln
(

1− Ī1−3
Jm

)
µ
2

(

1− Ī1−3
Jm

)−1
µ
2Jm

(

1− Ī1−3
Jm

)−2

Khajehsaeid et al. (2013)
A
[
1
a
ea(Ī1−3)− 1

a
− b

A
[
ea(Ī1−3)− b ln(Ī1−2)

]
A
[
aea(Ī1−3)− b(Ī1−2)

−1
]

+b(Ī1−2)
(
1− ln(Ī1−2)

)
]

Demiray (1971) c
(

eβ(Ī1−3)−1
)

cβeβ(Ī1−3) cβ2eβ(Ī1−3)

Demiray et al. (1988) α
4
(Ī1−3)

2+
β
4c

[

ec(Ī1−3)
2
−1
]

1
2
(Ī1−3)

[

α+βec(Ī1−3)
2
]

1
2

{

α+βec(Ī1−3)
2
[1+2c(Ī1−3)

2]
}

Da Silva Soares (2008)
µ1e
−(Ī1−3)(Ī1−3) µ1e

−(Ī1−3)(4− Ī1) −µ1e
−(Ī1−3)(5− Ī1)

+µ2 ln[1+a(Ī1−3)] +µ2a[1+a(Ī1−3)]
−1 −µ2a

2[1+a(Ī1−3)]
−2

Knowles (1977)
µ
2b

{[

1+ b
n
(Ī1−3)

]n
−1
}

µ
2

[

1+ b
n
(Ī1−3)

]n−1 µb(n−1)
2n

[

1+ b
n
(Ī1−3)

]n−2

Material Energy potential U(J) 1st derivative ∂U/∂J 2nd derivative ∂2U/∂J2

Sussman and Bathe (1987) 1
D1
(J−1)2 2

D1
(J−1) 2

D1

Simo and Taylor (1982) 1
D1
[(J−1)2+(lnJ)2] 2

D1

(

J+ lnJ
J
−1
)

2
D1

[

1+ 1
J2
(1− lnJ)

]



Finite element implementation of slightly compressible and incompressible... 797

Fig. 2. Uniaxial defomation of a single element: (a) distribution of the displacement,
(b) boundary conditions

Fig. 3. Uniaxial tension for various hyperelastic models; comparison of analytical and FE results:
(a) Demiray (1972), (b) Demiray et al. (1988), (c) Exp-Ln, (d) Gent, (e) Jemioło–Lopez-Pamies,

(f) Da Silva Soares

5.2. Simple shear

In the case of simple shear of an incompressible rectangular block in the X1 −X2 plane
(Fig. 4), the deformation is defined by the set of equations

x1 =X1+γX2 x2 =X2 x3 =X3 (5.4)

where γ > 0. The first invariant of the right C-G tensor is given as

I1 = γ
2+3 (5.5)

which yields the following components of the Lagrange stress tensor

T3×3 =
2

3

∂W

∂I1






−γ2 3γ 0
γ(γ2+3) −γ2 0
0 0 −γ2




 (5.6)



798 C. Suchocki

Fig. 4. Shear deformation of a single element: (a) distribution of the displacement,
(b) boundary conditions

The analytical formula for T12 given by Eq. (5.6) was utilized to verify the results of FE
calculations. In numerical simulation, a 15mm×15mm×15mm block was undergoing a simple
shear (Fig. 4). Again, the analysis was carried out using the penaltymethodwith a single C3D8
element (D1 =33E-9MPa

−1) andwas subsequently repeated for a hybrid element C3D8H. The
comparison of the numerical results and the analytical solution for the incompressible material
can be seen in Fig. 5. The FE simulations were later performed for the block meshed with 125
elements with exactly the same results.

Fig. 5. Simple shear for various hyperelastic models; comparison of analytical and FE results:
(a) Demiray (1972), (b) Demiray et al. (1988), (c) Exp-Ln, (d) Gent, (e) Jemioło–Lopez-Pamies,

(f) Da Silva Soares

6. Conclusions

In this paper, the FE implementation of slightly compressible, first invariant-based, isotropic
hyperelastic constitutive equations is discussed. Special attention is paid to the newly developed
models for polymers and some of the stored energy functions used in the soft tissue biome-
chanics. A user subroutine UMAT code is attached, which enables the implementation of the
aformentionedmodels intoAbaqus andSalome-MecaFEpackages. Theperformanceof this code
has been verified using some exemplary problems and an excellent agreement was found with



Finite element implementation of slightly compressible and incompressible... 799

the analytical solutions. It should be emphasized that the stress-stretch (or stress-amount of
shear) relation which yields from the potential function developed by Demiray et al. (1988) is
characterized by a very flat slope in the small strain domain (cf. Figs. 3b and 5b). Thus, for this
particularmodel, a considerably small strain increment should be used initially in order to avoid
convergence problems.ThepresentedUMATcode can be furthermodified in order to define any
constitutive theory that would be an extension of the slightly compressible, first invariant-based
hyperelasticity.

References

1. Da Silva Soares J.F., 2008, Constitutive modeling for biodegradable polymers for application
in endovascular stents, PhD thesis, Texas A&MUniversity

2. Demiray H., 1972, A note on the elasticity of soft biological tissues, Journal of Biomechanics, 5,
3, 309-311

3. Demiray H., Weizsäcker H.W., Pascale K., Erbay H.A., 1988, A stress-strain relation for
a rat abdominal aorta, Journal of Biomechanics, 21, 5, 369-374

4. Gent A.N., 1996, A new constitutive relation for rubber,Rubber Chemistry and Technology, 69,
59-61

5. Hibbit B., KarlssonB., SorensenP., 2008,ABAQUSTheoryManual, Hibbit, Karlsson&So-
rensen Inc.

6. Jemioło S., 2002,A study on the hyperelastic properties of isotropicmaterials (inPolish),Scienti-
fic Surveys ofWarsawUniversity of Technology, 140,WarsawUniversity ofTechnologyPublishing
House,Warsaw

7. Jemioło S., 2016, Constitutive relations of hyperelasticity (in Polish), Studies in the Field of
Engineering,94,TheCommittee onCivilEngineeringandHydroengineeringof thePolishAcademy
of Sciences,Warsaw

8. JemiołoS.,GajewskiM., 2014,Hyperelastoplasticity (inPolish),Monographs of theDepartment
of Strength of Materials, Theory of Elasticity and Plasticity, 4, Warsaw University of Technology
Publishing House,Warsaw

9. Khajehsaeid H., Arghavani J., Naghdabadi R., 2013, A hyperelastic constitutive model for
rubber-like materials,European Journal of Mechanics A/Solids, 38, 2, 144-151

10. Knowles J.K., 1977, The finite anti-plane shear field near the tip of a crack for a class of incom-
pressible elastic solids, International Journal of Fracture, 13, 5, 611-639

11. LiuC.H.,HofstetterG.,MangH.A., 1994, 3Dfinite element analysis of rubber-likematerials
at finite strains,Engineering Computations, 11, 111-128

12. Lopez-Pamies O., 2010,A new I1-based hyperelasticmodel for rubber elasticmaterials,Comptes
Rendus Mecanique, 338, 1, 3-11

13. Simo J.C., Taylor R.L., 1982, Penalty function formulations for incompressible nonlinear ela-
stostatics,Computer Methods in Applied Mechanics and Engineering, 35, 107-118

14. Stein E., SagarG., 2008, Convergence behavior of 3D finite elements for neo-Hookeanmaterial,
Engineering Computations: International Journal for Computer-Aided Engineering and Software,
25, 3, 220-232

15. Suchocki C., 2011, A finite element implementation of Knowles stored-energy function: theory,
coding and applications,Archive of Mechanical Engineering, 58, 3, 319-346

16. Suchocki C., 2013, A quasi-linear viscoelastic rheological model for thermoplastics and resins,
Journal of Theoretical and Applied Mechanics, 51, 1, 117-129



800 C. Suchocki

17. SussmanT., Bathe K.J., 1987,A finite element formulation for nonlinear incompressible hyper-
elastic and inelastic analysis,Computers and Structures, 26, 357-409

18. Young J.M., Yao J., Ramasubramanian A., Taber L.A., Perucchio R., 2010, Automatic
generationofusermaterial subroutines forbiomechanicalgrowthanalysis,Journal of Biomechanical
Engineering, 132, 10, doi: 10.1115/1.4002375

Manuscript received April 13, 2016; accepted for print January 16, 2017