ICABB-2010, Venice, Italy
October 14-16, 2010
Modeling and parameter estimation of rheological objects for
simultaneous reproduction of force and deformation
Zhongkui Wang and Shinichi Hirai
Department of Robotics, Ritsumeikan University
Noji-higashi 1-1-1, Kusatsu, 525-8577, Japan
Abstract— Many deformable objects in our living life demon- application of robotics and automation in food industry is a
strate rheological behaviors, such as human organs and tissues, lack of understanding of the food product properties as an
clays and various food products. Rheological objects include “engineering” material for handling operations . To date,
both elastic and plastic properties. Due to the presence of
residual (permanent) deformation, it is difﬁcult to model rheo- modeling of 2D/3D rheological objects have not been well
logical objects, especially to reproduce both force and residual developed and an effective method for determining physical
deformation simultaneously. In this paper, a series of physical properties of rheological object as an “engineering” material
models was investigated for simulating rheological behaviors. has not yet been established.
Generalized formulations of the constitutive laws were derived Generally, modeling a rheological object is more difﬁcult
for serial and parallel physical models, respectively. We found
that the serial models are appropriate for formulating strain, than modeling an elastic or plastic object due to the presence
whereas the parallel models allow a convenient calculation of of residual deformation. Early work on the modeling of
stress. Analytical expressions of force and residual deformation rheological objects dates back to Terzopoulos et al. , who
were then derived for generalized parallel models. Theoretical have proposed several physical models to describe inelastic
discussions suggested the difﬁculty to reproduce both force objects and a Burgers model was introduced to describe
and deformation simultaneously using linear physical models.
A 2D FE (ﬁnite element) model was then formulated and rheological behaviors. However, it is only a conceptual
an efﬁcient method for estimating physical parameters were description and they did not give any simulation results and
proposed by taking the advantages of analytical force expres- information of parameter determination. A particle-based (or
sions. Experimental results with commercial clay and Japanese mass-spring-damper) model was employed to model a food
sweets material were presented to validate our modeling and dough and a method for calibrating its physical parameters
parameter estimation methods. A dual-moduli viscous element
was also introduced to improve our FE model for reproducing was proposed based on a genetic algorithm (GA) optimiza-
rheological force and deformation simultaneously. tion , . The particle-based model has an advantage
Index Terms— Modeling; Parameter estimation; Rheological of less computation costs , but the formulation was not
object; Finite element method; Simulation. based on continuum mechanics and the simulation accuracy
is quite limited. A two-layered Maxwell model  and
I. I NTRODUCTION Fung’s viscoelastic model  has been used respectively
Modeling and simulation of soft deformable objects has to reproduce the force response of a “Norimaki-sushi” when
been studied for over 20 years ,  and has played an grasped by a robot hand. Good approximations of force
important role in many ﬁelds, including computer aided behaviors were obtained. Unfortunately, both models are 1D
surgery, food automation, and robotics. Deformable objects cases. In addition, the ISU exoskeleton technique has been
can be roughly divided into three categories based on their used in modeling clay to simulate an interaction between
deformation behavior: elastic object, in which the defor- virtual clay and a human ﬁnger .
mation is completely reversible; plastic object, in which Interestingly, above-mentioned work has focused on either
the deformation is completely maintained; and rheological reproduction of deformation alone ,  or reproduction
object, in which the deformation is partly reversible. Previous of force alone , . None of them considered both.
work on modeling of deformable objects has mainly focused However, in some application, e.g., surgical simulation with
on elastic and visco-elastic objects, especially in surgical haptic feedback or virtual manipulation of food products,
applications since most biological organs and tissues seem accurate results of both force and deformation are neces-
to be recoverable after loading-unloading operations. Some sary. Both force and deformation behaviors of organs and
organs and tissues, however, may fail to completely recover tissues were considered simultaneously in an ex vivo and
from the deformation. Hrapko et al. found that porcine an in vivo property characterization of porcine livers ,
brain tissue did not recover completely after a loading- . However to our knowledge, in current surgical related
unloading cycle . In vivo experimental results showed that applications, organs and tissues were mostly supposed to be
residual deformation may also present in human liver . In elastic or visco-elastic and no residual deformation was taken
addition, many other deformable objects, such as clay and into consideration.
various food products, demonstrate rheological behaviors. In our previous work, we have developed FE models
Chua et al. stated that the most critical barrier against the for simulating various deformable objects . A three-
E E1 c1
element  and a four-element  models were employed
respectively to describe rheological behaviors. A method c1
for identifying physical parameters was proposed based on E1 E1 c1
iterative optimization . We found that the four-element c1 E2 c2
model was more appropriate than the three-element model E1 c1
for reproducing rheological forces. In addition, experimental E2 c2
results with commercial clay indicated that a contradiction c3 c3
between the reproductions of force and residual deformation
(a) Serial models (b) Parallel models
may exist due to the linearity of the physical model.
In this paper, we summarized the physical models which Fig. 2. 2 groups of rheological models: (a) serial, and (b) parallel models.
could be used to represent rheological behaviors. According
to the conﬁguration among elements, we divided the phys-
ical models into two groups: serial and parallel models. A series. Note that the deformation generated in an elastic or a
general constitutive law for each conﬁguration was formu- Kelvin element is completely recoverable. Therefore, a serial
lated. To analyze the reason of the contradiction between rheological model must include a viscous element connected
reproductions of force and residual deformation, analytical in series, which causes the residual (permanent) deformation.
expressions of force and residual deformation were derived According to the presence of elastic element, serial models
and discussions were done. Two-dimensional FE dynamic can be divided into two types, as shown in Fig. 3. Let us
model was then given with simpler formulation comparing take the type 1 model as an example to show the derivation
with the previous one. Based on the analytical expressions procedure of the constitutive law.
of rheological forces, a more efﬁcient method for estimating Note that the constitutive law of four basic elements were
the physical parameters was proposed. Experimental results formulated as:
with commercial clay and Japanese sweets material were
presented to validate our discussions and methods. Finally, a Elastic element : σ = Eε,
dual-moduli viscous element was introduced to improve our Viscous element : ε
σ = c˙ ,
current model for simultaneously reproducing both rheolog- E (1)
ical force and deformation. Maxwell element : σ + σ = Eε,
Kelvin element : ε
σ = E ε + c˙ ,
II. S UMMARIZATION OF P HYSICAL M ODELS
Physical models are often employed to describe de- where vector σ and ε are the stress and strain generated
formable materials and objects, e.g., an elastic element (Fig. in basic elements. Constants E and c are the Young’s mod-
1a) and a viscous element (Fig. 1b) represent ideal elastic ulus and viscous modulus of elastic and viscous elements,
and viscous material, respectively. An elastic and a viscous respectively.
elements connected in series was called a Maxwell element Let εi and εn+1 be the strain at the i-th Kelvin element
(Fig. 1c), which denotes a simplest rheological material. An and the (n + 1)-th viscous element, respectively. Let Ei and
elastic and a viscous elements connected in parallel was ci be the Young’s modulus and viscous modulus of the i-th
called a Kelvin (sometimes called Kelvin-Voigt) element elastic and viscous elements, respectively. Due to the serial
(Fig. 1d), which denotes a visco-elastic material. We refer connections among these basic elements, the total stress at
the above four elements as basic elements. By connecting the serial model is equal to the stress at each basic element
some basic elements in different conﬁgurations, many phys- and the total strain at the serial model is equal to the sum
ical models could be obtained for simulating rheological of the strains at all basic elements. That is,
behaviors. We divided such models into two groups: serial
σ = Ei εi + ci εi ,
˙ 1 ≤ i ≤ n,
and parallel models (Fig. 2).
σ = cn+1 εn+1 ,
A. Generalized Serial Model n+1
A serial rheological model consists of numbers of Kelvin ε= ∑ εi .
elements and a viscous or a Maxwell element connected in
E1 E2 En
E c cn+1
(a) (b) c1 c2 cn
E (a) Type 1
E c E1 E2 En
c En+1 cn+1
(c) (d) c1 c2 cn
(b) Type 2
Fig. 1. The basic elements for describing deformable materials: (a) the
elastic; (b) the viscous; (c) the Maxwell; and (d) the Kelvin elements. Fig. 3. Generalized serial models: (a) type 1; and (b) type 2.
Taking Laplace transform of the above equations, we have is no constant term in the coefﬁcients of strain polynomial
(the left side of Eq. 12).
σ (s) = Ei εi (s) + ci sεi (s), 1 ≤ i ≤ n,
Following the same derivation procedure, we could obtain
σ (s) = cn+1 sεn+1 (s), the constitutive law of serial model of type 2 as follows:
ε (s) = ∑ εi (s). n+1
∑ Ai−1 ∂ t i
∑ Bs2 ∂ t j ,
Eliminating εi (s) from the above equations, we then have
1 1 1 1
ε (s) = ∑ s + ri ci + s cn+1 σ (s), (4)
, Bs2 = ∑+
, · · · , Bs2 =
i=1 En+1 i=1 ci cn+1 En+1 cn+1
where ri = Ei /ci . Let us deﬁne a polynomial as below: (14)
n Equation 13 suggested that the highest-order derivative of
∏(s + ri ) = An sn + An−1sn−1 + · · · + A1s + A0. (5) strain ε is equal to the highest-order of stress σ . Note that
i=1 the left side of Eq. 13 has the same form with the left side
The coefﬁcients of the above polynomial have the forms of: of Eq. 12.
n n n n
B. Generalized Parallel Models
An = 1, An−1 = ∑ ri , An−2 = ∑ ∑ ri r j , · · · , A0 = ∏ ri .
i=1 i=1 j=1 i=1 Two kinds of parallel rheological models were shown in
j=i Fig. 4. Due to the parallel connections among basic elements,
the total strain at the parallel model is equal to the strain at
Multiplying Eq. 4 by Eq. 5, we have
each basic element and the total stress at the parallel model
n n n
1 1 1 1 is equal to the sum of the stress at all basic elements. For
∏(s + ri )ε (s) = ∏(s + ri) ∑ s + ri ci +
parallel model of type 1, we therefore have
i=1 i=1 i=1
(s + r j )
n n s+r
j 1 Ei
= ∑∏ ci
σ (s). σi +
σi = Ei ε , 1 ≤ i ≤ n,
i=1 j=1 j=1
j=i σn+1 = cn+1 ε ,
We then ﬁnd the following equation: σ= ∑ σi .
∏ (s + r j ) = (s + r1) · · · (s + ri−1 )(s + ri+1) · · · (s + rn) Following the same derivation with serial models, we could
j=1 (8) end up with the constitutive law of parallel model of type 1
(Fig. 4a) as follows:
= sn−1 + Bi,1 sn−2 + · · · + Bi,n−2s + Bi,n−1,
∂ iσ n+1
where ∑ Ai ∂ t i = ∑ B p1 ∂ t j ,
n n n n i=0 j=1
Bi,1 = ∑ r j, Bi,2 = ∑ ∑ r j rk , · · · , Bi,n−1 = ∏ r j . (9) where
j=1 j=1 k=1 j=1
j=i j=i k=i j=i n
k= j Bn+1 = cn+1 , Bn = ∑ Ei + An−1cn+1 ,
Substituting Eqs. 5 and 8 into Eq. 7, we have the following i=1
Laplace transform equation: ··· , p1
B1 = ∑ Bi,n−1 Ei + A0cn+1 .
(An sn+1 + An−1sn + · · · + A0 s)ε (s) i=1
(10) Correspondingly, the constitutive law of parallel model of
= (Bs1 sn + Bs1 sn−1 + · · · + Bs1s + Bs1 )σ (s),
n n−1 1 0
type 2 (Fig. 4b) could be formulated as:
∂ iσ n
∑ ci + cn+1 , Bs1 = ∑ ci + cn+1 , · · · , Bs1 = cn+1 .
A0 ∑ Ai ∂ t i = ∑ B p2 ∂ t j ,
Applying the inverse Laplace transform to Eq. 10 yields the E1 c1 E1 c1
constitutive law of serial model of type 1 as follows: E2 c2 E2 c2
∂ iε n
∂ jσ En cn Ei ci
∑ Ai−1 ∂ t i = ∑ Bs1 ∂ t j .
j (12) cn+1 En M cn
(a) Type 1 (b) Type 2
Note that the highest-order derivative of strain ε is one order
larger than the highest-order of stress σ . In addition, there Fig. 4. Generalized parallel models: (a) type 1; and (b) type 2.
n n n
Bn = ∑ Ei , Bn−1 = ∑ Bi,1 Ei , · · · , B1 = ∑ Bi,n−1 Ei . (19)
p2 p2 p2
i=1 i=1 i=1
According to Eqs. 12, 13, 16, and 18, we found that the t
constitutive law of serial model of type 1 (Eq. 12) has
the identical form with the parallel model of type 1 (Eq. (a-1) Rheological force (a-2) Deformed shapes
(a) Experimental results with commercial available clay
16) except some coefﬁcients having different formulations.
Correspondingly, the constitutive laws of serial model of type
2 (Eq. 13) also has the same form with the parallel model
of type 2 (Eq. 18) by replacing the summation limit n + 1
in Eq. 13 by n. Note that same constitutive laws yield same
deformation behaviors. Therefore, for simulating a certain t
behavior, we could use either a serial model or a parallel
model. Actually, for any type of physical model, which (b-1) Rheological force (b-2) Deformed shapes
(b) Experimental results with Japanese sweets material
consists of any numbers of basic elements connected in any
conﬁguration, we are always able to ﬁnd a corresponding Fig. 5. Typical rheological behaviors of: (a) commercial available clay,
serial or parallel model which yields the same behaviors. and (b) Japanese sweets material.
This allows us to investigate only one representative model
instead of analyzing all models for understanding the ability
of physical models. In this paper, we investigated the parallel B. Analytical Expression of Rheological Force
models. In addition, according to Eq. 2, if the total stress We took the parallel model of type 1 as an example to
at the serial model was available, we could easily calculate show the derivation. In the pushing phase, we suppose the
the total strain by summing up the individual strain at each strain rate is constant, i.e., ε = p. Solving Eq. 15, we have the
element due to the independence among these strains. On the analytical force expression in the pushing phase as below:
other hand, Eq. 15 indicated that the convenient calculation n E
− ci t
of rheological force could be achieved by using the parallel σ (t) = ∑ ci p 1 − e i + cn+1 p, (0 ≤ t ≤ t p ). (20)
III. A NALYSIS OF PARALLEL M ODEL In the holding phase, solving Eq. (15) with ε = 0 and initial
A. Experimental Rheological Behaviors condition of σi (t p ), we can formulate the analytical force
expression in this phase as:
Typical rheological behaviors (force and deformed shapes)
n E E
of commercial available clay and Japanese sweets material − c i tp − c i (t−t p )
σ (t) = ∑ ci p 1 − e i e i , (t p ≤ t ≤ t p + th ).
were shown in Fig. 5. The sweets materials were provided i=1
by OIMATU, a sweets company in Kyoto. Two ﬂat squared (21)
objects made of the above-mentioned two materials were
pushed respectively by a motorized stage with a constant C. Analytical Expression of Residual Deformation
velocity during the pushing phase (0 ≤ t ≤ t p ). Before unload- After unloading, we intuitively considered to solve the
ing, the deformed shapes were maintained for a while, which constitutive law Eq. 16 with σ = 0 to formulate the strain
was called holding phase (t p ≤ t ≤ t p + th ). The rheological recovering curve. However, when the order of time derivative
forces (Fig. 5a-1 and 5b-1) were measured by a tactile of strain ε exceeds two, it becomes impossible to solve Eq.
sensor and deformed shapes (initial, held, and ﬁnal shapes 16 because we have no information about the initial condition
as shown Fig. 5a-2 and 5b-2) were recorded by a calibrated of strain derivative. We therefore turn our focusing to each
camera. During the pushing phase, the rheological forces viscous element which dominates the residual strain. Let
were increasing with the deformation increasing. During εiela (t) and εivis (t) be the strain at each elastic and viscous
the holding phase, however, the deformed shapes were kept element, respectively. Note that the stress at a Maxwell
unchanged while the rheological forces were decreasing element is equal to the stress at the elastic element and the
(force relaxation). After unloading, the rheological forces viscous element as well. Thus, total stress after unloading
went to zero and the deformed shapes were partially recov- could be formulated as:
ered. Figure 5 showed that rheological behaviors of different n+1 n
materials are quite different. Comparing with clay, the force σ (t) = ∑ σi (t) = ∑ ci εivis (t) + cn+1ε (t) = 0.
˙ ˙ (22)
relaxation behavior of sweets material was slower and the i=1 i=1
recovered deformation was smaller. Our target is to ﬁnd an Integrating the above equation from time t p + th to time
appropriate model to reproduce both rheological force and inﬁnity, we have
deformation behaviors simultaneously. Let us now investigate
n ∞ ∞
the ability of generalized parallel model for achieving this
∑ ci t p +th
εivis (t)dt + cn+1
t p +th
ε (t)dt = 0,
σ n +1 = c n +1 &
n σ (t p )
∑ ci εivis (∞) − εivis (t p + th ) + cn+1 ε (∞) − ε (t p + th ) = 0. Pushing
σ (t p + t h ) Pushing
σ (t p + t h )
σ (t p )
(24) ε (∞ ) tp tp+th ε (∞ )
It is important to note that the residual strain at every
viscous element in a parallel model should be the same and
equal to the total residual strain when time goes to inﬁnity,
i.e., ε1 (∞)=ε2 (∞)=· · · =εn (∞)=ε (∞), because all elastic
vis vis vis
(a) Parallel model type 1 (b) Parallel model type 2
elements completely recovered from the deformation. Thus,
Eq. 24 yields Fig. 6. Typical simulation results of rheological behaviors by using: (a)
5-element model, and (b) 2-layered Maxwell model.
ci εivis (t p + th ) cn+1 ε (t p + th )
ε (∞) = ∑ + . (25)
i=1 ∑n+1 c j
j=1 ∑n+1 c j
pushing phase for 5-element model (parallel type 1). This
Each viscous element has its own constitutive law as sudden drop was denoted by σ = cn+1 ε . Figure 6b showed
σi =ci εivis . Integrating this law through time 0 to time t p + th ,
˙ that the 2-layered Maxwell model (parallel type 2) results in
we have attenuated vibrations in both stress and strain after unloading
1 t p +th
εivis (t p + th) = σi (t)dt. (26) (Fig. 6b). Based on the above discussion, we could say that
ci 0 the parallel model has the ability of reproducing rheological
Substituting Eq. 26 into Eq. 25 and considering σ (t) = force behaviors. Our previous work  has shown good
∑n+1 σi (t), we ﬁnally end up with the expression of total
i=1 reproduction of rheological force for commercial clay. How-
residual strain: ever, we failed to reproduce the recovered shape.
1 t p +th Equation 27 implies that the residual strain is dominated
ε (∞) = n+1 σ (t)dt. (27) by the sum of viscous moduli ci . On the other hand, param-
∑i=1 ci 0
eters ci also strongly affect force amplitude as formulated in
This equation indicates that the ﬁnal residual deformation of Eqs. 20, 21, and 28. This caused a contradiction between the
a parallel rheological model depended on the sum of viscous reproduction of rheological force and recovered shape. For
moduli and the integration of force through the pushing and example, if we determine the parameters ci from force data,
holding phase. the sum of ci will therefore yields a certain recovered shape.
For the parallel model of type 2, we could obtain the We are unable to change this shape to our desired one. On
same formulations of force expression in the holding phase the contrary, if we ﬁrstly calculated the sum of ci based on
and the same formulation of ﬁnal residual deformation with Eq. 27, we have an upper limit (∑n+1 ci ) for each modulus
the summation limit n + 1 replaced by n in Eq. 27. The ci and we have to keep each ci under this limit during the
only difference of the parallel model of type 2 is the force reproduction of force behaviors. For some materials, we may
expression in the pushing phase, which is given by be able to achieve a good reproduction of force with ci
− ci t
under the limit. For most materials, however, this limit was
σ (t) = ∑ ci p 1 − e i , (0 ≤ t ≤ t p ). (28) always broken in order to well capture the force curve. The
above discussion suggests that parallel model do not have the
D. Discussions ability to reproduce both rheological force and deformation
Typical simulation results of rheological stress and strain simultaneously. Section 5 will give quantitative discussion
were shown in Fig. 6 by using a 5-element model (the last of this problem. The discussion in Section 2 indicated that
row of Fig. 2b) and a 2-layered Maxwell model (the middle we could always ﬁnd a corresponding parallel model for
row of Fig. 2b). According to Eqs. 20, 21, and 28, we found arbitrary physical model. We can therefore conclude that
that the stress curve was determined by viscous moduli ci physical models are not able to reproduce both rheological
and time coefﬁcients Ei /ci of exponential functions. The force and deformation simultaneously for most rheological
coefﬁcients Ei /ci determine the changing tendency in stress objects.
during the holding phase, as formulated in Eq. 21. In order
IV. FE DYNAMIC M ODEL AND PARAMETER E STIMATION
to obtain similar force relaxation curve as shown in Fig. 5, at
least two exponential terms are needed, one with large value A. Formulation of FE Dynamic Model
of Ei /ci and another one with small Ei /ci . The large Ei /ci The FE method has proven to be a powerful tool for
describes the rapid decreasing in force and the small one simulating complex behaviors of deformable objects. In FE
denotes the slow decreasing. For example, the stress curves in simulation, an object is described by a set of elements
Fig. 6 were both simulated with coefﬁcients of E1 /c1 = 0.005 (e.g., triangles in 2D and tetrahedrons in 3D). The dynamic
and E2 /c2 = 0.5. After determining coefﬁcients Ei /ci , Eq. behaviors of the object are then determined by analyzing the
20 could be used to determine the viscous moduli ci , which behaviors of individual elements. In this paper, we employ
dominated stress amplitude σ (t) in the pushing phase. Note a 2-layered Maxwell model as an example to present 2D FE
that there is a sudden drop in stress (Fig. 6a) at the end of dynamic model. We constructed a 2D object with triangle
n E E
− c i tp − c i (t−t p )
mesh (Fig. 5a-2 and 5b-2) and attached the two-layered F(t) = ∑ ci 1 − e i e i Mγ vPush , (t p ≤ t ≤ t p + th ),
Maxwell model onto each triangle. Let vectors uN and vN i=1
be a set of displacements and velocities of nodal points in (32)
2D mesh. Let vectors Frheo and Fext be a set of rheological where
forces and external forces on nodal points. Following the γ 1
Mγ = γλ Jλ + γμ Jμ = J + Jμ .
same formulation procedure presented in our previous work (1 + γ )(1 − 2γ ) λ 2(1 + γ )
, , , we could end up with an FE dynamic model
Vector vPush could be obtained from the simulation results
consists of following differential equations:
in the pushing phase with the known Poisson’s ratio γ .
uN = vN ,
˙ An objective function is deﬁned to be:
M˙ N − Aλ = −Frheo + Fext ,
v 2D N
−A vN = A (2ω vN + ω uN ),
˙ T 2 e(θ ) = ∑ fexp − fcal (t, θ ) 2 ,
i i (33)
F1 = − F1 + (λ1 Jλ + μ1 Jμ )vN ,
˙ ela ela
where N is the number of sampling points, vector fexp i
E2 and fcal (t, θ ) are experimental forces and calculated forces,
F2 = − F2 + (λ2 Jλ + μ2 Jμ )vN ,
˙ ela ela
respectively, at the i-th sampling time, and vector θ is the
parameters to be estimated. The optimization is terminated
where the last two equations were obtained from the ﬁrst when the tolerance on e(θ ) or the tolerance on θ less than
equation of Eq. 15 by changing 1D stress-strain relationship a threshold, 1 × 10−6. The basic estimation algorithm is as
to 2D force-displacement relationship. Two-dimensional rhe- follows:
ological force could be therefore formulated as
1) Initial setting of material parameters.
Frheo = F1 + F2 .
2D (30) 2) Using Eqs. 31 and 32 to calculate force.
3) Using Eq. 33 to calculate the objective function.
Let γ be the Poisson’s ratio and then variables λ1 , μ1 ,
4) If the terminate conditions are not satisﬁed, adjust
ela , and μ ela could be calculated as follows:
parameters and repeat steps 2–3.
Ei γ Ei
λiela = , μiela = , (i = 1, 2). This method is much faster than our previous simulation-
(1 + γ )(1 − 2γ ) 2(1 + γ ) based method. We could obtain a solution with only several
The detailed deﬁnition and description of other variables and seconds by using MATLAB optimization toolbox.
matrixes can be found in our previous work , , . In addition, extending Eq. 27 to 2D formulation, we have
Comparing with our previous FE dynamic model, this t p +th
model has a simpler formulation. This should thank to the Mγ uN (∞) = F(t)dt. (34)
advantage of parallel physical model, in which each layer of ∑n+1 ci
Maxwell element has an independent stress. Using Eq. 29, Note that residual deformation uN (∞) and force data F(t)
we can simulate rheological behaviors and typical results are available from experimental measurements. Matrix Mγ
were shown in Fig. 6b. can be prepared in advance since it only depends on initial
B. Parameter Estimation geometrical coordinates and Poisson’s ratio γ . This equation
allows us to calculate the sum of viscous moduli ∑n+1 ci .
Before simulating any real objects, their physical parame-
ters have to be available in advance. In the above FE model, V. VALIDATION AND N ONLINEAR M ODEL
there are a total of 5 unknown parameters: E1 , E2 , c1 , c2 , A. Experimental Validation
and γ . Our previous work  has proposed an approach to
Experimental data shown in Fig. 5 were taken as ex-
identify these parameters based on iterative FE simulation
amples to quantitatively demonstrate the contradiction in
and nonlinear optimization. The idea is to iterate FE simula-
reproduction of rheological force and residual deformation.
tion with updated physical parameters until the difference
We employed the two-layered Maxwell model and Eq. 29 to
between simulation and experiment becomes minimal. In
simulate the behaviors of these two objects. Using the pa-
, we found that we could estimate the Poisson’s ratio
rameter estimation methods proposed in the last section, we
γ separately because only γ affected the held shape. We
obtained two sets of physical parameters for each material,
could then identify other parameters by minimizing force
as given in Table I. The ﬁrst set of parameters was estimated
difference. This method is simple and works well but it
by minimizing the force difference without considering the
is time-consuming. For force optimization, it took hours or
ﬁnal residual deformation. In other words, we used four
days (depends on initial parameter setting) to obtain results.
unknown parameters (E1 , E2 , c1 , and c2 ) to optimize Eq.
However, according to the analysis in the previous section,
33. During the estimation of the second set of parameters,
we can take the advantage of analytical force expression to
however, we took the residual deformation into consideration
develop a more efﬁcient method.
by using Eq. 34 to calculate c1 + c2 , as given in the second
Extending Eqs. 28 and 21 to 2D formulations, we have
column of Table I. The other three parameters were then
− ci t
F(t) = ∑ ci 1 − e i Mγ vPush, (0 ≤ t ≤ t p ),
N (31) estimated by minimizing force difference. Table I suggested
i=1 that large values of c1 + c2 were needed to achieve good
E STIMATED PARAMETERS FOR BOTH MATERIALS
c1 + c2 E1 E2 c1 c2 e(θ )
(Pa·s) (Pa) (Pa) (Pa·s) (Pa·s) (N2 )
1.3988 × 107 3.1753 × 104 7.2168 × 104 1.3291 × 107 6.9737 × 105 4.0794
9.6961 × 106 3.7731 × 104 8.0952 × 104 9.2023 × 106 4.9380 × 105 24.515
1.3328 × 107 1.0553 × 104 3.7340 × 104 1.3212 × 107 1.1602 × 105 0.8388
1.6849 × 106 2.8909 × 103 1.1612 × 104 3.3585 × 105 1.3491 × 106 186.48
reproduction of rheological forces for both materials, where Experiment
clay has c1 + c2 = 1.3988 × 107 Pa·s and sweets material
has c1 + c2 = 1.3328 × 107 Pa·s. However, in order to well
reproduce the ﬁnal recovered shape, small values of c1 + c2
are necessary, e.g., c1 + c2 = 9.6961 × 106Pa·s for clay and (a-1) Rheological force (a-2) Held shape (a-3) Final shape
c1 + c2 = 1.6849 × 106 Pa·s for sweets material. For both (a) Using the parameters estimated by force optimization alone
materials, the values of c1 + c2 for force reproduction are Experiment Experiment Simulation
larger than the values for deformation reproduction. We Simulation
failed to satisfy the requirement of value c1 + c2 for both
force and deformation reproduction. Two values of c1 + c2
for clay are quite close but very different for sweets material.
It indicated that it is more difﬁcult for sweets material to (b-1) Rheological force (b-2) Held shape (b-3) Final shape
(b) Using the parameters estimated by force optimization with fixed c1+c2
reproduce force and deformation behaviors simultaneously.
Using the estimated parameters listed in Table I, we Fig. 8. Validation results of Japanese sweets material.
could simulate the rheological behaviors of both materials.
Comparison between simulation results and experimental
measurements were shown in Figs. 7 and 8, where Figs.
residual strain. Figures 7 and 8 also showed that the held
7a and 8a were simulated with the ﬁrst set of parameters
shapes simulated with different sets of parameters have
(Table I) for both materials and Figs. 7b and 8b were with
the identical results and quite close to the experimental
the second set of parameters. We found that clay has more
measurements. This is because the held shape was only
coincident results than sweets material using two sets of
affected by the Poisson’s ratio γ  and our FE model
parameters. This again attributes to the value of c1 + c2 with
can reproduce both force and held shape simultaneously. In
close values of c1 + c2 resulting in more coincident results,
addition, it is necessary to mention that the Poisson’s ratios
and vice versa. Figure 8 suggested that it is impossible for
γ listed in Table I were estimated in advance by minimizing
sweets material to reproduce both force and deformation
the difference of the held shape between simulation and
behaviors simultaneously by using the two-layered Maxwell
model. This impossibility will not be changed by adding
more elements on the physical model or by changing the con-
ﬁguration among basic elements. This impossibility comes B. Introduction of Dual-Moduli Viscous Element
from the linearity of the basic elements, especially the linear
viscous element which dominates both stress amplitude and According to the above discussions, we found that the
parameter c1 + c2 must take different values through the
simulation to capture both force and residual deformation.
Experiment Simulation Experiment Simulation
We therefore formulated a viscous element as
σ (t) = (κα + c)˙ (t), (35)
(a-1) Rheological force (a-2) Held shape (a-3) Final shape
where scalars α and c were parameters to be determined.
(a) Using the parameters estimated by force optimization alone Switch function κ takes the following value:
Experiment Experiment Simulation Experiment Simulation
1 t ≤ t p + th ;
−1 t > t p + th .
This dual-moduli viscous element has an ability to switch
(b-1) Rheological force (b-2) Held shape (b-3) Final shape
(b) Using the parameters estimated by force optimization with fixed c1+c2 the parameters from one set to another during simulation.
Introducing Eq. 35 into our FE model, we could modify our
Fig. 7. Validation results of commercial available clay. FE model by replacing the last two differential equations of
simultaneously for sweets materials and we believe that the
E STIMATED PARAMETERS FOR SWEETS WITH NONLINEAR MODEL
linear viscous elements caused the failure. Therefore, a dual-
Parameters Estimated results Parameters Estimated results moduli viscous element was then introduced to improve our
E1 (Pa) 1.0553 × 104 E2 (Pa) 3.7340 × 104 FE model and simultaneous reproduction of rheological force
c1 (Pa·s) 6.6093 × 106 c2 (Pa·s) 7.8618 × 104 and deformation were ﬁnally achieved.
α1 (Pa·s) 6.6027 × 106 α2 (Pa·s) 3.7403 × 104 In the future, nonlinear modeling of rheological objects
will be investigated except switching parameters. Biological
Experiment Experiment Simulation organs and tissues, such as porcine livers and brain tissues,
Simulation will be employed to validate our modeling and parameter
This research is supported in part by Grant in Aid for
Scientiﬁc Research (No. 20246049) and R-GIRO program
(a) Rheological force (b) Final shape
of Ritsumeikan University.
Fig. 9. Validation results of sweets material with the nonlinear model. R EFERENCES
 D. Terzopoulos, J. Platt, A. Barr, and K. Fleischer, “Elastically
deformable models”, Proc. 14th Annual Conference on Computer
Eq. 29 with the following two equations: Graphics and Interactive Techniques (SIGGRAPH ’87), pp. 205–214,
F1 = −
˙ F1 + (λ1 Jλ + μ1 Jμ )vN ,
ela ela  D. Terzopoulos and K. Fleischer, “Modeling inelatic deformation:
κα1 + c1 viscoelasticity, plasticity, fracture”, SIGGRAPH ’88, pp. 269–278,
(37) Atlanta, 1988.
F2 = −
˙ F2 + (λ2 Jλ + μ2 Jμ )vN .
ela ela  M. Hrapko, J.A.W. van Dommelen, G.W.M. Peters, and J.S.H.M.
κα2 + c2 Wismans, “The mechanical behaviour of brain tissue: Large strain
response and constitutive modelling”, Biorheology, vol. 43, pp. 623–
Using the parameter estimation methods proposed in the 636, 2006.
last section, we could determine all the parameters for  A. Nava, E. Mazza, M. Furrer, P. Villiger, and W.H. Reinhart, “In vivo
Japanese sweets material based on the new model, as given mechanical characterization of human liver”, Med. Image Anal., vol.
12, pp. 203–216, 2008.
in Table II. Simulation results compared with experimental  P.Y. Chua, T. Ilschner, and D.G. Caldwell, “Robotic manipulation of
measurements were shown in Fig. 9. We could now achieve food products—a review”, Ind. Robot, vol. 30, pp. 345–354, 2003.
good reproduction of both rheological force and residual  H. Yoshida, Y. Murata, and H. Noborio, “A Smart Rheologic MSD
Model Pushed/Calibrated/Evaluated by Experimental Impulses,” Proc.
deformation simultaneously. IEEE/RSJ International Conference on Intelligent Robots and Systems
(IROS ’05), pp. 269–276, Edmonton, 2005.
VI. C ONCLUSIONS AND F UTURE W ORK  T. Ikawa and H. Noborio, “On the Precision and Efﬁciency of
Hierarchical Rheology MSD Model,” IROS ’07, pp. 376–383, San
In this paper, we summarized physical models for sim- Diego, 2007.
ulating rheological objects. We formulated the generalized e
 B. A. Lioyd, G. Sz´ kely, and M. Harders, “Identiﬁcation of spring pa-
rameters for deformable object simulation,” IEEE Trans. Vis. Comput.
constitutive laws for both serial and parallel models. We Graph., vol.13, no.5, pp. 1081–1094, Sept./Oct., 2007.
found that serial and parallel models have the identical  N. Sakamoto, M. Higashimori, T. Tsuji, and M. Kaneko, “An Optimum
constitutive laws and they could be replaced by each other Design of Robotic Hand for Handling a Visco-elastic Object Based on
Maxwell Model,” Proc. IEEE International Conference on Robotics
with small changes of some coefﬁcients. We also found that and Automation (ICRA ’07), pp. 1219–1225, Roma, 2007.
the serial models are appropriate for calculating strain while  C.-H.D. Tsai, I. Kao, N. Sakamoto, M. Higashimori, and M. Kaneko,
the parallel models are convenient for calculating stress. “Applying Viscoelastic Contact Modeling to Grasping Task: An Ex-
perimental Case Study,” IROS ’08, pp. 1790–1795, Nice, 2008.
Analytical expressions of rheological forces and residual  Y.-H. Chai, G. R. Luecke, and J. C. Edwards, “Virtual Clay Modeling
deformation for generalized parallel models were derived. Using the ISU Exoskeleton,” Proc. IEEE Virtual Reality Annual
Theoretical discussions showed the difﬁculty for linear mod- International Symposium (VRAIS ’98), pp. 76–80, Atlanta, 1998.
 E. Samur, M. Sedef, C. Basdogan, L. Avtan, and O. Duzgun, “A
els to reproduce rheological forces and residual deformation robotic indenter for minimally invasive measurement and characteriza-
simultaneously. A simpler 2D FE dynamic model and a more tion of soft tissue response,” Med. Image Anal., vol.11, pp. 361–373,
efﬁcient method for parameter estimation were presented by Aug., 2007.
 B. Ahn and J. Kim, “Measurement and characterization of soft tissue
taking the advantages of the analytical force expressions. behavior with surface deformation and force response under large
Instead of iterative FE simulations, the parameter estima- deformations,” Med. Image Anal., vol.14, pp. 138–148, 2010.
tion method proposed in this paper only involved direct  J. Muramatsu, T. Ikuta, S. Hirai, and S. Morikawa, “Validation of
FE deformation models using ultrasonic and MR images,” in Proc.
calculations of force expressions and we could obtain an 9th International Conference on Control, Automation, Robotics and
optimal solution within only several seconds. To validate our Vision (ICARCV’06), pp. 1–6, Singapore, 2006.
modeling and parameter estimation methods, experiments  Z. Wang, K. Namima, and S. Hirai, “Physical parameter identiﬁcation
of rheological object based on measurement of deformation and force,”
with commercial clay and Japanese sweets materials were ICRA ’09, pp. 1238–1243, Kobe, 2009.
performed. We found that the value of c1 + c2 strongly  Z. Wang and S. Hirai, “Modeling and parameter identiﬁcation of
affected both force amplitudes and the residual deformation. rheological object based on FE method and nonlinear optimization,”
IROS ’09, pp. 1968–1973, St. Louis, 2009.
We failed to reproduce both force and recovered shape