The exponentiated Hencky energy:
Anisotropic extension and case studies
Jörg Schröder1 and Markus von Hoegen1 and Patrizio Neff2
1

Institut für Mechanik, Fakultät für Ingenieurwissenschaften / Abtl. Bauwissenschaften
Universität Duisburg-Essen, 45141 Essen, Universitätsstr. 15, Germany

arXiv:1702.00394v2 [] 27 Sep 2017

e-mail: j.schroeder@uni-due.de, phone: +49 201 183 2708, fax: +49 201 183 2680
e-mail: markus.von-hoegen@uni-due.de, phone: +49 201 183 3091, fax: +49 201 183 2680
2

Lehrstuhl für Nichtlineare Analysis und Modellierung, Fakultät für Mathematik
Universität Duisburg-Essen, 45127 Essen, Thea-Leymann-Straße 9, Germany
e-mail: patrizio.neff@uni-due.de, phone: +49 201 183 4243, fax: +49 201 183 4394

Abstract
In this paper we propose an anisotropic extension of the isotropic exponentiated Hencky energy, based on logarithmic strain invariants. Unlike other elastic formulations, the isotropic
exponentiated Hencky elastic energy has been derived solely on differential geometric grounds,
involving the geodesic distance of the deformation gradient F to the group of rotations. We
formally extend this approach towards anisotropy by defining additional anisotropic logarithmic strain invariants with the help of suitable structural tensors and consider our findings for
selected case studies.

Keywords: Hencky energy, logarithmic strain tensor, anisotropy, strain invariants,
biomechanics

1 . Introduction
In this article we consider a novel Hencky-type hyperelasticity model, the exponential Hencky-logarithmic strain energy proposed by Neff et al. [46], Neff et al. [47] and
Neff and Ghiba [42]. Here, we focus on an extension to anisotropy in a coordinate
invariant setting. Therefore, we apply the concept of structural tensors and introduce
additional mixed invariants. The flexibility of the proposed formulation is demonstrated
by identifying the linearized fourth-order elasticity tensor with the well-known coordinate
dependent representations. Furthermore, we propose an anisotropic exponential Hencky
model suitable for the description of soft biological tissues. The performance of this
model is demonstrated by the analysis of a patient specific artery.
The modeling of anisotropic soft tissues in the framework of nonlinear elasticity has made
considerable progress in the last decades. From the mathematical side, the polyconvexity
condition introduced by John Ball in his seminal paper Ball [4] is a strong mathematical
requirement which implies Legendre-Hadamard ellipticity (rank-one convexity) at all
deformation gradients F . In the early time after its introduction, polyconvexity was
exclusively used in the isotropic setting and it was unclear how to extend the framework
to anisotropy. In Hartmann and Neff [21] a large variety of isotropic strain energy
functions have been discussed. Two of the present authors have been able to solve
one of Ball’s major open problems, see Ball [5], namely the meaningful application of

The exponentiated Hencky energy: Anisotropic extension

2

polyconvexity to anisotropic materials. In a series of papers, Schröder and Neff [55],
Schröder et al. [57], Balzani et al. [6], Schröder et al. [58], Ebbing et al. [13], the concept
has been elaborated in detail, with papers from other authors following in due course,
see e.g. Itskov et al. [32], Ehret and Itskov [14]. It is, however, clear that polyconvexity
(or ellipticity) alone is not sufficient to characterize physically reasonable material
response: there exist polyconvex strain energies with unacceptable non-monotone
pressure-compression relation. Moreover, the identification of material parameters in the
above proposed anisotropic extensions heavily relies on nonlinear optimization strategies
after which the physical meaning of the obtained parameters is doubtful: a completely
different set of material parameters may equally well fit the available experimental data.
This already occurs for the isotropic Ogden-model, see Ogden [49]. The situation for the
anisotropic response can only be worse in general. Therefore, the need is to construct
strain-energy functions whose possibly few parameters have a clear physical meaning and
which are uniquely and easily identified from experiments. At the same time the proposed
strain energy should be Legendre-Hadamard elliptic at least in that range of deformation
which is typically encountered in the applications. In this paper we numerically explore
such a formulation based on the well known logarithmic Hencky strain.
In 1928 Heinrich Hencky [22] proposed the after him named strain-energy function WH for
finite isotropic elasticity. He replaced the small strain tensor ε in classical linear isotropic
energy by the Hencky or logarithmic strain measure log U , with the right stretch tensor U . For moderate deformations, this simple function WH with the two classical Lamé
constants is useful for a wide class of materials, see Anand [1; 2] and Bruhns et al. [11].
However, Hencky’s energy function is not rank-one convex, i.e. it does not fulfill the well
known Legendre-Hadamard, or ellipticity, condition. Bruhns et al. [12] derived necessary
and sufficient conditions for ellipticity in terms of principal stretches and computed the
largest common ellipticity region. They showed, in the case of positive Lamé constants,
that WH is elliptic whenever every principal stretch is in the range [0.21162..., 1.39561...].
Furthermore, Hencky’s strain-energy automatically satisfies the Baker-Ericksen inequality (Baker and Ericksen [3]) and Hill’s inequality (Hill [25; 26]), see Bruhns et al. [12],
Ghiba et al. [16].
1.1

The exponentiated Hencky energy

The exponentiated Hencky-logarithmic model was recently introduced by Neff et al. [46].
It is induced by the exponentiated Hencky strain energy
n
o


κ
µ
exp k̂ [tr(log U )]2
WeH (F ) = exp k kdevn log U k2 +
k
2#k̂
"
2
h
i
µ
U
κ
2
= exp k log
+
exp
k̂
(log
det
U
)
k
(det U )1/n
2k̂
where µ > 0 is the (infinitesimal) shear modulus, κ > 0 is
√ the bulk modulus, k and
k̂ are additional dimensionless material parameters, U = F T F is the right stretch
tensor corresponding to the deformation gradient F , log denotes the principal matrix
logarithm on the set of positive definite symmetric matrices, devn X = X − n1 trX1 and
p
kXk = tr(X T X) are the deviatoric part and the Frobenius matrix norm of an n × nmatrix X, respectively, and tr denotes the trace operator.

Schröder, von Hoegen, Neff

3

The exponentiated Hencky energy is based on the so-called volumetric and isochoric
logarithmic strain measures ωiso = kdevn log U k and ωvol = |tr log U | = | log det U | ,
respectively which have recently been given a purely geometric characterization not shared
by any other set of isotropic invariants (Neff et al. [48]): consider the general linear group
GL(n) endowed with the canonical left-invariant Riemannian metric g, which for A ∈
GL(n) and X, Y ∈ gl(n) = TA GL(n) ∼
= Rn×n is given by Martin and Neff [37]
gA (X, Y ) = hA−1 X, A−1 Y i
where hX, Y i = tr(XY T ) is the canonical inner product on the space of real n × nmatrices. Then the logarithmic strain measures are the geodesic distance of the isochoric
part detFF 1/n and the volumetric part (det F )1/n 1 to 1 of the deformation gradient to the
special orthogonal group SO(n), respectively, see Neff et al. [48] (Theorem 3.7):


F
, SO(n) ,
kdevn log U k = distgeod
(det F )1/n

|log det F | = distgeod (det F )1/n · 1, SO(n) .
(1)
These two quantities are thereby identified as the “natural” measures of strain in any
deformation, an observation which strongly suggests that an idealized elastic strain energy
function may depend on these quantities alone.1 An important example of such an energy
function is the classical quadratic Hencky energy
WH (F ) = µ kdevn log U k2 +
= µ log

U
(det U )1/n

κ
[tr(log U )]2
2
2
κ
+ (log det U )2 ,
2

(2)

which was introduced by Heinrich Hencky in 1929 (Hencky [23], Neff et al. [45]). While
the elasticity model induced by the Hencky energy is in very good agreement with
experimental observations for up to moderate strains for a large number of materials
Neff et al. [46], Anand [1], there are some major shortcomings of this model. For example,
the qualitative behavior of materials under very large deformations is not modeled accurately, and since the energy function is neither polyconvex nor quasiconvex or rank-one
convex (Neff [41], Ghiba et al. [16]), no known methods are available to ensure the existence of energy minimizers for general boundary value problems. Moreover, the pressurecompression relation is not monotone.
In order to alleviate some of these shortcomings, Neff et al. introduced the exponentiated Hencky energy WeH in a series of articles Neff et al. [46; 47], Neff and Ghiba [42],
Ghiba et al. [17]. This energy function closely approximates the classical quadratic
Hencky energy for small deformations, but aims to provide a more accurate model
for large deformations as well as an improvement in terms of common constitutive requirements; for example, WeH is polyconvex in the two-dimensional case Neff et al. [47],
and in the three-dimensional case the rank-one convexity domain contains the extremely large set {F ∈ GL+ (3)| kdev3 log U k ≤ 6}. Moreover, the induced mapping
1

Note that not every objective and isotropic energy function can be expressed in terms of the logarithmic strain measures alone, see Neff et al. [48], whereas every such energy can be expressed in terms
of the logarithmic strain tensor log U .

The exponentiated Hencky energy: Anisotropic extension

4

B 7→ σ of the Finger tensor B = F F T to the Cauchy stress tensor σ is invertible
Neff and Mihai [43], Mihai and Neff [38; 39], Jog and Patil [34], as is the case for suitable
variants of the Neo-Hooke and Mooney-Rivlin energies for slightly compressible materials
like rubber.
The low number of additional material parameters in the exponentiated Hencky model also
suggests that a good material fitting could be possible even without extensive experimental
measurements. Additionally, the exponentiated Hencky energy allows for the modeling of
3 κ−2 µ
in the finite strain regime: if the additional
a zero apparent Poisson’s modulus ν = 2(3
κ+µ)
parameters k, k̂ are chosen such that k = 23 k̂, then WeH can be written as





E
E
2
1
2
2
,
exp k kdevn log U k +
exp k (log det U )
2k 1 + ν
2 (1 − 2 ν)
3
κµ
where E = 39κ+µ
is Young’s modulus, and for ν = 0 we obtain a model with zero lateral
contraction under finite strains Neff et al. [46].

A variant of the exponentiated Hencky energy has previously been applied to so-called
tire derived materials and was found to be in good agreement with experimental data, see
Montella et al. [40]. In particular, the highly nonlinear equation of state (EOS) relating
pressure to purely volumetric deformations has been captured extraordinarily well. The
extra appearing non-dimensional parameters k and k̂ have an intuitive meaning: larger
k, k̂ lead to monotonically increased strain hardening. In principal these parameters can
be fitted independent of the shear and bulk modulus. Next, we extend the exponential
Hencky energy to the anisotropic case.

2 . Theoretical framework
2.1

Kinematics

For a better overview, the continuum-mechanical kinematic and constitutive quantities
are listed in Table 1. Let B ⊂ IR3 be the body of interest in the reference placement,
parametrized in X, and let S be the body in the current placement, parametrized in
x. The boundary ∂B of B is decomposed in ∂Bu and ∂Bt with ∂Bu ∪ ∂Bt = ∂B and
∂Bu ∩∂Bt = ∅. The nonlinear deformation map is given by x = ϕ(X). As basic kinematical
quantities we define the deformation gradient and the right Cauchy-Green tensor
F = Grad ϕ(X)

and

C = F TF = U2 ,

(3)

respectively. Here, 1 denotes the second-order identity tensor. The Jacobian of the deformation gradient has to satisfy J := det F > 0. The deformation gradient may be split
into
F = RU = V R,
(4)
where R ∈ SO(3) denotes a pure rotation tensor and U and V are the right and left
stretch tensors, respectively. In order to fulfill the principle of material frame indifference
a priori, we formulate the generic free energy function ψ in terms of the right CauchyGreen tensor, i.e. ψ = ψ(C). In spectral decomposition the right Cauchy-Green tensor C
and the left Cauchy-Green tensor B may be written as
C=

3
X
k=1

bk N k ⊗ N k ,
λ

B=

3
X
k=1

bk nk ⊗ nk
λ

(5)

Schröder, von Hoegen, Neff

5

Table 1: Kinematic and constitutive quantities.
Symbol

Continuum mechanical description

u
F = 1 + Gradu
C = FTF
B=√
FFT
U = √F T F
V = FFT
log U
log V

displacement vector
deformation gradient
right Cauchy-Green tensor
left Cauchy-Green tensor
right stretch tensor
left stretch tensor
right Hencky strain tensor
left Hencky strain tensor

ψ
WH
WeH

generic elastic energy
isotropic Hencky energy
exponentiated Hencky energy

S = 2∂C ψ(C)
τ = ∂log V ψ(log V )
σ = (det F )−1 τ

second Piola-Kirchoff stress tensor
Kirchoff stress tensor (see Vallée [61])
Cauchy stress tensor

bk denote the eigenvalues of C and B. The eigenvectors are expressed through N k
where λ
and nk associated to C and B, respectively. Therefore, we obtain the tensor functions
3
3 q
X
X
√
k
k
b
λk N k ⊗ N k ,
U= C=
λk N ⊗ N =

V =

√
B=

k=1
3 q
X

k=1
3
X

k=1

k=1

bk nk ⊗ nk =
λ

√

3

(6)

λk nk ⊗ nk ,

(7)
3

X1
X
1
!
bk N k ⊗ N k =!
log λk N k ⊗ N k ,
log U = log( C) = log C =
log λ
2
2
k=1
k=1
√

3

(8)

3

X1
X
1
!
bk nk ⊗ nk =!
log V = log( B) = log B =
log λk nk ⊗ nk ,
log λ
2
2
k=1
k=1

(9)

where λk denote the eigenvalues of U and V . The tensor log U is called right Hencky
strain tensor.
2.2

Stress measures

b
Let ∂log U ψ(log
U ) be the stress measure work conjugate to log U = 21 log C, then the
transformation rule for the second Piola-Kirchhoff stress tensor S is given by
b
b
b
∂ ψ(log
U)
∂ ψ(log
U ) ∂log U
∂ ψ(log
U)
=2
:
=
: PH
∂C
∂ log U
∂C
∂ log U
∂ log U
with PH = 2
.
∂C
S=2

(10)

The exponentiated Hencky energy: Anisotropic extension

6

The fourth-order tensor PH can only be derived with the help of the spectral decomposition
and yields

PH =

3 X
3
X

Pkkjj N k ⊗ N k ⊗ N j ⊗ N j

k=1 j=1

+2

3 X
3
X

Pkjkj (N k ⊗ N j ) ⊗ (N k ⊗ N j + N j ⊗ N k )

k=1 k6=j

Pkkjj
Pkjkj

(11)

b−1
= δkj λ
 1 kb 1 b
logλ − logλj

 2 λbk −λb2
j
k
 −1
1
=
bk − 1 log λ
bj
log λ
1
2
2
b
bk

(
=:
∂
log
λ
)
=
2λ
lim
bk 2
k
b b
b k −λ
bj
λ
λ
λk →λj

bk =
bj
for λ
6 λ
bk = λ
bj
for λ

see also Ogden [50] and Simo [59]. The first part of PH is related to the derivative of the
eigenvalues of log U with respect to C, while the second part is related to ∂C (N k ⊗ N k ).
In the isotropic case the following relations regarding the Kirchoff stress
e
b
∂ ψ(log
U) T
∂ψ # (U )
∂ ψ(log
V)
=R
R =R
U RT
∂log V
∂log U
∂U
+
∂ ψ̄(C) T
∂ψ (B)
∂ψ(F ) T
= 2F
F = 2B
=
F
∂C
∂B
∂F

τ =

(12)

hold true. But if anisotropic behavior is considered only the relations
τ = 2F

∂ ψ̄(C) T
∂ψ(F ) T
F =
F
∂C
∂F

(13)

remain valid. For the derivation of the above mentioned relations the reader is referred
to the appendix. The expression TBiot = ∂U ψ # (U ) is also known as Biot-stress. For the
linearization of the weak form, we need the tangent moduli
C=4

b
b
U)
∂ 2 ψ(log
∂ ψ(log
U)
= PH : CH : PH +
:K
∂C∂C
∂log U

(14)

with
∂ 2 ψ̂(log U )
∂PH
∂ 2 log U
and K = 2
=4
.
(15)
∂ log U ∂ log U
∂C
∂C∂C
The multiplicative volumetric isochoric decomposition of the deformation gradient
CH =

Fe = J −1/3 F

e = J −2/3 C
and C

(16)

was first proposed by Hans Richter [52], see also Flory [15]. In doing so we can express
the volumetric Hencky strain tensor with help of an additive split according to
log U = dev(log U ) +

1
tr(log U )1 ,
3

with tr(dev log U ) = 0 .

(17)

For the numerical treatment of an energy function ψ ∗ (dev log U ) we need the derivative
∂dev log U
1
= 1  1 − 1 ⊗ 1 = IP ,
∂log U
3

(18)

Schröder, von Hoegen, Neff

7

where  denotes the Kronecker product of second-order tensors. Let G and H denote
two second-order tensors and g and h two first-order tensors, then the operator is defined
by (G  H) : (g ⊗ h) = (Gg) ⊗ (Hh). Formulating a strain energy in dev log U , before
projecting the stress tensor and tangent moduli on C we first need to project them on
the Hencky strain log U . Therefore, we define
∂ψ ∗ (dev log U ) ∂dev log U
∂ψ ∗ (dev log U )
∂ψ ∗ (dev log U )
=
:
=
: IP ,
∂log U
∂dev log U
∂log U
∂dev log U

(19)

and for the linearization of the weak form it follows
CH =

∂ 2 ψ ∗ (dev log U )
e H : IP ,
= IP : C
∂ log U ∂ log U

(20)

∂ 2 ψ ∗ (dev log U )
.
∂ dev log U ∂ dev log U

(21)

with
eH =
C

The corresponding tensors are to be inserted in Eq. (10) and Eq. (14).
2.3

Isotropic and anisotropic invariants

The principal isotropic invariants of the right Cauchy-Green tensor C are given by
I1C = tr[C] = kF k2 ,
I3C

I2C = tr[Cof C] = kCof F k2
2

and

2

= det C = (det F ) = J .

(22)

Further, we introduce the basic invariants of the Hencky strain tensor log U
J1H = log(det U ) = tr(log U ) ,

J2H = klog U k2 = tr[(log U )2 ] and

J3H = tr[(log U )3 ] ,

(23)

already used by Richter [52]. Let A, with kAk = 1, be the preferred direction of the
transversely isotropic material, then the material symmetry group is defined by
Gti := {±1; Q(α, A) | 0 < α < 2π} ,

(24)

where Q(α, A) are all rotations along the A-axis. The structural tensor M whose invariance group preserves the material symmetry group Gti is given by the rank-one tensor
M = A ⊗ A,

(25)

see Boehler [8] and Boehler [9] regarding the concept of structural tensors. Based on the
structural tensor we define the mixed invariants
i

J5C = hCof C, M i ,

i

J5H = hlog(Cof U ), M i ,

I4C = hC i , M i ,
I4H = h(log U )i , M i ,

(26)

where i ∈ N, i > 0, denotes an exponent. Note that the cofactor Cof(log U ) has no
physical meaning and that for U ∈ Sym+ (3), Cof U is also positive definite. Because of

The exponentiated Hencky energy: Anisotropic extension

8

that we instead consider the logarithmic cofactor function
√ !
3
X


det
C
log(Cof U ) = log (detU )U −1 =
log
Nk ⊗ Nk
1/2
b
λk
k=1


3
X
√
1 b
k
k
=
log(det C) − logλ
k N ⊗N
2
k=1
1 b
1
b3 )N 1 ⊗ N 1 + ( 1 log λ
b1 + 1 log λ
b3 )N 2 ⊗ N 2
= ( logλ
log λ
2+
2
2
2
2
1
b1 + 1 log λ
b2 )N 3 ⊗ N 3
+ ( log λ
2
2
= tr(log U )1 − log U

(27)

and finally we observe the following
1

J5H = hlog(Cof U ), M i = tr(log U ) hM , 1i −hlog U , M i = J1H − I4H .
| {z }

(28)

1

3

Isotropic strain energy functions

3.1. Isotropic Hencky Energy
The isotropic Hencky energy was introduced in Hencky [23]. It measures the geodesic
distance of the deformation gradient to the special orthogonal group SO(n), as it was
discovered in Neff et al. [44]. The Hencky strain energy
WH (log U ) = µ kdev log U k2 +

κ
λ
[tr(log U )]2 = µk log U k2 + [tr(log U )]2
2
2

can be reformulated in principal logarithmic strains
"
2 
2 
2 #
1
1
1
b1 +
b2 +
b3
WH (log U ) = µ
log λ
log λ
log λ
2
2
2

2
1
1
λ 1
b
b
b
+
log λ1 + log λ2 + log λ3
2 2
2
2

(29)

(30)

based on the eigenvalues λ̂ of C, where the Lamé parameters λ and µ as well as the bulk
modulus κ are used. Note that κ = 3λ+2µ
and
3


s

2



s

2



s

2 

b1
b
b
1 
λ
 + log λ1  + log λ2  
log
b2
b3
b3
3
λ
λ
λ
"
2 
2 
2 #
2
1
1
1
b1 +
b2 +
b3
=
log λ
log λ
log λ
3
2
2
2


2 1
1
1
b
b
b
b
b
b
−
log λ1 log λ2 + log λ1 log λ3 + log λ2 log λ3 .
3 4
4
4

kdev (log U )k2 =

(31)

Schröder, von Hoegen, Neff

9

The function WH is not polyconvex, not quasiconvex, not coercive and not rank-oneelliptic, even for every admissible deformation state, see Neff et al. [46]. However, it holds
that WH (F ) = WH (F −1 ). The first and second derivative with respect to the Hencky
strain yield
∂WH
= 2µ dev (log U ) + κ tr(log U )1 ,
∂ log U
∂ 2 WH
CH =
= 2µ IP + κ 1 ⊗ 1 .
∂ log U ∂ log U

(32)
(33)

In the reference configuration with C = 1 the final material tangent C, in Voigt-notation2
denoted as CV , according to Eq. (14) simplifies to


κ + 34 µ κ − 23 µ κ − 23 µ 0 0 0
κ − 2 µ κ + 4 µ κ − 2 µ 0 0 0 
3
3
3


2
κ − 2 µ κ − 2 µ κ + 4 µ 0 0 0 
∂ WH
V
3
3
3

.
C |C=1 = 4
=
(34)

0
0
0
µ
0
0
∂C∂C C=1 


 0
0
0
0 µ 0
0
0
0
0 0 µ
3.2. Exponentiated Hencky energy
The exponentiated Hencky energy
i
h

κ
µ
2
2
WeH = exp k kdev (log U )k + exp k̂ (tr log U ) ,
k
2k̂
"
" 

 2#
2 #
1
1
µ
κ
1
1
= exp k dev
+ exp k̂ 1, log C
log C
, k > , k̂ >
(35)
k
2
2
3
8
2k̂
was introduced and described in Neff et al. [46]. It is still volumetric-isochoric decoupled
and polyconvex in 2D if k > 31 and k̂ > 18 , cf. Neff et al. [47]. Rank-one convexity is not
preserved in 3D, see Neff et al. [46]. However, numerical calculations show that the ellipticity domain contains the extremely large set {F ∈ GL+ (3)| kdev3 log U k ≤ 6}. In the
small strain regime for principal stretches λi ∈ (0.7, 1.4) it approximates the aforementioned isotropic Hencky energy quite well.
Reformulation in terms of the Lamé parameters µ and λ yields
"
" 

 2#
2 #
µ
1
(3λ + 2µ)/3
1
log C
+
exp k̂ 1, log C
.
WeH = exp k dev
k
2
2
2k̂

(36)

The derivatives with respect to the Hencky strain yield


∂WeH
= 2µ exp k kdev (log U )k2 dev (log U )
∂ log U
h
i
2
+ κ exp k̂ (tr log U ) (tr log U ) 1 ,
2

(37)

In the contracted notation the tensorial indices are allocated to the matrix indexes as follows
{11, 22, 33, 12, 23, 13} → {1, 2, 3, 4, 5, 6}.

The exponentiated Hencky energy: Anisotropic extension

10

"

 2#




2
∂
W
1
1
1
eH
H
C =
= 4µ exp k dev
log C
dev
log C ⊗ dev
log C
∂ log U ∂ log U
2
2
2
"
 2#

1
log C
IP
+ 2µ exp k dev
2
" 
2 #
1
+ κ exp k̂ 1, log C
1⊗1
2
" 
2 # 
2
1
1
+ 2κk̂ exp k̂ 1, log C
1, log C 1 ⊗ 1 .
(38)
2
2
In the reference configuration with C = 1 and log U = 0 the above equation simplifies to
∂ 2 WeH
∂ 2 log U

= 2µ IP + κ 1 ⊗ 1

(39)

C=1

and the final tangent according to Eq. (14) becomes


κ + 34 µ κ − 23 µ κ − 23 µ 0 0 0
κ − 2 µ κ + 4 µ κ − 2 µ 0 0 0 
3
3
3


κ − 2 µ κ − 2 µ κ + 4 µ 0 0 0 
∂ 2 WeH
V
3
3
3
,
C |C=1 = 4
=
0
0
µ 0 0
∂C∂C C=1 
 0

 0
0
0
0 µ 0
0
0
0
0 0 µ

(40)

which is identical to the elasticity tensor of WH , provided in Eq. (34).

4 . Anisotropic extension
4.1. Transverse isotropic Hencky and exponentiated Hencky models
1

2

In a first step we aim to investigate a strain energy function W1 (J1H , kdev log U k, I4H , I4H )
which basically extends the classical Hencky-strain energy by the basic mixed invariants,
introduced in chapter 2.3:
κ
W1 (log U ) = µT kdev log U k2 + [tr(log U )]2 + αhlog U , M i(tr log U )
2
1
+ 2(µL − µT )hM , (log U )2 i + βhM , log U i2 .
(41)
2
The parameters are chosen in analogy to Spencer [60]. Here, µT and µL are associated
to the shear moduli in the transverse isotropy plane and perpendicular to that, κ is
associated to the bulk modulus, β is associated to the stiffness in fiber direction. In a
small strain framework, i.e. replacing log U by ε = 21 (Grad u + GradT u), the above given
energy function would refer to transversely isotropic linear elasticity. In an exponential
framework the exponentiated transversely isotropic strain-energy function
κ
µT
W2 (log U ) =
exp[k1 kdev log U k2 ] +
exp[k2 (tr log U )2 ]
k1
2k2
α
2(µL − µT )
+ exp[k3 hlog U , M i(tr log U )] +
exp[k4 hM , (log U )2 i]
k3
k4
β
+
exp[k5 hM , log U i2 ]
(42)
2k5

Schröder, von Hoegen, Neff

11

will result in the same elasticity tensor C|C=1 at the identity, whereby ki > 0 are further
non-dimensional parameters. We also note that while the isotropic invariants have the
proposed differential geometric meaning, the novel exponential terms are formulated on
an ad hoc basis.
Due to the non-linearity of the above given equations we aim to identify the general
material parameters in the reference configuration with C = 1 and log U = 0. In case
of linearized transversely isotropic materials the elasticity tensor C may be formulated
in terms of five material parameters, see Eq. (43). If we choose the isotropic-plane to be
spanned by the X1 and X2 axis and the preferred direction to coincide with the coordinate
axis X3 , perpendicular to the isotropic plane, and M = diag(0, 0, 1) we may write
V
V
CV
11 C12 C13
V
CV
CV
11 C13
 12
V
V
C13 C13 CV
33
CV = 
 0
0
0

 0
0
0
0
0
0




0
0
0
0
0
0 

0
0
0 

V
1
(CV
0
0 
11 − C12 )

2

CV
0
0
44
V
0
0 C44

(43)

in Voigt-notation. In the reference configuration, the formulated transversely isotropic
strain energy functions W1 and W2 , both yield


λ + 2µT
λ
λ+α
0 0 0
 λ
λ + 2µT
λ+α
0 0 0


 λ+α

λ
+
α
λ
−
2µ
+
2α
+
4µ
+
β
0
0
0
T
L
V
,
C |C=1 = 
 0

0
0
µT 0


 0
0
0
0 µL 0 
0
0
0
0 0 µL
where the conversion λ = (3κ − 2µT )/3 was used. The components of the above presented
scheme are related to the five classical components in Eq. (43) through
µL = CV
44
1
V
µT = (CV
11 − C12 )
2
λ = CV
12

,

(44)

V
α = CV
13 − C12
V
V
V
β = CV
11 + C33 − 2C13 − 4C44

in analogy to Schröder and Gross [54]. The anisotropic characteristic of the strain energy
function W2 in Eq. (42) will be outlined on a number of numerical examples, where
different material parameters according to Table 2 are used.
Further, the sets 4, 5 and 6 are chosen such that only one term involving a structural tensor in Eq. (42) is active. Therefore, Set 4 is directly associated with the
term hlog U , M i(tr log U ), Set 5 with hlog U , M i and Set 6 with h(log U )2 , M i =
kdev(log U ) M k2 . In contrast, Set 3 will serve as the isotropic reference case. In order
to obtain comparable results a similar level of distinct anisotropy is chosen for each set.

The exponentiated Hencky energy: Anisotropic extension
Set
Set 1
Set 2
Set 3
Set 4
Set 5
Set 6

λ
µT
α
β
µL
k1
1000 175
10
10
375
1
5.64 2.64 1.27 0.29
5.66
1
5.5
2.5 0.00 0.00
2.5
1
5.5
14 40.75 0.00
14
1
5.5
2.5 0.00 104.5
2.5
1
5.5
2.5 0.00 0.00 28.625 1

k2
1
1
1
1
1
1

k3
1
1
75
75
75
75

k4
1
1
25
25
25
25

12
k5
1
1
45
45
45
45

Table 2: Different parameter sets for the numerical examples.

To achieve this, CV
33 is the same for all three sets regarding the reference configuration.
All parameter sets have been checked to be positive definite for CV |C=1 . Note that in
case of Set 4 the parameter µT needed to be increased in order to guarantee the positive
definiteness of CV |C=1 . The corresponding matrices are listed below:


33.5
 5.5

46.25
V
Set 4 : C |C=1 = 
 0

 0
0

10.5
 5.5

 5.5
V
Set 5 : C |C=1 = 
 0

 0
0

10.5
 5.5

 5.5
V
Set 6 : C |C=1 = 
 0

 0
0


5.5 46.25 0 0 0
33.5 46.25 0 0 0 

46.25 115 0 0 0 

0
0
14 0 0 

0
0
0 14 0 
0
0
0 0 14

5.5 5.5 0
0
0
10.5 5.5 0
0
0

5.5 115 0
0
0

0
0 2.5 0
0

0
0
0 2.5 0 
0
0
0
0 2.5

5.5 5.5 0
0
0
10.5 5.5 0
0
0 

5.5 115 0
0
0 

0
0 2.5
0
0 

0
0
0 28.625
0 
0
0
0
0
28.625

The implementation in a finite element framework in this work was done according to
the formulation in Schröder et al. [56] and Löblein et al. [36]. The weak form of balance
of momentum
Z
G(u, δu) := hDiv(∂F ψ(F )) + ρ0 (b − ẍ), δui dV
(45)
B

required for the finite element code may be reformulated such that we obtain


Z
Z
Z
G(u, δu) = hτ , ∇(δu)i dV −  hf , δui dA + hρ0 (b − ẍ), δui dV  = 0 .
|B

∂Bt

{z

Gint

}

|

(46)

B

{z

Gext

}

Here, the body force in the reference configuration is denoted by f , δu is the variation of
the displacement field and ẍ the acceleration. For the solution scheme in a finite element

Schröder, von Hoegen, Neff

13

framework a Newton iteration is required. Therefore, the linearization
LinG(ū, δu, ∆u) := G(ū, δu) + ∆G(ū, δu, ∆u)

(47)

at u = ū is required where the increment ∆G is defined through
Z
Z
∆G = h∇sym (δu), c : ∇sym (∆u)i dV + h∇(∆u) τ , ∇(δu)i dV ,
B

(48)

B

where the Eulerian tangent moduli c is obtained by the push-forward of the Lagrangian
tangent moduli C, i.e.
c = (F  F ) : C : (F T  F T )
(49)
and ∇sym (•) = 1/2 [grad(•) + gradT (•)]. For the numerical treatment the weak formulations of the aforementioned balance equation has been implemented in the finite element
analysis program FEAP of R.L. Taylor, University of California. For all of the following
examples quadratic triangular elements with six nodes per element were used.
Tension test: In a first example the transversely isotropic material behavior is to be
explained on the basis of a tensile test under plane strain conditions, see Fig. 1a). Material
parameter Set 1 was chosen and the preferred direction A, defined with help of the fiber
angle βf , was varied. The computed displacements of the nodes 1, 2 and 3 over the fiber
angle are plotted in Fig. 1b).The displacements δh1 , δh2 and δv3 are symmetric concerning
βf = 90◦ , while the vertical displacements δv1 and δv2 are antisymmetric. Surprisingly, the
horizontal displacements don’t reach their maximum value for βf = 90◦ . The load p0 has
been chosen such that large deformations are present.

2

A
βf

4
δhi |i=1,2 , δvi |i=1,2,3

10

3

1 p0 = 200
a)

10

δh2

δh1

3
2
1

δv1

0

δv2

-1

δv3

-2
0
b)

30

60

90
βf

120 150 180

Figure 1: Example 1: a) Tension test with one preferred direction and varying fiber orientation. b) Nodal displacements depending on the fiber orientations are plotted for parameter
Set 1.

Cooks Membrane: In a second example we consider the Cooks Membrane problem, as
depicted in Fig. 2a) which is dominated by non-homogenous stress distributions. Again
the fiber direction is to be varied and the body will undergo large deformations during
loading. In Fig. 2b), the vertical displacements δV of the node at the top right are plotted
for different fiber angles. Parameter Set 2 was considered. The anisotropic effect due to
the different fiber orientation clearly becomes apparent.

The exponentiated Hencky energy: Anisotropic extension

14

δV
1

16

0.8

A
βf
44

0.6

||F ||

F

βf = 0◦
βf = 30◦
βf = 60◦
βf = 90◦

0.4
0.2
0

a)

48

b)

0 2 4 6 8 10 12 14 16 18
δv

Figure 2: Example 2: a) Cooks Membrane with one preferred direction and varying fiber
orientation. b) Plot of the vertical displacement for different fiber orientations and parameter
Set 2.

Perforated Plate: Lastly, we want to consider three perforated plates, again under
consideration of plane strain conditions. The first two plates are squares and share the
same geometry. The first one, referred to as Example 3a, is subject to uniaxial tensile loads,
see Fig. 3a), while the second one, referred to as example 3b, is subject to biaxial tensile
loads, see Fig. 4a). In example 3c we consider a circular disk, depicted in Fig. 5a), which is
expanded on the inner ring, i.e. we apply a radial displacement of 3.5. All simulations are
displacement driven, i.e. only boundary conditions of Dirichlet type are present. Further,
all three bodies have one preferred direction A with an angle of 45◦ to the horizontal
axis. The deformed bodies for the parameter Sets 3, 4, 5 and 6 are plotted next to the
boundary value problems in Fig. 3, Fig. 4 and Fig. 5. The contour plots of the squares
show the horizontal displacements u1 . In order to highlight
the anisotropic characteristic
p
of the circular disk, the circumferential stretch λφ = hC, Nϕ ⊗ Nϕ i, where Nϕ denotes
the circumferential direction, is plotted for the third plate. The black lines in each plot
indicate the shape of the body in the undeformed configuration.
The comparison of the different parameter sets are intended to demonstrate different
anisotropic characteristics of different strain measures. Therefore, the parameter sets were
chosen such that only one anisotropic part of the strain energy in Eq. (42) is active. First
of all from Fig. 3 and Fig. 4 it becomes apparent that for the anisotropic Sets 4, 5 and
6 the displacements u1 are not symmetric with respect to the horizontal axis, which is
different for the isotropic Set 3. The different shapes of the holes very well emphasize
that the usage of different anisotropic invariants may lead to very different deformations,
even if the preferred direction is identical. In example 3c the inner ring of the disk is
exposed to predefined deformations, which is why the shapes of the holes are the same
for each parameter set. Nevertheless, the distribution of the circumferential stretch λϕ is
considerably different. Considering a polar coordinate system the circumferential stretch
is only depending on the radius, but independent of the polar angle for the isotropic Set
3. For Set 4 and Set 6 it can be seen that λϕ is smallest in the regions were the preferred
direction A and the circumferential direction Nϕ coincide and largest in the regions were
A and Nϕ are perpendicular. For Set 5 only the first of these two observations holds true.

Schröder, von Hoegen, Neff

15

We conclude that a considerable stiffening effect in the preferred direction is visible for
Set 4, 5 and 6.
3.75

3.75

3.75

u1

u1

1.25

A

Set 3

Set 4

3.75

45◦

b)

c)

a)
u1
Set 5

d)

Set 6

e)

Figure 3: Example 3a: a) Boundary conditions and b)-e) deformed bodies of a perforated
plate under uniaxial tension for parameter Set 3, 4, 5 and 6.

3.75

3.75

u2 = u1
3.75

u1

A

a)

u1

1.25

Set 3

Set 4

3.75

45◦

u2 = u1

b)

c)

u1
Set 5

d)

Set 6

e)

Figure 4: Example 3b: a) Boundary conditions and b)-e) deformed bodies of a perforated
plate under biaxial tension for parameter Set 3, 4, 5 and 6.

The exponentiated Hencky energy: Anisotropic extension

5.0

2.5 2.5

5.0

ur
A

a)

16

Set 3

Set 4

45◦

b)

c)

λϕ
Set 5

d)

Set 6

e)

Figure 5: Example 3c: a) Boundary conditions and b)-e) deformed bodies of a circular
perforated plate for parameter Set 3, 4, 5 and 6.

4.2

Orthotropic Hencky and exponentiated Hencky models

Orthotropic material behavior is symmetric regarding three orthogonal planes. These
three planes are spanned by the three preferred directions A1 , A2 and A3 , which are unit
vectors. The material symmetry group is then defined by
Go := {±1; R1 , R2 , R3 } ,

(50)

where R1 , R2 and R3 are the reflections with respect to the planes spanned by (A2 , A3 ),
(A1 , A3 ) and (A1 , A2 ), respectively. The three preferred directions Ai |i = 1, 2, 3 are orthogonal and form the three structural tensors Mi = Ai ⊗ Ai |i = 1, 2, 3, which satisfy the
orthogonality condition hMi , Mj i = δij . However, it is sufficient to formulate orthotropic
strain energies with two structural tensors M1 and M2 and additional isotropic principal
or main invariants. In the same manner as in the previous section we may formulate the
orthotropic energy function
κ
[tr(log U )]2
2
+ α1 hlog U , M1 i(tr log U ) + α2 hlog U , M2 i(tr log U )
+ 2µ1 hM1 , (log U )2 i + 2µ2 hM2 , (log U )2 i
1
1
1
+ β1 hM1 , log U i2 + β2 hM2 , log U i2 + β3 hM1 , log U ihM2 , log U i (51)
2
2
2

W3 (log U ) = µ kdev log U k2 +

Schröder, von Hoegen, Neff

17

and the corresponding exponentiated version
κ
µ
exp[k2 (tr log U )2 ]
W4 (log U ) = exp[k1 kdev log U k2 ] +
k1
2k2
α1
α1
+
exp[k3 hlog U , M1 i(tr log U )] +
exp[k4 hlog U , M2 i(tr log U )]
k3
k4
2µ2
2µ1
exp[k5 hM1 , (log U )2 i] +
exp[k6 hM2 , (log U )2 i]
+
k5
k6
β1
β
2
+
exp[k7 hM1 , log U i2 ] +
exp[k8 hM2 , log U i2 ]
2k7
2k8
β3
+
exp[k9 hM1 , log U ihM2 , log U i] .
(52)
2k9
The orthotropic elasticity tensor CV has nine independent variables. Choosing the structural tensors M1 = diag(1, 0, 0), M2 = diag(0, 1, 0) we obtain the general form
 V

V
C11 CV
0
0
0
12 C13
V
CV
CV
0
0
0 
22 C23
 12

V
V
V


C
C
C
0
0
0
V
13
23
33
.
(53)
C =
V
 0

0
0
C
0
0
44


 0
0
0
0 CV
0 
55
0
0
0
0
0 CV
66
The material tangent of both functions, W3 and

2µ + λ+
 2α1 + 4µ1 + β1 λ + α1 + α2 + β3


2µ + λ+
λ + α1 + α2 + β3

2α2 + 4µ2 + β2
CV |C=1 = 

λ
+
α
λ + α2
1


0
0


0
0
0
0

W4 both yield

λ + α1
λ + α2

0
0

0
0

0
0

2µ + λ
0
0
0
0
µ + µ2
0
0
0
0
µ + µ1
0
0
0
0
µ + µ1 + µ2













in the reference configuration, with λ = (3κ − 2µ)/3 . Following the same scheme as in
the transversely isotropic case the parameter identification gives the following relations
V
V
µ = CV
44 + C55 − C66
V
µ1 = C V
66 − C44
V
µ2 = C V
66 − C55
V
V
V
λ = CV
33 + 2(C66 − C44 − C55 )
V
V
V
V
α1 = CV
13 − C33 − 2(C66 − C44 − C55 )

α2 =
β1 =
β2 =
β3 =

CV
23
CV
11
CV
22
CV
12

−
+
+
−

CV
33
V
C33
CV
33
CV
13

−
−
−
−

V
V
2(CV
66 − C44 − C55 )
V
2CV
13 − 4C55
V
2CV
23 − 4C44
V
V
CV
23 + C33 + 2(C66 −

.

(54)

V
CV
44 − C55 )

A viscoelastic, orthotropic material model based on finite logarithmic strains has been
recently proposed by Latorre and Montáns [35].

The exponentiated Hencky energy: Anisotropic extension

4.3

18

Case study of the transversely isotropic model in logarithmic strain space

In order to study the anisotropic properties under compression we further introduce the
strain energy functions







µ1 

ti
i
2
ψ¶ =
exp k1 (hC , M i −1)  − 1 ,
(55)
C
| {z }

2k1 


i
I4C







µ1 
2

ti
i
(56)
ψ· =
exp k1 h(log U ) , M i  − 1 .
H
{z
}
|

2k1 


i
I4H

i

In the following we aim to investigate the performance of the anisotropic invariants I4C ,
i
I4H , see Eq. (26), respectively. The evolution of the invariants, as well as the stress response
i
i
of the transversely isotropic strain energy functions ψ ti (I4C ) and ψ ti (I4H ), are plotted for
different loading scenarios. The examples are evaluated such that the results are independent of any chosen isotropic strain energy function, since only the anisotropic stress
response will be plotted. The case distinction for compression and tension included in
the energy functions will be neglected, i.e. the fibers are allowed to induce stresses under
compression and the energy functions according to Eq. (55) and Eq. (56) will be used.
The parameters µ1 and k1 are set to one and the plotted evolution of the Cauchy stress
σ

aniso

∂ψ ti T
1
F
= 2F
J
∂C

(57)

will be normalized by the occurring maximum stress at the final deformation state, to
allow for a better comparison.
This study is restricted to classical homogenous deformation states, i.e. uniaxial tension
and compression, simple shear and biaxial loading conditions. During uniaxial tension,
uniaxial compression and biaxial loading the angles θk will remain constant and the body
is free of rotations, i.e. F = U . Only during the shear test they will change with a change
in the deformation.
4.3.1 Uniaxial tension and compression The considered problem is depicted in
Fig. 6. In this case the fiber direction is aligned with the loading direction. The reference
configuration as well as the deformed configurations under tension and compression are
shown on the right. The component F11 of the deformation gradient refers to the stretch
in fiber direction.
Since the body is considered to be incompressible we find that F22 =
√
F33 = 1/ F11 .
ti
The results for the transversely anisotropic Hencky function ψ·
are displayed in Fig. 7.
H
As already discussed in the previous section, in Fig. 7a) it becomes apparent that for
i
even exponents of i in I4H the values of the invariant are also positive under compression,
i.e. the sign of the invariant is not the right choice to distinguish between tensile and
compressive stretches. Nevertheless, the stress response seems to be adequate from a
1
physical point of view for each of the considered invariants. For I4H = hlog U , M i we
obtain a perfectly linear material behavior. Due to the logarithmic framework the stress

Schröder, von Hoegen, Neff

19

function is generally more sensitive to compression than to tension, see Fig. 7b). Note that
aniso
the stress function was normalized with the corresponding highest stress max(kσ11
)k
to allow for a better comparison. All stress-strain responses show the potential to exhibit
significant strain stiffening for both tension and compression. This is especially remarkable
for the compressive case. Hyun and Nakajima [31] for example found that porous copper
fabricated by unidirectional solidification behaves strongly anisotropic under compression
and exhibits considerable stiffening under large strains up to 80% due to the alignment
of the pores. Classical anisotropic material laws based on the invariant hC i , M i struggle
to reproduce this effect. The normalized stresses under compression according to Eq. (55)
are plotted in Fig. 8a). The plot reveals a strain softening behavior. Moreover, after a
certain point the stresses will begin to increase although the body is further compressed
which is strictly unphysical. Also classical polynomial laws of the form
µ1
(hC i , M i − 1)k1
2k1

ti
ψ¸
=
C

(58)

suffer from this effect, see Fig. 8b). Here, µ1 was set to one and k1 to two. In addition the
domain of definition of the above energy is restricted to even values of k1 in compression
which significantly limits the parameter fitting properties.
Z

A

X

comp.
Y

ref.
tension

Figure 6: Uniaxial compression and tension test, where the preferred direction and the
loading are aligned. The body is assumed to be incompressible with Vref = Vtension = Vcomp .

0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8

I4H

aniso
aniso
σ11
/max(||σ11
||)

I4H

i

The exponentiated Hencky energy: Anisotropic extension

i

i=1
i=2
i=3
i=4
0.6

0.8

a)

1
F11

1.2

0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

20

i

I4H

i=1
i=2
i=3
i=4

1.4

0.6

0.8

b)

1
F11

1.2

1.4

Figure 7: Monotonicity of the Cauchy stress σ11 as function of uniaxial stretch. Evolution
i
of a) the anisotropic Hencky-type invariants I4H for i = 1, 2, 3, 4 and b) corresponding
anisotropic stress response under uniaxial tension and compression for ψ ti . The stresses
·H
are normalized by the absolute value of the maximum arising compressive stress at F11 = 0.5.

0

i=1
i=2
i=3
i=4

-0.2
-0.4

aniso
aniso
σ11
/max(||σ11
||)

aniso
aniso
σ11
/max(||σ11
||)

0

-0.6
-0.8

I4C

i

-1

-0.2
-0.4

I4C

-0.6

i

-0.8
-1

0.5
a)

i=1
i=2
i=3
i=4

0.6

0.7 0.8
F11

0.9

1

0.5
b)

0.6

0.7 0.8
F11

0.9

Figure 8: Pathological non-monotonicity of the Cauchy stress σ11 as function of uniaxial
stretch. Normalized anisotropic stress response under uniaxial compression. In a) the exponential strain energy function ψ ti and in b) the polynomial function ψ ti is plotted.

¶C

¸C

1

Schröder, von Hoegen, Neff

21

4.3.2 Simple shear Next we investigate the behavior for simple shear, according to
Fig. 9. The shear direction will be aligned with the fiber direction and the amount of
shear
u
(59)
γ=
L
is defined as the quotient of the displacements by the length. Note that in this example
1
the fibers are not elongated at all, i.e. I4C = 1. The results for the transversely anisotropic
Hencky function are displayed in Fig. 10 and Fig. 11. Again the invariants of even and
odd powers take a different sign. The stress quantities which are not plotted in Fig. 11a)
aniso
are equal to zero. That means for even powers i, σ11
will be equal to zero.
Z

ref.
A
u

X

L
sheared

a)

0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5

1

i=1
i=2
i=3
i=4

I4H

i

aniso
aniso
σ13
/max(||σ13
||)

I4H

i

Figure 9: Shear test, where the preferred direction and the shear direction are aligned.
Simple shear deformation is incompressible with Vref = Vsheared .

i=1
0.8 i = 2
i=3
0.6 i = 4
0.4

i

I4H

0.2
0

0 0.2 0.4 0.6 0.8 1 1.2 1.4
amount of shear γ

b)

0 0.2 0.4 0.6 0.8 1 1.2 1.4
amount of shear γ
i

Figure 10: Evolution of the anisotropic Hencky-type invariants I4H for i = 1, 2, 3, 4 and b)
corresponding Cauchy stresses σ13 under simple shear. The stresses are normalized by the
absolute value of the maximum arising shear stress at γ = 1.5.

The exponentiated Hencky energy: Anisotropic extension
0

-0.2
-0.4

I4H

i

-0.6

-0.8 i = 1
i=3
-1
0 0.2 0.4 0.6 0.8 1 1.2 1.4
amount of shear γ
a)

aniso
aniso
σ33
/max(||σ33
||)

aniso
aniso
σ11
/max(||σ11
||)

0

22

b)

-0.2
-0.4

i

I4H

-0.6 i = 1
i=2
-0.8 i = 3
i=4
-1
0 0.2 0.4 0.6 0.8 1 1.2 1.4
amount of shear γ

Figure 11: Evolution of the Cauchy stresses a) σ11 and b) σ33 corresponding to the
anisotropic strain invariants plotted in Fig. 10a), under simple shear. Note that for i = 2
and i = 4, σ11 is equal to zero.

4.3.3 Biaxial tension In this case the body is exposed to biaxial tensile displacements, see Fig. 12. We consider two fiber families which are orientated symmetric regarding
the X-axis. Further, as indicated in the figure three different stretch ratios will be compared. The stretch ratio is here defined as the quotient F11 /F22 . For the three different
aniso
aniso
stretch ratios the evolution of the quotient σ22
/σ11
is plotted in Fig. 13, 14 and 15.
In each of the figures we find the results for the anisotropic invariants of C i on the left
hand side and the results for (log U )i on the right hand side. For the equi-biaxial test in
Fig. 13 the stress ratio is the same for each model, at each time. The stress ratio directly
aniso
aniso
= tan2 30◦ = 1/3.
/σ11
follows from the fiber angle with respect to the x-axis to be σ22
1
This ratio remains exactly the same for the other stretch ratios only if I4H is used, which
seems to be unreasonable. Regarding the invariants of C we see that the starting point
aniso
aniso
= 1/3, which is different for
/σ11
at nearly zero deformation is always defined by σ22
the Hencky-type strain measures. Further, the change in the stress ratio with increasing
deformation is less pronounced in the case that Hencky strains are used. This behavior can
be explained, when taking a look on Eq. (64). Since the angles θl | l = 1, 2, 3 are constant
only the logarithmic stretches are of interest. Due to the logarithmic function the slope is
decreasing when the stretch is increasing, i.e. for higher strains the slope is smaller than
for lower strains which is also reflected by the shown stress ratios. Generally the slope of
the stress ratios seem to have the opposite sign, regarding the basic strain measure. But
all stresses appear to have the same sign, independent of the stress measure.
1

The stress ratios of the I4C model can be exactly reproduced by the computation of the
fiber angle


hF A, ex i
,
(60)
βact = arccos
kF Ak kex k
aniso
aniso
where ex denotes the direction of the X-axis. Then σ22
/σ11
is equal to tan2 (βact ).

Schröder, von Hoegen, Neff

23

Z

A(2)
A

(1)

60◦

ref.
X

Y

1.35:1.7
1.7:1.35

1.7:1.7

0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

I4C

aniso
aniso
/σ11
σ22

aniso
aniso
/σ11
σ22

Figure 12: Biaxial tension test with two fiber families. The body is assumed to be incompressible with F33 = 1/(F11 F22 )

i

i=1
i=2
i=3
i=4
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

a)

0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

i

I4H

i=1
i=2
i=3
i=4
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

b)

a)

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

i
I4C

i=1
i=2
i=3
i=4
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

aniso
aniso
σ22
/σ11

aniso
aniso
σ22
/σ11

Figure 13: Evolution of the stress ratio for the equi-biaxial test with the stretch ratio
i
i
1.7 : 1.7 for a) the invariants I4C and b) the Hencky-type invariants I4H with the exponent
i = 1, 2, 3, 4.

b)

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

i

I4H

i=1
i=2
i=3
i=4

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

Figure 14: Evolution of the stress ratio for the biaxial test with the stretch ratio 1.7 : 1.35 for
i
i
a) the invariants I4C and b) the Hencky-type invariants I4H with the exponent i = 1, 2, 3, 4.

a)

3.5
3
2.5
2
1.5
1
0.5
0

i=1
i=2
i=3
i=4

aniso
aniso
σ22
/σ11

aniso
aniso
σ22
/σ11

The exponentiated Hencky energy: Anisotropic extension

i
I4C

1 1.05 1.1 1.15 1.2 1.25 1.3 1.35
F11

b)

3.5
3
2.5
2
1.5
1
0.5
0

24

i=1
i=2
i=3
i=4
i

I4H

1 1.05 1.1 1.15 1.2 1.25 1.3 1.35
F11

Figure 15: Evolution of the stress ratio for the biaxial test with the stretch ratio 1.35 : 1.7 for
i
i
a) the invariants I4C and b) the Hencky-type invariants I4H with the exponent i = 1, 2, 3, 4.

4.3.4 Biaxial-tension-compression In order to complete the study of classical mechanical loading scenarios we consider a biaxial combined tension and compression test,
assuming incompressibility. The problem is depicted in Fig. 16 and the fiber arrangement
is identical to the previous example in section 4.3.3. While the stretches λx and λy are
displacement driven, the principal stretch results form the incompressibility condition
λx λy λz = 1.
The fiber invariants, see Fig. 17, are identical for both fiber families and the evolution
is comparable to the uniaxial case. In Fig. 18, again the evolution of the stress ratio is
i
plotted, for both invariant sets. They appear to be quite different. For I4C the stress ratio
1
is approaching zero for infinite strains. For I4H we obtain a constant line as it was already
i
the case in the biaxial tension test. If i in I4H is an even number σ22 becomes negative,
thus the stress ratio becomes negative.
Z

A(2)
A(1)

60◦

ref.
X

Y

1.7 : 0.765
Figure 16: Biaxial combined compression and tension test, where the preferred direction and the loading are aligned. The body is assumed to be incompressible with F33 =
1/(F11 F22 ) and F33 = F22 .

14
12
10
8
6
4
2
0

25

0.25

i

i=1
i=2
i=3
i=4
I4H

I4C

i

Schröder, von Hoegen, Neff

I4C

i=1
0.2 i = 2
i=3
0.15 i = 4
0.1

I4H

i

i

0.05
0
-0.05

1
a)

1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

b)
i

i

Figure 17: Evolution of the anisotropic invariants a) I4C and b) I4H for the exponent
i = 1, 2, 3, 4 under biaxial-tension-compression.

aniso
aniso
/σ11
σ22

0.2
0.1
0
I4C

0.3
aniso
aniso
/σ11
σ22

i=1
i=2
i=3
i=4

0.3

i

0.2
0.1

i=1
i=2
i=3
i=4

I4H

i

0
-0.1

-0.1

-0.2
a)

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

b)

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
F11

Figure 18: Evolution of the stress ratio for the biaxial-tension-compression test for a) the
i
i
invariants I4C and b) the Hencky-type invariants I4H with the exponent i = 1, 2, 3, 4.

4.4

Adaption to biological soft tissues

Soft biological tissues as they occur in arterial walls have an anisotropic material behavior. In biomechanical applications we often assume that the material behaves nearlyincompressible.
The widely used anisotropic strain energy introduced in Holzapfel et al. [30] is given by



 
2 








 µ1 exp k hC, M i −1  − 1
if hC, M i = kF Ak2 ≥ 1




1
ti
|
{z
}

ψHGO
= 2k1 
(61)



C1

I
4


0
if hC, M i = kF Ak2 < 1 ,
1

with the material parameters µ1 and k1 . Recall from Eq. (26), that I4C = hC, M i. The
exponential function incorporated in the latter equation captures the material stiffening
in the high strain domain, caused by the fiber elongation. The case distinction for the
quadratic fiber elongation hC, M i = kF Ak2 in Eq. (61) aims to prevent the fibers from
inducing stiffness under compression. For further use below it is possible to rewrite the

The exponentiated Hencky energy: Anisotropic extension

26

switching criterion for which the anisotropic fiber contribution is neglected, as
hC − 1, M i = hC, M i − 1 < 0 .

(62)

When using anisotropic logarithmic invariants it seems natural to use criteria to switch off
the compression regime which are themselves defined in terms of logarithmic invariants.
However, as will be shown in the following, the case distinction for different anisotropic
invariants (introduced in Eq. (26)) will lead to considerable differences.
Let us consider the right Cauchy-Green tensor
C=

3
X

λ̂k Nk ⊗ Nk ,

with

k=1

1
and
λ̂1 = 0.9 ,
λ̂2 = 1.65 ,
λ̂3 =
λ̂1 λ̂2
 
 
 
1
0
0





N1 = 0 ,
N2 = 1 ,
N3 = 0 ,
0
0
1
represented in the spectral decomposition. The considered deformation is incompressible,
i.e. det C = 1. The set of all possible preferred directions A may be expressed with help
of the spherical coordinates (r, φ, θ) and
x = r sin θ cos φ ,
y = r sin θ sin φ ,
z = r cos φ .

(63)

Here, r denotes the radius, φ ∈ [−π, π] the polar angle and θ ∈ [0, π] the azimuthal angle.
Since kAk = 1, we choose r = 1 and consequently


sin θ cos φ
A =  sin θ sin φ  .
cos φ
Then the anisotropic invariants can be computed from
bi cos2 θ1 + λ
bi cos2 θ2 + λ
bi cos2 θ3
hC i , M i = λ
and
1
2
3

i

i

i
1
1
1
i
2
2
b1 cos θ1 +
b2 cos θ2 +
b3 cos2 θ3 ,
h(log U ) , M i =
log λ
log λ
log λ
2
2
2

(64)

where
cos θ1 = hA, N1 i

cos θ2 = hA, N2 i

cos θ3 = hA, N3 i .

The invariants, depending on the fiber orientation A(x, y, z), are plotted in Fig. 19. While
the principal distributions are similar for different exponents of C, it appears that the
distributions are very different for even and odd exponents i for the logarithmic invariants,
depending on log U . Please note that a change in the eigenvectors Nk | k = 1, 2, 3 would
merely lead to a rotation of the plotted sphere around the eigenvector base. Therefore,
in this scheme the eigenvalues remain as the only predefined variables. As our main goal

Schröder, von Hoegen, Neff

27

is to exclude the compression state from the anisotropic material response to induce
any stiffening we consider the sign of the invariants h(log U )i , M i and hC i , M i − 1 as
the determining criterion. The corresponding distributions are plotted in Fig. 20 over
the azimuthal and polar angle, which are sufficient in order to uniquely define the fiber
orientation. In that sense the black area representing negative values labels the fiber
directions for which the fiber response will be switched off. On the other hand the red
areas of positive values cover the fiber directions for which the anisotropic strain energy
function is switched on. For the invariants depending on C i the area of positive values will
increase with the exponent i and become more elliptic. The plots for sign(hC, M i − 1)
and sign(hlog U , M i) are generally similar. However, with help of Eq. (64) it appears that

hC, M i − 1 − hlog C, M i = hC − 1, M i − hlog C, M i
=

3
X
k

b − 1 − log λ
bk ) cos2 θk ≥ 0
(λ
{z
} | {z }
|k
≥0

≥0

and therefore
hC, M i − 1 ≥ hlog C, M i .

(65)

Note that cos2 θj is the same for each invariant, because the eigenvectors Nj for each
considered strain measure are the same. It immediately follows that

1
[hC, M i − 1] ≥ hlog C, M i = hlog U , M i for hC, M i ≥ 1 .
2

(66)

That means there exists a transition zone, where hC, M i−1 > 0, but h(log U ), M i < 0. In
other words, one may conclude that it is possible that the fiber direction may be stretched
and still the criterion h(log U ), M i will switch off the anisotropic response. This effect
may be favorably used when it is assumed that initially crimped fibers don’t exhibit
significant stiffness until they are straightened out. However, whenever the fiber direction
is compressed, anisotropic material response is precluded for the discussed logarithmic
transversely isotropic strain measure based on the criterion hlog U , M i < 0. Thus, there
is no anisotropic stiffening under compression.
For h(log U )2 , M i and h(log U )4 , M i (and also any other even exponent) the invariants
will always have a positive sign, since both, the structural tensor M as well as (log U )2
and (log U )4 are positive semidefinite. Consequently, these invariants may not serve as a
switching criterion.
In the light of the previous considerations, we may formulate the generalized strain energy

The exponentiated Hencky energy: Anisotropic extension

z

z

y

x
z

N2

y
z

hC 2 , M i − 1

N3

N1

z
hC 3 , M i − 1

N2

z

N2

N1

hC 4 , M i − 1

N2

x

h(log U , M i

x

y
z

N1

x
N3

x
N3

y

N2

N1

y

N3

y

x
N3

x

N1

N2

N1

h(log U )2 , M i

y

N3

h(log U )3 , M i

N2

N1

hC, M i − 1

N3

N3

N2

N1

y

x

Figure 19: Contour plot of different transversely isotropic invariants for a specific choice
of eigenvalues and eigenvectors. The coordinates (x, y, z) define the preferred direction A.

h(log U )4 , M i

z

28

φ

φ

θ
sign(hC 3 , M i − 1)

θ

φ

φ

θ
sign(hC 4 , M i − 1)

θ

φ

θ

sign(h(log U )3 , M i)

φ

sign(h(log U )2 , M i)

θ
sign(hC 2 , M i − 1)

θ

sign(h(log U )4 , M i)

φ

29

sign(h(log U ), M i)

sign(hC, M i − 1)

Schröder, von Hoegen, Neff

φ

θ

Figure 20: Plot of the sign of the invariants plotted in Fig. 19 over the spherical coordinates.
The boundary between positive and negative values marks the transition zone of a possible
case distinction.

The exponentiated Hencky energy: Anisotropic extension

30

functions





 µ1
ti
ψ¹ = 2k1
C




0





 µ1
ti
ψº = 2k1
H




0










i
2
exp k1 (hC , M i −1)  − 1
| {z }




i




if hC i , M i ≥ 1

(67)

I4C

if hC i , M i < 1 ,








2

ε
i
if hlog U , M i ≥ 0
exp k1 hlog U , M i h(log U ) , M i  − 1
{z
} |
{z
}
|




1
i
I4H

I4H

if hlog U , M i < 0 ,
(68)
i

i

in terms of the anisotropic invariants I4C |i = 1, 2, 3, 4 and I4H |i = 1, 2, 3, 4, defined in
Eq. (26). In Eq. (68) the value of ε is to be chosen positive and close to zero. The
incorporation of hlog U , M iε in the free energy function prevents jumps in the stresses at
the switch over point. Alternatively one may also incorporate the criterion (hC, M i − 1)ε
and formulate the case distinction based on the quadratic fiber stretch hC, M i. As it is
shown below a continuous material tangent is also ensured.
Case distinctions, like in Eq. (67) and Eq. (68) may generally lead to discontinuous functions. In order to avoid this it is evident that the stresses and the tangent must become
zero at each switch-point of the chosen criterion. Indeed one can show that
ti
∂ψ¹

∂C
ti
∂ψº

C

=0

if

hC i , M i = 1 ,

H

=0

if

hlog U , M i = 0 ,

∂C

(69)

ti
such that no jumps in the stresses at the switchover points are possible. Only ψ¹
of the
C
introduced strain energy calsses misses continuity, since
ti
∂ 2 ψ¹

C

∂C∂C
ti
∂ 2 ψº
H

∂C∂C

6= 0

if

hC i , M i = 1 ,

=0

if

hlog U , M i = 0 ,

(70)

ti
The material tangent of the function ψ¹
at the point C = 1 in the reference configuration
C
with a structural tensor M = diag(0, 0, 1),


0 0
0
0 0 0
0 0
0
0 0 0


0 0 4i2 µ1 0 0 0
ti
Ci
,
C|C=1 (ψ¹ (I4 |i = 1, 2, 3, 4)) = 
0 0
C
0
0 0 0


0 0
0
0 0 0
0 0
0
0 0 0

is already different from 0 and therefore violates the continuity requirement. For the class
ti
ψº
we find that C|C=1 = 0.
H

Schröder, von Hoegen, Neff

31

5 . Parameter adjustment
5.1. Parameter identification for soft biological tissues
In the following the proposed transversely isotropic Hencky models are adjusted to the
test data provided in Holzapfel [29]. There, a human, abdominal aorta from a human
cadaver was tested. The donor was female, 80 years old and suffering from congestive
cardiomyopathy. Arterial stripes were excised for two material layers, as indicated in
Fig. 21a). Therefore, we introduce the orthonormal coordinate system depending on the
circumferential direction Nϕ , the axial direction Nz and the radial direction Nr . For each
of the layers the tissue was stretched in either circumferential (see Fig. 21b)) or axial (see
Fig. 21c)) direction.
In an incompressible uniaxial tension test with two fiber families orientated in the Nϕ −Nz
plane, we may write the tensors
q



b
 λ1 q0

b2
F =U = 0
λ

0



log U = 

0
0
q

0

1
b1 )
log(λ
2

0
0

1
b1 λ
b2
λ

0
1
b2 )
log(λ
2
0



,




b1 0
λ

b2
C =0 λ
0 0


0

0 

and

1
b1 λ
b2
λ


0

0
.
1
1
log( λb λb )
2
1 2

Because of the incompressibility we have tr(log U ) = log(det U ) = 0 and the second part
of WH and WeH becomes automatically zero. Moreover, the structural tensors are
 2

 2

c
−cs 0
c cs 0
M(1) = −cs s2 0 and M(2) = cs s2 0
(71)
0
0 0
0 0 0
with c = cosβf and s = sinβf and therefore we have
i
I4H

I4C

i

i

i
1
1
2
b
b
log(λ1 ) (cos βf ) +
log(λ2 ) (sin βf )2
=
2
2
bi (cos βf )2 + λ
bi (sin βf )2 .
=λ


1

and

2

The angle βf denotes the angle between each fiber and the local circumferential direction,
while the angle between both fibers follows to 2βf . The second Piola-Kirchoff stresses in
this case may be written as
S11 = 2

∂ψ
p
+2 ,
b1
b1
∂λ
λ

S22 = 2

∂ψ
p
+2 ,
b2
b2
∂λ
λ

S33 = 2

∂ψ
∂ λb 1λb
1 2

+2

p
.
b1 λ
b2
λ

The Lagrange multiplier p is introduced in order to enforce the incompressibility and can
directly be calculated with help of the requirement that S22 must be equal to zero. The
b11 is known from the experiments and the remaining unknown λ
b22 is iterated
stretch λ

The exponentiated Hencky energy: Anisotropic extension

32

with help of Newton’s method, making use of the requirement that also S33 must be equal
to zero:
b2 ) =! 0
S33 (λ

bn −
bn+1 = λ
⇒λ
2
2

bn )
S33 (λ
2
,
bn )
LinS33 (λ
2

bn
bn ) = ∂S33 (λ2 ) .
with LinS33 (λ
2
bn
∂λ
2

(72)

The parameter fitting was performed with help of a Sequential Quadratic Programming
(SQP) algorithm for nonlinear numerical constrained optimization problems. The gradient
needed for the optimization procedure is calculated based on a finite difference scheme in
conjunction with the above described Newton iteration. The objective function
v
!2
nexp u
nmp
exp bi
u 1 X
sim (λ
X
bi , α)
σ
(
λ
)
−
σ
e
t
k
k
f obj (α) :=
(73)
exp
n
max(σ
)
e
mp i=1
e=1
is utilized as the optimization criterion. Here nexp and nmp denote the number of experiments to be fitted and the number of specific measuring points to be evaluated. The
bi . The
predefined amount of stretch associated to each measuring point i is labeled with λ
k
sim
obj
simulated Cauchy stresses σ
and the error f
are dependent on the chosen material
parameter set contained in the field α.
Nr
A(2) Nϕ
A(1)

Nz

2β Med

A(2)
A(1)

Nr

A(2)

f

A(1)

Nz

Nϕ

ref

act.

2β Adv
f

b)

A(2)

A(1)

Nr
Nϕ

Nz

ref
act.

Nz
Adven
titia

a)

ϕ
Media

Nr
c)

Figure 21: a) Excised tissue samples, b) uniaxial tension test in circumferential and c) in
axial direction.

The values of the shear modulus µ for the isotropic response will be directly estimated
from the experiments as the average of the circumferential and axial stiffness. This is
possible since the initial tangent C|C=1 of all adjusted strain energies is solely determined
by µ and will be independent of the anisotropic response. Then considering that
1
∆S = C|C=1 : ∆C ,
2
and C being a function only of µ in the reference configuration, the value of µ can be
directly computed for the first measuring point C(λ̂1k ). This is done separately for both

Schröder, von Hoegen, Neff

33

directions and the mean value is taken as the estimated parameter. However, µ could
as well be considered as a parameter to be included in the parameter set α in Eq. (73),
which is to be optimized. Due to the assumed incompressibility it is impossible to estimate
the bulk modulus κ. From the numerical point of view the bulk modulus may act like a
penalty constraint to enforce quasi incompressibility in more complex computations.

5.2

Media

The results of the parameter adjustment of the Media are given in Table 3 for 8 different
models and the corresponding stress-strain curves are plotted in Fig. 22. These models
result from the combination of the two isotropic strain energy functions in Eq. (29) and
Eq. (35) and the proposed transversely isotropic functions in Eq. (68). The exponent ε
for the switchover criterion is set to 0.1 for all functions. The fiber angle βf between
the fiber direction and the circumferential direction was part of the optimization, while
the parameters κ and k̂ were excluded from the optimization, due to the above assumed
quasi-incompressibility. The isotropic shear modulus was directly estimated from the experiments and was also not optimized.
Apparently, model 1 and 5 fail to accurately fit the experimental data which is evident in
the large values of the objective function f obj . When considering the exponentiated Hencky
energy WeH instead of the classical isotropic energy WH the fit quality is substantially
improved for models 2 and 6. Nevertheless, one has to admit that the associated fiber angle
becomes unsatisfactory small, which does not seem to be reasonable. The small fiber angle
is also accompanied by the effect that the isotropic material response in this case is highly
nonlinear. The functions which are based on invariants with even exponents generally
seem to perform better, independent of the chosen isotropic strain energy function. In
Holzapfel [29], a mean angle of 37.5◦ was reported. The optimized fiber angles of model 3,
4, 7 and 8 are of this order of magnitude. When WeH is used the error generally becomes
slightly smaller, since one additional parameter is available.

1
2
3
4
5
6
7
8

ψ in P2
1
ti
WH + a=1 ψ¹
(I4H(a) )
H
P
1
ti
WeH + 2a=1 ψ¹
(I4H(a) )
H
P
2
ti
WH + 2a=1 ψ¹
(I4H(a) )
H
P
2
ti
WeH + 2a=1 ψ¹
(I4H(a) )
H
P
3
ti
WH + 2a=1 ψ¹
(I4H(a) )
H
P
3
ti
WeH + 2a=1 ψ¹
(I4H(a) )
H
P
4
ti
WH + 2a=1 ψ¹
(I4H(a) )
H
P
4
ti
WeH + 2a=1 ψ¹
(I4H(a) )

µ in kPa
31.16

k in -

µ1 in kPa
0.0001

k1 in 948.81

βf in ◦
25.36

f obj
0.426

31.16

10.54

0.50

107.94

0.73

0.071

31.16

-

1204.86

1599.53

41.24

0.046

31.16

3.38

726.09

1848.66

40.68

0.044

31.16

-

11677.63

3112.51

0.0001

0.386

31.16

10.54

5033.61

17685.18

28.74

0.071

31.16

-

591428.36

51778.23

38.49

0.100

31.16

7.56

232287.68 174224.46

36.86

0.052

H

Table 3: Adjusted parameter sets of the Media.

circumf
axial.
model 1
model 2

1 1.05 1.1 1.15p
1.2 1.25 1.3 1.35

Cauchy sress σ in kPa

Stretch

140
120
100
80
60
40
20
0

140
120
100
80
60
40
20
0

1 1.05 1.1 1.15p
1.2 1.25 1.3 1.35

circumf
axial.
model 3
model 4

Stretch

b in λ

34

1 1.05 1.1 1.15p
1.2 1.25 1.3 1.35

b in λ

circumf
axial.
model 5
model 6

Stretch

Cauchy sress σ in kPa

140
120
100
80
60
40
20
0

Cauchy sress σ in kPa

Cauchy sress σ in kPa

The exponentiated Hencky energy: Anisotropic extension

140
120
100
80
60
40
20
0

b in λ

circumf
axial.
model 7
model 8

1 1.05 1.1 1.15p
1.2 1.25 1.3 1.35

Stretch

b in λ

Figure 22: Adjusted stress-strain curves for the media for model 1-8.

5.3

Adventitia

The case distinction in the isotropic response for the parameter fitting is neglected for the
Adventitia. Only four models are tested in the following. The optimized parameters are
shown in Table 4 and the corresponding stress-strain curves are given in Fig. 23. Please
note that in this case the constraint βf ≤ 75◦ was added to the optimization procedure,
which is in the range of experimentally obtained mean angles, see Holzapfel [29]. Again,
only model 4 and 8 are able to reproduce the experimental curves.

2 WeH +
4 WeH +
6 WeH +
8 WeH +

ψ in P2
ti
H1
a=1 ψ¹H (I4(a) )
P2
ti
H2
a=1 ψ¹H (I4(a) )
P2
ti
H3
a=1 ψ¹H (I4(a) )
P2
ti
H4
a=1 ψ¹ (I4(a) )

µ in kPa
15.90

k in 3.50

µ1 in kPa
7.4e-06

k1 in 246.90

βf in ◦
75.00

f obj
0.127

15.90

1.0e-08

0.05

3707.99

50.30

0.045

15.90

3.58

6.00

17045.18

63.28

0.128

15.90

1.0e-08

3973.68

56653.51

54.74

0.071

H

Table 4: Adjusted parameter sets of the Adventitia.

100

Cauchy sress σ in kPa

Cauchy sress σ in kPa

Schröder, von Hoegen, Neff

circumf
axial.
model 2
model 4

80
60
40
20
0
1

1.1

1.2 p1.3

Stretch

1.4

1.5

b in λ

100

35

circumf
axial.
model 6
model 8

80
60
40
20
0
1

1.1

1.2 p1.3

Stretch

1.4

1.5

b in λ

Figure 23: Evolution of the stress ratio induced by different anisotropic measures under
biaxial-tension-compression.

5.4

Artery

In this section the calibrated transversely isotropic model is applied to three-dimensional
simulations of coronary patient-specific arterial walls in order to prove the robustness of
the novel material formulation. Since the parameter fitting was based on an diseased abdominal artery, the results may merely be discussed in a qualitative context and conclusion
related to clinical interpretation are highly limited. An arterial geometry is reconstructed
based on two sequenced two-dimensional virtual histology (VH) intravascular ultrasound
(IVUS) images. For a detailed description of the three-dimensional reconstruction the
reader is referred to Balzani et al. [7].
The considered artery consists of two layers, see Fig. 24. The outer layer is the
Adventitia, the inner layer the Media. This artery was loaded in a finite element simulation with an inner pressure of 16 kPa. An augmented Lagrange strategy was applied to enforce quasi-incompressibility with an allowed tolerance of
1% in change of volume. For additional information the reader is referred to
Hestenes [24], Powell [51], Glowinski and Le Tallec [18], Glowinski and Le Tallec [19] and
Glowinski and Le Tallec [20]. Tetrahedron finite elements with ten nodes and quadratic
shape functions are used.
2

In a first simulation model 4, based on I4H , was used with the estimated parameters from the adjustment, see Table 3 and Table 4. These results are to be compared
with a second simulation, where the strain energy function and parameters according to
Schröder and Brinkhues [53] were used. There a Mooney-Rivlin model was used for the
isotropic part and the function from Holzapfel et al. [30], given in Eq. (61), for the superimposed transversely isotropic parts. Both parts refer to the strain measure C. The
parameters were adjusted to the same experimental data. The deformed configurations
for different stresses are shown in Fig. 25. While the general stress distributions are comparable the increase in the volume of the lumen is significantly larger for the Hencky
model.
The same effect for different strain energy functions, which were adjusted to the same
data and then used for numerical simulations of arterial segments was noticed in
Brands et al. [10].

The exponentiated Hencky energy: Anisotropic extension

Y

Y
X

Z

36

X

Y

Z

X

Z

Figure 24: Reconstructed artery consisting of Adventitia (outer layer) and Media (inner
layer).

Y
X

Z

Y

σ11

X

Z

Y

σ22

200
180
160
140
120
100
80
60
40
20
0

X

σ33

Z

80
70
60
50
40
30
20
10
0

300
270
240
210
180
150
120
90
60
30
0

a)
Y
X

Z

Y

σ11
200
180
160
140
120
100
80
60
40
20
0

b)

X

Z

Y

σ22

X

Z

σ33
80
70
60
50
40
30
20
10
0

300
270
240
210
180
150
120
90
60
30
0

in kPa
Figure 25: Deformed configuration and stress distributions in an arterial segment, a) sim2
ulated with model 4, based on I4H and b) the material model and parameters used in
1
Schröder and Brinkhues [53] based on I4C for an internal pressure p of 16 kPa.

6

Open problems

While the ellipticity domain of the isotropic exponentiated Hencky energy is reasonably
well understood (it is an extremely large cone in stretch space) the same is not clear for
the anisotropic logarithmic energy terms. The problem which has to be faced is due to
the analytical difficulties in treating derivatives of the right Hencky strain tensor; we
need to calculate for example ∂F2 (exphlog U , M i2 ).(ξ ⊗ η, ξ ⊗ η) which is manageable
along the lines of the presented algorithmic treatment in this paper; but it remains to
identify a suitable ellipticity region. The experience with the isotropic exponentiated
Hencky energy suggests, however, that the ellipticity domain largely contains the physical
range of arteries, i.e. principal stretches in a suitable λk ∈ [1/2, 2]. Therefore, our model
proposal is elliptic in the physiological range of arteries and this is all that must be
required on mathematical grounds.
Since our extension of the exponentiated Hencky energy towards anisotropy is done on an

Schröder, von Hoegen, Neff

37

ad hoc basis, it remains to study whether the differential-geometric program presented
in Neff et al. [48] may be extended to the anisotropic case. Of course, major technical
difficulties have to be solved. The benefit, however, cannot be overemphasized: There
would result finite strain anisotropic energies having a clear physical meaning. We will
pursue this issue in the future.
Acknowledgments The first two authors gratefully acknowledge support by the
Deutsche Forschungsgemeinschaft in the Priority Program 1748 under the project “Novel
finite elements for anisotropic media at finite strain” (SCHR 570/23-1), section 1-4.3. Further, the first two authors would like to acknowledge support by the Deutsche Forschungsgemeinschaft within the framework of the project “Domain-decomposition-based fluid
structure interaction algorithms for highly nonlinear and anisotropic elastic arterial wall
models in 3D“ (SCHR 570/15-2) under the D-A-CH agreement, section 4.4-5.

The exponentiated Hencky energy: Anisotropic extension

38

References
[1] L. Anand. On H. Henckys approximate strain energy function for moderate deformations. Journal of Applied Mechanics, 46:78–82, 1979.
[2] L. Anand. Moderate deformations in extension - torsion of incompressible isotropic
elastic materials. Journal of the Mechanics and Physics of Solids, 34:293 304, 1986.
[3] M. Baker and J. L. Ericksen. Inequalities restricting the form of the stressdeformation relations for isotropic elastic solids and Reiner-Rivlin fluids. Journal
of the Washington Academy of Sciences, 44:33–35, 1954.
[4] J. M. Ball. Convexity conditions and existence theorems in non-linear elasticity.
Archive of Rational Mechanics and Analysis, 63:337–403, 1977.
[5] J. M. Ball. Some open problems in elasticity. In P. Newton, P. Holmes, and A. Weinstein, editors, Geometry, mechanics, and dynamics, pages 3–59. Springer, New-York,
2002.
[6] D. Balzani, P. Neff, J. Schröder, and G. A. Holzapfel. A polyconvex framework for
soft biological tissues. Adjustment to experimental data. International Journal of
Solids and Structures, 43:6052–6070, 2006.
[7] D. Balzani, D. Böse, D. Brands, R. Erbel, A. Klawonn, O. Rheinbach, and J. Schröder.
Parallel simulation of patient-specific atherosclerotic arteries for the enhancement of
intravascular ultrasound diagnostics. Engineering Computations, 29:888–906, 2012.
[8] J. P. Boehler. Lois de comportement anisotrope des milieux continus. Journal de
Mécanique, 17:153–190, 1978.
[9] J. P. Boehler. A simple derivation of representations for non-polynomial constitutive
equations in some cases of anisotropy. Zeitschrift für angewandte Mathematik und
Mechanik, 59:157–167, 1979.
[10] D. Brands, A. Klawonn, O. Rheinbach, and J. Schröder. Modelling and convergence in
arterial wall simulations using a parallel FETI solution strategy. Computer Methods
in Biomechanics and Biomedical Engineering, 11:569–583, 2008.
[11] O. T. Bruhns, H. Xiao, and A. Meyers. Hencky’s elasticity model with the logarithmic
strain measure: A study on Poynting effect and stress response in torsion of tubes
and rods. Archives of Mechanics, 52:489–509, 2000.
[12] O. T. Bruhns, H. Xiao, and A. Meyers. Constitutive inequalities for an isotropic elastic strain-energy function based on Hencky’s logarithmic strain tensor. Proceedings
of the Royal Society of London, Series A : Mathematical, Physical and Engineering
Science, 457:2207–2226, 2001.
[13] V. Ebbing, J. Schröder, and P. Neff. Approximation of anisotropic elasticity tensors
at the reference state with polyconvex energies. Archive of Applied Mechanics, 79:
651–657, 2009.

Schröder, von Hoegen, Neff

39

[14] A. E. Ehret and M. Itskov. A polyconvex hyperelastic model for fiber-reinforced
materials in application to soft tissues. Journal of Material Science, 42:8853–8863,
2007.
[15] P. J. Flory. Thermodynamic relations for high elastic materials. Trans. Faraday Soc.,
57:829–838, 1961.
[16] I. D. Ghiba, P. Neff, and R. Martin. An ellipticity domain for the distortional Henckylogarithmic strain energy. Proceedings of the Royal Society of London, Series A :
Mathematical, Physical and Engineering Science, 471:20150510, 2015.
[17] I. D. Ghiba, P. Neff, and M. Šilhavỳ. The exponentiated Hencky-logarithmic strain
energy. Improvement of planar polyconvexity. International Journal of Non-Linear
Mechanics, 71:48–51, 2015.
[18] R. Glowinski and P. Le Tallec. Finite element analysis in nonlinear incompressible
elasticity. In J. Oden and G. Carey, editors, Finite elements, Vol V: Special Problems
in Solid Mechanics. Prentice-Hall, Englewood Cliffs, NH, 1984.
[19] R. Glowinski and P. Le Tallec. Augmented Lagrangian methods for the solution of
variational problems. Springer, Berlin, 1988.
[20] R. Glowinski and P. Le Tallec. Augmented Lagrangian and operator-splitting methods
in nonlinear mechanics, volume 9 of SIAM Studies in Applied Mathematics. Society
for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1989.
[21] S. Hartmann and P. Neff. Polyconvexity of generalized polynomial-type hyperelastic
strain energy functions for near-incompressibility. International Journal of Solids
and Structures, 40:2767–2791, 2003.
[22] H. Hencky. Über die Form des Elastizitätsgesetzes bei ideal elastischen Stoffen.
Zeitschrift für technische Physik, 9:215–220, 1928.
[23] H. Hencky. Welche Umstände bedingen die Verfestigung bei der bildsamen Verformung von festen isotropen Körpern? Zeitschrift für Physik, 55:145–155, 1929.
[24] M. Hestenes. Multiplier and gradient methods. Journal of Optimization Theory and
Applications, 4:303–320, 1969.
[25] R. Hill. On constitutive inequalities for simple materials. Journal of the Mechanics
and Physics of Solids, 16:229–242, 1968.
[26] R. Hill. Constitutive inequalities for isotropic elastic solids under finite strain. Proceedings of the Royal Society of London, Series A : Mathematical, Physical and Engineering Science, 314:457–472, 1970.
[27] R. Hill. Aspects of invariance in solid mechanics. Advances in Applied Mechanics,
18:1–75, 1978.
[28] A. Hoger. The stress conjugate to logarithmic strain. International Journal of Solids
and Structures, 23:1645–1656, 1987.

The exponentiated Hencky energy: Anisotropic extension

40

[29] G. A. Holzapfel. Determination of material models for arterial walls from uniaxial
extension tests and histological structure. Journal of Theoretical Biology, 238:290–
302, 2006.
[30] G. A. Holzapfel, T. C. Gasser, and R. W. Ogden. A new constitutive framework
for arterial wall mechanics and a comparative study of material models. Journal of
Elasticity, 61:1–48, 2000.
[31] S. K. Hyun and H. Nakajima. Anisotropic compressive properties of porous copper
produced by unidirectional solidification. Material Science and Engineering, A340:
258–264, 2003.
[32] M. Itskov, A. E. Ehret, and D. Mavrilas. A polyconvex anisotropic strain-energy
function for soft collagenous tissues. Biomechanics and Modelling in Mechanobiology,
5:17–26, 2006.
[33] C. S. Jog. Derivatives of the stretch, rotation and exponential tensors in ndimensional vector spaces. Journal of Elasticity, 82:175–192, 2006.
[34] C. S. Jog and K. D. Patil. Conditions for the onset of elastic and material instabilities
in hyperelastic materials. Archive of Applied Mechanics, 83(5):661–684, 2013.
[35] M. Latorre and F. J. Montáns. Anisotropic finite strain viscoelasticity based on
the Sidoroff multiplicative decomposition and logarithmic strains. Computational
Mechanics, 56:503–531, 2015.
[36] J. Löblein, J. Schröder, and F. Gruttmann. Application of generalized measures to
an orthotropic finite elasto-plasticity model. Computational Materials Science, 28:
696–703, 2003.
[37] R. J. Martin and P. Neff. Minimal geodesics on GL(n) for left-invariant, right-O(n)invariant Riemannian metrics. Journal of Geometric Mechanics, 8(3):323–357, 2016.
[38] L. A. Mihai and P. Neff. Hyperelastic bodies under homogeneous Cauchy stress
induced by three-dimensional non-homogeneous deformations. to appear in Mathematics and Mechanics of Solids, 2017.
[39] L. A. Mihai and P. Neff. Hyperelastic bodies under homogeneous Cauchy stress
induced by non-homogeneous finite deformations. International Journal of NonLinear Mechanics, 89:93–100, 2017.
[40] G. Montella, S. Govindjee, and P. Neff. The exponentiated Hencky strain energy
in modelling tire derived material for moderately large deformations. Journal of
Engineering Materials and Technology, 138:031008, 2016.
[41] P. Neff. Mathematische Analyse multiplikativer Viskoplastizität. Ph.D. Thesis, Technische Universität Darmstadt. Shaker Verlag, Aachen, 2000.
[42] P. Neff and I. D. Ghiba. The exponentiated Hencky-logarithmic strain energy. Part
III: Coupling with idealized isotropic finite strain plasticity. Continuum Mechanics
and Thermodynamics, 28:477–487, 2016.

Schröder, von Hoegen, Neff

41

[43] P. Neff and L. A. Mihai. Injectivity of the Cauchy-stress tensor along rank-one
connected lines under strict rank-one convexity condition. to appear in Journal of
Elasticity, 2016.
[44] P. Neff, B. Eidel, F. Osterbrink, and R. J. Martin. The Hencky strain energy ||log U ||2
measures the geodesic distance of the deformation gradient to SO(3) in the canonical
left-invariant Riemannian metric on GL(3). Proceedings of Applied Mathematics and
Mechanics, 13:369–370, 2013.
[45] P. Neff, B. Eidel, and R. J. Martin. The axiomatic deduction of the quadratic Hencky
strain energy by Heinrich Hencky. arXiv:1402.4027, 2014.
[46] P. Neff, I. Ghiba, and J. Lankeit. The exponentiated Hencky-logarithmic strain
energy. Part I: Constitutive issues and rank-one convexity. Journal of Elasticity, 121:
143–234, 2015.
[47] P. Neff, J. Lankeit, I. D. Ghiba, R. J. Martin, and D. J. Steigmann. The exponentiated Hencky-logarithmic strain energy. Part II: Coercivity, planar polyconvexity
and existence of minimizers. Zeitschrift für angewandte Mathematik und Physik, 66:
1671–1693, 2015.
[48] P. Neff, B. Eidel, and R. J. Martin. Geometry of logarithmic strain measures in solid
mechanics. Archive of Rational Mechanics and Analysis, 222:507–572, 2016.
[49] R. W. Ogden. Large deformation isotropic elasticity - On the correlation of theory
and experiment for incompressible rubberlike solids,. Proceedings of the Royal Society
London A, 326:565584, 1972.
[50] R. W. Ogden. Non-Linear elastic deformations. Courier Corporation, 1997.
[51] M. Powell. A method for nonlinear constraints in minimization problems. In
R. Fletcher, editor, Optimization, pages 283–298. Academic Press, New York, 1969.
[52] H. Richter. Das isotrope Elastizitätsgesetz. Zeitschrift für angewandte Mathematik
und Mechanik, 28:205–209, 1948.
[53] J. Schröder and S. Brinkhues. A novel scheme for the computation of residual stresses
in arterial walls. Archive of Applied Mechanics, 84:881–898, 2014.
[54] J. Schröder and D. Gross. Invariant formulation of the electromechanical enthalpy
function of transversely isotropic piezoelectric materials. Archive of Applied Mechanics, 73:533–552, 2004.
[55] J. Schröder and P. Neff. Invariant formulation of hyperelastic transverse isotropy
based on polyconvex free energy functions. International Journal of Solids and Structures, 40:401–445, 2003.
[56] J. Schröder, F. Gruttmann, and J. Löblein. A simple orthotropic finite elastoplasticity
model based on generalized stress-strain measures. Computational Mechanics, 30:48–
64, 2002.

The exponentiated Hencky energy: Anisotropic extension

42

[57] J. Schröder, P. Neff, and D. Balzani. A variational approach for materially stable
anisotropic hyperelasticity. International Journal of Solids and Structures, 42(15):
4352–4371, 2005.
[58] J. Schröder, P. Neff, and V. Ebbing. Anisotropic polyconvex energies on the basis of
crystallographic motivated structural tensors. Journal of the Mechanics and Physics
of Solids, 56:3486–3506, 2008.
[59] J. C. Simo. Numerical Analysis and Simulation of Plasticity, Vol. 6 of Handbook of
Numerical Analysis. Elsevier Science, 1998.
[60] A. J. M. Spencer. Kinematic constraints, constitutive equations and failure rules for
anisotropic materials. In J. P. Boehler, editor, Applications of Tensor Functions in
Solid Mechanics, volume 292 of CISM Courses and Lectures, pages 187–197. Springer,
1987.
[61] C. Vallée. Lois de comportement élastique isotropes en grandes déformations. International Journal of Engineering Science, 16:451–457, 1978.

Schröder, von Hoegen, Neff

A
A.1

43

Appendix
Notes on the Hencky tensor

The Hencky strain tensor is defined through
log U =

1
log C
2

(74)

and
tr(log U ) = log(detU ) .

(75)

The symmetric right Cauchy-Green tensor C in spectral decomposition is given by

C=

3
X

bk N k ⊗ N k =
λ

k=1

3
X

bk P k ,
λ

with P k = N k ⊗ N k

(76)

k=1

and for the Hencky strain we obtain

log U =

3
X
1
k=1

2

bk )N k ⊗ N k =
log(λ

3
X
1
k=1

2

bk )Pk ,
log(λ

with Pk = N k ⊗ N k .

(77)

The first derivative of log U with respect to C can be computed as
3

bk 1
∂ 12 log λ
∂log U X
bk ∂Pk
=
Pk ⊗
+ log λ
∂C
∂C
2
∂C
k=1
=

=

3
X
k=1
3
X
k=1

Pk ⊗

bk ∂ λ
bk 1
∂ 12 log λ
bk ∂Pk
+ log λ
bk
∂C
2
∂C
∂λ

bk 1
∂λ
1 b−1
bk ∂Pk
λk Pk ⊗
+ log λ
2
∂C
2
∂C

(78)

Considering that
bk
∂λ
= Pk
∂C

3
X
Pk  PjT + Pj  PkT
∂Pk
=
,
bk − λ
bj
∂C
λ
j=1,j6=k

and

(79)

see for instance Jog [33], we find that
3
3 X
3 1
bk ) − 1 (log λ
bj )
X
X
(log λ
∂ log U
1 b−1
2
2
λk Pk ⊗ Pk +
(Pk  Pj + Pj  Pk ) . (80)
=
bk − λ
bj
∂C
2
λ
k=1
k=1 j=1
k6=j

The exponentiated Hencky energy: Anisotropic extension

44

The second derivative for the linearization is given by
"
#
3546
3
b−1 1
T
T
∂ 21 λ
∂ 2 log U X
∂P
∂P
k
k
k
b−1
=
+ λ
⊗ Pk
Pk ⊗ Pk ⊗
+ Pk ⊗
k
∂C∂C
∂C
2
∂C
∂C
k=1
1 b−1
2
1
bk ) ∂ Pk + ∂Pk ⊗ ∂ 2 λk
+ (logλ
2
∂C∂C
∂C
∂C"
#

3546
3
X 1
T T
∂P
1
∂P
k
k
−2
−1
b Pk ⊗ Pk ⊗ Pk + λ
b
=
⊗ Pk
− λ
+ Pk ⊗
k
k
2
2
∂C
∂C
k=1
2
1
b k ) ∂ Pk + 1 λ
b−1 ∂Pk ⊗ Pk ,
+ (logλ
2
∂C∂C 2 k ∂C

(81)

where
2



3
X

∂ Pk
1

T
=
 Pj ⊗
bk − λ
bj
∂C∂C j=1 λ



∂Pk
∂C

12 !23
T
T


23
∂Pj T
+ Pk ⊗
∂C

k6=j

+

+

PkT ⊗



1
bk − λ
bj
λ

∂Pj
∂C

12 !23
T
T



23
∂Pk T 
+ Pj ⊗

∂C

2 (Pk  Pj + Pj  Pk ) ⊗ (Pj − Pk ) .

(82)

The exponent and the logarithm of an arbitrary symmetric tensor may also be be expressed
with help of a Taylor expansion of the form
exp(•) = 1 +
log(•) =

∞
X
1
(•)m
m!
m=1

∞
X
−1k−1
k=1

k

and

[(•) − 1]k ,

(83)
(84)

where the latter is convergent in a neighborhood of 1.
A.2. Conjugate stress tensors
The following considerations are adapted from Ogden [50].
The constitutive equation for the stresses are derived form the (isothermal) entropy inequality
ψ̇ − hP , Ḟ i ≥ 0 ,
(85)
From the latter we deduce the constitutive relation P = ∂F ψ. Let the generalized Lagrangean strain measures
(
1
(U m − 1)
m=
6 0
E (m) = 21
(86)
log U
m=0
2

Schröder, von Hoegen, Neff

45

and Eulerian strain measures
(
K (m) =

1
m

(V m − 1)
logV

m 6= 0
m=0

(87)

be given, we aim to find the corresponding constitutive equations. The so called stress
power may be written as
hP , Ḟ i = hτ , Di = hS, Ėi = h∂E (m) ψ, Ė (m) i = h∂K (m) ψ, K̇ (m) i ,
(88)


where D = 21 (L + LT ) and L = grad ẋ. Considering that Ė = 21 Ḟ T F + F T Ḟ , we
obtain the relations
P = F S = τ F −T .
(89)
The pairs in Eq. (88) are said to be work
 conjugate.
 By making use of the fact that
RT Ṙ = −ṘT R, we may rewrite Ė = 12 U U̇ + U̇ U and we are able to reformulate
D

S,

E
E D 1
1
U U̇ + U̇ U
=
(SU + U S), U̇ ,
2
{z
}
|2

(90)

TBiot

such that we directly obtain the Biot stress TBiot = ∂U ψ # (U ), work conjugate to U from
the entropy inequality. With U = RT V R we are able to relate E (m) and K (m) and the
corresponding time derivatives as follows:
E (m) = RT K (m) R and Ė (m) = RT K̇ (m) R + E (m) RT Ṙ − RT ṘE (m) .

(91)

Regarding the generalized stress-power it follows


h∂E (m) ψ, Ė (m) i = hR(∂E (m) ψ)RT , K̇ (m) i + h (∂E (m) ψ)E (m) − E (m) (∂E (m) ψ) , RT Ṙi (92)
Only if E (m) ∂E (m) ψ = ∂E (m) ψE (m) , i.e. ∂E (m) ψ is coaxial with E (m) , it immediately
follows that the constitutive law results in
∂ψ
= R(∂E (m) ψ)RT
(93)
∂K (m)
and the stress power is expressed through hK̇ (m) , R(∂E (m) ψ)RT i, which is identical to
Eq. (88).
The case that ∂E (m) ψ is coaxial with E (m) implies that also ∂E (m) ψ and U are coaxial.
Under this assumption one may show that
∂E (m) ψ = TBiot U −(m−1) = U −(m−1) TBiot ,

(94)

and it follows that
∂E (m) ψ = TBiot U −m+1

and ∂E (0) ψ = ∂log U ψ = TBiot U = RT τ R ,

(95)

for isotropic materials. Inserting the latter result in Eq. (93) we obtain the relation
τ =

∂ψ(logV )
.
∂logV

(96)

Regarding the conjugate stress to log U the reader is also referred to Hill [26; 27] and
Hoger [28].

