Continuous modeling of a multi-link flexible transmission

The problem of dynamic, infinite dimension, modeling of a transmission is considered. An accurate Laplace transfer function matrix of the system that consists of flexible shafts connected by gears that are either rigid or flexible is found. The first step is deriving a set of single input, infinite dimension, transfer functions for a single uniform link. The building blocks of those transfer functions are time delays, representing the wavemotion, and low order rational expressions, representing the boundary phenomena. The next step is combining these individual transfer functions into an overall model of the transmission, by means of the link reaction approach that makes use of the geometric relationships and reaction moments between neighboring links. The outcome is a generalized dynamic model with the moments in the gear pairs as the generalized state vector. The explicit and highly structured form of the transfer functions allows physical insight into the system, exact calculation of natural frequencies and the construction of exact simulation schemes built from standard blocks that are available in multi-purpose simulation software.


Introduction
Modeling of the angular motion of flexible shafts and gears transmissions, shown schematically in Fig. 1, is a well known problem. The dynamics of such systems is governed by both shaft and gear flexibility. A model for the gear flexibility, first suggested by Tuplin [1], consists of an infinitesimal gear mesh spring-damper element that is positioned in series with the gear transmission error excitation [2,3]. This model is widely used for gear modeling [2][3][4][5][6] and is adopted in this work as well.
An accurate model of the system requires modeling of the flexible shafts and their supports (boundary conditions) in addition to the gears flexibility. The dynamics of these elements is governed by the wave equation [4,7,8]. In [4] a continuous model for the gears and shafts, based on modal approach, was suggested. Since the system has infinitely many modes, the solution is given in terms of infinite series. In order to use such a model, approximations, such as truncating the high frequency modes, are necessary. Another method is using the Finite Element Method (FEM) to model the continuous system [8][9][10], with accuracy that depends upon the number of elements used. In [10] an accurate model for a uniform shaft with multiple lumped elements is presented, and used to numerically find the natural frequencies and mode shapes of the shaft.
The use of transfer functions to model infinite dimension systems is not a new idea, and many works used that approach, for example to model heat transfer systems [11]. The approach was also applied to flexible structures and results for some cases have been reported in [12][13][14][15][16][17]. The problem is often expressed in terms of Green's function of the system, as the transfer function is a Laplace transform of that function [14]. In some works, e.g. [15], the results are in the spirit of the modal approach and the transfer functions are given as an infinite series of modal terms. In others, e.g. [13], closed form expressions where derived only for a free-free system. A general method, based on a matrix exponent in the s domain and applicable to a wide range of systems, is presented in a series of works of Yang and colleagues [16,17], that are also summarized in [18].
In [19][20][21][22] a method for Laplace transfer function modeling of flexible structures governed by the wave equations was presented. The expressions are accurate, closed form and compact, consisting on second order rational expressions and delays. These transfer functions provide insight into the mechanism of the response, resulting in dedicated control laws with special properties. Since the transfer functions accurately model the system at all frequencies, control laws that are based on them, e.g. [21], do not suffer from spillover, resulting from unmodeled dynamics [23,24]. The transfer functions are also easy to use for simulation and analysis purposes without the need to refer to finite dimension approximations. From a physical point of view, the transfer functions correspond to the traveling wave approach [7,8,12].
First attempts to increase the approach to multi-link structures were made in [20] where the individual infinite dimension, accurate, Laplace transfer functions were collected in a systematic way to model a system of multiple shafts and flywheels. In this paper, the use of these Laplace transfer functions is extended to model flexible transmissions. In addition to the added complexity by the flexible gears, the paper extends the 'feedback approach' in [20] from two links to a general structure with any number of links. The modified approach is based in this paper on matrix representation rather then equivalent single-input single-output closed loop (feedback) representations. Therefore the term 'feedback' is no longer valid and the modified method is denoted by 'link reaction approach'. Preliminary and partial results on modeling flexible gear systems were given in [22].
An advantage of the transfer function modeling approach is that damping in the bearings enters the transfer function via the boundary conditions. Unlike the modal approach, the incorporation of damping at the boundary conditions is a natural, integral part of the model and does not require any further assumptions (except linearity). As structural damping is often negligible, damping in the boundary conditions, and in the gears, is most often the dominant damping in the system.
The derivation begins by presenting some results for a single shaft system. Next these results are extended to multi-link flexible transmissions, using the link reaction approach. The obtained model is an infinite dimension transfer function relating the angle at any point along the transmission to the input and load external torques. All other quantities such as angular velocity and internal moments, stress and strain are immediately available from the angle information.
The resulting model structure allows the development of a time simulation scheme of the multi-link transmissions which is based on generic link blocks. These blocks may be realized using standard ODE solvers and simulation tools such as Matlab-Simulink, without the need to use more complicated PDE solvers, thus allowing the simulation of constrained and free response of both damped and undamped systems.
The paper is organized as follows. The transfer function of a single link is introduced in Section 2. Section 3 contains the main results of the paper, namely the derivation of the transfer function of the entire transmission. Time domain simulation is discussed in Section 4. Examples and comparison with lumped and finite element modeling, are given in Section 5. The results are summarized and discussed in Section 6.

Transfer function of a single flexible shaft
Consider the single flexible link (shaft) of length L subjected to a lumped moment M in (t), at x = x 0 , shown in Fig. 2. The dynamics of the link is governed by the wave equation [7].
where I p is the shaft polar moment of inertia, ρ is the link density, G is the shear elasticity modulus and the wave propagation velocity in the link is c = (G/ρ) 1/2 . Two inertias, J 1 and J 2 are connected to the link at the ends and K 1 , K 2 , D 1 and D 2 , if exist, are the springs and dampers connecting the link to a fixed frame (skyhooks). The boundary conditions are therefore  This general setting includes all linear boundary conditions of interest. For example, in a free end all the coefficients are identically zero and a fixed end is obtained when K → ∞. Using Laplace transform on the partial differential Eq. (1) converts it to an ordinary differential equation in x. Solving this equation leads to the general transfer function from the input moment M in (s), acting at x = x 0 , to the torsion angle at any point x along the link [19,20].
where Table 1 lists some specific transfer functions that are used in the next sections. Notice the reciprocity existing in these expressions, e.g. the transfer function from a moment acting at x 0 = 0 to the torsion angle at x 1 = L is identical to that for x 0 = L and x 1 = 0. The boundary conditions are represented by R i , which are dynamic reflection coefficients [21]. For example, for a free end R i = 1 and for a fixed end R i = −1. The velocity transfer function is similar, without the s in the denominator.  Fig. 3. Natural feedback between links 1 and 2.

General theory
We move now to deriving the overall model of the system in Fig. 1, where the flexible transmissions are assumed to be free of backlash and shaft bending. One way to extend the results from single to multiple links is the link reaction approach that makes use of the natural feedback that exists in the system. Consider the first two links in Fig. 1. The first link rotates the inertia J 21 in Fig. 1 by an angle θ 2 (0) and as a result, the second link applies a restoring torque M 2 on the first link and so on, see Fig. 3. This should not be confused with the dynamics within the link, where a wave generated by a torque acting at any point along the link is reflected from the boundary conditions to the actuation point.
The first step in deriving the system Laplace transfer function is to separate each link from the system and replace its interaction with its neighbors by the reactions at its ends. Then, using superposition of the responses to reactions and external (input) moments in the link, the Laplace transfer function for any point along each link is computed. Finally, the transfer function of the whole system is obtained.
Consider a multi-link flexible transmission consisting of n flexible shafts, connected by flexible transmissions, as in Fig. 1. A free body diagram of the i-th link of such a system is shown in Fig. 4. Each link may also be connected to the environment by linear, or linearized, dampers and springs, representing the effect of the bearings damping and flexibility, and, at least theoretically, by skyhook springs.
We make the standard assumption that the only flexibility of the gear pair is in the meshing teeth [2][3][4][5][6], and that the transmission is free of backlash. It is further assumed that there is no shaft bending. Following the approach in [3,4], the nonlinear flexibility of the meshing teeth and the time varying behavior are modeled by the input of a displacement excitation e i (t), as shown in Fig. 5. Other periodic disturbances which are not included in the pure torsion model, may be included in e i (t) as well. The specific choice of e i (t) is discussed in [3]. k i and c i in Fig. 5 are the time average of the mesh stiffness and damping.
Using Newton's second law, with uniform positive direction for all the angles, we calculate the force in the springdamper element as a result of the gears rotation. All the quantities are Laplace transforms and the variable 's' is omitted for the ease of notation. This force generates moments acting on the gears inertia as follows Where the subscripts 1 and 2 denote the left (input) and right (output) ends of the shaft respectively.
is the transmission ratio between links i and i -1. The angle along the i-th link is a superposition of the responses to the reactions acting at its ends.
Defining M i ≡ M i,1 and combining Eqs (6) and (8) we have the moment equation The arguments of the transfer function that multiplies M i+1 in Eq. (9) were reversed, using the reciprocity property. The various transfer functions G i (x 0i , x 1i ) in Eqs (8)- (12) are listed in Table 1, and R i,1 and R i,2 represent the inertia of the gears connected to the shaft and the bearing damping.
Having the individual transfer functions and the connectivity between them, we are now in a position to derive the Laplace transfer function of the full system in Fig. 1. Equations (8)-(12) may be written in a general form as where the input, u ∈ 2 , the generalized state vector, z ∈ n−1 , the displacement excitation, e ∈ n−1 , the output, y ∈ n and the system matrices, are given by 0 · · · 0 a n−2Ḡn−1 a n−1 0 0 · · · 0 a n−1Ḡn Equations (13)- (14) are in a familiar form that often appears in modeling of dynamic systems. Equation (13) is the generalized state equation, representing the system dynamics. P (s) models the internal dynamics and Q(s) and E(s) are the dynamic input matrices. Notice that all of these quantities are independent of the locations of the measurements x i . Equation (14) is the dynamic output equation. If only some pre-specified angles are of interest, R(x,s) and V (x,s) reduce to the relevant rows of the full matrices in Eqs (22)- (23). The transfer functions, from u(s) and e(s) to y(s), are given by Equation (24) gives the explicit, accurate, infinite dimension transfer function of the system. An immediate application of it is the calculation of the Frequency Response Function (FRF) of the system. The first advantage is accuracy at all frequency ranges. This is in contrast with FEM models whose range of acceptable accuracy depends on their size and are always inaccurate at higher frequencies. High frequency behavior is important for example in feedback control design, where the effect of unmodeled high frequencies, known as spillover [23,24], is an important issue. Having the full FRF, known in controls context as Bode plot, one can see clearly the effect of the control law on all natural frequencies. The second advantage is the computational load. Obtaining the FRF from Eq. (24) requires inversion of an (n − 1) × (n − 1) matrix. In FEM based models the FRF requires inversion, actually equivalent yet more efficient algorithms, of a matrix with the order of the number of dof's (twice for non-proportional damping), which can be many orders of magnitude larger. Equation (24) indicates that the system poles are the solution of the characteristic equation det (P (s)) = 0 For conservative systems (D i,{1,2} = 0) all the poles are on the imaginary axis, and the natural frequencies are calculated from det(P (iω)) = 0. Here again we have the advantage that all the natural frequencies can be calculated accurately, without the degradation at the high frequencies that always exists in finite dimension approximations. The starting points for the numerical algorithm, e.g. Newton-Raphson or bi-section, can be obtained by simply observing the plot of the determinant vs. ω, or an automated version of that. The efficient algorithms for solving the finite dimension eigenvalue problem, at least for conservative systems, might be of comparable computation effort for small to medium size models. However as the order increases they inevitably require larger computation effort than solving Eq. (25) whose dimension is just the number of gear pairs.
Remark. The method of assembling the individual transfer functions into an overall model is called link reaction approach because the basic unit is the full link (shaft and two gears). In [20], it was denoted as the 'Feedback Approach' since for two links one may draw feedback block diagrams and use standard feedback formulas to obtain the overall transfer function. A different formulation, denoted as 'Dynamic Stiffness Approach' [20], uses the inertias between the shafts as building blocks. It results in a model with a general structure similar to Eqs (13)- (14) but with different variables and matrices.

Other configurations
The derivation so far can be used, with slight modifications, to obtain models for other configurations that may exist in a flexible gear system.

Rigid gears
The transfer function for rigid gears is obtained by substituting k i → ∞ and e i = 0 into Eq. (10) that becomes where a i is still as given in Eq. (12), andḠ i is replaced bȳ The combined model Eqs (13) x 1

Inertia in the middle of the shaft
Consider the system shown in Fig. 6(a) where J 2 is an inertia other than a gear, e.g. flywheel, with possible skyhook damper and spring. Since the entire derivation is based on shafts with uniform cross-section, a change in diameter is also modeled as the system in Fig. 6, with J 2 = 0. The side inertias J 11 and J 22 are the gears connecting this sub-structure to the rest of the system.
Splitting artificially J 2 into J 12 and J 21 as shown of Fig. 6(b) we arrive at a situation similar to a rigid transmission with N 2 = −1 since there is no change in the rotation direction. The model is then obtained as in the previous case. It is straightforward to show that the overall model is independent of the particular partitioning of the lumped inertia and damping. Any partition of the form and similar expressions for D 1,2 , D 2,1 and K 1,2 , K 2,1 , even with w J = w D = w K , results in the same transfer functions to the angles of the two flexible shafts.

Flexible joint
A configuration of a flexible joint is shown in Fig. 7, where k i and c i are now torsional spring and damper. A derivation, similar to the flexible gears case, shows that the model for such a system is obtained by substituting in Eqs (10)-(12) e i = 0, N i = −1 (no change in the rotation direction), and replacing the block

Input moment to load side angle
The most useful transfer function seems to be the one from M in , the input moment acting on the first link, to θ n,2 , the angle at the load side, for e = 0 and no load moment (i.e. M load = 0). The output matrix R(x,s) reduces in this case to its last row, where only the last component is nonzero. Consequently, only the (1, n − 1) component of P −1 (s) has to be calculated. From the structure of P (s) it follows that the transfer function is given as
where N = ΠN i is the overall transmission ratio. The 'closed loop' transfer function Eq. (29) has the forward path in its numerator, and a sum of combined paths in the denominator, which is in accordance with Mason's rule for multi-feedback schemes.

Two link system
Consider the simplest transmission, consisting of two shafts connected by a single gear pair. We wish to find explicit expressions for the angles of both shafts for the flexible and rigid gear case. In this simplified case, the model matrices reduce to the scalars P (s) =Ḡ 2 , Q(s) = [−a 1 G 2 (l 2 , 0)] (the first and last rows coincide) and E = −1/r 2,1 . Hence, Eq. (13) may be written as As was already explained, the solution for the rigid gears case is recovered by substituting e 2 = 0,and k 2 → ∞.
The result in this case is the removal of the last terms in the numerator and the denominator of both expressions in Eq. (31).

Simulation
One of the main advantages of the transfer function modeling approach is that it enables the construction of simple and efficient simulation schemes for the accurate, infinite dimension, model of the entire transmission. Starting at a single link level, notice that the denominator of the transfer function Eq. (3) can be treated as a result of an internal feedback in a system consisting of rational transfer functions (based on R i (s)) and delays. This is demonstrated in Fig. 8 for G(0,l). Since rational transfer functions and delays are standard blocks in numerical solvers such as Matlab -Simulink, the scheme can be easily implemented.  Fig. 9. Generic, single link block diagram.
The implicit form of Eqs (9)-(10) is convenient for analysis and leads to the generalized state Eq. (13). For simulation purposes however, an explicit expression for the moment M i is required, which is immediate from Eq. (9) as Equation (32), together with its output Eq. (8), may be represented by the block diagram in Fig. 9. Notice that G i (0,0) is a bi-causal (delay wise) transfer function hence its inverse can be realized. To represent a system of n links, n − 1 single-link models are connected in that manner to each other. For the first link M 1 = M in , and for the last one M n+1 = M load . A block diagram for a four link system is shown in Fig. 10, where each grey block represents the grey area in Fig. 9. The feedback nature of the aggregated system is emphasized by Fig. 10, where each link is subjected to torques from the neighboring links. The block diagram may be adapted for the rigid transmission case by setting e i = 0 and k i → ∞, which is equivalent to removing the green (dashed) lines and blocks from Figs 9 and 10. Thus the scheme is applicable to any combination of flexible gears, rigid gears, shafts with middle inertia and flexible joints. Fig. 10. Block diagram of a multi-link system.

Example 1
Consider the gear system in Fig. 11, consisting of three flexible shafts, with lengths 0.21, 0.14, 0.14 m, diameters 0.025, 0.05, 0.05 m, density 7800 kg/m 3 and shear modulus 7.95·10 10 Pa. The two gear pairs are identical with d 1 = 0.047 m, d 2 = 0.094 m, J 1 = 1.15·10 −3 kg·m 2 and J 2 = 2.3·10 −3 kg·m 2 . The actuator inertia is 5.75 · 10 −3 kg·m 2 and load inertia is 2.3·10 −2 kg·m 2 . The gears flexibility was modeled by k 2 = 200 ·10 6 N/m and k 3 = 300·10 6 N/m. The first 10 natural frequencies, calculated from Eq. (25), were compared to lumped approximation and FEM with linear bar elements. The lumped model is made of 6 inertias, which include the shafts inertias, connected by springs representing the shafts and gears flexibility. The results are displayed in Table 2. The lumped model gives a fairly good approximation of the first 2 natural frequencies and an acceptable approximation of the next 3. As expected, the accuracy of the FEM model increases with the number of elements and decreases for the higher frequencies. The transfer function model suggested in this paper gives any required number of natural frequencies, in any practical frequency range, with absolute accuracy up to numerical issues. It is worth noting that the transfer function model gives also the exact full frequency response function (FRF) of the system with no additional computational effort.
The same procedure was repeated when all the shafts are 5 times longer, i.e. 1.05, 0.7, 0.7 m. The results of that case are given in Table 3. The lumped model completely misses the 4-th natural frequency and the accuracy for the next frequencies is poor. Similar to the lumped model, FEM with the same number of elements provides less accurate estimates for the longer shafts.
The system response to an input moment was simulated twice: using the transfer function model, as derived in Section 4, and a lumped approximation such as that of the previous example.  The simulation results are shown in Fig. 12. While the lumped model captures the general behavior of the system response, it misses the high frequency content of the response. The transfer function model results show the complex behavior of the response. As expected, the effect of the pure time delay may only be observed when using this model. Figure 13 shows the force acting on the gears teeth. The same conclusions, regarding accuracy and delay, may be drawn here as well. It should be noted that accurate modeling of the high frequency modes is important in predicting the gears wear. Figure 14 repeats the simulation for the case where all the damping elements are reduced by 50 percent.
The simulation results in Figs 12-14 clearly show the effect of damping on the system. The vibrations energy is absorbed in the dampers and the response is smoother than with smaller or no damping. At steady state the input moment equals the equivalent moment applied by the dampers and the system moves at a constant speed with no vibration. The transient parts of the time response in these figures are analogous to the free response of the system, e.g. when the input moment is turned off.

Conclusion
A new method for a simple and exact continuous (infinite dimension) modeling of a flexible transmission is proposed. The method applies to any number of gear pairs and is built in a modular fashion from the single link transfer functions. The models were derived by combining, in a systematic manner, the full, infinite dimension models of the individual shafts. The basic units in building the equations of motion of the overall system are the links, consisting of a shaft connected to two inertias. Hence the name 'link reaction approach'. The overall model is given by generalized state equations where the generalized state vector consists of the reaction moments between the links, and therefore has a dimension of n − 1, for an n is links system.   The gears were assumed to be rigid inertias except for the meshing teeth that were represented by spring-damper and deflection excitation. By increasing the gear flexibility to infinity, the models for flexible gears converge to those for rigid gears. The proposed model is general and can include several other cases such as a change in the shaft diameter or a spring damper connection instead of gears, etc. Mixed kinds of connections between the links may also be represented by connecting the appropriate single link models.
An immediate application of the transfer function model is the calculation of the natural frequencies and frequency response function of the system. The FRF is exact at all frequency ranges and the required computation is much The compact form and the natural feedback built into the model makes it also a useful tool for time domain simulation. As demonstrated in Section 4, generic link model may be created and the system is then modeled by connecting several such sub-models. This enables the use of standard software, suitable for ODE's solutions with time delay option, instead of the "heavier" PDE's solvers.