Electric vehicles have become the most important segment of the automobile market, and the continuous improvement of customer expectation has imposed increasing demands on the noise, vibration and harshness behaviour of the electric powertrain. The reducer, as the key transmission component of an electric vehicle, transforms the high-speed and high-torque output of the driving motor into a suitable operating range for the wheels. The vibration behaviour of the reducer is therefore one of the most decisive attributes of the overall vehicle comfort. In the present contribution, I investigate the influence of microscopic tooth surface errors, particularly surface roughness and waviness, on the mesh characteristics and on the system-level vibration response of an electric vehicle reducer. The work is based on the combination of elastohydrodynamic lubrication theory, potential energy equivalence, slicing techniques, dynamic substructure methods and design-of-experiment sensitivity analysis. Throughout the paper, the term “electric vehicle” is used to emphasize the target application, and all excitation and response analyses are discussed in the context of high-speed electric powertrains.
My research starts from the fact that in precision manufacturing of gears, micro-scale geometric deviations are inevitable. Surface roughness is mainly related to the random material removal process, while waviness is generated by periodic relative vibration between the grinding wheel and the workpiece. These deviations are often treated in existing simulation tools as random disturbances or even ignored. In practice, however, they produce systematic modifications of the instantaneous contact condition and of the friction force, which can either amplify or suppress the dynamic response of an electric vehicle reducer. The objective of this work is to establish an accurate and efficient computational framework that links the microscopic tooth surface error parameters to the final vibration level of the housing. The entire research is built on four interlocking steps: (i) topological modelling of roughness and waviness; (ii) analytical calculation of the time-varying mesh stiffness and static transmission error under non-ideal tooth surfaces; (iii) construction and validation of a rigid–flexible coupled dynamic model of the reducer; and (iv) application of global sensitivity analysis to identify the dominant surface error parameters in different speed ranges. All the numerical models are validated by test-bench measurements on a single-stage helical gear reducer intended for an electric vehicle.
1. Interpretation and mathematical representation of tooth surface errors
In this section, I explain how the tooth surface morphology is decomposed into roughness and waviness components. For a ground gear flank, the measured profile is usually a superposition of short-wavelength roughness and medium-wavelength waviness. The cut-off between roughness and waviness depends on the sampling length, and I therefore selected the profile arithmetic mean deviation \(R_a\) as the primary roughness index. This parameter is defined in the sampling length \(l\) as
\[
R_a = \frac{1}{l}\int_0^l \left|\Delta(x)\right| \mathrm{d}x,
\tag{1}
\]
where \(\Delta(x)\) is the departure of the measured profile from the least-squares mean line. In order to align with international practice, the evaluation length is taken as five times the sampling length, as recommended in the standard. Table 1 presents the recommended sampling lengths that I used for the subsequent calculation of the roughness-dependent friction coefficient.
| \(R_a\) (\(\mu\)m) | Sampling length \(l\) (mm) | Evaluation length \(l_n\) (mm) |
|---|---|---|
| 0.008 – 0.02 | 0.08 | 0.4 |
| 0.02 – 0.1 | 0.25 | 1.25 |
| 0.1 – 2 | 0.8 | 4 |
| 2 – 10 | 2.5 | 12.5 |
Table 1. Sampling and evaluation lengths recommended for roughness evaluation.
In order to introduce the roughness effect into the dynamics, I adopted an elastohydrodynamic lubrication approach to compute the instantaneous friction coefficient of the gear pair. The friction coefficient depends not only on roughness but also on the entrainment velocity \(V_e\), the slide-to-roll ratio \(\mathrm{SR}\), the equivalent curvature radius \(R\), the maximum Hertzian pressure \(P_h\) and the absolute viscosity \(\eta_0\). I used a semi-empirical relation in the form
\[
\mu(t) = e^{f(t)} b_0 \, [P_h(t)]^{b_1} |\mathrm{SR}(t)|^{b_2} [V_e(t)]^{b_3} \eta_0^{b_4} [R(t)]^{b_5},
\tag{2}
\]
where \(b_0,\dots,b_5\) are constants identified from numerical regression and \(f(t)\) is an auxiliary function depending on the instantaneous surface roughness and lubrication regime:
\[
f(t) = b_6 + b_7 \ln\left(\frac{h(t)}{R_a}\right) + b_8 e^{-|\mathrm{SR}(t)|P_h(t)\eta_0} + b_9 \ln(\eta_0 V_e(t)).
\tag{3}
\]
The friction coefficient changes continuously as the meshing point moves from the root to the tip. In particular, the relative sliding velocity reverses its sign at the pitch point; therefore, the friction force changes direction and becomes zero exactly at the pitch point. This time-varying characteristic makes the friction force a non-negligible source of parametric excitation inside an electric vehicle reducer.
2. Topological model of the wavy tooth surface
The second type of microscopic error investigated in this work is surface waviness. In cylindrical grinding, the periodic vibration between the workpiece and the grinding wheel causes a periodic variation of the material removal rate. If the workpiece experiences \(q\) radial oscillations per revolution, \(q\) waviness lobes are formed on the surface. The waviness signal is intrinsically periodic, and I represented it by a Fourier series. The complete waviness of a gear flank cannot be measured on one tooth alone because the measurement is performed tooth by tooth. To reconstruct the full topology, I introduced the gear rotation angle \(\varphi\). For a helical gear with helix angle \(\beta\), the axial separation between the two ends of a tooth is \(\Delta Z\), and the corresponding rotation angle increment is
\[
\Delta\varphi = \frac{2 \Delta Z \tan\beta}{D},
\tag{4}
\]
where \(D\) is the pitch diameter. By mapping every measured point onto a continuous rotation coordinate, the waviness of all teeth can be connected into a single closed curve. The resulting curve was then fitted by a sinusoidal compensation function. The waviness angle \(\theta\) is related to the phase variation in the profile direction, \(\mathrm{Ph\_p}\), and in the tooth direction, \(\mathrm{Ph\_h}\), through
\[
\tan\theta = \frac{\mathrm{Ph\_p}}{\mathrm{Ph\_h}}.
\tag{5}
\]
If the tooth is discretized into \(n_r\) thin slices in the face-width direction and the total number of teeth is \(Z\), the phase shift between adjacent slices of the same tooth is
\[
\Delta \psi_{\text{row}} = \frac{\mathrm{Ph\_h}}{n_r Z}.
\tag{6}
\]
With these definitions, I constructed the full waviness matrix of the whole gear. The value at a discrete location of the unfolded tooth surface is given by
\[
\mathrm{Matrix\_val}(x) = \frac{\mathrm{PP\_value}}{2} \sin\left[2\pi \frac{\mathrm{waves} \cdot r}{\mathrm{len}} \left(x + \Delta\psi_{\text{row}} (j-1)\right)\right],
\tag{7}
\]
where \(\mathrm{PP\_value}\) is the peak-to-valley waviness amplitude, \(\mathrm{waves}\) is the number of waviness cycles over the total developed length, \(r\) is the slice index and \(\mathrm{len}\) is the developed length. From this global matrix, the waviness of every tooth is extracted by picking the corresponding diagonal sub-matrix. The topology model described by equations (4)–(7) provides the geometric deviation input for the subsequent mesh stiffness and transmission-error calculation.
3. Time-varying mesh stiffness of the gear pair under tooth surface errors
The computation of the mesh stiffness is the first step to determine the internal dynamic excitation. For the reducer of an electric vehicle, I considered a single-stage helical gear pair. The basic geometry of the gear pair is listed in Table 2, and the main material parameters are listed in Table 3.
| Parameter | Symbol | Value |
|---|---|---|
| Number of pinion teeth | \(Z_1\) | 23 |
| Number of gear teeth | \(Z_2\) | 62 |
| Normal module (mm) | \(m_n\) | 1.77 |
| Active face width (mm) | \(b\) | 20 |
| Normal pressure angle (deg) | \(\alpha\) | 18.5 |
| Helix angle (deg) | \(\beta\) | 30 |
| Profile roughness (\(\mu\)m) | \(R_a\) | 0.8 |
Table 2. Basic parameters of the helical gear pair.
| Parameter | Unit | Value |
|---|---|---|
| Young’s modulus | \(\mathrm{N/m^2}\) | \(2.06\times10^{11}\) |
| Poisson’s ratio | – | 0.3 |
Table 3. Material properties of the gears.
For a spur gear slice, the potential energy method assumes that the tooth is a cantilever beam fixed at the root. Four deformation contributions were considered: Hertzian contact deformation, bending deformation, shear deformation and axial compression deformation. The corresponding potential energies \(U_h\), \(U_b\), \(U_s\) and \(U_a\) are
\[
U_h = \frac{F^2}{2k_h}, \quad U_b = \frac{F^2}{2k_b}, \quad U_s = \frac{F^2}{2k_s}, \quad U_a = \frac{F^2}{2k_a},
\tag{8}
\]
with \(k_h\), \(k_b\), \(k_s\) and \(k_a\) being the Hertzian, bending, shear and axial-compressive stiffness components respectively. The Hertzian contact stiffness is
\[
k_h = \frac{\pi E b}{4(1-\nu^2)},
\tag{9}
\]
where \(E\), \(b\) and \(\nu\) are the elastic modulus, face width and Poisson’s ratio. The bending, shear and axial stiffness components are expressed by integrals over the effective path from the root to the contact point. The classical expressions are
\[
k_b = \sum_{i=1}^{n} \frac{1}{\int_{\alpha_1}^{\alpha_2} \frac{3 \left\{1+\cos\alpha_1\left[(\alpha_2-\alpha)\sin\alpha – \cos\alpha\right]\right\}^2 (\alpha_2-\alpha)\cos\alpha}{2 E L \left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]^3} \,\mathrm{d}\alpha},
\tag{10}
\]
\[
k_s = \sum_{i=1}^{n} \frac{1}{\int_{\alpha_1}^{\alpha_2} \frac{1.2(1+\nu)(\alpha_2-\alpha)\cos\alpha \cos^2\alpha_1}{E L \left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]} \,\mathrm{d}\alpha},
\tag{11}
\]
\[
k_a = \sum_{i=1}^{n} \frac{1}{\int_{\alpha_1}^{\alpha_2} \frac{(\alpha_2-\alpha)\cos\alpha \sin^2\alpha_1}{2 E L \left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]} \,\mathrm{d}\alpha}.
\tag{12}
\]
In the above expressions, \(L\) represents the tooth width of the slice and \(n\) is the number of slices. Because friction modifies the internal force distribution on the cantilever, I introduced the friction coefficient \(\mu\) into the force decomposition. The normal mesh force \(F\) is resolved into two components: one perpendicular to the tooth centre-line and one parallel to it. When the friction force is added, additional bending moments and additional shear forces appear. Combining these effects, the modified bending, shear and axial stiffness terms become
\[
k_{b,f} = \frac{1}{\int_{\alpha_1}^{\alpha_2} \frac{3\left\{1+\cos\alpha_1\left[(\alpha_2-\alpha)\sin\alpha-\cos\alpha\right]\right\}^2(\alpha_2-\alpha)\cos\alpha}{2 E L \left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]^3} \left(1 \pm \mu \frac{\mu_{b}}{\mu} \right) \,\mathrm{d}\alpha},
\tag{13}
\]
\[
k_{s,f} = \frac{1}{\int_{\alpha_1}^{\alpha_2} \frac{1.2(1+\nu)(\alpha_2-\alpha)\cos\alpha \cos^2\alpha_1}{E L \left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]} \left(1 \pm \mu \frac{\mu_{s}}{\mu} \right) \,\mathrm{d}\alpha},
\tag{14}
\]
\[
k_{a,f} = \frac{1}{\int_{\alpha_1}^{\alpha_2} \frac{(\alpha_2-\alpha)\cos\alpha \sin^2\alpha_1}{2 E L \left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]} \left(1 \pm \mu \frac{\mu_{a}}{\mu} \right) \,\mathrm{d}\alpha},
\tag{15}
\]
where the \(\mu_b\), \(\mu_s\) and \(\mu_a\) coefficients arise from the decomposition of the friction vector in the bending, shear and axial directions respectively. The summation symbol indicates that all the slices contributing at the same angular position are taken into account. In addition to the tooth stiffness, the flexibility of the gear body is not negligible. The gear body stiffness \(k_f\) was estimated by the Sainsot–Velex formula:
\[
\frac{1}{k_f} = \frac{\cos^2\alpha_1}{E L} \left[ L^*\left(\frac{u_f}{S_f}\right)^2 + M^*\left(\frac{u_f}{S_f}\right) + P^*\left(1+Q^* \tan^2\alpha_1\right)\right],
\tag{16}
\]
in which \(u_f\) is the distance between the load application point and the root circle, \(S_f\) is the tooth-root arc length, and \(L^*\), \(M^*\), \(P^*\), \(Q^*\) are polynomial coefficients depending on the ratio \(h_f/r_f\). Finally, the single-pair mesh stiffness of a spur gear slice is obtained by considering all components in series:
\[
\frac{1}{k_{\text{single}}} = \frac{1}{k_h} + \frac{1}{k_{b,f}} + \frac{1}{k_{s,f}} + \frac{1}{k_{a,f}} + \frac{1}{k_f}.
\tag{17}
\]
Equation (17) is valid for a straight slice of the gear. To extend it to a helical gear, I applied the slice method. The helical gear is divided into \(n\) thin straight slices in the face-width direction. The contact line of a helical gear does not appear simultaneously over the whole width; it grows from zero to a maximum value and then decreases. The dimensionless mesh cycle was divided into two regimes according to the relationship between the transverse contact ratio \(\varepsilon_\alpha\) and the axial contact ratio \(\varepsilon_\beta\). The maximum contact line length is
\[
L_{\max} = \begin{cases}
\dfrac{B}{\cos\beta_b}, & \text{if } \varepsilon_\alpha > \varepsilon_\beta,\\[0.5em]
\dfrac{B}{\tan\beta_b}, & \text{if } \varepsilon_\alpha < \varepsilon_\beta,
\end{cases}
\tag{18}
\]
where \(\beta_b\) is the base helix angle and \(B\) is the total face width. Let \(\varepsilon_1=\min(\varepsilon_\alpha,\varepsilon_\beta)\) and \(\varepsilon_2=\max(\varepsilon_\alpha,\varepsilon_\beta)\). In a normalized mesh cycle of period \(T\), the time-varying contact line length \(L(t)\) can be expressed as
\[
L(t)=\begin{cases}
L_{\max}\dfrac{t}{\varepsilon_1 T}, & 0\le t\le \varepsilon_1 T,\\[0.5em]
L_{\max}, & \varepsilon_1 T \le t\le \varepsilon_2 T,\\[0.5em]
L_{\max}\left(1-\dfrac{t-\varepsilon_2 T}{\varepsilon_1 T}\right), & \varepsilon_2 T\le t\le (\varepsilon_1+\varepsilon_2)T.
\end{cases}
\tag{19}
\]
The total mesh stiffness of the gear pair at any instant is obtained by superimposing the contributions of all tooth pairs lying in the contact zone. Since the tooth pitch corresponds to one base-pitch spacing, the single-pair stiffness curve of successive tooth pairs is shifted in time by one base-pitch interval. If \(z_{\mathrm{act}}\) tooth pairs mesh simultaneously, the combined time-varying mesh stiffness can be written as
\[
k_m(t) = \sum_{j=1}^{z_{\mathrm{act}}} k_{\mathrm{single}}\left(t – (j-1)\Delta T\right),
\tag{20}
\]
where \(\Delta T\) is the base-pitch period. I implemented the above formulation for the gear pair of Table 2. Figure 1 (not referenced) shows a comparison between the stiffness computed with and without including the surface-friction correction. The inclusion of friction reduces the effective mesh stiffness because friction introduces additional deformation energy; nevertheless, the time-varying shape remains qualitatively similar.
4. Calculation of static transmission error under non-ideal tooth surfaces
Transmission error is the most direct measure of the deviation of a gear pair from perfect conjugacy. In the time domain, the transmission error along the line of action is defined as
\[
\Delta E = r_{b2}\theta_2 – r_{b1}\theta_1,
\tag{21}
\]
where \(r_{b1}\) and \(r_{b2}\) are the base radii of the pinion and the gear, and \(\theta_1\), \(\theta_2\) are the rotation angles of the two gears. For an ideal and rigid gear pair, this quantity is zero. In a real electric vehicle reducer, the elastic deformation and the tooth surface waviness produce a non-zero transmission error. Under an applied torque \(T_{\mathrm{in}}\), the static transmission error may be approximated initially by
\[
\mathrm{TE}_0 = \frac{P}{k_m}, \qquad P = \frac{T_{\mathrm{in}}}{r_{b1}},
\tag{22}
\]
where \(P\) is the nominal transmitted load. If a waviness amplitude \(\varepsilon\) exists at a particular contact point, the elastic deformation \(\delta\) and the transmission error \(x\) are linked by
\[
x = \delta + \varepsilon \quad \Rightarrow \quad \delta = x – \varepsilon.
\tag{23}
\]
When the contact is not established, \(\delta < 0\), the corresponding slice does not transmit load. Otherwise, the load on slice \(j\) of tooth \(i\) is determined from
\[
P_{i,j} = k_{i,j}\,\delta_{i,j} = k_{i,j}(x-\varepsilon_{i,j}),
\tag{24}
\]
and the total transmitted load over all active slices and tooth pairs must equal the external load \(P\):
\[
\sum_{i=1}^{m}\sum_{j=1}^{n} P_{i,j} = P.
\tag{25}
\]
Equations (23)–(25) must be solved iteratively because the unknown transmission error appears in the left-hand side through equation (24). At every mesh position, I start with the elastic value \(\mathrm{TE}_0\), compute the deformation of all slices, sum the individual loads, and compare the result with the applied load. If the computed load is lower than \(P\), the transmission error is increased by a small step \(\delta_{it}\); otherwise it is decreased. The iteration stops when the equilibrium condition is satisfied within a prescribed tolerance. The calculation flow is summarized as follows:
- Define the gear macro-geometry, load and surface error matrix.
- Compute the friction coefficient and the time-varying mesh stiffness.
- Initialize transmission error using the ideal formula (22).
- Determine, for every rotation angle, all active slices and their waviness error \(\varepsilon_{i,j}\).
- Compute \(P_{i,j}=k_{i,j}(x-\varepsilon_{i,j})\) and sum over all active slices.
- Compare the computed total load with the nominal load \(P\).
- Adjust \(x\) by \(\pm \delta_{it}\) until convergence.
- Record the converged value as the static transmission error at that angular position.
- Move to the next angular position and repeat the loop.
The algorithm yields the static transmission error curve over a complete mesh cycle. In order to validate the accuracy of the calculation, I carried out measurements on a reducer static transmission-error test stand. The test rig comprises a driving motor with a maximum speed of 60 rpm, a loading motor with a maximum torque of 1200 Nm, high-resolution angle encoders with an accuracy of \(\pm1.94”\), and the data-acquisition system. The reducer input shaft and output shaft were equipped with angle encoders, and tests were performed under two loading conditions listed in Table 4. The measured angular transmission error was converted into a length-equivalent value by multiplication with the base radius.
| Input speed (rpm) | Input torque (Nm) | Pinion acquisition cycles |
|---|---|---|
| 10 | 100 ±2% | 60 |
| 10 | 300 ±2% | 60 |
Table 4. Conditions used in the static transmission-error test.
The simulation was run with the same input speed and torque values. Figures of the resulting comparison are omitted here for brevity, but the mean relative error was used as a quantitative indicator:
\[
\mathrm{MRE} = \frac{1}{N}\sum_{i=1}^{N} \left|\frac{y_{\mathrm{sim},i}-y_{\mathrm{exp},i}}{y_{\mathrm{exp},i}}\right|\times 100\%,
\tag{26}
\]
where \(y_{\mathrm{sim},i}\) and \(y_{\mathrm{exp},i}\) are the simulated and measured transmission errors at the \(i\)-th sample, respectively. Under the 100 Nm load, the mean relative error obtained by the model including tooth surface errors was about 15.39%, while the ideal-surface model provided a relative error as large as 36.16%. Under the 300 Nm load, the error of the corrected model was 12.38%, compared with 29.19% for the uncorrected model. This comparison demonstrates that the iterative transmission-error algorithm, which accounts for waviness and friction, accurately extracts the internal displacement excitation of an electric vehicle reducer.
5. Dynamic substructure modelling of the reducer
The static transmission error is not sufficient to represent the vibration of the complete reducer because the rotating shafts, bearings, and housing all contribute to the dynamic behaviour. In this section I describe the rigid–flexible coupled dynamic model that I built for the electric-vehicle reducer under investigation. The main component bodies were discretized by the finite-element method and then reduced by dynamic substructuring. The physical model comprises the input shaft, the output shaft, the gear pair, four rolling bearings and the housing. The housing is modelled by a reduced stiffness matrix with four master nodes located at the bearing seats. The complete system is represented by 180 degrees of freedom: 13 nodes for the input shaft, 13 nodes for the output shaft and 4 nodes for the housing.
5.1 Shaft element
Each shaft was modelled with Timoshenko beam elements that include transverse shear deformation and rotary inertia. The nodes of the mesh are placed at the bearing centres and at the gear centres. Each shaft node carries six degrees of freedom: three translational displacements \(x,y,z\) and three rotations \(\theta_x,\theta_y,\theta_z\). The element displacement vector for a two-node shaft element is
\[
\mathbf{q}_e = \left[x_1,\ y_1,\ z_1,\ \theta_{x1},\ \theta_{y1},\ \theta_{z1},\ x_2,\ y_2,\ z_2,\ \theta_{x2},\ \theta_{y2},\ \theta_{z2}\right]^T.
\tag{27}
\]
The element stiffness matrix \(\mathbf{K}_s\) and the consistent mass matrix \(\mathbf{M}_s\) were constructed for a circular cross-section, using the elementary expressions for axial, torsional and bending deformation. For a hollow shaft of external diameter \(D\) and internal diameter \(d\), the area moment of inertia and polar moment of inertia are respectively
\[
I = \frac{\pi}{64}\left(D^4-d^4\right), \qquad J = \frac{\pi}{32}\left(D^4-d^4\right).
\tag{28}
\]
Damping in the shafts is represented by Rayleigh damping:
\[
\mathbf{C}_s = \alpha \mathbf{M}_s + \beta \mathbf{K}_s,
\tag{29}
\]
where the coefficients \(\alpha\) and \(\beta\) are selected according to the modal damping ratios in the frequency range of interest.
5.2 Gear-pair mesh element
The gear mesh element connects the pinion and gear nodes through a spring–damper set acting along the line of action. The projection vector \(\mathbf{V}\) is derived from the base helix angle \(\beta_b\), the pressure angle and the mounting phase angle \(\varphi\). The resulting relative deformation along the line of action is
\[
\delta = \mathbf{V} \mathbf{q}_m,
\tag{30}
\]
where \(\mathbf{q}_m\) is the vector containing the six degrees of freedom of both gears. The instantaneous mesh stiffness \(k_m\) computed in Section 3 is projected onto the physical degrees of freedom by the dyadic product
\[
\mathbf{K}_m = k_m \mathbf{V}\mathbf{V}^T.
\tag{31}
\]
The velocity-proportional mesh damping is determined by the mean mesh stiffness \(k_{\mathrm{mean}}\) and the equivalent masses of the gear pair:
\[
c_m = 2\zeta_m \sqrt{k_{\mathrm{mean}} \frac{m_{eq,p}m_{eq,g}}{m_{eq,p}+m_{eq,g}}},
\tag{32}
\]
where \(\zeta_m\) is a mesh damping ratio that usually lies between 0.03 and 0.17. The external torque is converted into a generalized nodal force vector
\[
\mathbf{F} = \mathbf{V} F,
\tag{33}
\]
with \(F\) being the instantaneous mesh force caused by the transmitted torque. The dynamic equation of the gear-pair element is:
\[
\mathbf{M}_m \ddot{\mathbf{q}}_m + \mathbf{C}_m \dot{\mathbf{q}}_m + \mathbf{K}_m \mathbf{q}_m = \mathbf{F}_m.
\tag{34}
\]
5.3 Bearing and housing elements
Each rolling bearing is modelled as a linear spring system between a shaft node and a housing master node. The bearing stiffness matrix \(\mathbf{K}_b\) is diagonal in the principal directions:
\[
\mathbf{K}_b = \mathrm{diag}\left(k_{xx}, k_{yy}, k_{zz}, k_{\theta_x\theta_x}, k_{\theta_y\theta_y}, k_{rr}\right).
\tag{35}
\]
The infinite-dimensional housing was first discretized with solid finite elements. Only four master nodes at the bearing seats were retained by static condensation. This operation produces an equivalent stiffness matrix \(\mathbf{K}_h\) of size \(24\times24\), coupling the four housing nodes:
\[
\mathbf{K}_h = \begin{bmatrix}
\mathbf{K}_{11} & \mathbf{K}_{12} & \mathbf{K}_{13} & \mathbf{K}_{14}\\
\mathbf{K}_{21} & \mathbf{K}_{22} & \mathbf{K}_{23} & \mathbf{K}_{24}\\
\mathbf{K}_{31} & \mathbf{K}_{32} & \mathbf{K}_{33} & \mathbf{K}_{34}\\
\mathbf{K}_{41} & \mathbf{K}_{42} & \mathbf{K}_{43} & \mathbf{K}_{44}
\end{bmatrix}.
\tag{36}
\]
The housing mass matrix \(\mathbf{M}_h\) was condensed similarly. Since the bearing mass is negligible compared with the shaft, gear and housing masses, I neglected its inertia contribution. The vector of global degrees of freedom is denoted by
\[
\mathbf{q} = \left[\mathbf{q}_{\mathrm{in}}^T,\ \mathbf{q}_{\mathrm{out}}^T,\ \mathbf{q}_{h}^T\right]^T,
\tag{37}
\]
where \(\mathbf{q}_{\mathrm{in}}\), \(\mathbf{q}_{\mathrm{out}}\) and \(\mathbf{q}_{h}\) refer to the input-shaft, output-shaft and housing master nodes. After assembling all element stiffness, damping and mass matrices, the global equations of motion of the electric vehicle reducer become
\[
\mathbf{M}\ddot{\mathbf{q}} + \mathbf{C}\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} = \mathbf{F}_s + \mathbf{F}_{\mathrm{mesh}},
\tag{38}
\]
where \(\mathbf{F}_s\) contains the externally applied motor torque and load torque, while \(\mathbf{F}_{\mathrm{mesh}}\) contains the internal mesh forces due to the time-varying mesh stiffness and static transmission error.
5.4 Modal analysis and time integration
Before solving the forced vibration, I performed a free-vibration analysis using the average mesh stiffness. Solving the eigenvalue problem
\[
\left(\mathbf{K} – \omega_i^2 \mathbf{M}\right)\boldsymbol{\phi}_i = \mathbf{0}
\tag{39}
\]
provides the natural frequencies \(f_i=\omega_i/(2\pi)\) and the mode shapes \(\boldsymbol{\phi}_i\). The first twenty natural frequencies are summarized in Table 5.
| Order | Frequency (Hz) | Order | Frequency (Hz) |
|---|---|---|---|
| 1 | 713.14 | 11 | 3351.69 |
| 2 | 785.24 | 12 | 3485.63 |
| 3 | 1018.40 | 13 | 3748.29 |
| 4 | 1257.43 | 14 | 3815.07 |
| 5 | 1694.24 | 15 | 4026.50 |
| 6 | 1796.95 | 16 | 4361.34 |
| 7 | 2207.31 | 17 | 4813.60 |
| 8 | 2225.33 | 18 | 5168.45 |
| 9 | 2580.78 | 19 | 5562.34 |
| 10 | 2814.37 | 20 | 6475.35 |
Table 5. Natural frequencies of the reducer model.
Because the time-varying stiffness and the mesh forces are periodic functions of the shaft rotation, I solved equation (38) directly in the time domain. To avoid a prohibitive computation in the physical coordinates, I transformed the equations into modal coordinates:
\[
\mathbf{q} = \boldsymbol{\Phi}\boldsymbol{\eta},
\tag{40}
\]
where \(\boldsymbol{\Phi}\) is the truncated modal matrix. Premultiplying equation (38) by \(\boldsymbol{\Phi}^T\), and using the orthogonality conditions, each modal coordinate satisfies a single-degree-of-freedom equation:
\[
M_i \ddot{\eta}_i + C_i \dot{\eta}_i + K_i \eta_i = F_i(t),
\tag{41}
\]
with the modal parameters \(M_i, C_i, K_i\) and the modal force \(F_i = \boldsymbol{\phi}_i^T \mathbf{F}\). The modal forces contain two contributions: the external torque projected onto the modal space and the internal mesh-force vector produced by the static transmission error \(\mathrm{TE}\). The internal part was represented as an equivalent displacement excitation:
\[
\mathbf{F}_{\mathrm{mesh}}(t) = k_m(t)\, \mathrm{TE}(t)\, \mathbf{V},
\tag{42}
\]
where \(\mathbf{V}\) is the projection vector expanded to the global dimension. In this way, every instantaneous value of the non-ideal transmission error is converted into a distributed mesh-force vector. The second-order modal equations were solved by the Newmark-\(\beta\) time integration scheme, with a constant time step small enough to resolve the highest frequency of interest. After all modal coordinates were obtained, the physical response was recovered through equation (40).
6. Vibration response analysis and experimental verification
During the operation of an electric vehicle reducer, the motor speed changes continuously; therefore, conventional constant-speed frequency analysis is not sufficient. I used order-tracking analysis to separate the rotationally synchronous components from structural resonances. In this technique, the vibration signal is resampled in the angular domain so that a component whose frequency is proportional to the shaft speed appears at a fixed order. The order \(O\), the frequency \(f\) and the shaft speed \(n\) in rpm are related by
\[
O = \frac{60 f}{n}.
\tag{43}
\]
Since the pinion of the studied reducer has 23 teeth, the gear meshing order is 23. I extracted the 23rd-order component from both the measured and simulated housing acceleration signals under an input torque of 50 Nm and a speed sweep from 600 to 6000 rpm. The experiments were performed on the same test bench used for the transmission-error measurements, with the reducer connected to a driving motor and a loading motor. Three orthogonal accelerometers were installed on the housing near the bearing seat of the output side. The accelerometer characteristics used in the tests are listed in Table 6.
| Parameter | Unit | Value |
|---|---|---|
| Range | \(g\) | \(\pm500\) |
| X-sensitivity | mV/g | 9.61 |
| Y-sensitivity | mV/g | 9.74 |
| Z-sensitivity | mV/g | 9.57 |
| Frequency range | Hz | 2 – 15000 |
| Resonant frequency | Hz | 30000 |
Table 6. Specifications of the tri-axial accelerometers used in the vibration test.
The predicted vibration level in the X, Y and Z directions was compared to the measured value at the bearing-seat location. The comparison showed that the main resonance intervals and the amplitudes of the 23rd-order component follow the same trend. The detailed amplitudes at the input-shaft and output-shaft bearing nodes are provided in Tables 7 and 8.
| Bearing position | Direction | Speed at maximum amplitude (rpm) | Maximum amplitude (m/s²) |
|---|---|---|---|
| Input shaft left bearing | X | 5800 | 7.32 |
| Input shaft left bearing | Y | 5800 | 13.16 |
| Input shaft left bearing | Z | 5800 | 1.37 |
| Input shaft right bearing | X | 4400 | 13.75 |
| Input shaft right bearing | Y | 4400 | 14.68 |
| Input shaft right bearing | Z | 5800 | 12.01 |
Table 7. Predicted response of the input-shaft bearing nodes at the 23rd mesh order.
| Bearing position | Direction | Speed at maximum amplitude (rpm) | Maximum amplitude (m/s²) |
|---|---|---|---|
| Output shaft left bearing | X | 5800 | 11.22 |
| Output shaft left bearing | Y | 5800 | 12.19 |
| Output shaft left bearing | Z | 5800 | 36.69 |
| Output shaft right bearing | X | 5800 | 16.61 |
| Output shaft right bearing | Y | 4400 | 14.75 |
| Output shaft right bearing | Z | 5800 | 55.69 |
Table 8. Predicted response of the output-shaft bearing nodes at the 23rd mesh order.
Both the numerical and experimental results reveal three important resonance regimes. The first regime appears around 1800–2200 rpm, where the mesh frequency approaches the first two structural modes (713 and 785 Hz). The second regime appears between 2700 and 3200 rpm, where the mesh frequency crosses the third and fourth modes (1018 and 1257 Hz). The strongest resonance occurs at about 4400 rpm and above. At 4400 rpm, the mesh frequency is 1685.87 Hz, which is very close to the fifth mode at 1694.24 Hz; therefore, the lateral vibrations of the right bearings reach their peak values. At 5800 rpm, the mesh frequency is 2222.28 Hz, simultaneously exciting the seventh and eighth modes, which produces the largest housing vibration levels, especially in the Z direction at the output-shaft right bearing, reaching an order-component amplitude of 55.69 m/s² in the numerical result. The test data confirmed the same dominant speed ranges, although the absolute values show slight differences due to measurement location and structural damping assumptions.
The good qualitative agreement between the simulation and the test indicates that the rigid–flexible coupled model can be used to investigate the influence of manufacturing error parameters on the electric-vehicle-reducer vibration. A comparison of the X-direction, Y-direction and Z-direction waterfall slices showed that the peak amplitudes and their speed positions are faithfully reproduced by the model. Therefore, the model is considered suitable for the subsequent parameter study.
7. Influence of surface roughness and waviness on the vibration of the electric vehicle reducer
After validating the dynamic model, I used it to study the effect of the microscopic gear-surface characteristics on the vibration response. The analysis was performed at an input torque of 50 Nm and a speed of 4400 rpm, which corresponds to a typical high-speed cruising condition of the electric vehicle reducer. The three independent parameters were the surface roughness \(R_a\), the waviness amplitude (peak-to-valley value) \(\mathrm{PP\_value}\) and the waviness angle \(\theta\).
7.1 Effect of surface roughness
Four values of the roughness, \(R_a = 0.4, 0.8, 1.2, 1.6\ \mu\mathrm{m}\), were simulated while keeping the waviness parameters constant. The root-mean-square and peak-to-peak values of the dynamic transmission error increase with roughness, although the change is relatively small. In the considered range, the root-mean-square of the dynamic transmission error increases by about 0.15%, while the peak-to-peak value increases by about 0.73%. The perceived vibration level is more affected than the static level because the roughness changes the local friction force, hence the high-frequency excitation. The housing acceleration in the X direction shows the largest sensitivity to roughness; the peak-to-peak value rises by about 1.32% and the RMS value by 1.41% when \(R_a\) is changed from 0.4 to 1.6 \(\mu\mathrm{m}\). The corresponding increases in the Y direction are 1.12% and 1.16%, while the Z-direction increase remains below 0.6%. This direction-dependent effect is explained by the fact that friction modifies the axial and radial projections of the mesh force in a complex manner.
7.2 Effect of waviness amplitude
The waviness amplitude is a direct measure of the periodic deviation of the tooth flank from its ideal involute shape. I varied the peak-to-valley waviness amplitude from 2.5 to 10 \(\mu\mathrm{m}\). The dynamic transmission error exhibits a nearly linear growth of its peak-to-peak value, from about 2.0 \(\mu\mathrm{m}\) at 2.5 \(\mu\mathrm{m}\) to about 7.0 \(\mu\mathrm{m}\) at 10 \(\mu\mathrm{m}\), corresponding to an increase of approximately 247%. The root-mean-square value of the transmission error varies only weakly, because waviness primarily causes oscillations around a mean elastic deformation. The housing acceleration amplitudes increase significantly in all three directions. In the Z direction, the peak-to-peak acceleration grows by about 208% and the RMS by 156% when the waviness amplitude is raised from 2.5 to 10 \(\mu\mathrm{m}\). This indicates that the waviness amplitude strongly excites the axial mode of the reducer.
7.3 Effect of waviness angle
The waviness angle is defined as the inclination of the waviness lines with respect to the tooth trace. For the studied gear, the base helix angle is approximately \(28^\circ\). When the waviness angle approaches the base helix angle, the contact lines periodically cross the waviness peaks and valleys, producing a kind of “overlap resonance” that amplifies the transmission-error fluctuation. I tested angles \(10^\circ, 20^\circ, 30^\circ\) and \(40^\circ\). The largest transmission-error fluctuation occurs at \(30^\circ\), with a peak-to-peak value of 0.833 \(\mu\mathrm{m}\), compared to values of 0.415, 0.450 and 0.452 \(\mu\mathrm{m}\) at \(10^\circ\), \(20^\circ\) and \(40^\circ\), respectively. The same angular alignment leads to the highest housing vibration, especially in the Z direction, where the peak-to-peak acceleration at \(30^\circ\) is about 156% larger than the minimum value at \(20^\circ\). Thus, the waviness angle is an important parameter for tuning the meshing excitation in high-speed electric vehicle applications.
8. Global sensitivity analysis of the tooth-surface-error parameters
The previous single-parameter analysis provides qualitative trends, but in a real manufactured gear all three parameters vary simultaneously. To quantify the relative importance of each parameter, I integrated the dynamic model into a design-of-experiment framework using Latin hypercube sampling. The input variables are \(R_a\), \(\mathrm{PP\_value}\) and \(\theta\), and the output variables are the peak-to-peak value and the root-mean-square value of the housing vibration acceleration. The Latin hypercube method was chosen because it provides a good space-filling property with a limited number of simulations. The contribution of each input variable was evaluated by means of analysis of variance and visualized in Pareto charts.
8.1 Results at medium speed (4400 rpm)
At 4400 rpm, the contributions of the three surface-error parameters to the peak-to-peak acceleration are listed in Table 9. The waviness angle is the most influential parameter, with a contribution of 44.80%. The waviness amplitude and surface roughness contribute 33.26% and 21.94%, respectively. Therefore, if the vibration peak amplitude is excessive, the waviness angle should be controlled first. In contrast, the RMS acceleration at 4400 rpm is dominated by the surface roughness, whose contribution is 39.62%, followed by the waviness angle (32.06%) and the waviness amplitude (28.32%). This shows that the surface roughness mainly controls the vibration energy, while the waviness angle mainly controls the vibration fluctuation.
| Speed (rpm) | Response metric | \(R_a\) contribution (%) | Waviness amplitude contribution (%) | Waviness angle contribution (%) |
|---|---|---|---|---|
| 4400 | Peak-to-peak | 21.94 | 33.26 | 44.80 |
| 4400 | RMS | 39.62 | 28.32 | 32.06 |
| 8000 | Peak-to-peak | 7.72 | 54.03 | 38.24 |
| 8000 | RMS | 24.74 | 53.80 | 21.46 |
Table 9. Contribution rates of the tooth-surface-error parameters to the housing acceleration vibration indicators at different speed conditions.
The contrast between the two indicators demonstrates that the same geometric error cannot be optimized simultaneously for both peak fluctuations and total energy. For example, reducing the surface roughness is an effective way to lower the RMS acceleration, but it has a weaker effect on the peak acceleration; conversely, optimizing the waviness angle is a priority when the main problem is a strong periodic beating of the housing.
8.2 Results at high speed (8000 rpm)
Because modern electric vehicle reducers often run at extremely high speeds, I repeated the DOE analysis at 8000 rpm. The results, also included in Table 9, reveal a dramatic change in the ranking of the parameters. At 8000 rpm, the waviness amplitude becomes the most important parameter for both the peak-to-peak and RMS acceleration, with contributions of 54.03% and 53.80%, respectively. The waviness angle is the second most important parameter for the peak-to-peak value (38.24%), while the surface roughness is more important than the waviness angle for the RMS value (24.74% vs. 21.46%). The contribution of surface roughness to the peak-to-peak acceleration drops to only 7.72%, indicating that friction-induced excitation is less amplified than the geometrical waviness excitation at very high mesh frequencies. This trend can be explained by the frequency content of the excitations: waviness with a large amplitude creates a stronger displacement excitation at multiples of the mesh frequency, and these components eventually coincide with the structural modes of the electric vehicle reducer. The high-speed operation therefore shifts the dominating error source from the waviness direction to the waviness magnitude.
8.3 Comparison and design implication
Table 10 summarizes the dominant parameters in the two speed ranges. In the medium-speed range, the vibration of the housing is mainly governed by the orientation of the waviness for peak signals and by the frictional influence of the surface roughness for the RMS signal. At high speed, the waviness amplitude becomes the most critical factor for both indicators. This finding implies that a dedicated manufacturing specification should be used for high-speed electric vehicle reducers: an upper limit for the waviness amplitude should be prioritized, while the control of the surface roughness may be relaxed if the operation is dominated by high-speed conditions.
| Operating condition | Dominant parameter for peak-to-peak | Dominant parameter for RMS |
|---|---|---|
| 4400 rpm | Waviness angle | Surface roughness |
| 8000 rpm | Waviness amplitude | Waviness amplitude |
Table 10. Dominant parameters at medium and high speeds.
9. Concluding remarks and outlook
In this work I have systematically investigated the influence of microscopic tooth-surface geometry on the vibration of an electric vehicle reducer. The main contributions can be summarized as follows:
- I established a complete mathematical description of the tooth surface including the roughness and the waviness. The roughness was introduced through the arithmetic mean deviation \(R_a\), and the waviness was represented by a Fourier-series-based full-topology matrix. The proposed model allows direct simulation of measured tooth-surface error patterns without introducing artificial random disturbances.
- I developed an analytical approach for the mesh characteristics that combines the potential energy method, the slice method and a friction-force decomposition. The corrected time-varying mesh stiffness reflects the additional compliance caused by friction. The iterative transmission-error algorithm, accounting for the local waviness deformation, produces accurate predictions of the static transmission error as confirmed by test-bench data.
- I constructed a rigid–flexible coupled dynamic model of the electric vehicle reducer with 180 degrees of freedom. The model includes Timoshenko shafts, a six-degree-of-freedom gear-pair element, bearing stiffness, and a condensed housing model. The Newmark-\(\beta\) integration combined with modal superposition provided stable solutions for the full-speed vibration response. Order-tracking analysis of the 23rd-order component reproduced the measured resonance regions and amplitude trends.
- I carried out a global sensitivity analysis using the design-of-experiment method. The result demonstrates that the ranking of the error parameters is strongly speed-dependent. At medium speed, the waviness angle governs the vibration peaks, while the surface roughness governs the vibration RMS. At high speed, the waviness amplitude becomes the most critical parameter for both indicators.
The findings provide clear guidance for the low-vibration manufacturing of electric vehicle reducers. In particular, when an electric vehicle reducer is intended to operate at high rotational speeds, the waviness amplitude should be strictly controlled, while the surface roughness requirement can be balanced with the manufacturing cost. Conversely, the waviness angle should be controlled near the base helix angle to avoid the overlap-resonance effect that amplifies transmission error and housing vibration.
There are further research directions that I plan to explore in the future. One important extension is the inclusion of wear evolution because the roughness and waviness change during operation. Another direction is the development of machine-learning surrogates derived from the current dynamic model to enable multi-objective optimization of tooth-surface error parameters. A digital-twin approach could combine the dynamic model with real-time sensor information to predict the state of the electric vehicle reducer over its entire lifetime. These future studies will improve the fidelity and efficiency of the simulation framework and will bring the virtual design of electric-vehicle reducers closer to industrial application.