A new method to estimate effective elastic torsional compliance of single-stage Cycloidal drives

Mechanism and Machine Theory 105 (2016) 185–198

Contents lists available at ScienceDirect

Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmt

A new method to estimate effective elastic torsional compliance of single-stage Cycloidal drives Naren Kumar ⁎, Vladis Kosse 1, Adekunle Oloyede 2 Science and Engineering Faculty, School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, 2, George Street, GPO Box 2434, Brisbane, Queensland 4000, Australia

a r t i c l e

i n f o

Article history: Received 22 October 2015 Received in revised form 18 May 2016 Accepted 25 June 2016 Available online xxxx Keywords: Cycloidal drives Torsional stiffness Speed reducer Mechanical power transmission KHV Planetary gears

a b s t r a c t Standard single-stage tooth-to-pin-contact Cycloidal drives deliver various benefits such as high-ratio speed reduction, highly efficient torque delivery, compact physical structure etc. They lack torsional rigidity on their own because of inherent “lost-motion” and “backlash”, and hence inappropriate for precision-motion mechanical systems. However from dynamics point of view, this is beneficial for non-precision-motion systems, as it reduces the drive-train's shock factor. Currently there are no design standards for Cycloidal drives owing to their complicated component stiffness behaviour, increased tooth-load-sharing at overloads etc., which further obscures the analytical estimation of torsional rigidity of a given configuration. This paper presents a novel method comprising both analytical and numerical techniques for the effective determination of the elastic torsional compliance of single-stage Cycloidal drives based on static experimental results conducted on a commercially available gear-drive. We establish a unique key parameter − ηOP, ‘torque transfer efficiency’, of output-shaft-pins from mechanism-kinematics to be included in the system's dynamic model. Applying the techniques and outcomes presented here in a lumped mass/inertia dynamic model yielded agreeable natural frequency of torsional oscillations in comparison to experimental results obtained under the same loading conditions. This can lead to dynamically optimised designs and hence their standardisation. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Tooth-to-pin contact single-stage Cycloidal drives (invented in the 1930s [1]), are currently becoming increasingly popular for mechanical power transmission – because of their high-efficiency, high-ratio torque delivery in a compact physical form and high shock endurance factor. They are mainly classified as “KHV drives”, which are derived from the family of planetary gearing systems [2]. The working principle of Cycloidal drives is well described in most available literature [1–7]. Several design variants of the Cycloidal drive are available at present intended for various applications [8–20]. However, in today's industrial practice, two types of Cycloidal drives are commonly used – those without precision output shaft rotation, and those with precision output; as presented schematically in Fig. 1(a,b) respectively. Single-stage low-precision output Cycloidal drives (classified as KHV type transmissions), are best suited for general applications where accurate output shaft rotation is not critical, e.g.: conveyor belt drives, drum agitators etc., (see Fig. 1a). Contrarily, for transmissions used in the field of robotics and

⁎ Corresponding author. E-mail addresses: [email protected], [email protected] (N. Kumar), [email protected] (V. Kosse), [email protected], [email protected] (A. Oloyede). 1 This Author's permanent address is 117, Jones Road, Carina Heights, Queensland, Australia 4152. 2 This author's current address is Elizade University, P.M.B. 002, Ilara-Mokin, Ondo State, Nigeria, Africa.

http://dx.doi.org/10.1016/j.mechmachtheory.2016.06.023 0094-114X/© 2016 Elsevier Ltd. All rights reserved.

186

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

Fig. 1. Schematic construction of Cycloidal drives; (a) Single-stage KHV type Cycloidal drive; (b) Two-stage 2KV type Cycloidal drive with primary planetary reduction.

servo systems where precision output motion and repeatability are vital, a two-stage high precision-motion Cycloidal drive, known as the 2KV type [21] (see Fig. 1b), is preferred. This two-stage gear transmission has an initial low-ratio planetary gear reducer prior to the cycloidal stage. The planetary stage greatly contributes to the overall torsional rigidity of the drive-train, which is not so pronounced in the case of single-stage KHV variants. The large input torsional compliance in KHV Cycloidal drives acts as a buffer to shock loads (torsional shock absorber), which are distributed among multiple load-bearing surfaces of drive-train components almost instantaneously [22]. Hence they are most suitable for applications with moderate to high operating shocks and frequent rotation direction reversals, as they excel in performance compared to conventional equivalents. Despite the increasing popularity, there are currently no available design standards for Cycloidal drives classified under mechanical transmission methods [23]. One of the prime literature start points of Cycloidal drive research is the original patent of probably the first tooth-to-pin contact single-stage Cycloidal drive namely “Gear Transmission” (1928 US Patent No. 1,694,031), by Lorenz K. Braren [1]. He has thoroughly described the construction of his patented design with geometric and kinematic groundwork. However very little information is available on the effect of drive-train dynamics or component design methods to endure the effects of dynamic forces. Other early patents of similar variants also provide valuable information although the list is not exhaustive [1,24–36]. Among the non-patent literature resources available on Cycloidal drives – Botsiber and Kingston [5], have provided the construction method of the Cycloidal disc (gear) starting from basic parameters and speed ratio. They have highlighted that the operating as well as the predominant surface stresses in the system are not the same as in conventional (spur/helical) counterparts. They state that the loading conditions are analogous to those of rolling element bearings or “free-wheeling clutches”, since the prime components of the drive-train are also made of high strength bearing-grade steel. Malhotra and Parameswaran [6], have presented the procedure to calculate the forces acting on the drive-train elements and in addition, have provided methods to derive the theoretical efficiency, determine the thickness of the Cycloidal disc, described the effects of design parameters on contact stresses and the force system in the drive-train. Yang and Blanche [37] have presented the drive selection method (similar to selection of rolling element bearings) for particular applications based on transmission ratio and operating conditions. They also have presented two empirical formulae to determine backlash and torque ripple. Litvin, along with Fuentes [38,39] have clearly explained methods to write rigid body equations of loci of points of interest on conventional as well as Cycloidal tooth surfaces using coordinate transformation. Later researchers have used this technique to develop the tooth profile as well as generation of tooth meshing equations for Cycloidal drives. Hwang and Hsieh have applied Litvin's approach to gerotor type hydraulic pump mechanisms besides Cycloidal drives, to develop mathematical models of the teeth-meshing process that aids in its design [40]. Focusing on design of gear geometry, others [41,42] have developed formulae to avoid undercutting in Cycloidal drives. Cycloidal drive mechanisms can have single or multiple ‘teeth differences’ (between the numbers of teeth on the ring gear and the disc respectively). Most researchers including Meng et al. [43] have presented analysis results of mathematical models of the cycloidal mechanism with one tooth difference, while several others [44–48] have analysed the teeth-mesh-forces using models and computer simulations as this is a computation-intensive process. Sensinger and other researchers [49–51] have provided unified equations for calculating stress, efficiency and optimisation of the mechanism synthesis, besides comparing them with high ratio single-stage harmonic drives. So far, information about the general design and construction of both KHV and 2KV drives are available in several literature sources [1,5,32,52,53]. However, only a few researchers have vaguely shed light on the dynamics and stiffness characteristics of the drive-train components [54–58]. Time-dependent dynamic forces such as shocks or from the high starting torque of the motor/driver can severely load the components of the drive-train of the transmission system. This necessitates a design shock factor to be incorporated into the drive-train system design. In tooth-to-pin contact single-stage Cycloidal drives large input torsional compliance is inherent by design. Then the method to include a shock factor for design safety remains as a question. Experimentally, dynamic analyses of drive-trains are usually performed using special purpose-built test rigs that support both load and speed transition inputs while the subjected drive-train is transmitting power [7,59]. Alternatively a predictive rigid multi-body model with discreet flexible element-joints can be used to analytically determine the dynamic responses to sudden load transitions. This can be applied to any generic type of drive-train as is currently done for wind-turbine planetary spur/helical transmissions [60,61]. Such models are currently unavailable for high-transmission ratio speed reducers especially – Cycloidal

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

187

drives. Owing to the complexity of their kinematic linkage, all parameters necessary to build a dynamic model are quite challenging to obtain solely from analytical methods. Large torsional compliance and backlash at the input end result in dead zones, which further complicate the dynamic model synthesis. Computational methods such as Finite Element Analyses (FEA) can also be used along with analytical techniques to determine the stiffness and inertia properties of all element-joints of the drive-train. This paper attempts to fill this gap by presenting a new modelling technique using both analytical and computational approaches to determine the effective torsional compliance of single-stage KHV type Cycloidal drives validated by experimental research conducted on a commercially available reducer. The principle of this technique can be applied to generic high-ratio Cycloidal drives (any KHV variant), to determine the torsional compliance and critical joint-torque transfer efficiency necessary to synthesize a multi-body dynamic model of the drive-train. 2. Analysis of a typical single stage Cycloidal drive Mechanically, ‘torsional stiffness’ of a transmission system is the ratio of the applied torque to the resultant elastic torsional deformation of the loaded drive-train. Mathematically torsional compliance is the reciprocal of torsional stiffness. The total lost motion of a Cycloidal drive system is made up of individual torsional compliances of all drive-train components. Practically, manufacturing and assembly tolerances inevitably add up to the theoretical torsional compliance, resulting in the total actual lost motion. In this paper, as part of an investigation of Cycloidal drive-train dynamics, techniques to develop a component-level predictive dynamic model of a typical Cycloidal drive-train using analytical and computational methods (FEA) are presented, which are validated by system level experimentation performed on an actual transmission system. 2.1. System-level experimental analysis At the Queensland University of Technology, Brisbane, Australia (QUT), pilot system-level studies were conducted by Kosse V., [62] through static torsional stiffness study of an output-shaft-locked single-stage Cycloidal drive with speed ratio 87:1 (model CHH413087, manufactured by Sumitomo Heavy Industries Ltd., Japan [3]). The same experimental setup was used to assess the linear part of the torsional stiffness curve and its related hysteresis in this research. The output (slow-speed) shaft of the drive was locked and secured to the housing using a cap and high-strength bolts. The input end showed large free-play comprising both lost motion and backlash totalling to about 56° for this particular model; i.e., about 28° to either side (clockwise or counter-clockwise) from a middle datum. The setup is as shown in Fig. 2 (photograph on the left). The Cycloidal drive unit was clamped to the workbench. A large pulley carrying a 1/2° protractor dial on its rim is attached at the input (high-speed) shaft end of the drive system. A steel rope is passed around the pulley-system, whose ends are connected to an aluminium crossbar and a weight hanger. Using this arrangement, a static force couple can be applied to the input shaft by placing a weight on the hanger. Since the output shaft of the drive is fixed in this case, with no rotation, the drive train components move using the available relative degrees of freedom within the system. Cycloidal tooth profiles are generally wider than the pins they mesh with, and also clearances are present between mating parts to allow for kinematic errors (if any), assembly tolerances, thermal expansion etc. – resulting in a large angular free play at the input shaft. The pulley could be turned freely by hand until it reaches a dead-stop, in either direction of rotation. Additional weights added to the hanger below further increase the input-end load, forcing the drive-train components to deform elastically and redistribute the forces according to their acquired new positions. The angle of twist was read directly from the protractor. This speed reducer has a manufacturer's rated output torque of 585 Nm for an input power of 1.14 kW at 1500 RPM [3], from which the nominal input torque TN was calculated as the ratio of input power to the angular velocity – which evaluated to be 7.2575 Nm in this case. Torque was

Fig. 2. Static loading of output shaft-locked single-stage Cycloidal drive; (left) Experiment setup; (right) loading and unloading responses of the drive-train up to 4.3 TN (TN = 7.2575 Nm - nominal input torque).

188

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

calculated for each load step of dead weight applied on the weight hanger, while the corresponding angle of elastic deformation was read directly from the protractor dial. All torque values are plotted as multiples of the nominal input torque TN on the ordinate corresponding to their angular deformation θexp values (in degrees) on the abscissa of the plot shown in Fig. 2. Only the steady clockwise loading and unloading responses are shown in this plot as angle of turn vs. multiples of rated input torque load. The loading curve begins from +28° (after the dead zone), of input shaft rotation from the middle datum and rises linearly up to about 2.5 TN resulting in an input shaft turn angle of 75°–80°, and then becomes non-linear for higher loads. The unloading path differs greatly from the loading, thus depicting the hysteresis of the system. Cycloidal drives exhibit “Hereditary” characteristics in their loading-unloading hysteresis loops: meaning – previous cycles of loading and unloading affect the current and future loading and unloading cycles. This is mainly because of clearance among the mating components and the large torsional compliance of the drive-system. Harmonic drives, which are also well-known for high ratio speed reduction in a compact single stage, also depict hereditary behaviour [63]. Literature or research on this aspect of all high ratio transmission systems in general, particularly Cycloidal drives, is unavailable at present. In hereditary systems the initially set datum point (the selected zero load point before applying the load to the system for the first time), is likely to drift away from their original position over successive loading and unloading cycles as is shown in Fig. 3. In Fig. 3, the plots of hysteresis loops of four successive trials of loading up to 1 TN before fully unloading the drive-train are shown. Trial runs comprising of trial names “Loading 1” and “Unloading 1” form a completed hysteresis loop. It is clearly noticeable that the initial hysteresis loop named “Loading 1” and “Unloading 1”, moved left to a new position “Loading 2” and “Unloading 2” for the second run, before shifting towards the right over the successive runs (3 & 4) again. Also, each run depicts slightly varying nonlinearity in both loading and unloading responses – indicating that the initially selected middle datum (+28°) was not constant throughout and; the current and successive loaded or unloaded positions are dependent on effects of previous loading cycles (hereditary hysteresis behaviour). To prevent this, before every new experiment run, the system had to be ‘reset’ by turning the input shaft in both clockwise and anticlockwise directions twisting slightly beyond both dead-centres to bring the datum back to the initially set zero point. This resets the system and re-aligns the internal components so that a trial response could be read starting from +28° mark. If the resetting procedure was not performed, the middle datum, which was initially at 0°, would move over several runs up to 1° either in the positive or the negative direction, for successive loading-unloading trials leading to a misconception of a progressively increasing dead zone. 2.2. Component-level FE analysis Response to the static load cycle shown in the plot of Fig. 2 can be better understood if the same loading cycle was applied in a component-level study such as a 3D Finite Element model simulation. Hence ANSYS Workbench FEA code (version 14.5.7 – ANSYS Inc. USA), was used on a 3D Finite Element model of a simplified drive-train model to get a clearer insight of component-interactions. The 3D CAD model used for FEA was prepared based on accurate measurements taken off components of the actual drive system. Since geometric and positional tolerance information of drive-train components are not available they have been ignored in this study. Individual output shaft pin deformations, which result due to the movement and spin of the Cycloidal discs as they are carried around the periphery of the fixed annulus because of the tooth mesh, would contribute to the overall torsional compliance of the

Fig. 3. Hysteresis loops of Cycloidal drives exhibit hereditary characteristics as the original datum (start point) drifts away due to successive loading cycles.

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

189

system. As the arrangement of the components is symmetric, only one Cycloidal disc and its associated components will be best suited for FEA scenarios. A simplified 3D CAD model of a single-disc Cycloidal drive-train (see Fig. 4), was used in this study, mainly to reduce the FEA problem size so that a practical solution could be achieved. In an actual system however, two or more opposed discs are usually used to dynamically balance each other's eccentric rotating masses. As is seen in Fig. 4, the simplified drive-train has some merged components and only the output-shaft-flange with bushes fused to their pins was considered, ignoring its remaining portion. Input torque was applied to the simplified eccentric cam while the output flange and the annulus (ring gear), were held fixed. The Cycloidal disc was positioned such that its top most tooth space was in full contact with the top most pin of the ring gear. Subsequent tooth spaces will be in partial contact with their conjugate pins. All contacting teeth form the tooth mesh zone of this disc. There are 2 Cycloidal discs and 8 output shaft pins in this actual Cycloidal drive configuration. Only 4 out of the 8 output shaft pins actively transmit torque per disc in this arrangement. The two discs are mutually transposed with a phase difference of 180° for dynamic balance. Hence when one disc is meshing with the Annulus pins at the top, the other will be interacting with those at the bottom. Kinematically each disc will be actively transmitting torque through only half the total number of output pins that are located on the opposite half of the disc relative to its tooth mesh zone at any given time. The other half will receive power from the second Cycloidal disc at the same time. Hence at the output end of this FEA model, contacts were considered only for four (of the total eight), output-shaft pins located at the left-half of the single Cycloidal disc. This is because only these four pins can actively transmit torque for the given tooth mesh zone position. Angles of twist determined from the hysteresis experiment (θexp from Fig. 2), were used as inputs in individual static structural FEA scenarios to obtain resultant moment reactions as the outputs. ‘Rough’ contacts were applied to one half of the ring-gear's fused pin surfaces and the disc's teeth surfaces they mesh with, in the direction of disc travel. ‘Rough contact’ in ANSYS Workbench means: normal separation of the mating parts is permitted, but not sliding [64]. Frictional contacts for machined surfaces (friction coefficient 0.05), were provided to output pin-bushes and the oversized-hole-surfaces of the disc they touch. A ‘fine swept mesh’ comprising hexahedral elements was prepared with 25% global ‘relevance’ for this model, and a ‘face sizing’ of 1 mm was applied to contact surfaces of oversized-holes in the disc with ‘mapped mesh faces’ to get uniform sized regular hexahedral mesh-geometry, complete with 958,773 nodes in 249,795 elements altogether. The calculated output torque from the hysteresis experiment was compared with that obtained from FEA as is shown in Fig. 5. The FEA results also showed a linear trend for simulated four small loads only, as this was a computation-intensive problem involving complex geometry. Higher ANSYS results deviated by about 15% in comparison to those from the hysteresis experiment. This would probably be due to the simplified drive-train geometry used in the FEA model where the entire output shaft was not considered. The difference could also be because of the pure torsional loads used in FEA, whereas in the experiment, the vertically hanging weight imparted a combined load among mating drive-train components. 3. Analysis FEA results showed bending deformations of output shaft pins, loaded like individual cantilevers. Fig. 6 shows exaggerated bending of pins for clarity. From FEA results, it was discovered that even for light loads, some of the output-pins bent radially outwards instead of bending in the tangential direction caused by the torque-producing tooth-mesh force component. This meant that at any given instant

Fig. 4. Simplified 3D CAD model of the Cycloidal drive-train.

190

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

Fig. 5. FEA results also showed a linear trend matching closely with those of the hysteresis experiment.

in the tooth-mesh cycle, some pins will be dissipating energy by bending radially outwards and back – hence not contributing to torquetransfer while other pins do. This means that torque transferred to the output shaft is unevenly delivered by its pins, proving that not only the radial and tangential force components are irregular over the tooth-mesh cycle, but also their magnitudes and directions depend on tooth-mesh positions at any given instant. This led to further investigation wherein a map of active pins was developed to estimate the number of output pins actively transferring torque at any instant over the tooth-mesh cycle per Cycloidal disc. 4. Map of active pins Fig. 7 shows a front view of the Cycloidal disc assembled within the annulus (only annulus pin-sockets are shown here for clarity), whose top most tooth is in full contact with the corresponding (top most) ‘Pin 0’ of the annulus. The Cycloidal disc is carried by the eccentric bearing (only the inner bearing race and one roller is shown here for clarity), at the central axis. The Cycloidal disc follows the eccentric bearing, as it rotates clockwise for example, since their races are mutually concentric. Due to the annulus-pin mesh, the disc is forced to spin about its centroidal axis at the speed ratio because of the difference between the numbers of annulus pins and disc teeth. In this particular drive-train example, 4 annulus-pins (ring gear pins) were kinematically in mesh with the disc teeth at any given time by observation; and one tooth surface was in full contact mesh with its conjugate pin. Let this geometric condition where the middle tooth-space is in full mesh and the neighbouring pins in partial meshing position be defined as the “tooth-mesh zone”. For convenience, let the topmost annulus-pin be designated Pin 0 and successively Pins 1 and 2 etc., are located clockwise relative to central axis (see Fig. 7). As the eccentric bearing race (input) rotates clockwise, the disc changes position and simultaneously begins to spin about its own axis due to the pin-tooth mesh, thus moving the tooth-mesh zone (which was initially under Pin 0), to be in full mesh successively

Fig. 6. Irregular bending of output shaft pins revealed in FEA.

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

191

Fig. 7. Map of active pins on the left half of the disc for a full clockwise rotation of the tooth mesh zone.

with Pins 1, 2, and so on. At these disc-orientations, the positions of the oversized-holes of the disc would also change accordingly. There are 8 equispaced output shaft pins on the flange in this drive-train design. Output pin-centres are located on a theoretical pin-pitch circle on the flange. The number of pins is chosen by the manufacturer to fulfil the design for strength to distribute the large output torque (585 Nm in this case), among multiple pins. The number of outputpins employed can vary in various Cycloidal/KHV drives of similar construction, depending on their sizes and power capacities. A 2D kinematic position analysis was performed using SolidWorks CAD software (version 2012, Dassault Systèmes, USA), wherein the loci of the oversized-hole centres on the left half of the disc were plotted for a full clockwise travel of the tooth mesh zone. In Fig. 7, let H1, H2, H3 and H4 be designated for holes on the left-half of the disc. Due to the compound motion, the centres of all oversized-holes trace hypotrochoids as is seen in Fig. 8(b). The part of the hypotrochoid locus which is outside the boundary of the pins' pitch circle on the flange, forces the output-pin to bend radially outwards. Since there is no tangential component of the reaction force, torque is not produced at this position of the locus. The centre of the oversized-hole in the disc has to be on the wider portion of the hypotrochoid locus (path shown in green colour in Fig. 8c), to transmit torque by pushing the pin in the tangential direction relative to pin-pitch circle – which is determined by the tooth mesh zone location. Since movement of the output shaft pins is governed by the over-sized hole-positions on the Cycloidal disc, each output pin has an active torque transmission zone over the full tooth-pin mesh cycle. In other words, each output shaft pin will actively transmit torque to the output shaft only when the tooth mesh zone of the Cycloidal disc is within a certain angular range somewhere around the central axis of the annulus. Thus it is understandable that torque outflow per pin is not continuous but pulsating, due to uneven load sharing among output pins. This necessitated a further study through the generation of a map of all torque transmitting zones of the holes for one full rotation of the tooth mesh zone to estimate the number of active pins at any given instant. Using 2D geometry analysis in SolidWorks software

Fig. 8. Locus of oversized-hole centres of the moving Cycloidal disc trace hypotrochoids, which cross the pins' pitch circle; (a) partial front view; (b) enlarged view of oversized-hole and the pin-bush at H1; (c) hypotrochoid locus of oversized-hole centre.

192

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

exact tooth mesh positions of the pin at H1 were determined, where torque transmission would start and cease over the mesh cycle – viz. the orange zone in Fig. 7. In the present case the tooth-mesh zone has to move almost 66° from Pin 0, before the pin at H1 begins to transmit torque to the output shaft. Similar maps were prepared for all other hole-positions (as is shown in Fig. 9 – coloured individual sectors marked H1-H4); when all plots are superimposed, a map of the active zones of all holes located on the left-half of the disc results (refer Fig. 7). Each coloured sector in Fig. 9 represents one individual output pin's active torque transmission zone when the tooth mesh zone moves clockwise from 0 to 180° along the annulus periphery. As the colours overlap, the numbers of pins within that position in the range add up – viz. transmit torque together in parallel. It can be seen that all four pins will be actively transmitting torque only for a short duration near the 90° position of tooth mesh zone. From 180° to 360° positions of the tooth mesh zone, the pins on the right half of the disc will become active in the same sense but inverted. Hence, a copy of the colour-superimposed map in Fig. 9 is inverted and superimposed again onto itself and hence the pin numbers are further added where colours overlap – to result in a full map showing the active pins per disc over a full 360° rotation of the tooth mesh zone (see third image from the left in Fig. 10). As there are two Cycloidal discs in this version, the number of pins in the map (third from the left in Fig. 10), are doubled to get the full map of active pins in the drive-system shown on the far right of Fig. 10. Here an uneven torque distribution is clearly evident (only 4 pins are active in some positions; see highlighted portion on the far right image of Fig. 10). An average of all pin numbers in the far right image of Fig. 10, is 6.75 – indicating that only 6.75 pins (statistically speaking), out of the total 8, are active at any given instant. Then the overall output torque transmitted by one pin can be expressed as:


TO ¼


0  TO NOP  N OP NOP

ð1Þ

TO⁎ is the average output torque due to active pins; TO is the rated output torque of the Cycloidal drive; NOP is the total number of output-pins present in the system and NOP′ is the number of active output-pins (pins which participate in torque transfer). The term within parenthesis in Eq. (1) is the “Output torque transfer efficiency” of a pin ηOP: ηOP ¼


0  NOP NOP

ð2Þ

In this example, ηOP ¼

6:75 ¼ 0:84375 ¼ 84:375% 8

ð3Þ

5. Discussion McIver in his 1950 patent [32] of a similar high-ratio speed reducer, states that the number of output pins to be chosen for any KHV design has to be a minimum of 9 or any greater feasible odd number. His application was to reduce speed from a turbine shaft running at very high speeds. For having nine or any higher odd number of output pins, he supports his claim with experimental evidence and shows that the output shaft's pulsating torque is kept to a minimum. From the results of experiments and FEA, the research outcome described in this paper has established that not all pins will be actively transmitting torque at any given tooth mesh position throughout its cycle; and has also derived the output pins' torque transfer efficiency ηOP. With NOP being 9 or any higher odd number, the active zones will overlap more closely (as is shown in ′ and NOP, thereby proving that output torque of a Cycloidal drive of this Figs. 9 and 10), thus reducing the difference between NOP type is not constant but pulsating – emphasizing the significance of drive-train dynamics.

Fig. 9. Maps of holes H1, H2, H3 and H4 overlapped show the number of active pins at various positions.

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

193

Fig. 10. Map (far right) showing active number of pins at all tooth mesh zone positions for two Cycloidal discs.

Torsional stiffness of the Cycloidal drive-train plays an important role in system dynamics, which can be used to predict the transient response of the transmission system in a lumped inertia model of the configuration. The output shaft pins are of cantilever-type in this design, whose bending deflections explains the linear trend in loading curves of the experiment as well as the FEA scenario mentioned before – highlighting the vital role of the output shaft in Cycloidal drive design and functionality. To validate ηOP, a theoretical stiffness model of the full output shaft alone was considered. Fig. 11 shows the profile view of the output shaft with many steps. Analytical torsional compliance formulae for commonly used mechanical components are given in several Russian literatures [65–70]. Torsional compliance formulae (in radians per newton-meter), for keyway, shaft-steps and the frustum of cone attaching shaft to the flange shown in Fig. 11, are given by [67,68]. The stepped output shaft acts as an assembly of linear torsional springs connected in series from the free end up to the flange. Onto the right of the flange are the output-pins, which transfer torque by bending elastically. Hence their bending deflections have to be converted to torsional deflection appropriately. The reciprocal of effective torsional stiffness of multiple springs connected in series is equal to the sum reciprocals of individual spring stiffness. But the effective torsional compliance of series connected springs simply add up, because compliance is the reciprocal of stiffness [71]. Effective torsional compliance of the above shaft comprising several steps up to the flange from the free end, is given by summing individual torsional compliances of steps using the analytical formula given by [67,68]: t

C OS ¼ ckeyway þ c1;2 þ c2;3 þ c3;4 þ c4;5 þ c5;6

ð4Þ

t is the torsional compliance of the output shaft up to the flange from the free end (i.e. starting from the far left end of the shaft in where COS Fig. 11 and moving towards the flange on the right end); ckeyway is the torsional compliance of the keyway; c1,2 , c2,3 , c3,4 , c4,5 are the respective torsional compliances of shaft-steps and c5,6 is that of the frustum of cone joining all the steps to the flange. ckeyway given by [67,68] is:

ckeyway

0 1 10hl 32 Bl þ D C ¼ @ A πG D4

ð5Þ

Fig. 11. Profile view of the output shaft of a typical Cycloidal drive-train.

194

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

where G is the torsional rigidity of the shaft material; h is the height and l is the length of the keyway respectively; D is the shaft diameter at the keyway. Step : c1,2, c2,3, c3,4, and c4,5 are given by a general step equation cOS Step

cOS ¼

" 32 l1 l 1 þ 2 þ πG d4 D4 8d3

3

1−

d D3

!# ð6Þ

where l1, l2, d, and D are lengths and diameters of smaller and larger parts of individual steps respectively according to [67,68] and:

c5;6

" !# 2 32 l d d d ¼ 1þ þ 2 πG d4 3D D D

ð7Þ

where d and D represent the smaller and larger diameters of the frustum respectively, separated by distance l. The pins on the flange deform elastically by bending while transmitting torque. Hence their bending equivalent torsional compliance has to be calculated. In Fig. 12, δ (in μm), is the elastic pin bending deflection; φ is the equivalent angular deflection (in radians), of one pin situated on the pin-pitch circle of radius ROP caused by the tangential force component Fn (in newtons); which is given by: −1

φ ¼ tan




δ ROP

 ð8Þ

and Fn ¼

T
O ROP

ð9Þ

from Geometry and Statics respectively. Here in Eq. (9), TO⁎ involves ηOP – see Eqs. (1) and (2). Also, in this particular design, the span to depth ratio of pins was less than 10 classifying them as short beams for which the normal beam bending theory is invalid [72], and Timoshenko beam theory has to be applied instead [73]. The elastic beam bending deflection δ was calculated in FEA using loads caused by multiples of applied input torque: 1 TN, 2 TN, 3 TN, 4 TN and 4.5 TN (see Fig. 13). The bending equivalent elastic torsional compliance of the shaft contributed by one pin is therefore given by: C OP ¼

φ ð F n  ROP Þ

ð10Þ

Fig. 12. Bending equivalent angular deformation is used to calculate the effective torsional compliance.

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

195

Fig. 13. Linear bending stiffness of the short-beam type pin was determined by static loading in FEA.

Then the bending equivalent torsional compliance of all pins together with the torsional compliance of the free end of the shaft (to the left of the flange in Fig. 11), is given by: equiv:

C OS

t

¼ C OS þ

C OP NOP

ð11Þ

Since the pins bend separately, they act as springs connected in parallel, hence COP has to be divided by NOP to get the equivalent compliance due to all pins. The total equivalent torsional compliance of the output shaft w.r.t the input shaft of the drive-system is usually preferred in system dynamics, because it is based on one reference member from a kinetic energy point of view [74]. Hence the equivalent reflected torsional compliance of the output shaft (w.r.t. the input shaft), is given by multiplying Eq. (11) by the square of the overall drive-train reduction ratio i, (which in this case is 87):


C OS ¼


   C t 2 rad C OS þ OP i NOP Nm

ð12Þ

The reciprocal of COS⁎ is KOS⁎, the reflected linear elastic torsional stiffness of the output shaft in newton-meters per radian. In this research, the theoretically determined COS⁎ was used to verify the actual large input shaft torsional compliance (described in the hysteresis experiment of Section 2.1), by assuming that only the output shaft was torsionally flexible and the rest of the drive-train components

Fig. 14. Comparing experimental input shaft turn angle θexp with the developed theoretical angle θ′considering only the output shaft to be flexible.

196

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

being rigid. Angle θ′ (in degrees), the input shaft would turn due to elastic torsional compliance of output shaft alone, was calculated as: 0

θ ½ ° ¼ ηOP


 ðC OS  nT N Þ  180 ° þγ° π

ð13Þ

where n = (1, 2, 3, 4, 4.5) – viz. the multiples of nominal input torque; γ° is the angle of input shaft free play in degrees – which in this case is 28° (see Section 2.1). θ′ was compared with θexp values obtained from the hysteresis experiment (in Fig. 2), by plotting both against nTN as is shown in Fig. 14. Theoretical estimate of the input shaft turn angle due to the torsional flexibility of the output shaft alone, agreed with the actual experimental results with a difference of about 8.5% at 4 TN as is seen in Fig. 14. This confirms that the output shaft's flexibility plays a major role in the drive-train dynamics and hence the design of the Cycloidal transmission system. Although this is a static test-based evaluation, dynamic loads such as shocks experienced in direction reversals (under load), or sudden torque fluctuations due to a prime-mover or a driven mechanism/machine, can be a higher multiple of TN for a short instant of time causing severe component deformations due to peak forces among the components. The techniques presented here can be applied to any type of high-ratio speed reducers with linear stiffness behaviour, such as planetary or cycloidal drive systems (KHV, 2KV or equivalent), to determine applicable equivalent torsional stiffness of the drive-system. Torsional stiffness of components of any drive-train is vital for dynamic analyses using lumped mass/inertia mathematical models, which can be used to characterise the system's transient response to time-dependent loads. 6. Conclusions In this research, effective torsional compliance of a single stage Cycloidal drive-train system was investigated starting from a system-level study using a commercially available actual model. The model was subjected to torsional loads at its input shaft while holding the output shaft immobile. The linearity of torque to twist angle relationship was confirmed even in a component-level FEA study. Upon observing the bending deformation of the output shaft pins in FEA results, it was discovered that not all slow-speed (output) pins bend in the same direction at any given time. Some output-pins transfer torque by bending tangentially to their pitch circle, while others bend radially depending on their location relative to the tooth-mesh zone. Maps were hence prepared using overlapping coloured areas, to determine the number of active output-pins at any given instant. This developed to the introduction of a new term, “torque transmission efficiency” of the output shaft pins ηOP, which is the ratio of the number of active output-pins to the total number of output-pins in the transmission. For the particular Cycloidal drive design used in this research, the value of ηOP evaluated to be 84.375%. From first principles, the reflected equivalent torsional compliance (and stiffness) of the output shaft were developed using ηOP; from which theoretical input shaft turn angles were calculated and compared with the actual (hysteresis experiment) agreeing with each other to be within a 8.5% margin up to four times the nominal input torque (4 TN). This highlights the role of the output shaft in the design and drive-train dynamics of Cycloidal drives while validating the theory developed. The theoretical torsional compliance method presented in this paper (ignoring geometric and positional tolerances of components), can be directly applied to planetary or KHV drives in general, particularly for single-stage Cycloidal drives, provided – the drive system has linear stiffness behaviour and all parameters are known. In case of simple geometries, analytical torsional stiffness (or compliance), formulae can be available from various literature sources for their prescribed static load range. For unknown parameters (such as short-beam bending stiffness of output shaft pins in this case), FEA or experimental or other computational techniques can be used as aids to determine the theoretical values where required. As is always advisable, in general, theoretical formulations have to be practically verified through experimentation before standardising the application. Using these techniques, effective torsional stiffness of a mechanical transmission can be determined which is required for lumped mass/inertia model analyses to predict transient responses of drive-trains under time-dependent (shock) loads accurately enough. Transient response analyses of mechanical drive-trains help optimise their system dynamic (shock) factor, thus redefining the design philosophy of single-stage high-ratio speed reduction drive-trains from a dynamics point of view. Acknowledgements The authors wish to thank all personnel at Queensland University of Technology, Brisbane, Australia, who provided support for this research. References [1] Braren, L.K., (1928), Gear Transmission, U.S. Patent, # 1,694,031, [2] D. Yu, KHV planetary gearing, Gear Technology, Gear Technology, USA Nov/Dec 1987, pp. 21–31 (48). [3] SM-CYCLO, SM-CYCLO Speed Reducers and Gearmotors Catalogue [990341, 12/97], Sumitomo Cyclo Europe, Sumitomo Heavy Industries Cyclo Drive (Europe) GmbH, 1997. [4] Shimpo Circulute 3000 Catalogue, Shimpo Drives, Inc., USA, 2007. [5] D.W. Botsiber, L. Kingston, Design and performance of cycloid speed reducer, Mach. Des. 28 (13) (1956) 65–69. [6] S.K. Malhotra, M.A. Parameswaran, Analysis of a cycloid speed reducer, Mech. Mach. Theory 18 (6) (1983) 491–499.

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

197

[7] N. Kumar, V. Kosse, A. Oloyede, A novel testing system for a cycloidal drive, 12th International Conference on Intelligent Systems Design and Applications (ISDA), Cochin University of Science And Technology, Kochi, India, 2012 (Institute of Electrical and Electronics Engineers (IEEE)). [8] I. Hayashi, N. Iwatsuki, M. Kawai, J. Shibata, T. Kitagawa, Development of a piezoelectric cycloid motor, Mechatronics 2 (5) (1992) 433–444. [9] I. Gu, Design of antibacklash pin-gearing, KSME Int. J. 11 (6) (1997) 611–619. [10] T. Guan, L. Lei, Y. Sun, Y. Ma, L. Zhang, Inverse design used to new FA series cycloid drive, China Mech. Eng. 14 (Copyright 2003, IEE) (2003) 68–71. [11] L. Xin, H. Weidong, L. Lixing, L.C. Schmidt, A new cycloid drive with high-load capacity and high efficiency, Trans. ASME J. Mech. Des. 126 (Copyright 2004, IEE) (2004) 683–686. [12] T.-S. Lai, Geometric design of roller drives with cylindrical meshing elements, Mech. Mach. Theory 40 (1) (2005) 55–67. [13] T.-S. Lai, W.-H. Hsieh, G.-T. Chen, Geometrical design of roller drives with two-tooth difference, J. Chin. Soc. Mech. Eng. 28 (6) (2007) 641–648. [14] C. Gorla, P. Davoli, F. Rosa, C. Longoni, F. Chiozzi, A. Samarani, Theoretical and experimental analysis of a cycloidal speed reducer, J. Mech. Des. Trans. ASME 130 (Compendex) (2008) 1126041–1126048. [15] Z. Caichao, L. Mingyong, D. Xuesong, X. Ning, Z. Bin, Analysis on transmission characteristics of new axis-fixed cycloid gear, Adv. Mater. Res. 97-101 (Copyright 2011, The Institution of Engineering and Technology) (2010) 60–63. [16] W.-K. Nam, S.H. Oh, A Design of Speed Reducer with Trapezoidal Tooth Profile for Robot Manipulator, Journal of Mechanical Science and Technology, Vol. 25(1 (2011)), Springer, 2010 171–176. [17] C.C. Zhu, Z.X. Luo, M.Y. Liu, N. Xiao, Analysis on new gyratory fixed-spindle cycloid drive, Chongqing Daxue Xuebao/J. Chongqing Univ. 33 (3) (2010) 7–12. [18] M. Blagojevic, N. Marjanovic, Z. Djordjevic, B. Stojanovic, A. Disic, A new design of a two-stage cycloidal speed reducer, J. Mech. Des. 133 (8) (2011) 085001–085007. [19] Y. Zhang, W. He, X. Yang, Dynamical formulation and analysis of double crank ring-plate-type pin-cycloidal gear planetary drive, 2010 International Conference on Frontiers of Manufacturing and Design Science, ICFMD2010, December 11, 2010–December 12, 2010, Trans Tech Publications, Chongqing, China, 2011. [20] J. Liu, B. Chen, S. Matsumura, C. Li, Design of a novel cycloid drive with a cycloid-arc gear and analysis of its meshing characteristic, J. Adv. Mech. Des. Syst. Manuf. 6 (2) (2012) 310–322. [21] B. Zhu, W. Qin, J. Liu, Y. Fu, Simulation and analysis of dynamical transmission precision of 2K-V cycloid pin gear reducerbased on multi-body system dynamics, 2011 International Conference on Advanced Design and Manufacturing Engineering, ADME 2011, September 16, 2011–September 18, 2011, Trans Tech Publications, Guangzhou, China, 2011. [22] N. Kumar, Investigation of drive-train dynamics of mechanical transmissions incorporating Cycloidal drives(PhD Thesis) Science and Engineering Faculty, School of Chemistry Physics and Mechanical Engineering, Queensland University of Technology, Brisbane, Australia 2015, p. 316. [23] S. Kuzmanovic, M. Rackov, K. Rafa, Benefits of the application of cycloidal backlash gear reducers with more eccentric shafts, Balk. J. Mech. Transm. 2 (1) (2012) 33–38. [24] Regan, S., Power-transmitter, 1895, Google Patents [US Patent # 546,249]. [25] Williams, H.D., Power-transmission gearing, 1909, Google Patents [US Patent # 908,529]. [26] Harrison, W.L., Mechanism for transmitting rotary motion, 1910, Google Patents [US Patent # 978,371]. [27] Hatlee, D.H., Speed-reducer, 1916, Google Patents [US Patent # 1,192,627]. [28] Braren, L.K., (1931), Production of Cycloidal Curves, U.S. Patent, # 1,817,405. [29] Braren, L.K., (1932), Gear Transmission, U.S. Patent, # 1,867,492. [30] Fliesberg, E.A.M., Transmission gear, 1936, Google Patents [US Patent # 2,049,696]. [31] Jackson, J.A., (1949), Reduction Gear, U.S. Patent, # 2,475,504. [32] McIver, W.K., (1950), Gear Transmission, U.S. Patent, # 2,508,121. [33] Menge Sr, T.L., (1961), Speed Reducer, U.S. Patent, # 2,972,910. [34] Sundt, E.V., (1962), Differential Gear Reducer, U.S. Patent, # 3,037,400. [35] Braren, R., (1963), Planetary Gear, U.S. Patent, # 3,073,184. [36] Lee, R., (1964), Speed Reducers, U.S. Patent, # 3,129,611. [37] J.G. Blanche, D.C.H. Yang, Cycloid drives with machining tolerances, J. Mech. Transm. Autom. Des. 111 (Compendex) (1989) 337–344. [38] F.L. Litvin, Theory of Gearing, National Aeronautics and Space Administration, Scientific and Technical Information Division, 1989. [39] F.L. Litvin, A. Fuentes, Gear Geometry and Applied Theory, second ed. Cambridge University Press, New York, 2004 xvi (800 pp.). [40] Y.-W. Hwang, C.-F. Hsieh, Determination of surface singularities of a cycloidal gear drive with inner meshing, Math. Comput. Model. 45 (Compendex) (2007) 340–354. [41] Z. Ye, W. Zhang, Q. Huang, C. Chen, Simple explicit formulae for calculating limit dimensions to avoid undercutting in the rotor of a cycloid rotor pump, Mech. Mach. Theory 41 (2006) 405–414. [42] B. Chen, T.T. Fang, C.Y. Li, S.Y. Wang, Gear geometry of cycloid drives, Sci. China Ser. E Technol. Sci. 51 (5) (2008) 598–610. [43] Y. Meng, C. Wu, L. Ling, Mathematical modeling of the transmission performance of 2K-H pin cycloid planetary mechanism, Mech. Mach. Theory 42 (7) (2007) 776–790. [44] Y.-s. Sun, T.-m. Guan, Theory analysis of forces on two teeth difference cycloid gear drive, J. Anshan Steel Inst. 24 (5) (2001) 346–349. [45] Y.-s. Sun, T.-m. Guan, G.-s. Wang, L. Liu, Analysis of the many teeth difference cycloid gear drive, J. Dalian Railw. Inst. 22 (2) (2001) 23–26. [46] C. Wan, W. Zhao, L. Li, Optimization of the curve parameters within the pin gear housing in two teeth difference cycloid speed-reducer, Jixie Gongcheng Xuebao/ Chin. J. Mech. Eng. 39 (6) (2003) 124–127 (+134). [47] S. Yingshi, G. Tianmin, The modeling and simulation method to calculate force in the equivalent substitution flank profile two tooth difference cycloid pin gear reducer cycloid gear, Digital Manufacturing and Automation (ICDMA), 2010 International Conference on, 2010. [48] Y. Sun, T. Guan, Modeling and mechanical simulation of cycloid gear for two teeth difference cycloid pin reducer, Appl. Mech. Mater. 43 (9) (2010) 9–13 (Advance in Mechatronics Technology). [49] J.W. Sensinger, Unified approach to cycloid drive profile, stress, and efficiency optimization, J. Mech. Des. 132 (2) (2010) 024503. [50] J.W. Sensinger, J.H. Lipsey, Cycloid vs. harmonic drives for use in high ratio, single stage robotic transmissions, Robotics and Automation (ICRA), 2012 IEEE International Conference on, 2012. [51] J.W. Sensinger, Efficiency of high-sensitivity gear trains, such as cycloidal drives, J. Mech. Des. 135 (7) (2013) 0710061–0710069. [52] E.P. Pollitt, Some applications of the cycloid machine design, ASME J. Eng. Ind. 1960 (November, 1960) 407–414. [53] B. Borislavov, I. Borisov, V. Panchev, Design of a planetary-cyclo-drive speed reducer cycloid stage, geometry, element analyses, Mechanical Engineering, Vol. 90, Linnaeus University, Växjö, Sweden, 2012. [54] M. Blagojević, V. Nikolić-Stanojević, N. Marjanović, L. Veljović, Analysis of Cycloid Drive Dynamic Behavior, in Scientific Technical Review, Vojnotehnički Institut Faculty of Mechanical Engineering, University of Kragujevac, Belgrade, Serbia, 2009 52–56. [55] K.-H. Kim, C.-S. Lee, H.-J. Ahn, Torsional rigidity of a two-stage cycloid drive, Trans. Korean Soc. Mech. Eng. A 33 (Compendex) (2009) 1217–1224. [56] Y. Zhang, J. Xiao, W. He, Dynamical Formulation and Analysis of RV Reducer, Engineering Computation, 2009. ICEC '09. International Conference on, 2009. [57] Y. Zhang, W. He, J. Xiao, Dynamical model of RV reducer and key influence of stiffness to the nature character, Information and Computing (ICIC), 2010 Third International Conference on Engineering Computation, 2010. [58] L. Lei, T. Guan, J. Li, Rigidity analysis of FA pin-cycloid drive, 2011 International Conference on Materials and Products Manufacturing Technology, ICMPMT 2011, October 28, 2011 - October 30, 2011, Trans Tech Publications, Chengdu, China, 2011. [59] A. Mihailidis, I. Nerantzis, a new system for testing gears under variable torque and speed, Recent Pat. Mech. Eng. 2 (3) (2009) 179–192. J. [60] Peeters, D. Vandepitte, P. Sas, Analysis of internal drive train dynamics in a wind turbine, 7th German Wind Energy Conference, Wiley Interscience, Wilhelmshaven, Germany, 2004. [61] J. Peeters, Simulation of dynamic drive train loads in a wind turbine, Faculteit Ingenieurswetenschappen Departement Werktuigkunde Afdeling Productietechnieken Machinebouw en Automatisering, Katholieke Universiteit Leuven, Leuven (Heverlee), Belgium 2006, p. 336.

198

N. Kumar et al. / Mechanism and Machine Theory 105 (2016) 185–198

[62] V. Kosse, Experimental investigation of torsional stiffness of mechanical drives with large reduction ratio - part 1: cyclodrives, in: J. Mathew, et al., (Eds.), Proceedings of the 10th Asia-Pacific Vibration Conference, Queensland University of Technology, Brisbane, Queensland 2003, pp. 318–322. [63] R. Dhaouadi, F.H. Ghorbel, P.S. Gandhi, A new dynamic model of hysteresis in harmonic drives, IEEE Trans. Ind. Electron. 50 (6) (2003) 1165–1171. [64] H.-H. Lee, Finite Element Simulations with ANSYS Workbench 12, Schroff Development Corporation (SDC Publications), 2010 586. [65] Z.M. Levina, D.H. Reshetov, Contact Stiffness of Machines, Mashinostroyenie, Moscow, 1971 (in Russian). [66] G.S. Maslov, Calculations of Shaft Vibrations, Mashinostroyenie, Moscow, 1980 (in Russian). [67] B.K. Chistiakov, Dynamics of Piston and Combined Internal Combustion Engines, Mashinostroyenie, Moscow, 1989 (in Russian). [68] O.P. Mihaylov, Dynamics of Electromechanical Drives of Machine Tools, Mashinostroyenie, Moscow, 1989 (in Russian). [69] I.A. Birger, Y.G. Panovko, in: I.A. Birger, Y.G. Panovko (Eds.), Strength, Stability, Vibrations, (Reference Book in 3 Volumes), Mashinostroyenie, Moscow, 1968 (in Russian). [70] G.V. Kreinina, in: G.V. Kreinina (Ed.), Dynamics and Control of Machines, Mashinostroyenie, Moscow, 1988 (in Russian). [71] S.S. Rao, Mechanical Vibrations, second ed. Addison-Wesley Publishing Company, 1990. [72] W.C. Young, R.G. Budynas, Roark's Formulas for Stress and Strain, seventh ed., General Engineering Series, McGraw-Hill International Edition, 2002. [73] S.P. Timoshenko, On the correction factor for shear of the differential equation for transverse vibrations of bars of uniform cross-section, Pilosophical Magazine 1921, p. 744. [74] W.J. Palm III, System Dynamics, second ed. McGraw Hill Higher Education, 2010.