A method for quantitatively adjusting contact marks of an arc tooth bevel gear

Abstract

The application discloses a kind of contact imprint quantization adjustment methods of spiral bevel gear, it is related to imprint quantization technical field, comprising: using incremental loading and gradually iterative solution obtains bearing contact state distribution data set and bearing transmission error curve;Obtain the topological mapping model of modified tooth surface contact stiffness;Frequency-dependent dynamic sensitivity matrix is constructed;Based on dimension reduction bias field and frequency-dependent dynamic sensitivity matrix, the optimal processing parameter correction is obtained by genetic algorithm, and the optimal processing parameter correction is written into the global variable register of numerical control system by parameterized programming interface.This application describes the transmission influence of the unit change of each processing parameter on the normal contact force and transmission error of any contact node pair, and continuously weights the participation coefficient of boundary contact node, which converts the original mutation and non-differentiable contact boundary movement into a smooth continuous process, solving the problem of unstable boundary sensitivity and difficult convergence of optimization.

Description

Technical Field

[0001] This invention relates to the field of imprint quantification technology, specifically a method for quantifying and adjusting contact imprints of spiral bevel gears. Background Technology

[0002] Spiral bevel gears are key components for power transmission between intersecting or staggered shafts and are widely used in transmission systems. The position, shape, and size of their contact marks directly determine the transmission smoothness, noise, and load-bearing life of the gear pair. In actual production, due to systematic errors in milling parameters, tooling errors, and heat treatment deformation, there will always be deviations between the tooth surface geometry and the theoretical design. Therefore, after gear pairing, contact mark tests on a rolling inspection machine are necessary to repeatedly adjust machining parameters or mounting distances to correct the marks to the ideal position.

[0003] In existing technologies, adjusting contact marks relies heavily on the operator's experience, making it a typical trial-and-error process. Technicians observe the scratched or stained marks on the tooth surface and intuitively judge the direction and magnitude of adjustments needed to parameters such as radial tool position and rolling ratio. They then re-grind the tool, adjust the machine tool, or perform another rolling inspection.

[0004] However, the mapping relationship between minute changes in machining parameters and normal clearance at various points on the tooth surface, and how these changes nonlinearly affect the pressure distribution and transmission error on the final contact imprint, is unclear. Although traditional tooth surface contact analysis and load-bearing contact analysis can perform forward problem calculations, the computational load is large and highly nonlinear, making it impossible to directly use them to inversely calculate corrections. Furthermore, when actual gears operate at high speeds, machining parameter errors are transformed into periodic geometric displacement excitations related to the meshing frequency and its harmonics, causing significant fluctuations in dynamic contact force and transmission error. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for quantitatively adjusting the contact imprint of spiral bevel gears, comprising the following specific steps: S1. Based on the machining parameters of the large and small wheels fixed in the global installation coordinate system, the discretized grid points of the tooth surface are generated and the initial normal clearance is calculated. Based on the local contact stiffness of the tooth surface, the mechanical equilibrium equations of the flexible contact of the tooth surface are constructed. The load-bearing contact state distribution dataset and the load-bearing transmission error curve are obtained by incremental loading and step-by-step iterative solution. S2. Based on the final state in the load-bearing contact state distribution dataset as the reference point for linearization expansion, a first-order Taylor expansion of the nonlinear mechanical relationship is performed to obtain a linear relationship. The influence of machining parameter perturbation on the geometric boundary is introduced, and a mapping model is formed by combining global equations. Based on the preset transition zone threshold, the critical potential contact nodes identified in the load-bearing contact state distribution dataset are weighted and corrected to obtain the corrected tooth surface contact stiffness topology mapping model. S3. The changes in each processing parameter are transformed into periodic geometric displacement excitation functions related to the meshing frequency and harmonics. The numerical extension solution is performed based on the nonlinear time-domain dynamic equation. The obtained additional normal displacement compensation is superimposed point by point onto the initial normal clearance of all potential contact node pairs used in the load contact state distribution dataset. Based on the modified tooth surface contact stiffness topology mapping model, the updated load contact state and transmission error are obtained, and the frequency-related dynamic sensitivity matrix is ​​constructed. S4. Quantitative feature extraction is performed based on the fluorescent contact imprint image set to obtain the mechanical feature deviation vector of the bearing contact imprint. Then, the tooth surface morphology is measured and interpolated to form the tooth surface normal deviation data matrix. Then, the tooth surface normal deviation data matrix is ​​decomposed into intrinsic orthogonal modes. The obtained dominant intrinsic orthogonal decomposition modes and mode coefficients are linearly superimposed to reconstruct the reconstructed dimensionless deviation field. Based on the dimensionless deviation field and the frequency-related dynamic sensitivity matrix, the optimal machining parameter correction amount is obtained through a genetic algorithm. The optimal machining parameter correction amount is written into the global variable register of the CNC system through the parametric programming interface.

[0006] Preferably, S1 establishes the tooth surface equations of the large gear and the small gear respectively through a homogeneous coordinate transformation matrix based on the machining parameters of the large gear and the small gear. The machining parameters include radial tool position, angular tool position, vertical gear position, bed position correction, gear position correction, blank mounting angle, rolling ratio, tool tip radius, and cutter head tooth profile angle. The homogeneous coordinate transformation matrix unifies the tooth surfaces of the large gear and the small gear from spatial coordinates to a fixed global mounting coordinate system. This spatial coordinate transformation process strictly follows the design axis intersection angle and design offset distance of the gear pair. The design axis intersection angle refers to the design angle between the axis of the large gear and the axis of the small gear in space, and the design offset distance refers to the shortest vertical distance between the axis of the large gear and the axis of the small gear. The mesh is divided along the tooth length and tooth height directions under a fixed global installation coordinate system. Specifically, in the tooth tip and tooth root edge regions where the absolute value of the tooth surface curvature gradient is greater than the overall median of the tooth surface curvature gradient, the mesh spacing is set to a small proportion of the tooth surface working height. In the gentle middle region of the tooth surface where the absolute value of the tooth surface curvature gradient is less than and / or equal to the overall median of the tooth surface curvature gradient, the mesh spacing is set to a large proportion of the tooth surface working height, so that the mesh density in the edge region is greater than the mesh density in the middle region. The generated one-to-one correspondence of the mesh points on the large gear tooth surface and the mesh points on the small gear tooth surface, when the large gear and the small gear are at the theoretical meshing angle position, the mesh points in any row and any column on the large gear tooth surface and the mesh points in the same row and the same column on the small gear tooth surface are most close in spatial position, forming a potential contact node pair; For each potential contact node pair generated, calculate the coordinate difference between the grid point on the large gear tooth surface and the grid point on the small gear tooth surface along the normal direction of that point on the large gear tooth surface. Use the absolute value of the coordinate difference as the initial normal clearance of the potential contact node pair. Record all potential contact node pairs and their initial normal clearances as a set of full-tooth-surface potential contact node pairs and their initial normal clearances.

[0007] Preferably, in step S1, the normal contact force required to produce a unit normal elastic deformation per unit area for a pair of contact nodes on the large and small wheels is taken as the local contact stiffness of the tooth surface. The value of the local contact stiffness of the tooth surface is calculated in advance for each node pair by Hertz contact theory, combined with the elastic modulus and Poisson's ratio of the gear material measured by standard tensile test. Based on the full tooth surface potential contact node pairs and the initial normal clearance set, the criteria for mutual contact and the definition of mechanical relationship for each potential contact node pair are defined. Specifically, when the sum of the relative elastic deformation of the mesh points of the large wheel and the small wheel under the action of external load exceeds the initial normal clearance, it is determined that they have entered the contact state. The difference between the normal contact force and the sum of the relative elastic deformation of the large wheel and the small wheel in the contact state minus the initial normal clearance is expressed by the nonlinear increasing function relationship determined by Hertz contact theory. The nonlinear increasing function relationships of all potential contact node pairs are combined and added to the macroscopic torque balance equation of the gear pair. The macroscopic torque balance equation of the gear pair means that the sum of the torques of the normal contact forces of all contact node pairs on the gear rotation axis is equal to the input torque, forming a global nonlinear mechanical balance equation system. The unknowns of this equation system are the normal contact forces borne by each potential contact node pair and the macroscopic load transmission error of the gear pair caused by the elastic deformation of the tooth surface. The load transmission error is defined as the deviation of the actual output angle of the large gear from the theoretical output angle of the large gear under the condition of constant input torque of the small gear. The numerical method based on incremental loading and step-by-step iteration solves the global nonlinear mechanical equilibrium equations. The total working torque is divided into several load increment steps starting from zero. The division rule is as follows: in the light load range where the contact state changes gently, a preset constant large step size with the total working torque as the reference is used for division; in the heavy load range where the contact state changes drastically, the step size is divided by successively decreasing. The decreasing factor makes the subsequent step size decrease in a geometric series. In each load increment step, the initial set of contact node pairs for the current step is obtained first based on the bearing contact state obtained by convergence at the end of the previous load increment step. The global nonlinear mechanical equilibrium equations are expanded using a first-order Taylor series under the initial set of contact node pairs to obtain linearized compliance equations. These compliance equations are then solved using a sparse matrix direct solver based on multi-wavefront LU decomposition to obtain trial solutions for the normal contact force and load transmission error of each contact node in the current step. The initial set of contact node pairs is then checked and updated based on these trial solutions. Initial contact node pairs with normal contact forces less than or equal to zero are removed, and new contact node pairs with closed normal gaps are added. All incremental steps are completed sequentially, and the records are merged to obtain the load contact state distribution dataset and the load transmission error curve under the full working torque. The load contact state distribution dataset contains the spatial coordinates of all contact node pairs in the final step, the normal contact force values, and the spatial boundary coordinates of the load contact area defined by these coordinates.

[0008] Preferably, S2 uses the final state obtained by convergence under the final load increment step, based on the concentrated load contact state distribution data, as the linearization reference point. For each potential contact node pair, the nonlinear functional relationship between the normal contact force, the sum of relative elastic deformation, and the initial normal gap is expanded using a first-order Taylor series to obtain a linear relationship. It should be noted that the linear relationship means that the small change in normal contact force is equal to the small change in the local contact stiffness of the node multiplied by the sum of relative elastic deformation, minus the small change in the local contact stiffness multiplied by the initial normal gap. Then, the initial normal of each potential contact node pair is obtained based on the homogeneous coordinate transformation matrix. The partial derivatives of the clearance with respect to each machining parameter are substituted into the linear relationship to directly express the change in normal contact force as a linear combination of the small perturbations of each machining parameter. The linearized forms of the linear relationship with respect to all potential contact nodes and the macroscopic torque balance equation are combined to form a global linear equation system, which forms the tooth surface contact stiffness topology mapping model. In this application, the coefficient matrix of this equation system is the tooth surface contact stiffness topology mapping model. Each element quantitatively describes the change in normal contact force or macroscopic load transmission error caused by the unit perturbation of a specific machining parameter by changing the initial normal clearance of the entire tooth surface.

[0009] Preferably, in step S2, all potential contact node pairs located on the spatial boundary of the bearing contact zone in the bearing contact state distribution dataset are extracted, and the statistical distribution of the ratio of normal contact force to local contact stiffness of all potential contact node pairs is calculated. The lower quantile and upper quantile of the ratio are used as the lower limit and upper limit of the transition zone threshold, respectively. In some other specific embodiments, the lower quantile and upper quantile are selected based on the attenuation characteristics of the probability density function of the statistical distribution at the tail. For each pair of critical potential contact nodes located at the spatial boundary of the bearing contact zone in the bearing contact state distribution dataset, the logistic function is used as a smooth step function. Its steepness parameter is calibrated through contact imprint rolling test of the same batch of gears. The goal is to minimize the overall mean of the sum of squares of the difference between the boundary movement predicted by the tooth surface contact stiffness topology mapping model and the measured movement. The result is determined by single variable scanning optimization. When solving, if the sum of the normal relative position and deformation of the critical potential contact node pair is less than the initial normal clearance and the difference is greater than the upper limit of the transition zone, the result is zero. If the sum of the normal relative position and deformation of the node pair is greater than the initial normal clearance and the difference is greater than the upper limit, the result is one. The transition is smooth within the transition zone to obtain the contact participation coefficient. Based on the contact participation coefficient as a weight multiplied by the local contact stiffness of the corresponding critical potential contact node pair, the diagonal sub-blocks related to the contact stiffness in the coefficient matrix of the global linear equation system of the tooth surface contact stiffness topology mapping model are weighted and corrected to obtain the tooth surface contact stiffness topology mapping model weighted and corrected by the contact participation coefficient. Thus, the correction in this application makes the contact boundary movement a continuously differentiable process.

[0010] Preferably, step S3 transforms the changes in each processing parameter into a periodic geometric displacement excitation function related to the meshing frequency and harmonics, so that the steady-state vibration response of the large wheel and the small wheel is expanded into a Fourier series. The highest order of expansion is determined by gradually increasing it from the fundamental wave until the amplitude of the newly added harmonics decays to the fundamental wave amplitude. The vibration response in the form of the Fourier series is substituted into the nonlinear time-domain dynamic equation. The nonlinear time-domain dynamic equation is transformed into a set of coupled nonlinear algebraic equations using the harmonic balance principle. The unknowns are the amplitude and phase of each harmonic. The equations are solved based on the numerical extension algorithm. Starting from the steady-state solution without perturbation, the perturbation amplitude of the processing parameter is gradually increased from zero by a preset proportion of the maximum allowable perturbation of the processing parameter as a fixed step size. Each step uses the convergent solution of the previous step as the initial value of the Newton iteration. The iteration continues until the relative change of the Euclidean norm of the solution vector is less than a minimum threshold, thus obtaining the harmonic amplitude of the vibration response under each extension step. Based on the obtained fluctuation amplitude, it is converted into periodic displacement fluctuation along the meshing line, which is used as an additional normal displacement compensation amount. It is then superimposed point by point onto the initial normal clearance of all potential contact node pairs used in the load contact state distribution dataset. It is substituted into the tooth surface contact stiffness topology mapping model after the contact participation coefficient is weighted and corrected, and the incremental loading and step-by-step iterative solution of S1 are repeated to obtain the updated load contact state and transmission error. After repeating the operation for each processing parameter, the coordinates of the center position of the bearing contact area before and after, the non-uniformity of the normal force distribution, and the change in the amplitude of the specific order harmonic of the bearing transmission error are compared to construct a two-dimensional numerical matrix as a frequency-related dynamic sensitivity matrix. Its rows correspond to the processing parameters, the columns correspond to the evaluation index of the bearing contact dynamic response under the current working condition, and the elements are the quantitative changes in the index caused by the perturbation of a unit processing parameter.

[0011] Preferably, step S4 involves uniformly spraying a dry powder fluorescent detection agent onto the tooth surfaces of the large and small gears using a rolling inspection machine. The dry powder fluorescent detection agent is calibrated by a laser thickness gauge to the median of the recommended optimal development thickness range for the fluorescent detection agent. The same preset total working torque used in S1 is applied to the gear pair for rolling inspection. The number of rotations is set to the minimum number of rotations required to ensure sufficient stability of the contact imprint. The tooth surface is excited by a narrow-band filtered ultraviolet light source with a center wavelength in the ultraviolet band. A high-resolution monochrome industrial camera is used in a darkroom environment to acquire several images from multiple angles within a preset symmetrical angle range with respect to the tooth surface normal. After averaging and noise reduction, a fluorescent contact imprint image set is formed. For each image in the fluorescent contact imprint image set, an adaptive local thresholding method is used for binarization. The neighborhood window size is taken as a proportion of the length of the minor axis of the theoretical equivalent ellipse of the contact imprint predicted by S1. The pixel threshold is taken as a multiple of the Gaussian weighted gray average value in the neighborhood minus the standard deviation of the gray values ​​in the neighborhood. An active contour model based on gradient vector flow is applied. The initial contour is set as the convex hull of the gray-level abrupt change region in the binary image. When the contour convergence criterion is satisfied, it converges to the sub-pixel accuracy boundary. The contour convergence criterion is satisfied based on the overall change rate of the area enclosed by the contour in the previous two iterations. The area, equivalent ellipse major axis, minor axis, center coordinates and tilt angle are calculated based on the sub-pixel accuracy boundary. Based on the pre-calibrated fluorescence intensity-normal contact pressure proportional relationship using pressure-sensitive film, the grayscale of pixels inside the sub-pixel precision boundary is converted and integrated point by point to reconstruct the relative normal force distribution map. The relative normal force distribution map is then compared with the ideal uniform force distribution constructed by the linearized reference point to calculate the center point coordinate offset vector and the non-uniformity statistics of the force distribution, thus forming the mechanical characteristic deviation vector of the bearing contact imprint.

[0012] Preferably, step S4 integrates the line laser profile scanning probe to the side of the rolling inspection machine, and performs reciprocating scanning of each tooth groove of the large and small wheels. In the specific implementation process, the scanning line spacing is set to a value that balances accuracy and efficiency based on the optical resolution of the probe and the curvature of the tooth surface. Two-dimensional profile lines are acquired in real time, and all profile lines are registered to a fixed global installation coordinate system through an iterative nearest point algorithm. In each iteration, the nearest corresponding point is searched using a KD tree. The optimal rigid body transformation matrix is ​​solved by the singular value decomposition method to minimize the sum of squared distances between point pairs. The process continues until the average distance change between adjacent iterations is less than a preset proportion of the nominal accuracy of the line laser profile scanning probe, at which point convergence occurs. The installation eccentricity error is compensated synchronously. The registered measured point cloud is then matched with the theoretical tooth surface for best fit. The deviation value of each measured point along the tooth surface normal is calculated and mapped onto the tooth surface discretization grid through Kriging space interpolation to form a tooth surface normal deviation data matrix.

[0013] Preferably, in step S4, the tooth surface normal deviation data matrix is ​​decomposed into an intrinsic orthogonal matrix. The tooth surface normal deviation data matrix is ​​straightened into a high-dimensional column vector. A deviation sample matrix is ​​constructed based on the deviation column vectors of the current measurement sample and multiple samples accumulated from historical measurements. After calculating the covariance matrix of the deviation sample matrix, the dimension of the covariance matrix is ​​equal to the number of tooth surface grid points. Using the implicit restart Arnoldi iterative algorithm suitable for solving the eigenvalues ​​of large sparse matrices, all eigenvalues ​​of the covariance matrix and the eigenvectors corresponding to each eigenvalue are calculated. Each eigenvector obtained is rearranged according to the row and column order of the original deviation data matrix to obtain a two-dimensional spatial mode diagram with the same dimension as the original deviation data matrix. This mode is used as the intrinsic orthogonal decomposition mode. The magnitude of the eigenvalues ​​corresponding to this mode characterizes the deviation energy contained in this mode. Among all spatially orthogonal intrinsic orthogonal decomposition modes obtained, the eigenvalues ​​are arranged in descending order. The first few modes whose cumulative energy of the eigenvalues ​​reaches a dominant proportion of the total energy of all eigenvalues ​​are selected as the dominant intrinsic orthogonal decomposition modes. It should be noted that this dominant proportion is determined by analyzing the inflection point of the decay curve of the eigenvalue sequence. Modes after the inflection point are considered to be mainly caused by noise. The selected dominant intrinsic orthogonal decomposition modes and their corresponding modal coefficients are linearly superimposed to reconstruct a dimension-reduced deviation field that reflects the characteristics of tooth surface deviation. This dimension-reduced deviation field is used to compress and reduce noise in the original measurement data.

[0014] Preferably, step S4 involves constructing an optimization objective function and obtaining the optimal processing parameter correction amount through a genetic algorithm based on the dimensionality-reduced bias field and the frequency-related dynamic sensitivity matrix, specifically: By minimizing the arithmetic mean of the absolute values ​​of the normal deviations of the reconstructed dimensionless bias field at all grid points within the region; It should be noted that the bandwidth of the region is taken as the final bearing contact area width of S1 multiplied by a ratio greater than one, in order to ensure coverage of the bearing contact area and its directly adjacent tooth surface area. By minimizing the Euclidean norm of the deviation of the center position of the bearing contact area, characterized by the mechanical characteristic deviation vector of the bearing contact imprint, and the weighted sum of the non-uniformity of the normal contact force distribution inside the bearing contact area; It should be noted that the weighting coefficient is set by taking the ratio of the Euclidean norm of the deviation of the center position of the bearing contact area to the reciprocal of the standard deviation of the normal contact force inside the bearing contact area in the historical convergence data, so as to eliminate the influence of the two due to the difference in the dimensions and scale of change. By minimizing the amplitude of the first harmonic of the load-bearing transmission error curve; It should be noted that the above three objective optimization functions set constraints for the multi-objective optimization problem, including: the absolute value of each machining parameter after correction cannot exceed the upper and lower limits of the physical adjustment stroke of each motion axis of the gear milling machine, which are directly obtained from the technical specifications provided by the machine tool manufacturer; the correction amount of each machining parameter cannot exceed the maximum allowable correction amount calculated from the cutting parameter life curve provided by the tool supplier, under the premise of ensuring a reasonable tool life. This forms the definition of a multi-objective optimization problem containing three objective functions and constraints. The system performs global optimization of the constructed objective function based on a genetic algorithm and a frequency-related dynamic sensitivity matrix, outputting a frontier solution set containing multiple non-dominated solutions. Based on minimizing the ideal point distance, the system calculates the weighted Euclidean distance from each non-dominated solution in the frontier solution set to the ideal point composed of the ideal values ​​of each objective function. The solution with the smallest weighted Euclidean distance is selected as the final set of optimal machining parameter corrections, which are then written into the global variable register of the CNC system through a parametric programming interface.

[0015] Beneficial effects This invention provides a quantitative adjustment method for contact imprints of spiral bevel gears, which has the following beneficial effects: by describing the influence of the unit change of each processing parameter on the transmission of normal contact force and transmission error at any contact node, and by continuously weighting and correcting the participation coefficient of the boundary contact node, the original abrupt and non-differentiable contact boundary movement is transformed into a smooth and continuous process, thus solving the problems of unstable boundary sensitivity and difficulty in convergence of optimization. Furthermore, by solving harmonic balance and numerical extension, the additional normal displacement compensation under dynamic conditions is obtained and superimposed on the initial normal clearance to obtain the updated contact state and transmission error considering dynamic effects. On this basis, the influence of machining parameter perturbation on dynamic contact center offset, dynamic force distribution non-uniformity and transmission error characteristic order harmonic amplitude is quantitatively evaluated. The adjustment method is no longer limited to meeting the static rolling inspection mark qualification, but can actively predict and suppress the dynamic excitation caused by machining error, ensuring that the gear pair can still maintain excellent contact stability and low vibration and noise level under high speed and heavy load conditions. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating a method for quantitative adjustment of contact marks on spiral bevel gears according to the present invention. Detailed Implementation

[0017] The present invention will now be further described with reference to the accompanying drawings.

[0018] Please see Figure 1 This invention provides a method for quantitative adjustment of contact imprints of spiral bevel gears, including the following steps: S1, based on the machining parameters of the large and small gears fixed in the global installation coordinate system, generate discretized grid points on the tooth surface and calculate the initial normal clearance, and construct a set of mechanical equilibrium equations for flexible contact of the tooth surface based on the local contact stiffness of the tooth surface, and obtain the load-bearing contact state distribution dataset and the load-bearing transmission error curve by incremental loading and step-by-step iterative solution; S1 establishes the tooth surface equations of the large and small gears respectively through homogeneous coordinate transformation matrices based on the machining parameters of the large and small gears. The machining parameters include radial tool position, angular tool position, vertical gear position, bed position correction, gear position correction, blank mounting angle, rolling ratio, tool tip radius, and cutter head tooth profile angle. The homogeneous coordinate transformation matrix unifies the tooth surfaces of the large and small gears from spatial coordinates to a fixed global mounting coordinate system. This spatial coordinate transformation process strictly follows the design axis intersection angle and design offset distance of the gear pair. The design axis intersection angle refers to the design angle between the axis of the large gear and the axis of the small gear in space, and the design offset distance refers to the shortest vertical distance between the axis of the large gear and the axis of the small gear. The mesh is divided along the tooth length and tooth height directions under a fixed global installation coordinate system. Specifically, in the tooth tip and tooth root edge regions where the absolute value of the tooth surface curvature gradient is greater than the overall median of the tooth surface curvature gradient, the mesh spacing is set to a small proportion of the tooth surface working height. In the gentle middle region of the tooth surface where the absolute value of the tooth surface curvature gradient is less than and / or equal to the overall median of the tooth surface curvature gradient, the mesh spacing is set to a large proportion of the tooth surface working height, so that the mesh density in the edge region is greater than the mesh density in the middle region. The generated one-to-one correspondence of the mesh points on the large gear tooth surface and the mesh points on the small gear tooth surface, when the large gear and the small gear are at the theoretical meshing angle position, the mesh points in any row and any column on the large gear tooth surface and the mesh points in the same row and the same column on the small gear tooth surface are most close in spatial position, forming a potential contact node pair; For each potential contact node pair generated, calculate the coordinate difference between the grid point on the large gear tooth surface and the grid point on the small gear tooth surface along the normal direction of that point on the large gear tooth surface. Use the absolute value of the coordinate difference as the initial normal clearance of the potential contact node pair. Record all potential contact node pairs and their initial normal clearances as a set of full-tooth-surface potential contact node pairs and their initial normal clearances.

[0019] S1 uses the normal contact force required to produce a unit normal elastic deformation per unit area for a pair of contact nodes on the large and small wheels as the local contact stiffness of the tooth surface. The value of the local contact stiffness of the tooth surface is calculated in advance for each node pair by Hertz contact theory, combined with the elastic modulus and Poisson's ratio of the gear material measured by standard tensile tests. Based on the full tooth surface potential contact node pairs and the initial normal clearance set, the criteria for mutual contact and the definition of mechanical relationship for each potential contact node pair are defined. Specifically, when the sum of the relative elastic deformation of the mesh points of the large wheel and the small wheel under the action of external load exceeds the initial normal clearance, it is determined that they have entered the contact state. The difference between the normal contact force and the sum of the relative elastic deformation of the large wheel and the small wheel in the contact state minus the initial normal clearance is expressed by the nonlinear increasing function relationship determined by Hertz contact theory. The nonlinear increasing function relationships of all potential contact node pairs are combined and added to the macroscopic torque balance equation of the gear pair. The macroscopic torque balance equation of the gear pair means that the sum of the torques of the normal contact forces of all contact node pairs on the gear rotation axis is equal to the input torque, forming a global nonlinear mechanical balance equation system. The unknowns of this equation system are the normal contact forces borne by each potential contact node pair and the macroscopic load transmission error of the gear pair caused by the elastic deformation of the tooth surface. The load transmission error is defined as the deviation of the actual output angle of the large gear from the theoretical output angle of the large gear under the condition of constant input torque of the small gear. The numerical method based on incremental loading and step-by-step iteration solves the global nonlinear mechanical equilibrium equations. The total working torque is divided into several load increment steps starting from zero. The division rule is as follows: in the light load range where the contact state changes gently, a preset constant large step size with the total working torque as the reference is used for division; in the heavy load range where the contact state changes drastically, the step size is divided by successively decreasing. The decreasing factor makes the subsequent step size decrease in a geometric series. In each load increment step, the initial set of contact node pairs for the current step is obtained first based on the bearing contact state obtained by convergence at the end of the previous load increment step. The global nonlinear mechanical equilibrium equations are expanded using a first-order Taylor series to obtain a linearized compliance equation. This compliance equation is then solved using a sparse matrix direct solver based on multi-wavefront LU decomposition, yielding trial solutions for the normal contact force and load transmission error of each node in the current step. Next, the contact node pair set is checked and updated based on these trial solutions: Each potential contact node pair in the set is checked to ensure that the normal contact force is greater than zero; node pairs with normal contact forces less than or equal to zero are removed from the set. Simultaneously, it is checked whether there are potential contact node pairs outside the set whose total relative elastic deformation exceeds their initial normal gap; if so, they are added to the set. This prediction-checking-update iteration is repeated until the contact node pair set in the current step no longer changes, and the relationship between the normal contact force and the total relative elastic deformation of all node pairs entering the contact state in this step satisfies a preset convergence tolerance, which is a minimal proportion of the overall average normal contact force of all contact node pairs in the current step. At this point, the current increment step converges, thus obtaining the load contact force distribution and load transmission error value under that increment step. The calculations for all increment steps are completed sequentially, and the load contact force distributions obtained under each increment step are merged and recorded, ultimately yielding the load contact state distribution dataset under full working torque (i.e., when the last load increment step is completed).

[0020] This dataset contains the spatial coordinates of all discrete contact node pairs on the tooth surface determined to be in contact at the final load increment step, the numerical values ​​of the bearing normal contact force they each bear, and the coordinate information of the spatial boundary of the bearing contact area defined by these spatial coordinates. Simultaneously, the bearing transmission error values ​​obtained at each increment step are connected to form a curve, resulting in the bearing transmission error curve during the full working torque loading process. The contact state at the completion of the final load increment step is the final state mentioned in subsequent steps.

[0021] S2. Based on the final state in the load-bearing contact state distribution dataset as the reference point for linearization expansion, a first-order Taylor expansion of the nonlinear mechanical relationship is performed to obtain a linear relationship. The influence of machining parameter perturbation on the geometric boundary is introduced, and a mapping model is formed by combining global equations. Based on the preset transition zone threshold, the critical potential contact nodes identified in the load-bearing contact state distribution dataset are weighted and corrected to obtain the corrected tooth surface contact stiffness topology mapping model. S2 uses the final state obtained by converging the obtained load contact state distribution data at the final load increment step as the reference point for linearization expansion. In this final state, the normal contact force, the sum of relative elastic deformation, and the load transmission error of all potential contact node pairs in contact are known and definite values. For each potential contact node pair, the nonlinear function relationship between the associated normal contact force, the sum of relative elastic deformation, and the initial normal gap, defined in sub-step 1.3, is expanded at the reference point of this final state using a first-order Taylor expansion. Higher-order infinitesimal terms are discarded to obtain a linear relationship. This linear relationship is expressed as follows: the small change in the normal contact force of a node pair is equal to the local contact stiffness of the node pair multiplied by the small change in the sum of relative elastic deformation of the node pair, minus the local contact stiffness multiplied by the small change in the initial normal gap of the node pair, where the local contact stiffness refers to the derivative of the nonlinear function at the reference point. Then, by establishing the homogeneous coordinate transformation matrix used in the gear geometric model, the partial derivatives of the initial normal clearance of each potential contact node pair with respect to each machining parameter are obtained, that is, the sensitivity coefficient of the small change in the initial normal clearance caused by the small change in the machining parameter. Substituting this sensitivity relationship into the above linear relationship, eliminating the small change in the initial normal gap, the small change in the normal contact force of the node is directly expressed as a linear combination of the small perturbations of each processing parameter. The combination coefficient is the sensitivity of the node to the normal contact force of each processing parameter. Finally, by combining the above linear relationships of all potential contact node pairs with the linearized form of the macroscopic torque balance equation of the gear pair, a global linear equation system is formed. The input vector of this equation system is the small perturbation of all machining parameters, and the output vector is the small change in the normal contact force of all potential contact node pairs and the small change in the load transmission error. The coefficient matrix of this global linear equation system constitutes the tooth surface contact stiffness topological mapping model. Each element of this coefficient matrix quantitatively describes the unit perturbation of a specific machining parameter. By changing the initial normal clearance of each node pair on the entire tooth surface, it causes a change in the normal contact force or macroscopic load transmission error of another specific node pair, thus characterizing the topological relationship between machining parameters, tooth surface micromorphology, and load contact force distribution.

[0022] This application introduces a new continuous state description variable, named the contact participation coefficient, for each pair of critical potential contact nodes that are near the spatial boundary of the bearing contact area in the bearing contact state distribution dataset and that are not in contact in the current state but are adjacent to the boundary. To determine which node pairs belong to the critical potential contact node pairs, a transition zone threshold needs to be determined; the method for determining the transition zone threshold is as follows: Extract all contact node pairs on the spatial boundary of the final bearing contact area in step 1, calculate the ratio of the normal contact force of each node pair to the local contact stiffness of the node pair, and form a statistical distribution of these ratios; take the lower quantile of the statistical distribution as the lower limit of the transition zone threshold, and take the upper quantile as the upper limit of the transition zone threshold. The lower and upper quantiles are selected based on the decay characteristics of the probability density function of the ratio distribution at the tail, so as to cover the main area of ​​the contact state transition near the boundary. For each pair of critical nodes falling within the range defined by the threshold of this transition zone, the rules for determining the value of its contact participation coefficient are designed: the logistic function is used as the smooth step function, and its steepness parameter is determined by the contact imprint rolling inspection calibration test of gears of the same model and batch. The optimization objective is to minimize the overall mean of the sum of squares of the difference between the contact zone boundary movement predicted by the mapping model and the measured movement. The optimal steepness parameter is determined by a single variable scanning optimization method and used as a preset fixed parameter for subsequent calculations. Based on the logistic function, the contact participation coefficient continuously ranges from 0 to 1: when the sum of the relative elastic deformations of the critical node pair is less than its initial normal clearance, and the difference between the two is greater than the upper limit of the transition zone threshold, the contact participation coefficient is 0; when the sum of the relative elastic deformations of the node pair is greater than its initial normal clearance, and the difference between the two is greater than the upper limit of the transition zone threshold, the contact participation coefficient is 1; when the absolute value of the difference between the two is within the transition zone threshold, the contact participation coefficient smoothly transitions from zero to 1 according to the logistic function. Subsequently, this contact participation coefficient is used as a weight, multiplied by the local contact stiffness of the corresponding critical node pair, and the diagonal sub-blocks related to the contact stiffness in the coefficient matrix of the global linear equation system of the tooth surface contact stiffness topology mapping model are weighted and corrected to obtain the tooth surface contact stiffness topology mapping model weighted by the contact participation coefficient. This correction makes the mechanical response near the contact boundary change from a non-differentiable abrupt change to a differentiable smooth change.

[0023] S3. The changes in each processing parameter are transformed into periodic geometric displacement excitation functions related to the meshing frequency and harmonics. The numerical extension solution is performed based on the nonlinear time-domain dynamic equation. The obtained additional normal displacement compensation is superimposed point by point onto the initial normal clearance of all potential contact node pairs used in the load contact state distribution dataset. Based on the modified tooth surface contact stiffness topology mapping model, the updated load contact state and transmission error are obtained, and the frequency-related dynamic sensitivity matrix is ​​constructed. Furthermore, this application transforms the minute changes in each processing parameter into a set of periodic geometric displacement excitation functions applied to the tooth surface, which are related to the meshing frequency and its integer multiples, based on the multiple relationship between the corresponding gear rotation frequency and the meshing frequency. The steady-state vibration response of the gear system under this periodic excitation can be expanded into a Fourier series. This series includes a constant term, a fundamental term with a frequency equal to the meshing frequency, and harmonic terms with frequencies that are integer multiples of the meshing frequency. The Fourier coefficients of each term, including amplitude and phase, are unknown.

[0024] The highest expansion order is determined through trial calculations, starting from the fundamental wave and gradually increasing the order until the amplitude of the newly added harmonics decays to a very small proportion of the fundamental wave amplitude. The vibration response expression containing these unknown Fourier coefficients is substituted into the nonlinear time-domain dynamic equations of the system. Using the harmonic balance principle, it is required that after performing a Fourier transform on both sides of the equations, the harmonic component coefficients of the same frequency are equal. This transforms the time-domain differential equations into a set of coupled nonlinear algebraic equations, the unknowns of which are the amplitudes and phases of the harmonic components of the vibration response. A numerical extension algorithm is used to solve this set of nonlinear algebraic equations under conditions of no machining parameters. The steady-state solution is used as the starting point for solving the problem, and is assigned to the system of equations as the initial solution. Then, the amplitude of the perturbation of the processing parameter is used as the extension parameter. Starting from zero, the parameter is gradually increased by a fixed step size according to a preset proportion of the maximum allowable perturbation of the processing parameter. At each increment step of the extension parameter, the solution obtained by the convergence of the previous parameter step is used as the initial value of the Newton iteration method under the current parameter step. By repeatedly solving the Jacobian matrix and updating the variables, the method converges when the relative change of the Euclidean norm of the solution vectors of two adjacent iterations is less than a preset minimum positive threshold. The amplitude and phase of each harmonic of the vibration response under the current extension parameter value are solved in turn.

[0025] Subsequently, the dynamic transmission error fluctuation amplitude obtained from the solution and related to the meshing frequency is converted into periodic displacement fluctuation along the meshing line direction, and this periodic displacement fluctuation is used as an additional normal displacement compensation amount and is added point by point to the initial normal clearance of all potential contact node pairs used in the load contact state calculation in S1. Substitute this updated clearance into the tooth surface contact stiffness topology mapping model that has been corrected by the contact participation coefficient, and re-execute the incremental loading and step-by-step iterative calculation process to obtain the updated load-bearing contact state distribution dataset and the updated load-bearing transmission error curve after considering dynamic effects. Apply a unit perturbation to each processing parameter and repeat the above process. By comparing the coordinates of the center position of the bearing contact area before and after the perturbation, the non-uniformity of the normal contact force distribution inside the bearing contact area, and the change in the amplitude of a specific order harmonic of the bearing transmission error curve, a two-dimensional numerical matrix is ​​constructed, namely the frequency-related dynamic sensitivity matrix. Each row of the matrix corresponds to a processing parameter, and each column corresponds to a load contact dynamic response evaluation index under the current working speed and working torque. The matrix elements are the quantitative changes in the evaluation index caused by the perturbation of a unit processing parameter.

[0026] S4. Quantitative feature extraction is performed based on the fluorescent contact imprint image set to obtain the mechanical feature deviation vector of the bearing contact imprint. Then, the tooth surface morphology is measured and interpolated to form the tooth surface normal deviation data matrix. Then, the tooth surface normal deviation data matrix is ​​decomposed into intrinsic orthogonal modes. The obtained dominant intrinsic orthogonal decomposition modes and mode coefficients are linearly superimposed to reconstruct the reconstructed dimensionless deviation field. Based on the dimensionless deviation field and the frequency-related dynamic sensitivity matrix, the optimal machining parameter correction amount is obtained through a genetic algorithm. The optimal machining parameter correction amount is written into the global variable register of the CNC system through the parametric programming interface.

[0027] S4 uses a rolling inspection machine to uniformly spray dry powder fluorescent detection agent onto the tooth surfaces of large and small gears. The dry powder fluorescent detection agent is calibrated by a laser thickness gauge to the median value of the recommended optimal development thickness range for the fluorescent detection agent. The same preset total working torque used in S1 is applied to the gear pair for rolling inspection. The number of rotations is set to the minimum number of rotations required to ensure sufficient stability of the contact imprint. The tooth surface is excited by a narrow-band filtered ultraviolet light source with the center wavelength in the ultraviolet band. A high-resolution monochrome industrial camera is used in a dark room environment to acquire several images from multiple angles within a preset symmetrical angle range with the tooth surface normal. After mean superposition and noise reduction, a fluorescent contact imprint image set is formed. For each image in the fluorescent contact imprint image set, an adaptive local thresholding method is used for binarization. The neighborhood window size is taken as a proportion of the length of the minor axis of the theoretical equivalent ellipse of the contact imprint predicted by S1. The pixel threshold is taken as a multiple of the Gaussian weighted gray average value in the neighborhood minus the standard deviation of the gray values ​​in the neighborhood. An active contour model based on gradient vector flow is applied. The initial contour is set as the convex hull of the gray-level abrupt change region in the binary image. When the contour convergence criterion is satisfied, it converges to the sub-pixel accuracy boundary. The contour convergence criterion is satisfied based on the overall change rate of the area enclosed by the contour in the previous two iterations. The area, equivalent ellipse major axis, minor axis, center coordinates and tilt angle are calculated based on the sub-pixel accuracy boundary. Based on the pre-calibrated fluorescence intensity-normal contact pressure proportional relationship using pressure-sensitive film, the grayscale of pixels inside the sub-pixel precision boundary is converted and integrated point by point to reconstruct the relative normal force distribution map. The relative normal force distribution map is then compared with the ideal uniform force distribution constructed by the linearized reference point to calculate the center point coordinate offset vector and the non-uniformity statistics of the force distribution, thus forming the mechanical characteristic deviation vector of the bearing contact imprint.

[0028] The S4 integrates a line laser profile scanning probe next to the rolling inspection machine, reciprocatingly scanning each tooth groove of the large and small wheels. In the specific implementation process, the scanning line spacing is set to a value that balances accuracy and efficiency based on the optical resolution of the probe and the curvature of the tooth surface. Two-dimensional profile lines are acquired in real time, and all profile lines are registered to a fixed global installation coordinate system through an iterative nearest point algorithm. In each iteration, the nearest corresponding point is searched using a KD tree. The optimal rigid body transformation matrix is ​​solved by the singular value decomposition method to minimize the sum of squared distances between point pairs. The process converges when the average distance change between adjacent iterations is less than the preset proportion of the nominal accuracy of the line laser profile scanning probe. The installation eccentricity error is compensated synchronously. The registered measured point cloud is then best-fitted with the theoretical tooth surface, and the deviation value of each measured point along the tooth surface normal is calculated. This deviation value is then mapped onto the tooth surface discretization grid through Kriging space interpolation to form a tooth surface normal deviation data matrix.

[0029] S4 performs intrinsic orthogonal decomposition on the tooth surface normal deviation data matrix, straightening the tooth surface normal deviation data matrix into a high-dimensional column vector. Based on the deviation column vectors of the current measurement sample and multiple samples accumulated from historical measurements, a deviation sample matrix is ​​constructed. After calculating the covariance matrix of the deviation sample matrix, the dimension of the covariance matrix is ​​equal to the number of tooth surface grid points. Using the implicit restart Arnoldi iterative algorithm suitable for solving eigenvalues ​​of large sparse matrices, all eigenvalues ​​of the covariance matrix and the eigenvector corresponding to each eigenvalue are calculated. Each eigenvector obtained is rearranged according to the row and column order of the original deviation data matrix to obtain a two-dimensional spatial mode diagram with the same dimension as the original deviation data matrix, which is used as the intrinsic orthogonal decomposition mode. The magnitude of the eigenvalue corresponding to this mode characterizes the deviation energy contained in this mode. Among all spatially orthogonal intrinsic orthogonal decomposition modes obtained, the eigenvalues ​​are arranged in descending order. The first few modes whose cumulative energy of the eigenvalues ​​reaches a dominant proportion of the total energy of all eigenvalues ​​are selected as the dominant intrinsic orthogonal decomposition modes. It should be noted that this dominant proportion is determined by analyzing the inflection point of the decay curve of the eigenvalue sequence. Modes after the inflection point are considered to be mainly caused by noise. The selected dominant intrinsic orthogonal decomposition modes and their corresponding modal coefficients are linearly superimposed to reconstruct a dimension-reduced deviation field that reflects the characteristics of tooth surface deviation. This dimension-reduced deviation field is used to compress and reduce noise in the original measurement data.

[0030] S4 constructs an optimization objective function and obtains the optimal processing parameter correction amount through a genetic algorithm based on the dimensionality-reduced bias field and the frequency-related dynamic sensitivity matrix, specifically: By minimizing the arithmetic mean of the absolute values ​​of the normal deviations of the reconstructed dimensionless bias field at all grid points within the region; It should be noted that the bandwidth of the region is taken as the final bearing contact area width of S1 multiplied by a ratio greater than one, in order to ensure coverage of the bearing contact area and its directly adjacent tooth surface area. By minimizing the Euclidean norm of the deviation of the center position of the bearing contact area, characterized by the mechanical characteristic deviation vector of the bearing contact imprint, and the weighted sum of the non-uniformity of the normal contact force distribution inside the bearing contact area; It should be noted that the weighting coefficient is set by taking the ratio of the Euclidean norm of the deviation of the center position of the bearing contact area to the reciprocal of the standard deviation of the normal contact force inside the bearing contact area in the historical convergence data, so as to eliminate the influence of the two due to the difference in the dimensions and scale of change. By minimizing the amplitude of the first harmonic of the load-bearing transmission error curve; It should be noted that the above three objective optimization functions set constraints for the multi-objective optimization problem, including: the absolute value of each machining parameter after correction cannot exceed the upper and lower limits of the physical adjustment stroke of each motion axis of the gear milling machine, which are directly obtained from the technical specifications provided by the machine tool manufacturer; the correction amount of each machining parameter cannot exceed the maximum allowable correction amount calculated from the cutting parameter life curve provided by the tool supplier, under the premise of ensuring a reasonable tool life. This forms the definition of a multi-objective optimization problem containing three objective functions and constraints. The system performs global optimization of the constructed objective function based on a genetic algorithm and a frequency-related dynamic sensitivity matrix, outputting a frontier solution set containing multiple non-dominated solutions. Based on minimizing the ideal point distance, the system calculates the weighted Euclidean distance from each non-dominated solution in the frontier solution set to the ideal point composed of the ideal values ​​of each objective function. The solution with the smallest weighted Euclidean distance is selected as the final set of optimal machining parameter corrections, which are then written into the global variable register of the CNC system through a parametric programming interface.

[0031] In some other specific implementations, a non-dominated sorting genetic algorithm with an elitist strategy is used as the global optimizer. Algorithm parameters, such as crossover distribution index, mutation distribution index, and mutation probability, are set according to the size of the optimization variables and the complexity of the problem, along with recommended values. In each iteration, after generating new offspring individuals with different parameter combinations through crossover and mutation operations, the frequency-related dynamic sensitivity matrix is ​​first used to calculate the product of the processing parameter perturbation vector represented by the offspring individual and the sensitivity matrix using matrix multiplication. This yields the predicted values ​​of the increments in the three optimization objective functions corresponding to that individual. Based on this prediction result, the generated offspring individuals are initially screened. The screening rules are as follows: if an offspring individual is predicted to have no deterioration in any of the three objective function values, i.e., the change is no greater than zero, it is retained; if it is predicted to have at least one objective function value improved, i.e. decreased, and the deterioration of the other objective function values ​​does not exceed an allowable proportion of the optimal objective value of its respective parent, it is also retained; otherwise, the offspring individual is directly eliminated. Only the surviving offspring individuals are evaluated using a precise finite element or digital twin simulation based on the tooth surface contact stiffness topology mapping model to obtain the true objective function value.

Claims

1. A method for quantitatively adjusting the contact imprint of spiral bevel gears, characterized in that: The specific steps include the following: S1. Based on the machining parameters of the large and small wheels fixed in the global installation coordinate system, the discretized grid points of the tooth surface are generated and the initial normal clearance is calculated. Based on the local contact stiffness of the tooth surface, the mechanical equilibrium equations of the flexible contact of the tooth surface are constructed. The load-bearing contact state distribution dataset and the load-bearing transmission error curve are obtained by incremental loading and step-by-step iterative solution. S2. Based on the final state in the load-bearing contact state distribution dataset as the reference point for linearization expansion, a first-order Taylor expansion of the nonlinear mechanical relationship is performed to obtain a linear relationship. The influence of machining parameter perturbation on the geometric boundary is introduced, and a mapping model is formed by combining global equations. Based on the preset transition zone threshold, the critical potential contact nodes identified in the load-bearing contact state distribution dataset are weighted and corrected to obtain the corrected tooth surface contact stiffness topology mapping model. S3. The changes in each processing parameter are transformed into periodic geometric displacement excitation functions related to the meshing frequency and harmonics. The numerical extension solution is performed based on the nonlinear time-domain dynamic equation. The obtained additional normal displacement compensation is superimposed point by point onto the initial normal clearance of all potential contact node pairs used in the load contact state distribution dataset. Based on the modified tooth surface contact stiffness topology mapping model, the updated load contact state and transmission error are obtained, and the frequency-related dynamic sensitivity matrix is ​​constructed. S4. Quantitative feature extraction is performed based on the fluorescent contact imprint image set to obtain the mechanical feature deviation vector of the bearing contact imprint. Then, the tooth surface morphology is measured and interpolated to form the tooth surface normal deviation data matrix. Then, the tooth surface normal deviation data matrix is ​​decomposed into intrinsic orthogonal modes. The obtained dominant intrinsic orthogonal decomposition modes and mode coefficients are linearly superimposed to reconstruct the reconstructed dimensionless deviation field. Based on the dimensionless deviation field and the frequency-related dynamic sensitivity matrix, the optimal machining parameter correction amount is obtained through a genetic algorithm. The optimal machining parameter correction amount is written into the global variable register of the CNC system through the parametric programming interface.

2. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 1, characterized in that: S1 establishes the tooth surface equations of the large wheel and the small wheel respectively through the homogeneous coordinate transformation matrix based on the machining parameters of the large wheel and the small wheel. The homogeneous coordinate transformation matrix unifies the tooth surface of the large wheel and the tooth surface of the small wheel from spatial coordinate transformation to a fixed global installation coordinate system. The mesh is divided along the tooth length and tooth height directions under a fixed global installation coordinate system. Specifically, in the tooth tip and tooth root edge regions where the absolute value of the tooth surface curvature gradient is greater than the overall median of the tooth surface curvature gradient, the mesh spacing is set to a small proportion of the tooth surface working height. In the gentle middle region of the tooth surface where the absolute value of the tooth surface curvature gradient is less than and / or equal to the overall median of the tooth surface curvature gradient, the mesh spacing is set to a large proportion of the tooth surface working height. The generated one-to-one correspondence of the mesh points on the large gear tooth surface and the mesh points on the small gear tooth surface, when the large gear and the small gear are at the theoretical meshing angle position, the mesh points in any row and any column on the large gear tooth surface and the mesh points in the same row and the same column on the small gear tooth surface are most close in spatial position, forming a potential contact node pair; For each potential contact node pair generated, calculate the coordinate difference between the grid point on the large gear tooth surface and the grid point on the small gear tooth surface along the normal direction of that point on the large gear tooth surface. Use the absolute value of the coordinate difference as the initial normal clearance of the potential contact node pair. Record all potential contact node pairs and their initial normal clearances as a set of full-tooth-surface potential contact node pairs and their initial normal clearances.

3. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 1, characterized in that: S1 uses the normal contact force required to produce a unit normal elastic deformation per unit area on a pair of contact nodes on the large and small wheels as the local contact stiffness of the tooth surface. Based on the full tooth surface potential contact node pairs and the initial normal clearance set, the criteria for mutual contact and the definition of mechanical relationship for each potential contact node pair are defined. Specifically, when the sum of the relative elastic deformation of the mesh points of the large wheel and the small wheel under the action of external load exceeds the initial normal clearance, it is determined that they have entered the contact state. The difference between the normal contact force and the sum of the relative elastic deformation of the large wheel and the small wheel in the contact state minus the initial normal clearance is expressed by the nonlinear increasing function relationship determined by Hertz contact theory. The nonlinear increasing function relationships of all potential contact node pairs are combined and added to the macroscopic torque balance equation of the gear pair to form a global nonlinear mechanical balance equation system. The numerical method based on incremental loading and step-by-step iteration solves the global nonlinear mechanical equilibrium equations. The total working torque is divided into several load increment steps starting from zero. Within each load increment step, the initial set of contact node pairs for the current step is obtained based on the bearing contact state obtained at the end of the previous load increment step. (The above is completely copied.) The global nonlinear mechanical equilibrium equations are expanded using a first-order Taylor series under the initial set of contact node pairs to obtain linearized compliance equations. The compliance equations are then solved using a sparse matrix direct solver based on multi-wavefront LU decomposition to obtain trial solutions for the normal contact force and load transmission error of each contact node in the current step. The initial set of contact node pairs is checked and updated based on the trial solutions. Initial contact node pairs with normal contact forces less than or equal to zero are removed, and new contact node pairs with closed normal gaps are added. All incremental steps are completed sequentially, and the records are merged to obtain the load contact state distribution dataset and load transmission error curve under full working torque.

4. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 1, characterized in that: The S2 method uses the final state obtained by convergence under the final load increment step, based on the concentrated data of the bearing contact state distribution, as the linearization reference point. For each potential contact node pair, the nonlinear functional relationship between the normal contact force, the sum of the relative elastic deformation, and the initial normal clearance is expanded by a first-order Taylor expansion to obtain a linear relationship. Then, based on the homogeneous coordinate transformation matrix, the partial derivatives of the initial normal clearance of each potential contact node pair with respect to each processing parameter are obtained. The obtained partial derivatives are substituted into the linear relationship, and the linearized form of the linear relationship of all potential contact node pairs and the macroscopic torque balance equation are combined to form a global linear equation system, forming a tooth surface contact stiffness topology mapping model.

5. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 4, characterized in that: S2 extracts all potential contact node pairs located on the spatial boundary of the bearing contact zone from the bearing contact state distribution dataset, calculates the statistical distribution of the ratio of normal contact force to local contact stiffness of all potential contact node pairs, and uses the lower quantile and upper quantile of the ratio as the lower limit and upper limit of the transition zone threshold, respectively. For each pair of critical potential contact nodes located at the spatial boundary of the bearing contact zone in the bearing contact state distribution dataset, the logistic function is used as a smooth step function. The objective is to minimize the overall mean of the sum of squares of the differences between the boundary movement predicted by the tooth surface contact stiffness topology mapping model and the measured movement. This is determined by single-variable scanning optimization. When solving, if the sum of the normal relative position and deformation of the critical potential contact node pair is less than the initial normal gap and the difference is greater than the upper limit of the transition zone, the value is zero. If the sum of the normal relative position and deformation of the node pair is greater than the initial normal gap and the difference is greater than the upper limit, the value is one. This yields the contact participation coefficient. Based on the contact participation coefficient as a weight multiplied by the local contact stiffness of the corresponding critical potential contact node pair, the diagonal sub-blocks related to the contact stiffness in the coefficient matrix of the global linear equation system of the tooth surface contact stiffness topology mapping model are weighted and corrected to obtain the tooth surface contact stiffness topology mapping model weighted by the contact participation coefficient.

6. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 5, characterized in that: The S3 process transforms the changes in each processing parameter into a periodic geometric displacement excitation function related to the meshing frequency and harmonics, thereby expanding the steady-state vibration response of the large and small wheels into a Fourier series. The highest expansion order is determined by gradually increasing the order from the fundamental wave until the amplitude of the newly added harmonics decays to the fundamental wave amplitude. The vibration response in the form of the Fourier series is substituted into the nonlinear time-domain dynamic equation. The nonlinear time-domain dynamic equation is transformed into a set of coupled nonlinear algebraic equations using the harmonic balance principle. The equations are then solved based on the numerical extension algorithm to obtain the harmonic amplitude of the vibration response under each extension step. Based on the obtained fluctuation amplitude, it is converted into periodic displacement fluctuation along the meshing line, which is used as an additional normal displacement compensation amount. It is then superimposed point by point onto the initial normal clearance of all potential contact node pairs used in the load contact state distribution dataset. It is substituted into the tooth surface contact stiffness topology mapping model after the contact participation coefficient is weighted and corrected, and the incremental loading and step-by-step iterative solution of S1 are repeated to obtain the updated load contact state and transmission error. After repeating the operation on each processing parameter, the coordinates of the center position of the bearing contact area before and after, the non-uniformity of the normal force distribution, and the change in the amplitude of a specific order harmonic of the bearing transmission error are compared to construct a two-dimensional numerical matrix, which serves as the frequency-related dynamic sensitivity matrix.

7. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 1, characterized in that: The S4 is based on the rolling inspection machine to uniformly spray dry powder fluorescent detection agent on the tooth surfaces of large and small gears, and excites the tooth surface based on a narrow-band filtered ultraviolet light source with the center wavelength in the ultraviolet band. A high-resolution monochrome industrial camera is used to collect several images from multiple angles within a preset symmetrical angle range with the tooth surface normal. After average superposition and noise reduction, a fluorescent contact imprint image set is formed. Based on each image in the fluorescent contact imprint image set, an adaptive local thresholding method is used for binarization. An active contour model based on gradient vector flow is applied. The initial contour is set as the convex hull of the gray-level abrupt change region in the binary image. The contour is iteratively evolved until the contour convergence criterion is satisfied, and then converges to the sub-pixel accuracy boundary. The area, equivalent ellipse major axis, minor axis, center coordinates and tilt angle are calculated based on the sub-pixel accuracy boundary. Based on the pre-calibrated fluorescence intensity-normal contact pressure proportional relationship using pressure-sensitive film, the grayscale of pixels inside the sub-pixel precision boundary is converted and integrated point by point to reconstruct the relative normal force distribution map. The relative normal force distribution map is then compared with the ideal uniform force distribution constructed by the linearized reference point to calculate the center point coordinate offset vector and the non-uniformity statistics of the force distribution, thus forming the mechanical characteristic deviation vector of the bearing contact imprint.

8. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 7, characterized in that: The S4 integrates the line laser profile scanning probe to the side of the rolling inspection machine, and performs reciprocating scanning of each tooth groove of the large and small wheels. It acquires two-dimensional profile lines in real time, and registers all profile lines to a fixed global installation coordinate system through an iterative nearest point algorithm. After solving the optimal rigid body transformation matrix by the singular value decomposition method to minimize the sum of squared distances between point pairs, it simultaneously compensates for installation eccentricity errors. Then, it matches the registered measured point cloud with the theoretical tooth surface for best fit, calculates the deviation value of each measured point along the tooth surface normal, and maps it to the tooth surface discretization grid through Kriging space interpolation to form a tooth surface normal deviation data matrix.

9. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 8, characterized in that: S4 performs intrinsic orthogonal decomposition on the tooth surface normal deviation data matrix, straightens the tooth surface normal deviation data matrix into a high-dimensional column vector, and constructs a deviation sample matrix based on the deviation column vectors of the current measurement sample and multiple samples accumulated from historical measurements. After calculating the covariance matrix of the deviation sample matrix, it calculates all eigenvalues ​​of the covariance matrix and the eigenvector corresponding to each eigenvalue. Each eigenvector obtained is rearranged according to the row and column order of the original deviation data matrix to obtain a two-dimensional spatial mode diagram with the same dimension as the original deviation data matrix, which is used as the intrinsic orthogonal decomposition mode. Among all the spatially orthogonal eigenorthogonal decomposition modes obtained, the eigenvalues ​​are arranged in descending order, and the top few modes whose cumulative energy of the eigenvalues ​​reaches a dominant proportion of the sum of all eigenvalues ​​are selected as the dominant eigenorthogonal decomposition modes. The selected dominant intrinsic orthogonal decomposition modes and their corresponding modal coefficients are linearly superimposed to reconstruct a dimensionless deviation field that reflects the characteristics of tooth surface deviation, which is then used as the reconstructed dimensionless deviation field.

10. The method for quantitative adjustment of contact imprints of spiral bevel gears according to claim 9, characterized in that: The S4 method constructs an optimization objective function and obtains the optimal processing parameter correction amount through a genetic algorithm based on the dimensionality-reduced bias field and the frequency-related dynamic sensitivity matrix. Specifically: By minimizing the arithmetic mean of the absolute values ​​of the normal deviations of the reconstructed dimensionless bias field at all grid points within the region; By minimizing the Euclidean norm of the deviation of the center position of the bearing contact area, characterized by the mechanical characteristic deviation vector of the bearing contact imprint, and the weighted sum of the non-uniformity of the normal contact force distribution inside the bearing contact area; By minimizing the amplitude of the first harmonic of the load-bearing transmission error curve; The system performs global optimization of the constructed objective function based on a genetic algorithm and a frequency-related dynamic sensitivity matrix, outputting a frontier solution set containing multiple non-dominated solutions. Based on minimizing the ideal point distance, the system calculates the weighted Euclidean distance from each non-dominated solution in the frontier solution set to the ideal point composed of the ideal values ​​of each objective function. The solution with the smallest weighted Euclidean distance is selected as the final set of optimal machining parameter corrections, which are then written into the global variable register of the CNC system through a parametric programming interface.

Images