Finite Element Modeling and Physical Property Estimation of Rheological Food Objects

The purpose of this study is to accurately simulate the rheological behaviors of food objects undergoing a loading-unloading operation using finite element (FE) model. Due to the presence of residual deformation, it is difficult to model rheological behaviors. Especially, it is hard to accurately reproduce both rheological force and residual deformation simultaneously. In this study, objects made of food materials were tested. Force and deformation measurements were recorded for parameter estimation. Constitutive models were investigated for describing rheological behaviors. A parallel five-element model including two dual-moduli viscous elements was proposed to accurately predict both rheological force and residual deformation simultaneously. 2D/3D FE model was formulated for simulating rheological behaviors. To estimate the parameters, an effective four-step method was established based on nonlinear optimization which aimed at minimizing the differences of forces and deformation between simulation and experiments. The proposed FE model and parameter estimation method were validated in both 2D and 3D cases and good agreements were achieved in both rheological forces and deformation between numerically simulated and experimentally measured data.


Introduction
The modeling and simulation of deformable objects, such as biological organs and tissues, cloth, clay, and various kinds of food products, has been studied more than twenty years.Many important applications have been involved, including computer graphics [Terzopoulos & Fleischer (1988)], surgical simulation [Bro-Nielsen (1998), Cotin et al. (1999)], robot manipulation [Inoue & Hirai (2006)] and food engineering [Liu & Scanlon (2003), Martins (2006)].Three important issues should be considered during the modeling and simulation of deformable objects: 1) the constitutive model, which normally is used to govern the physical behaviors of the material, 2) the modeling method, which is used to construct the geometry of the object and formulate the dynamic equations, and 3) parameter estimation, which is dedicated to find appropriate physical parameters involving in the model.Among various constitutive models, the ones consisting of a number of elastic and viscous elements connected in a certain configuration were widely used in literatures.Such models have explicit physical meaning, simple formulation, and are easy to be constructed, but they are linear models.Conventional modeling methods are the mass-spring-damper method (MSD) [Waters (1987)], the finite difference method (FDM) [Terzopoulos et al. (1987)], the boundary element method (BEM) [James & Pai (1999)], and the finite element method (FEM) [Beo-Nielsen & Cotin (1996)], in increasing order of computation cost and simulation accuracy.Most recently, meshfree particle methods, such as smoothed particle hydrodynamics (SPH), have also been used to model solid deformable objects [Zhu et al. (2010)].A combination of a constitutive model and a modeling method can be employed to construct a two-dimensional (2D) or three-dimensional (3D) dynamic model.In addition, to simulate real-world deformable objects, important physical parameters need to be estimated in advance.So far, the optimization-based methods have been the most popular ones for the purpose of parameter estimation.
In our definition, we roughly divided deformable objects into three categories: elastic, plastic, and rheological objects, depending on their behaviors after a loading-unloading operation.Rheological objects have both elastic and plastic properties and always yield residual deformation after unloading.Residual deformation is important for some applications, such as manufacture of pottery and food products, where a desired final shape is required during forming process.Therefore, it is necessary to investigate the modeling and parameter estimation of such rheological objects.
Early work on the modeling of rheological objects dates back to Terzopoulos et al. (1987) and Terzopoulos & Fleischer (1988), who have proposed a four-element model to describe rheological behaviors.Detailed studies on the modeling and parameter estimation of rheological objects have been investigated by Noborio et al., who employed a three-element model to describe rheological behaviors and an MSD modeling method to construct food dough.The authors explored the lattice [Noborio et al. (2003), Nogami et al. (2004a)], the truss [Nogami et al. (2004b)], and the hierarchical [Ikawa & Noborio (2007)] structures, with decreased MSD elements connected between nodes to reduce the computation cost.Two optimization methods, modified randomized algorithm [Noborio et al. (2003)] and genetic algorithm [Yoshida et al. (2005)], were used to estimate the physical parameters.The authors successfully captured the deformation behaviors, but failed to reproduce the force responses at the same time [Yoshida et al. (2005)].The MSD modeling method has advantages of simple formulation and relatively low computation cost, but the formulation is not based on continuum mechanics and the geometrical topology significantly affects the simulation accuracy.
Two layered Maxwell model [Sakamoto et al. (2007)] and Fung's viscoelastic model [Tsai et al. (2008)] have been used respectively to simulate the force response of 'Norimaki-sushi' when it was grasped by a robot hand.Good agreements of rheological forces were achieved.Unfortunately, both models are 1D cases and the residual deformation was not considered during the modeling and parameter estimation.FE method has been utilized to model and simulate the indentation of bread crumbs [Liu & Scanlon (2003)].It has also been used to evaluate food quality and safety losses during processing, storage, and distribution [Martins (2006)].In addition, FE simulation has been employed to investigate the dependence of temperature and water content on processing time during meat cooking [Purlis and Salvadori (2005)].It has been reported that the most critical barrier against the application of robotics and automation in food industry is a lack of understanding of the food product properties as an "engineering" material for handling operations [Chua et al. (2003)].Therefore, it is necessary to investigate the methods for modeling and parameter estimation of rheological objects.
In this study, three kinds of rheological food materials: commercial available clay (food-like), Japanese sweets, and bacon, were tested with a loading-unloading operation.Rheological force and deformation, including the residual deformation, were recorded.To find an appropriate model for describing these behaviors, the constitutive laws of generalized models were formulated and theoretical analysis was done.A parallel five-element model with two dual-moduli viscous elements was proposed to accurately predict both rheological force and residual deformation simultaneously.The 2D/3D FE model was then formulated by imposing a five-element model onto each triangle (2D) or tetrahedron (3D) to govern its behaviors.Based on the analytical expressions of rheological forces, an efficient approach was proposed to estimate the physical parameters involving in the FE model.This method aimed at minimizing the difference in both force and deformation between simulation results and experimental measurements.The proposed FE model and parameter estimation method were then evaluated through a series of comparisons of experimental measurements and predicted results from FE simulation.
There are some potential applications of this research in food industry.Using the proposed model and method, we are able to estimate the physical properties of food products and then establish a relationship between the physical parameters and the taste of the food product.As we know, the taste of a food product includes not only the flavor but also the chew feeling, such as the hardness and the viscosity of the product.We can use the estimated parameters to evaluate the taste of food products.The second potential application is to simulate and predict the deformed shape of food product.For example, we may need to perform a pick-and-place operation in a sushi manufacturing line.In order to grasp the sushi stably, a certain grasping deformation is required.On the other hand, large deformation may affect the final appearance of the product.There is a compromise between enough grasping deformation and appropriate final appearance.We can determine an appropriate grasping deformation by performing a series of experiments which may be time-consuming and troublesome.Instead, we are able to use the FE model of the sushi product to simulate such grasping operation to find an optimal grasping deformation.In addition, the developed FE model can be also utilized to establish a virtual scenario of food making process with haptic feedback since good reproduction of force response can be achieved using the proposed parameter estimation method.Such a virtual scenario can serve as a digital demonstration to show the skillful making process of some traditional food products.

Materials and Methods
The goal of this study is to develop an FE model and establish a parameter estimation method to accurately predict rheological force and deformation behaviors simultaneously.To this end, indentation experiments were performed using typical rheological materials.Experimental measurements of force and deformation were used to estimate the physical parameters.The FE model and estimated parameters were then used to predict rheological behaviors of same materials.

Experimental materials
Three typical rheological materials, i.e., commercial available clay, Japanese sweets, and bacon, as shown in Figure 1, were used in our indentation experiments.Clay is made of flour, salt, and water and supposed to be played by children over 3 years old.Three kinds of sweets materials were provided by OIMATU (a sweets company in Kyoto) and used to make traditional Japanese sweets products.Each sweets material is made of flour, water, and one kind of bean powder mixed at a specific ratio.Bacon was bought from a supermarket and normally was cut in pieces and eaten with bread for breakfast.Several object samples made from each material were prepared and indented by a linear stage.

Indentation experiments
For each material, two kinds of indentation experiments were performed.The first kind is for parameter estimation.The sample object was indented from one side using a motorized linear stage (KX1250C-L, SURUGA SEIKI Co.), as shown in Figure 2a.A push-keep-release displacement function, as shown in Figure 2b, was used to deform the object.The input deformation is assumed to be plane stress.The force response on the other side of the object was recorded using a tactile sensor (I-SCAN100L, NITTA Co.).Three static images, which denote the initial, keep, and final shapes respectively, were recorded using a digital camera (EOS Kiss X2, Canon Inc.).The deformation during keep phase and residual deformation were then calculated based on these images.As an example, Figure 3 shows the measured rheological behaviors of object made by clay material.These measurements of force and deformation were used to estimate the physical parameters.
The purpose of the second kind of experiments is to evaluate the performance of estimated parameters for predicting the behaviors of each material with different indentation operations, non-homogeneous, and irregular-shaped objects respectively.For clay material, sample objects were indented from the center part of one side (Figure 4a) instead of the entire side.Such an operation may often be used in the grasping of an object.In terms of Japanese sweets materials, non-homogeneous three-layered objects were deformed from both entire and center part of the side (Figure 4b).Irregular shaped objects (Figure 4c) were used to evaluate the parameters of bacon material.A push-keep-release procedure was also used in these experiments.Detailed information of all experimental trials is given in Table 1.In addition, it is worth to be mentioned that we performed the deformation measurements in 2D scenario and ignored the buckling or bulging behaviors which slightly happened during the deformation.

Constitutive models
The constitutive models, such as Maxwell and Voigt models, were widely used to describe the behaviors of deformable objects.In such models, we commonly utilize an elastic element (denoted by Young's modulus E) and a viscous element (denoted by viscous modulus c) to describe the elastic and viscous behaviors, respectively.By connecting several such elements into different configurations, we can construct various models for simulating rheological behaviors, such as the three-element and four-element models used in [Wang et al. (2009), Wang & Hirai (2009)].In order to investigate the capability of the constitutive models and to establish a criterion for selecting appropriate models, we have summarized such models into two categories: serial and parallel models, as shown in Figure 5.The constitutive laws of generalized serial and parallel models were derived in [Wang & Hirai (2010a)] and summarized in Table 2.We found that the constitutive laws of serial and parallel models have the same forms and therefore yield the same behaviors.In other words, for a certain serial model A, we are always able to find a corresponding parallel model B, which has the same constitutive law and yields the same behaviors with model A, and vice versa.In addition, we can always obtain analytical expressions of strain (displacement) over time by solving the constitutive laws of serial models due to the serial connections between elements.On the other hand, parallel models always result in straightforward calculation of analytical stress (force) expressions due to the parallel connection.This suggests us the way to select appropriate one from serial and parallel models.If we are interested in the deformation, it is better to choose a serial model.On the contrary, we should go with parallel model if the force is of concern.
In our case, we are more interested in rheological force since deformation measurements only include static images.We have therefore chosen parallel models.In [Wang & Hirai (2009)], we also found that at least two exponential terms in force expressions are necessary to accurately reproduce force relaxation behavior during keep phase (Figure 3a).Consequently, we need at least two Maxwell elements in the constitutive model, e.g., a parallel five-element model (the bottom one in Figure 5b with n = 2).By using such a model, we are able to accurately reproduce the rheological force behaviors.However, we failed to predict the residual deformation at the same time.We found a contradiction between the reproductions of rheological force and residual deformation simultaneously because the linear viscous elements (denoted by parameter c i in Figure 5b) affected both force and residual deformation [Wang & Hirai (2010a)].One single set of parameter c i can guarantee a good reproduction of either force or residual deformation.However, it cannot cover both at the same time.We have therefore proposed a dual-moduli viscous element, as shown in Figure 6a, to solve this problem.The governing equation of this element can be formulated as where dual  and dual  are stress and strain vectors generated on the element, scalars  and c are parameters to be determined, switch variable  takes the following values The dual-moduli viscous element works as a switch to change the viscous coefficient from one value to another during simulation whenever a criterion is satisfied.The physical meaning of this element can be explained as the property change of a material during loading and unloading operations.For example, elastic materials experience a hysteresis phenomenon during which the material properties are slightly changed.Some metal materials also demonstrate strain hardening behavior when they are strained beyond the yield point [Shames & Cozzareli (1992)].The properties of the materials also changed after strain hardening.In rheological materials, both hysteresis and strain hardening may also happen and may be in a stronger way.In other words, the physical parameters of rheological object may be continuously changing during loading and reach another set of values once loading operation finishes.Unfortunately, continuous change of parameters brings troubles in parameter estimation.Therefore, we suppose that the parameters take one set of values during loading operation and switch to another set once the operation finishes.
The criterion used in Eq. 2 has different options depending on applications.If the operation time is available, the simulation time can be a criterion.In some applications, the simulation time may be not available in advance.Fortunately in most cases, an interaction often happens between the object and external instruments.This interaction can also serve as a criterion for switching parameters, as presented in [Wang & Hirai (2010b)].By introducing the dual-moduli viscous elements, we finally end up with an effective model (Figure 6b) for predicting both rheological force and residual deformation simultaneously.

FE dynamic model
The FEM is the most successful method for numerical computation of object deformation.In FE modeling, an object is described by a set of elements (e.g., triangles in 2D and tetrahedra in 3D cases).Dynamic behaviors of the object are then determined by analyzing the behaviors of individual element.In this paper, we assume that the behaviors of individual element are governed by the model given in Figure 6b.Note that the constitutive law of a Maxwell model with a dual-moduli viscous element located at the first row of the five-element model (Figure 6b) is formulated as By performing a series of replacements as done in [Wang & Hirai (2011a)], Eq. 3 can be transformed to 2D/3D force-displacement relationship as where 1 F is the force vector generated on the first row of the five-element model, vector N v consists of velocity components of all nodes in the FE mesh,  J and  J are referred to as connection matrices that are constant matrices and only depend on geometric quantities, say, initial coordinates of nodes, scalars 1 ela  and 1 ela  denote Lamé's constants and can be calculated as: where parameter E 1 is Young's modulus and  is Poisson's ratio which denotes the compressibility of the material.For an incompressible material, parameter  usually takes a value smaller than but very close to 0.5, e.g., 0.499 . For compressible materials, parameter γ should take a value much smaller than 0.5.Similarly, the force-displacement relationships at the second and third row of five-element model (Figure 6b) can be formulated as where parameters 3 vis  and 3 vis  describe the model's viscosity and are defined as Due to the parallel connections, the total force rheo F generated on the model can be calculated by summing up the contributions of each row as Supposing that one side of the object is fixed on the ground and the other side is pushed downward with a displacement function of ( ) t d .Two constraints on the nodes of both sides are formulated as follows using the constraint stabilization method (CSM) [Baumgarte (1972)] where matrices A and B denote the nodes to be constrained on each side, respectively, scalar  is a predetermined angular frequency and is set to 2000 for both constraints.
Let M be an inertia matrix and 1 λ and 2 λ be the Lagrange multipliers that denote a set of constraint forces corresponding to both geometric constraints (Eq.9).Applying Lagrange equations of motion, a set of dynamic equations of nodes is formulated as Combining Eqs. 4, 6, 8, 9, 10 and considering N N   v u , we have a set of differential equations for simulating the dynamic behaviors of a rheological object.Note that the above formulations are applicable in both 2D and 3D scenarios.The main difference between the 2D and 3D models is the calculation of connection matrices  J and  J , which depend on the geometrical meshes.

Parameter estimation
In order to accurately predict the behaviors of real-world objects, physical parameters involved in the FE model need to be determined.The FE model presented in the previous section includes the following eight physical parameters: Young's moduli E 1 , E 2 , viscous moduli 1  , 2  , c 1 , c 2 , c 3 and Poisson's ratio  .In [Wang & Hirai (2010a)], we found that parameter c 3 accounts for attenuating the vibration in both force and deformation during object recovery.Simulation showed that a small value of c 3 is sufficient to attenuate the vibration and without significantly affecting the amplitudes of force and deformation.Therefore, we preassigned a value of 10 2 Pa•s to c 3 .It is quite small comparing with c 1 and c 2 which normally have a magnitude around 10 6 Pa•s.After determining parameter c 3 , we have seven physical parameters to be estimated.We therefore propose an estimation method based on inverse FE optimization, which aims at minimizing the difference in force and deformation between simulation or theoretical calculation and experimental measurements.The measurements used here include force data during loading (includes push and keep phases) and three static images denoting the initial, keep, and final shapes of the object.The parameter estimation method proposed in this paper is summarized into the following four steps: .

Estimation of Poisson's ratio 
In [Wang & Hirai (2009)], we found that only Poisson's ratio  affects the keep-shape and all the other parameters do not have contributions to keep-shape.This coincides with the definition of Poisson's ratio which indicates a ratio of lateral strain over normal one and it happens in all constitutive models constructed by linear elements shown in Figure 5. Therefore, we are able to estimate Poisson's ratio  separately by minimizing the difference in keep-shapes.The objective function is formulated as where Vector push N v consists of velocities of all nodes in push phase.We assume that it consists of constant values if the object was deformed with a small and constant velocity.Therefore, it can be easily calculated using the displacements of nodal points divided by time period p t .
Similarly by solving Eq. 6, we can formulate the force expressions at the second and third row as Combining the forces of three rows, the rheological force generated on the five-element model during the push phase ( In Using Eqs. 14 and 15, we are able to calculate the forces in both push and keep phases.The difference between the calculated forces and experimentally measured ones is then minimized to estimate the parameters.The objective function is formulated as where vector consists of the parameters to be determined.Vector exp i

F
is the force measurements from experiments at the i-th sampling time and vector The termination threshold is the tolerance on   E θ or θ less than 6 1 10   .Instead of using the calculated forces, we can also obtain a set of rheological forces by running FE simulations during the optimization process.However, the optimization will become very time-consuming since four unknown parameters are involved and the FE simulation has to be iterated many times.In [Wang & Hirai (2010a)], we found that it is hard to accurately predict both rheological force and the final-shape simultaneously by using a linear FE model.There is a contradiction between the reproduction of force and the approximation of final-shape.We have therefore introduced a dual-moduli viscous element (Figure 6a) to deal with this issue.This element switches the viscous moduli from one value ( and 2 rel c need to be estimated.Since force is out of concern after releasing, we only focus on the deformation and we can determine these two parameters by minimizing the difference in final-shapes between simulation and experiment.The objective function is formulated as

Estimation of parameters
For some applications, we may only concern about the forces but not residual deformation.In such cases, it is not necessary to include the dual-moduli viscous elements in the constitutive model and we should also skip the third and fourth steps during parameter estimation.Note that the parameter estimation method proposed in this section can be applied in both 2D and 3D models.

Parameter estimation results
Generally, the material property of an object will not differ even though the object may be subjected to different operations or has different shapes or sizes.This feature allows us to use regular-shaped objects with simple loading operations to estimate their parameters.As a result, the estimated parameters should be able to simulate arbitrary shaped objects with any operations.In our experiments, we used flat-squared objects pushed from one side with constant velocity to estimate the parameters.Using the method proposed in Section 2.5, we estimated the physical parameters for the experimental trials given by Table 1 with 'Est-' in the trial names.Estimation results are listed in Table 3.
Note that all objects made of clay material are supposed to have the same properties since they were bought in the same pack and at the same time.The tendency of estimated parameters relative to different compressing velocities is shown in Figure 7.We found that there is no clear tendency of parameters  , E 1 , and E 2 relative to compressing velocities.On the other hand, parameters c 1 and c 2 are slightly increased along with the decreasing of compressing velocities.This is reasonable since the viscous moduli are normally used to describe the velocity related behaviors.
Because three kinds of sweets materials are made of different bean powders and with different mixture ratios, the estimated parameters are also quite different among one another.For the bacon material, due to the non-homogeneous and anisotropic properties, different compressing directions affected the physical parameters, which can be seen from the estimated parameters of bacon trials listed in Table 3. 3D mesh information and time cost for estimating parameters of two bacon trials are given in Table 4.We can see that time cost in the second step is very low comparing with the first and third steps by taking the advantage of analytical expressions of rheological force.In the first and third steps, the FE simulation was iterated several times and it costs longer time to reach an optimal solution.Comparing with the first step, two variables ( 1 rel c , 2 rel c ) are involved in the third step and yield higher cost in parameter estimation.The estimated parameters listed in Table 3 will be evaluated in the next subsection by predicting unknown behaviors of each material under different operations and with different shaped object.

Evaluation results of clay material
Evaluation trials of clay materials were performed using white colored objects with compressing from the top-center area and different compressing velocities.Detailed experimental information of these trials is given in Table 1 with 'Eva-' in the beginning of the trial names.The average values of estimated parameters are used to predict the rheological behaviors of these trials.Simulation results comparing with experimental measurements are shown in Figure 8.We found that both keep-shapes and final-shapes are pretty well matched between simulations and experiments but the predictions of rheological forces showed some errors in all trials.

Evaluation results of sweets material
To evaluate the estimated parameters of sweets materials, several non-homogeneous layered objects were tested.Each of them consists of three layers with two different kinds of sweets materials.Two kinds of compressing operations were performed on these objects.One is from entire top side of the objects and the other one is from the center area of the top side.The FE modeling of non-homogeneous rheological object was presented in [Wang & Hirai (2011b)].The idea is to virtually separate a non-homogeneous object into several homogeneous ones with their own physical properties.The behaviors of the non-homogeneous object were then simulated by combining the behaviors of individual homogeneous ones using a set of boundary constraints.The estimated parameters for each sweets material listed in Table 3 are used to simulate the non-homogeneous objects.
Simulation results comparing with experimental measurements are shown in Figures 9 and 10. Figure 9 shows that both deformed shapes and forces are pretty well matched between simulation and experiment.In Figure 10, however, the force profiles experience certain errors comparing with the results in Figure 9.This is because the compressing operation.In Figure 9, the objects were compressed from the entire top sides which are the same with the operation used for parameter estimation.On the other hand, the compressing operation in Figure 10 is different from the one used in parameter estimation.Therefore, we can conclude that similar operations used for parameter estimation will yield better performance in prediction.Even with different operations as shown in Figures 8 and 10, we still can obtain acceptable prediction results, which denote that our FE model and parameter estimation method are able to predict both rheological force and deformation behaviors simultaneously, even for non-homogeneous object with different operations.

Evaluation results of bacon material
Evaluation results presented above are all based on 2D FE model.To validate our model and method in 3D cases, we performed 3D parameter estimation and evaluation using bacon material.Two flat-squared bacon slices were firstly tested with a push-keep-release procedure for estimating the parameters.Experimental data consist of force and static images, which are same with 2D case.However during parameter estimation, 3D FE model was used instead of the 2D one.3D tetrahedra meshes are generated using an open-source mesh generator named 'TetGen ' [Si (2009)], which creates a 3D mesh based on the information of the object boundary.The estimated parameters are listed in Table 3.
Two experimental trials were performed with irregular shaped (trapezoid and diamond) bacon slices for evaluating the estimated parameters.The muscle fiber textures of trial 'Eva-B-V-T' and 'Eva-B-P-D' are vertical and parallel to the compressing direction respectively.To show the effects of the muscle fiber texture, both sets of estimated parameters are used to simulate the two evaluation trials.Simulation results comparing with experimental ones are shown in Figures 11 and 12, respectively.
Figure 11 shows that the simulated deformation behaviors matched the experimental ones quite well with both sets of estimated parameters.Unfortunately, the force experiences certain errors with both sets of parameters.We can see that the simulated force profiles are quite similar for both sets of parameters but the force amplitudes are quite different.The parameter set 'Est-B-V' yields larger force amplitude (Figure 11b-1) comparing with 'Est-B-P'.From Table 3 we also found that the parameter set 'Est-B-V' has larger values than parameter set 'Est-B-P'.This means that we feel harder if we compress the bacon slice in a direction vertically to the muscle fiber texture.From Figure 11, it is hard to say which set of parameters is better for predicting the force profile.
On the other hand, from Figure 12 we can easily tell that the parameter set 'Est-B-P', whose muscle texture direction is the same with the evaluation trial 'Eva-B-P-D', yields much better predictions of force profile.This tells us that the compressing direction affects the parameter estimation results for anisotropy objects and further affects the prediction results as well.Therefore, it would be better if we could perform the experiments for estimation using the same operation, especially the compressing direction, with the one using in prediction.
The simulated deformation behaviors shown in Figures 11 and 12 are 2D surface view projected from 3D simulations for the convenience of comparison.Figure 13 shows the 3D simulation snapshots for both evaluation trials.The 3D views of keep-shapes for both trials are shown in Figures 13a and 13b.Figures 13c and 13d are used to demonstrate the change in object thickness of both trials.To accurately approximate the experimental condition (Figure 2), one surface of the objects (left sides in Figures 13c and 13d) is constrained to make sure no deformation happens during simulation and only the other side of the objects is allowed to deform freely.Evaluation results shown in Figures 11, 12, and 13 suggested that our FE dynamic model and parameter estimation method are applicable in 3D cases as well.

Conclusions
In this paper, we presented 2D/3D FE dynamic model and a parameter estimation method for accurately simulating the behaviors of rheological objects.Due to the presence of residual deformation, modeling of rheological object is more difficult than doing an elastic one.Especially, it is hard to accurately predict both rheological force and residual deformation simultaneously.Therefore in this paper, we proposed a parallel five-element model with two dual-moduli viscous elements to deal with this issue.The dual-moduli viscous element has an ability to switch the viscous coefficient from one value to the other during simulation.This makes the object behave differently under loading and unloading operations, which can be physically explained as property change caused by the deformation.By imposing a parallel 5-element model onto each triangle or tetrahedron to govern its behaviors, 2D/3D FE dynamic model was formulated to simulate homogeneous objects and it can be easily extended to simulate non-homogeneous objects as well.There are 7 unknown physical parameters in this model.To determine these parameters, a four-step method for parameter estimation was proposed.This method aims at minimizing the difference in rheological force and deformation between simulation and experiments.By taking the advantage of parallel configuration of the constitutive model, we can analytically solve the force expressions.Consequently, the efficiency of the estimation method was improved significantly.It can reach an optimal solution in the second step within several seconds.This method can be also applied in other FE models as long as a parallel model is used.
Three typical rheological materials: clay, Japanese sweets, and bacon, were tested to validate the proposed FE model and parameter estimation method.Several square-flat shaped objects made of each material were indented and the measurements of force and deformation were recorded.Parameter estimations were performed with 2D FE model for clay and Japanese sweets materials and with 3D model for bacon material.Estimated parameters show that the indentation displacement and velocity do not affect the estimation results much but the indentation directions for an anisotropic object do affect the estimated parameters.To evaluate the estimated parameters, indentation experiments were performed for each material with different compressing operations, non-homogeneous, and irregular-shaped objects, respectively.Evaluation results show that we can predict the deformation quite well but the predictions of rheological forces experience certain errors for some trials.Using same compressing operations and same indentation directions in both estimation and evaluation always yield better prediction results.Due to the difference in object shapes, properties, and compressing operations, we think the prediction errors in rheological forces shown in Figures 8,9,10,11,and 12 are acceptable and the prediction results are accurate enough for some applications, such as haptic display of object manipulation.Evaluation results finally suggested us that our FE model and parameter estimation method are suitable for simulating rheological objects and can successfully reproduce both rheological force and deformation behaviors simultaneously, even for non-homogeneous, anisotropic, and 3D objects.
This paper is a study on methodology of object modeling but not on particular application.The proposed FE model and parameter estimation method can be easily applied in various applications as long as the measurement data of force and deformation field are available.They are also not limited to model rheological objects but can be used to model elastic, viscoelastic, and plastic objects as well with slight change in the constitutive model.Note that the dual-moduli viscous element is not necessary to be used in modeling elastic or viscoelastic objects since there is no residual deformation happens in such objects.Accordingly, the third and fourth steps in parameter estimation method can be skipped for dealing with elastic or viscoelastic objects.
In this paper, only compression loadings and isotropic properties were considered.In the future, tensile loading, orthotropic, and anisotropic properties will be considered as well to better understand the behaviors of rheological objects.The first term in the trial name denotes that the trial is used to estimate the parameters ('Est-') or evaluate the parameters ('Eva-'), respectively.The second term denotes different materials with '-C-' standing for clay, '-S-' for sweets, and '-B-' for bacon materials, respectively.The third term in clay and sweets trials denotes the object color with '-R-' for red, '-B-' for blue, '-W-' for white, '-Y-' for yellow, and '-P-' for purple, respectively.The third term in bacon trials denotes that the pushing direction is parallel ('-P-') or vertical ('-V-') to the texture of muscle fiber.The remaining terms in different trials are some parameters used to distinguish between one another, for example, the numbers ('-06', '-08', and '-10') in the clay trials for parameter estimation denote rough indentation displacements during experiments, the letters ('-D' and '-T') in last two trials of bacon denote that the objects have irregular shapes of diamond and trapezoid, respectively.
In serial models, 'free elastic' indicates that an elastic element appears without a viscous element connected in parallel.In parallel models, 'free viscous' indicates that a viscous element appears without an elastic element connected in serial.

Sweets 2D
Est-S-W 0.3746 1.3468×10 4 2.4695×10 4 1.4820×10 7 5.3855×10 4 1.4811×10 7 1.8527×10 4 Est-S-Y 0.3353 1.0553×10 4 3.7276×10 4 6.6096×10 6 7.8271×10 4 6.6034×10 6 3.7659×10 4 Est-S-P 0.3267 9.1565×10 3 5.0802×10 4 4.0958×10 6 8.2198×10 4 4.0851×10 6 5.2072×10 4 Bacon 3D Est-B-P 0.4928 5.8250×10 4 1.1640×10 4 2.8537×10 7 1.8828×10 4 2.8513×10 7 1.7652×10 4 Est-B-V 0.4923 6.7607×10 4 1.2988×10 4 3.7292×10 7 2.4118×10 4 3.7243×10 7 2.2872×10 4 Parameters of clay and sweets materials were estimated using 2D FE model and parameters of bacon were estimated using 2D measurements but 3D FE model.The object is pushed from time 0 to t p by the linear stage with a constant velocity and this time period is called push phase.During time t p to t p +t k (keep phase), the linear stage is stopped and keep the deformation unchanged.Accordingly, the deformed shape during this period is called keep-shape.After time t p +t k , the linear stage is moved back to the original position and the deformation generated inside the object is allowed to recover freely and finally reach a permanent shape, which is called final-shape in this paper.During push phase, the rheological force is increasing in an almost linear manner.During the keep phase, the deformation is kept unchanged (keep-shape) but the force response is decreasing in a nonlinear manner, which is referred to as force relaxation.After releasing, the rheological force goes to zero but the deformation starts to recover and finally reach a permanent shape called final-shape.
Estimation of parameter  by optimizing the keep-shapes, .⑵ Estimation of parameters E 1 , E 2 vectors in the keep phase from simulation and experiment, respectively.Scalar m=2N with N being the total number of nodes in the FE mesh.The optimization is terminated when the tolerance on function value   Considering the force response during loading, let variable  in Eq. 2 take 1 and we have Thanks to the parallel configuration of the constitutive model, we are able to derive the analytical expression of rheological force during loading.By solving Eq. 4 in push phase ( p 0 t t   ), we have the force expression at the first row of the five-element model as To accurately reproduce the final-shape, two more parameters 1 rel c

Figure 1 .
Figure 1.Objects made by three typical rheological materials used in experimentsThe black dots on the surface of the objects were manually drawn before experiments for the convenience of making FE meshes and clearly capturing the deformation.

Figure 3 .
Figure 3. Experimental measurements of rheological behaviors

Figure 4 .
Figure 4. Different pushing operations and different shaped objects used in evaluation experiments.In clay case, the square-shaped objects were deformed from the center part of the top sides instead of entire top sides in estimation trials.In the case of Japanese sweets, non-homogeneous layered objects were used instead of homogeneous objects in estimation trials.In bacon case, irregular-shaped objects were used instead of regular-shaped objects in estimation trials.

Figure 6 .
Figure 6.The dual-moduli viscous element (a) and the parallel five-element model (b) with two dual-moduli viscous elements, which is used in this paper to model rheological objects

Figure 8 .
Figure 8. Evaluation results of clay objects with a compressing velocity of (a) 0.5m/s, (b) 0.2m/s, and (c) 0.1m/s, respectivelyBecause the deformation behaviors are complicated (especially in the contact corners), a triangular mesh with 17×17 nodal resolution was used in simulations, but only quadrilateral meshes (each quadrilateral consists of eight triangles) were demonstrated for the convenience of comparisons.

Figure 9 .
Figure 9. Evaluation results of non-homogeneous sweets objects with a top-side compression

Figure 11 .Figure 13 .
Figure 11.3D evaluation results of bacon object with a trapezoid shape (Eva-B-V-T).The estimated parameters used in both (a) and (b) are from 3D estimation results of trial Est-B-P and Est-B-V, respectively 

Table 1 .
Detailed information of all experiments trials

Table 2 .
The constitutive laws of generalized models

Table 3 .
Estimated parameters for trials with 'Est-' in the trial names

Table 4 .
3D mesh information and time cost for parameter estimation of bacon trials