Shaft-hub press fit subjected to bending couples: Analytical evaluation of the shaft-hub detachment couple

Applied Mathematical Modelling 50 (2017) 135–160

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Shaft-hub press fit subjected to bending couples: Analytical evaluation of the shaft-hub detachment couple E. Radi a, L. Lanzoni b, A. Strozzi b, E. Bertocchi b,∗ a b

DISMI-Dipartimento di Scienze e Metodi dell’Ingegneria, Università di Modena e Reggio Emilia, 42122 Reggio Emilia, Italy DIEF-Dipartimento di Ingegneria Enzo Ferrari, Università di Modena e Reggio Emilia, 41125 Modena, Italy

a r t i c l e

i n f o

Article history: Received 28 October 2016 Revised 4 April 2017 Accepted 8 May 2017 Available online 13 May 2017 Keywords: Press-fit Shaft-hub incipient detachment Bending couple Beam model Winkler foundation Analytical solution

a b s t r a c t A mathematical modelling of a shaft-hub press-fit subjected to bending couples applied to the shaft extremities is developed, and the value of the bending couple inducing an undesired shaft-hub incipient detachment is analytically determined. The shaft-hub contact is modelled in terms of two elastic Timoshenko beams connected by a distributed elastic spring, whose stiffness is analytically evaluated. Two models of the distributed spring are considered. The first model expresses the combined deformability of both the shaft and the hub cross sections. The second model accounts for the stiffening effect exerted by the shaft portion protruding from the hub on the adjacent shaft part that is in contact with the hub, and, consequently, it assumes only a rigid body motion of the shaft cross section, thus neglecting its deformability. Based upon this beam-like model, the bending couple producing the incipient detachment between the shaft and the hub is theoretically determined in term of the shaft-hub geometry, of the initial shaft-hub interference, and of the elastic constants. Comparisons with selected Finite Element (FE) forecasts indicate that the first modelling produces an incipient detachment couple that appreciably overrates the FE forecasts, whereas the second modelling lowers the error down to technically acceptable predictions. © 2017 Elsevier Inc. All rights reserved.

1. Introduction This paper presents a theoretical elastic analysis of a shaft-hub press-fit in the presence of a bending couple C applied to the shaft extremities protruding from the hub, as shown in Fig. 1, and it focuses upon the determination of the bending couple that locally annihilates the initial shaft-hub press-fit, and begins the detachment between the shaft and the hub. The shaft loading by a bending couple in a shaft-hub press-fit is considered in the standards DIN 7190, see [1, Fig. 11.13] and [2, Section 4]. Typical applications of interference fits are to be found in railway rolling stock, in steam turbines, in built-up crankshafts, in drive and powertrain engineering, and in gear and bearing mounting. More generally, assembly by interference may be an effective method of attaching bored members to shafts.

∗

Corresponding author. E-mail addresses: [email protected] (E. Radi), [email protected] (L. Lanzoni), [email protected] (A. Strozzi), [email protected] (E. Bertocchi). http://dx.doi.org/10.1016/j.apm.2017.05.018 0307-904X/© 2017 Elsevier Inc. All rights reserved.

136

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Nomenclature a a(r) b b(θ ) c ci ci fi (θ ) gi (r) l ld p ppf qi r x Ai A0 A1 B0 C C0 D0 E E0 E1 Fi Fr (i) Fθ (i) F0 F1 G Gi Hi I Ii K L Mi Q Qi R S Ue V y z

α γ γ rθ δ ε εr εθ θ κ λ

shaft radius, hub inner radius function hub outer radius function constant constant coefficients function function hub half length shaft-hub detachment axial length contact pressure contact pressure due to the press-fit alone contact force radial coordinate coordinate area of the cross section constant constant constant bending couple constant constant Young’s modulus constant constant constant body force body force radial and circumferential body force radial and circumferential body force shear modulus constant constant diametral interference moment of inertia Winkler foundation coefficient hub length bending moment shear force shear force radius of the bore filleted edge variable elastic strain energy relative displacement vertical coordinate axial coordinate a/b parameter shear strain radial interference constant radial strain hoop strain angular coordinate Kolosov’s constant parameter

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

ν vi

σa σr σθ τ rz (i) τ θ z (i) τ rθ (i) ϕi χi ω  i ∇2

137

Poisson’s ratio deflection axial stress radial stress hoop stress shear stress shear stress shear stress rotation of the beam cross section shear coefficient root relative rotation Airy stress function Laplacian operator

The title problem falls within the general heading of the possible separation between two initially contacting surfaces, as a result of an externally applied load. The detachment between two mating surfaces often constitutes an undesired event in mechanics, since it promotes tribological damage, e.g. [3,4]. The evaluation of the loading that causes separation between two contacting bodies has been the subject of several papers. For instance, in [5,6] a pin-lug contact is considered, and the possible bush-lug separation is investigated. Similarly, in [7] the loosening mechanism is examined of the bush press-fitted in the small end of a connecting rod. In [8] a shrink-fit contact separation problem is considered in the thermo-plastic regime within the assumption of axisymmetry. A problem closely connected to separation is the undesired slip between two contacting bodies. As for the separation problem, the evaluation of the slippage incipient condition requires a detailed knowledge of the contact pressure distribution. For instance, in pressed-up crankshafts, relative slippage between the crankpin and the crankweb occurs if the applied torque exceeds a threshold. This slippage problem has been investigated essentially experimentally in [9], with a simplified model in [10], and it has been reconsidered in [11] by adopting a rigorous approach based upon bipolar coordinates. Going back to the title problem, the analysis of the stresses induced by a mechanical connection achieved via press-fitting has been the subject of a conspicuous number of studies spread over various decades. A perfunctory literature review on the stress field exclusively induced by the press-fit is presented in the following, whereas a more exhaustive bibliographic survey may be traced for a solid shaft in [12] and for its extension to a hollow shaft in [13]. The first useful result is the plane, axisymmetric solution of a shaft-hub press-fit, achieved in terms of the Lamé equations [14]. An inadequacy of the above plane, polar solution is its inability to describe the often relevant contact pressure peaks at the shaft-hub contact extremities. Such pressure bumps are extremely localized, the axial length on which each pressure peak insists being often in the region of the fillet radius R of the hub bore, Fig. 1 and [12]. This inadequacy has been overcome in [12,13], where, by extending the normalization developed in [15] for a pin-lug connection, a non dimensional factor has been proposed, useful for expressing the stress concentrations at the contact rounded extremities in terms of the radius of the hub bore rounded edge, of the shaft and hub radii, of the imposed press-fit interference, and of the elastic constants. In the previous analysis addressing the press-fit stress field, the interference-fit problem has been modelled as cylindrical, since the stress field is axisymmetric, but it generally varies along the shaft axial direction, as evidenced by the above discussed presence of the lateral pressure bumps. In this cylindrical modelling, the effect of friction is usually neglected

Fig. 1. Shaft-hub press-fit under bending.

138

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

for simplicity [12,16], since the coefficient of friction for steel is relatively low, e.g. [17], although it is known that friction may substantially alter the whole contact pressure [18], often mitigating the stress concentrations at the shaft-hub contact extremities [12]. The presence of torque is generally neglected, e.g. [17], p. 4, [19], since it does not alter the shaft-hub contact pressure distribution due to the press-fit alone. In addition, in the previously examined references addressing a cylindrical modelling of the press-fit, the commonly encountered situation in which the shaft undergoes bending has been ignored for simplicity. The consequences on the press-fit contact pressure of the application of a bending moment to the shaft are addressed in the following. The bending moment confers a curvature to the shaft, which thus exhibits a concave and a convex side, Fig. 3 below; it is conjectured in [20,21], analytically proved in [19], and numerically confirmed in [17], p. 35 and in [22] that the shaft curvature modifies the stress field with respect to its counterpart addressing the press-fit alone, by increasing (lowering) the lateral contact pressure peaks in the shaft concave (convex) side. (An opposite trend holds for the central contact pressure.) Consequently, the stress field is no longer axisymmetric, but it becomes three-dimensional, e.g. [23]; as a result of the shaft rotation (the bending moment is supposed not to rotate with the shaft), the stresses are no longer steady, but they are subjected to a fatigue cycle. The above three-dimensional problem remains linear to loads that are insufficient to cause separation between the shaft and the hub [24]. As the bending couple applied to the shaft extremities is increased, the lateral contact pressure peaks diminish until the shaft detaches from the hub, e.g. [23]. The occurrence of a shaft-hub partial detachment constitutes a symptom of poor mechanical design, since it promotes undesired phenomena of fretting fatigue [25]. The above mentioned determination of the loading causing shaft-hub detachment must account for the geometry of the press-fit, for the intensity of the applied loading, and for the degree of interference. Two main categories of press-fits subject to bending may be proposed, see [26, Figs. 6 and 10]. The first group collects problems in which the shaft laterally protrudes from both the hub sides, and the shaft protruding portions are subjected to bending, as in Fig. 1. Bearing mounting often belongs to this category. The second category embraces situations where the shaft protrudes from only one side of the hub, whereas the remaining shaft extremity is often aligned with one of the two hub lateral faces; in this case, the loading is applied to the shaft protrusion and to the hub. The crankpin press-fitting into a crankweb belongs to this second group. (The above proposed classification does not exhaust all the possible press-fit configurations. For instance, in [20, Fig. 3] the geometry is examined of a split shaft, fitted into the hub in correspondence of its separation.) Two markedly different stress fields are expected to take place for the two above categories. In the first group, the shaft cross sections along the press-fit axial length are subject to bending, whereas in the second group the bending moment is null at the shaft extremity aligned with the corresponding hub lateral face. This difference is expected to cause a variation in the distribution of the shaft-hub reaction force in the shaft axial direction, see Fig. 9 below, which in turn will modify the value of the loading initiating the shaft-hub separation. Comparison between the two experimental deflection curves 3 and 6 in diagram 3 of [20], referring to a split and to a continuous shaft, respectively, evidences the different mechanical response of the two fits. It may however be surmised on physical grounds that the detachment loading becomes comparable for the two groups if the hub axial length is very high. Papers examining the shaft-hub detachment and belonging to the first group include [22,23,26]. Investigations addressing the shaft-hub detachment and belonging to the second group comprise [17,19,20,26–28]. This paper exclusively considers the loading illustrated in Fig. 1, belonging to the first group. In addition to the effect of the macro-geometry of the press-fit, the geometrical detail of the fillet radius R of the hub bore edge, Fig. 1, is equally relevant in the determination of the contact pressure distribution and of the detachment loading. In fact, the numerical forecasts of [23] indicate that, if the hub bore edges are filleted, the couple that induces the shaft-hub detachment is higher than its counterpart for sharp hub bore edges; in practical situations, this increase may be of, say, 20 percent. A thorough modelling of the shaft-hub detachment problem should therefore incorporate the hub bore fillet radius effect. Papers addressing the shaft-hub separation problem are reviewed below, following a logical rather than a chronological order. It has previously been mentioned that, to simulate the bending couple effect, Fig. 1, the cylindrical model of the shafthub press-fit has to be abandoned in favour of a three-dimensional approach. In [20] a beam-like model has been employed to mimic such three-dimensional aspects. A simplified pressure distribution between the shaft and the hub, linear in the axial direction and cosinusoidal in the circumferential direction, has been assumed, and an analytical expression of the bending moment initiating the shaft-hub detachment has been derived. Friction effects have been included, but the above commented press-fit pressure peaks at the shaft-hub contact extremities have been neglected. This simplified modelling is analogous to the approximate analysis developed in [29] to redistribute the contact pressure between the stem and the head of a femoral prosthesis when the applied load becomes inclined. In [19] a considerably more advanced beam-like model has been developed, that accounts for both the shaft bending and the elasticity of the shaft and hub cross sections. For simplicity, both the shaft shear deformation and the hub flexibility have been neglected, having in mind applications to hollow shafts and to narrow press-fit widths. An elastic foundation has been introduced to model the deformability of the shaft-hub cross section. Finally, friction has been disregarded. A closedform solution for the shaft-hub reaction force has been achieved. The loading initiating the shaft-hub detachment has been analytically evaluated. In [22] a three dimensional FE analysis has been performed for a shaft-hub press-fit, where the hub bore edges have been assumed as sharp. Friction has been neglected. In [23] a 3D FE study of the possible shaft-hub detachment has been

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

139

performed in the presence of a bending couple, although in the absence of friction, and a preliminary design diagram expressing the normalized value of the couple applied to the shaft that initiates the shaft-hub detachment has been compiled for a selection of shaft-hub aspect ratios. The effect of the fillet radius of the hub bore edges has been explored. Derived from numerical interpolation, an approximate formula expressing the bending couple that begins the shaft-hub detachment has also been proposed for a limited range of shaft-hub geometries. In the comprehensive PhD thesis [17] a complete coverage of the state of the art is provided with regard to the salient design aspects of cylindrical and conical press-fits. The stress fatigue cycles, the fretting fatigue aspects, and the occurrence of a shaft-hub detachment, have been included. In particular, various literature formulae expressing the shaft-hub detachment loading have been ranked. The press-fits addressed in [17] belong to the above defined second group. Going back to the present paper, the geometry and loading of the press-fit described in Fig. 1 have been considered owing to their relative simplicity. A rigorous three-dimensional analytical elastic study of the title problem described in Fig. 1 being prohibitively complex, it was decided to develop the analytical beam-like solution of [19] by including the hub flexural and shear deformation, as well as the shaft shear deformation. In an effort to increase its accuracy, the beam model has been endowed with a foundation model, as in [19]. In detail, the interaction between the two Timoshenko beams describing the shaft and the hub has been mimicked by introducing a distributed linear spring between them, as in the Winkler foundation model [30]. In this paper the Winkler constant is exactly evaluated by considering the in-plane deformability of the shaft and hub cross sections, due to the shear load exchanged between them. To estimate the deformability of the shaft-hub contact, reference has been made to the plane model of an elastic circular section surrounded by an elastic annulus, both loaded by a body force distribution equal to the shear stress distribution caused by the shear force, according to the de Saint–Venant (dSV) beam theory [31,32]. A simple interpolating formula of the Winkler foundation coefficient is provided. Following [19], the detachment initiation couple is analytically determined by imposing that the compressive contact pressure at the hub extremities and caused by the press-fit alone, be equal and opposite to the tensile pressure due to the bending couple alone, where a tensile pressure is mathematically justified by the employment of a bilateral modelling of the shaft-hub contact. Three simplifications have been adopted: (a) a distributed linear spring has been employed, as in [19], although it is known that it is unable to correctly model the contact zones characterized by a high gradient of the reaction force, e.g. [33]; (b) friction has been disregarded, as in [19]; (c) the effect of the fillet radius of the hub bore edge has been ignored, because it cannot be modelled with the beam theory. The errors ascribable to the several approximations adopted−application of the beam model to non sufficiently slender bodies, plane analytical modelling of the Winkler foundation constant, employment of an axially constant contact pressure distribution according to the plane axisymmetric press-fit solution, sinusoidal distribution of the contact pressure in the circumferential direction – are cumulatively quantified by comparing the analytical solution developed in this paper to selected FE forecasts. As a result of the above approximations, it is expected that the beam-like approach favoured in this study may only predict an indicative value of the bending couple initiating the shaft-hub detachment rather than providing its accurate threshold. It is however believed that the knowledge of the indicative value of the detachment bending couple for this complex contact problem still constitutes a useful information from a technical viewpoint. Extensions of the solution developed in this paper to more complex loadings are possible. Aspects dealing with tribological damage, fretting fatigue, and fatigue notch factor, fall beyond the scope of the present paper. The organization of this paper is as follows. In Section 2 the press-fit problem between the shaft and the hub is addressed. In Section 3 the model is developed of two Timoshenko beams, simulating the shaft and the hub, respectively, connected by a continuous distribution of elastic springs according to the Winkler foundation model, and loaded by two opposite couples applied to the shaft extremities. The spring stiffness, namely the Winkler foundation constant, is then analytically evaluated in Section 4. The bending couple that initiates the shaft-hub detachment is determined in Section 5, and the analytical forecasts are compared to available FE predictions and to an approximate model in Section 6. In Section 7 an improved solution is developed, that is based upon an alternative evaluation of the Winkler constant. Finally, in Section 8 an assessment is proposed on the applicability of this elastic approach.

2. Press-fit contact pressure Fig. 1 displays a shaft-hub press-fit in the presence of a bending couple C applied to the shaft extremities protruding from the hub, and it clarifies the meaning of various symbols. The inner and outer radii are denoted by a and b, and the hub length is L, whereas the hub semi-length is l. The radius of the filleted edge of the hub bore is R. The applied couple is C. The axial, z, vertical, y, radial, r, and angular, θ , coordinates are also illustrated. The Lamé plane axisymmetric (polar) contact pressure due to the press fit alone, ppf (the index pf denotes “press-fit”), in the shaft-hub central contact zone sufficiently far from the hub lateral extremities, and in the absence of a bending couple, is [14]:

ppf =

 Eδ  1 − α2 2a

(1)

140

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Fig. 2. Normalized contact pressure versus normalized axial coordinate from the hub sides, for α = 0.5, a/L = 1/4, δ = 0.01 a and for the normalised fillet radius R/a = 0.05.

where a denotes the nominal value of both the shaft radius (the shaft is assumed to be solid) and the hub inner radius, b represents the hub outer radius, Fig. 1, and α = a/b is the ratio between the inner and outer radii of the hub. The shaft and hub are modelled as linear elastic and isotropic bodies, where E denotes the Young’s modulus of the materials of both the shaft and hub, and δ is the radial interference. It has been observed in the Introduction that, in the absence of a bending couple, the actual axisymmetric contact pressure profile in the axial direction is camel backed, that is, it exhibits a flattish central part, whose value is correctly forecast by formula (1), and two lateral pressure bumps, occurring where the hub bore edges are rounded [12]. Unfortunately, the plane solution of formula (1) does not mimic such lateral pressure bumps. (If the hub bore edges are sharp and not rounded, such elastic pressure peaks become unbounded. If the shaft extremity is aligned with the hub lateral face, the pressure peak is absent, e.g. [19].) It is also known that the shaft-hub detachment due to the bending couple C initiates from the hub extremities e.g. [23], Fig. 3 below, where the contact pressure peaks in the absence of the couple C are considerably higher than the pressure central values. Consequently, at first sight it seems essential to account for such pressure peaks in the theory under development, the Lamé plane axisymmetric modelling of formula (1) appearing totally inadequate. Instead, the following observations show that the presence of the lateral pressure spikes may be neglected in an approximate determination of the detachment couple, and that the Lamé plane axisymmetric formula may be appropriate. In fact, a specific FE analysis has been carried out to clarify the distribution of the contact force when the shaft has started to move apart from the hub. Fig. 2 depicts the FE shaft-hub contact pressure in the axial direction z for a hub with filleted bore. More precisely, the contact pressure p is reported for a value of the bending couple C that has just begun the shaft-hub detachment along the flat portion of the hub bore axial profile, see Fig. 3(a) below. The contact pressure refers to an angular position for which the applied bending couple C causes the initiation of the shaft-hub separation. This angular position corresponds to the shaft maximum tensile bending stress; Fig. 1 clarifies that this angular position is defined by θ = 0. As a result of the press-fit alone, i.e., for C = 0, only a limited portion of the rounded edge of the hub bore axial profile is in contact with the shaft; in this zone the contact pressure exhibits pronounced bumps. Consequently, for C = 0 the contact pressure is null along a normalized axial length ld /R just lower than unity (ld is measured from the hub lateral surfaces, see Fig. 3 below). When, as a consequence of the application of the couple C, the shaft begins to detach from the hub along the initially flat portion of the hub bore axial profile, the contact pressure along the shaft axial direction is null along the normalized axial length ld /R = 1, and, more important, it no longer exhibits lateral pressure bumps, but it smoothly decays from its flattish central value to zero (Fig. 2). In addition, in [23, Fig. 6] the FE curve expressing the normalized axial extent of the detachment zone in terms of the normalized surplus of the applied couple C beyond its critical value that initiates the shaft-hub detachment in a sharply edged hub, is reported for a sharply edged hub bore, and for two differently rounded edge bores. Two aspects emerge from the above quoted figure. First, when the detachment axial extent is sufficiently higher than the fillet radius, i.e., ld >> R, the detachment extent for rounded bore edges, Fig. 3a, is very similar to that referring to a sharply edged hub bore (Fig. 3(b)). In other words, provided that ld >> R, the presence of rounded edges moderately affects the axial detachment extent for an imposed value of the bending couple beyond its detachment value. Secondly, a comparison may be carried out between the values of the bending couple that initiates the detachment between the shaft and the hub,

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

141

Fig. 3. Detachment axial length ld for a hub bore rounded edge and for a sharp edge.

for sharp and for rounded edges of the hub bore, respectively. It appears from [23, Fig. 6] that the value of the bending couple that, for filleted hub bore edges, begins the detachment of the initially flat portion of the axial profile of the hub bore, is perceivably higher than its counterpart valid for sharp edges. For the two practically relevant fillet radii R defined by the ratios R/a = 0.025 and R/a = 0.05, the detachment couple is higher than its sharp edge analogue by about 10 and 20 percent, respectively. From the previous discussion it may be concluded that (a) in computing the bending couple that initiates the shafthub detachment, it is acceptable to refer to the flattish central contact pressure of formula (1), and to neglect the lateral pressure bumps; (b) the bending couple that initiates the shaft-hub detachment in the presence of filleted hub bore edges is higher than that referring to sharp edges, but it is of the same order of magnitude; therefore, the presence of rounded edges may be neglected for simplicity in the evaluation of an indicative value of the bending couple beginning the shaft-hub detachment, rather than of its accurate threshold, and reference may be made to sharp edges. In formula (24) of [19] and in formula (7) of [20] the presence of the lateral pressure peaks has been disregarded.

3. Contact pressure due to the bending couple Following [19], the frictionless contact pressure due to the bending couple alone, and acting between the shaft and the hub, is examined in the presence of a perfect fit, i.e., in the absence of interference or clearance. The main aim is to evaluate the contact pressure profile especially at the shaft-hub contact extremities, where the shaft begins to detach from the hub (Fig. 3). A bilateral modelling of the shaft-hub contact is postulated, an assumption, this, which mathematically justifies the possible occurrence of a tensile contact pressure. (The total contact pressure, evaluated as the sum of the pressure due to the press-fit, and of its counterpart imputable to bending alone, must be positive.) Both the shaft and the hub are described in terms of the Timoshenko beam theory; various tests have shown that the introduction of the shear deformation produces a more realistic contact pressure between the shaft and the hub, see Fig. 7 below. With this idealization, the title problem is described by a fourth order differential equation with constant coefficients, whose solution may be expressed in terms of exponential functions [31,34]. The bilateral interaction between the two components is modelled by a continuous distribution of linear elastic springs, whose stiffness is determined by analysing the deformation of the shaft and hub cross sections under the action of two equal and opposite shear forces exchanged between these two components; the shear stress distribution over their cross section is provided by the classical solution of the dSV problem [35]. From a mathematical standpoint, the introduction of a Winkler foundation model in the description of the contact pressure is physically advantageous, since it avoids the outcome of stress singularities at the contact extremities [35,36], which would preclude the possibility of a comparison with the bounded Lamé contact pressure (1) due to the press-fit alone.

3.1. Governing equations A brief description is reported in the following of the mathematical modelling of the frictionless shaft-hub contact problem in terms of the two Timoshenko beams describing the hub and the shaft. The bending couple C is applied to the shaft extremities, Figs. 1 and 4.

142

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Fig. 4. The beam model of the shaft-hub contact problem.

Fig. 5. Cross sections of the hub (0) and shaft (1) under the shear loading Q.

The hub is denoted by index 0, and the shaft by index 1, Fig. 5 below. Then, the governing equations for the Timoshenko beams (0) and (1) are:

d Qi = −qi , dz d ϕi Mi = , dz E Ii

d Mi = Qi , dz dv χQ ϕi + i = i i , dz GAi

(2)

( i = 1, 0 )

where y and z are the vertical and axial variables, respectively, Fig. 1, Mi denotes the bending moment, Qi is the shear force, vi designates the deflection of the axes of the shaft and hub along the y-axis direction, ϕ i is the rotation of the cross sections of the beam, positive if anticlockwise, and qi is the distributed load applied along the beam resulting from the contact interaction between shaft and hub, positive if applied in the y-axis direction. In addition, Ai and Ii denote the area and the moment of inertia of the cross sections of the two beams, respectively, G is the shear modulus, and χ i is the shear

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

coefficient, namely [37,38]

χ0 =

7 + 34 α 2 + 7 α 4 + 6 (1 +



ν

2

1+ν 2 2

2

(1 − α 2 )

α )

,

χ1 =

7 + 14 ν + 8 ν 2 6 (1 + ν )2

≈ 1.176,

143

(3)

where ν is the Poisson’s ratio (note that χ 1 can be obtained as the limit of χ 0 for α → 0). The numerical value in Eq. (3) refers to ν = 0.3. Note that continuity of tractions exchanged between the shaft and the hub requires

q1 + q0 = 0.

(4)

We assume a linear relation between the distributed load and the relative displacement between the shaft and the hub along the y-axis direction, as a result of the cross section deformability. In fact, only incipient separation between the shaft and the hub is investigated, and a relative detachment, which would involve nonlinearity, see the Introduction, is not modelled in this paper. For simplicity, the contact between the shaft and the hub is described as frictionless, see the comments in the Introduction and in [19]. By introducing the relative displacement V in the y direction, and the relative rotation 

V = v1 − v0 ,

 = ϕ1 − ϕ0 ,

(5)

the relation between the distributed contact load, q0 and q1 , and the total relative displacements between shaft and hub, resulting from the cross section deformability due to the shear stress and the normal bending stress arising in the two components, is expressed as

q0 = −q1 = K V,

(6)

where the balance condition (4) has been considered, and K is the Winkler constant, which will be evaluated in Section 4 in terms of the cross section geometry and elastic constants. From Eqs. (2) and (6) one obtains:

d Q0 = −K V, dz

d Q1 = K V. dz

(7)

By calculating the difference between Eq. (2) evaluated for i = 1 and i = 0, respectively, one obtains:

dV χ1 Q1 χ0 Q0 = − , dz GA1 GA0

+

M0 d M1 = − , dz E I1 E I0

(8)

Q0 d2  Q1 = − . E I1 E I0 d z2 By differentiating once the first and the last Eq. (8) and using Eq. (7), two coupled ODEs for V and  may be derived:

d d 2V 4γ 2 + = 2 V, dz d z2 l d3  4λ4 = 4 V, d z3 l

(9)

where the following non-dimensional parameters have been introduced

λ4 =





K l4 1 1 + , 4 E I1 E I0

γ2 =





K l 2 χ1 χ0 + . 4 GA1 GA0

(10)

From Eq. (9)2 , the expression of V in terms of  is determined:

V =

l 4 d3  . 4λ4 dz3

(11)

In addition, the above Eq. (9) may combined to obtain the following ordinary differential equation (ODE) in the unknown function  alone:

4ε λ2 d2  d4  4λ4 − + 4  = 0, 4 2 2 dz l dz l

(12)

where the non-dimensional constant ε has been introduced

   γ 2  χ1 χ0  K E I0 I1 1+ν K 1 + α2  ε= 2 = + = (1 − α 2 ) χ1 + α 2 χ0 . 2 GA1 GA0 4 (I0 + I1 ) 2 πE 1−α λ

(13)

144

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Since the Winkler foundation coefficient K is proportional to E, then ε is a purely geometrical parameter independent of the Young’s modulus. A different combination of Eq. (9) would produce an alternative formulation of the governing Eq. (12) in terms of the relative displacement V, which is physically more sound in this shaft-hub detachment problem. However, the boundary conditions are more naturally expressed in terms of the relative rotation , and, consequently, the formulation of the governing Eq. (12) in terms of  has been preferred. The solution of the fourth order linear ODE (12) may be expressed in terms of the sum of exponential functions [19,31]. To this aim, it is set:

(z ) = c eω z/l

(14)

where c is a constant. The fourth degree characteristic equation becomes:

ω4 − 4ε λ2 ω2 + 4λ4 = 0

(15)

whose four solutions are:

  ω = ±λ 1 + ε ± i 1 − ε ,

(16)

where the expressions for λ and ε are given in Eqs. (10) and (13), respectively, and i = Then, the general solution of the ODE (12) in terms of  becomes:

(z ) =

√

−1.

      z z z λ 1 − ε + c2 cos λ 1 − ε cosh λ 1 + ε l l l        z z z + c3 sin λ 1 − ε + c4 cos λ 1 − ε sinh λ 1 + ε , Cl E I1



c1 sin

l

l

(17)

l

where C is the couple applied to the shaft extremities (Figs. 1 and 4) and ci (i = 1, 2, 3, 4) are non-dimensional constants to be determined by imposing the boundary conditions at z = 0, l. 3.2. Imposition of the boundary conditions The unknown constants ci of (17) are evaluated by imposing four boundary conditions derived from Fig. 1. It is required that, due to symmetry, the rotation and the shear force vanish at the shaft-hub contact midpoint, namely at z = 0, for both beams:

ϕ0 ( 0 ) = ϕ1 ( 0 ) = 0,

d 2 ϕ0 ( 0 ) d 2 ϕ1 ( 0 ) = = 0. 2 dz d z2

(18)

The bending moment at the hub extremity, namely at z = l, is set equal to the applied moment C for the shaft, whereas it is null for the hub. Moreover, the shear force must vanish therein for both beams:

d ϕ1 ( l ) C = , dz E I1

d ϕ0 ( l ) = 0, dz

d 2 ϕ1 ( l ) = 0, d z2

d 2 ϕ0 ( l ) = 0. d z2

(19)

The boundary conditions (18) and (19) provide the four conditions for the function (z)

d2  ( 0 ) = 0, d z2

(0 ) = 0,

d C (l ) = , dz E I1

d2  ( l ) = 0. d z2

(20)

By imposing the symmetry conditions, namely the first and second boundary conditions (20) to expression (17) for , the coefficients c2 and c3 of expression (17) are null, and the solution in terms of  simplifies to

(z ) =



Cl c1 sin EI1

        z z z z λ 1 − ε cosh λ 1 + ε + c4 cos λ 1 − ε sinh λ 1 + ε . l

l

l

l

(21)

Then, the third and fourth boundary conditions (20) supply the following expressions for the constants c1 and c4

 √   √   √   √  √ 1 − ε 2 sin λ 1 − ε cosh λ 1 + ε − ε cos λ 1 − ε sinh λ 1 + ε  √  √  √  c1 = , √ λ 1 + ε sin 2λ 1 − ε + 1 − ε sinh 2λ 1 + ε  √   √   √   √  √ 2 1 − ε 2 cos λ 1 − ε sinh λ 1 + ε + ε sin λ 1 − ε cosh λ 1 + ε  √  √  √  c4 = . √ λ 1 + ε sin 2λ 1 − ε + 1 − ε sinh 2λ 1 + ε 2

(22)

(23)

The introduction of constants (22) and (23) into Eq. (17) provides the exact solution for the function . Then, the function V can be determined from  by using Eq. (11). According to Eq. (6), the distributed load attains the following magnitude at the hub extremity

q0 (l ) = KV (l ) =

K l 4 d3 (l ) Cl = K l S(λ, ε ), E I1 4λ4 dz3

(24)

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

145

where the function S(λ,ε ) is evaluated by differentiating  of Eq. (21) three times with respect to z, and by computing its value for z = l. One obtains

S(λ, ε ) =

√ √ √ √ 1 1 + ε sin(2λ 1 − ε ) − 1 − ε sinh(2λ 1 + ε ) . √ √ √ √ 2λ2 1 + ε sin(2λ 1 − ε ) + 1 − ε sinh(2λ 1 + ε )

(25)

In Fig. 1 and in [19, Formula (19)], boundary conditions different from those expressed in Fig. 1 of the present paper are employed. Such boundary conditions describe a shaft-hub press-fit in which the shaft protrusion from the hub is subject to a couple and to a shear force, and where the remaining shaft extremity is aligned with one of the web lateral faces; the equilibrating loads are applied to the hub. These boundary conditions may be applied to solution (17) of the present paper. 4. Evaluation of the Winkler foundation coefficient It has been observed in Section 3 that the introduction of a Winkler foundation model in the description of the contact pressure between the shaft and the hub is physically advantageous, since it avoids the outcome of stress singularities at the contact extremities. In this section the available information on the evaluation of the Winkler constant, see formula (6), is reviewed, and an improved model is developed. 4.1. Available Winkler foundation approaches In [19] the deformability of the shaft-hub cross section has been estimated by adopting a plane model formed by a solid disk surrounded by a ring (Fig. 5). The variation of the shear force in the axial direction produces a corresponding variation of the shear stresses distributed over the lateral faces of the disk and ring, whose effect is accounted for by applying to the fit cylindrical surface a sinusoidal shear stress τ rθ . A closed form solution is obtained for the elastic stress and displacement field. The Winkler constant K is assumed in [19] as the ratio between the contact force per unit length and the variation of the two segments AB and CO of Fig. 5. For the reference geometry of the shaft-hub cross section defined by α = 1/2, the normalized Winkler constant K/(2 G) is 7.29. Alternatively, by referring to the variation of the vertical distance between the two points D and O, K/(2 G) becomes 2.51. In [24] the compliance of the circular cross section (not surrounded by a ring) of a beam is analytically examined. It is noted that, since the circular slice representing the shaft cross is taken to be thin, the differences between the values of τ rz and τ θ z acting on the two plane faces of the slice may be replaced by body forces. Since the exact distribution of the shear stresses τ rz and τ θ z due to a shear force is analytically known for both a solid disk and a ring, [38, p. 335], the above model has here been developed for a solid disk surrounded by a ring, describing a shaft-hub cross section, and the Winkler constant has been evaluated with an energy approach. This model is presented in detail in the following section. For the above reference geometry α = 1/2, the normalized Winkler constant K/(2 G) is 1.90. In Section 7 it is argued that the shaft cross section may be treated as rigid. With this assumption, K/(2 G) becomes 2.31. This increase is of about 20 percent. Alternatively, the body forces representative of the shear stresses may be assumed for simplicity to be uniformly distributed according to a gravitational field within both the solid disk and the ring, see Fig. 2(a) of [24,39], p. 94, and [40]. Details of the corresponding calculations are omitted for brevity. With this model, K/(2 G) is 3.68. If the shaft cross section is treated as rigid, then K/(2 G) becomes 4.81. This increase is of about 30 percent. The above numerical values evidence that the Winkler constant noticeably varies according to the approach employed. The present authors believe that the evaluation of the Winkler constant based upon the exact distribution of the shear stresses, developed in the following section, provides the most reliable value of the Winkler constant, but it is equally felt that only a three-dimensional accurate assessment of the various Winkler models may classify their reliability. Such critical comparison is beyond the scope of this paper. 4.2. Accurate Winkler foundation model In the following, the Winkler constant K appearing in Eq. (6) is determined for a section describing the shaft-hub pressfit, composed by a ring surrounding a sliding solid disk (Fig. 5) and subjected to two equal and opposite shear forces Q0 and Q1 , being Q0 +Q1 = 0 due to balance conditions. To this aim, the deformation of the shaft-hub cross section due to the shear stresses induced by the shear forces Q0 and Q1 is evaluated, and the constant K expressing the relative displacement between the two cross sections due to the distributed load is analytically computed using the Clapeyron theorem. The elastic springs connecting the two components are modelled as linear (i.e., the Winkler coefficient is assumed to be independent of the deflection), since the ring-disk interface is supposed to behave bilaterally due to the press-fit pressure, and only incipient separation between the shaft and the hub is investigated, whereas a detachment of finite extent, which would involve nonlinearity, is not considered. 4.2.1. Plane problem Here we investigate the deformability of the two cross sections, originating from the shear stress distribution and from the corresponding in-plane stress field arising from the frictionless contact condition. Following [24], we consider the plane

146

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

strain problem for both cross sections, subjected to body forces displaying the same distribution of the shear stresses provided by the classical solution of the shear bending problem for a beam with solid circular (1) and annular (0) cross section (Fig. 5). According to [37,39,41], the shear stress distribution in an annular cross section subjected to the upwards shear force Q0 acting along the y-axis direction is:



τrz

Q0 3 + 2ν = 8I0 1 + ν

τθ(z0)

Q0 3 + 2ν =− 8I0 1 + ν

(0 )

a2 b2 a + b − 2 − r2 r 2





2

cos θ ,

1 − 2ν 2 a2 b2 a +b + 2 − r 3 + 2ν r 2

2



(26) sin θ .

In Eq. (26), I0 is the moment of inertia of the ring, and a and b are its inner and outer radii, respectively. In addition, the angular coordinate θ is taken from the direction of Q0 (Fig. 5). Similarly, the shear stress distribution in a solid disk of radius a (the solid disk radius nominally coincides with the ring inner radius), subjected to the shear force Q1 , is:

Q1 3 + 2ν 2 (a − r2 ) cos θ , 8I1 1 + ν   Q 3 + 2ν 2 1 − 2ν 2 τθ(z1) = − 1 a − r sin θ , 8I1 1 + ν 3 + 2ν

τrz(1) =

(27)

where I1 is the moment of inertia of the disk. 4.2.2. Airy stress function formulation The three-dimensional equilibrium equations for the stress components in polar coordinates write

∂ σr 1 ∂ τrθ σr − σθ ∂ τrz + + + = 0, ∂r r ∂θ r ∂z ∂ τr θ 1 ∂ σ θ τ ∂ τθ z + + 2 rθ + = 0. ∂r r ∂θ r ∂z

(28)

Following [24], the last terms in Eq. (28) may be interpreted as the radial and circumferential components of body forces F( 0) and F( 1) acting on the ring and disk, respectively. Therefore, in the following, the determination is carried out of the in-plane stress components σ r , σ θ , and τ r θ induced by the body forces

Fr(i ) =

∂τθ(zi) ∂τrz(i) , Fθ(i ) = , ( i = 0, 1 ) ∂z ∂z

(29)

whose expressions in terms of q0 and q1 can be derived from (26) and (27) taking into consideration the first relation in (2). The in-plane stress components and the body force components must fulfil the following compatibility Eq. [42]

∇ ( σr 2

4 + σθ ) = − 1+κ



 ∂ Fr Fr 1 ∂ Fθ + + , ∂r r r ∂θ

(30)

where the Laplacian operator ∇ 2 in plane polar coordinates is



∇

2

1 ∂ ∂ = r r ∂r ∂r



+

1 r2

∂2 , ∂θ2

(31)

and where κ = 3 −4ν for plane strain, whereas κ = (3 −ν )/(1 +ν ) for plane stress [39]. It is noted that the body force components Fr and Fθ are not conservative (i.e. (curlF)z = 0) [39] since (i ) ∂ Fr(i) ∂ (r Fθ ) = (i = 0, 1). ∂θ ∂r

(32)

Following [43], two potentials φ and ψ are introduced for representing the non-conservative, non uniform body force distribution, such that

Fr = −

∂φ φ − ψ − , ∂r r

Fθ = −

1 ∂ψ . r ∂θ

(33)

According to (33), the integral representations of the two potentials φ and ψ in terms of the body force components Fr and Fθ are

ψ = −r φ=−



Fθ dθ + a(r ),     1 r

r Fr +



Fθ dθ dr −





a(r ) dr +

b( θ ) , r

(34)

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

147

where a(r) and b(θ ) are two unknown functions that may be neglected, since they provide vanishing body forces. The introduction of the Airy stress function  [39] allows the in-plane stress components to be written in the form [19,39,42,43]

σr =

1 ∂ 1 + 2 r ∂r r

∂ 2 + φ, ∂θ2

∂ 2 + ψ, ∂ r2   ∂ 1 ∂ τr θ = − . ∂ r r ∂θ σθ =

(35)

The above expressions (35) for the in-plane stresses fulfil the equilibrium Eq. (28), provided that the body force components Fr and Fθ are expressed in terms of the two potentials φ and ψ according to Eq. (33). The formulation of the compatibility Eq. (30) in terms of the Airy stress function  and of the two potentials φ and ψ is considered in the sequel. Based upon Eq. (33), the following relation is obtained

1 ∂ 2ψ ∂ Fr Fr 1 ∂ Fθ ∂ 2 φ 1 ∂ ( ψ − 2φ ) + + =− 2 − + . 2 ∂r r r ∂θ r ∂r r ∂θ ∂ r2

(36)

In addition, from Eq. (35), the following equality holds

∇ 2 (σr + σ θ ) = ∇ 4  + ∇ 2 (φ + ψ ).

(37)

Consequently, the compatibility Eq. (30) becomes

∇ 4 =

4 1+κ



 ∂ 2 ψ ∂ 2 φ 1 ∂ ( ψ − 2φ ) + − − ∇ 2 (φ + ψ ). r ∂r ∂θ2 ∂ r2

1 r2

(38)

The body forces Fr and Fθ of Eq. (29) are expressed in terms of the shear stresses τ rz and τ θ z of Eqs. (26) and (27); they are then introduced in Eq. (34) to provide the following expressions of the potentials φ and ψ for the rings





φ0 =

r3 q0 3 + 2ν (a2 + b2 ) r − 8I0 1 + ν 3 + 2ν

ψ0 =

1 − 2ν 3 q0 3 + 2ν a2 b2 (a2 + b2 ) r + − r 8I0 1 + ν r 3 + 2ν

cos θ ,



and for the disk



φ1 ψ1

1 − 2ν 3 q1 3 + 2ν 2 = a r− r 8I1 1 + ν 3 + 2ν

cos θ ,

2



(39) cos θ ,



q1 3 + 2ν = 8I1 1 + ν

r3 a r− 3 + 2ν





(40) cos θ .

Note that expressions (39) and (40) already provide the correct expressions of the body forces (33) according to Eq. (34), without needing to introduce integration constants, which would provide vanishing body forces and, in turn, no effects on the stress field. From Eqs. (39) and (40) the following relations useful for evaluating the compatibility Eq. (38) are derived:

∇ 2 ( φi + ψ i ) = −

2 ( 1 − ν ) qi r cos θ 1 + ν Ii

( i = 0, 1 ),

(41)

and

1 r2

qi ∂ 2 ψ i ∂ 2 φi 1 ∂ ( ψ i − 2 φi ) + − = − r cos θ (i = 0, 1 ). r ∂r Ii ∂θ2 ∂ r2

(42)

With the aid of Eqs. (41) and (42), the compatibility Eq. (38) for the ring and the disk becomes

∇ 4 i = −2

 2 1+κ

−

1−ν 1+ν

q

i

Ii

r cos θ ,

( i = 0, 1 ).

(43)

148

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

The integration of Eq. (43) yields the following Airy stress functions 0 for the ring and 1 for the disk, respectively

q 0 = 0 I0





r5 1−ν r 1 2 A0 r + B0 r log + C0 − − a r 96 1 + κ 1+ν 3





r5 1−ν q 2 1 = 1 A 1 r 3 − − I1 96 1 + κ 1+ν







cos θ + D0 r θ sin θ , (44)

cos θ ,

where A0 , B0 , C0 , D0 and A1 are constants to be determined by imposing the boundary conditions at r = a, b; in addition, the result

∇ 4 (r5 cos θ ) = 192 r cos θ ,

(45)

has been used. Comparison between the expressions of 0 and 1 of Eq. (44) and the expression of φ of [19, Eq. (5)] clarifies the degree of approximation of [19]. 4.2.3. In-plane stress, strain and displacement fields The in-plane stress components both in the disk and ring are computed from Eqs. (35) and (44) in terms of the unknown constants. The stress field in the ring (0) becomes



σr

(0 )

 3 + 2ν r 3 2 C0 q0 B0 + 2 D0 r 2 = 2A0 r + − 3 + a + b2 − I0 r 8 1+ν 24 r 

τr θ ( 0 ) = − σθ(0) = −

q0 2 C0 r3 B0 2A0 r + − 3 − I0 r 24 r



q0 5 r3 B0 2 C0 6A0 r + + 3 − I0 r 24 r

 2

1+κ

 2

1+κ

−

1−ν 1+ν

−

1−ν 1+ν



sin



+

 2 1+κ

2+ν + 1+ν

 cos

θ

θ 

3 + 2ν 1 − 2ν a2 b2 + r 2 (a2 + b2 ) − r 4 8 r (1 + ν ) 3 + 2ν

 cos

θ, (46)

and stress field in the disk (1) becomes







q1 r 3 + 2ν 2 2+ν 48A1 + 3a2 − r2 + cos θ 24 I1 1+ν 1+κ 1+ν    1−ν q r 2 − sin θ τrθ (1) = 1 48A1 − r2 24 I1 1+κ 1+ν    2+ν q r 3 + 2ν 10 − r2 − cos θ . σθ(1) = 1 144A1 + 3a2 24 I1 1+ν 1+κ 1+ν

σr ( 1 ) =

(47)

The vanishing shear stress τ r θ (1) at the border of the disk, i.e. at r = a, as required by the frictionless contact condition, provides the constant





1−ν a2 2 − . 48 1 + κ 1+ν

A1 =

(48)

The displacement fields in the disk and ring follow from integration of the strain components provided by the linear elastic constitutive relation under plane stress or plane strain conditions [39, p. 43]

 1  (1 + κ ) σr(i) − (3 − κ ) σθ(i) , 8G

εr(i) =

εθ(i) =

 1  (1 + κ ) σθ(i) − (3 − κ ) σr(i) , 8G

γr(θi) =

τr(θi) G

,

(49)

where G is the elastic shear modulus, and the expressions of κ for plane stress and plane strain have been reported after formula (31). Then, the displacements ur ( i ) and uθ ( i ) for both the disk and the ring are computed by integration [39], as

ur(i ) =



εr(i) dr +

d f i (θ ) , dθ

uθ(i ) =



(r εθ(i) − ur(i) ) dθ + gi (r ),

(50)

where fi (θ ) and gi (r) are two unknown functions that are determined by employing the definition of the shear strain γr(θi ) [39], namely from:

τr(θi) G

=

1 ∂ ur(i ) + r ∂θ

∂ uθ(i) uθ(i) − . ∂r r

(51)

The following equation in the unknown functions fi (θ ) and gi (r) is obtained

−r

d gi ( r ) d 2 f i (θ ) + gi ( r ) = + fi (θ ). dr dθ 2

(52)

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

149

Therefore, both the right-hand and left-hand sides must be constant. The expressions of f1 (θ ) and g1 (r) thus are

fi (θ ) = Ei sin θ + Fi cos θ + Hi ,

gi (r ) = Gi r + Hi ,

(53)

where the constants Fi , Gi , and Hi are set to zero to preserve symmetry and to ignore a rigid body rotational motion [39]. Note that the constant Ei is instead responsible for an admissible rigid body translation along the direction of the shear force. Correspondingly, the radial displacement in the ring is

ur(0) =



B q0 C D r4 (κ − 2 )A0 r2 − 0 (1 − κ ) ln r + 20 + 0 (1 + κ ) ln r + 2 G I0 2 2 96 r  3 + 2ν  2 2 2 2 2 − r (a + b )(1 − κ ) + a b (3 − κ ) ln r cos θ , 32 (1 + ν )

and the radial displacement in the disk is

ur(1) =

  10

q1 r4 192 G I1

1+κ

−

5 + 4ν 1+ν





+ r 2 a2 8 κ − 6

3+κ 5−κ + 1+κ 1+ν



 10 1+κ

−

5 + 4ν 1+ν

 (54)



− E1 cos θ ,

(55)

where E1 is an additional constant responsible for a relative rigid body translation along the direction of the shear force, and E0 has been set equal to 0 without loss of generality. Six integration constants occur in Eqs. (47), (48), (54) and (55). The constant A1 has already been computed in Eq. (49). The four constants A0 , B0 , C0 , D0 , are evaluated by imposing the vanishing of the shear stress and the continuity of radial stress along the contact region between the disk and the ring, i.e. at r = a, and the vanishing of the shear and radial stress at the outer boundary of the ring, i.e. at r = b, together with the accomplishment of Eq. (52):

A0 = B0 =



b2 (α 4 + α 2 + 1 ) 48 (1 + α 2 ) 4

b 8

κ −1 κ +1



2α 2 a2 b4 C0 = 3 − + 1+ν 48(1 + α 2 ) D0 = −





 2

1+κ

−

1−ν 1+ν

α 2 (3 + κ ) − 6 1+κ



+3

1−κ 1+κ





(56)

a2 b2 1 2 2+ + 2 . 16 1+ν α

Finally, the remaining constant E1 is determined by imposing the continuity of the radial displacement between the shaft and the ring at r = a:

E1 =



  κ − 1 − 4 log α 2(7 + 12κ log α ) 4 ( 3 + 2 ν ) 3 + 4 ( 2κ − 1 ) − + α 1+ν 1+κ (1 + α 2 )(1 − α 4 )   46 12[3 + 2ν + κ (5 + 4ν )] 4−κ 5 − 2κ + + α 2 10(1 + κ ) − + − log α . 1+ν 1+κ 1+ν (1 + ν )(1 + κ ) a4

(57)

4.2.4. Determination of the Winkler constant The in-plane stresses and strains allow the elastic strain energy to be evaluated and, consequently, the Winkler constant K to be quantified by using the Clapeyron theorem, namely by comparing the elastic strain energy to the work made by the two opposite shear forces Q. To this aim, we calculate the elastic strain energy stored in both components, written in terms of stresses 1  1  Ue = 16 G Ai i=0



 2  2   2  i (i ) (i ) (i ) (i ) + σθ + 2 (κ − 3 ) σr σθ + 8 τr(θ ) dA . ( κ + 1 ) σr

(58)

Then, Clapeyron theorem and Eq. (6) provide

Ue =

1 1 q20 q0 V = . 2 2 K

The Winkler constant K therefore is



1   2 2 2 1 q2o K= = 8Gq20 {(κ + 1 )[(σr(i) ) + (σθ(i) ) ] + 2 (κ − 3 ) σr(i) σθ(i) + 8 (τr(θi) ) } dA 2 Ue Ai i=0

(59)

 −1

,

(60)

150

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Fig. 6. The Winkler constant K normalized over 2 G versus the aspect ratio α = a/b, for ν = 0.3.

which is independent of the contact force q0 . The complete expression of the Winkler constant K under plane strain conditions is



K = 576π G(1 − ν )(1 + ν )2 1 − α 4



2 







1 + α 2 /{2 1 − α 2 [59 − 32ν − 227ν 2 + 92ν 3 − 272ν 4





+2α 2 137 + 16ν − 553ν 2 − 16ν 3 + 368ν 4 + 3α 4 177 − 128ν − 473ν 2 + 132ν 3 + 304ν 4







+12α 6 (1 − ν )(3 + 2ν ) 7 − 3ν − 16ν 2 ]



−72 1 + α 2



   (1 + ν )2 (3 − 4ν ) + 4α 2 1 − ν 2 (3 + 2ν ) + 2α 4 (1 − ν )(3 + 2ν )2 lnα},

and for ν = 0.3 it simplifies to



K = Gπ 1 − α 4



2 

(61)



1 + α 2 /[0.0988 + 0.4670α 2 + 0.3329α 4 − 0.4850α 6 − 0.4136α 8



− 0.3214 + 1.7060α 2 + 3.3018α 4 + 1.9172α 6 lnα ] .

(62)

Fig. 6 reports the normalized Winkler constant K/(2 G) versus the aspect ratio α = a/b, for ν = 0.3 and for plane stress, plane strain, and generalized plane strain [44]. The plane stress and plane strain curves in Fig. 6 have been derived from the analytical solution; the generalized plane strain curve is numerical; numerical solutions have also been employed to validate the plane stress and plane strain analytical curves. The three curves are very similar; plane strain is usually adopted [24]. The following interpolation formula for the normalized Winkler constant K/(2 G) in plane strain and for ν = 0.3, is valid for the interval 0.2 ≤ α ≤ 0.8; its maximum error is 1.6 percent, and it occurs for α = 0.8.

K = 7.27 α 3 − 15.24 α 2 + 6.93 α + 1.34 . 2G

(63)

In [19, Formula (12)] an analytical expression of the Winkler coefficient is reported. Recast with the symbols of the present paper, and with reference to a solid shaft, this formula becomes K/(2 G) = π (1 +ν )(1 +α 2 )/(1-ν ). For α = 0.5 and ν = 0.3, K/(2 G) = 7.29, see Section 4.1, whereas Fig. 6 indicates for plane strain K/(2 G) = 1.90. This disagreement may partially derive from the differences between expressions of 0 and 1 of Eq. (44) and the expression of φ of [19, Eq. (5)] (0 and 1 include r5 terms), and partially from the different approaches employed for defining the relative displacement of the cross sections, since in [19] the radial strain has been integrated along the thicknesses of the hub and the shaft, see formulae (11) of [19], whereas in the present paper an energy definition has been employed, see formula (59). 5. Incipient detachment condition To determine the incipient detachment condition, the contact pressure (1) due to the press-fit alone is compared to the contact pressure exclusively due to the bending couple at the hub extremities, Eq. (24), as in [19]. Consistent with the plane

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

151

Fig. 7. Analytical and FE normalized contact pressure profiles versus the normalized axial coordinate z/L, for α = 0.5, a/L = 1/8, and C = π a3 E/4.

modelling of the shaft and hub contact interaction developed in Section 4, the radial stress components (46)1 and (47)1 display a cosinusoidal profile along the shaft-hub circumferential contact region, namely at r = a (it is recalled that the shaft-hub contact is assumed as bilateral, see the comments in the Introduction). (A review of different pressure profiles employed in practical applications is presented in [45].) Consequently, by employing Eq. (24), the maximum contact pressure at the hub extremity, namely at z = l, due to the bending couple C is

p( l ) =

q0 ( l ) Cl K l = S(λ, ε ), πa E I1 π a

(64)

where l is the hub axial half-length (Fig. 1), and where the function S(λ,ε ) is explicitly reported in Eq. (25). Therefore, the incipient detachment condition becomes

 Cl Kl Eδ  S (λ, ε ) = 1 − α2 . E I1 π a 2a

(65)

By adopting the regrouping of the variables employed in [23], the normalized detachment couple C is (it is recalled that L = 2 l, see Fig. 1)

 2 π 1 − α 2  , S(λ, ε )

Ca2 E a = 2δ EI1 K L

(66)

where 2δ is the diametral interference. Fig. 7 reports the analytical shaft-hub contact pressure profile induced by the bending couple alone, normalized over E. This contact pressure is displayed along the normalized axial direction z/L, and for an angular position for which the shaft bending stress is most tensile, i.e. for θ = 0 of Fig. 1. Two analytical pressure curves are displayed, referring to a purely flexural (Euler-Bernoulli) beam model, and to a shear-deformable (Timoshenko) model, respectively. The behaviour of the purely flexural beam is recovered from the shear-deformable beam for χ 0 = χ 1 = 0, see Eqs. (10) and (12). The FE pressure forecasts for a sharp hub bore edge are included for comparison. The following values have been adopted: α = 0.5, a/L = 1/8, and C = Eπ a3 /4, so that the shaft maximum bending stress is unity. Since in the analytical solution a bilateral contact is assumed between the shaft and the hub bore, the possible occurrence of negative contact pressures is justified. The contact pressure curves of Fig. 7 exhibit a bathtub-shaped profile; the analytical shear-deformable model more accurately replicates the FE forecasts than its purely flexural counterpart, and, therefore, it has been preferred in this study, see the Introduction and Section 3. For low values of z/L, i.e., at the shaft-hub central contact region, the contact pressure remains reasonably constant, and the degree of agreement between the analytical and the FE forecasts is good. At the

152

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Fig. 8. The superposition of the contact pressure due to the press-fit alone and to the bending couple alone.

shaft-hub contact extremities, where z/L = 0.5, the contact pressure exhibits a finite peak in the analytical solution, and a (virtually) infinite peak in the FE solution. Consequently, an unavoidable deviation is appreciable between the analytical and the numerical curves especially at the shaft-hub contact extremity, see the remarks at the beginning of Section 3. Fig. 8 graphically illustrates the approximate approach favoured in this study for the determination of the detachment couple. The contact pressure induced by the press-fit alone is assumed to be constant in the axial direction, Fig. 8(a), since the lateral pressure bumps may be neglected, as commented in Section 2. The contact pressure induced by the application of a bending couple C, being based on a bilateral modelling of the shaft-hub contact, is partially tensile and partially compressive (Fig. 8(b)). As discussed at the end of Section 3, the modelling of the shaft-hub contact deformability in terms of the Winkler foundation excludes the presence of a concentrated force in the description of the contact pressure, and it produces a finite pressure peak at the contact extremities. Finally, the superposition of the two pressure profiles generates a wholly compressive contact pressure that, for a particular value of the couple C, becomes null at the contact extremities, thus fulfilling the condition of incipient shaft-hub detachment (Fig. 8(c)). Comments are presented in the following concerning the validity of the assumption of a linear distribution in the axial direction for the contact force qi , Fig. 4, due to the shaft bending alone and not to the press-fit. A linear approximation is favoured in [20]; in [19] it is observed that for the loading of [19, Fig. 1] the analytical distribution of the reaction force becomes linear for a hub of small axial length. If a linear force distribution is employed for the problem of Fig. 1, as a result of symmetry the reaction force due to the shaft bending alone would be constant, whereas the presence of lateral reaction force peaks is expected on physical grounds. In fact, in the problem of Fig. 1, pertaining to the first group discussed in the Introduction, a shaft is laterally bent and centrally surrounded by a hub (Fig. 9(a)). By adopting a beam-like, purely flexural model for the shaft and the hub, and by neglecting the deformations of the shaft and hub cross sections (i.e., by excluding a Winkler foundation), the reaction force between the shaft and the hub is a pure couple applied at both the hub extremities (Fig. 9(b)). (The values C/2 of the reaction couples drawn in Fig. 9(b) refer to a shaft-hub geometry characterized by I0 = I1 .) Instead, if the shear and the Winker compliances are included, as in the present paper, the concentrated couples are dampened to become two adjacent contact force peaks, whose areas have been shaded, and whose resultant forces have been represented in Fig. 9(c), see Fig. 7. It clearly emerges from Fig. 9(c) that for the problem of Fig. 1 the reaction force cannot be acceptably modelled in terms of a linear distribution. Conversely, in the problem of Fig. 9(d), pertaining to the second group, the bending couple C is applied to the shaft protruding extremity and to a hub side. Fig. 9(e) shows the mutual loading for a purely flexural model, where the contact reaction is constituted by concentrated couples C/2. Following [17, Figs. 4–8], Fig. 9(f) reports a reaction force distribution based on FE, but also valid for a shear deformable beam model equipped with a Winkler foundation; since the reaction force distribution is skew-symmetric in character, it may be acceptably modelled in terms of a linear approximation. In conclusion, it is believed that the linear approximation of the contact force is plausible for the press-fits belonging to the second group defined in the Introduction, whereas the linear approximation becomes unrealistic for press-fits belonging to the first category. This observation confirm the different mechanical response of the two groups examined in the Introduction. In [17], p. 54, an approximate shaft-hub pressure distribution is proposed, that is sufficiently adaptable to mimic the pressure peaks appearing at the shaft-hub contact extremities. This curve depends on two constants, that are calibrated versus the FE forecasts. The employment of a pressure curve exhibiting lateral peaks confirms that the linear pressure curve of [20] is often an unsuitable approximation. 6. Comparison with FE predictions In this section, the accuracy is explored of Eq. (66) in determining the bending couple that causes the shaft-hub incipient detachment for technically significant shaft-hub geometries. The errors ascribable to the several approximations adopted, discussed in the Introduction, are cumulatively quantified by comparing the analytical solution based upon formula (66) to selected FE predictions extracted from [23], referring to a sharply edged hub bore. The analytical solution here developed confirms the regrouping of the variables defining the normalized detachment bending couple employed in [23].

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

153

Fig. 9. (a) The shaft loading with a couple C applied to both shaft extremities; (b) the appearance of contact reactions in terms of couples in a purely flexural model; (c) the dampening of the couples as the shear and Winkler compliances are included; (d) the shaft loading with a couple C applied to the shaft protruding extremity; (e) the mutual contact reactions for a purely flexural model; (f) the qualitative distributed contact force.

The commercial FE program MSC Marc 2013 has been employed in this study. The 3D mesh was formed by about 60 0,0 0 0 nodes, whose element size grades smoothly from 2.2 × 10−5 to 0.1 times the shaft radius a. The nodes in contact between the shaft and the hub rounded edges in the axial direction were more than ten [23]. The symmetry of this problem was exploited to contain the number of nodes. For selected test cases an even thinner mesh was employed to assess the numerical convergence. Fig. 10 presents a typical mesh adopted in the vicinity of the shaft-hub contact extremities. The shaft-hub separation condition was recently reconsidered in [46], where the nonlinear progressive contact problem was dealt with by likening it to an auxiliary problem extracted from the fracture mechanics realm, see also [47]. The results in terms of detachment couple retrieved with this completely different approach correlate favourably with the FE predictions of [23]; the maximum error of 1.7 percent for the 25 cases examined testifies to the validity of the detachment couple forecasts of [23], assumed as reference values in this paper. Conversely, a certain disagreement between FE and beam-like analytical forecasts is unavoidable, since he beam model adopts simplifying assumptions (e.g. conservation of plane sections in bending), absent in the 3D FE model. In Fig. 11 the normalized detachment bending couple Ca2 /(IEI1 ), where I is the diametral interference (I = 2δ ) and I1 is the moment of inertia of the shaft, is reported versus the aspect ratio α = a/b for two curves referring to a/L = 1/3 and 1/8, respectively. The abscissa variable covers the interval 0.3 ≤ α ≤ 0.7. Two sets of curves are reported, namely the FE and the analytical shear-deformable beam-like forecasts. The analytical solution considerably overrates the FE predictions. In particular, for the common value α = 0.5, and for a/L = 1/8, the error with respect to the FE predictions is 53 percent, whereas for α = 0.5 and for a/L = 1/3, the error becomes 89 percent. This overprediction is partially imputable to the fact, thoroughly discussed in Section 2, that the contact pressure smoothly decays from its flattish central value to zero, Fig. 2, and, consequently, the lateral pressure bumps disappears, provided that the shaft-hub detachment reaches the flat portion of the hub bore axial profile. Therefore, expression (66) of the detachment couple, neglecting the lateral pressure bumps, describes the occurrence of a detachment of finite size along the flat portion of the hub bore profile, rather than the occurrence of the incipient detachment condition, see Section 2, thus unavoidably overrating the incipient detachment couple.

154

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Fig. 10. The FE mesh in the vicinity of the shaft-hub contact extremities. The presence of an initial gap between the shaft and the hub, due to the hub filleted edge, is modelled with a gap function.

Fig. 11. Normalized detachment couple Ca2 /(IEI1 ) versus a/b for a/L = 1/3 and 1/8, in the situation of incipient shaft-hub detachment, for ν = 0.3. In the FE study the hub edges are sharp.

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

155

Fig. 12. Shaft with protrusion (a); shaft without protrusion (b).

The effect of the radius of the hub bore rounded edges has not been included into this analytical solution. It is frankly admitted that it is complex to incorporate the rounded edge effect into a beam-like theory. 7. Improved solution As it appears from Section 6, the analytical solution achieved in this paper significantly overrates the FE forecasts, the error ranging between 53 percent and 89 percent for the geometries examined. The following comments describe a different physical interpretation of the Winkler model that yields a modified value of its constant, and it enables the above error of the analytical solution to be reduced to more acceptable values. It is recalled that the axial distribution of the shaft-hub contact force due to the bending alone and not to the interference fit remains reasonably constant in the central portion of the contact, and that it exhibits a pronounced peak at the contact extremities (Fig. 7). Since along the majority of the shaft-hub contact extent the contact force remains relatively constant, it seems appropriate to compute the Winkler constant by adopting a plane modelling of this contact, and to employ this Winkler constant along the whole shaft-hub contact length, even in the lateral regions where the contact force peak occurs. In evaluating the Winkler constant according to a plane modelling, the in-plane deformability of both the shaft and the hub is accounted for, as done in Sections 4 and 6. Conversely, the most relevant aspect describing this shaft-hub incipient detachment problem is the correct evaluation of the separation couple, which in turn involves an accurate description of the shaft-hub contact zones especially in the vicinity of the contact extremities, where the above commented contact force peaks take place, i.e., where the title contact problem mostly deviates from a plane modelling. Consequently, to get an accurate prediction of the shaft-hub detachment couple, it seems vital to thoroughly model the title contact problem especially at the contact extremities. The following observations address the approximate evaluation of the Winkler constant at the shaft-hub contact ends. Since the shaft portion protruding from the hub extremities is not directly loaded by the shaft-hub contact force, but it is merely bent by the applied couple C, Fig. 1, its cross section remains virtually undeformed apart from the Poisson’s effect associated to bending [33,34]. Consequently, the in-plane deformation of the shaft section aligned with the hub lateral faces is hindered by the axially contiguous, virtually undeformed sections of the shaft protruding portion. If, conversely, the shaft does not protrude from the hub lateral faces, this restraining effect is missing, see also the comments in Section 2 regarding the absence of the lateral pressure peaks in an axisymmetric press-fit. Fig. 12 shows two shaft-hub couplings. In Fig. 12(a) the shaft protrudes from the hub lateral face, and it is loaded by the couple C, whereas in Fig. 12(b) the shaft extremal face is aligned with the hub lateral face, so that no shaft protrusion occurs. Fig. 13 displays the corresponding FE deformations of the two shafts, in the presence, Fig. 13 left, and in the absence, Fig. 13 right, of shaft protrusion. This FE study refers to a situation of incipient shaft-hub detachment, and to α = 0.5 and to a/L = 1/4; it shows that the contraction of the shaft vertical diameter due to the shaft-hub contact force in the absence of a shaft protrusion is about 1.5 times its analogue in the presence of a shaft protrusion. This result suggests that at the shaft-hub contact extremities the Winkler constant may be evaluated by adopting the extremal assumption of neglecting the in-plane deformation of the shaft cross section, i.e. by assuming the shaft cross section as rigid, and by accounting for the hub in-plane deformability alone. As a result of the above observations, the best possible modelling of the title problem would probably consist in employing an axially variable Winkler constant, e.g. [48]; in the contact central portion, the Winkler constant might be computed according to a plane modelling that accounts for the in-plane deformability of both the shaft and the hub, see Section 4; at the shaft-hub contact extremities, the Winkler constant might be evaluated by considering the hub in-plane deformability alone. For α = 0.5, the increase of the Winkler constant due to the assumption of rigid shaft is about 20 percent, see also Section 4.1.

156

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

Fig. 13. FE vertical displacement of a protruding (left) and unprotruding (right) shaft, in the situation of incipient shaft-hub detachment, and for α = 0.5 and a/L = 1/4.

However, the two following observations suggest that a simpler approach may be preferable in practice, in which a single value of the Winkler constant is adopted for the whole contact length, as in [19]. First, an axially variable Winkler constant could preclude the possibility of achieving a closed-form solution of the shaft-hub contact force distribution. Secondly, the analytical solution of Section 4 shows that an increase of the Winkler constant of 20 percent produces a reasonably comparable raise of the lateral contact force peaks of 11 percent, and an appreciably lower diminution of the central contact force, of 3 percent. So, the desired raise of the lateral contact force peak is not accompanied by an undesired modification of the central distribution of the contact force, where its axial gradient is low. It may be concluded that, since the lateral contact pressure peaks are more sensitive than the central contact pressure to variations of the Winkler constant, and since the correct evaluation of the lateral contact pressure peaks in determining the detachment couple is more relevant than the determination of the central contact force, it is important that the single value adopted for the Winkler constant be realistic at the shaft-hub contact extremities more than at the contact central portion. In summary, the best possible modelling of the title contact problem that may be achieved by employing a single value of the Winkler constant, is to adopt, along the whole shaft-hub contact region and even along the central portions far from the contact extremities, the value of the Winkler constant that is correct at the shaft-hub extremities, i.e., the value computed by neglecting the shaft in-plane deformability. In Fig. 14 the normalized detachment bending couple Ca2 /(IEI1 ) computed for the improved value of the Winkler constant is reported versus the aspect ratio α = a/b for two curves referring to a/L = 1/3 and 1/8, respectively. The abscissa variable covers the interval 0.3 ≤ α ≤ 0.7. Two sets of curves are reported, namely the FE and the analytical shear-deformable forecasts. The theoretical modelling favoured in this section considerably lowers the error of the analytical forecasts with respect to those of Fig. 10; the error now ranges between 34 percent and 59 percent for the geometries examined. It is expected that the error increases as the hub axial length decreases. Analytical detachment couple forecasts in better agreement with the FE predictions would be achieved by employing the Winkler constant provided by a uniform distribution of the shear stresses across the hub and shaft sections, according to a gravitational field, Section 4.1. Considering the deformability of both the shaft and the hub cross sections, for α = 0.5 and a/L = 1/8, then K/(2 G) is 3.68, and Ca2 /(IEI1 ) = 0.41; considering only the deformability of the hub cross section, for α = 0.5 and a/L = 1/8, then K/(2 G) is 4.81, and Ca2 /(IEI1 ) = 0.36. The errors of the analytical forecast versus their FE counterparts become less than 10 percent and less than 5 percent, respectively. It is however felt that this model, based upon an unrealistically uniform distribution of a shear stress, is physically less justified than its analogue based on the exact distribution of the stresses due to a shear force, developed in Section 4.2.

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

157

Fig. 14. Normalized detachment couple Ca2 /(IEI1 ) versus a/b for a/L = 1/3 and 1/8, in the situation of incipient shaft-hub detachment, for ν = 0.3. In the FE study the hub edges are sharp. The improved Winkler foundation model has been employed.

In conclusion, this simplified beam-like model may only predict an indicative value of the incipient shaft-hub detachment couple C, rather than its accurate threshold. It is however believed that this analytical solution still retains its technical usefulness. It is recalled that, for the specific loading depicted in Fig. 1, numerical values of the incipient shaft-hub detachment couple appear to be available only in [23]; in fact, in the rest of the pertinent literature only numerical values are available that refer to the loadings pertaining to the second group defined in the Introduction. It seems therefore impossible to assess the analytical forecasts of this paper versus additional reference values. However, to appraise the order of magnitude of the detachment value of the applied couple, it was decided to carry out a comparison with [20, Formula (7)], dealing with a shaft protruding from only one side of the hub, see [20, Fig. 1(c)]. In fact, as noted in the Introduction, if the hub axial length is sufficiently high, it may be surmised on physical grounds that the detachment loading of [20, Fig. 1(c)] becomes comparable to that of the title problem. Recast with the symbols of the present paper, the normalized incipient detachment couple Ca2 /(IEI1 ) according to [20, Formula (7)] in the absence of friction becomes

Ca2 = IEI1

  π 1 − α 2 a2 l 2 48I1

 2

=

 l 1  1 − α2 12 a

,

(67)

where formula (1) has been employed for the evaluation of the press-fit pressure; in addition, the axial length of the shaft portion in contact with the hub of [20, Fig. 1(c)] has been assumed as half the hub total axial length L of Fig. 1, since the geometry of Fig. 1 may be interpreted as a mirroring of the geometry of [20, Fig. 1(c)] with respect to a vertical axis passing through the shaft extremity not protruding from the hub. For α = 1/2 and l/a = 4, the normalized incipient detachment couple is 1, whereas the FE forecast is about 0.4. Formula (67) therefore provides a value 2.5 times higher than the FE forecast. One source of inaccuracy is the adoption in [20] of a linear axial contact force distribution between the shaft and the hub, which, especially for long hubs, fails to model the lateral contact force peak, Fig. 7, and, therefore, it overestimates the detachment couple. It is finally briefly mentioned that the influence on the detachment couple of the cross section deformation caused by the axial bending stresses and due to the Poisson’s ratio effect has been investigated with the aid of the beam modelling.

158

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

It was however found that, for the geometries examined, this contribution to the evaluation of the detachment couple is negligible, being lower than 1 percent. 8. Possible occurrence of yielding The following observations aim at assessing whether the incipient detachment situation produces high enough stresses to cause yielding, an event, this, that would jeopardize the reliability of an elastic approach. It is preliminarily observed that, as a result of the press-fit alone, the bore surface in the central part of the hub is subjected to a markedly deviatoric stress state (the radial and hoop stresses exhibit opposite signs), whereas the shaft stress state is more hydrostatic (both the radial and hoop stresses are compressive). Moving to the shaft-hub contact extremities, the hoop stress in the hub becomes compressive, especially for small fillet radii of the hub bore [12]. The photoelastic readings of [41] experimentally confirm the above sign of the hoop stress in the vicinity of the contact extremities for a comparable fit problem. It is known that a hydrostatic stress field restricts the tendency to yielding [49, p. 21]. An assessment on the applicability of the elastic solution may therefore be restricted to the central portion of the hub bore surface. An experimental support to the above conclusion is provided by the axial profile of the web bore of a built-up crankshaft, after journal dismantling. In [9, Fig. 34] a greater permanent set was recorded at half length of the shrink fit rather than at its extremities, thus corroborating the occurrence of a more pronounced yielding at the press-fit centre. (However, it must be frankly admitted that the measurements presented in Fig. 16 of the same paper provided contrasting results.) Since the contact force remains reasonably constant in the shaft-hub central zone, the hub stress field may be evaluated for simplicity at the shaft-hub contact centre. To compute the hub stresses at the centre of its bore surface and due to the shaft-hub bending alone, the relative displacement V is determined from Eq. (11) by employing formula (21) of  and the expressions of the constants c1 and c4 from Eqs. (22) and (23). The value of V is then computed for the hub bore centre by imposing z = 0, and the distributed load q0 is derived from V with the aid of the Winkler foundation constant, see Eq. (6). Then, the in-plane stresses σ r and σ θ at the hub bore centre, exclusively due to the shaft-hub bending, are evaluated from Eq. (46) by employing the values provided in (56) for the constants A0 , B0 , C0 , D0 . (The stress component τ rθ is null at the shaft-hub interface, since the contact has been assumed as frictionless.) The out of plane, bending stress σ a may be computed as M0 a/I0 , where M0 is evaluated at the press-fit centre by combining the second Eq. (8) with the equilibrium equation M0 + M1 = C, to provide



M0 =

C − ddz  1

z=0 + II10

E I1

.

(68)

Moving to the hoop stresses at the press-fit centre, at the hub inner border, and due to the press-fit alone, the radial stress may be computed from formula (1) as σ r = –ppf , whereas the Lamé solution provides for the hub hoop stress the value of σ θ . Superposition of the two above stress fields, referring to the shaft-hub bending and to the press-fit, respectively, provides preliminary indications about the possible occurrence of yielding. A numerical example is carried out in the following for a shaft-hub geometry defined by α = 1/2 and a/L = 1/8. It is assumed that the shaft bends exhibiting concavity upwards. The hub stresses due to the bending alone are a radial stress (due to the shaft-hub contact as a result of the shaft bending) and an axial bending stress. The radial stress σ r is evaluated at the hub central portion and at the points B and E of Fig. 5, and its maximum value in the central part of the shafthub contact may be estimated from Fig. 7, where p/E = –σ r /E in the hub central zone ≈ –0.05, and considering that the pressure curves are proportional to the variable 4C/(π a3 E). Therefore σ r = ±0.05 × E × 4C/(π a3 E) = ±0.063 C/a3 , where the sign + (–) is adopted for σ r at the point B (E) of Fig. 5. The hub axial stress σ a due to bending is computed by evaluating the fraction of the total couple C that bends the hub. This fraction of C is denoted M0 , and its value is found from formula (68) to be 0.89 C. (An approximate approach for computing M0 is to impose equal curvature to the shaft and to the hub. From this approximate imposition, M0 = CI0 /(I0 +I1 ). For the reference geometry, M0 = 0.94C.) Then, σ a = ±0.89 Ca/I0 = ±0.89 Ca/[π (b4 -a4 )/4] = ±0.08 C/a3 , where the sign + (–) holds for the point E (B) of Fig. 5. Reference is made to the value of C of incipient detachment for the reference configuration. From Fig. 14, Ca2 /(IEI1 ) ≈ 0.4, from which C ≈ 0.4 IEI1 /a2 = 0.32 IEa2 . By substituting this value for C into the expressions of σ r and σ a due to bending alone, one obtains σ r = ±0.005 IE/a, σ a = ±0.03 IE/a. Moving to the hub stresses due to the press-fit alone, they are a radial stress and a circumferential stress, and they attain equal values at the points B and E, respectively. The radial stress σ r is computed from Eq. (1), σ r = –ppf = –EI(1–α 2 )/(4a) = – 0.20 EI/a; the circumferential stress σ c is determined from the Lamé solution: σ c = ppf × (1 +α 2 )/(1–α 2 ), σ c = EI(1 +α 2 )/(4a) = 0.31EI/a. At the point B, the static Tresca equivalent stress is [0.31–(–0.2 + 0.005)]IE/a = 0.51 IE/a. At the point E, the equivalent stress is [0.31–(–0.2–0.005)]IE/a = 0.52 IE/a. Typical values of IE/a adopted for steel press-fits are of the order of 40 × 10−3 × 210,0 0 0/20 = 420 MPa. For this value, the equivalent stress at the point E is less than 218 MPa, a relatively low stress value. It may be concluded that shaft-hub detachment may occur without implying the occurrence of distributed yielding, and, therefore, an elastic modelling is suitable for predicting the shaft-hub detachment.

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

159

9. Conclusions A mathematical modelling of a shaft-hub press-fit subjected to a bending couple applied to the shaft extremities has been developed, and the value of the bending couple beginning the shaft-hub detachment, which produces undesired wear and fretting fatigue, has been determined analytically. The beam-like model presented in [19] has been improved by accounting for the hub flexural and shear deformation, the shaft shear deformation, and by developing a more realistic Winkler-type foundation model. As in [19], friction has been neglected. The two contact pressures due to the press-fit alone and to the application of the bending couple alone have been modelled separately. The contact pressure imputable to the press-fit has been forecast in terms of the plane Lamé solution. The analytical modelling of the shaft-hub contact has been developed in terms of two shear-deformable Timoshenko beams, mutually connected by a Winkler foundation; the Winkler constant has been theoretically determined. Based on this modelling, the bending couple producing the incipient detachment between the shaft and the hub has been analytically determined in term of the shaft-hub geometry, of the initial interference, and of the elastic constants. Two models of the Winkler foundation have been considered. In the first modelling, a plane contact between the shaft and the hub has been assumed, that expresses the combined deformability of both the shaft and the hub. In the second modelling, the stiffening effect exerted by the shaft portion protruding from the hub on the adjacent shaft part that is in contact with the hub has been accounted for, and, therefore, only the in-plane hub deformability has been considered. Comparisons with selected FE forecasts have shown that the first modelling produces an incipient detachment couple that appreciably overrates the FE forecasts, whereas the second modelling lowers the error down to technically acceptable predictions. This beam-like modelling is unable to account for the influence of the hub bore edge radius. Possible improvements of the analytical model consist in introducing the frictional effects, in employing a more advanced contact model, e.g. a two parameter elastic interface, [30], and in modelling the effect of the hub bore edge radius by likening it to an auxiliary problem extracted from the fracture mechanics realm [46,47]. Acknowledgements The authors are grateful to an anonymous referee for drawing their attention to Ref. [17] and related bibliography. References [1] GWJ Technology GmbH, Benutzerhandbuch zur webbasierten Berechnungs so_ware eAssistant, http://www.eassistant.eu/_leadmin/dokumente/ eassistant/pdf/Hilfe/Handbuch/eAssistantHandb.pdf, 52, (2013). [2] http://www.mitcalc.cz/doc/shaftconf/help/en/shaftconftxt.htm [3] E. Leidich, J. Vidner, B. Bruzek, Integral approach for endurance limit evaluation on shrink-fitted assemblies, Bull. Appl. Mech. 18 (2009) 44–49. [4] E. Leidich, B. Bruzek, J. Vidner, Investigation on the shrink-fitted assemblies with a circumferential groove, in: Proceedings of the ICMFF9, 2013, pp. 761–768. [5] N. Antoni, F. Gaisne, Analytical modelling for static stress analysis of pin-loaded lugs with bush fitting, Appl. Math. Model. 35 (2011) 1–21. [6] N. Antoni, A study of contact non-linearities in pin-loaded lugs: separation, clearance and frictional slipping effects, Int. J. Non-Linear Mech. 58 (2014) 258–282. [7] L. Marmorini, A. Baldini, E. Bertocchi, M. Giacopini, R. Rosi, A. Strozzi, On the loosening mechanism of a bush press-fitted in the small end of a connecting rod, Proc. Int. Mech. Eng. D: J. Autom. Eng. 226 (2012) 312–324. [8] N. Antoni, Contact separation and failure analysis of a rotating thermo-elastoplastic shrink-fit assembly, Appl. Math. Model. 37 (2013) 2352–2363. [9] A.S.T. Thomson, A.W. Scott, C.M. Moir, Shrink-fit Investigations on simple rings and on full-scale crankshaft webs, in: Proceedings of the Institution of Mechanical Engineers, 168, 1954, pp. 797–830. [10] A. Strozzi, P. Vaccari, On the press fit of a crankpin into a circular web in pressed-up crankshafts, J. Strain Anal. Eng. Des. 38 (2003) 189–199. [11] E. Radi, A. Strozzi, Jeffery solution for an elastic disk containing a sliding eccentric circular inclusion assembled by interference fit, Int. J. Solids Struct. 46 (2009) 4515–4526. [12] A. Strozzi, A. Baldini, M. Giacopini, E. Bertocchi, L. Bertocchi, Normalization of the stress concentrations at the rounded edges of a shaft-hub interference fit, J. Strain Anal. Eng. Des. 46 (2011) 478–491. [13] D. Croccolo, M. De Agostinis, N. Vincenzi, Normalization of the stress concentrations at the rounded edges of a shaft–hub interference fit: extension to the case of a hollow shaft, J. Strain Anal. Eng. Des. 47 (2012) 131–139. [14] D. Croccolo, M.N. Vincenzi, Stress concentration factors in compression—fit couplings, Proc. Ints. Mech. Eng. C: J. Mech. Eng. Sci. 224 (2010) 1143–1152. [15] M. Ciavarella, A. Baldini, J.R. Barber, A. Strozzi, Reduced dependence on loading parameters in almost conforming contacts, Int. J. Mech. Sci. 48 (2006) 917–925. [16] E. Radi, D. Bigoni, A. Tralli, On uniqueness for frictional contact rate problems, J. Mech. Phys. Solids 47 (1999) 275–296. [17] T. Smetana, Untersuchungen zum Übertragungsverhalten biegebelasteter Kegel-und Zylinderpressverbindungen, Shaker-Verlag, 2001. [18] M. Ast, H. Rösle, R. Schenk, FEM-Analysereibschlüssiger Welle-Nabe-Verbindungen, in: Proceedings of VDI-Berichte 1384, ‘Welle-Nabe-Verbingunen System-komponenten im Wandel, Fulda, Germany, 1998, pp. 227–244. [19] G. Lundberg, Spannungen in Pressverbänden bei Belastung, 3, Die Kugellager Zeitschrift SKF, 1958, pp. 55–63. [20] A.S.T. Thomson, A.W. Scott, W. Ferguson, Strength of shrink fits in bending and combined bending and torsion, Engineer 213 (1962) 178–180. [21] D.J. White, J. Humpherson, Finite-element analysis of stresses in shafts due to interference-fit hubs, J. Strain Anal. 2 (1969) 105–114. [22] M.D. Garnett, T.R. Grimm, Finite Element analysis of interference-fitted shafts subjected to bending, Comput. Engng ASME 2 (1989) 273–280. [23] A. Strozzi, E. Bertocchi, A. Baldini, S. Mantovani, Normalisation of the stress concentrations at the rounded edges of an interference fit between a solid shaft subjected to bending and a hub, Mech. Based Des. Struc. Mach. 44 (2015) 405–425. [24] J. Castillo, J.R. Barber, Lateral contact of slender prismatic bodies, Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 453 (1997) 2397–2412. [25] M. Ciavarella, D.A. Hills, G. Monno, The influence of rounded edges on indentation by a flat punch, Proc. Inst. Mech. Eng. C: J. Mech. Eng. Sci. 212 (1998) 212–319. [26] W. Peppler: Preßverbindungen. VDI-Berichte, p. 60–72. Düsseldorf: VDI-Verlag 1956. [27] T. Schwämmle, Betriebsverhalten von konventionellen und fugendruckhomogenisierten Pressverbänden unter Biegelast, (Dissertation) Stuttgart University (2010).

160

E. Radi et al. / Applied Mathematical Modelling 50 (2017) 135–160

[28] J.F. Heydt, Untersuchungen zum dynamischen Verhalten von topologisch optimierten Pressverbänden bei Umlaufbiegung, (Dissertation) Stuttgart University (2012). [29] H. Fessler, D.C. Fricker, A study of stresses in alumina universal heads of femoral prostheses, Proc. Inst. Mech. Eng. Part H J. Eng. Med. 203 (1989) 15–34. [30] A.D. Kerr, Elastic and viscoelastic foundation models, J. Appl. Mech. 31 (1964) 491–498. [31] S. Timoshenko - Strength of Materials . Vol 2 , Van Nostrand , Canada, 1956. [32] L. Lanzoni, E. Radi, A loaded Timoshenko beam bonded to an elastic half plane, Int. J. Solids Struct. 92 (2016) 76–90. [33] A. Strozzi, Theoretical analysis of a segmented locking ring for shell-bottom connection in pressure vessels, J. Strain Anal. Eng. Des. 19 (1984) 153–160. [34] F. Essenburg, Shear deformation in beams on elastic foundations, J. Appl. Mech. 29 (1962) 313–317. [35] E. Dragoni, A. Strozzi, A. Analysis of a split ring inserted into a circular housing, J. Strain Anal. Eng. Des. 21 (1986) 59–70. [36] A. Strozzi, A. Baldini, M. Giacopini, R. Rosi, E. Bertocchi, Contact stresses within a split ring inserted into a circular housing, J. Strain Anal. Eng. Des. 44 (2009) 671–688. [37] J.D. Renton, A note on the form of the shear coefficient, Int. J. Solids Struct. 34 (1997) 1681–1685. [38] A.E.H. Love, A Treatise On the Mathematical Theory of Elasticity, Dover Publications, New York, 1944 (Vol. 1). [39] J.R. Barber, Elasticity, third ed, Kluver, Dordrecht, 2010. [40] K.T. Chau, X.X. Wei, Stress concentration reduction at a reinforced hole loaded by a bonded circular inclusion, J. Appl. Mech. 68 (2001) 405–411. [41] H. Fessler, H.B. Padgham, A contribution to the stress analysis of piston pins, J. Strain Anal. Eng. Des. 1 (1966) 422–428. [42] M.H. Sadd, Elasticity: Theory, Applications, and Numerics, Academic Press, 2009. [43] K.C. Ho, K.T. Chau, An infinite plane loaded by a rivet of a different material, Int. J. Solids Struct. 34 (1997) 2477–2496. [44] A.D. Cheng, J.J. Rencis, Y. Abousleiman, Generalized plane strain elasticity problems, WIT Trans. Model. Simul. (1970) 10. [45] A. Strozzi, F. De Bona, Hoop stresses in the con-rod small end, Proc. Ints. Mech. Eng. D: J. Automob. Eng 219 (2005) 1331–1345. [46] E. Bertocchi, S. Mantovani, A. Strozzi, Contact stresses in shaft-hub press-fits with rounded bore edges in the presence of shaft bending, in: Proceedings of the 45th Meeting AIAS, Triest, Italy, 2016 7–10 September(In Italian). [47] A. Strozzi, E. Bertocchi, S. Mantovani, M. Giacopini, A. Baldini, Analytical evaluation of the peak contact pressure in a rectangular elastomeric seal with rounded edges, J. Strain Anal. Eng. Des. 51 (2016) 304–317. [48] P. Djondjorov, V. Vassilev, V. Dzhupanov, Dynamic stability of fluid conveying cantilevered pipes on elastic foundations, J. Sound Vibr. 247 (2001) 537–546. [49] K.L. Johnson, Contact Mechanics, Cambridge University Press, 1987.