ARTICLE INFO

ABSTRACT

1. INTRODUCTION
In the human gait analysis, three different foot-rocker mechanisms are known as heel rocker, ankle rocker and forefoot rocker, which are generated from pendulum dynamics with three fulcrum points of heel, ankle and toe [ 1 ]. In a gait cycle, there are two phases as swing stance phases and particularly in the human gait analysis, the stance phase is decomposed to the initial contact for touching with the ground from the heel, the loading response for shifting of the center of mass from back to front, the mid-stance for lifting the heel up, the terminal stance for transition to the swing phase, and then pre-swing [ 2 ]. Such as bipedal gait cycle can be simplified to a mechanical system based on a coupled pendulum as demonstrated by McGeer [ 3 ], called a passive dynamic walking robot without any electromechanical actuator to mimic three rockers. In this case, the compliant human leg dynamics has been studied with a series of contact points on the rolling surface [ 4 – 6 ] and with the effective roll-over geometry of leg motion [ 7–9 ]. The fact suggests that the utilization of the simple skeleton model provides a large benefit to compensate the limitation of experimental measurements in the gain analysis to realize the actual phenomenon of the ground reaction force, because there are evidences in the walking condition [10,11], such as demonstrated two peaks in the temporal evolution of the ground reaction force during the stance phase, while there is a lack of evidences in other gait patterns. In robotic applications inspired from biological mechanisms, an articulated leg with hoof was inspired from the horse [ 12 ] and a linkage mechanism to mimic the horse leg trajectory by Batbaatar and Wagatsuma [ 13 ] and sim- ilar approaches were found in closed-loop linkage mechanisms [ 14,15 ]. In the animal locomotion, there are various gait patterns known as walk, trot, gallop and so on [ 16 ]. Therefore, the estab- lishment of the detail analytical method for the ground reaction force is highly important to elucidate the principle to differ gaits in the sense of the energy consumption. Focusing on the ground con- tact phenomenon, the hoof part of the horse leg mechanism can be modelled in the form of the triangle structure to connect three rep- resentative points as the toe, heel and ankle, which reproduces the rocker mechanisms. If it is possible to provide kinematic/kinetics of the hoof-ground contact, the viscoelastic effect is also considerable to analyze how much the leg absorbs the impact force when touch- ing with the ground. It realizes the detail analysis of the amount of the energy storage discussed in the distal leg movement and landing motions [ 17,18 ]. The purpose of the present study is to establish the basement analytical method for the ground reaction forces in gaits and demonstrate the efficacy of the method in the simplified leg model to focus on the hoof-ground interaction. This paper is divided into following sections. Section 2 introduces the general Multibody Dynamics (MBD) formulation applied to the proposed leg mechanism as well as the viscoelastic contact force introduced to simulate realistic hoof-ground contact. Section 3 contains the kinematic and kinetic analysis focusing on the hoof–ground inter- action, while Section 4 summarizes the main results and discusses the potential improvements and limitation of the study.
2.  MATERIALS AND METHODS
2.1.  Multibody System Formulation for
Horse Leg Mechanism
The initial proposal of horse leg mechanism was introduced in the article by Batbaatar and Wagatsuma [ 13 ] and its kinematics/ dynamics formulated with MBD. Mathematical interpretation and initial state for horse leg mechanism can be written by generalized coordinates q as
The vector of q contains 42 (= 14 × 3) elements which are 14 rigid links and their center of position x , y , and orientation φ obtained in the present analysis. A set of kinematic algebraic con-
where the first 41 elements of the column matrix of kinematic con- straint equation are derived from the absolute constraints between body and fixed ground node. The last element defines the driver constraint (Crankshaft) of the proposed leg mechanism. The par- tial derivative of kinematic constraint equation respect to the gen- eralized absolute Cartesian coordinates q is Jacobian matrix Φ  is obtained as
where it allows us to investigate placement, velocity and accelera- tion analyses kinematically. The forward dynamics analysis intro- duces the mass matrix M = (42 × 42), and the generalized external force vector h (*) = (42 × 1), as follows:
where m, is the mass of rigid link to point i, J, = 21/3 is the polar moment of inertia of rigid link to point i , and g is the gravitational acceleration. The equation motion of the system for the computer system analysis can be is expressed in general matrix form as
= 0,
where ¨ is the acceleration vector, h ( α ) denotes the generalized external force vector, λ is the Lagrange multipliers with vector Φ ij = −(Φ q) x q − 2Φ n q − Φ n = γ . The parameter α and β are the feedback compensator in Baumgarte stabilization method [ 19 ] used to guarantee accuracy of the numerical solution. In the for- ward dynamic analysis, new coordinates and velocities require two arrays for ¶  and ¶ for the time step t + Δt as
ics analysis requires the initial state of the target mechanism in order to proceed numerical simulation in next time step.
2.2.  Contact Analysis Procedure in
Planar Mechanical System
Subject to the Impact
When describing the contact impact event between the hoof and ground, we need to determine the state of the contact and 80 accurately acquire the relative deformation and speed at the instant of contact. In Figure 1 , non-contact scenario and representation of viscoelastic ground contact model with variable damping if contact occurs in the system were shown in the kinematic model. When the leg contacts with the ground, the deformation or penetration is estimated as
where the P is the normal component of the position vector or the potential contact point as K = [ K , K ] (toe-first) or H = [ H , H ] for the case of heel first contact. The ground plane γ c is consistently levelled from the lowest point of the locomotive trajectory in the sense of the relative indentation as
In Equation (11) , the positive value of δ is that distance represent a separation, while negative values denote relative penetration of the contacting bodies. Therefore, the sign of penetration indicates the phase transition from swing to stance and vice versa for the case of hoof–ground contact. By using Equation (10) , it is clear to discrim- inate walking phases and contact forces as continuous function of deformation. Based on updated position and velocity of potential contact points on the hoof K = [ K , K ] or H = [ H , H ] from the multibody kinematic analysis, relative normal and tangential veloc- ity among the contact points are estimated as
where t is the unit vector in tangential direction obtained by rotat- ing unit normal vector n counter-clockwise by 90°. When bodies come into contact, the normal contact force can be estimated with contact force model introduced by Lankarani and Nikravesh [ 20 ]. In the damping factor model, the crucial parameter known as coef- ficient of restitution was introduced which defines the whether contact is fully elastic cr = 1 (or plastic) cr = 0, resulting impact force model can be written as
where cr denotes the coefficient of restitution which is defined as the ratio of relative approach velocity ♂ ( * ) and relative departing velocity of  δ (-) . The tangential friction force in contact is calcu- lated using modified Coulomb friction law [ 21 ]. The dynamic fric- tion forces F , in the presence of sliding [ 22 ] can be written as
where μ , the kinetic friction coefficient. To take the friction force effect into account, the direction of the tangential velocity vector c , = | ν , | with the target body velocity ν , which is opposite to direc- tion of motion. c , is the dynamic correction coefficient given by Ambrósio [ 21 ]. If the normal force F_ is obtained from a contact with a force-deformation [Equation (14)] of contact force model, a logical point-to-point spring-damper element [ 23 ] is only active during the period of contact in normal and tangential direction as shown in ground contact model in Figure 1 . Finally, absolute values of the force vector in normal and tangential defines the resultant ground  reaction  forces   GRF = exerted when the hoof F + F ( press against the ground. Numerical solution of equation of motion as well as contact force estimation based on state variables can be treated in one computa- tional loop with the iterative computational procedure as shown in Figure 2 .
3. RESULTS
According to the MBD formulation above, the constitutive con- tact force model is calculated numerically with kinematic and dynamic analyses, which allow to visualize temporal evolutions of the locomotive trajectory, resultant ground reaction force of the leg mechanism. Table 1 displays the parameters used in the numerical simula- tion for contact force analyses of horse leg mechanisms. Matlab- based numerical simulation was performed with a combination of the Euler method with the time step of 1 × 10 −3 s. For the contact analysis to prevent an unnatural rise of ground reaction force, the stiffness is no smaller than 2 × 10 4 N/m for the contacting bodies is considered according to simulation in Ristow [ 24 ]. It is simply that the ground deformation is equivalent to the resultant motion of the stiffness-damping system under the compression. Two tra- jectories were selected as typical trajectories to reproduce walking and running gaits according to the mechanics of the simple model, which are generated by changing the control parameters associated with driver constraint in the horse leg mechanism. The trajectories shown in Figure 3 , the 1st trajectory considered to be a walking in term of duty factor and step length, evaluation was done in recent study by Batbaatar and Wagatsuma [ 25 ]. The 2nd trajectory which has longer step length compared to the 1st trajectory considered to be running gait in which leg orientation controlled by swinging motion with respect to the body, which may potentially generate the propulsive ground reaction force due the intrinsic high rate of angular oscillation at the hoof. In the temporal evaluation of the hoof angle as shown in Figure 4 , walking and running gaits were generated in accordance with the first and second trajectories in
(a), running (b). Figure 3 . In the waling gait, the transition of sub-phases commonly occurred in the rolling motion of the triangle structure in both cases. In the walking gait, those sub-phases from heel-strike, flat and toe-off were gradually generated (Figure 4a top), which is con- sistent with the gradual change of the hoof angle ( Figure 4b top), while in the running gait, the transition started from the toe-strike that causes the roll-over sub-phase quickly and then reached to the toe-off sub-phase with a large gap from the ground rather than that in the walking gait. The result was numerically examined as the first contact occurs at instant time t, = 0.267 s and α, = −6.83° leaves the ground at t, = 0.569 s and α  = −20.9° for the 1st (walking) trajectory. In the case, the hoof rolling over the ground and during the mid-stance mul- tiple contact points were examined at α t = −0.23° in which hoof was almost parallel with the ground. When the hoof leaves the ground, angular variation during the contact phase was reached 26°. For the 2nd (running) trajectory, the contact phenomenon occurred at the instant time t j = 0.208 s and α j = 6.20° firstly, shifted to leaving from the ground at t , = 0.478 s and α , = 71.08°. Finally, the hoof angle was reached 64.8°, which is approximately 2.5 times larger than the result in the walking gait. Interestingly, the hoof–ground interaction was significantly influenced by the initial impact phenomenon, which differentiate the grounding part either heel or toe and it reflects to successive sub-phases too. Even in the simple model, the differentiation was clearly observed not only in the trajectory level as a posture and kinematic analy- sis but also in the kinetic level as the dynamics analysis, which is easily extended to the energy analysis. All these sub-phases are characterized as representative sub-phase to reconstruct the target gait and it can provide the parametric analysis with stiffness and damping, which is associated with an actual parameter from soft tissues in the distal limb. It is because that the compliant contact force model was theoretically described as shown in the method section, which allow to change those parameters related to reac- tion forces easily. In the Matlab-based numerical simulation, ground reaction forces were evaluated by using Lankarani and Nikravesh [ 20 ] model by changing restitution coefficients as cr = 1 for elastic and cr = 0 for inelastic contact cases. According to the analysis, the max- imum contact force was obtained as 453.65 N for the 1st (walking) trajectory in the inelastic contact condition and 175.34 N in the elastic case. Average normal force was 80.83 and 34.88 N respec- tively. For the 2nd (running) trajectory, the maximum GRF was 8.694 kN for the inelastic contact condition and 5.719 kN in the elastic case. Average normal force was 978 and 800 N respectively. Result showed that the significant reduction of ground reaction force with respect to the viscoelastic contact was estimated as 61.35% in 1st trajectory and 34.22% in 2nd trajectory. In the walk- ing trajectory, multiple point support occurring in which hoof loaded the body weight at mid-stance as shown in Figure 5a , the maximum pressure or force peak indicated in only when the leg first impact and leaves the ground which was consistent with the force pattern observed in the human walking gait [ 26 , 27 ]. The GRF vector by arrow in Figure 6 represents magnitude, direc- tion of resultant ground reaction force through the entire stance phase in temporal evaluation of Figure 3 for walking and running trajectory. It simply shows that braking force in walking trajectory because in during stance phase ground reaction force was directed to the against moving direction. On the other hand, slight rotation of force vector respect to the gait cycle in stance phase in running gait which indicates the propulsive directional force acting toward center of pressure of moving leg. However, it should be noted that in our simulation leg mechanism was fixed in sagittal plane, there- fore detail analysis of force direction could be estimated precisely when multiple leg considered to regenerate the different locomo- tive gaits [ 28 ].
4.  CONCLUSION AND DISCUSSION
For the sake of the establishment of the generalized analysis of the viscoelastic effect with respect to various gait patterns, we intro- duced a novel theoretical method as the integrative framework 82 of the MBD and contact force model and then the computational analysis clearly demonstrated results in four different conditions from the combination of walking/running and elastic/inelastic contact. In the simplification of the hoof mechanism by using a triangle structure to focus on three nodes as heel, ankle and toe, the result revealed not only the two peak generation in the temporal evolution of the ground reaction force during the stance phase of the walking gait, but also a single peak observation in the case of the running gait. The phenomenon was generated from that wheel- like rolling motion of the hoof support and stabilizes the body in mid-stance, and partially consistent with the observation in the human foot–ground interaction [ 1,10 ]. In spite of the existence of the complexity in the horse leg structure, our computational result suggested a possible way to analyze the spring-damper effect embedded in principle to the biological mechanism for absorb- ing the ground reaction force flexibly as a common tendency of the walking and running gaits. The analysis with damping factor model by Lankarani and Nikravesh [ 20 ] realized the fact that a soft grounding effectively reduces the ground reaction force. In the actual phenomenon due to the interaction in the musculoskeletal systems such as a vibration from the sensory feedback loop [ 29 , 30 ], the further detail modeling is crucial. Phenomenologically in the level of the force generation, we successfully observed a time delay in a peak in the force time profile specifically in running gaits, which indicates an intrinsic behavior due to the damping factor is highly important in the kicking motion of the hoof. The phe- nomenon is also need to be verified with more detail leg models and evidences from the biological system [ 31 , 32 ]. Further consid- erations are possible to compare in the energy consumption due to the elasticity [ 18 ] and robotic application [ 33 ]. Geometric char- acteristics of the hoof is of interest to researchers studying on the biological nature mechanical in the hoof shape and its influence in joint angles and contact timing [ 9,34 ].
CONFLICTS OF INTEREST
The authors declare they have no conflicts of interest.
ACKNOWLEDGMENTS
This work was supported in part by Mongolia-Japan higher Engineering Education Development (MJEED-JICA) through the joint research project "Development of industrial, service and intelligent systems based on advanced technologies" (J23A16), JSPS KAKENHI (16H01616, 17H06383) and the New Energy and Industrial Technology Development Organization (NEDO), and Project on Regional Revitalization Through Advanced Robotics (Kyushu Institute of Technology/Kitakyushu city, Japan). [1] D.J. Mayich, A. Novak, D. Vena, T.R. Daniels, J.W. Brodsky, Gait analysis in orthopedic foot and ankle surgery—topical review, part 1: principles and uses of gait analysis, Foot Ankle Int. 35 (2014), 80–90. [2] L.J. Dockery, Analysis of gait propulsion during the ankle and forefoot rocker sub-phases, PhD Thesis, University of Central Arkansas, Conway, AR, 2019. [3] T. McGeer, Passive dynamic walking, Int. J. Robot. Res. 9 (1990), 62–82. [4] P. Mahmoodi, R.S. Ransing, M.I. Friswell, Modelling the effect of 'heel to toe' roll-over contact on the walking dynamics of passive biped robots, Appl. Math. Modell. 37 (2013), 7352–7373. [5] D. Torricelli, J. Gonzalez, M. Weckx, R. Jiménez-Fabián, B. Vanderborght, M. Sartori, et al., Human-like compliant loco- motion: state of the art of robotic implementations, Bioinspir. Biomim. 11 (2016), 051002. [6] E. Corral, M.J. Gómez García, C. Castejon, J. Meneses, R. Gismeros, Dynamic modeling of the dissipative contact and friction forces of a passive biped-walking robot, Appl. Sci. 10 (2020), 2342. [7] A.H. Hansen, D.S. Childress, E.H. Knox, Roll-over shapes of human locomotor systems: effects of walking speed, Clin. Biomech. (Bristol, Avon) 19 (2004), 407–414. [8] A.H. Hansen, D.S. Childress, Investigations of roll-over shape: implications for design, alignment, and evaluation of ankle-foot prostheses and orthoses, Disabil. Rehabil. 32 (2010), 2201–2209. [9] A. Karatsidis, G. Bellusci, H. Martin Schepers, M. de Zee, M.S. Andersen, P.H. Veltink, Estimation of ground reaction forces and moments during gait using only inertial motion capture, Sensors (Basel) 17 (2017), 75. [10] A. Bonnefoy, S. Armand, Normal Gait. Orthopedic Management of Children with Cerebral Palsy: A Comprehensive Approach, Nova Science Publishers Inc., Hauppauge, NY, 2015, pp. 199–213. [11] A. Leardini, J.J. O'Connor, S. Giannini, Biomechanics of the nat- ural, arthritic, and replaced human ankle joint, J. Foot Ankle Res. 7 (2014), 8. [12] W. Lee and M. Raibert, Control of hoof rolling in an articulated leg, Proceedings of the 1991 IEEE International Conference on Robotics and Automation, IEEE, Sacramento, CA, USA, 1991, pp. 1386–1391. [13] D. Batbaatar, H. Wagatsuma, A proposal of the kinematic model of the horse leg musculoskeletal system by using closed link- ages, Proceedings of the 2019 IEEE International Conference on Robotics and Biomimetics (ROBIO), IEEE, Dali, China, 2019, pp. 869–874. T. Jansen, The Great Pretender: [works of art by], nai010 [14] Publishers, Rotterdam, 2007. [15] K. Komoda, H. Wagatsuma, Energy-efficacy comparisons and multibody dynamics analyses of legged robots with differ- ent closed-loop mechanisms, Multibody Syst. Dyn. 40 (2017), 123–153. [16] S.D. Starke, J.J. Robilliard, R. Weller, A.M. Wilson, T. Pfau, Walk– run classification of symmetrical gaits in the horse: a multidi- mensional approach, J. R. Soc. Interface 6 (2009), 335–342. [17] J.J. Thomason, M.L. Peterson, Biomechanical and mechanical investigations of the hoof-track interface in racing horses, Vet. Clin. Equine Pract. 24 (2008), 53–77. R. McNeill Alexander, Elastic energy stores in running verte- [18] brates, Am. Zool. 24 (1984), 85–94. [19] J. Baumgarte, Stabilization of constraints and integrals of motion in dynamical systems, Comput. Methods Appl. Mech. Eng. 1 (1972), 1–16. [20] H.M. Lankarani, P.E. Nikravesh, A contact force model with hysteresis damping for impact analysis of multibody systems, J. Mech. Des. 112 (1990), 369–376. [21] J.A.C. Ambrósio, Impact of rigid and flexible multibody systems: deformation description and contact models, in: W. Schiehlen, M. Valášek (Eds.), Virtual Nonlinear Multibody Systems. NATO ASI Series (Series II: Mathematics, Physics and Chemistry), vol 103, Springer, Dordrecht, 2003, pp. 57–81. [22] V.P. Chervinskii, Dissipative function of coulomb friction on elliptical contact area, J. Friction Wear 32 (2011), 386–393. [23] P.E. Nikravesh, Planar Multibody Dynamics: Formulation, Programming with MATLAB*, and Applications, 2nd ed., CRC Press, Boca Raton, 2018, p. 430. [24] G.H. Ristow, Simulating granular flow with molecular dynamics, J. Phys. I France 2 (1992), 649–662. [25] D. Batbaatar, H. Wagatsuma, An extension of kinematic model for the linkage mechanism to analyze the workspace and gait tra- jectory induced by the swing leg motion, Proceedings of the 2020 3rd International Conference on Control and Robots (ICCR), IEEE, Tokyo, Japan, 2020, pp. 57–64. [26] R. McNeill Alexander, Energetics and optimization of human walking and running: the 2000 Raymond Pearl memorial lecture, Am. J. Human Biol. 14 (2002), 641–648. He received his B. Eng. degree in mechanical engineering from Mongolian University of Science and Technology (MUST), Mongolia, in 2013, the M. Eng. degree in the field of mechatronics from MUST. Currently, he is a PhD degree can- didate in Graduate School of Life Science and Systems Engineering in Kyushu Institute of Technology, Japan. His research interests include mechanism design, computational non-linear dynamics, and mechatronics. [27] R.M. Alexander, A.S. Jayes, Fourier analysis of forces exerted in walking and running, J. Biomech. 13 (1980), 383–390. [28] H.M. Clayton, S.J. Hobbs, Ground reaction forces: the Sine Qua Non of legged locomotion, J. Equine Vet. Sci. 76 (2019), 25–35. [29] D.F.B. Haeufle, K. Stollenmaier, I. Heinrich, S. Schmitt, K. Ghazi- Zahedi, Morphological computation increases from lower- to higher-level of biological motor control hierarchy, Front. Robot. AI 7 (2020), 511265. [30] Z.H. Shen, J.E. Seipel, A fundamental mechanism of legged loco- motion with hip torque and leg damping, Bioinspir. Biomim. 7 (2012), 046010. [31] I. Abraham, Z. Shen, J. Seipel, A nonlinear leg damping model for the prediction of running forces and stability, J. Comput. Nonlinear Dynam. 10 (2015), 051008. [32] A.M. Wilson, M.P. McGuigan, A. Su, A.J. van den Bogert, Horses damp the spring in their step, Nature 414 (2001), 895–899. [33] E. Garcia, J.C. Arevalo, G. Muñoz, P. Gonzalez-de-Santos, On the biomimetic design of agile-robot legs, Sensors (Basel) 11 (2011), 11305–11334. [34] C. Johnston, W. Back, Hoof ground interaction: when biome- chanical stimuli challenge the tissues of the distal limb, Equine Vet. J. 38 (2006), 634–641. He received a M.S. in Mathematical Sciences in 1997 and a Dr. Sci. in Mathematical Sciences in 2005 from Tokyo Denki University. In 2000, he became a Special Postdoctoral Researcher at RIKEN for studying computational models focus- ing on the brain oscillation. From 2003 to 2008, he was a Research Scientist in the Laboratory for Dynamics of Emergent Intelligence at the RIKEN Brain Science Institute. His research background started from theoretical modeling of brain oscilla- tions, the memory integration process of experienced episodes, and the implementation of oscillatory neural networks into neurorobotics. He is currently an Associate Professor of Kyushu Institute of Technology (KYUTECH) in the department of Brain Science and Engineering, Graduate School of Life Science and Systems Engineering. His research interests have been extended to include Bio-medical Signal Processing and Sparse Coding, Sport Dynamics and Synergy analysis, Brain-Inspired Robotics and Neuroinformatics Institute of Technology, Japan. 84