Implicit Identification of Contact Parameters in a Continuous Chain Model

Accurate contact modeling is of great importance in the field of dynamic chain simulations. In this paper emphasis is on contact dynamics for a time-domain simulation model of large chains guided in a closed loop track. The chain model is based on theory for unconstrained rigid multibody dynamics where contact within the chain and with the track is defined through continuous point contacts using the contact indentation and rate as means. This paper presents an implicit method to determine contact parameters of the chain model through the use of none gradient optimization methods. The set of model parameters are estimated by minimizing the residual between simulated and measured results. The parameter identification is tested on four different formulations of the Hunt-Crossly hysteresis damping factor with the aim of recognizing a superior model.


Introduction
Contact models play a vital role in a variety of scientific fields.The importance of accurate, reliable and fast contact models is shown by the wide range of literature on the subject dating back to 1882 where (Hertz, 1882) formulated one of the first contact models.The complex mechanics of interacting bodies make contact difficult to model as kinematics, geometry, material and surface properties have to be taken into account.Therefore, most models presented in the literature are simplified through assumptions that quantify specific types of body interactions.Assumptions like, pure static contact, low impact velocity, pure impact, point contact, no plastic deformation and pure linear elastic deformation are often used as means to simplify the model and speed up the computational time.Today a widely used contact model is the Hunt-Crossly model.The Hunt-Crossly model evolves from impact theory to replace the linear spring and damper model, which in the contact period may introduce both an infinite force gradient and tensional force.The Hunt-Crossly model assumes contact of perfect convex geometries by which the reaction force is applied in a single point on the two contacting bodies.The Hunt-Crossly contact model is a continuous force model derived from the local indentation, δ and indentation rate, δ at the point of contact, see (Hunt and Crossley, 1975).The contact model evolves around the pure elastic contact model formulated by Hertz, as in which K and n are the stiffness coefficient and the power exponent computed from material and geometric properties in the local contact region of the contacting bodies.With the aim of achieving a continuous force contact Hunt and Crossly replaced the common linear viscous damper with a hysteresis damping factor λ and δ n , see Equation (2).A hysteresis shaped contact force is obtained, in this way by which the enclosed area of the contact force versus indentation corresponds to the energy dissipated in the period of contact.
Tensional force in the contact period is avoided by the Hunt-Crossly model as long as One main concern of the Hunt-Crossly model is the formulation of a proper hysteresis damping factor, λ in relation to the coefficient of restitution (CoR), e.A collection of different proposals for the hysteresis damping factor are presented in (Zhang and Shaft, 2009).The proposed formulations of λ both have approximate and exact relationship with e, leaving five different formulation of the Hunt-Crossly model.The presented models of the hysteresis damping factor have never been experimental validated as stated by (Zhang and Shaft, 2009), despite their extensive use in literature.Likewise, these models have never been compared to find the superior model.The identification of damping parameters in large models is often a complicated task as the dissipated energy within a mechanical system depends on the configurations of materials, structural design and internal contact.Therefore, damping in the contact model is often a joint term which covers several effects like internal heating, plastic deformation, viscoelastic effects and deformation wave.In some models damping is simply estimated by the trial-and-error method, see (Pedersen, 2006) and (Moreira et al., 2010).Although, several researchers have paid their attention to identification of damping parameters.In (Zhang and Shaft, 2009) and (Sondergaard et al., 1990) extensive experimental work on solid spheres bouncing off a flat plate is conducted by the use of high speed cameras to identify the pre and post impact velocity.Also, in (Labous et al., 1997) a more complex experimental setup of two colliding solid spheres is performed with emphasis on the frictional effects during impact.The resonance method is another popular method used to determine the system damping in more complicated structural designs, see (Ta and Lardiés, 2006).An experimental setup capable of isolating both stiffness and damping parameters from the system response in a microscale structure using a sophisticated experimental setup is presented in (Shi and Polycarpou, 2009).
Direct experimental methods for parameter identification are infeasible means in large and complex structures as the utilized equipment often has a limited range of use.Implicit parameter identification methods are therefore often used in large applications.An implicit method for model parameter identification of an induction motor is presented in (Ursem and Vadstrup, 2003).The exact parameters of the motor model are identified by minimizing the residual between measured and simulated current using a differential evolution algorithm as means.The use of parameter optimization also proves to be efficient in other well defined models like for example magnetorheological fluid dampers and magnetic hysteresis characteristics of construction steel, see (Kwok et al., 2006) and (Lederer et al., 1999).Implicit parameter identification method has also been used on contiguous multibody models.A gradient based optimization method for estimation of model parameters in a multibody vehicle model is presented in (Serban and Freeman, 2001).By minimizing the difference between measured and simulated accelerations, a reasonable estimation of both stiffness and damping parameters for the suspension system of a full scale multibody vehicle model is obtained.In this paper an implicit parameter identification method is presented.It is applied to a previously developed continuous model of large scale loop-sortingsystems (LLS), see (Sørensen et al., 2010(Sørensen et al., , 2011)).Damping parameters of the LSS model is estimated by minimizing the difference between simulated and measured data using optimization methods as means.The parameter identification is conducted on several LSS to account for any variation between the layouts.Furthermore, four of the models presented by (Zhang and Shaft, 2009) are tested to identify the superior hysteresis damping model of the chain dynamics.

Model of LSS
The main feature of the LSS is a closed chain of carts which runs in an enclosed track.Sorting is conducted by each cart in the chain which transports items from an induction point located along the track to one of several chutes also located along the track.Each cart in the chain consists of two components; a cart and a carrier.The cart interfaces to the track by four wheels and interfaces to the cart ahead and behind through a spherical plain bearing joint, see Figure 1.The carrier is mounted on top of the cart and consists of a tray which carries the item on top, see Figure 1.The length of the cart can be between 500mm and 1250mm depending on the size of the items.The closed track is designed by eleven different standard types, like for The only physical contact between the chain of carts and the track is via the four wheels mounted on the cart.The chain of carts is propelled by linear motors which are mounted to the straight tracks underneath the chain of carts, see Figure 2. The linear motors are distributed equally along the track forcing the chain of carts forward through a magnetic field applied to an iron core located at the bottom of the carts.The chain of carts normally operates at a speed between 1.5m/s and 3.0m/s.

The Multibody Chain Model
The LSS is modeled as a chain of kinematically decoupled bodies interacting with each other through a contact force formulation.Each cart is assumed to be rigid with a local body coordinate frame at each body's center of gravity, see Figure 3.Because of the spatial track layout a spatial chain model is used by which the cart i is defined in space by a position vector r i and an orientation vector p i .The enclosed track is assumed to be rigid with perfect aligned and tangent joints between each track element.The kinematically decoupled bodies lead to a set of independent equations of motions which are solved individually at each time step.The second order differential equations of motions are numerically solved utilizing the fifth order Cash-Karp-Runge-Kutta solver provided by (Press et al., 1996).Great emphasis is put in the development of a fast and robust chain simulation model.By means of which computation of the chain dynamics has become linearly dependent on the number of carts using only 4s to simulate the dynamics of a layout with 52 carts.

Points of Contact
As all deformations within the chain contacts occur in a local region they are assumed to happen at a single point.A local vector from the body reference frame defines seven contact points on each cart within the chain, see Figure 3.That is: • Front and rear joint, s f i and s r i .
• Left and right steering wheel, s lws i and s rsw i .
• Left and right driving wheel, s lrw i and s rrw i .
• Motor force, s m i .
The Hunt-Crossly model in Equation ( 2) is used in all contact points except from the force applied by the motor at point s m i which is computed by the speed controller.

Hysteresis Damping Models
The Hunt-Crossly contact model was developed in studies of impact contact to avoid the non-physical effects of the common linear spring and damper model, see (Hunt and Crossley, 1975).Effects like infinite force gradients and tensional forces in the initial and final part of the contact period are often found.Effects like infinite force gradients and tensional forces are avoided in the Hunt-Crossly model as damping is assumed be linearly dependent with the contact indentation.Energy dissipated within the contact period arises only from internal heating of the contacting materials.This is also observed from the model which does not take any permanent deformation of the interacting bodies into account.Studies of colliding spheres show that impacts of pure internal heating occurs only in a small range of e > 0.7 and v i < 0.5m/s, see (Zhang and Shaft, 2009), (Lankarani and Nikravesh, 1994) and (Lankarani and Nikravesh, 1990).Despite its defined limited range of use the Hunt-Crossly model has been used, with great success, in a wide variety of multibody models, see for example (Flores et al., 2006) and (Pedersen et al., 2004).Several different proposals for the Hunt-Crossly hysteresis damping factor may be found in literature.An overview of five different proposed formulations of the hysteresis damping factor is presented by (Zhang and Shaft, 2009), see Table 1.In which v i is the relative impact velocity between the colliding bodies.The contact models gathered by (Zhang and Shaft, 2009) all revolve around the formulation of λ in relation to e.In general two different approaches are used to define the hysteresis damping factor.The hysteresis damping factor in (Hunt and Crossley, 1975), (Lankarani and Nikravesh, 1994) and (Lankarani and Nikravesh, 1990) derive from an energy approach by comparing the work dissipated by the Hunt-Crossly model to the energy dissipated by the CoR.A relationship between λ and e is obtained in (Herbert and McWhannell, 1977), (Lee and Wang, 1983) and (Gonthier et al., 2004) by combining the equation of motion of the colliding bodies with the Hunt-Crossly contact model.Through different assumptions (Herbert and McWhannell, 1977) and (Lee and Wang, 1983) obtain a closed form solution whereas (Gonthier et al., 2004) presents the exact open form solution.The five hysteresis damping factors are all according to (Zhang and Shaft, 2009) derived scientifically cor-rect despite the different approaches used to obtain the relationship with the CoR.As stated in (Seifried et al., 2010) the CoR is a joined term based on different mechanisms of kinetic energy loss such as plastic deformation, viscoelastic material behavior and wave propagation in the contacting bodies.This further complicates the model as the CoR becomes a function of different variables like for example the indentation rate, body geometries and material properties.Different models of the CoR have been proposed in literature.Both linear and higher order models of the CoR in terms of the indentation rate are presented in (Seifried et al., 2010).A model for identifying the CoR taking the effects of wave propagation in different contacting structures into account are also presented in (Zhang and Shaft, 2009).The presented model requires information derived from a finite element model of the structures subjected to impact in order to deduct the CoR.
In the developed simulation model the CoR is defined as a constant and the formulation corresponding to Equation ( 4) to ( 7) can be utilized.The formulation in Equation ( 8) is avoided because of its iterative nature.HC, HW, LW and LN are used in the paper as abbreviated forms for the four hysteresis models in Equation ( 4) to (7), see Table 1.

The Parameter Identification Method
Most of the geometric and physical parameters of the chain model are simply found by using CAD software.
Wheel and joint stiffness are obtained through experiments, measuring the coherence between indentation and contact force.Coefficient K and n of the Hertz's contact model are obtained by using a non-linear regression scheme.Data of the rolling resistance of the wheels is provided by the wheel manufacturer.
To estimate other, more elusive, parameters an implicit optimization technique is presented.The elusive model parameters are considered as design variables in the implicit identification method from which the optimal solution is obtained by minimizing the difference between simulated and measured data.Different track layouts and cart configurations are considered to test the robustness of the parameter identification method.

The Implicit Parameter Identification Method
An estimate of the CoR is found by minimizing the difference between measured and simulated data.The parameter identification method is formulated in Equa-  Herbert and McWhannell (1977) λ Lee and Wang ( 1983) Lankarani and Nikravesh (1994) and Lankarani and Nikravesh ( 1990) tion ( 9) as an unconstrained optimization problem. Minimize Where f (x) is the implicit objective function of the design variables x and m is the number of design variables.The vectors a and b are the lower and upper explicit bounds of the design variables.The objective function of the parameter optimization problem is defined as the sum of squares of the normalized root mean square deviation (NRMSD) between measured and simulated data, as where u is the number of forces measured by the transducer cart and H j (x) is the NRMSD defined by in which I j (x) is the root mean square deviation (RMSD) and f m j is the vector of measured forces (subscript m).The RMSD is found by in which f m j,k is the measured force and f s j,k (x) is the simulated force and v is the number of points in the measured data serials.The use of experimental data and the highly non-linear chain model as objective function makes the chance of a continuously differentiable design space most unlikely.With the possibilities of a non-convex design space with several local minima it is chosen to utilize a non-gradient based optimization algorithm by which the genetic Complex method is used.The Complex method which is an evolutionary optimization algorithm was proposed by (Manetsch, 1990) and later on further improved by (Seifried et al., 2010).Tests in (Seifried et al., 2010) show that the complex method can handle up to around 36 design variables with a high chance of reaching the global optimum.The size of the complex population is recommended to m + 2 by (Seifried et al., 2010).

Design Variables of the Parameter Identification
The main purpose is to identify the best set of damping parameters of the contact model.This introduces four design variables as the CoR of the following contacts: • Cart joint -axial ( η ) direction, e JA .
• Driving wheel vs. track, e RW .
The explicit bounds of the CoR are set to 0.2 ≤ x ≤ 0.9 throughout all conducted optimizations.These design variables are independent of the chosen models, Equation ( 4) to (7).Besides the CoR there are also a number of layout dependent variables defined.They are: • Sorter speed reference value, v s .
• Pretension of the chain, d.
• Trigger point from which the data is logged, s M .
• Rolling resistance of the wheels, m d .
• Proportional gain of sorter speed PI controller, K P .
• Integral gain of sorter speed PI controller, K I .
The layout dependent design variables are individually defined by each measurement in order to get the exact time domain of the simulation model.Likewise the explicit bounds are individually defined.The wheel friction is included in the parameter optimization, as the individual track layout may be worn and cause changes to the rolling resistance.The speed control parameters of each layout are hard to extract from the real system due to the complex higher order system.A simplified model is then used in the continuous chain model by which the chain speed is controlled through the applied motor force.The vector of design variables are hereby defined as x = [e JA , e JR , e SW , e RW , v s , d, µ d , s M , K P , K I ], (13) in which ten design variables are used in the conducted optimization.

Test Layouts
Seven different test layouts are used in the parameter optimization, see Table 2.The test layouts all vary in track length in order to get a large representation of feasible LSS layouts.The number of level changes are also well represented.A top view of test layout 2 is illustrated in Figure 4.The two parallel closed lines illustrate the track width and the layout shape.Lines within the track illustrate level changes in the track layout and arrows along the track show the driving direction of the chain of carts.Figures along the track indicate the distance from the starting point of the conducted measurements.Similar drawings of the remaining six test layouts are found in Appendix.The configuration of the carts in the seven test layouts are differentiated.Three different wheel types have been used, denoted type A to C as shown in Table 2.The three wheel types all have different elastic properties which are accounted for in the chain simulations, see Figure 5. Likewise two different types of joints between the carts are utilized in the test layouts.Both joint type A and B are constructed as spherical plain bearings although they have significant different elastic properties, see Figure 6.
The different cart configurations also induce different configurations of the item carrier.This entails different mass properties of the carts which of course is accounted for in each simulation.Different pretension and speed of the chain is commenced to test the chain simulation model against the 0 0.5 1 1.5 2  2 shows the number of the conducted measurements on each layout.

Measurements of Chain Forces
A transducer cart has been developed in order to measure forces in the chain of a LSS, see Figure 7.The transducer cart is adapted with the same geometric interface and the same flexibility as a normal cart.
Eleven different sensors are fit into the transducer cart making it capable of measuring forces in the front joint in ξ-, η-and ζ-directions, force on the steering wheels in η-direction, force on the driving wheels in ζ-direction and the rotational speed of all four wheels.All force sensors measure shear forces using strain gauges to measure shear strain in a thin web design.Variations in temperature are compensated utilizing a full Wheatstone bridge on each web.Wheel speed is computed from inductive sensors providing 15 pulses per revolution.Data is collected using two Spider8 from Hottinger Baldwin Messtechnik GmbH.The data is continuously logged on an onboard laptop.All equipment is mounted on top of the cart along with a power supply which provides 12V DC to the two data logging units.Measurements are conducted by replacing one cart in the chain with the transducer cart.A photocell mounted on the cart is used to trigger start and stop data logging as it hits a well defined trigger point along the track layout.Experience from the different conducted measurements shows that a sampling fre-

Results
The presented parameter identification method is conducted for all 22 measurements for each of the four hysteresis contact models.
For clarity the presented results are divided into three parts.First results from the presented parameter identification method are presented.Second, results form the comparison of the four hysteresis contact models are presented.Finally, an analysis of variance of the CoR is presented.

Parameter Identification
Due to the extensive amount of results gained from the performed parameter identification only some selected results are highlighted in the paper.The best results are in general achieved in measurements with low velocities and high pretension of the chain.With high pretension of the chain the steering wheels of the carts are in constant contact with the track.Hereby, a loose and constant shifting contact between left and right steering wheel is avoided as this may introduce high impact forces in the chain simulation model.Results form the LW contact model are shown in Table 3.
Fractions of the initial guess and the final joint forces of the parameter identification is shown in Figure 8 and 9. Figure 8 shows the joint force in ξ-direction for measurement 1 with HC contact model and Figure 9 shows the joint force in ξ-direction for measurement 18 with contact model HM.The two dotted lines illustrate the initial guess of the chain dynamics and the optimum results.The initial set of design variables for the parameter optimization is defined through manual interpolation before commencing optimization.The manual interpolation induces that some initial guess of the design variables are very close to the final result of the parameter optimization.Minor improvement is achieved for measurement 1 as the NRMSD is reduced from 21% to 18%.Whereas improvements are more significant in measurement 18 as the NRMSD is reduced from 32% to 13%.The delay observed between the measured date and the optimal simulation in Figure 9 result from several circumstances.One reason is the small geometric variation within the chain and track which can introduce local delays along some parts of the track.Also, the formulation of the objective function emphasizes minimization of the largest residual between the seven forces.Minimizing the largest force residual may flavor a delay of some of the other forces.A common iteration history from the conducted optimizations is shown in Figure 10.The iteration history shows the reduction in the objective function on test layout 5, measurement 13 for the four contact models.Most improvements are gained in the first 100 iterations.Around 20% of the conducted optimizations stopped at the maximum number of 5000 iterations.As shown in Figure 10 no significant improvements over the last 2000 to 4000 iterations are obtained.A perfect match between simulation and measured data will in accordance with Equation (10) make the objective turn zero.As shown in Figure 10 the perfect match is not reached in the iteration history.It is believed that the optimum reached in the presented results is close to the global optimum and the reason why the perfect match has not been reached is the assumptions made in the multibody model and the lack of flexibility in the chosen design variables.Other optimization algorithms might prove superior in finding the global optimum although the Complex method is known for its robustness in handling non-convex problems.
The square term of the objective function in Equation (10) forces the optimization routine to reduce the largest residual of the seven forces.This is also shown from the initial and final size of the NRMSD which tend to reach the same level at the optimum, see Figure 11.
The fact that the parameter optimization finds an improved solution proves that an enhanced set of parameters of the chain simulation model can be identified.The implicit parameter identification method provides an estimate of the model parameters for the single LSS track layout by which it does not necessarily represent the parameters for all LSS layouts.An investigation of an overall set of parameter are further investigated by analyzing the variance of the CoR between the conducted measurements.

Comparison of Contact Models
The superior hysteresis damping models tested are identified by comparing the final objectives of the conducted parameter optimizations.A box-plot of the best objectives for the four contact models is shown in Figure 12.Both the medians, the 25th and 75th percentiles are very much alike in the four tested hysteresis damping models.Results from a one-way ANOVA test show that the means of the objectives for the four contact models are the same.Hence, none of the four contact models are statistically better than the other.Although, the objective function for the HW contact model is some magnitudes smaller than the other models.The median and the standard deviation are smaller for the HW model and the extreme data is located closer to each other.

Variance of the CoR
From the conducted parameter identifications an estimate of the average set of damping parameters is acquired.This is done by analyzing the results of the conducted parameter optimizations.Data from the conducted optimizations are divided into groups of the four contact models.Each group is then further divided according to joint or wheel types.The group with joint type B only contains 3 samples and is not included in further discussion.To further investigate the mean CoR for wheel type B additional experimental data should be provided to have sufficient statistical material.A box-plot of the CoR in the axial and radial joint contact is shown in Figure 13 and 14 for the four contact models.The dotted lines indicate the upper and lower explicit bound of the parameter optimizations.Each contact type contain 19 sample points on 6 test layouts, see Table 2.
Both axial and radial CoR vary insignificantly.Especially contact model LN and LW have large variation in the size of e with extreme point in both ends of the explicit bounds.The 25th and 75th percentiles of the HW and HC contact models are within a sufficient variation.A fine approximation of the wheels CoR is obtained.In the presented results wheel type B and C are not included since these groups only contain three samples each.A box-plot of both the driving and steering wheel CoR for wheel type A is shown in Figure 15 and 16.Each contact type contain 16 sample points on 5 test layouts, see Table 2.
The CoR median of the four contact models are almost the same around 0.47.Although, the variation of the CoR is large for contact model HC and LN.The 25th and 75th percentiles of LW are fine whereas HW is significantly smaller than the other three models.The variation of the CoR is in general smaller for the wheel contact than in the joint contact.The difference observed between the wheel and joint contact can not be explained although it indicates that further improvement of the chain model can be gained by improving the joint contact formulation.
The HW contact model proves to have the smallest variance of all four CoR.Especially the CoR of the two wheels have a small variance, see Figure 15 and 16.This shows that HW supplies the best accordance of the chain contact forces which is in accordance to the box-plot of the objective functions, see Figure 12.

Conclusion
An implicit method for identification of model parameters has been proposed.The parameter identification is defined as an optimization problem minimizing the residual between measured and simulated forces.The proposed parameter identification method has been tested on a complex continuous chain model with the aim of identifying contact damping parameters.The presented method has been tested on 22 measurements on 7 different LSS layouts.Apart from the four design variables dedicated the CoR six layout dependent design variables are utilized to synchronize the time variance.Reasonable results have been reached for the individual measurements with residuals less than 10% for some simulations.The fine individual accordance shows that the presented method is an efficient way to estimate complex model parameters in large multibody models.
Analyzing results between the 22 measurements show some variation in the CoR.The exact cause of this variation is not fully identified.Although, two plausible causes may be recognized: • The global solution has not always been reached due the non-convex design space of the optimization problem.
• The utilized continuous chain model is insufficient to simulate all variations in the complex dynamics of LSS.
Theses issues might be further treated to get an enhanced accuracy.Although it is believed that the presented results provide a sound estimate of the average CoR for a general purpose use.
The Hunt-Crossly model proves to be accurate and efficient as shown in previous work on the dynamic chain model, see (Sørensen et al., 2010(Sørensen et al., , 2011)).Although, from literature, other models of the hysteresis damping factor may hold a superior formulation.Four different proposals of the hysteresis damping factor has therefore been tested in the continuous chain model.Tests reveal that none of the proposed contact models have a significant superior accuracy.Although, from the commenced parameter identifications the HW hysteresis damping model indicates a minor variation within the objectives and CoR.

Figure 1 :Figure 2 :
Figure 1: Main components of the cart design.

Figure 3 :
Figure 3: Points of contact on each cart within the chain.

Figure 5 :
Figure 5: The elastic properties of wheel type A, B and C.

Figure 6 :
Figure 6: The radial elastic properties of joint types A and B

Figure 7 :
Figure 7: The left hand side shows the design of the transducer cart.To the right is shown the transducer cart as mounted into a chain of carts.

Figure 8 :
Figure 8: Joint forces in ξ-direction for measurement 1 with contact model HC.

Figure 9 :
Figure 9: Joint forces in ξ-direction for measurement 18 with contact model HW.

Figure 10 :
Figure 10: Iteration history of measurement 13 on test layout 5.

Figure 11 :
Figure 11: Initial and final result of measurement 8 with the LN hysteresis contact formulation.

Figure 12 :
Figure 12: Box-plot of the final objectives for the four contact models.

Figure 13 :Figure 14 :Figure 16 :
Figure 13: Box-plot of the CoR of the axial joint contact for the four contact models.

Table 1 :
Formulations of the hysteresis damping factor for the Hunt-Crossly contact model.

Table 2 :
Main data on the seven test layouts.

Table 3 :
CoR from parameter identification with contact model LW.
quency of 300Hz is sufficient in order to provide accurate results.