Multibody Dynamics Formulation for Modeling and Simulation of Roller Chain Using Spatial Operator

This paper proposes a new approach for modeling and simulating roller-chain drives using multibody formulations. The presented approach models the chain links as spatial decoupled dynamic bodies. The chain links are modeled using two different dynamic representations for the pin-link and the bushing-link. The bushing-link is modeled with two descendent dynamic bodies to represent the rollers. The adjacent pin and bushing links are connected by compliant bushing force elements. An efficient search algorithm is used to detect the contact between the roller and the sprocket teeth while a nonlinear force module is used to predict the contact force. A generalized sprocket representation is used to model the sprocket. The spatial motion of the chain links allows the out-of-plane vibrational motion of the links as well as simulating sprockets misalignment. Using the compliant connection between links avoids using the iterative calculation needed to satisfy the joint constraint equations leading to more efficient calculation scheme. The nonlinear dynamic equations of motion are solved using recursive approach. Complex roller chains drives, bicycle chain and conveyance systems can be easily modeled and analyzed using the proposed approach.


Introduction
Recent design advances in roller chain opened new areas of simulation-based-optimization to improve the product quality.Xi [1] utilized multibody dynamic simulation of roller chains to identify the noise sources in diesel engines.The effect of the chain impact on the noise was also investigated by Zheng [2] where the noise level was correlated to the rigid cylinder acceleration.The finite element approach was used and results were compared to experimental results.Ryabov [3] presented an optimization approach for selecting the sprocket teeth based on chain drive cost and service life.
Purdue [4] studied the effect of varying the bike configuration on improving the transmitted pedaling force using commercial program in his investigation.Assuming the chain motion is planar, the simulation software was used to study the effect of non-circular sprocket on minimizing the force exerted by the rider.Approximate method of modeling the chain as a travelling heavy link was used by Rivola [5] to investigate how the transmitted power is fluctuated as a result of the sprocket eccentricity and polygonal effect.Altered conjugated sprocket tooth profile was presented by Wang [6] with the objective of reducing the polygonal action and impact in high speed drive chain systems.Another investigation utilized multibody numerical results to study the impact of using asymmetric sprocket tooth.The effect of using new composite materials for links in conveyer chains was a Corresponding author: momar@taibahu.edu.sa,momar23@gmail.cominvestigated by Kadam [7].The contact models in cylindrical objects were reviewed by Pereira [8] and the application range of the nonlinear contact forces as function of the indentation were examined.
Choi [9] examined the roller chains with tensioners and the chain polygonal effect across the chain span length.Veikos [10] modeled the chain as lumped masses connected with simple spring-damper element to study the chain vibration.Modeling the link-sprocket contact was investigated by Liu [11] using impulsive function.Sheu [12] proposed a vector approach for analyzing the friction kinetostatics and the transmission efficiency.The forces in the chain links were measured experimentally by Kidd et al [13] and showed that sprocket misalignment leads to localized bending of the chain links.Marshek et al [14] showed that the misalignment leads to rubbing of the links with the sprockets.
Pin-Link Body (PLB) This review shows the importance and need to develop an efficient multibody simulation capability that should be able to address the out-of-plane link motion, sprocket misalignment, comprehensive algorithms for contact detection between links and sprockets, and generalized description of the sprocket tooth and profile.In earlier investigation, the author [15] presented a novel approach for modeling chain drives without taking the roller into account.This could lead to unrealistic high friction and contact forces between the sprocket and the links.This paper presents a modified formulation that accounts for the roller dynamics.The system nonlinear dynamic equations are modified accordingly and presented in details.
This paper is organized as follows.In the following section, the link kinematic description and the equations of motion are derived.In section 3, the compliant bushing force element is presented, the roller-bushing contact force module is also discussed.A brief discussion of the nonlinear link-sprocket contact force model is presented.Finally section 4 presents an example of the proposed approach.

Link equations of motion
The chain links will be divided into two distinct dynamic body types; the pin-link body (PLB) and the bushing-link body (BLB).The PLB is composed of two side plates press-fitted with two pins to form one rigid body.Each PLB is considered to be kinematically independent from other links and will be modeled as a rigid body with 6-DOF relative to the global frame of reference (GFOR).Similarly, the BLB is composed of two side plates press-fitted with two bushings and form one rigid body.The rollers are considered as dynamic bodies connected to BLB with planar joint with 3 DOF (two translational and one rotational DOF).In order to increase the computational efficiency, the two rollers and the BLB will be treated as on super body link (SBL).The BLB will be modeled with 6-DOF relative to the GFOR while the two rollers are modeled with 3-DOF each relative to the BLB.The SBL will have a total of 12-DOF.
The vector of spatial momentum of the PLB expressed in the GFOR could be differentiated to obtain the body equation motion.The spatial momentum vector of pin-link-body j, could be calculated as follows: where 0 ,0 j P is the spatial momentum vector of body j, 0 j M is the inertia matrix transformed to GFOR, m j is the body mass, I is 3x3 identity matrix, 0 j J is the 3x3 matrix representing the moment of inertia defined at a triad located at the body center of mass and transformed to GFOR, and 0 0 j v is the spatial velocity vector in GFOR, and 0 0 j r and 0 0 j Z are the translational and angular velocity vectors in GFOR.
Differentiating (1) leads to the BLB equation of motion as follows: , where where 0 j G is the spatial force vector including gravity defined in GFOR, 0 j W and 0 j g represents the external moments and forces acting on the body.The mass matrix in (2) could be inverted to calculate the body absolute accelerations which can be directly integrated to obtain the velocity and displacement of the PLB as follows: In the case of SLB, the rollers are considered as descendants of the BLB.The relative velocity of child bodies relative to the parent could be expressed in local and global coordinate systems using the joint variable velocity as follows [15]: h is a 6x3 joint influence coefficient matrix, as shown in Figure 1, k q is 3x1 vectors of joint relative velocities, and 0 0 k k X is the spatial transformation matrix between the body k and GFOR.The joint influence coefficient matrices could also be used to project the spatial force vectors, F k , into the joint subspace force vector, Q k , as follows: Using the recursive formulation, the absolute velocity, of the roller bodies, k and l, could be calculated using the BLB parent, j, and the relative joint velocities as follows: Differentiating velocity equations, the acceleration of the child bodies across the joint DOF could be obtained as follows: where k q is 3x1 vectors of joint relative velocity and acceleration, and 0 0k γ is the vector of quadratic velocities.
The velocity and accelerations of the SLB could be written in a compact form as follows: where T a is called the system topology matrix and H a is the system influence coefficient matrix.The equation of motion of the connected bodies j, k, and l could be driven from the free body diagram.For body j, the equation of motion in its local coordinate system could be written as: These equations could be arranged in a matrix form as follows: G vP a G v P (10) 01003-p.3 Equation ( 10) represent the SLB equation of motion that contains the bodies Cartesian accelerations and reaction forces between the rollers and the BLB.It could be written in a compact matrix form as: Using Equations ( 5) , ( 8) and (11) , the system equations could be rearranged as follows: The system coefficient matrix in the left hand side of equation ( 12) could be inverted to solve for the unknown joint variables accelerations as follows: The detailed block-entries of the inertia matrix and generalized force vector are written as follows: Noting that constructing the H matrices doesn't involve computations, rather, it is picking the right column from the spatial transformation matrix.This 12x12 matrix could be easily evaluated and inverted to compute the motion derivatives.Using Equation (3) and Equation ( 14) all the accelerations of the links in the chain could be calculated and then integrated very efficiently.The connections between the pin-links and the bushing-links are modeled using bushings force elements.The bushing stiffness accounts for the pin-link flexibility and the bushing-link flexibility.

Bushing, sprocket and contact force elementS
The repeated pattern of pin-link and bushing-link connection allows using the same bushing parameters over the chain.As shown in Figure 2, two markers are located in each link O pi1 and O pi2 , to define the two pin-centers relative to the body reference-frame, O pi .Similarly, two markers, O bi1 and O bi2 , are used to define the bushing center relative to the bushing-link reference frame, O bi .The bushing force element connects the pin center marker to the corresponding bushing center marker.At the beginning of the simulation when there is no deflection or relative rotation between the pin-link and bushing-link, the pin-center and bushing-center markers are considered to be coincident and parallel.To calculate the elastic forces of the bushing element, the translational deflection and relative orientation between the pin-center marker and the bushing-center marker need to be calculated every time step.
Due to the bushing elastic deformation and structural material damping, the bushing force element develop a linear elastic forces and torques proportional to the bushing deformation.Also, calculating the approximated viscous damping forces and torques proportional to the relative deformation velocity could be done.The elastic torque component about the bushing axis is assumed to be zero in order to allow relative rotation.The bushing forces and torques can be computed as follows: , where K x , K y , K z , K tx , K ty , and K tz , the stiffness coefficients associated with translational and rotational displacement respectively and C x , C y , C z , C tx , C ty , and C tz the damping coefficients associated with translational and rotational velocities respectively.
Estimating the contact force between the roller and the bushing is crucial in studying the durability and wear of the two parts.The contact between the bushing and the roller is modelled using the Hertz contact model.The roller has relative planar displacement and rotation with respect to the BL.The magnitude of penetration, G, could be calculated from the components of relative displacements.Employing Johnson's formula based on Hertz theory [16], the normal contact force, P, could be calculated as function of material properties * Using the relative velocities the penetration rate could be calculated and used to determine the damping force.The magnitude of the friction force and torque could be calculated from the normal force, coefficient of friction, and the relative velocity.
The contact between the sprocket teeth and chain links is the main motion transmission mechanism in chain drive system.The contact forces are modeled to prevent penetration of the chainrollers into the sprocket teeth.The contact force magnitude depends on tooth curvature at the contact point and on the localized mutual deformation of the roller and the tooth.Predicting the contact location requires accurate calculation of the roller center position relative to the tooth profile.The friction force due to roller sliding over the tooth requires accurate calculation of the relative velocity between the roller and sprocket tooth at the point of contact.In the current proposed approach, the sprocket tooth is modeled using five line segments to approximate the tangent lines in plane, shown in Figure 2 (b).More accurate representation of the tooth profile could utilize spline to smooth the tooth geometry.The kinematic model for contact point detection is shown in Figure 2 (c).
The relative position and relative velocity vectors between the roller center and the sprocket center, , s sb s r and , s sb s r , are calculated.Knowing the tangent line segment and the roller radius, the magnitude of penetration could be calculated and the contact forces are estimated.In order to calculate the normal and sliding velocity components, , s sb s r the relative velocity vector needs to be transformed using the inclination angle φ as follows:

Demonstrating example
This example represents a bicycle chain that contains 56 PLB and 56 SLB.The multibody model of the bike is composed of: front drive sprocket, wheel sprocket, the lower idler sprocket, and the upper idler sprocket.The bicycle chassis is fixed to the GFOR and all sprockets are attached to the chassis with revolute joints.Every PLB and SLB body is referenced to the ground with full 6-DOF joints.The drive sprocket attachment position is offset from the wheel sprocket attachment position by 30mm in the lateral direction.The chain link-pitch is 12.7mm, the bushing outer radius is 3mm, the pin radius is 2.1mm, and roller radius is 4mm.The sprocket tooth thickness is 3.2mm.A kinematic driver is applied to the driving sprocket ramping up the speed from zero to 1 rev/s in 3 seconds, as shown in Figure 3.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences

Figure 3
(a) Drive sprocket speed (b) Lateral acceleration of chain link (c) Bushing tension force Each link travels a distance of 1.422m around the sprockets to complete one cycle.Sample simulation results are shown in Figure 3.
Figure 3 (b) shows the links' lateral accelerations.The axial tension force in link 1 during the simulation is shown in Figure 3 (c).
17)The sliding velocity could be computed as: The nonlinear sliding friction force is a function of the normal elastic force and the friction coefficient.The friction coefficient can be considered as smooth function of the sliding velocity as follows: P is a predefined maximum sliding velocity at which the friction force reaches maximum value, S P is the friction force softening.Then contact forces are now applied in the sprocket reference frame and the link reference frame as: