Download ANSYS_NeoHook_usermat PDF

TitleANSYS_NeoHook_usermat
TagsElasticity (Physics) Building Engineering Materials Stress (Mechanics) Deformation (Mechanics)
File Size181.5 KB
Total Pages15
Table of Contents
                            Neo--Hookean hyperelastic material
	Stress
	Tangent stiffness
Corotational frame
	Polar decomposition
	Stress
	Tangent stiffness
Voigt notation
	Stress
	Tangent stiffness
	Change of basis
Implementation
	Sample simulation
		Input
		Output
	Code listing
                        
Document Text Contents
Page 1

A large deformation Neo–Hookean user material in ANSYS.

ANSYS, Inc.
275 Technology Drive, Canonsburg, PA 15317

December 2008

Abstract

The development and implementation of a 3–dimensional, large deformation, hyperelastic user ma-
terial routine for ANSYS is presented. The constitutive model development follows the well–known
Neo–Hookean hyperelastic formulation which yields the stress and material tangent in the spatial config-
uration. The model is recast in the corotated frame required by ANSYS. An algorithm for determination
of the corotated frame rotation via the polar decomposition of the deformation gradient is also presented
along with the conversion of the stress, material tangent and change of basis operations to Voigt nota-
tion. The Fortran code listing for the user material routine and a sample input file simulating the simple
shearing of a block are given.

1 Neo–Hookean hyperelastic material

The strain energy potential for a nearly incompressible hyperelastic material is

W ≡ Φ
(


)
+ φ (J) , (1-1)

where C̄ = F̄ T F̄ is the isochoric right Cauchy–Green deformation tensor, F̄ = J−1/3F is the isochoric
deformation gradient, F is the deformation gradient and J = det (F ). The 2nd Piola–Kirchhoff stress is

S = 2
∂W

∂C
= 2

∂Φ
∂C

+
∂φ

∂J
JC−1, (1-2)

and material stiffness tensor is

C = 2
∂S

∂C
= 4

∂2Φ
∂C∂C

+ 2


∂C

(
∂φ

∂J

)
⊗ JC−1 + 2

∂φ

∂J



∂C

(
JC−1

)
. (1-3)

For a Neo–Hookean material

Φ
(


)
=

µ

2
(
ĪC − 3

)
,

φ (J) =
K

2
(J − 1)2 , (1-4)

where ĪC = trace
(


)
, and µ and K are the Neo–Hookean shear and bulk moduli. The 1st order derivatives

in Eqs. (1-2) and (1-3) are

∂Φ
∂C

=
µ

2
J−2/3

(
I −

IC
3

C−1
)

,

∂φ

∂J
= K (J − 1) ,



∂C

(
JC−1

)
= J

(
C−1 ⊗C−1 − 2I

)
, (1-5)

where IC = trace (C), I is the 2nd order identity tensor, and I is a 4th order tensor with components
Iijkl = C−1ik C

−1
jl . The 2

nd order derivatives in Eqs. (1-2) and (1-3) are

∂2Φ
∂C∂C

=
µ

2
J−2/3

3

[
−C−1 ⊗ I +

IC
3

C−1 ⊗C−1 − I ⊗C−1 + ICI
]

,



∂C

(
∂φ

∂J

)
=

J

2
KC−1. (1-6)

Page 2

LARGE DEFORMATION USER MATERIAL IN ANSYS DECEMBER 2008

1.1 Stress

Substituting Eq. (1-5) into (1-2), the 2nd Piola–Kirchhoff stress is then

S = µJ−2/3
(

I −
IC
3

C−1
)

+ KJ (J − 1) C−1, (1-7)

and the Cauchy stress is

σ =
1
J

τ =
1
J

FSF T

=
1
J

(
µJ−2/3

(
b−

IC
3

I

)
+ KJ (J − 1) I

)
,

=
1
J

µ

(
b̄−

ĪC
3

I

)
+ K (J − 1) I, (1-8)

where τ is the Kirchhoff stress, b = FF T is the left Cauchy–Green deformation tensor and b̄ = F̄ F̄ T .

1.2 Tangent stiffness

Substituting Eqs. (1-5) and (1-6) into (1-3), the material tangent stiffness is

C =
2
3
µJ−2/3

[
−C−1 ⊗ I − I ⊗C−1 + ICI +

IC
3

C−1 ⊗C−1
]

+

J2KC−1 ⊗C−1 + KJ (J − 1)
(
C−1 ⊗C−1 − 2I

)
(1-9)

or, with some simplification,

C =
2
3
µJ−2/3

[
−C−1 ⊗ I − I ⊗C−1 + ICI +

IC
3

C−1 ⊗C−1
]

+

KJ
[
(2J − 1) C−1 ⊗C−1 − 2 (J − 1) I

]
. (1-10)

Using the Piola transform, the spatial tangent is

c =
2
3
µJ−1

[
−b̄⊗ I − I ⊗ b̄ + ĪCi +

ĪC
3

I ⊗ I
]

+ K [(2J − 1) I ⊗ I − 2 (J − 1) i] . (1-11)

where i is a 4th order identity tensor with components iijkl = IikIjl.

2 Corotational frame

Large deformation in ANSYS is formulated in a corotated frame given by the rotation, R, from the polar
decomposition of the deformation gradient

F = RU , (2-1)

where U is the stretch tensor. The rotation of the corotated frame is then

R = FU−1. (2-2)

The stress and tangent stiffness tensors are returned to ANSYS in the corotated frame.

©2008 ANSYS, Inc. 2

Page 7

LARGE DEFORMATION USER MATERIAL IN ANSYS DECEMBER 2008

3.3 Change of basis

For an arbitrary change of basis given by the orthogonal tensor T , the components of the 2nd order stress
tensor and the 4th order tangent stiffness tensors transform as

σ′ij = TimσmnTjn,
∇′
c ijkl = TimTjn


cmnop TkoTlp. (3-11)

And in Voigt notation,

σ′α = Qαξσξ,
∇′
c αβ = Qαξ


c ξη Qβη. (3-12)

Comparing Eqs. (3-11) and (3-12), the change of basis tensor Q is given by

Q [T ] =
[
[P1] 2 [P2]
[P3] [P4]

]
, (3-13)

where

[P1] =


T 211 T 212 T 213T 221 T 222 T 223

T 231 T
2
32 T

2
33


 , [P2] =


T11T12 T12T13 T11T13T21T22 T22T23 T21T23

T31T32 T32T33 T31T33


 ,

[P3] =


T11T21 T12T22 T13T23T21T31 T22T32 T23T33

T31T11 T32T12 T33T13


 , [P4] =


T11T22 + T12T21 T12T23 + T13T22 T11T23 + T13T21T21T32 + T22T31 T22T33 + T23T32 T21T33 + T23T31

T31T12 + T32T11 T32T13 + T33T12 T31T13 + T33T11


 .
(3-14)

For a change of basis from the spatial frame to the corotated frame, T = RT .

©2008 ANSYS, Inc. 7

Page 8

LARGE DEFORMATION USER MATERIAL IN ANSYS DECEMBER 2008

APPENDIX

A Implementation

Shown below in §A.2 is the usermat3d() subroutine that implements the Neo–Hookean Cauchy stress and
tangent stiffness calculations of Eqs. (1-8) and (3-10). Two auxiliary subroutines to perform the polar
decomposition and form the 6–dimensional change of basis matrix are also given. The coding is done in fixed
source form modern Fortran and has favored clarity over computational efficiency.

The implementation is 3-dimensional and, as presented above, a displacement only formulation. It
requires the Neo-Hookean Young’s Modulus and Poisson’s ratio as the first and second material properties
respectively and has no user state variables. The stored energy used by ANSYS as an output variable is
calculated as Eq. (1-1) with (1-4) and the dissipated energy is zero. The transverse shear stiffness terms are
arbitrarily set to the (6, 6) component of the material tangent stiffness in the corotated frame.

A sample simulation input file and solution output are given below in §A.1. It simulates a cycle of simple
shear of a block that has the initial dimensions 1× 1× 1. The block is discretized with a 4× 4× 4 array of
SOLID185, 8–node, trilinear, hexahedral elements. The Neo–Hookean Young’s Modulus and Poisson’s ratio
are 3× 107 and 0.3 respectively. The simulation uses a fixed 10 substeps for the loading and unloading steps
and it writes the time, effective stress and effective strain, for an element near the center of the block, to a
text file named stressstrain usr.log.

A.1 Sample simulation

A.1.1 Input

!*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
! Sample input file to demonstrate the large deformation Neo-Hookean
! user material implementation. Simulates a simple shear cycle of a
! block with initial dimensions 1x1x1.
!*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
!---------------------------------------------------------
! Preprocessing
!---------------------------------------------------------
/PREP7

! Geometry and mesh
Lx=1
Ly=1
Lz=1
BLOCK,0,Lx,0,Ly,0,Lz
ET,1,SOLID185
ESIZE,Lx/4.
MAT,1
VMESH,1

! Material

! User material parameters
E = 30E6 ! Young’s Modulus
nu = 0.3 ! Poisson’s ratio
TB,USER,1,1,2,
TBDATA,1,E,nu

! Neo-hookean hyperelastic parameters
!G = E/2/(1+nu) ! shear modulus

©2008 ANSYS, Inc. 8

Page 14

LARGE DEFORMATION USER MATERIAL IN ANSYS DECEMBER 2008

!right Cauchy-Green
C(:,:) = matmul(transpose(F),F)
CC(:,:) = matmul(C,C)

!invariants
IC = C(1,1) + C(2,2) + C(3,3)
IIC = 0.5_RK*(IC**2 - (CC(1,1) + CC(2,2) + CC(3,3)))
IIIC = C(1,1) * (C(2,2)*C(3,3) - C(2,3)*C(3,2))
& + C(1,2) * (C(2,3)*C(3,1) - C(2,1)*C(3,3))
& + C(1,3) * (C(2,1)*C(3,2) - C(2,2)*C(3,1))

! eigenvalues of sqrt(C)
p = IIC - (IC**2)/3._RK
q = -(2._RK/27._RK)*IC**3+IC*IIC/3._RK-IIIC
if(abs(p)<epsilon(1._RK))then
l1 = sqrt( abs(-abs(q)**(1._RK/3._RK) + IC/3._RK) )
l2 = l1
l3 = l2

else
m = 2._RK*sqrt(abs(p)/3._RK)
n = 3._RK*q/(m*p)
if(abs(n)>1._RK) n = sign(1._RK, n)
t = atan2(sqrt(1._RK-n**2),n)/3._RK
l1 = sqrt( abs(m*cos(t) + IC/3._RK) )
l2 = sqrt( abs(m*cos(t+2._RK/3._RK*pi) + IC/3._RK) )
l3 = sqrt( abs(m*cos(t+4._RK/3._RK*pi) + IC/3._RK) )

endif

!stretch and inverse
!invariants
IU = l1 + l2 + l3
IIU = l1*l2 + l1*l3 + l2*l3
IIIU = l1*l2*l3
D = IU*IIU-IIIU

U(:,:) = (-CC(:,:) + (IU**2-IIU)*C(:,:) + IU*IIIU*I(:,:))/D
Ui(:,:) = (C(:,:) - IU*U(:,:) + IIU*I(:,:))/IIIU

! Rotation
R(:,:) = matmul(F,Ui)

end subroutine polarRU
!----------------------------------------------------------------------

SUBROUTINE trans_matrx_6(Q,R)
!----------------------------------------------------------------------
! Construct the 6x6 Voigt rotation matrix
!----------------------------------------------------------------------

IMPLICIT NONE
INTEGER,PARAMETER :: RK=KIND(1.D0) ! real kind

! arguments
REAL(RK), INTENT(OUT) :: Q(6,6)
REAL(RK), INTENT(IN) :: R(3,3)

!----------------------------------------------------------------------
Q(1,1) = R(1,1)**2

©2008 ANSYS, Inc. 14

Page 15

LARGE DEFORMATION USER MATERIAL IN ANSYS DECEMBER 2008

Q(2,1) = R(2,1)**2
Q(3,1) = R(3,1)**2
Q(4,1) = R(1,1)*R(2,1)
Q(5,1) = R(2,1)*R(3,1)
Q(6,1) = R(3,1)*R(1,1)

Q(1,2) = R(1,2)**2
Q(2,2) = R(2,2)**2
Q(3,2) = R(3,2)**2
Q(4,2) = R(1,2)*R(2,2)
Q(5,2) = R(2,2)*R(3,2)
Q(6,2) = R(3,2)*R(1,2)

Q(1,3) = R(1,3)**2
Q(2,3) = R(2,3)**2
Q(3,3) = R(3,3)**2
Q(4,3) = R(1,3)*R(2,3)
Q(5,3) = R(2,3)*R(3,3)
Q(6,3) = R(3,3)*R(1,3)

Q(1,4) = R(1,1)*R(1,2)*2._RK
Q(2,4) = R(2,1)*R(2,2)*2._RK
Q(3,4) = R(3,1)*R(3,2)*2._RK
Q(4,4) = R(1,1)*R(2,2)+R(1,2)*R(2,1)
Q(5,4) = R(2,1)*R(3,2)+R(2,2)*R(3,1)
Q(6,4) = R(3,1)*R(1,2)+R(3,2)*R(1,1)

Q(1,5) = R(1,2)*R(1,3)*2._RK
Q(2,5) = R(2,2)*R(2,3)*2._RK
Q(3,5) = R(3,2)*R(3,3)*2._RK
Q(4,5) = R(1,2)*R(2,3)+R(1,3)*R(2,2)
Q(5,5) = R(2,2)*R(3,3)+R(2,3)*R(3,2)
Q(6,5) = R(3,2)*R(1,3)+R(3,3)*R(1,2)

Q(1,6) = R(1,1)*R(1,3)*2._RK
Q(2,6) = R(2,1)*R(2,3)*2._RK
Q(3,6) = R(3,1)*R(3,3)*2._RK
Q(4,6) = R(1,1)*R(2,3)+R(1,3)*R(2,1)
Q(5,6) = R(2,1)*R(3,3)+R(2,3)*R(3,1)
Q(6,6) = R(3,1)*R(1,3)+R(3,3)*R(1,1)
END SUBROUTINE trans_matrx_6

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

©2008 ANSYS, Inc. 15

Similer Documents