# EHL Analysis of Spiral Bevel Gear Pairs considering the Contact Point Migration due to Deformation under Load

Authors: Xiaoyu Sun; Yanping Liu; Yongqiang Zhao; Ming Liu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The actual contact point of a spiral bevel gear pair deviates from the theoretical contact point due to the gear deformation caused by the load. However, changes in meshing characteristics due to the migration of contact points are often ignored in previous studies on the elastohydrodynamic lubrication (EHL) analysis of spiral bevel gears. The purpose of this article is to analyze the impact of contact point migration on the results of EHL analysis. Loaded tooth contact analysis (LTCA) based on the finite element method is applied to determine the loaded contact point of the meshing tooth pair. Then, the osculating paraboloids at this point are extracted from the gear tooth surface geometry. The geometric and kinematic parameters for EHL simulation are determined according to the differential geometry theory. Numerical solutions to the Newtonian isothermal EHL of a spiral bevel gear pair at the migrated and theoretical contact points are compared to quantify the error involved in neglecting the contact point adjustment. The results show that under heavy-loaded conditions, the actual contact point of the deformed gear pair at a given pinion (gear) roll angle is different from the theoretical contact point considerably, and so do the meshing parameters. EHL analysis of spiral bevel gears under significant load using theoretical meshing parameters will result in obvious errors, especially in the prediction of film thickness.


---

## Body

## 1. Introduction

Spiral bevel gears, the most common form of gearing used for transmitting significant power between nonparallel shafts, have the potential to wear, pitting, and scuffing due to lubrication failure. Therefore, the lubrication behaviour of these gears has long been the concern of many research studies. However, fruitful results have only been yielded after the 21st century thanks to the significant improvement in tooth contact analysis technology, elastohydrodynamic lubrication (EHL) theory, numerical algorithms, and computer performance.Early studies by Simon [1] and Chao [2] analyzed the lubricating properties of spiral bevel and hypoid gears assuming smooth surface due to the insufficient computational resource and techniques at that time. Recently, Xu and Kahraman [3] predicted the transmission efficiency of hypoid gear pairs, where the friction coefficient is calculated by applying a line contact EHL model which is only reasonable for the theoretical conjugate tooth surfaces. In applications, however, it is imperative to modify the tooth profiles of the mating gears to make them mismatched, thereby being capable of tolerating a certain amount of misalignments caused by the deflection of the supporting shaft and housing, as well as installation and machining errors. Due to the mismatch, an elliptical contact footprint, instead of the theoretical line contact one, is formed between each pair of teeth in contact. Kolivand et al. [4] regressed the friction coefficient formulae for the entire lubrication range, while the line contact assumption was still adopted. Simon [5] reported a thermo-EHL analysis of spiral bevel and hypoid gears where although a more realistic elliptical point contact model is applied, the entrainment flow vector is assumed to be along the minor axis of the elliptical footprint. Similarly, Zhang et al. [6] also predicted the tooth surface friction coefficient and wear through thermo-EHL analysis with the assumption that the rolling velocity direction coincides with the minor axis of the contact ellipse. Mohammadpour et al. [7] provided an isothermal Newtonian EHL solution of hypoid gear pairs at high loads, considering angled entrainment flow of the lubricant. This work was then extended to take non-Newtonian shear of thin films into account [8]. Lately, the unified lubrication model covering all the lubrication regions [9] was recognized by more and more researchers in the field. Pu et al. [10] used it to conduct the mixed EHL analysis of spiral bevel and hypoid gears with machined roughness. Cao et al. [11] also employed this model to study the effect of contact path on the lubrication performance, friction, and fatigue life. More recently, Gan et al. [12] further applied this model to investigate the temperature behaviour under mixed lubrication condition.The data required for EHL analysis include the load acting on the meshing tooth pair and the geometric and kinematic meshing parameters at the contact points. In studies mentioned above, the load sharing among simultaneously engaged tooth pairs was evaluated using the loaded tooth contact analysis (LTCA), and the meshing parameters were calculated at the theoretical contact point (TCP) obtained through tooth contact analysis (TCA) which is, however, only valid under the no-load and light-loaded conditions. Many researchers such as Wang [13], Zhou et al. [14], and Sun et al. [15] have revealed that the contact path (which consists of all contact points over a tooth engagement cycle) of the mating spiral bevel gear pair under loaded condition is different from the theoretical one due to deformation of gear teeth and the body caused by the load. It turns out that the contact path changed by adjusting the mounting position of mating gears has a nonnegligible impact on lubrication performance [11]. However, changes in the contact path due to the load-induced deformation are usually ignored in EHL analysis.The present study proposes a method for calculating the meshing parameters at the loaded contact point (LCP) which is synthesized from the load distribution results obtained by LTCA. Once the LCP is obtained, the osculating paraboloids at that point are extracted from the geometry of pinion and gear tooth surfaces. Then, the geometric and kinematic parameters for EHL simulation are determined according to the differential geometry theory. The flowchart of the proposed methodology is shown in Figure1. Numerical solutions to the Newtonian isothermal EHL of a Gleason-type spiral bevel gear pair at the loaded and theoretical (unloaded) contact points are obtained and compared to quantify the error involved in neglecting the contact point adjustment due to deformation caused by the load.Figure 1
Flowchart of the methodology used to calculate the meshing parameters at loaded contact points.


## 2. Determination of Meshing Parameters

TCA algorithms reported in the literatures [16, 17] and many others can be applied to determine the TCPs of spiral bevel gear pairs under no-load and light-loaded conditions where the deformation caused by the load can be ignored. At high loads, the obtainment of LCPs needs LTCA. The LCP, in the present work, is defined as the contact point of the undeformed mating tooth surfaces of the pinion and the gear as shown in Figure 2. It is known that, for a pair of mismatched spiral bevel gears, an elliptical contact footprint is formed between the contacting tooth surfaces due to elastic contact deformation. The geometric center of the contact ellipse can be viewed as the projection point of the LCP on the deformed tooth flanks (Figure 2). Furthermore, the center of the contact pressure distribution, which is readily available from LTCA, is close to the center of the contact ellipse. Therefore, the center of the contact pressure distribution is utilized to determine the LCP in the present work. The center of the contact pressure is sometimes defined as the acting point of the effective meshing force due to contact pressure and named the “effective meshing point” (EMP) [13, 15], which will also be used in the following description.Figure 2
Definition of the loaded contact point.The EMPs at each meshing position over a tooth engagement cycle (from a pair of teeth coming into contact to exit contact) can be obtained through quasistatic LTCA based on either a conventional finite element method [14] or a finite element/contact mechanism method [18]. For the second method, readers may refer to references [13, 19], where a specialized LTCA tool [20] is utilized. In this work, the first method is adopted, and LTCA is performed with the aid of the commercial FEA software ABAQUS [21].The coordinate systems shown in Figure3 are established to describe the relative kinematics of a spiral bevel gear pair. It should be noted that the proposed method is also applicable to the hypoid gear pair having a shaft offset. The global fixed reference frame Oxyz whose origin is fixed to the intersection point of the pinion and gear rotation axes is established automatically in the simulation environment of ABAQUS, and the results from LTCA are given in this coordinate system. O1x1y1z1 and O2x2y2z2 are the local fixed reference frames used to describe the rotational positions of the gear and the pinion, respectively. The gear and the pinion rotate about the x-axis of their local fixed coordinates. γ is the angle between the gear and the pinion axis. Ogxgygzg and Opxpypzp are the body-fixed coordinate systems of the gear and the pinion in which their tooth surfaces are represented. The symbol Acb will be often used in this article where the subscript c indicates the coordinate system in which A is represented and the superscript b indicates the object to which A belongs.Figure 3
Coordinate systems of a spiral bevel gear pair.The procedures of using ABAQUS for LTCA have been introduced in detail by Zhou [14] and Sun et al. [15]. The position vector of EMP can be obtained conveniently as long as the history output XN (the center of the total force due to contact pressure [21]) of the contacting tooth pair is requested. Besides, the effective meshing force vectors that will be used to determine the position of tooth surfaces after deformation can be acquired from the history output CFN (the total force due to contact pressure [21]).The next focus of the present study is the calculation of the coordinates of LCPs corresponding to the EMPs. As shown in Figure4(a), it is assumed that Pm is the EMP at a given meshing position and its position vector represented in the fixed coordinate system Oxyz is rmxm,ym,zm. FmFmx,Fmy,Fmz is the vector of effective mesh force acting on the EMP. Pm can first be projected onto the projection plane P of the gear and pinion teeth as shown in Figure 4(c). As illustrated in Figure 4(b), the coordinates of the projected points for the pinion (PPp) and the gear (PPg) can be calculated by(1)xPp=xm,yPp=ym2+zm2,xPg=xm2+zm2,yPg=−ym.Figure 4
The EMP and the corresponding projection point of the gear and the pinion.
(a)(b)(c)The mathematical model of the tooth flank geometry depends on how the gear teeth are machined but can be generally represented as [16](2)rbs,θ,ϕ=Mbtϕrts,θ,fs,θ,ϕ=nt⋅vtb=0,where rts,θ indicates the cutting surface formed by the rotation of the cutting edge about the axis of the cutter head, s is the parameter describing the cutting edge geometry, θ is the cutter phase angle, and Mbtϕ is the transformation matrix representing the coordinate transformation from the coordinate system fixed to the cutter head to the coordinate system fixed to the gear blank. The form of the matrix Mbtϕ hinges on the motion of the machine tool. Readers can refer to relevant literatures (e.g., [16, 17]) for the definition of the matrix Mbtϕ. The cutting surface will generate a family of surfaces rbs,θ,ϕ in the coordinate system rigidly connected to the gear blank. ϕ is the parameter of the family of surfaces and can either be the blank phase angle or the cradle angle. According to the enveloping theorem, the tooth surface of the machined gear is the envelope surface of the family of surfaces rbs,θ,ϕ. The points on rbs,θ,ϕ which simultaneously satisfy the equation of meshing fs,θ,ϕ are those on the envelope surface (namely, the tooth surface). In the equation of meshing, nt denotes the unit normal vector of the cutting surface and vtb is the velocity in relative motion between the cutting surface and the tooth surface.Applying equation (2) and the coordinate relationship between the tooth surface point and the projection point (as depicted in Figures 3 and 4(c)), the LCPs of the pinion and the gear can be calculated, respectively, by(3)fpsp,θp,ϕp=nt,p⋅vtb,p=0,xPp=xppsp,θp,ϕp,yPp=yppsp,θp,ϕp2+zppsp,θp,ϕp2,fgsg,θg,ϕg=nt,g⋅vtb,g=0,yPg=xggsg,θg,ϕg,xPg=yggsg,θg,ϕg2+zggsg,θg,ϕg2,where xii, yii, and ziii=p,g are the three components of the position vector of the LCP riixii,yii,zii represented in the body-fixed frames of the pinion and the gear (Oixiyizi).The osculating paraboloids at LCP are extracted from the tooth surfaces of the pinion and the gear by the following method. In order to facilitate the following description, the equation of the gear tooth surface can be represented in a two-parameter form as follows:(4)riiθi,ϕi=riisiθi,ϕi,θi,ϕi,i=p,gwhere siθi,ϕi is derived from the equation of meshing fs,θ,ϕ.As mentioned previously, the LCP is different from the theoretical unloaded contact point calculated based on the principle of local conjugation [16]. That is to say, the contact conditions are not satisfied at the LCP for the undeformed gear set. As illustrated in Figure 5(a), the LCPs on the pinion (Pp) and gear (Pg) tooth surfaces do not coincide and their normals (np and ng) are not collinear. If the gear and pinion tooth surfaces are in contact at the LCP, it is most likely that the relative position of the gear and pinion tooth flanks changes due to the deformation of the gear set caused by the load. The total deformation of the gear set is mainly composed of the tooth surface contact deformation, the bending and shear deflection of the gear teeth, and the bending and torsional deformation of the gear body and the shaft, as well as the deformation of housing. The rotation of the gear and pinion teeth is one of the manifestations of the total deformation [22]. They are, respectively, simulated by the rotation angles of the gear (Δθg) and the pinion (Δθp) around their rotational shafts. According to the calculation process of the LCP described before, it is obvious that Pp and Pg will coincide after the rotation as illustrated in Figure 5(b). The gear and pinion tooth surfaces are now represented in the coordinate system Oxyz as follows:(5)rg=MOO1γMO1OgΔθgrgg,(6)rp=MOO2γMO2OpΔθprpp,and at the LCP, there is(7)rg=rp.Figure 5
Coordinate systems applied for the movement of the gear and pinion tooth surfaces.
(a)(b)(c)Usually,γ=π/2, and thus the transformation matrices in equations (5) and (6) are defined as(8)MOOgΔθg=0sinΔθgcosΔθg−1000−cosΔθgsinΔθg,MOOpΔθp=1000−sinΔθpcosΔθp0−cosΔθp−sinΔθp.Then, the rotation angle of the gear tooth (Δθg) and the pinion tooth (Δθp) are determined by simultaneously solving equations (5)–(7) at the LCP.At this time, the tooth surfaces of the gear and the pinion have coinciding position vectors at the LCP, while the contact condition is still not satisfied because the surface normal vectors at the LCP (ng and np) are not collinear as illustrated in Figure 5(b). It is difficult to describe the changing process of the tooth surface position due to the composite deformation of the gear set. Nevertheless, it is reasonable to assume that the deformation makes the normal vectors of the gear and pinion tooth surfaces at the LCP coincide with the direction of the effective mesh force obtained by LTCA. Based on this assumption, the next step is to eliminate the angle between the tooth surface normal vector and the direction vector of the effective mesh force. It is assumed that the to-be-eliminated angles for the gear and the pinion are Δγg and Δγp, respectively (Figure 5(c)). Taking the pinion as an example, the procedure of eliminating the angle Δγp is described. A local coordinate system OLpxLpyLpzLp which is rigidly attached to the pinion tooth surface at the point Pp is established (Figure 5(c)). The unit vectors along the axes of OLpxLpyLpzLp are defined as(9)xLp=npnp,yLp=np×Lmnp×Lm,zLp=xLp×yLpxLp×yLp,where the three components of the direction vector Lmnmx,nmy,nmz of the effective mesh force FmFmx,Fmy,Fmz are obtained by(10)nml=FmlFnet,l=x,y,zFnet=Fmx2+Fmy2+Fmz2.Before rotating, the pinion tooth surface should be represented in the local body-fixed coordinate systemOLxLpyLpzLp:(11)rLp=MOLOprp,and the coordinate transformation matrix from Oxyz to OLxLpyLpzLp is(12)MOLOp=xLp⋅xxLp⋅yxLp⋅z−xpxLp⋅x−ypxLp⋅y−zpxLp⋅zyLp⋅xyLp⋅yyLp⋅z−xpyLp⋅x−ypyLp⋅y−zpyLp⋅zzLp⋅xzLp⋅yzLp⋅z−xpzLp⋅x−ypzLp⋅y−zpzLp⋅z0001,where x, y, and z are the unit vectors along the axes of the coordinate system Oxyz and xp, yp, and zp are the three components of the position vector of Pp represented in Oxyz. Then, the pinion tooth surface along with the coordinate system OLxLpyLpzLp rotates about the axis yLp to the position shown by OL′xL′pyL′pzL′p. The rotation matrix MOL′OLp is(13)MOL′OLp=cosΔγp0−sinΔγp00100sinΔγp0cosΔγp00001.The normal vector of the pinion tooth surface at the LCP is now parallel with the direction vector of the effective mesh force. A similar approach can be used to rotate the gear tooth flank to make its normal vector collinear with the direction vector of the effective mesh force as well. At this point, the tooth surfaces of the gear and the pinion meet the contact conditions at the LCP, and locally, they are osculating paraboloids at the LCP.The tooth surfaces of the gear and the pinion should be represented in the same coordinate system, e.g.,Oxyz, and then the geometric characteristic parameters at the LCP can be determined according to the differential geometry principle [23]. Note that the representation of the tooth surface remains unchanged in the body-fixed coordinate system, and thereby rL′i=rLii=p,g. Therefore, the tooth surfaces of the gear and the pinion in Oxyz are defined by(14)ri=MOOL′irL′i=MOLOi−1MOL′OLi−1rL′i,i=p,gAs depicted in Figure6(a), it is assumed that the position vectors of the LCP for the pinion and gear tooth surfaces represented in Oxyz are rpθp,ϕp and rgθg,ϕg, and the corresponding normal vectors of tooth surfaces at the LCP are npθp,ϕp and ngθg,ϕg. The induced normal curvature of two conjugate surfaces, which is defined as the difference between the normal curvatures of the two surfaces in arbitrary direction, describes the closeness of the two surfaces along the specified direction near the contact point. The smaller the value of the induced normal curvature, the closer the two surfaces. Therefore, the direction in which the induced normal curvatures are minimum and maximum can be, respectively, defined as the direction of the major and minor axes of the instantaneous contact ellipse.Figure 6
(a) Contact relation of the gear and the pinion and (b) principal directions and the contact ellipse.
(a)(b)To obtain the minimum value of the induced normal curvature, the principal curvatures and principal directions at the LCP can be calculated firstly. The normal curvature in the arbitrary direction is given by(15)kn=Φ2Φ1=Ldθ2+2Mdθdϕ+Ndϕ2Edθ2+2Fdθdϕ+Gdϕ2,where Φ1 and Φ2 are the first and second fundamental forms of the surfaces, the coefficients (E, F, and G) of the first fundamental form are given as(16)E=rθ2=∂r∂θ2,F=rθ⋅rϕ=∂r∂θ⋅∂r∂ϕ,G=rϕ2=∂r∂ϕ2,and the coefficients (L, M, and N) of the second fundamental form are defined as(17)L=n⋅rθθ=n⋅∂2r∂θ2,M=n⋅rθϕ=n⋅∂2r∂θ∂ϕ,N=n⋅rϕϕ=n⋅∂2r∂ϕ2,where the unit normal to the surface is calculated by(18)n=rθ×rϕrθ×rϕ=∂r/∂θ×∂r/∂ϕ∂r/∂θ×∂r/∂ϕ.Equation (15) can be rewritten as follows:(19)knE−Lμ2+2knF−Mμ+knG−N=0,where μ=dθ/dϕ indicates the tangential direction of a surface corresponding to the normal curvature. Since the principal curvature is the maximum or minimum of the normal curvature, taking the derivative of equation (19) respect to μ and letting dkn/dμ=0, one can obtain the relationship between the principal curvatures and principal directions as(20)knE−Ldθ+knF−Mdϕ=0,knF−Mdθ+knG−Ndϕ=0.The equation of principal curvatures can be derived by eliminatingdθ/dϕ from equation (20):(21)EG−F2kn2−EN−2FM+GLkn+LN−M2=0.The calculated principal curvatures can be substituted into equation (19) to evaluate the corresponding principal directions. After getting the principal curvatures of both the gear (kg1 and kg2) and pinion (kp1 and kp2) tooth surfaces, as well as the corresponding principal directions (eg1, eg2, ep1, and ep2), the minimum value of the induced normal curvature can be calculated as follows: as shown in Figure 6(b), assume that the unit vectors in the direction of the major and minor axes of the instantaneous contact ellipse are w and u, respectively. w has an unknown angle α from the principal direction of the gear eg1. β is the angle between eg1 and ep1 and can be calculated by(22)β=arccoseg1⋅ep1eg1ep1.Utilizing the Euler equation, the gear and pinion normal curvatures alongw are given, respectively, as(23)kgw=kg1cos2α+kg2sin2α,(24)kpw=kp1cos2α+β+kp2sin2α+β.Hence, the induced normal curvature alongw is defined as(25)krw=kgw−kpw=kg1cos2α+kg2sin2α−kp1cos2α+β−kp2sin2α+β.Taking the extremum of the induced normal curvaturedkrw/dα=0, one can obtain two values of α. One of them is the angle between eg1 and w and the other is the angle between eg1 and u. The normal curvature of the gear (kgw and kgu) and the pinion (kpw and kpu) in the direction of w and u can then be obtained by substituting α into equations (23) and (24). The radii of curvature in the direction of w and u is the inverse of the corresponding normal curvature values:(26)Rgw=1kgw,Rgu=1kgu,Rpw=1kpw,Rpu=1kpu.In addition to the geometrical properties of the gear and pinion tooth surfaces at the LCP, the implementation of EHL analysis also requires the motion characteristics. As illustrated in Figure6(a), the velocity vectors of the gear and pinion at the ACP are given as(27)vg=ωg×rg,vp=ωp×rp,where ωg and ωp are the angular velocity vectors of the gear and pinion, respectively. The velocity vectors vg and vp can be resolved along the unit normal of the tooth surfaces and in the tangential plane T along the major and minor axes of the elliptical contact footprint as shown in Figure 7. These tangential components, which are used to calculate the entrainment and sliding velocities, can be obtained by the vector dot product:(28)vgw=vg⋅w,vgu=vg⋅u,vpw=vp⋅w,vpu=vp⋅u.Figure 7
Velocity relation.Therefore, the speed of entraining motion along the major and minor axes of the contact ellipse is defined as(29)vew=12vgw+vpw,veu=12vgu+vpu.The entrainment and sliding velocities can then be calculated, respectively, by(30)ve=veww+veuu,vs=vgT−vpT=vgww+vguu−vpww+vpuu.


## 3. EHL Basic Dimensionless Equations

The Reynolds equation describing the formation of hydrodynamic pressure for elliptical contacts with an angled entrainment flow in the steady-state condition is(31)εx∂∂Xρ¯H3η¯∂P∂X+εy∂∂Yρ¯H3η¯∂P∂Y=cosθ∂ρ¯H∂X+sinθk∂ρ¯H∂Y,with the following nondimensional variables: X=x/a, Y=y/b, H=h/a, P=p/pH, η¯=η/η0, and ρ¯=ρ/ρ0. Here, a and b are, respectively, the half-contact width (along the minor axis) and the half-contact length (along the major axis) of the Hertzian contact ellipse and can be calculated according to ref. [24]; pH is the maximum Hertzian pressure, and for elliptical contacts, pH=3w/2πab where w is the load carried by a pair of meshing teeth; η0 and ρ0 denote the ambient viscosity and density of the lubricant, respectively; θ is the angle of entrainment velocity relative to the minor axis of the contact ellipse. The dimensionless coefficients for Poiseuille flow εx and εy are given by(32)εx=apH12η0ve,εy=a3pH12b2η0ve=εxk2,where k is the contact ellipticity and defined by k=b/a. The dimensionless boundary conditions of the Reynolds equation are PX0,Y=PX,Y0=0 at the inlet and PXe,Y=∂PXe,Y/∂X=PX,Ye=∂PX,Ye/∂Y=0 at the outlet. The cavitation condition imposes PX,Y≥0.The non-Newtonian behaviour of the lubricant is not taken into account in the present work for simplicity. This simplification does not have a significant impact on the findings of this research because Newtonian models are suitable for the prediction of EHL film thickness, which are the index for evaluating lubrication performance used in this work, even under heavy-loaded conditions. It has been found in many previous studies that the entraining action that governs the formation of the lubricant film occurs in the inlet zone of the EHL contact where the pressure is relatively low, and thus the non-Newtonian behaviour of the lubricant is not apparent.For the purpose of investigating the impacts of changes in meshing parameters due to the contact point migration on analysis results, smooth surfaces are assumed to exclude the influence of surface topography. Therefore, the film thickness expression for the smooth elliptical contact is applied:(33)H=Hc0+aX22Reu+kbY22Rew+VX,Y,where Hc0 is the normal approach of the gear and pinion tooth surfaces and needs to be adjusted during the iterative process to make the calculation result meet the load balance equation and aX2/2Reu and kbY2/2Rew represent the undeformed contact geometry where Reu and Rew are, respectively, the equivalent radii along the minor (u) and major (w) axes of the contact ellipse. Applying the radii of curvature in the direction of u and w (calculated by equation (26)), one can obtain Reu and Rew by(34)1Reu=1Rgu+1Rpu,1Rew=1Rgw−1Rpw.It is important to note that “+” is used before 1/Rpu as the curvature directions of the gear and pinion tooth surfaces along u are opposite, whereas “−” is used before 1/Rpw since the curvature directions of the two surfaces along w are the same.The localized elastic contact deformationVX,Y is calculated by the Boussinesq integral:(35)VX,Y=2kpHπE′∬ΩPX′,Y′X−X′2+k2Y−Y′2dX′dY′,where E' is the reduced Young’s modulus of the two surfaces.Isothermal condition is assumed in this work, so the viscosity and density of lubricants are only influenced by pressure. Although Doolittle’s free-volume model has been proven to be the most accurate description of the real piezo-viscous behaviour at the typical EHL pressures [25], it is not convenient to use as several Doolittle–Tait parameters must be experimentally measured. In contrast, the Roelands equation is relatively easy to use and can give reasonable calculation results close to those from the free-volume model [26]. Therefore, the Roelands equation is adopted here to describe the pressure-viscosity relation for simplicity:(36)η¯=explnη0+9.67−1+1+5.1×10−9pHPZ,Z=αcplnη0+9.67,where Z is the pressure-viscosity index, cp=1.96MPa is a constant, η0 denotes the atmospheric viscosity of lubricants, and α is the pressure-viscosity coefficient. The variation of lubricant density with pressure is described by the Dowson–Higginson pressure-density relationship:(37)ρ¯=1+0.6×10−9pHP1+1.7×10−9pHP.The lubricant reaction force defined by the integral of the pressure over the whole solution domain should be balanced by the load acting on each pair of meshing teeth:(38)∬ΩPX,YdXdY=W=wabpH=2π3.In the above equations the unit of maximum pressurepH is Pa.


## 4. Numerical Method

The finite difference method is applied in the discretization of the Reynolds equation (AppendixA). In order to ensure the solution convergence and numerical stability under heavy-loaded conditions, the semisystem approach described in ref. [27] is adopted. In accordance with this method, the coefficient matrix of the discretized Reynolds equation is constructed utilizing both the pressure flow and the entraining flow terms, and thereby the diagonal dominance of the coefficient matrix can be guaranteed even when the pressure flow terms almost vanish under extremely heavy load. The discrete Reynolds equation is solved by Gauss–Seidel line relaxation. The full multigrid algorithm [28] is applied to accelerate the solution process and improve numerical accuracy. Surface contact deformation is calculated using the multilevel multi-integration technique [28] (Appendix B). The domain of computation is defined as −2.5≤X≤1.5 and −2≤Y≤2. Four layers of the grid are used, with the finest grid consisting of 257×257 nodes equally spaced. The next coarser grid has 128+1×128+1 nodes, etc. The finest grid has the dimensionless mesh spacing of ΔX=ΔY=0.015625, which meets the requirement of most engineering applications. The first target grid on which the full multigrid process is started should be sufficiently fine [28], and thus the grid having 64+1×64+1 nodes is selected. The pressure and load convergence criteria are used to determine whether the result converges on the target grid:(39)∑i∑jPi,jn−Pi,jn−1Pi,jn≤εp,2π/3−ΔXΔY∑i∑jPi,jn2π/3≤εw,where Pi,jn and Pi,jn−1, respectively, denote the new and old approximation of pressure in each grid point and the limit error εp and εw are set to be 10−5 and 10−3, respectively.Figure8 shows the flowchart of the present numerical procedure. Note that the normal approach Hc0 will also be adjusted on the coarsest grid according to the multigrid method [28], which is not illustrated in Figure 8.Figure 8
The flowchart of the numerical procedure.


## 5. Results and Discussion

A pair of mismatched spiral bevel gears manufactured by the Gleason face milling approach is treated. The pinion is the driving member. The drive sides of the gear pair, namely, the concave side of the pinion tooth and the mating convex side of the gear tooth, are the objects of interest. The basic design parameters are listed in Table1, with which the machine tool settings used to derive tooth geometry are obtained by the Gleason CAGE4win. The material properties of the gear pair and the lubricant properties are shown in Table 2.Table 1
Design parameters of the example spiral bevel gear pair.
ParametersPinion (left hand)Gear (right hand)Number of tooth1134Helical angle (rad)0.646Pressure angle (rad)0.349Pitch angle (rad)0.3131.257Pitch diameter (mm)71.5221Pitch cone distance (mm)116.139Width of tooth (mm)35Outer transverse module6.5Shaft angle (deg)90Table 2
Properties of contacting solids and the lubricant.
ParametersValueEffective elasticity modulus (E′)226 (GPa)Poisson’s ratio0.3Ambient atmosphere (p0)1.96 (MPa)Inlet density (ρ0)846 (kg/m3)Atmospheric viscosity (η0)0.135 (Pa⋅s)Pressure-viscosity coefficient (α)2.6e − 8 (Pa−1)A series of quasistatic isothermal EHL solutions are given at each meshing position over a tooth engagement cycle. The meshing position represented by the pinion roll angle commences at the beginning of the tooth engagement. To investigate the effect of torque loads on contact geometry and lubrication performance, we performed EHL simulations under two typical torque loads of 50 Nm and 300 Nm, respectively, representing lightly loaded and heavily loaded conditions. The rated speed of the driving pinion is 500 r/min. TCA and LTCA were, respectively, performed to obtain the theoretical unloaded contact path and loaded contact paths at 50 Nm and 300 Nm torque load. LTCA also provides variations of contact load (shown in Figure9), the coordinates of EMPs, and the direction cosine of the effective mesh force vector (listed in Tables 3 and 4) at each meshing position over the tooth engagement. The support of the example gear set is assumed to be rigid in LTCA, so the deformations of bearing and housing are not considered.Figure 9
Contact load variations over a tooth engagement cycle.Table 3
Results from LTCA at 50 Nm torque load.
Pinion roll angle (rad)EMP coordinatexm (mm)EMP coordinateym (mm)EMP coordinatezm (mm)Direction cosinenmxDirection cosinenmyDirection cosinenmz0.24399.7576−30.5724−3.6266−0.6314−0.16160.75840.27998.5905−30.3121−3.4352−0.6292−0.16120.76040.31597.5790−30.1337−3.0758−0.6287−0.16010.76100.35196.8110−30.0895−2.6075−0.6300−0.15850.76020.38796.1453−30.1450−2.0752−0.6325−0.15670.75850.42395.5419−30.2295−1.5013−0.6355−0.15510.75640.45994.9332−30.2986−0.9084−0.6384−0.15370.75420.49594.2516−30.3383−0.3472−0.6411−0.15240.75210.53193.5125−30.31300.0996−0.6430−0.15170.75070.56792.5259−30.13820.3517−0.6429−0.15180.75080.60391.2030−29.73990.3843−0.6401−0.15300.75290.63989.6923−29.15380.2597−0.6352−0.15500.75660.67587.9616−28.53990.1086−0.6299−0.15690.76060.71186.6346−28.07200.0663−0.6265−0.15770.76330.74785.5889−27.81060.2067−0.6257−0.15740.76400.78384.7482−27.69270.5181−0.6271−0.15650.76310.81984.0397−27.63410.9342−0.6294−0.15550.76130.85583.3518−27.61261.3925−0.6323−0.15460.75920.89182.8174−27.64941.8746−0.6357−0.15350.75650.92782.3046−27.74192.3997−0.6400−0.15210.75320.96381.8595−27.84712.9825−0.6445−0.15050.7496Table 4
Results from LTCA at 300 Nm torque load.
Pinion roll angle (rad)EMP coordinatexm (mm)EMP coordinateym (mm)EMP coordinatezm (mm)Direction cosinenmxDirection cosinenmyDirection cosinenmz0.171101.923−31.0176−4.3476−0.6337−0.16990.75470.207100.879−30.8057−4.0301−0.6324−0.16690.75640.24399.9242−30.6475−3.6532−0.6320−0.16410.75740.27999.0651−30.5566−3.2164−0.6324−0.16150.75760.31598.2943−30.5335−2.7238−0.6338−0.15920.75690.35197.5905−30.5640−2.1880−0.6358−0.15720.75560.38796.9235−30.6219−1.6305−0.6383−0.15550.75390.42396.2570−30.6751−1.0807−0.6408−0.15410.75210.45995.5506−30.6918−0.5697−0.6430−0.15310.75050.49594.7658−30.6475−0.1247−0.6444−0.15260.74930.53193.8746−30.52920.2380−0.6451−0.15240.74880.56792.8700−30.33740.5163−0.6449−0.15260.74890.60391.7699−30.08420.7222−0.6440−0.15300.74960.63990.6137−29.78990.8781−0.6426−0.15350.75060.67589.4492−29.47831.0113−0.6411−0.15380.75180.71188.3185−29.17201.1476−0.6398−0.15390.75300.74787.2473−28.88941.3072−0.6389−0.15370.75380.78386.2443−28.64361.5036−0.6386−0.15340.75410.81985.3090−28.44261.7442−0.6390−0.15290.75390.85584.4419−28.29052.0331−0.6400−0.15240.75310.89183.6497−28.18862.3731−0.6418−0.15170.75170.92782.9418−28.13622.7653−0.6443−0.15080.74970.96382.3190−28.12933.2080−0.6474−0.14980.74730.99981.7642−28.16013.6929−0.6509−0.14840.74451.03581.2408−28.21504.2013−0.6547−0.14690.74141.07180.7052−28.27484.7043−0.6586−0.14530.73841.10780.1340−28.31695.1663−0.6621−0.14400.7355Figures10(a) and 10(b), respectively, compare the contact paths (composed of the EMPs at each meshing position in the projection plane) at 50 Nm and 300 Nm torque loads with the theoretical contact path under no-load condition on the concave side of the pinion tooth and the convex side of the gear tooth. It is shown that the EMPs for both the pinion and the gear move from the heel to the toe of teeth throughout the engagement. Meanwhile, for the pinion, the EMP moves from the tooth root to the tooth top, but for the gear, on the contrary, it shifts from the tooth top to the tooth root. It can also be observed that the increase in load results in the extension of the EMP footprint toward the heel and the toe of teeth. As shown in Figure 10, the contact paths at 50 Nm and 300 Nm torque loads are longer than the theoretical unloaded contact path, and among them, the contact path at 300 Nm is the longest. It indicates that under loaded conditions, the tooth pair comes into meshing earlier and secedes meshing later than the unloaded case due to the deformation.Figure 10
Contact paths on (a) the concave side of the pinion tooth and (b) the convex side of the gear tooth.
(a)(b)Utilizing the data in Tables3 and 4, one can calculate the geometrical and kinematical parameters for EHL analysis according to the previously described approaches. Figure 11 displays the major axis direction of the contact ellipse at the LCP of each meshing position calculated by the proposed method. The contact pattern, which is the envelope of the ends of the contact ellipse predicted from LTCA shown in Figure 12, is also drawn in Figure 11 for comparison. It can be seen that the obtained major axis directions of each elliptical contact footprint in the contact pattern are in good agreement with those from LTCA. This indicates that the proposed method is capable of evaluating the contact geometry at the LCP while ensuring that the calculated contact footprints are consistent with the footprints predicted from LTCA.Figure 11
Comparison of contact footprints predicted by the proposed method and those from LTCA at (a) 50 Nm and (b) 300 Nm torque load.
(a)(b)Figure 12
Instantaneous contact ellipses from LTCA.Figures13 and 14, respectively, show the variations of the radii of curvature and entraining velocities in the direction of the major and minor axes of the contact ellipse over one tooth engagement. They are the input geometric and kinematic meshing parameters of EHL analysis. It can be seen that since the LCPs do not significantly deviate from the TCPs at light load, as shown in Figure 10, the variation trend of meshing parameters over one tooth engagement under light-loaded conditions is similar to their theoretical change trend. In contrast, at heavy loads, a more noticeable difference can be found. Besides, in different meshing positions shown in Figure 8, the deviation direction of the LCP from the TCP is not consistent, so no obvious change rule of meshing parameters with load can be found from Figures 13 and 14.Figure 13
Variations of the radii of curvature of the mating gear and pinion tooth surfaces along the major and minor axes of the elliptical contact footprint.
(a)(b)(c)(d)Figure 14
Variations of the entraining velocity along the (a) major and (b) minor axes of the elliptical contact footprint.
(a)(b)Both the meshing parameters obtained at the LCPs and the corresponding values calculated at the TCPs are used in the EHL simulation to observe the error involved by ignoring the adjustment of contact points due to deformation. Note that in the areas where the meshing parameters are only presented under loaded conditions as shown in Figures13 and 14, if the gear pair is considered to be a rigid body, the tooth pair of interest is out of contact, while the adjacent tooth pair is in contact. Therefore, only the calculation results in regions where the tooth pair of interest is in contact under both loaded and unloaded conditions are compared. Figure 15 shows the comparison of the central and minimum film thickness. It demonstrates that, at 50 Nm torque load, the results without considering the effect of deformation on the contact point are similar to those calculated at the LCPs. The maximum differences, which occur at the pinion roll angle of 0.711 rad, are about 6.4% for the central film thickness and 6.5% for the minimum film thickness. However, these differences become significant under the torque load of 300 Nm, reaching to 13.2% for the central film thickness and 15.2% for the minimum film thickness. Figure 16 compares the pressure distributions and film thickness profiles along the major and minor axes of the contact ellipse at the pinion roll angle of 0.711 rad. Though there is a clear difference between the film thickness predicted at the LCP and the TCP, the difference in the corresponding primary pressure peak is relatively small. Recall that although the gear tooth contact geometry at the LCP is different from that at the TCP, the loads carried by the meshing tooth pair at these two positions are the same (751.2 N and 4416 N at the applied torque of 50 Nm and 300 Nm, respectively, shown in Figure 9). The half-contact width a and half-contact length b of the Hertzian contact ellipse in the four cases are 0.171 mm and 3.220 mm (LCP, 50 Nm), 0.168 mm and 3.014 mm (TCP, 50 Nm), 0.313 mm and 6.166 mm (LCP, 300 Nm), and 0.304 mm and 5.758 mm (TCP, 300 Nm). Hence, the corresponding maximum Hertzian pressures pH are 0.65 GPa (LCP, 50 Nm), 0.7 GPa (TCP, 50 Nm), 1.1 GPa (LCP, 300 Nm), and 1.2 GPa (TCP, 300 Nm). It can be found that the difference in the maximum Hertzian contact pressure caused by the difference in the tooth surface geometry between the LCP and the TCP is not very large. This is why the primary pressure peak calculated at the LCP does not differ from that obtained at the TCP significantly. However, in addition to the tooth contact geometry, the film thickness also largely depends on the entrainment velocity in the direction of the minor axis of the contact ellipse, which has been found in many previous studies. As can be seen in Figure 14, the entrainment velocity along the minor axis of the contact ellipse calculated at the LCP is higher than the one obtained at the TCP. This may be the dominant reason why the film thickness calculated at the LCP is greater than the oil film thickness calculated at the TCP. Overall, the results indicate that the contact point migration, which may not significantly affect pressure distribution, has nonnegligible effect on the prediction of film thickness. An interesting result is that if the effect of gear deformation on the position of the contact point is taken into account, the lubrication performance of the example gear pair is improved to some extent after loading. For example, at the pinion roll angle of 0.711 rad where the load carried by the meshing tooth pair reaches to the maximum value as shown in Figure 9, the film thickness predicted at the LCP is larger than the value evaluated at the TCP. This indicates that owing to the deformation, the actual contact point moves to a position that is favorable for lubrication.Figure 15
Comparison of central and minimum film thickness at (a) 50 Nm torque load and (b) 300 Nm torque load.
(a)(b)Figure 16
Comparison of pressure distributions and oil film thickness profiles along the minor and major axes of contact ellipses at the pinion roll angle of 0.711 rad under (a, b) 50 Nm and (c, d) 300 Nm torque load.
(a)(b)(c)(d)


## 6. Conclusion

The present study provides a solution of the isothermal elliptical contact EHL for spiral bevel gears, considering the migration of contact points of the meshing tooth pair due to load-induced deformation. The loaded contact points at each meshing position during the tooth engagement are determined by the loaded tooth contact analysis. A method for calculating the geometric and kinematic meshing parameters at the loaded contact points is proposed.At high loads, the loaded contact points of the deformed gear pair deviate significantly from the theoretical contact points of the corresponding rigid body. If the effect of deformation on the position of the contact point is not taken into account, it may cause a significant error in the EHL simulation results, especially the prediction of film thickness. Therefore, under tremendous loads, it is recommended to use the meshing parameters at the loaded contact points for EHL simulation. Changes in contact points due to gear deformation may improve or worsen lubrication performance. Thus, obtaining the lubrication performance at the loaded contact points under the actual working conditions will help to guide the tooth surface modification.



---

*Source: 2047876-2020-06-22.xml*