Chain Drive Simulation Using Spatial Multibody Dynamics

Abstract

This paper presents an efficient approach for modeling chain derives using multibody dynamics formulation based on the spatial algebra. The recursive nonlinear dynamic equations of motion are formulated using spatial Cartesian coordinates and joint variables to form an augmented set of differential-algebraic equations. The spatial algebra is used to express the kinematic and dynamic equations leading to consistent and compact set of equations. The connectivity graph is used to derive the system connectivity matrix based on the system topological relations. The connectivity matrix is used to eliminate the Cartesian quantities and to project the forces and inertia into the joint subspace. This approach will result in a minimum set of equation and can avoid iteratively solving the system of differential and algebraic equations to satisfy the constraint equations. In order to accurately capture the full dynamics of the chain links, each link in the chain is modeled as rigid body with full 6 degrees of freedom. To avoid singularities in closed loop configurations, the chain drive is considered a kinematically decoupled subsystem and the interaction between the links and other system components is modeled using force elements. The out-of-plane misalignment between the sprockets can be easily modeled using a compliant force element to model the joints between each two adjacent links. The nonlinear three dimensional contact forces between the chain links and the sprockets are modeled using elastic spring-damper element and accounts for the sliding friction. The proposed approach can be used to model complex drive chain, bicycle chain as well as conveyance systems. Results show that realistic behavior of the chain as well as out-of-plane vibration can be easily captured using the presented approach. The proposed approach for chain drive subsystem could be easily appended to any other multibody simulation system.

1. Introduction

Chain drives are very common reliable mechanisms that have been used over centuries for power transmission and conveyance. The design of chain drive system may vary based on the application, the rated power, and the operating conditions. Valve trains in internal combustion engines, bicycle and motorcycles, conveyer belts, over-head conveyers, and hoist mechanisms are typical applications of the chain drive systems. Majority of the commonly used chain drives are made of steel and the links are considered rigid. The roller chain systems are usually composed of articulated links forming a closed kinematic chain. A typical power transmission roller chain is composed of outer links (pin links), inner links (bushing links), pins, bushings, and rollers, as shown in Figure 1. The pins are typically press-fitted to the two outer plates to form the pin link while the bushing is press-fitted to the two inner plates to form the bushing link. The bushing can rotate freely around the pins while the rollers are freely assembled into the bushing.

Figure 1:

Typical roller chain system.

In the simplest chain drive system, the chain is driven by a driving sprocket converting the driving torque into a tensile force that will drive the output sprocket and/or move the load. Figure 2 shows a typical arrangement of the chain drive system used in bicycles. Variations from this configuration may include the following: adding idle sprockets, using chain support guides and chain tensioners, and moving loads. Tensioners are usually introduced to eliminate the chain slag, reduce noise and vibration, and improve the efficiency. Chain guides are used to reduce the in-plane vibration and eliminate the out-of-plan chain motion. In many cases where guides are used, the chain motion could be accurately modeled by planar motion. Out-of-plane motion may be excited in other cases as follows: when the axes of drive and driven sprockets are not perfectly parallel, if the designer introduced lateral shift between the drive and driven sprocket planes, in the case of high speed sprockets with large slack, when excitation is transmitted from other machine components, and so forth. Also, during operation, multiple contacts between the sprocket teeth and the chain rollers lead to high impact forces and sudden variation of link velocities and create lateral and out-of-plane vibration. In the above mentioned cases, planar dynamic simulation will not be sufficient to model and capture the chain dynamics. The joint connection between the roller and bushing typically has clearance leading to more complex dynamic contact.

Figure 2:

Typical bicycle chain arrangement.

In automotive applications, the chain drive systems are designed to eliminate out-of-plane motion using chain guides and tensioners. In such cases of planar motion, 2D simulation of the chain drive system could be sufficient to study the dynamic behavior and extract the necessary information about the velocities, accelerations, contact forces, and friction. In other cases where chain guides cannot be introduced in the design, the chain drive links can experience lateral vibration. In the case of bicycle chain, for example, the chain must undergo out-of-plan motion in order to allow the driver to shift between gears. Due to the lateral misalignment between the drive sprocket plane and driven sprocket plane, the chain motion will not be planar and the links will be subjected to bending and side contact with the sprockets. The chain engagement-disengagement with the sprockets during shifting are very complicated processes and involve sever friction and rubbing between the chain plates and sprockets' faces. If the chain is stiff, shifting the chain from one sprocket to another may cause glitches resulting in rough acceleration that will cause discomfort to the rider. For this reason, the chain could be designed with lateral flexibility or induced elasticity. Polygonal action is another phenomenon that can cause meshing noise and lead to wear and premature chain fatigue life [1–5].

In general the roller chain systems exhibit very complex dynamics behavior [2]. For this reason it was impossible to find general analytical procedures that can describe thoroughly the chain drive dynamic problems. In the above mentioned cases, 2D dynamic model is not sufficient to model the transitional dynamic during speed shift. Due to the large number of bodies in the chain, the number of degrees of freedom will become significantly high and the dynamic simulation will require huge computational resources. The objective of this research is to device a formulation that will simulate the spatial chain dynamic motion while reducing the computational effort. This research aims to develop an algorithm for modeling the chain drive system that can be appended to any multibody dynamics platform.

In this proposed approach, links, sprockets, and rollers are modeled as rigid bodies. The spatial motion of each link will be modeled using full 6 degrees of freedom (DOF). As the number of chain links and the associated DOF increase, the computational effort required to factor the inertia matrix will significantly increase [6, 7]. In order to avoid increasing the computational effort, the proposed approach is designed to decouple the chain links equations of motion from the main system equation of motion as will be discussed in the following sections. In the following section, some of the existing chain modeling approaches will be reviewed. Section 3 summarizes the joint-based multibody dynamics formulation approach and the system dynamic equations will be formulated using spatial algebra. The proposed approach to model the chain links will be presented and discussed. The interaction between the rollers and the sprocket teeth is represented as contact force model and will be discussed. The compliant model for the connection between the links will be also presented. The recursive solution algorithm is presented in Section 4. Section 5 present the dynamic model of a simple bicycle chain drive system and the simulation results will be presented and discussed. Section 6 summaries the work and presents the concluding remarks.

2. Chain Modeling Approaches

Hippmann et al. [8] presented a 2D chain drive model based on the relative coordinates. The proposed approach in [8] was developed to study the gross motion of the chain along its trajectory. The model was reported to be reliable and efficient in simulating highly complex dynamical behavior of chain drives without numerical damping. The model also embedded substructure techniques for the chain drive model in multibody system models of combustion engines. To solve the equation of motion, implicit integration approach was adapted.

A mathematical model for motorcycle drive train was presented by Gamba [1]. The model can be implemented in a multibody motorcycle model for handling and maneuverability studies and long time-scale analysis. It was reported that this model is computationally fast. Pedersen [2] presented a novel formulation for the simulation of the dynamics of roller chain drives using a continuous contact force method. Modeling the contact surface between the rollers and sprocket has shown to be an important issue regarding the numerical stability of the simulation program. The proposed model concluded that using the real tooth profile leads to superior results over other simplified methods.

Pereira et al. [3, 4] presented a multibody methodology to address the kinematic and dynamic effects of roller chain drives. In his proposed approach the chain was modeled as a collection of rigid bodies, connected to each other by revolute clearance joints. Each clearance revolute joint, representing the connection between pair of links, is made up of the pin link/bushing link plus the bushing link/roller pairs, if the chain is a roller chain.

Although, the above mentioned models were reported to be efficient, they lack the full dynamics of the multibody systems. When each link is modeled as rigid body with 6-DOF the computational effort will significantly increase. Tracked vehicle dynamics analysis algorithms are very similar to chain drive algorithms and track modeling approaches and could be adapted to model chain drive systems. Rubinstein and Hitron [9] created a model which incorporates a detailed description of the track links and the dynamic interaction with other components. Each track shoe was considered as a rigid body and is connected to its neighboring track shoes via a kinematic revolute joint constraint. Ryu et al. [10] created a three-dimensional multibody model with compliant track-shoe connections. The revolute joints which connect track shoes with their neighbors are replaced with compliant force elements which are described by stiffness and damping values. Ryu et al. [8, 10, 11] built on the methodology proposed in [10] by further developing the contact force models that include link contact with the sprocket.

3. Kinematic and Dynamic Equations of Motion

The connectivity graph is one way to graphically represent the kinematic model and the topological connectivity between the different bodies in the multibody system. Figure 3 shows a simple bike model and the corresponding connectivity graph. In the connectivity graph, each body is connected to one parent through an arc joint that can allow one or more DOF. The root body is connected/referenced to the ground and the system under consideration has no kinematic loops. While each body can have only one parent it could have one or more descendent (child bodies). In the connectivity graph, the root body is numbered as 0 while the other bodies are numbered from 1 to n_b consecutively. The body and the joint connecting it to its parent are given the same number [6, 7]. Bodies (links) not connected to the ground or to any parent with arc joint are considered floating body (FB) with 6-DOF joint. Based on the joint connectivity, a parent child list could be developed and stored to be used later in the recursive calculations [6, 7, 12–14].

Figure 3:

Simple connectivity graph diagram for a bike.

The following development refers frequently to the triad or marker. The triad (marker) is composed of three orthogonal unit vectors used to measure relative and absolute motion between different points. The triad position is defined by the position of its origin while the orientation is defined by a set of rotational parameters. The child body is joined to its parent through two triads; one triad is fixed in the parent side while the other triad is on the child side, as shown in Figure 4. The triad in the child side is considered as the output triad and is used as the reference frame of this body. The triad in the parent side is considered as the input triad to the child body, connector triad. The position and orientation of the connectors are measured relative to the parent reference triad through the spatial transformation matrix. The joint DOF or joint variables are defined as the allowed relative motion between the two connecting triads. Relative motion between the triads could be translation or rotation along one of the connector triad axes. Massless triads could be inserted between the two bodies to represent joints with more than one DOF. The joint displacement, velocities, and accelerations are used as the body state variables. The Cartesian displacement, velocity, acceleration vectors, and the joint reaction forces are considered as augmented algebraic variables.

Figure 4:

Kinematic connection and the joint influence coefficient matrix.

3.1. Kinematic Equations

The spatial velocity vector of a body can be represented by a combination of two 3 × 1 vectors $\underset{\_}{\omega }$ that represents the body angular velocity and ${\underset{\_}{\dot{r}}}_o$ represents the linear velocity of the body-fixed point that is passing through an arbitrary given point, O, in the motion space [6, 11, 15–18]. The spatial velocity vector can be written as:

${\mathbf{v}}^i={\left[\begin{array}{cc}\underset{\_}{\omega }& {\underset{\_}{\dot{r}}}_o\end{array}\right]}^T.$ (1)

The spatial velocity expressed in coordinate system (CS) at j is transformed to a CS i through the spatial transformation matrix in motion space as follows:

${\mathbf{v}}^i={}_{}{}^m{X}_{ij}^{ij}{\mathbf{v}}^j,$ (2)

where ${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{ij}^{ij}$ is the spatial transformation matrix, the left superscript m refers to motion space, i and j indicate the coordinate systems. The spatial velocity vector ${\mathbf{v}}_{ij,k}^i$ is defined in a specific coordinate system i (superscript i) measuring the velocity difference between two triads i and j (subscript ij) and is defined with respect to a certain observer at k (subscript after comma k). The transformation matrix ${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{ij}^{ij}$ is the result of rotation followed by a translation as follows:

${}^m{X}_{ij}^{ij}=L_{ij}^{ij}{R}^{ij}=\left[\begin{array}{cc}I& \mathrm{0}\\ {\underset{\_}{\tilde{r}}}_{ij}^i& I\end{array}\right]\left[\begin{array}{cc}{\underset{\_}{R}}^{ij}& \mathrm{0}\\ \mathrm{0}& {\underset{\_}{R}}^{ij}\end{array}\right]=\left[\begin{array}{cc}{\underset{\_}{R}}^{ij}& \mathrm{0}\\ {\underset{\_}{R}}^{ij}{\underset{\_}{\tilde{r}}}_{ij}^i& {\underset{\_}{R}}^{ij}\end{array}\right],$ (3)

where L _ij ^ij is the 6 × 6 spatial translation matrix, ${\underset{\_}{r}}_{ij}^i$ is the 3 × 1 translational displacement vectors of CS j relative to CS i defined in the CS i, ${\underset{\_}{\tilde{r}}}_{ij}^i$ is 3 × 3 skew symmetric matrix representing the cross-product operation of ${\underset{\_}{r}}_{ij}^i$ , R ^ij is the 6 × 6 spatial rotation matrix, ${\underset{\_}{R}}^{ij}$ is the 3 × 3 transformation matrix relating the coordinate frames i and j. It should be noted that the general spatial transformation matrix is not orthonormal. Displacements, velocities and accelerations can be expressed in motion space while forces, torques, momentums and inertia can be expressed in the dual or function space [19]. Vectors in the function space can be transformed using ${}_{}{}^f{X}_{ij}^{ij}$ , [6, 11, 15, 16]. The spatial transformation matrix in function space ${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^f{X}_{ij}^{ij}$ can be written as follows:

${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^f{X}_{ij}^{ij}={\left({}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{ij}^{ij}\right)}^{-T}={\left(L_{ij}^{ij}{R}^{ij}\right)}^{-T}.$ (4)

The relative spatial velocity vector between body i and body j defined in CS of body j observed at the origin of j can be written as follows:

${\mathbf{v}}_{ij,j}^j=\left[\begin{array}{c}{\underset{\_}{\omega }}_{ij}^j\\ {\underset{\_}{\dot{r}}}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}ij,j}^j\end{array}\right],$ (5)

where ${\underset{\_}{\omega }}_{ij}^j$ is the time derivative of the orientation variables between CS i and CS j defined in CS j and ${\underset{\_}{\dot{r}}}_{ij,j}^j$ is the translational velocity between CS i and CS j. The relative spatial velocity vector across primitive joint jbetween bodies i and j, v ^ij , can be evaluated using as follows [6, 14, 20]:

${\mathbf{v}}_{ij,j}^j={\mathbf{h}}_{j,j}^j{\dot{\mathbf{q}}}_j,$ (6)

where ${\mathbf{h}}_{j,j}^j$ is a 6 × n _d matrix representing the joint influence coefficient matrix defined in the CS of body j, n_d is the number of joint degrees of freedom, and ${\dot{\mathbf{q}}}_j$ is an n_d vector of the velocities of the joint variables. ${\mathbf{h}}_{j,j}^j$ is composed of the columns of the identity matrix corresponding to the primitive joint DOF. The spatial velocity vector of body j can be referenced to the inertial frame of reference through its parent body i as follows:

${\mathbf{v}}_{0j,j}^j={}_{}{}^m{X}_{ji}^{ji}{\mathbf{v}}_{0i,i}^i+{\mathbf{v}}_{ij,j}^j={\mathbf{v}}_{0i,j}^j+{\mathbf{h}}_{j,j}^j{\dot{\mathbf{q}}}_j,$ (7) ${\mathbf{v}}_{0j,0}^0={}^m{X}_{0j}^{0j}{\mathbf{v}}_{0j,j}^j.$ (8)

The velocities of the bodies in the kinematic tree can be calculated recursively according to the parent child list [6, 14] as will be shown later.

The influence coefficient matrix could be transformed using the spatial transformation matrix as follows:

${\mathbf{h}}_{j,i}^i=X_{ij}^{ij}{\mathbf{h}}_{j,j}^j.$ (9)

Since ${\mathbf{h}}_{j,j}^j$ contains columns of the identity matrix, ${\mathbf{h}}_{j,i}^i$ will contain the corresponding columns from the transformation matrix X _ij ^ij . Note that the joint influence coefficient matrix could be written as

${\mathbf{h}}_{j,0}^0={}_{}{}^m{X}_{0j}^{0j}{\mathbf{h}}_{j,j}^j.$ (10)

Since each joint could have six DOF, ${\mathbf{h}}_{j,j}^j$ could be considered as a projection operator that selects only the free joint DOF out of the possible six DOF to contribute to the body velocity. In the same context, a rejection operator ${\stackrel{-}{\mathbf{h}}}_{j,j}^j$ could be used to filter out the locked joint DOF. ${\stackrel{-}{\mathbf{h}}}_{j,j}^j$ is a matrix that contains the remaining 6-n_d columns of the identity matrix after extracting the joint influence coefficients ${\mathbf{h}}_{j,j}^j$ . A permuted joint matrix could be defined as ${\overline{\mathbf{H}}}_{j,j}^j=\left[\begin{array}{cc}{\mathbf{h}}_{j,j}^j& {\stackrel{-}{\mathbf{h}}}_{j,j}^j\end{array}\right]$ . Similarly, the permuted joint influence matrix could be used to extract the working and nonworking reaction forces. Based on this permutation, the joint influence coefficient matrix could be used to project the Cartesian forces into the free joint subspace as follows:

$Q_j={\mathbf{h}}_{j,j}^j{\mathbf{F}}_{j,j}^j,$ (11)

where ${\mathbf{F}}_{j,j}^j$ is a 6 × 1 force vector in the Cartesian space and Q _j is n _d × 1 force project into the free DOF joint subspace.

The relative spatial acceleration vector could be obtained in the global coordinate system by differentiating the velocity vector in (8) as follows:

${\mathbf{a}}_{ij,\mathrm{0}}^{\mathrm{0}}={}_{}{}^m{X}_{0j}^{0j}{\mathbf{h}}_{j,j}^j{\ddot{\mathbf{q}}}_j+{}_{}{}^m{\dot{X}}_{0j}^{0j}{\mathbf{h}}_{j,j}^j{\dot{\mathbf{q}}}_j,$ (12)

where ${\ddot{\mathbf{q}}}_j$ is joint variables acceleration vector and ${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{\dot{X}}_{\mathrm{0}j}^{\mathrm{0}j}$ is the derivative of the spatial transformation matrix that could be written as

${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{\dot{X}}_{\mathrm{0}j}^{\mathrm{0}j}={\tilde{\mathbf{v}}}_{\mathrm{0}j,\mathrm{0}}^{\mathrm{0}}{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{\mathrm{0}j}^{\mathrm{0}j}={}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{\mathrm{0}j}^{\mathrm{0}j}{\tilde{\mathbf{v}}}_{\mathrm{0}j,j}^j.$ (13)

Substituting (13) and (10) into (12) we get

${\mathbf{a}}_{ij,\mathrm{0}}^{\mathrm{0}}={\mathbf{h}}_{j,0}^0{\ddot{\mathbf{q}}}_j+{\tilde{\mathbf{v}}}_{0j,0}^0{\mathbf{h}}_{j,0}^0{\dot{\mathbf{q}}}_j.$ (14)

The child acceleration is related to its parent acceleration by ${\mathbf{a}}_{\mathrm{0}j,\mathrm{0}}^{\mathrm{0}}={\mathbf{a}}_{\mathrm{0}i,\mathrm{0}}^{\mathrm{0}}+{\mathbf{a}}_{ij,\mathrm{0}}^{\mathrm{0}}$ , which can be rewritten using (14) as follows:

${\mathbf{a}}_{\mathrm{0}j,\mathrm{0}}^{\mathrm{0}}={\mathbf{a}}_{\mathrm{0}i,\mathrm{0}}^{\mathrm{0}}+{\mathbf{h}}_{j,0}^0{\ddot{\mathbf{q}}}_j+{\tilde{\mathbf{v}}}_{0j,0}^0{\mathbf{h}}_{j,0}^0{\dot{\mathbf{q}}}_j={\mathbf{a}}_{0i,0}^0+{\mathbf{h}}_{j,0}^0{\ddot{\mathbf{q}}}_j+{\mathbf{\gamma }}_{0j,0}^0,$ (15)

where ${\mathbf{h}}_{j,\mathrm{0}}^{\mathrm{0}}$ are selected columns of the spatial transformation matrix and ${\mathbf{\gamma }}_{0j,0}^0={\tilde{\mathbf{v}}}_{0j,0}^0{\mathbf{h}}_{j,0}^0{\dot{\mathbf{q}}}_j={\tilde{\mathbf{v}}}_{0j,0}^0{\mathbf{v}}_{0j,0}^0$ is the quadratic velocity term. The acceleration could be expressed as follows:

${\mathbf{a}}_{ij,j}^j={\mathbf{h}}_{j,j}^j{\ddot{\mathbf{q}}}_j+{\mathbf{\gamma }}_{0j,j}^j,\hspace{1em}\hspace{1em}{\mathbf{a}}_{0j,j}^j={\mathbf{a}}_{0i,j}^j+{\mathbf{h}}_{j,j}^j{\ddot{\mathbf{q}}}_j+{\mathbf{\gamma }}_{0j,j}^j.$ (16)

Rearranging the terms in (6), (7), and (8) we can write the velocity equation as follows:

$\left[\begin{array}{cc}I& \mathrm{0}\\ -{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{ji}^{ji}& I\end{array}\right]\left[\begin{array}{c}{\mathbf{v}}_{\mathrm{0}i,i}^i\\ {\mathbf{v}}_{\mathrm{0}j,j}^j\end{array}\right]=\left[\begin{array}{cc}{\mathbf{h}}_{i,i}^i& \mathrm{0}\\ \mathrm{0}& {\mathbf{h}}_{j,j}^j\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{q}}}_i^{}\\ {\dot{\mathbf{q}}}_j^{}\end{array}\right].$ (17)

Considering the bike topology of the system with multiple bodies, shown in Figure 3, and using a parent child list, we can calculate velocities as follows:

$\left[\begin{array}{cccccc}I& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{\mathrm{21}}^{\mathrm{21}}& I& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{\mathrm{32}}^{\mathrm{32}}& \mathrm{0}& I& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{\mathrm{42}}^{\mathrm{42}}& \mathrm{0}& \mathrm{0}& I& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{\mathrm{52}}^{\mathrm{52}}& \mathrm{0}& \mathrm{0}& I& \mathrm{0}\\ \mathrm{0}& -{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^m{X}_{\mathrm{62}}^{\mathrm{62}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& I\end{array}\right]\left[\begin{array}{c}{\mathbf{v}}_{\mathrm{01,1}}^1\\ {\mathbf{v}}_{\mathrm{02,2}}^2\\ {\mathbf{v}}_{\mathrm{03,3}}^3\\ {\mathbf{v}}_{\mathrm{04,4}}^4\\ {\mathbf{v}}_{\mathrm{05,5}}^5\\ {\mathbf{v}}_{\mathrm{06,6}}^6\end{array}\right]=\left[\begin{array}{cccccc}{\mathbf{h}}_{\mathrm{1,1}}^1& 0& 0& 0& 0& 0\\ 0& {\mathbf{h}}_{\mathrm{2,2}}^2& 0& 0& 0& 0\\ 0& 0& {\mathbf{h}}_{\mathrm{3,3}}^3& 0& 0& 0\\ 0& 0& 0& {\mathbf{h}}_{\mathrm{4,4}}^4& 0& 0\\ 0& 0& 0& 0& {\mathbf{h}}_{\mathrm{5,5}}^5& 0\\ 0& 0& 0& 0& 0& {\mathbf{h}}_{\mathrm{6,6}}^6\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{q}}}_1^{}\\ {\dot{\mathbf{q}}}_2^{}\\ {\dot{\mathbf{q}}}_3^{}\\ {\dot{\mathbf{q}}}_4^{}\\ {\dot{\mathbf{q}}}_5^{}\\ {\dot{\mathbf{q}}}_6^{}\end{array}\right].$ (18)

The equation can be written in a compact form as follows:

$T_a^{\ell \ell }{\mathbf{v}}_a^{\ell }={\mathbf{H}}_a^{\ell }{\dot{\mathbf{q}}}_a={\mathbf{v}}_r^{\ell },$ (19)

where T _a ^ℓℓ is the topological assembly matrix and it could be easily constructed using the connectivity graph and the parent child list, the superscript ℓ refers to local matrices, the superscript a refers to the system assembly matrix, and ${\mathbf{H}}_a^{\ell }$ is the influence coefficient matrices grouped by joint block [6, 21]. Similarly, the acceleration could be calculated as follows:

$T_a^{\ell \ell }{\mathbf{a}}_a^{\ell }={\mathbf{H}}_a^{\ell }{\ddot{\mathbf{q}}}_a+{\dot{\mathbf{H}}}_a^{\ell }{\dot{\mathbf{q}}}_a={\mathbf{H}}_a^{\ell }{\ddot{\mathbf{q}}}_a+{\tilde{\mathbf{v}}}_r^{\ell }{\mathbf{v}}_r^{\ell }={\mathbf{H}}_a^{\ell }{\ddot{\mathbf{q}}}_a+{\mathbf{\gamma }}_a^{\ell }.$ (20)

Knowing the velocity and accelerations of the root body, (19) and (20) could be used to recursively calculate the velocity and acceleration of the descendent bodies marching outward from the root to the leaves. Similarly, the system topological matrix could be used to extract the joint force as follows:

$Q_a={\left(H_a^l\right)}^T{\mathbf{F}}_a^l.$ (21)

3.2. Dynamic Equations of Motion

The equation of motion of any link/body could be obtained by differentiating the vector of spatial momentum with respect to time. The spatial inertia matrix ${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^f{M}_{j,GG}^{jj}$ of body j that is defined at the body's center of mass marker G, can be written as follows:

${}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^f{M}_{j,GG}^{jj}=\left[\begin{array}{cc}{\underset{\_}{J}}_{j,GG}^{jj}& \mathrm{0}\\ \mathrm{0}& m_j\underset{\_}{I}\end{array}\right],$ (22)

where m_j is the body mass, $\underset{\_}{I}$ is 3 × 3 identity matrix, and ${\underset{\_}{J}}_{j,GG}^{jj}$ is the 3 × 3 matrix representing the body j moment of inertia defined at a triad located at the body center of mass. Note that the inertia matrix is defined in the function space. The spatial momentum vector ${\mathbf{P}}_{j,G}^j$ could be calculated as follows:

${\mathbf{P}}_{j,G}^j=M_{j,GG}^{jj}{\mathbf{v}}_{0j,G}^j.$ (23)

In order to obtain the equation of motion, the spatial momentum should be expressed in the global coordinate system as follows:

$\begin{array}{l}{\mathbf{P}}_{j,0}^0=M_{j,00}^{00}{\mathbf{v}}_{0j,0}^0,\\ \text{where}\mathrm{\hspace{0.17em}\hspace{0.17em}}{M}_{j,\mathrm{00}}^{\mathrm{00}}=X_{\mathrm{0}j}^{\mathrm{0}j-T}{M}_{j,jj}^{jj}{X}_{\mathrm{0}j}^{\mathrm{0}j-\mathrm{1}}=X_{j\mathrm{0}}^{j\mathrm{0}T}{M}_{j,jj}^{jj}{X}_{j\mathrm{0}}^{j\mathrm{0}}.\end{array}$ (24)

Differentiating the momentum equation will lead to the body equation of motion as follows:

$M_{j,\mathrm{00}}^{\mathrm{00}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathbf{a}}_{\mathrm{0}j,\mathrm{0}}^{\mathrm{0}}-{\tilde{\mathbf{v}}}_{0j,0}^{0T}{M}_{j,00}^{00}{\mathbf{v}}_{0j,0}^0={\mathbf{G}}_{j,0}^0,$ (25)

where ${\mathbf{G}}_{j,\mathrm{0}}^{\mathrm{0}}={\left[\begin{array}{cc}{\underset{\_}{\tau }}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}j,\mathrm{0}}^{\mathrm{0}}& {\underset{\_}{g}}_j^{\mathrm{0}}\end{array}\right]}^T$ is the spatial force vector including the gravitational forces, ${\underset{\_}{\tau }}_{j,\mathrm{0}}^{\mathrm{0}}$ is the sum of all external torques acting on the body, and ${\underset{\_}{g}}_j^{\mathrm{0}}$ is the sum of all external forces acting on the body.

The equation of motion of any two connected bodies i and j could be driven from the free body diagram as follows. For body j, the equation of motion in its local coordinate system could be written as

$M_{j,jj}^{jj}{\mathbf{a}}_{\mathrm{0}j,j}^j-{\mathbf{F}}_{j,j}^j={\mathbf{G}}_{j,j}^j+{\tilde{\mathbf{v}}}_{0j,j}^{jT}{\mathbf{P}}_{j,j}^j.$ (26)

While for body i, the equation of motion in its local coordinate system could be written as

$M_{i,ii}^{ii}{\mathbf{a}}_{\mathrm{0}i,i}^i-{\mathbf{F}}_{i,i}^i+{}_{}{}^f{X}_{ji}^{jiT}{\mathbf{F}}_{j,j}^j={\mathbf{G}}_{i,i}^i+{\tilde{\mathbf{v}}}_{0i,i}^{iT}{\mathbf{P}}_{i,i}^i.$ (27)

The equations could be rearranged in matrix form as follows:

$\left[\begin{array}{cc}{M}_{i,ii}^{ii}& \mathrm{0}\\ \mathrm{0}& M_{j,jj}^{jj}\end{array}\right]\left[\begin{array}{c}{\mathbf{a}}_{\mathrm{0}i,i}^i\\ {\mathbf{a}}_{\mathrm{0}j,j}^j\end{array}\right]+\left[\begin{array}{cc}I& -{}_{\mathrm{\hspace{0.17em}\hspace{0.17em}}}{}^f{X}_{ji}^{jiT}\\ \mathrm{0}& I\end{array}\right]\left[\begin{array}{c}-{\mathbf{F}}_{i,i}^i\\ -{\mathbf{F}}_{j,j}^j\end{array}\right]=\left[\begin{array}{c}{\mathbf{G}}_{i,i}^i+{\tilde{\mathbf{v}}}_{0i,i}^{iT}{P}_{i,i}^i\\ {\mathbf{G}}_{j,j}^j+{\tilde{\mathbf{v}}}_{0j,j}^{jT}{\mathbf{P}}_{j,j}^j\end{array}\right].$ (28)

The system equation of motion could be written in a more compact form as follows:

$M_{}^{ll}{\mathbf{a}}_a^l-{\left(T_a^{ll}\right)}^T{\mathbf{F}}_a^l={\mathbf{G}}_a^l+{\left({\tilde{\mathbf{v}}}_a^l\right)}^T\mathbf{P}.$ (29)

The system of equations in (29) contains the unknown Cartesian acceleration of the bodies as well as the joint reaction forces. In order to be able to solve the equation of motion, the kinematic relations from (20), the constraint equations from (21), and the system dynamic equation (29) are rearranged in matrix form as follows:

$\left[\begin{array}{ccc}{M}_{}^{\ell \ell }& 0& {\left(T_a^{\ell \ell }\right)}^T\\ 0& 0& -{\mathbf{H}}_a^{\ell }\\ T_a^{\ell \ell }& -{\left({\mathbf{H}}_a^{\ell }\right)}^T& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{a}}_a^{\ell }\\ {\ddot{\mathbf{q}}}_a\\ -{\mathbf{F}}_a^{\ell }\end{array}\right]=\left[\begin{array}{c}{\mathbf{G}}_a^{\ell }+{\left({\tilde{\mathbf{v}}}_a^{\ell }\right)}^T{\mathbf{P}}_a^{\ell }\\ Q_a\\ {\mathbf{\gamma }}_a^{\ell }\end{array}\right].$ (30)

Since the inertia matrix M^ℓℓ is invertible, (30) could be solved to get the unknown joint accelerations ${\ddot{\mathbf{q}}}_a$ which can used to calculate the joint Cartesian accelerations ${\mathbf{a}}_a^{\ell }$ , the joint reaction forces ${\mathbf{F}}_a^{\ell }$ as follows:

$\left({\left({\mathbf{H}}_{\mathbf{a}}^{\mathcal{l}}\right)}^T{\overline{M}}_a{\mathbf{H}}_a^{\ell }\right){\ddot{\mathbf{q}}}_a=Q_a+{\left({\mathbf{H}}_a^{\ell }\right)}^T{\left(T_a^{\ell \ell }\right)}^{-T}{\overline{G}}_a.$ (31)

Equation (31) has been simplified to understand the significance of each term as follows: ${\overline{M}}_a={(T_a^{\ell \ell })}^{-T}{M}^{\ell \ell }{(T_a^{\ell \ell })}^{-\mathrm{1}}$ represents the system apparent inertia matrix. This term could be calculated by backward marching through the kinematic tree from the children to the root and adding the inertia matrix of the decedent bodies into their parent inertia matrix. The term ${\overline{G}}_a=({\mathbf{G}}_a^{\ell }+{({\tilde{\mathbf{v}}}_a^{\ell })}^T{\mathbf{P}}_a^{\ell })-M^{\ell \ell }{(T_a^{\ell \ell })}^{-\mathrm{1}}{\mathbf{\gamma }}_a^{\ell }$ accumulates the quadratic velocity term to the Coriolis forces for each body and adds it to the vector of external forces. Multiplying these accumulated forces by ${(T_a^{\ell \ell })}^{-T}$ implies backward marching from the children to the root and adding this force vector to the parent body. The accumulated forces projected them into the joint subspace through multiplying by ${({\mathbf{H}}_a^{\ell })}^T{(T_a^{\ell \ell })}^{-T}{\overline{G}}_a$ . The double product ${({\mathbf{H}}_a^{\ell })}^T{\overline{M}}_a{\mathbf{H}}_a^{\ell }$ projects the accumulated apparent inertia matrix of each body into the body joint subspace. The equation of motion could be written in a simple form as follows:

$\widehat{M_q}{\ddot{\mathbf{q}}}_a=\widehat{Q_a},$ (32)

where $\widehat{M_q}={({\mathbf{H}}_a^{\ell })}^T{\overline{M}}_a{\mathbf{H}}_a^{\ell }$ in the projection of the apparent inertia matrix into the joint coordinate subspace and $\widehat{Q_a}=Q_a+{({\mathbf{H}}_a^{\ell })}^T{(T_a^{\ell \ell })}^{-T}{\overline{G}}_a$ is the projection of the cartesian forces into the joint coordinate subspace. The inertia matrix in (32) could be inverted and the joint accelerations could be computed. Featherstone [6] demonstrated that the number of projection across the kinematic tree is highly dependent on the number of joint degrees of freedom between the child and root body, n_list. The number of operation is in the order of ${(n_{\text{list}})}^3$ . For this reason, the depth of the kinematic tree should be minimized in order to improve the computational efficiency.

3.3. Floating Body Equation of Motion

In order to reduce the problem dimensionality and develop a general purpose independent algorithm, the chain links will be decoupled from the kinematic tree of the main system. All the chain links have neither parents nor children and can be considered as floating bodies. This means they will not project any inertia into their parents subspace and their inertia will not be affected by any decedent bodies. FB representing chain links will be interacting with other system components through compliant forces: contact with sprocket, contact with guide rails, and connection to other links. These contact forces can be calculated in the global coordinate system, as will be shown in the following sections, and then will be transformed into the local link reference CS. The equation of motion of FB could be written in local CS in order to improve the computational efficiency. This approach will significantly improve the computational efficiency by reducing the depth of the kinematic tree. Since all external force vectors are calculated and accumulated in the global CS, it can be then transformed once to the FB local CS. Also, there will be no need to transform the FB mass matrix to the inertial CS each time step.

The equation of motion of FB in local CS can be written as follows:

$M_{j,jj}^{jj}{\mathbf{a}}_{0j,j}^j={\mathbf{G}}_{j,j}^j+{\tilde{\mathbf{v}}}_{0j,j}^{jT}{M}_{j,jj}^{jj}{\mathbf{v}}_{0j,j}^j.$ (33)

The link mass matrix can be easily inverted and the link accelerations could be directly integrated as follows:

${\mathbf{a}}_{\mathrm{0}j,j}^j={\left(M_{j,jj}^{jj}\right)}^{-\mathrm{1}}\left({\mathbf{G}}_{j,j}^j+{\tilde{\mathbf{v}}}_{0j,j}^{jT}{M}_{j,jj}^{jj}{\mathbf{v}}_{0j,j}^j\right).$ (34)

Since each FB is referenced with 6-DOF, the rotational transformation matrix need to be updated using the body rotational parameters. In the proposed approach, Euler parameters are used describe the orientation of the FB and that requires using one additional constraint equation for each body. It was reported by many researchers [7, 15, 20, 22] that in the case of long chain or track systems, the Euler parameters norm in the integrator may drift from unity causing phase shift. In order to improve the computational efficiency, the constraint equation could be eliminated and replaced with error penalty feedback within the integrator [23, 24]. Assume that Euler parameters are defined as $E={\left[\begin{array}{cccc}{e}_{\mathrm{0}}& e_{\mathrm{1}}& e_{\mathrm{2}}& e_{\mathrm{3}}\end{array}\right]}^T$ and $E_s={\left[\begin{array}{ccc}{e}_{\mathrm{1}}& e_{\mathrm{2}}& e_{\mathrm{3}}\end{array}\right]}^T$ . The normalized Euler parameters can be calculated as follows:

${\underset{\_}{E}}_s=\frac{E_s}{\sqrt{E_s^T{E}_s}}.$ (35)

The error in calculating Euler parameters can be approximated as:

$E_{\mathrm{err}}=\mathrm{1}-\left|{\underset{\_}{E}}_s\right|.$ (36)

This error may be fed back into the Euler parameter differential equations as

${\underset{\_}{\dot{E}}}_s⟸{\underset{\_}{\dot{E}}}_s-\mathrm{2}{E}_{\mathrm{err}}\hspace{0.17em}{\underset{\_}{E}}_s.$ (37)

Equation (37) shows that a small negative correction, proportional to error and the Euler parameters, is fed back to the differential equations. Implementing the feedback term into (37) will normalize the Euler parameters without introducing constraint equations to the system.

3.4. System Interaction

In order to complete the equation of motion of the chain drive system within a multibody system, the interaction between the chain links and the other components in the system should be modeled. In the following sections, the contact between the chain links and the sprockets will be presented. The compliant bushing element that will be used to model the connection between two consecutive links will be discussed. It should be noted that the roller will be ignored in this analysis.

3.4.1. Bushing Force

In the proposed approach, a compliant bushing force element is used to model the connection between two consecutive links. The chain is considered to be composed of two repeated links, the pin link and the bushing link as shown in Figure 5. Each pin link has two triads, O_pi1 and O_pi2, to locate the center of the two pins with respect to the pin link reference frame, O _pi . Similarly, each bushing link has two triads, O_bi1 and O_bi2, to locate the center of the two bushing centers with respect to the bushing link reference frame, O _bi . Each pin triad of the pin link is connected to the bushing triad of the bushing link via a bushing force element. The pin and bushing triads connected by the bushing element are considered to be coincident and parallel if there is no deflection or relative rotation between links, as shown in Figure 5. In order to calculate the bushing forces, the translational deflection and relative orientation between the pin and bushing triad should be calculated and updated every time step.

Figure 5:

Kinematics of the compliant bushing force element.

During the dynamic simulation, the position and velocity vectors of each body's reference frame are calculated. Also, the pin and bushing global position vectors, ${\underset{\_}{r}}_{\mathrm{0}p}^{\mathrm{0}}$ and ${\underset{\_}{r}}_{\mathrm{0}b}^{\mathrm{0}}$ , and the global transformation matrices, ${\underset{\_}{R}}^{\mathrm{0}p}$ and ${\underset{\_}{R}}^{\mathrm{0}b}$ , can be calculated using (3) and (7). The relative displacement vector between the two triads and the relative orientation matrix could be calculated in the pin coordinate system as follows.

$\begin{array}{l}{\underset{\_}{r}}_{pb}^{\mathrm{0}}={\underset{\_}{r}}_{\mathrm{0}b}^{\mathrm{0}}-{\underset{\_}{r}}_{\mathrm{0}p}^{\mathrm{0}},\hspace{1em}\hspace{1em}{\underset{\_}{r}}_{pb}^p={\left({\underset{\_}{R}}^{\mathrm{0}p}\right)}^{-\mathrm{1}}{\underset{\_}{r}}_{pb}^{\mathrm{0}},\\ {\underset{\_}{R}}^{pb}={\left({\underset{\_}{R}}^{\mathrm{0}p}\right)}^{-\mathrm{1}}{\underset{\_}{R}}^{\mathrm{0}b}.\end{array}$ (38)

Once transformed into the pin coordinate system, the bushing deformation is represented by the relative translational displacement vector ${\underset{\_}{r}}_{pb}^p$ and the relative orientation matrix ${\underset{\_}{R}}^{pb}$ . The bushing rotational deformation can be calculated in two steps from orientation matrix as follows. Using the relative orientation matrix ${\underset{\_}{R}}^{pb}$ Euler parameters can be calculated [25] by expressing the transformation matrix in terms of Euler parameters as follows:

${\underset{\_}{R}}^{pb}=\underset{\_}{I}+\mathrm{2}{\widetilde{\underset{\_}{E}}}_s\left(e_{\mathrm{0}}^{}\underset{\_}{I}+{\widetilde{\underset{\_}{E}}}_s\right),\hspace{1em}\hspace{1em}{\underset{\_}{E}}_s={\left[\begin{array}{ccc}{e}_{\mathrm{1}}& e_{\mathrm{2}}& e_{\mathrm{3}}\end{array}\right]}^T.$ (39)

Then Euler parameters can be converted into a rotational angle θ around a vector ${\underset{\_}{u}}^r$ as follows [25]:

$\begin{array}{c}{e}_{\mathrm{0}}=\cos \left(\frac{\theta }{2}\right),\hspace{1em}\hspace{1em}{e}_{\mathrm{1}}=u_{\mathrm{1}}^r\sin \left(\frac{\theta }{2}\right),\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}\\ e_{\mathrm{2}}=u_{\mathrm{2}}^r\sin \left(\frac{\theta }{2}\right),\hspace{1em}\hspace{1em}{e}_{\mathrm{3}}=u_{\mathrm{3}}^r\sin \left(\frac{\theta }{2}\right).\end{array}$ (40)

The relative velocity vector between the pin and bushing triads can be calculated and transformed into the pin coordinate frame as follows:

$\begin{array}{c}{\mathbf{v}}_{0b,p}^0=\left[\begin{array}{cc}I& \mathrm{0}\\ {\underset{\_}{\tilde{r}}}_{\mathrm{0}p}^{\mathrm{0}}-{\underset{\_}{\tilde{r}}}_{\mathrm{0}b}^{\mathrm{0}}& I\end{array}\right]\left[\begin{array}{c}{\omega }_{\mathrm{0}b}^{\mathrm{0}}\\ {\underset{\_}{\dot{t}}}_{\mathrm{0}b,b}^{\mathrm{0}}\end{array}\right],\\ {\mathbf{v}}_{pb,p}^p={\left(X_{0p}^{0p}\right)}^{-1}\left({\mathbf{v}}_{0b,p}^0-{\mathbf{v}}_{0p,p}^0\right),\\ {\mathbf{v}}_{pb,p}^p=\left[\begin{array}{c}{\omega }_{pb}^p\\ {\underset{\_}{\dot{r}}}_{pb,p}^p\end{array}\right]=\left[\begin{array}{cc}I& \mathrm{0}\\ {\underset{\_}{\tilde{r}}}_{pb}^p& I\end{array}\right]\left[\begin{array}{c}{\omega }_{pb}^p\\ {\underset{\_}{\dot{r}}}_{pb,p}^p\end{array}\right].\end{array}$ (41)

It is assumed that the bushing element produces three force components F _x , F _y , and F _z that are functions of the relative displacements and velocities between the pin and bushing triads and two torque components T _x and T _y that are functions of the relative angular rotations and angular velocities between the triads. The torque component T _z is assumed to be zero in order to allow relative rotation. The bushing forces and torques can be computed in the pin triad coordinate system as follows:

$F_{p,p}^p=\left[\begin{array}{c}-\left(\theta K_{tx}\right)u_{\mathrm{1}}^r-{\omega }_x{C}_{tx}\\ -\left(\theta K_{ty}\right)u_{\mathrm{2}}^r-{\omega }_y{C}_{ty}\\ -\left(\theta K_{tz}\right)u_{\mathrm{3}}^r-{\omega }_z{C}_{tz}\\ -t_x{K}_x-{\dot{t}}_x{C}_x\\ -t_y{K}_y-{\dot{t}}_y{C}_y\\ -t_z{K}_z-{\dot{t}}_z{C}_z\end{array}\right],$ (42)

where K _tx , K _ty , K _tz , C _tx , C _ty , and C _tz represent the rotational stiffness and damping coefficients K _x , K _y , K _z , C _x , C _y , and C _z represent the translational stiffness and damping coefficients.

The forces and torques could now be transformed into the bushing triad coordinate system and applied to the bushing link body using (4). The forces from the bushing force element could be then applied to the link equation of motion (24). The bushing element stiffness coefficients can be calculated experimentally or by using a finite element model for the pin-link and bushing link while the damping coefficient could be chosen to represent the structural damping.

3.4.2. Link Sprocket Contact

In chain drive systems, motion is transmitted between sprockets through contact between the chain bushing/rollers and the sprocket teeth. In the case of planar chains dynamic simulation, the code accounts for the planar contact between the surface of the sprocket teeth and the bushing cylinder; however contact may occur between the sprocket and the side plates in bicycle simulation. During speed shift and when the rear sprockets are not in the same plane of the drive sprocket, side contact may occur. The sprocket tooth could be described by detailed geometric features, x-y spline curve, or approximated using tangent lines. The first two approaches are computationally intensive while the last approach leads to reasonably accurate results faster [20].

In this analysis, the tooth profile is approximated using six points defined in the x-y plane relative to the sprocket center, as shown in Figure 6. The tooth profile will be replicated around z-axis according to the number of sprocket teeth. The tooth thickness could be specified to fully define the sprocket geometry and to detect the side contact with the link plates. Five possible contact surfaces could be defined using the profile points and the sprocket thickness. Prior to simulation, the profile tangent and normal vectors could be determined for each surface segment and saved for repeated calculations later in order to reduce the simulation run time.

Figure 6:

Approximated sprocket tooth profile.

The contact detection algorithm loops over all the sprockets in the system to identify the links that are within the prespecified proximity range of each sprocket. A typical proximity range would be the sprocket outer radius plus one link pitch. Referring to Figure 7, the relative position vector of the chain link reference frame with respect to the sprocket reference could be defined in the sprocket frame as follows:

$\begin{array}{l}{\underset{\_}{r}}_{sl}^{\mathrm{0}}={\underset{\_}{r}}_{\mathrm{0}l}^{\mathrm{0}}-{\underset{\_}{r}}_{\mathrm{0}s}^{\mathrm{0}},\hspace{1em}\hspace{1em}{\underset{\_}{t}}_{sl}^s={\left({\underset{\_}{R}}^{\mathrm{0}s}\right)}^{-\mathrm{1}}{\underset{\_}{t}}_{sl}^{\mathrm{0}},\\ {\underset{\_}{R}}^{sl}={\left({\underset{\_}{R}}^{\mathrm{0}s}\right)}^{-\mathrm{1}}{\underset{\_}{R}}^{\mathrm{0}l},\end{array}$ (43)

where ${\underset{\_}{r}}_{\mathrm{0}s}^{\mathrm{0}}$ and ${\underset{\_}{r}}_{\mathrm{0}l}^{\mathrm{0}}$ are the global position vector of the sprocket and the link reference frames, ${\underset{\_}{r}}_{sl}^{\mathrm{0}}$ and ${\underset{\_}{r}}_{sl}^s={\left[\begin{array}{ccc}{r}_x& r_y& r_z\end{array}\right]}^T$ are the relative position vector in the global and in the sprocket frame, ${\underset{\_}{R}}^{\mathrm{0}s}$ and ${\underset{\_}{R}}^{\mathrm{0}l}$ are the global orientation matrix of the sprocket and bushing frames, ${\underset{\_}{R}}^{sl}$ is the relative orientation matrix between the sprocket and the link. The distance d between the sprocket center and the link reference frame is computed as

$d=\sqrt{\left(r_x^{\mathrm{2}}+r_y^{\mathrm{2}}\right)}.$ (44)

If the distance d is less than the prespecified proximity, δ_prox, then contact might occur between this link and the sprocket teeth profile. In order to exactly identify the contact point, the position vector of the bushing center is calculated and transformed into the sprocket plane as follows:

$\begin{array}{l}{\underset{\_}{r}}_{sb}^{\mathrm{0}}={\underset{\_}{r}}_{\mathrm{0}b}^{\mathrm{0}}-{\underset{\_}{r}}_{\mathrm{0}s}^{\mathrm{0}},\hspace{1em}\hspace{1em}{\underset{\_}{r}}_{sb}^s={\left({\underset{\_}{R}}^{\mathrm{0}s}\right)}^{-\mathrm{1}}{\underset{\_}{r}}_{sb}^{\mathrm{0}},\\ {\underset{\_}{R}}^{sb}={\left({\underset{\_}{R}}^{\mathrm{0}s}\right)}^{-\mathrm{1}}{\underset{\_}{R}}^{\mathrm{0}b},\end{array}$ (45)

where ${\underset{\_}{r}}_{\mathrm{0}s}^{\mathrm{0}}$ and ${\underset{\_}{r}}_{\mathrm{0}b}^{\mathrm{0}}$ are the global position vector of the sprocket and the bushing triad, ${\underset{\_}{r}}_{sb}^{\mathrm{0}}$ and ${\underset{\_}{r}}_{sb}^s={\left[\begin{array}{ccc}{r}_x& r_y& r_z\end{array}\right]}^T$ are the relative position vector in global and in the sprocket frame, ${\underset{\_}{R}}^{\mathrm{0}s}$ and ${\underset{\_}{R}}^{\mathrm{0}b}$ are the global orientation matrix of the sprocket and bushing frames, and ${\underset{\_}{R}}^{sb}$ is the relative orientation matrix between the sprocket and the bushing. The components of the relative position vector could be used to determine the angle ψ between the relative position vector and the vertical y-axis of the sprocket reference frame and hence the contact tooth index could be identified. The relative position between the profile points to bushing center is determined and lines L_1–6 are defined, as shown in Figure 7. Contact can occur between the bushing circle and the profile line segment if the projection of the tangent vector passes through the center of the circle or if either segment end-point falls inside the circle. The normal distance between the bushing center and the contact profile segment could be easily determined and compared to the bushing radius, R _b , to determine if penetration occurred. The position vector of the exact contact point, p_c, relative to the sprocket reference frame could be determined as ${\underset{\_}{r}}_{sp}^s$ , and the angle φ could be calculated from the vector coordinates.

Figure 7:

Bushing sprocket detection.

The global velocities of the bushing link and sprocket, ${\mathbf{v}}_{0b,b}^0$ and ${\mathbf{v}}_{0s,s}^0$ , are computed every time step. The relative velocity between bushing link and sprocket can be computed as follows. The bushing velocity is computed at the sprocket reference frame:

${\mathbf{v}}_{0b,s}^0=\left[\begin{array}{cc}I& \mathrm{0}\\ {\underset{\_}{\tilde{r}}}_{\mathrm{0}b}^{\mathrm{0}}-{\underset{\_}{\tilde{r}}}_{\mathrm{0}s}^{\mathrm{0}}& I\end{array}\right]\left[\begin{array}{c}{\omega }_{\mathrm{0}b}^{\mathrm{0}}\\ {\underset{\_}{\dot{r}}}_{\mathrm{0}b,b}^{\mathrm{0}}\end{array}\right].$ (46)

The relative velocity can be computed in global at the sprocket reference frame as follows:

${\mathbf{v}}_{sb,s}^s={\left(X_{0s}^{0s}\right)}^{-1}\left({\mathbf{v}}_{0b,s}^0-{\mathbf{v}}_{0s,s}^0\right).$ (47)

The relative velocity between the bushing and the sprocket at the contact point could be computed as follows:

${\mathbf{v}}_{sb,p}^s=\left[\begin{array}{c}{\omega }_{sb}^s\\ {\underset{\_}{\dot{r}}}_{sb,s}^s\end{array}\right]=\left[\begin{array}{cc}I& \mathrm{0}\\ {\underset{\_}{\tilde{r}}}_{sp}^s& I\end{array}\right]\left[\begin{array}{c}{\omega }_{sb}^s\\ {\underset{\_}{\dot{r}}}_{sb,s}^s\end{array}\right].$ (48)

In order to calculate the normal and sliding velocity components, ${\underset{\_}{\dot{r}}}_{sb,s}^s={\left[\begin{array}{ccc}{v}_{rx}& v_{ry}& v_{rz}\end{array}\right]}^T$ needs to be transformed using the φ as follows:

$\left[\begin{array}{c}{v}_{x-\text{slid}}\\ v_{\mathrm{norm}}\\ v_{z-\text{slid}}\end{array}\right]=\left[\begin{array}{ccc}\cos \phi & \sin \phi & \mathrm{0}\\ -\sin \phi & \cos \phi & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{c}{v}_{rx}\\ v_{ry}\\ v_{rz}\end{array}\right].$ (49)

The sliding velocity could be computed as

$\begin{array}{c}{\underset{\_}{v}}_s=\left|{\underset{\_}{v}}_s\right|{\underset{\_}{u}}_{vs},\hspace{1em}\text{where}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left|{\underset{\_}{v}}_s\right|=\sqrt{v_{x-\text{slid}}^{\mathrm{2}}+v_{z-\text{slid}}^{\mathrm{2}}},\\ {\underset{\_}{u}}_{vs}=\frac{{\underset{\_}{v}}_s}{\left|{\underset{\_}{v}}_s\right|}.\end{array}$ (50)

The normal contact force results from elastic spring force and damping force can be evaluated as follows:

$F_n=K\delta +C{v}_{\mathrm{norm}}.$ (51)

The sliding friction force is a function of the normal elastic force and the friction coefficient. Due to impact between the bushing and the sprocket, the normal force component may be significantly high while the sliding velocities are significantly small. The large normal force can lead to very large friction forces causing the simulation to fail. To reduce the effect of large friction force, the friction coefficient can be considered as smooth function of the sliding velocity as follows:

$\begin{array}{l}{\underset{\_}{F}}_f=\mu S_{\mu }\left(K\delta \right){\underset{\_}{u}}_{vs},\\ \text{where}\mathrm{\hspace{0.17em}\hspace{0.17em}}{S}_{\mu }=\min \left(\mathrm{1,2}\xi -{\xi }^{\mathrm{2}}\right),\hspace{1em}\hspace{1em}\xi =\frac{\left|{\underset{\_}{v}}_s\right|}{v_{\mu }},\end{array}$ (52)

where v_μ is a predefined maximum sliding velocity at which the friction force reaches maximum value. The friction force softening parameter S_μ is shown in Figure 8. Other softening functions could be effectively used [25–27]. The contact and friction forces vectors are rotated back to the sprocket reference frame as follows:

${\underset{\_}{F}}^c=\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]=\left[\begin{array}{ccc}\cos \phi & -\sin \phi & \mathrm{0}\\ \sin \phi & \cos \phi & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{c}{F}_{fx}\\ F_{\mathrm{norm}}\\ F_{fz}\end{array}\right].$ (53)

Then contact forces are now translated to the sprocket reference frame and the link reference frame as follows:

$F_{s,s}^s={\left[\begin{array}{cc}I& \mathrm{0}\\ {\underset{\_}{\tilde{r}}}_{sp}^s& I\end{array}\right]}^{-T}\left[\begin{array}{c}\mathrm{0}\\ {\underset{\_}{F}}^c\end{array}\right],\hspace{1em}\hspace{1em}{F}_{l,l}^l={\left(X_{sl}^{sl}\right)}^{-T}{F}_{s,s}^s.$ (54)

Similarly the contact between the sprocket tooth and the bushing link side plates could be detected as shown in Figure 9. The tooth profile points are projected into the bushing plate surface plane. The angle α could be calculated from the relative orientation matrix between the bushing link and the sprocket, ${\underset{\_}{R}}^{sb}$ . The normal distance between the profile point and the project point is calculated, δ _s . Contact can occur if δ_s-pen = δ _s – δ _b <0. The relative velocity at the contact point could be calculated using (49) to (51). The relative velocity and the penetration δ_s-pen could be used to calculate the contact and friction forces between the sprocket tooth and bushing plates, (52) to (54).

Figure 8:

The friction force smoothing parameter, S_μ.

Figure 9:

Contact between the sprocket tooth and bushing link side plates.

It should be noted that contact can occur between the sprocket tooth and both the bushing cylinder and the side plate, as shown in Figure 9. The contact algorithm should be able to detect both contacts and add the resulting forces and moments. In reality, the sprocket teeth are usually chamfered to reduce the side contact especially in applications that may undergo out-of-plane motion like bicycle chains. It should be noted that the bushing cylinder is assumed to be rigid and will not deform under the effect of the contact forces.

4. Recursive Solution Procedure

The proposed approach was developed using general purpose joint-based recursive approaches [17, 18, 20, 25–29]. Prior to solving the equation of motion, the system must be assembled and the chain must be wrapped around the sprocket without penetrating into the sprocket. The initial conditions are determined to minimize the constraint violation and eliminate any interference between the metals especially contact between links and sprocket. A static equilibrium could be performed to obtain excellent initial conditions that could be fed into the dynamic simulation. At the beginning of the dynamic simulation, the structure of the equation of motion is determined based on a preliminary analysis of the system topology. The dependent and independent variable sets are determined using the generalized coordinate partitioning approach. The solver integrates both dependent and independent variable sets and the kinematic constraints are enforced. The input states to the integrator represent the first and second time derivative of the joint variables of the bodies in the main system and the FB. The output states are the joint displacements and the joint velocities of the bodies in the main system and the FB. The following summarizes the structure of multibody dynamic simulation recursive algorithm including the chain drives.

1. Using the joint variables returned from the integrator, the Cartesian displacements, velocities, and accelerations are calculated. Forward evaluation scheme is utilized (starting from the root body to the descendent/branch bodies).

2. The Cartesian displacements, velocities, and accelerations of the chain links are calculated.

3. The generalized coordinate portioning technique is applied to determine the sets of independent and dependent variables.

4. The internal and external forces are calculated and applied to the different bodies in the main system.

5. Process of the chain links:

6. Transform of the inertia matrices into global coordinate system.

7. Calculation of the inertia forces, centrifugal and Coriolis forces in Cartesian space.

8. For all the chain links, the inertia matrix, the centrifugal and Coriolis forces, and the external forces are transformed into the global local coordinate system.

9. Propagation of the external and inertia forces from the descendent bodies into their parents.

10. Projection of the Cartesian forces into the body joint space.

11. Factoring the mass matrix based on the independent and dependent variables sets.

12. Calculation of second derivatives of the joint variables.

13. Inversion of the chain links (floating body) mass matrix and calculation of the body second derivatives.

14. Sending the states to the integrator.

15. Appending the chain links states to the main system states and sending it to the integrator.

This algorithm utilizes a predictor-corrector explicit integrator with variable-order interpolation polynomials and variable time step [28, 30]. This explicit integrator insures the stability of the solution and the ability to capture the high-speed impacts between the different machine parts.

5. Numerical Example

The following example demonstrates the results of the proposed approach. The example represents a bicycle chain with 56 pin links and 56 roller links. The bike model is comprised of front drive sprocket (DS), wheel sprocket (WS), the lower idler sprocket (LIS), and the upper idler sprocket (UIS). The connectivity graph of the system is shown in Figure 10. To demonstrate the chain dynamics, the bicycle chassis is assumed to be fixed at the ground inertial frame of reference. The DS, WS, and the derailleur are attached to the chassis with revolute joints. The lower and upper idler sprockets are attached to the derailleur with revolute joints. All the chain links are referenced to the ground with full 6-DOF joints (floating bodies). The details of the sprockets design and locations are listed in Table 1. It should be noted that the DS attachment position is offset from the WS attachment position by 30 mm in the lateral z-direction as shown in Figure 2. The chain link pitch is 12.7 mm, the bushing radius is 4 mm, and the pin radius is 2 mm. The sprocket tooth thickness is 3.2 mm.

Table 1:

Characteristics of the bicycle chain sprockets system.

	N-Teeth	Pitch radius	Outer radius	X	Y	Z
DS	23	46.64	50.01	0	0	0.03
WS	16	32.55	35.73	−486.2	0	0
UIS	10	20.55	23.35	−474.1	−88.9	0
LIS	10	20.55	23.35	−486.2	−177.8	0

Figure 10:

Connectivity graph of the bicycle chain and the drive system.

A kinematic driver is applied to the driving sprocket ramping up the speed from zero to 1 revolution per second in 3 seconds, as shown in Figure 11. A spline curve is used to smoothly ramp up the sprocket velocity and avoid high jerks. This resulting bicycle speed can be considered as slow speed for a beginner rider who can achieve approximately 10 km/hr on a 26″ bike.

Figure 11:

Speed of the drive sprocket.

In order for each link to complete one cycle, it will have to travel a distance of 1.422 m around the sprockets. After the drive sprocket reaches its maximum speed, the link forward speed will be averaged around 0.293 m/s. Each link will require about 5 seconds to complete one revolution. Figures 12 and 13 represent the acceleration and forces results of link 1, shown in Figure 2, as it moves two complete revolutions from its initial position. Figure 12(a) shows the link lateral acceleration as the sprocket speed increases from zero to maximum and after reaching a steady state. It could be seen from the Figure that the link vibration occurs in repeated pattern with time period of about 5 seconds. Figure 12(b) shows FFT of the link lateral acceleration after reaching the maximum speed. It can be seen that the dominant amplitudes peaks occur at multiples of 23 which is the number of teeth in the driving sprocket. The main reason for that can be explained by the lateral position shift between the front and rear sprocket. In the rare section of the mechanism, the three sprockets are in one plane and this configuration forces the chain links to move in the plane containing these sprockets. As the links pass this zone and approach the front driving sprocket, they engage with the driving sprocket teeth. The links then align themselves in the driving sprocket's plane. This action generates high bending moment around the vertical axis causing out-of-plane vibration. This vibration wave travels across the chain each time a link engages with the sprocket leading to the peaks that correspond to the number of front sprocket teeth.

Figure 12:

Lateral acceleration of link 1.

Figure 13:

Bushing forces in the chain as the driving sprocket speed is ramping up from zero to maximum.

The axial tension force in link 1 during the simulation is shown in Figure 13(a). It could be noted that during the startup period, the tension force is high. After the bicycle speed reaches the steady state, the tension force settles around lower amplitude and becomes periodic as shown in Figure 13. It should be noted that the tensile forces are always positive. The lateral shear force in the link bushing element is shown in Figure 13(b). It should be also noted the sign of the force changes depending on the link location along the chain as shown in Figure 13.

Figure 14 shows the bushing force vectors between the chain links at different instants of time after the sprocket reached the maximum speed. It should be noted that the dominant component of force vector is the tensile force along the chain length except where the chain is in contact with the sprockets.

Figure 14:

Changes in the bushing force vectors in the X-Y plan.

Figure 15 shows the bushing force vectors in space after the sprocket reached the maximum speed. It should be noted that the dominant components of force vectors do not change significantly with time. However, the minor changes in the force vector components are due to induced out-of-plane vibration. The out-of-plan vibration could be damped out by increasing the artificial damping in the bushing parameters.

Figure 15:

Changes in the bushing spatial force vectors.

6. Summary and Conclusions

This paper proposed a simple approach to include large chain dynamics into joint-based multibody system formulation. Each chain link is represented by a body with full 6-DOF to capture the exact link dynamics and allow out-of-plane motion. As the chain links are not kinematically coupled to the other bodies in the multibody system, it was noticed that there will be no inertia projections associated with their coordinates. For this reason, it was easy to separate those bodies from the system kinematic chain. The equation of motion of the chain links was formulated in the link coordinate system without including them in the system spanning tree. This approach resulted in significant reduction in the depth of the system spanning tree and problem dimensionality. The resulting mass matrix was significantly smaller. Euler parameters were used to represent the link orientation. A velocity feedback was used to eliminate the violation in the Euler parameters constraints without introducing algebraic constraint equations.

A complaint model was used to represent the connection between the adjacent chain links. The joint forces can be easily calculated without using Lagrange multipliers. The out-of-plane misalignment between the different sprockets can be easily modeled. Both in-plane and out-of-plane vibration of the chain can be easily captured. The contact between the chain links was modeled using a full geometric description of the sprocket tooth. The proposed approach can allow as many sprockets as the designer needs.

The proposed approach show significant numerical stability and computational efficiency. Through the application of the proposed approach, it was shown that the significant dynamics of the chain components easily captured and the contact problem is characterized. Ride and shift quality, contact forces, wear prediction, and noise emission are typical problems that can be easily investigated using the proposed modeling approach.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.