A new analytical model to predict the transversal deﬂection under load of stepped shafts

Shaft deﬂection is a common phenomenon in machine design that has an important inﬂuence on the behavior of many transmission elements supported by the shafts, like gears, pulleys, sprockets, etc. This deﬂection can be estimated eﬃciently by using any beam theory but machine shafts are usually stepped shafts and it has been demonstrated that beam theories do not predict accurately the deﬂection of this type of shafts. Thus, in this work, an equivalent model of the shaft to be used in conjunction with the Timoshenko beam theory has been proposed to improve the accuracy of the computed deﬂections. The new model substitutes the steps of the shaft by a linear variation of the diameter of the cross-section, removing the discontinuities caused by these steps. The slope of the linear variation provides new variables (angles) that are optimized locally and globally for the best coincidence with the deﬂection obtained with a realistic ﬁnite element model of the shaft. Then, two diﬀerent approaches (with locally and globally optimized angles) are proposed and their accuracy is investigated with 56 cases of study and with the analysis of a realistic gear drive. The results demonstrate the high accuracy of the proposed model and the improvement with respect to the direct application of the Timoshenko beam theory.


Introduction (revised)
Computing accurately the transversal deflection of shafts is of major importance in many mechanical applications.An example of this is the stress analysis of gear transmissions.Gears are extremely sensitive to variations of the load distribution along its teeth [1], as small changes in the tensional state may result in substantial shorter lifespans and premature failure, as well as vibrations, noise, wear and an increase of heat generation.And this load distribution is highly dependent on the alignment status of the gears [2][3][4][5], which is influenced by the deformation of all structural parts of the transmission.For this reason, a stress analysis of gear transmissions should consider the flexibility of gears, shafts, rolling bearings and even the box that contains them all.This can be done extensively by using a realistic 3D finite element (FE) model of the whole gear drive, but this approach is laborious and computationally expensive, slowing down the design process.
A faster solution can be obtained with simplified models, where the flexibility of some elements of the transmission is considered.Among them, the model including only gears and shafts (assuming the last ones as simply supported beams) is frequently used [6][7][8] to have a first approach to the contact and structural problem taking into account the natural misalignment of the gears during the transmission of power.In this model, the use of beam theories [9] to compute the transversal deflection of shafts increases the efficiency of the approach, constituting an interesting alternative to realistic FE models.
The classical and most widely used beam theories are the ones proposed by Euler-Bernoulli, Saint-Venant and Timoshenko.The first two (EBT and SVBT) do not consider transverse shear deformations and are frequently discarded in short and thick beams because the relative contribution of the shear in the deflection of the beam can be significant.The Timoshenko beam theory (TBT) considers first order shear deformation along the cross-section of the beam, constituting a better approximation to the deflection of short and thick beams, but it requires the use of a shear correction factor in order to compensate for the error caused by assuming a constant transverse shear stress distribution along the beam depth.In this line, several research works [10][11][12][13] were directed to improve the global response of beam theories based on the TBT through the use of the appropriate shear correction factors.In other works [14,15] this research line was extended by computing shear correction factors for some particular cases like arbitrary shaped cross-sections; torsional and flexural shearing stresses in prismatic beams; and wide, thin-walled, and bridge-like structures.
Apart from the shear issue, there are many other effects that are not considered by the classical theories, such as warping, out-plane and in-plane deformations, torsion-bending coupling, localized boundary conditions, etc., that are usually due to small slenderness ratios, thin walls, and the anisotropy of the materials.Thus, a number of methods have been proposed [16] to overcome these limitations and to extend the applicability of the beam theories.
Within these methods, El Fatmi [17][18][19] presented a beam theory with a non-uniform warping including the effects of torsion and shear forces, and valid for any homogeneous cross-section made of isotropic elastic material.Berdichevsky [20] developed a variationally and asymptotically consistent theory in order to derive the governing equations of anisotropic thin-walled beams with closed sections.The theory is based on an asymptotic analysis of two-dimensional shell theory and represents the starting point of an alternative approach to constructing refined beam theories where a characteristic parameter (e.g., the cross-sectional thickness of a beam) is exploited to build an asymptotic series.Other valuable contributions on asymptotic methods are those related to variational asymptotic beam sectional (VABS) models developed by Volovoi, Hodges, Popescu, Yu and other collaborators [21][22][23][24][25][26][27].
On the other hand, generalized beam theory (GBT) was developed from conventional beam theory to consider cross-section distortion.It was originally proposed by Schardt [28], extended by Davies [29] and developed extensively by Silvetre and Camotim and their co-workers [30,31] during the last years.As a consequence of this development, GBTs can be applied without restrictions to any prismatic thin-walled structural member exhibiting straight or curved axial axis (any loading, any cross-section geometry, any boundary conditions).Apart from GBT, many other higher-order theories which are based on enhanced displacement fields over the beam cross-section have been introduced to include non-classical effects.
Beam theories can be applied to estimate the deflection of machine shafts under load, but shafts have an important particularity that may condition the accuracy of the results: they usually have steps to provide shoulders to support rings, rolling bearings, etc. and have elements (gears, pulleys, etc.) that generate additional steps.These steps constitute abrupt changes in the cross-section that lead to discontinuities along the shaft, generating two main problems related to the application of beam theories: • Integration problem: the discontinuities make difficult to obtain a unique closed-form expression for the elastic curve of the whole shaft from the integration of the differential equations.
• Lack of accuracy: as it will be proved in section 2, the classical beam theories have a poor accuracy when estimating the deflections of stepped shafts.The reason is that stress concentration in the steps modifies the strain field along the shaft and it is not properly described by the mathematical model of these beam theories.
The first problem has been investigated by researchers in order to find a solution to certain cases.Romano [32] solved the integration problem obtaining a closed-form solution by applying the EBT to beams with circular cross-section and linear and parabolic variation of the diameter.Then [33], he provided a similar solution for the TBT and beams with rectangular cross-section and linear and parabolic variations of width and depth.Cueva [34] also solved the integration problem for the EBT by using Macaulay functions to deal with the singularities generated by the steps.Yavari [35] applied the distribution theory to obtain a single displacement function and a single rotation function for the EBT and the TBT in beams with jump discontinuities in slope, deflection, flexural stiffness and mechanical properties, identifying the cases where this is not possible.And Rencis [36] applied the classical method of segments to provide a solution procedure for axially and torsionally loaded stepped beams.
Regarding the lack of accuracy, the vibratory behavior of stepped beams have also been a recent topic for several researchers.Lee [37] proposed a Chebyshev-tau method based on both Euler-Bernoulli and Timoshenko beam theories for the free vibration analyses of stepped beams, Mao [38,39] employed the Adomian decomposition method to investigate the free vibrations of the Euler-Bernoulli beams with multiple cross-section steps and Naguleswaran [40] proposed an analytical method to calculate the frequencies of beams with up to three step changes in cross-section; to cite a few.
But it was Sanderson [41][42][43] one of the first researchers that investigated this lack of accuracy, proposing a new model of equivalent shaft with a variable second moment of area (inertia) that can be used to obtain the elastic curve of the shaft by means of the EBT.The study was performed by comparing finite element models with experimental tests and, as a result, a graph with different empiric curves (for different geometries of the step) representing the evolution of the inertia around the region of the step to build the equivalent shaft was proposed.However, the model was designed not to have shear load in the steps and the complexity of the resulting diagrams of equivalent inertia makes difficult to automate the computation for a general case of stepped shaft.
The lack of accuracy of the classical beam theories in stepped shafts can be tackled with two different strategies: (i) the use or development of beam theories with higher order mathematical models able to describe the complex strain field in stepped shafts and (ii) the use of simpler beam theories with an adapted physical model that represents the true stiffness of each part of the shaft.The research of Sanderson inspired the work presented here that is in line with this second strategy.Thus, the objective is to develop an equivalent physical model of the stepped shaft that, combined with the Timoshenko beam theory, allows to compute the elastic curve of the shaft with a much higher level of accuracy than the one obtained when the classical beam theories are applied to a realistic model of the shaft.The main advantage of using the TBT is the simplicity of the formulation while considering bending and shear effects (being both important in short stepped shafts).Furthermore, the developed physical model pretends to be simple and general, to make easy the automated computation of the elastic curve in any stepped shaft.

Modified Timoshenko beam model for stepped shafts
Unlike the EBT, the TBT does not assume that the cross-sections remain normal to the neutral axis after deformation.Thus, the cross-section may rotate an angle φ and a constant transverse shear strain distribution along the thickness coordinate of the beam is considered.As a result of this assumption, a shear correction factor must be applied in order to compensate the error caused.This factor depends on the material and the geometry of the beam, as well as the loading and boundary conditions [9].
Considering a beam with a coordinate system where X-axis is placed along the thickness of the beam, Y -axis is located along its width and Z-axis along its length, the differential equations governing the bending moment M x and the shear force Q y provided by the TBT are: where ω is the transverse deflection of a point in the neutral axis, E is the elastic modulus, I x is the second moment of area with respect to the X-axis, G yz is the shear modulus in the Y Z plane and A is the area of the cross-section.K s represents the shear correction factor required to compensate the error caused by the assumed transverse shear stress distribution.For circular cross-section beams a value of K s = 0.89 is generally adopted, as suggested by Cowper [10].
If the beam is loaded with a distributed load q in the negative direction of the Y -axis, the equilibrium equations of the beam in plane YZ are: Subtituting M x and Q y from Eqs. ( 1) and (2) in Eqs. ( 3) and (4), a system of two differential equations are obtained.These equations can be integrated to obtain φ first and ω next.Assuming that E and G yz are both constant, the resulting expression of the deflection is: where C 1 , C 2 , C 3 , and C 4 are constants of integration that can be obtained from the boundary conditions of the beam and from connectivity conditions between segments of the beam.Despite the TBT predicts very well the deflection of beams with uniform circular cross-sections, some investigations demonstrated that the direct application of the TBT to stepped shafts does not provide accurate results [43,44].To illustrate this limitation, a stepped shaft was compared to a solid cylindric shaft, both made of steel, highly loaded and equally long (Fig. 1).In both cases, the elastic curve of the shafts was calculated by means of the TBT and also by solving the structural problem from a realistic finite element (FE) model (details of this model are provided in section 3.2).It is important to remark that the Timoshenko beam model (TBM) takes into account the area and inertia of the cross-sections of the actual stepped shaft, reproducing accurately the geometry of the steps.And, here, the application of the Timoshenko beam theory to the Timoshenko beam model is called Timoshenko beam approach (TB approach).The results confirm the influence of the step in the accuracy of the elastic curve obtained by the TB approach.In the case of a solid cylinder (Fig. 1a), where there is no presence of steps in the shaft, both approaches TB and FE have a very high degree of coincidence.However, in the case of the stepped shaft (Fig. 1b) an error is induced which makes the TBM significantly more rigid than the realistic FE model.And it has been observed that this difference becomes more important as the diameter variation of the step increases.
Because of the reasons stated before, a new beam model based on the Timoshenko beam theory (called here modified Timoshenko beam model, mTBM) has been developed regarding its application to stepped shafts.The new model emerged from analyzing the stress field of a stepped shaft under bending load (Fig. 2) and after stating that the axial stress vanishes in the corners of each abrupt change.This means that these regions do not contribute significantly to the stiffness of the shaft but they are fully considered in the original TBM, increasing the error in stepped shafts.Taking this into account, the new model deals with abrupt changes in the cross-section by discarding these corners and, thus, considering a progressive linear change of the section (conic segment) in each step.By doing this, a continuous transition between the properties of the sections before and after the step is modeled and the discontinuity resulting from the geometry is avoided.These new intermediate conic segments are characterized by a cone angle α that is initially unknown, and the determination of its optimal value is one of the main objectives of this work.In a stepped shaft, each step can be characterized by the ratios b/d and D/d, where b, d and D are defined in Fig. 3.Then, the new model is built from the geometry of the shaft and, depending on the ratios, different situations may arise.The most common situations are the ones presented in Fig. 3, where the geometry of the proposed mTBM (shaded) is compared to the actual geometry of the shaft (outline), that coincides with the TBM.The description of these situations is the following: • When the step is characterized by high b/d and low D/d ratios (Fig. 3a), the higher diameter segment of the shaft is defined by two conic segments and a middle cylindrical one.
• When the step is characterized by low b/d and high D/d ratios (Fig. 3b), the higher diameter segment of the shaft is defined by two opposed conic segments.
• When the lower diameters at both sides of the step are different (Fig. 3c), different angles α are used and the conic segments end at the higher diameter, like in Fig. 3a, or when they intersect each other, like in Fig. 3b.The variation of the diameter in the new model as a function of the z-coordinate can be easily obtained by means of a simple trigonometric analysis.In a generic step (Fig. 4), the z-coordinate of the theoretical intersection of both cones (z i ) can be obtained by means of Eq. ( 6) and the diameter can be expressed as a piecewise function as shown in Eq. (7).When there are nested steps, this analysis can be repeated to obtain the full piecewise function of the diameter.In order to apply the theory of Timoshenko to the conic segments of the model, let us consider a segment from a tapered beam with a circular cross-section as the one illustrated in Fig. 5, with L being its length, q(z) a generic distributed load applied along its length and d(z) its diameter.
The cross-section of this conic segment is defined by Eq. ( 8), where d(z) is the diameter as a function of z-coordinate; d 0 is starting diameter and α is the slope angle of the conic segment.
Consequently, the cross-section area A and moment of inertia I x will also gradually change along the length of the conic segment, through the diameter: But the influence of both area and moment of inertia over the transversal deflection of the model can be different and it is directly reflected in the shear load and bending moment applied, respectively.Depending on the length of the shaft, the bending moment and the shear load have different effects: in short beams the deformation caused by the shear load is usually significant compared to the deformation caused by the bending moment, whereas in long beams the deformation caused by the shear load is usually negligible.Taking this into consideration, two separate angles are used to define this distinct behavior in the conic section: α A to describe the variation of the cross-section area and α I that of the moment of inertia.This is equivalent to considering two different cone segments in the model, one to describe the variation of the area and another to describe the variation of the inertia of the cross-section.Thus, Eq. ( 9) and ( 10) can be written in terms of these angles: These expressions have been implemented in Eq. ( 5), whose terms have been integrated analytically for a constant distributed load q(z) = q, obtaining the equation of the deflection for the conic segment with the constants of integration as unknowns.For cylindric segments, the obtained expression of the deflection could also be applied by imposing α A = α I = 0, but it generates a mathematical singularity.Authors like Romano [32,33] noticed this issue and attempted to obtain a formulation compatible with both conic and cylindric segments by applying McLaurin series.However, in this work, Eq. ( 5) has been solved again for cylindric segments (where I x and A are constant) to avoid the singularity.Thus, different equations for the deflection of conic and cylindric segments are used.
In the most general case where the distributed load q(z) is not constant, Eq. ( 5) cannot be solved analytically but its terms can be integrated numerically to obtain the deflection of the segment.The computational cost is slightly higher, but the approach can be applied in the same way.
In summary, the new proposed model is composed of several cylindric and conic segments.To evaluate the transversal deflection of this beam model, the method of segments [36] is used.In this method, the deflection of the different segments is evaluated separately and the constants of integration are determined by establishing continuity conditions between adjacent segments, as well as boundary conditions that simulate the supports of the shafts.

Model optimization
As described in the previous section, the values of α A and α I are the unknown variables of the new mTBM and they define the regions of the stepped shaft that are not considered by the model when computing the elastic curve because they have a low contribution to the stiffness of the shaft.It is possible that the optimal values of these angles (those that provide the most realistic elastic curve) are dependent on the ratios that define the step (b/d and D/d in Fig. 3).This possibility has been investigated by calibrating the model with a set of cases, where each case is composed of a simple stepped shaft with a distributed load and only the geometry of the step changes from one case to others.Thus, for each case, the following steps have been performed: • A realistic finite element model of the stepped shaft has been built, the structural problem has been solved and the resulting elastic curve has been obtained.This curve has been considered as the reference elastic curve.
• An mTBM of the stepped shaft has been built following the instructions of the proposed approach that was described in section 2. This model depends on α A and α I (explanatory set of variables).It has been solved with the Timoshenko beam theory and the approximated elastic curve of the shaft (dependent variable) has been obtained.
• An optimization process (calibration) has been performed to compute the optimal value of the explanatory variables for the approximated elastic curve to best fit the reference elastic curve, quantified by a comparison error that is minimized.
Thus, the result of this investigation is the optimal value of α A and α I for each case, that will provide information about the variation of the optimal angles with the geometry of the step.

Set of cases
The optimization has been performed with 56 different cases representing typical step geometries that are commonly found in machine shafts, in general, and in gear shafts, in particular.The geometric and load characteristics of these cases are based on the same model of the shaft, shown in Fig. 6.This base model consists of a symmetrical shaft made of steel of total length L = 270 mm which is simply supported at both ends.The shaft is made up of two external segments with a circular cross-section of diameter d = 50 mm and a central segment with a constant cross-section of length b and diameter D, being D > d.A high load is convenient to have significant deflections and avoid numerical problems.Thus, a uniform distributed load F = 50 kN is applied at the central segment.Finally, each case is defined by a unique combination of D/d and b/d ratios.However, not all steps of this space are equally found in real shafts.After a brief analysis, it has been observed that most steps are included in two regions of the studied space.Region 1 includes steps that are wide (b/d ≥ 2) and with a small change of diameter (D/d ≤ 1.5).These steps are very common to provide shoulders for bearings, rings and other elements.On the other hand, region 2 includes steps with a moderate ratio b/d (0.5 ≤ b/d ≤ 1.5) but with a wide range of values in D/d, that are common in gears, belt pulleys and other elements.Finally, the steps outside of these regions are considered uncommon and no special attention is paid to them in this study.

Reference FE model
To obtain the reference elastic curve, a realistic three-dimensional finite element model has been created for each case of Fig. 7.The characteristics of this model are the following: • The model (Fig. 8) has been generated with second order (i.e., quadratic) tetrahedral finite elements using generic meshing algorithms to accurately reproduce the geometry of the stepped shaft.The global size of the element has been specified as 3 mm and the fillets have been refined with 1 mm elements.Furthermore, meshing restrictions have been included to ensure that there is a line of nodes in the neutral axis of the shaft.These nodes have been used to obtain the resulting (discrete) elastic curve of the deformed shaft.
• Rigid surfaces have been created at each shaft end.The movement of these surfaces, which contain the corresponding nodes, has been coupled to the movement of the reference nodes A and B (Fig. 8).
• Simply supported boundary conditions have been established by restricting the proper degrees of freedom of the reference nodes.
• A total load F has been introduced as a distributed load to avoid stress concentration, as presented in Fig. 8.
After generating the model, the FE problem has been solved by using the generic FE solver Abaqus and the reference elastic curve has been obtained from the displacements of the nodes of the neutral axis of the shaft.

Error measurement
With the analysis of each case, two elastic curves are obtained: the reference elastic curve from the FE model and the approximated elastic curve from the application of the new proposed approach, being this last one dependent on α I and α A , as explained in section 2. When these two curves must be compared for a specific purpose, there are several possibilities and, in this work, two different error ratios are proposed.The first one estimates the error of the approximated elastic curve compared to the reference one as follows: where A approx is the area enclosed between the approximated elastic curve and the undeformed neutral axis, A ref is the corresponding area in the reference elastic curve and L is the length of the shaft (Fig. 9).Both areas are computed numerically (trapezoidal rule method) from the results of the compared approaches.This error can be described as the area enclosed between the two elastic curves divided by the length of the shaft.Thus, the resulting ratio has length units and represents the average deflection error of the approximated elastic curve.Consequently, this error is dependent on the load applied to the shaft.A different way to measure the error of the approximated elastic curve when compared to the reference elastic curve is the following: This error represents the area between both elastic curves relative to the area under the reference elastic curve, it is dimensionless and can be expressed as percent.

Optimization results and discussion
In a first step, a local optimization of the mTBM has been performed.According to this, the unknown variables of the model (α I and α A ) have been determined to minimize the error ε def l individually for each  case of the selected set by using a Nelder-Mead simplex algorithm [45].The algorithm converged in all cases and the obtained optimal values are shown in Fig. 10, where it can be observed that, within the scope of this study, α I presents a variation from 30 • to 70 • and α A is 90 • in a wide area of the space but lowers up to 0 • for low values of b/d.Considering the regions of interest (Fig. 7), in region 1 (cyan) α I has a high gradient and α A has a uniform value of 90 • .On the other hand, in region 2 (green), α I changes from 30.2 • to 34.8 • and α A has two significant areas: one with a uniform value of 90 • and another with a moderate gradient lowering up to 13.7 • .In a second step, a global optimization process has been achieved.In this case, the variables α I and α A have been optimized to minimize the sum of errors (ε def l ) of all cases in the selected set.As a result, the unique global values of 34.872 • and 90 • have been obtained for α IG and α AG (where subscript "G" stands for "Globally optimized angle"), respectively.
With these results, when a practical case of a stepped shaft must be solved with the proposed mTBM, it is possible to use two approaches.The first one (called here modified Timoshenko beam model with local angles approach or mTBM-LA approach) consists in computing the parameters D/d and b/d of each step of the shaft, interpolating the optimal values of α I and α A from the graphs in Fig. 10 and using these values to build the model of the shaft.The second strategy (called modified Timoshenko beam model with global angles approach or mTBM-GA approach) consist in using uniformly the values of α IG and α AG obtained with the global optimization in all steps of the shaft to build the model.After that, the elastic curve is computed by means of the TBT in both cases.
To compare these two approaches, each case of the set presented in section 3.1 has been solved with both of them and the obtained elastic curve has been compared with the reference one by computing the relative error (ε rel ).The results of this comparison are presented in Fig. 11.Additionally, all selected cases have also been solved with the original TB approach and, similarly, the obtained elastic curve has been compared with the reference elastic curve obtaining ε rel in all cases.The result of this comparison can be observed in Fig. 12.Finally, a summary of these results for the regions of interest are shown in Tab. 1, where the minimum, maximum and average values of ε rel can be observed.The first conclusion that can be obtained (Tab. 1) is that, as expected, the approach that provides the elastic curve best fitting the reference curve is mTBM-LA, followed by mTBM-GA.On the other hand, the TB approach proves to have results with important errors.
Another important conclusion is that the use of α IG and α AG in all steps of the shaft (mTBM-GA approach) provides errors almost as low as when locally optimized values for these angles are used (mTBM-LA approach), even in areas of the space where the globally optimized angles are very different from the locally optimized angles.For example, in the case where D/d = 1.2 and b/d = 3, the optimal angles (Fig. 10) are α I = 70.1 • and α A = 90 • and, with these angles, the mTBM-LA approach provides an elastic curve with a relative error ε rel = 0.36%.On the other hand, the mTBM-GA approach (that uses the angles α IG = 34.872• and α AG = 90 • ) produces an error as low as ε rel = 1.23%.The reason for the low error of the mTBM-GA approach despite the important difference in α I is that the elastic curve obtained with the mTBM is not equally sensitive to variations of α I (and α A ) in all cases of the considered set.This is justified in Fig. 13, where three cases have been analyzed.The figure shows, in the upper part, the actual geometry of each case that coincides with the original TBM.In the lower part, the effective longitudinal sections obtained by the mTBM-GA (gray) and mTBM-LA (black) approaches are overlapped.
It can be observed that, when D/d is low and b/d is high (Fig. 13a), a change of the considered angle α implies a small difference in the longitudinal section of the model and, consequently, in its global stiffness.And this sensitivity of the section to the angle α is lower for high values of this angle.As a consequence of this, the obtained values of ε rel with the mTBM-LA, mTBM-GA and TB approaches are small and very similar, despite the values of α I and α A used by these models being quite different.
On the other hand, when D/d is high and b/d is low (Fig. 13b), the moderate change of the angle α I between mTBM-LA (α I = 44.6 • ) and mTBM-GA (α IG = 34.872• ) implies a small variation in the longitudinal section.Thus, the values of ε rel obtained with these models (0.03% and 0.81%, respectively) are both similar.However, it is easy to observe that the sensitivity of the longitudinal section to α I (resp.α A ) increases when the value of this angle also increases.Thus, in the TBM (where α I = 90 • ) the longitudinal section changes significantly compared with the mTBM-LA model, and so does the obtained error (ε rel = 6.8%).
Finally, when both D/d and b/d are high (Fig. 13c), the longitudinal section of the mTBM proves to be very sensitive to angle α, so the resulting error obtained with the mTBM-GA approach (ε rel = 4.35%) is higher than the one obtained with the mTBM-LA (ε rel = 0.83%), but the error of the TB approach (ε rel = 31.0%) is much higher.
In conclusion, although the use of locally optimized values of α I and α A for each step (mTBM-LA approach) is the best strategy to obtain an accurate elastic curve of the stepped shaft, the use of the globally optimized values of the angles for all steps (mTBM-GA approach) constitutes a simpler approach that provides accurate results improving significantly the results of the classic Timoshenko beam model in stepped shafts.

Application to the gear drive transmission design
The full contact analysis of a gear transmission considering the flexibility of the shafts is a coupled nonlinear problem where the deformation of the shafts produces the misalignment of the gears, which changes the tooth contact pressure distribution, which, in turn, changes the deformation of the shafts.
This problem can be solved accurately with the finite element method, but this approach is slow and computationally expensive.Alternatively, a first approximation can be obtained by means of an iterative algorithm with simplified models (Fig. 14), under certain hypotheses, where both coupled phenomena are computed separately [7,8,46].In this algorithm, step 1 constitutes a structural problem with a static model of a flexible shaft and non-deformable gears that is used to compute the elastic curve of the shaft.On the other hand, in step 4, a contact model of a misaligned flexible gear pair (without shafts) is considered to compute the contact pressure distribution associated with the transmission of power.Solving these problems, the algorithm converges to a realistic solution, ending when no change is produced in the misalignment of the gears (step 3) according to certain tolerance.
Since the contact pattern of a gear pair is extremely sensitive to the misalignment of the gears and since this misalignment is highly dependent on the deformation of shafts, obtaining realistic elastic curves of the shafts is of greatest importance for the design of the gear transmission.Taking into consideration that gear shafts are usually stepped shafts, the model proposed in this work is very convenient for this problem.Thus, to demonstrate the benefits of the new model compared to the classical TBM, a new test with a realistic gear transmission has been performed.The selected transmission is shown in Fig. 15, where the upper part shows a drawing of the gear unit and the lower part shows the geometric information of gears and shafts considered in the problem.In addition, the specifications of the gears and the transmitted torque from input shaft to output shaft are presented in Tab. 2.
The test performed here consists in computing step 1 for the first iteration of the algorithm.According to this, the shafts have been analyzed separately and a uniform contact pressure distribution has been applied to the gears at the tooth contact surfaces.This initial contact pressure has been obtained from the contact forces in each gear pair (Tab.2), that have been calculated from the transmitted torque and the geometry of the gears.
Under these conditions, a realistic FE model of each shaft (similar to the model described in section 3.2) has been created.This model has been used to obtain a realistic elastic curve that will be used as reference.Similarly, the TB, mTBM-LA and mTBM-GA approaches have been used to compute the elastic curve of 0. Assume uniform pressure distributions in gear tooth contacts.
1. Consider separately each shaft with its gears (stepped shaft model).Apply the pressure distribution to gears in the shaft and compute the elastic curve of the shaft.
2. For each gear pair, compute the gear mesh misalignment from the elastic curves of the shafts.4. For each gear pair, using a model of the misaligned gears in contact (without shafts), perform a contact analysis to obtain the new contact pressure distribution.each shaft and, finally, a comparison between the approximated elastic curves with the reference one of each shaft has been performed to test the accuracy of the approaches.
To build the mTBM-LA model, the optimal angles α I and α A had to be estimated by interpolation (and extrapolation) for each step of the shafts (labeled in Fig. 15).These angles have been obtained from the graphs in figure 10  In the problem raised here, since all gears are spur gears and since the case corresponds to the first iteration of the algorithm (where the gears are perfectly aligned), the tooth contact line of each gear pair is parallel to the axis of the gear and the contact pressure is assumed to be uniform.Therefore, the external load of shaft 1 is planar and contained in the plane of action of the gear set 1-2 (Fig. 16).Choosing conveniently the axis Y S1 , the elastic curve will be contained in plane Z S1 -Y S1 .For the same reason, the elastic curve of the shaft 3 will be contained in plane Z S3 -Y S3 .But in the case of shaft 2, it contains two gears with different planes of actions, so the external load for this shaft is not planar and, as a consequence, the resulting elastic curve will be a spatial curve.Thus, an arbitrary set of axis (X S2 , Y S2 ) has been selected to represent the elastic curve of this shaft.
After solving the problem, the resulting elastic curves in the described planes are presented in Fig. 17, where the better approximation of the proposed models (mTBM-LA and mTBM-GA) compared to the TBM can be observed.To perform a quantitative comparison, the relative error ( rel ) has been computed (Tab.4).
From the results of this example, several conclusions can be extracted: • The error values of Tab. 4 corroborate that, as it was commented in previous sections, the best results are obtained with the mTBM-LA approach, followed closely by the mTBM-GA approach and far away by the TB approach.
• The mTBM-LA approach provides more accurate results than the mTBM-GA approach in all shafts, as expected.But the error difference is small (lower than 2% in all shafts of this example), indicating that both approaches provide very similar results.
• The mTBM-GA approach tends to slightly overestimate the stiffness of the shaft compared to the  mTBM-LA approach.The reason is that in steps where the stiffness of the shaft is very sensitive to the angles α I and α A (steps with a high ratio D/d, as gear-type steps in Tab.3), the locally optimized angles are usually lower than the globally optimized values.
• The TB approach demonstrates to overestimate significantly the stiffness of the stepped shafts, particularly when there are steps with a high ratio D/d (like the ones generated by gears).And this is the reason why this approach provides results far from reality.
• The mTBM-LA approach, that uses locally optimized angles, has the best results in all shafts, but   it slightly differs from the reference elastic curve in some cases.This is due to the simplicity of the proposed model and, consequently, the low number of degrees of freedom in the optimization process.For example, in shaft 3 (Fig. 17d), there are only 4 degrees of freedom (two steps and two angles per step), but they are enough to lower the relative error from 20.81% (TB approach) to 2.04%, as it is shown in Tab. 4.
In conclusion, it has been demonstrated that the proposed approaches constitute a much better approximation than the direct application of the Timoshenko beam theory to compute the elastic curve of a stepped shaft.Furthermore, it is clear that the mTBM-LA approach provides the highest accuracy compared to reference elastic curve of the shaft, but the mTBM-GA approach is a simpler method (where there is no need to interpolate the cone angles for each step) that provides results of similar accuracy.Thus, the mTBM-GA approach has the best ratio between accuracy and simplicity and constitutes a good approximation to the computation of the elastic curves of stepped shafts.

Conclusions
In this work, it has been demonstrated that the Timoshenko beam theory cannot always provide accurate results when computing the elastic curve of a stepped shaft.For this reason, a new model of the shaft based on the Timoshenko beam theory but modified with a linear (instead of abrupt) variation of the diameter in the steps of the shaft has been proposed.Each linear variation of diameter generates a conic segment that is defined by a cone angle (α).In the new model, two angles, α I and α A , have been considered to compute the deformations caused by bending moment and shear load.Then, the model has been optimized for a number of stepped shaft cases in a design space that is characterized by two parameters: D/d and b/d, being D, d and b the higher diameter, lower diameter and width of the step, respectively.In this optimization, the elastic curve obtained with the proposed model has been compared with the one obtained with a realistic finite element model of the stepped shaft.As a result, the optimal values of α I and α A have been obtained for each case of the design space.
Regarding the optimal values of α I , three different regions can be distinguished in the design space: a wide region where there is a small variation (low gradient) from 30 • to 38 • and two regions where there are high gradients up to 70 • .But it has been demonstrated that the stiffness of the shaft cases in the regions with high gradients have a very low sensitivity to the value of α I , so good approximations of the elastic curve can be obtained even when using values of α I different from the optimal values.In the case of α A , the optimal value is 90 • in practically all cases of the design space where b/d ≥ 1.When b/d ≤ 1 there is a

Figure 2 :
Figure 2: Distribution of σz in plane YZ of a stepped shaft.

1 2Figure 3 :
Figure 3: Actual geometry and proposed physical model of stepped shafts.

Figure 5 :
Figure 5: Definition of a tapered beam with a generic distributed load applied.

Figure 6 :
Figure 6: Defining parameters of the analyzed cases.

Figure 7 :
Figure 7: Study cases considered for optimization.

Figure 8 :
Figure 8: Tridimensional finite element model of a stepped shaft.

Figure 9 :
Figure 9: Areas considered when comparing the approximated elastic curve with the reference elastic curve.

igure 10 :
Local optimization results for angles α I and α A .
of GA approach (a) Relative error of LA approach

Figure 12 :
Figure 12: Relative error obtained with TB approach.

Figure 13 :
Figure 13: Influence of α I in different cases.

Figure 14 :
Figure 14: Algorithm for contact analysis of gear drives considering the flexibility of shafts.
starting from the values of D/d and b/d of the step.A summary of the obtained values is presented in Tab. 3. On the other hand, regarding the mTBM-GA model, the global values of α IG = 34.872• and α AG = 90 • have been used in all steps of all shafts.

Figure 17 :
Figure 17: Deflection curves of the of the gear box.

Table 1 :
Comparison of errors of mTBM-LA, mTBM-GA and TB approaches.

Table 2 :
Design data and theoretical working conditions of the gear drive.

Table 3 :
Cone angles used by mTBM-LA approach.