DEFORMATION OF PEACHES SUBMITTED TO CYCLIC LOADING USING THE DISCRETE ELEMENT METHOD

This study aimed to use and validate a simple contact force model that describes the behavior of peaches stored in a wooden box, loaded cyclically, by using the discrete element method. The Kelvin-Voigt contact force model was used and the cyclic loading in the simulation was performed with an amplitude of 1 mm and frequency of 12 Hz for time intervals of 3600, 5400, and 7200 s. The laboratory test was performed with a peach box attached to a vibrating table by applying the same excitation conditions used in the simulation. For validation, model results were compared with those from the laboratory test by using deformation values of fruits in the contact regions. To measure this deformation, we adapted the methodology for determining the radius of curvature. The Kelvin-Voigt linear contact model, despite its simplicity in representing the viscoelastic behavior, adequately estimated the final deformation in peaches submitted to cyclic loading.


INTRODUCTION
To understand the behavior of fruits submitted to mechanical stresses it is essential to know their response to these stresses. Mechanical damage in fruits and vegetables is the rupture of the cellular structure caused by mechanical interactions such as impact, cutting, abrasion, and compression. Many studies found in the literature (Li, 2013;Tezotto et al., 2011;Sousa et al., 2017;Komarnicki et al., 2016;Stropek & Golacki, 2013;Fadiji et al., 2016) carried out analyses in the laboratory and related the physiological behavior observed with the applied force or tension as the final objective.
In general, studies associate damage caused by vertical vibration (transportation) with box position in the pile and pile position inside the truck (Barchi et al., 2002;ZHOU et al., 2007;Dantas et al., 2013). The study of Van Zeebroeck (2007) was the only one found in the literature that predicted the damage extension in spherical fruits conditioned in the box and submitted to cyclic loading, which evidences the lack of knowledge about fruit behavior at the time vibration occurs.
The discrete element method (DEM) is based on particle-particle or surface-particle dynamic interaction analysis, reducing the problem to integrating Newton's translational motion equation and Euler's rotational motion equation (Haug, 1992). Interactions between particles require a force model that describes how the interactions occur, being important to determine the time interval between successive integrations since it is directly related to the simulation time accuracy. DEM has application in solving different problems involving particle behavior from different areas such as simulation of concrete failure (Skarzynski et al., 2015) and concrete flow (Remond & Pizette, 2014), failure propagation in rocks (Trivino & Mohanty, 2015) and study of soil movement (Asaf et al., 2007). This method is also used in studies of grain post-harvest (Lu et al., 1997;Sakaguchi & Favier, 2000;Raji & Favier, 2004;Coetzee & Lombard, 2013;Hou et al., 2014;Leblicq et al., 2016).
In this study, we propose to use and validate a linear contact force model that describes the behavior of fruits conditioned in a box, loaded cyclically, by means of DEM.

MATERIAL AND METHODS
For the simulation, we used a wooden box containing 24 fruits arranged in two symmetrical layers of 12 fruits (Figure 1). Peaches of the variety BRS Rubimel, commercially classified as caliber 4 (diameters of 56 to 61 mm), was used in the experiment. In the vertical direction, three contact points were considered for assessing the simulation: PCT1, the contact point subject to the mass of both layers; PCT2, the contact point subject to the mass of one layer; and PCT3, the contact point of the upper layer ( Figure 1).
The software YADE ® was used for the simulation with DEM since it is a free software. Some characteristics of fruits and box, such as dimension, density, coefficient of friction, and parameters for the contact force model, were necessary to describe the bodies. These characteristics, as well as the type of model to be used, interactions that must occur, iteration time, and excitation to be applied in the set, were declared in the YADE ® input algorithm.
For the physical characterization of the wooden box used in the validation test (a virtual box with the same characteristics was used in the simulation), its density, calculated in 514.18 kg m −3 , was determined by using a body of proof with dimensions of 27.5 × 10.75 × 8.4 mm. Since the rigidity of the box material is much greater than peach rigidity, we considered the box material as purely rigid and its coefficient of rigidity was determined in a uniaxial test, being found a value of 9840 N m −1 . The coefficient of friction of 0.837 was obtained between two wood surfaces through a characterization table, using their densities and coefficients of rigidity as inputs.
Wood-wood coefficient of friction was determined because the software applies a coefficient of friction for each material involved in the interaction. Thus, it is necessary to know the wood-wood and peach-peach coefficients of friction so that the software simulates the interaction between materials in a correct manner.
In the simulation, spheres represented the fruits. Laboratory tests were performed in triplicate for time intervals of 3600, 5400, and 7200 s. Sphere diameter was determined from the average diameter of peaches that occupied the same position at the same time interval.
Peach density was calculated from the average values obtained by using an analytical balance (precision of 0.0001g) and the water volume displaced (precision of 1 mm) when individually emerging each fruit in a beaker.
The device used for defining the coefficient of friction consisted of a rectangular structure with dimensions of 200 × 180 mm and 60 mm high with bearings supported by a second rectangular structure of 200 × 180 mm and 65 mm high. Both structures were constructed in acrylic and arranged in such a way as to allow the shearing of their material. The force required to move the product inside the acrylic structures was determined from both rectangular structures filled with peaches and the addition of mass in a container at the end of the cable. The mass of peaches and that added to the container was measured. Based on the Coulomb's Law, the coefficient of friction was determined by [eq. (1)]. (1) Where, F is the force of friction (N); µ is the coefficient of friction, and N is the normal force.
The Kelvin-Voigt model (Equation 2), which is a simple linear viscoelastic model, was used as a model of contact force.
(2) Where, F is the force (N); k is the coefficient of rigidity (N m −1 ); b is the damping coefficient (N s m −1 ); δ is the deformation, and is the strain rate (m s −1 ).
In order to determine the coefficients of rigidity and damping, which are fruit mechanical parameters represented by the Kelvin-Voigt model, a fluency test was carried out in a cylindrical specimen extracted from peach pulp measuring 12 mm in height and 10 mm in diameter with normal loading of 2.5 N. The data obtained from this test were used as input to the software MATLAB ® , as well as the equation of the model found in Mohsenin (1986) (Equation 3) for curve fitting and coefficient obtaining. These values were used as input data in the software YADE ® . For calibration, the fluency curve obtained in the software YADE ® was compared with the curve constructed in the software MATLAB ® with the test data. The input values from YADE ® were corrected until both programs presented overlapped curves.
The characteristics found for fruits of peach were the sphere diameter, which varied according to the product diameter, a density of 132 kg m −3 , a coefficient of rigidity of 1044 N m −1 , a damping coefficient of 38750 N s m −1 , and a coefficient of friction of 0.512.
The software requires that the interactions occurred in the simulation be declared in the algorithm. In this case, the interactions are between the sphere-surface of the box (horizontal and vertical directions) and sphere-sphere.
In this study, time interval was 8 × 10 −6 s because it was the longest time the software could reproduce the behavior exhibited by the fruit material in the laboratory tests.
A vertical excitation was applied in the set with a frequency of 12 Hz and an amplitude of 1 mm for periods of 3600, 5400, and 7200 s (t60, t90, and t120, respectively).
To obtain sphere deformations in the specified contact regions, the final positions of spheres were registered and, by difference from the initial positions, deformations at points of contact were calculated. The deformation values obtained in the simulation were compared to the corresponding values obtained in the laboratory test.
For the validation test, we used a vibrating table with vertical oscillation, fixed amplitude of 1 mm, and variable frequency. A wooden box with fruits arranged similarly to the configuration used in the simulation was placed on this table.
Peaches were arranged to avoid contact points in the apical and peduncle insertion regions. In order to assess fruit rotations, a line was painted in the middle part of each fruit, representing the intersection with vertical planes passing through the center of each peach, being used to align the fruits (Figure 3).
The tests were carried out randomly, in triplicate, at a frequency of 12 Hz and duration of 3600, 5400, and 7200 s.
After the test, fruits were numbered according to the layer (layer 1 at the bottom and layer 2 at the top) and their position in the layer (Figure 2) to facilitate measurement. We waited 24 hours to carry out the measurement in order to visualize the darkening by the release of enzymes and facilitate the measurement. For deformation measurement, we adapted the Mohsenin (1986) methodology for determining the radius of curvature. Based on the use of a clock comparator and fixed rods equally spaced from clock rod and knowing these distances (points A, B, and C of Figure 4) and having only the measure of the clock comparator as a variable (point D of Figure 3), the radius of curvature of the product can be determined. We measured at pre-and postdamage point of contact (points D and D' of Figure 3) and calculated the difference of measurements to determine fruit deformation (D in Figure 3). The difficulty of this procedure lies in performing the pre-damage measure without knowing the point of contact. Thus, a study of diameter variation was conducted in 45 randomly selected peaches from the test batch. For this, the diameter of each peach was measured in its equatorial region, not considering the suture, making three measurements spaced 10 mm in diameter and calculating the mean and coefficient of variation. The average diameter found was 59.11 mm and the coefficient of variation was 3%, which is considered low. For this reason, we assumed that the diameter value at the exact location of the point of contact before the test is identical to the diameter value at the edge of the deformed region due to cyclic excitation. Thus, measurements were taken at the center of the deformed area (post-damage measure) and at the edge of the deformed region (equivalent to the pre-damage measure). The difference between pre-damage and post-damage measures, obtained in the clock comparator, provided the product deformation from the test. The measuring device was built in the Laboratory of Prototypes of the Faculty of Agricultural Engineering of UNICAMP.
The simulation validation was performed by comparing the deformation values obtained from the simulation with the deformation values obtained from the test.

RESULTS AND DISCUSSION
A simple linear model was selected since we verified the lack of knowledge of the response of these models applied in spherical fruits loaded cyclically, in addition to being easy models to determine the parameters.
Moreover, Van Zeebroeck (2007) used the model of Kuwabara and Kono and had a great difficulty in obtaining the parameters due to the model complexity.
Simulation validation was performed by comparing experimental and simulated deformation values in order to identify whether the response of the points of contact obtained with the simulation followed the same trend of the real points of contact. Van Zeebroeck (2007) also assessed the trend of simulated behavior in comparison to the real, i.e. how close the simulated behavior reached the real.
In general, peaches presented no significant damage after the tests, except for some fruits with a higher degree of ripeness. This behavior suggests that peaches submitted to the applied frequency value, along with the test duration times, would not be visually affected, being suitable for commercialization.
For each point of contact in both layers, a pair of deformation values was obtained, one of them experimentally and other from the simulation. The results show how close are the simulated deformation values from those obtained experimentally, regardless of the position occupied by the spheres or peaches in the box. In addition, graphs were generated with the 12 events, i.e. 12 experimental-simulated pairs obtained for each point of   Both experimental and simulated deformation values for PCT1 presented the same order of magnitude. The experimental deformation values varied from 0.29 to 1.16 mm at t60, 0.14 to 1.23 mm at t90, and 0.36 to 1.42 mm at t120, while in the simulation this variation ranged from 0.56 to 0.96 mm for all time intervals.
Among the pairs that presented the lowest deformation values (events from 0 to 6), the highest difference between test and simulation occurred for t90, with a value of 0.42 mm, while among the highest deformation values (events from 7 to 12), the highest difference was observed for t120, which presented a value of 0.5 mm. Because the highest difference between them was 0.5 mm, and considering an average diameter of the sample peaches of 59.11 mm, this difference corresponds to 0.84% of the diameter, which is low and evidence that the deformation values predicted by the model satisfactorily represented those values obtained in the test for PCT1. Percentage differences in this order of magnitude (less than 1% of the diameter) can be considered as irrelevant when considering a biological material and its intrinsic variability.
Similarly, Figures 7 to 9 show the results for PCT2, the point of contact submitted to the mass of one of the layers.  In this case, test values varied from 0.23 to 1.22 mm at t60, 0.24mm to 1.55 mm at t90, and 0.21 to 1.96 mm at t120. In the simulation, variation ranged from 0.42 to 0.73 mm. At the lowest deformation values, the highest difference was 0.21 mm. Thus, we can consider that for the lowest values, the simulation was able to reproduce the behavior of the test. When analyzing the events from 7 to 12, three points obtained differences higher than 0.59 mm (0.82 mm at t90 and 0.64 and 1.23 mm at t120), exceeding 1% of the average diameter, which could be considered as an error. However, this value is still lower than the 3% variation of the sample diameter. PCT2 presented the highest differences in the events from 7 to 12 when compared to PCT1. On the other hand, the events from 0 to 6 presented closer values in PCT2 when compared to PCT1. Considering that only three points exceeded 1% of the diameter, the model was able to predict satisfactorily the tests in PCT2.   In PCT3, the values obtained in the test ranged from 0.15 to 0.87 mm at t60, 0.10 to 0.93 mm at t90, and 0.21 to 1.83 mm at t120 whereas, in the simulation, these values ranged from 0.28 to 0.5 mm. As for PCT2, in PCT3, for events from 0 to 6, the model satisfactorily predicted the observed experimental behavior since the value of the highest difference between test and simulation was 0.18 mm. For events from 7 to 12 at t60 and t90, no differences higher than 0.5 mm was observed between test and simulation. However, at t120, the differences in events from 10 to 12 reached a range from 1 to 2.5% of the average diameter (0.73, 1.07, and 1.33 mm), but lower than the diameter variation, i.e. they present small differences considering acceptable variations from 15 to 25% for biological materials (Pimentel-Gomes, 2009). Thus, the model was able to reproduce adequately the deformation response reached in the test, as at the other points of contact.
Considering the events shown in Figures 4 to 12, we can affirm that 67.6% of them presented a difference of up to 0.2 mm when comparing experimental and simulated deformation values. This difference represents 0.34% of the average fruit diameter, which is a value considered as acceptable once that damages in fruits with deformations in this order of magnitude will hardly be perceived. These results, in addition to the natural variations of the biological material (peach pulp), show that simulation, in general, presented a good result, with values within the expected and close to the deformation values found in the test, validating the model. Another way to assess the similarity between simulated and real results is through the coefficient of determination (R 2 ). Considering the results obtained in the laboratory test as a base (x-axis), we determined the R 2 of the simulated results (Table 1), reinforcing the satisfactory values presented in the simulation. When stipulating differences higher than 0.2 mm as underestimated and overestimated between the deformation values tested and those simulated, the linear model used, in general, overestimated 9.3% and underestimated 23% of the events. In this case, the event 12 shown in Figure 12 was underestimated in 72.7% or 1.33 mm. Van Zeebroeck (2007) used the more complex Kuwabara-Kono model and obtained values overestimated in 2.5 mm, which is equivalent to variations of 100 to 300% in apples and tomatoes. Considering that the model is an adequate predictor, the spatial distribution of deformations is presented considering the position of spheres and peaches inside the box for a time interval of 7200 s.
Due to the different cycling loadings in the points of contact, we expect that PCT1 have the greatest deformations because it is subject to the mass of two fruit layers, PCT2 presents lower values due to the loading of a single layer, and PCT3 presents even lower values due to vertical free movement of the upper layer. This expectation is confirmed by the results found in the simulation shown in Figure 13, which represents the same arrangement shown in Figure 2 and numbered from 1 to 12 in Figure  13.
FIGURE 13. Spatial distribution of deformations obtained in the simulation for t120. Figure 14 shows the same spatial distribution of the deformation response obtained in the experimental test for the same time interval of 7200 s. The results shown in Figures 13 and 14 did not evidence any trend in the deformation values regarding the possible interactions between peaches or spheres between themselves or with the box sides. This could be attributed to a small influence of horizontal forces and friction at the used loading levels.
FIGURE 14. Distribution of deformations obtained in the test for t120.