CFD simulations of gearboxes: implementation of a mesh clustering algorithm for efficient simulations of complex system’s architectures

In the last decade, computer-aided engineering (CAE) tools have become a determinant factor in the analysis of engineering problems. In fact, they bring a clear reduction of time in the design phase of a new product thanks to parametrical studies based on virtual prototypes. The application of such tools to gearboxes allowed engineers to study the efficiency and lubrication inside transmissions. However, the difficulties of handling the computational domain are still a concern for complex system configurations. For this reason, the authors maintain that it is fundamental to introduce time efficient algorithms that enable the effective study of any kind of gear, e.g., helical and bevel configurations. In this work, a new mesh handling strategy specifically suited for this kind of studies is presented. The methodology is based on the Global Remeshing Approach with Mesh Clustering (GRAMC) process that drastically reduces the simulation time by minimizing the effort for updating the grids. This procedure was tested on spur, helical, and bevel gears, thus demonstrating the flexibility of the approach. The comparison with experimentally measured power losses highlighted the good accuracy of the strategy. The algorithm was implemented in the opensource software OpenFOAM®.


Introduction
Thanks to continuous technological developments in computer science, simulation software packages have gained their consideration not only among academic researchers, but also in industries. The advantages introduced by virtual prototyping, namely the reduction of time and costs with respect to real prototypes, have encouraged engineers to exploit these new opportunities. These technologies fall into the "CAE" category. For structural analysis the finite element analysis (FEA) is used, while the study of fluids is achieved through the computational fluid dynamics (CFD). The fluid-structure interaction (FSI) represents the coupling of the above-mentioned tools. While at the beginning of the twentyfirst century opensource software were not widespread and used mainly by the developers, nowadays they are likewise widespread in academia and industry. They show the same robustness and maturity of commercial software. One of the most appreciated aspects of opensource software is related to the possibility to use parallel calculus without additional licenses for it. Parallel computing can dramatically speed up the calculation time of the simulations, especially in the case of huge computational domains involving several physical phenomena at the same time. This aspect is fundamental especially in industrial applications for which the time to market is a key factor for the success. On the other hand, the main drawbacks of free simulation tools are represented by a steep learning curve, meaning the user must invest some time to master the code.
The application to gearboxes of CFD codes represents an important aspect that opens the way to the design of more and more efficient solutions or to the optimization of the existing ones. Indeed, important considerations on power losses and lubricant fluxes inside the transmission can be obtained numerically. Conversely, the experimental testing of these systems is not trivial since it requires skilled technician capable of designing such experiments and advanced equipment for measuring the resistant torque (torque sensors) and monitor the lubrication capabilities (particle image velocimetry-PIV measurements). For these reasons, the authors consider that specific studies exploiting CFD can be beneficial in the sector of tribology on one hand considering the more and more stringent requirements of energy efficient solutions and on the other hand for a better understanding of the physics involved. In this work, a timeefficient mesh handling procedure is used on configurations for which no previous analyses have been made. The opensource code OpenFOAM® (OpenFOAM n.d.) was used. In particular, while works dealing with systems in which extrusion algorithms can be adopted (e.g., spur gears, epicyclic spur gearboxes, cycloidal drives) have been already published (Gorla et al. 2013;Liu et al. 2018;Concli et al. 2019), particular gears' geometries such as helical and bevel ones, present in many of real applications, were poorly touched by numerical studies. This is mainly due to the difficulties of handling the remeshing of such domains ensuring a sufficient quality of the grid during the simulation and the convergence of the solution. Trying to fill this lack, the authors have therefore implemented a remeshing strategy that can be applied to any system configuration, thus overcoming the current limitations of the available approaches.
First, a comprehensive overview of the dynamic mesh methodologies that can treat interacting bodies as in gearboxes is presented. Advantages and disadvantages of each approach are illustrated. The second part of this manuscript is devoted to the description of the new mesh handling strategy that allowed a time efficient simulation of real application gearboxes. Good results in terms of power losses are obtained by comparing the numerical predictions with experimental data from other authors. The approach shows a very high computational efficiency.

Analysis of mesh moving strategies and their application to gearboxes
In order to discretize the gears inside a transmission, it is fundamental to rely on a robust mesh generator. Native applications of OpenFOAM® are blockMesh and snappyHexMesh. The former is a utility that builds full structured mesh providing the user with extremely high control on every aspect of the mesh generation process; however, it can be applied to rather simple domains and is not suitable for gearboxes. The latter is a utility that builds hexahedral dominant grids and can be applied also to complex domains. Nevertheless, it is quite a time-and memory-consuming tool that does not represent the optimal solution in these situations. Another alternative is represented by cfMesh, an external grid generator that in its base version can be installed freely as a library that communicates with OpenFOAM®. It exploits shared-memory parallelization (SMP) and is generally faster than snappyHexMesh, which is based on distributed-memory parallelization (DMP). However, the meshing process is not so easy to control especially near the gears where refinements are applied. Third party software (commercial or opensource) can be used as well to create the numerical mesh that can be successively converted to the OpenFOAM® format. In this regard, Salome (SALOME n.d.) is a powerful opensource preprocessor that can build polyhedral grids giving the user a good control over the meshing parameters. This software is highly appreciated because it comes with a developed Graphical User Interface (GUI) that makes the meshing process simpler.
Since the correct setup of the CFD simulation of gearboxes is already demanding from a modeling point of view (multiphase with dynamic mesh and turbulence must be considered simultaneously), it is fundamental to choose the meshing strategy that mostly reduces the complexity of the virtual model and, possibly, that requires the lowest interaction with the user, i.e., automatized process. Before explaining the remeshing strategy used in this work, a review of the possibilities currently available is provided, and the main advantages and limitations are discussed.

Overset mesh approach
The overset mesh approach (also called overlapping grids method or chimera framework) was implemented by (Benek et al. 1944;Benek et al. 1985). The strategy is based on the creation of a fixed background mesh that contains the meshes of the moving features, which are not conformal among them, but overlap themselves. By doing so, the grid is rigidly moved during the calculations, thus the overall mesh is not deformed. During the simulation the computational cells can be identified by the solver as holes (that represent the patches and do not require any solution since no fluid is passing through), calculated (that represent the regions of the domain where no overlapping occurs), and interpolated (that are the cells in the overset regions in which the solution of the equations are computed by interpolating the results between the background mesh and the overset meshes). This approach is very effective since there is no need to create a single mesh that contains all the features, but all the grids corresponding to the moving domains can be created separately and then "overlapped" to the background mesh. This implies that a high-quality global mesh can be obtained due to the decomposition in smaller domains. However, the main drawback of this approach is related to the computational effort in gearboxes simulation derived from the continuous interpolation process that takes place every timestep. Proprietary codes such as STAR-CCM+ (STAR-CCM+ n.d.), ANSYS FLUENT (ANSYS FLUENT n.d.), and FLOW-3D (FLOW-3D n.d.) are implemented with this technique which comes from years of usage and validation. On the contrary, the overlapping grids method was implemented in OpenFOAM® only in some distributions, meaning that further developments are needed to reach a maturity comparable to commercial software. Application of this technique to geared transmissions can be found in (Arisawa et al. 2017;Arisawa et al. 2009;Saegusa and Kawai 2014;Renjith et al. 2015;Deshpande et al. 2018;Cho et al. 2019;Patel and Anto 2019;Hu et al. 2019;Renjith et al. 2017). The results of the simulations allowed to calculate the efficiency of the transmissions and to verify the lubricant distribution in several operating conditions. Studies on gearboxes exploiting the overset mesh technique have been done only in commercial codes.

Immersed boundary method
The immersed boundary method (IBM) was initially used by (Peskin 1972) to study the flow around a heart valve. This method does not require the grid to follow the geometrical features; a source term in the momentum equation is added allowing the adaptation of the flow to the considered solid boundaries rather than an actual mesh motion. The complete description on the implementation of this method in OpenFOAM® can be found in (Jasak and Tukovic n.d.). The main drawbacks of this method are related to the low resolution of the computational cells near the boundary layer of the walls and to the complex implementation of the boundary conditions at the immersed boundary. This approach was used by (Burberi et al. 2016) to study two mating spur gears. They concluded that the IBM lead to an overestimation of the losses, but can be a valid approach to initialize the flow fields. This can be explained through the non-alignment of the edges to the wall. Therefore, this method is usually considered not appropriate for the simulation of gearboxes.

Local remeshing Approach
One of the most advanced techniques to handle mesh deformations is the local remeshing approach (LRA).
This method requires to have tetrahedrons in the remeshing region and exploits the spring-based smoothing algorithm. Being the elements' edges considered as linked springs, the movement of the grid in the dynamic mesh region causes a force on each node that connects the springs that can be formulated exploiting the Hook's law: where k ij is the spring constant between nodes i and j, Δ x j * and Δ x * i represent the displacements of the node j and i, and ni stands for the nodes linked to node i. The term k ij is calculated as: The solution of the equilibrium of forces is achieved iteratively from: This equation can be solved by updating the positions of the nodes (which are known from the displacements) at the next timestep: Based on the quality of the cells in the remeshing domain, the too distorted elements are substituted with new valid cells and the fields are interpolated. By doing so, the movement of the boundaries is accomplished. The application of this method to gearboxes led to significant advancements in the analysis of such systems, allowing the study of spur gears (Gorla et al. 2013), spiral bevel gears (Hu et al. 2018), and planetary gearboxes (Liu et al. 2018) in commercial software.
OpenFOAM® handles small mesh deformation applying a Laplacian smoothing. However, when the deformation is excessive, the mesh should be (locally or globally) updated/substituted. Some commercial software are already implemented with the local remeshing approach, which has proven to be effective, but the mesh updating procedure is left to the software and not to the user. In this regard, a global mesh substitution could introduce significant benefits, which are discussed in the next paragraph.

Global remeshing approach
In order to model intersecting solid bodies as mating gears, a full automatic procedure (Concli et al. 2016) has been implemented in OpenFOAM®: the global remeshing approach (GRA). The algorithm exploits Salome and Python (Python n.d.) to automate the geometry and meshing procedure, and a script in Bash (Bash n.d.) that controls all the dictionaries and the solution in Open-FOAM®. The Laplace smoothing equation and the pseudo solid equation are used to describe the motion (z) of the boundaries: where γ is the diffusivity, x is the grid's nodes position and t is the current timestep. The rotation of the wheels is accomplished through the creation of several grids corresponding to subsequent gears' positions. The results are interpolated from mesh to mesh until the regime condition is reached. The interpolation procedure follows an inverse distance weighting 2 nd order method.
This modeling technique has the great advantage of giving the user a perfect control on the elements' size, which are kept constant throughout the simulation, thus allowing a stable timestep while maintaining the Courant number lower than 1 and avoiding convergence issues: where δt is the timestep, |u| is velocity magnitude in a cell and δx represents the length of a computational cell in the flow direction. The computational advantage given by the remeshing process is enhanced by the mesh generation itself. In fact, for spur gears it is sufficient to mesh only a surface representing the 2D representation of the box and the wheels' partition. This plane is then appropriately extruded creating the whole 3D domain. The element type is therefore a triangular prism. This kind of cells exhibits outstanding quality parameters, making it possible to reach a maximum nonorthogonality below 35 and a maximum skewness below 1. When the mesh deforms as a consequence of the gears' rotation, these quality indicators worsen. The allowed maximum non-orthogonality and skewness after the deformation is imposed to be 70 and 2 respectively. If one of these thresholds are exceeded, a new mesh is provided as input and the results are mapped onto the new valid grid. The non-orthogonality (angle between the line connecting two cell centers and the normal of their common face) and the skewness (distance between the intersection of the line connecting two cell centers with their common face and the center of that face) are fundamental parameters that indicate the quality of a mesh. Especially in cases as this one (where the domain deforms as a consequence of the boundaries rotation), it is determinant for the success of the simulation to build high-quality grids. Running simulations with poor quality parameters will make it difficult to reach stability and convergence, introducing diffusion errors in the simulation. A maximum non-orthogonality of 75 and a maximum skewness of 4 are considered as threshold to run a simulation. Since the simulation of gearboxes involves also dynamic domains, these limits have been further lowered to 70 and to 2 to try to reach maximum stability. Works dealing with this methodology can be found in (Mastrone et al. 2020;Concli and Gorla 2017;Lucchini et al. 2015;Montenegro et al. 2014).
The main limitation of the GRA is represented by the fact that the procedure is effective as far as an extrusion algorithm can be exploited (as in spur gears). Geometries that do not fall into this category cannot benefit from the fast remeshing process. For this reason, the authors have implemented a strategy that can be used to study also more complex geometries with a dramatic reduction of the computational effort. The detailed explanation of this procedure is explained in the next paragraph.

Global remeshing approach with mesh clustering
In order to be able to simulate domains which cannot be meshed with extrusions' algorithms (that were shown by the authors in previous studies), a new strategy that exploits the repetitiveness of the wheels' profiles after one engagement is implemented. In fact, in the GRA explained in paragraph 2.4, the mesh updating algorithm is built on a grid-quality basis rather than on a geometrical basis, meaning that each grid is exploited until its validity limit without specific constraints on the wheels position. This is related to the fact that the mesh generation requires just few seconds on a good hardware. Conversely, domains involving complex gears' geometry must be meshed with tetrahedrons: in fact, these elements adapt well to the geometrical features almost automatically. The generation time of this type of grids is consistently higher with respect to prismatic meshes, forcing the user to find a different strategy in order to contain the computational costs.
In this regard a new algorithm called global remeshing approach with mesh clustering (GRA MC ) was implemented. As suggested by the name, this strategy is based on the creation of a set of numerical meshes that covers one engagement. After the first engagement occurred, the wheels find themselves in the exact same position of the first mesh. This means that the grids corresponding to the position of the gears during the first engagement can be recursively reused to describe the whole rotation of the gears.
Despite being a rather simple consideration, its implementation is not trivial because the set of mesh necessary to cover one engagement is dependent from case to case, and must be calculated following a specific procedure based not only on the mesh quality indicators, but also on the design parameters of the gearbox (e.g., the transmission ratio). Moreover, all the libraries (packages of code) that control the time must be updated automatically at each mesh substitution.
The first step is the creation of a high-quality mesh. Indeed, OpenFOAM® is very sensitive to the mesh quality: a value of maximum non-orthogonality of 70 is considered a limit. In the authors experience, the value of 70 should not be overtaken in this kind of simulations. Values between 70 and 75 could also be accepted, if some adjustments in the numerical schemes and the solution's algorithms are implemented. Once the first mesh is computed, it is necessary to understand the rotation that can be imposed to the wheels before the mesh quality decreases unacceptably, as previously mentioned in paragraph 2.4. To accomplish this objective, a very useful utility of OpenFOAM® is exploited. This utility, called moveDynamicMesh, allows the user to move the boundaries according to the specified motion in the dictionary that controls the boundaries' movement without solving the flow equations, but just accounting for the topological changes of the domain. The advantages of this utility are related not only to the possibility to verify if the boundaries move correctly, but also to quickly check the mesh quality parameters after each iteration. For this purpose, a generic rotational speed is imposed to the gears. When the mesh quality exceeds the imposed limits, the simulation is stopped and the last valid timestep is taken as a reference. At this point, the angle between two successive teeth (α) is calculated for each gear and the time needed for completing 1 engagement (T eng ) can be obtained: The following step is to identify the updating time of each mesh (t up ) that allows to reach T eng with a set of predefined meshes. This is done by detecting a suitable time bounded in the range obtained with the moveDyna-micMesh utility. This should be an integer divisor of T eng . The number of meshes needed to cover 1 engagement (N mesh ) is then calculated as: The set of grids is implemented in Salome by imposing successive rotations (θ) to the gear and to the pinion according to the transmission ratio: At this point, the set of mesh is ready and can be used to describe the whole rotation of the gears. The governing equations are solved subsequently for each grid; the results are saved and mapped onto the following mesh that is immediately provided as input. A consistent mapping process is used, thus ensuring that all the fields are patched between conformal domains. Furthermore, thanks to the very good control over the mesh generation parameters, the fields' interpolation occurs between very similar grids in terms of elements' size and, therefore, the errors are significantly reduced. Having this control over the mesh creation, mass conservation issues have been avoided. When the regime condition is reached, the results are postprocessed. The regime condition was monitored throughout the simulation by observing the variation of a physical quantity, namely the power losses. The regime condition was considered reached when the variation of the power losses did not differ from the mean value among a sufficiently long period by more the 5%. The average value was then taken as reference to evaluate the power loss for a specific operating condition. While the updating of the time libraries is performed at the end of each t up , the choice of which mesh to use at every step is managed by an if-else statement inside the procedure loop itself. The general workflow of the solution algorithm, which is completely automated in a Bash script, is summarized in Fig. 1.

Analyzed systems
Three transmissions were implemented with the GRA MC : spur, helical, and spiral bevel gearboxes. The spur and helical gearbox are FZG type with dimensions 210 × 175 × 56 mm 3 , while the gearbox used for the spiral bevel gears 300 × 300 × 300 mm 3 . Being able to model these systems is fundamental for a deeper understanding of lubrication in many real applications. With respect to spur gears, helical gears are characterized by a gradual engagement; therefore, the teeth do not undergo an impact load as in spur gears but a smoother loading. This has significative implications on the reduction of noise and vibrations, and the consequent increase in lifetime. On the other side, spiral bevel gearboxes are employed to transmit power between non-parallel shafts. Also in this case, the curved profile of the teeth allows a gradual engagement during operation, thus making them widespread in high load applications, e.g. helicopters, trucks and vehicles in general.
Considering the amount of different applications in which these mechanisms are involved, it is fundamental to dispose of accurate modeling techniques capable of correctly predicting the efficiency and the lubrication of these transmissions. The detailed explanation of how to apply the GRA MC to these specific cases and its enormous advantages in the reduction of the simulation time are reported in the next paragraphs.

Spur gears
Even if it is already possible to simulate spur gears with GRA very efficiently, the GRA MC was applied in order to show the potential further reduction of the computational effort. The gear pair characteristics (taken from (Höhn et al. 2011) for validation) are reported in Table 1.  Exploiting the geometrical characteristics of gearbox, an extruded mesh consisting in triangular prisms was generated. A local refinement of 1 mm was imposed on the gears' edges, while a global size of 5 mm for the rest of the domain was used. The computational domain of the wheels is shown in Fig. 2.
Only half of the domain was modeled by taking advantage of the symmetry of the system. The quality parameters of the set of the numerical grids are listed in Table  2 including the initial and the deformed mesh.
The final set of mesh that respect these indicators is made of 15 meshes, which are recursively reused at each t eng . At the 16 th mesh substitution, the algorithm foresees to reuse the 1 st mesh. In Fig. 3, different wheels positions within 1 engagement are reported.

Helical gears
The modeling of helical gears is based on the experimental tests performed at FZG (Höhn et al. 2011), which were used for validation. An electric motor enables the rotation of the helical wheels whose geometrical characteristics are reported in Table 3.
The GRA MC was adopted to mesh the complex domain involving helical gears. The first step is the creation of a high-quality tetrahedral mesh. A local size of 1 mm was used to refine the mesh on the wheels' edges, while a global size of 5 mm for the rest of the domain was used. A growth rate defining the difference between two adjacent elements of 20% was imposed. The number of loops to optimize the surface and volume mesh was set to 5. The computational domain of the wheels is shown in Fig. 4.
The quality parameters of the set of the numerical grids are listed in Table 4 including the initial and the deformed mesh.
The final set of mesh that respect these indicators is made of 10 meshes, which are recursively reused at each t eng . In Fig. 5 different wheels positions within 1 engagement are reported.

Spiral bevel gears
A spiral bevel gearbox based on the geometry by (Hu et al. 2018) is modeled in OpenFOAM® exploiting the GRA MC . The characteristics of the gear pair are summarized in Table 5.
The computational domain of the wheels is shown in Fig. 6.
The same procedure used to mesh the helical gears was followed. The quality parameters of the set of the numerical grids are listed in Table 6 including the initial and the deformed mesh.
The set of mesh for the spiral bevel gears is made of 10 meshes, which are recursively reused at each t eng .
In Fig. 7, different wheels positions within 1 engagement are reported.

Computational enhancement
The set of mesh for each gearbox was computed on a single core 1.7 GHz CPU (6.8 GFLOPS) machine. A comparison between the GRA MC and the GRA is provided in terms of computational effort to generate the numerical grids and to run the full simulations. The meshing time normalized to one rotation of the bigger gear for each of the three systems is reported in Table 7.
For sake of clarity, these considerations are done considering a serial mesh generation process. This means that additional time reduction can be gained by parallelizing the process. It should be noticed that the adoption of tetrahedral grids (helical and spiral bevel gears) implies a remeshing time of a higher order of magnitude (for comparable gears) with respect to the extruded approach (spur gears). In any case, the proportion of the computational gain is similar for the three systems (since they have almost the same number of teeth): in fact, in all cases, the effort of the meshing process is reduced by about 95% by exploiting the GRA MC . This is related to the fact that this approach foresees the generation of a limited number of meshes with the wheels in predefined locations (according to what explained in paragraph 3). All the dictionaries that control the time are set such that the gears pass through the control positions defined by the imposed rotations. This has positive impact on the simulation itself which is drastically reduced. In Fig.  8, the relative time gain vs. the number of gears rotation is shown for gears with different number of teeth.
Generally, the higher the rotational velocity, the longer it will take to reach the regime condition because of chaotic oil splashing around. The aim of the plot is to highlight that, with the increasing number of gear's rotation, the time due to the mesh handling decreases continuously, thus demonstrating this procedure being an effective solution for the reduction of the computational effort in the simulation of gearboxes in various operating conditions.

Numerical approach
In the following paragraphs, the description of the numerical models used in the simulations is provided. These models describe the mathematical formulations implemented in CFD codes to solve the conservation laws. However, these usually refer to monophasic conditions, which are not realistic in gearboxes due to the simultaneous presence of more fluids. Therefore, the mass and momentum conservation equations are not enough to describe the involved physics even in the complete filling condition. As it will be explained, phase changes (cavitation) and interaction with air (dip lubrication) occur continuously, thus needing additional equations to be considered.

Mathematical description
CFD software are based on the solution of three governing equations: mass, momentum, and energy equation. In this study, the problem was modeled as isothermal. Therefore, the energy equation was not activated. In this way, the solution is limited to the mass and momentum equations which can be written as: These equations are valid only in simulations involving one phase. This condition is not realistic in gearboxes.  In fact, phenomena like cavitation and aeration, which imply continuous phases' transitions, are not included in the previous equations. Therefore, it is necessary to introduce a balance equation that considers the presence of two or more phases. This quantity is represented by the volumetric fraction, which represents the percentage of one fluid in every discretized cell. The VOF (Volume of Fluid) model (Hirt and Nichols 1981) is used. The equation of the volumetric fraction can be expressed as follows: The MULES (Multidimensional Universal Limiter with Explicit Solution) (Rusche 2002) correction was added in the solver algorithm in order to obtain a more stable and bounded solution of the volumetric fraction. This is achieved by adding a dummy velocity field (u c ) in the conservation equation of the volumetric fraction: An additional source term (S g ) must be added to the equation to account for a possible phase change: To calculate the source term, a phase change model must be introduced. The most used in CFD are those by Kunz (Kunz et al. n.d.), Merkl (Merkle et al. n.d.), and Saurer (Saurer 2000). In this work, the Kunz model was used to model one operating condition (to account for cavitation, as it will be explained in detail later). The great advantage of this formulation is related to the fact the source term is independent from the pressure. The vaporization (ṁ v ) is modeled as proportional to the liquid fraction (α) and to the quantity of pressure under the saturation pressure (p sat ). The condensation (ṁ c ) is modeled analogously: With the above-mentioned equations, it is possible to model all the operating conditions object of this work.

Power losses calculation
The mathematical formulation that enables the calculation of the power losses was implemented as a utility that reads the fields at each timestep and calculate inertial ( F p i ) and viscous ( F τ i ) contribution of the losses. The initial transitory is characterized by a fluctuation of the forces which does not allow a correct evaluation of the efficiency of the gearbox. Therefore, it is necessary to observe the variation of these quantities in time. Once these do not oscillate too much and stabilize around a mean value within a certain percentage limit, the solution reached the regime condition and the fields could be postprocessed to compute the losses, which are calculated as:

Numerical schemes
The PIMPLE (merged PISO-SIMPLE) algorithm was used in all simulations. This algorithm allows a better control in transient simulations. In fact, it is possible to tune the correctors of the conservation equations within one timestep in order to reach the best compromise between computational effort and stability of the simulation. A convergence criterion of 1e-6 was imposed to all field's variables. The generalized geometric-algebraic multi-grid (GAMG) solver was used for the pressure, since, in the authors experience, it is much faster than other solvers (like the PCG-preconditioned conjugate gradient) in this kind of simulations. The velocity was solved with a smooth solver (Gauss-Seidel). The Courant number was limited to 1 to ensure the stability of the simulations. The first-order implicit Euler scheme was used for time-stepping. A total variation diminishing (TVD) scheme using the vanLeer limiter was adopted for the convection of the volumetric fraction. The convective fluxes in the turbulence equations were discretized using second order schemes.

Simulations' parameters
In Table 8, the investigated operating conditions are summarized. Three solvers were used: interDyMFoam for the multiphase condition; pimpleDyMFoam for the completely filled pressurized condition; interPhase-ChangeDyMFoam for the completely filled nonpressurized conditions (in order to account for the phase change). The properties of the lubricants for the two configurations are reported in Table 9.

Results
Three operating conditions were analyzed numerically: partial immersion (interDyMFoam), complete immersion in oil pressurized at 6 bar (pimpleDyMFoam), and complete immersion in oil non-pressurized (interPhase-ChangeDyMFoam). These conditions were simulated with the models mentioned previously. The validation of the simulations is done through the comparison with experimental data. Despite being different gear configurations, the experimental equipment to measure the noload losses was similar for the different systems. An electric motor supplies the power (that is completely dissipated) to the test gearbox. A sensor mounted on the shaft measures the resistant torque. For the spur and helical gearbox (FZG), an intermediate gearbox is also present. The no-load losses of the bearings and of other auxiliaries were subtracted from the total losses in order to isolate those given by the gears. In Fig. 9, a simplified scheme of the experimental apparatus is reported.

Spur gears
For the case of complete oil immersion non-pressurized, the cavitation model was introduced in order to better explain from a physical point of view what observed experimentally by (Höhn et al. 2011). The pressurization of the gearbox has been applied to ensure the complete wetting of the wheels also during operation. However, being the lubricant considered as incompressible, the only way the gears can result non-completely wet by the oil is by admitting a phase change from liquid to vapor (cavitation). Indeed, without cavitation, the pressurization should not affect the resistant torque (being the pressure hydrostatic), as instead it was observed experimentally. On the other side, cavitation makes it possible to account for the vaporization pressure during the wheels' rotation which in turn allows to account for a phase change of the lubricant from the liquid state to the vapor state. Consequently, a compensation of the negative pressures on the teeth flanks takes place: in fact, the pressure could increase/decrease its value without constraints if cavitation is not considered. The negative pressure would contribute to the total losses with a suction effect. Instead, considering cavitation, the pressure cannot reach values lower than the vaporization pressure, thus reducing the resistant torque and justifying what observed in experimental tests. Therefore, the local phase changes are responsible for the difference of the losses between the pressurized and non-pressurized condition. Figure 9 shows the comparison between experimental and numerical values in terms of non-dimensional power losses vs. tangential velocities for the three operating conditions. The dotted lines represent a thirdorder polynomial interpolation of the experimental data with an error bar of 10%.
As it can be noticed, the power losses increase dramatically in the complete filling pressurized condition, while in the dip lubrication they are considerably lower. The case in which cavitation occurs shows slightly higher power losses with respect to the partial immersion condition. The solvers showed good accuracy for all the three operating conditions, with almost all points in the 10% error range, suggesting numerical simulations to be an effective tool for the investigation of different scenarios in the design and optimization of new gearboxes. Moreover, thanks to the reduced computational time of the proposed mesh handling technique, it is possible to gather a large amount of information for different configurations in a reasonable amount of time overcoming the current limitations in terms of computational effort. Figure 10 shows the evolution of oil fraction for the dip lubrication case. The oil tends to reach the top of the gearbox creating a cavity below the mating region. At high rotational speeds the oil flow reaches the walls of the box and then falls down causing the entrapment of air in the lubricant. In Fig. 11, it is possible to appreciate the axial gradient in the meshing region confirming the presence of squeezing/pocketing effects. These are related to the pressure increase in the meshing region. When the teeth leave the contact point, the pressure decreases, thus originating a suction effect which is highlighted by the axial velocity streamlines.

Helical gears
The same operating conditions with respect to spur gears were analyzed for the helical gearbox. Figure 12 shows the comparison between experimental and numerical values in terms of non-dimensional power losses vs tangential velocities for the three operating conditions. The dotted lines represent a 3rd order polynomial interpolation of the experimental data with an error bar of 10%.
Similar considerations with respect to spur gears can be made for the analyzed helical gearbox regarding the ratio of the power losses among the three different situations. However, it can be noted that helical gears exhibit lower power losses with respect to spur gears. This is related to the design characteristics of helical gears that allow a better drainage of the fluid thanks to the helix angle. This would suggest that helical gearboxes have higher efficiency with respect to spur gears. Nevertheless, they are more expensive and more difficult to produce than spur gears. Moreover, they are subjected to radial as well as thrust forces, and so they need bearings than can withstand both forces. This points out the importance for engineers to find the best compromise among all these aspects and to choose the best option for the specific application of interest. Figure 13 reports the comparison between pressurized and non-pressurized condition in terms of pressure distribution on the gear. The characteristic peaks associated with the windage phenomenon on both the front and the posterior flanks in the pressurized condition can be noticed, while in the non-pressurized simulation that posterior flank region is cavitating (being the pressure limited to the vaporization pressure) and therefore no peaks are present. This is confirmed by the analysis of the volumetric fraction that reports values lower than one on the teeth flanks in the cavitation case (Fig. 14).

Spiral bevel gears
The spiral bevel gearbox was simulated for dip lubrication conditions at different rotational speeds. In Fig. 15, the results of the simulations are compared to experimental data. The dotted lines represent a third-order polynomial interpolation of the experimental data with an error bar of 10%.
This system allowed to study the influence of the oil fill level, which has a huge impact on the power losses. As it can be intuited, the high number of variables that can influence the efficiency of geared transmissions put some limitations to experimental campaigns in terms of time and costs. On the other hand, numerical tools allow engineers to estimate the power losses of gearboxes by implementing parametric models which allows a rapid examination of different operating conditions reducing the costs associated with physical prototypes.
In Fig. 16, the velocity vectors are illustrated. A recirculation flow around the wheels can be noticed (Fig.  16a). The flow is sucked in the mating region in the upper part of the domain (Fig. 16b) and is ejected in the lower part (Fig. 16c).
Moreover, additional flows are expelled with a 45°a ngle with respect to the rotational axis of each gear. The squeezing effects are highlighted by the oil fraction field near the meshing region (Figs. 17 and 18).

Discussion
Differently from previous analyses, complex domains involving helical and spiral bevel gearboxes were simulated in an opensource environment with a general procedure that drastically reduce the computational effort of the simulations. The algorithm is based on the clustering of a set of grids that covers one engagement. This set is recursively reused at each engagement, thus limiting the complexity of the remeshing process. This algorithm was introduced to overcome the high computational effort required for such simulations, which was substantially reduced due to the efficient management of the meshes. A specific procedure based on mesh quality and design parameters was implemented to determine the optimum number of meshes that allows to complete the gears' rotation. While in the past only spur gears have been simulated with the GRA, a modified version of it, the GRA MC , has been implemented and applied to gearboxes whose geometry has been poorly touched numerically and only in commercial software. The general purpose of this approach makes it ideal in the design and optimization of gearboxes, helping engineers to investigate the influence of the parameters as lubricant level and rotational speed on the efficiency of the transmission. In fact, very important information on the physical phenomena that occur in gearboxes can be obtained numerically. Moreover, information on the transmissions' efficiency in various operating conditions can be gathered fast: indeed, once the case is setup, it sufficient to coherently change the input parameters in the scripts that manage the simulation process to investigate the desired  operating conditions. The methodology can be applied with every solver (i.e., operating condition). The procedure is fully automatic and is implemented in opensource software only. The approach is completely general; it can be implemented in each programming language and coupled transversally with any CFD solver. The specific choice of implementing this algorithm within Open-FOAM® is related to its simplicity: OpenFOAM® is opensource and already written in C++ which promotes a complete integration with the actual code. This paper is aimed to show the computational advantages of the new implemented method without which such simulations would take days/weeks. All the examples consider a single gear stage only and simple housing geometries. Industrial gearboxes have often more than one stage. Moreover, in the presence of high reduction ratios, the slow wheel has a speed of order of magnitude lower than the input pinion. This implies that, in order to reach the physical regime condition, many revolutions should be simulated. In such context, even a moderate time saving at each timestep could lead to tremendous time gains on the global simulation. The implemented algorithm fit in this context to reduce the computational effort of CFD gearboxes simulations.

Conclusions
The application of simulation tools in mechanical design has increased significantly over the last years. This is related to the technological advancements in computer science that have characterized the last decade. The analysis of lubrication and efficiency of gearboxes has benefitted from these improvements. In particular, the application of CFD allowed engineers to reduce time and costs associated to design new transmissions. However, the management of dynamic domains involving multiphase conditions represents still a concern for the complexity of the model's implementation requirements. This is supported by the fact that most of the studies in literature are related to spur gears, which can benefit from faster meshing processes related to extrusion algorithms. For this reason, the authors have implemented an efficient remeshing process, named GRA MC , for the simulation of complex geometries. The algorithm was validated with spur, helical, and spiral bevel gearboxes demonstrating its feasibility for different gears' type and its accuracy in the power losses trend prediction (maximum 20% error with respect to experimental data). In this framework, this development is even more significant for real industrial multi-stage gearboxes, where the current available finite volume methods would suffer from severe computational costs. With the new implemented approach, even the simulation of more complex gearboxes can be approached bringing clear computational enhancements with respect to standard remeshing approach.
Considering that all cases have been simulated by parallelizing the computations among only 4 cores (Intel(R) Xeon(R) Gold 6154 CPU @ 3.00 GHz), it is clearly noticeable that, on a more performant hardware, with this modeling technique it would be possible to approach complex multi-stage gearboxes in a comparable computational time with respect to the analyzed cases. While the aim of this study was to show the capabilities of the new approach and its transversal applicability to many geometries and lubrication conditions, from the testing simulations some interesting considerations emerged. This kind of data could be very interesting to understand the physics behind experimental observations. The simulation process is fully automatized in a Bash script and requires as input from the user only the set of mesh. The opensource software OpenFOAM® was used because of its high flexibility and customization capabilities.