A Bimodular Polyconvex
Anisotropic Strain Energy
Function for Articular Cartilage
Stephen M. Klisch
A strain energy function for finite deformations is developed that has the capability to
Associate Professor
describe the nonlinear, anisotropic, and asymmetric mechanical response that is typical
Mechanical Engineering Department,
of articular cartilage. In particular, the bimodular feature is employed by including strain
California Polytechnic State University,
energy terms that are only mechanically active when the corresponding fiber directions
San Luis Obispo, CA 93407
are in tension. Furthermore, the strain energy function is a polyconvex function of the
e-mail: sklisch@calpoly.edu
deformation gradient tensor so that it meets material stability criteria. A novel feature of
the model is the use of bimodular and polyconvex “strong interaction terms” for the
strain invariants of orthotropic materials. Several regression analyses are performed
using a hypothetical experimental dataset that captures the anisotropic and asymmetric
behavior of articular cartilage. The results suggest that the main advantage of a model
employing the strong interaction terms is to provide the capability for modeling aniso-
tropic and asymmetric Poisson’s ratios, as well as axial stress–axial strain responses, in
tension and compression for finite deformations. DOI: 10.1115/1.2486225
Introduction
The extracellular solid matrix of articular cartilage contains
proteoglycans and a crosslinked collagen network. The proteogly-
cans are negatively charged molecules that primarily resist com-
pressive loads 1,2 while the collagen network primarily resists
tensile and shear loads 3,4. Due in part to its complex molecular
structure, articular cartilage typically behaves as an anisotropic
material with substantial tension-compression asymmetry 5–10
and likely experiences finite, multi-dimensional strains when sub-
ject to typical loads 11,12. In particular, both the Young’s modu-
lus and Poisson’s ratio
1
are anisotropic and strain dependent, and
can be approximately two orders of magnitude greater in tension
than in compression 9,10,13–19. Consequently, the development
of accurate finite deformation models of the equilibrium elastic
response is challenging.
Bimodular elastic and biphasic models have been developed
that can model the asymmetric tensile and compressive mechani-
cal properties for infinitesimal strains 8,10. Those models were
based on a general bimodular theory for infinitesimal strains 20
in which the material constants may be discontinuous or jump
across a surface of discontinuity in strain space, provided that
stress continuity conditions are satisfied at the surface. Several
exponential models for finite deformations allowing for different
mechanical properties in tension and compression have been used
for the arterial wall 21 and the annulus fibrosus 22. However,
those models have not employed a general bimodular theory that
ensures stress continuity across the surface of discontinuity. Re-
cently, a general bimodular theory employing second-order and
exponential stress–strain equations was shown to be capable of
modeling the anisotropy and asymmetry in Young’s modulus for
finite deformations 23. Those results suggested that, when using
the bimodular feature, second-order models might provide a ma-
terial description as accurate as those provided by exponential
models. However, the models studied in Ref. 23 were not ca-
1
In this paper, the terms “Young’s modulus” and “Poisson’s ratio” will be used to
refer to strain-dependent functions because a finite deformation theory is used.
Contributed by the Bioengineering Division of ASME for publication in the J
OUR-
NAL OF
BIOMECHANICAL ENGINEERING. Manuscript received May 1, 2006; final manu-
script received September 15, 2006. Review conducted by Clark T. Hung.
pable of modeling the anisotropy and asymmetry in Poisson’s ra-
tio for finite deformations and were not appropriate for use in
computational solutions.
2
The overall goal of this study is to develop an elastic strain
energy function for finite deformations of the articular cartilage
solid matrix that meets several criteria. First, it should be capable
of modeling the nonlinearity, anisotropy, and asymmetry in
Young’s modulus and Poisson’s ratio. Although the desired accu-
racy of the stress–strain equation may not be the same for all
applications, the level of accuracy sought here is likely to be
crucial in continuum growth analysis.
3
In order to meet this crite-
rion, the bimodular feature is employed. Second, it should satisfy
stability criteria so that numerical stability of computational solu-
tions can be expected. In order to meet this criterion, a polycon-
vex strain energy function is developed; polyconvexity guarantees
the existence of local minimizers of the strain energy function
when subject to boundary conditions 24 while not sharing the
limitations of convexity with respect to the violation of invariance
requirements and global uniqueness. Third, it should use a rela-
tively low number of parameters needed to model the desired
elastic response, so that the material constants are based on a
model that is not over-parameterized
4
and can be determined from
a combination of several common experimental protocols.
In a preliminary study, a bimodular polyconvex strain energy
function was developed for articular cartilage based on the strain
invariants for an orthotropic material 25; however, that model,
nor the earlier second-order and exponential models 23, were
capable of modeling the anisotropy and asymmetry of Poisson’s
ratio. In that polyconvex model 25, there were no strong inter-
action or coupling terms for the orthotropic strain invariants. Al-
though recent studies have proposed 26 or used 27 strong in-
teraction terms for orthotropic strain invariants, preliminary
studies for this work were not successful in using those terms with
the bimodular feature.
5
Also, those studies have not discussed
2
The second-order model in terms of the first Piola-Kirchhoff stress developed in
Ref. 23 was shown to satisfy stability criteria; however, the corresponding Cauchy
stress was not.
3
Applications are presented in the “Discussion.”
4
See the “Discussion” for comments on “over-parameterization” in the context of
the nonlinear regression analysis used here.
5
In particular, strong interaction terms that satisfy the bimodular stress–strain
continuity conditions stated in Eq. 7 were not found for an orthotropic material.
250 / Vol. 129, APRIL 2007 Copyright © 2007 by ASME Transactions of the ASME
possible microstructural causes. An aim of this work is to derive
bimodular strong interaction terms that are “simple” enough to
allow the experimenter to investigate possible microstructural
mechanisms.
A recent development of Ref. 28 used two mechanically
equivalent secondary fiber families, in addition to primary fiber
families, to represent the phenomena of collagen crosslinking in
the annulus fibrosus tissue; that model was capable of producing
tensile Poisson’s ratios that are an order of magnitude greater than
those of our earlier studies 23,25. Since the secondary fiber
families introduced in Ref. 28 basically serve as strong interac-
tion or coupling terms for the strain invariants related to the pri-
mary fiber families, here it was hypothesized that the introduction
of strain invariants generated by secondary fiber families will al-
low a more accurate description of tensile Poisson’s ratios for
articular cartilage. In contrast to Ref. 28, this development is
incorporated into a bimodular polyconvex strain energy function.
The specific objectives are to: 1 adapt the bimodular theory
for finite deformations to the present application; 2 develop a
bimodular polyconvex anisotropic strain energy function using
primary fibers and strong interaction terms generated by second-
ary fibers; and 3 compare the predictive capability of models
with and without the strong interaction terms using experimental
data gathered from the literature. The results suggest that using
both the bimodular feature and the strong interaction terms facili-
tates the accurate description of the anisotropic and asymmetric
mechanical properties of articular cartilage in finite deformations.
Methods
Background. The right Cauchy–Green deformation tensor C is
defined as
C = F
T
F 1
where F is the deformation gradient tensor and the superscript T
signifies the transpose operator. The Cauchy, first Piola–Kirchhoff,
and second Piola–Kirchhoff stress tensors denoted as T, P, and S,
respectively are related by
JT = PF
T
= FSF
T
2
where J is the determinant of F. The stress constitutive equations
for a Green-elastic material may be expressed as
n
W
W
I
i
S =2 =2
3
C
C
i=1
I
i
where W = W
ˆ
I
i
is a scalar strain energy function that depends on
a set of invariants I
i
corresponding to the material symmetry
group. The fourth-order elasticity tensor is defined as
S
C = 4
C
Bimodular Elasticity for Finite Deformations. Due to the ob-
served tension–compression asymmetry of the articular cartilage
solid matrix, a bimodular theory is used. Earlier models 8,10,23
were based on a bimodular theory 20 in which the material
constants may be discontinuous or jump across a surface of dis-
continuity in strain space, provided that stress continuity condi-
tions are satisfied at the surface. The bimodular theory of Ref. 20
was developed in terms of the second Piola–Kirchhoff stress and
Lagrange strain tensors. Here, that theory is reformulated to use C
instead of the Lagrange strain tensor.
A scalar valued function of C identifying a surface of discon-
tinuity in the six-dimensional strain space of C is defined as
gC =0 5
and is restricted to be a function of the invariants corresponding to
the material symmetry group. Different strain energy functions
may be specified on either side of a surface of discontinuity; i.e.
Fig. 1 Schematic of the coordinate system and experimental
specimen orientations in relation to anatomical directions. The
unit vector E
1
is parallel to the local split-line direction, the unit
vector E
3
is perpendicular to the articular surface, and the unit
vector E
2
is perpendicular to the split-line direction and parallel
to the surface. The cylinders labeled P
11
, P
22
, and P
33
represent
specimens loaded in tension or compression along the E
1
,E
2
,
and E
3
directions, respectively.
W = W
+
if gC 0, W = W
−
if gC 0 6
In a similar fashion, different stress and elasticity tensors may be
specified on either side of a surface of discontinuity; i.e., as
S
+
,S
−
, C
+
,C
−
In Ref. 20, a theorem was proved establishing necessary and
sufficient conditions for stress continuity across the surface of
discontinuity. Introducing a slight modification in Lemma 3.2 of
Ref. 20, one obtains the following necessary and sufficient con-
ditions for stress continuity across the surface of discontinuity
g
g
S = S
+
= S
−
, C = C
+
− C
−
= sC �
7
C
C
for all C that satisfy gC= 0, where C represents the jump in
the elasticity tensor, sC is a scalar valued function of C, and
�
is the tensor dyadic product.
Structural anisotropy. Spencer 29 proposed a general theory
capable of modeling an anisotropic material as a composite mate-
rial consisting of an isotropic matrix reinforced with fiber fami-
lies. That theory has been used to develop strain energy functions
for cartilaginous tissues 22,28,30–32. For example, in Ref. 30
two mechanically equivalent fiber families were used to model the
annulus fibrosus in finite deformations; in Ref. 28 that model
was generalized to include two mechanically equivalent fiber
families representing crosslinking phenomena. Also, in Ref. 21
two families of fibers were used to model arterial tissue, but these
fiber families were not assumed to be mechanically equivalent
because a bimodular feature was used.
6
Here, secondary fibers are
used as in Ref. 28 without assuming that the fiber families are
mechanically equivalent so that the bimodular feature can be used
as in Ref. 21.
First, three fiber families are introduced that are parallel to three
mutually orthogonal basis vectors E
1
, E
2
, E
3
in a stress-free ref-
erence configuration; these will be referred to as “principal fi-
bers.” As seen below in Eq. 13, the principal fibers generate
strain invariants for orthotropic materials. Structural tensors
M
1
,M
2
,M
3
are defined as in the case of an orthotropic material
M
1
= E
1
� E
1
, M
2
= E
2
� E
2
, M
3
= E
3
� E
3
8
The unit vectors used to form these structural tensors correspond
to the following anatomical directions: E
1
is parallel to the local
split-line direction, E
3
is perpendicular to the articular surface,
and E
2
is perpendicular to the split-line direction and parallel to
the surface Fig. 1.
Second, two fiber families are introduced in each of the three
6
See the comment below following Eq. 21.
Journal of Biomechanical Engineering APRIL 2007, Vol. 129 / 251
Fig. 2 Schematic of the principal and secondary fiber orienta-
tions in relation to anatomical directions in the 1-2 plane. The
two principal fiber directions are parallel to the unit vectors E
1
and E
2
and the two secondary fiber directions, denoted as E
±12
,
are oriented at angles of ±
12
to the E
1
direction. The weights
of the line elements represent the relative strength of the fiber
directions as predicted by regression analysis; i.e., the princi-
pal fibers along the E
1
direction are the strongest while the
secondary fibers are the weakest.
planes formed by the basis vectors E
1
,E
2
,E
3
; these will be re-
ferred to as “secondary fibers.” As seen below in Eq. 14, the
secondary fibers generate strain invariants that represent strong
interaction terms between the orthotropic strain invariants. Con-
sider the 1-2 plane, which contains the E
1
,E
2
unit vectors Fig.
2. The secondary fiber families are defined to lie at angles of
±
12
to the E
1
direction. These secondary fiber directions are
denoted as E
+12
,E
−12
and are expressed as
E
±12
= cos
12
E
1
± sin
12
E
2
9
Corresponding structural tensors M
+12
,M
−12
are defined as in
Eq. 8
M
±12
= E
±12
� E
±12
= cos
2
12
E
1
� E
1
+ sin
2
12
E
1
� E
1
± cos
12
sin
12
E
1
� E
2
+ E
2
� E
1
10
The secondary fiber directions introduced in the 1-3 and 2-3
planes are denoted as E
+13
,E
−13
and E
+23
,E
−23
, respectively,
are expressed as
E
±13
= cos
13
E
1
± sin
13
E
3
, E
±23
= cos
23
E
2
± sin
23
E
3
11
Corresponding structural tensors M
+13
,M
−13
,M
+23
, M
−23
are
defined as in Eq. 10
M
±13
= cos
2
13
E
1
� E
1
+ sin
2
13
E
3
� E
3
± cos
13
sin
13
E
1
� E
3
+ E
3
� E
1
M
±23
= cos
2
23
E
2
� E
2
+ sin
2
23
E
3
� E
3
± cos
23
sin
23
E
2
� E
3
+ E
3
� E
2
12
Following Ref. 29, the strain energy function W is assumed to
be an isotropic function of C and the nine structural tensors intro-
duced above. In Ref. 33, a procedure is outlined for obtaining
minimal lists of irreducible scalar invariants for an arbitrary finite
number of symmetric structural tensors; however, that procedure
was only employed for up to six symmetric structural tensors.
Here, only the decoupled first-order scalar invariants are used in
an attempt to obtain a relatively low number of material constants
to prevent the model from becoming overparameterized. Invari-
ants associated with the primary fibers include
2 2
M
1
· C,M
2
· C,M
3
· C = C
11
,C
22
,C
33
=
1
,
2
,
2
3
13
where
A
2
represents the square of the stretch of the material line
element initially oriented along the principal fiber direction E
A
.
Invariants associated with the secondary fibers include
M
±12
· C,M
±13
· C,M
±23
· C
= C
11
cos
2
12
+ C
22
sin
2
12
±2C
12
cos
12
sin
12
,C
11
cos
2
13
+ C
33
sin
2
13
±2C
13
cos
13
sin
13
,C
22
cos
2
23
+ C
33
sin
2
23
±2C
23
cos
23
sin
23
2 2 2
=
±12
,
±13
,
±23
14
where
+
2
AB
and
−
2
AB
represent the squares of the stretches of the
material line elements initially oriented along the secondary fiber
directions E
+AB
and E
−AB
, respectively. Note that these invariants
represent strong interaction or coupling terms for the orthotropic
strain invariants as discussed in Refs. 26,34; for example, the
invariant
+12
2
is a function of the invariants C
11
=
1
2
and C
22
=
2
2
, thereby coupling these invariants. It is important to empha-
size that this approach adopts a lesser level of material symmetry
than orthotropy because the invariants in Eq. 14 can easily be
shown to not be invariant under transformations due to reflections
about three orthogonal planes.
7
Including the three principal invariants of C, the strain energy
function for the model proposed here can be expressed as a func-
tion of 12 invariants
2 2 2 2 2 2
W = W
ˆ
tr C,tradj C,det C,
1
,
2
,
3
,
±12
,
±13
,
±23
15
where tr is the trace operator, adj C= det CC
−1
is the adjugate of
C, and det is the determinant operator.
Bimodular Polyconvex Strain Energy functions. In recent
years, polyconvex strain energy functions have been proposed for
anisotropic materials 26,27,34,36; discussion of the rationale for
using polyconvex strain energy functions is in the “Introduction.”
A sufficient condition for polyconvexity is as follows 26:ifthe
strain energy function WF satisfies the additive decomposition
WF = W
1
F + W
2
adj F + W
3
det F 16
and each of the functions W
1
F,W
2
adj F,W
3
det F is a con-
vex function of F , adj F , det F , respectively, then WF is poly-
convex. Furthermore, addition of two or more polyconvex func-
tions results in a polyconvex function.
Here, W is additively decomposed into two terms W
O
and W
BIM
representing nonbimodular and bimodular contributions, respec-
tively. In general, W
O
can be anisotropic; a general polynomial
form is proposed in Ref. 36. Here, a simple isotropic function is
adopted from Ref. 36 for W
O
1
tr C −3 + tradj C −3 −3 lndet C 17W
O
=
2
This term is polyconvex if
is positive and contributes a stress
term as follows 36
S
O
=
I − det CC
−2
+ det Ctr C
−1
−3C
−1
18
Then, it is assumed that W
BIM
represents the collagen network
molecules that account for all of the tissue anisotropy. To model
tension–compression asymmetry, it is assumed that all fiber fami-
lies can only support tensile stresses; consequently, a total of nine
surfaces of discontinuity are used
7
If the secondary fibers are not bimodular and assumed mechanically equivalent,
then the symmetry reduces to orthotropy as in 30,35.
252 / Vol. 129, APRIL 2007 Transactions of the ASME
g
1
= M
1
· C −1=0, g
2
= M
2
· C −1=0, g
3
= M
3
· C −1=0
g
±12
= M
±12
· C −1=0, g
±13
= M
±13
· C −1=0
g
±23
= M
±23
· C −1=0 19
For example, the surface g
1
=M
1
·C −1=
1
2
−1=0 defines a five-
dimensional hyperplane that divides the C space into two half-
spaces corresponding to tensile and compressive strains in the
principal fiber direction E
1
. The following bimodular form is used
1 1 1
1
1
1
2
−1
3
+
2
2
2
2
−1
3
+
3
3
3
2
−1
3
W
BIM
=
6 6 6
1
2
1
2
+ −1
3
+ −1
3
6
±12
±12
±12
6
±13
±13
±13
+
1
2
−1
3
20
6
±23
±23
±23
where
1
,
2
,
3
,
±12
,
±13
,
±23
are six material constants that
represent bimodular terms via the definitions
1
0if
1
1
1
1
=
0
if
1
1
+12
0if
+12
1
+12
+12
= , etc 21
0
if
+12
1
and the angles
12
,
13
,
23
that appear in Eq. 14 can be re-
garded as three additional material constants. In this general for-
mulation, the two secondary fiber families in any of the three
planes will have the same stiffness if both are active i.e.,
+12
=
−12
, but are not assumed to be mechanically equivalent as
defined in Ref. 29 because in some shearing deformations one
fiber family may be in tension while the other may be in compres-
sion. Considering Eq. 21, it is evident that each of the material
constants are related to one of the surfaces of discontinuity de-
fined in Eq. 19. For example, the material constant
1
defines a
strain energy term that can jump across the surface g
1
=M
1
·C
−1=
1
−1=0.
In the Appendix, the proposed strain-energy function W
BIM
is
shown to satisfy both the bimodular stress continuity and polycon-
vexity conditions when the material constants are defined as in
Eq. 21. Since preliminary statistical results suggested that the
model defined by Eq. 20 was overparameterized given the ex-
perimental dataset,
8
it is further assumed that the material con-
stants associated with the secondary fibers are equal when ac-
tive; i.e.,
±AB
=
. This reduced strain energy function contributes
a stress term as
2
S
BIM
=
1
1
1
2
−1
2
M
1
+
2
2
2
2
−1
2
M
2
+
3
3
3
2 2
−1
2
M
3
+
±12
±12
−1
2
M
±12
+
±13
±13
2
−1
2
M
±13
+
±23
±23
−1
2
M
±23
22
so that the stress constitutive equation is defined by Eqs. 18 and
22 as
S = S
O
+ S
BIM
23
with a total of eight material constants
,
1
,
2
,
3
,
,
12
,
13
,
23
.
Experimental Data. A hypothetical experimental dataset was
developed that approximates the equilibrium elastic response of
the solid matrix of adult human cartilage in the surface region
Tables 1 and 2, assuming homogeneous tissue composition and
elastic properties of test specimens that are free of residual stress.
In order to construct enough data to prevent the models presented
here from being overparameterized, it was necessary to use data
Table 1 Values of tangent Young’s modulus „MPa… in tension
at 0% strain „E
+0
… and 16% strain „E
+0.16
… and in compression at
0% strain „E
−0
… and 16% strain „E
−0.16
… in the1,2,and3direc-
tions for the experimental dataset used.
Direction
Parameter 1 2 3
E
+0
E
+0.16
E
−0
E
−0.16
7.8
42.8
0.18
0.26
5.9
26.3
0.18
0.26
1.2
9.0
0.18
0.26
from several studies representing different anatomic locations,
species, etc.
Based on how mechanical properties were calculated in the
studies used here, the first Piola–Kirchhoff stress and Biot strain
tensors are used. In particular, the first Piola–Kirchhoff stress nor-
malizes load by original cross-sectional area. Also, the Biot strain
B B
tensor has principal strain components e.g., E
11
B
that,E
22
,E
33
correspond to the definition of the infinitesimal strain tensor
B
e.g., E
11
=
11
=
1
−1, etc.. Consequently, Poisson’s ratios de-
fined in terms of the Biot strain tensor correspond to the Poisson’s
ratios defined in terms of
used in the studies mentioned below.
For example, the Poisson’s ratio
12
is defined here as
12
=− E
22
B
/E
11
B
=−
2
−1/
1
−1 =−
22
/
11
24
The data used corresponds to uniaxial tension UT and uncon-
fined compression UCC experiments along three directions: 1
=parallel to the split-line; 2=perpendicular to the split line and
parallel to the surface; and 3=perpendicular to the surface Fig.
1. Exponential functions were used to generate axial i.e., along
the direction of applied loading stress–axial strain data and linear
functions were used to generate transverse strain-axial strain data
from 0% to 20% strain in 2% increments. UT axial stress–axial
strain data were adopted from Refs. 17,18. UT Poisson’s ratios
were assumed based on the results of several studies
13,15,16,18,37. UCC axial stress–axial strain data were adopted
from Refs. 9,14,18 and assumed to be the same in all three
directions.
9
UCC Poisson’s ratios were assumed based on the re-
sults of Refs. 10,19,38,39.
It is important to note that this hypothetical dataset includes not
only substantial anisotropy and asymmetry in the axial stress-
strain response Table 1, but also substantial anisotropy and
asymmetry in the Poisson’s ratios Table 2. In particular, the
Poisson’s ratios in UT can be approximately two orders of mag-
nitude greater than those in UCC, and in both UT and UCC the
Poisson’s ratios
13
and
23
have been measured to be greater than
those in other directions; see Refs. 13,37 for UT and Refs.
10,19 for UCC.
Regression Analysis. A simultaneous nonlinear regression al-
gorithm was performed in Mathematica Wolfram, V5.0 based on
an approach developed in Refs. 23,30. The Levenberg–
Marquardt method is used to minimize an error term representing
the sum of squared differences between theoretical and experi-
mental stress values. Although additional models were studied,
10
only results of three regression analyses are presented Table 3.
An eight-parameter model 8-PAR defined by Eqs. 18 and 22
was studied. To provide a comparison with a model that does not
use the strong interaction terms, a four-parameter model 4-PAR
was studied, for which the material constant
was set equal to
zero. For the 8- and 4-PAR models, the assumed Poisson’s ratios
were used to prescribe the transverse strains for UT and UCC.
9
This limitation is addressed in the “Discussion.”
8
These results are summarized in the “Discussion.”
10
The results of other models are summarized in the Discussion.
Journal of Biomechanical Engineering APRIL 2007, Vol. 129 / 253
Table 2 Numerical values of Poisson’s ratios in tension „
+ij
… and compression „
−ij
… obtained
from regression analysis; i =loading direction, j=direction of transverse strain component. The
range given corresponds to values at 0% and 20% strain. The assumed values are given for
comparison in the column labeled Range; tensile values were assumed to be linear functions
of strain and the compressive values were assumed to be constant.
Model
Parameter Range 8-PAR 4-PAR 8-PAR-B
+12
0.5– 1.0 0.44– 0.83 0.25– 0.26 0.79– 1.21
+13
1.0– 2.0 0.75– 1.59 0.25– 0.26 1.22– 1.73
+21
0.5– 1.0 0.46– 0.86 0.25– 0.26 0.96– 1.54
+23
1.0– 2.0 0.75– 1.59 0.25– 0.26 1.08– 1.82
+31
0.5– 1.0 0.32– 0.45 0.25– 0.26 0.63– 0.81
+32
0.5– 1.0 0.32– 0.45 0.25– 0.26 0.71– 0.92
−12
0.1
0.10– 0.03 0.11– 0.06 0.01–−0.06
−13
0.2
0.19– 0.09 0.60– 0.10 0.08– 0.04
−21
0.1
0.09– 0.03 0.12–−0.06 0.01–−0.06
−23
0.2
0.20– 0.09 0.21– 0.11 0.12– 0.12
−31
0.1
0.11– 0.04 0.13– 0.05 0.05– 0.02
−32
0.1
0.13– 0.05 0.15– 0.06 0.03– 0.01
Then, a composite function representing a total of 18 equations
was derived: six axial stress–axial strain equations three each in
UT and UCC, and 12 transverse stress–axial strain equations cor-
responding to the traction-free boundary conditions six each in
UT and UCC. To provide a comparison with a model that does
not explicitly include the 12 traction-free boundary condition
equations, an additional regression with the eight-parameter
model was performed that included only the six axial stress–axial
strain equations obtained after prescribing the transverse strains
8-PAR-B. In all cases, the UCC stress values were weighted by
multiplying each stress value by 100, since the UT stress response
is two orders of magnitude greater than the UCC stress response.
After the nonlinear regression analysis was performed, the de-
termined model parameters were used to derive numerical solu-
tions to the UT and UCC boundary-value problems, including
theoretical predictions of Poisson’s ratios.
Results
The numerical values for the material constants are presented in
Table 3. The nonlinear regression analyses always converged to
results consistent with the stability criteria; i.e.,
,
1
,
2
,
3
,
were all positive. The calculated error terms were 0.773 and 0.876
for the 8-PAR and 4-PAR models, respectively, and 0.206 for the
8-PAR-B model. It is important to note that this latter error term
cannot be directly compared to the others because fewer equations
were used in the nonlinear regression.
Table 3 Numerical values for the material parameters ob-
tained from regression analysis. The constants „
,
1
,
2
,
3
,
…
are in MPa and the constants „
12
,
13
,
23
… are in degrees.
Model
Parameter 8-PAR 4-PAR 8-PAR-B
0.035 0.035 0.035
23.13 23.88 21.92
1
14.83 15.61 1.78
2
3.86 4.51 3.26
3
12.07 — 436.8
45 — 46
12
35 — 40
13
35 — 42
23
The predictions of the 8-PAR and 4-PAR models Figs. 3 and
4, as well as the 8-PAR-B model for axial stresses were qualita-
tively similar, with one exception. For the 8-PAR-B model that
did not explicitly include the traction-free boundary conditions,
Fig. 3 Predictions of the eight-parameter model „8-PAR… for
the uniaxial tension „UT… response in the 1, 2, and 3 directions
and the unconfined compression „UCC… response in the 1 di-
rection. The theoretical UCC curves in the 2 and 3 directions
are within 1% of the curve shown. UCC stress and strain val-
ues, although negative by definition, are plotted as positive
numbers.
Fig. 4 Predictions of the four-parameter model „4-PAR… for the
uniaxial tension „UT… response in the 1, 2, and 3 directions and
the unconfined compression „UCC… response in the 1 direction.
The theoretical UCC curves in the 2 and 3 directions are within
3% of the curve shown. UCC stress and strain values, although
negative by definition, are plotted as positive numbers.
254 / Vol. 129, APRIL 2007 Transactions of the ASME