Document Sample

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 gr046074@ed.ritsumei.ac.jp, hirai@se.ritsumei.ac.jp 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 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. [2], 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 [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 ﬁ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 [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 c2 respectively to describe rheological behaviors. A method c1 c2 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 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ε, ˙ ˙ c 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 (2) A serial rheological model consists of numbers of Kelvin ε= ∑ εi . i=1 elements and a viscous or a Maxwell element connected in E1 E2 En E c cn+1 L (a) (b) c1 c2 cn E (a) Type 1 E c E1 E2 En c En+1 cn+1 L (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: n+1 (3) ε (s) = ∑ εi (s). n+1 ∑ Ai−1 ∂ t i ∂ iε = n+1 ∂ jσ ∑ Bs2 ∂ t j , j (13) i=1 i=1 j=0 Eliminating εi (s) from the above equations, we then have where n 1 1 1 1 ε (s) = ∑ s + ri ci + s cn+1 σ (s), (4) Bs2 = n+1 1 n 1 , Bs2 = ∑+ An + n An−1 , · · · , Bs2 = 0 A0 . 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, (6) 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 n (s + r j ) n n s+r j 1 Ei = ∑∏ ci +∏ s cn+1 σ (s). σi + ˙ ci σi = Ei ε , 1 ≤ i ≤ n, ˙ i=1 j=1 j=1 j=i σn+1 = cn+1 ε , ˙ (15) (7) n+1 We then ﬁnd the following equation: σ= ∑ σi . i=1 n ∏ (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 j=i (Fig. 4a) as follows: = sn−1 + Bi,1 sn−2 + · · · + Bi,n−2s + Bi,n−1, n ∂ 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 = n n 1 An n−1 n Bi,1 An−1 ∑ ci + cn+1 , Bs1 = ∑ ci + cn+1 , · · · , Bs1 = cn+1 . 0 A0 ∑ Ai ∂ t i = ∑ B p2 ∂ t j , j (18) i=0 j=1 i=1 i=1 (11) 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 n+1 ∂ 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. where 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 t 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 t 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) 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 ﬂ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 target. ∑ ci t p +th εivis (t)dt + cn+1 ˙ t p +th ε (t)dt = 0, ˙ (23) i=1 σ n +1 = c n +1 & ε and thus n σ (t p ) ∑ ci εivis (∞) − εivis (t p + th ) + cn+1 ε (∞) − ε (t p + th ) = 0. Pushing Holding σ (t p + t h ) Pushing Holding σ (t p + t h ) σ (t p ) i=1 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 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. n ci εivis (t p + th ) cn+1 ε (t p + th ) ε (∞) = ∑ + . (25) i=1 ∑n+1 c j j=1 ∑n+1 c j j=1 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 ﬁ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 i=1 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 i=1 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 ), N 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 2D 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 N 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 ˙ T 2 e(θ ) = ∑ fexp − fcal (t, θ ) 2 , i i (33) (29) i=1 E1 F1 = − F1 + (λ1 Jλ + μ1 Jμ )vN , ˙ ela ela where N is the number of sampling points, vector fexp i c1 E2 and fcal (t, θ ) are experimental forces and calculated forces, i F2 = − F2 + (λ2 Jλ + μ2 Jμ )vN , ˙ ela ela respectively, at the i-th sampling time, and vector θ is the c2 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 , ela ela 4) If the terminate conditions are not satisﬁed, 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 deﬁnition 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 1 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 . i=1 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 ﬁ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 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 Simulation 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 ﬁ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 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 γ [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 experiment. 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 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 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 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 estimation methods. ACKNOWLEDGMENTS 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 [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. E1 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. E2 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 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 [8] 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 [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 coefﬁcients. 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 difﬁculty 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, efﬁcient 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 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 [16] 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

DOCUMENT INFO

Shared By:

Categories:

Tags:

Stats:

views: | 4 |

posted: | 3/14/2012 |

language: | English |

pages: | 8 |

OTHER DOCS BY zhouwenjuan

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.