A Detailed Investigation of Gear Body-Induced Tooth Deﬂections and Development of an Improved Analytical Solution

: Many researchers have developed analytical methods to evaluate gear meshing sti ﬀ ness. Some of these methods ignored the e ﬀ ect of the gear body while others used a simpliﬁed model to consider its e ﬀ ect. Until now, a detailed investigation of gear body-induced tooth deﬂections has been rare, especially for the double-tooth-pair meshing period. In this study, we present a detailed investigation of gear body-induced tooth deﬂections. To be speciﬁc, we will discuss how to accurately evaluate gear body-induced tooth deﬂections using the ﬁnite element analysis, and what are the e ﬀ ect of parameters such as loading and gear parameters on gear body-induced tooth deﬂections. Then, an improved solution is developed for evaluating the body-induced tooth deﬂection. In the single-tooth-pair meshing period, the improved formula is developed based on a popular formula proposed by Sainsot and Velex. This is achieved by optimizing the coe ﬃ cients used in their formula to make the formula more accurate to evaluate gear body-induced tooth deﬂections. Meanwhile, we introduce a new term called a ﬃ liated body sti ﬀ ness to evaluate the body-induced tooth deﬂections in the double-tooth-pair meshing period. The improved method can give higher accuracy in evaluating gear body-induced tooth deﬂections of spur gears with a pressure angle of 20 ◦ .


Introduction
The magnitude and variation of tooth pair compliance affect tooth loading and gear dynamics significantly [1]. For an elastic gear, when a force F is applied on a gear tooth, a tooth deflection δ is generated along the line of force as illustrated in Figure 1. The tooth deflection can be calculated using the finite element method (FEM) or analytical method [2,3]. With the development of computer science technology, the FEM can give very accurate results if the boundary conditions and the degree of discretization are properly defined [4]. Meanwhile, many researchers have been developing analytical formulas to achieve the fast and accurate calculation of gear tooth deflections. Analytical formulas are easy to use and time-saving in computation [5]. In Ref. [1], gear teeth were assumed to be rigid in evaluating gear-body induced tooth deflections in their finite element model. This assumption will generate some errors. We will remove this assumption in our analysis and give a detailed explanation of how to develop a finite element model to evaluate gear-body induced tooth deflections.
In the existing analytical methods, tooth deflection is considered from two parts, the deflection from gear tooth (assuming rigid gear body) and the deflection caused by gear body. One popular analytical method to estimate the deflection from the first part is the potential energy method [6,7], in which the gear tooth is modeled as a cantilever beam of variable cross-section and the energy stored in the beam is considered as the summation of Hertzian contact energy, bending energy, shear energy, and axial compressive energy. Many researchers have applied this method to evaluate the time-varying meshing stiffness of gears [8][9][10][11]. The deflection induced by a gear body is also called tooth fillet/foundation deflection [1]. Weber [12], Attia [13], Cornell [14], and Sainsot and Velex [1] developed analytical formulas to calculate the gear body-induced tooth deflections. Their analytical formulas all have the same format as Equation (1)  cos 1 tan where δ is tooth deflection under a force F as shown in Figure 1, E represents Young's modulus of gear body material, α denotes the gear pressure angle, b is the tooth face width, u represents the distance between the root circle and the intersection of the tooth center-line and the gear meshing action line (see Figure 1), h denotes the ratio of gear root circle radius (Rf) and gear bore radius (Ra), θf represents the angle between tooth center-line and the junction with the root circle, and Sf is the arc length corresponding to the angle θf. Table 1 gives a summary of the L, M, P, and Q values used in the four papers [1,[12][13][14]. Authors in Refs. [12][13][14] used constant L, M, P, and Q values in their formulas while Sainsot and Velex [1] considered L, M, P, and Q as functions of h and θf. The formulas reported in Ref. [1] were developed based on the theory of Muskhelishvili [15] applied to circular elastic rings with the assumption of linear and constant stress variation at the root circle. It has been widely used in the gear meshing stiffness evaluation even for gears with a tooth crack [16][17][18][19]. Cirelli et al. [20] studied non-linear dynamics of spur gears using a multi-body model. In their model, the gear body induced tooth deflection was calculated using the model proposed by Sainsot et al. [1]. Vivet et al. [21] gave one analytical formula to evaluate gear body induced tooth deflections. But, they did not give details how the formula was derived. Cappellini et al. [22] developed a combined analytic-numerical contact In Ref. [1], gear teeth were assumed to be rigid in evaluating gear-body induced tooth deflections in their finite element model. This assumption will generate some errors. We will remove this assumption in our analysis and give a detailed explanation of how to develop a finite element model to evaluate gear-body induced tooth deflections.
In the existing analytical methods, tooth deflection is considered from two parts, the deflection from gear tooth (assuming rigid gear body) and the deflection caused by gear body. One popular analytical method to estimate the deflection from the first part is the potential energy method [6,7], in which the gear tooth is modeled as a cantilever beam of variable cross-section and the energy stored in the beam is considered as the summation of Hertzian contact energy, bending energy, shear energy, and axial compressive energy. Many researchers have applied this method to evaluate the time-varying meshing stiffness of gears [8][9][10][11]. The deflection induced by a gear body is also called tooth fillet/foundation deflection [1]. Weber [12], Attia [13], Cornell [14], and Sainsot and Velex [1] developed analytical formulas to calculate the gear body-induced tooth deflections. Their analytical formulas all have the same format as Equation (1) but with different L, M, P, and Q values.
where δ is tooth deflection under a force F as shown in Figure 1, E represents Young's modulus of gear body material, α denotes the gear pressure angle, b is the tooth face width, u represents the distance between the root circle and the intersection of the tooth center-line and the gear meshing action line (see Figure 1), h denotes the ratio of gear root circle radius (R f ) and gear bore radius (R a ), θ f represents the angle between tooth center-line and the junction with the root circle, and S f is the arc length corresponding to the angle θ f . Table 1 gives a summary of the L, M, P, and Q values used in the four papers [1,[12][13][14]. Authors in Refs. [12][13][14] used constant L, M, P, and Q values in their formulas while Sainsot and Velex [1] considered L, M, P, and Q as functions of h and θ f . The formulas reported in Ref. [1] were developed based on the theory of Muskhelishvili [15] applied to circular elastic rings with the assumption of linear and constant stress variation at the root circle. It has been widely used in the gear meshing stiffness evaluation even for gears with a tooth crack [16][17][18][19]. Cirelli et al. [20] studied non-linear dynamics of spur gears using a multi-body model. In their model, the gear body induced tooth deflection was calculated using the model proposed by Sainsot et al. [1]. Vivet et al. [21] gave one analytical formula to evaluate gear body induced tooth deflections. But, they did not give details how the formula was derived. Cappellini et al. [22] developed a combined analytic-numerical contact model to simulate gear dynamics, but they did not specifically investigate gear body induced tooth deflections.  [13] 0.32 Cornell [14] 5.306 1.4 (plane stress) 1.534 0.32 1.14 (plane strain) Sainsot and Velex [1] L However, the analytical equations derived for L, M, P, and Q in Ref. [1] are very complicated and hard to use. To address this problem, Sainsot and Velex [1] numerically calculated the values of L, M, P and Q based on the plane strain assumption over a realistic range of values (h between 1.4 and 7, θ f between 0.01 and 0.12) and fitted these values using polynomial functions of the form as shown in where the constants A i , B i , C i , D i , E i , and F i are given in Table 2. Table 2. Coefficients of the polynomial functions [1]. For a pair of spur gears with a contact ratio between one and two, single-tooth-pair meshing and double-tooth-pair meshing take place alternatively as shown in Figure 2. The formula reported in Ref. [1] can only be used in the single-tooth-pair meshing duration but not in the double-tooth-pair meshing duration as their model only considered a single loading force on a gear. To address this shortcoming, Xie et al. [23] introduced a correction factor to consider the load dependence in the double-tooth-pair meshing duration and to account for the errors of using Equation (1) in the double-tooth-pair meshing duration. Meanwhile, Xie et al. [24] derived new formulas for gear body induced deflections considering cubic and parabolic stress distribution around the tooth root and also load dependence between teeth. Chen et al. [25] extended the formula reported in Ref. [1] by considering the effect of tooth root crack on gear body-induced tooth deflections.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 18 model to simulate gear dynamics, but they did not specifically investigate gear body induced tooth deflections.  [13] 0.32 Cornell [14] 5.306 1.4 (plane stress) 1.534 0.32 1.14 (plane strain) Sainsot and Velex [1] However, the analytical equations derived for L, M, P, and Q in Ref. [1] are very complicated and hard to use. To address this problem, Sainsot and Velex [1] numerically calculated the values of L, M, P and Q based on the plane strain assumption over a realistic range of values (h between 1.4 and 7, θf between 0.01 and 0.12) and fitted these values using polynomial functions of the form as shown in Equation (2) θ θ θ θ where the constants Ai, Bi, Ci, Di, Ei, and Fi are given in Table 2. For a pair of spur gears with a contact ratio between one and two, single-tooth-pair meshing and double-tooth-pair meshing take place alternatively as shown in Figure 2. The formula reported in Ref. [1] can only be used in the single-tooth-pair meshing duration but not in the double-tooth-pair meshing duration as their model only considered a single loading force on a gear. To address this shortcoming, Xie et al. [23] introduced a correction factor to consider the load dependence in the double-tooth-pair meshing duration and to account for the errors of using Equation (1) in the doubletooth-pair meshing duration. Meanwhile, Xie et al. [24] derived new formulas for gear body induced deflections considering cubic and parabolic stress distribution around the tooth root and also load dependence between teeth. Chen et al. [25] extended the formula reported in Ref. [1] by considering the effect of tooth root crack on gear body-induced tooth deflections.  Based on our literature review, the formula reported in Ref. [1] is most widely used by researchers to calculate gear body-induced deflections. The formula reported in Ref. [24] is an improved version of the formula reported in Ref. [1]. However, these methods still generate a large error in calculating gear body-induced deflections. There are mainly two reasons for the errors. First, the formulas were derived considering the gear body as an elastic ring without considering the effect of teeth on body deflections. Second, the stress distribution assumed in deriving these formulas is not proper. Linear and constant stress variation around the tooth root circle was assumed in Ref. [1] while cubic and parabolic stress distributions were assumed in Ref. [24]. Based on our finite element analysis, these assumptions are far from real situations. Detailed descriptions will be given in this study.
In this paper, we will discuss how to develop an accurate finite element model, and analyze the effect of different parameters such as gear number of teeth, location of loads, and gear geometry on gear body-induced deflections. Then, an improved analytical method will be proposed to deliver a more accurate estimation of gear body-induced deflections. The remaining part of this paper is organized as follows. Section 2 describes our finite element model and investigates the effect of gear parameters on gear body-induced deflections. Section 3 presents the proposed method to derive an improved formula to evaluate gear body-induced tooth deflections. Section 4 gives case studies to validate the accuracy of the proposed method. Section 5 draws a conclusion.

Finite Element Analysis of Gear Body-Induced Tooth Deflections
In this section, we will investigate the effect of important factors such as loading location and gear geometry such as the number of teeth and gear bore size on gear body-induced tooth deflections. A finite element model is developed to conduct the analysis. The knowledge generated by this analysis can also be used to develop an analytical or numerical solution to evaluate gear body-induced tooth deflections.

Comparisons between Rigid Teeth and Elastic Teeth Finite Element Models
In [1,24], finite element models were developed for gears with the assumption of the elastic body and rigid teeth to evaluate gear body-induced deflections. In these models, a contact force F was applied along the action line and the tooth deflection δ was considered as the deflection caused by gear body along the action line. The rigid teeth in these finite element models generate large errors because these rigid teeth increase gear body stiffness. In our analysis, we will use elastic teeth and body to demonstrate the difference.
Through our analysis, we find there are two ways to accurately obtain gear body-induced deflections through FEM. We will describe them one by one in the following two paragraphs.
Method 1: Build an elastic gear (elastic teeth and body) using FEM and apply a contact force F on a tooth (see Figure 1). Extract acting forces generated by the force F on the tooth root (see the acting forces in Figure 3). Then, apply the same acting forces on the root of the gear rather than apply the force F on a rigid tooth, and record gear tooth deflection δ b along the action line. In this way, δ b is the gear body-induced tooth deflection (the deflection induced by the tooth is not included) corresponding to the contact force F. Method 2: Build an elastic gear (elastic teeth and body) using FEM and apply a contact force F on a tooth (see Figure 1). Record the tooth deflection δ along the action line that is the deflection induced by both gear teeth and gear body. Then, in the finite element model, change the gear body to be rigid. Apply the same contact force F on the same location of the tooth and record the tooth deflection δ t along the action line that is the deflection caused by the gear tooth. The gear body-induced tooth deflection corresponding to the contact force F can be calculated as Method 2: Build an elastic gear (elastic teeth and body) using FEM and apply a contact force F on a tooth (see Figure 1). Record the tooth deflection δ along the action line that is the deflection induced by both gear teeth and gear body. Then, in the finite element model, change the gear body to be rigid. Apply the same contact force F on the same location of the tooth and record the tooth deflection δt along the action line that is the deflection caused by the gear tooth. The gear bodyinduced tooth deflection corresponding to the contact force F can be calculated as δb = δ − δt.
A finite element model reported in Ref. [11] is modified to test the performance of our two methods and the method used in Refs. [1,24]. The gears are modeled using the element type SOLID185 and the element shape of mapped hexahedral. Six teeth are refined when the middle one or two teeth are in meshing. In our model, meshing nodes of the pinion and gear are defined theoretically and then coupled. Gear parameters are given in Table 3 and Figure 3 shows the finite element model (a two-dimensional plane strain model) developed in Ansys.  Figure 4 gives a comparison of three methods: the method reported in [1,24] and our methods 1 and 2. A constant force F (contact force per unit of tooth width F/b = 10 N/m) is applied on a single gear tooth from the tooth tip to the tooth root. The force F mimics the meshing force between a pair of gears with the transmission ration 1:1. The x-axis is the gear angular displacement from tooth tip to tooth root. The tooth tip corresponds to the angle 0 at the very left of the x-axis while the tooth root corresponds to the angle at the very right of the x-axis. The y-axis is gear-body induced deflection along the action line. From the comparison, we can see that the rigid teeth in the finite element model will lead to less tooth deflection. The rigid teeth make the gear more rigid, which generates errors. Our two methods give almost the same result and they are more accurate than the model reported in [1,24]. A finite element model reported in Ref. [11] is modified to test the performance of our two methods and the method used in Refs. [1,24]. The gears are modeled using the element type SOLID185 and the element shape of mapped hexahedral. Six teeth are refined when the middle one or two teeth are in meshing. In our model, meshing nodes of the pinion and gear are defined theoretically and then coupled. Gear parameters are given in Table 3 and Figure 3 shows the finite element model (a two-dimensional plane strain model) developed in Ansys. Table 3. A pair of gears with the transmission ratio 1:1.

Number of Teeth
Gear Module Bore Radius Pressure Angle 19 3.175 17.5 20 • Figure 4 gives a comparison of three methods: the method reported in [1,24] and our methods 1 and 2. A constant force F (contact force per unit of tooth width F/b = 10 N/m) is applied on a single gear tooth from the tooth tip to the tooth root. The force F mimics the meshing force between a pair of gears with the transmission ration 1:1. The x-axis is the gear angular displacement from tooth tip to tooth root. The tooth tip corresponds to the angle 0 at the very left of the x-axis while the tooth root corresponds to the angle at the very right of the x-axis. The y-axis is gear-body induced deflection along the action line. From the comparison, we can see that the rigid teeth in the finite element model will lead to less tooth deflection. The rigid teeth make the gear more rigid, which generates errors. Our two methods give almost the same result and they are more accurate than the model reported in [1,24]. Figure 5 gives deformation and stress information (obtained from our FEM model) along a tooth root when a meshing force (10 N/m) is applied on the tooth tip. In Figure 5, the x-axis represents the range of tooth angle on the root circle with the angle 0 corresponding to the tooth central line.
The very left and right angles in the x-axis of Figure 5 correspond to tooth fillet. When the teeth are modeled as rigid, radial deformation changes linearly along the circumferential direction of the tooth root, and the circumferential deformation is almost constant. When the teeth are elastic, near a sinusoidal radial deformation is observed along the circumferential direction of the tooth root, and the circumferential deformation is almost a constant except for linear variation near the tooth fillet. In the middle part, the normal stress in the radial direction varies linearly and the shearing stress along the circumferential direction is almost constant, but, large stress variation occurs nears tooth fillets. In addition, the change trend is quite different between a rigid tooth and an elastic tooth. The big differences in the deformation and stress demonstrate that it is not proper to use rigid teeth in FEM to obtain gear body-induced deformations.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 18  Figure 5 gives deformation and stress information (obtained from our FEM model) along a tooth root when a meshing force (10 N/m) is applied on the tooth tip. In Figure 5, the x-axis represents the range of tooth angle on the root circle with the angle 0 corresponding to the tooth central line. The very left and right angles in the x-axis of Figure 5 correspond to tooth fillet. When the teeth are modeled as rigid, radial deformation changes linearly along the circumferential direction of the tooth root, and the circumferential deformation is almost constant. When the teeth are elastic, near a sinusoidal radial deformation is observed along the circumferential direction of the tooth root, and the circumferential deformation is almost a constant except for linear variation near the tooth fillet. In the middle part, the normal stress in the radial direction varies linearly and the shearing stress along the circumferential direction is almost constant, but, large stress variation occurs nears tooth fillets. In addition, the change trend is quite different between a rigid tooth and an elastic tooth. The big differences in the deformation and stress demonstrate that it is not proper to use rigid teeth in FEM to obtain gear body-induced deformations.     Figure 5 gives deformation and stress information (obtained from our FEM model) along a tooth root when a meshing force (10 N/m) is applied on the tooth tip. In Figure 5, the x-axis represents the range of tooth angle on the root circle with the angle 0 corresponding to the tooth central line. The very left and right angles in the x-axis of Figure 5 correspond to tooth fillet. When the teeth are modeled as rigid, radial deformation changes linearly along the circumferential direction of the tooth root, and the circumferential deformation is almost constant. When the teeth are elastic, near a sinusoidal radial deformation is observed along the circumferential direction of the tooth root, and the circumferential deformation is almost a constant except for linear variation near the tooth fillet. In the middle part, the normal stress in the radial direction varies linearly and the shearing stress along the circumferential direction is almost constant, but, large stress variation occurs nears tooth fillets. In addition, the change trend is quite different between a rigid tooth and an elastic tooth. The big differences in the deformation and stress demonstrate that it is not proper to use rigid teeth in FEM to obtain gear body-induced deformations.

Effect of Loading and Gear Parameters on Tooth Root Stress
Linear and constant stress variation around tooth root circle was assumed in Ref. [1] while cubic and parabolic stress distribution was used in Ref. [24]. These assumptions may not be accurate. This sub-section will reveal more realistic stress distribution around tooth root circle and evaluate the effect of loading location and gear parameters on stress distribution.

Effect of Loading Location
A force F (contact force per unit of tooth width F/b = 10 N/m) is applied on four locations of a tooth along the tooth action line (see Figure 6) to test the effect of loading locations on tooth root stress. Locations p 1 and p 2 lie in the double-meshing period while locations p 3 and p 4 lie in the single-meshing period. These four locations are specified using the values given in Table 4. effect of loading location and gear parameters on stress distribution.

Effect of Loading Location
A force F (contact force per unit of tooth width F/b = 10 N/m) is applied on four locations of a tooth along the tooth action line (see Figure 6) to test the effect of loading locations on tooth root stress. Locations p1 and p2 lie in the double-meshing period while locations p3 and p4 lie in the singlemeshing period. These four locations are specified using the values given in Table 4.  Figure 1: An elastic tooth subjected to a force F) 6.4977 3.4018 3.1860 2.1688   Figure 8 gives shearing stress distribution. As the loading location moves from tooth tip to tooth root, the normal stress and shearing stress both decrease. The normal stress distribution is like a sinusoidal curve while the shearing stress distribution looks like a mountain valley. The tooth root stress distribution is not linear or cubic or parabolic.

Effect of Number of Teeth on Stress Distribution
This section investigates the effect of number of teeth on gear tooth root stress distribution. Gear parameters are given in Table 3 except we test 3 different teeth. A force F (contact force per unit of tooth width F/b = 10 N/m) is applied on the tooth tip.
With the increase of the number of teeth, the tooth root normal stress and shearing stress decrease slightly, and the circumferential angle range covered by a single tooth will become smaller, as shown in Figure 9.

Effect of Number of Teeth on Stress Distribution
This section investigates the effect of number of teeth on gear tooth root stress distribution. Gear parameters are given in Table 3 except we test 3 different teeth. A force F (contact force per unit of tooth width F/b = 10 N/m) is applied on the tooth tip.
With the increase of the number of teeth, the tooth root normal stress and shearing stress decrease slightly, and the circumferential angle range covered by a single tooth will become smaller, as shown in Figure 9.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 18 Figure 9. Effect of number of teeth on tooth root stress distribution.

Effect of Gear Bore Size on Tooth Root Stress Distribution
In this section, we only change the gear bore size to achieve different h values: 1.5, 2.1, and 3.5. A larger h value corresponds to a smaller gear bore size. From our analysis results, we find the effect of gear bore size on tooth root stress distribution is very small, as shown in Figure 10.

Effect of Gear Bore Size on Tooth Root Stress Distribution
In this section, we only change the gear bore size to achieve different h values: 1.5, 2.1, and 3.5. A larger h value corresponds to a smaller gear bore size. From our analysis results, we find the effect of gear bore size on tooth root stress distribution is very small, as shown in Figure 10. Figure 9. Effect of number of teeth on tooth root stress distribution.

Effect of Gear Bore Size on Tooth Root Stress Distribution
In this section, we only change the gear bore size to achieve different h values: 1.5, 2.1, and 3.5. A larger h value corresponds to a smaller gear bore size. From our analysis results, we find the effect of gear bore size on tooth root stress distribution is very small, as shown in Figure 10.  Similarly, we also analyzed the effect of gear modulus on tooth root stress distribution. We find its effect is negligible.
Overall, the gear tooth root stress is largely affected by the loading location, and slightly affected by the number of teeth and gear bore size. The effect of gear modulus on tooth root stress is negligible.

An Improved Formula to Evaluate Gear-Body Induced Deflections
In this section, we develop an improved formula to evaluate the gear body-induced tooth deflections. For the single-tooth-pair meshing duration, Equations (1) and (2) are still used following Ref.
[1] but we will update the coefficients A i , B i , C i , D i , E i , and F i of Equation (2) by optimizing them according to our finite element results. For the double-tooth-pair meshing duration, Equation (1) is not proper to use directly. We will give a new formula.

Single-Tooth-Pair Meshing Duration
This sub-section aims to optimize the coefficients A i , B i , C i , D i , E i , and F i of Equation (2) in order to get higher accuracy when we evaluate gear body induced tooth deflection in the single-tooth-pair meshing duration. The objective function of this optimization problem is to minimize the difference between the results generated by Equation (1) and a well-established finite element model. A finite element model reported in Ref. [11] is modified to evaluate the gear body-induced tooth deflections of 72 standard involute spur gears. The gears are modeled using the element type SOLID185 and the element shape of mapped hexahedral. Six teeth are refined when the middle one or two teeth are in meshing. In our model, meshing nodes of the pinion and gear are coupled. ANSYS is used to implement the finite element analysis. These gears (72 = 9 × 8) are a combination of of 9 teeth and 8 different h values, as shown in Table 5. These gears cover a realistic range of θ f and h values: θ f between 0.03 and 0.15 (radian), and h between 2.1 and 7 (remember that L, M, P and Q are functions of θ f and h). All the gears have a module of 3.175 and a pressure angle of α = 20 • . Gear body-induced tooth deflections of these 72 gears are evaluated using the finite element model to optimize the coefficients in Equation (2). Based on our finite element analysis, gear body-induced tooth deflections are independent of the module of gears. This conclusion can also be confirmed by the analytical formulas reported in Ref. [1]. We can see that the module of gears is not a variable in Equation (1). Therefore, our equations derived in this paper can be used for a spur gear with any modulus. However, since we only used the pressure angle of α = 20 • in our finite element analysis, the optimized coefficients may not give accurate results for other pressure angles of a spur gear.
In the reported finite element models [1,23], gear teeth were assumed as rigid in evaluating the gear body-induced tooth deflections. However, this assumption causes errors of about 10% in the gear body-induced tooth deflections based on our analysis. This assumption is released in our finite element analysis. When a force F is applied on the tooth profile (see Figure 3), we first extract acting forces generated by the force F on the tooth root. Then, we apply the same acting forces on the root of the gear and record gear tooth deflection δ along the force action line. In this way, δ is the gear body-induced tooth deflection and it is obtained from an elastic gear (elastic gear teeth and body).
A two-dimensional (2-D) plane strain model is used in our finite element analysis with element type PLANE183 as it can give similar accuracy in meshing stiffness evaluation as a 3-D finite element model. We use the 2-D model because it requires much lower computational cost. To demonstrate that it is acceptable to use the plane strain conditions, we did tests on gears with the same parameters (number of teeth: 19, gear module: 3.175, bore radius: 17.5 mm, meshing force per unit of tooth width F/b = 10 N/m, transmission ratio 1:1, pressure angle α = 20 • ) except for the tooth face width. Figure 11 shows the comparison results. The x-axis of Figure 11 is the gear angular displacement corresponding to the meshing force moving from the tooth tip to root. The y-axis is the gear body stiffness (some researchers call it body stiffness) per unit of tooth width which is expressed as k b = F/b δ , where F is the meshing force, b represents tooth width, and δ denotes the gear body-induced tooth deflection. In the whole process, F is applied on a single gear tooth since we focus on the single-tooth-pair meshing duration in this sub-section. The force F mimics the meshing force between a pair of gears with the transmission ration 1:1. From Figure 11, we can see that when the tooth width is above 20 mm, the error generated by the 2-D plane strain model is negligible. Even when the tooth width is 10 mm, the error is less than 10%. It is acceptable to use the 2-D plane strain model in our finite element analysis. When the tooth width is much smaller than 10 mm, we do not recommend using the tooth deflection equations derived in this paper.
Using the above described 2-D plane strain finite element model, we evaluated the gear body stiffness for 72 standard involute spur gears as given in Table 5. Then the Genetic Algorithm was used to optimize the 24 coefficients that are required in Equation (2). The objective of the optimization is to let the gear body-induced tooth deflections evaluated using the Equation (1) closest to the results obtained from the finite element analysis. Table 6 gives the updated coefficients obtained from our optimization.  tooth deflection. In the whole process, F is applied on a single gear tooth since we focus on the singletooth-pair meshing duration in this sub-section. The force F mimics the meshing force between a pair of gears with the transmission ration 1:1. From Figure 11, we can see that when the tooth width is above 20 mm, the error generated by the 2-D plane strain model is negligible. Even when the tooth width is 10 mm, the error is less than 10%. It is acceptable to use the 2-D plane strain model in our finite element analysis. When the tooth width is much smaller than 10 mm, we do not recommend using the tooth deflection equations derived in this paper. Figure 11. Comparison between 2-D and 3-D models.
Using the above described 2-D plane strain finite element model, we evaluated the gear body stiffness for 72 standard involute spur gears as given in Table 5. Then the Genetic Algorithm was used to optimize the 24 coefficients that are required in Equation (2). The objective of the optimization is to let the gear body-induced tooth deflections evaluated using the Equation (1) closest to the results obtained from the finite element analysis. Table 6 gives the updated coefficients obtained from our optimization. A comparison is conducted between our proposed coefficients and the coefficients reported in Ref. [1]. We use these coefficients to obtain the gear body-induced tooth deflections analytically, respectively. The y-axis of Figure 12 gives the maximum relative difference of the gear body stiffness per unit of tooth width between the analytical method and the 2-D plane strain finite element model. The maximum error generated using our proposed coefficients is below 10% whereas the coefficients reported in Ref. [1] can give the error as high as 20%. In 50 out of the total 72 cases (see Table 5), the errors generated using our coefficients are lower than those in Ref. [1]. In the other 22 cases, the method with our coefficients performs a little worse (the average difference between the two methods  A comparison is conducted between our proposed coefficients and the coefficients reported in Ref. [1]. We use these coefficients to obtain the gear body-induced tooth deflections analytically, respectively. The y-axis of Figure 12 gives the maximum relative difference of the gear body stiffness per unit of tooth width between the analytical method and the 2-D plane strain finite element model. The maximum error generated using our proposed coefficients is below 10% whereas the coefficients reported in Ref. [1] can give the error as high as 20%. In 50 out of the total 72 cases (see Table 5), the errors generated using our coefficients are lower than those in Ref. [1]. In the other 22 cases, the method with our coefficients performs a little worse (the average difference between the two methods is smaller than 1.5%). Overall, our proposed coefficients give more accurate results in evaluating the gear body-induced tooth deflections in the single-tooth-pair meshing duration.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 18 is smaller than 1.5%). Overall, our proposed coefficients give more accurate results in evaluating the gear body-induced tooth deflections in the single-tooth-pair meshing duration.

Double-Tooth-Pair Meshing Duration
In Section 2.1, we optimized the 24 coefficients used in Equation (1) that was derived by Sainsot et al. [1]. Our coefficients give better performance in the evaluation of gear body-induced tooth deflections during the single-tooth-pair meshing period. However, this analytical formula is not proper to use directly for the double-tooth-pair meshing duration. This section aims to derive an

Double-Tooth-Pair Meshing Duration
In Section 2.1, we optimized the 24 coefficients used in Equation (1) that was derived by Sainsot et al. [1]. Our coefficients give better performance in the evaluation of gear body-induced tooth deflections during the single-tooth-pair meshing period. However, this analytical formula is not proper to use directly for the double-tooth-pair meshing duration. This section aims to derive an improved formula for the double-tooth-pair meshing duration.
In the double-tooth-pair meshing duration, there are two pairs of teeth meshing simultaneously. Figure 13 is a simple illustration of gear meshing forces and their resulting defections. In this illustration, we assume gear teeth do not have tooth deflection since our focus is only the gear body deflections. In the figure, F 1 and F 2 represent the gear meshing forces of the first and the second tooth pairs, respectively. And P 1 and P 2 denote the gear meshing points of the first and second gear tooth pairs, respectively. Imagine if the first tooth pair starts to mesh first and, correspondingly, a deflection ∆ 11 is generated in the mesh point P 1 and a deflection ∆ 21 is generated in the meshing point P 2 . Later, when the second tooth pair starts to mesh (the first tooth pair is still in meshing), further deflections ∆ 12 and ∆ 22 are generated in the meshing points P 1 and P 2 , respectively. Therefore, the total deflection is ∆ 11 + ∆ 12 (δ 1 = ∆ 11 + ∆ 12 ) at the meshing point P 1 and ∆ 21 + ∆ 22 (δ 2 = ∆ 21 + ∆ 22 ) at the meshing point P 2 . The values of ∆ 11 and ∆ 22 can be obtained from the analytical equation derived for the single-tooth-pair meshing duration. Then, our remaining problem is to get the values of ∆ 12 and ∆ 21 .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 18 Figure 13. Illustration of meshing forces and the resulting deflections.
As we described in the previous paragraph, if we apply the force F1 followed by F2, the energy stored in the system can be obtained as follows: Similarly, if we apply the force F2 followed by F1, the energy stored in the system can be expressed as Equation (4) based on one assumption that the deflection is proportional to the load. This assumption is generally valid for elastic materials with small deflection.
From Equations (3) and (4), we can conclude: Further, we can get: where  denotes the load sharing ratio of a pair of spur gears with the expression Figure 13. Illustration of meshing forces and the resulting deflections.
As we described in the previous paragraph, if we apply the force F 1 followed by F 2 , the energy stored in the system can be obtained as follows: Similarly, if we apply the force F 2 followed by F 1 , the energy stored in the system can be expressed as Equation (4) based on one assumption that the deflection is proportional to the load. This assumption is generally valid for elastic materials with small deflection.
From Equations (3) and (4), we can conclude: Further, we can get: where λ denotes the load sharing ratio of a pair of spur gears with the expression λ = F 1 F 1 +F 2 . Equation (6) gives the relationship between ∆ 12 and ∆ 21 (see Figure 13 for details about these two deflections).
For better illustration, ∆ 21 can be expressed as follows where k A 21 is the ratio of F 1 and ∆ 21 . We call k A 21 affiliated body stiffness in this paper. Similarly, we can define the other affiliated body stiffness k A 12 with the expression k A 12 = F 2 ∆ 12 . From Equation (5), we can find that k A 21 is equal to k A 12 . From Equation (7), we can see that the deflection ∆ 21 can be obtained if the affiliated body stiffness is known. Then, from Equation (6), the deflection of ∆ 12 is available if ∆ 21 and the load sharing ratio λ are known.
We used the above-mentioned 2-D plane strain finite element model to evaluate the affiliated body stiffness for the 72 gears with the parameters given in Table 5. Figure 14 gives an example to illustrate the affiliated body stiffness of two gear pairs. Parameters of the two gear pairs are given in Table 7. The affiliated body stiffness is not constant while it is almost linearly proportional to the gear rotation angle (we know gear meshing point moves as the gear rotates). But, we find the deviation of the affiliated body stiffness values from their average value is small. The first gear pair has a maximum deviation of −3.18% while the second gear pair owns a maximum deviation of 3.23% (see Figure 14). For all other gears, the maximum deviation is around or below 3%. Therefore, we believe it is reasonable to use the average value of the affiliated body stiffness values of gear to represent this gear's affiliated body stiffness. Some errors will be generated, but in most cases, the error is around or below 5%. Table 8 gives the average of the affiliated body stiffness values per unit of tooth width for 72 gears. Notably, the affiliated body stiffness is independent of one gear's module, the same as the gear body stiffness.  Notably, the affiliated body stiffness is independent of one gear's module, the same as the gear body stiffness.

Meshing Stiffness of a Pair of Gears
According to the potential energy method [2], the meshing stiffness of a pair of spur gears can be evaluated as follows: , double tooth pairs in meshing sin gle tooth pair in meshing (8) where F 1 and F 2 are gear meshing forces, and δ 1 , δ 1 , δ 2 , and δ 2 are tooth deflections along the action line corresponding to the driving gear tooth of the first tooth pair, the driven gear tooth of the first tooth pair, the driving gear tooth of the second tooth pair, and the driven gear tooth of the second tooth pair, respectively. Each of the tooth deflections (δ 1 , δ 1 , δ 2 , and δ 2 ) can be considered as the summation of two parts: the deflection from gear tooth (assuming rigid gear body) and the deflection caused by gear body. The deflection from gear tooth (assuming rigid gear body) can be evaluated using the potential energy method [3] or the FEM [11]. The deflections from the gear body can be evaluated using the method proposed in this paper. In the single-tooth-pair meshing duration, gear body induced tooth deflection can be evaluated using Equation (1) with the parameters given in Table 6. In the double-tooth-pair meshing duration, the gear body induced tooth deflection needs to be considered as the summation of two parts. Part 1 can be evaluated using Equation (1) with the parameters given in Table 6. Part 2 needs to be evaluated using Equations (6) and (7) with the affiliated body stiffness given Table 8. In this process, the load sharing ratio of a pair of gears needs to be known.

Validation
Two pairs of involute spur gears are used to validate the performance of our proposed method. The parameters of the gears are given in Table 9. The torque per unit of tooth width applied on the driving gear is 10 N.m/m.  Figure 15 gives the comparisons between the proposed method and the FEM in evaluating the meshing stiffness of two gear pairs. The 2-D finite element model described in Section 2 is used. The red lines in Figure 15 are obtained from the 2-D finite element model. The blue lines are generated considering two parts. Part 1 is the gear body induced deflection that is evaluated using the analytical method proposed in Section 2. Part 2 is the deflection from gear tooth (assuming rigid gear body) that is also obtained from the 2-D finite element model. We use the finite element model to evaluate the tooth induced deflections because our focus is to show the performance of the analytical method for evaluating gear body induced deflections.
From the comparison, we can see that the difference between the two methods in the single-tooth-pair meshing duration is very small. The difference in the double-tooth-pair meshing duration is larger but it is still below 10%. We think a 10% difference is still acceptable. A big reason for the difference may because we used the average value of the affiliated body stiffness in the meshing stiffness calculation. Actually, it is time-varying. A time-varying function may give results that are more accurate in gear meshing stiffness evaluation.
Our study also has limitations in that we did not give affiliated body stiffness values when gears are rotating fast with transmission errors. This will be part of our future work. This is also a challenging work. After this challenge is addressed, we will further investigate gear dynamics such as transmission errors in the future.
From the comparison, we can see that the difference between the two methods in the singletooth-pair meshing duration is very small. The difference in the double-tooth-pair meshing duration is larger but it is still below 10%. We think a 10% difference is still acceptable. A big reason for the difference may because we used the average value of the affiliated body stiffness in the meshing stiffness calculation. Actually, it is time-varying. A time-varying function may give results that are more accurate in gear meshing stiffness evaluation.

Conclusions
This study gives a detailed investigation of gear body-induced tooth deflections and presents an improved analytical solution to evaluate gear body-induced tooth deflections. First, we find it is not proper to assume rigid teeth in finite element models to evaluate gear body-induced tooth deflections.
The assumption generates a large error. We also proposed two ways to accurately evaluate gear body-induced tooth deflections using the finite element method. Then, through finite element analysis, we find that loading location has a large effect on gear body induced tooth root stress distribution, the effect from the number of teeth and the bore size is small, and the effect from gear modulus is negligible. Gear tooth root stress distribution is not linear or cubic or parabolic. The gear tooth normal stress caused by a gear body is kind of a sinusoidal curve while the shearing stress looks like a mountain valley. Last but not least, an improved analytical solution is developed to accurately evaluate gear body-induced deflections. In the single-tooth-pair meshing duration, we optimized the coefficients of the formula reported in Ref. [1]. In the double-tooth-pair meshing duration, we introduced the concept of affiliated body stiffness and provided a table to obtain the affiliated body stiffness. The improved analytical solution is more accurate than existing analytical methods to evaluate the gear body-induced tooth deflections, which can also be used to incorporate the effect of gear body in gear meshing stiffness evaluation. However, our study is limited to spur gears with a pressure angle of 20 • . In the future, we will investigate gear body-induced tooth deflections for spur gear with other pressure angles and other types of gears such as helical gear and bevel gear.