ORIGINAL PAPER
E. Mazza Æ O. Papes Æ M. B. Rubin Æ S. R. Bodner
N. S. Binur
Nonlinear elastic-viscoplastic constitutive equations for aging facial
tissues
Received: 17 November 2004 / Accepted: 19 April 2005 / Published online: 12 August 2005
Springer-Verlag Berlin Heidelberg 2005
Abstract This paper reports on the initial stages of a
project to simulate the nonlinear mechanical behavior of
an aging human face. A cross-section of the facial
structure is considered to consist of a multilayered
composite of tissues with differing mechanical behavior.
The constitutive properties of these tissues are incorpo-
rated into a finite element model of the three-dimen-
sional facial geometry. Relatively short time (elastic-
viscoplastic) behavior is governed by equations previ-
ously developed which are consistent with mechanical
tests. The long time response is controlled by the aging
elastic components of the tissues. An aging function is
introduced which, in a simplified manner, captures the
observed loss of stiffness of these aging elastic compo-
nents due to the history of straining as well as other
physiological and environmental influences. Calcula-
tions have been performed for 30 years of exposure to
gravitational forces. Progressive gravimetric soft tissue
descent is simulated, which is regarded as the main
indication of facial aging. Results are presented for the
deformations and stress distributions in the layers of the
soft tissues.
1 Introduction
Biomedical research is being influenced by mechanics in
a number of important ways. Examples are the use of
computer methods to simulate trauma (Snedeker et al.
2002), and surgical planning and surgical training using
virtual reality (Brett et al. 1995; Burdea 1996; Koch et al.
1996; Avis 2000; Szekely 2003). In most cases the goal of
simulations rela ted to the fa ce has been to produce
realistic animations of facial expressions (Parke 1982;
Lee et al. 1995; Koch et al. 1996; Zhang et al. 2004).
The objective of the present paper is to focus atten-
tion on the important role of realistic modeling of the
mechanical response of facial tissues to loads. Previous
research (Har-Shai et al. 1996, 199 7) has considered
mechanical experiments to study the viscoplastic prop-
erties of the skin and the underlying supportive tissue
SMAS (superficial musculoaponeurotic system).
Mechanical constitutive equations were developed in a
simple one-dimensional form (Rubin et al. 1998) and in
a more general three-dimensional form (Rubin and
Bodner 2002). More specifically, the constitutive equa-
tions discussed there are in the class of elastic-visco-
plastic phenomenological equations which can model
the time-dependent, nonreversible mat erial response
observed in most solids. In elastic-viscoplastic theory,
which is akin to nonlinear visoelasticity, the stress re-
sponse is determined by history-dependent variables
which are introduced through evolution equations ra-
ther than from hereditary integrals commonly used in
viscoelasticity.
With regard to the modeling of facial tissues, the
history-dependent variables and their associated evolu-
tion equations are introduced to model the main mac-
roscopic response caused by the microscopic
morphology of the tissues and relevant biological pro-
cesses. In general, the equations and the associated
material parameters are determined by comparing
predictions of the theory with macroscopic experimen-
tal data. A direct connection between microscopic
E. Mazza (&) Æ O. Papes
Institute of Mechanical Systems ETH,
8092 Zurich, Switzerland
E-mail: edoardo.mazza@imes.mavt.ethz.ch
Tel.: +41-1-6325574
Fax: +41-1-6321145
M. B. Rubin Æ S. R. Bodner
Faculty of Mechanical Engineering,
Technion—Israel Institute of Technology,
32000 Haifa, Israel
N. S. Binur
The Cosmetic Surgery Center of SE Texas Port Arthur,
Texas, 77642, USA
Biomechan Model Mechanobiol (2005) 4: 178–189
DOI 10.1007/s10237-005-0074-y
biological processes and the material parameters of this
constitutive model remains illusive at this time. Never-
theless, a reasonable understanding of the mechanical
behavior of the composite tissue structure should be
helpful to guide new developments in clinical practice.
Regarding the aging of facial tissues, LaTrenta (2004,
pp 46–47) states that: ‘‘The most commonly held theory
is that facial aging is the result of progressive gravimetric
soft tissue descent . Over time, the soft tissues of the face
simply sag off the bones of the face, forming the dis-
tinctive wrinkles, furrows, folds, and eventual tissue
redundancy of the aged face. Gravimetric soft tissue
descent is complex, how ever, and enc ompasses several
distinct processes. One of the most important processes
is actinic damage or solar elastosis.... Wrinkles become
apparent in a woman’s skin in her mid-30s as estrogen
levels begin to decline from their peak. The dermis
begins to lose collagen and elastin... Fat, unlike muscle,
is supported solely by facial ligaments. After years of
being pulled and stretched, these facial ligaments never
regain their tautness.’’
In the proposed constitutive model, aging is charac-
terized by a reduction of stiffness of facial tissues. This
modeling approach is justified at the histological level by
tissue degradation processes. The dermis becomes atro-
phied during aging, with a reduction of the volume
fraction of glycosaminoglicans (specifica lly, the hyal-
uronic acid) and collagen fibers of types I, III and VII
(Craven et al. 1997; Fleischmajer et al. 1972). In par-
ticular, the dermal elastic fiber network (oxitalan fibers)
decreases significantly with age (Cotta-Pereira et al.
1978). In the face, these processes are accelerated by
damage due to sunlight exposure (photoaging) and are
complemented by solar elastosis, an accumulation of
truncated, disorganized elastic fibers in the dermis
(Craven et al. 1997). Progressive atrophy of superficial
fat occurs in the face at distinct locations (Donofrio
2000) which might contribute to the reduction of stiff-
ness of the SMAS (Har-Shai et al. 1998).
Age-dependent gravimetric descent is also related to
the increase in volume of deep fat and the accumulation
of fat around the eyes, in the cheeks and under the chin,
as well as to progressive lengthening of the musculature,
laxity of the ligaments and skeletal resorption (LaTrenta
2004). These processes were not considered in detail in
the present model. One other typical indication of face
aging that could not be simulated with the present model
is the formation of deep wrinkles and skin folds. Wrin-
kles evolve from mimetic lines such as those around the
eyes and the mouth. This proce ss is related to localized
histological modifications, in particular actinic elastosis
and disappearance of oxitalan fibers (Contet-Audo n-
neau et al. 1999), which might be influenced by the
concentration of stress and strain at these locations.
In this paper, the gravimetric descent of the facial
tissues is modeled by implementing the nonlinear con-
stitutive equations into the commercial finite element
computer code ABAQU S (Hibbit et al. 2002). Specifi-
cally, the constitutive equations in Rubin et al. (1998)
and Rubin and Bodner (2002), which were developed for
a relatively short time response to cyclic loading and
relaxation tests, are generalized to include an aging
function which captures the main effect of tissue deg-
radation (loss of stiffness) exhibited in long time
behavior. In particular, an aging (damage) quantity x is
proposed as a nonlinear function of an auxiliary aging
parameter a (see Fig. 1), which itself is determined by
integrating a rather simple evolution equation for its
time rate of change.
The three-dimensional finite element model of the
face is based on a range laser scan (Vannier et al. 1991)
of a young man. Kinematic boundary conditions are
defined from the description of face anatomy by Barton
(2001) and LaTrenta (2004). The facial tissue is modeled
as having four layers (skin, SMAS and superficial fat,
deep fat and mucosa, see Fig. 2). In the present exercise,
calculations have been performed for 30 years of expo-
sure to gravitational forces from the onset of aging (age
assumed to be about 30). This exposure period would
correspond to about 45 years of living (assuming the
face is in an erect position for 16 h per day). Such
exposure would cause progressive gravimetric soft tissue
Fig. 1 Plot of the assumed stiffness reduction factor due to aging as
a function of the aging parameter a
Fig. 2 Simple model of the tissue layers in the face
179
descent, which is regarded as the main indication of
facial aging (LaTrenta 2004, p. 46). Results are pre-
sented for the deformations and stress dis tributions in
the layers of the soft tissues.
An outline of the paper is as follows. Section 2
describes the nonlinear constitutive equations and the
aging function which are used to model each of the facial
tissues. Sec tion 3 discusses the procedure for determin-
ing the material constants so that the model matches
available mechanical experimental data. Section 4 briefly
discusses the numerical model of the face. Section 5
describes the main resul ts and Sect. 6 presents conclu-
sions and possible directions for future res earch.
2 Constitutive equations
Rubin and Bodner (2002) developed nonlinear three-
dimensional constitutive equations for facial tissues
which are valid for arbitrarily large defo rmations and
which produce reasonable agreement with the experi-
mental data of Har-Shai et al. (1996). In that work, the
tissue was modeled as a composite material with a fully
elastic compone nt and a dissipative component which
contains both elastic and viscous elements (similar to a
Maxwell model in viscoelasticity theory). The specific
(per unit mass) strain energy function w was specified in
the form
q
0
w ¼
l
0
2q
expðqgÞ1½; ð1Þ
where q
0
is the constant reference mass density, l
0
is a
constant shear modulus and q controls nonlinearity.
Moreover, the function g was specified in an additive
form
g ¼ g
1
þ g
2
þ g
3
þ g
4
; ð2Þ
with {g
1
, g
2
, g
3
} characterizing the response of the fully
elastic component and with g
4
associated with the dis-
sipative compon ent. Specifically, g
1
characterizes the
fully elastic response to dilatation (volumetric changes),
g
2
characterizes the fully elastic response to distortional
deformations and g
3
characterizes the fully elastic
response to extension of specific fibers. The function g
4
characterizes the respon se to elastic distortional defor-
mation of the dissipative component. It is noted that
the dissipative component is modeled using elastic-vi-
scoplastic constitutive equations whic h for low-stress
levels can produce nearly nondissipative elastic re-
sponse and for high-stress levels can produce dissipative
viscoplastic response. Consequently, the word ‘‘dissi-
pative’’ is used to describe the dissipative component
because that component has the potential to exhibit
dissipation even though it may respond essentially
elastically in some ranges of stress. For all stress levels,
the stress tensor is a function of elastic deform ation
measures associated with both the fully elastic and the
dissipative components.
Here, the effects of specific fibers are neglected
(g
3
=0) due to lack of experimental data and the func-
tion g
2
is modified (relative to Rubin and Bodner 2002)
to model the reduction in elastic stiffness caused by
distortional deformations of the elastic component using
an aging term x. Due to this modification, the elastic
component of the present model is referred to as an
‘‘aging elastic’’ component rather than a ‘‘fully elastic’’
component. Except for this modification, most of the
details of the following developments can be found in
Rubin and Bodner (2002). Specifically, the functions {g
1
,
g
2
, g
3
, g
4
} are taken in the forms
g
1
¼ 2m
1
ðJ 1Þ1nðJÞ½; g
2
¼ð1 xÞm
2
ðb
1
3Þ;
g
3
¼ 0; g
4
¼ a
1
3; ð3Þ
where J is the dilatation, b
1
and a
1
are measures of the
elastic distortions of the elastic and dissipative compo-
nents, and m
1
, m
2
are material constants. Details of the
independent variables in the theory will be presented in
Eqs. 7–19 after presenting the constitutive equations for
the stresses. Next, the constitutive equation for the
Cauchy stress T is hype relastic with T obtained by a
derivative of w and T being separated ad ditively into
three components
T ¼ T
ð1Þ
þ T
ð2Þ
þ T
ð4Þ
; T
ð1Þ
¼m
1
l
1
J
1
I; ð4a; bÞ
T
ð2Þ
¼ð1 xÞm
2
lJ
1
B0
1
3
ðB0IÞI
;
T
ð4Þ
¼ lJ
1
B0
de
1
3
ðB0
de
IÞI
; (4c, d)
where l is a nonlinear shear modulus, and B¢ and B¢
de
are tensorial measures of elastic distortion of the elastic
and dissipative components. In these equations, T
(1)
characterizes the elastic response to dilatation, T
(2)
characterizes the aging elastic response to distor-
tional deformations, and T
(4)
characterizes the dissipa-
tive response to distortional deformat ions. Also,
AÆB=tr(AB
T
) denotes the inner product between two
second-order tensors A and B.
From Eq. 4a it can be seen that the pressure p becomes
p ¼
1
3
T I ¼ m
1
l
1
J
1
; ð5Þ
where the bulk modulus m
1
l is defined in terms of the
nonlinear shear modulus l
l ¼ l
0
expðqgÞ: ð6Þ
The total dilation J is determined by the standard
evolution equation
_
J ¼ JD I; ð7Þ
where a superposed dot denotes material time differen-
tiation, and D is the symmetric part of the velocity
gradient L=ƶv/ƶx
180
D ¼
1
2
ðL þL
T
Þ: ð8Þ
In Eq. 4c, B¢ is a unimodular second order tensor
(Flory 1961)
det ðB
0
Þ¼1; ð9Þ
which is a measure of total distortional deformation,
and is determined by the evolution equation
_
B
0
¼ LB0þB
0
L
T
2
3
ðD IÞB
0
: ð10Þ
This equation is purely kinematical and can be
obtained by differentiating the expression
B
0
¼ J
2=3
B ð11Þ
where B is the left Cauchy-Green deformation tensor.
Also, the scalar b
1
in (2.3) is a measure of total distor-
tional deformation given by
b
1
¼ B
0
I: ð12Þ
In Eq. 4d, the dissipative component is modeled using
equations similar to those for viscoplasticity. Speci fi-
cally, B¢
de
is a unimodular second order tensor
detðB
0
de
Þ¼1; ð13Þ
which is a measure of the elastic distortional deforma-
tion associated with the dissipative compon ent, and is
determined by the evolution equation
_
B
0
de
¼ LB
0
de
þ B
0
de
L
T
2
3
ðD IÞB
0
de
CA
d
ð14Þ
Also, the scalar a
1
in Eq. 3 is given by
a
1
¼ B
0
de
I ð15Þ
Equation 14 is a generalization of Eq.10 with the rate
of inelastic deformation being controlled by the term C
Ad. In particul ar, when C A
d
vanishes, Eq. 14 has the
same form as Eq. 10 and the dissipative component re-
sponds elastically. Moreover, the tensor A
d
controls the
direction of inelastic deformation rate and is specified by
the form
A
d
¼ B
0
de
3
B
0
de
1
I
I ð16Þ
which ensures that B¢
de
remains unimodular Eq. 13.
The scalar function C in the evolution Eq. 14 requires
a constitutive equation which is specified by
C ¼ C
1
þ C
2
_
e½exp
1
2
b
b
de
2n
"#
; ð17Þ
similar to that of Bodner and Partom (1975). In
this equation, {C
1
, C
2
, n} are non-negative material
constants, b
de
is the magnitude of the deviatoric tensor
B
00
de
b
de
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
3
2
B
00
de
B
00
de
r
; B
00
de
¼ B
0
de
1
3
ðB
0
de
IÞI; ð18Þ
e is the equivalent total distortional strai n determined by
the evolution equation
_
e ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2
3
D
0
D
0
r
; D
0
¼ D
1
3
ðD IÞI; ð19Þ
and b is a measure of ha rdening of the dissipative
component characterized by the evolution equation
_
b ¼
r
1
r
3
þ r
2
_
e
r
3
þ
_
e
Cb
de
r
4
b
r
5
: ð20Þ
with {r
1
,r
2
,r
3
,r
4
,r
5
} being material constants. Using the
discussion in Rubin (1994), it can be shown that these
constitutive equations are properly invarian t under
superposed rigid body motions. Also, it is noted that the
term C
2
_
e in Eq. 17 was introduced in Rubin and Bodner
(2002) to cause the constitutive response to be consistent
with the observations on soft tissues which indicate that
hysteresis loops to the same stress levels are nearly
independent of strai n rate over a wide range of rates
(Fung 1993, p. 281).
For small values of b
de
relative to b, it follows from
Eq. 17 that the exponential term causes the rate of
inelastic deformation to be vanishingly small so the
dissipative component responds essentially elastically.
On the other hand, when b
de
attains values on the order
of b, then the magnitude of the inelastic response is
controlled by the constants {C
1
, C
2
}. The sharpness of
the transition between elastic and dissipative response is
controlled by the constant n.
The evolution Eq. 20 for the hardening parameter b
models two main effects. The first term on the right-
hand-side of Eq. 20 causes an increase in b due to
inelastic deformation rate. The second term on the right-
hand-side of Eq. 20 causes b to recover to the value zero.
This overall response attempts to model effects of fluid
flow in the tissue, the fluid being expelled when the tissue
is deformed (hardening) and the fluid being imbi bed by
the tissue over time as the tissue returns to its nearly
unstressed state (recovery of hardening). Specifically, the
constant r
1
controls the rate of hardening during relax-
ation tests ð
_
e ¼ 0Þ; the constant r
2
controls the rate of
hardening during loading (large values of
_
eÞ; and the
constant r
3
controls the value of strain rate
_
e associated
with the transition be tween these two responses. In this
regard, it shou ld be emphasized that this particular
functional form for hardening is quite simplistic and can
be modified when additional experimental data is
available. Also, the constants r
4
,r
5
control the rate and
shape of recovery of hardening.
The functional dependence on the total strain rate
_
e in
Eq. 17 and Eq. 20 attempts to capture differences in the
observed responses to loading and relaxation tests .
Specifically, for loading with large values of
_
e (i.e.
C
2
_
e C
1
and
_
e r
3
Þ the evolution Eq.14 characterizes
nearly rate-independent response. On the other hand,
181
during relaxation tests (with
_
e ¼ 0Þ the evolution Eq.14
characterizes viscoplastic rate-dependent response.
Moreover, the transitions between these two types of
response depend mainly on the constants {C
2
,r
3
}.
The main modification in Eq. 3, relative to the
constitutive equations of Rubin and Bodner (2002), is
the presence of the aging (damage) term x which causes
the stiffness of the elastic response to distortional
deformation (see g
2
in Eq. 3) to decrease with increas-
ing x. Consistent with information in the medical lit-
erature, damage due to aging of the tissues is
considered to be nonreversible. Most models for dam-
age of materials propose an evolution equation for the
time rate of change of a damage variable like x which
is a highly nonlinear function of x and stress. Here, an
alternative procedure is proposed which improves the
stability of the numerical integr ation of the aging
model. Specifically, the aging quantity x is taken to be
a function of an auxiliary aging parameter a which
produces a generic S-like curve (see Fig. 1)
x ¼ xðaÞ¼ð1 a
1
Þ
a
2
ð3a=2Þ
a
3
1 a
2
þ a
2
ð3a=2Þ
a
3
;
0 a
1
1; 0 a
2
1; a
3
0;
ð21Þ
where {a
1
,a
2
,a
3
} are material constants and a is
determined by an evolution equation. This functional
form has been chosen so that x is bounded by zero and
the value (1–a
1
)
0 x 1 a
1
; xð0Þ¼0; xð1Þ ¼ 1 a
1
: ð22Þ
Also, it can be shown that
dx
da
¼
3ð1 a
1
Þa
2
a
3
ð3a=2Þ
a
3
1
2½1 a
2
þ a
2
ð3a=2Þ
a
3
2
0; ð23Þ
which indicates that x is a monotonically increasing
function of a. Next, the aging parameter a is determ ined
by the relatively simple evolution equation
_
a ¼ C
3
e þC
4
; ð24Þ
where {C
3
, C
4
} are additional material constants, with
C
3
controlling the dependence on strain and C
4
con-
trolling the dependence on other physiological and
environmental effects of aging. Also, the factor 3/2 in
Eq. 21 was chosen for convenience and could be ad-
justed by changing the values of {C
3
, C
4
}.
With regard to more standard formulations of dam-
age (e.g. Bodner and Chan 1986) it is noted that Eqs. 21,
23 and 24 could be combined to obtain an evolution
equation directly for the damage parameter x which is
independent of the auxiliary parameter a. However, the
resulting equation would be a highly nonlinear function
of x which could cause difficulties in numerical inte-
gration. In contrast, the procedure used here embeds
most of the nonlinearity in the functional form Eq. 21
for x(a) and leaves a rather simple evolution equation
Eq. 24 for a. Thus, this procedure of specifying x as a
function of an auxiliary parameter a may be useful for
other constitutive equations that include continuum
damage parameters.
Next, using the above constitutive equations and the
conservation of mass
qJ ¼ q
0
; ð25Þ
which gives an expression for the current mass density q,
it can be shown that the rate of material dissipation D
D ¼ T D q
_
w 0; ð26Þ
requires
D ¼ T D q
_
w
¼
1
2
lJ
1
CA
d
I þ m
2
ðb
1
3Þ
_
x½0; ð27Þ
which is satisfied for all processes since {l
0
,m
2
, C} are
nonnegative.
The material constants {l
0
,q,m
1
,m
2
} control the
elastic response of the tissue. More specifically, q con-
trols nonlinear elastic effects through the strength of the
exponential function and it has insignificant influence
for small deformations (with l l
0
). Consequently, the
remaining elastic constants can be identified by consid-
ering small deformations. The constant l
0
controls the
elastic shear modulus of the dissipative component and
the elastic moduli of the other components have been
normalized by l
0
. Thus, the small deformation bulk
modulus k
1
of the tissue is given by
k
1
¼ m
1
l
0
: ð28Þ
The small deformation shear modulus of the aging
elastic component is given by (1x)m
2
l
0
. Thus, the
aging function x tends to reduce the magnitude of this
shear modulus as x and a increase (see Fig. 1). The
material constants {a
1
,a
2
,a
3
} control the shape of the
aging function Eq. 21. From Eq. 22 and Eq. 23 it is
clear that x smoothly changes from the value 0 to its
maximum value 1a
1
as the aging parameter a in-
creases. This means that the small deformation shear
modulus of the fully aged elastic component becomes
a
1
m
2
l
0
, which is reduced from the young tissue value of
m
2
l
0
by the constant a
1
. The constant a
2
controls the
extent of this transition when a=2/3, since
x ¼ð1 a
1
Þa
2
for a ¼
2
3
: ð29Þ
Also, the transition becomes sharper as the value of
a
3
is increased. As noted earlier, the value of a=2/3 in
Eq. 21 was chosen for convenience but could be adjusted
by changing the factor 3/2 in Eq. 21. Furthermore, the
evolution Eq. 24 for the aging parameter a models the
effect of straining on the inhomogeneous evolution of
aging by the constan t C
3
. It is well known that exposure
to sun causes aging of skin but quantif ying the cumu-
lative effect of a time-dependent history of sun exposure
182
