A Multi Domain Modeling Approach for the CFD Simulation of Multi-Stage Gearboxes

: The application of Computer-Aided Engineering (CAE) tools in mechanical design has consistently increased over the last decades. The beneﬁts introduced by virtual models in terms of time and cost reductions are the main drivers for their exploitation in industry as well as for research purposes in academia. In this regard, Computational Fluid Dynamics (CFD) can be exploited to study lubrication and efﬁciency of gears. However, the mesh handling complexities deriving from the boundary motion is still a concern for its application to multi-stage gearboxes. In this work, an innovative multi domain partitioning method for the simulation of a two-stage industrial speed reducer is presented. The implemented solution foresees the combination of two remeshing strategies, namely GRA (Global Remeshing Approach) and GRA MC (GRA with Mesh Clustering), and resulted in a computationally effective performance. The results were compared with experimental data obtained with measurements on the real system, providing a good agreement in the power losses prediction. Considering the complexity of obtaining such results experimentally, the proposed numerical algorithm can offer substantial beneﬁts for an estimation of the transmissions’ efﬁciency in various operating conditions. The numerical model was built in the open-source environment OpenFOAM ® .


Introduction
The reduction of wear in mechanical components as well as the more stringent requirements in terms of energy savings require the proper design of mechanical transmissions (i.e., gearboxes) to ensure a proper lubrication supply. Indeed, the lubricant is a fundamental design factor to guarantee a long functioning time and to reduce the downtime periods of the system. By establishing a thin film between the surfaces in contact, it decreases the frictional (load-dependent) power losses. Lubricant is also involved in heat dissipation. The determination of the correct lubricant quantity is, therefore, a key aspect to achieve eco-friendly solutions. At the same time, due to the mutual interaction between the lubricant and the system's components, additional power losses originate. These losses are called load-independent (or no-load) and can be subdivided in three main contributions: churning, squeezing and windage. The churning phenomenon is typical of dip lubrication conditions, where the oil resistance to the mechanical components' movement causes high pressure and viscous effects. Windage refers to the losses that arise from the interaction of mechanical components and one fluid only (e.g., gears completely immersed in oil). Lastly, the pressure gradients in the engagement region generate axial fluxes that promote the squeezing power losses (also called pocketing).
The first analyses aimed at proposing empirical equations for the estimation of the power losses were performed by [1][2][3], who derived analytical equations for the calculation of the churning losses of a smooth disk in dip lubrication conditions. The first work on gears was done by [4], who proposed an equation that considers the oil viscosity, the tangential velocity of the gear, the tip diameter and the face width. One of the milestones in the study of the power losses in gears is represented by the work of Mauz [5], who, for the first time, considered the mutual interaction between mating gears in the analytical formulations. Despite being one of the most complete models available in the literature, it has some important limitations, among which the fact that it can be used only for spur gears with reduction ratios between 1 and 2. Moreover, the equations are valid for tangential velocities between 10 and 60 m/s only. Specific studies focusing on the gears' design parameters (module, diameter, face width, etc.) and on the lubricant characteristics (immersion depth, viscosity, density, etc.) were done by [6,7]. The experiments were performed varying the rotational speed between 100 rpm and 3000 rpm, the module from 2 mm to 8 mm and the oil viscosity from 20 mm 2 /s to 2000 mm 2 /s. The influence gearbox dimensions and the radial wall distance were also investigated. An increase of the churning losses with the immersion depth and the oil viscosity and a digressive growth with the angular velocity of the gears was observed. Equations based on dimensionless numbers for the calculation of the churning power losses were proposed by [8,9]. The scholars noticed that the module has little influence on the churning losses and that the wheel's diameter becomes an important factor only at high rotational speed. Moreover, they observed that the total losses generated by two mating gears are not equal to the individual contributions of each one of them when measured separately. Specific experimental analyses on a gear pair were performed by [10], who analyzed various operating conditions including the rotation direction of the wheels (clockwise and counter clockwise). While all the above-mentioned studies have focused on spur gears, recent formulations for the calculation of the churning power losses of spiral bevel gears were proposed by [11,12], who concluded that the immersion depth has a strong influence on the churning power losses of this type of gears.
Nowadays, it is possible to study the complex phenomena involved in gears' lubrication exploiting numerical solutions. CFD offers the possibility to model the desired system without the limitations on the operating conditions typical of the empirical models. Two possibilities are available to build digital twins of geared systems: meshless and meshbased approaches. Meshless approaches are based on the discretization of the fluid with discrete elements (contrarily to mesh-based approaches that discretize the spatial domain). The great advantage of these techniques is related to the fact that the issues associated with the topological changes are avoided. Meshless models are capable of predicting the main lubricant fluxes, while the power losses calculation are not particularly accurate [13]. Complex domains involving worm and bevel gears as well as multi-stage gearboxes were studied exploiting meshless methods as the SPH (Smoothed Particle Hydrodynamics) and the MPS (Moving Particle Semi-Implicit) [14][15][16][17]. These works highlight the effectiveness of meshless approaches in the prediction of the lubricant fluxes even for complex kinematics. On the other side, mesh-based methods foresee the discretization of the "empty space" between the moving objects, i.e., the regions where the oil can flow. Simplified numerical models exploiting sectorial symmetry were presented by [18][19][20] for the prediction of the windage power losses of a single gear. The assumption of a 2D model led to an underestimation of the losses as compared with measurement. Dai et al. [21] exploited the MRF (Moving Reference Frame) approach to study the windage power losses of a spur gear, showing good agreement with experimental data, while Cavotta et al. [22] studied the churning power losses of a disk. Following the same approach, Bianchini et al. [23] analyzed the load independent power losses of an epicyclic gearbox, concluding that the trend of dissipated power with rotational speed is less than cubic. Extensive analyses based on a dynamic mesh approach were done by [24][25][26], who validated the numerical models of a FZG gearbox (ISO 14635-1 [27]) with experimental data in several operating conditions, demonstrating the accuracy of the implemented virtual models.
Thanks to recent advancements in computer science, more complex domains with respect to spur gears were approached with numerical tools. Peng et al. [28][29][30] built a CFD model of an hypoid gear pair. The results were compared with experimental observations in terms of fluid flow and power losses. Lu et al. [31,32] analyzed the churning power losses and the temperature profile of an intermediate gearbox of an helicopter. The simulations were massively parallelized to reduce the computational time. Oil jet lubrication was studied by [33,34] in order to find the position of the nozzle that guarantees the optimal lubrication of the teeth. A numerical model for the prediction of churning and windage power losses of a single spiral bevel gear was presented by the authors [35]. The simulations could predict the trend of the losses justifying from ta numerical point of view the experimental observations.
As can be noticed, all the works deal with single-stage gearboxes of various gear types (spur, helical, bevel), but not with multi-stage industrial gearboxes. On the other side, works dealing with complex domains have been found only meshless approaches, in which the issues of mesh deformation are avoided.
In this work, a Finite Volume Method (FVM) model for the study of a two-stage gearbox has been implemented in the open-source environment OpenFOAM ® [36]. The numerical model is based on the exploitation of two remeshing strategies previously developed by the authors (GRA-Global Remeshing Approach and GRA MC -Global Remeshing Approach with Mesh Clustering). These have been combined together in order to be able to handle the rotation of multiple gear pairs. By doing so, a CFD model for the efficient simulation of an industrial gearbox has been developed. Several operating conditions have been investigated, including different rotational speeds and temperatures. The numerical results have been compared with the experimental data. The aim of the proposed approach is to show how, by exploiting an efficient strategy for the remeshing process, it is possible to simulate a multi-stage gearbox efficiently even with a mesh-based approach.

Materials
Multi-stage gearboxes are characterized by a compact design with high transmission ratios and the possibility to combine different gears' types. In the current analysis, an industrial gearbox produced by Comer Industries was considered. It is mainly employed in combined harvesters and self-propelled machines, and in agricultural machinery and garden equipment. Its main geometrical characteristics are reported in Table 1. The analyzed solution consists in a two-stage gearbox, with a total transmission ratio of 4.76. The power enters from the input shaft (on which the first gear is mounted) and is transmitted through the gear train to the output shaft (where the fourth gear is mounted). Each of the three shafts is supported by a pair of ball bearings. The cross-section of the gearbox is illustrated in Figure 1. The gearbox was tested in no-load operating conditions. An electric motor supplies the power that is completely dissipated in terms of heat by the mechanical components. The current testing configuration allows a direct measurement of the load-independent power losses. It is based on the principle of power circulation. The transmitted torque is applied by a load clutch and measured by a load torque meter. The electric engine only supplies the total power loss, which is measured by a loss torque meter between the electric engine and the power circle. In no-load operating conditions, only negligible torque is applied to the load clutch, so that the load-independent power loss can be evaluated directly. A Shell Spirax S4 TXM oil, whose properties are listed in Table 2, is used. This lubricant exhibits excellent characteristics for what regards stability against oxidation and resistance to wear, corrosion, and foaming effects. It is mainly used for hydraulic systems and auxiliary systems of agricultural machines and tractors. Viscosity index -138 Flash point • C 220 Pour point • C −42

Mathematical Description
The numerical model has been implemented in OpenFOAM ® . CFD codes are based on the mass, momentum, and energy conservation equations. These are second-order partial differential equations. Except for few simple cases in which the system of equations can be solved analytically, numerical methods should be introduced. In the current analysis, the thermal effects have not been considered. Therefore, the energy equation was not activated.
The FVM method [37] is adopted. This implies that the computational domain is subdivided into finite volumes. The integration of the conservation equations on each finite volume leads to the following: where ρ is the density, U is the velocity, P is the pressure, = τ is the viscous stress tensor, G represents the external forces, dA is the differential area, and dV is the differential volume. These equations are approximated to linear equations through the Gauss quadrature. Generally, higher interpolation orders and finer mesh size lead to more accurate numerical results. In order to consider the simultaneous presence of more than one fluid, the Volume of Fluid (VOF) method [38] was adopted. This model foresees the introduction of the volumetric fraction, representing the percentage of each fluid in a computational cell: where α is the volumetric fraction (0 for pure air and 1 for pure oil, in this case). The physical properties of the equivalent fluid (ϕ eq ), calculated as a combination of the ones of the two fluids (ϕ 1 and ϕ 2 ), are defined as: A two-phase solver for incompressible and isothermal conditions, and immiscible fluids was used. The properties of the lubricant have been set a priori based on the following relationships [39]: where ρ is the density, ν is the viscosity, and θ is the operating temperature. Each gear is assigned with its rotational speed according to the transmission ratio, while the pressure and the volumetric fraction must be calculated (Neumann boundary condition). The box, on the contrary, is fixed in space and has a zero-velocity boundary condition, while the pressure and the volumetric fraction must be calculated (Neumann boundary condition). The PIMPLE (merged PISO-SIMPLE) algorithm was adopted to solve the transient simulations. The SIMPLE algorithm allows to reach the solution's convergence very fast, but it can be used only for steady-state cases. The effectiveness of this algorithm is limited by the fact that it does not contain temporal information. Conversely, the PISO algorithm is time conservative and can be used in transient simulations. However, the solution's convergence is associated with a reduction of the maximum allowable time step. The PIMPLE combines the advantages of the two schemes: it allows to use big time steps with relaxation (from the SIMPLE) and, at the same time, it assures the time consistency of the solution (from the PISO). Three correctors of the pressure equation were imposed to maintain the stability of the solution. A convergence criterion of 1 × 10 −5 was imposed to all field's variables. The PCG (Preconditioned Conjugate Gradient) solver was used for the pressure. The velocity was solved with the PBiCG (Stabilized preconditioned bi-conjugate gradient). The timestep was set to achieve a maximum Courant number of 1. The time derivative was discretized with the first-order implicit Euler scheme, while the velocity with the second-order linearUpwind scheme and the convection of the volumetric fraction with the vanLeer scheme.

Numerical Approach
In order to simulate the deformation of the mesh arising from the engaging wheels, a dynamic mesh approach must be introduced. Indeed, the gears' rotation requires to update the numerical grid every little angular step so to ensure that low-quality elements do not impact the solution's convergence. The mutual engagement of the wheels cannot be approached with simple modeling techniques as the Moving Reference Frame or the Sliding Mesh approaches. Therefore, an advanced mesh handling strategy has been employed. Generally, two remeshing methods are used to simulate multiple gears: Local Remeshing Approach (LRA) and Global Remeshing Approach (GRA). The LRA is implemented in some commercial software to handle dynamic domains. With this technique, the distorted elements are substituted with new valid ones based on the imposed quality criteria. This approach is very effective and allows to simulate not only spur gears, but also more complex systems such as helical and bevel gears. In OpenFOAM ® , this method is not available, and a different procedure should be followed, namely the GRA (implemented by the authors). The main advantage of the GRA is that the user has a direct control on the elements' size, which can be kept constant throughout the simulation without the risk of very small elements that may impact the time step and, consequently, the computational effort required. On the other side, the GRA is an efficient methodology as far as the mesh generation is performant, e.g., for spur gears, where it is sufficient to mesh only a 2D plane and then extrude it properly to create the full 3D swept mesh. For additional details on the GRA, refer to [40]. The authors have implemented another remeshing strategy, named GRA MC , in order to efficiently simulate also gears' geometries that do not belong to the previous category. The latter approach is based on the recursive use of a set of meshes that covers one engagement. By doing so, with only few numerical grids it is possible to simulate the whole wheels' rotations. Exploiting this algorithm, it was possible to also model helical and bevel gears [41,42].
The system under investigation is a two-stage helical gearbox. A straight application of the GRA MC is not possible. In fact, this method requires to coherently set the time libraries and the rotational speeds based on the transmission ratio. In this case, the transmission ratio is an irrational number. This means that it is not possible to synchronize the two stages of the gearbox so only a limited number of meshes can be used. Indeed, the positions of each of the two gear pairs are not periodic, and therefore, another solution must be considered. In order to contain the computational cost of the remeshing process deriving from a pure GRA implementation, a multi-domain approach that exploits a combination of GRA and GRA MC has been adopted. In particular, the GRA MC was used for the input gear pair (1&2), while the GRA was used for the output gear pair (3&4). By doing so, the whole domain was decomposed in two sub-regions, each one assigned with a different remeshing strategy. Salome [43] and Python [44] scripts were used to create the numerical grids. The Netgen [45] algorithm was used. In a first step, the corners' points are computed. Then, the edges are meshed in segments. Finally, an advancing front surface algorithm discretize the faces. The elements are generated by a Delaunay triangulation, which subdivides the edges, seeds each surface with internal points, and connects them together. The final elements' number of the assembled partitions was about 600k tetrahedral cells. The numerical connection between the two domains was achieved via the introduction of an interface that links the non-conformal grids (Figure 2). The input gear pair was handled with the GRA MC . One complete engagement required the computation of 12 mesh. At the 13th mesh substitution, the first mesh could be reused. On the other side, the output gear pair was managed with the GRA. This means that the grids are created and substituted based on the rotation of the input gear pair. However, since the output gear rotates at a much slower angular velocity, the mesh deformation is lower. Therefore, there was no need to update the mesh of the output domain every time the mesh of the input domain was substituted. In the current analysis, the output gear pair mesh could withstand a deformation equal to 3 mesh substitutions of the input gear pair. Afterwards, a new valid mesh with the wheels in the new positions was provided and the simulation could restart. In Figure 3, the flowchart of the described mesh handling strategy is reported. The results were mapped from mesh to mesh until a quasi-stationary condition of the fluid dynamics forces was reached. In Figure 4, the scheme of the simulation procedure is shown. In order to ensure the solution convergence and stability, a very strict control over the mesh quality parameters has been implemented. The maximum allowable nonorthogonality and skewness were set to 75 and 2, respectively. A Bash [46] script that controls the quality of each generated mesh was coded. If the imposed quality parameters were satisfied, the mesh was kept and saved for the simulation process; otherwise, a new mesh was generated for the same gear position with slightly different meshing parameters until the criteria were satisfied. In Figure 5, the scheme of the meshing process is reported.

Results
The validation of the virtual model was made by means of the comparison of the power losses with experimental measurements. The numerical calculation considers the inertial and viscous contribution given by each gear: where F pi are the inertial forces, F τ i are the viscous forces, r i is the distance of each face center from the specified axis of rotation, and ω is the rotational speed. The individual contributions are then summed, providing the total power dissipation. The measured power losses include the bearing (P LB ) and the seal contributions (P LD ). In order to isolate the gears' losses (P LG ), the P LB and the P LD contributions have been calculated according to the ISO 14179 [47] and then subtracted from the total losses. The power loss for a specific operating condition is taken as a mean value over subsequent engaging periods. The simulation was considered stabilized when less than 5% different was observed between two subsequent averages. The gearbox was tested in a vertical position and filled with 0.8 L (corresponding to half of the output gear filling level). In Figure 6, the experimental apparatus is shown. An electric motor supplies the input power to the big wheel, from which the motion of the system is initialized. The power is completely dissipated in terms of heat by the mechanical components. Three rotational velocities (3000 ÷ 9000 rpm) and three temperatures (40-60-80 • C) for a total of nine operating points were analyzed (3 × 3 full factorial Design of Experiments). In Figure 7, the comparison between the experimental data and the numerical predictions in terms of non-dimensional quantities are reported. Indeed, because of confidential reasons, the results have been divided by the maximum predicted power loss, such that the maximum value is 1. The best agreement was found for 3000 rpm. The results at 6000 rpm agree with experimental data with a maximum error of 14%. For the 9000 rpm cases, the numerical results exhibit an error about 30%. Considering that the prediction of the efficiency's trend can already give important design information, the implemented approach seems to be an effective tool even for complex geometries as multi-stage gearboxes. The effect of temperature on the power losses in the range of investigated velocities can be clearly identified: the higher the temperature, the lower the losses. The other variable (the velocity) being the same, this trend can be explained through the temperaturedensity/viscosity relationships. Indeed, the increase of the temperature is associated with a decrease of the density and viscosity of the lubricant. This translates into a reduction of the power dissipation, but also into poorer lubrication capabilities and lower system reliability. During the experimental tests, it was observed that at high rotational speeds, the temperature effects became not-negligible. Consequently, the physical properties of the lubricant change during operation. This could be a possible explanation for the higher discrepancy between the numerical results and the experimental data at 9000 rpm. In Figure 8, the thermal measurement performed with a thermal camera are reported. A maximum of 95 • C can be noted on the input seal, while the mean temperature on the box is about 83 • C. In Figure 9, the percentage contributions of the three sources of losses (gears, bearings, and seals) are reported. It can be observed that the gears losses represent the majority of the total losses. The bearings' contributions were significantly lower with respect to the gears' contributions. The power dissipation given by the seals gave a non-negligible contribution, especially at low speeds. In Figure 10, the percentage contribution of the four wheels to the total load independent power losses is shown. It can be noticed that for all the analyzed configurations the bigger gear gives the highest contribution. The losses are generally proportional to the rotational speed. This effect is evident comparing the losses of gears 1 and 2. Gear 2 and gear 3 have the same rotational speed, but different pitch diameters, i.e., different tangential velocities: this reflects the lower power dissipation given by the latter. The high losses exhibited by gear 4, which has the same pitch velocity of gear 3, is related to the different loss mechanism (oil churning vs air windage). Indeed, the lubricant density and viscosity are of some order of magnitudes higher than the ones of air. This has a direct impact on the individual loss contributions. This effect could be clearly observed in Figure 11, where the shear of viscous contribution is negligible for gear 4 (oil churning). For the other three gears, the total losses depend on both the inertial and the viscous contribution (windage). Most of the oil is involved in the lubrication of the output gear pair, thus generating power losses due to churning, while the input gear pair is less wetted and the windage effects are predominant ( Figure 12). In Figure 13, it is possible to see the streamlines of the lubricant fluxes inside the gearbox. The majority of the droplets of lubricant escape radially from the wheels.

Discussion
Differently from previous analyses, for the first time a mesh-based approach was successfully applied for the CFD simulation of an industrial multi-stage gearbox. This was possible thanks to the combination of two remeshing strategies (multi-domain approach), each one assigned to one stage of the gearbox. The proposed approach has a general validity and can be applied to study any multi-stage gearbox. Indeed, once a gear pair is selected for the remeshing with the GRA MC , the remaining stages can be meshed with a simple GRA and the domains merged together through the introduction of interfaces. The simulations were parallelized among four cores (INTEL Xeon ® Gold 6154 CPU, 4 Cores, 3 GHz−48 GFLOPS) and required about 1 day to reach the regime condition where the efficiency was evaluated (given that all the meshes have been calculated in advance-1 /2 day of computation-and are immediately available for the interpolation process). The simulation results showed a good agreement with experimental measurements, despite some discrepancies at 9000 rpm. The assumptions that have been made (bearing and seals are not modeled and the box shape is simplified) can be a possible explanation for this mismatch. Moreover, experimentally, it was observed that at higher velocities, the hypothesis of isothermal fluid loses its validity; the inclusion of a thermal solver would probably increase the accuracy of the numerical results. Furthermore, complex physical interactions between the lubricant and the air were neglected. Recent findings by the authors [48] have shown that, especially for fluids having a high energetic content (namely for high rotational speeds), aeration could reduce the losses with respect to the ones calculated for completely immiscible fluids. This could be another explanation for the discrepancies at high speeds.
The analyzed gearbox includes small vanes and complex gears' geometries. Additionally, the complexity of the box shape is higher with respect to the standard test rigs that are used to validate CFD models of gears. Despite the adopted geometrical simplifications, the computational domain is still of another level of difficulty in terms of meshing procedure and topological management. Without an efficient mesh handling strategy, the simulation of such a system would become unacceptably expensive from a computational point of view. This work highlights the maturity of CFD codes as well as their possible application in industrial practice, where such a model can be exploited for an initial estimation of the efficiency trend in various operating conditions.

Conclusions
The introduction of virtual prototyping in mechanical design has been continuously increasing. Indeed, the advancements in hardware and software technology led to the possibility to study specific physical phenomena numerically (structural-fluid dynamicsmultibody-problems). These numerical tools allowed a better understanding of the physical phenomena and to preliminary evaluations of the systems' performance. In this work, a CFD model of an industrial two-stage gearbox was implemented. The decomposition of the domain in two sub-domains allowed to exploit a limited set of mesh for the input gear pair, and to create the mesh of the output gear pair at occurrence, thus limiting the computational costs of the remeshing process. The virtual model allowed to study the efficiency of the gearbox in different operating conditions. The best accuracy was achieved at low (3000 rpm) and medium (6000 rpm) rotational speed, while at 9000 rpm, some discrepancies were noticed. The trend of the power losses at different operating temperatures could be in any case identified. Moreover, some interesting considerations on the wetted surfaces and velocity streamlines could be gathered. Considering the amount of information that such a model can provide in the design phase, the benefits of the introduction of numerical tools in industrial practice can be clearly understood.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, [F.C.], upon reasonable request.