Modeling and parameter estimation of rheological objects for

Document Sample
Modeling and parameter estimation of rheological objects for Powered By Docstoc
					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 [5]. To date,
both elastic and plastic properties. Due to the presence of
residual (permanent) deformation, it is difficult 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 difficult
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. [2], who
were then derived for generalized parallel models. Theoretical     have proposed several physical models to describe inelastic
discussions suggested the difficulty to reproduce both force        objects and a Burgers model was introduced to describe
and deformation simultaneously using linear physical models.
A 2D FE (finite element) model was then formulated and              rheological behaviors. However, it is only a conceptual
an efficient 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 [6], [7]. The particle-based model has an advantage
   Index Terms— Modeling; Parameter estimation; Rheological        of less computation costs [8], 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 [9] and
                     I. I NTRODUCTION                              Fung’s viscoelastic model [10] 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 [1], [2] and has played an          grasped by a robot hand. Good approximations of force
important role in many fields, 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 finger [11].
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 [6], [7] or reproduction
object, in which the deformation is partly reversible. Previous    of force alone [9], [10]. 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 [12],
brain tissue did not recover completely after a loading-           [13]. However to our knowledge, in current surgical related
unloading cycle [3]. In vivo experimental results showed that      applications, organs and tissues were mostly supposed to be
residual deformation may also present in human liver [4]. 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 [14]. A three-
                                                                                                       E                                                  E1          c1
element [15] and a four-element [16] models were employed
respectively to describe rheological behaviors. A method                                               c1
for identifying physical parameters was proposed based on                                        E1                                                       E1         c1
                                                                                                             E2             c2
iterative optimization [16]. We found that the four-element                                      c1                                                       E2         c2

model was more appropriate than the three-element model                                                                                                    E1        c1
                                                                                            E1               E2
for reproducing rheological forces. In addition, experimental                                                                                             E2          c2
                                                                                            c1               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 configuration 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 configuration 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 efficient 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˙ ,
   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 configurations, 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 coefficients 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
                                                                                                                                                   ∂ iε
                                                                                                                                                                        ∂ jσ
                                                                                                                                                                ∑ Bs2 ∂ t j ,
                                                                                                                                                                   j                    (13)
                                                                                                                                i=1                             j=0
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 = ∑+
                                                                                                                                                   , · · · , Bs2 =
                                   i=1                                                                           En+1         i=1 ci   cn+1 En+1                   cn+1
where ri = Ei /ci . Let us define 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 coefficients 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                                       +
                                                                                   s cn+1
                                                                                             σ (s)
                                                                                                        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 cn+1
                                                              σ (s).                                                           σi +
                                                                                                                                               σi = Ei ε , 1 ≤ i ≤ n,
                               i=1 j=1             j=1
                                         j=i                                                                                                σn+1 = cn+1 ε ,
                                                                                                                                                          ˙                             (15)
                                                                                                  (7)                                                     n+1
We then find 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
                                                                                                                                                                        ∂ jε
where                                                                                                                                ∑ Ai ∂ t i          =     ∑ B p1 ∂ t j ,
                                                                                                                                                                   j                    (16)
                 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 ,
                                                                                                                       p1            p1

Substituting Eqs. 5 and 8 into Eq. 7, we have the following                                                                                                     i=1
                                                                                                                                                          n                             (17)
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:
where                                                                                                                                 n
                                                                                                                                                  ∂ iσ          n
                                                                                                                                                                        ∂ jε
Bs1 =
          1     An
                                Bi,1 An−1
      ∑ ci + cn+1 , Bs1 = ∑ ci + cn+1 , · · · , Bs1 = cn+1 .
                                                       A0                                                                            ∑ Ai ∂ t i           =    ∑ B p2 ∂ t j ,
                                                                                                                                                                   j                    (18)
                                                                                                                                     i=0                       j=1
      i=1                   i=1
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

                                                                                                                                      M                                        M
                                              ∂ iε           n
                                                                     ∂   jσ                                                    En            cn                          Ei        ci

                           ∑ Ai−1 ∂ t i              =   ∑ Bs1 ∂ t j .
                                                            j                                    (12)                                cn+1                                En    M   cn
                           i=1                           j=0
                                                                                                                              (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 coefficients 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
configuration, we are always able to find 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)
models.                                                                           i=1

              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 flat 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 final 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 find an         Integrating the above equation from time t p + th to time
appropriate model to reproduce both rheological force and          infinity, 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,
                                                                                                                                                             ˙                     (23)
                                                                                           σ n +1 = c n +1 &
and thus
  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 )

                                                                                      tp                       tp+th
                                                        (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 infinity,
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 [16] has shown good
∑n+1 σi (t), we finally 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 final 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 firstly calculated the sum of ci based on
and the same formulation of final 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
                     n                E
                                    − 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 find 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 coefficients Ei /ci of exponential functions. The                  force and deformation simultaneously for most rheological
coefficients 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 coefficients of E1 /c1 = 0.005              (e.g., triangles in 2D and tetrahedrons in 3D). The dynamic
and E2 /c2 = 0.5. After determining coefficients 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 + γ )
[14], [15], [16], 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 defined 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)
                                                        (29)                                            i=1
              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 first                 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 ,
                                                       ela   ela
                                                                            4) If the terminate conditions are not satisfied, adjust
 ela , and μ ela could be calculated as follows:
λ2          2
                                                                                 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 definition and description of other variables and            seconds by using MATLAB optimization toolbox.
matrixes can be found in our previous work [14], [15], [16].                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
                                                                                                               i=1             0

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 [16] 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-
[16], 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 first 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
                                                                         final 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 efficient 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
              n                E
                             − 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
                                                                                                      TABLE I
                                                                       E STIMATED PARAMETERS FOR BOTH MATERIALS

                                      c1 + c2                                      E1                          E2                     c1                        c2                 e(θ )
             Material                                              γ
                                       (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
                Clay                                         0.2902
                                   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
              Sweets                                         0.3353
                                   1.6849 × 106                               2.8909 × 103                1.1612 × 104           3.3585 × 105              1.3491 × 106           186.48

                                                                                                                                                                                    Experiment    Simulation
reproduction of rheological forces for both materials, where                                                                      Experiment
                                                                                                                                                       Experiment    Simulation

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 final 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
                                                                                                                                                       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 difficult 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 first 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 γ [16] 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-
figuration 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 ;
                                                                                                                                        κ=                                                               (36)
                                                                                                                                                    −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
                                     TABLE II
                                                                                       simultaneously for sweets materials and we believe that the
                                                                                       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 finally 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
                                                                                       estimation methods.
                                                                                         This research is supported in part by Grant in Aid for
                                                                                       Scientific 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
                                                                                        [1] 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,
                                                                                            Anaheim, 1987.
           F1 = −
           ˙               F1 + (λ1 Jλ + μ1 Jμ )vN ,
                                  ela     ela                                           [2] 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                                           [3] 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                                 [4] 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                              [5] 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                                [6] 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                           [7] T. Ikawa and H. Noborio, “On the Precision and Efficiency 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
                                                                                        [8] B. A. Lioyd, G. Sz´ kely, and M. Harders, “Identification 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                                [9] 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 coefficients. We also found that                                  and Automation (ICRA ’07), pp. 1219–1225, Roma, 2007.
the serial models are appropriate for calculating strain while                         [10] 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                              [11] 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 difficulty for linear mod-                                International Symposium (VRAIS ’98), pp. 76–80, Atlanta, 1998.
                                                                                       [12] 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,
efficient method for parameter estimation were presented by                                  Aug., 2007.
                                                                                       [13] 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                                [14] 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                                 [15] Z. Wang, K. Namima, and S. Hirai, “Physical parameter identification
                                                                                            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                                 [16] Z. Wang and S. Hirai, “Modeling and parameter identification 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

Shared By: