Exploration of trade-offs between NVH and efficiency in bevel gear design

Minimizing NVH and friction-induced power losses is becoming paramount in the design of geared transmissions. The aim of this paper is to present an automatic methodology to explore Pareto-optimal designs of bevel gears when minimization of noise and frictional losses is essential. In the first part, a semi-empirical model to estimate frictional power losses under elasto-hydrodynamic lubrication is described. The model has been validated against experimental data available in the literature in previous works by the authors. The efficiency calculation is coupled with a state-of-the-art loaded tooth contact analysis (LTCA) tool to obtain accurate predictions of the instantaneous load shared by the mating tooth pairs during the meshing cycle. In the second part, an automatic framework based on multi-objective optimization (MOO) is presented where the tooth micro-geometry is systematically designed. The design variables are represented by few coefficients of a polynomial basis that embodies the tooth flank ease-off topography. To ensure manufacturability, the polynomial modifications are projected onto the feasible set of the machine-tool envelopes. This step is achieved through a state-of-the-art identification algorithm that the authors have developed in previous work. Frictional losses are estimated with the aforementioned model, whereas the NVH level is measured by the loaded transmission error (LTE), directly available from the simulation tool. The maximum contact pressures are limited by the material properties, thus proper nonlinear constraints are prescribed. Application to a test case involving the design of a spiral bevel gearset reveals that the methodology presented allows the designer to obtain Pareto-optimal solutions in a systematic and automatic manner.


Introduction
Spiral bevel gears are extensively utilized in diverse industrial applications owing to their capability to efficiently transmit high torque, operate smoothly, and offer high reliability in power transmission.Nevertheless, the deployment of spiral bevel gears in high-end applications is confronted with some design challenges, such as minimization of efficiency losses, noise, vibration and harshness (NVH).
The efficiency of spiral bevel gears may be impacted by various factors, which can be categorized into two types: load-independent and load-dependent losses.The former encompass churning and windage losses, which refer to the lubricant being pumped between the mating members and its splashing caused by inertial effects.The latter pertains to rolling and sliding frictional losses, with the former being negligible in comparison to the latter, particularly in geared transmissions.Our study specifically concentrates on sliding frictional losses, which are primarily influenced by the sliding motion between the teeth, gear mesh stiffness, and tooth modifications.
The core aspects involved in correctly modeling the mechanical losses are: 1. loaded tooth contact analysis (LTCA) to predict load sharing between the mating teeth and the pressure distribution over the contact zones; 2. calculation of the geometric and kinematic properties of the contacting zones (i.e., entraining velocity of the lubricant, sliding velocity, and surface curvatures);

estimation of the elasto-hydrodynamic lubrication (EHL) friction coefficient;
Numerous models and methods exist in the literature for evaluating friction-induced losses.The primary differences between these models and methods lie in the semi-empirical formulas used to estimate the EHL friction coefficient and/or the tools used to perform the LTCA.
Li et al. [1] proposed an EHL friction model for helical gears.In this model, the contacting bodies are sliced along the entire contact zone, and a line contact model is postulated on each slice.The authors noted a conflicting trend on tooth flank modifications between mechanical efficiency and other transmission specifications, such as maximum pressure and loaded transmission error (LTE).
Xu et al. [2] applied a similar EHL friction model to estimate efficiency losses in hypoid gears.In this model, a large set of results from EHL analyses were fitted to obtain a semi-empirical coefficient formula, while LTCA was carried out using the same simulation tool employed in this paper (Calyx, [3]).However, their EHL model did not account for mixed and boundary lubrication conditions.
Kolivand et al. [4] started from the same EHL friction model and improved it to also account for mixed and boundary lubrication, with the aim of investigating efficiency losses in hypoid gears.However, the LTCA tool was replaced by the Hypoid Analysis Program (HAP) [5], which uses a semi-analytical shell theory to compute the tooth compliance and an ease-off-based approach to perform LTCA analyses.If compared to the former, HAP K trades computational efficiency for accuracy, especially at high loads.
Cao et al. [6] evaluated the lubrication performance under different possible contact paths on bevel gears.Mohammadpour et al. [7] carried out an EHL analysis on hypoid gears under relatively high loads, while Pu [8] analyzed the lubrication of hypoid gears taking into account the threedimensional surface roughness.However, despite their accurate tribological analyses, many of the cited contributions fall short in providing reliable contact analysis results.In fact, an accurate LTCA tool has a paramount importance for the evaluation of the contact pressures and the subsequent estimation of the friction coefficient.Moreover, most of the cited references do not offer precise estimations of the lubricant's entraining velocity, often relying instead on calculating it as the average of the absolute velocities of the mating pairs at their contact points.
Ziegltrum et al. [9] provided an accurate estimation of the local friction coefficient of spur gears under transient thermo-elasto-hydrodynamic lubrication (TEHL) conditions using FE-based multiphysics software.Their results were also validated against experimental data for different lubricants.
Paschold et al. [10] mathematically modeled the heat transfer during meshing of worm gears in order to accurately investigate their overall efficiency losses.
NVH issues in spiral bevel gears can largely be attributed to gear meshing errors, manufacturing errors, and misalignments.These problems may result in increased noise levels, vibration, and reduced reliability of the gear system.Consequently, understanding the underlying causes of NVH problems and efficiency losses in spiral bevel gears and developing strategies to alleviate them is crucial to enhance their performance and reliability.In this study, improvements in NVH are dealt with by minimizing the harmonic amplitudes of the LTE through appropriate ease-off modifications.
The relationship between LTE and micro-geometry has been extensively demonstrated in the literature.Simon [11] demonstrated how adjusting higher-order motions associated with the cradle radial setting and the modified roll could decrease LTE.Similarly, Su et al. [12] proposed a methodology to design a seventh-order transmission ratio function based on rigid TCA assumptions.A numerical test case demonstrated that this method could provide advantages in terms of transmission error under load.However, the effective values of the parameters representing the seventh-order function are assumed to be known, and no guidelines are provided for deriving them.Along the same lines, [13] proposed an optimization framework to minimize LTE.Here, higher-order motions associated with the coefficients of the modified vertical and helical motions are identified during the optimization process.
Several studies in the literature have proposed automatic optimization methods to improve the performance of spiral bevel gears during meshing.Mermoz et al. [14] used accurate finite element analyses to carry out LTCA and determine the contact pressure, which then served as the objective function minimized by the BOSS/Quattro commercial optimization software.Similarly, Astoul et al. [15] also used this methodology to minimize quasi-static transmission errors.However, both studies directly optimize a limited subset of machine settings, specifically the modifiedroll coefficients, which may not be ideal when starting from a basic design that exhibits poor contact patterns.
Likewise, a recent study by Li et al. [16] uses finite element analyses to perform multi-objective optimization, wherein efficiency, contact pressure, and LTE are simultaneously minimized.Machine settings are selected as the optimization variables following a sensitivity analysis.
However, it should be noted that the presented optimization test case is somewhat questionable, as the LTCA is carried out with just 10 Nm pinion torque (a very lightly loaded tooth contact analysis).
Notably, the Particle Swarm Optimization-Gravitational Search Algorithm (PSOGSA) [17] and the Fast Elitist Nondominated Sorting Genetic Algorithm (NSGA-II) [18] have been effectively utilized for optimizing and exploring Pareto-optimal solutions in spiral bevel and hypoid gears.However, it is widely recognized that these algorithms are computationally expensive, especially in terms of number of function evaluations to convergence (and the number of iterations to convergence is often set a priori).Therefore, to obtain solutions within reasonable time frames, these solvers still need to be coupled with relatively simple and approximate contact models.
The present study involves the automated optimization of tooth flank topography using ease-off polynomial functions to explore potential trade-offs between NVH and efficiency.Initially, an EHL friction model is briefly presented, where the entraining velocity is rigorously evaluated to estimate the shear strain rate within the film.This model has already been introduced and validated against experimental data in [19,20], specifically for a hypoid gearset with a 75W90 oil, but further improvements are described in the present work.Other losses, such as windage, churning, and rolling resistance, are excluded due to their limited dependence on the tooth flank topography.The study uses the MIL-L-23699 lubricant, which is a neopentyl polyol ester base oil commonly used in aerospace applications such as gas turbine engines and aircraft gearboxes.
Subsequently, the multi-objective optimization framework previously described in [21] is briefly reviewed.To ensure accurate LTCA analyses, the Calyx contact solver is employed, although the methodology proposed in this work may well be complemented by other contact analysis tools (i.e.Becal [22] or other general-purpose FEM software).The optimizations are carried out by a relatively simple direct-search algorithm (MATLAB's patternsearch).An important characteristic of this innovative framework is that, despite utilizing ideal polynomial functions to steer the optimization process, a tailored identification algorithm is applied at each iteration to project the ease-off into the set of deviations that can be actually manufactured, as presented in [21].This approach offers the benefit of driving the expensive optimization loop with a minimal number of variables, while determining, with the embedded identification algorithm, the most appropriate manufacturable flank modifications from a relatively extensive set of machinetool settings.
To determine the boundaries of the Pareto front, singleobjective optimizations are initially carried out.Then, the Pareto front is explored by solving additional multi-objective optimization problems.The contact patterns for all points sampled on the Pareto front are also shown.

EHL friction
There are a number of factors that influence friction, such as the lubrication regime (full film, mixed, boundary), the behaviour of the lubricant with varying operating conditions (temperature, pressure, shear rate) and the surrounding environment (the boundary conditions of the lubricated contact).
An EHL friction coefficient model similar to the one described in [19] is proposed here.The model is briefly reintroduced in this section, since there are some differences compared to the previous implementation, particularly in the load sharing function and in the rheological model of the lubricant.
The friction coefficient f is determined based on a load sharing factor between boundary and full-film lubrication conditions, which is described by the function X. /.It is well known in the literature that f can be calculated using the following equation: Here f h and f b are the friction coefficients related to the hydrodynamic (full fluid) and the boundary lubrication conditions, respectively.is the ratio between the film thick-ness h and the equivalent root-mean-square roughness of the contacting body surfaces: Usually, the central film thickness h c is used, and R q1 and R q2 are the root-mean-square roughnesses of the two surfaces.
Several expressions can be found for the load sharing factor X. / in the literature.Differently from what has been described in [19], the recent findings by Taylor and Sherrington [23] are employed in this work.Their study suggests a reverse S-shaped curve for X. /, as represented by the following formula: According to [23], and based on experimental data from several lubricant oils, the recommended values for k and a here employed are 1.453 and 1.32, respectively.The coefficient of boundary friction f b is considered constant, and it has to be experimentally derived.Based on data available in the literature [23], its value ranges from 0.1 to 0.14, depending on the employed lubricant oil.In this work, a worst-case value of f b = 0.14 is assumed.
The mean value of f h can be obtained by dividing the shear stress, , by the mean contact pressure p m : In our framework, p m is obtained directly from contact analysis tool for each contact zone.The shear stress is estimated according to the Bair-Winer rheological model: Here, L is the limiting shear stress, Á is the dynamic viscosity, and the shear strain rate is approximated by P = v s h c , with u being the sliding velocity between the two mating teeth at their nominal contact point.
The limiting shear stress L is derived for specific pressure and temperature conditions according to the following exponential relation [24]:

K
The dynamic viscosity Á variation with temperature T and pressure p is described by the Roelands and modified Barus formulas: The viscosity-pressure coefficient ˛can be estimated as follows: To estimate the remaining coefficients in Eqs. ( 8), ( 7) and ( 9), empirical data specific to the lubricant being used is necessary.Table 1 contains the relevant data for the MIL-L-23699 lubricant.Different expressions for the evaluation of the central film thickness can be found in the literature.Four different lubrication regimes can be observed for nonconformal contacts, depending on the elastic deformation of the bodies and the variation of viscosity with pressure.The four regimes are usually indicated as isoviscous-rigid IR, piezoviscous-rigid PR, isoviscous-elastic IE and piezoviscouselastic PE [25].Bassani et al. [26] provided semi-empirical formulas for the central film thickness for the more general case in which the entraining velocity is not collinear with any principal direction of the curvatures: Here, E is the equivalent elastic modulus, u is the entraining velocity, and R e and R s are the equivalent radii of curvature, respectively, parallel and perpendicular to u.The calculation of those quantities is detailed in the next subsection.The proper lubrication regime is determined in a practical way by taking the highest of the values given by Eqs. ( 10)- (13).
Thermal effects are included using a reduced value for the central film thickness.This is obtained via scaling it by a dimensionless reduction factor ˚, for which several models are available in the literature.The one described in [27] is employed in our model: where S = u=u and L = ˇÁ0 u 2 k is the dimensionless thermal loading parameter.Their values can be elicited from Table1.
It should be emphasized that the equations introduced within this section serve as preliminary estimations of the lubricant film thickness and coefficient of friction.Of noteworthy consideration, the thermal effects linked to high sliding speeds are currently unaccounted for in the given formulas.

Kinematic relationships at the contact point
In order to use the formulas presented in the previous subsection, it is necessary to calculate the entraining velocity u and the slide-to-roll ratio S. The entraining velocity of the lubricant is the relative velocity of the surface with respect to the contact point, or vice versa.The overall entraining velocity of the lubricant is the average of the entraining velocities of the two mating surfaces.
Previous studies [20,28] have presented formulas for the entraining velocity that presume the availability of a local orthogonal parametrization of surface coordinates [29].This is usually too strong of an hypothesis and certainly not very convenient when dealing with spiral bevel (and hypoid) gears.Here, we set out to generalize the calculations of the entraining velocity.To this sake, we first introduce the general equations that describe the contact point Given c 1 = c 1 .˛1/and c 2 = c 2 .˛2/ in tangent contact, the following two equations are satisfied: Here, ˛1; ˛2 2 R 2 represent the surface coordinates and d 1 and d 2 mark the position of a generic point on the axis of each mating surface.The first equation expresses coincidence of the two surfaces at the contact point, while the second is the equation of meshing.The symbol n 1 denotes the unit normal vector of surface 1 at the contact point, while b ! 1 and b ! 2 are the spatial angular velocities (with respect to the fixed frame S 0 ) in their skew-symmetric matrix representation.
To determine the contact point velocity, we differentiate the system of Eqs.(15) with respect to time, which results in the following linear system of equations in P ˛1 and P ˛2: where and  16) allows us to obtain the velocity of the contact point as it moves over the two surfaces as follows: The overall entraining velocity then becomes: The surface contact points c j , tangential vectors c j;˛j , and normals n j (j = 1,2) required in Eq. ( 16) can be obtained at each time step and contact zone through post-processing of the LTCA data.Additionally, the principal curvatures K x and K y and their corresponding directions x and y for the equivalent contact can also be extracted.

K
Finally, to obtain the radius of curvature R e along the entraining direction and along the side-leakage direction R s , the angle Â between x and u e can be determined by: which can then be used to transform the principal radii of curvature:

Frictional power losses estimation
The friction coefficient computed using Eq. ( 1) represents the mean coefficient of friction for a mating tooth pair over the instantaneous contact zone.Although for the determination of such friction coefficient the sliding and entraining velocities are evaluated at the mean contact point, power loss estimation requires that the sliding velocity field is considered over the entire contact zone.To achieve this, sliding velocities and normal loads applied to each cell of the computational grid are extracted from the simulation post-process data.The computational grid is the local surface mesh employed by the FEA solver to compute contact pressures (through the Boussinesq solution for an elastic half-space), as shown in Fig. with a friction coefficient f ik , is then estimated as the sum of the products of the sliding velocities v s ikl and the normal loads F ikl for all N k cells within the contact zone:

Optimization and Pareto-front exploration
This section introduces the framework for optimizing the micro-geometry of face-milled spiral bevel gears.The optimization pipeline for the micro-geometry is depicted in Fig. 4, outlining the steps presented in the next sections.
The framework utilizes MATLAB interfaced with a FEAbased contact solver to launch an analysis and subsequently post-process the results, as initially described in [21].
The main components of the pipeline are the following: 1. definition of the initial geometry and the optimization variables; 2. identification of the machine-tool settings which generate the closest flank modifications described by the polynomial ease-off; 3. description of the interface between the optimizer and the LTCA simulation; 4. post-processing of the LTCA simulation results for the calculation of objective functions and constraints.
The LTCA solver, by being finite element-based, is known to be computationally expensive, especially when employed in automatic optimizations.Furthermore, numerical noise caused by the discretization of the underlying problem would hamper its direct use in gradient-based optimizers.
To overcome both limitations at once, the (direct search) patternsearch algorithm, from the MATLAB's Optimization Toolbox, has been employed with parallel calls to LTCA simulations.At each iteration, a complete poll is performed in parallel, exploring the variable space along all its directions, thereby significantly reducing the optimization time.It is important to note, however, that patternsearch is designed as a local minimizer.During the exploration of the Pareto-front within a multiobjective optimization scenario (presented in the numerical section), our expertise guided the choice of initial guesses to ensure that a satisfactory solution is obtained by the optimizer, particularly for the initial single-objective minimizations.Subsequently, the initial guesses were provided by averaging the solutions found during the Pareto-front exploration.

Initial gear geometry generation and optimization variables
The surface geometry generation is a well-established topic in the gear literature.Pioneered by Litvin [30], mathematical tools to simulate the envelope surface between the grinding tool and the working blank, borrowed from differential geometry, have been extensively described and applied.Nonetheless, to make this step more systematic and elegant, the Lie Group method borrowed from robotics [31] has been employed to parameterize the face-milling gear generation.
In order to clearly describe the machine-tool setting identification algorithm, it is necessary to reintroduce the fundamental mathematical concepts involved in describing the of tool motions and the well-known equation of meshing.Fig. 5 shows the cross section profile of the grinding wheel.The tool surface points are defined by p t .xt ; ; Â/ and the associated normals by n t .xt ; ; Â/, where .; Â/ are the surface variables and x t are the tool settings that define its shape.The envelope family of tool surfaces (and its normals) during motion can be expressed as follows1 : p g .x;/ = g gt .xm ; '/p t .xt ; ; Â/ n g .x;/ = g gt .xm ; '/n t .xt ; ; Â/ ; The homogeneous transformation matrix g gt .xm ; '/ 2 SE.3/ maps the tool reference frame fT g to the gear blank frame fGg, according to the Gleason UMC hypoid generator shown in Fig. 6, and ' represents the cradle rotation.
The vector x m includes all the machine settings.The points of the tool family belong to the envelope surface if the equation of meshing is satisfied: f .x;/ = n t .xt ; ; Â/ T b V t gt .xm ; '/p t .xt ; ; Â/ = 0.
Here, b V t gt .xm ; '/ 2 se.3/ is the homogeneous form of the body twist that describes the rigid-body velocity of fT g with respect to fGg, with components expressed in frame fT g.With = Â ' T , we compactly represent the enveloping triplet.
To optimize the micro-geometry, ease-off modification surfaces are described by polynomials p e .u;vI c/ expressed as: p e .u;vI c/ = S.u; v/ T c: (28) Specifically, shape functions S.u; v/ [32] are used to construct a quartic polynomial.A vector c is constructed from nine coefficients, each corresponding to a node, as shown Polynomial ease-off used to drive the optimization, represented in the zR axial plane of the blank.The optimization variables c represent the nodal values of the quartic element employed to describe the ease-off in Fig. 7.This vector is used to define the shape of the tooth flank modification, and it is employed as the vector of optimization variables.However, only eight of the nine coefficients are actually used as optimization variables: the 9-th node is kept fixed (at zero) to avoid tooth thickness variations.
After sampling the polynomial over a 11 22 grid of N = 242 points, the ease-off values are mapped onto the zR axial plane of the tooth flank, as in Fig. 7.
The ease-off h i at the i-th point of the grid then defines a target grid point according to the following equation: p i = p g .x;i / + h i n g .x;i /; .i= 1; :::; N / (29) This target grid is then used to identify the machine-tool settings x that generate the best-fit envelope surface, as described in the following section.

Machine-Tool Settings Identification
In a previous study [21], we demonstrated the benefits of using a sparse formulation variant of [33] to describe the nonlinear programming (NLP) problem of identifying the machine-tool settings required to obtain a given tooth flank modification.This approach offered benefits in terms of both robustness and computational efficiency.The computational efficiency was so remarkable that it allowed the identification step to be integrated directly into the optimization loop, immediately preceding the LTCA simulation.
We briefly introduce here the algorithm formulation.The identification is framed as the following constrained nonlinear least-squares problem (NLP):

<
:::: p. i ; x/ + n. i ; x/h i − p i = 0 f .i ; x/ = 0 ::: .i = 1; :::; N / (30) where: x = x t x m T is the vector of machine-tool settings; of the pinion tooth surface; h is the vector of the residuals (deviations) between the basic flank surface and the target surface; h i defines the i-th component (positive for material removal); p. i ; x/ and n. i ; x/ are the generic position and normal vectors of the family of the tool surfaces during motion; f .i ; x/ is the equation of meshing associated with the i-th point.
The CasADi framework [34] is used to create efficient symbolic expressions for the NLP problem and its gradients.
The interior-point algorithm IPOPT [35], directly interfaced with CasADi, is employed as back-end solver.

MOO problem
The micro-geometry design is formulated as a multi-objective optimization problem according to the method described in [36]: where f = f 1 f 2 T is the vector of the two nonlinear objectives that need to be simultaneously minimized, and c are the polynomial coefficients introduced previously.In the present paper, the trade-offs between NVH performances and efficiency losses are sought.Therefore: x .c///,root mean square of the first three harmonics of the loaded transmission error (LTE) (in rad); f 2 .x.c//= e L .x .c//,efficiency loss (in percentage points).
The following nonlinear constraint needs to be satisfied: x .c//= max.proot ; p tip ; p heel ; p toe / − p edge Ä 0 The constraints ensure that edge-loading is avoided by specifying that the contact pressure values at the tooth edges (root, tip, heel, and toe) must not exceed a certain threshold p edge .The maximum contact pressure is not accounted for in the objective vector.However, to take account of the material limitations in terms of pitting load carrying capacity, its maximum value is restrained by the tolerance p tol in the second nonlinear constraint of Eq. (32).
To guarantee manufacturability of the tooth flank modifications, the previously described identification is performed before the LTCA step.This is highlighted by the dependency of the objectives and constraints on x .c/:these are the identified machine-tool settings x associated with the ease-off represented by specific c values.As a matter of fact, while the optimization is guided by the polynomial coefficients c, only the identified machine-tool settings x are fed back to the LTCA software to re-generate the geometry and to carry out the analysis, at each iteration.
The approach described in [36] is adopted to reformulate the problem as a single-objective optimization problem via a scalarizing achievement function s W R 3 R 3 !R and to enforce the constraints through an exact penalty formulation: min c2B s.f.x .c//Ie f/ + 2 X k=1 max.g k .x.c//;0/ (33) where E e f is the selected reference point and is the penalty coefficient.The design variables are bounded to a reasonable set B to better guide the optimization solver.The following non-differentiable achievement function is employed: In Eq. ( 34), f (id) i and f (nad) i represent the ideal and the nadir values of the objectives.They are respectively the lower and the upper bounds of the Pareto front and can be obtained by minimizing the two objectives individually.Importantly, they allow to properly normalize the objectives.A small augmentation parameter 1 (usually 10 −4 ) is used to avoid obtaining weakly Pareto-optimal solutions.
Compared to the more traditional weighting methods used in similar optimizations, Eq. ( 34) allows non-convex portions of the Pareto front to be explored.For further details, the reader is referred to [36].
It is important to remark that the optimization methods proposed here revolve on simulating only the nominal operating conditions in terms of speed, torque and misalignments.Different load cases may significantly steer the corresponding misalignments which can result in sub-optimal meshing conditions.To properly account for different operating loads, further steps are necessary, as highlighted for example in [37].

Numerical test case
The spiral bevel gear set used in this study is described in Table 2, which provides the main data about its macrogeometry.The optimization focuses on the operating conditions of the drive side and, as such, the ease-off functions are utilized to modify the pinion's concave flank.
Before proceeding to the optimization, the machine-tool settings are derived to define a virtual conjugate gearset, which operates under the specified misalignments, according to the methods outlined in [21].The machine-tool settings that are associated with the concave flank of the conjugate pinion are listed in Table 3.The optimization variables are the c coefficients described previously.The following machine-tool settings are selected to identify manufacturable tooth flank modifications prior to the LTCA step: all the tool settings excluding: edge radius, Toprem angle, Flankrem angle (see Fig. 5); the following zero-th degree coefficients: R 0 , q 0 , E 0 , S 0 and D 0 ; modified roll coefficients up to the 4-th order: m, C , D, E; modified radial motion coefficients up to the 4-th order: R 1 , R 2 , R 3 , R 4 ; a b c Fig. 8 Optimal solution corresponding to the minimum of e L .a Contact pressure pattern: the sliding velocity pattern is overlayed to highlight the concentration of contacts near the instantaneous axis of rotation.b Optimal manufacturable ease-off (achieved via new machinetool settings).c Residual error after identification of the polynomial ease-off modified helical motion coefficients up to the 4-th order: Since the ease-off function used for optimization is a quartic polynomial, using machine settings associated with higherdegree motions may not be worthwhile.
The constraints are defined to enforce the maximum contact pressure to be below 1600 MPa (i.e., p tol = 1600 MPa) and below 450 MPa at the tooth edges (i.e., p edge = 450 MPa).The requirement of minimizing the LTE does not map distinctively into a predictable contact pattern shape.b Optimal manufacturable ease-off (achieved via new machine-tool settings).c Residual error after identification of the polynomial ease-off

Single objective minimizations
The two single-objective minimizations are carried out first.This step is fundamental when dealing with objectives with different units of measurements and/or different scales.
The solution obtained by minimizing the efficiency loss e L , shown in Fig. 8, features a contact pattern (Fig. 8a) that is localized around the trace of the instantaneous axis of rotation on the tooth surface.Further localization of the contact patterns is limited by the maximum allowable contact pressure p tol = 1600 MPa.It is worth noting that a peak value p max = 1595 MPa was registered, which indeed complies with the introduced pressure constraint.
Fig. 8b shows a contour plot of the optimal manufacturable ease-off, defined as x .c/,at the solution.In Fig. 8c, the difference between optimal ease-off associated with c and optimal manufacturable ease-off x .c/(after machinetool setting identification has been completed) is illustrated.The picture remarks the importance of the embedded identification step (see Fig. 4).If the identification process had been performed post-optimization as a secondary step, the discrepancy of the registered error in the range of 5 to 15 micrometers could have significantly affected the resultant contact pattern, due to the potential deviation of the identified parameters from the optimal values obtained through the optimization process.
It is worth also noting that, despite minimizing energy losses, the e L solution is associated with a higher friction coefficient (see the lower values in Fig. 11) and possibly with a higher risk of pitting caused by an increased chance of contacts among surface asperities.
On the contrary, as shown in Fig. 9a, the solution minimizing LTE is spread over the whole tooth surface.In this case, the requirement of minimizing LTE does not map distinctively into a predictable contact pattern shape.This brings about the chance that multiple contact patterns share the same performance in terms of LTE minimization with the associated drawback of multiple local minima.Fig. 9b and c illustrate again the importance of the embedded identification step to minimize the risk of potential deviations when performing machine-tool setting identification as an afterthought.
The significant disparity in frictional losses between two single-objective solutions lies in the different averaged sliding velocities, which are calculated by taking a weighted average with respect to the local contact pressure over the instantaneous contact zones.This is highlighted in Fig. 10.

Pareto front exploration
As a last step, we demonstrate that the proposed framework allows to systematically explore distinct multi-objective optimal solutions by simply selecting different values K Fig. 10 Comparison of the averaged sliding velocities of the contact zones over the meshing cycle between the two single-objective solutions.for the reference point e f.Fig. 12 displays the explored Pareto front, which includes the solutions from the previous single-objective optimizations, labeled P 1 and P 6 , as well as four additional Pareto-optimal solutions, namely the points P 2 ; :::; P 5 .Quantitative values of the solutions are provided in Table 4.
Notably, Pareto solutions P 2 and P 3 , despite belonging to a concave region of the front, have been correctly determined by our procedures.This clearly indicates the advantage of the employed scalarized achievement function defined in Eq. ( 34) rather than the typical weighting method, which would have been unable to obtain such solutions.Moreover, solution P 5 is nearly dominant over P 6 , with a significant reduction of 0.135% in efficiency losses (from 0.890% to 0.755%) and only a negligible increase in RMS.H LTE / (from 0.76 rad to 1.23 rad).
It is worth noting that simultaneous minimization of maximum contact pressure (p max ) and efficiency losses (e L ) would have resulted in a qualitatively similar Pareto front, indeed confirming their conflicting nature.Intuitively, contact pattern localization to minimize efficiency losses leads inevitably to higher maximum contact pressure.However, the same cannot be concluded when considering RMS.H LTE / and p max .Their correlation appears to be more subtle and unpredictable as the minimum of the RMS.H LTE / is not necessarily associated with conjugate mating surfaces, since elastic tooth contact and bending deflections play a significant role.
If preventing gear pitting is a major concern, P 6 could be the best solution, although it presents significant efficiency losses.Solutions P 2 to P 5 provide a good compromise between efficiency, LTE, and contact pressures.However, P 1 may pose an increased risk of gear pitting, but it may be the optimal solution for gears with very good surface finish and advanced heat treatments, since it maximizes efficiency.

Conclusions
In this paper, an automatic methodology to explore Paretooptimal designs of bevel gears has been presented.Concurrent minimization of two conflicting objectives, namely NVH level and frictional losses, was demonstrated.In the first part, a semi-empirical model to estimate frictional power losses under elasto-hydrodynamic lubrication was described, which has been coupled to a state-of-the-art loaded tooth contact analysis tool to obtain accurate predictions of the instantaneous load during meshing.The automatic design relied on a multi-objective optimization approach, where the design variables are represented by few coefficients of a polynomial basis employed to define the tooth flank ease-off topography.Manufacturability constraints were considered by a properly devised mechanism that allows to project the polynomial modifications onto the feasible set of the machine-tool envelopes.The frictional losses were evaluated with the aforementioned model, whereas the NVH level was measured by the loaded transmission error, directly available from the simulation tool.The contact pressures, limited by the material properties, were properly included as nonlinear constraints in the optimization problem from the outset.Application to a test case involving the design of a spiral bevel gearset has revealed that our methodology allows the designer to find multiple Pareto-optimal solutions in a systematic and automatic manner.We believe that the availability of a tool that can present a range of solutions, enabling the designer to choose the most appropriate one based on the unique requirements of the given application, represents a significant advancement in the design of robust, low-noise, and highly efficient gears.Future research will be devoted to obtaining robust solutions which also account for multiple operating loads and corresponding misalignments.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.

Fig. 1
Fig. 1 Schematic representation of two contacting bodies

Fig. 2
Fig. 2 Reference frames of two contacting bevel gears with misalignments .E; P; G; ˛/.Symbol ˙denotes the shaft angle

Fig. 9
Fig.9 Optimal solution associated to the minimum of RMS .H LTE /. a The contact pressure pattern is spread over the whole tooth surface.The requirement of minimizing the LTE does not map distinctively into a predictable contact pattern shape.b Optimal manufacturable ease-off (achieved via new machine-tool settings).c Residual error after identification of the polynomial ease-off

Fig. 11
Fig. 11 Comparison of the lambda ratio over the meshing cycle between the two single-objective solutions

Acknowledgements
Financed by the European Union -NextGener-ationEU (National Sustainable Mobility Center CN00000023, Italian Ministry of University and Research Decree n. 1033 -17/06/2022, Spoke 11 -Innovative Materials & Lightweighting).The opinions expressed are those of the authors only and should not be considered as representative of the European Union or the European Commission's official position.Neither the European Union nor the European Commission can be held responsible for them.Funding Open access funding provided by Università di Pisa within the CRUI-CARE Agreement.

Table 1
Coefficients of the MIL-L-23699 oil

Table 2
Macro-geometry data for the spiral bevel gearset under study

Table 3
Initial machine-tool settings for the conjugate pinion

Table 4
Detailed properties of the found Pareto-optimal solutions