Contact stresses and bending stresses analysis of curvilinear cylindrical gears generated by a face-milling cutter with parabolic profile

<!-- metadata commented in wiki content


<span id='_Hlk41161685'></span><span id='_Hlk41225231'></span><span id='OLE_LINK2'></span><span id='_Hlk38552328'></span><span id='_Hlk42204287'></span><span id='_Hlk41222095'></span><span id='_Hlk41126098'></span><big>'''Contact stresses and bending stresses analysis of curvilinear cylindrical gears generated by a face-milling cutter with parabolic profile'''</big>

<span id='OLE_LINK8'></span><span id='OLE_LINK9'></span>'''Zhang Xuegang <sup>a</sup> '''(Corresponding author)

a. School of Intelligent Manufacturing, Panzhihua University, Panzhihua 617000, Sichuan, China.

Email: [mailto:xuegang.zhang@yahoo.com xuegang.zhang@yahoo.com].

==Huang Xiaobin<sup>a</sup> ==

a. School of Intelligent Manufacturing, Panzhihua University, Panzhihua 617000, Sichuan, China.

==Zhang Qi<sup>a</sup> ==

a. School of Intelligent Manufacturing, Panzhihua University, Panzhihua 617000, Sichuan, China.

==Gan Bin<sup>b</sup>==

b. Sichuan Provincial Key Lab of Process Equipment and Control, Sichuan University of Science & Engineering, Zigong 643002, Sichuan, China.

<span id='_Hlk41216933'></span>
-->
==Abstract:==

Curvilinear cylindrical gear drives are a type of parallel-axis gear set, and the longitudinal geometrical shape of the gear teeth is a part of the circular arc. This study aims to develop a curvilinear cylindrical gear drive with superior transmission performance.

<span id='_Hlk41225915'></span><span id='_Hlk41226138'></span>This study proposes a curvilinear cylindrical gear drive with a parabolic profile. The gear tooth surfaces are generated by a face-milling cutter with a straight-line profile, and the matching gear tooth surfaces are generated by a face-milling cutter with the parabolic profile. Initially, mathematical models of the gear set are described based on the theory of gearing. The finite element method and tooth contact analysis algorithms were applied to contact and bending stress analyses. Finally, four numerical examples are provided to illustrate the transmission performance of the curvilinear cylindrical gear drive with the parabolic profile. The study results show that curvilinear cylindrical gear drives have superior transmission performance compared to those with a straight-line profile.

==Keywords:==

Curvilinear cylindrical gears; Parabolic profile; Tooth contact analysis; FEM; Stress analysis

==<br/>1. Introduction==

<span id='_Hlk41307839'></span><span id='_Hlk38401023'></span><span id='_Hlk41310383'></span>The curvilinear cylindrical gear is a kind of novel gear which is the same as ‘S’ gear [1], it was first proposed in 1914 [2]. Until the 1980s, curvilinear cylindrical gears were industrially applied in china steel plants, aluminum rolling mills, and cement equipment plants [3]. Owing to the geometrical shape of the gear teeth as a part of a circular arc along the longitudinal direction of the gears, lubricating oil will be retained within the concave tooth surface during operation. Therefore, curvilinear cylindrical gear drives have superior lubrication performance. In addition, the axial thrust forces will be avoided during transmission. Nevertheless, curvilinear cylindrical gear drives cannot be employed for internal gear transmission systems. Therefore, gear drives have not been widely applied in the industry.

<span id='_Hlk41310563'></span>According to the longitudinal pressure angle of gear teeth, curvilinear cylindrical gear drives are to be divided into two categories: curvilinear cylindrical gears with teeth of constant pressure angle (CC-gears) and curvilinear cylindrical gears with teeth of variable pressure angle (CV-gears) [4,5]. According to the contact pattern on the gear tooth surfaces, they have divided into two categories: one type is the curvilinear cylindrical gear drive with line contact, and another type is the curvilinear cylindrical gear drive with point contact.

Regarding curvilinear cylindrical gear drives with line contact, it has been investigated in the aspect of geometric model [6,7], stress analysis [8,9], processing method [10,11], and elastohydrodynamic lubrication [12,13],the  et al. In practice, the slightly angular misalignments will bring out edge contact at the edge of the gear tooth surfaces. Therefore, this type of gear drive is sensitive to assembly errors of gear sets [14]. Usually, the gear tooth surfaces need to be modified. Thereby, the instantaneous line contact of tooth surfaces will turn into the instantaneous point contact in meshing.

<span id='_Hlk38488204'></span><span id='_Hlk38491508'></span><span id='_Hlk41316562'></span><span id='_Hlk38532763'></span><span id='_Hlk38533859'></span><span id='_Hlk41318075'></span><span id='_Hlk41317889'></span><span id='_Hlk41318034'></span>Regarding curvilinear cylindrical gear drives with point contact, Tsay and his team have investigated the generation of the gear sets by two face-milling cutters with straight-line profiles [15,16] and hob cutters [17,18], respectively. Meanwhile, the mathematical model, undercutting conditions, and contact characteristics of these gear sets are investigated by the theory of gearing and tooth contact analysis (TCA). Wu [19] and Chen [20] et al. utilized two complemented circular-arc rack cutters to generate the gear and matching gear, respectively. Chen [21] et al. employed two same circular-arc rack cutters to generate the gear sets. Alfonso et al. [22] proposed that two independent fixed-setting cutters (FSC) are used to generate the concave and convex sides of the pinion tooth surfaces, respectively, and the wheel tooth surfaces are generated by a spread-blade cutter. Additionally, the contact patterns and stresses of the gear sets were compared with gear sets manufactured by two spread-blade cutters (SBC). In [23], four kinds of gear drives (FSC-gears, SBC-gears, spur, and helical gears) were compared in terms of contact stresses and bending stresses.

This paper proposed a kind of curvilinear cylindrical gear generated by a face-milling cutter with the parabolic profile. For this type of gear drive, the main research contents are as follows:

<span id='_Hlk41327331'></span>(1) Generation of the geometry of the wheel tooth surfaces by a face-milling cutter with the straight-line profile. Nevertheless, the generation of the geometry of the pinion tooth surfaces by a face-milling cutter with the parabolic profile.

(2) Under different assembly error conditions, parabola parameters, and cutter mean pitch radius, The TCA algorithm is utilized to investigate the gear set's initial contact position. Meanwhile, contact stresses and bending stresses are also analyzed by the finite element method.

(3) Four numerical group examples are used to expose the transmission performances of the proposed curvilinear cylindrical gear drives.

<span id='_Hlk38137371'></span><span id='_Hlk38136944'></span>
=2. Geometry of curvilinear cylindrical gears=

<span id='sec2_1'></span><span id='_Hlk38137224'></span>''2.1. Generating principle of curvilinear cylindrical gears by face-milling cutters''

<span id='_Hlk39784102'></span><span id='_Hlk39838220'></span><span id='_Hlk39770907'></span><span id='_Hlk39831173'></span><span id='_Hlk39924172'></span>In order to generate the gear tooth surfaces, the first task is the establishment of the virtual face-milling cutter surfaces. A brief generating mechanism of the gear tooth surfaces is shown in Fig. 1. The outer and inner blades of the face-milling cutter are employed for the generation of concave and convex sides of the gear tooth surfaces, respectively. Generating process of the cutter blade surfaces is as follows:

<span id='_Hlk40632260'></span>(1) The outer blade profile (including straight line and fillet) can be represented in coordinate system ''S<sub>ob</sub>'', and the inner blade profile can be represented in coordinate system ''S<sub>ib</sub>''. When the outer and inner blade profiles of the face-milling cutter are represented in coordinate system ''S<sub>b</sub>'', coordinate transformation is necessary from the coordinate system ''S<sub>ob</sub>'' (''S<sub>ib</sub>'') to ''S<sub>b</sub>''. Coordinate transformation matrices are expressed by '''M''' ''<sub>b,ob</sub>'' and '''M''' ''<sub>b,ib</sub>'', which are given in section 2.2.

<span id='_Hlk39924252'></span><span id='_Hlk43228639'></span>(2) When blade profiles of the face-milling cutter are represented in coordinate system ''S<sub>c </sub>'', the virtual face-milling cutter surfaces are expressed by a coordinate transformation matrix '''M''' ''<sub>c,b </sub>''(''θ''). Matrix '''M''' ''<sub>c,b </sub>''(''θ'') is from coordinate system ''S<sub>b</sub>'' to ''S<sub>c </sub>'':

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq1'></span> <math display="inline">\boldsymbol{M}_{c,b}(\theta )=\left[\begin{array}{cccc}
cos\theta  & 0 & -sin\theta  & -Rcos\theta \\
0 & 1 & 0 & 0\\
sin\theta  & 0 & cos\theta  & -Rsin\theta \\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}


Where, symbol ''θ'' expresses the rotation angle of the face-milling cutter during gear cutting. Symbol ''R'' expresses the mean pitch radius of the face-milling cutter.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image2-c.png|600px]] </div>

<div id="fig1" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 1'''. A brief generating mechanism of curvilinear cylindrical gears</div>

After a series of coordinate transformations generate the virtual face-milling cutter surfaces, the gear tooth surfaces can be generated by the curvilinear gear machining principle. The curvilinear gear machining principle is as follows:

<span id='_Hlk39853596'></span><span id='_Hlk39858046'></span><span id='_Hlk39858659'></span><span id='_Hlk39858748'></span>(1) Before processing curvilinear cylindrical gears, the cutter pitch plane always remains tangent to the pitch circle of the gear blank. Meanwhile, the rotation axis of the face-milling cutter has to retain one translational and one rotational degree of freedom, and the rotation axis of the gear blank only retains one rotational degree of freedom. Here, the direction of the rotation axis of the face-milling cutter is the same as axis ''y<sub>c</sub>'', and located at the geometric center of the face-milling cutter, the rotation axis of the gear blank is the same as axis ''z<sub>p</sub>'', and located at the geometric center of the gear blank.

(2) In the process of machining the gear, the face-milling cutter continuously rotates around axis ''y<sub>c</sub>'' with angular velocity ''ω<sub>c</sub>'', which does not effect on the generated geometry of the gear tooth surfaces. Furthermore, the face-milling cutter reciprocates at lineal velocity ''v<sub>c</sub>'' along the direction orthogonal to the gear axis ''z<sub>p</sub>''. Meanwhile, the gear blank is rotated alternately clockwise and counterclockwise by a specific angle'' ψ ''under the angular velocity ''ω<sub>p</sub>''. Namely, the cutter always does pure rolling relative to the gear blank. Then, the moving velocity of the cutter and the angular velocity of the gear blank must satisfy the following equation:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq2'></span> <math display="inline">v_c=r{\omega }_p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}


Here, ''r'' expresses the pitch radius of the to-be-generated curvilinear cylindrical gear.

<span id='_Hlk43228848'></span><span id='_Hlk40429911'></span><span id='_Hlk40375028'></span>(3) When the gear blank moves with the cutter at the same time based on the method in the previous step, the gear tooth surfaces are generated as the envelope to the family of positions of the virtual cutter surfaces. The family of positions of the virtual cutter surfaces can be represented in coordinate system ''S<sub>p</sub>''. The coordinate transformation from coordinate system ''S<sub>c</sub>'' to ''S<sub>p</sub>'' is represented by the following matrix

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq3'></span> <math display="inline">\boldsymbol{M}_{p,c}(\psi )=\left[\begin{array}{cccc}
cos\psi  & sin\psi  & 0 & r(sin\psi -\psi cos\psi )+Rcos\psi \\
-sin\psi  & cos\psi  & 0 & r(cos\psi +\psi sin\psi )-Rsin\psi \\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}


After one tooth space is generated by the face-milling cutter, the gear blank is indexed by one tooth pitch angle, and the generating mechanism of the gear continues to generate the next tooth space, until all tooth spaces have been formed.

<span id='sec2_2'></span><span id='_Hlk41335196'></span><span id='_Hlk41335359'></span>''2.2. Derivation of mathematical model of the wheel tooth surfaces''

<span id='_Hlk40031870'></span><span id='_Hlk40020610'></span><span id='_Hlk40020643'></span><span id='_Hlk40610880'></span>Fig. 1 (Right) shows the straight-line profiles of the face-milling cutter for the generation of the wheel tooth surfaces. The inner and outer blade profiles are divided into two parts: the straight-line profile, which will be used for generation of the working part of the wheel tooth surfaces, and the circular arc profile, which will be used for generation of the fillet part of the wheel tooth surfaces. The variable ''u<sub>g </sub>''represents the distance between arbitrary points on the straight-line profile and the origin of coordinate system ''S<sub>ib</sub>'' or ''S<sub>ob</sub>'', and the variable ''λ<sub>g</sub>'' represents the polar angle of arbitrary points on the circular arc profile.

The working and fillet parts of the wheel tooth surfaces will be derived, respectively.

<span id='_Hlk40612456'></span><span id='_Hlk40001787'></span><span id='_Hlk40031278'></span><span id='_Hlk40031291'></span>The position vectors of a point on the straight-line profiles can be expressed in the coordinate systems ''S<sub>ib</sub>'' and ''S<sub>ob</sub>'', as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq4'></span> <math display="inline">\boldsymbol{r}_{}^{\left(G\right)}(u_g)=\boldsymbol{r}_{ib}^{\left(G_{ib}\right)}(u_g)=\boldsymbol{r}_{ob}^{\left(G_{ob}\right)}(u_g)=\left[\begin{array}{c}
0\\
u_g\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}


<span id='_Hlk40386595'></span>The coordinate transformations from ''S<sub>ib</sub>'' and ''S<sub>ob</sub>'' to ''S<sub>b</sub>'' are represented by the following matrixes

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq5'></span> <math display="inline">\boldsymbol{M}_{b,ib}^{\left(G_{ib}\right)}=\left[\begin{array}{cccc}
cos{\alpha }_n & sin{\alpha }_n & 0 & \frac{\pi m}{4}\\
-sin{\alpha }_n & cos{\alpha }_n & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq6'></span> <math display="inline">\boldsymbol{M}_{b,ob}^{\left(G_{ob}\right)}=\left[\begin{array}{cccc}
cos{\alpha }_n & -sin{\alpha }_n & 0 & -\frac{\pi m}{4}\\
sin{\alpha }_n & cos{\alpha }_n & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}


<span id='_Hlk40387120'></span>Then, the inner and outer blade profiles are expressed in coordinate system ''S<sub>b</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq7'></span> <math display="inline">\boldsymbol{r}_b^{\left(G_{ib}\right)}(u_g)=\boldsymbol{M}_{b,ib}^{\left(G_{ib}\right)}\boldsymbol{r}_{}^{\left(G\right)}(u_g)=\left[\begin{array}{c}
u_gsin{\alpha }_n+\frac{\pi m}{4}\\
u_gcos{\alpha }_n\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq8'></span> <math display="inline">\boldsymbol{r}_b^{\left(G_{ob}\right)}(u_g)=\boldsymbol{M}_{b,ob}^{\left(G_{ob}\right)}\boldsymbol{r}_{}^{\left(G\right)}(u_g)=\left[\begin{array}{c}
-u_gsin{\alpha }_n-\frac{\pi m}{4}\\
u_gcos{\alpha }_n\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}


<span id='_Hlk41337238'></span><span id='_Hlk40387366'></span><span id='_Hlk41336060'></span>The virtual face-milling cutter surface is represented in coordinate system ''S<sub>c</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq9'></span><math>\boldsymbol{r}_c^{\left(G\right)}({\theta }_g,u_g)=\boldsymbol{M}_{c,b}({\theta }_g)\boldsymbol{r}_b^{\left(G\right)}(u_g)=\left[\begin{array}{c}
(\pm \mbox{ }\mbox{ }u_gsin{\alpha }_n-R\, \pm \, \frac{\pi m}{4})cos{\theta }_g\\
u_gcos{\alpha }_n\\
(\pm \mbox{ }u_gsin{\alpha }_n-\, R\mbox{ }\pm \frac{\pi m}{4})sin{\theta }_g\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}


<span id='_Hlk40630336'></span>From now on, if both upper and lower signs appear in matrixes. Then, the upper sign represents the inner blade or convex side, while the lower sign represents the outer blade or concave side.

<span id='_Hlk40387216'></span>Regarding the fillet parts of the cutter, the position vectors of a point on the circular arc profiles can be expressed in coordinate systems ''S<sub>fib</sub>'' and ''S<sub>fob</sub>'', as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq10'></span> <math display="inline">\boldsymbol{r}_{fib}^{\left(G_{ib}\right)}({\lambda }_g)=\left[\begin{array}{c}
\rho mcos{\lambda }_g\\
-\rho msin{\lambda }_g\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq11'></span> <math>\boldsymbol{r}_{fob}^{\left(G_{ob}\right)}({\lambda }_g)=\left[\begin{array}{c}
-\rho mcos{\lambda }_g\\
-\rho msin{\lambda }_g\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}


<span id='_Hlk40630456'></span>The coordinate transformations from ''S<sub>fib</sub>'' and ''S<sub>fob</sub>'' to ''S<sub>b</sub>'' are represented by the following matrixes

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq12'></span> <math display="inline">\boldsymbol{M}_{b,fib}^{\left(G_{ib}\right)}=\left[\begin{array}{cccc}
1 & 0 & 0 & \frac{\pi m}{4}-tan{\alpha }_n(\rho msin{\alpha }_n+bm-\rho m)-\rho mcos{\alpha }_n\\
0 & 1 & 0 & -(bm-\rho m)\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq13'></span> <math display="inline">\boldsymbol{M}_{b,fob}^{\left(G_{ob}\right)}=\left[\begin{array}{cccc}
1 & 0 & 0 & -(\frac{\pi m}{4}-tan{\alpha }_n(\rho msin{\alpha }_n+bm-\rho m)-\rho mcos{\alpha }_n)\\
0 & 1 & 0 & -(bm-\rho m)\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}


<span id='_Hlk40630502'></span>Then, the inner and outer circular arc profiles are expressed in coordinate system ''S<sub>b</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq14'></span> <math>\begin{array}{c}
\boldsymbol{r}_b^{\left(G_{ib}\right)}({\lambda }_g)=\boldsymbol{M}_{b,fib}^{\left(G_{ib}\right)}\boldsymbol{r}_{fib}^{\left(G_{ib}\right)}({\lambda }_g)
\mbox{=}\left[\begin{array}{c}\mbox{ }\frac{\pi m}{4}-\mbox{tan}{\alpha }_n\mbox{(}bm-\rho m\mbox{ + }\rho m\mbox{sin}{\alpha }_n\mbox{)}-\rho m(\mbox{cos}{\alpha }_n-\mbox{cos}{\lambda }_g)\\
\rho m-bm-\rho m\mbox{sin}{\lambda }_g\\
0\\
1
\end{array}\right]
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq15'></span> <math display="inline">\begin{aligned}
\mathbf{r}_b^{\left(G_{o b}\right)}\left(\lambda_g\right) & =\mathbf{M}_{b, f o b}^{\left(G_{o b}\right)} \mathbf{r}_{f o b}^{\left(G_{o b}\right)}\left(\lambda_g\right) \\
& =\left[\begin{array}{c}
-\frac{\pi m}{4}+\tan \alpha_n\left(b m-\rho m+\rho m \sin \alpha_n\right)+\rho m\left(\cos \alpha_n-\cos \lambda_g\right) \\
\rho m-b m-\rho m \sin \lambda_g \\
0 \\
1
\end{array}\right]
\end{aligned}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}


<span id='_Hlk41337192'></span>The virtual fillet parts of the face-milling cutter are represented in coordinate system ''S<sub>c</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq16'></span> <math>\begin{aligned}
\mathbf{r}_c^{(G)}\left(\lambda_g, \theta_g\right) & =\mathbf{M}_{c, b}\left(\theta_g\right) \mathbf{r}_b^{(G)}\left(\lambda_g\right) \\
& =\left[\begin{array}{c}
\left( \pm \rho m \cos \lambda_g \mp \tan \alpha_n\left(b m-\rho m+\rho m \sin \alpha_n\right) \pm \frac{\pi m}{4} \mp \rho m \cos \alpha_n-R\right) \cos \theta_g \\
\rho m-b m-\rho m \sin \lambda_g \\
\left( \pm \rho m \cos \lambda_g \mp\left(\tan \alpha_n\left(b m-\rho m+\rho m \sin \alpha_n\right) \pm \frac{\pi m}{4} \mp \rho m \cos \alpha_n-R\right) \sin \theta_g\right. \\
1
\end{array}\right]
\end{aligned}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}


<span id='_Hlk40630886'></span>According to the proposed methodology of generation of the gear tooth surfaces in section 2.1, the parametric representations of the family of positions of the cutter surfaces in coordinate system ''S<sub>1</sub>'', are given by Eqs. (17) and (18).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq17'></span> <math>\begin{aligned}
\mathbf{r}_1^{(G)}\left(u_g, \theta_g, \psi_g\right) & =\mathbf{M}_{p, c}\left(\psi_g\right) \mathbf{r}_c^{(G)}\left(\theta_g, u_g\right) \\
& =\left[\begin{array}{c}
u_g\left(\cos \alpha_n \sin \psi_g \pm \cos \psi_g \sin \alpha_n \cos \theta_g\right)+R \cos \psi_g \\
+r_g\left(\sin \psi_g-\psi_g \cos \psi_g\right)-\left(R \pm \frac{\pi m}{4}\right) \cos \psi_g \cos \theta_g \\
u_g\left(\cos \alpha_n \cos \psi_g \mp \sin \alpha_n \cos \theta_g \sin \psi_g\right)-R \sin \psi_g \\
+r_g\left(\cos \psi_g+\psi_g \sin \psi_g\right)+\left(R \mp \frac{\pi m}{4}\right) \sin \psi_g \cos \theta_g \\
\left( \pm \frac{\pi m}{4}-R \pm u_g \sin \alpha_n\right) \sin \theta_g \\
1
\end{array}\right]
\end{aligned}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq18'></span> <math>\left.\begin{array}{rl}
\mathbf{r}_1^{(G)}\left(\lambda_g, \theta_g, \psi_g\right)= & \mathbf{M}_{p, c}\left(\psi_g\right) \mathbf{r}_c^{(G)}\left(\lambda_g, \theta_g\right) \\
= & {\left[\begin{array}{c}
\left(r_g-b m+\rho m-\rho m \sin \lambda_g\right) \sin \psi_g \\
-\left(r_g \psi_g-R\right) \cos \psi_g-\left(R \mp \rho m \cos \lambda_g \pm S\right) \cos \psi_g \cos \theta_g \\
\left(r_g-b m+\rho m-\rho m \sin \lambda_g\right) \cos \psi_g \cdots \\
+\left(r_g \psi_g-R\right) \sin \psi_g+\left(R \mp \rho m \cos \lambda_g \pm S\right) \sin \psi_g \cos \theta_g \\
\left( \pm \rho m \cos \lambda_g \mp S-R\right) \sin \theta_g \\
1
\end{array}\right]} \\
S= & \tan \alpha_n\left(b m-\rho m+\rho m \sin \alpha_n\right)-\frac{\pi m}{4}+\rho m \cos \alpha_n
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}


<span id='_Hlk40631882'></span><span id='_Hlk41316617'></span>According to the theory of gearing in [12], the equations of meshing may be derived as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq19'></span> <math display="inline">f(u_g,{\theta }_g,{\psi }_g)=(\frac{\partial \boldsymbol{r}_1^{\left(G\right)}}{\partial {\theta }_g}\times \frac{\partial \boldsymbol{r}_1^{\left(G\right)}}{\partial u_g})\cdot \frac{\partial \boldsymbol{r}_1^{\left(G\right)}}{\partial {\psi }_g}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq20'></span> <math display="inline">f({\lambda }_g,{\theta }_g,{\psi }_g)=(\frac{\partial \boldsymbol{r}_1^{\left(G\right)}}{\partial {\lambda }_g}\times \frac{\partial \boldsymbol{r}_1^{\left(G\right)}}{\partial {\theta }_g})\cdot \frac{\partial \boldsymbol{r}_1^{\left(G\right)}}{\partial {\psi }_g}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}


Therefore, the generated working and fillet surfaces can be obtained as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq21'></span> <math>\lbrace \begin{array}{c}
\boldsymbol{r}_1^{\left(G\right)}(u_g,{\theta }_g,{\psi }_g)=\boldsymbol{M}_{p,c}({\psi }_g)\boldsymbol{r}_c^{\left(G\right)}({\theta }_g,u_g)\\
f(u_g,{\theta }_g,{\psi }_g)=\mbox{0}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq22'></span> <math>\lbrace \begin{array}{c}
\boldsymbol{r}_1^{\left(G\right)}({\lambda }_g,{\theta }_g,{\psi }_g)=\boldsymbol{M}_{p,c}({\psi }_g)\boldsymbol{r}_c^{\left(G\right)}({\lambda }_g,{\theta }_g)\\
f({\lambda }_g,{\theta }_g,{\psi }_g)=\mbox{0}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}


<span id='_Hlk40787685'></span>In Eq. (21), ''u<sub>g</sub>'' and ''θ<sub>g</sub>'' are the parameters of the working tooth surfaces. In Eq. (22), ''λ<sub>g</sub>'' and ''θ<sub>g</sub>'' are the parameters of the fillet surfaces, and ''ψ<sub>g</sub>'' is the generalized parameter of movement. Arbitrary point coordinates on the wheel tooth surfaces will be solved by numerical method.

<span id='sec2_3'></span>''2.3. Derivation of mathematical model of pinion tooth surfaces''

<span id='_Hlk40441175'></span><span id='_Hlk40471207'></span><span id='_Hlk40471968'></span><span id='_Hlk40883215'></span><span id='_Hlk40645047'></span>The generation of a mathematical model of the pinion tooth surfaces is similar to that of the wheel tooth surfaces. They distinguish only in different cutter blade profiles. Here, the face-milling cutter with the parabolic blade profiles is used to generate the pinion tooth surfaces. Fig. 2 shows the parabolic profiles of the face-milling cutter. The parabolic profiles of the cutter retain tangent to the standard blade with straight-line profiles. The point of tangency is located at arbitrary positions on the straight-line profiles. The distance between the point of tangency and the cutter pitch plane is Δ''y<sub>tp</sub>'' along the vertical direction. The distance between the center points of circular arc profiles and the origin of the coordinate system ''S<sub>b</sub>'', is ''x<sub>fc</sub>'' alongthe horizontal direction. The variable ''u<sub>p</sub>'' represents the distance between arbitrary points on the parabolic profile and the origin of the coordinate system ''S<sub>ib</sub>'' or ''S<sub>ob</sub>''. The variable ''λ<sub>p</sub>'' represents the polar angle of arbitrary points on the circular arc profiles.

The working and fillet parts of the pinion tooth surfaces will be derived, respectively.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image24.png|234px]] </div>

<div id="fig2" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 2'''. The parabolic profiles of cutter for generation of the pinion tooth surfaces</div>

The position vectors of a point on the parabolic profiles can be expressed in the coordinate systems ''S<sub>ib</sub>'' and ''S<sub>ob</sub>'', as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq23'></span> <math display="inline">\boldsymbol{r}_{ib}^{\left(P_{ib}\right)}(u_p)=\left[\begin{array}{c}
a_{pf}u_p{}^2\\
u_p\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq24'></span> <math display="inline">\boldsymbol{r}_{ob}^{\left(P_{ob}\right)}(u_p)=\left[\begin{array}{c}
-a_{pf}u_p{}^2\\
u_p\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}


The coordinate transformations from ''S<sub>ib</sub>'' and ''S<sub>ob</sub>'' to ''S<sub>b</sub>'' are represented by the following matrixes

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq25'></span> <math display="inline">\boldsymbol{M}_{b,ib}^{\left(P_{ib}\right)}=\left[\begin{array}{cccc}
cos{\alpha }_n & sin{\alpha }_n & 0 & \frac{\pi m}{4}+\Delta y_{tp}tan{\alpha }_n\\
-sin{\alpha }_n & cos{\alpha }_n & 0 & \Delta y_{tp}\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq26'></span> <math display="inline">\boldsymbol{M}_{b,ob}^{\left(P_{ob}\right)}=\left[\begin{array}{cccc}
cos{\alpha }_n & -sin{\alpha }_n & 0 & -(\frac{\pi m}{4}+\Delta y_{tp}tan{\alpha }_n)\\
sin{\alpha }_n & cos{\alpha }_n & 0 & \Delta y_{tp}\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}


Then, the inner and outer blade profiles are expressed in the coordinate system ''S<sub>b</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq27'></span> <math display="inline">\mathbf{r}_b^{\left(P_{i b}\right)}\left(u_p\right)=\mathbf{M}_{b, i b}^{\left(P_{i b}\right)} \mathbf{r}_{i b}^{\left(P_{i b}\right)}\left(u_p\right)=\left[\begin{array}{c}
a_{p f} u_p{ }^2 \cos \alpha_n+u_p \sin \alpha_n+\left(\frac{\pi m}{4}+\Delta y_{t p} \tan \alpha_n\right) \\
-a_{p f} u_p{ }^2 \sin \alpha_n+u_p \cos \alpha_n+\Delta y_{t p} \\
0 \\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq28'></span> <math display="inline">\mathbf{r}_b^{\left(P_{o b}\right)}\left(u_p\right)=\mathbf{M}_{b, o b}^{\left(P_{o b}\right)} \mathbf{r}_{o b}^{\left(P_{o b}\right)}\left(u_p\right)=\left[\begin{array}{c}
-a_{p f} u_p{ }^2 \cos \alpha_n-u_p \sin \alpha_n-\left(\frac{\pi m}{4}+\Delta y_{t p} \tan \alpha_n\right) \\
-a_{p f} u_p{ }^2 \sin \alpha_n+u_p \cos \alpha_n+\Delta y_{t p} \\
0 \\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}


The virtual cutter surfaces are represented in the coordinate system ''S<sub>c</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq29'></span> <math display="inline">\begin{aligned}
\mathbf{r}_c^{(P)}\left(\theta_p, u_p\right)= & \mathbf{M}_{c, b}\left(\theta_p\right) \mathbf{r}_b^{(P)}\left(u_p\right) \\
& =\left[\begin{array}{c}
\left( \pm a_{p f} u_p{ }^2 \cos \alpha_n \pm u_p \sin \alpha_n \pm \frac{\pi m}{4} \pm \Delta y_{t p} \tan \alpha_n\right) \cos \theta_p-R \sin \theta_p \\
-a_{p f} u_p{ }^2 \sin \alpha_n+u_p \cos \alpha_n+\Delta y_{t p} \\
\left( \pm a_{p f} u_p{ }^2 \cos \alpha_n \pm u_p \sin \alpha_n \pm \frac{\pi m}{4} \pm \Delta y_{t p} \tan \alpha_n\right) \sin \theta_p-R \sin \theta_p \\
1
\end{array}\right]
\end{aligned}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}


Regarding the fillet parts of the cutter, the position vectors of a point on the circular arc profiles can be expressed in the coordinate systems ''S<sub>fib</sub>'' and ''S<sub>fob</sub>'', as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq30'></span> <math display="inline">\boldsymbol{r}_{fib}^{\left(P_{ib}\right)}({\lambda }_p)=\left[\begin{array}{c}
\rho mcos{\lambda }_p\\
-\rho msin{\lambda }_p\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq31'></span> <math>\boldsymbol{r}_{fob}^{\left(P_{ob}\right)}({\lambda }_p)=\left[\begin{array}{c}
-\rho mcos{\lambda }_p\\
-\rho msin{\lambda }_p\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}


The coordinate transformations from ''S<sub>fib</sub>'' and ''S<sub>fob</sub>'' to ''S<sub>b</sub>'' are represented by the following matrixes

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq32'></span> <math>\boldsymbol{M}_{b,fib}^{\left(P_{ib}\right)}(x_{fc})=\left[\begin{array}{cccc}
1 & 0 & 0 & x_{fc}\\
0 & 1 & 0 & -(bm-\rho m)\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq33'></span> <math>\boldsymbol{M}_{b,fob}^{\left(P_{ob}\right)}(x_{fc})=\left[\begin{array}{cccc}
1 & 0 & 0 & -x_{fc}\\
0 & 1 & 0 & -(bm-\rho m)\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}


Then, the inner and outer circular arc profiles are expressed in the coordinate system ''S<sub>b</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq34'></span> <math>\boldsymbol{r}_c^{\left(P\right)}({\lambda }_p,{\theta }_p,x_{fc})=\boldsymbol{M}_{b,fib}^{\left(P_{ib}\right)}(x_{fc})\boldsymbol{r}_{fib}^{\left(P_{ib}\right)}({\lambda }_p)=\left[\begin{array}{c}
x_{fc}\mbox{+}\rho m\mbox{cos}{\lambda }_p\\
\mbox{ }\rho m-bm-\rho m\mbox{sin}{\lambda }_p\\
0\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq35'></span> <math display="inline">\mathbf{r}_b^{\left(P_{o b}\right)}\left(\lambda_p, x_{f c}\right)=\mathbf{M}_{b, f o b}^{\left(P_{o b}\right)}\left(x_{f c}\right) \mathbf{r}_{f o b}^{\left(P_{o b}\right)}\left(\lambda_p\right)=\left[\begin{array}{c}
-\left(x_{f c}+\rho m \cos \lambda_p\right) \\
\rho m-b m-\rho m \sin \lambda_p \\
0 \\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}


The virtual cutter fillet surfaces are represented in the coordinate system ''S<sub>c</sub>'' as

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq36'></span> <math display="inline">\begin{aligned}
\mathbf{r}_c^{(P)}\left(\theta_p, u_p\right)= & \mathbf{M}_{c, b}\left(\theta_p\right) \mathbf{r}_b^{(P)}\left(u_p\right) \\
& =\left[\begin{array}{c}
\left( \pm a_{p f} u_p{ }^2 \cos \alpha_n \pm u_p \sin \alpha_n \pm \frac{\pi m}{4} \pm \Delta y_{t p} \tan \alpha_n\right) \cos \theta_p-R \sin \theta_p \\
-a_{p f} u_p{ }^2 \sin \alpha_n+u_p \cos \alpha_n+\Delta y_{t p} \\
\left( \pm a_{p f} u_p{ }^2 \cos \alpha_n \pm u_p \sin \alpha_n \pm \frac{\pi m}{4} \pm \Delta y_{t p} \tan \alpha_n\right) \sin \theta_p-R \sin \theta_p
\end{array}\right]
\end{aligned}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq37'></span> <math>\boldsymbol{r}_c^{\left(P\right)}({\lambda }_p,{\theta }_p,x_{fc})=\boldsymbol{M}_{c,b}({\theta }_p)\boldsymbol{r}_b^{\left(P\right)}({\lambda }_p,x_{fc})=\left[\begin{array}{c}
(\pm x_{fc}-R\pm \rho m\mbox{cos}{\lambda }_p\mbox{)cos}{\theta }_p\\
\rho m-bm-\rho m\mbox{sin}{\lambda }_p\\
(\pm x_{fc}-R\pm \rho m\mbox{cos}{\lambda }_p\mbox{)sin}{\theta }_p\\
1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}


<span id='_Hlk38221162'></span>According to the proposed methodology of generation of the gear tooth surfaces in section 2.1, the parametric representations of the family of positions of the cutter blade surfaces in coordinate system ''S<sub>1</sub>'', is given by Eqs. (38) and (39).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq38'></span> <math>\begin{aligned}
\mathbf{r}_1^{(P)}\left(u_p, \theta_p, \psi_p\right)= & \mathbf{M}_{p, c}\left(\psi_p\right) \mathbf{r}_c^{(P)}\left(\theta_p, u_p\right) \\
& {\left[\begin{array}{l}
-a_{p f} u_p^2\left(\sin \alpha_n \sin \psi_p \mp \cos \alpha_n \cos \psi_p \cos \theta_p\right) \cdots \\
+u_p\left(\cos \alpha_n \sin \psi_p \pm \cos \psi_p \sin \alpha_n \cos \theta_p\right)+R \cos \psi_p+\Delta y_{t p} \sin \psi_p \\
+r_p\left(\sin \psi_p-\psi_p \cos \psi_p\right) \pm \cos \theta_p \cos \psi_p\left(\frac{\pi m}{4}+\Delta y_{t p} \tan \alpha_n-R\right) \\
-a_{p f} u_p^2\left(\cos \psi_p \sin \alpha_n \pm \cos \alpha_n \cos \theta_p \sin \psi_p\right) \cdots \\
+u_p\left(\cos \alpha_n \cos \psi_p \mp \sin \alpha_n \cos \theta_p \sin \psi_p\right)-R \sin \psi_p+\Delta y_{t p} \cos \psi_p \cdots \\
+r_p\left(\cos \psi_p+\psi_p \sin \psi_p\right) \mp \cos \theta_p \sin \psi_p\left(\frac{\pi m}{4}+\Delta y_{t p} \tan \alpha_n+R\right) \\
\left( \pm a_{p f} u_p^2 \cos \alpha_n \pm u_p \sin \alpha_n \pm \frac{\pi m}{4}+\Delta y_{t p} \tan \alpha_n-R\right) \sin \theta_p \\
1
\end{array}\right] }
\end{aligned}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq39'></span> <math>\left.\begin{array}{rl}
\mathbf{r}_1^{(P)}\left(\lambda_p, \theta_p, \psi_p, x_{f c}\right)= & \mathbf{M}_{p, c}\left(\psi_p\right) \mathbf{r}_c^{(P)}\left(\lambda_p, \theta_p, x_{f c}\right) \\
& =\left[\begin{array}{c}
\left(r_p-b m+\rho m-\rho m \sin \lambda_p\right) \sin \psi_p-S\left(\lambda_p, \psi_p\right) \cos \psi_p \\
\left(r_p-b m+\rho m-\rho m \sin \lambda_p\right) \cos \psi_p+S\left(\lambda_p, \psi_p\right) \sin \psi_p \\
-S\left(\lambda_p\right) \sin \theta_p \\
1
\end{array}\right] \\
&S\left(\lambda_p, \psi_p\right)=r_p \psi_p-R+\left(R \mp x_{f c} \mp \rho m \cos \lambda_p\right) \cos \theta_p
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}


The equations of meshing for the pinion tooth surfaces may be derived as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq40'></span> <math display="inline">f({\theta }_p,u_p,{\psi }_p)=(\frac{\partial \boldsymbol{r}_1^{\left(P\right)}}{\partial {\theta }_p}\times \frac{\partial \boldsymbol{r}_1^{\left(P\right)}}{\partial u_p})\cdot \frac{\partial \boldsymbol{r}_1^{\left(P\right)}}{\partial {\psi }_p}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq41'></span> <math display="inline">f({\lambda }_p,{\theta }_p,{\psi }_p,x_{fc})=(\frac{\partial \boldsymbol{r}_1^{\left(P\right)}}{\partial {\lambda }_p}\times \frac{\partial \boldsymbol{r}_1^{\left(P\right)}}{\partial {\theta }_p})\cdot \frac{\partial \boldsymbol{r}_1^{\left(P\right)}}{\partial {\psi }_p}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}


Therefore, the generated working and fillet surfaces can be obtained as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq42'></span> <math>\lbrace \begin{array}{c}
\boldsymbol{r}_1^{\left(P\right)}(u_p,{\theta }_p,{\psi }_p)=\boldsymbol{M}_{p,c}({\psi }_p)\boldsymbol{r}_c^{\left(P\right)}({\theta }_p,u_p)\\
f(u_p,{\theta }_p,{\psi }_p)=\mbox{0}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq43'></span> <math>\lbrace \begin{array}{c}
\boldsymbol{r}_1^{\left(P\right)}({\lambda }_p,{\theta }_p,{\psi }_p,x_{fc})=\boldsymbol{M}_{p,c}({\psi }_p)\boldsymbol{r}_c^{\left(P\right)}({\lambda }_p,{\theta }_p,x_{fc})\\
f({\lambda }_g,{\theta }_g,{\psi }_g,x_{fc})=\mbox{0}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}


In Eq. (42), ''u<sub>p</sub>'' and ''θ<sub>p</sub>'' are the parameters of the working tooth surfaces. In Eq. (43), ''λ<sub>p</sub>'' and ''θ<sub>p</sub>'' are the parameters of the fillet surfaces. ''ψ<sub>p</sub>'' is the generalized parameter of movement.

<span id='_Hlk40650762'></span><span id='_Hlk40650626'></span><span id='_Hlk40782851'></span>Regarding the determination of parameter ''x<sub>fc</sub>''. Due to the edge radius of the cutter blades and parabola coefficient being constant. When the parabolic profiles tangent to the circular arc profiles, unit normal vectors of the parabolic profiles are equal to unit normal vectors of the circular arc profiles at the point of tangency. Simultaneously, the position vectors of the parabolic profiles are equal to the position vectors of the circular arc profiles at the point of tangency. Therefore, the parameter ''x<sub>fc</sub>'' can be determined by Eq. (44).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq44'></span> <math display="inline">\lbrace \begin{array}{c}
\boldsymbol{r}_b^{\left(P\right)}(u_p)=\boldsymbol{r}_b^{\left(P\right)}({\lambda }_p,x_{fc})\\
\boldsymbol{n}_b^{\left(P\right)}(u_p)=\boldsymbol{n}_b^{\left(P\right)}({\lambda }_p,x_{fc})
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}


Where,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq45'></span> <math display="inline">\boldsymbol{n}_b^{\left(P\right)}(u_p)=\frac{\partial \boldsymbol{r}_b^{\left(P\right)}(u_p)/\partial u_p}{\vert \partial \boldsymbol{r}_b^{\left(P\right)}(u_p)/\partial u_p\vert }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq46'></span> <math>\boldsymbol{n}_b^{\left(P\right)}({\lambda }_p,x_{fc})=\frac{\partial \boldsymbol{r}_b^{\left(P\right)}({\lambda }_p,x_{fc})/\partial {\lambda }_p}{\vert \partial \boldsymbol{r}_b^{\left(P\right)}({\lambda }_p,x_{fc})/\partial {\lambda }_p\vert }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}


<span id='_Hlk41337680'></span>In Eq. (44), there are three independent scalar equations with three unknown parameters (''u<sub>p</sub>'', ''λ<sub>p</sub>'' and ''x<sub>fc</sub>''). Therefore, the three unknown parameters are easy to solve.

=3. Finite element model for stress analysis=

In order to create a precise finite element model (FEM), point clouds on the tooth surfaces are employed to generate nodes and hexahedral elements of FEM. Here, the point cloud coordinates come from the tooth surface equations '''r'''(''P'') 1and '''r'''(''G'') 1. This method of establishing nodes and elements of a gear set comes from Litvin et al [12]. Accurately defining the initial contact position and boundary conditions of the gear set is a key step for establishing a FEM of gear sets.

''3.1. Initial contact position of gear sets''

The definition of initial contact position is essential for finite element analysis of the gear sets. It takes the position where the gear sets with assembly errors just contact as the solution target. That is, when the tooth surfaces of the gear sets are just in contact, the initial rotation angle of the two matched gears is calculated.

<span id='_Hlk40808455'></span>The assumption, the pinion is fixed to coordinate system ''S<sub>1</sub>'', and it rotates around the axis ''z<sub>f</sub>'' of coordinate system ''S<sub>f</sub>'' which is rigidly fixed to the gear transmission frame. Moreover, the wheel is fixed to the coordinate system ''S<sub>2</sub>'' which is separated with respect to the coordinate system ''S<sub>f</sub>'' by a distance ''C''+Δ''C'' along the axis ''y<sub>f</sub>'' direction. The relative position of the wheel with respect to the pinion has four assembly errors: axial displacement of the wheel with respect to the pinion (Δ''A''), the center distance error (Δ''C''), intersecting shaft angle error (Δ''V''), crossing shaft angle error (Δ''H''). (see Fig. 3)

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image49.png|240px]] </div>

<span id='fig3'></span><div id="_Hlk40960466" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 3'''. Definition of coordinate systems under gear assembly errors</div>

<span id='_Hlk40882869'></span><span id='_Hlk40869782'></span>To solve the initial contact points on the gear tooth surfaces when the wheel is in meshing with the pinion. Both the wheel and pinion tooth surfaces have to be expressed in the same coordinate system ''S<sub>f</sub>'' . The pinion tooth surfaces are represented in the coordinate system ''S<sub>f</sub>'' as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq47'></span> <math>\boldsymbol{r}_f^{\left(P\right)}(u_p,{\theta }_p,{\phi }_p)=\boldsymbol{M}_{f,1}({\phi }_p)\boldsymbol{r}_1^{\left(P\right)}(u_p,{\theta }_p,{\psi }_p(u_p,{\theta }_p))</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}


Where,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq48'></span> <math display="inline">\boldsymbol{M}_{f,1}({\phi }_p)=\left[\begin{array}{cccc}
cos{\phi }_p & -sin{\phi }_p & 0 & 0\\
-sin{\phi }_p & cos{\phi }_p & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}


<span id='_Hlk40872497'></span>In Eq. (47), the function ''ψ<sub>p</sub>'' (''u<sub>p</sub>'', ''θ<sub>p</sub>'') represents that the parameter ''ψ<sub>p</sub>'' in Eq.(40) is expressed by the parameter ''u<sub>p</sub>'' and ''θ<sub>p</sub>''.

Introducing four assembly errors, the wheel tooth surfaces are represented in the coordinate system ''S<sub>f</sub>'' as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq49'></span> <math display="inline">\boldsymbol{r}_f^{\left(g\right)}(u_g,{\theta }_g,{\phi }_g)\mbox{=}\boldsymbol{M}_{f,k}\boldsymbol{M}_{k,l}\boldsymbol{M}_{l,m}\boldsymbol{M}_{m,n}\boldsymbol{M}_{n,2}({\phi }_g)\boldsymbol{r}_1^{\left(g\right)}(u_g,{\theta }_g,{\psi }_g(u_g,{\theta }_g))</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}


Where,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq50'></span> <math display="inline">\boldsymbol{M}_{f,k}=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & C+\Delta C\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq51'></span> <math display="inline">\boldsymbol{M}_{k,l}=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 1 & \Delta A\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq52'></span> <math display="inline">\boldsymbol{M}_{l,m}=\left[\begin{array}{cccc}
cos(\Delta H) & 0 & sin(\Delta H) & 0\\
0 & 1 & 0 & 0\\
-sin(\Delta H) & 0 & cos(\Delta H) & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq53'></span> <math display="inline">\boldsymbol{M}_{m,n}=\left[\begin{array}{cccc}
1 & 0 & 0 & 0\\
0 & cos(\Delta V) & -sin(\Delta V) & 0\\
0 & sin(\Delta V) & cos(\Delta V) & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq54'></span> <math display="inline">\boldsymbol{M}_{n,2}({\phi }_g)=\left[\begin{array}{cccc}
cos{\phi }_g & -sin{\phi }_g & 0 & 0\\
sin{\phi }_g & cos{\phi }_g & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}


In Eq. (49), the function ''ψ<sub>g</sub>'' (''u<sub>g</sub>'', ''θ<sub>g</sub>'') represents that the parameter ''ψ<sub>g</sub>'' in Eq.(19) is expressed by the parameter ''u<sub>g</sub>'' and ''θ<sub>g</sub>''.

When the wheel is in meshing with the pinion, Both the wheel and pinion tooth surfaces must be in continuous tangency, their position vectors and unit normal vectors coincide at arbitrary points of tangency. Thus,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <span id='eq55'></span> <math>\lbrace \begin{array}{c}
\boldsymbol{r}_f^{\left(P\right)}(u_p,{\theta }_p,{\phi }_p)=\boldsymbol{r}_f^{\left(g\right)}(u_g,{\theta }_g,{\phi }_g)\\
\boldsymbol{n}_f^{\left(P\right)}(u_p,{\theta }_p,{\phi }_p)=\boldsymbol{n}_f^{\left(g\right)}(u_g,{\theta }_g,{\phi }_g)
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}


In Eq.(55), there are five independent scalar equations with six unknowns. If the initial rotation angle of the pinion ''φ<sub>p</sub>'' is given, the position of the initial contact point on the gear tooth surfaces, and the initial rotation angle of the wheel ''φ<sub>g</sub>'' will be solved by numerical method. The initial rotation angle (''φ<sub>p</sub>'','' φ<sub>g</sub>'') and four assembly errors (Δ''A'', Δ''C'', Δ''V'', Δ''H'') will be used for the assembly of the gear set in the finite element model.

''3.2. Definition of boundary conditions''

Finite element method is one of the basic approaches of stress analysis, and it can be used for the simulation of contact and bending stresses of gear drives [24-26]. Boundary conditions are a significant part of FEM. The following is a detailed creation method of FEM of the gear set.

<span id='_Hlk40975243'></span><span id='_Hlk40977079'></span><span id='_Hlk40977095'></span><span id='OLE_LINK3'></span><span id='OLE_LINK4'></span><span id='_Hlk40992386'></span>Fig. 4 shows the schematic diagram of the three-tooth finite element model for stress analysis. The number of nodes and elements in the three-tooth model is simplified and does not represent the actual number of nodes and elements. Generally, to ensure enough boundary conditions for precise stress analysis, the three-tooth models are employed for finite element analysis of the gear drives. The number of elements of each gear model is 73500 with 88995 nodes. The element type is C3D8I, the hexahedral first-order three-dimensional solid element. The material of the gear set is steel with specific material properties (see Table 1). In Fig. 4(c), the gear ring surfaces and two gear rim surfaces are rigidly constrained to the reference points on the gear axes. All degrees of freedom of the reference point (RP-1) on the pinion axis are constrained, but the rotation degree of freedom of the reference point (RP-2) on the wheel axis is only released. In addition, a resisting torque is applied to the ''z<sub>2</sub>'' axis of the wheel. Account of no remarkable advantage has been found when the concave or convex sides of the gear tooth surface as driving surfaces [20]. In this paper, the concave sides of the pinion tooth surfaces are defined as master surfaces (Fig. 4(a)), while the convex sides of the wheel tooth surfaces are defined as slave surfaces (Fig. 4(b)). When the stress is solved for a contact position, the pinion rotates ''φ<sub>pi</sub>'' (''φ<sub>pi</sub> ''= (2π/''N''<sub>1</sub>)/21) around the ''z<sub>1</sub>'' axis, and the stress solution is performed for the next con tact position. Until stress analysis of 21 contact positions along a meshing cycles is entirely completed. Finally, collecting the contact and bending stresses at all contact positions.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image59.png|336px]] </div>

<span id='fig4'></span><div id="_Hlk41058038" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 4'''. Definition of three-tooth finite element model</div>

=4. Numerical examples=

<span id='_Hlk41052970'></span><span id='_Hlk41060064'></span><span id='_Hlk41047663'></span>In this paper, four group examples are employed for investigation of contact stresses, bending stresses of the proposed gear drives, and all examples of curvilinear cylindrical gear drives remain free of undercutting and interference. The basic design data of the curvilinear cylindrical gear drives are shown in Table 1 for all examples.

<span id='Table1'></span><span id='_Hlk41052948'></span>'''Table 1''' Common basic design data of the curvilinear cylindrical gear drives

{| style="width: 39%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Parameter
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Values
|-
|  style="text-align: center;"|Number of teeth, ''N''
|  style="text-align: center;"|31/40
|-
|  style="text-align: center;"|Module,'' m'' (mm)
|  style="text-align: center;"|4
|-
|  style="text-align: center;"|Normal pressure angle, ''α<sub>n</sub>'' (deg)
|  style="text-align: center;"|20
|-
|  style="text-align: center;"|Addendum coefficient,'' a''
|  style="text-align: center;"|1
|-
|  style="text-align: center;"|Dedendum coefficient,'' b''
|  style="text-align: center;"|1.25
|-
|  style="text-align: center;"|Face width,'' W'' (mm)
|  style="text-align: center;"|30
|-
|  style="text-align: center;"|Edge radius coefficient, ''ρ''
|  style="text-align: center;"|0.38
|-
|  style="text-align: center;"|Young’s Modulus,'' E'' (MPa)
|  style="text-align: center;"|210000
|-
|  style="text-align: center;"|Poisson’s ratio, ''ν''
|  style="text-align: center;"|0.3
|-
|  style="border-bottom: 1pt solid black;text-align: center;"|Resisting torque, ''T'' (Nm)
|  style="border-bottom: 1pt solid black;text-align: center;"|400
|}


<span id='sec5_1'></span>''4.1. Example 1''

<span id='_Hlk41129104'></span><span id='_Hlk41139534'></span><span id='_Hlk41129162'></span>In this group example, four assembly errors are determined as independent experiment variable, respectively. The main purpose is to investigate the influence of assembly errors on contact stresses, bending stresses. Here, the cutter mean pitch radius ''R'' = 40 mm, the parabola coefficient ''a<sub>pf</sub>'' = 0.0015, and the offset distance of parabolic profile ''∆y<sub>tp</sub>'' = 0 mm. The design data of the assembly errors are shown in Table 2. There are four kinds of assembly errors, each assembly error is composed of 4 cases.

<span id='Table2'></span><span id='_Hlk41061998'></span>'''Table 2''' Design data of four kinds of assembly errors

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;"|<span id='_Hlk38053475'></span>Type of assembly errors
|  colspan='4'  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|∆''V'' (°)
|  colspan='5'  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|∆''H'' (°)
|  colspan='5'  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|∆''A'' (mm)
|  colspan='6'  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|∆''C'' (mm)
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;"|Error values
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.05
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.1
|  colspan='2'  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.2
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.05
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.1
|  colspan='2'  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.2
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.1
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.2
|  colspan='2'  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.5
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.1
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.2
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.5
|}


<span id='_Hlk41080114'></span><span id='_Hlk41078428'></span><span id='_Hlk41139911'></span><span id='_Hlk41138911'></span><span id='OLE_LINK5'></span>Fig. 5 shows the evolution of contact stresses for the gear drives designed for example 1. The intersecting shaft (Δ''V'') and crossing shaft angle errors (Δ''H'') have no effect on magnitudes of contact stresses, but the contact position is slightly shifted toward the gear edge with increasing assembly error (Fig. 6(a) and (b)). The red curves represent the centerline of the gear along the tooth width in Fig. 6, and the red curves in stress contours that appear later in this article also represent the centerline of the gear. The axial displacement error (Δ''A'') and center distance errors (Δ''C'') have no effect on contact stresses unless the values of axial displacement error and center distance error are too large. If Δ''A'' and Δ''C ''are too large, severe contact stresses appear on the addendum edge of the gear at contact positions 10 and 20. With the gear set rotates, the pinion rotates exactly one cycle from contact position 10 to 20, so the stresses are almost equal at these two contact positions. Among them, the stress contours of contact position 10 is shown in Fig.6(c) and (d). In Fig.6(c), it can be seen that the offset of the contact patterns increases with the increase of the magnitude of Δ''A''. In Fig.6(d), there is no offset of the contact patterns with the increase of the magnitude of Δ''C''.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image60.jpeg|492px]] </div>

<span id='fig8'></span><span id='fig5'></span><span id='_Hlk41134141'></span><span id='_Hlk41067486'></span><div id="OLE_LINK1" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 5'''. Evolution of contact stresses for gear drives designed for example 1</div>
[[File:Draft_Xuegang_614753900-image61.png|centre|600x600px]]
[[File:Draft_Xuegang_614753900-image62.png|centre|600x600px]]
[[File:Draft_Xuegang_614753900-image63.png|centre|600x600px]]
[[File:Draft_Xuegang_614753900-image64.png|centre|600x600px]]

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 6'''. Stress contours at specific contact positions under four assembly errors</div>

<span id='_Hlk43241972'></span>Fig. 7 shows the evolution of bending stresses for the gear drives designed for example 1. The intersecting shaft error (Δ''V'') and crossing shaft angle error (Δ''H'') have no effect on bending stresses. The axial displacement (Δ''A'') error and center distance error (Δ''C'') have no effect on bending stress unless the values of axial displacement error are too large. If Δ''A ''is too large, the mean bending stresses are increased. To sum up, the suitable assembly error is not sensitive to the contact stress and bending stress of the curvilinear gear sets.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image65.jpeg|492px]] </div>

<span id='fig9'></span><div id="fig6" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 7'''. Evolution of bending stresses for gear drives designed for example 1</div>

''4.2. Example 2''

<span id='_Hlk41145372'></span><span id='_Hlk41159564'></span><span id='_Hlk41135990'></span>In section 4.1, it has been confirmed that assembly errors benefit the reduction of transmission errors for the gear drives with parabolic profile. In this section, assume that the assembly error of the gear drives is defined as ''∆H ''= 0.2°, the cutter mean pitch radius ''R'' = 40 mm, the offset distance of parabolic profile ''∆y<sub>tp</sub>'' = 0 mm. The parabola coefficient is determined as the independent experiment variable. The design data of parabola coefficients are shown in Table 3.

<span id='Table3'></span>'''Table 3''' Design data of different parabola coefficients

{| style="width: 55%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Case
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|1
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|2
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|3
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|4
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|Parabola coefficient, ''a<sub>pf</sub>''(1/mm)
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.0005
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.001
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0.0015
|}


<span id='_Hlk41138762'></span><span id='_Hlk41146502'></span>Fig. 8 shows the evolution of contact stresses for the gear drives designed for example 2. Maximum contact stresses during the meshing cycle occur at contact positions 10 and 20, and the contact stresses decrease with a higher parabola coefficient at contact positions 10 and 20. It can be seen from Fig. 9 that the maximum contact stress at contact position 10 occurs at the tooth tip, which is caused by the insufficient modification of the tooth tip. Therefore, the parabolic modification can eliminate the edge contact at the addendum.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image66.jpeg|246px]] </div>

<div id="_Hlk41140305" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 8'''. Evolution of contact stresses for gear drives designed for example 2</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image67.png|600px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 9'''. Stress contours at contact position 10 for example 2</div>

Fig. 10 shows the evolution of bending stresses for the gear drives designed for example 2. In the meshing cycle, the maximum bending stress magnitude increases with the parabola coefficient increment. In the first half cycle (Contact position 1-11), with the increase of the parabolic coefficient, the bending stress first decreases and then increases; in the second half cycle (Contact position 11-21), it is just the opposite.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image68.jpeg|246px]] </div>

<span id='fig13'></span><div id="_Hlk41140319" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 10'''. Evolution of bending stresses for gear drives designed for example 2</div>

<span id='_Hlk40992703'></span>''4.3. Example 3''

<span id='_Hlk41140543'></span>In this group example, the cutter mean pitch radius is determined as the independent experiment variable. The axial displacement error ''∆A ''= 0.2 mm, the parabola coefficient ''a<sub>pf</sub>'' = 0.0015 1/mm, and the offset distance of parabolic profile ''∆y<sub>tp</sub>'' = 0 mm. The design data of the cutter mean pitch radius are shown in Table 4.

<span id='Table4'></span><span id='_Hlk41140960'></span>'''Table 4''' Design data of the cutter mean pitch radius

{| style="width: 41%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Case
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|1
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|2
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|3
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|4
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Mean radius of cutter, ''R ''(mm)
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|40
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|60
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|80
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|100
|}


<span id='_Hlk41143385'></span><span id='_Hlk41141627'></span><span id='_Hlk41141608'></span>Fig. 11 shows the evolution of contact stresses for the gear drives designed for example 3. In the meshing cycle, the maximum contact stresses at any contact position are reduced with the increment of the cutter mean pitch radius. It can be seen from Fig. 12 that as the cutter radius increases, the width of the contact pattern also increases, at the same time, the increase in the width of the contact pattern leads to a decrease in the contact stress. However, the center point of the contact pattern is more seriously offset under the same assembly error.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image69.jpeg|246px]] </div>

<div id="fig16" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 11'''. Evolution of contact stresses for gear drives designed for example 3</div>
[[File:Draft_Xuegang_614753900-image70.png|centre|600x600px]]

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 12'''. Stress contours at contact position 10 for example 3</div>

<span id='_Hlk41142233'></span><span id='_Hlk41164607'></span>Fig. 13 shows the evolution of bending stresses for the gear drives designed for example 3. In the meshing cycle, maximum bending stresses are less affected by the cutter mean pitch radius. To sum up the above, the results show that a higher cutter mean pitch radius can improve the meshing characteristics of the gear drives.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image71.jpeg|246px]] </div>

<div id="fig17" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 13'''. Evolution of bending stresses for gear drives designed for example 3</div>

''4.4. Example 4''

In this group example, the offset distance of the parabolic profiles is determined as the independent experiment variable under ideal assembly conditions. The cutter mean pitch radius ''R'' = 40 mm, and the parabola coefficient ''a<sub>pf</sub>'' = 0.0005/mm. The design data of the offset distance of the parabolic profiles are shown in Table 5.

<span id='Table5'></span><span id='_Hlk41145289'></span><span id='_Hlk41160164'></span><span id='_Hlk41147159'></span>'''Table 5''' Design data of the offset distance of parabolic profiles

{| style="width: 70%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|<span id='_Hlk38054711'></span>Case
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|1
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|2
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|3
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|4
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|5
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|<span id='_Hlk38059192'></span>Offset distance of parabolic profile, ''∆y<sub>tp</sub> ''(mm)
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|0
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|-1
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|-2
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|-3
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|-4
|}


<span id='_Hlk41147437'></span><span id='_Hlk41149792'></span><span id='_Hlk41147109'></span><span id='_Hlk41148022'></span>Fig. 14 shows the evolution of contact stresses for the gear drives for example 4. In the meshing cycle, the maximum contact stress occurs at contact position 10 and 20, it reduces with higher the offset distance towards tooth root at these contact position. The contact stresses increase with higher the offset distance except contact position 10 and 20. It can be seen from Fig. 15 that the reason for the reduction of the contact stress at the contact position 10 is that the addendum edge contact is gradually eliminated in the process of increasing the offset amount. Therefore, an appropriate offset amount can reduce the maximum contact stress during the meshing cycle.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image72.jpeg|246px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 14'''. Evolution of contact stresses for gear drives designed for example 4</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image73.png|600px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 15'''. Stress contours at contact position 10 for example 4</div>

<span id='_Hlk41148276'></span><span id='_Hlk41148933'></span><span id='_Hlk41164751'></span>Fig. 16 shows the evolution of bending stresses for the gear drives designed for example 4. In the meshing cycle, maximum bending stresses are reduced with the increment of the offset distance of the parabolic profile. Bending stresses are increased with the increment of the offset distance during the first half cycle (Contact position 1-11). In addition, bending stresses are reduced with the increment of the offset distance during the second half cycle (Contact position 11-21). To sum up the above, the obtained results show higher offset distance toward the tooth root will bring about worse mechanical behavior, higher vibration, and noise for the gear drives without edge contact.

<div id="_GoBack" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Xuegang_614753900-image74.jpeg|246px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Fig. 16'''. Evolution of bending stresses for gear drives designed for example 4</div>

=5. Conclusions=

Contact and bending stresses of the curvilinear cylindrical gear drives are investigated by four groups of examples, respectively. Based on the research work, the following conclusions can be drawn:

(1) The assembly error is not sensitive to contact stress and bending stress of gear drives. The curvilinear cylindrical gear drives are more suitable for industrial applications.

(2) Tooth flank modification can eliminate tooth tip-edge contact, but it will increase bending stress at the tooth root.

(3) The offset of the modified parabola can also eliminate the tip-edge contact, which reduces the bending stress at the tooth root.

(4) If the cutter pitch mean radius is increased, then the width of the contact pattern is increased, and contact stresses on gear tooth surfaces is decreased. Therefore, a higher cutter mean pitch radius can improve the mechanical performances of the curvilinear cylindrical gear drives.

==Acknowledgements==

<span id='_Hlk49700416'></span><span id='OLE_LINK6'></span><span id='OLE_LINK7'></span>The authors express their deep gratitude to the financial support of Sichuan Provincial Key Lab of Process Equipment and Control (No. GK202211, No. GK201907), Science & Technology Plan Project of Panzhihua City (No. 2021ZD-G-1), Panzhihua Key Laboratory of Advanced Manufacturing Technology (No. 2022XJZD05), Key Program of Panzhihua University (2020ZD008)

==Conflicts of interests==

The authors declared that they have no conflicts of interest to this work.

==Nomenclature==

{| style="width: 81%;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;vertical-align: top;"|''R''
|  style="border-top: 1pt solid black;vertical-align: top;"|<span id='_Hlk41060195'></span>Cutter mean pitch radius
|-
|  style="vertical-align: top;"|''a''
|  style="vertical-align: top;"|Addendum coefficient
|-
|  style="vertical-align: top;"|''b''
|  style="vertical-align: top;"|Dedendum coefficient
|-
|  style="vertical-align: top;"|''W''
|  style="vertical-align: top;"|Face width of gear
|-
|  style="vertical-align: top;"|''N''
|  style="vertical-align: top;"|Number of teeth
|-
|  style="vertical-align: top;"|''m''
|  style="vertical-align: top;"|Module
|-
|  style="vertical-align: top;"|''ρ''
|  style="vertical-align: top;"|Edge radius coefficient
|-
|  style="vertical-align: top;"|''r''
|  style="vertical-align: top;"|Pitch radius of the gear
|-
|  style="vertical-align: top;"|''α<sub>n</sub>''
|  style="vertical-align: top;"|Normal pressure angle
|-
|  style="vertical-align: top;"|''S<sub>i</sub> ''
|  style="vertical-align: top;"|Coordinate system ''i'' with three mutual perpendicular axes ''x<sub>i</sub>'', ''y<sub>i</sub>'' and ''z<sub>i</sub>''
|-
|  style="vertical-align: top;"|'''M''' ''<sub>j,i</sub>''
|  style="vertical-align: top;"|Homogeneous matrix of coordinate transformation from ''S<sub>i</sub>'' to ''S<sub>j</sub>''
|-
|  style="vertical-align: top;"|'''r'''
|  style="vertical-align: top;"|Position vector of gear tooth surface
|-
|  style="vertical-align: top;"|'''n'''
|  style="vertical-align: top;"|Unit normal vector of gear tooth surface
|-
|  style="vertical-align: top;"|''ψ''
|  style="vertical-align: top;"|Generalized parameter of motion during generating process of gear
|-
|  style="vertical-align: top;"|''C''
|  style="vertical-align: top;"|Center distance
|-
|  style="vertical-align: top;"|''ΔC''
|  style="vertical-align: top;"|Center distance error
|-
|  style="vertical-align: top;"|''ΔA''
|  style="vertical-align: top;"|Axial displacement error
|-
|  style="vertical-align: top;"|''ΔH''
|  style="vertical-align: top;"|Crossing shaft angle error 
|-
|  style="vertical-align: top;"|''ΔV''
|  style="vertical-align: top;"|Intersecting shaft angle error
|-
|  style="vertical-align: top;"|''a<sub>pf</sub>''
|  style="vertical-align: top;"|<span id='_Hlk41060845'></span><span id='_Hlk41325039'></span>Parabola coefficient
|-
|  style="vertical-align: top;"|Δ''y<sub>tp</sub>''
|  style="vertical-align: top;"|<span id='_Hlk41143678'></span>Offset distance of the parabolic profile
|-
|  style="vertical-align: top;"|'''Superscripts:'''
|  style="vertical-align: top;"|
|-
|  style="vertical-align: top;"|''p'','' g''
|  style="vertical-align: top;"|Pinion and wheel of gear set
|-
|  style="vertical-align: top;"|'''Subscripts:'''
|  style="vertical-align: top;"|
|-
|  style="vertical-align: top;"|''P<sub>ib </sub>'','' G<sub>ib</sub>''
|  style="vertical-align: top;"|Convex side of the pinion and wheel
|-
|  style="border-bottom: 1pt solid black;vertical-align: top;"|''P<sub>ob </sub>'','' G<sub>ob</sub>''
|  style="border-bottom: 1pt solid black;vertical-align: top;"|Concave side of the pinion and wheel
|}

==References==

<span id='ref1'></span>1. Yang H-C, Tsai P-Y (2021) Mathematical model of S-shaped gear surface. Journal of Mechanical Science and Technology 35(7): 2841-2850. doi:10.1007/s12206-021-0607-z

2. Arafa HA (2016) C-Gears: Geometry and Machining. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 219 (7):709-726. doi:10.1243/095440605x31481

<span id='ref2'></span>3. Liu S-T (1988) Curvilinear cylindrical gears. Gear Technology 5 (3):8-12

<span id='ref3'></span>4. HA A, M B C-gears in rotorcraft transmissions: a novel design paradigm. In: International powered lift conference, Philadelphia, PA, USA, October 5-7, 2010 2010. pp 93-100

<span id='ref4'></span>5. HA A, M B Manufacturability and Viability of Different C-Gear Types: A Comparative Study

In: ASME (ed) Proceedings of the ASME Design Engineering Technical Conferences, Chicago, IL, USA, August 12-15, 2012 2012. ASME, pp 381-391

<span id='ref5'></span>6. Jia F, Hou L, Wei Y, Li B, You Y (2015) Modelling and bending strength analysis of cylindrical gears with arcuate tooth trace. Australian Journal of Mechanical Engineering 13 (2):77-86. doi:10.7158/m13-068.2015.13.2

<span id='ref6'></span>7. Tang R, Hou L, Zhou B (2015) Modelling and strength analysis for manufacture-oriented cylindrical gear with arcuate tooth trace. Australian Journal of Mechanical Engineering 14 (2):88-94. doi:10.1080/14484846.2015.1093217

<span id='ref7'></span>8. Zhijun S, Li H, Jinge W, Wei L, Qinglin C (2014) Contact strength analysis of circular-arc-tooth-trace cylindrical gear. Journal of the Brazilian Society of Mechanical Sciences and Engineering 38 (3):999-1005. doi:10.1007/s40430-014-0272-6

<span id='ref8'></span>9. Zhang X, Xie Y, Tan X (2016) Design, Meshing Characteristics and Stress Analysis of Cylindrical Gears with Curvilinear Tooth Profile. Transactions of FAMENA 40 (4):27-44

<span id='ref9'></span>10. Qinglin C, Hou L (2016) Parallel Translating Mechanism Process-Oriented Mathematical Model and 3-D Model for Cylindrical Gears with Curvilinear Shaped Teeth. Jordan Journal of Mechanical and Industrial Engineering 10 (3):171-177

<span id='ref10'></span>11. Zhang Q, Hou L, Tang R, Wen G (2016) Method of Processing and an Analysis of Meshing and Contact of Circular Arc Tooth Trace Cylindrical Gears. Transactions of FAMENA 40 (4):11-24. doi:10.21278/tof.40402

<span id='ref11'></span>12. Yongqiao W, Shuhong Y, Qi Z, Yulin W, Li H (2018) Numerical Analysis of Isothermal Elastohydrodynamic Lubrication of Cylindrical Gears with Variable Hyperbolic Circular Arc and Tooth Trace. Transactions of FAMENA 42 (1):61-72. doi:10.21278/tof.42106

13. Wei Y, Yang Y, Guo R, et al. (2022) Integrated wear prediction model for cylindrical gear with variable hyperbolic circular arc tooth trace under mixed elastohydrodynamic lubrication. Journal of Mechanical Science and Technology 36(8): 4053-4065. doi:10.1007/s12206-022-0726-1

<span id='ref12'></span>14. Litvin FL, Fuentes A (2004) Gear Geometry and Applied Theory. Cambridge University Press, Cambridge. doi:10.1016/0094-114X(95)90001-4

<span id='ref13'></span>15. Tseng R-T, Tsay C-B (2001) Mathematical model and undercutting of cylindrical gears with curvilinear shaped teeth. Mechanism and Machine Theory 36 (11):1189-1202. doi:10.1016/S0094-114X(01)00049-0

<span id='ref14'></span>16. Tseng R-T, Tsay C-B (2004) Contact characteristics of cylindrical gears with curvilinear shaped teeth. Mechanism and Machine Theory 39 (9):905-919. doi:10.1016/j.mechmachtheory.2004.04.006

<span id='ref15'></span>17. Tseng J-T, Tsay C-B (2005) Mathematical Model and Surface Deviation of Cylindrical Gears With Curvilinear Shaped Teeth Cut by a Hob Cutter. Journal of Mechanical Design 127 (5). doi:10.1115/1.1876437

<span id='ref16'></span>18. Tseng J-T, Tsay C-B (2006) Undercutting and Contact Characteristics of Cylindrical Gears With Curvilinear Shaped Teeth Generated by Hobbing. Journal of Mechanical Design 128 (3). doi:10.1115/1.2181605

<span id='ref17'></span>19. Wu Y-C, Chen K-Y, Tsay C-B, Ariga Y (2009) Contact Characteristics of Circular-Arc Curvilinear Tooth Gear Drives. Journal of Mechanical Design 131 (8):081003. doi:10.1115/1.3151805

<span id='ref18'></span>20. Chen Y-C, Li Z-W, Lo C-C, Wang Z-G (2017) A Study on The Improvement of Meshing Characteristics of A Complementary Curvilinear Gear Set Generated by Complementary Rack Cutters. Transactions of The Canadian Society for Mechanical Engineering 41 (2):281-291. doi:10.1139/tcsme-2017-1019

<span id='ref19'></span>21. Chen Y-C, Lo C-C (2014) Contact stress and transmission errors under load of a modified curvilinear gear set based on finite element analysis. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 229 (2):191-204. doi:10.1177/0954406214532907

<span id='ref20'></span>22. Fuentes A, Ruiz-Orzaez R, Gonzalez-Perez I (2014) Computerized design, simulation of meshing, and finite element analysis of two types of geometry of curvilinear cylindrical gears. Computer Methods in Applied Mechanics and Engineering 272:321-339. doi:10.1016/j.cma.2013.12.017

<span id='ref21'></span>23. Fuentes-Aznar A, Ruiz-Orzaez R, Gonzalez-Perez I (2016) Comparison of spur, helical and curvilinear gear drives by means of stress and tooth contact analyses. Meccanica 52 (7):1721-1738. doi:10.1007/s11012-016-0515-y

24. Wang G, Luo Q, Zou S (2021) Time-varying meshing stiffness calculation of an internal gear pair with small tooth number difference by considering the multi-tooth contact problem. Journal of Mechanical Science and Technology 35(9): 4073-4083. doi:10.1007/s12206-021-0819-2

25. Zhao Z, Yang Y, Ma H, et al. (2023) Meshing characteristics of spur gear pairs with tooth modification under different assembly errors and sensitivity analysis for impact factors. Journal of Mechanical Science and Technology 37(1): 149-162. doi:10.1007/s12206-022-1215-2

26. Li D, Liu Y, Gong J, et al.(2022) A novel method for longitudinal modification and tooth contact analysis of non-circular cylindrical gears. Journal of Mechanical Science and Technology 36(12): 6157-6170. doi:10.1007/s12206-022-1130-6