Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

53, 4, pp. 895-910, Warsaw 2015
DOI: 10.15632/jtam-pl.53.4.895

DETERMINATION OF MATERIAL PARAMETERS OF ISOTROPIC
AND ANISOTROPIC HYPER-ELASTIC MATERIALS USING

BOUNDARY MEASURED DATA

Maedeh Hajhashemkhani, Mohammad R. Hematiyan
Shiraz University, Department of Mechanical Engineering, Shiraz, Iran

e-mail: maedeh.hajhashemkhani@gmail.com; mhemat@shirazu.ac.ir

Identification of mechanical properties of isotropic and anisotropic materials that demon-
strate non-linear elastic behavior, such as rubbers and soft tissues of human body, is critical
formany industrial andmedical purposes. In this paper, amethod is presented to obtain the
mechanical constants of Mooney-Rivlin and Holzapfel hyper-elastic material models which
are employed to describe the behavior of isotropic and anisotropic hyper-elastic materials,
respectively. By using boundarymeasured data from a samplewith non-standard geometry,
and by using an iterative inverse analysis technique, the material constants are obtained.
The method uses the results of different experiments simultaneously to obtain the mate-
rial parameters more accurately. The effectiveness of the proposedmethod is demonstrated
through three examples. In the two first examples, the simulated measured data are used,
while in the third example, the experimental data obtained from a polyvinyl alcohol sample
are used.

Keywords: hyper-elstic material, inverse analysis, anisotropic, polyvinyl alcohol

1. Introduction

Many materials such as rubbers, elastomers and tissues of human body experience large defor-
mation under small loads. Further, some materials such as wood, fiber-reinforced composites
and some body tissues like arterial walls, show anisotropic behavior in addition to large elastic
deformation due to the presence of preferred directions in their structure. The large strains in
the response of these materials clarify that their mechanical behavior is nonlinear, and getting
back to the reference configuration after load removal demonstrates their elastic response. Be-
cause of the wide use of these materials in industry and the crucial role of different tissues in
human body, presenting a model with known parameters that can predict the mechanical be-
havior of these materials is essential. When these materials experience small deformations (less
than 2 to 5 percent), their mechanical behavior can be modeled using common linear elastic
models (Czabanowski, 2010), but under large deformations, their mechanical response must be
represented by nonlinear models such as hyper-elastic material models. Therefore, many efforts
have beenmade to present constitutive laws thatmodel the behavior of thesematerials proper-
ly. Pamidi and Advani (1987) proposed nonlinear constitutive relations for human brain tissue.
Fung (1993) showed that the elastic properties of rabbits’ mesentery could be simply modeled
as an exponential function. Holzapfel andGasser (2000) considered arterial walls as thin-walled
cylinders and presented a constitutive equation for describing their behavior. They obtained the
numerical values of the constants of theirmodel using experimental data. A constitutive law for
arterial layers with distributed collagen fiber orientation was presented in the work of Gasser et
al. (2006).
Many researchers used hyper-elastic models to predict material behavior especially for hu-

man body tissues. Moulton et al. (1995) used a combination of finite element model, nonlinear



896 M. Hajhashemkhani, M.R. Hematiyan

optimization algorithm, and a set of strains obtained by magnetic resonance imaging (MRI) to
determine passivemyocardial material properties.Miller andChinzei (1997) conducted pressure
experiments on brain tissue to model the destructions that may happen during brain surgery.
They obtained force-displacement diagrams from the experiments. Then, they computed the
unknown constants for different hyper-elastic models using the least squaresmethod. Krouskop
et al. (1998) conducted some pressure experiments on different parts of breast and prostate
tissue, such as fat, and obtained Young’s modules for each part.
Hartman (2001) conducted some tension and torsion tests on different samples made of

rubber to identify their properties based on incompressible Rivlin’s hyper-elastic model. Ogden
et al. (2004) used a non-linear least squares optimization method to obtain the constants of
hyper-elastic models for incompressible materials by fitting experimental data to their model.
In the work of Hu andDesai (2004), a tissue indentation test was conducted on liver to identify
its biomechanical properties.
Balaraman et al. (2005) conducted different experiments on 19 samples of muscle tissue and

used the experimental data in a finite element analysis. They obtained mechanical properties
of the muscle by an inverse method using Taguchi’s approach. Mehrabian (2008) estimated the
constants of an incompressible hyper-elastic material based on three different energy functions.
Some samples made of polyvinyl alcohol were used in pressure experiments. The tests were
modeled in the finite element analysis and the unknown constants were computed using the
least square method.
Ahn et al. (2008) conducted biomechanical experiments on macro and micro liver samples

to characterize theirmechanical behavior using a nonlinear least squaresmethod and an inverse
finite element (FE) method based on a parameter estimation algorithm. In the work of Mesa-
Múnera et al. (2012), several hyper-elastic models were used to predict and model mechanical
behavior of a silicone rubber. They conducted a uniaxial compression test on a phantom with
mechanical properties close to brain tissue. They could obtain force-displacement diagrams for
compression tests on the phantom.
Rauchs et al. (2010) used a depth-sensing spherical indentation technique to determine the

parameters of different rubber materials by using an inverse method based on a gradient-
-based numerical optimization method. Czabanowski (2010) conducted some pressure expe-
riments on machinery elastomers to obtain their material properties. Czabanowski obtained
force-displacement diagrams first. Then, by using a function in ABAQUS finite element softwa-
re, he determined the constants of Mooney-Rivlin model for the materials.
Abyaneh et al. (2013) used an inverse method to obtain viscoelastic material properties of

porcine cornea by performing tension and indentation tests. They used a sample with possible
maximum length for the tension test. Using the results of the tension test, material properties
for an isotropic viscoelastic material were obtained using a curve fitting procedure. Due to
insufficiency of this model, to match the indentation test results, an anisotropic model with
the results of the indentation test was also used in an inverse algorithm to obtain anisotropic
parameters. Parameters of the isotropicmodel were the initial guesses for the inverse algorithm.
Baker andShrot (2013) proposedanewmethod for inverse identification ofmaterial parameters.
They used auxiliary quantities to describematerial behavior.
In all of the works reviewed above, the sample that was used to obtain the measurement

data was cut to a cubic or cylindrical shape to provide standard test specimens. Therefore,
these methods cannot be applicable to find the material constants of a hyper-elastic material
with a non-standard shape nondestructively. The inverse method presented in this paper uses
the measured displacement data from a sufficient number of sampling points on the original
hyperelastic body to obtain unknown material constants. If the properties of the material are
dependent on size and geometry, the properties, which are obtained by thismethod, are suitable
for that special shape.Unlike previous researches, themethodpresented in this research uses the



Determination of material parameters of isotropic and anisotropic... 897

results of different experiments simultaneously to obtain the material parameters more accura-
tely. This inverse problem is highly nonlinear and encounters various difficulties. Simultaneous
use ofmeasured data from several experimentswith different load cases reduces the ill-posedness
of the inverse problem. It is observed that using the measured data from several experiments
results in a better solution than the case where only the measured data from one experiment
is used. The Mooney-Rivlin and Holzapfel models are considered for isotropic and anisotropic
materials, respectively. The Tikhonov regularization method is used in the inverse analysis. A
method to determine the initial guesses for the unknown constants is also presented. The inverse
methodneeds a sensitivity analysis, which is carried out using the finite elementmethod (FEM).

2. Materials modeling

2.1. Hyper-elastic materials

For hyper-elasticmaterials, the stress-strain relations are determined using the strain energy
density function ψ which is defined in terms of a deformation gradient or strain tensor. The
derivative of ψ with respect to a component of strain gives its corresponding stress component
(Holzapfel, 2000)

Sij =
∂ψ

∂εij
(2.1)

whereS is the secondPiola-Kirchhoff stress tensor and ε is the Lagrangian strain tensor defined
as follows

εij =
1
2
(Cij − δij) (2.2)

In Eq. (2.2),C is the right Cauchy-Green deformation tensor (C=FTF).F is the deformation
gradient, which can be expressed in terms of the displacement vector u (F=∇u+ I).
Local shape deformation of amaterial element is described byF. The strain energy function

is usually a function of deformation gradientF. Indeed, the function ψ is expressed in terms of
the rightCauchy-Green deformation tensor, i.e.,ψ = ψ(C).The following relation holds between
the second Piola-Kirchhoff and Cauchy stress tensors (Holzapfel, 2000)

S= JF−1σF−T (2.3)

where J =detF. Therefore, according to Eqs. (2.1), (2.2) and (2.3) the Cauchy stress tensor is
defined as follows

σ =2J−1F
∂ψ

∂C
FT (2.4)

For an isotropic material, ψ is dependent onC based on its invariants. The invariants ofC are

I1 = trC I2 =
1
2
[(trC)2− trC2] I3 =detC (2.5)

For an isotropic material ψ is just dependent on I1, I2 and I3 (Holzapfel andOgden, 2010) and
the Cauchy stress equation is expanded as follows

σ =2J−1
(
F
∂ψ

∂I1

∂I1
∂C
FT+F

∂ψ

∂I2

∂I2
∂C
FT+F

∂ψ

∂I3

∂I3
∂C
FT
)

(2.6)



898 M. Hajhashemkhani, M.R. Hematiyan

By obtaining the derivatives of the invariants with respect toC (Holzapfel, 2000) and knowing
thatB=FFT Eq. (2.6) yields to

σ =2J−1[ψ1B+ψ2(I1B−B2)+ I3ψ3I] (2.7)

where ψi = ∂ψ/∂Ii. For an incompressible isotropic material I3 = detF = 1 and the Cauchy
stress tensor is modified as follows (Holzapfel, 2000)

σ =−pI+2F
∂ψ

∂C
FT (2.8)

where p is a scalar identified as hydrostatic pressure. Therefore, the Cauchy stress tensor for an
incompressible material with respect toB is expressed as follows

σ =−pI+2ψ1B+2ψ2(I1B−B2) (2.9)

2.2. Isotropic hyper-elastic materials

The polynomial hyper-elastic model is a proper model for rubbers and other soft mate-
rials with isotropic behavior. In this research, the Mooney-Rivlin hyper-elastic model, which is
obtained from the polynomial model, is used to predict the mechanical behavior of isotropic
materials. In this model, density of strain energy for an incompressible material is a function of
the first and second invariants of the left Cauchy-Green deformation tensor. The strain energy
function for the polynomial model is defined as follows (Rivlin and Saunder, 1951)

ψ =
n∑

i,j=0

Aij(I1−3)i(I2−3)j (2.10)

where Aij are the constants of the material and A00 = 0. The Mooney-Rivlin strain energy
function is obtained by considering n =1 and A11 =0 in equation (2.10), i.e.

ψ = A10(I1−3)+A01(I2−3) (2.11)

By considering relations (2.9) and (2.11), theCauchy stress for theMooney-Rivlin hyper-elastic
material is obtained as follows

σ =−pI+2A10B+2A01(I1B−B2) (2.12)

2.3. Anisotropic hyper-elastic materials

In the structure of human soft tissues, the presence of collagen fibers (Unnikrishnan, 2012)
causes thematerial to have one ormore preferred direction. This preferred direction is presented
here by M. In this case, the strain energy density is a function of both C and M. Two more
dependent pseudo invariants are defined for these materials (Holzapfel andGasser, 2000)

I4 =M(CM) I5 =M(C
2M) (2.13)

where, for example,CM represents the action of the second order tensorC on the vector M.
For an incompressible material reinforced with one family of fibers, ψ is dependent on I1,

I2, I4 and I5. In this case, the Cauchy stress has two additional terms which show the effect
of anisotropy. Therefore, the Cauchy stress tensor can be expressed as follows (Holzapfel and
Ogden, 2010)

σ =−pI+2ψ1B+2ψ2(I1B−B)2+2ψ4m⊗m+2ψ5[m⊗Bm+Bm⊗m] (2.14)



Determination of material parameters of isotropic and anisotropic... 899

where ⊗ denotes the dyadic product of two vectors and m=FM is the deformed form of the
vectorM in current configuration. For some tissues like arterial walls, two families of fiberswith
different directions can be detected within the tissue. M′ is considered to be the unit vector in
the direction of the second family of fibers. In this case, three more invariants are considered,
which are expressed as follows (Holzapfel and Ogden, 2010)

I6 =M
′(CM′) I7 =M

′(C2M′) I8 = [M(CM
′)](MM′) (2.15)

In this case, the Cauchy stress is expressed as follows (Holzapfel andOgden, 2010)

σ =−pI+2ψ1B+2ψ2(I1B−B)2+2ψ4m⊗m+2ψ5[m⊗Bm+Bm⊗m]
+2ψ6m

′⊗m′+2ψ7[m′⊗Bm′+Bm′⊗m′]+2ψ8(M⊗M′)(m⊗m′+m′⊗m)
(2.16)

In this research, the Holzapfel’s strain energy function for anisotropic materials is used. This
function is described as follows (Holzapfel and Ogden, 2010)

ψ = R10(I1−3)+
k1
k2
{exp[k2(I∗4 −1)2]−1} (2.17)

R10, k1, k2 and κ are material constants and I∗4 is defined as follows

I∗4 = κI1+(1−3κ)I4 (2.18)

In this model, it is assumed that the direction of each family of fibers is dispersed around a
mean direction. This dispersion is shown by κ (0¬ κ ¬ 1/3) (Holzapfel and Ogden, 2010). By
using relations (2.16), (2.17) and (2.18), the Cauchy stress tensor can be obtained.

3. Computation of unknown material parameters of hyper-elastic materials

In this paper, an inversemethod is used to obtain the constants of hyper-elasticmaterialmodels.
The inverse analysis is used to convert themeasureddata to some informationabout thematerial
or system under study. In the inverse analysis, the outputs of the systemmay be available from
an experiment but loading parameters, material properties, geometry of structure, boundary
conditions or a combination of these factors have to be determined (Liu and Han, 2003). In
the present study, unknowns are constants of hyperelastic materials. These unknowns are found
using the measured displacements at some sampling points on the boundary of the body. The
measured displacements are collected from a few load cases (for example, 3 cases). Then by
using the inverse analysis in an iterative process, the constants are obtained. A proper initial
guess is required for the iterative process.

3.1. Inverse analysis

In this section, the inverse formulation for evaluation of constants of isotropic andanisotropic
hyper-elastic material models is presented.
For an incompressible material, the Mooney-Rivlin hyper-elastic model constants are A10

and A01. Therefore, the vector of unknowns for the problem is defined as follows

x=
[
x1 x2

]T
(3.1)

where x1 = A10 and x2 = A01. For the Holzapfel hyper-elastic model, these constants are R10,
k1, k2 and κ. The vector of unknowns for this case is defined as follows

x=
[
x1 x2 x3 x4

]T
=
[
R10 k1 k2 κ

]T
(3.2)



900 M. Hajhashemkhani, M.R. Hematiyan

In order to obtain material constants, a few experiments (load cases) are performed on the
material. These experiments must be different with each other to provide suitable data for the
inverse analysis. The number of these experiments is 2, 3 ormore. In this work, 3 load cases are
considered for describing the formulation.
In each experiment, the displacement at some boundary points is measured and restored in

vectors Y(1),Y(2) andY(3). If they are N1, N2 and N3 measured data in load cases 1, 2 and 3,
respectively, Y(i) (i =1,2,3) are defined as follows

Y(i) =
[
Y
(i)
1 Y

(i)
2 · · · Y

(i)
Ni

]
i =1,2,3 (3.3)

By solving the problem numerically using the estimated constants, the displacement vector at
sampling points is obtained and stored inD(1),D(2) andD(3) which are defined as follows

D(i) =
[
D
(i)
1 D

(i)
2 · · · D

(i)
Ni

]
i =1,2,3 (3.4)

In order to find the constants of the material, the Tikhonov regularization method is used. In
this method, the objective function Π is defined as follows

Π =(Y−D)T(Y−D)+αXTX (3.5)

where the vectors Y andD are defined as follows

Y=
[
Y(1) Y(2) Y(3)

]T
D=

[
D(1) D(2) D(3)

]T
(3.6)

In equation (3.5), α is the regularization parameter. The objective function Π should be mi-
nimized to find the unknown vector X. By taking the derivative of Π with respect to X and
equating to zero, the following equation is obtained

∂Π

∂X
=−2ST(Y−D)+2αX=0 (3.7)

In Eq. (3.7), S is the sensitivity matrix and is defined as follows

S(L) =




S
(L)
11 S

(L)
12 · · · S

(L)
1q

S
(L)
21 S

(L)
22 · · · S

(L)
2q

...
...
...

...
S
(L)
NL1

S
(L)
NL2

· · · S(L)NLq




(3.8)

where L is the number of the load case and q is the number of unknowns or the number of
components ofX.
The components of the sensitivity matrix are defined as follows

S
(L)
ij =

∂D
(L)
i

∂Xj
(3.9)

In order to compute the constants in an iterative process and calculate the sensitivity matrix,
a program is written in python for the finite element software ABAQUS. This program uses
the finite difference method to calculate the sensitivity matrix. For this purpose, each material
constant is changed by a very small value and displacements at sampling points are obtained.
The differences between these displacements are used to compute the sensitivity coefficients.
The sensitivity coefficient S

(L)
ij is computed through the following finite difference equation

S
(L)
ij =

D
(L)
i (X1, . . . ,Xj +µXj, . . . ,Xq)−D

(L)
i (X1, . . . ,Xq)

µXj
(3.10)

where µ is a small value such as 0.001.



Determination of material parameters of isotropic and anisotropic... 901

For example, if there are fourunknownmaterial constants, the sensitivity of thedisplacement
at sampling point 3 with respect to material constant 2 is computed as follows

S
(2)
32 =

D
(2)
3 (X1,X2+µX2,X3,X4)−D

(2)
3 (X1,X2,X3,X4)

µX2
(3.11)

It is possible to obtain the unknownvectorX through equation (3.7). Let the vector of unknown
constants in the current step of the iterative process is represented by Xc and the correspon-
ding displacement vectors for load cases 1, 2, and 3 are represented by D

(1)
c , D

(2)
c and D

(3)
c ,

respectively. The total displacement vector Dc is defined as follows

Dc =
[
D
(1)
c D

(2)
c D

(3)
c

]T
(3.12)

The total displacement vector in the new step, i.e. D, can be expressed as follows

D=Dc +S(X−Xc) (3.13)

whereX is the vector of unknown constants in the new step. By using Eqs. (3.13) and (3.7) the
following equation is obtained

X= [STS+αI]−1[ST(Y−Dc)+STSXc] (3.14)

Equation (3.14) is used iteratively through the following equation to obtain the constants

Xk+1 = [(Sk)TSk +αkI]−1[(Sk)T(Y−Dk)+(Sk)TSkXk] (3.15)

k and k+1 are the numbers of iterations, and the convergence rule is defined as follows

‖Xk+1−Xk‖¬ e (3.16)

where e is the specified tolerance. In the cases where a large number of sampling points and/or
many load cases are considered, selecting a small value or even zero for α results in stable
solutions.
Inequality constraints on material constants (for example 0¬ κ ¬ 1/3) are not considered

in the inverse formulation of this study. In other words, the inverse formulation of this research
uses an unconstrained optimization formulation. However, the obtained results must satisfy all
the constraints. Otherwise, the problem has to be solved with different initial guesses.

3.2. Initial guesses

Inorder to start the iterative process andobtain thematerial constants, proper initial guesses
are required. Selecting appropriate initial guesses reduces the number of iterations and increases
the convergence rate. To find proper initial guesses for the unknown parameters of the hyper-
-elastic material, we try to find unknown Young’s modulus E and Poisson’s ratio ν for an
isotropic linear elastic material, which can approximately reconstruct the measured data for
the original hyper-elastic material. For this purpose, themeasured data Y are used in a simple
inverse analysis to find appropriate values for the two parameters of the pseudo linear elastic
material. Then, the calculated parameters are used to compute initial guesses for the original
material parameters. A similar approach has been used byHematiyan et al. (2012) for providing
initial guesses for identification of material constants of linear elastic anisotropic materials.
The pseudomaterial is assumed incompressible with ν =0.5. The displacement vector obta-

ined fromtheexperiments on thehyper-elasticmaterial, i.e.Y, is assumedtobethedisplacement
vector for the pseudomaterial with unknownE. Theunknownparameter E can be simply found
as follows.



902 M. Hajhashemkhani, M.R. Hematiyan

The displacement vector at the sampling points of the pseudomaterial can be computed as

D=
D

E
(3.17)

where D is the displacement vector at the same sampling points for a linear elastic material
with E =1Pa and ν =0.5. The vectorD can be computed using a direct analysis by the FEM.
To find a suitable value of E, weminimize the difference betweenD andY. This can be carried
out byminimizing the following expression

F =
(
Y−D

E

)T(
Y−D

E

)
(3.18)

which gives

E =
D
T
D

D
T
Y

(3.19)

Now, we have found elastic constants of a pseudo linear elastic isotropicmaterial for the inverse
problem. After obtaining the Youngmodulus from equation (3.19), by selecting different values
for ε, a set of stress-strain data is provided using the relation σ = Eε. This set of data is used
to find initial guesses for the unknown constants of the original hyper-elastic material.
Consider a case of uniaxial loading with the axial stress σ which causes the axial strain ε in

direction 1 for the pseudomaterial. The deformation of thematerial can be expressed as follows

x1 = λ1X1 x2 = λ2X2 x3 = λ3X3 (3.20)

where (X1,X2,X3) are rectangular Cartesian coordinates that identify material particles in
some unstressed reference configuration, (x1,x2,x3) are the corresponding coordinates after
deformationwith respect to the sameaxes, and the coefficients (λ1,λ2,λ3) are principal stretches
of the deformation. Therefore, the deformation gradient tensor Fij = ∂xi/∂Xj is expressed as
follows

F=




λ1 0 0
0 λ2 0
0 0 λ3



 (3.21)

Substituting Bij = FikFjk in equation (2.12), the following equation is obtained

σ =2
(
λ21−

1
λ1

)(
A10+

1
λ1
A01

)
(3.22)

On the other hand, for the pseudomaterial we have

ε =
∂(x1−X1)

∂X1
=1−λ1 (3.23)

Substituting λ1 =1+ε into Eq. (3.22) results in

σ =2
[
(1+ε)2−

1
1+ε

](
A10+

1
1+ε

A01
)

(3.24)

Equation (3.24) can be written in the following form

Y = A10X +A01 (3.25)



Determination of material parameters of isotropic and anisotropic... 903

where Y and X are

Y =
σ(1+ε)

2
[
(1+ε)2− 1

1+ε

] X =1+ε (3.26)

By using the previously generated set of stress-strain data, a set of X-Y data is computed using
Eq. (3.26). Then a line is fitted through the X-Y data by a simple linear regression. Considering
Eq. (3.25) and by using the coefficients of the fitted line, the values of A10 and A01 are obtained.
These values are considered as initial guesses for the inverse analysis.
A method to obtain the initial guesses for constants of anisotropic hyper-elastic materials

is also presented here. These constants are R10, k1, k2, and κ, see Eqs (2.17) and (2.18). Since
0¬κ ¬ 1/3, we simply select the value of 0.25 as an initial guess for κ. To simplify the process
of obtaining the initial guesses we also set k1 = k2. Therefore, only the initial guesses for R10
and k1 should be determined. Substituting Eqs. (2.17) and (2.18) into Eq. (2.16), the Cauchy
stress is expressed as follows

σ =2R10B+k1
(I1
4
+
I4
4
−1

)
exp

[
k1
(I1
4
+
I4
4
−1

)2]
(B+m⊗m) (3.27)

As discussed for the isotropic case, by using Eq. (3.19) Young modulus is obtained for the
pseudo incompressible linear elasticmaterial. Using the calculatedYoung’smodulusandHooke’s
law, the components of the stress tensor σ and the strain tensor ε at two different points are
computed. The strain tensor is used in Eq. (2.2) to obtain the right Cauchy-Green deformation
tensorC fromwhich thedeformationgradient tensorFand the tensorB canbeobtained.Having
F and M, the vector m=FM is also known and the unknowns in Eq. (3.27), i.e. R10 and k1,
can be found. These computed values are considered as initial guesses for the corresponding
parameters.

4. Numerical examples

In this Section, three examples are presented to show the efficiency of the presented method.
In the two first examples, numerically simulated measured data are used; however, in the third
example, the experimental data obtained from a polyvinyl alcohol sample are used for the
identification of material parameters.

4.1. Example 1. Isotropic hyper-elastic material

In this example, a rectangular plate with a non-circular hole made of an isotropic hyper-
-elasticmaterial is considered.As shown inFig. 1, the body is subjected to three load cases. The
loading location and directions vary in each case to obtain more useful measurement data. The
measurement data were numerically simulated. For this purpose, three direct problems corre-
sponding to the three load cases were solved by the exact material constants and displacement
at sampling points where obtained. Then, some errors with a normal distributionwere added to
the displacement data to account for practical inaccuracies.
The exact material constants of the Mooney-Rivlin hyper-elastic model for the material

are assumed 80Pa and 20Pa for A10 and A01, respectively. The initial guesses for the ma-
terial constants obtained from the method presented in Section 3.2, are A10 = 90.34Pa and
A01 =3.497Pa. These initial guesses are used in the developed program, and by using equation
(3.15) the unknowns are updated in each step until they satisfy convergence rule.
At first, the problemwas solved without anymeasurement error. Table 1 shows the number

of iterations and obtained values for the constants when the results of the load cases where
used separately (1, 2 or 3) or together (1+2 or 1+2+3) for the inverse analysis. The tolerance



904 M. Hajhashemkhani, M.R. Hematiyan

Fig. 1. Load cases applied to the material

for the convergence rule was assumed 0.01. Then the problem was solved by applying 1 and 5
percent error in the measurement data. These errors are applied in the exact solution of the
problem using the normal distribution function. The results of these cases are shown in Table 1
as well. The convergence process for the case with 5% measurement error is shown in Fig. 2.
From the results of Table 1, it is clear that when there is no error in the measurements, the
solution for all load cases converge easily to thematerial constants and the number of iterations
is almost the same whether the load cases are solved together or separately. Further, when the
error percentage inmeasurements increases, the number of iterations and errors in the obtained
results increase too. In addition, the number of iterations is less for the cases that the problem
is solved by using 2 or 3 load cases. Moreover, it is seen that the obtained results are more
accurate whenmore than one load case is considered in the inverse analysis.

Table 1.Results formaterial constants with different values of measurement errors, Example 1

Load case 1 2 3 1+2 1+2+3

without
measu-
rement
error

No. of iterations 6 6 6 5 5

Constants [Pa]
(Error [%])

A10=80 A10=80 A10=80 A10=80 A10=80
(0) (0) (0) (0) (0)

A01=20 A01=20 A01=20 A01=20 A01=20
(0) (0) (0) (0) (0)

1%
measu-
rement
error

No. of iterations 8 12 8 6 7

Constants [Pa]
(Error [%])

A10=77.17 A10=80.36 A10=77.66 A10=81.15 A10=80.81
(3.5) (0.4) (2.9) (1.4) (1)

A01=23.18 A01=19.52 A01=22.15 A01=18.74 A01=19.04
(15.9) (2.4) (10.7) (6.3) (4.8)

5%
measu-
rement
error

No. of iterations 10 11 17 4 5

Constants [Pa]
(Error [%])

A10=70.058 A10=72.97 A10=75.66 A10=86.43 A10=80.53
(12.4) (8.7) (5.4) (8) (0.6)

A01=31.09 A01=25.93 A01=25.18 A01=12.60 A01=19.39
(55.4) (29.6) (25.9) (36.9) (3)

4.2. Example 2. Anisotropic hyper-elastic material

In this example, the same geometry of the previous example but with an anisotropic hyper-
-elasticmaterial is considered.The exactmaterial constants of theHolzapfel hyper-elasticmodel
for this material are assumed 7.64Pa, 996.6, 524.6 and 0.226 for R10, k1, k2 and κ, respectively.
The direction of the fibers in a hyper-elasticmaterial that determines theM andM′ vectors can
be obtained from two-dimensional images of soft biological tissues (Schriefl et al., 2012). In this



Determination of material parameters of isotropic and anisotropic... 905

Fig. 2. Convergence of (a) A10 and (b) A01 for different load cases with 5% error in measurements

example, these vectors are assumed tobe known (M=0.766e1−0.643e2,M′ =0.5e1−0.866e2).
The initial guesses for the material constant, which are obtained from the method presented in
Section 3.2, are R10 =10.65Pa, k1 = k2 =144.2 and κ =0.25. At first, the problemwas solved
without anymeasurement error. Table 2 shows the number of iterations and obtained constants
when the results of different load cases were used separately or together. The tolerance for the
convergence rulewas assumed0.01.Then, the problemwas solved by considering 1 and5percent
error in themeasurement data. The results of these cases are shown inTable 2. The convergence
process of constants for the case with 1% error in measurements is shown in Fig. 3. From the
results given in Table 2, it is seen that the inverse solution diverges for load cases 2 and 3 even
without any measurement error. Further, it is seen that when an error exists in measurements,
the solution cannot be obtained by only one load case. It is also observed that the number of
iterations for the solution with 3 load cases is smaller than the solution with 2 load cases.

Table 2.Results for material constants with different values ofmeasurement errors, Example 2

Load case 1 2 3 1+2 1+2+3

without
measu-
rement
error

No. of iterat. 7

D
iv
er
ge
d

D
iv
er
ge
d

7 8

Constants
(Error [%])

R10=7.640Pa (0) R10=7.640Pa (0) R10=7.640Pa (0)
k1 =996.6 (0) k1 =996.6 (0) k1 =996.6 (0)
k2 =524.6 (0) k2 =524.6 (0) k2 =524.6 (0)
κ =0.226 (0) κ =0.226 (0) κ =0.226 (0)

1%
measu-
rement
error

No. of iterat.

Diverged

D
iv
er
ge
d

D
iv
er
ge
d

18 7

Constants
(Error [%])

R10=7.645Pa(0.06) R10=7.645Pa(0.06)
k1 =1015 (1.8) k1 =1020 (2.3)
k2 =528.8 (0.8) k2 =527.3 (0.5)
κ =0.227 (0.4) κ =0.227 (0.4)

5%
measu-
rement
error

No. of iterat.

Diverged

D
iv
er
ge
d

D
iv
er
ge
d

37 15

Constants
(Error [%])

R10=7.607Pa (0.4) R10=7.650Pa (0.1)
k1 =1079 (8.2) k1 =926.0 (7)
k2 =509.8 (2.8) k2 =572.8 (9.1)
κ =0.228 (0.8) κ =0.225 (0.4)

4.3. Example 3. A rectangular plate with a hole made of polyvinyl alcohol

In order to verify the proposed method for obtaining material constants of hyper-elastic
materials, a sample was made of polyvinyl alcohol which exhibits hyperelastic mechanical be-
havior. Polyvinyl alcohol is usually used to simulate soft tissues of human body (Hebden et al.,



906 M. Hajhashemkhani, M.R. Hematiyan

Fig. 3. Convergence of (a) R10, (b) k1, (c) k2 and (d) κ for different load cases with 1% error in
measurements

Fig. 4. The polyvinyl alcohol sample used in the experiments

2006). The sample used in the experiments (Fig. 4) is a rectangular plate with dimensions of
80.04mm×10.02 containing a 33.28mm diameter hole.
Two load cases are considered for the polyvinyl alcohol sample. In the first case, the sample

is placed horizontally on the desk and by using a forcemeter, a load is applied to one endwhile
the other end is fixed. This load case is shown in Fig. 5. In this case, a load of 1.2N is applied
to the sample.
In the second case, the sample is hanged by its weight and by using a hook, a load is applied

to the sample. This case is shown in Fig. 6. A load of 1.5N is applied in this case. It should
be mentioned that, in this case, there exists a body force in the same direction of loading too.
Therefore, the deflection of the member in the second case is much more than the deflection in
the first load case.
For measuring the displacement at the sampling points, these points were marked before

the test. By using themeasured displacements, thematerial constants are obtained for a simple
cylindrical sample and for the rectangular sample with a hole. The finite element model (con-



Determination of material parameters of isotropic and anisotropic... 907

Fig. 5. First load case for the rectangular sample

Fig. 6. Second load case for the rectangular sample

Fig. 7. Load application in FEA for the two cases

structed in ABAQUS) and displacement for the two load cases are shown in Figs. 7 and 8. As
it is shown in Fig. 7, the load is not directly applied to thematerial. A rigid plate is used in the
FEAmodel for load application in order to simulate the experiment more accurately.

Thematerial constants obtained from loading the cylindrical sample arepresented inTable 3.
Thematerial constants obtained for the rectangular sample are presented in Table 4, where the
results of the first and second cases are used separately and together. FromTable 4 it is seen that
when each test is used separately to obtain the material constants, the number of iterations is
larger and the constants deviatemore from those obtained for the cylindrical sample inTable 3.
Further, when two tests are used together to obtain the constants, the deviation from constants
of Table 3 and number of iterations are smaller.



908 M. Hajhashemkhani, M.R. Hematiyan

Fig. 8. Displacement contours in the sample for two load cases

Table 3. Material constants obtained for the cylindrical sample

C10 [Pa] C01 [Pa] Number of iterations

4.964 0.541 19

Table 4. Material constants obtained for the rectangular sample with a hole

Test C10 [Pa] C01 [Pa] Number of iterations

1 6.231 0.756 28
2 5.843 0.641 43
1+2 5.006 0.524 13

5. Conclusions

Amethod to obtain the constants of isotropic and anisotropic incompressible hyper-elastic ma-
terials is presented. In the proposed method, by using the measured displacements at some
sampling points on the boundary of a member with a non-standard shape, the unknownmate-
rial constants of the material are computed. The measured data are obtained from more than
one test with different load cases. It is shown that themethod is more efficient when the results
of two or three load cases are simultaneously used to solve the problem.
From the results of the numerical examples, it is concluded that when the measurement

error is zero for isotropic materials, the solution process corresponding to all load cases easily
converge to the exact solution and the number of iterations is almost the samewhether the load
cases are solved together or separately. However, for anisotropic hyper-elastic materials, even in
the cases without any measurement error, the solution corresponding to some load cases may
diverge. Since the number of constants to be determined in the anisotropic case is more than
the isotropic case, more load cases are needed to find the unknown constants of hyper-elastic
anisotropic materials.
For both isotropic and anisotropic hyper-elastic materials, when the measurement error in-

creases, the number of iterations and the error in the obtained results increase too. Further, it
can be seen that usually the number of iterations is smaller for the cases that the problem is
solved by using the measured data from 2 or 3 load cases and the results obtained are more
accurate.
In this research, the preferred direction vectors are considered to be known. However, if one

considers these parameters as unknowns, the number of the unknowns of the inverse problem
increases. This will complicate the problem more and more, and the need for conducting more
experiments increases as well.



Determination of material parameters of isotropic and anisotropic... 909

References

1. Abyaneh M.H., Wildman R.D., Ashcroft I.A., Ruiz P.D., 2013, A hybrid approach to
determining cornea mechanical properties in-vivo using a combination of nano-indentation and
inverse finite element analysis, Journal of the Mechanical Behavior of Biomedical Materials, 27,
239-248

2. AhnB., Kim Y., Kim J., 2008, Biomechanical characterizationwith inverse FEmodel parameter
estimation: Macro and Micro applications, International Conference on Control, Automation and
Systems, 1769-1772

3. BakerM., ShrotA., 2013, Inverse parameter identificationwith finite element simulation using
knowledge-based descriptors,Computational Materials Science, 69, 128-136

4. Balaraman K., Mukherjee S., Chawla A., Malhotra R., 2005, Inverse finite element cha-
racterization of soft tissue using impact experiments and Taguchi method, SAE Internationals,
2005-01-4044

5. Czabanowski R., 2010, Experimental identification of hyperelastic material parameters for cal-
culations by the finite elementmethod, Journal of KONES, 17, 1, 87-92

6. Fung Y.C., 1993,Biomechanics: Mechanical Properties of Living Tissues, Springer, NewYork

7. Gasser T.C., Ogden R.W., Holzapfel G.A., 2006, Hyperelastic modeling of arterial layers
with distributed collagen fiber orientations, Journal of Royal Society Interface, 3, 6, 15-35

8. Hartman S., 2001, Numerical studies on the identification of the material parameters of Rivlin’s
hyper-elasticity using tension-torsion tests,Acta Mechanica, 148, 129-155

9. Hebden J.C., Price B.D., Gibson A.P., Royle G., 2006, A soft deformable tissue equivalent
phantom for diffuse optical tomography,Physics in Medicine and Biology, 51, 21, 5581-5590

10. Hematiyan M.R., Khosravifard A., Siah Y.C., Tan C.L., 2012, Identification of material
parameters of two-dimensional anisotropic bodies using an inversemulti loading boundary element
technique,Computer Modeling in Engineering and Science, 87, 1, 55-76

11. Holzapfel G.A., Gasser T.C., 2000, A new constitutive framework for arterial wall mechanics
and a comparative study of material models, Journal of Elasticity, 61, 1-48

12. Holzapfel G.A., 2000,Nonlinear Solid Mechanics, Technical University Graz, Austria

13. HolzapfelG.A.,OgdenR.W., 2010,Constitutivemodelingofarteries,Journal ofRoyal Society,
466, 2118, 1551-1597

14. HuT., Desai J.P., 2004,Characterization of Soft-TissueMaterial Properties: Large Deformation
Analysis, Drexel University, Philadelphia

15. Krouskop T.A., Wheeler T.M., Kallel F., Garra B.S. Hall T., 1998, Elastic moduli of
breast and prostate tissues under compression,Ultrasonic Imaging, 20, 260-274

16. Liu G.R., Han X., 2003,Computational Inverse Techniques in Nondestructive Evaluation, CRC
Press, Boca Raton, London, NewYork,Washington, D.C.

17. Miller K., Chinzei K., 1997, Constitutive modeling of brain tissue: experiment and theory,
Journal of Biomechanics, 30, 1115-1121

18. Mehrabian H., 2008, Soft tissue hyper-elastic parameter reconstruction for breast cancer asses-
sment, Mastered Engineering Science Thesis, University of Toronto

19. Mesa-MúneraE.,Raḿirez-SalazarJ.F.,BoulangerP.,BranchJ.W., 2012, Inverse-FEM
characterization of a brain tissue phantom to simulate compression and indentation, Ingeniería y
Cienc, 8, 11-36

20. Moulton M.J., Creswell L.L., Actis R.L., Myers K.W., Vannier M.W., Szabo B.A.,
Pasque M.K., 1995, An inverse approach to determiningmyocardialmaterial properties, Journal
of Biomechanics, 28, 8, 935-948



910 M. Hajhashemkhani, M.R. Hematiyan

21. Ogden R.W., Saccomandi G., Sgura I., 2004, Fitting hyper-elastic models to experimental
data,Computational Mechanics, 34, 6, 484-502

22. PamidiM.R.,Advani S.H., 1978,Nonlinear constitutive relations for humanbrain tissue,Journal
of Biomechanical Engineering, 100, 1, 44-48

23. RauchsG.,BardonJ.,GeorgesD., 2010, Identification of thematerial parameters of a viscous
hyper-elastic constitutive law from spherical indentation tests of rubber and validation by tensile
tests,Mechanic of Materials, 42, 961-973

24. Rivlin R.S., Saunder D.W., 1951, Large elastic deformations of isotropic materials VII. Expe-
riments on the deformation of rubber, Journal of Royal Society, 243, 865, 251-288

25. Schriefl A.J., Reinisch A.J., Sankaran S., Pierce D.M., Holzapfel G.A., 2012, Quan-
titative assessment of collagen fibre orientations from two-dimensional images of soft biological
tissues, Journal of Royal Society, 9, 76, 3081-3093

26. Unnikrishnan V.U., Unnikrishnan G.U., Reddy J.N., 2012, Biomechanics of breast tumor:
effect of collagen and tissue density, International Journal of Mechanics and Materials in Design,
8, 3, 257-267

Manuscript received August 25, 2014; accepted for print May 8, 2015